text:
# =================================================================================================
# =                                         Hyland83.jl                                           =
# =================================================================================================
@testset "Hyland83" begin
# Appendix of ref. [1]. First table
@testset "Specific volume of saturated ice" begin
@test 1.0768e-3 ≈ MoistAir.volumeice(173.15) atol=0.0001e-3
@test 1.0829e-3 ≈ MoistAir.volumeice(223.15) atol=0.0001e-3
@test 1.0909e-3 ≈ MoistAir.volumeice(273.16) atol=0.0001e-3
end
@testset "Specific volume of saturated liquid water" begin
@test 1.00021e-3 ≈ MoistAir.volumewater(273.16) atol=0.00001e-3
@test 1.01215e-3 ≈ MoistAir.volumewater(323.15) atol=0.00001e-3
@test 1.04346e-3 ≈ MoistAir.volumewater(373.15) atol=0.00001e-3 # In the paper there is a typo in the power of 10: 1.04346e10
@test 1.09050e-3 ≈ MoistAir.volumewater(423.15) atol=0.00001e-3
@test 1.15653e-3 ≈ MoistAir.volumewater(473.15) atol=0.00001e-3
end

@testset "Saturation enthalpy of saturated ice" begin
@test -507.215e3 ≈ MoistAir.enthalpyice(173.15) atol=0.001e3
@test -429.413e3 ≈ MoistAir.enthalpyice(223.15) atol=0.001e3
@test -333.429e3 ≈ MoistAir.enthalpyice(273.15) atol=0.001e3  # The table presents -333.409. Typo???
end

@testset "Saturation enthalpy of saturated liquid water" begin
@test 0.0 ≈ MoistAir.enthalpywater(273.16) atol=1.0
@test 209.330e3 ≈ MoistAir.enthalpywater(323.15) atol=5
@test 419.158e3 ≈ MoistAir.enthalpywater(373.15) atol=5
@test 632.210e3 ≈ MoistAir.enthalpywater(423.15) atol=9
@test 852.329e3 ≈ MoistAir.enthalpywater(473.15) atol=15
end

# Table 2 of reference [2]
@testset "Enhancement factor" begin
P = 0.1e6
@test 1.0105 ≈ MoistAir.efactor(173.15, P) atol=0.0001
@test 1.0039 ≈ MoistAir.efactor(273.15, P) atol=0.0001
@test 1.0039 ≈ MoistAir.efactor(363.15, P) atol=0.0001
P = 0.5e6
@test 1.054 ≈ MoistAir.efactor(173.15, P) atol=0.001
@test 1.0177 ≈ MoistAir.efactor(273.15, P) atol=0.0006  # I think the table has a typo!!?
@test 1.0180 ≈ MoistAir.efactor(363.15, P) atol=0.0001
@test 1.0188 ≈ MoistAir.efactor(373.15, P) atol=0.0001
@test 1.0022 ≈ MoistAir.efactor(423.15, P) atol=0.0001
P = 1.0e6
@test 1.113 ≈ MoistAir.efactor(173.15, P) atol=0.001
@test 1.0353 ≈ MoistAir.efactor(273.15, P) atol=0.0015  # I think the table has a typo!!?
@test 1.0284 ≈ MoistAir.efactor(363.15, P) atol=0.0001
@test 1.0295 ≈ MoistAir.efactor(373.15, P) atol=0.0001
@test 1.0288 ≈ MoistAir.efactor(423.15, P) atol=0.0001
@test 1.0235 ≈ MoistAir.efactor(433.15, P) atol=0.0001
@test 1.0142 ≈ MoistAir.efactor(443.15, P) atol=0.0001

P = 5.0e6
@test 1.820 ≈ MoistAir.efactor(173.15, P) atol=0.001
@test 1.191 ≈ MoistAir.efactor(273.15, P) atol=0.008  # I think the table has a typo!!?
@test 1.1102 ≈ MoistAir.efactor(363.15, P) atol=0.0001
@test 1.1082 ≈ MoistAir.efactor(373.15, P) atol=0.0001
@test 1.111 ≈ MoistAir.efactor(423.15, P) atol=0.001
@test 1.114 ≈ MoistAir.efactor(433.15, P) atol=0.001
@test 1.116 ≈ MoistAir.efactor(443.15, P) atol=0.001
@test 1.117 ≈ MoistAir.efactor(453.15, P) atol=0.001
@test 1.116 ≈ MoistAir.efactor(473.15, P) atol=0.001
end

@testset "Saturation pressure of water vapor over ice" begin
# Second table in the appendix of reference [1]
@test 1.40510e-3 ≈ MoistAir.Pws_s(173.15) atol=0.00001e-3
@test 6.11153e2 ≈ MoistAir.Pws_s(273.15) atol=0.00001e2
end

@testset "Saturation pressure of water vapor over liquid water" begin
# Second table in the appendix of reference [1]
@test 6.11213e2 ≈ MoistAir.Pws_l(273.15) atol=0.00001e2
@test 1.01419e5 ≈ MoistAir.Pws_l(373.15) atol=0.00001e5
@test 1.55507e6 ≈ MoistAir.Pws_l(473.15) atol=0.00001e6
end

@testset " saturation temperature function" begin
@test MoistAir.Tws(MoistAir.Pws_s(173.15)) ≈ 173.15 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_s(223.15)) ≈ 223.15 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_s(273.15)) ≈ 273.15 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_s(273.16)) ≈ 273.16 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_l(273.16)) ≈ 273.16 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_l(323.15)) ≈ 323.15 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_l(373.15)) ≈ 373.15 atol=1e-5
@test MoistAir.Tws(MoistAir.Pws_l(473.15)) ≈ 473.15 atol=1e-5
end

@testset "First virial coefficient of saturated vapor B'" begin
# Second table in the appendix of reference [1]
@test -3.2939e-5 ≈ MoistAir.Blin(173.15) atol=0.0001e-5
@test -8.3497e-7 ≈ MoistAir.Blin(273.15) atol=0.0001e-7
@test -1.4658e-7 ≈ MoistAir.Blin(373.15) atol=0.0001e-5
@test -5.0508e-8 ≈ MoistAir.Blin(473.15) atol=0.0001e-5
end

@testset "Second virial coefficient of saturated vapor C'" begin
# Second table in the appendix of reference [1]
@test -4.6563e-9 ≈ MoistAir.Clin(173.15) atol=0.0001e-9
@test -2.0928e-12 ≈ MoistAir.Clin(273.15) atol=0.0001e-12
@test -5.7548e-14 ≈ MoistAir.Clin(373.15) atol=0.0001e-14
@test -6.3933e-15 ≈ MoistAir.Clin(473.15) atol=0.0001e-15
end

@testset "Specific enthalpy of saturated water vapor" begin
# First table of the appendix of ref. [1]
@test 2315.87 ≈ MoistAir.enthalpyvapor(173.15)/1000 atol=0.01
@test 2408.41 ≈ MoistAir.enthalpyvapor(223.15)/1000 atol=0.01
@test 2500.81 ≈ MoistAir.enthalpyvapor(273.16)/1000 atol=0.01
@test 2591.29 ≈ MoistAir.enthalpyvapor(323.15)/1000 atol=0.01
@test 2675.46 ≈ MoistAir.enthalpyvapor(373.15)/1000 atol=0.01
@test 2746.15 ≈ MoistAir.enthalpyvapor(423.15)/1000 atol=0.01
@test 2793.11 ≈ MoistAir.enthalpyvapor(473.15)/1000 atol=0.01
end

@testset "Specific volume of saturated water vapor" begin
# First table of the appendix of ref. [1]
@test 5.6873e7 ≈ MoistAir.volumevapor(173.15) atol=0.0001e7
@test 2.6146e4 ≈ MoistAir.volumevapor(223.15) atol=0.0001e4
@test 2.0601e2 ≈ MoistAir.volumevapor(273.16) atol=0.0001e2
@test 1.2030e1 ≈ MoistAir.volumevapor(323.15) atol=0.0001e1
@test 1.6718e0 ≈ MoistAir.volumevapor(373.15) atol=0.0001e0
@test 3.9253e-1 ≈ MoistAir.volumevapor(423.15) atol=0.0001e-1
@test 1.2722e-1 ≈ MoistAir.volumevapor(473.15) atol=0.0001e0
end

@testset "Specific entropy of saturated water vapor" begin
# First table of the appendix of ref. [1]
@test 14.30387 ≈ MoistAir.entropyvapor(173.15)/1000 atol=0.00005
@test 11.10965 ≈ MoistAir.entropyvapor(223.15)/1000 atol=0.00005
@test  9.15510 ≈ MoistAir.entropyvapor(273.16)/1000 atol=0.00005
@test  8.07477 ≈ MoistAir.entropyvapor(323.15)/1000 atol=0.00005
@test  7.35365 ≈ MoistAir.entropyvapor(373.15)/1000 atol=0.00005
@test  6.83731 ≈ MoistAir.entropyvapor(423.15)/1000 atol=0.00005
@test  6.43218 ≈ MoistAir.entropyvapor(473.15)/1000 atol=0.00005
end
end


text:
# # Introduciton to Lotka-Volterra Competition

using Parameters
using DifferentialEquations
using ForwardDiff
using Plots
pyplot();

# # Introduction/Overview
# This notebook is simply looking at the two species L-V competiti9on model. Here, we are exploring
# the geometry and dynamics of this model.
#
# # Model and Equilibria
# The model we will use is the classical LV competition model:

include("lv_model.jl")

# # Evaluationg the model -- time series

let
    u0 = [0.5, 0.5]
    t_span = (0.0, 100.0)
    p = LVPar()

    prob = ODEProblem(lv_comp, u0, t_span, p)
    sol = solve(prob, reltol = 1e-8)

    plot(sol)
    xlabel!("time")
    ylabel!("Density")
end

# # Equilibria and Isoclines
# From the above equations we can generate the Jacobian:

# Make a numerical (not symbolic version using the ForwardDiff.jl library.
# This is what you will want to do most of the time)
jac(u, p) = ForwardDiff.jacobian(u -> lv_comp(u, p, NaN), u)

# Here we have the jacobian evaluated at the point `u = [0.5, 0.5]`, with the parameter set `par`
jac([0.5, 0.5], LVPar())

# First solve for $f_1$ and $f_2$ to determine when functions are equal to 0 (i.e., when $du_1/dt$
# and $du_2/dt = 0$). These are known as the isoclines, or nullclines, and describe the set of
# solutions when $u_1$ and $u_2$ do not change.

include("lv_symbolic.jl")

# We want to solve our equations for $u_1$ so that we can plot the isoclines as a function of  $u_1$
# ($u_1$ on the y-axis), but here I am keeping all possible solutions, with respect to both $u_1$
# and $u_2$. You should be able to see that the first two solutions ($u_1 = K_1 - \alpha_{12}u_2$
# and $u_2 = (K_1 - u_1)/\alpha_{12}$) are equivalent.
#
# The solutions $u_1 = K_1 - \alpha_{12}u_2$ and $u_1 = (K_2 - u_2)/\alpha_{21}$ are the equations
# for our isoclines for $du_1/dt = 0$ and $du_2/dt = 0$, respectively. We can find the interior
# equilibrium where these lines intersect (i.e., both $u_1$ and $u_2$ do not change). In this case
# we get one interior equilibrium, but note that in more complicated models we can get more interior
# equilibria.

# # Now plot the isoclines
# We will again use the manipulate function to see how our parameters change the isoclines. This can
# immediately tell us a lot about our equilibria and stability.

function iso1(u2, p)
    @unpack α12, K1 = p
    return K1 - α12 * u2 / K1
end
#-
function iso2(u2, p)
    @unpack α21, K2 = p
    return (K2 - u2) / α21
end
#-
let
    p = LVPar()
    u2s = range(0, 5, length = 100)

    fig = plot(u2s, [iso1(u2, p) for u2 in u2s], label = "u1 = 0")
    plot!(u2s, [iso2(u2, p) for u2 in u2s], label = "u2 = 0")
    xlims!(0, 3)
    ylims!(0, 5)
    xlabel!("u2")
    ylabel!("u1")
end

ext: jl

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
# | Read a .2dm (SMS Aquaveo) mesh-file and construct a @Mesh@
function load(st::Stream{format"2DM"}; facetype=GLTriangleFace, pointtype=Point3f)
    faces = facetype[]
    vertices = pointtype[]
    io = stream(st)
    for line = readlines(io)
        if !isempty(line) && !all(iscntrl, line)
            line = chomp(line)
            w = split(line)
            if w[1] == "ND"
                push!(vertices, Point3f(parse.(Float32, w[3:end])))
            elseif w[1] == "E3T"
                push!(faces, GLTriangleFace(parse.(Cuint, w[3:5])))
            elseif w[1] == "E4Q"
                push!(faces, convert_simplex(facetype, QuadFace{Cuint}(parse.(Cuint, w[3:6])))...)
            else
                continue
            end
        end
    end
    return Mesh(vertices, faces)
end

function print_face(io::IO, i::Int, f::TriangleFace)
    println(io, "E3T $i ", join(Int.(value.(f)), " "))
end

function print_face(io::IO, i::Int, f::QuadFace)
    println(io, "E4Q $i ", join(Int.(value.(f)), " "))
end

# | Write @Mesh@ to an IOStream
function save(st::Stream{format"2DM"}, m::AbstractMesh)
    io = stream(st)
    println(io, "MESH2D")
    for (i, f) in enumerate(faces(m))
        print_face(io, i, f)
    end
    for (i, v) in enumerate(coordinates(m))
        println(io, "ND $i ", join(v, " "))
    end
    return
end

show(io::IO, ::MIME"model/2dm", mesh::AbstractMesh) = save(io, mesh)


text:
"""
    Monomials(n) = [ x-> x^k for k = 0:n-1 ]
"""
struct Monomials <: AbstractBasis
    n::Int
end

Base.length(B::Monomials) = B.n

approxpoints(B::Monomials) = exp.(2π*im/length(B).*(0:length(B)-1))
approxtransform(B::Monomials,f) = fft(convert(Array,f),1)./length(B)

function iterate_basis(B::Monomials, x)
    length(B) == 0 && return nothing
    p = one(x)
    return p,(2,p)
end
function iterate_basis(B::Monomials, x, (i,p))
    i > length(B) && return nothing
    p = x*p
    return p,(i+1,p)
end

evaluate_basis(B::Monomials,i,x) = x^(i-1)


text:
# ------------------------------------------------------------------
# It is advantageous to have the coeffs as
# one single vector in order to make use of BLAS
# and related libraries when defining operators on
# our space of functions
#
# In this script, we "vectorize" our functions,
# not in the regular meaning of the word, but in the
# sense that dictionaries will all be recelld by vectors
#
# The main difficulty with this is that, unlike with dictionaries
# there is no easy way to go from (level, cell, mode) to a
# corresponding index in a 1-D vector.
#
# For this, we have the reference functions:
# full/sparse_referenceV2D/D2V
# D2V generates a dict that, upon input of a level, cell, mode
# gives the corresponding index in a vector
# V2D generates a vector, with row i having the three numbers
# level, cell, mode corresponding to index i in the vector
#
# These methods work in all dimensions, and there are ones for both
# full and sparse grids.
#
# The entire script culminates in a final result: a matrix
# representation of the derivative operator
# (both in full and sparse bases)
# ------------------------------------------------------------------

# Calculates the exact dimension of interpolating
# basis functions in the full or sparse grid schemes
# using the degree k-1 Galerkin polynomials
# with multi-resolution up to n
function get_size(D::Int, k::Int, n::Int; scheme="sparse")
	cutoff	= get_cutoff(scheme, D, n)
	ls		= ntuple(i->(n+1), D)
	size	= 0
	for level in CartesianRange(ls)
		cutoff(level) && continue

		ks = ntuple(q -> 1<<max(0, level[q]-2), D)
		size += prod(ks)*k^D
	end
	return size
end

function D2V{d,T<:Real}(D::Int, k::Int, n::Int,
						coeffs::Dict{CartesianIndex{d}, Array{Array{T,d},d}};
						scheme="sparse")
	(d != D) && throw(TypeError(:coeffs))
	cutoff		= get_cutoff(scheme, D, n)
	size		= get_size(D, k, n; scheme=scheme)
	vect		= Array{T}(size)
	modes	= ntuple(i-> k, D)
	ls			= ntuple(i->(n+1), D)
	j = 1
	for level in CartesianRange(ls) #This really goes from 0 to l_i for each i
		cutoff(level) && continue

		ks = ntuple(q -> 1<<max(0, level[q]-2), D)  #This sets up a specific k+1 vector
		for cell in CartesianRange(ks)
			for mode in CartesianRange(modes)
				vect[j] = coeffs[level][cell][mode]
				j += 1
			end
		end
	end
	return vect
end

function V2D{T<:Real}(D::Int, k::Int, n::Int, vect::Array{T}; scheme="sparse")
	cutoff		= get_cutoff(scheme, D, n)
	coeffs		= Dict{CartesianIndex{D}, Array{Array{T,D},D}}()
	modes		= ntuple(q-> k, D)
	ls			= ntuple(i->(n+1), D)
	j = 1
	for level in CartesianRange(ls) #This really goes from 0 to l_i for each i
		cutoff(level) && continue

		ks = ntuple(q -> 1<<max(0, level[q]-2), D)  #This sets up a specific k+1 vector
		level_coeffs = Array{Array{T}}(ks) #all the coefficients at this level
		for cell in CartesianRange(ks)
			cell_coeffs = Array{T}(modes)
			for mode in CartesianRange(modes)
				cell_coeffs[mode] = vect[j]
				j += 1
			end
			level_coeffs[cell] = cell_coeffs
		end
		coeffs[level] = level_coeffs
	end
	return coeffs
end

function D2Vref(D::Int, k::Int, n::Int; scheme="sparse")
	cutoff		= get_cutoff(scheme, D, n)
	size		= get_size(D, k, n; scheme=scheme)
	dict		= Dict{NTuple{3,CartesianIndex{D}}, Int}()
	modes		= ntuple(q-> k, D)
	ls			= ntuple(i->(n+1), D)
	j = 1
	for level in CartesianRange(ls)
		cutoff(level) && continue

		ks = ntuple(q -> 1<<max(0, level[q]-2), D)  #This sets up a specific k+1 vector
		lvl = ntuple(i -> level[i]-1,D)
		for cell in CartesianRange(ks)
			for mode in CartesianRange(modes)
				dict[(level, cell, mode)] = j
				j += 1
			end
		end
	end
	return dict
end

function V2Dref(D::Int, k::Int, n::Int; scheme = "sparse")
	cutoff		= get_cutoff(scheme, D, n)
	size		= get_size(D, k, n; scheme=scheme)
	vect		= Array{NTuple{3,CartesianIndex{D}}}(size)
	modes	= ntuple(q-> k, D)
	ls			= ntuple(i->(n+1), D)
	j = 1
	for level in CartesianRange(ls)
		cutoff(level) && continue

		ks = ntuple(q -> 1<<max(0, level[q]-2), D)  #This sets up a specific k+1 vector
		lvl = ntuple(i -> level[i]-1,D)
		for cell in CartesianRange(ks)
			for mode in CartesianRange(modes)
				vect[j] = (level, cell, mode)
				j += 1
			end
		end
	end
	return vect
end

# -----------------------------------------------------
# Let's now make the coefficient operators work on
# and return vectors
# -----------------------------------------------------
function vcoeffs_DG(D::Int, k::Int, n::Int, f::Function;
								rel_tol=REL_TOL, abs_tol=ABS_TOL,
								max_evals=MAX_EVALS,
								scheme="sparse")
	cutoff		= get_cutoff(scheme, D, n)
	len			= get_size(D, k, n; scheme=scheme)
	coeffs		= Array{Float64}(len)
	modes		= ntuple(i-> k, D)
	ls			= ntuple(i-> (n+1), D)
	j = 1
	for level in CartesianRange(ls)     # This really goes from 0 to l_i for each i,
		cutoff(level) && continue

		ks = ntuple(i -> 1<<max(0, level[i]-2), D)  #This sets up a specific k+1 vector
		lvl = ntuple(i -> level[i]-1,D)
		for cell in CartesianRange(ks)
			for mode in CartesianRange(modes)
				coeffs[j] = get_coefficient_DG(k, lvl, cell, mode, f;
												rel_tol=rel_tol, abs_tol=abs_tol,
												max_evals=max_evals)
				j += 1
			end
		end
	end
	return coeffs
end


text:
using Base.Test
include("utils.jl")

@testset "WorkUnit Tests" begin

@testset "WorkUnit constructor tests"  begin
    @test_throws ErrorException WorkUnit(-1, 2)
    @test_throws ErrorException WorkUnit(2, 0)
    @test_throws ErrorException WorkUnit(0, 0)
    @test typeof(WorkUnit(2, 3)) == WorkUnit{typeof(2), typeof(3)}
end

@testset "split_work tests" begin
    testwork = WorkUnit(10, 3)
    testwork2 = WorkUnit(555, 4)
    testwork3 = WorkUnit(1, 1)

    @testset "split_work domain tests" begin
        @test_throws DomainError split_work(testwork, -0.5)
        @test_throws DomainError split_work(testwork2, 1.63)
        @test_throws DomainError split_work(testwork3, .5)
    end

    @testset "split_work 50% test" begin
        u1, v1 = split_work(testwork, 0.5)
        u2, v2 = split_work(testwork2, 0.5)
        @test u1.units == 5
        @test v1.units == 5
        @test u2.units == 278
        @test v2.units == 277
    end

    @testset "split_work 0 + ϵ test" begin
        u1, v1 = split_work(testwork, eps(Float64))
        u2, v2 = split_work(testwork2, eps(Float64))
        @test u1.units == 1
        @test v1.units == 9
        @test u2.units == 1
        @test v2.units == 554
    end

    @testset "split_work 1 - ϵ test" begin
        u1, v1 = split_work(testwork, 1 - eps(Float64))
        u2, v2 = split_work(testwork2, 1 - eps(Float64))
        @test u1.units == 9
        @test v1.units == 1
        @test u2.units == 554
        @test v2.units == 1
    end

end

end

@testset "Utility Tests" begin

    @testset "Single process" begin
        @test w_idx(myid()) == myid()
        let msg_chls = create_msg_chls(16)
            @test length(msg_chls) == 1
            @test msg_chls[1].whence == myid()
            @test msg_chls[1].where == myid()
            @test msg_chls[w_idx(myid())] == msg_chls[1]
        end
    end

    @testset "Multiprocess" begin
        @testset "1 Controller, 1 Worker Tests" begin
            addprocs(1)
            @spawnat 2 include("utils.jl")
            @test w_idx(myid()) == 0
            let msg_chls = create_msg_chls(16)
                @test length(msg_chls) == 1
                @test msg_chls[1].whence == myid()
                @test msg_chls[1].where == workers()[1]
                @test_throws BoundsError msg_chls[w_idx(myid())]
                #
                @test remotecall_fetch(length, workers()[1], msg_chls) == 1
                @test_broken remotecall_fetch((y)->(x = get_msg_chl(myid(), y); x.where),
                                       workers()[1],
                                       msg_chls) == workers()[1]
            end
            rmprocs(workers()[1])
        end

    end

end

max_stars_repo_name: dillondaudert/jlance

max_issues_repo_name: dillondaudert/jlance

max_forks_repo_name: dillondaudert/jlance


text:
"""
The `mode` type parameter must be either `:inflate` or `:deflate`.
"""
type Sink{mode,T<:BufferedOutputStream}
    output::T
    zstream::ZStream
    state::State
end

# inflate sink constructors
# -------------------------

function InflateSink{T<:BufferedOutputStream}(output::T, gzip::Bool)
    zstream = init_inflate_zstream(gzip)
    zstream.next_out = pointer(output)
    zstream.avail_out = BufferedStreams.available_bytes(output)
    return Sink{:inflate,T}(output, zstream, initialized)
end

function InflateSink(output::BufferedOutputStream, bufsize::Integer, gzip::Bool)
    return InflateSink(output, gzip)
end

function InflateSink(output::IO, bufsize::Integer, gzip::Bool)
    output_stream = BufferedOutputStream(output, bufsize)
    return InflateSink(output_stream, gzip)
end

function InflateSink(output::Vector{UInt8}, bufsize::Integer, gzip::Bool)
    return InflateSink(BufferedOutputStream(output), gzip)
end

"""
    ZlibInflateOutputStream(output[; <keyword arguments>])

Construct a zlib inflate output stream to decompress gzip/zlib data.

# Arguments
* `output`: a byte vector, IO object, or BufferedInputStream to which decompressed data should be written.
* `bufsize::Integer=8192`: input and output buffer size.
* `gzip::Bool=true`: if true, data is gzip compressed; if false, zlib compressed.
"""
function ZlibInflateOutputStream(output; bufsize::Integer=8192, gzip::Bool=true)
    return BufferedOutputStream(InflateSink(output, bufsize, gzip), bufsize)
end

# deflate sink constructors
# -------------------------

function DeflateSink{T<:BufferedOutputStream}(
        output::T, gzip::Bool, level::Integer, mem_level::Integer, strategy)
    zstream = init_deflate_zstream(gzip, level, mem_level, strategy)
    zstream.next_out = pointer(output)
    zstream.avail_out = BufferedStreams.available_bytes(output)
    return Sink{:deflate,T}(output, zstream, initialized)
end

function DeflateSink(output::BufferedOutputStream, bufsize::Integer, gzip::Bool,
                     level::Integer, mem_level::Integer, strategy)
    return DeflateSink(output, gzip, level, mem_level, strategy)
end

function DeflateSink(output::IO, bufsize::Integer, gzip::Bool, level::Integer,
                     mem_level::Integer, strategy)
    output_stream = BufferedOutputStream(output, bufsize)
    return DeflateSink(output_stream, gzip, level, mem_level, strategy)
end

function DeflateSink(output::Vector{UInt8}, bufsize::Integer, gzip::Bool,
                     level::Integer, mem_level::Integer, strategy)
    return DeflateSink(BufferedOutputStream(output), gzip, level, mem_level, strategy)
end

"""
    ZlibDeflateOutputStream(output[; <keyword arguments>])

Construct a zlib deflate output stream to compress gzip/zlib data.

# Arguments
* `output`: a byte vector, IO object, or BufferedInputStream to which compressed data should be written.
* `bufsize::Integer=8192`: input and output buffer size.
* `gzip::Bool=true`: if true, data is gzip compressed; if false, zlib compressed.
* `level::Integer=6`: compression level in 1-9.
* `mem_level::Integer=8`: memory to use for compression in 1-9.
* `strategy=Z_DEFAULT_STRATEGY`: compression strategy; see zlib documentation.
"""
function ZlibDeflateOutputStream(output;
                                 bufsize::Integer=8192,
                                 gzip::Bool=true,
                                 level::Integer=6,
                                 mem_level::Integer=8,
                                 strategy=Z_DEFAULT_STRATEGY)
    return BufferedOutputStream(
        DeflateSink(output, bufsize, gzip, level, mem_level, strategy),
        bufsize)
end

"""
    writebytes(sink, buffer, n, eof)

Write some bytes from a given buffer. Satisfies the BufferedStreams sink
interface.
"""
function BufferedStreams.writebytes{mode}(
        sink::Sink{mode},
        buffer::Vector{UInt8},
        n::Int, eof::Bool)
    if sink.state == finalized
        return 0
    elseif sink.state == finished
        reset!(sink)
    end

    @trans sink (
        initialized => inprogress,
        inprogress  => inprogress
    )

    BufferedStreams.flushbuffer!(sink.output)
    sink.zstream.next_in = pointer(buffer)
    sink.zstream.avail_in = n
    n_in, _ = process(sink, mode == :deflate && eof ? Z_FINISH : Z_NO_FLUSH)
    return n_in
end

function Base.flush(sink::Sink)
    if sink.state == finalized
        return
    elseif sink.state == inprogress
        process(sink, Z_FINISH)
    end
    flush(sink.output)
    return
end

function process{mode}(sink::Sink{mode}, flush)
    @assert sink.state == inprogress
    # counter of processed input/output bytes
    n_in = n_out = 0
    output = sink.output
    zstream = sink.zstream

    #println("--- Sink{", mode, "} ---")
    @label process
    zstream.next_out = pointer(output)
    zstream.avail_out = BufferedStreams.available_bytes(output)
    old_avail_in = zstream.avail_in
    old_avail_out = zstream.avail_out
    if mode == :inflate
        ret = inflate!(zstream, flush)
    else
        ret = deflate!(zstream, flush)
    end
    n_in += old_avail_in - zstream.avail_in
    n_out += old_avail_out - zstream.avail_out
    output.position += old_avail_out - zstream.avail_out

    if ret == Z_OK
        if zstream.avail_out == 0
            BufferedStreams.flushbuffer!(output)
            @goto process
        end
    elseif ret == Z_STREAM_END
        @trans sink inprogress => finished
    elseif ret == Z_BUF_ERROR
        # could not consume more input or produce more output
    elseif ret < 0
        zerror(zstream, ret)
    else
        @assert false
    end

    return n_in, n_out
end

function Base.close{mode}(sink::Sink{mode})
    if sink.state == finalized
        isopen(sink.output) && close(sink.output)
        return
    end
    if mode == :inflate
        @zcheck end_inflate!(sink.zstream)
    else
        @zcheck end_deflate!(sink.zstream)
    end
    @trans sink (
        initialized => finalized,
        inprogress  => finalized,
        finished    => finalized
    )
    close(sink.output)
    return
end

function reset!{mode}(sink::Sink{mode})
    if mode == :inflate
        @zcheck reset_inflate!(sink.zstream)
    else
        @zcheck reset_deflate!(sink.zstream)
    end
    @trans sink (
        initialized => initialized,
        inprogress  => initialized,
        finished    => initialized,
        finalized   => initialized
    )
    return sink
end


text:
export AtmosFilterPerturbations

import ..Mesh.Filters:
    AbstractFilterTarget,
    vars_state_filtered,
    compute_filter_argument!,
    compute_filter_result!

struct AtmosFilterPerturbations{M} <: AbstractFilterTarget
    atmos::M
end

vars_state_filtered(target::AtmosFilterPerturbations, FT) =
    vars_state_conservative(target.atmos, FT)

function compute_filter_argument!(
    ::AtmosFilterPerturbations,
    filter_state::Vars,
    state::Vars,
    aux::Vars,
)
    # copy the whole state
    parent(filter_state) .= parent(state)
    # remove reference state
    filter_state.ρ -= aux.ref_state.ρ
    filter_state.ρe -= aux.ref_state.ρe
end
function compute_filter_result!(
    ::AtmosFilterPerturbations,
    state::Vars,
    filter_state::Vars,
    aux::Vars,
)
    # copy the whole filter state
    parent(state) .= parent(filter_state)
    # add reference state
    state.ρ += aux.ref_state.ρ
    state.ρe += aux.ref_state.ρe
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
#TODO: Document interface to support static size
# Define maybestaticsize, maybestaticlength, and maybestaticfirstindex

@inline static(::Val{N}) where {N} = StaticInt{N}()
@inline static(::Nothing) = nothing
@generated function static_sizeof(::Type{T}) where {T}
    Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:curly, :StaticInt, sizeof(T))))
end

# @inline static_last(::Type{T}) where {T} = static(known_last(T))
# @inline static_first(::Type{T}) where {T} = static(known_first(T))
# @inline static_length(::Type{T}) where {T} = static(known_length(T))

@inline maybestaticfirst(a) = static_first(a)
@inline maybestaticlast(a) = static_last(a)
@inline maybestaticlength(a) = static_length(a)

# @inline _maybestaticfirst(a, ::Nothing) = first(a)
# @inline _maybestaticfirst(::Any, L) = StaticInt{L}()
# @inline maybestaticfirst(a::T) where {T} = _maybestaticfirst(a, known_first(T))

# @inline _maybestaticlast(a, ::Nothing) = last(a)
# @inline _maybestaticlast(::Any, L) = StaticInt{L}()
# @inline maybestaticlast(a::T) where {T} = _maybestaticlast(a, known_last(T))

# @inline _maybestaticlength(a, L::Int, U::Int) = StaticInt{U-L}()
# @inline _maybestaticlength(a, ::Any, ::Any) = length(a)
# @inline _maybestaticlength(a, ::Nothing) = _maybestaticlength(a, known_first(a), known_last(a))
# @inline _maybestaticlength(::Any, L::Int) = StaticInt{L}()
# @inline maybestaticlength(a::T) where {T} = _maybestaticlength(a, known_length(T))

@inline maybestaticrange(r::Base.OneTo{T}) where {T} = ArrayInterface.OptionallyStaticUnitRange(StaticInt{1}(), last(r))
@inline maybestaticrange(r::UnitRange) = r
# @inline maybestaticrange(r::AbstractStaticIntUnitRange) = r
@inline maybestaticrange(r) = maybestaticfirst(r):maybestaticlast(r)
# @inline maybestaticaxis(A::AbstractArray, ::Val{I}) where {I} = maybestaticfirstindex(A, Val{I}()):maybestaticsize(A, Val{I}())

@inline maybestaticsize(::NTuple{N}, ::Val{1}) where {N} = StaticInt{N}() # should we assert that i == 1?
@inline maybestaticlength(::NTuple{N}) where {N} = StaticInt{N}()
@inline maybestaticsize(::LinearAlgebra.Adjoint{T,V}, ::Val{1}) where {T,V<:AbstractVector{T}} = One()
@inline maybestaticsize(::LinearAlgebra.Transpose{T,V}, ::Val{1}) where {T,V<:AbstractVector{T}} = One()
@inline maybestaticlength(B::LinearAlgebra.Adjoint) = maybestaticlength(parent(B))
@inline maybestaticlength(B::LinearAlgebra.Transpose) = maybestaticlength(parent(B))
@inline maybestaticsize(A, ::Val{N}) where {N} = ArrayInterface.size(A)[N]
# @inline maybestaticsize(B::LinearAlgebra.Adjoint{T,A}, ::Val{1}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{2}())
# @inline maybestaticsize(B::LinearAlgebra.Adjoint{T,A}, ::Val{2}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{1}())
# @inline maybestaticsize(B::LinearAlgebra.Transpose{T,A}, ::Val{1}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{2}())
# @inline maybestaticsize(B::LinearAlgebra.Transpose{T,A}, ::Val{2}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{1}())

# These have versions that may allow for more optimizations, so we override base methods with a single `StaticInt` argument.
for (f,ff) ∈ [(:(Base.:+),:vadd), (:(Base.:-),:vsub), (:(Base.:*),:vmul), (:(Base.:<<),:vshl), (:(Base.:÷),:vdiv), (:(Base.:%), :vrem), (:(Base.:>>>),:vashr)]
    @eval begin
        # @inline $f(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{$f(M, N)}()
    # If `M` and `N` are known at compile time, there's no need to add nsw/nuw flags.
        @inline $ff(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{$f(M, N)}()
        # @inline $f(::StaticInt{M}, x) where {M} = $ff(M, x)
        @inline $ff(::StaticInt{M}, x) where {M} = $ff(M, x)
        # @inline $f(x, ::StaticInt{M}) where {M} = $ff(x, M)
        @inline $ff(x, ::StaticInt{M}) where {M} = $ff(x, M)
    end
end
for f ∈ [:vadd, :vsub, :vmul]
    @eval @inline $f(::StaticInt{M}, n::Number) where {M} = $f(M, n)
    @eval @inline $f(m::Number, ::StaticInt{N}) where {N} = $f(m, N)
end

@inline vsub(::StaticInt{N}, ::Zero) where {N} = StaticInt{N}()
# @inline vsub(::Zero, ::StaticInt{N}) where {N} = StaticInt{-N}()
@inline vsub(::Zero, ::Zero) = Zero()
@inline vsub(a::Number, ::Zero) = a
@inline vsub(a, ::Zero) = a

@inline vadd(::StaticInt{N}, ::Zero) where {N} = StaticInt{N}()
@inline vadd(::Zero, ::StaticInt{N}) where {N} = StaticInt{N}()
@inline vadd(::Zero, ::Zero) = Zero()
@inline vadd(a::Number, ::Zero) = a
@inline vadd(::Zero, a::Number) = a

@inline vmul(::StaticInt{N}, ::Zero) where {N} = Zero()
@inline vmul(::Zero, ::StaticInt{N}) where {N} = Zero()
@inline vmul(::Zero, ::Zero) = Zero()
@inline vmul(::StaticInt{N}, ::One) where {N} = StaticInt{N}()
@inline vmul(::One, ::StaticInt{N}) where {N} = StaticInt{N}()
@inline vmul(::One, ::One) = One()
@inline vmul(a::Number, ::One) = a
@inline vmul(::One, a::Number) = a
@inline vmul(::Zero, ::One) = Zero()
@inline vmul(::One, ::Zero) = Zero()
@inline vmul(i::MM{W,X}, ::StaticInt{N}) where {W,X,N} = MM{W}(vmul(data(i), StaticInt{N}()), StaticInt{X}() * StaticInt{N}())
@inline vmul(i::MM{W,X}, ::StaticInt{1}) where {W,X} = i
@inline vmul(::StaticInt{N}, i::MM{W,X}) where {W,X,N} = MM{W}(vmul(data(i), StaticInt{N}()), StaticInt{X}() * StaticInt{N}())
@inline vmul(::StaticInt{1}, i::MM{W,X}) where {W,X} = i

@generated staticp1(::StaticInt{N}) where {N} = StaticInt{N+1}()
@inline staticp1(N) = vadd(N, one(N))
@inline staticp1(i::Tuple{}) = tuple()
@inline staticp1(i::Tuple{I}) where {I} = @inbounds (staticp1(i[1]),)
@inline staticp1(i::Tuple{I1,I2}) where {I1,I2} = @inbounds (staticp1(i[1]), staticp1(i[2]))
@inline staticp1(i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = @inbounds (staticp1(i[1]), staticp1(Base.tail(i))...)
@generated staticm1(::StaticInt{N}) where {N} = StaticInt{N-1}()
@inline staticm1(N) = vsub(N, one(N))
@inline staticm1(i::Tuple{}) = tuple()
@inline staticm1(i::Tuple{I}) where {I} = @inbounds (staticm1(i[1]),)
@inline staticm1(i::Tuple{I1,I2}) where {I1,I2} = @inbounds (staticm1(i[1]), staticm1(i[2]))
@inline staticm1(i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = @inbounds (staticm1(i[1]), staticm1(Base.tail(i))...)
@generated staticmul(::Type{T}, ::StaticInt{N}) where {T,N} = sizeof(T) * N
@generated staticmul(::Type{T}, ::Val{N}) where {T,N} = sizeof(T) * N
@inline staticmul(::Type{T}, N) where {T} = vmul(N, sizeof(T))
@inline staticmul(::Type{T}, i::Tuple{}) where {T} = tuple()
@inline staticmul(::Type{T}, i::Tuple{I}) where {T,I} = @inbounds (vmul(i[1], sizeof(T)),)
@inline staticmul(::Type{T}, i::Tuple{I1,I2}) where {T,I1,I2} = @inbounds (vmul(sizeof(T), i[1]), vmul(sizeof(T), i[2]))
@inline staticmul(::Type{T}, i::Tuple{I1,I2,I3,Vararg}) where {T,I1,I2,I3} = @inbounds (vmul(sizeof(T), i[1]), staticmul(T, Base.tail(i))...)
for T ∈ [:VecUnroll, :Mask, :MM]
    @eval begin
        @inline Base.:(+)(x::$T, ::Zero) = x
        @inline Base.:(+)(::Zero, x::$T) = x
        @inline Base.:(-)(x::$T, ::Zero) = x
        @inline Base.:(*)(x::$T, ::One) = x
        @inline Base.:(*)(::One, x::$T) = x
        @inline Base.:(*)(::$T, ::Zero) = Zero()
        @inline Base.:(*)(::Zero, ::$T) = Zero()
    end
end
@inline Base.:(+)(m::Mask{W}, ::StaticInt{N}) where {N,W} = m + vbroadcast(Val{W}(), N)
@inline Base.:(+)(::StaticInt{N}, m::Mask{W}) where {N,W} = vbroadcast(Val{W}(), N) + m
# @inline Base.:(*)(::StaticInt{N}, m::Mask{W}) where {N,W} = vbroadcast(Val{W}(), N) * m
@inline vadd(x::VecUnroll, ::Zero) = x
@inline vadd(::Zero, x::VecUnroll) = x
@inline vsub(x::VecUnroll, ::Zero) = x
@inline vmul(x::VecUnroll, ::One) = x
@inline vmul(::One, x::VecUnroll) = x
@inline vmul(::VecUnroll, ::Zero) = Zero()
@inline vmul(::Zero, ::VecUnroll) = Zero()


text:
# ---
# title: 972. Equal Rational Numbers
# id: problem972
# author: Indigo
# date: 2021-02-20
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/equal-rational-numbers/description/>
# hidden: true
# ---
#
# Given two strings `S` and `T`, each of which represents a non-negative
# rational number, return **True** if and only if they represent the same
# number. The strings may use parentheses to denote the repeating part of the
# rational number.
#
# In general a rational number can be represented using up to three parts: an
# _integer part_ , a  _non-repeating part,_ and a  _repeating part_. The number
# will be represented in one of the following three ways:
#
#   * `<IntegerPart>` (e.g. 0, 12, 123)
#   * `<IntegerPart><.><NonRepeatingPart>`  (e.g. 0.5, 1., 2.12, 2.0001)
#   * `<IntegerPart><.><NonRepeatingPart><(><RepeatingPart><)>` (e.g. 0.1(6), 0.9(9), 0.00(1212))
#
# The repeating portion of a decimal expansion is conventionally denoted within
# a pair of round brackets.  For example:
#
# 1 / 6 = 0.16666666... = 0.1(6) = 0.1666(6) = 0.166(66)
#
# Both 0.1(6) or 0.1666(6) or 0.166(66) are correct representations of 1 / 6.
#
#
#
# **Example 1:**
#
#
#
#     Input: S = "0.(52)", T = "0.5(25)"
#     Output: true
#     Explanation: Because "0.(52)" represents 0.52525252..., and "0.5(25)" represents 0.52525252525..... , the strings represent the same number.
#
#
# **Example 2:**
#
#
#
#     Input: S = "0.1666(6)", T = "0.166(66)"
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: S = "0.9(9)", T = "1."
#     Output: true
#     Explanation:
#     "0.9(9)" represents 0.999999999... repeated forever, which equals 1.  [[See this link for an explanation.](https://en.wikipedia.org/wiki/0.999...)]
#     "1." represents the number 1, which is formed correctly: (IntegerPart) = "1" and (NonRepeatingPart) = "".
#
#
#
# **Note:**
#
#   1. Each part consists only of digits.
#   2. The `<IntegerPart>` will not begin with 2 or more zeros.  (There is no other restriction on the digits of each part.)
#   3. `1 <= <IntegerPart>.length <= 4 `
#   4. `0 <= <NonRepeatingPart>.length <= 4 `
#   5. `1 <= <RepeatingPart>.length <= 4`
#
#
## @lc code=start
using LeetCode

function is_rational_equal(s::String, t::String)
    function to_rational(s::String)::Rational
        pos1, pos2 = findfirst('.', s), findfirst('(', s)
        (pos1 === nothing) && return parse(Int, s)
        (pos1 == length(s)) && return parse(Int, s[1:end-1])
        pos2 === nothing && return parse(Int, s[1:pos1-1]) + parse(Int, s[pos1 + 1:end]) // 10 ^ (length(s) - pos1)
        pos3 = findlast(')', s)
        int_part = parse(Int, s[1:pos1-1]) |> Rational
        nr_part = pos1 == pos2 - 1 ? 0 : parse(Int, s[pos1+1:pos2-1]) // 10 ^ (pos2 - pos1 - 1)
        rp_part = parse(Int, s[pos2+1:pos3-1]) // (10 ^ (pos3 - pos2 - 1) - 1) // 10 ^ (pos2 - pos1 - 1)
        int_part + nr_part + rp_part
    end
    to_rational(s) == to_rational(t)
end
## @lc code=end
is_rational_equal("0.9(9)", "1.")


text:
module Coordinates

using Base: @_propagate_inbounds_meta, @_inline_meta
using Base.Cartesian: @ntuple

export Coordinate, coordinateaxes

"""
    Coordinate(axes...)
    Coordinate{N}(axis)

Construct `Coordinate` from `axes`.
If only single `axis` is given with number of dimensions `N`,
the `axis` is used for all dimensions.

# Examples
```jldoctest
julia> Coordinate(1:3, 2:4)
3×3 Coordinate{2,Tuple{Int64,Int64},Tuple{UnitRange{Int64},UnitRange{Int64}}}:
 (1, 2)  (1, 3)  (1, 4)
 (2, 2)  (2, 3)  (2, 4)
 (3, 2)  (3, 3)  (3, 4)

julia> Coordinate(1:2, [3.0,4.0], "abc")
2×2×3 Coordinate{3,Tuple{Int64,Float64,Char},Tuple{UnitRange{Int64},Array{Float64,1},String}}:
[:, :, 1] =
 (1, 3.0, 'a')  (1, 4.0, 'a')
 (2, 3.0, 'a')  (2, 4.0, 'a')

[:, :, 2] =
 (1, 3.0, 'b')  (1, 4.0, 'b')
 (2, 3.0, 'b')  (2, 4.0, 'b')

[:, :, 3] =
 (1, 3.0, 'c')  (1, 4.0, 'c')
 (2, 3.0, 'c')  (2, 4.0, 'c')

julia> Coordinate{2}(0:3)
4×4 Coordinate{2,Tuple{Int64,Int64},Tuple{UnitRange{Int64},UnitRange{Int64}}}:
 (0, 0)  (0, 1)  (0, 2)  (0, 3)
 (1, 0)  (1, 1)  (1, 2)  (1, 3)
 (2, 0)  (2, 1)  (2, 2)  (2, 3)
 (3, 0)  (3, 1)  (3, 2)  (3, 3)
```
"""
struct Coordinate{N, ElType, Axes <: Tuple{Vararg{Any, N}}} <: AbstractArray{ElType, N}
    axes::Axes
    function Coordinate{N, ElType, Axes}(axes::Axes) where {N, ElType, Axes}
        check_coordinate_parameter(Val(N), ElType, Axes)
        new{N, ElType, Axes}(axes)
    end
end

@generated function check_coordinate_parameter(::Val{N}, ::Type{ElType}, ::Type{Axes}) where {N, ElType, Axes}
    if Tuple{map(eltype, Axes.parameters)...} != ElType
        return :(throw(ArgumentError("`ElType` and `eltype`s of given axes must match.")))
    end
end

# constructors
@generated function Coordinate(axes::Tuple{Vararg{Any, N}}) where {N}
    ElType = Tuple{map(eltype, axes.parameters)...}
    quote
        Coordinate{N, $ElType, $axes}(axes)
    end
end
Coordinate(axes...) = Coordinate(axes)
Coordinate{N}(ax) where {N} = Coordinate(ntuple(i -> ax, Val(N)))

coordinateaxes(C::Coordinate) = C.axes
coordinateaxes(C::Coordinate, i::Int) = (@_propagate_inbounds_meta; coordinateaxes(C)[i])

Base.size(A::Coordinate) = map(length, coordinateaxes(A))

_eachindex(x) = firstindex(x):lastindex(x)
_eachindex(x::Union{AbstractArray, Tuple}) = eachindex(x)
Base.axes(A::Coordinate) = map(_eachindex, coordinateaxes(A))

# getindex
@generated function Base.getindex(A::Coordinate{N}, I::Vararg{Int, N}) where {N}
    quote
        @_inline_meta
        @boundscheck checkbounds(A, I...)
        @inbounds @ntuple $N i -> coordinateaxes(A, i)[I[i]]
    end
end
@generated function Base.getindex(A::Coordinate{N}, I::Vararg{Union{Int, AbstractUnitRange, Colon}, N}) where {N}
    quote
        @_inline_meta
        @boundscheck checkbounds(A, I...)
        @inbounds Coordinate(@ntuple $N i -> coordinateaxes(A, i)[I[i]])
    end
end
@inline function Base.getindex(A::Coordinate{N}, I::CartesianIndices{N}) where {N}
    @boundscheck checkbounds(A, I)
    @inbounds A[I.indices...]
end

Base.transpose(c::Coordinate{2}) = Coordinate(coordinateaxes(c, 2), coordinateaxes(c, 1))
Base.adjoint(c::Coordinate{2}) = transpose(c)

end # module


text:
using Test
using Quiqbox

@testset "TwoBody.jl tests" begin
    errT = 1e-10
    nucs = ["H", "H"]
    cens = [[-0.7, 0.0, 0.0], [ 0.7, 0.0, 0.0]]
    bs = genBasisFunc.(cens, fill("6-31G", 2)) |> flatten
    bsiz = basisSize(bs) |> sum
    vec = [[1.1401448653048112 0.6114018855091738 0.38206058180946545 0.4478418840074861;
            0.6114018855091738 0.626369765524018 0.34091363711214673 0.4981842896930728;
            0.38206058180946545 0.34091363711214673 0.6748720646269153 0.41372200192631736;
            0.4478418840074861 0.4981842896930728 0.41372200192631736 0.5200082440216406],
           [0.6114018855091738 0.34515151828171575 0.21634301778438714 0.253673516795947;
            0.34515151828171575 0.3803285686123761 0.20500253723707199 0.30420736880783067;
            0.21634301778438714 0.20500253723707199 0.4137220019263174 0.25314595728085076;
            0.253673516795947 0.30420736880783067 0.25314595728085076 0.3221534573785314],
           [0.38206058180946545 0.21634301778438714 0.17120945350064398 0.16981096699715015;
            0.21634301778438714 0.238582175642025 0.16981096699715015 0.20918696249164206;
            0.17120945350064398 0.16981096699715015 0.38206058180946545 0.21634301778438708;
            0.16981096699715015 0.20918696249164206 0.21634301778438708 0.2385821756420249],
           [0.4478418840074861 0.253673516795947 0.16981096699715015 0.19017743024985695;
            0.253673516795947 0.28111489545152707 0.1642082324755541 0.2316959448024114;
            0.16981096699715015 0.1642082324755541 0.34091363711214673 0.20500253723707199;
            0.19017743024985695 0.2316959448024114 0.20500253723707199 0.25160103976844717],
           [0.6114018855091738 0.34515151828171575 0.21634301778438714 0.253673516795947;
            0.34515151828171575 0.3803285686123761 0.20500253723707199 0.30420736880783067;
            0.21634301778438714 0.20500253723707199 0.4137220019263174 0.25314595728085076;
            0.253673516795947 0.30420736880783067 0.25314595728085076 0.3221534573785314],
           [0.626369765524018 0.3803285686123761 0.238582175642025 0.28111489545152707;
            0.3803285686123761 0.4829562753692258 0.25160103976844717 0.3918319845719349;
            0.238582175642025 0.25160103976844717 0.5200082440216407 0.32215345737853146;
            0.28111489545152707 0.3918319845719349 0.32215345737853146 0.43088268760281656],
           [0.34091363711214673 0.20500253723707199 0.16981096699715015 0.1642082324755541;
            0.20500253723707199 0.25160103976844717 0.19017743024985695 0.2316959448024114;
            0.16981096699715015 0.19017743024985695 0.4478418840074862 0.253673516795947;
            0.1642082324755541 0.2316959448024114 0.253673516795947 0.281114895451527],
           [0.4981842896930728 0.30420736880783067 0.20918696249164206 0.2316959448024114;
            0.30420736880783067 0.3918319845719349 0.2316959448024114 0.33726555786356255;
            0.20918696249164206 0.2316959448024114 0.4981842896930728 0.30420736880783067;
            0.2316959448024114 0.33726555786356255 0.30420736880783067 0.3918319845719349],
           [0.38206058180946545 0.21634301778438714 0.17120945350064398 0.16981096699715015;
            0.21634301778438714 0.238582175642025 0.16981096699715015 0.20918696249164206;
            0.17120945350064398 0.16981096699715015 0.38206058180946545 0.21634301778438708;
            0.16981096699715015 0.20918696249164206 0.21634301778438708 0.2385821756420249],
           [0.34091363711214673 0.20500253723707199 0.16981096699715015 0.1642082324755541;
            0.20500253723707199 0.25160103976844717 0.19017743024985695 0.2316959448024114;
            0.16981096699715015 0.19017743024985695 0.4478418840074862 0.253673516795947;
            0.1642082324755541 0.2316959448024114 0.253673516795947 0.281114895451527],
           [0.6748720646269153 0.4137220019263174 0.38206058180946545 0.34091363711214673;
            0.4137220019263174 0.5200082440216407 0.4478418840074862 0.4981842896930728;
            0.38206058180946545 0.4478418840074862 1.1401448653048112 0.6114018855091738;
            0.34091363711214673 0.4981842896930728 0.6114018855091738 0.626369765524018],
           [0.41372200192631736 0.25314595728085076 0.21634301778438708 0.20500253723707199;
            0.25314595728085076 0.32215345737853146 0.253673516795947 0.30420736880783067;
            0.21634301778438708 0.253673516795947 0.6114018855091738 0.34515151828171575;
            0.20500253723707199 0.30420736880783067 0.34515151828171575 0.3803285686123761],
           [0.4478418840074861 0.253673516795947 0.16981096699715015 0.19017743024985695;
            0.253673516795947 0.28111489545152707 0.1642082324755541 0.2316959448024114;
            0.16981096699715015 0.1642082324755541 0.34091363711214673 0.20500253723707199;
            0.19017743024985695 0.2316959448024114 0.20500253723707199 0.25160103976844717],
           [0.4981842896930728 0.30420736880783067 0.20918696249164206 0.2316959448024114;
            0.30420736880783067 0.3918319845719349 0.2316959448024114 0.33726555786356255;
            0.20918696249164206 0.2316959448024114 0.4981842896930728 0.30420736880783067;
            0.2316959448024114 0.33726555786356255 0.30420736880783067 0.3918319845719349],
           [0.41372200192631736 0.25314595728085076 0.21634301778438708 0.20500253723707199;
            0.25314595728085076 0.32215345737853146 0.253673516795947 0.30420736880783067;
            0.21634301778438708 0.253673516795947 0.6114018855091738 0.34515151828171575;
            0.20500253723707199 0.30420736880783067 0.34515151828171575 0.3803285686123761],
           [0.5200082440216406 0.3221534573785314 0.2385821756420249 0.25160103976844717;
            0.3221534573785314 0.43088268760281656 0.281114895451527 0.3918319845719349;
            0.2385821756420249 0.281114895451527 0.626369765524018 0.3803285686123761;
            0.25160103976844717 0.3918319845719349 0.3803285686123761 0.4829562753692258]]
    eeI = reshape(vec |> flatten, 4,4,4,4,1)

    randIdx = rand(1:4, 4)
    randVal = eeI[randIdx...,1] * ones(1,1,1,1,1)
    @test isapprox(eeInteraction(bs[randIdx]...), randVal, atol=errT)

    eeI1, uniqueIdx1 = Quiqbox.eeInteractionsCore(bs, outputUniqueIndices=true)
    isapprox(eeInteractions(bs), eeI1[:,:,:,:,1], atol=1e-15)
    @test isapprox(eeI1, eeI, atol=errT)
    uniqueInts1 = [eeI1[i...] for i in uniqueIdx1]
    @test unique(eeI1) ⊆ uniqueInts1
    @test isapprox(sort([eeI[i,j,k,l,1] for i = 1:size(eeI,1)
                                        for j = 1:i
                                        for k = 1:i
                                        for l = 1:(k==i ? j : k)]),
                   sort(uniqueInts1), atol=errT)

    eeI2, uniqueIdx2 = Quiqbox.eeInteractionsCore(bs[[1,3,2,4]], outputUniqueIndices=true)
    uniqueInts2 = [eeI2[i...] for i in uniqueIdx2]
    @test isapprox(sort(uniqueInts1), sort(uniqueInts2), atol=errT)
end


text:
function selfplay(env::GameEnv, nn::NeuralNet, num_ro::Int = 800)
  #= Plays out a self-play match, returning
  - the final position
  - the n x N² tensor of floats representing the mcts search probabilities
  - the n-ary tensor of floats representing the original value-net estimate
  where n is the number of moves in the game =#

  # Disable resign in 5% of games
  if rand() < 0.05
    resign_threshold = -1.0
  else
    resign_threshold = -0.9
  end

  player = MCTSPlayer(env, nn, num_readouts = num_ro, resign_threshold = resign_threshold)

  readouts = player.num_readouts
  initialize_game!(player)

  # Must run this once at the start, so that noise injection actually
  # affects the first move of the game.
  first_node = select_leaf(player.root)
  prob, val = nn(first_node.position)
  incorporate_results!(first_node, prob.tracker.data, val.tracker.data, first_node)

  while true
    inject_noise!(player.root)
    current_readouts = N(player.root)

    # we want to do "X additional readouts", rather than "up to X readouts".
    while N(player.root) < current_readouts + readouts
      tree_search!(player)
    end

    if should_resign(player)
      set_result!(player, -player.root.position.to_play, true)
      break
    end
    move = pick_move(player)
    play_move!(player, move)

    #println(player.root.position)
    if is_done(player.root)
      set_result!(player, result(player.root.position), false)
      break
    end
  end
  return player
end


text:
export CIFAR10
module CIFAR10
    using DataDeps
    using BinDeps
    using ColorTypes
    using FixedPointNumbers
    using ..MLDatasets: bytes_to_type, datafile, download_dep, download_docstring,
                        _colorview, _channelview

    export

        classnames,

        traintensor,
        testtensor,

        trainlabels,
        testlabels,

        traindata,
        testdata,

        convert2image,

        download

    @deprecate convert2features reshape

    const DEPNAME = "CIFAR10"
    const NCHUNKS = 5

    filename_for_chunk(file_index::Int) =
        joinpath("cifar-10-batches-bin", "data_batch_$(file_index).bin")

    const TESTSET_FILENAME =
        joinpath("cifar-10-batches-bin", "test_batch.bin")

    const CLASSES = [
        "airplane",
        "automobile",
        "bird",
        "cat",
        "deer",
        "dog",
        "frog",
        "horse",
        "ship",
        "truck",
    ]

    """
        download([dir]; [i_accept_the_terms_of_use])

    Trigger the (interactive) download of the full dataset into
    "`dir`". If no `dir` is provided the dataset will be
    downloaded into "~/.julia/datadeps/$DEPNAME".

    This function will display an interactive dialog unless
    either the keyword parameter `i_accept_the_terms_of_use` or
    the environment variable `DATADEPS_ALWAYS_ACCEPT` is set to
    `true`. Note that using the data responsibly and respecting
    copyright/terms-of-use remains your responsibility.
    """
    download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

    include(joinpath("Reader","Reader.jl"))
    include("interface.jl")
    include("utils.jl")

    function __init__()
        register(DataDep(
            DEPNAME,
            """
            Dataset: The CIFAR-10 dataset
            Authors: Alex Krizhevsky, Vinod Nair, Geoffrey Hinton
            Website: https://www.cs.toronto.edu/~kriz/cifar.html
            Reference: https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf

            [Krizhevsky, 2009]
                Alex Krizhevsky.
                "Learning Multiple Layers of Features from Tiny Images",
                Tech Report, 2009.

            The CIFAR-10 dataset is a labeled subsets of the 80
            million tiny images dataset. It consists of 60000
            32x32 colour images in 10 classes, with 6000 images
            per class.

            The compressed archive file that contains the
            complete dataset is available for download at the
            offical website linked above; specifically the binary
            version for C programs. Note that using the data
            responsibly and respecting copyright remains your
            responsibility. The authors of CIFAR-10 aren't really
            explicit about any terms of use, so please read the
            website to make sure you want to download the
            dataset.
            """,
            "https://www.cs.toronto.edu/~kriz/cifar-10-binary.tar.gz",
            "c4a38c50a1bc5f3a1c5537f2155ab9d68f9f25eb1ed8d9ddda3db29a59bca1dd",
            post_fetch_method = file -> (run(BinDeps.unpack_cmd(file,dirname(file), ".gz", ".tar")); rm(file))
        ))
    end
end


text:
#   Unit "JoB.jl" of the Diagonalization.jl Package for Julia language
#
#   MIT License
#   Copyright (c) 2019,
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements a very general form of the OJoB algoithm
#   (Congedo et al., 2012)

# """
#  Orthogonal Joint Blond Source Separation (OJoB) iterative algorithm and
#  Nonorthogonal Joint Blond Source Separation (NoJoB) iterative algorithm:
#  'Orthogonal and Non-Orthogonal Joint Blind Source Separation in the
#  Least-Squares Sense'
#  by Marco Congedo, Ronald Phlypo, Jonas Chatel-Goldman, EUSIPCO 2012
#  https://hal.archives-ouvertes.fr/hal-00737835/document
#
# if `algo` is :OJoB, runs the OJoB algorithm
# if `algo` is :NoJoB, runs the NoJoB algorithm

#  These algorithms handle several use cases corresponding to several
#  combinations of the number of datasets (m) and number of observations (k):
#  1) m=1, k>2
#  2) m>1, k=1
#  3) m>1, k>1
#
#  They take as input either the data matrices or the covariance and cross-
#  covariance matrices. Argument `input` tells the algos to take as input
#  data matrices (`input=:d`) or cov/cross-cov matrices (`input=:d`)
#
#                                   INPUT
# ---------------------------------------------------------------------------
#             Data Matrices                  Cov/cross-cov matrices
# ---------------------------------------------------------------------------
#  1) a vector of k data matrices     an array (k, 1, 1) of covariance matrices
#  2) a vector of m data matrices     an array (1, m, m) of covariance matrices
#  3) k vectors of m data matrices    an array (k, m, m) of covariance matrices
#
#  If data matrices are given as input, the covariance matrices 𝒞 are computed
#  with optional keyword arguments `covEst`, `dims`, `meanX`, `trace1` and `w`.
#
#  OJoB seeks m orthogonal matrices 𝐔=U_1,...,U_m such that the sum of the
#  squares of the diagonal part of matrices U_i'𝒞(κij)U_j
#  is maximized for all i≠j=1:m and κ=1:k.
#  The covariance is maximized rotating each data matrix as: U_i' X_κi,
#  for all i=1:m and κ=1:k.
#
#  NoJoB seeks m non-singular matrices 𝐔=U_1,...,U_m with vectors
#  u_qi, for all i=1:n scaled such that the sum u_qi'(𝒞(κj)u_qju_qj^H𝒞(κj)^H)u_qi=1
#  where the sum is over κ=1:k and j=1:m, separetedly for all q=1:m
#  The covariance is maximized transforming each data matrix as: U_i' X_κi,
#  for all i=1:m and κ=1:k.
#
#  if `fullModel=true` is passed as optional keyword argument,
#  then also the sum of the squares of the diagonal part of matrices
#  U_i'𝒞(κij)U_j is maximized for all ``i=1:m`` and ``κ=1:k``.
#  This amounts to require the intra-covariance to be diagonalized as well.
#
#  if `whitening` = true is passed as optional keyword argument,
#  for OJoB:
#       the 𝐔=U_1,...,U_m are no longer constrained to be orthogonal;
#       instead they will be invertible and a scaling constraint is imposed on them
#       so as to satisfy mean_κ=1:k U_i'𝒞(κii)U_j=I for all i=1:m.
#  for NoJoB:
#       The same constraint is imposed
#
#  For both algorithms, given m whitening matrices W_i,
#  the covariance matrices will be transformed
#  as 𝒢(κij) = W_i'𝒞(κij)W_j, for all i,j=1:m and κ=1:k.
#  Then, the OJoB or NoJob algos will be run on the set 𝒢(κij)
#  and the final solutions will be given by U_iW_i, for all i=1:m
#  Using OjOB, this amounts to requiring a generalized Canonical Correlation
#  Analysis instead of a generalized Maximum Covariance Analysis.
#  By this pre-whitening the dimension of the problem can be reduced
#  requesting matrices W_i to be rectangular. This is achieved using
#  argument `eVar`. If eVar is an integer p, the matrices in 𝒢(κij) will all
#  be pxp. If it is equal to nothing the explained variance is set to
#  the default level (0.999), otherwise if eVar is a float in [0..1], this
#  will be taken as explained variance.
#  if m=1 the eigenvalues of mean_κ=1^k 𝒞(κ11) is considered, otherwise
#         the eigenvalues of mean_κ=1^k, i=1^m 𝒞(κii) is considered.
#  These eigenvalues are normalized to unit sum and the cumulated sum is
#  considered. Finally p is found on this accumulated regularized ev vector
#  using the `eVarMethod` function.
#
#  if sort=true (default) the column vectors of the 𝐔 matrices are reordered
#  so as to allow positive and sorted in descending order
#  diagonal elements of the products ``U_i'𝒞(κij)U_j`` as much as possible.
#  If the NoJoB algorithm is used the column vectors of the matrices U_i
#  are normalized to unit norm.
#  Note that if `whitening` is true the output matrices U_i will not have
#  unit norm columns as they are multiplied by the whiteners after being
#  scaled and sorted and before being returned.

#  By passing a matrix if k=1 or a vector of k matrices if k>1 as `init`,
#  you can smartly initialize `𝐔`, for example, when excluding some subjects
#  and running the OJoB again.
#
#  `tol` is the convergence to be attained.
#
#  `maxiter` is the maximum number of iterations allowed.
#
#  if `verbose`=true, the convergence attained at each iteration and other
#  information will be printed.
#
#  if `threaded`=true (default) and n>x and x>1, where x is the number
#  of threads Julia is instructed to use (the output of Threads.nthreads())
#  the algorithms run in multithreaded mode paralellising several comptations
#  over n.
#
#  Besides optionally multi-threaded, these algorithms heavely use BLAS.
#  Before running this function you may want to set the number of threades
#  Julia is instructed to use to the number of logical CPUs of your machine
#  and set `BLAS.set_num_threads(Sys.CPU_THREADS)`. See:
#  https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Threads-1
# """
function JoB(𝐗::AbstractArray, m::Int, k::Int, input::Symbol, algo::Symbol, type;
              covEst   :: StatsBase.CovarianceEstimator=SCM,
              dims     :: Int64 = 1,
              meanX    :: Tmean = 0,
              trace1   :: Bool = false,
              w        :: Twf = ○,
          fullModel :: Bool = false,
          preWhite  :: Bool = false,
          sort      :: Bool = true,
          init      :: Union{Matrix, Nothing} = ○,
          tol       :: Real = 0.,
          maxiter   :: Int  = 1000,
          verbose   :: Bool = false,
          threaded  :: Bool = true,
      eVar     :: TeVaro = ○,
      eVarMeth :: Function = searchsortedfirst)

    k<3 && m<2 && throw(ArgumentError("Either `k` must be equal to at least 3 or `m` must be equal to at least 2"))
    #if m==1 fullModel=false end

    # if input ≠:d the input is supposed to be the cross-covariance matrices
    input==:d ? 𝒞=_crossCov(𝐗, m, k;
                            covEst=covEst, dims=dims, meanX=meanX) : 𝒞=𝐗

    𝒢=deepcopy(𝒞)
    if trace1 || w ≠ ○ _normalize!(𝒢, m, k, trace1, w) end

    if eVar isa Int && size(𝒢[1, 1, 1], 1) ≤ eVar ≤ 0
        eVar=size(𝒢[1, 1, 1], 1)
    end

    iter, conv, oldconv, 😋, conv_ = 1, 0., 1.01, false, 0.
    tol==0. ? tolerance = √eps(real(type)) : tolerance = tol

    # pre-whitening
    if preWhite
        # for the moment being computed on the average across k and m
        # if eVar is not an Int.
        # This way for all i=1:m the dimensionality reduction is fixed
        # for the case m=1 this is the classical whitening
        if !(eVar isa Int)
            G=ℍ(𝛍(𝒢[κ, i, i] for κ=1:k, i=1:m))
            p, arev = _getssd!(eVar, eigvals(G), size(𝒢[1, 1, 1], 1), eVarMeth)
            𝑾=[whitening(ℍ(𝛍(𝒢[κ, i, i] for κ=1:k)); eVar=p) for i=1:m]
        else
            𝑾=[whitening(ℍ(𝛍(𝒢[κ, i, i] for κ=1:k)); eVar=eVar) for i=1:m]
        end

        if m==1
            @inbounds for κ=1:k
                𝒢[κ, 1, 1] = 𝑾[1].F' * 𝒢[κ, 1, 1] * 𝑾[1].F
            end
        else
            @inbounds for κ=1:k, i=1:m-1, j=i+1:m
                𝒢[κ, i, j] = 𝑾[i].F' * 𝒢[κ, i, j] * 𝑾[j].F
                𝒢[κ, j, i] = 𝒢[κ, i, j]'
            end
            @inbounds for κ=1:k, i=1:m
                𝒢[κ, i, i] = 𝑾[i].F' * 𝒢[κ, i, i] * 𝑾[i].F
            end
        end
    end
    n=size(𝒢[1, 1, 1], 1)

    ⏩ = n>=Threads.nthreads() && Threads.nthreads()>1 && threaded

    # initialization of 𝐔[i], i=1:m, as the eigenvectors of sum_k,j(𝒢_k,i,j*𝒢_k,i,j')
    # see Eq. 17 of Congedo, Phlypo and Pham, 2011, or with the provided matrices
    # note: gemm supports complex matrices
    ggt(κ::Int, i::Int, j::Int) = BLAS.gemm('N', 'T', 𝒢[κ, i, j], 𝒢[κ, i, j])
    if m>1 # gmca, gcca, majd
        if init === nothing
            if algo==:NoJoB
                𝐔 = [Matrix{type}(I, n, n) for i=1:m]
            else
                if fullModel
                    𝐔 = [eigvecs(𝛍(ggt(κ, i, j) for j=1:m, κ=1:k)) for i=1:m]
                else
                    𝐔 = [eigvecs(𝛍(ggt(κ, i, j) for j=1:m, κ=1:k if j≠i)) for i=1:m]
                end
            end
        else
            𝐔 = init
        end
    else  # ajd
        if init === nothing
            if algo==:NoJoB
                𝐔 = [Matrix{type}(I, n, n)]
            else
                𝐔 = [eigvecs(𝛍(ggt(κ, 1, 1) for κ=1:k))]
            end
        else
            𝐔 = [init]
        end
    end

    function updateR!(η, i, j)  # 𝐑[η] += (𝒢[κ, i, j] * 𝐔[j][:, η]) times its transpose
        # don't use BLAS for complex data
        if ⏩ # if threaded don't share memory for Ω but use 𝛀
            if type<:Real
                @inbounds for κ=1:k
                    𝛀[η][:, κ] = BLAS.gemv('N', 𝒢[κ, i, j], 𝐔[j][:, η])
                end
                𝐑[η] += Hermitian(BLAS.gemm('N', 'T', 𝛀[η], 𝛀[η])) # (Ω * Ω')
            else
                @inbounds for κ=1:k
                    𝛀[η][:, κ] = 𝒢[κ, i, j] * 𝐔[j][:, η]
                end
                𝐑[η] += Hermitian(𝛀[η] * 𝛀[η]') # (Ω * Ω')
            end
        else
            if type<:Real
                @inbounds for κ=1:k
                    Ω[:, κ] = BLAS.gemv('N', 𝒢[κ, i, j], 𝐔[j][:, η])
                end
                𝐑[η] += Hermitian(BLAS.gemm('N', 'T', Ω, Ω)) # (Ω * Ω')
            else
                @inbounds for κ=1:k
                    Ω[:, κ] = 𝒢[κ, i, j] * 𝐔[j][:, η]
                end
                𝐑[η] += Hermitian(Ω * Ω') # (Ω * Ω')
            end
        end
    end

    function update!(m, n, i)

        function udR1!(η) # case m=1
            fill!(𝐑[η], type(0))
            updateR!(η, 1, 1)
        end

        function udRm!(η, i)  # case m>1
            fill!(𝐑[η], type(0))
            for j=1:m i≠j ? updateR!(η, i, j) : nothing end # j ≠ i
            fullModel ? updateR!(η, i, i) : nothing         # j = i
        end

        if m==1
            ⏩ ? (@threads for η=1:n udR1!(η) end) : (for η=1:n udR1!(η) end)
        else
            ⏩ ? (@threads for η=1:n udRm!(η, i) end) : (for η=1:n udRm!(η, i) end)
        end
    end

    # pre-allocate memory
    𝐑 = HermitianVector([Hermitian(zeros(type, n, n)) for η=1:n])
    ⏩ ? (𝛀 = [Matrix{type}(undef, n, k) for η=1:n]) :
         (Ω = Matrix{type}(undef, n, k))

    # # # # # here starts the algorithm

    if algo==:OJoB
        verbose && @info("Iterating OJoB algorithm...")
        while true
            conv_ =0.
            @inbounds for i=1:m # m optimizations for updating 𝐔[1]...𝐔[m]

                update!(m, n, i)

                # do 1 power iteration, not worth threading here
                for η=1:n 𝐔[i][:, η] = 𝐑[η] * 𝐔[i][:, η] end

                conv_ += PosDefManifold.ss(𝐔[i]) # square of the norms of power iteration vectors

                # Lodwin Orthogonalization and update 𝐔[i]<-UV', with svd(𝐔[i])=UwV'
                𝐔[i] = PosDefManifold.nearestOrth(𝐔[i])
            end

            conv_ =sqrt(conv_ /(n^2*m))
            iter==1 ? conv=1. : conv = abs((conv_-oldconv)/oldconv)  # relative change

            verbose && println("iteration: ", iter, "; convergence: ", conv)
            (diverging = conv < 0) && verbose && @warn("OJoB diverged at:", iter)
            (overRun = iter == maxiter) && @warn("OJoB: reached the max number of iterations before convergence:", iter)
            (😋 = 0. <= conv <= tolerance) || overRun==true ? break : nothing
            oldconv=conv_
            iter += 1
        end # while
    else
        verbose && @info("Iterating NoJoB algorithm...")

        # solve Lx=𝐑[η]*𝐔[i][:, η] for x and L'y=x for y
        # and scale the ηth column as 𝐔[i][:, η] <- y/sqrt(y'𝐑[η]t)
        function triS!(cho, η,  i)
            y=cho.U\(cho.L\(𝐑[η]*𝐔[i][:, η]))
            𝐔[i][:, η]=y*inv(sqrt(PosDefManifold.qf(y, 𝐑[η])))
        end

        # thread sum of 𝐑[η] and the solutions to triangular systems
        # only if n is at least Threads.nthreads()*2 (no worth otherwise)
        ⏩x = ⏩ && n≥Threads.nthreads()*2

        while true
            conv_ =0.
            @inbounds for e2=1:2, i=1:m # m optimizations for updating 𝐔[1]...𝐔[m]
                                        # double loop to avoid oscillating convergence
                update!(m, n, i)

                # do 1 power iteration
                cho=cholesky(⏩x ? fVec(sum, 𝐑) : sum(𝐑)) # Cholesky LL'of 𝐑[1]+...+𝐑[n]

                ⏩x ? (@threads for η=1:n triS!(cho, η,  i) end) :
                            (for η=1:n triS!(cho, η,  i) end)

                conv_+=PosDefManifold.ss(𝐔[i])
            end

            conv_=conv_/(2*n^2*m) # 2 because of the double loop
            iter==1 ? conv=1. : conv = abs((conv_-oldconv)/oldconv)  # relative change

            verbose && println("iteration: ", iter, "; convergence: ", conv)
            (diverging = conv < 0) && verbose && @warn("NoJoB diverged at:", iter)
            (overRun = iter == maxiter) && @warn("NoJoB: reached the max number of iterations before convergence:", iter)
            (😋 = 0. <= conv <= tolerance) || overRun==true ? break : nothing
            oldconv=conv_
            iter += 1
        end # while
    end

    verbose && @info("Convergence has "*(😋 ? "" : "not ")*"been attained.\n\n")

    # scale and permute the vectors of U_1,...,U_m
    if sort
        algo==:NoJoB ? for i=1:m normalizeCol!(𝐔[i], 1:size(𝐔[i], 2)) end : nothing
        λ = m==1 ? _permute!(𝐔[1], 𝒢, k, :c) :
                   _flipAndPermute!(𝐔, 𝒢, m, k, :c)
    else
        λ = m==1 ? 𝛍([𝐔[1][:, η]'*𝒢[l, 1, 1]*𝐔[1][:, η] for η=1:n] for l=1:k) :
                   𝛍([𝐔[i][:, η]'*𝒢[l, i, j]*𝐔[j][:, η] for η=1:n] for l=1:k, j=1:m, i=1:m if i≠j)
    end

    if preWhite
        algo==:OJoB ? 𝐕=[𝐔[i]'*𝑾[i].iF for i=1:m] :
                      𝐕=[pinv(𝐔[i])*𝑾[i].iF for i=1:m] # algo==:NoJoB
        for i=1:m 𝐔[i] = 𝑾[i].F * 𝐔[i] end
    else
        algo==:OJoB ? 𝐕=[Matrix(𝐔[i]') for i=1:m] :
                      𝐕=[pinv(𝐔[i]) for i=1:m] # algo==:NoJoB
    end

    return m>1 ? (𝐔, 𝐕, λ, iter, conv) : (𝐔[1], 𝐕[1], λ, iter, conv)
end


text:
using LinearAlgebra
using SparseArrays
using Plots

"This routine solves the heat equation using Backwards Euler in time and finite differences
in space. The output is saved as a .gif file (which should play automatically when this file is run)."

m = 100 # number of points
dt = .01 # time-step
T = 2.0

u0(x) = 0.0 # initial condition
f(x) = 5*(Float64((x > -.5) && (x <= 0.0)) - Float64((x < .5) && (x >= 0.0)))
f(x,t) = f(x)*exp(-t)
α(t) = 1.0
β(t) = sin(2*pi*t) # du/dx

# define spatial grid
x = LinRange(-1,1,m+2)
xint = x[2:end-1]
h = x[2]-x[1]
A = (1/h^2) * spdiagm(0=>2*ones(m),-1=>-ones(m-1),1=>-ones(m-1))

function F(xint,t)
    b = f.(xint,t)
    b[1] += α(t)/h^2
    b[end] += β(t)/h^2
    return b
end

u = u0.(xint)
Nsteps = ceil(Int,T/dt)
dt = T / Nsteps

interval = 2
@gif for i = 1:Nsteps
    t = (i+1)*dt
    u .= (I + dt*A)\(u + dt*F(xint,t))
    if i % interval == 0 # % = mod function
        plot(xint,u,linewidth=2,label="Solution",ylims=(-1.0,3.0))
        println("on timestep $i out of $Nsteps.")
    end
end every interval


text:
ENV["GKSwstype"] = "nul"  # For StatsPlots
println("pid: ", getpid())

# Read command line args.
if length(ARGS) > 1
  resultsdir = ARGS[1]
  awsbucket = ARGS[2]
  snum = parse(Int, ARGS[3])
  println("ARGS: ", ARGS)
else
  resultsdir = "results/test/"
  awsbucket = nothing
  snum = 4
end
flush(stdout)

# Create results/images dir if needed
imgdir = "$(resultsdir)/img"
mkpath(imgdir)

import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))

using CytofDensityEstimation
const CDE = CytofDensityEstimation
using LaTeXStrings
using Distributions
using StatsPlots
using HypothesisTests
import Random
using RCall
using BSON
import CytofDensityEstimation.Model.default_ygrid
using StatsFuns

include("scenarios.jl")

plotsize = (400, 400)

# Simulate data.
Random.seed!(5);
simdata = scenarios(snum)
yC, yT = simdata[:yC], simdata[:yT]

# Plot true data density.
function plot_trude_data_density()
  plot(default_ygrid(), pdf.(simdata[:mmC], default_ygrid()), lw=3,
       label=L"y_C", color=:blue)
  plot!(default_ygrid(), pdf.(simdata[:mmT], default_ygrid()), lw=3,
        label=L"y_T", color=:red)
end
plot_trude_data_density()
plot!(size=plotsize)
savefig(joinpath(imgdir, "data-true-density.pdf"))
closeall()

# Plot observed data density.
function plot_observed_data_density()
  density(yC[isfinite.(yC)],  lw=3, label=L"y_C", color=:blue)
  density!(yT[isfinite.(yT)], lw=3, label=L"y_T", color=:red)
end
plot_observed_data_density()
plot!(size=plotsize)
savefig(joinpath(imgdir, "data-kde.pdf"))
closeall()

# Define data, prior, and initial state.
Random.seed!(7)
K = 3
data = CDE.Model.Data(yC, yT)
prior_mu = let
  yfinite = [data.yC_finite; data.yT_finite]
  Normal(mean(yfinite), std(yfinite))
end
prior = CDE.Model.Prior(K, mu=prior_mu, nu=LogNormal(1.6, 0.4), p=Beta(100, 100),
                        omega=InverseGamma(.1, .1), psi=Normal(-1, 1))
state = CDE.Model.State(data, prior)
tuners = CDE.Model.Tuners(K, 0.1)

println("Priors:")
for fn in fieldnames(CDE.Model.Prior)
  println(fn, " => ", getfield(prior, fn))
end

# Specify parameters to keep constant.
fix = Symbol[:beta, :p]

# Run chain.
state.p = 0.5

state.beta = false
@time chain0, laststate0, summarystats0 = CDE.Model.fit(
    state, data, prior, deepcopy(tuners), fix=fix,
    nsamps=[2000], nburn=2000, thin=1, rep_aux=10)

state.beta = true
@time chain1, laststate1, summarystats1 = CDE.Model.fit(
    state, data, prior, deepcopy(tuners), fix=fix,
    nsamps=[2000], nburn=2000, thin=1, rep_aux=10)

# Posterior Odds (in favor of treatment effect).
@time post_odds = CDE.Model.posterior_odds(data, chain0[1], chain1[1])
println("Posterior Odds: ", post_odds); flush(stdout)

imgdir0 = "$(resultsdir)/bf0/img"
mkpath(imgdir0)
@time CDE.Model.plotpostsummary(chain0, summarystats0, data.yC, data.yT, imgdir0)
flush(stdout)

imgdir1 = "$(resultsdir)/bf1/img"
mkpath(imgdir1)
@time CDE.Model.plotpostsummary(chain1, summarystats1, data.yC, data.yT, imgdir1)
flush(stdout)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# Copyright (c) 2021 Oscar Dowson and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Performance tips

# This tutorial explains how to write performant JuMP code.

# !!! tip
#     Read the [Performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html)
#     section of the Julia manual. The most important rule is to avoid global
#     variables!

using JuMP  # hide

# ## Use macros to build expressions

# ### What

# Use JuMP's macros (or [`add_to_expression!`](@ref) to build expressions.
# Avoid constructing expressions outside the macros.

# ### Why

# Constructing an expression outside the macro results in intermediate copies of
# the expression. For example,
# ```julia
# x[1] + x[2] + x[3]
# ```
# is equivalent to
# ```julia
# a = x[1]
# b = a + x[2]
# c = b + x[3]
# ```
# Since we only care about `c`, the `a` and `b` expressions are not needed and
# constructing them slows the program down!

# JuMP's macros rewrite the expressions to operate in-place and avoid these
# extra copies. Because they allocate less memory, they are faster, particularly
# for large expressions.

# ### Example

model = Model()
@variable(model, x[1:3])

# Here's what happens if we construct the expression outside the macro:

@allocated x[1] + x[2] + x[3]

# !!! info
#     The `@allocated` measures how many bytes were allocated during the
#     evaluation of an expression. Fewer is better.

# If we use the [`@expression`](@ref) macro, we get many fewer allocations:

@allocated @expression(model, x[1] + x[2] + x[3])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Example models to simulate from...
abstract type EBayesSimulation end

struct NormalEBayesSimulationResult{Sim<:EBayesSimulation, NS<:NormalSamples}
    ss::NS
    true_μs::Vector{Float64}
    sim::Sim
end

function NormalSamples(nsim::NormalEBayesSimulationResult)
   nsim.ss
end

#predict(method, sim_res::NormalEBayesSimulationResult) = predict(method, sim_res.ss)

struct NormalNormalSimulation <: EBayesSimulation
   n::Int64
   prior_μ::Float64
   A::Float64
   σ::Float64
end

function NormalNormalSimulation(;n=1000,
                                 prior_μ=0.0,
                                 A=1.0,
                                 σ=1.0)
   NormalNormalSimulation(n, prior_μ, A, σ)
end

function rand(x::NormalNormalSimulation)
   n = x.n
   A = x.A
   σ = x.σ
   prior_μ = x.prior_μ

   μs = randn(n).*sqrt(A) .+ prior_μ
   Zs = randn(n).*σ .+ μs
   σs = fill(σ, n)

   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end
#----------------------------------------------------------------------
# Xie-Kou-Brown 2012 include 6 examples for their simulation benchmark
#----------------------------------------------------------------------

struct XieKouBrownExample1 <: EBayesSimulation
    n::Int64
end

function rand(x::XieKouBrownExample1)
   n = x.n
   As = rand(Uniform(0.1,1), n)
   σs = sqrt.(As)
   μs = randn(n)
   Zs = randn(n).*σs .+ μs
   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end

struct XieKouBrownExample2 <: EBayesSimulation
    n::Int64
end

function rand(x::XieKouBrownExample2)
   n = x.n
   As = rand(Uniform(0.1,1), n)
   σs = sqrt.(As)
   μs = rand(n)
   Zs = randn(n).*σs .+ μs
   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end

struct XieKouBrownExample3 <: EBayesSimulation
    n::Int64
end

function rand(x::XieKouBrownExample3)
   n = x.n
   As = rand(Uniform(0.1,1), n)
   σs = sqrt.(As)
   μs = As
   Zs = randn(n).*σs .+ μs
   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end

struct XieKouBrownExample4 <: EBayesSimulation
    n::Int64
end

function rand(x::XieKouBrownExample4)
   n = x.n
   As = rand(InverseGamma(5, 1/2), n)
   σs = sqrt.(As)
   μs = As
   Zs = randn(n).*σs .+ μs
   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end

struct XieKouBrownExample5 <: EBayesSimulation
    n::Int64
end

function rand(x::XieKouBrownExample5)
   n = x.n
   groups = rand(Bernoulli(1/2), n)
   As = [ (g==0) ? 0.1 : 0.5 for g in groups]
   σs = sqrt.(As)
   μs = [ (g==0) ? rand(Normal(2.0,sqrt(0.1))) : rand(Normal(0.0, sqrt(0.5))) for g in groups]
   Zs = randn(n).*σs .+ μs
   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end

struct XieKouBrownExample6 <: EBayesSimulation
    n::Int64
end

function rand(x::XieKouBrownExample6)
   n = x.n
   As = rand(Uniform(0.1,1), n)
   σs = sqrt.(As)
   μs = As
   Zs =  μs .+ rand.(Uniform.( -sqrt.(3 .* As), sqrt.(3 .* As)))
   NormalEBayesSimulationResult(NormalSamples(Zs,σs), μs, x)
end

#-----------------------------------------
# Simulations with covariates
#-----------------------------------------

struct CovariateNormalEBayesSimulationResult{Sim<:EBayesSimulation,
                                    Tbl,
                                    NS<:NormalSamples}
    ss::NS
    X::Tbl
    true_ms::Vector{Float64}
    true_μs::Vector{Float64}
    sim::Sim
end

function NormalSamples(nsim::CovariateNormalEBayesSimulationResult)
   nsim.ss
end

struct FriedmanFayHerriotSimulation{T<:Number} <: EBayesSimulation
    n::Int
    p::Int
    A_sqrt::T
    σ::T
end

friedman_reg(x) = 10*sin(π*x[1]*x[2]) + 20*(x[3] - 1/2)^2 + 10*x[4] +5*x[5]

function rand(sim::FriedmanFayHerriotSimulation)
   p = sim.p
   n = sim.n
   σ = sim.σ
   A_sqrt = sim.A_sqrt
   X = rand(Uniform(0,1), n, p)
   ms = map(friedman_reg, eachrow(X))
   μs = ms .+ rand(Normal(0, A_sqrt), n)
   Zs = μs .+ rand(Normal(0, σ), n)
   X = MLJBase.table(X)
   CovariateNormalEBayesSimulationResult(NormalSamples(Zs,fill(σ, n)), X, ms, μs, sim)
end


text:
# test ScatterAlignPose2

using Test
using Images
using Caesar
using Distributions
using Manifolds

# test plotting helper functions
using Gadfly
using Random

import Rotations as _Rot

##
@testset "Test ScatterAlignPose2" begin
##

x = -10:0.1:10;
y = -10:0.1:10;

x = -10:0.1:10;

σ = 0.1

g = (x,y)->pdf(MvNormal([3.;0],[σ;σ]),[x;y]) + pdf(MvNormal([8.;0.0],[σ;σ]),[x;y]) + pdf(MvNormal([0;5.0],[σ;σ]),[x;y])

bIM1 = zeros(length(x),length(y))
bIM2 = zeros(length(x),length(y))

oT = [2.; 0]
oΨ =  pi/6

M = SpecialEuclidean(2)
e0 = identity_element(M)
pCq = [oT;oΨ]
qGp = inv(M, exp(M, e0, hat(M, e0, pCq)))
qTp = affine_matrix(M, qGp )

qCp = vee(M, e0, log(M, e0, qGp))

##

for (i,x_) in enumerate(x), (j,y_) in enumerate(y)
  bIM1[i,j] = g(x_,y_)
  v = qTp*[x_;y_;1.0]
  _x_, _y_ = v[1], v[2]
  bIM2[i,j] = g(_x_, _y_)
end

sap = ScatterAlignPose2(bIM1, bIM2, (x,y); sample_count=100, bw=1.0, cvt=(im)->im)

# requires IIF at least v0.25.6
@test sample(sap.cloud1,1) isa Tuple
@test sample(sap.cloud2,10)[1] isa AbstractArray

## test plotting function

snt = overlayScatterMutate(sap; sample_count=50, bw=1., user_coords=[0.;0;0*pi/6]); # , user_offset=[0.;0;0.]);
plotScatterAlign(snt; title="\npCq=$(round.(pCq,digits=2))")

##

# inverse for q --> p
@test isapprox( pCq[1:2], snt.best_coords[1:2]; atol=0.3 )
@test isapprox( pCq[3], snt.best_coords[3]; atol=0.2 )

## check packing and unpacking

psap = convert(PackedScatterAlignPose2, sap);
sap_ = convert(ScatterAlignPose2, psap);

@test sap.gridscale == sap_.gridscale
@test sap.sample_count == sap_.sample_count
@test sap.bw == sap_.bw

@test isapprox(sap.cloud1, sap_.cloud1, mmd_tol=1e-2)
@test isapprox(sap.cloud2, sap_.cloud2, mmd_tol=1e-2)

## check that optimize works (using the same tfg)

tfg = initfg()
getSolverParams(tfg).attemptGradients = false
M = getManifold(sap)
e0 = identity_element(M)
meas = sample(sap.cloud1,100)[1], [ProductRepr(sample(sap.cloud2,1)[1][1],[1 0; 0 1.]) for _ in 1:100], M
δ(x) = calcFactorResidualTemporary(sap, (Pose2,Pose2), meas, (e0,ProductRepr(x[1:2],_Rot.RotMatrix(x[3]))) , tfg=tfg)[1]

@show δ([0;0;0.]);
@show δ([1.;0;0.]);

@test isapprox(δ([0;0;0.]), δ([0;0;0.]); atol=1e-6)
@test isapprox(δ([10;0;0.]), δ([10;0;0.]); atol=1e-6)
@test !isapprox( δ([0;0;0.]), δ([0.1;0;0.]), atol=1e-6 )
@test !isapprox( δ([0;0;0.]), δ([0;0;0.1]), atol=1e-6 )

## build into graph

fg = initfg()
getSolverParams(fg).inflateCycles=1

addVariable!(fg, :x0, Pose2)
addVariable!(fg, :x1, Pose2)

addFactor!(fg, [:x0;], PriorPose2(MvNormal([0;0;0.],[0.01;0.01;0.01])))
addFactor!(fg, [:x0;:x1], sap, inflation=0.0)

## use in graph

X1 = approxConvBelief(fg, :x0x1f1, :x1)

c1 = AMP.makeCoordsFromPoint(getManifold(Pose2), mean(X1))
@show c1

##

@test isapprox( pCq[1:2], c1[1:2], atol=0.5 )
@test isapprox( pCq[3], c1[3],   atol=0.3 )

## check save and load of sap

saveDFG("/tmp/caesar/test_sap", fg)
fg_ = loadDFG("/tmp/caesar/test_sap")

##

sap = getFactorType(fg, :x0x1f1)
sap_ = getFactorType(fg_, :x0x1f1)

@test isapprox( sap.cloud1, sap_.cloud1)
@test isapprox( sap.cloud2, sap_.cloud2)

@test sap.gridscale == sap_.gridscale
@test sap.sample_count == sap_.sample_count
@test sap.bw == sap_.bw

##
end

@testset "test ScatterAlignPose2 with MKD direct" begin
##

# setup

oT = [2.; 0]
oΨ =  pi/6

M = SpecialEuclidean(2)
e0 = identity_element(M)
pCq = [oT;oΨ]
qTp = affine_matrix(M, exp(M, e0, hat(M, e0, pCq)))

##

# Points in XY only

# g = (x,y)->pdf(MvNormal([3.;0],[σ;σ]),[x;y]) + pdf(MvNormal([8.;0.0],[σ;σ]),[x;y]) + pdf(MvNormal([0;5.0],[σ;σ]),[x;y])
p1 = vcat([0.1*randn(2).+[3;0.] for i in 1:50], [0.1*randn(2)+[8.;0] for i in 1:50], [0.1*randn(2)+[0;5.] for i in 1:50])
# foreach(pt->(pt[1] += 100), p1)
shuffle!(p1)
P1 = manikde!(getManifold(Point2), p1)

p2 = vcat([0.1*randn(2).+[3;0.] for i in 1:50], [0.1*randn(2)+[8.;0] for i in 1:50], [0.1*randn(2)+[0;5.] for i in 1:50])
# foreach(pt->(pt[1] += 100), p2)
# adjust points
for (i,pt) in enumerate(p2)
  v = qTp*[pt;1.0]
  pt[1:2] .= v[1:2]
end
shuffle!(p2)
P2 = manikde!(getManifold(Point2), p2)

sap = ScatterAlignPose2(;cloud1=P1, cloud2=P2, sample_count=100, bw=1.0)

## test plotting function

snt = overlayScatterMutate(sap; sample_count=100, bw=2.0, user_coords=[0.;0;0*pi/6]);
plotScatterAlign(snt; title="\npCq=$(round.(pCq,digits=2))")

##

# inverse for q --> p
@test isapprox( oT, snt.best_coords[1:2]; atol=1.0 )
@test isapprox( oΨ,   snt.best_coords[3]; atol=0.5 )

##

fg = initfg()
getSolverParams(fg).inflateCycles=1

addVariable!(fg, :x0, Pose2)
addVariable!(fg, :x1, Pose2)

addFactor!(fg, [:x0], PriorPose2(MvNormal([0.01;0.01;0.01])))
addFactor!(fg, [:x0;:x1], sap, inflation=0.0)

## check residual calculation

# see #1415
meas = sample(P1,100)[1], [ProductRepr([0;0.],[1 0; 0 1.]) for _ in 1:100], M
δ1 = calcFactorResidualTemporary(sap, (Pose2,Pose2), meas, (e0,e0))

meas = sample(P1,100)[1], [ProductRepr(sample(P2,1)[1][1],[1 0; 0 1.]) for _ in 1:100], M
δ2 = calcFactorResidualTemporary(sap, (Pose2,Pose2), meas, (e0,e0))

# check different cloud samplings produce different residual values
@test !isapprox(δ1, δ2,  atol=1e-4)

## check that optimize works (using the same tfg)

tfg = initfg()
# meas = sample(P1,100)[1], [ProductRepr(sample(P2,1)[1][1],[1 0; 0 1.]) for _ in 1:100], M
meas = sampleFactor(fg, :x0x1f1)[1]
δ(x) = calcFactorResidualTemporary(sap, (Pose2,Pose2), meas, (e0,ProductRepr(x[1:2],_Rot.RotMatrix(x[3]))), tfg=tfg)[1]

@show δ([0;0;0.])
@show δ([1.;0;0.])

@test !isapprox( δ([0;0;0.]), δ([1.;0;0.]), atol=1e-6 )

##

X1 = approxConvBelief(fg, :x0x1f1, :x1)
c1 = AMP.makeCoordsFromPoint(getManifold(Pose2), mean(X1))
@show c1

##

@test isapprox( pCq[1:2], c1[1:2], atol=0.5 )
@test isapprox( pCq[3], c1[3],   atol=0.3 )

##
end

#


text:
abstract type SourceSymbol end

struct DCSymbol <: SourceSymbol
    bits::UInt8
end

struct ACSymbol <: SourceSymbol
    bits::UInt8
    ACSymbol(bits::UInt8) = new(bits)
    ACSymbol(run::Integer, amp::Integer) = new(run << 0x4 | (amp & 0x0F))
end
const ZRL = ACSymbol(0xF0)
const EOB = ACSymbol(0x00)
const EOB0 =  ACSymbol(0x00)
const EOB1 =  ACSymbol(0x10)
const EOB2 =  ACSymbol(0x20)
const EOB3 =  ACSymbol(0x30)
const EOB4 =  ACSymbol(0x40)
const EOB5 =  ACSymbol(0x50)
const EOB6 =  ACSymbol(0x60)
const EOB7 =  ACSymbol(0x70)
const EOB8 =  ACSymbol(0x80)
const EOB9 =  ACSymbol(0x90)
const EOB10 = ACSymbol(0xA0)
const EOB11 = ACSymbol(0xB0)
const EOB12 = ACSymbol(0xC0)
const EOB13 = ACSymbol(0xD0)
const EOB14 = ACSymbol(0xE0)

amplitude(s::DCSymbol) = s.bits
amplitude(s::ACSymbol) = s.bits & 0x0F
runlength(s::ACSymbol) = s.bits >> 0x4
iszrl(s::ACSymbol) = s.bits === 0xF0
iseob(s::ACSymbol) = s.bits === 0x00
iseobn(s::ACSymbol) = (s.bits & 0x0F === 0x00) & (s.bits < 0xF0)

function Base.show(io::IO, s::T) where T <: Union{DCSymbol, ACSymbol}
    if get(io, :typeinfo, Any) === T
        print(io, "0x", hex(s.bits))
    else
        print(io, T, "(0x", hex(s.bits), ")")
    end
end

Base.convert(T::Type{<:SourceSymbol}, bits::UInt8) = T(bits)

"""
    DCTCoeff{T}

A type for the intermediate representation of quantized DCT coefficients.
See also [`DCCoeff`](@ref) and [`ACCoeff`](@ref).
"""
struct DCTCoeff{T<:SourceSymbol}
    symbol::T
    data::CodeBits
end

const DCCoeff = DCTCoeff{DCSymbol}
const ACCoeff = DCTCoeff{ACSymbol}

"""
    DCCoeff(v)

Create a DC coefficient with value of `v`.
"""
function DCCoeff(v::Integer)
    c = encode_coeff(v)
    DCCoeff(DCSymbol(length(c)), c)
end

"""
    ACCoeff(v, run::Integer=0)
    ACCoeff(v, s::ACSymbol)

Create an AC coefficient with value of `v`.
"""
function ACCoeff(v::Integer, run::Integer=0)
    c = encode_coeff(v)
    ACCoeff(ACSymbol(run, length(c)), c)
end
function ACCoeff(v::Integer, s::ACSymbol)
    if iseobn(s)
        c = encode_eoblength(v)
        ACCoeff(ACSymbol(length(c), 0), c)
    else
        c = encode_coeff(v)
        ACCoeff(ACSymbol(runlength(s), length(c)), c)
    end
end

# 0 => length=0
# -1, 1 => length=1
# -3:-2, 2:3 => length=2
# -7:-4, 4:7 => length=3
const N_DIGITS = SVector{256, Int}([ceil(log2(i)) for i = 1:256])

function encode_coeff(v::Integer)
    n = abs(v)
    @inbounds len = n < 256 ? N_DIGITS[n + 1] : N_DIGITS[n >> 0x8 + 1] + 8
    bits = unsafe_trunc(UInt16, v < 0 ? v - 1 : v)
    CodeBits(len, bits)
end
function decode_coeff(c::CodeBits)
    len = length(c)
    len == 0 && return zero(Int16)
    mask = (0xFFFF >> UInt8(len) << UInt8(len)) + 0x001
    signed(word(c) + ifelse(c[1], 0x0000, mask))
end
function encode_eoblength(v::Integer)
    n = max(v, oneunit(v))
    @inbounds len = n < 256 ? N_DIGITS[n + 1] : N_DIGITS[n >> 0x8 + 1] + 8
    CodeBits(len - 1,  unsafe_trunc(UInt16, v))
end
function decode_eoblength(c::CodeBits)
    signed(0x0001 << UInt8(length(c)) + word(c))
end

"""
    length(coeff::ACCoeff)

Return the length of coefficients with the run-length of zero coefficients.

!!! note
    Since EOB and EOBn have no finite length of coefficients, this method
    returns `0` and does not throw errors.
    If you want to get the length of bands which have an end-of-band, convert
    an ACCoeff with EOBn to `Signed`, e.g. ```Int(coeff)```.
"""
function Base.length(coeff::ACCoeff)
    iseobn(coeff.symbol) && return 0
    iszrl(coeff.symbol) && return 15
    runlength(coeff.symbol) + 1
end

function Base.getindex(coeff::ACCoeff, i::Integer)
    length(coeff) == i && return Int(coeff)
    zero(i)
end
function Base.iterate(coeff::ACCoeff, state::Int=0)
    state >= length(coeff) && return nothing
    state += 1
    (coeff[state], state)
end

function Base.convert(::Type{T}, coeff::DCCoeff) where T <: Signed
    convert(T, decode_coeff(coeff.data))
end
function Base.convert(::Type{T}, coeff::ACCoeff) where T <: Signed
    if iseobn(coeff.symbol)
        convert(T, decode_eoblength(coeff.data))
    else
        convert(T, decode_coeff(coeff.data))
    end
end

(::Type{T})(coeff::DCTCoeff) where T <: Signed = convert(T, coeff)

function Base.show(io::IO, coeff::DCTCoeff)
    show(IOContext(io, :compact=>true), MIME("text/plain"), coeff)
end

function Base.show(io::IO, ::MIME"text/plain", coeff::DCCoeff)
    if get(io, :typeinfo, Any) === DCCoeff
        print(io, Int(coeff))
    else
        compact = get(io, :compact, false)
        name = compact ? "DCCoeff" : string(typeof(coeff))
        print(io, name, "(", Int(coeff), ")")
    end
end

function Base.show(io::IO, ::MIME"text/plain", coeff::ACCoeff)
    s = coeff.symbol
    compact = get(io, :compact, false)
    name = compact ? "ACCoeff" : string(typeof(coeff))
    if iseob(s)
        print(io, name, "(1, EOB)")
    elseif iseobn(s)
        print(io, name, "(", Int(coeff), ", EOB", runlength(s), ")")
    elseif length(coeff) == 1
        print(io, name, "(", Int(coeff), ")")
    elseif compact
        print(io, name, "(", Int(coeff), ", ", runlength(s), ")")
    else
        print(io, length(coeff), "-element ", typeof(coeff), ":")
        for c in coeff
            println(io)
            print(io, " ", c)
        end
    end
end


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastLoop
using GeometricExamples

const Δt     = 5.
const ntime  = 10000
const nloop  = 200
const nsave  = 10
const nplot  = 10
const run_id = "poincare_1st_fast_vprk_dt5"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_1st.jl")


text:
function count_nucleotides(strand::AbstractString)
    counter = Dict('A' => 0, 'C' => 0, 'G' => 0, 'T' => 0)
    for sym in strand
        sym in keys(counter) || throw(DomainError(sym, "only A, C, G and T are valid nucleotides"))
        counter[sym] += 1
    end
    return counter
end


text:
## sample space/domain

"""
`F <: VariateForm` specifies the form of the variate or
dimension of a sample, univariate (scalar), multivariate (vector), matrix-variate (matrix).
"""
abstract type VariateForm end
struct Univariate    <: VariateForm end
struct Multivariate  <: VariateForm end
struct Matrixvariate <: VariateForm end

"""
`S <: ValueSupport` specifies the support of sample elements,
either discrete or continuous.
"""
abstract type ValueSupport end
struct Discrete   <: ValueSupport end
struct Continuous <: ValueSupport end

## Sampleable

"""
    Sampleable{F<:VariateForm,S<:ValueSupport}

`Sampleable` is any type able to produce random values.
Parametrized by a `VariateForm` defining the dimension of samples
and a `ValueSupport` defining the domain of possibly sampled values.
Any `Sampleable` implements the `Base.rand` method.
"""
abstract type Sampleable{F<:VariateForm,S<:ValueSupport} end

"""
    length(s::Sampleable)

The length of each sample. Always returns `1` when `s` is univariate.
"""
Base.length(s::Sampleable) = prod(size(s))
Base.length(::Sampleable{Univariate}) = 1
Base.length(s::Sampleable{Multivariate}) = throw(MethodError(length, (s,)))

"""
    size(s::Sampleable)

The size (i.e. shape) of each sample. Always returns `()` when `s` is univariate, and
`(length(s),)` when `s` is multivariate.
"""
Base.size(s::Sampleable)
Base.size(s::Sampleable{Univariate}) = ()
Base.size(s::Sampleable{Multivariate}) = (length(s),)

"""
    eltype(s::Sampleable)

The default element type of a sample. This is the type of elements of the samples generated
by the `rand` method. However, one can provide an array of different element types to
store the samples using `rand!`.
"""
Base.eltype(s::Sampleable{F,Discrete}) where {F} = Int
Base.eltype(s::Sampleable{F,Continuous}) where {F} = Float64

"""
    nsamples(s::Sampleable)

The number of values contained in one sample of `s`. Multiple samples are often organized
into an array, depending on the variate form.
"""
nsamples(t::Type{Sampleable}, x::Any)
nsamples(::Type{D}, x::Number) where {D<:Sampleable{Univariate}} = 1
nsamples(::Type{D}, x::AbstractArray) where {D<:Sampleable{Univariate}} = length(x)
nsamples(::Type{D}, x::AbstractVector) where {D<:Sampleable{Multivariate}} = 1
nsamples(::Type{D}, x::AbstractMatrix) where {D<:Sampleable{Multivariate}} = size(x, 2)
nsamples(::Type{D}, x::Number) where {D<:Sampleable{Matrixvariate}} = 1
nsamples(::Type{D}, x::Array{Matrix{T}}) where {D<:Sampleable{Matrixvariate},T<:Number} = length(x)

"""
    Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S}

`Distribution` is a `Sampleable` generating random values from a probability
distribution. Distributions define a Probability Distribution Function (PDF)
to implement with `pdf` and a Cumulated Distribution Function (CDF) to implement
with `cdf`.
"""
abstract type Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S} end

const UnivariateDistribution{S<:ValueSupport}   = Distribution{Univariate,S}
const MultivariateDistribution{S<:ValueSupport} = Distribution{Multivariate,S}
const MatrixDistribution{S<:ValueSupport}       = Distribution{Matrixvariate,S}
const NonMatrixDistribution = Union{UnivariateDistribution, MultivariateDistribution}

const DiscreteDistribution{F<:VariateForm}   = Distribution{F,Discrete}
const ContinuousDistribution{F<:VariateForm} = Distribution{F,Continuous}

const DiscreteUnivariateDistribution     = Distribution{Univariate,    Discrete}
const ContinuousUnivariateDistribution   = Distribution{Univariate,    Continuous}
const DiscreteMultivariateDistribution   = Distribution{Multivariate,  Discrete}
const ContinuousMultivariateDistribution = Distribution{Multivariate,  Continuous}
const DiscreteMatrixDistribution         = Distribution{Matrixvariate, Discrete}
const ContinuousMatrixDistribution       = Distribution{Matrixvariate, Continuous}

variate_form(::Type{Distribution{VF,VS}}) where {VF<:VariateForm,VS<:ValueSupport} = VF
variate_form(::Type{T}) where {T<:Distribution} = variate_form(supertype(T))

value_support(::Type{Distribution{VF,VS}}) where {VF<:VariateForm,VS<:ValueSupport} = VS
value_support(::Type{T}) where {T<:Distribution} = value_support(supertype(T))

# allow broadcasting over distribution objects
# to be decided: how to handle multivariate/matrixvariate distributions?
Broadcast.broadcastable(d::UnivariateDistribution) = Ref(d)

## TODO: the following types need to be improved
abstract type SufficientStats end
abstract type IncompleteDistribution end

const DistributionType{D<:Distribution} = Type{D}
const IncompleteFormulation = Union{DistributionType,IncompleteDistribution}

"""
    succprob(d::DiscreteUnivariateDistribution)

Get the probability of success.
"""
succprob(d::DiscreteUnivariateDistribution)

"""
    failprob(d::DiscreteUnivariateDistribution)

Get the probability of failure.
"""
failprob(d::DiscreteUnivariateDistribution)

# Temporary fix to handle RFunctions dependencies
"""
    @rand_rdist(::Distribution)

Mark a `Distribution` subtype as requiring RFunction calls. Since these calls
cannot accept an arbitrary random number generator as an input, this macro
creates new `rand(::Distribution, n::Int)` and
`rand!(::Distribution, X::AbstractArray)` functions that call the relevant
RFunction. Calls using another random number generator still work, but rely on
a quantile function to operate.
"""
macro rand_rdist(D)
    esc(quote
        function rand(d::$D, n::Int)
            [rand(d) for i in Base.OneTo(n)]
        end
        function rand!(d::$D, X::AbstractArray)
            for i in eachindex(X)
                X[i] = rand(d)
            end
            return X
        end
    end)
end


text:
struct UDPHeader
    src_port::UInt16
    dst_port::UInt16
    length::UInt16
    checksum::UInt16
end

struct UDPPacket
    header::UDPHeader
    payload::UnsafeArray{UInt8, 1}
end

function decode_udp(data::DenseVector{UInt8})
    p = Base.unsafe_convert(Ptr{UInt8}, data)
    GC.@preserve data begin
        h = unsafe_load(convert(Ptr{UDPHeader}, p))
    end
    h = UDPHeader(
        ntoh(h.src_port),
        ntoh(h.dst_port),
        ntoh(h.length),
        ntoh(h.checksum))
    UDPPacket(h, UnsafeArray{UInt8, 1}(p + 8, (Int(h.length - 8),)))
end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols