text:
using CUDA, CuYao
using DelimitedFiles
device!(parse(Int, ARGS[1]))

include("dataload.jl")

function load_and_contract(::Type{T}, n::Int; seed::Int, usecuda=false, datafile="ising.hdf5") where T
	loadeddata = HDF5.h5open(TropicalTensors.project_relative_path("data", datafile), "r")
    instance = loadeddata["n$n"]["seed$seed"]
    arrays = read(instance, "tensors")
    entros = read(instance, "entros")
    order = read(instance, "order")
    neighbors = read(instance, "labels")
    close(loadeddata)

    # build tensor network
	labels = tensor_labels(neighbors)
	tensors = map(1:n) do i
		arr = T.(permutedims(arrays[:,:,:,i], (3,2,1)))
		degen = T.(round.(exp.(permutedims(entros[:,:,:,i], (3,2,1)))))
		LabeledTensor(CountingTropical{T}.(arr, degen), labels[:,i])
	end
	# circle layout
	metas = [TensorMeta((0.5+0.5*cos(i/n*2π), 0.5+0.5*sin(i/n*2π)), string(i)) for i=1:n]
    tn = TensorNetwork(tensors, metas=metas)
    if usecuda
        tn = TropicalTensors.togpu(tn)
    end
    tree = build_tree(order)
    Array(TropicalTensors.contract(tn, tree).array)[]
end

function run(::Type{T}, n::Int; dataset) where T
    saveto = joinpath(@__DIR__, "$(dataset)_n$(n)_elsl.dat")
    elsl = zeros(T, 3, 100)
    datafile = dataset*".hdf5"
	loadeddata = HDF5.h5open(TropicalTensors.project_relative_path("data", datafile), "r")
    println("number of seeds = $(length(read(loadeddata["n$n"])))")
    t = @elapsed res = load_and_contract(T, n; seed=10, usecuda=true, datafile=dataset*".hdf5")
    for seed = 1:100
        try
            t = @elapsed res = load_and_contract(T, n; seed=seed, usecuda=true, datafile=datafile)
            println("seed = $seed, Δt = $t, maximum energy = $(res.n), degeneracy = $(res.c)")
            elsl[1,seed] = res.n/n
            elsl[2,seed] = log(res.c)/n
            elsl[3,seed] = t
        catch e
            println("Fail on seed $seed. $e")
        end
    end
    writedlm(saveto, elsl)
    return elsl
end

const n = parse(Int, ARGS[3])
@time run(Float32, n; dataset=ARGS[2])


text:
# utilities to handle scalars as uniform arrays
# (instead writing multiple implementations of functions, proxes, gradients)

get_kth_elem(n::N, k) where {N <: Number} = n
get_kth_elem(n::T, k) where {T} = n[k]

ext: jl


text:
# Wrap builders with logging
fns =  [:add!, :alloca!, :and!, :array_alloca!, :ashr!, :bitcast!, :br!, :call!, :dispose, :extract_value!,
        :fadd!, :fcmp!, :fdiv!, :fmul!, :fneg!, :fpext!, :fptosi!, :fptoui!, :fptrunc!, :frem!, :fsub!,
        :gep!, :globalstring_ptr!, :icmp!, :inttoptr!, :load!, :lshr!, :mul!, :neg!, :not!, :or!, :position,
        :position!, :ptrtoint!,
        :ret!, :sdiv!, :select!, :sext!, :shl!, :sitofp!, :srem!, :store!, :struct_gep!, :sub!,
        :trunc!, :udiv!, :uitofp!, :unreachable!, :urem!, :xor!, :zext!]
for f in fns
    @eval LLVM.$f(x::LoggingBuilder, args...) = (Base.@_inline_meta; O = LLVM.$f(x.builder, args...); @debug "LLVM output for $(x.name):" O; O)
end

export basedump
basedump(@nospecialize(fun), @nospecialize(argtypes)) =
    print(Base._dump_function(fun, argtypes, false, true, false, true, :att, false))

# Temporary workaround for an LLVM bug
function LLVM.ConstantInt(typ::LLVM.IntegerType, val::UInt64)
    bits = reinterpret(Culonglong, val)
    return LLVM.ConstantInt(LLVM.API.LLVMConstInt(LLVM.ref(typ), bits, convert(LLVM.Bool, false)))
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 86a927d5-0b05-4a7a-83c4-0f98fb15509e
begin
	using PlutoUI
	using Plots
end

# ╔═╡ 5d15bae3-d408-4f74-98a3-8b5e53562770
md"""
# Chapter 4
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ 435ece74-d26a-11eb-1bd1-5b65d058f3a1
md"""
## 4.1
トリボナッチ数列の第$N$項の値を求める再起関数を設計してください.
"""

# ╔═╡ 4a01190f-1e0e-4371-9847-e66eeab849b9
function f₁(n)
	if n == 0
		return 0
	elseif n == 1
		return 0
	elseif n == 2
		return 1
	end

	return f₁(n - 1) + f₁(n - 2) + f₁(n - 3)
end

# ╔═╡ 0224c51b-3028-45f0-bf6c-7423aac643a6
with_terminal() do
	for n ∈ 0:9
		println("n = $n: $(f₁(n))")
	end
end

# ╔═╡ aa51aa81-b8be-4e65-b259-e816655d01c9
plot(n -> (@timed f₁(n)).time, 1:33,
	title = "Computational Time of Naive Tribonacci",
	xlabel = "n",
	ylabel = "Seconds",
	legend = false
)

# ╔═╡ 5d8bedcf-6f20-4be3-961b-3b81ef486209
md"""
## 4.2
問題4.1で設計した再帰関数をメモ化によって効率化してください.
また, メモ化を実施後の計算量を評価してください.
"""

# ╔═╡ 1f806382-2995-4da2-b33d-c0403c530163
function f₂(n, memo = -ones(Int, n))
	if n == 0
		return 0
	elseif n == 1
		return 0
	elseif n == 2
		return 1
	end

	if memo[n] != -1
		return memo[n]
	end

	res = f₂(n - 1, memo) + f₂(n - 2, memo) + f₂(n - 3, memo)
	memo[n] = res
	return res
end

# ╔═╡ 645cdd0a-97c5-4a90-a70e-09a3571142db
with_terminal() do
	for n ∈ 0:9
		println("n = $n: $(f₂(n))")
	end
end

# ╔═╡ d65730a4-ec5f-436f-876b-fc98eb203107
md"""
計算量は$T_0, \dots, T_{N-1}$までメモ化されている時, $T_N$の計算量が$O(1)$であることから, $O(N)$.
"""

# ╔═╡ 3ac6e8f6-a3e9-4d81-a271-5457b20477bc
plot(n -> median((@timed f₂(n)).time for _ in 1:100), 1:3000,
	title = "Computational Time of Memorized Tribonacci",
	xlabel = "n",
	ylabel = "Seconds",
	legend = false
)

# ╔═╡ f6e73830-ad09-4718-868d-26f6a2c11b16
md"""
## 4.3
フィボナッチ数列の一般項が
$F_N = \frac{1}{\sqrt{5}}\left[\left(\frac{1 + \sqrt{5}}{2}\right)^N -
\left(\frac{1 - \sqrt{5}}{2}\right)^N \right]$
で表されることを示してください.
"""

# ╔═╡ 65503b27-14d4-42c8-abec-70060988d465
md"""
**答**
$F_0 = 1, F_1 = 1$は自明. 一般項を$F_{N-1}$まで仮定して,

$\begin{aligned}
\sqrt{5}(F_{N-1} + F_{N-2}) &= \left(1 + \frac{2}{1 + \sqrt{5}}\right)
\left(\frac{1 + \sqrt{5}}{2}\right)^{N-1} -
\left(1 + \frac{2}{1 - \sqrt{5}}\right)
\left(\frac{1 - \sqrt{5}}{2}\right)^{N-1} \\
&= \left(\frac{1 + \sqrt{5}}{2}\right)^N -
\left(\frac{1 - \sqrt{5}}{2}\right)^N.
\end{aligned}$
よって,
$F_N = F_{N - 1} + F_{N - 2}$.
帰納法から, 命題は示された.
"""

# ╔═╡ ae54b3e7-e04c-457f-839d-8aeecb6eba0d
md"""
## 4.4
コード 4.5で示したアルゴリズムの計算量が $O\left(\left(\frac{1 + \sqrt{5}}{2}\right)^N\right)$ で与えられることを示してください.
"""

# ╔═╡ e0cd89ce-9a49-4d51-9c36-1b467b3e4add
md"""
**答** コード4.5のアルゴリズムの計算量が$O(f(N))$であると仮定する.
定義から, $f(0) = f(1) = 1$は自明. また, Big Oの計算規則とコードの定義から,

$f(N) = f(N - 1) + f(N - 2).$

これはフィボナッチ数列と一致し, 一般項は問4.3で示された通りである.
支配項は$\left(\frac{1 + \sqrt{5}}{2}\right)^N$であるから, 命題は示された.
"""

# ╔═╡ 0303b656-d949-4365-b884-bf7bbc026a6d
md"""
## 4.6
部分和問題に対する再帰関数を用いる計算量$O(2^N)$のコード4.9に対してメモ化して, $O(NW)$ の計算量で動作するようにしてください.
"""

# ╔═╡ 14bcda24-9431-468a-ac5b-4bf6feb1e8e8
function f₆(as, w, i = length(as), memo = -ones(Int, length(as), w))
	if w < 0
		return false
	elseif w == 0
		return true
	elseif i == 0
		return false
	end

	if memo[i, w] == 1
		return true
	elseif memo[i, w] == 0
		return false
	end

	if f₆(as, w, i - 1, memo)
		memo[i, w] = 1
		return true
	elseif f₆(as, w - as[i], i - 1, memo)
		memo[i, w] = 1
		return true
	else
		memo[i, w] = 0
		return false
	end

end

# ╔═╡ 501824d0-017f-4d05-a477-cd7f4cef1df8
with_terminal() do
	println("Given a = [2, 4, 8]")
	for w ∈ 1:16
		println("w = $w: $(f₆([2, 4, 8], w))")
	end
end

# ╔═╡ Cell order:
# ╟─5d15bae3-d408-4f74-98a3-8b5e53562770
# ╠═86a927d5-0b05-4a7a-83c4-0f98fb15509e
# ╟─435ece74-d26a-11eb-1bd1-5b65d058f3a1
# ╠═4a01190f-1e0e-4371-9847-e66eeab849b9
# ╟─0224c51b-3028-45f0-bf6c-7423aac643a6
# ╟─aa51aa81-b8be-4e65-b259-e816655d01c9
# ╟─5d8bedcf-6f20-4be3-961b-3b81ef486209
# ╠═1f806382-2995-4da2-b33d-c0403c530163
# ╟─645cdd0a-97c5-4a90-a70e-09a3571142db
# ╟─d65730a4-ec5f-436f-876b-fc98eb203107
# ╟─3ac6e8f6-a3e9-4d81-a271-5457b20477bc
# ╟─f6e73830-ad09-4718-868d-26f6a2c11b16
# ╟─65503b27-14d4-42c8-abec-70060988d465
# ╟─ae54b3e7-e04c-457f-839d-8aeecb6eba0d
# ╟─e0cd89ce-9a49-4d51-9c36-1b467b3e4add
# ╟─0303b656-d949-4365-b884-bf7bbc026a6d
# ╠═14bcda24-9431-468a-ac5b-4bf6feb1e8e8
# ╟─501824d0-017f-4d05-a477-cd7f4cef1df8


text:
using CompScienceMeshes
using Test

m = meshrectangle(1.0,1.0,1.0,3)
p = chart(m, cells(m,i))

# bouding box of an array of points
c, s = boundingbox(p.vertices)
@test isa(c, eltype(p.vertices))
@test isa(s, Number)

# bouding box of a point
c, s = boundingbox(p.vertices[1])
@test c == p.vertices[1]
@test s == 0


text:
import Base: sum, push!

"""
This is a binary tree where the leaves are values
and the nodes are sums of those values. It is meant to make it
easier to find the leaf such that the sum of it and all previous
leaves is greater than a given value.
"""
struct PrefixSearchTree{T<:Real}
	array::Array{T,1}
	depth::Int
	offset::Int
end

"""
    PrefixSearchTree(ElementType::DataType, count::Int)

Constructor for an empty tree of type ElementType and length count.
"""
function PrefixSearchTree(ElementType::DataType, count::Int)
	new(zeros(ElementType, count), 0, 0)
end

"""
    PrefixSearchTree(summables)

Constructor of a prefix search tree from an iterable list of real numbers.
"""
function PrefixSearchTree(summables)
	a = Array(summables)
	@assert length(a) > 0
	depth = Int(ceil(log(2, length(a)))) + 1
	offset = 2^(depth - 1)
	b = zeros(eltype(summables), 2^depth - 1)
	b[offset:(length(a) + offset - 1)] = a
	#print("depth $depth offset $offset length $(length(b))\n")
	pst = PrefixSearchTree{eltype(summables)}(b, depth, offset)
	calculate_prefix!(pst)
	pst
end

"""
    choose(pst::PrefixSearchTree, value)

Find the minimum index such that the prefix is greater than the given
value.

Precondition: The value must be strictly less than the total for
the tree.
"""
function choose(pst::PrefixSearchTree, value)
	if value >= pst.array[1]
		error("Value $value not less than total $(pst.array[1])")
	end

	requested_value = value
    index = 1
    for level = 0:(pst.depth - 2)
        left_child = 2 * index
        if pst.array[left_child] > value
            index = left_child
        else
            index = left_child + 1
            value -= pst.array[left_child]
        end
    end
    index - pst.offset + 1, pst.array[index] - value
end

sum(pst::PrefixSearchTree) = pst.array[1]

"""
If there are multiple values to enter, then present them
at once as pairs of tuples, (index, value).
"""
function push!(pst::PrefixSearchTree, pairs)
	modify=Set{Int}()
	for (pos, value) in pairs
		index = pos + pst.offset - 1
		pst.array[index] = value
		push!(modify, div(index, 2))
	end

	# everything at depth-1 is correct, and changes are in modify.
    for depth = (pst.depth - 2):-1:0
        parents = Set{Int}()
        for node_idx in modify
            pst.array[node_idx]=pst.array[2 * node_idx] + pst.array[2 * node_idx + 1]
            push!(parents, div(node_idx, 2))
        end
        modify = parents
    end
end

function push!(pst::PrefixSearchTree{T}, index::Int, value::T) where T
    push!(pst, [(index, value)])
end

# Private
function calculate_prefix!(pst::PrefixSearchTree)
	for d = (pst.depth - 2):-1:0
		for node_idx = (2^d):(2^(d + 1) - 1)
			pst.array[node_idx] = pst.array[2 * node_idx] + pst.array[2 * node_idx + 1]
		end
	end
end


text:
module PoincareInvariant2ndTest

    using GeometricIntegrators
    using PoincareInvariants
    using SymPy

    const Δt = 10.
    const B₀ = 1.
    const r₀ = 0.5
    const z₀ = 0.0
    const z₁ = 0.1
    const u₀ = 5E-1
    const u₁ = 5E-2
    const kz = 2
    const ku = 4

    function B(t, q)
        B₀ * (1 + q[1]^2 + q[2]^2)
    end

    function gcϑ(t, q, ϑ)
        ϑ[1] =-q[2] * (1 + B(t,q)) / 4
        ϑ[2] =+q[1] * (1 + B(t,q)) / 4
        ϑ[3] = q[4]
        ϑ[4] = 0
        nothing
    end

    function gcϑ(t, q)
        ϑ = zero(q)
        gcϑ(t, q, ϑ)
        return ϑ
    end

    function gcω(t, q, Ω)
        Ω[1,1] = 0
        Ω[1,2] =-B(t,q)
        Ω[1,3] = 0
        Ω[1,4] = 0

        Ω[2,1] =+B(t,q)
        Ω[2,2] = 0
        Ω[2,3] = 0
        Ω[2,4] = 0

        Ω[3,1] = 0
        Ω[3,2] = 0
        Ω[3,3] = 0
        Ω[3,4] =+q[4]

        Ω[4,1] = 0
        Ω[4,2] = 0
        Ω[4,3] =-q[4]
        Ω[4,4] = 0

        nothing
    end

    function gc_surface_q(s,t)
        x  = r₀*(s-0.5)
        y  = r₀*(t-0.5)
        z  = z₀ + z₁ * cos(2π*kz*s) * cos(2π*kz*t)
        u  = u₀ + u₁ * sin(2π*ku*s) * sin(2π*ku*t)

        [x, y, z, u]
    end

    function gc_surface_p(s,t)
        gcϑ(0, gc_surface_q(s,t))
    end

    function gc_dummy(t, a, b, c)
        nothing
    end

    function gc_dummy_pode(q₀, p₀)
        PODE(gc_dummy, gc_dummy, q₀, p₀)
    end

    function gc_dummy_iode(q₀)
        p₀ = [gcϑ(0,q) for q in q₀]
        IODE(gc_dummy, gc_dummy, gc_dummy, q₀, p₀; v̄=gc_dummy)
    end

    function compute_canonical_invariant_approxfun(nx, ny, nt=0)
        pinv = PoincareInvariant2ndCanonical(gc_dummy_pode, gc_surface_q, gc_surface_p, Δt, 4, nx, ny, nt)
        sol  = Solution(pinv.equ, Δt, nt)

        I = evaluate_poincare_invariant(pinv, sol)

        return I[0]
    end

    function compute_invariant_approxfun(nx, ny, nt=0)
        pinv = PoincareInvariant2ndApproxFun(gc_dummy_iode, gc_surface_q, gcω, Δt, 4, nx, ny, nt)
        sol  = Solution(pinv.equ, Δt, nt)

        I, J, ΔI, ΔJ = evaluate_poincare_invariant(pinv, sol)

        return I[0], J[0]
    end

    function compute_invariant_cop(nx, ny, nt=0)
        pinv = PoincareInvariant2ndCOP(gc_dummy_iode, gc_surface_q, gcω, Δt, 4, nx, ny, nt)
        sol  = Solution(pinv.equ, Δt, nt)

        I, J, ΔI, ΔJ = evaluate_poincare_invariant(pinv, sol)

        return I[0], J[0]
    end

    function compute_invariant_trapezoidal(nx, ny, nt=0)
        pinv = PoincareInvariant2ndTrapezoidal(gc_dummy_iode, gc_surface_q, gcω, Δt, 4, nx, ny, nt)
        sol  = Solution(pinv.equ, Δt, nt)

        I, J, ΔI, ΔJ = evaluate_poincare_invariant(pinv, sol)

        return I[0], J[0]
    end

    function compute_invariant_analytical()
        local verbosity = get_config(:verbosity)

        if get_config(:verbosity) > 1
            println()
            println("Second Euler-Poincaré Integral Invariant (Analytical)")
            println("=====================================================")
            println()
        end

        s, t = Sym("s, t")
        x, y, z, u = gc_surface_q(s,t)

        S = simplify(B₀ * (1 + x^2 + y^2) * ( diff(x, s) * diff(y, t) + diff(x, t) * diff(y, s) )
                   + u * ( diff(z, s) * diff(u, t) + diff(z, t) * diff(u, s) ) )
        I = N(SymPy.integrate(S, (s, 0, 1), (t, 0, 1)), 16)

        verbosity > 1 ? println("           I = ", I) : nothing

        return I
    end

    export compute_canonical_invariant_approxfun,
           compute_invariant_approxfun,
           compute_invariant_cop,
           compute_invariant_trapezoidal,
           compute_invariant_analytical

end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames: DataFrame
using CUDA
using Random
using Suppressor

include("../helper/gpu.jl")

@testset "Marginals" begin
    EPS = 1e-5

    prob_circuit = zoo_psdd("little_4var.psdd");
    @test prob_circuit(false, false, false, missing) ≈ -1.0498221

    data_marg = DataFrame([false false false false;
                      false true true false;
                      false false true true;
                      false false false missing;
                      missing true false missing;
                      missing missing missing missing;
                      false missing missing missing])
    true_prob = [0.07; 0.03; 0.13999999999999999;
                    0.3499999999999; 0.1; 1.0; 0.8]

    calc_prob = exp.(MAR(prob_circuit, data_marg))
    @test true_prob ≈ calc_prob atol=EPS
    @test marginal_log_likelihood_avg(prob_circuit, data_marg) ≈ sum(log.(true_prob))/7
    marginal_all_result = marginal_all(prob_circuit, data_marg);

    marginal_all_true_answer = [0.0 -Inf -Inf -Inf -Inf -Inf 0.0 0.0 0.0 0.0 -0.356675 -2.30259 -2.65926
        0.0 -Inf -Inf 0.0 0.0 -Inf 0.0 -Inf -Inf 0.0 -2.30259 -1.20397 -3.50656
        0.0 -Inf -Inf -Inf 0.0 0.0 0.0 0.0 -Inf -Inf -0.356675 -1.60944 -1.96611
        0.0 -Inf -Inf -Inf -Inf 0.0 0.0 0.0 0.0 0.0 -0.356675 -0.693147 -1.04982
        0.0 -Inf 0.0 0.0 -Inf 0.0  0.0 -Inf 0.0 0.0 -1.60944 -0.693147 -2.30259
        0.0 -Inf 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.98023f-8 -7.45058f-9 -3.72529f-8
        0.0 -Inf -Inf 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.223144 -7.45058f-9 -0.223144];

    for pair in zip(marginal_all_result, marginal_all_true_answer)
        @test pair[1] ≈ pair[2]  atol=EPS
    end

    cpu_gpu_agree_approx(data_marg) do d
        marginal_all(prob_circuit, d)
    end

    function test_complete_mar(circ, data, weights = nothing, atol = 1e-6)
        r1 = isnothing(weights) ? EVI(circ, data) : EVI(circ, data, weights)
        @test isgpu(data) == isgpu(r1)
        r2 = isnothing(weights) ? MAR(circ, data) : MAR(circ, data, weights)
        @test isgpu(data) == isgpu(r2)
        @test r1 ≈ r2 atol = atol
    end

    data_full = generate_data_all(num_variables(prob_circuit))

    test_complete_mar(prob_circuit, data_full)
    CUDA.functional() && test_complete_mar(prob_circuit, to_gpu(data_full))

    cpu_gpu_agree_approx(data_full) do d
        marginal_all(prob_circuit, d)
    end

    # make sure log-likelihoods are -Inf when the input is not satisfied
    data = DataFrame([false true false missing;
                      false true true false;
                      missing missing missing false])
    alltrue = multiply(pos_literals(ProbCircuit,4))
    @test all(MAR(alltrue, data) .== -Inf)

    cpu_gpu_agree(data) do d
        MAR(alltrue, d)
    end

    # Strudel Marginal Flow Test
    rng = MersenneTwister(100003); # Fix the seed
    samples, _ = sample(prob_circuit, 100000; rng)
    mix, weights, _ = learn_strudel(DataFrame(convert(BitArray, samples)); num_mix = 10,
                                    init_maxiter = 20, em_maxiter = 100, verbose = false)
    mix_calc_prob = exp.(MAR(mix, data_marg, weights))
    for mix_pair in zip(true_prob, mix_calc_prob)
        @test mix_pair[1] ≈ mix_pair[2]  atol=0.1
    end

    test_complete_mar(mix, data_full, weights, 0.1)
end

@testset "Marginals batch" begin
    prob_circuit = zoo_psdd("little_4var.psdd");
    @test prob_circuit(false, false, false, missing) ≈ -1.0498221

    data_marg = DataFrame([false false false false;
                      false true true false;
                      false false true true;
                      false false false missing;
                      missing true false missing;
                      missing missing missing missing;
                      false missing missing missing])
    batched_data_marg = batch(data_marg, 1; shuffle = false)

    data_marg = DataFrame([false false false false;
                      false true true false;
                      false false true true;
                      false false false missing;
                      missing true false missing;
                      missing missing missing missing;
                      false missing missing missing])
    weights = [0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6]
    data_marg_w = weigh_samples(data_marg, weights)
    batched_data_marg_w = batch(data_marg_w, 1)

    true_prob = [0.07; 0.03; 0.13999999999999999;
                    0.3499999999999; 0.1; 1.0; 0.8]

    data_marg = DataFrame([false false false false;
                      false true true false;
                      false false true true;
                      false false false missing;
                      missing true false missing;
                      missing missing missing missing;
                      false missing missing missing])

    calc_prob = exp.(MAR(prob_circuit, data_marg))
    @test true_prob ≈ calc_prob atol=1e-7
    @test marginal_log_likelihood(prob_circuit, batched_data_marg) ≈ sum(log.(true_prob))
    @test marginal_log_likelihood_avg(prob_circuit, batched_data_marg) ≈ sum(log.(true_prob)) / 7
    @test marginal_log_likelihood(prob_circuit, batched_data_marg_w) ≈ sum(log.(true_prob)) * 0.6
    @test marginal_log_likelihood_avg(prob_circuit, batched_data_marg_w) ≈ sum(log.(true_prob)) / 7

    mar = marginal(prob_circuit, batched_data_marg)
    @test true_prob ≈ exp.(mar) atol=1e-7
end

@testset "Marginal flows" begin

    prob_circuit = zoo_psdd("little_4var.psdd");

    function test_flows(data)
        # Comparing with down pass with fully observed data

        data_f = CUDA.@allowscalar Float64.(data)

        _, f1, _ = satisfies_flows(prob_circuit, data_f)
        _, f2, _ = marginal_flows(prob_circuit, data)

        # note: while downward pass flows should be the same,
        # the upward pass is *not* supposed to be the same (parameters used vs not)

        f1 = to_cpu(f1[:,3:end]) # ignore true and false leaf
        f2 = to_cpu(f2[:,3:end]) # ignore true and false leaf

        @test f1 ≈ exp.(f2) atol=1e-6
    end

    data_full = generate_data_all(num_variables(prob_circuit))

    test_flows(data_full)
    CUDA.functional() && test_flows(to_gpu(data_full))

    cpu_gpu_agree_approx(data_full) do d
        _, f, _ = marginal_flows(prob_circuit, d)
        f[:,3:end] # ignore true and false leaf
    end

    # Validating one example with missing features done by hand
    data_partial = DataFrame([missing true missing true;])
    prob_circuit = zoo_psdd("little_4var.psdd");
    _, f, _ = marginal_flows(prob_circuit, data_partial)
    f = exp.(f)

    @test f[end] ≈ 1.0
    @test f[end-1] ≈ 1.0
    @test f[end-2] ≈ 1.0
    @test f[end-4] ≈ 2/3
    @test f[end-5] ≈ 0.0 atol=1e-7
    @test f[end-6] ≈ 1/2
    @test f[end-7] ≈ 1.0
    @test f[end-8] ≈ 1/3
    @test f[end-9] ≈ 1
    @test f[end-10] ≈ 1/2

    # correctness on gpu by transitivy with above test

end


text:
##################################################
# Structure which gathers the quantities used in computing K and dK
##################################################
"""
Structure which gathers all the necessary arrays needed to compute the geometrical factors (Knn and its derivatives).

# Summary

struct IntTable <: Any

# Fields
- `tabR      :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table for the radius of the first particle.
- `tabg      :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table for the prefactor of the first particle.
- `tabRp     :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table for the radius of the second particle.
- `tabgp     :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table for the prefactor of the second particle.
- `tabw      :: MArray{Tuple{nbK+2},Int64,1,nbK+2}`: Table for the sorting of the arrays.
- `tabh1     :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table for the h1 of the first particle.
- `tabh2     :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table for the h2 of the second particle.
- `tabdeltaP :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table of the partial sum deltaP.
- `tabdeltaQ :: MArray{Tuple{nbK},Float64,1,nbK}`  : Table of the partial sum deltaQ.

"""
struct IntTable
    tabR::MArray{Tuple{nbK},Float64,1,nbK}	# Table for the radius of the first particle
    tabg::MArray{Tuple{nbK},Float64,1,nbK}	# Table for the prefactor of the first particle
    tabRp::MArray{Tuple{nbK},Float64,1,nbK}	# Table for the radius of the second particle
    tabgp::MArray{Tuple{nbK},Float64,1,nbK}	# Table for the prefactor of the second particle
    tabw::MArray{Tuple{nbK+2},Int64,1,nbK+2}	# Table for the sorting of the arrays
    tabh1::MArray{Tuple{nbK},Float64,1,nbK}	# Table for the h1 of the first particle
    tabh2::MArray{Tuple{nbK},Float64,1,nbK}	# Table for the h2 of the second particle
    tabdeltaP::MArray{Tuple{nbK},Float64,1,nbK}	# Table of the partial sum deltaP
    tabdeltaQ::MArray{Tuple{nbK},Float64,1,nbK} # Table of the partial sum deltaQ
end
##################################################
# Creation of an IntTable structure
##################################################

"""
    IntTable_create!(t1=tab1, t2=tab2)

Creates a IntTable structure with all fields set as zeros-arrays.
"""
function IntTable_create!()
    return IntTable(@MVector(zeros(Float64,nbK)),@MVector(zeros(Float64,nbK)),@MVector(zeros(Float64,nbK)),
           @MVector(zeros(Float64,nbK)),@MVector( zeros(Int64,nbK+2)),@MVector(zeros(Float64,nbK)),
           @MVector(zeros(Float64,nbK)),@MVector(zeros(Float64,nbK)),@MVector(zeros(Float64,nbK)))
end
##################################################
# Initialization of a structure
##################################################

"""
    IntTable_init!(Table::IntTable, <keyword arguments>)

Initializes the IntTable structure `Table` at the given arguments.

# Arguments
- `n  ::Int64`  : resonance number of the first orbit.
- `np ::Int64`  : resonance number of the second orbit.
- `a  ::Float64`: semi-major axis of the first orbit.
- `j  ::Float64`: reduced angular momentum of the first orbit.
- `ap ::Float64`: semi-major axis of the decond orbit.
- `jp ::Float64`: reduced angular momentum of the second orbit.

"""
function IntTable_init!(Table::IntTable,n::Int64,np::Int64,a::Float64,j::Float64,ap::Float64,jp::Float64)
    for k=1:nbK # Loop over the sample points
        cosf = tabcosf_Klnnp[k] # Extracting the value of cos(f)
        cosnf, cosnpf = tabcosnf_Klnnp[abs(n),k], tabcosnf_Klnnp[abs(np),k] # Extracting the value of cos(n*f) and cos(np*f) !! ATTENTION, we use (abs(n),abs(np)) for the array access to be ok.
        R, Rp = a*(j^(2))/(1.0+sqrt(1.0-j^(2))*cosf), ap*(jp^(2))/(1.0+sqrt(1.0-jp^(2))*cosf) # Radii of the first and second particles
        dMdf, dMpdfp = (R/a)^(2)/(j), (Rp/ap)^(2)/(jp) # Computing the Jacobian dM/df and dMp/dfp
        dM = ddMdfdj(j,cosf) # Computing j-derivative of the Jacobian
        rr = reduced_rad(j,cosf) # Computing reduced radius
        drr = Dreduced_rad(j,cosf) # Computing j-derivative of reduced radius
        #####
        Table.tabR[k] , Table.tabg[k]  = R , dMdf   * cosnf
        Table.tabRp[k], Table.tabgp[k] = Rp, dMpdfp * cosnpf
        Table.tabh1[k], Table.tabh2[k] = dM*cosnf, rr*drr/j *cosnf #h1, h2
    end
    Table.tabw[1] = 0 # Initial term. !! ATTENTION, indices are shifted by one
    icount = 0 # Temporary counter for filling on tabw
    for k=1:nbK
        while (icount < nbK && Table.tabR[icount+1] <= Table.tabRp[k]) # !! ATTENTION, to the order of the tests to avoid segmentation faults
            icount += 1 # Updating the counter
        end
        Table.tabw[k+1] = icount # ATTENTION, indices are shifted by one
    end
    Table.tabw[nbK+2] = nbK # Last term. !! ATTENTION, indices are shifted by one
end


text:
"""
    basis = basis_rotation(Ω, i, ω)
    basis = basis_rotation(xvec, yvec)
    basis_rotation() = one(RotMatrix{3, Float64})

From the RAAN `Ω`, inclination `i`, and argument of the periapsis `ω`, computes the rotation matrix that transforms
an orbit's perifocal plane to the inertial frame. For this matrix, the first column is a unit vector that points toward
the periapsis of the orbit, and the third column is a unit vector that points in the normal direction.

If two vectors `xvec` and `yvec` are provided, a rotation matrix is constructed that describes the plane the two vectors
define, where `normalize(xvec)` is the first basis vector.
"""
function basis_rotation(Ω, i, ω)
    # Ω, ω, i = Float64(Ω), Float64(ω), Float64(i)
    # R₁ = [cos(Ω) -sin(Ω)  0.;    # rotate around z-axis to match Ω
    #       sin(Ω)  cos(Ω)  0.;
    #       0.      0.      1.]
    # R₂ = [1.      0.      0.;    # rotate around x-axis to match i
    #       0.  cos(i) -sin(i);
    #       0.  sin(i)  cos(i)]
    # R₃ = [cos(ω) -sin(ω)  0.;    # rotate around z-axis to match ω
    #       sin(ω)  cos(ω)  0.;
    #       0.      0.      1.]
    # return R₃*R₂*R₁
    return RotZXZ(Ω, i, ω)
end

function basis_rotation(xvec::AbstractVector{<:Real}, yvec::AbstractVector{<:Real})
    @assert length(xvec) == length(yvec) == 3
    xbasisvec = normalize(xvec)
    ybasisvec = normalize(yvec - xbasisvec .* dot(xbasisvec, yvec))
    zbasisvec = cross(xbasisvec, ybasisvec)
    return RotMatrix(hcat(xbasisvec,ybasisvec,zbasisvec))
end

basis_rotation() = one(RotMatrix{3, Float64})

# methods for switching orbital reference frames

"""
    pvec = inertial_to_perifocal_bases(vec, basis)
    pvec = inertial_to_perifocal_bases(vec, orb)

Transforms a vector `vec` from the inertial frame to an orbit's perifocal frame,
defined by `basis`.
"""
inertial_to_perifocal_bases(vec::AbstractVector{<:Real}, basis::Rotation{3,Float64}) = basis\vec
inertial_to_perifocal_bases(vec::AbstractVector{<:Real}, orb::Orbit) = inertial_to_perifocal_bases(vec, orb.basis)

"""
    ivec = perifocal_to_inertial_bases(vec, basis)
    ivec = perifocal_to_inertial_bases(vec, orb)

Transforms a vector `vec` from an orbit's perifocal frame, defined by `basis`,
to the inertial frame.
"""
perifocal_to_inertial_bases(vec::AbstractVector{<:Real}, basis::Rotation{3,Float64}) = basis*vec
perifocal_to_inertial_bases(vec::AbstractVector{<:Real}, orb::Orbit) = perifocal_to_inertial_bases(vec, orb.basis)

"""
    rotmat = align_vectors(vec, vecref)

Computes the rotation matrix to align `vec` with `vecref`.
"""
function align_vectors(vec::AbstractVector{<:Real}, vecref::AbstractVector{<:Real})
    nvec, nvecref = normalize(vec), normalize(vecref)
    axis = normalize(cross(nvec, nvecref))
    angle = wrap_acos(dot(nvec, nvecref))
    if angle == 0
        return one(RotMatrix{3, Float64})
    end
    return AngleAxis(angle, axis...)
end


text:
using FinancialMonteCarlo

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 300;
sigma = 0.2;
variate_ = FinancialMonteCarlo.ControlVariates(Forward(T), MonteCarloConfiguration(1000, 100))
variate_2 = FinancialMonteCarlo.ControlVariates(AsianFloatingStrikeOption(T), MonteCarloConfiguration(1000, 100))
mc = MonteCarloConfiguration(Nsim, Nstep, variate_);
mc2 = MonteCarloConfiguration(Nsim, Nstep, variate_2);
mc1 = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;

spotData1 = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

AsianPrice1 = pricer(Model, spotData1, mc, AsianFloatingStrikeData);
AsianPrice2 = pricer(Model, spotData1, mc, AsianFixedStrikeData);
AsianPrice2 = pricer(Model, spotData1, mc2, AsianFixedStrikeData);


text:
include("consts.jl")
include("parseutils.jl")

struct URLParsingError <: Exception
    msg::String
end
Base.show(io::IO, p::URLParsingError) = println(io, "HTTP.URLParsingError: ", p.msg)

@enum(http_host_state,
    s_http_host_dead,
    s_http_userinfo_start,
    s_http_userinfo,
    s_http_host_start,
    s_http_host_v6_start,
    s_http_host,
    s_http_host_v6,
    s_http_host_v6_end,
    s_http_host_v6_zone_start,
    s_http_host_v6_zone,
    s_http_host_port_start,
    s_http_host_port,
)

const blank = SubString("blank", 1, 0)

# url parsing
function parseurlchar(s, ch::Char, strict::Bool)
    @anyeq(ch, ' ', '\r', '\n') && return s_dead
    strict && (ch == '\t' || ch == '\f') && return s_dead

    if s == s_req_spaces_before_target || s == s_req_target
        (ch == '/') && return s_req_path
        isalpha(ch) && return s_req_schema
    elseif s == s_req_schema
        isalphanum(ch) && return s
        ch == ':' && return s_req_schema_slash
    elseif s == s_req_schema_slash
        ch == '/' && return s_req_schema_slash_slash
        isurlchar(ch) && return s_req_path
    elseif s == s_req_schema_slash_slash
        ch == '/' && return s_req_server_start
        isurlchar(ch) && return s_req_path
    elseif s == s_req_server_with_at
        ch == '@' && return s_dead
        ch == '/' && return s_req_path
        ch == '?' && return s_req_query_string_start
        (isuserinfochar(ch) || ch == '[' || ch == ']') && return s_req_server
    elseif s == s_req_server_start || s == s_req_server
        ch == '/' && return s_req_path
        ch == '?' && return s_req_query_string_start
        ch == '@' && return s_req_server_with_at
        (isuserinfochar(ch) || ch == '[' || ch == ']') && return s_req_server
    elseif s == s_req_path
        (isurlchar(ch) || ch == '@') && return s
        ch == '?' && return s_req_query_string_start
        ch == '#' && return s_req_fragment_start
    elseif s == s_req_query_string_start || s == s_req_query_string
        isurlchar(ch) && return s_req_query_string
        ch == '?' && return s_req_query_string
        ch == '#' && return s_req_fragment_start
    elseif s == s_req_fragment_start
        isurlchar(ch) && return s_req_fragment
        ch == '?' && return s_req_fragment
        ch == '#' && return s_req_fragment
    elseif s == s_req_fragment
        isurlchar(ch) && return s
        (ch == '?' || ch == '#') && return s
    end
    #= We should never fall out of the switch above unless there's an error =#
    return s_dead
end

function http_parse_host_char(s::http_host_state, ch)
    if s == s_http_userinfo || s == s_http_userinfo_start
        ch == '@' && return s_http_host_start
        isuserinfochar(ch) && return s_http_userinfo
    elseif s == s_http_host_start
        ch == '[' && return s_http_host_v6_start
        ishostchar(ch) && return s_http_host
    elseif s == s_http_host
        ishostchar(ch) && return s_http_host
        ch == ':' && return s_http_host_port_start
    elseif s == s_http_host_v6_end
        ch == ':' && return s_http_host_port_start
    elseif s == s_http_host_v6
        ch == ']' && return s_http_host_v6_end
        (ishex(ch) || ch == ':' || ch == '.') && return s_http_host_v6
        s == s_http_host_v6 && ch == '%' && return s_http_host_v6_zone_start
    elseif s == s_http_host_v6_start
        (ishex(ch) || ch == ':' || ch == '.') && return s_http_host_v6
        s == s_http_host_v6 && ch == '%' && return s_http_host_v6_zone_start
    elseif s == s_http_host_v6_zone
        ch == ']' && return s_http_host_v6_end
        (isalphanum(ch) || @anyeq(ch, '%', '.', '-', '_', '~')) && return s_http_host_v6_zone
    elseif s == s_http_host_v6_zone_start
        (isalphanum(ch) || @anyeq(ch, '%', '.', '-', '_', '~')) && return s_http_host_v6_zone
    elseif s == s_http_host_port || s == s_http_host_port_start
        isnum(ch) && return s_http_host_port
    end
    return s_http_host_dead
end

function http_parse_host(host::SubString, foundat=false)

    host1 = port1 = userinfo1 = 1
    host2 = port2 = userinfo2 = 0
    s = ifelse(foundat, s_http_userinfo_start, s_http_host_start)

    for i in eachindex(host)
        @inbounds p = host[i]

        new_s = http_parse_host_char(s, p)
        if new_s == s_http_host_dead
            throw(URLParsingError("encountered invalid host character: \n" *
                                  "$host\n$(lpad("", i-1, "-"))^"))
        end
        if new_s == s_http_host
            if s != s_http_host
                host1 = i
            end
            host2 = i

        elseif new_s == s_http_host_v6
            if s != s_http_host_v6
                host1 = i
            end
            host2 = i

        elseif new_s == s_http_host_v6_zone_start ||
               new_s == s_http_host_v6_zone
            host2 = i

        elseif new_s == s_http_host_port
            if s != s_http_host_port
                port1 = i
            end
            port2 = i

        elseif new_s == s_http_userinfo
            if s != s_http_userinfo
                userinfo1 = i
            end
            userinfo2 = i
        end
        s = new_s
    end
    if @anyeq(s, s_http_host_start, s_http_host_v6_start, s_http_host_v6,
                 s_http_host_v6_zone_start, s_http_host_v6_zone,
                 s_http_host_port_start, s_http_userinfo, s_http_userinfo_start)
        throw(URLParsingError("ended in unexpected parsing state: $s"))
    end

    return SubString(host, host1, host2),
           SubString(host, port1, port2),
           SubString(host, userinfo1, userinfo2)
end

http_parser_parse_url(url::AbstractString) = http_parser_parse_url(String(url))

function http_parser_parse_url(url::String)

    s = s_req_spaces_before_target

    uf = -1
    old_uf = -1
    off1 = off2 = 0
    foundat = false

    empty = SubString(url, 1, 0)
    scheme = userinfo = host = port = path = query = fragment = empty

    mask = 0x00
    end_i = endof(url)
    for i in eachindex(url)
        @inbounds p = url[i]
        olds = s
        s = parseurlchar(s, p, false)
        if s == s_dead
            throw(URLParsingError(
                "encountered invalid url character for parsing state = " *
                "$(ParsingStateCode(olds)):\n$url)\n$(lpad("", i-1, "-"))^"))
        elseif @anyeq(s, s_req_schema_slash,
                         s_req_schema_slash_slash,
                         s_req_server_start,
                         s_req_query_string_start,
                         s_req_fragment_start)
            if s == s_req_server_start
                host = blank
            end
            if i != end_i
                continue
            else
                uf = s
            end
        elseif s == s_req_server_with_at
            foundat = true
            uf = s_req_server
        elseif @anyeq(s, s_req_schema,
                         s_req_server_with_at,
                         s_req_server,
                         s_req_path,
                         s_req_query_string,
                         s_req_fragment)
            uf = s
        else
            throw(URLParsingError("ended in unexpected parsing state: $s\n$url"))
        end
        if uf == old_uf
            off2 = i
            if i != end_i
                continue
            end
        end

        @label save_part
        if old_uf != -1
            part = SubString(url, off1, off2)
            old_uf == s_req_schema             && (scheme = part)
            old_uf == s_req_server             && (host = part)
            old_uf == s_req_path               && (path = part)
            old_uf == s_req_query_string       && (query = part)
            old_uf == s_req_fragment           && (fragment = part)
            old_uf == s_req_query_string_start && (query = blank)
            old_uf == s_req_fragment_start     && (fragment = blank)
        end

        off1 = i
        off2 = i
        if i == end_i && uf != old_uf
            old_uf = uf
            @goto save_part
        end
        old_uf = uf
    end
    if !isempty(scheme) && isempty(host) && isempty(path)
        throw(URLParsingError("URI must include host or path with scheme\n$url"))
    end
    if !isempty(host)
        host, port, userinfo = http_parse_host(host, foundat)
        if foundat && isempty(userinfo)
            userinfo = blank
        end
    end
    return URI(url, scheme, userinfo, host, port, path, query, fragment)
end

const normal_url_char = Bool[
#=   0 nul    1 soh    2 stx    3 etx    4 eot    5 enq    6 ack    7 bel  =#
        false,   false,   false,   false,   false,   false,   false,   false,
#=   8 bs     9 ht    10 nl    11 vt    12 np    13 cr    14 so    15 si   =#
        false,   true,   false,   false,   true,   false,   false,   false,
#=  16 dle   17 dc1   18 dc2   19 dc3   20 dc4   21 nak   22 syn   23 etb =#
        false,   false,   false,   false,   false,   false,   false,   false,
#=  24 can   25 em    26 sub   27 esc   28 fs    29 gs    30 rs    31 us  =#
        false,   false,   false,   false,   false,   false,   false,   false,
#=  32 sp    33  !    34  "    35  #    36  $    37  %    38  &    39  '  =#
        false,   true,   true,   false,   true,   true,   true,  true,
#=  40  (    41  )    42  *    43  +    44  ,    45  -    46  .    47  /  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#=  48  0    49  1    50  2    51  3    52  4    53  5    54  6    55  7  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#=  56  8    57  9    58  :    59  ;    60  <    61  =    62  >    63  ?  =#
        true,   true,   true,   true,   true,   true,   true,   false,
#=  64  @    65  A    66  B    67  C    68  D    69  E    70  F    71  G  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#=  72  H    73  I    74  J    75  K    76  L    77  M    78  N    79  O  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#=  80  P    81  Q    82  R    83  S    84  T    85  U    86  V    87  W  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#=  88  X    89  Y    90  Z    91  [    92  \    93  ]    94  ^    95  _  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#=  96  `    97  a    98  b    99  c   100  d   101  e   102  f   103  g  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#= 104  h   105  i   106  j   107  k   108  l   109  m   110  n   111  o  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#= 112  p   113  q   114  r   115  s   116  t   117  u   118  v   119  w  =#
        true,   true,   true,   true,   true,   true,   true,  true,
#= 120  x   121  y   122  z   123  {   124,   125  }   126  ~   127 del =#
        true,   true,   true,   true,   true,   true,   true,   false,
]

@inline isurlchar(c) =  c > '\u80' ? true : normal_url_char[Int(c) + 1]


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

mutable struct PacketCapture <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, PacketCaptureParameters } # spec name: properties

    function PacketCapture(;properties=nothing)
        o = new()
        validate_property(PacketCapture, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type PacketCapture

const _property_map_PacketCapture = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"))
const _property_types_PacketCapture = Dict{Symbol,String}(Symbol("properties")=>"PacketCaptureParameters")
Base.propertynames(::Type{ PacketCapture }) = collect(keys(_property_map_PacketCapture))
Swagger.property_type(::Type{ PacketCapture }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_PacketCapture[name]))}
Swagger.field_name(::Type{ PacketCapture }, property_name::Symbol) =  _property_map_PacketCapture[property_name]

function check_required(o::PacketCapture)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

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


text:
# Cross-check against data published in Frost 2017 PRB.

@testset "FrostPolaronMobility2017" begin

T=300

#Ts,Kμs, Hμs, FHIPμs, ks, Ms, As, Bs, Cs, Fs, Taus
#effectivemass=0.12 # the bare-electron band effective-mass.
# --> 0.12 for electrons and 0.15 for holes, in MAPI. See 2014 PRB.
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
MAPIe=polaronmobility(T, 4.5, 24.1, 2.25, 0.12)
MAPIh=polaronmobility(T, 4.5, 24.1, 2.25, 0.15)

# Hellwarth mobility
@test MAPIe.Hμ[1] ≈ 136.42 rtol=0.02
# Test variational parameters
@test MAPIe.v[1] ≈ 19.86 rtol=0.02
@test MAPIe.w[1] ≈ 16.96 rtol=0.02

# Same for the MAPI holes @ 300 K
@test MAPIh.Hμ[1] ≈ 94.15 rtol=0.02
@test MAPIh.v[1] ≈ 20.09  rtol=0.02
@test MAPIh.w[1] ≈ 16.81  rtol=0.02

end


text:
using Measurements
using DustExtinction: aa_to_invum,
                      g03_invum,
                      g16_invum,
                      G03_SMCBar,
                      G16

@testset "G16" begin

    function get_fA_1_ref()
        # From Fitzpatrick (1999) Table 3
        x_invum = [0.377, 0.820, 1.667, 1.828, 2.141, 2.433, 3.704, 3.846]
        ref_values = [0.265, 0.829, 2.688, 3.055, 3.806, 4.315, 6.265, 6.591] ./ 3.1
        return (x_invum, ref_values)
    end

    # Test Rv=3.1, fA=1.0 mixture
    x_invum, ref_values = get_fA_1_ref()
    model = G16(Rv=3.1, f_A=1.0)
    wave = 1e4 ./ x_invum
    @test model.(wave) ≈ ref_values rtol = 2e-3

    # Test RV=2.74, f_A=0 mixture
    tmodel = G16(Rv=2.74, f_A=0.0)
    gmodel = DustExtinction.G03_SMCBar()
    x_invum = gmodel.obsdata_x
    wave = 1e4 ./ x_invum
    ref_values = gmodel.obsdata_axav
    tolerance = gmodel.obsdata_tolerance
    @test collect(tmodel.(wave)) ≈ collect(ref_values) rtol = tolerance

    # Test out of bounds
    model = G16(Rv=3.1, f_A=1.0)
    bad_waves = [100, 4e4]
    @test model.(bad_waves) == zeros(length(bad_waves))
    @test_throws ErrorException g16_invum(aa_to_invum(bad_waves[1]), 3.1, 1.0)
    @test_throws ErrorException g16_invum(aa_to_invum(bad_waves[2]), 3.1, 1.0)

    # uncertainties
    x_invum, ref_values = get_fA_1_ref()
    wave = 1e4 ./ x_invum
    model = G16(Rv=3.1, f_A=1.0)
    noise = rand(length(wave)) .* 0.01
    wave_unc = wave .± noise
    reddening = map(w -> @uncertain(model(w)), wave_unc)
    @test Measurements.value.(reddening) ≈ ref_values rtol = 1e-3

    # Unitful
    model = G16(Rv=3.1, f_A=1.0)
    x_invum, ref_values = get_fA_1_ref()
    wave = 1e4 ./ x_invum
    wave_u = wave * u"angstrom"
    reddening = @inferred map(model, wave_u)
    @test eltype(reddening) <: Gain
    @test ustrip.(reddening) ≈ ref_values rtol = 0.016

end

@testset "G03_SMCBar" begin

    # Test output
    model = G03_SMCBar()
    x, y, tolerance = model.obsdata_x, model.obsdata_axav, model.obsdata_tolerance
    w = 1e4 ./ x
    @test model.(collect(w)) ≈ collect(y) rtol = tolerance

    # Test out of bounds
    model = G03_SMCBar(Rv=2.74)
    bad_waves = [100, 4e4]
    @test model.(bad_waves) == zeros(length(bad_waves))
    @test_throws ErrorException g03_invum(aa_to_invum(bad_waves[1]), 3.1)
    @test_throws ErrorException g03_invum(aa_to_invum(bad_waves[2]), 3.1)

    # uncertainties
    model = G03_SMCBar(Rv=2.74)
    x, y, tolerance = model.obsdata_x, model.obsdata_axav, model.obsdata_tolerance
    wave = 1e4 ./ x
    noise = rand(length(wave)) .* 0.01
    wave_unc = wave .± noise
    reddening = map(w -> @uncertain(model(w)), wave_unc)
    @test Measurements.value.(reddening) ≈ collect(y) rtol = tolerance

    # Unitful
    model = G03_SMCBar(Rv=2.74)
    x, y, tolerance = model.obsdata_x, model.obsdata_axav, model.obsdata_tolerance
    wave = 1e4 ./ x
    wave_u = collect(wave) * u"angstrom"
    reddening = @inferred map(model, wave_u)
    @test eltype(reddening) <: Gain
    @test ustrip.(reddening) ≈ collect(y) rtol = tolerance
end


text:
# this is the 2nd attempt; in this attempt I will try use MCTS and gradients
using StatsBase, DataFrames, Flux, Plots, JLD2, FileIO, BSON
include("2048_fn.jl")

struct IdentitySkip
   inner
   activation
end

(m::IdentitySkip)(x) = m.activation.(m.inner(x) .+ x)
IdentitySkip(Dense(128, 128), relu)

policy = Flux.Chain(
  grid_to_feature,
  Dense(16*14, 128, relu),
  IdentitySkip(Dense(128, 128), relu),
  Dense(128, 32, relu),
  IdentitySkip(Dense(32, 32), relu),
  Dense(32, 4),
  softmax
  )

loss(x, (move, reward)) = reward*log(policy(init)[move .== DIRS][1] + eps(Float64))
opt=ADAM(Flux.params(policy))

function simulate_game_on_policy(policy)
  agrid = init_game()
  init = copy(agrid)
  # play according to policy
  vm = valid_moves(agrid)
  rewards = Int[]
  seq = Symbol[]
  carts = CartesianIndex{2}[]
  tfs = Int[]

  while length(vm) > 0
    sample_weight = softmax(Tracker.data(policy(agrid))[indexin(vm, DIRS)])
    sampled_move = sample(vm, Weights(sample_weight))
    (agrid, cart, tf, reward) = simulate_one!(agrid, sampled_move)
    push!(seq, sampled_move)
    push!(rewards, reward)
    push!(carts, cart)
    push!(tfs, tf)
    vm = valid_moves(agrid)
  end
  (init, agrid, seq, carts, tfs, rewards)
end

@time (init, fnl, seq, cart, tf, reward) = simulate_game_on_policy(policy)

function convert_to_training_sample(init, seq, carts, tfs, rewards)
  reward_disc_rate = 0.99
  agrid = copy(init)
  discr = discounted_rewards(rewards, reward_disc_rate)
  #discr = (discr .- StatsBase.mean(discr))./StatsBase.std(discr)

  res_vec = Tuple{Array{Int,2}, Tuple{Symbol, Float64}}[]
  for (move, cart, tf, reward) in zip(seq, carts, tfs, discr)
    (rotated_grids, rotated_moves) = rotations(copy(agrid), move)
    mean_random_reward = avg_disc_reward(agrid, 20, reward_disc_rate)
    repeat_reward = repeat([reward-mean_random_reward], inner=8)
    append!(res_vec, zip(rotated_grids, zip(rotated_moves, repeat_reward)) |> collect)
    replay_move!(agrid, move, cart, tf)
  end
  res_vec
end

@time res = convert_to_training_sample(init, seq, cart, tf, reward)

function train_sample_batch(policy, batch_size)
  res_vec = Tuple{Array{Int,2}, Tuple{Symbol, Float64}}[]
  while length(res_vec) < batch_size
    (init, _, seq, cart, tf, reward) = simulate_game_on_policy(policy)
    res = convert_to_training_sample(init, seq, cart, tf, reward)
    append!(res_vec, res)
  end
  res_vec
end

max_prob = 0.25

while true
  global max_prob
  res_vec = train_sample_batch(policy, 20480)

  # pb2 = findmax([policy(r[1]) |> Tracker.data |> maximum for r in res_vec])
  # best_one = res_vec[pb2[2]][1]
  # vm = valid_moves(best_one)
  # best_move = DIRS[findmax(policy(best_one) |> Tracker.data)[2]]
  # if any(best_move .== vm) #& (pb2[1] > max_prob)
  #   pb = findmax(policy(best_one) |> Tracker.data)
  #   display("*********************************")
  #   display(best_one)
  #   display((pb[1], DIRS[pb[2]]))
  #   display(Dict(zip(DIRS, policy(best_one) |> Tracker.data)))
  #   display(avg_rewards_next_move(best_one, 100))
  #   max_prob = pb[1]
  # end
  @time Flux.train!(loss, res_vec, opt)
  @save "policy.jld2" policy
  # if any(best_move .== vm) #& (pb2[1] > max_prob)
  #   display(Dict(zip(DIRS, policy(best_one) |> Tracker.data)))
  #   #display((pb[1], DIRS[pb[2]]))
  #   max_prob = pb[1]
  # end
  display(simulate_game_on_policy(policy)[2])
end

function abc()
  init = init_game()
  (init, fnl, seq, cart, tf, reward) = simulate_game_on_policy(policy)
  maximum(fnl)
end

policy_max = [abc() for i=1:10000]
random_max = [maximum(simulate_game!(init_game())[2]) for i = 1:10000]

countmap(policy_max)
countmap(random_max)

#nback = Int(round(length(seq)/2, digits=0))
for i = 1:length(seq)-1
  bgrid = sim_seq!(copy(init), seq[1:end-i], cart[1:end-i], tf[1:end-i])
  prob = round.(policy(bgrid) |> Tracker.data, digits=4)
  if all(prob .!= 0.25)
    display(bgrid)
    display(DataFrame(dir = DIRS, prob = prob))
  end
end

Flux.train!(loss, res_vec, opt)
softmax(policy(bgrid))

loss(res_vec[1][1], res_vec[2][2])

loss(init, (:left, 1))

loss2(grid,())reward*Flux.crossentropy(softmax(policy(move)), move .== DIRS)

loss1(grid,(move, reward)) = (grid, move, reward)
loss1(res_vec[1]...)

display(agrid)
display(transpose(agrid))

:left => :up
:right => :down
:up => :left
:down => :right

(agrid, cart, tf, reward) = simulate_one!(agrid, :left)

typeof(cart)

function train1!(policy, critic, policy_loss, critic_loss, opt, opt_critic, n;
  γ = 0.99 # next critic moves downward rate
  )
  for i = 1:n
    agrid = init_game()
    vm = valid_moves(agrid)
    while length(vm) > 0
      # use policy to generate searches
      p = Tracker.data(policy(agrid))[indexin(vm, DIRS)]

      sampled_move = sample(vm, Weights(p))
      (new_grid, _, _, reward) = simulate_one!(copy(agrid), sampled_move)

      # sample a new move for the new_grid
      new_vm = valid_moves(new_grid)

      if length(new_vm) > 0
        p = Tracker.data(policy(new_grid))[indexin(new_vm, DIRS)]
        sampled_next_action = sample(new_vm, Weights(p))

        # evaluate critic at the new location
        discounted_critic_value = γ*Tracker.data(critic((new_grid, sampled_next_action)))[1]
      else
        discounted_critic_value = 0
      end

      # update the critic
      Flux.train!(critic_loss, [((agrid, sampled_move), reward + discounted_critic_value)], opt_critic)

      # now use the updated critic to update to the policy
      best_move_according_to_new_critic = findmax(Dict{Symbol, Float64}((vm1, Tracker.data(critic((agrid, vm1)))[1]) for vm1 in vm))[2]
      Flux.train!(policy_loss, [(agrid, Flux.onehot(best_move_according_to_new_critic, DIRS))], opt)

      agrid = new_grid
      vm = new_vm
    end
    i += 1
    if i % 1000 == 0
      @save "actor_critic_model.jld2" policy critic
      display(agrid)
    end
  end
end

init = init_game()

policy(init)
[critic((init,vm)) for vm in DIRS]
@time train1!(policy, critic, policy_loss, critic_loss, opt, opt_critic, 1_000_000)
policy(init)
critic((init,:left))

(init, fnl, seq, cart, tf, rewards) = simulate_game!(init)

nmovesback = 1
last_grid = sim_seq(init, seq[1:end-nmovesback], cart[1:end-nmovesback], tf[1:end-nmovesback])
display(last_grid)
display(policy(last_grid))
display([(vm, Tracker.data(critic((last_grid, vm)))[1]) for vm in valid_moves(last_grid)])

#manually train this
agrid = copy(last_grid)
vm = valid_moves(agrid)

# use policy to generate searches
p = Tracker.data(policy(agrid))[indexin(vm, DIRS)]

sampled_move = sample(vm, Weights(p))
(new_grid, _, _, reward) = simulate_one!(copy(agrid), sampled_move)

# sample a new move for the new_grid
new_vm = valid_moves(new_grid)

if length(new_vm) > 0
  p = Tracker.data(policy(new_grid))[indexin(new_vm, DIRS)]
  sampled_next_action = sample(new_vm, Weights(p))

  # evaluate critic at the new location
  discounted_critic_value = γ*Tracker.data(critic((new_grid, sampled_next_action)))[1]
else
  discounted_critic_value = 0
end

critic((agrid, :right))
# update the critic
Flux.train!(critic_loss, [((agrid, sampled_move), reward + discounted_critic_value)], opt_critic)
display((critic((agrid, :right))|> Tracker.data, critic((agrid, :left))|>Tracker.data))

# now use the updated critic to update to the policy
best_move_according_to_new_critic = findmax(Dict{Symbol, Float64}((vm1, Tracker.data(critic((agrid, vm1)))[1]) for vm1 in vm))[2]
Flux.train!(policy_loss, [(agrid, Flux.onehot(best_move_according_to_new_critic, DIRS))], opt)

agrid = new_grid
old_grid = copy(new_grid)
display(agrid)
vm = new_vm

critic((init,:left))

agrid = last_grid

vm = valid_moves(agrid)

(init, fnl, seq, cart, tf, reward) = simulate_game!(agrid)

# obtain the last move
last_move = 2
before_last = sim_seq(init, seq[1:end-last_move], cart[1:end-last_move], tf[1:end-last_move])

vm = valid_moves(before_last)
@time avg_rewards_next_move(before_last, 100)

function train_n_games(n, model, opt)
  for i=1:n
    agrid = init_game()

    vm = valid_moves(agrid)
    while length(vm) > 0
      # compute average rewards
      avgr = avg_rewards_next_move(agrid, 100)
      grid_as_features = grid_to_feature(agrid)
      #model(grid_as_features)
      Flux.train!(loss, [(grid_as_features, findmax(avgr)[2])], opt)

      # now choose an action based on the updated policy
      themove = sample(vm, Weights(Tracker.data(model(grid_as_features))[indexin([:left,:up], DIRS)]))
      simulate_one!(agrid, themove)
      vm = valid_moves(agrid)
    end
    println(i)
    display(agrid)

    @save "policy_gradient_model.jld2" model

  end
end

@time train_n_games(1, model, opt)

agrid = init_game()
(init,_,seq,cart, tf,_) = simulate_game!(agrid)

new_grid = sim_seq(init, seq[1:end-2], cart[1:end-2], tf[1:end-2])

model(grid_to_feature(new_grid))
avgr = avg_rewards_next_move(new_grid, 100)
best_move = findmax(avgr)[2]
Flux.@epochs 10000 Flux.train!(loss, [(grid_to_feature(new_grid), Flux.onehot(best_move, DIRS))], opt)
model(grid_to_feature(new_grid))

Float64.(Flux.onehot(best_move, DIRS))

crossentropy(model(grid_to_feature(new_grid)), Float64.(Flux.onehot(best_move, DIRS)))

logloss(model(grid_to_feature(new_grid)), Float64.(Flux.onehot(:up, DIRS)))

grid_to_feature(before_last) |>
  model |>
  Tracker.data

grid_to_feature(before_last) |>
  model |>
  Tracker.data

function assess_model(model)
  init_grid = init_game()
  @time (init_grid, fnl, seq, cart, tf, reward) = simulate_game!(init_grid, model, 0.0)

  # obtain the discounted rewards
  y = discounted_rewards(reward, 0.99)[1:end-1];

  # score every move using  model
  updated_model = [Tracker.data(model(grid_move_to_feature(sim_seq(init_grid, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]
  display(sim_seq(init_grid, seq, cart, tf))
  g = plot(updated_model, label="model")
  plot!(y, label="discounted rewards")
  (init_grid, fnl, seq, cart, tf, reward, g)
end

ϵ = 0.5

model = nothing
@load "no_neg_model.jld2"
model

@time for i=1:10000
  global ϵ
  ϵ -= 0.0001
  simulate_one_game_and_train!(model, opt, ϵ)
  if i % 100 == 0
    @save "no_neg_model.jld2" model;
    println("epsilon: $ϵ")
    init_grid = init_game()
    res = simulate_game!(init_grid, model, 0)
    display(res[2])
    @save "no_neg_model.jld2" model;
  end
end

@time (init, fnl, seq, cart, tf, reward, g) = assess_model(model);
g

last_n_move = length(seq)-20
grid_new = sim_seq(init, seq[1:end-last_n_move], cart[1:end-last_n_move], tf[1:end - last_n_move])
grid_try = copy(grid_new)
display(grid_new)
vm = valid_moves(grid_new)
x = vec(Tracker.data(model(hcat(grid_move_to_feature.([grid_new], vm)...))))
y = x .- minimum(x)

using DataFrames
DataFrame(vm=vm, chosen=vm.== seq[end-last_n_move], x=vec(x),  prob_exp = round.(exp.(y)/sum(exp.(y)),digits=2), prob_square =((y.+1).^2)./sum((y.+1).^2))


text:
# Getting length
print("Emter length: ")
length = readline()
length = parse(Float64, length) # Converting string to Float

# Getting width
print("Emter width: ")
width = readline()
width = parse(Float64, width) # Converting string to Float

# Calculating Area
area = length*width

println("Area of Rectangle is ", area)


text:
using HDF5
using CRC32c
using Test

@testset "plain" begin

# Create a new file
fn = tempname()
f = h5open(fn, "w")
@test isopen(f)
# Write scalars
f["Float64"] = 3.2
f["Int16"] = Int16(4)
# compression of empty array (issue #246)
f["compressedempty", "shuffle", (), "compress", 4] = Int64[]
# compression of zero-dimensional array (pull request #445)
f["compressed_zerodim", "shuffle", (), "compress", 4] = fill(Int32(42), ())
f["bloscempty", "blosc", 4] = Int64[]
# Create arrays of different types
A = randn(3, 5)
write(f, "Afloat64", convert(Matrix{Float64}, A))
write(f, "Afloat32", convert(Matrix{Float32}, A))
Ai = rand(1:20, 2, 4)
write(f, "Aint8", convert(Matrix{Int8}, Ai))
f["Aint16"] = convert(Matrix{Int16}, Ai)
write(f, "Aint32", convert(Matrix{Int32}, Ai))
write(f, "Aint64", convert(Matrix{Int64}, Ai))
write(f, "Auint8", convert(Matrix{UInt8}, Ai))
write(f, "Auint16", convert(Matrix{UInt16}, Ai))
write(f, "Auint32", convert(Matrix{UInt32}, Ai))
write(f, "Auint64", convert(Matrix{UInt64}, Ai))
# Arrays of bools (pull request #540)
Abool = [false, true, false]
write(f, "Abool", Abool)

salut = "Hi there"
ucode = "uniçº∂e"
write(f, "salut", salut)
write(f, "ucode", ucode)
# Manually write a variable-length string (issue #187)
let
    dtype = HDF5Datatype(HDF5.h5t_copy(HDF5.H5T_C_S1))
    HDF5.h5t_set_size(dtype.id, HDF5.H5T_VARIABLE)
    HDF5.h5t_set_cset(dtype.id, HDF5.cset(typeof(salut)))
    dspace = HDF5.dataspace(salut)
    dset = HDF5.d_create(f, "salut-vlen", dtype, dspace)
    HDF5.h5d_write(dset, dtype, HDF5.H5S_ALL, HDF5.H5S_ALL, HDF5.H5P_DEFAULT, [pointer(salut)])
end
# Arrays of strings
salut_split = ["Hi", "there"]
write(f, "salut_split", salut_split)
salut_2d = ["Hi" "there"; "Salut" "friend"]
write(f, "salut_2d", salut_2d)
# Arrays of strings as vlen
vlen = HDF5Vlen(salut_split)
d_write(f, "salut_vlen", vlen)
# Empty arrays
empty = UInt32[]
write(f, "empty", empty)
# Empty strings
empty_string = ""
write(f, "empty_string", empty_string)
# Empty array of strings
empty_string_array = String[]
write(f, "empty_string_array", empty_string_array)
# Array of empty string
empty_array_of_strings = [""]
write(f, "empty_array_of_strings", empty_array_of_strings)
# Attributes
species = [["N", "C"]; ["A", "B"]]
attrs(f)["species"] = species
C∞ = 42
attrs(f)["C∞"] = C∞
dset = f["salut"]
@test !isempty(dset)
label = "This is a string"
attrs(dset)["typeinfo"] = label
close(dset)
# Scalar reference values in attributes
attrs(f)["ref_test"] = HDF5.HDF5ReferenceObj(f, "empty_array_of_strings")
@test read(attrs(f)["ref_test"]) === HDF5.HDF5ReferenceObj(f, "empty_array_of_strings")
# Group
g = g_create(f, "mygroup")
# Test dataset with compression
R = rand(1:20, 20, 40);
g["CompressedA", "chunk", (5, 6), "shuffle", (), "compress", 9] = R
g["BloscA", "chunk", (5, 6), "shuffle", (), "blosc", 9] = R
close(g)
# Copy group containing dataset
o_copy(f, "mygroup", f, "mygroup2")
# Copy dataset
g = g_create(f, "mygroup3")
o_copy(f["mygroup/CompressedA"], g, "CompressedA")
o_copy(f["mygroup/BloscA"], g, "BloscA")
close(g)
# Writing hyperslabs
dset = d_create(f, "slab", datatype(Float64), dataspace(20, 20, 5), "chunk", (5, 5, 1))
Xslab = randn(20, 20, 5)
for i = 1:5
    dset[:,:,i] = Xslab[:,:,i]
end
# More complex hyperslab and assignment with "incorrect" types (issue #34)
d = d_create(f, "slab2", datatype(Float64), ((10, 20), (100, 200)), "chunk", (1, 1))
d[:,:] = 5
d[1,1] = 4
# 1d indexing
d = d_create(f, "slab3", datatype(Int), ((10,), (-1,)), "chunk", (5,))
@test d[:] == zeros(Int, 10)
d[3:5] = 3:5
# Create a dataset designed to be deleted
f["deleteme"] = 17.2
close(f)
@test !isopen(f)
# Test the h5read/write interface, with attributes
W = copy(reshape(1:120, 15, 8))
Wa = Dict("a" => 1, "b" => 2)
h5write(fn, "newgroup/W", W)
h5writeattr(fn, "newgroup/W", Wa)

# Read the file back in
fr = h5open(fn)
x = read(fr, "Float64")
@test x == 3.2 && isa(x, Float64)
y = read(fr, "Int16")
@test y == 4 && isa(y, Int16)
zerodim = read(fr, "compressed_zerodim")
@test zerodim == 42 && isa(zerodim, Int32)
bloscempty = read(fr, "bloscempty")
@test bloscempty == Int64[] && isa(bloscempty, Vector{Int64})
Af32 = read(fr, "Afloat32")
@test convert(Matrix{Float32}, A) == Af32
@test eltype(Af32) == Float32
Af64 = read(fr, "Afloat64")
@test convert(Matrix{Float64}, A) == Af64
@test eltype(Af64) == Float64
@test eltype(fr["Afloat64"]) == Float64  # issue 167
Ai8 = read(fr, "Aint8")
@test Ai == Ai8
@test eltype(Ai8) == Int8
Ai16 = read(fr, "Aint16")
@test Ai == Ai16
@test eltype(Ai16) == Int16
Ai32 = read(fr, "Aint32")
@test Ai == Ai32
@test eltype(Ai32) == Int32
Ai64 = read(fr, "Aint64")
@test Ai == Ai64
@test eltype(Ai64) == Int64
Ai8 = read(fr, "Auint8")
@test Ai == Ai8
@test eltype(Ai8) == UInt8
Ai16 = read(fr, "Auint16")
@test Ai == Ai16
@test eltype(Ai16) == UInt16
Ai32 = read(fr, "Auint32")
@test Ai == Ai32
@test eltype(Ai32) == UInt32
Ai64 = read(fr, "Auint64")
@test Ai == Ai64
@test eltype(Ai64) == UInt64

Abool_read = read(fr, "Abool")
@test Abool_read == Abool
@test eltype(Abool_read) == Bool

salutr = read(fr, "salut")
@test salut == salutr
salutr = read(fr, "salut-vlen")
@test salut == salutr
ucoder = read(fr, "ucode")
@test ucode == ucoder
salut_splitr = read(fr, "salut_split")
@test salut_splitr == salut_split
salut_2dr = read(fr, "salut_2d")
@test salut_2d == salut_2dr
salut_vlenr = read(fr, "salut_vlen")
@test salut_vlenr == salut_split
Rr = read(fr, "mygroup/CompressedA")
@test Rr == R
Rr2 = read(fr, "mygroup2/CompressedA")
@test Rr2 == R
Rr3 = read(fr, "mygroup3/CompressedA")
@test Rr3 == R
Rr4 = read(fr, "mygroup/BloscA")
@test Rr4 == R
Rr5 = read(fr, "mygroup2/BloscA")
@test Rr5 == R
Rr6 = read(fr, "mygroup3/BloscA")
@test Rr6 == R
dset = fr["mygroup/CompressedA"]
@test get_chunk(dset) == (5, 6)
@test name(dset) == "/mygroup/CompressedA"
dset2 = fr["mygroup/BloscA"]
@test get_chunk(dset2) == (5, 6)
@test name(dset2) == "/mygroup/BloscA"
Xslabr = read(fr, "slab")
@test Xslabr == Xslab
Xslabr = h5read(fn, "slab", (:, :, :))  # issue #87
@test Xslabr == Xslab
Xslab2r = read(fr, "slab2")
target = fill(5, 10, 20)
target[1] = 4
@test Xslab2r == target
dset = fr["slab3"]
@test dset[3:5] == [3:5;]
emptyr = read(fr, "empty")
@test isempty(emptyr)
empty_stringr = read(fr, "empty_string")
@test empty_stringr == empty_string
empty_string_arrayr = read(fr, "empty_string_array")
@test empty_string_arrayr == empty_string_array
empty_array_of_stringsr = read(fr, "empty_array_of_strings")
@test empty_array_of_stringsr == empty_array_of_strings
@test a_read(fr, "species") == species
@test a_read(fr, "C∞") == C∞
dset = fr["salut"]
@test a_read(dset, "typeinfo") == label
close(dset)
# Test ref-based reading
Aref = fr["Afloat64"]
sel = (2:3, 1:2:5)
Asub = Aref[sel...]
@test Asub == A[sel...]
close(Aref)
# Test iteration, name, and parent
for obj in fr
    @test filename(obj) == fn
    n = name(obj)
    p = parent(obj)
end
# Test reading multiple vars at once
z = read(fr, "Float64", "Int16")
@test z == (3.2, 4)
@test typeof(z) == Tuple{Float64,Int16}
# Test function syntax
read(fr, "Float64") do x
    @test x == 3.2
end
read(fr, "Float64", "Int16") do x, y
    @test x == 3.2
    @test y == 4
end
# Test reading entire file at once
z = read(fr)
@test z["Float64"] == 3.2
close(fr)

# Test object deletion
fr = h5open(fn, "r+")
@test exists(fr, "deleteme")
o_delete(fr, "deleteme")
@test !exists(fr, "deleteme")
close(fr)

# Test the h5read interface
Wr = h5read(fn, "newgroup/W")
@test Wr == W
rng = (2:3:15, 3:5)
Wr = h5read(fn, "newgroup/W", rng)
@test Wr == W[rng...]
War = h5readattr(fn, "newgroup/W")
@test War == Wa

# more do syntax
h5open(fn, "w") do fid
    g_create(fid, "mygroup") do g
        write(g, "x", 3.2)
    end
end
fid = h5open(fn, "r")
@test names(fid) == ["mygroup"]
g = fid["mygroup"]
@test names(g) == ["x"]
close(g)
close(fid)
rm(fn)

# more do syntax: atomic rename version
tmpdir = mktempdir()
outfile = joinpath(tmpdir, "test.h5")

# create a new file
h5rewrite(outfile) do fid
    g_create(fid, "mygroup") do g
        write(g, "x", 3.3)
    end
end
@test length(readdir(tmpdir)) == 1
h5open(outfile, "r") do fid
    @test names(fid) == ["mygroup"]
    @test names(fid["mygroup"]) == ["x"]
end

# fail to overwrite
@test_throws ErrorException h5rewrite(outfile) do fid
    g_create(fid, "mygroup") do g
        write(g, "oops", 3.3)
    end
    error("failed")
end
@test length(readdir(tmpdir)) == 1
h5open(outfile, "r") do fid
    @test names(fid) == ["mygroup"]
    @test names(fid["mygroup"]) == ["x"]
end

# overwrite
h5rewrite(outfile) do fid
    g_create(fid, "mygroup") do g
        write(g, "y", 3.3)
    end
end
@test length(readdir(tmpdir)) == 1
h5open(outfile, "r") do fid
    @test names(fid) == ["mygroup"]
    @test names(fid["mygroup"]) == ["y"]
end
rm(tmpdir, recursive=true)

test_files = joinpath(@__DIR__, "test_files")

d = h5read(joinpath(test_files, "compound.h5"), "/data")
@test typeof(d[1]) === HDF5.HDF5Compound{4}
@test length(d) == 2
dtypes = [typeof(x) for x in d[1].data]
@test dtypes == [Float64, Vector{Float64}, Vector{Float64}, Float64]
@test length(d[1].data[2]) == 3
@test d[1].membername == ("wgt", "xyz", "uvw", "E")

# get-datasets
fn = tempname()
fd = h5open(fn, "w")
fd["level_0"] = [1,2,3]
grp = g_create(fd, "mygroup")
fd["mygroup/level_1"] = [4, 5]
grp2 = g_create(grp, "deep_group")
fd["mygroup/deep_group/level_2"] = [6.0, 7.0]
datasets = get_datasets(fd)
@test sort(map(name, datasets)) ==  sort(["/level_0", "/mygroup/deep_group/level_2", "/mygroup/level_1"])
close(fd)
rm(fn)

# File creation and access property lists
cpl = p_create(HDF5.H5P_FILE_CREATE)
cpl["userblock"] = 1024
apl = p_create(HDF5.H5P_FILE_ACCESS)
apl["libver_bounds"] = (HDF5.H5F_LIBVER_EARLIEST, HDF5.H5F_LIBVER_LATEST)
h5open(fn, false, true, true, true, false, cpl, apl) do fid
    write(fid, "intarray", [1, 2, 3])
end
h5open(fn, "r", "libver_bounds",
    (HDF5.H5F_LIBVER_EARLIEST, HDF5.H5F_LIBVER_LATEST)) do fid
    intarray = read(fid, "intarray")
    @test intarray == [1, 2, 3]
end

# Test null terminated ASCII string (e.g. exported by h5py) #332
h5open(joinpath(test_files, "nullterm_ascii.h5"), "r") do fid
    str = read(fid["test"])
    @test str == "Hello World"
end

@test HDF5.unpad(UInt8[0x43, 0x43, 0x41], 1) == "CCA"

# Test the h5read/write interface with a filename as a first argument, when
# the file does not exist
rm(fn)
h5write(fn, "newgroup/W", W)
Wr = h5read(fn, "newgroup/W")
@test Wr == W
close(f)
rm(fn)

if !isempty(HDF5.libhdf5_hl)
    # Test direct chunk writing
    h5open(fn, "w") do f
      d = d_create(f, "dataset", datatype(Int), dataspace(4, 4), "chunk", (2, 2))
      raw = HDF5ChunkStorage(d)
      raw[1,1] = 0, collect(reinterpret(UInt8, [1,2,5,6]))
      raw[3,1] = 0, collect(reinterpret(UInt8, [3,4,7,8]))
      raw[1,3] = 0, collect(reinterpret(UInt8, [9,10,13,14]))
      raw[3,3] = 0, collect(reinterpret(UInt8, [11,12,15,16]))
    end

    @test h5open(fn, "r") do f
      vec(f["dataset"][:,:])
    end == collect(1:16)

    close(f)
    rm(fn)
end

# Test that switching time tracking off results in identical files
h5open("tt1.h5", "w") do f
    f["x", "track_times", false] = [1, 2, 3]
end
h5open("tt2.h5", "w") do f
    f["x", "track_times", false] = [1, 2, 3]
end

@test open(crc32c, "tt1.h5") == open(crc32c, "tt2.h5")

end # testset plain

@testset "complex" begin
  HDF5.enable_complex_support()

  fn = tempname()
  f = h5open(fn, "w")

  f["ComplexF64"] = 1.0 + 2.0im
  attrs(f["ComplexF64"])["ComplexInt64"] = 1im

  Acmplx = rand(ComplexF64, 3, 5)
  write(f, "Acmplx64", convert(Matrix{ComplexF64}, Acmplx))
  write(f, "Acmplx32", convert(Matrix{ComplexF32}, Acmplx))

  HDF5.disable_complex_support()
  @test_throws ErrorException f["_ComplexF64"] = 1.0 + 2.0im
  @test_throws ErrorException write(f, "_Acmplx64", convert(Matrix{ComplexF64}, Acmplx))
  @test_throws ErrorException write(f, "_Acmplx32", convert(Matrix{ComplexF32}, Acmplx))
  HDF5.enable_complex_support()

  close(f)

  fr = h5open(fn)
  z = read(fr, "ComplexF64")
  @test z == 1.0 + 2.0im && isa(z, ComplexF64)
  z_attrs = attrs(fr["ComplexF64"])
  @test read(z_attrs["ComplexInt64"]) == 1im

  Acmplx32 = read(fr, "Acmplx32")
  @test convert(Matrix{ComplexF32}, Acmplx) == Acmplx32
  @test eltype(Acmplx32) == ComplexF32
  Acmplx64 = read(fr, "Acmplx64")
  @test convert(Matrix{ComplexF64}, Acmplx) == Acmplx64
  @test eltype(Acmplx64) == ComplexF64

  HDF5.disable_complex_support()
  z = read(fr, "ComplexF64")
  @test isa(z, HDF5.HDF5Compound{2})

  Acmplx32 = read(fr, "Acmplx32")
  @test eltype(Acmplx32) == HDF5.HDF5Compound{2}
  Acmplx64 = read(fr, "Acmplx64")
  @test eltype(Acmplx64) == HDF5.HDF5Compound{2}

  close(fr)

  HDF5.enable_complex_support()
end

# test strings with null and undefined references
@testset "undefined and null" begin
fn = tempname()
f = h5open(fn, "w")

# don't silently truncate data
@test_throws ArgumentError write(f, "test", ["hello","there","\0"])
@test_throws ArgumentError write(f, "trunc1", "\0")
@test_throws ArgumentError write(f, "trunc2", "trunc\0ateme")

# test writing uninitialized string arrays
undefstrarr = similar(Vector(1:3), String) # strs = String[#undef, #undef, #undef]
@test_throws UndefRefError write(f, "undef", undefstrarr)

close(f)
rm(fn)

end # testset null and undefined

# test writing abstract arrays
@testset "abstract arrays" begin

# test writing reinterpreted data
fn = tempname()
try
    h5open(fn, "w") do f
        data = reinterpret(UInt8, [true, false, false])
        write(f, "reinterpret array", data)
    end

    @test h5open(fn, "r") do f
        read(f, "reinterpret array")
    end == UInt8[0x01, 0x00, 0x00]
finally
    rm(fn)
end

# don't silently fail for arrays with a different stride
fn = tempname()
try
    data = rand(UInt16, 2, 3);
    pdv_data = PermutedDimsArray(data, (2, 1))

    @test_throws ArgumentError h5write(fn, "pdv_data", pdv_data)
finally
    rm(fn)
end

end # writing abstract arrays


text:
using DataDeps

register(DataDep(
  "magic-gamma-telescope",
  "http://archive.ics.uci.edu/ml/datasets/MAGIC+Gamma+Telescope",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/magic/magic04.data",
  "e9314b7ebd4b4b59a3b3d65f7316663963777b16a46786877651dbbaa640b36a",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=11,
        feature_indices=1:10,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Test
using myJuliaUtils
using Statistics

T = 100;
K = rand(1:20);
x = abs.(randn(T, K))
@testset "testing dataTransf" begin

    a = dataTransf(ones(T), 5)
    @test a.y == zeros(T)
    @test a.n_loss == 0

    @test dataTransf(x, 1) == dataTransf(x, ones(Int, K))

    transfId = rand(1:8)
    (a1, b1) = dataTransf(x, transfId)
    (a2, b2) = dataTransf(x, transfId * ones(Int, K))

    @test b1 == b2
    @test a1[b1 + 1:end,:] == a2[b2 + 1:end,:]
end

T = 100;
K = rand(1:20);
x = abs.(randn(T, K))

@testset "testing detrend functions" begin
    @test detrend(x, power = 0) ≈ x .- mean(x, dims = 1)
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
"""
    refco(phpa, tc, rh, wl)

Determine the constants A and B in the atmospheric refraction model
dZ = A tan Z + B tan^3 Z.

Z is the "observed" zenith distance (i.e. affected by refraction)
and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo)
zenith distance.

### Given ###

- `phpa`: Pressure at the observer (hPa = millibar)
- `tc`: Ambient temperature at the observer (deg C)
- `rh`: Relative humidity at the observer (range 0-1)
- `wl`: Wavelength (micrometers)

### Returned ###

- `refa`: tan Z coefficient (radians)
- `refb`: tan^3 Z coefficient (radians)

### Notes ###

1. The model balances speed and accuracy to give good results in
   applications where performance at low altitudes is not paramount.
   Performance is maintained across a range of conditions, and
   applies to both optical/IR and radio.

2. The model omits the effects of (i) height above sea level (apart
   from the reduced pressure itself), (ii) latitude (i.e. the
   flattening of the Earth), (iii) variations in tropospheric lapse
   rate and (iv) dispersive effects in the radio.

   The model was tested using the following range of conditions:

     lapse rates 0.0055, 0.0065, 0.0075 deg/meter
     latitudes 0, 25, 50, 75 degrees
     heights 0, 2500, 5000 meters ASL
     pressures mean for height -10% to +5% in steps of 5%
     temperatures -10 deg to +20 deg with respect to 280 deg at SL
     relative humidity 0, 0.5, 1
     wavelengths 0.4, 0.6, ... 2 micron, + radio
     zenith distances 15, 45, 75 degrees

   The accuracy with respect to raytracing through a model
   atmosphere was as follows:

                          worst         RMS

     optical/IR           62 mas       8 mas
     radio               319 mas      49 mas

   For this particular set of conditions:

     lapse rate 0.0065 K/meter
     latitude 50 degrees
     sea level
     pressure 1005 mb
     temperature 280.15 K
     humidity 80%
     wavelength 5740 Angstroms

   the results were as follows:

     ZD       raytrace     eraRefco   Saastamoinen

     10         10.27        10.27        10.27
     20         21.19        21.20        21.19
     30         33.61        33.61        33.60
     40         48.82        48.83        48.81
     45         58.16        58.18        58.16
     50         69.28        69.30        69.27
     55         82.97        82.99        82.95
     60        100.51       100.54       100.50
     65        124.23       124.26       124.20
     70        158.63       158.68       158.61
     72        177.32       177.37       177.31
     74        200.35       200.38       200.32
     76        229.45       229.43       229.42
     78        267.44       267.29       267.41
     80        319.13       318.55       319.10

    deg        arcsec       arcsec       arcsec

   The values for Saastamoinen's formula (which includes terms
   up to tan^5) are taken from Hohenkerk and Sinclair (1985).

3. A wl value in the range 0-100 selects the optical/IR case and is
   wavelength in micrometers.  Any value outside this range selects
   the radio case.

4. Outlandish input parameters are silently limited to
   mathematically safe values.  Zero pressure is permissible, and
   causes zeroes to be returned.

5. The algorithm draws on several sources, as follows:

   a) The formula for the saturation vapour pressure of water as
      a function of temperature and temperature is taken from
      Equations (A4.5-A4.7) of Gill (1982).

   b) The formula for the water vapour pressure, given the
      saturation pressure and the relative humidity, is from
      Crane (1976), Equation (2.5.5).

   c) The refractivity of air is a function of temperature,
      total pressure, water-vapour pressure and, in the case
      of optical/IR, wavelength.  The formulae for the two cases are
      developed from Hohenkerk & Sinclair (1985) and Rueger (2002).

   d) The formula for beta, the ratio of the scale height of the
      atmosphere to the geocentric distance of the observer, is
      an adaption of Equation (9) from Stone (1996).  The
      adaptations, arrived at empirically, consist of (i) a small
      adjustment to the coefficient and (ii) a humidity term for the
      radio case only.

   e) The formulae for the refraction constants as a function of
      n-1 and beta are from Green (1987), Equation (4.31).

### References ###

- Crane, R.K., Meeks, M.L. (ed), "Refraction Effects in the Neutral
    Atmosphere", Methods of Experimental Physics: Astrophysics 12B,
    Academic Press, 1976.

- Gill, Adrian E., "Atmosphere-Ocean Dynamics", Academic Press,
    1982.

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
    1987.

- Hohenkerk, C.Y., & Sinclair, A.T., NAO Technical Note No. 63,
    1985.

- Rueger, J.M., "Refractive Index Formulae for Electronic Distance
    Measurement with Radio and Millimetre Waves", in Unisurv Report
    S-68, School of Surveying and Spatial Information Systems,
    University of New South Wales, Sydney, Australia, 2002.

- Stone, Ronald C., P.A.S.P. 108, 1051-1058, 1996.

"""
function refco(phpa, tk, rh, wl)
    refa = Ref(0.0)
    refb = Ref(0.0)
    ccall((:eraRefco, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          phpa, tk, rh, wl, refa, refb)
    refa[], refb[]
end

"""
    rm2v(r)

Express an r-matrix as an r-vector.

### Given ###

- `r`: Rotation matrix

### Returned ###

- `w`: Rotation vector (Note 1)

### Notes ###

1. A rotation matrix describes a rotation through some angle about
   some arbitrary axis called the Euler axis.  The "rotation vector"
   returned by this function has the same direction as the Euler axis,
   and its magnitude is the angle in radians.  (The magnitude and
   direction can be separated by means of the function eraPn.)

2. If r is null, so is the result.  If r is not a rotation matrix
   the result is undefined;  r must be proper (i.e. have a positive
   determinant) and real orthogonal (inverse = transpose).

3. The reference frame rotates clockwise as seen looking along
   the rotation vector from the origin.

"""
function rm2v(r)
    w = zeros(3)
    ccall((:eraRm2v, liberfa), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}),
          r, w)
    w
end

"""
    rv2m(w)

Form the r-matrix corresponding to a given r-vector.

### Given ###

- `w`: Rotation vector (Note 1)

### Returned ###

- `r`: Rotation matrix

### Notes ###

1. A rotation matrix describes a rotation through some angle about
   some arbitrary axis called the Euler axis.  The "rotation vector"
   supplied to This function has the same direction as the Euler
   axis, and its magnitude is the angle in radians.

2. If w is null, the unit matrix is returned.

3. The reference frame rotates clockwise as seen looking along the
   rotation vector from the origin.

"""
function rv2m(w)
    r = zeros((3, 3))
    ccall((:eraRv2m, liberfa), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}),
          w, r)
    r
end

"""
    rxr(a, b)

Multiply two r-matrices.

### Given ###

- `a`: First r-matrix
- `b`: Second r-matrix

### Returned ###

- `atb`: A * b

### Note ###

   It is permissible to re-use the same array for any of the
   arguments.

### Called ###

- `eraCr`: copy r-matrix

"""
function rxr(a, b)
    atb = zeros((3, 3))
    ccall((:eraRxr, liberfa),
          Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          a, b, atb)
    atb
end

"""
    rx(phi, r)

Rotate an r-matrix about the x-axis.

### Given ###

- `phi`: Angle (radians)

### Given and returned ###

- `r`: r-matrix, rotated

### Notes ###

1. Calling this function with positive phi incorporates in the
   supplied r-matrix r an additional rotation, about the x-axis,
   anticlockwise as seen looking towards the origin from positive x.

2. The additional rotation can be represented by this matrix:

       (  1        0            0      )
       (                               )
       (  0   + cos(phi)   + sin(phi)  )
       (                               )
       (  0   - sin(phi)   + cos(phi)  )

"""
rx

"""
    ry(phi, r)

Rotate an r-matrix about the y-axis.

### Given ###

- `theta`: Angle (radians)

### Given and returned ###

- `r`: r-matrix, rotated

### Notes ###

1. Calling this function with positive theta incorporates in the
   supplied r-matrix r an additional rotation, about the y-axis,
   anticlockwise as seen looking towards the origin from positive y.

2. The additional rotation can be represented by this matrix:

       (  + cos(theta)     0      - sin(theta)  )
       (                                        )
       (       0           1           0        )
       (                                        )
       (  + sin(theta)     0      + cos(theta)  )

"""
ry

"""
    rz(phi, r)

Rotate an r-matrix about the z-axis.

### Given ###

- `psi`: Angle (radians)

### Given and returned ###

- `r`: r-matrix, rotated

### Notes ###

1. Calling this function with positive psi incorporates in the
   supplied r-matrix r an additional rotation, about the z-axis,
   anticlockwise as seen looking towards the origin from positive z.

2. The additional rotation can be represented by this matrix:

       (  + cos(psi)   + sin(psi)     0  )
       (                                 )
       (  - sin(psi)   + cos(psi)     0  )
       (                                 )
       (       0            0         1  )

"""
rz

for name in ("rx",
             "ry",
             "rz")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(a, r)
            ccall(($fc, liberfa), Cvoid,
                  (Cdouble, Ptr{Cdouble}),
                  a, r)
            r
        end
    end
end

"""
    rxpv(r, pv)

Multiply a pv-vector by an r-matrix.

### Given ###

- `r`: R-matrix
- `pv`: Pv-vector

### Returned ###

- `rpv`: R * pv

### Note ###

   It is permissible for pv and rpv to be the same array.

### Called ###

- `eraRxp`: product of r-matrix and p-vector

"""
function rxpv(r, p)
    rp = zeros((2, 3))
    ccall((:eraRxpv, liberfa), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            r, p, rp)
    rp
end

"""
    rxp(r, p)

Multiply a p-vector by an r-matrix.

### Given ###

- `r`: R-matrix
- `p`: P-vector

### Returned ###

- `rp`: R * p

### Note ###

   It is permissible for p and rp to be the same array.

### Called ###

- `eraCp`: copy p-vector

"""
function rxp(r, p)
    rp = zeros(3)
    ccall((:eraRxp, liberfa), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            r, p, rp)
    rp
end


text:
@testset "4.6 Cosecant" begin
    include("4.6.0 (a csc)^m (b trg)^n.jl")
    include("4.6.1.2 (d csc)^n (a+b csc)^m.jl")
    include("4.6.1.3 (d cos)^n (a+b csc)^m.jl")
    include("4.6.1.4 (d cot)^n (a+b csc)^m.jl")
    include("4.6.11 (e x)^m (a+b csc(c+d x^n))^p.jl")
    include("4.6.3.1 (a+b csc)^m (d csc)^n (A+B csc).jl")
    include("4.6.4.2 (a+b csc)^m (d csc)^n (A+B csc+C csc^2).jl")
    include("4.6.7 (d trig)^m (a+b (c csc)^n)^p.jl")
end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/.jl

max_issues_repo_path: test/4 Trig functions/4.6 Cosecant/.jl

max_forks_repo_path: test/4 Trig functions/4.6 Cosecant/.jl


text:
# using Revise

##

using RoME

#  Do some plotting
# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)

function driveHex(fgl::G, posecount::Int) where G <: AbstractDFG
    # Drive around in a hexagon
    for i in (posecount-1):(posecount-1+5)
        psym = Symbol("x$i")
        posecount += 1
        nsym = Symbol("x$(i+1)")
        addVariable!(fgl, nsym, Pose2)
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
        addFactor!(fgl, [psym;nsym], pp, graphinit=false )
    end

    return posecount
end

function offsetHexLeg(fgl::G, posecount::Int; direction=:right) where G <: AbstractDFG
    psym = Symbol("x$(posecount-1)")
    nsym = Symbol("x$(posecount)")
    posecount += 1
    addVariable!(fgl, nsym, Pose2)
    pp = nothing
    if direction == :right
        pp = Pose2Pose2(MvNormal([10.0;0;-pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    elseif direction == :left
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    end
    addFactor!(fgl, [psym; nsym], pp, graphinit=false )
    return posecount
end

## start with an empty factor graph object
fg = initfg()
getSolverParams(fg).isfixedlag = true
getSolverParams(fg).qfl = 15
posecount = 0

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)
posecount += 1

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0.0; 0.0; 0.0],
                                           Matrix(Diagonal([0.1;0.1;0.05].^2))) ), graphinit=false )

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x0; :l1], p2br, graphinit=false )

## hex 1

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, graphinit=false )

# writeGraphPdf(fg, show=true)

# debugging options
getSolverParams(fg).drawtree = true
getSolverParams(fg).multiproc = false
getSolverParams(fg).dbg = true
getSolverParams(fg).limititers = 50
# getSolverParams(fg).downsolve = false
getSolverParams(fg).async = true

# do inference over factor graph
tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg))

fetchCliqTaskHistoryAll!(smt, hist)
assignTreeHistory!(tree, hist)

## test first cliq for inferdim

printCliqHistorySummary(tree, :x0)
# cliq6 = getClique(tree, :x0)

sfg = hist[6][9][4].cliqSubFg

isInitialized(sfg,:x0)
getVariableInferredDim(sfg,:x0)
isInitialized(fg,:x0)
getVariableInferredDim(fg,:x0)

isInitialized(sfg,:x1)
getVariableInferredDim(sfg,:x1)
isInitialized(fg,:x1)
getVariableInferredDim(fg,:x1)

isInitialized(sfg,:l1)
getVariableInferredDim(sfg,:l1)
isInitialized(fg,:l1)
getVariableInferredDim(fg,:l1)

## check the up message sent to

sfg = hist[6][9][4].cliqSubFg
upmsg = prepCliqInitMsgsUp(sfg, getClique(tree, :x0))
upMsg(tree, :x0)

# # FIXME these are deprecated
# dwnmsg = getDwnMsgs(tree, :x3)
# dwnmsg = getDwnMsgs(tree, :x1)
# dwnmsg = getDwnMsgs(tree, :x0)

printCliqHistorySummary(tree, :x1)

## test second inference/init clique for correct inferred dimensions

printCliqHistorySummary(tree, :x2)
getClique(tree, :x2)

sfg = hist[5][17][4].cliqSubFg

isInitialized(sfg,:x1)
getVariableInferredDim(sfg,:x1)
isInitialized(fg,:x1)
getVariableInferredDim(fg,:x1)

isInitialized(sfg,:x3)
getVariableInferredDim(sfg,:x3)
isInitialized(fg,:x3)
getVariableInferredDim(fg,:x3)

isInitialized(sfg,:x2)
getVariableInferredDim(sfg,:x2)
isInitialized(fg,:x2)
getVariableInferredDim(fg,:x2)

sfg = hist[5][17][4].cliqSubFg
upmsg = prepCliqInitMsgsUp(sfg, getClique(tree, :x2))
upMsg(tree, :x2)

## test root inference/init clique for correct inferred dimensions

printCliqHistorySummary(tree, :x3)
getClique(tree, :x3)

sfg = hist[1][10][4].cliqSubFg

isInitialized(sfg,:x3)
getVariableInferredDim(sfg,:x3)
isInitialized(fg,:x3)
getVariableInferredDim(fg,:x3)

isInitialized(sfg,:x5)
getVariableInferredDim(sfg,:x5)
isInitialized(fg,:x5)
getVariableInferredDim(fg,:x5)

isInitialized(sfg,:l1)
getVariableInferredDim(sfg,:l1)
isInitialized(fg,:l1)
getVariableInferredDim(fg,:l1)

sfg = hist[1][10][4].cliqSubFg
upmsg = prepCliqInitMsgsUp(sfg, getClique(tree, :x3))
upMsg(tree, :x3)

## test root inference/init clique for correct inferred dimensions

printCliqHistorySummary(tree, :x1)
getClique(tree, :x1)

sfg = hist[2][16][4].cliqSubFg
getVariableInferredDim(sfg,:x1)

sfg = hist[2][17][4].cliqSubFg
getVariableInferredDim(sfg,:x1)

stuff = sandboxCliqResolveStep(tree,:x1,16)

isInitialized(fg,:x1)
getVariableInferredDim(fg,:x1)

isInitialized(sfg,:l1)
getVariableInferredDim(sfg,:l1)
isInitialized(fg,:l1)
getVariableInferredDim(fg,:l1)

isInitialized(sfg,:x3)
getVariableInferredDim(sfg,:x3)
isInitialized(fg,:x3)
getVariableInferredDim(fg,:x3)

sfg = hist[2][16][4].cliqSubFg
upmsg = prepCliqInitMsgsUp(sfg, getClique(tree, :x1))
getUpMsgs(tree, :x3)

# make sure all infer dims are properly set

isInitialized(fg, :x0)
getVariableInferredDim(fg, :x0)

isInitialized(fg, :x1)
getVariableInferredDim(fg, :x1)

isInitialized(fg, :x2)
getVariableInferredDim(fg, :x2)

isInitialized(fg, :x3)
getVariableInferredDim(fg, :x3)

isInitialized(fg, :x4)
getVariableInferredDim(fg, :x4)

isInitialized(fg, :x5)
getVariableInferredDim(fg, :x5)

isInitialized(fg, :x6)
getVariableInferredDim(fg, :x6)

isInitialized(fg, :l1)
getVariableInferredDim(fg, :l1)

## Hex 2

posecount = offsetHexLeg(fg, posecount, direction=:right)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l2, Point2, tags=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l2], p2br, graphinit=false )

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l2], p2br2, graphinit=false )

# writeGraphPdf(fg,show=true)

getSolverParams(fg).downsolve = false

tree, smt, hist = solveTree!(fg, tree, recordcliqs=ls(fg))

## debugging

# fetchCliqTaskHistoryAll!(smt, hist)
# assignTreeHistory!(tree, hist)
# printCliqHistorySummary(hist, tree, :x10)
# makeCsmMovie(fg, tree)

# down msg to :x10
dwinmsgs = prepCliqInitMsgsDown!(fg, tree, getClique(tree, :x9), getClique(tree,:x11))

# up msg to :x6
getClique(tree, :x6)
printCliqHistorySummary(hist, tree, :x6)
sfg = hist[11][9][4].cliqSubFg
# writeGraphPdf(sfg, show=true)
upmsg = prepCliqInitMsgsUp(sfg, getClique(tree, :x6))

## PROBLEM: initializes and solves, but inferdim not populated
# deeper issue is that neither x5 or x6 have valid inferdim values
sfg = hist[11][8][4].cliqSubFg
isInitialized(sfg,:x7)
getVariableInferredDim(sfg,:x7)
isInitialized(sfg,:x5)
getVariableInferredDim(sfg,:x5)
sfg = hist[11][9][4].cliqSubFg
isInitialized(sfg,:x7)
getVariableInferredDim(sfg,:x7)

#


text:
# ========== Transforms On A Group Of Signals ==========
# ----- ACDWT on a set of signals -----
"""
    acdwtall(x, wt[, L])

Computes the autocorrelation discrete wavelet transform (ACDWT) on each slice of signal.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input `N-1`-D signals, where each signal is sliced
  at dimension `N`.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `minimum(size(xw)[1:end-1]) |> maxtransformlevels`) Number of
  levels of wavelet transforms.

# Returns
- `::Array{T}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACDWT on all signals in x
xw = acdwtall(x, wt)
```

**See also:** [`acdwt`](@ref)
"""
function acdwtall(x::AbstractArray{T},
                  wt::OrthoFilter,
                  L::Integer = minimum(size(x)[1:end-1]) |> maxtransformlevels) where
                  T<:Number
    @assert 2 ≤ ndims(x) ≤ 3
    # Allocate space for transforms
    sz = size(x)[1:(end-1)]
    k = ndims(x)==2 ? L+1 : 3*L+1
    N = size(x)[end]
    xw = Array{T}(undef, (sz...,k,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            acdwt!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iacdwtall(xw[, wt])

Computes the inverse autocorrelation discrete wavelet transform (IACDWT) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: ACDWT-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
- `::Array{T}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACDWT on all signals in x
xw = acdwtall(x, wt)

# IACDWT on all signals in xw
x̂ = iacdwtall(xw)
```

**See also:** [`iacdwt`](@ref)
"""
function iacdwtall(xw::AbstractArray{T}, wt::Union{OrthoFilter,Nothing} = nothing) where
                   T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iacdwt!(xᵢ, xwᵢ)
        end
    end
    return x
end

# ----- ACWPT on a set of signals -----
"""
    acwptall(x, wt[, L])

Computes the autocorrelation wavelet packet transform (ACWPT) on each slice of signal.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input `N-1`-D signals, where each signal is sliced
  at dimension `N`.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `minimum(size(xw)[1:end-1]) |> maxtransformlevels`) Number of
  levels of wavelet transforms.

# Returns
- `::Array{T}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPT on all signals in x
xw = acwptall(x, wt)
```

**See also:** [`acwpt`](@ref)
"""
function acwptall(x::AbstractArray{T},
                  wt::OrthoFilter,
                  L::Integer = minimum(size(x)[1:end-1]) |> maxtransformlevels) where
                  T<:Number
    @assert 2 ≤ ndims(x) ≤ 3
    # Allocate space for transforms
    sz = size(x)[1:(end-1)]
    k = ndims(x)==2 ? 1<<L : Int(4^L)
    N = size(x)[end]
    xw = Array{T}(undef, (sz...,k,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            acwpt!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iacwptall(xw[, wt])

Computes the inverse autocorrelation wavelet packet transform (IACWPT) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: ACWPT-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
- `::Array{T}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPT on all signals in x
xw = acwptall(x, wt)

# IACWPT on all signals in xw
x̂ = iacwptall(xw)
```

**See also:** [`iacwpt`](@ref)
"""
function iacwptall(xw::AbstractArray{T}, wt::Union{OrthoFilter,Nothing} = nothing) where
                   T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iacwpt!(xᵢ, xwᵢ)
        end
    end
    return x
end

# ----- ACWPD on a set of signals -----
"""
    acwpdall(x, wt[, L])

Computes the autocorrelation wavelet packet decomposition (ACWPD) on each slice of signal.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input `N-1`-D signals, where each signal is sliced
  at dimension `N`.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `minimum(size(xw)[1:end-1]) |> maxtransformlevels`) Number of
  levels of wavelet transforms.

# Returns
- `::Array{T}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPD on all signals in x
xw = acwpdall(x, wt)
```

**See also:** [`acwpd`](@ref)
"""
function acwpdall(x::AbstractArray{T},
                  wt::OrthoFilter,
                  L::Integer = minimum(size(x)[1:end-1]) |> maxtransformlevels) where
                  T<:Number
    @assert 2 ≤ ndims(x) ≤ 3
    # Allocate space for transforms
    sz = size(x)[1:(end-1)]
    k = ndims(x)==2 ? 1<<(L+1)-1 : sum(4 .^(0:L))
    N = size(x)[end]
    xw = Array{T}(undef, (sz...,k,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            acwpd!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iacwpdall(xw[, wt, L])
    iacwpdall(xw, L)
    iacwpdall(xw, wt, tree)
    iacwpdall(xw, tree)

Computes the inverse autocorrelation wavelet packet decomposition (IACWPD) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: ACWPD-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.
- `L::Integer`: (Default: `minimum(size(xw)[1:end-2]) |> maxtransformlevels`) Number of
  levels of wavelet transforms.
- `tree::BitVector`: Binary tree for inverse transform to be computed accordingly.

# Returns
- `::Array{T}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPD on all signals in x
xw = acwpdall(x, wt)

# IACWPD on all signals in xw
x̂ = iacwpdall(xw)
x̂ = iacwpdalll(xw, maketree(x))
x̂ = iacwpdall(xw, 5)
```

**See also:** [`iacwpd`](@ref)
"""
function iacwpdall(xw::AbstractArray{T},
                   wt::Union{OrthoFilter,Nothing} = nothing,
                   L::Integer = minimum(size(xw)[1:end-2]) |> maxtransformlevels) where
                   T<:Number
    return iacwpdall(xw, L)
end

function iacwpdall(xw::AbstractArray{T}, L::Integer) where T<:Number
    return iacwpdall(xw, maketree(size(xw)[1:(end-2)]..., L))
end

function iacwpdall(xw::AbstractArray{T},
                   wt::Union{OrthoFilter, Nothing},
                   tree::BitVector) where T<:Number
    return iacwpdall(xw, tree)
end

function iacwpdall(xw::AbstractArray{T}, tree::BitVector) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iacwpd!(xᵢ, xwᵢ, tree)
        end
    end
    return x
end