text:
"""
    ReifiedInSet(x::AbstractIntVar, s::IntSetVar, b::BoolVar, trailer::Trailer)

ReifiedInSet contrainst, states that b ⟺ x ∈ s
"""
struct ReifiedInSet <: Constraint
    x       ::AbstractIntVar
    s       ::IntSetVar
    b       ::AbstractBoolVar
    active  ::StateObject{Bool}
    function ReifiedInSet(x::AbstractIntVar, s::IntSetVar, b::AbstractBoolVar, trailer::Trailer)
        constraint = new(x, s, b, StateObject{Bool}(true, trailer))
        addOnDomainChange!(x, constraint)
        addOnDomainChange!(s, constraint)
        addOnDomainChange!(b, constraint)
        return constraint
    end
end

"""
    propagate!(constraint::ReifiedInSet, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`ReifiedInSet` propagation function.
"""
function propagate!(constraint::ReifiedInSet, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    surely_in_set = all(is_required(constraint.s.domain, v) for v in constraint.x.domain)
    surely_notin_set = all(!is_possible(constraint.s.domain, v) for v in constraint.x.domain)

    if isbound(constraint.b)
        # Feasibility
        if assignedValue(constraint.b) && surely_notin_set
            return false
        end
        if !assignedValue(constraint.b) && surely_in_set
            return false
        end

        # Filtering s & x
        if assignedValue(constraint.b) # Make sure x ∈ s
            # Filter x
            x_values = collect(constraint.x.domain)
            for v in x_values
                if !is_possible(constraint.s.domain, v)
                    remove!(constraint.x.domain, v)
                    addToPrunedDomains!(prunedDomains, constraint.x, [v])
                    triggerDomainChange!(toPropagate, constraint.x)
                end
            end

            # Filter s
            if isbound(constraint.x)
                if !is_required(constraint.s.domain, assignedValue(constraint.x))
                    require!(constraint.s.domain, assignedValue(constraint.x))
                    addToPrunedDomains!(prunedDomains, constraint.s, SetModification(;excluded=[assignedValue(constraint.x)]))
                    triggerDomainChange!(toPropagate, constraint.s)
                end
            end
        else # Make sure x ∉ s
            # Filter x
            x_values = collect(constraint.x.domain)
            for v in x_values
                if is_required(constraint.s.domain, v)
                    remove!(constraint.x.domain, v)
                    addToPrunedDomains!(prunedDomains, constraint.x, [v])
                    triggerDomainChange!(toPropagate, constraint.x)
                end
            end

            # Filter s
            if isbound(constraint.x)
                if is_possible(constraint.s.domain, assignedValue(constraint.x))
                    exclude!(constraint.s.domain, assignedValue(constraint.x))
                    addToPrunedDomains!(prunedDomains, constraint.s, SetModification(;excluded=[assignedValue(constraint.x)]))
                    triggerDomainChange!(toPropagate, constraint.s)
                end
            end
        end
    else
        # No feasibility test because if b not assigned, it is always feasible

        # Filtering b
        if surely_in_set
            removed = assign!(constraint.b, true)
            addToPrunedDomains!(prunedDomains, constraint.b, removed)
            triggerDomainChange!(toPropagate, constraint.b)
        end
        if surely_notin_set
            removed = assign!(constraint.b, false)
            addToPrunedDomains!(prunedDomains, constraint.b, removed)
            triggerDomainChange!(toPropagate, constraint.b)
        end
    end

    surely_in_set = all(is_required(constraint.s.domain, v) for v in constraint.x.domain)
    surely_notin_set = all(!is_possible(constraint.s.domain, v) for v in constraint.x.domain)

    # Deactivation
    if isbound(constraint.b) && ((assignedValue(constraint.b) && surely_in_set) || (!assignedValue(constraint.b) && surely_notin_set))
        setValue!(constraint.active, false)
    end

    return true
end

variablesArray(constraint::ReifiedInSet) = [constraint.x, constraint.s, constraint.b]

function Base.show(io::IO, ::MIME"text/plain", con::ReifiedInSet)
    println(io, typeof(con), ": ", con.b.id, " ≡ ", con.x.id, " ∈ ", con.s.id, ", active = ", con.active)
    println(io, "   ", con.b)
    println(io, "   ", con.x)
    print(io, "   ", con.s)
end

function Base.show(io::IO, con::ReifiedInSet)
    print(io, typeof(con), ": ", con.b.id, " ≡ ", con.x.id, " ∈ ", con.s.id)
end


text:
using Clang.Generators
using Clang.LibClang.Clang_jll

cd(@__DIR__)

options = load_options(joinpath(@__DIR__, "generator.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()

headers = ["./rtmidi/rtmidi_c.h"]
# there is also an experimental `detect_headers` function for auto-detecting top-level headers in the directory
# headers = detect_headers(clang_dir, args)

# create context
ctx = create_context(headers, String[], options)

# run generator
build!(ctx)

ext: jl


text:
#=
Created on Tuesday 11 August 2020
Last update: Tuesday 18 August 2020

@author: Michiel Stock
michielfmstock@gmail.com

Simple demonstration of linear regression with zero-norm minimization.
=#

using BioLizard.ZeroNorm
using Plots, StatsBase

# generate some data

n, p, σ = 150, 50, 1

X = randn(n, p)
strue = rand(p) .< 0.2  # use 20% of the features
wtrue = randn(p)
y = X * (wtrue .* strue) .+ σ * randn(n)

# simple squared loss, easy as pie
l(w) = mean((X * w .- y).^2)

# initial guess for w
w = zeros(p)
λ = 0.05

# use the method to find the norm
s, w = findzeronorm(w, l, λ, SimulatedAnnealing(), verbose=true)

# can also use an initial s
#s0 = ones(Bool, p)
#s, obj_vals = findzeronorm!(w, l, λ, s0, verbose=true)

# get confusion table
tp = sum(strue[s])  # true pos
fp = sum(.!strue[s])
tn = sum(.!strue[.!s])
fn = sum(strue[.!s])

println("""
TP = $tp
FP = $fp
TN = $tn
FN = $fn
""")

scatter(w[strue], wtrue[strue], color=1s[strue])
xlabel!("w fitted")
ylabel!("w true")


text:
@testset "spatial" begin

@testset "no units, no time" begin
    A = NamedAxisArray(reshape(1:12, 3, 4), x = 1:3, y = 1:4);
    @test_throws ArgumentError time_axis(A)
    @test !has_timedim(A)
    @test_throws ArgumentError timedim(A)
    @test ntime(A) == 1
    @test @inferred(pixel_spacing(A)) === (1,1)
    @test @inferred(spatial_directions(A)) === ((1,0),(0,1))
    @test @inferred(spatialdims(A)) === (1, 2)
    @test @inferred(spatial_order(A)) === (:x, :y)
    @test @inferred(spatial_size(A)) === (3,4)
    @test @inferred(spatial_indices(A)) === (Base.OneTo(3), Base.OneTo(4))
    assert_timedim_last(A)
end

@testset "units, no time" begin
    mm = u"mm"     # in real use these should be global consts
    m = u"m"
    A = NamedAxisArray(reshape(1:12, 3, 4), x = 1mm:1mm:3mm, y = 1m:2m:7m)
    @test_throws ArgumentError time_axis(A)
    @test !has_timedim(A)
    @test_throws ArgumentError timedim(A)
    @test ntime(A) == 1
    @test @inferred(pixel_spacing(A)) === (1mm,2m)
    @test @inferred(spatial_directions(A)) === ((1mm,0m),(0mm,2m))
    @test @inferred(spatialdims(A)) === (1,2)
    @test @inferred(spatial_order(A)) === (:x,:y)
    @test @inferred(spatial_size(A)) === (3,4)
    @test @inferred(spatial_indices(A)) === (Base.OneTo(3),Base.OneTo(4))
    assert_timedim_last(A)
end

@testset "units, time" begin
    s = u"s" # again, global const
    A = NamedAxisArray(reshape(1:12, 3, 4), x = 1:3, time = 1s:1s:4s)
    @test @inferred(time_keys(A)) == 1s:1s:4s
    @test has_timedim(A)
    @test timedim(A) == 2
    @test ntime(A) == 4
    @test @inferred(pixel_spacing(A)) === (1,)
    @test @inferred(spatial_directions(A)) === ((1,),)
    @test @inferred(coords_spatial(A)) === (1,)
    @test @inferred(spatialorder(A)) === (:x,)
    @test @inferred(size_spatial(A)) === (3,)
    @test @inferred(spatial_indices(A)) === (Base.OneTo(3),)
    assert_timedim_last(A)
end

@testset "units, time first" begin
    s = u"s" # global const
    A = NamedAxisArray(reshape(1:12, 4, 3), time = 1s:1s:4s, x = 1:3)
    @test @inferred(time_keys(A)) === 1s:1s:4s
    @test has_timedim(A)
    @test timedim(A) == 1
    @test ntime(A) == 4
    @test @inferred(pixel_spacing(A)) === (1,)
    @test @inferred(spatial_directions(A)) === ((1,),)
    @test @inferred(coords_spatial(A)) === (2,)
    @test @inferred(spatial_order(A)) === (:x,)
    @test @inferred(spatial_size(A)) === (3,)
    @test @inferred(spatial_indices(A)) === (Base.OneTo(3),)
    @test_throws ErrorException assert_timedim_last(A)
end

end


text:
function acqr_kfilter(y,a,c,q,r,x10,P10)
	#
	# Kalman filter for model
	#
	# x_{t+1} = a*x_{t} + w_{t}
	# y_{t}   = c*x_{t} + v_{t}
	#
	# cov(w_{t},v_{t}) = [q 0;0 r]
	#
	# javier.cara@upm.es, 2020-05
	#

	nt = length(y)

	# allocation
	xtt = zeros(nt)
	Ptt = zeros(nt)
	xtt1 = zeros(nt+1)
	Ptt1 = zeros(nt+1)
	et = zeros(nt)
	St = zeros(nt)
	Kt = zeros(nt)
	loglik = 0.0

	# Filter
	xtt1[1] = x10
	Ptt1[1] = P10
	for t in 1:nt
		#  innovations
		et[t] = y[t] - c*xtt1[t]
		St[t] = c^2*Ptt1[t] + r # et variance

		# Kalman gain
		Kt[t] = Ptt1[t]*c/St[t] # kalman gain

		# filtered values
		xtt[t] = xtt1[t] + Kt[t]*et[t]
		Ptt[t] = (1 - Kt[t]*c)*Ptt1[t]

		# one-step ahead prediction
		xtt1[t+1] = a*xtt[t]
		Ptt1[t+1] = a^2*Ptt[t] + q

		# likelihood
		#l0 = et[t]^2/St[t] # typeof(l0) = Array{Float64,1}
		loglik = loglik + log(St[t]) + et[t]^2/St[t]
	end

	loglik =  - nt/2*log(2*pi) - 0.5*loglik

	return xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik

end


text:
# Useful Functions

"""
    max(vec::Vector{Atom}, para)

Do this will return the maximum data of `para` in `vec_atom`.

# Example
```julia-repl
max_data = max(vec_atom, "coord")
```
"""
function max(vec::Vector{Atom}, para)
    fields = fieldnames(typeof(vec[1]))
    para = Meta.parse(para)
    if !in(para, fields)
        error(join(["Data_Type Atom doesn't have field: ", string(para)]))
    end
    data_para = getfield(vec[1], para)
    len = length(data_para)
    if len != 1
        max = Vector{typeof(data_para).parameters[1]}(undef, len)
        max[:] = data_para[:]
        for i = 2 : length(vec)
            data_now = getfield(vec[i], para)
            for j = 1 : len
                if max[j] <= data_now[j]
                    max[j] = data_now[j]
                end
            end
        end
    else
        max = getfield(vec[1], para)
        for i = 2 : length(vec)
            data_now = getfield(vec[i], para)
            if max <= data_now
                max = data_now
            end
        end
    end
    max
end

function max(array_01::AbstractArray, array_02::AbstractArray)
    if typeof(array_01).parameters[2] == 1
        num_cols = 1
        num_rows = length(array_01)
    else
        num_rows, num_cols = size(array_01)
    end
    [array_01[row, col]>array_02[row, col] ? array_01[row, col] : array_02[row, col] for row = 1:num_rows, col = 1:num_cols]
end
"""
    min(vec::Vector{Atom}, para)

Do this will return the minimum data of `para` in `vec_atom`.

# Example
```julia-repl
in_data = min(vec_atom, "coord")
```
"""
function min(vec::Vector{Atom}, para)
    fields = fieldnames(typeof(vec[1]))
    para = Meta.parse(para)
    if !in(para, fields)
        error(join(["Data_Type Atom doesn't have field: ", string(para)]))
    end
    data_para = getfield(vec[1], para)
    len = length(data_para)
    if len != 1
        min = Vector{typeof(data_para).parameters[1]}(undef, len)
        min[:] = data_para[:]
        for i = 2 : length(vec)
            data_now = getfield(vec[i], para)
            for j = 1 : len
                if min[j] >= data_now[j]
                    min[j] = data_now[j]
                end
            end
        end
    else
        min = getfield(vec[1], para)
        for i = 2 : length(vec)
            data_now = getfield(vec[i], para)
            if min >= data_now
                min = data_now
            end
        end
    end
    min
end

function min(array_01::AbstractArray, array_02::AbstractArray)
    if typeof(array_01).parameters[2] == 1
        num_cols = 1
        num_rows = length(array_01)
    else
        num_rows, num_cols = size(array_01)
    end
    [array_01[row, col]<array_02[row, col] ? array_01[row, col] : array_02[row, col] for row = 1:num_rows, col = 1:num_cols]
end

"""
    get_data(vec_unit::Vector{T}, para::AbstractString) where T <: Data_Unit
Do this will return a vector of variable `para` for all elements in `vec_unit`.
# Example
```julia-repl
str = Si3N4()
cell = genr_cell([10, 10, 10])
data = genr(cell, str)
get_data(data.vec_angle, "atom")
```
"""
function get_data(vec_unit::Vector{T}, para::AbstractString) where T <: Unit
    fields = fieldnames(typeof(vec_unit[1]))
    para = Meta.parse(para)
    if !in(para, fields)
        error(join(["Data_Type Atom doesn't have field: ", string(para)]))
    end
    num_units = length(vec_unit)
    data = getfield(vec_unit[1], para)
    res = zeros(typeof(data[1]), num_units, length(data))

    for unit = 1 : num_units
        res[unit, :] .= getfield(vec_unit[unit], para) # try a=zeros(3, 1) a[1, :] = 1
    end
    return res
end

"""
    add(vec_unit::Vector{T}, tilt, para::AbstractString) where T <: Data_Unit

Do this will add `tilt` to variable `para` for all elements in `vec_unit`

# Example
```julia-repl
str = Si3N4()
cell = genr_cell([10, 10, 10])
data = genr(cell, str)
add!(data.vec_bond, 3, "typ")
```
"""
function add!(vec_unit::Vector{T}, tilt, para::AbstractString) where T <: Unit
    fields = fieldnames(typeof(vec_unit[1]))
    para = Meta.parse(para)
    if !in(para, fields)
        error(join(["Data_Type Atom doesn't have field: ", string(para)]))
    end
    num_units = length(vec_unit)

    for unit = 1 : num_units
        data = getfield(vec_unit[unit], para)
        len = length(data)
        if len == 1
            data += tilt
        else
            data .+= tilt
        end
        setfield!(vec_unit[unit], para, data)
    end
end

"""
    set!(vec_unit::Vector{T}, tilt, para::AbstractString) where T <: Unit

Do this will change variable `para` of all elements in `vec_unit` to `tilt`

# Example
```julia-repl
str = Si3N4()
cell = genr_cell([10, 10, 10])
data = genr(cell, str)
change(data.vec_bond, 3, "typ")
```
"""
function set!(vec_unit::Vector{T}, goal, para::AbstractString) where T <: Unit
    fields = fieldnames(typeof(vec_unit[1]))
    para = Meta.parse(para)
    if !in(para, fields)
        error(join(["Data_Type Atom doesn't have field: ", string(para)]))
    end
    num_units = length(vec_unit)

    for unit = 1 : num_units
        data = getfield(vec_unit[unit], para)
        data = goal
        setfield!(vec_unit[unit], para, data)
    end
end

function set!(vec_unit::Vector{T}, goal::AbstractArray, para::AbstractString) where T <: Unit
    fields = fieldnames(typeof(vec_unit[1]))
    para = Meta.parse(para)
    dims = size(goal)[2]
    if !in(para, fields)
        error(join(["Data_Type Atom doesn't have field: ", string(para)]))
    end
    num_units = length(vec_unit)

    for unit = 1 : num_units
        data = getfield(vec_unit[unit], para)
        if dims == 1
            data = goal[unit]
        else
            data = goal[unit, :]
        end
        setfield!(vec_unit[unit], para, data)
    end
end

function diag(vec::AbstractArray)
    len = length(vec)
    result = zeros(len, len)
    for i = 1 : len
        result[i, i] = vec[i]
    end
    result
end

function conv(vec::Vector, dim)
    len = length(vec)
    para = typeof(vec).parameters
    if dim == 1
        result = convert.(para[1], zeros(1, len))
    elseif dim == 2
        result = convert.(para[1], zeros(len, 1))
    else
        error("dim should be 1 or 2, representing Row or Column vector respectively")
    end
    for i in 1 : len
        result[i] = vec[i]
    end
    result
end

function conv(mat::Matrix, redu=0)
    para = typeof(mat).parameters
    if redu == 1
        (row, col) = size(mat)
        len = row>col ? row : col
        result = convert.(para[1], zeros(len))
    elseif redu == 0
        result = convert.(para[1], zeros(size(mat)[2], size(mat)[1]))
    else
        error("redu should be 0 or 1, representing wether Martix will be reduced to Vector.")
    end
    result[:] = mat[:]
end

"""
    dist(atom::Atom, org::Array)

Do this will return the distance between coord of `atom` and coord of `org`

# Example
```julia-repl
r = dist(atom, [0 0 0])
```
"""
function dist(atom::Atom, org::Array)
    r = 0
    for dim = 1 : 3
        r += (org[dim] - atom.coord[dim])^2
    end
    sqrt(r)
end

"""
    dist(atom::Atom, org::Array, dim)

Do this will return the 2-d distance between coord of `atom` and coord of `org` in plane vertical to axies `dim`.
`dim = 1, 2, 3` represents axies x, y, z respectively.

# Example
```julia-repl
r = dist(atom, [0 0 0], 1)
```
"""
function dist(atom::Atom, org::Array, dim)
    r = 0
    dim_range = [1 2 3]
    dim_range = dim_range[1:end .!= dim]
    for i in dim_range
        r += (org[i] - atom.coord[i])^2
    end
    sqrt(r)
end

"""
    dist(pos::Array, org::Array)

Do this will return the distance between coord of `pos` and coord of `org`

# Example
```julia-repl
r = dist([1 2 3], [0 0 0], 1)
```
"""
function dist(pos::Array, org::Array)
    r = 0
    for dim = 1 : 3
    r += (org[dim] - pos[dim])^2
    end
    sqrt(r)
end

"""
    dist(pos::Array, org::Array, dim)

Do this will return the 2-d distance between coord of `pos` and coord of `org` in plane vertical to axies `dim`.
`dim = 1, 2, 3` represents axies x, y, z respectively.

# Example
```julia-repl
r = dist([1 2 3], [0 0 0], 1)
```
"""
function dist(pos::Array, org::Array, dim)
    r = 0
    dim_range = [1 2 3]
    dim_range = dim_range[1:end .!= dim]
    for i in dim_range
        r += (org[i] - pos[i])^2
    end
    sqrt(r)
end

"""
    norm_vec(vec::Vector)
Do this will return a normalized vector along the direction of `vec`
"""
function norm_vec(vec::Vector)
    return vec ./ norm(vec)
end

"""
    function rot_mat(θx, θy, θz)
This will return a rotation matrix, coord times which will have a rotation along x, y, and z axies with angle θx, θy, θz respectively.

More details in [Wiki](https://zh.wikipedia.org/wiki/%E6%97%8B%E8%BD%AC%E7%9F%A9%E9%98%B5)
"""
function rot_mat(θx, θy, θz)
    x_mat = [
        1  0        0
        0  cos(θx)  -sin(θx)
        0  sin(θx)  cos(θx)
    ]

    y_mat = [
        cos(θy)  0  sin(θy)
        0        1  0
        -sin(θy) 0  cos(θy)
    ]

    z_mat = [
        cos(θz)  -sin(θz) 0
        sin(θz)  cos(θz)  0
        0        0        1
    ]

    return x_mat * y_mat * z_mat
end

"""
    function central_point_atom(data::Data)
This will return a vector of the coordinate of the center points, based on the atom infomation.
**Recommand** to use when model is small.
"""
function central_point_atom(data::Data)
	coord_vec = hcat([data.vec_atom[n].coord for n = 1:data.data_basic.num_atoms]...)
	return mean(coord_vec, dims=2)
end

"""
    function central_point_atom(data::Data)
This will return a vector of the coordinate of the center points, based on the box infomation.
**Recommand** to use when model is large.
"""
function central_point_box(data::Data)
	box_size = data.data_basic.box_size
	return (box_size[:, 2] .- box_size[:, 1])/2
end

"""
    genr_box_diag(data::Data)
This will return the matrix of box size information.
"""
function genr_box_diag(data::Data)
	box_info = data.data_basic.box_size
	box_info = box_info[:, 2] - box_info[:, 1]
	box_diag = diag(box_info)
	try
		data.data_basic.box_tilt
	catch
		return box_diag
	end
	box_diag[2, 1] = data.data_basic.box_tilt[1]
	box_diag[3, 1] = data.data_basic.box_tilt[2]
	box_diag[3, 2] = data.data_basic.box_tilt[3]
	return box_diag
end

"""
    genr_box_inv(data::Data)
This will return the inverse matrix of box size information.
"""
function genr_box_inv(data::Data)
	box_diag = genr_box_diag(data)
	return inv(box_diag)
end

function copy_array(goal)
    res = zeros(typeof(goal).parameters[1], size(goal))
    res[:] = goal[:]
    return res
end

"""
    function grid(x, y)
This will generate two array from two vector for contour.
"""
function grid(x, y)
    len_x = length(x)
    len_y = length(y)

    x_mat = zeros(len_y, len_x)
    y_mat = zeros(len_y, len_x)
    z_mat = zeros(len_y, len_x)

    for i = 1:len_y
        for j = 1:len_x
            x_mat[i, j] = x[j]
            y_mat[i, j] = y[i]
        end
    end

    return x_mat, y_mat
end

function in_str(data_str::Str, vec_str::Vector{T}) where T <: Str
    for str in vec_str
        if str.atom_name == data_str.atom_name
            return true
        end
    end
    return false
end


text:
"Mixture of Factor Analyzers"
function fit_mm(::Type{FactorAnalysis}, X::AbstractMatrix{T},
                k::Integer;            # the number of factor analyzers to use
                factors::Integer = 1,  # the number of factors in each analyzer
                tol::Real=1.0e-6,      # convergence tolerance
                maxiter::Integer=1000, # number of iterations
                μs::Union{AbstractMatrix{T}, Nothing} = nothing,
                Σs::Union{AbstractArray{T,3}, Nothing} = nothing
            ) where {T<:AbstractFloat}

    d, n = size(X)
    CMVN = (2π)^(-d/2)
    mX=mean(X, dims=2)
    cX=cov(X, dims=2)
    sc=det(cX)^(1/d)
    l = factors

    # initialize parameters
    πⱼ = fill(one(T)/k, k)
    Wⱼ = zeros(T,d,l,k)
    μⱼ = zeros(T,d,k)
    Ψⱼ = zeros(T,d,k)
    rn = 1:min(d,l)
    for j in 1:k
        if Σs !== nothing
            ev = eigvecs(Σs[:,:,j], sortby=-)
            Wⱼ[:,rn,j] = ev[:,rn].*sqrt(det(Σs[:,:,j])^(1/d))
        else
            Wⱼ[:,:,j] = randn(d,l) * sqrt(sc/l)
        end
        μⱼ[:,j] = μs !== nothing ? μs[:,j] : (vec(randn(d)' * sqrt(cX)) .+ mX)
        Ψⱼ[:,j] = diag(Σs !== nothing ? Σs[:,:,j] : cX) .+ eps()
    end

    hᵢⱼ = zeros(T,k,n)
    πᵢⱼ = similar(hᵢⱼ)
    E₍zzᵀₗxᵢωⱼ₎ = zeros(T,l,l,k)
    Σ⁻¹ⱼ = zeros(T,d,d,k)
    ZIᵀ=zeros(T,n,l+1,k)
    ZIᵀ[:,end,:] .= one(T)
    Σ⁻¹Y = similar(X)
    Y = similar(X)

    E₍zₗxᵢωⱼ₎ = view(ZIᵀ,:,1:l,:)

    L_old = 0.0
    for itr in 1:maxiter
        # E Step
        # 1) Calculate hᵢⱼ, E₍zₗxᵢωⱼ₎
        for j in 1:k
            Ψ⁻¹ = Diagonal(1 ./ Ψⱼ[:, j])
            W = @view Wⱼ[:,:,j]
            Σ⁻¹ = @view Σ⁻¹ⱼ[:,:,j]
            # V⁻¹ = inv(I + WᵀΨ⁻¹*W)
            WᵀΨ⁻¹ = W'*Ψ⁻¹
            # Σ⁻¹ = Ψ⁻¹ - Ψ⁻¹*W*V⁻¹*WᵀΨ⁻¹
            Σ⁻¹[:] = Ψ⁻¹ - Ψ⁻¹*W*inv(I + WᵀΨ⁻¹*W)*WᵀΨ⁻¹
            detΣ⁻¹ = sqrt(det(Σ⁻¹))
            # Y = X .- view(μⱼ,:,j)
            broadcast!(-,Y,X,view(μⱼ,:,j))
            Σ⁻¹Y .= Σ⁻¹*Y
            # view(hᵢⱼ,j,:) .= πⱼ[j]*(CMVN*detΣ⁻¹).*exp.(-0.5.*vec(sum(Y.*Σ⁻¹Y, dims=1)))
            broadcast!(*,Y,Y,Σ⁻¹Y)
            sum!(view(hᵢⱼ,j,:)', Y)
            c = πⱼ[j]*(CMVN*detΣ⁻¹)
            broadcast!(e->c*exp(-0.5*e),view(hᵢⱼ,j,:),view(hᵢⱼ,j,:))
            # Ez[:,:,j] = V⁻¹*WᵀΨ⁻¹*Y # W'*Σ⁻¹*Y
            view(E₍zₗxᵢωⱼ₎,:,:,j) .= (W'Σ⁻¹Y)'
        end

        # 2) Calculate log(L)
        hᵢⱼ[hᵢⱼ .< eps()] .= eps()
        hᵢ = sum(hᵢⱼ, dims=1)
        L = sum(log.(hᵢ))

        # 3) Calculate πⱼ
        # πᵢⱼ = hᵢⱼ./hᵢ
        broadcast!(/,πᵢⱼ,hᵢⱼ,hᵢ)
        πⱼ .= sum(πᵢⱼ, dims=2) |> vec

        # Check exit conditions
        if itr > 1
            chg = L - L_old
            if chg < 0
                @warn "log likelihood decreased" itr=itr change=chg
            end
            @debug "log(L)=$L" itr=itr increment=chg πⱼ=πⱼ' # Ψⱼ=Ψⱼ
            if (abs(chg) < tol)
                break
            end
        end
        L_old = L

        # 4) Calculate E₍zzᵀₗxᵢωⱼ₎
        for j in 1:k
            W = @view Wⱼ[:,:,j]
            Σ⁻¹ = @view Σ⁻¹ⱼ[:,:,j]
            broadcast!(-,Y,X,view(μⱼ,:,j))
            βⱼ = W'*Σ⁻¹
            Sᵢ = (view(πᵢⱼ, j, :)'.*Y)*Y'./πⱼ[j]
            view(E₍zzᵀₗxᵢωⱼ₎,:,:,j) .= I - βⱼ*W + βⱼ*Sᵢ*transpose(βⱼ)
        end

        # M Step
        for j in 1:k
            πᵢ = @view πᵢⱼ[j, :]
            broadcast!(*,Y,πᵢ',X) # πᵢ'.*X
            QWnum = Y*view(ZIᵀ,:,:,j)
            QWden1 = πⱼ[j]*view(E₍zzᵀₗxᵢωⱼ₎,:,:,j)
            QWden23 = view(E₍zₗxᵢωⱼ₎,:,:,j)'*πᵢ
            QWden = [QWden1 QWden23; QWden23' πⱼ[j]]
            R = QWnum/QWden

            Wⱼ[:,:,j] .= R[:,1:end-1]
            μⱼ[:,j] .= R[:,end]

            Ψⱼ[:,j] .= diag(Y*X' - R*QWnum')./πⱼ[j]
        end
    end
    return FinateMixtureModel([FactorAnalysis(μⱼ[:,j], Wⱼ[:,:,j], Ψⱼ[:,j]) for j in 1:k], πⱼ./n)
end

fit_mm(::Type{FactorAnalysis}, X::AbstractMatrix{T}; m=2, k=1, kwargs...) where {T<:AbstractFloat} =
    fit_mm(FactorAnalysis, X, m; factors=k, kwargs...)


text:
# Counting
# ========

@generated function _temporal_counts!(counts::ContingencyTable{T,N,A}, weights,
                           seqs::Vararg{AbstractVector{Residue},N}) where {T,N,A}
    quote
        @assert N == length(seqs) "Number of residue arrays and table dimension doesn't match."
        # seq_1 = seqs[1]
        # seq_2 = ...
        @nextract $N seq d -> seqs[d]
        len = length(seq_1)
        # @assert len == length(seq_1) "Residue arrays have different lengths"
        # @assert len == length(seq_2) ...
        @nexprs $N d -> @assert len == length(seq_d) "Residue arrays have different lengths."
        if isa(weights, AbstractArray)
            @assert len == length(weights) "Residue array and weights sizes doesn't match."
        end
        _cleanup_temporal!(counts)
        temporal = counts.temporal
        @inbounds @simd for index in 1:length(seq_1)
            # temporal[Int(seq_1[index]), Int(seq_2... += getweight(weights, index)
            @nref($N, temporal, d -> Int(seq_d[index])) += getweight(weights, index)
        end
        counts
    end
end

"""
It populates a `ContingencyTable` (first argument) using the frequencies in the sequences
(last positional arguments). The dimension of the table must match the number of sequences
and all the sequences must have the same length. You must indicate the used weights and
pseudocounts as second and third positional arguments respectively. You can use
`NoPseudofrequencies()` and `NoClustering()` to avoid the use of sequence weighting and
pseudocounts, respectively.
"""
function count!(table::ContingencyTable{T,N,A},
                weights,
                pseudocounts::Pseudocount,
                seqs::Vararg{AbstractVector{Residue},N}) where {T,N,A}
    _temporal_counts!(table, weights, seqs...)
    apply_pseudocount!(table, pseudocounts)
    _update!(table)
    table
end

function count!(table::Counts{T,N,A}, args...) where {T,N,A}
    count!(getcontingencytable(table), args...)
    table
end

# Default counters
# ================

function _count(alphabet::A,
                weights, pseudocounts,
                seqs::Vararg{AbstractVector{Residue},N}
                )::ContingencyTable{Float64,N,A} where {N,A <: ResidueAlphabet}
    table = ContingencyTable(Float64, Val{N}, alphabet)::ContingencyTable{Float64,N,A}
    count!(table, weights, pseudocounts, seqs...)
    table
end

"""
It returns a `ContingencyTable` wrapped in a `Counts` type with the frequencies of residues
in the sequences that takes as arguments. The dimension of the table is equal to the number
of sequences. You can use the keyword arguments `alphabet`, `weights` and `pseudocounts`
to indicate the alphabet of the table (default to `UngappedAlphabet()`), a clustering
result (default to `NoClustering()`) and the pseudocounts (default to `NoPseudocount()`)
to be used during the estimation of the frequencies.
"""
function Base.count(seqs::Vararg{AbstractVector{Residue},N};
                    alphabet::ResidueAlphabet = UngappedAlphabet(),
                    weights = NoClustering(),
                    pseudocounts::Pseudocount = NoPseudocount()) where N
    Counts(_count(alphabet, weights, pseudocounts, seqs...))
end

# Probabilities
# =============

"""
It populates a `ContingencyTable` (first argument) using the probabilities in the sequences
(last positional arguments). The dimension of the table must match the number of sequences
and all the sequences must have the same length. You must indicate the used weights,
pseudocounts and pseudofrequencies as second, third and fourth positional arguments
respectively. You can use `NoClustering()`, `NoPseudocount()` and `NoPseudofrequencies()`
to avoid the use of sequence weighting, pseudocounts and pseudofrequencies, respectively.
"""
function probabilities!(table::ContingencyTable{T,N,A},
                        weights,
                        pseudocounts::Pseudocount,
                        pseudofrequencies::Pseudofrequencies,
                        seqs::Vararg{AbstractVector{Residue},N}) where {T,N,A}
    count!(table, weights, pseudocounts, seqs...)
    normalize!(table)
    apply_pseudofrequencies!(table, pseudofrequencies)
    table
end

function probabilities!(table::Probabilities{T,N,A}, args...) where {T,N,A}
    probabilities!(getcontingencytable(table), args...)
    table
end

# Default probabilities
# =====================

function _probabilities(alphabet::A,
                        weights, pseudocounts, pseudofrequencies,
                        seqs::Vararg{AbstractVector{Residue},N}
                        )::ContingencyTable{Float64,N,A} where {N,A <: ResidueAlphabet}
    table = ContingencyTable(Float64, Val{N}, alphabet)::ContingencyTable{Float64,N,A}
    probabilities!(table, weights, pseudocounts, pseudofrequencies, seqs...)
    table
end

"""
It returns a `ContingencyTable` wrapped in a `Probabilities` type with the frequencies of
residues in the sequences that takes as arguments. The dimension of the table is equal to
the number of sequences. You can use the keyword arguments `alphabet`, `weights`,
`pseudocounts` and `pseudofrequencies` to indicate the alphabet of the table
(default to `UngappedAlphabet()`), a clustering result (default to `NoClustering()`),
the pseudocounts (default to `NoPseudocount()`) and the pseudofrequencies
(default to `NoPseudofrequencies()`) to be used during the estimation of the probabilities.
"""
function probabilities(seqs::Vararg{AbstractVector{Residue},N};
                       alphabet::ResidueAlphabet = UngappedAlphabet(),
                       weights = NoClustering(),
                       pseudocounts::Pseudocount = NoPseudocount(),
                       pseudofrequencies::Pseudofrequencies = NoPseudofrequencies()) where N
    Probabilities(_probabilities(alphabet,weights,pseudocounts,pseudofrequencies,seqs...))
end


text:
#=
These examples use ArbNumerics v0.4.4, no earlier version supports `exp(matrix)`.
So, either first do `pkg> up` or do `pkg> rm ArbNumerics` then `pkg> add ArbNumerics`.
Sometimes, rebuilding with `pkg> build ArbNumerics` cleans up after itself.
All the matrix handling is new code, the examples here are working as of now.
If you find anything unclear or something awry, post an issue and
please use slack's direct message to draw my attention.
=#

#=
Using Float64 matrices with entries from `rand` in 0.0..1.0:
with a   250x250 matrix of Float64s, working with 350 bits used  10 seconds, trace ~ 1.0e+53
with a   500x500 matrix of Float64s, working with 350 bits used  50 seconds, trace ~ 1.0e+108
with a   750x750 matrix of Float64s, working with 350 bits used 140 seconds, trace ~ 1.0e+162
with a 1000x1000 matrix of Float64s, working with 350 bits used 380 seconds, trace ~ 1.0e+215
with a   250x250 matrix of Float64s, working with 700 bits used  25 seconds
with a   500x500 matrix of Float64s, working with 350 bits used  25 seconds
with a   750x750 matrix of Float64s, working with 175 bits used  25 seconds
=#

using ArbNumerics # GenericLinearAlgebra, LinearAlgebra, Readables

#=
ArbNumerics' interactive users see the "viewable precision"; computations utilize "working precision".
By definition, _viewable precision <= working_precision_. They may share the same number; then the
results shown are precisely the computed numerical results. Should they differ, there are "extra bits".
Julia should be restarted just before these precisions are altered.  When the precision is being altered
in the course of iterative refinement, it is the precision of successive values which moves, not their
bedrock relationship.  Individual values may be constructed more or less precisely than their sources.
Arithmetic over values of differing precision yields a result that adopts the lesser precision.
When the precision is lower than that which permits crisp numerical results, the interval known
to contain the actual result eclipses the magnitude of that value.  It appears so: `[+/- 12.1]`.
=#

setextrabits(0)  #  obtain the raw results
setextrabits(24) #  more accurately, less precision is observed
#  This is not a means of protection from numerical misadventure.
#  This may be adaptive accomodation to the partially [un]known.

                                    # Use one of these perspectives on bits that count.
setworkingprecision(ArbFloat, 125)  # This focus is on available computational precision.
setprecision(ArbReal, 100)          # Or prefer less precision with more accuracy.

                                    # the most recent `set_` gives the operative precision
                                    # for constructing new Arb values using type name only
                                    # the setting is shared by all of the Arb types
                                    # and that makes it easier to characterize accuracy

#=
#   During the course of computation, ArbFloat, ArbReal, ArbComplex values
#   may be constructed at various precisions.  Once given, precision stays.
#   Change of precision occurs through construction of that changed value.
=#

a₂  = ArbFloat(1, digits=100, base=2);
a₁₀ = ArbFloat(1, digits=29, base=10);
workingprecision.((a₂, a₁₀))
(100, 98)

b₂  = ArbFloat(1, digits=3000, base=2);
b₁₀ = ArbFloat(1, digits=902, base=10);
workingprecision.((b₂, b₁₀))
(3000, 2999)

workingprecision(ArbReal(1)) == 128
workingprecision(ArbFloat(1, bits=1024) + ArbFloat(1, bits=256)) == 256

#>  ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...  <#

#=
    To regard a matrix of extended precision ArbFloats or ArbReals,
    where the elements are in the Float64 domain: `Float64.(matrix)`.
    To characterize or classify a large matrix of the same nature:
    `approx = Float64.(matrix); result = describe(process(approx))`
#
const Arb{P} = Union{ArbFloat{P}, ArbReal{P}} where {P}
#=
#   Obtaining the exponential of a matrix
#
#     (a) provision a square matrix
#     (b) generate a similar matrix of eltype ArbReal
#     (c) receive each of (b) carried from (a)
#     (d) request exp(c) of ArbNumerics                (faster than exp(Matrix{ArbFloat}))
#     (e) generate a similar matrix given (a)
#     (f) receive each of (e) carried from (d)
=#
function expmatrix(m::Array{T,2}) where {T<:Real}
    a = Array{ArbReal,2}(undef, size(m))
    a[:] = m[:]
    z = exp(a)
    result = similar(m)
    result[:] = z[:]
    return result
end


text:
export CiteSeer

"""
    CiteSeer

The CiteSeer citation network dataset from Ref. [1].
Nodes represent documents and edges represent citation links.
The dataset is designed for the node classification task.
The task is to predict the category of certain paper.
The dataset is retrieved from Ref. [2].

# Interface

- [`CiteSeer.dataset`](@ref)

# References

[1]: [Deep Gaussian Embedding of Graphs: Unsupervised Inductive Learning via Ranking](https://arxiv.org/abs/1707.03815)
[2]: [Planetoid](https://github.com/kimiyoung/planetoid)
"""
module CiteSeer

using DataDeps
using ..MLDatasets: datafile, read_planetoid_data
using DelimitedFiles: readdlm

using PyCall

const DEPNAME = "CiteSeer"
const LINK = "https://github.com/kimiyoung/planetoid/raw/master/data"
const DOCS = "https://github.com/kimiyoung/planetoid"
const DATA = "ind.citeseer." .* ["x", "y", "tx", "allx", "ty", "ally", "graph", "test.index"]

function __init__()
    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $DOCS
        """,
        map(x -> "$LINK/$x", DATA),
        "7f7ec4df97215c573eee316de35754d89382011dfd9fb2b954a4a491057e3eb3",  # if checksum omitted, will be generated by DataDeps
        # post_fetch_method = unpack
    ))
end

"""
    dataset(; dir=nothing, reverse_edges=true)

Retrieve the CiteSeer dataset. The output is a named tuple with fields
```juliarepl
julia> keys(CiteSeer.dataset())
(:node_features, :node_labels, :adjacency_list, :train_indices, :val_indices, :test_indices, :num_classes, :num_nodes, :num_edges, :directed)
```

In particular, `adjacency_list` is a vector of vector,
where `adjacency_list[i]` will contain the neighbors of node `i`
through outgoing edges.

If `reverse_edges=true`, the graph will contain
the reverse of each edge and the graph will be undirected.

See also [`CiteSeer`](@ref).

# Usage Examples

```julia
using MLDatasets: CiteSeer
data = CiteSeer.dataset()
train_labels = data.node_labels[data.train_indices]
```
"""
dataset(; dir=nothing, reverse_edges=true) =
    read_planetoid_data(DEPNAME, dir=dir, reverse_edges=reverse_edges)

end #module


text:
"""
    MRCHeader

An MRC header.

    MRCHeader()
    MRCHeader(; kwargs...)

Construct a header with suitable defaults, optionally substituting any provided entries.

See the [MRC2014 format spec](https://www.ccpem.ac.uk/mrc_format/mrc2014.php) for a
description of the header fields.
"""
mutable struct MRCHeader
    nx::Int32
    ny::Int32
    nz::Int32
    mode::Int32
    nxstart::Int32
    nystart::Int32
    nzstart::Int32
    mx::Int32
    my::Int32
    mz::Int32
    cella_x::Float32
    cella_y::Float32
    cella_z::Float32
    cellb_alpha::Float32
    cellb_beta::Float32
    cellb_gamma::Float32
    mapc::Int32
    mapr::Int32
    maps::Int32
    dmin::Float32
    dmax::Float32
    dmean::Float32
    ispg::Int32
    nsymbt::Int32
    extra1::NTuple{8,UInt8}
    exttyp::String
    nversion::Int32
    extra2::NTuple{84,UInt8}
    origin_x::Float32
    origin_y::Float32
    origin_z::Float32
    map::String
    machst::NTuple{4,UInt8}
    rms::Float32
    nlabl::Int32
    label::NTuple{10,String}
end
function MRCHeader(; kwargs...)
    h = MRCHeader(
        ntuple(_ -> 0, Val(13))..., # nx, ny, nz, mode, nxstart, nystart, nzstart, mx, my, mz, cella_x, cella_y, cella_z
        90, # cellb_alpha
        90, # cellb_beta
        90, # cellb_gamma
        1, # mapc
        2, # mapr
        3, # maps
        0, # dmin
        -1, # dmax < dmin means both are undetermined
        -2, # dmean < min(dmin, dmax) means dmean is undetermined
        1, # ispg, is a volume
        0, # nsymbt, empty extended header
        ntuple(_ -> UInt8(0), Val(8)), # extra
        "", # exttyp
        MRC2014_VERSION, # nversion
        ntuple(_ -> UInt8(0), Val(84)), # extra
        ntuple(_ -> 0, Val(3))..., # origin
        MAP_NAME, # map
        NTuple{4,UInt8}(machstfrombyteorder()), # machst
        -1, # rms < 0 means rms is undetermined
        1, # nlabl
        (
            "Created by MRC.jl at $(Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))",
            ntuple(_ -> "", Val(9))...,
        ), # label
    )
    for (k, v) in pairs(kwargs)
        setproperty!(h, k, v)
    end
    return h
end

Base.sizeof(::Type{MRCHeader}) = 1024

function Base.show(io::IO, ::MIME"text/plain", h::MRCHeader)
    print(io, "MRCHeader(")
    for f in fieldnames(typeof(h))
        v = sprint(show, getfield(h, f); context=io)
        print(io, "\n    $(f) = $(v),")
    end
    print(io, "\n)")
    return nothing
end

function Base.copyto!(hto::MRCHeader, hfrom::MRCHeader)
    hto === hfrom && return hto
    for f in fieldnames(typeof(hto))
        setfield!(hto, f, getfield(hfrom, f))
    end
    return hto
end

function sizeoffield(name, type)
    name === :exttyp && return 4
    name === :map && return 4
    name === :label && return 800
    return sizeof(type)
end

function bytestoentry(name, type, pointer)
    name === :exttyp && return _unsafe_string(pointer, 4)
    name === :map && return _unsafe_string(pointer, 4)
    name === :label && return ntuple(Val(10)) do i
        return strip(_unsafe_string(pointer + (i - 1) * 80, 80), ' ')
    end
    return unsafe_load(convert(Ptr{type}, pointer))
end

function entrytobytes(name, value)
    name === :exttyp && return padtruncto!(Vector{UInt8}(value), 4)
    name === :map && return padtruncto!(Vector{UInt8}(value), 4)
    name === :label && return vcat(ntuple(Val(10)) do i
        return padtruncto!(Vector{UInt8}(value[i]), 80)
    end...)
    return reinterpret(UInt8, [value])
end

@generated function fieldoffsets(::Type{T}) where {T}
    sizes = map(sizeoffield, fieldnames(T), T.types)
    offsets = cumsum(sizes)
    pop!(offsets)
    pushfirst!(offsets, 0)
    return offsets
end

Base.size(h::MRCHeader) = (h.nx, h.ny, h.nz)

Base.length(h::MRCHeader) = prod(size(h))

function Base.ndims(h::MRCHeader)
    mrcdtype = mrcdatatype(h)
    return if mrcdtype == Image
        2
    elseif mrcdtype in (ImageStack, Volume, Unknown)
        3
    else
        error("Data type of $(mrcdtype) is not currently supported.")
    end
end

"""
    cellangles(h::MRCHeader) -> (α, β, γ)

Get cell angles in degrees.
"""
cellangles(h::MRCHeader) = (h.cellb_alpha, h.cellb_beta, h.cellb_gamma)

"""
    cellangles!(h::MRCHeader, (α, β, γ))

Set cell angles in degrees.
"""
function cellangles!(h::MRCHeader, cellb)
    h.cellb_alpha = cellb[1]
    h.cellb_beta = cellb[2]
    h.cellb_gamma = cellb[3]
    return cellb
end
cellangles!(h::MRCHeader, cellb::Number) = cellangles!(h, ntuple(_ -> cellb, Val(3)))

"""
    cellsize(h::MRCHeader) -> (x, y, z)

Get cell dimensions in angstroms.
"""
cellsize(h::MRCHeader) = (h.cella_x, h.cella_y, h.cella_z)

"""
    cellsize!(h::MRCHeader, (x, y, z))

Set cell dimensions in angstroms.
"""
function cellsize!(h::MRCHeader, cella)
    h.cella_x = cella[1]
    h.cella_y = cella[2]
    h.cella_z = cella[3]
    return cella
end
cellsize!(h::MRCHeader, cella::Number) = cellsize!(h, ntuple(_ -> cella, Val(3)))

"""
    gridsize(h::MRCHeader) -> (x, y, z)

Get size of sampling grid in angstroms.
"""
gridsize(h::MRCHeader) = (h.mx, h.my, h.mz)

"""
    gridsize!(h::MRCHeader, (x, y, z))

Set size of sampling grid in angstroms.
"""
function gridsize!(h::MRCHeader, m)
    h.mx = m[1]
    h.my = m[2]
    h.mz = m[3]
    return m
end
gridsize!(h::MRCHeader, m::Number) = gridsize!(h, ntuple(_ -> m, Val(3)))

"""
    origin(h::MRCHeader) -> (x, y, z)

Get phase origin in pixels or origin of subvolume in angstroms.
"""
origin(h::MRCHeader) = (h.origin_x, h.origin_y, h.origin_z)

"""
    origin!(h::MRCHeader, (x, y, z))

Set phase origin in pixels or origin of subvolume in angstroms.
"""
function origin!(h::MRCHeader, origin)
    h.origin_x = origin[1]
    h.origin_y = origin[2]
    h.origin_z = origin[3]
    return origin
end
origin!(h::MRCHeader, origin::Number) = origin!(h, ntuple(_ -> origin, Val(3)))

"""
    start(h::MRCHeader) -> (nx, ny, nz)

Get location of first column, first row, and first section in unit cell.
"""
start(h::MRCHeader) = (h.nxstart, h.nystart, h.nzstart)

"""
    start!(h::MRCHeader, (nx, ny, nz))

Set location of first column, first row, and first section in unit cell.
"""
function start!(h::MRCHeader, nstart)
    h.nxstart = nstart[1]
    h.nystart = nstart[2]
    h.nzstart = nstart[3]
    return nstart
end
start!(h::MRCHeader, nstart::Number) = start!(h, ntuple(_ -> nstart, Val(3)))

"""
    voxelaxes(h::MRCHeader, i) -> StepRangeLen
    voxelaxes(h::MRCHeader) -> NTuple{3,StepRangeLen}

Get range of voxel positions along axis `i` or all axes.
"""
voxelaxes(h::MRCHeader, i) = StepRangeLen(origin(h)[i], voxelsize(h)[i], size(h)[i])
voxelaxes(h::MRCHeader) = ntuple(i -> voxelaxes(h, i), Val(3))

"""
    voxelsize(h::MRCHeader, i)
    voxelsize(h::MRCHeader)

Get size of dimension `i` of voxel in angstroms.
"""
function voxelsize(h::MRCHeader, i)
    cs = cellsize(h)[i]
    return iszero(cs) ? cs : cs / gridsize(h)[i]
end
voxelsize(h::MRCHeader) = ntuple(i -> voxelsize(h, i), Val(3))

"""
    voxelsize!(h::MRCHeader, s, i)
    voxelsize!(h::MRCHeaders, (x, y, z))

Set size of dimension `i` of voxel in angstroms.
"""
function voxelsize!(h::MRCHeader, s, i)
    setproperty!(h, (:cella_x, :cella_y, :cella_z)[i], s * getfield(h, (:mx, :my, :mz)[i]))
    return s
end
function voxelsize!(h::MRCHeader, s)
    voxelsize!(h, s[1], 1)
    voxelsize!(h, s[2], 2)
    voxelsize!(h, s[3], 3)
    return s
end

Base.minimum(h::MRCHeader) = h.dmin

Base.maximum(h::MRCHeader) = h.dmax

Base.extrema(h::MRCHeader) = (minimum(h), maximum(h))

Statistics.mean(h::MRCHeader) = h.dmean

function mrcdatatype(nz, mz, ispg)
    return if ispg == 0
        ifelse(mz == 1, Image, ImageStack)
    elseif ispg == 1 && mz == nz
        Volume
    elseif ispg ∈ 401:630
        VolumeStack
    else
        Unknown
    end
end
mrcdatatype(header) = mrcdatatype(header.nz, header.mz, header.ispg)

function datatype(mode)
    if !haskey(MODE_TO_TYPE, mode)
        error("Mode $(mode) does not correspond to a known data type.")
    end
    return MODE_TO_TYPE[mode]
end
datatype(header::MRCHeader) = datatype(header.mode)


text:
using ContactStructures
using Plots

f = (x,y,z) -> cos(2π*x)
g = (x,y,z) -> sin(2π*x)
giroux_form = f * dz + g * dy

plotly()
plot_along_xy_grid(0:0.125:1,0:0.125:1,0,giroux_form;  markersize = 0,ε = 0.05)
plot!(xlims=(-0.1, 1.1), ylims=(-0.1,1.1), zlims=(-0.6,0.6),legend=false)


text:
# test reservation table
@time @testset "ReservationTable" begin
let
    table = ReservationTable{Int}(10)
    @test is_valid(table)
    for (interval,val) in [
        ((2,4), true),
        ((0,1), true),
        ((5,6), true),
        ]
        @test is_reserved(table,ResourceReservation(1,2,interval)) == !val
        @test reserve!(table,ResourceReservation(1,2,interval)) == val
        @test is_reserved(table,ResourceReservation(1,2,interval))
    end
    @test is_valid(table)
end
let
    mapf = init_mapf_1()
    env = build_env(mapf,initialize_root_node(mapf),1)
    table = ReservationTable{Float64}(num_states(mapf)+num_actions(mapf))
    S = state_type(mapf)
    A = action_type(mapf)
    n = PathNode(S(1,0),A(Edge(1,2),1),S(2,1))
    create_reservations(env,n)
    @test reserve!(table,env,get_s(n),get_a(n),get_sp(n))
    s2 = get_sp(n)
    a2 = CRCBS.wait(env,get_sp(n))
    sp2 = get_next_state(env,s2,a2)
    n2 = PathNode(s2,a2,sp2)
    @test reserve!(table,env,get_s(n2),get_a(n2),get_sp(n2))
end
let
    mapf = init_mapf_1()
    env = MetaAgentCBS.construct_team_env(
        [build_env(mapf,initialize_root_node(mapf),1)],
        [1]
        )
    S = state_type(mapf)
    A = action_type(mapf)
    n = PathNode(
        MetaAgentCBS.State([S(1,0)]),
        MetaAgentCBS.Action([A(Edge(1,2),1)]),
        MetaAgentCBS.State([S(2,1)])
        )
    s2 = get_sp(n)
    a2 = CRCBS.wait(env,get_sp(n))
    sp2 = get_next_state(env,s2,a2)
    n2 = PathNode(s2,a2,sp2)
    let
        table = ReservationTable{Float64}(num_states(mapf)+num_actions(mapf))
        @test reserve!(table,env,get_s(n),get_a(n),get_sp(n))
        @test !reserve!(table,env,get_s(n),get_a(n),get_sp(n))
        @test reserve!(table,env,get_s(n2),get_a(n2),get_sp(n2))
        @test !reserve!(table,env,get_s(n2),get_a(n2),get_sp(n2))
    end
    let
        table = PIBTReservationTable(num_states(mapf),num_actions(mapf))
        @test reserve!(table,env,get_s(n),get_a(n),get_sp(n))
        @test !reserve!(table,env,get_s(n),get_a(n),get_sp(n))
    end
end
let
    table = ReservationTable{Int}(10)
    @test is_valid(table)
    for (interval,val) in [
            ((1,2), true),
            ((1,2), false),
            ((0,1), true),
            ((0,1), false),
            ((2,3), true),
            ((2,3), false),
            ((4,7), true),
            ((5,8), false),
            ((3,6), false),
            ((5,6), false),
            ((2,6), false),
        ]
        @test is_reserved(table,ResourceReservation(1,2,interval)) == !val
        @test reserve!(table,ResourceReservation(1,2,interval)) == val
        @test is_reserved(table,ResourceReservation(1,2,interval))
        # @test is_available(table,ResourceReservation(1,2,interval))
        @test 2 in reserved_by(table,ResourceReservation(1,2,interval))
    end
    @test is_valid(table)
end
let
    table = ReservationTable{Float64}(10)
    @test is_valid(table)
    for (interval, val) in [
        ((1.0, 2.0), true),
        ((1.0, 2.0), false),
        ((0.0, 1.0), true),
        ((0.0, 1.0), false),
        ((2.0, 3.0), true),
        ((2.0, 3.0), false),
        ((4.0, 5.0), true),
        ((3.5, 4.5), false),
        ((4.5, 5.5), false),
        ((4.5, 4.6), false),
        ((3.5, 5.5), false),
    ]
        @test is_reserved(table, ResourceReservation(1, 2, interval)) == !val
        @test reserve!(table, ResourceReservation(1, 2, interval)) == val
        @test is_reserved(table, ResourceReservation(1, 2, interval))
    end
    @test is_valid(table)
end
end # end testset

# For testing Priority Inheritance with Backtracking (PIBT)
@time @testset "PIBT" begin
let
    mapf = init_mapf_1()
    solver = PIBTPlanner{Float64}()
    cache = CRCBS.pibt_init_cache(solver,mapf)
end
# test PIBT on mapf test problems
let
    solver = PIBTPlanner{Float64}()
    set_iteration_limit!(solver,20)
    set_verbosity!(solver,0)
    for f in mapf_test_problems()
        mapf = f()
        reset_solver!(solver)
        solution, valid = pibt!(solver, mapf)
        @test valid
        reset_solver!(solver)
        solution, cost = solve!(solver, mapf)
        @test cost != typemax(cost_type(mapf))
    end
end
# test PIBT with ragged plans on mapf test problems
let
    solver = PIBTPlanner{Float64}()
    set_iteration_limit!(solver,20)
    set_verbosity!(solver,0)
    for f in mapf_test_problems()
        mapf = f()
        reset_solver!(solver)
        solution, valid = pibt!(solver, mapf)
        @test valid
        # remove the path of robot 1 and replan
        p1 = get_paths(solution)[1]
        @show length(p1)
        # for i in reverse(2:length(p1))
        #     deleteat!(p1.path_nodes,i)
        # end
        empty!(p1.path_nodes)
        set_cost!(p1,get_initial_cost(mapf))
        p2 = deepcopy(get_paths(solution)[2])
        # replan with the ragged plan
        reset_solver!(solver)
        # cache = CRCBS.pibt_init_cache(solver,mapf,solution)
        # @show map(string, cache.states)
        # @show map(string, cache.actions)
        # @show cache.active_countdowns
        # CRCBS.pibt_update_cache!(solver,mapf,cache)
        # @show cache.active_countdowns
        # @show map(p->length(p), get_paths(cache.solution))
        solution, valid = pibt!(solver, mapf)
        @test valid
        p2b = get_paths(solution)[2]
        @test length(p2) == length(p2b)
        @test get_cost(p2) == get_cost(p2b)
        for (n2a,n2b) in zip(p2.path_nodes,p2b.path_nodes)
            @test get_s(n2a) == get_s(n2b)
            @test get_a(n2a) == get_a(n2b)
            @test get_sp(n2a) == get_sp(n2b)
        end
    end
end
# Test PIBT on problem from paper
let
    solver = PIBTPlanner{Float64}()
    mapf = init_mapf_5()
    set_max_iterations!(solver,50)
    # set_verbosity!(solver,4)
    solution, status = pibt!(solver,mapf)
end

end # end testset


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

# Note: Most (if not all) of these types are temporary, to be replaced later
# on by generic types from RadiationDetectorSignals.jl (when the latter is
# mature enough).

# Internal units should be SI units
const internal_length_unit  = u"m"
const internal_angle_unit   = u"rad"
const internal_time_unit    = u"s"
const internal_voltage_unit = u"V"
const internal_efield_unit  = u"V / m"
const internal_energy_unit  = u"eV"

to_internal_units(u_internal::Unitful.Units, x::Real) = x
to_internal_units(u_internal::Unitful.Units, x::Quantity) = ustrip(uconvert(u_internal, x))

to_internal_units(u_internal::Unitful.Units, x::AbstractArray{<:Real}) = x
to_internal_units(u_internal::Unitful.Units, x::AbstractArray{<:Quantity}) = ustrip.(uconvert.(u_internal, x))

from_internal_units(u_external::typeof(Unitful.NoUnits), u_internal::Unitful.Units, x::Real) where {T<:Real} = x
from_internal_units(u_external::Unitful.Units, u_internal::Unitful.Units, x::Real) where {T<:Quantity} = uconvert(u_external, x * u_internal)

from_internal_units(u_external::typeof(Unitful.NoUnits), u_internal::Unitful.Units, x::AbstractArray{<:Real}) where {T<:Real} = x
from_internal_units(u_external::Unitful.Units, u_internal::Unitful.Units, x::AbstractArray{<:Real}) where {T<:Quantity} = uconvert.(u_external, x * u_internal)

const MaybeWithUnits{T} = Union{T, Quantity{<:T}}
const RealQuantity = MaybeWithUnits{<:Real}

DetectorHitEvents = TypedTables.Table{
    <:NamedTuple{
        (:evtno, :detno, :thit, :edep, :pos),
        <:Tuple{
            Union{Integer, AbstractVector{<:Integer}},
            Union{Integer, AbstractVector{<:Integer}},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:AbstractVector{<:RealQuantity}}
        }
    }
}


text:
# This file is translated from `aobench.c` by Syoyo Fujita
module AOBench

export Vect, main

using LinearAlgebra, UnPack, Setfield, MuladdMacro, Images, FileIO
using Base.Cartesian

struct Vect{T} <: AbstractVector{T}
    x::T
    y::T
    z::T
end
Vect(x, y, z) = Vect(promote(x, y, z)...,)
Base.getindex(v::Vect, i::Int) = getfield(v, i)
Base.size(::Vect) = (3,)
Base.Tuple(v::Vect) = v.x, v.y, v.z
Vect(xx::Tuple) = Vect(xx...,)
Vect(x::Number) = Vect(@ntuple 3 i->x)

struct Intersect{T,V}
    t::T
    p::V
    n::V
    hit::Int
end

struct Sphere{V,T}
    center::V
    radius::T
end

struct Plane{V}
    p::V
    n::V
end

struct Ray{V}
    org::V
    dir::V
end

struct State{T,V}
    spheres::NTuple{3,Sphere{V,T}}
    plane::Plane{V}
    width      ::Int
    height     ::Int
    nsubsamples::Int
    nao_samples::Int
end

@muladd LinearAlgebra.dot(v0::Vect, v1::Vect) = v0.x * v1.x + v0.y * v1.y + v0.z * v1.z
@muladd LinearAlgebra.cross(v0::Vect, v1::Vect) = Vect(
    v0.y * v1.z - v0.z * v1.y,
    v0.z * v1.x - v0.x * v1.z,
    v0.x * v1.y - v0.y * v1.x,
   )
LinearAlgebra.norm(v::Vect) = sqrt(dot(v, v))
LinearAlgebra.normalize(v::Vect) = v * inv(norm(v))
for op in [:+, :-, :*, :/]
    @eval Base.$op(v0::Vect, v1::Vect) = Vect(@ntuple 3 i -> $op(v0[i], v1[i]))
    @eval Base.$op(v0::Vect, s::Number) = Vect(@ntuple 3 i -> $op(v0[i], s))
    if op !== :/
        @eval Base.$op(s::Number, v0::Vect) = Vect(@ntuple 3 i -> $op(s, v0[i]))
    end
end
Base.:(-)(v::Vect{T}) where T = Vect(@ntuple 3 i -> -v[i])
Base.muladd(v1::Vect, v2::Vect, v3::Vect) = Vect(@ntuple 3 i -> muladd(v1[i], v2[i], v3[i]))
Base.zero(::Type{Vect{T}}) where T = Vect(@ntuple 3 _->zero(T))
Base.zero(::T) where {T<:Vect} = zero(T)

@fastmath @muladd function ray_intersect(isect::Intersect, ray::Ray, sphere::Sphere)
    rs = ray.org - sphere.center
    B = rs'ray.dir
    C = rs'rs - sphere.radius^2
    D = B^2 - C
    if D > zero(D)
        t = -B - sqrt(D)
        if t > zero(t) && t < isect.t
            @set! isect.t = t
            @set! isect.hit = 1

            @set! isect.p = ray.org + ray.dir * t
            @set! isect.n = normalize(isect.p - sphere.center)
        end
    end
    return isect
end

@muladd function ray_intersect(isect::Intersect, ray::Ray, plane::Plane)
    d = -plane.p'plane.n
    v = ray.dir'plane.n

    abs(v) < 1.0e-17 && return isect

    t = -(ray.org'plane.n + d) / v

    if t > zero(t) && t < isect.t
        @set! isect.t = t
        @set! isect.hit = 1

        @set! isect.p = ray.org + ray.dir * t
        @set! isect.n = plane.n
    end
    return isect
end

function orthobasis(n::Vect)
    b1 = zero(n)
    b2 = n
    if -0.6 < n.x < 0.6
        @set! b1.x = 1
    elseif -0.6 < n.y < 0.6
        @set! b1.y = 1;
    elseif -0.6 < n.z < 0.6
        @set! b1.z = 1;
    else
        @set! b1.x = 1
    end
    b0 = normalize(b1 × b2)
    b1 = normalize(b2 × b0)
    return b0, b1, b2
end

@fastmath @muladd function ambient_occlusion(isect::Intersect{T,V}, scene::State) where {T,V}
    @unpack nao_samples = scene
    ntheta = nao_samples
    nphi   = nao_samples
    ε = 0.0001
    p = isect.p + ε * isect.n
    b0, b1, b2 = orthobasis(isect.n)
    occlusion = zero(T)
    for j in 1:ntheta, i in 1:nphi
        theta = sqrt(rand())
        phi = 2pi * rand()

        x = cos(phi) * theta
        y = sin(phi) * theta
        z = sqrt(1 - theta^2)

        dir = x * b0 + y * b1 + z * b2
        ray = Ray(p, dir)

        occisect = Intersect(1e17, zero(V), zero(V), 0)
        occisect = ray_intersect(occisect, ray, scene.spheres[1])
        occisect = ray_intersect(occisect, ray, scene.spheres[2])
        occisect = ray_intersect(occisect, ray, scene.spheres[3])
        occisect = ray_intersect(occisect, ray, scene.plane)
        iszero(occisect.hit) || (occlusion += 1.0)
    end
    occlusion = (ntheta * nphi - occlusion) / (ntheta * nphi)
    return occlusion
end

function render!(img, h::Int, w::Int, nsubsamples::Int, scene::State)
    ss = inv(nsubsamples^2)
    fill!(img, 0)

    for y in 0:h-1, x in 0:w-1
        for v in 0:nsubsamples-1, u in 0:nsubsamples-1
            px = (x + (u / nsubsamples) - (w / 2.0)) / (w / 2.0)
            py = -(y + (v / nsubsamples) - (h / 2.0)) / (h / 2.0)
            dir = normalize(Vect(px, py, -1))
            ray = Ray(zero(dir), dir)
            isect = Intersect(1e17, zero(dir), zero(dir), 0)
            isect = ray_intersect(isect, ray, scene.spheres[1])
            isect = ray_intersect(isect, ray, scene.spheres[2])
            isect = ray_intersect(isect, ray, scene.spheres[3])
            isect = ray_intersect(isect, ray, scene.plane)

            if !iszero(isect.hit)
                col = ambient_occlusion(isect, scene)
                img[y+1, x+1] += Gray(col)
            end
        end

        img[y+1, x+1] *= ss
    end
    return img
end

function init_scene(;width=256, height=256, nsubsamples=2, nao_samples=8)
    spheres = (
               Sphere(Vect(-2.0, 0.0, -3.5), 0.5),
               Sphere(Vect(-0.5, 0.0, -3.0), 0.5),
               Sphere(Vect(1.0, 0.0, -2.2), 0.5),
              )
    plane = Plane(Vect(0.0, -0.5, 0.0), Vect(0.0, 1.0, 0.0))
    return State(spheres, plane, width, height, nsubsamples, nao_samples)
end

function main(;width=256, height=256, nsubsamples=2, nao_samples=8)
    scene = init_scene(;
                        width=width,
                        height=height,
                        nsubsamples=nsubsamples,
                        nao_samples=nao_samples,
                      )
    img = Array{RGB{Float64}}(undef, scene.height, scene.width)
    render!(img, scene.height, scene.width, scene.nsubsamples, scene)
    save(File(format"PNG", "ao.png"), img)
    return img
end

#=
function saveppm(fname::String, w::Int, h::Int, img::Array{UInt8,3})
    open(fname, "w+") do f
        println(f, "P6")
        println(f, "$w $h")
        println(f, "255")
        write(f, vec(UInt8.(img)))
    end
    return fname
end
=#

end # module

max_stars_repo_name: YingboMa/AOBench

max_issues_repo_name: YingboMa/AOBench

max_forks_repo_name: YingboMa/AOBench


text:
# An helper private function
#
# Role: modify src index taken into account that dest component are removed
#
# Algo: for each dest index inferior to a given src index, src index
# must be reduced by one unit.
#
function src_indices_after_dest_removed(src::AbstractVector{Int},
                                        dest::AbstractVector{Int})
    @assert issorted(dest) # only dest has to be sorted (reason: the
                           # insertion algo requires an ordered dest)
    @assert length(src)==length(dest)

    n = length(src)

    count =  zeros(Int,n)
    i_count = 1

    # note this works even if src, dest are not sorted but this has
    # n^2 complexity
    for src_i in src
        for dest_i in dest
            if dest_i < src_i
                count[i_count] += 1
            else
                @assert dest_i != src_i "this makes no sense"
            end
        end
        i_count += 1
    end
    src - count
end

@doc raw"""

Define a map ``g`` of the form
```math
\theta = g_{\hat{\theta}_f}(\hat{\theta}_m)
```

Where the resulting ``\theta`` vector is computed as follows:

- compute ``\tau=f_{\hat{\theta}_f}(\hat{\theta}_m[\text{src}])`` an
  intermediate vector ``\tau`` using the wrapped map ``f`` action on a
  subset of the ``\hat{\theta}_m`` vector. This subset being defined
  by an index array `src`.

- create the resulting ``\theta`` by inserting components of ``\tau``
  into the ``\hat{\theta}_m`` vector at indices stored in an index
  array `dest` (`src` and `dest` have the same length).

This is certainly more comprehensible with an example.

# Example

```jldoctest
using NLS_Fit

f(θ̂m_src,θ̂f) = θ̂f[1] .+ θ̂f[2] * θ̂m_src

f_map = Map_From_VectFunc(2,f)
θ̂f    = Float64[0,-1] # change sign

# insert -θ̂m[4], at position 2
# insert -θ̂m[1], at position 3
src  = [4,1]
dest = [2,3]

# delete elements that are going to be computed & inserted
#
# we see that elements of `src` indices are kept,
# whereas elements of `dest` indices are removed
#
θ̂m = Float64[1:5;]       #  [1,2,3,4,5]
θ̂m = deleteat!(θ̂m,dest)  #  [1,4,5]

map_s_d = NLS_Fit.Transform_Src_Insert_Dest_Map(f_map,src=>dest)

θ = eval_map(map_s_d,θ̂m,θ̂f)

# output
5-element Vector{Float64}:
  1.0
 -4.0
 -1.0
  4.0
  5.0
```

# Constructor

```julia
Transform_Src_Insert_Dest_Map(f_map, src=>dest)
```
- `f_map` is the stored [`Abstract_Map`](@ref) map,

- `src=>dest` is a pair of two `Int` vectors used to define
  indices. Note that we must have `src ∩ dest = ∅`. The reason is that
  as dest indices are removed, there would no more associated src
  component to apply ``f``

# Also see

This `map` can be directly used by [`Model2Fit_Transformed_Parameters`](@ref).

"""
struct Transform_Src_Insert_Dest_Map{MAP<:Abstract_Map} <: Abstract_Map
    _f_map::MAP
    _src_dest::Matrix{Int} # src=col(1), dest=col(2)
end

function Transform_Src_Insert_Dest_Map(f_map::Abstract_Map, src_dest::Pair{SOURCE,DEST}) where {SOURCE<:AbstractVector{Int},DEST<:AbstractVector{Int}}

    (src, dest) = src_dest

    @assert length(src) == length(dest)

    # Form [src,dest] matrix & sort rows wrt dest
    src_dest = hcat(src, dest)
    src_dest = sortslices(src_dest, dims = 1, by = x -> last(x))

    # Update src indices
    src_dest[:, 1] .= src_indices_after_dest_removed(src_dest[:, 1], src_dest[:, 2])

    Transform_Src_Insert_Dest_Map(f_map, src_dest)
end

parameter_size(tm::Transform_Src_Insert_Dest_Map) = parameter_size(tm._f_map)

function eval_map(tm::Transform_Src_Insert_Dest_Map,hat_model_θ::AbstractVector,hat_map_θ::AbstractVector)
    @assert parameter_size(tm) == length(hat_map_θ)

    src  = @view tm._src_dest[:,1]
    dest = @view tm._src_dest[:,2]

    θf = eval_map(tm._f_map,(@view hat_model_θ[src]),hat_map_θ)

    θ = insert_some_elements(hat_model_θ, dest, θf)

    θ
end


text:
include("./exportData.jl")
include("./fenkar_pde.jl")
using .exportData, .fenkar_pde, DifferentialEquations, DiffEqParamEstim, PyPlot, BlackBoxOptim, Optim, DelimitedFiles

#solve and take data
prob = ODEProblem(fenkar!, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=tspan[begin]:5.0:tspan[end])

x = collect(0:(N-1)).*0.015;
plt.pcolormesh(sol.t, x, sol[:,1,:], rasterized=true, snap=true, vmin=0.0, vmax=1.0)
plt.colorbar()
plt.savefig("./pde/bbo/baseline.pdf")
plt.close()

# getExpData
ind=100
#pp = readdlm("./pde/bbo/78_params.txt")
#p .= pp[end,:]
#p[1:13] .= [29.0, 7.0, 7.0, 10.0, 60.0, 250.0, 0.25, 12.0, 100.0, 10.0, 0.13, 0.04, 0.5]
t, data, p[15] = getExpData(ind; tInds=1:250); save_idxs = Int(round(N//2)):Int(round(N//2))

# remake problem
prob = remake(prob, tspan=(0.0, t[end]), p=p)

# use the save_idxs to select the mid-point of the u-variable
sol = solve(prob, Tsit5(); saveat=t, save_idxs=save_idxs)

fig = plt.figure(constrained_layout=true)
plt.plot(t, data[:], "-", label="Exp. Sample")
plt.plot(sol.t, sol[1,:], "--", label="Initial Parameters")
plt.xlabel("\$ t \$")

# define the cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t,data), maxiters=100000, verbose_opt=true, verbose_steps=5, saveat=t, save_idxs=Int(round(N//2)):Int(round(N//2)));

# define limits for parameters
#          tsi,  tv1m,   tv2m,   tvp,   twm,    twp,  td,    to,    tr,   xk,  uc,    uv, ucsi,   t0,       TI,  IA,   θ
#   p = [ 33.4,  25.7, 1162.6,   4.00, 26.7,  662.8, 0.26, 18.7,  29.0, 10.0, 0.26, 0.13,  0.41,  1.0,     BCL, 0.1, 0.0]
lower = [  1.0,  10.0,  400.0,   2.33, 10.0,  100.0, 0.10, 10.0,  10.0,  9.0, 0.10, 0.01,  0.10,  0.0,    40.0, 0.0, -pi]
upper = [ 50.0,  50.0, 2000.0,   4.33, 50.0, 1000.0, 0.50, 50.0,  50.0, 11.0, 0.90, 0.90,  0.90,  t[1], 1000.0, 1.0, +pi]

#@assert all(p .>= lower) && all(p .<= upper)

bound = [(lower[n], upper[n]) for n in 1:length(p)]

# fit
result = bboptimize(cost_function; SearchRange = bound, MaxSteps = 5e4)#, NThreads=Threads.nthreads()-1)
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.archive_output.best_candidate) \n"
prob = remake(prob, p=result.archive_output.best_candidate)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=save_idxs)
plt.plot(sol.t, sol[1,:], "--", label="BBO Parameters")

result = optimize(cost_function, lower, upper, result.archive_output.best_candidate, Fminbox(NelderMead()))
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.minimizer) \n"
prob = remake(prob, p=result.minimizer)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=save_idxs)
plt.plot(sol.t, sol[1,:], "--", label="Nelder-Mead Parameters")

plt.savefig("./pde/bbo/$(ind)_fittings.pdf")
plt.close()

open("./pde/bbo/$(ind)_params.txt", "w") do io
	writedlm(io, ["tsi" "tv1m" "tv2m" "tvp" "twm" "twp" "td" "to" "tr" "xk" "uc" "uv" "ucsi" "t0" "TI" "IA"])
	writedlm(io, "\n")
	writedlm(io, transpose(p))
	writedlm(io, "\n")
	writedlm(io, transpose(result.minimizer))
end


text:
#println("Works")

function reetertr() #args...)
#println("Test Works ", args...)
#return 1234567890
end

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
################################################################################
#                          TEST: DCM <=> Euler Angles
################################################################################

for k = 1:samples
    for rot_seq in rot_seq_array
        # Sample three angles form a uniform distribution [-pi,pi].
        θx       = -pi + 2*pi*rand()
        θy       = -pi + 2*pi*rand()
        θz       = -pi + 2*pi*rand()
        eulerang = EulerAngles(θx, θy, θz, rot_seq)

        # Get the error matrix related to the conversion from DCM => Euler
        # Angles => DCM.
        error1 = angle2rot(DCM, dcm2angle(angle2rot(eulerang),rot_seq)) -
                 angle2rot(DCM, θx, θy, θz, rot_seq)
        error2 = angle2rot(dcm2angle(angle2rot(θx, θy, θz, rot_seq),rot_seq)) -
                 angle2rot(eulerang)

        # If everything is fine, the norm of the matrix error should be small.
        @test norm(error1) < 1e-10
        @test norm(error2) < 1e-10
        @test error1 ≈ error2
    end
end

for k = 1:samples
    # Sample three angles form a uniform distribution [-0.0001,0.0001].
    θx       = -0.0001 + 0.0002*pi*rand()
    θy       = -0.0001 + 0.0002*pi*rand()
    θz       = -0.0001 + 0.0002*pi*rand()
    eulerang = EulerAngles(θx, θy, θz, :XYZ)

    # Get the error between the exact rotation and the small angle
    # approximation.
    error1 = angle2rot(eulerang) -
             smallangle2rot(DCM, eulerang.a1, eulerang.a2, eulerang.a3)
    error2 = angle2rot(θx, θy, θz, :XYZ) -
             smallangle2rot(eulerang.a1, eulerang.a2, eulerang.a3)

    # If everything is fine, the norm of the matrix error should be small.
    #
    # TODO: 2018-06-21: Some tests were failing with the tolerance `5e-7` in
    # nightly builds in Windows. Those failures were not seen in Linux or macOS.
    @test norm(error1) < 7e-7
    @test norm(error2) < 7e-7
    @test error1 ≈ error2
end

# Test exceptions.
@test_throws ArgumentError angle2dcm(0,0,0,:xyz)
@test_throws ArgumentError angle2dcm(0,0,0,:zyx)
@test_throws ArgumentError angle2dcm(0,0,0,:xyx)
@test_throws ArgumentError angle2dcm(0,0,0,:abc)
@test_throws ArgumentError angle2quat(0,0,0,:xyz)
@test_throws ArgumentError angle2quat(0,0,0,:zyx)
@test_throws ArgumentError angle2quat(0,0,0,:xyx)
@test_throws ArgumentError angle2quat(0,0,0,:abc)


text:
import Core: CodeInfo, svec, ReturnNode, SlotNumber, SSAValue, LineInfoNode

function forward_callable(@nospecialize f)
    @assert isdefined(Base, nameof(f))
    for m in methods(f)
        forward_method(m)
    end
end

function forward_method(m::Method)
    local src
    try
        src = make_forward_method(m)
        isnothing(src) && return
        Core.eval(@__MODULE__, Expr(:thunk, src))
    catch err
        if !@isdefined(src)
            @error "failed to create a method to forward `Core.Compiler` types, inspect `$(@__MODULE__).m`"
            Core.eval(@__MODULE__, :(m = $m))
        else
            @error "failed to define a method to forward `Core.Compiler` types, inspect `$(@__MODULE__).m` and `$(@__MODULE__).src`"
            Core.eval(@__MODULE__, :(m = $m; src = $(QuoteNode(src))))
        end
        rethrow(err)
    end
end

function make_forward_method(m::Method)
    # step 0: check if the exact same method already exists
    # the method might be defined already, because:
    # - some are defined manually within `Base.IRShow`
    # - some methods work only with `Core` types
    f = getfield(Base, m.name)
    let
        tt1 = Base.unwrap_unionall(m.sig).parameters[2:end]
        for m in methods(f)
            tt2 = Base.unwrap_unionall(m.sig).parameters[2:end]
            if tt1 == tt2
                # this method is defined already, because:
                # - some are defined manually within `Base.IRShow`
                # - some methods work only with `Core` types
                return nothing
            end
        end
    end

    # step 1. change the method signature to work with `Core.Compiler` types
    msig = m.sig
    sparams = TypeVar[]
    while isa(msig, UnionAll)
        push!(sparams, msig.var)
        msig = msig.body
    end
    msig = msig::DataType
    atype = collect(msig.parameters)
    atype[1] = typeof(f)
    # step 0 is not enough to catch all cases, mostly because of typevars with different identities
    # and so here we make an extra effort to avoid such methods by checking if the signature
    # actually interacts with any of `Core.Compiler` types
    local shouldforward = false
    for t in atype
        if is_corecompiler_type(t)
            shouldforward |= true
            break
        end
    end
    shouldforward || return nothing
    atype = svec(atype...)
    sparams = svec(sparams...)
    loc = LineNumberNode(@__LINE__, @__FILE__)
    sig = svec(atype, sparams, loc)

    # step 2. create a method body to forward the call to the corresponding `Core.Compiler` definition
    slotnames = ccall(:jl_uncompress_argnames, Vector{Symbol}, (Any,), m.slot_syms)
    body = let body = Any[]
        local bodyidx::Int = 0
        function insertbody!(@nospecialize x)
            push!(body, x)
            bodyidx += 1
        end
        if m.nospecialize ≠ 0
            nospecialize = Expr(:meta, :nospecialize)
            for i in 1:m.nargs
                if m.nospecialize & (1 << (i - 1)) ≠ 0
                    push!(nospecialize.args, SlotNumber(i))
                end
            end
            insertbody!(nospecialize)
        end
        if m.isva
            if m.nargs-1 > 1
                tpl = Expr(:call, GlobalRef(Core, :tuple))
                for i in 2:(m.nargs-1)
                    push!(tpl.args, SlotNumber(i))
                end
                insertbody!(tpl)
                call = Expr(:call, GlobalRef(Core, :_apply_iterate),
                    GlobalRef(Core.Compiler, :iterate), GlobalRef(Core.Compiler, m.name),
                    SSAValue(bodyidx), SlotNumber(Int(m.nargs)))
                insertbody!(call)
            else
                call = Expr(:call, GlobalRef(Core, :_apply_iterate),
                    GlobalRef(Core.Compiler, :iterate), GlobalRef(Core.Compiler, m.name),
                    SlotNumber(Int(m.nargs)))
                insertbody!(call)
            end
        else
            call = Expr(:call, GlobalRef(Core.Compiler, m.name))
            for i in 2:m.nargs
                push!(call.args, SlotNumber(i))
            end
            insertbody!(call)
        end
        insertbody!(ReturnNode(SSAValue(bodyidx)))
        make_codeinfo(body;
                      slotnames,
                      inlineable=true, # force inlining and eliminate the (possible) overhead of the forwarding
                      )
    end

    # step 3. make a definition thunk
    def = Any[]
    local defidx::Int = 0
    function insertdef!(@nospecialize x)
        push!(def, x)
        defidx += 1
    end
    insertdef!(Expr(:method, nothing, sig, body))
    insertdef!(ReturnNode(SSAValue(defidx)))

    return make_codeinfo(def)
end

function make_codeinfo(code::Vector{Any};
                       slotnames::Vector{Symbol} = Symbol[],
                       inlineable::Bool = false,
                       )
    codeinfo = ccall(:jl_new_code_info_uninit, Ref{CodeInfo}, ())
    codeinfo.code = code
    n = length(codeinfo.code)
    codeinfo.codelocs = Int32[1 for _ in 1:n]
    codeinfo.linetable = Any[LineInfoNode(@__MODULE__, :none, Symbol(@__FILE__), @__LINE__, 0)]
    codeinfo.ssavaluetypes = n
    codeinfo.ssaflags = UInt8[0x00 for _ in 1:n]
    codeinfo.slotnames = slotnames
    codeinfo.slotflags = UInt8[0x00 for _ in 1:length(slotnames)]
    if inlineable
        codeinfo.inlineable = inlineable
        for i in 1:length(codeinfo.ssaflags)
            codeinfo.ssaflags[i] |= Core.Compiler.IR_FLAG_INLINE
        end
    end
    return codeinfo
end

let
    mods = Module[]
    for name in names(Core.Compiler; all=true)
        x = getfield(Core.Compiler, name)
        if isa(x, Module)
            push!(mods, x)
        end
    end

    global function is_corecompiler_type(@nospecialize t)
        if t === Union{}
            return false
        elseif Base.isvarargtype(t)
            return is_corecompiler_type(t.T)
        elseif isa(t, Union)
            return is_corecompiler_type(t.a) || is_corecompiler_type(t.b)
        elseif isa(t, TypeVar)
            return is_corecompiler_type(t.lb) || is_corecompiler_type(t.ub)
        elseif isa(t, UnionAll)
            return is_corecompiler_type(Base.unwrap_unionall(t))
        end
        return parentmodule(t) in mods
    end
end

max_stars_repo_name: aviatesk/CompilerBridge

max_issues_repo_name: aviatesk/CompilerBridge

max_forks_repo_name: aviatesk/CompilerBridge


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;

using Test
using Documenter
using SafeTestsets

@time @safetestset "Вспомогательные                                  " begin include("utils.jl") end

@time @safetestset "Якобиан прямой задачи на статической сетке.      " begin include("direct_jacobian.jl") end
@time @safetestset "Якобиан прямой задачи на динамической сетке.     " begin include("direct_jacobian_nonuniform.jl") end
@time @safetestset "Модельная невзяка для прямой задачи              " begin include("direct_check.jl") end

@time @safetestset "Якобиан сопряженной задачи на статической сетке. " begin include("adjoint_jacobian.jl") end
@time @safetestset "Модельная невязка для сопряженной задачи         " begin include("adjoint_check.jl") end

@time @safetestset "Неоднородность на точных данных                  " begin include("heterogeneity_check.jl") end
@time @safetestset "Новая аппроксимация дельта-функции               " begin include("delta_check.jl") end
@time @safetestset "Вырожденые корни                                 " begin include("degenerated_check.jl") end

@time @safetestset "Градиент, функционал                             " begin include("functional.jl") end


text:
module MooreJuliaMultigrid

using LinearAlgebra
using SparseArrays

"""
Run the Weighted Jacobi method on matrix system Ax=b for a total of numiters
    iterations, starting with an initial guess of x0
"""
function Jacobi(A,b,x0; numiters=1,weight=1)
    wDinv = weight*inv(Diagonal(diag(A)))
    R = UniformScaling(1-weight) + wDinv*(A - Diagonal(diag(A)))
    for i = 1:numiters
        x0 = wDinv*b - R*x0
    end
    return x0
end

"""
V-cycle method to solve the matrix system Ax=b
Inputs:
    A, b, x0 : A matrix, RHS vector, initial guess
    relaxFunc : Function to complete the relaxation
      -- argument form : relaxFunc(A,b,x0; named options)
    preRelaxOpts : Options for relaxFunc before coarse-grid transfer
      -- named tuple with named options for relaxFunc
    postRelaxOpts : Options for relaxFunc after coarse-grid correction
      -- named tuple with named options for relaxFunc
    prolongOp : Prolongation (interpolation) operator
      -- matrix expected - may change
    restrictOp : Restriction operator
      -- matrix expected - may change
"""
function Vcycle(A,b,x0,relaxFunc,preRelaxOpts, postRelaxOpts, prolongOp, restrictOp)
    # First relax numprerelax times
    if preRelaxOpts[:numiters] > 0
        x0 = relaxFunc(A,b,x0;preRelaxOpts...)
    end

    # Transfer residual to coarser grid
    r = restrictOp*(b - A*x0)

    # Transfer A to coarser grid
    Acoarse = restrictOp*A*prolongOp

    # Solve the coarse problem
    ecourse = Acoarse \ r

    # Interp ecourse to fine grid
    e = prolongOp*ecourse

    # Correct approximation with the residual
    x0 = x0 + e

    # Relax numpostrelax times
    if postRelaxOpts[:numiters] > 0
        x0 = relaxFunc(A,b,x0;postRelaxOpts...)
    end
    return x0
end

"""
BuildLinearInterpOp builds the linear interpoltion operator for transferring
    a vector of length fineSize to a vector of length (fineSize+1)/2 - 1.

    Currently expects only an odd length vector

    TODO: Add support for even length vectors
"""
function BuildLinearInterpOp(fineSize)
    numCols::Int = (fineSize+1)/2 - 1
    I = Vector{Int}(undef,3*numCols)
    J = Vector{Int}(undef,3*numCols)
    V = Vector{Float64}(undef,3*numCols)
    for i = 1:numCols
        I[3i-2:3i] = [2i-1; 2i; 2i+1]
        J[3i-2:3i] = i*ones(3,1)
        V[3i-2:3i] = [1; 2; 1]
    end
    return 0.5*sparse(I,J,V)
end

export Jacobi, Vcycle, BuildLinearInterpOp

end # module


text:
include("distribution.jl")

export TE

ϵ = 0.0000001

function TE_ML(data::Matrix{Int64}; base=2) # D_KL(p(x|y,z)||p(x|y))

  x2         = data[2:end,1]
  x1         = data[1:end-1,1]
  y1         = data[1:end-1,2]
  x_max      = maximum(data[:,1])
  y_max      = maximum(data[:,2])
  px2x1y1    = Shannon.fe3p(hcat(x2,x1,y1))
  px1y1      = sum(px2x1y1, 1)
  px2x1      = sum(px2x1y1, 3)
  px1        = sum(px2x1y1, (1,3))

  px2_c_x1y1 = zeros(x_max, x_max, y_max)
  px2_c_x1   = zeros(x_max, x_max)

  for x2_index = 1:x_max
    for x1_index = 1:x_max
      px2_c_x1[x2_index, x1_index] = px2x1[x2_index, x1_index, 1] / px1[1, x1_index, 1]
      for y1_index = 1:y_max
        px2_c_x1y1[x2_index, x1_index, y1_index] = px2x1y1[x2_index, x1_index, y1_index] / px1y1[1, x1_index, y1_index]
      end
    end
  end

  r = 0
  for x2_index = 1:x_max
    for x1_index = 1:x_max
      for y1_index = 1:y_max
        if abs(px2x1y1[x2_index, x1_index, y1_index]) > ϵ && abs(px2_c_x1y1[x2_index, x1_index, y1_index]) > ϵ && abs(px2_c_x1[x2_index, x1_index]) > ϵ
          r = r + px2x1y1[x2_index, x1_index, y1_index] * (log(base, px2_c_x1y1[x2_index, x1_index, y1_index]) - log(base, px2_c_x1[x2_index, x1_index]))
        end
      end
    end
  end
  return r
end

  # Transfer Entropy
  function TE(data::Matrix{Int64}; base=2, mode="ML", pseudocount=0)
    @assert (size(data)[2] == 2) "TE required two columns (X,Y) to calcualte TE(X->Y)."

    modes      = ["ML", "Maximum Likelihood"]
    umodes     = map(x->uppercase(x), modes)
    known_mode = uppercase(mode) in umodes
    pmodes     = map(x->", $x",modes)
    pmodes     = foldl(*, pmodes[1][3:end], pmodes[2:end])
    @assert known_mode "Mode may be any of the following: [$pmodes]."

    if uppercase(mode) in umodes[1:2]
      return TE_ML(data, base)
    end
  end


text:
using FFTW
using MAT

@testset "calculateSNR" begin
    # MATLAB output
    test = matread("resources/calculateSNR.mat")

    # low pass
    u = 0.5 * sin.(2*pi*85/8192 * (0:8191))
    ntf = synthesizeNTF(5, 32, 1)
    v, = simulateDSM(u, ntf)
    spec = fft(v .* ds_hann(8192)/2048)
    snr = calculateSNR(spec[3:129], 83)
    @test isapprox(test["snr_low"], snr)

    # band pass
    ntf = synthesizeNTF(8, 64, 1, 1.5, 0.125)
    u = 0.5 * sin.(2*pi*1045/8192 * (0:8191))
    v, = simulateDSM(u, ntf)
    spec = fft(v .* ds_hann(8192)/2048)
    snr = calculateSNR(spec[992:1056], 54)
    @test isapprox(test["snr_band"], snr)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
## Load the packages.
using LoopVectorization
using Interpolations
using OffsetArrays
using Statistics
using BenchmarkTools
using PyPlot

## Create the refine functions.
function refine_field_int!(hi, hi_tmp, lo, n_hi, n_lo)
    if n_hi ÷ n_lo == 2
        @tturbo for k in 0:size(lo, 3)-3, j in 0:size(lo, 2)-3, i in 0:size(lo, 1)-3
            i_lo = i + 2
            j_lo = j + 2
            k_lo = k + 2
            i_hi = i*2 + 2
            j_hi = j*2 + 2
            k_hi = k*2 + 2
            for kk in 0:1, jj in 0:1, ii in 0:1
                hi_tmp[i_hi+ii, j_hi+jj, k_hi+kk] = lo[i_lo, j_lo, k_lo]
                # hi[i_hi+ii, j_hi+jj, k_hi+kk] = lo[i_lo, j_lo, k_lo]
            end
        end

        # @tturbo for k in 2:size(hi, 3)-1, j in 2:size(hi, 2)-1, i in 2:size(hi, 1)-1
        #     hi[i, j, k] = hi_tmp[i, j, k]
        # end

        @tturbo a_hi_tmp[1, :, :] = a_hi_tmp[end-1, :, :]
        @tturbo a_hi_tmp[end, :, :] = a_hi_tmp[2, :, :]
        @tturbo a_hi_tmp[:, 1, :] = a_hi_tmp[:, end-1, :]
        @tturbo a_hi_tmp[:, end, :] = a_hi_tmp[:, 2, :]
        @tturbo a_hi_tmp[:, :, 1] = a_hi_tmp[:, :, end-1:]
        @tturbo a_hi_tmp[:, :, end] = a_hi_tmp[:, :, 2]

        coef = OffsetArray(zeros(3, 3, 3), -1:1, -1:1, -1:1)
        for kk in -1:1, jj in -1:1, ii in -1:1
            coef[ii, jj, kk] = (1//2)^(3 + abs(ii) + abs(jj) + abs(kk))
        end

        @tturbo for k in 2:size(hi, 3)-1, j in 2:size(hi, 2)-1, i in 2:size(hi, 1)-1
            hi[i, j, k] = 0
            for kk in -1:1, jj in -1:1, ii in -1:1
                hi[i, j, k] += coef[ii, jj, kk] * hi_tmp[i+ii, j+jj, k+kk]
            end
        end
    elseif n_hi ÷ n_lo == 3
        @tturbo for k in 0:size(lo, 3)-3, j in 0:size(lo, 2)-3, i in 0:size(lo, 1)-3
            i_lo = i + 2
            j_lo = j + 2
            k_lo = k + 2
            i_hi = i*3 + 2
            j_hi = j*3 + 2
            k_hi = k*3 + 2
            for kk in 0:2, jj in 0:2, ii in 0:2
                hi_tmp[i_hi+ii, j_hi+jj, k_hi+kk] = lo[i_lo, j_lo, k_lo]
                # hi[i_hi+ii, j_hi+jj, k_hi+kk] = lo[i_lo, j_lo, k_lo]
            end
        end

        # @tturbo for k in 2:size(hi, 3)-1, j in 2:size(hi, 2)-1, i in 2:size(hi, 1)-1
        #     hi[i, j, k] = hi_tmp[i, j, k]
        # end

        @tturbo a_hi_tmp[1, :, :] = a_hi_tmp[end-1, :, :]
        @tturbo a_hi_tmp[end, :, :] = a_hi_tmp[2, :, :]
        @tturbo a_hi_tmp[:, 1, :] = a_hi_tmp[:, end-1, :]
        @tturbo a_hi_tmp[:, end, :] = a_hi_tmp[:, 2, :]
        @tturbo a_hi_tmp[:, :, 1] = a_hi_tmp[:, :, end-1]
        @tturbo a_hi_tmp[:, :, end] = a_hi_tmp[:, :, 2]

        @tturbo for k in 2:size(hi, 3)-1, j in 2:size(hi, 2)-1, i in 2:size(hi, 1)-1
            hi[i, j, k] = 0
            for kk in -1:1, jj in -1:1, ii in -1:1
                hi[i, j, k] += 1//27*hi_tmp[i+ii, j+jj, k+kk]
            end
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

function refine_field_int_ref!(hi, hi_tmp, lo, x_hi, x_lo, n_hi, n_lo)
    # interp = LinearInterpolation((x_lo, x_lo, x_lo), lo)
    interp = interpolate((x_lo, x_lo, x_lo), lo, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    # interp = interpolate((x_lo, x_lo, x_lo), lo, (Gridded(Constant()), Gridded(Constant()), Gridded(Constant())))
    # dx_lo = x_lo[2] - x_lo[1]
    # interp = scale(interpolate(lo, (BSpline(Linear()), BSpline(Linear()), BSpline(Linear()))),
    #                -dx_lo/2:dx_lo:1+dx_lo, -dx_lo/2:dx_lo:1+dx_lo, -dx_lo/2:dx_lo:1+dx_lo)
    hi[2:end-1, 2:end-1, 2:end-1] = interp(x_hi[2:end-1], x_hi[2:end-1], x_hi[2:end-1])
end

## Set up the grids.
n_hi = 384
n_lo = n_hi ÷ 2

a_lo = rand(n_lo, n_lo, n_lo)
a_hi = zeros(n_hi, n_hi, n_hi)

a_hi_gc = zeros(n_hi+2, n_hi+2, n_hi+2)
a_lo_gc = zeros(n_lo+2, n_lo+2, n_lo+2)

a_hi_gc[2:n_hi+1, 2:n_hi+1, 2:n_hi+1] = a_hi[:, :, :]
a_lo_gc[2:n_lo+1, 2:n_lo+1, 2:n_lo+1] = a_lo[:, :, :]

## Cyclic BCs in all 3 directions.
a_hi_gc[1, :, :] = a_hi_gc[end-1, :, :]
a_hi_gc[end, :, :] = a_hi_gc[2, :, :]
a_hi_gc[:, 1, :] = a_hi_gc[:, end-1, :]
a_hi_gc[:, end, :] = a_hi_gc[:, 2, :]
a_hi_gc[:, :, 1] = a_hi_gc[:, :, end-1]
a_hi_gc[:, :, end] = a_hi_gc[:, :, 2]

a_lo_gc[1, :, :] = a_lo_gc[end-1, :, :]
a_lo_gc[end, :, :] = a_lo_gc[2, :, :]
a_lo_gc[:, 1, :] = a_lo_gc[:, end-1, :]
a_lo_gc[:, end, :] = a_lo_gc[:, 2, :]
a_lo_gc[:, :, 1] = a_lo_gc[:, :, end-1]
a_lo_gc[:, :, end] = a_lo_gc[:, :, 2]

dx_hi = 1/n_hi; dx_lo = 1/n_lo
x_hi = dx_hi/2:dx_hi:1 |> collect
x_lo = dx_lo/2:dx_lo:1 |> collect
xh_hi = 0:dx_hi:1 |> collect
xh_lo = 0:dx_lo:1 |> collect

x_hi_gc = -dx_hi/2:dx_hi:1+dx_hi |> collect
x_lo_gc = -dx_lo/2:dx_lo:1+dx_lo |> collect

## Compute.
a_hi_tmp = zeros(size(a_hi_gc))
@btime refine_field_int!(a_hi_gc, a_hi_tmp, a_lo_gc, n_hi, n_lo)
a_hi_int = a_hi_gc[2:end-1, 2:end-1, 2:end-1]

## Compute a reference using Interpolations.jl
a_hi_tmp = zeros(size(a_hi_gc))
a_hi_gc[:, :, :] .= 0.
@btime refine_field_int_ref!(a_hi_gc, a_hi_tmp, a_lo_gc, x_hi_gc, x_lo_gc, n_hi, n_lo)
a_hi_int_ref = a_hi_gc[2:end-1, 2:end-1, 2:end-1]

println("Mean equal to lo: ", mean(a_lo) ≈ mean(a_hi_int))
println("Mean equal to ref: ", mean(a_hi_int_ref) ≈ mean(a_hi_int))
println("Interpolations equal: ", a_hi_int ≈ a_hi_int_ref)

## Plot the output.
figure(figsize=(10, 5))
subplot(121)
pcolormesh(xh_lo, xh_lo, a_lo[:, :, 1], vmin=0, vmax=1)
subplot(122)
pcolormesh(xh_hi, xh_hi, a_hi_int[:, :, 1], vmin=0, vmax=1)
tight_layout()
display(gcf())
show()


text:
@symbol_func function q_95(cur_reactor::AbstractReactor)
  cur_q = 5.0

  cur_q *= cur_reactor.B_0

  cur_q *= a(cur_reactor) ^ 2

  cur_q *= _f_q(cur_reactor)

  cur_q /= cur_reactor.R_0

  cur_q /= cur_reactor.I_P

  cur_q
end


text:
# coding:utf-8

@windows_only println("unix環境推奨")

""" @swap! """
macro swap!(x,y)
  quote
    $(esc(x)) = $(eval(y))
    $(esc(y)) = $(eval(x))
    return
  end
end

""" @cd """
macro cd(obj)
  :(cd($obj))
end

""" @pwd """
macro pwd()
  :(pwd())
end

""" @head """
macro head(filename)
  quote
    $(run(`head $filename -n 3`))
  end
end

""" @dir """
macro dir(pass...)
  return readdir(pass...)
end

""" @replace! """
macro replace!(string, pat, r)
  quote
    $(esc(string)) = $(eval(replace(eval(string),eval(pat), eval(r))))
    return
  end
end

""" @assign! """
macro assign!(a, b)
  quote
    $(esc(a)) = $(eval(b))
    return
  end
end

""" @toCmd """
macro toCmd(string)
  return `$string`
end

""" @run """
macro run(cmd)
  quote
    $(run(eval(cmd)))
  end
end

""" @split! """
macro split!(string, chars)
  quote
    $(esc(string)) = $(eval(split(eval(string),eval(chars))))
    return
  end
end

""" @rstrip! """
macro rstrip!(string)
  quote
    $(esc(string)) = $(eval(rstrip(eval(string))))
    return
  end
end

""" @println """
macro println(obj...)
  :(println($obj...))
end

""" @print """
macro print(obj...)
  :(print($obj...))
end

""" @inc! """
macro inc!(x)
   :($x += 1)
end

""" @decr! """
macro decr!(x)
  :($x -= 1)
end

""" iconv """
function iconv(inputfile, outputfile, from, to)
  return readall(pipeline(`iconv -f $from -t $to $inputfile`, outfile))
end

# convert any character encoding to utf8
""" readall_utf8 """
function readall_utf8(filename)
  return readall(pipeline(`iconv -t UTF-8 $filename`))
end

function readall_utf8(filename, from)
  return readall(pipeline(`iconv -f $from -t UTF-8 $filename`))
end

""" cabocha_parser """
function cabocha_parser(sentence, option)
  return readall(pipeline(`echo $sentence`, `cabocha -$option`))
end

function cabocha_parser(sentence)
  sentence = replace(sentence, "\n", "")
  result = readall(pipeline(`echo $sentence`, `cabocha -f1`))
  d = Dict(:chunk_id => Array{Int, 1}(), :chunk => Array{UTF8String,1}(), :tok_surface => Array{Array{UTF8String,1},1}(),
           :tok_feature => Array{Array{UTF8String,1},1}(), :link => Array{Int, 1}(), :head => Array{Int, 1}(), :func => Array{Int, 1}())
  for line in split(result, "\n")
    line == "EOS" && break
    if ismatch(r"^\* \d+", line)
     tmp = split(line)
     push!(d[:chunk_id], parse(tmp[2])+1)
     push!(d[:chunk], "")
     push!(d[:link], parse(tmp[3][1:end-1])+1)
     head, func = map(parse, split(tmp[4], "/"))
     push!(d[:head], head+1)
     push!(d[:func], func+1)
     push!(d[:tok_surface], Array{UTF8String,1}())
     push!(d[:tok_feature], Array{UTF8String,1}())
    else
     surface, feature = split(line, "\t")
     d[:chunk][end] = d[:chunk][end]*surface
      push!(d[:tok_surface][end], surface)
      push!(d[:tok_feature][end], feature)
    end
  end
  return d
end


text:
function Mdot(m, v::Array{Float64, 1})
  w = ones(size(m, 1), size(m, 2))
  for i in 1:size(m, 1)
    for j in 1:size(m, 2)
      w[i, j] = (m[i, j, :]'*v)[1]
    end
  end
  return w
end


text:
# BASED ON TIKZPICTURES.JL TESTS

using TikzCDs
using TikzCDs: Styles
using Test

for file in ["testCD.pdf", "testCDDoc.pdf", "testCD.tex", "triangle.pdf", "triangle.svg"]
  if isfile(file)
    rm(file)
  end
end

# check that the TEX file contains the desired environments
function has_environment(content::String, environment::String)
    has_begin = occursin("\\begin{$environment}", content)
    has_end = occursin("\\end{$environment}", content)
    if has_begin && has_end
        return true # has both
    elseif !has_begin && !has_end
        return false # has neither
    else
        error("\\begin{$environment} and \\end{$environment} do not match")
    end
end

# Test tikz-cd

data = L"""
A \arrow[rd] \arrow[r, "\phi"] & B \\
                               & C
"""
tp = TikzCD(data, options="scale=0.25")

td = TikzDocument()
push!(td, tp, caption="hello")

save(TEX("testCD"), tp)
@test isfile("testCD.tex")

filecontent = join(readlines("testCD.tex", keep=true)) # read with line breaks
@test occursin(strip(data[2:end-1]), filecontent) # also check that the data is contained
@test has_environment(filecontent, "tikzcd")
@test has_environment(filecontent, "document")

if success(`lualatex -v`)
    save(PDF("testCD"), tp)
    @test isfile("testCD.pdf")

    save(PDF("testCDDoc"), td)
    @test isfile("testCDDoc.pdf")
else
    @warn "lualatex is missing; can not test compilation"
end

# Test q.uiver.app support
function testrender(diagram, name::String)
    if success(`lualatex -v`)
        save(PDF(name), diagram)
        @test isfile("$name.pdf")
        save(SVG(name), diagram)
    else
        @warn "lualatex is missing; can not test compilation"
    end
end

triangle = L"""        A &&& Q \\
               \\
               &&& P
               \arrow["h", from=1-4, to=3-4]
               \arrow["f", from=1-1, to=1-4]
               \arrow["g"', from=1-1, to=3-4]
            """
tridiagram = TikzCD(triangle, preamble=Styles.Quiver)
testrender(tridiagram, "triangle")


text:
using SimilaritySearch
import StatsBase: fit
using Random
export dnet

"""
A `k`-net is a set of points `M` such that each object in `X` can be:
- It is in `M`
- It is in the knn set of an object in `M` (defined with the distance function `dist`)

The size of `M` is determined by ``\\leftceil|X|/k\\rightceil``

The dnet function uses the `callback` function as an output mechanism. This function is called on each center as `callback(centerId, res)` where
res is a `KnnResult` object (from SimilaritySearch.jl).

"""
function dnet(callback::Function, dist::Function, X::Vector{T}, k) where {T}
    N = length(X)
    metadist = (a::Int, b::Int) -> dist(X[a], X[b])
    I = fit(Sequential, shuffle!(collect(1:N)))
    res = KnnResult(k)
    i = 0

    while length(I.db) > 0
        i += 1
        empty!(res)
        c = pop!(I.db)
        search(I, metadist, c, res)
        callback(c, res)
        @info "computing dnet point $i, dmax: $(covrad(res))"

        j = 0
        for p in res
            I.db[p.objID] = I.db[end-j]
            j += 1
        end

        for p in res
            pop!(I.db)
        end
    end
end