text:
using Primes

# Divide and multiply integers and rationals
# like a CAS does. Always return exact result (no floats).
# Return Int if possible. We want to avoid infection Rationals
# and putting explicit conversions everywhere. Examples:

#  12/6 --> 2, not 2.0
#  13/6 --> 13//6
#  (13//6) * 6 --> 13
rat_to_int(r::Rational{T}) where {T<:Integer} = r.den == 1 ? r.num : r

@doc """
    mmul(x,y), mpow(x,y), mplus(x,y), mminus(x), mdiv(x,y)
arithmetic funtions for Symata. These are similar to `*`, `^`, `+`, etc. in Julia. But there are important differences.

- rational numbers with denominator equal to `1` are converted to an integer type.
- complex numbers with imaginary part equal to zero are converted to a real type.
- for symbolic arguments, or mixed symbolic-numeric arguments a symata expression of type `Mxpr` is constructed.
- integer factors are extracted from roots. (overzealously at the moment).
""" arithmetic

@doc (@doc arithmetic) mmul
@doc (@doc arithmetic) mplus
@doc (@doc arithmetic) mpow
@doc (@doc arithmetic) mdiv
@doc (@doc arithmetic) mminus

## FIXME: Int or Integer ?
mmul(x::Int, y::Rational{T}) where {T<:Integer} = (res = x * y; return res.den == 1 ? res.num : res )

mmul(x::Complex{V}, y::Rational{T}) where {T<:Integer,V<:Integer} = (res = x * y; return (res.im.den == 1 && res.re.den == 1) ? Complex(res.re.num,res.im.num) : res )
mmul(x::Rational{V}, y::Complex{T}) where {T<:Integer,V<:Integer} = (res = x * y; return (res.im.den == 1 && res.re.den == 1) ? Complex(res.re.num,res.im.num) : res )

## ???
## mmul{T<:Integer,U<:Rational}(x::Complex{T}, y::Complex{U}) = mmul(y,x)

function mmul(x::Complex{U}, y::Complex{T}) where {T<:Integer,U<:Rational}
    res = x * y
    isa(res,Complex) && return res
    res.den == 1 && return res.num
    res
end

## Julia returns different types depending on the order of x and y
function mmul(x::Complex{T}, y::Complex{U}) where {T<:Integer,U<:Rational}
    res = x * y
    if isa(res,Complex)
        imag(res).num != 0 && return res
        real(res).den == 1 && return real(res).num
    end
    res
end

mmul(x::Rational{T}, y::Int) where {T<:Integer} = (res = x * y; return res.den == 1 ? res.num : res )

# Oct 2016. Added generic mxpr(:Times,x,y). Does not break tests.
# Why ? This allows Compiled (ie. Julia) functions to take far more arguments.
# FIXME. Use mmul in Symata code instead of *. We want to remove undesired methods for *
mmul(x::T,y::V) where {T<:Number, V<:Number} = x * y
mmul(x,y) = mxpr(:Times, x, y)

mplus(x::Rational{T}, y::Rational{V}) where {T<:Integer,V<:Integer} = rat_to_int(x+y)

mplus(x::T,y::V) where {T<:Number, V<:Number} = x + y
mplus(x,y) = mxpr(:Plus, x, y)

mminus(x) = mxpr(:Times,-1,x)
mminus(x::Number) = -x

mminus(x,y) = mxpr(:Plus, x, mxpr(:Times,-1,y))
mminus(x::Number,y::Number) = x-y

# mdiv is apparently used in do_Rational, which seems to be the only entry point to mdiv.
mdiv(x::Integer, y::Integer) =  y == 0 ? ComplexInfinity : rem(x,y) == 0 ? div(x,y) : x // y
_mdiv_float(x, y) = y == 0 ? ComplexInfinity : div(x,y)
mdiv(x::Integer, y::AbstractFloat) = _mdiv_float(x, y)
mdiv(x::AbstractFloat, y::Integer) = _mdiv_float(x, y)
mdiv(x::AbstractFloat, y::AbstractFloat) = _mdiv_float(x, y)

mdiv(x::Integer, y::Rational) = (res = x / y; return res.den == 1 ? res.num : res )
mdiv(x::Number,y::Number) = x/y
# Should be a fixed Mxpr ?
mdiv(x,y) = mxpr(:Times, x, mxpr(:Power,y,-1))

function mpow(x,y)
#    @show (x,y)
    _mpow(x,y)
end
function _mpow(x::Number, y::Number)
    try
        x^y
    catch e
        println("Error computing $x^$y")
        rethrow(e)
    end
end
function _mpow(x::Rational,y::Integer)
    if x > 0 && y < 0
        my = -y
        (denominator(x)^(my))//(numerator(x)^(my))
    else
        x^y  # will fail
    end
end
function _mpow(x::Irrational, n::Number)
    n >= 0 && return x^n
    (x1,n1) = promote(x,n)
    x1^n1
end

_mpow(x,y) =  mxpr(:Power, x, y)

function _mpow(x::Integer, y::V) where V<:AbstractFloat
    x >= 0 && return convert(V,x)^y
    res = convert(Complex{V},x)^y
    imag(res) == 0 && return real(res)
    res
end

_mpow(x::Complex{T},y::V) where {T<:Integer,V<:AbstractFloat} = convert(Complex{V},x)^y
_mpow(x::Integer, y::Integer) = y >= 0 ? x^y : x == 0 ? ComplexInfinity : 1//(x^(-y))
_mpow(x::AbstractFloat, y::Integer) = y >= 0 ? x^y : x == 0 ? ComplexInfinity : x^y
_mpow(x::Symata.Mxpr{:DirectedInfinity}, y::T) where {T<:Real} = y > 0 ? x : 0
_mpow(x::Symata.Mxpr{:DirectedInfinity}, y::Complex{T}) where {T<:Real} = y.re > 0 ? x : 0
# could be more efficient: cospi(exp) + im * sinpi(exp)
_mpow(b::T,exp::V) where {T<:AbstractFloat,V<:Number} = b < 0 ? (cospi(exp) + im * sinpi(exp)) * abs(b)^exp : b^exp

# FIXME: This code and function do_Power{T<:Integer}(mx::Mxpr{:Power},b::T,e::Rational)
# in apprules are in conflict.
# Maybe Fixed. We are using this routine.
# Note: do_Power{T<:Integer}(mx::Mxpr{:Power},b::T,e::Rational) is disabled.
# We could research how to do this instead of rolling our own. But, this seems to
# work.
# function _mpow(x::T,y::Rational{V}) where {T<:Integer, V<:Integer}
function _mpow(x::Integer, y::Rational{V}) where {V<:Integer}
    x == 1 && return x
    local gotneg::Bool
    if x < 0
        x = -x
        gotneg = true
    else
        gotneg = false
    end
    facs = factor(x)
    newfacs = newargs()
    for (fac,mul) in facs
        (q,r) = divrem(mul,y.den)
        if r == 0
            nf = fac^(q*abs(y.num))
            nf != 1 ? y.num > 0 ? push!(newfacs,nf) : push!(newfacs,1//nf) : nothing
            continue
        else
            nf = fac^q
            nf != 1 ? y.num > 0 ? push!(newfacs,nf) : push!(newfacs,1//nf) : nothing
        end
        (nf1,r1) = divrem(r*y.num,y.den)
        if nf1 != 0
            newfac = mpow(fac,nf1)
            push!(newfacs,newfac)
        end
        if r1 != 0
            newexp = r1//(y.den)
            push!(newfacs,mxpr(:Power,fac,newexp)) # Optimize. at some point, do mxprfc and test
        end
    end
    if gotneg == true  # This at least not wrong, but it's not like Mma for y < 0
        (n,r) = divrem(y.num,y.den)
        fac = y.den == 2 ? I : mxpr(:Power, -1, r//y.den)
        if iseven(n)
            push!(newfacs, fac)
        else
            push!(newfacs,mxpr(:Times, -1, fac))
        end
    end
    newfacs = mpow_test_for_integer_factors(newfacs)
    length(newfacs) == 1 && return newfacs[1]
    # We may need to doeval this to use in Julia functions. Or doeval in the function. Eg.
    # pattern to test for perfect squares
    # t = _:?(:( (x) -> typeof(mpow(x,1//2)) <: Integer ))
    expr = mxpr(:Times, newfacs...)
    expr
end

# If all factors are integers, just multiply them
function mpow_test_for_integer_factors(facs)
    allints = true
    for x in facs
        typeof(x) <: Integer && continue
        allints = false
        break
    end
    return allints ? [*(facs...)] : facs
end

# The following expression will call this method:  (-27/64)^(2/3) == 9/16*((-1)^(2/3))
# Also called by Norm([1,I/2])
_mpow(x::Rational{T}, y::Rational{V}) where {T<:Integer, V<:Integer} =  mxpr(:Times,mpow(x.num,y), mxpr(:Power,mpow(x.den,y), -1))

### Code below handles z^n for z complex and n real.

# Is there not an automatic way ?
_divide_pow(x::T,y::T) where {T<:AbstractFloat} = x/y
_divide_pow(x::Complex{T},y::T) where {T<:AbstractFloat} = x/y
_divide_pow(x::T,y::T) where {T<:Integer} = x//y
_divide_pow(x::Complex{T},y::T) where {T<:Integer} = x//y
_divide_pow(x::Complex{Rational{T}},y::Rational{T}) where {T<:Integer} = x//y

function _mpow(base::Complex{T}, expt::V) where {T<:Real,V<:Union{AbstractFloat, Integer}}
    expt >= 0 && return base^expt
    if expt == -one(expt)
        if base == complex(zero(real(base)),one(real(base)))
            return complex(zero(real(base)),-one(real(base)))
        else
            return _divide_pow(conj(base),abs2(base))
        end
    else
        res = _divide_pow((conj(base)^(-expt)) , abs2(base)^(-expt))
        if imag(res) == zero(imag(res))
            return real(res)
        else
            return res
        end
    end
end

#### mabs

mabs(x::Number) = abs(x)

mabs(x) = mxpr(:Abs, x)

# TODO.  T Rational
function mabs(x::Complex{T}) where T<:Integer
    r,i = reim(x)
    sq = r*r + i*i
    x = isqrt(sq)
    if x * x == sq
        return x
    else
        res = mxpr(:Power,sq,1//2)
        setfixed(res)
        mergesyms(res,:nothing)
        return res
    end
end

# mpow for Rational is wrong often wrong.
# so this will be broken, too.
function mabs(x::Complex{Rational{T}}) where T<:Integer
    r,i = reim(x)
    sq = mplus(mmul(r,r),mmul(i,i))
    return mpow(sq,1//2)
end

## copied from base/operators, a great trick
for op = (:mplus, :mmul)
    @eval begin
        # note: these definitions must not cause a dispatch loop when +(a,b) is
        # not defined, and must only try to call 2-argument definitions, so
        # that defining +(a,b) is sufficient for full functionality.

        # ($op){T<:Number}(a::T, b::T, c::T)        = ($op)(($op)(a,b),c)
        # ($op){T<:Number}(a::T, b::T, c::T, xs::T...) = ($op)(($op)(($op)(a,b),c), xs...)

        # Oct 2016 Compiled Symata translates mxpr(:Times,x,y,z,...) to mmul(x,y,z,...).
        # In order for this to work, we allow any combination of numbers.
        ($op)(a::Number, b::Number, c::Number)        = ($op)(($op)(a,b),c)
        ($op)(a::Number, b::Number, c::Number, xs::Number...) = ($op)(($op)(($op)(a,b),c), xs...)
        # a further concern is that it's easy for a type like (Int,Int...)
        # to match many definitions, so we need to keep the number of
        # definitions down to avoid losing type information.
    end
end

## `Compile` and other functions may generate `mmul(x,y,z)` where some of the args are not numbers.
## These catch cases where some of args are not numbers.
## This may want to implment more logic for efficiency
mmul(args...) = mxpr(:Times, args...)
mplus(args...) = mxpr(:Plus, args...)


text:
abstract type Basis end

# Assumes Unit vectors
struct CubicBasis <: Basis
    positions::Vector{Point3f0}
end

"""
    sc(pos = Point3f0(0.))
"""
sc(pos::Point3f0 = Point3f0(0.)) = CubicBasis([pos])

"""
    bcc(pos = Point3f0(0.))
"""
bcc(pos::Point3f0 = Point3f0(0.)) = CubicBasis([pos, pos + Point3f0(0.5)])

"""
    fcc(pos = Point3f0(0.))
"""
function fcc(pos::Point3f0 = Point3f0(0.))
    CubicBasis([
        pos,
        pos + Point3f0(0.5, 0.5, 0.0),
        pos + Point3f0(0.0, 0.5, 0.5),
        pos + Point3f0(0.5, 0.0, 0.5)
    ])
end

"""
    diamond(pos = Point3f0(0.))
"""
function diamond(pos::Point3f0 = Point3f0(0.))
    CubicBasis([
        pos,
        pos + Point3f0(0.5, 0.5, 0.0),
        pos + Point3f0(0.0, 0.5, 0.5),
        pos + Point3f0(0.5, 0.0, 0.5),
        pos + Point3f0(0.25, 0.25, 0.25),
        pos + Point3f0(0.75, 0.75, 0.25),
        pos + Point3f0(0.25, 0.75, 0.75),
        pos + Point3f0(0.75, 0.25, 0.75)
    ])
end

"""
    merge(cubics::CubicBasis...)
"""
function merge(cubics::CubicBasis...)
    positions = Point3f0[]
    for c in cubics
        for p in c.positions
            if !(p in positions)
                push!(positions, p)
            end
        end
    end

    Cubic(positions)
end

function center!(cubic::CubicBasis; epsilon::Float32 = 1f-10)
    es = [Point3f0(1, 0, 0), Point3f0(0, 1, 0), Point3f0(0, 0, 1)]
    for i in eachindex(cubic.positions)
        for j in eachindex(cubic.positions[i])
            while cubic.positions[i][j] > 0.5f0 + epsilon
                cubic.positions[i] = cubic.positions[i] - es[j]
            end
            while cubic.positions[i][j] < -0.5f0 - epsilon
                cubic.positions[i] = cubic.positions[i] + es[j]
            end
        end
    end
    nothing
end

function expand(
        cubic::CubicBasis;
        min::Point3f0 = Point3f0(-1.0),
        max::Point3f0 = Point3f0(1.0),
        epsilon::Float32 = 1f-10
    )

    center!(cubic)
    xrange = floor(Int64, min[1]):ceil(Int64, max[1])
    yrange = floor(Int64, min[2]):ceil(Int64, max[2])
    zrange = floor(Int64, min[3]):ceil(Int64, max[3])

    positions = Point3f0[]
    for x in xrange, y in yrange, z in zrange
        r = Point3f0(x, y, z)
        for p in cubic.positions
            v = p + r
            if all(map((a, x, b) -> a-epsilon <= x <= b+epsilon, min, v, max))
                push!(positions, p + r)
            end
        end
    end

    positions
end

ext: jl


text:
using SpatialFields
using Base.Test

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

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

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

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

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

@testset "linear fit" begin
    include("linear_fit.jl")
end

@testset "benchmarks" begin
    include("../perf/runbenchmarks.jl")
end


text:
@testset "circular_turn_buffer" begin
    @testset "scalar state" begin
        b = circular_RTSA_buffer(; capacity = 2, state_eltype = Int, state_size = ())

        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true

        push!(b; reward = 0.0, terminal = true, state = 1, action = 1)

        @test length(b) == 0  # a turn is not finished yet
        @test isempty(b) == false

        push!(b; reward = 0.3, terminal = true, state = 2, action = 2)
        push!(b; reward = 0.6, terminal = false, state = 3, action = 3)

        @test length(b) == 2
        @test isempty(b) == false
        @test isfull(b) == true

        # make sure state/action is not overwritten yet
        @test b[1] == (
            state = 1,
            action = 1,
            reward = 0.3f0,
            terminal = true,
            next_state = 2,
            next_action = 2,
        )
        @test b[end] == (
            state = 2,
            action = 2,
            reward = 0.6f0,
            terminal = false,
            next_state = 3,
            next_action = 3,
        )

        push!(b; reward = 0.9, terminal = true, state = 4, action = 4)

        # make sure that old elements are removed
        @test b[1] == (
            state = 2,
            action = 2,
            reward = 0.6f0,
            terminal = false,
            next_state = 3,
            next_action = 3,
        )
        @test b[end] == (
            state = 3,
            action = 3,
            reward = 0.9f0,
            terminal = true,
            next_state = 4,
            next_action = 4,
        )
    end

    @testset "2d state" begin
        b = circular_RTSA_buffer(
            ;
            capacity = 3,
            state_eltype = Int,
            state_size = (2, 2),
        )

        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true

        push!(b; reward = 1.0, terminal = false, state = [1 1; 1 1], action = 1)
        push!(b; reward = 1.0, terminal = false, state = [2 2; 2 2], action = 2)

        @test length(b) == 1
        @test isfull(b) == false
        @test isempty(b) == false

        push!(b; reward = 2.0, terminal = false, state = [3 3; 3 3], action = 3)
        push!(b; reward = 3.0, terminal = false, state = [4 4; 4 4], action = 4)

        @test length(b) == 3
        @test isfull(b) == true
        @test isempty(b) == false
        @test b[1] == (
            state = [1 1; 1 1],
            action = 1,
            reward = 1.0f0,
            terminal = false,
            next_state = [2 2; 2 2],
            next_action = 2,
        )
        @test b[end] == (
            state = [3 3; 3 3],
            action = 3,
            reward = 3.0f0,
            terminal = false,
            next_state = [4 4; 4 4],
            next_action = 4,
        )

        push!(b; reward = 4.0, terminal = false, state = [5 5; 5 5], action = 5)

        # old experience should be removed
        @test length(b) == 3
        @test isfull(b) == true
        @test b[1] == (
            state = [2 2; 2 2],
            action = 2,
            reward = 2.0f0,
            terminal = false,
            next_state = [3 3; 3 3],
            next_action = 3,
        )
        @test b[end] == (
            state = [4 4; 4 4],
            action = 4,
            reward = 4.0f0,
            terminal = false,
            next_state = [5 5; 5 5],
            next_action = 5,
        )
    end
end


text:
@testset "build_single_discrimiantor" begin
    noise_size = (64, 64, 3, 10)
    noise = randn(Float32, noise_size)
    dscr = SinGAN.build_single_discriminator(5, 128)
    @info dscr
    @test size(dscr(noise)) == (64, 64, 1, 10)
end

@testset "build_single_gen_layers" begin
    noise_size = (64, 64, 3, 10)
    noise = randn(Float32, noise_size)
    gen = SinGAN.build_single_gen_layers(5, 128)
    @info gen
    @test size(gen(noise)) == (64, 64, 3, 10)
end

@testset "DiscriminatorPyramid" begin
    @info DiscriminatorPyramid(4, 5)
end

@testset "NoiseConnection" begin
    nc = NoiseConnection(Conv((3, 3), 3 => 3, pad = 1), 5)
    @info nc
    prev = randn(Float32, 42, 42, 3, 1)
    noise = randn(Float32, 42, 42, 3, 1)
    @test size(nc(prev, noise)) == (32, 32, 3, 1)
    # @code_warntype nc(prev, noise)
end

# @testset "NoiseConnectionInferType" begin
#     nc = NoiseConnection(Conv((3, 3), 3 => 3, pad = 1), 5)
#     prev = randn(Float32, 128, 128, 3, 1)
#     noise = randn(Float32, 128, 128, 3, 1)
#     @time for i in 1:100 nc(prev, noise) end
# end

@testset "GeneratorPyramid" begin
    image_shapes = [(32, 32), (44, 44)]
    genp = GeneratorPyramid(image_shapes, 5)
    @info genp
    @test genp.noise_shapes == [(42, 42), (54, 54)]
    xs = [randn(Float32, 42, 42, 3, 1), randn(Float32, 54, 54, 3, 1)] |> gpu
    @test size(genp(xs, 0, false)) == (32, 32, 3, 1)
    @test size(genp(xs, 0, true)) == (42, 42, 3, 1)
    @test size(genp(xs, 1, false)) == (32, 32, 3, 1)
    @test size(genp(xs, 1, true)) == (54, 54, 3, 1)
    @test size(genp(xs, 2, false)) == (44, 44, 3, 1)
    # @code_warntype genp(xs, 2, false)
end

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
using Pkg
Pkg.activate(@__DIR__)

# For linear algebra
using StaticArrays: SVector
using LinearAlgebra
using SparseArrays: SparseMatrixCSC, sparse
# Fast Cholesky
using KoLesky

# multiple dispatch
import IterativeSolvers: cg!
import LinearAlgebra: mul!, ldiv!
import Base: size

# parser
using ArgParse

# logging
using Logging

using JLD

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--alpha"
            help = "α"
            arg_type = Float64
            default = 1.0
        "--m"
            help = "m"
            arg_type = Int
            default = 3
        "--kernel"
            arg_type = String
            default = "Matern7half"
            # default = "Gaussian"
        "--sigma"
            help = "lengthscale"
            arg_type = Float64
            default = 0.3
        "--h"
            help = "grid size"
            arg_type = Float64
            default = 0.02
        "--nugget"
            arg_type = Float64
            default = 1e-10
        "--GNsteps"
            arg_type = Int
            default = 3
        "--rho_big"
            arg_type = Float64
            default = 3.0
        "--rho_small"
            arg_type = Float64
            default = 3.0
        "--k_neighbors"
            arg_type = Int
            default = 3
        "--compare_exact"
            arg_type = Bool
            default = false
    end
    return parse_args(s)
end

## PDEs type
abstract type AbstractPDEs end
struct NonlinElliptic2d{Tα,Tm,TΩ} <: AbstractPDEs
    # eqn: -Δu + α*u^m = f in [Ω[1,1],Ω[2,1]]*[Ω[1,2],Ω[2,2]]
    α::Tα
    m::Tm
    Ω::TΩ
    bdy::Function
    rhs::Function
end

## sample points
function sample_points_rdm(eqn::NonlinElliptic2d, N_domain, N_boundary)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]

    X_domain = hcat(rand(Float64,(N_domain, 1))*(x1r-x1l).+x1l,rand(Float64,(N_domain, 1))*(x2r-x2l).+x2l)

    N_bd_each=convert(Int64, N_boundary/4)
    if N_boundary != 4 * N_bd_each
        println("[sample points] N_boundary not divided by 4, replaced by ", 4 * N_bd_each)
        N_boundary = 4 * N_bd_each
    end

    X_boundary = zeros((N_boundary, 2))
    # bottom face
    X_boundary[1:N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2l*ones(N_bd_each))
    # right face
    X_boundary[N_bd_each+1:2*N_bd_each, :] = hcat(x1r*ones(N_bd_each),(x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    # top face
    X_boundary[2*N_bd_each+1:3*N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2r*ones(N_bd_each))

    # left face
    X_boundary[3*N_bd_each+1:N_boundary, :] = hcat(x1l*ones(N_bd_each), (x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    return X_domain', X_boundary'
end
function sample_points_grid(eqn::NonlinElliptic2d, h_in, h_bd)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]
    x = x1l + h_in:h_in:x1r-h_in
    y = x2l + h_in:h_in:x2r-h_in
    X_domain = reduce(hcat,[[x[i], y[j]] for i in 1:length(x) for j in 1:length(x)])

    l = length(x1l:h_bd:x1r-h_bd)
    X_boundary = vcat([x1l:h_bd:x1r-h_bd x2l*ones(l)], [x1r*ones(l) x2l:h_bd:x2r-h_bd], [x1r:-h_bd:x1l+h_bd x2r*ones(l)], [x1l*ones(l) x2r:-h_bd:x1l+h_bd])
    return X_domain, X_boundary'
end

## exact algorithms
# assemby Gram matrices
function get_Gram_matrices(eqn::NonlinElliptic2d, cov::KoLesky.AbstractCovarianceFunction, X_domain, X_boundary, sol_now)
    d = 2
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    Δδ_coefs = -1.0
    δ_coefs_int = eqn.α*eqn.m*(sol_now.^(eqn.m-1))

    # get linearized PDEs correponding measurements
    meas_δ = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_Δδ = [KoLesky.ΔδPointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), Δδ_coefs, δ_coefs_int[i]) for i = 1:N_domain]
    meas_test_int = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]

    Theta_train = zeros(N_domain+N_boundary,N_domain+N_boundary)

    measurements = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,2)
    measurements[1] = meas_δ; measurements[2] = meas_Δδ
    cov(Theta_train, reduce(vcat,measurements))

    Theta_test = zeros(N_domain,N_domain+N_boundary)
    cov(view(Theta_test,1:N_domain,1:N_boundary), meas_test_int, meas_δ)
    cov(view(Theta_test,1:N_domain,N_boundary+1:N_domain+N_boundary), meas_test_int, meas_Δδ)
    return Theta_train, Theta_test

end
# iterative GPR
function iterGPR_exact(eqn, cov, X_domain, X_boundary, sol_init, nugget, GNsteps)
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    # get the rhs and bdy data
    rhs = [eqn.rhs(X_domain[:,i]) for i in 1:N_domain]
    bdy = [eqn.bdy(X_boundary[:,i]) for i in 1:N_boundary]
    v = sol_init

    for _ in 1:GNsteps
        Theta_train, Theta_test = get_Gram_matrices(eqn, cov, X_domain, X_boundary, v)
        rhs_now = vcat(bdy, rhs+eqn.α*(eqn.m-1)*v.^eqn.m)

        v = Theta_test*(((Theta_train+nugget*diagm(diag(Theta_train))))\rhs_now)
    end
    return v
end

## algorithm using Kolesky and pcg
# struct that stores the factor of Theta_train
abstract type implicit_mtx end
struct approx_Theta_train{Tv,Ti,Tmtx<:SparseMatrixCSC{Tv,Ti}} <: implicit_mtx
    P::Vector{Ti}
    U::Tmtx
    L::Tmtx
    δ_coefs::Vector{Tv}
    N_boundary::Ti
    N_domain::Ti
end
struct precond_Theta_train{Tv,Ti,Tmtx<:SparseMatrixCSC{Tv,Ti}} <: implicit_mtx
    P::Vector{Ti}
    U::Tmtx
    L::Tmtx
end

function size(A::approx_Theta_train, num)
    return size(A.U,num)
end

function mul!(x, Θtrain::approx_Theta_train, b)
    @views temp = vcat(b[1:Θtrain.N_boundary],Θtrain.δ_coefs.*b[Θtrain.N_boundary+1:end],b[Θtrain.N_boundary+1:end])
    temp[Θtrain.P] = Θtrain.L\(Θtrain.U\temp[Θtrain.P])

    @views x[1:Θtrain.N_boundary] = temp[1:Θtrain.N_boundary]
    @views x[Θtrain.N_boundary+1:end] .= Θtrain.δ_coefs.*(temp[Θtrain.N_boundary+1:Θtrain.N_boundary+Θtrain.N_domain]) .+ temp[Θtrain.N_boundary+Θtrain.N_domain+1:end]
end

function ldiv!(x, precond::precond_Theta_train, b)
    x[precond.P] = precond.U*(precond.L*b[precond.P])
end

function iterGPR_fast_pcg(eqn, cov, X_domain, X_boundary, sol_init, nugget, GNsteps; ρ_big =4.0, ρ_small=6.0, k_neighbors = 4, lambda = 1.5, alpha = 1.0)

    @info "[Algorithm]: iterative GPR + fast cholesky factorization + pcg"
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)

    # get the rhs and bdy data
    rhs = [eqn.rhs(X_domain[:,i]) for i in 1:N_domain]
    bdy = [eqn.bdy(X_boundary[:,i]) for i in 1:N_boundary]

    sol_now = sol_init

    # form the fast Cholesky part that can be used to compute mtx-vct mul for Theta_test
    d = 2
    Δδ_coefs = -1.0
    δ_coefs = 0.0
    meas_δ = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_δ_int = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]
    meas_Δδ = [KoLesky.ΔδPointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), Δδ_coefs, δ_coefs) for i = 1:N_domain]
    measurements = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,3)
    measurements[1] = meas_δ; measurements[2] = meas_δ_int
    measurements[3] = meas_Δδ

    @info "[Big Theta: implicit factorization] time"
    @time implicit_bigΘ = KoLesky.ImplicitKLFactorization_FollowDiracs(cov, measurements, ρ_big, k_neighbors; lambda = lambda, alpha = alpha)
    # @time implicit_bigΘ = KoLesky.ImplicitKLFactorization(cov, measurements, ρ_big, k_neighbors; lambda = lambda, alpha = alpha)
    # @time implicit_bigΘ = KoLesky.ImplicitKLFactorization_DiracsFirstThenUnifScale(cov, measurements, ρ_big, k_neighbors; lambda = lambda, alpha = alpha)

    @info "[Big Theta: explicit factorization] time"
    @time explicit_bigΘ = KoLesky.ExplicitKLFactorization(implicit_bigΘ; nugget = nugget)
    U_bigΘ = explicit_bigΘ.U
    L_bigΘ = sparse(U_bigΘ')
    P_bigΘ = explicit_bigΘ.P

    Θtrain = approx_Theta_train(P_bigΘ, U_bigΘ, L_bigΘ,zeros(N_domain),N_boundary,N_domain)

    implicit_factor = nothing
    for step in 1:GNsteps

        @info "[Current GN step] $step"
        # get Cholesky of Theta_train
        Δδ_coefs = -1.0
        δ_coefs_int = eqn.α*eqn.m*sol_now.^(eqn.m-1)

        meas_Δδ = [KoLesky.ΔδPointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), Δδ_coefs, δ_coefs_int[i]) for i = 1:N_domain]
        measurements = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,2)
        measurements[1] = meas_δ; measurements[2] = meas_Δδ

        @info "[Theta Train: implicit factorization] time"
        @time if implicit_factor === nothing
            implicit_factor = KoLesky.ImplicitKLFactorization(cov, measurements, ρ_small, k_neighbors; lambda = lambda, alpha = alpha)
        else
            implicit_factor.supernodes.measurements .= reduce(vcat, collect.(measurements))[implicit_factor.P]
        end
        @info "[Theta Train: explicit factorization] time"
        @time explicit_factor = KoLesky.ExplicitKLFactorization(implicit_factor; nugget = nugget)

        U = explicit_factor.U
        L = sparse(U')
        P = explicit_factor.P
        rhs_now = vcat(bdy, rhs.+eqn.α*(eqn.m-1)*sol_now.^eqn.m)

        # use the approximate solution as the initial point for the pCG iteration
        Θinv_rhs = U*(L*rhs_now[P])
        Θinv_rhs[P] = Θinv_rhs

        # pcg step for Theta_train\rhs
        Θtrain.δ_coefs .= δ_coefs_int
        precond = precond_Theta_train(P,U,L)

        @info "[pcg started]"
        @time Θinv_rhs, ch = cg!(Θinv_rhs, Θtrain, rhs_now; Pl = precond, log=true)
        @info "[pcg finished], $ch"

        # get approximated sol_now = Theta_test * Θinv_rhs
        tmp = vcat(Θinv_rhs[1:N_boundary], δ_coefs_int .* Θinv_rhs[N_boundary+1:end], Θinv_rhs[N_boundary+1:end])
        tmp[P_bigΘ] = L_bigΘ\(U_bigΘ\tmp[P_bigΘ])
        @views sol_now = tmp[N_boundary+1:N_domain+N_boundary]
    end
    @info "[solver finished]"
    return sol_now
end

function main(args)
    α = args.alpha::Float64
    m = args.m::Int
    Ω = [[0,1] [0,1]]
    # ground truth solution
    freq = 600
    s = 6
    function fun_u(x)
        ans = 0
        @inbounds for k = 1:freq
            ans += sin(pi*k*x[1])*sin(pi*k*x[2])/k^s
            # H^t norm squared is sum 1/k^{2s-2t}, so in H^{s-1/2}
        end
        return ans
    end

    # right hand side
    function fun_rhs(x)
        ans = 0
        @inbounds for k = 1:freq
            ans += (2*k^2*pi^2)*sin(pi*k*x[1])*sin(pi*k*x[2])/k^s
        end
        return ans + α*fun_u(x)^m
    end

    # boundary value
    function fun_bdy(x)
        return fun_u(x)
    end

    @info "[solver started] NonlinElliptic2d"
    @info "[equation] -Δu + $α u^$m = f"
    eqn = NonlinElliptic2d(α,m,Ω,fun_bdy,fun_rhs)

    h_in = args.h::Float64; h_bd = h_in
    X_domain, X_boundary = sample_points_grid(eqn, h_in, h_bd)
    # X_domain, X_boundary = sample_points_rdm(eqn, 10000, 400)
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    @info "[sample points] grid size $h_in"
    @info "[sample points] N_domain is $N_domain, N_boundary is $N_boundary"

    lengthscale = args.sigma
    if args.kernel == "Matern5half"
        cov = KoLesky.MaternCovariance5_2(lengthscale)
    elseif args.kernel == "Matern7half"
        cov = KoLesky.MaternCovariance7_2(lengthscale)
    elseif args.kernel == "Matern9half"
        cov = KoLesky.MaternCovariance9_2(lengthscale)
    elseif args.kernel == "Matern11half"
        cov = KoLesky.MaternCovariance11_2(lengthscale)
    elseif args.kernel == "Gaussian"
        cov = KoLesky.GaussianCovariance(lengthscale)
    end
    @info "[kernel] choose $(args.kernel), lengthscale $lengthscale\n"

    nugget = args.nugget::Float64
    @info "[nugget] $nugget"

    GNsteps_approximate = args.GNsteps::Int
    @info "[total GN steps] $GNsteps_approximate"

    ρ_big = args.rho_big::Float64
    ρ_small = args.rho_small::Float64
    k_neighbors = args.k_neighbors::Int
    @info "[Fast Cholesky] ρ_big = $ρ_big, ρ_small = $ρ_small, k_neighbors = $k_neighbors"

    sol_init = zeros(N_domain) # initial solution
    truth = [fun_u(X_domain[:,i]) for i in 1:N_domain]

    fast_solve() = iterGPR_fast_pcg(eqn, cov, X_domain, X_boundary, sol_init, nugget, GNsteps_approximate; ρ_big = ρ_big, ρ_small = ρ_small, k_neighbors=k_neighbors);
    sol = fast_solve()

    # time the second time, to avoid compilation time count
    time = @elapsed fast_solve()

    pts_accuracy = sqrt(sum((truth-sol).^2)/N_domain)
    @info "[L2 accuracy: pCG method] $pts_accuracy"
    pts_max_accuracy = maximum(abs.(truth-sol))
    @info "[Linf accuracy: pCG method] $pts_max_accuracy"

    # @btime @time iterGPR_fast_pcg($eqn, $cov, $X_domain, $X_boundary, $sol_init, $nugget, $GNsteps_approximate; ρ_big = $ρ_big, ρ_small = $ρ_small, k_neighbors=$k_neighbors);

    if args.compare_exact
        GNsteps_exact = 2
        @info "[comparison: exact method]"
        @time sol_exact = iterGPR_exact(eqn, cov, X_domain, X_boundary, sol_init, nugget, GNsteps_exact)
        pts_accuracy_exact = sqrt(sum((truth-sol_exact).^2)/N_domain)
        @info "[L2 accuracy: exact method] $pts_accuracy_exact"
        pts_max_accuracy_exact = maximum(abs.(truth-sol_exact))
        @info "[Linf accuracy: exact method] $pts_max_accuracy_exact"
    end

    return time, pts_accuracy, pts_max_accuracy
end

args = parse_commandline()
# args = (; (Symbol(k) => v for (k,v) in args)...) # Named tuple from dict
# main(args)
arr_kernel = ["Matern5half"]
# arr_kernel = ["Matern5half", "Matern7half", "Matern9half"]
# arr_h = [0.02,0.01]
arr_h = [0.02,0.01,0.005,0.0025]
arr_ρ = [2.0, 3.0, 4.0]

result = Dict()
for kernel in arr_kernel
    result[("kernel",kernel)] = Dict()
    args["kernel"] = kernel
    for h in arr_h
        args["h"] = h
        result[("kernel",kernel)][("h",h)] = Dict()

        for ρ in arr_ρ
            result[("kernel",kernel)][("h",h)][("rho",ρ)] = Dict()
            args["rho_big"] = ρ
            args["rho_small"] = ρ
            args_now = (; (Symbol(k) => v for (k,v) in args)...) # Named tuple from dict

            @info "-------------------------------"
            time, pts_accuracy, pts_max_accuracy = main(args_now)
            @info "-------------------------------"
            result[("kernel",kernel)][("h",h)][("rho",ρ)]["time"] = time
            result[("kernel",kernel)][("h",h)][("rho",ρ)]["L2"] = pts_accuracy
            result[("kernel",kernel)][("h",h)][("rho",ρ)]["Linf"] = pts_max_accuracy
        end
    end
end

save("NonlinElliptic2d_data.jld", "result", result)
# save results


text:
@testset "AWSException" begin
    code = "InvalidSignatureException"
    message = "Signature expired: ..."
    body = """
        {
           "__type": "$code",
           "message": "$message"
        }
        """
    headers = ["Content-Type" => "application/x-amz-json-1.1"]
    status_code = 400

    # This does not actually send a request, just creates the object to test with
    req = HTTP.Request("GET", "https://amazon.ca", headers, body)
    resp = HTTP.Response(status_code, headers; body=body, request=req)
    ex = AWSException(HTTP.StatusError(status_code, resp))

    @test ex.code == code
    @test ex.message == message
end


text:
"""

    simplex_split(;k::Int = 2, dim::Int = 3, orientations = false)

Instructions for splitting a ´dim´-dimensional simplex with a size
reducing factor `k`.
"""
function simplex_split(;k::Int = 2, dim::Int = 3, orientations = false)

    sequences::Array{Int, 2} = tensordecomp(k, dim)
    n_seq = size(sequences, 1)

    χ1 = sequences .* (dim + 1)
    χ2 = repeat(transpose(collect(1:dim)), n_seq, 1)
    χ::Array{Int, 2} = sort(χ1 .+ χ2, dims = 2)

    matrices_simplicial_subdivision::Array{Int, 2} = zeros(Int, (dim + 1) * k^dim, k)

    for a = 1:n_seq
        tmp = ones(Int, k * (dim + 1))
        for b = 1:(dim - 1)
            tmp[(χ[a, b] + 1):(χ[a, b+1])] = (b + 1) * ones(χ[a, b + 1] - χ[a, b], 1)
        end

        tmp[(χ[a, dim] + 1):size(tmp, 1)] = (dim + 1) * ones(size(tmp, 1) - χ[a, dim], 1)

        indices = ((a - 1) * (dim + 1) + 1):(a*(dim + 1))
        matrices_simplicial_subdivision[indices, :] = reshape(tmp, dim + 1, k)
    end

    simplex_orientations = zeros(Float64, k^dim, 1)
    Χ_values = zeros(Int, dim, 1)

    for i = 1:size(simplex_orientations, 1)
        indices::Vector{Int} = ((dim + 1) * (i - 1) + 1):((dim + 1) * i)
        Χi = matrices_simplicial_subdivision[indices, :]

        ΔΧi = (Χi[2:(dim + 1), :] - Χi[1:dim, :]) .* repeat(transpose(collect(1:k)), dim, 1)

        transition_ind = round.(Int, ΔΧi * ones(k, 1))

        M = zeros(dim, dim)

        for j = 1:dim
            Χ_values[j] =  Χi[j, transition_ind[j]]
            M[j, Χ_values[j]] = 1
        end

        simplex_orientations[i] = det(M) / k^dim
    end

    if orientations # Should the orientations of the simplices also be returned?
        return matrices_simplicial_subdivision, simplex_orientations
    else
        return matrices_simplicial_subdivision
    end

end

"""
    heaviside0(x)

Returns 0 if `x < 0` and 1 if `x >= 0` . If `x` is a vector, do the same, but element-wise.
"""
function heaviside0 end

heaviside0(x::Real) = ifelse(x < 0, zero(x), ifelse(x >= 0, one(x), oftype(x, 0)))
heaviside0(x::AbstractArray{T}) where {T<:Real} = heaviside0.(x)

"""
    simplex_splitting_rules(k, dim)

Generate generic rules for generating vertices when splitting a `dim`-dimensional
simplex with a size reducing factor of k.

## Returns
1. An array of dimension `r * k` (``R = (k + dim)!/(k! dim!) - dim - 1``), where each
    row contains the indices of the original vertices (to be linearly combined to
    form new vertices?); these vertices defining the new vertices of the simplicial
    subdivision. The number of columns in this array defines how many new vertices
    are generated each split.
2. An array of integers with dimension `(dim + 1)*k^dim x 1`. It consists of a
    stack of k^E columns of dimension E+1, each stack containing
    the indices of the vertices generating the corresponding simplex
    of the subdivision in the order given by
    `[rows(original_vertices); rows(new_vertices)]`
"""
function simplex_splitting_rules(k::Int, dim::Int)

    # Get simplex splitting rules, ignore orientations for now.
    # What are these instructions?
    splits = simplex_split(;k = k, dim = dim, orientations = false)

    # Retrieve integer labels from the tensor decomposition. A way of uniquely identifying
    # the rows of splits. Each integer corresponds to a row in
    # splits (but are not necessarily in order).
    integer_labels = (splits .- 1) * ((dim + 1) .^ collect(0:(k - 1)))
    uniquelabels = unique(integer_labels)
    n_rows = size(uniquelabels, 1)

    J = indexin(integer_labels, uniquelabels)

    # What does I contain?
    # IndicesOfUnqiueIntegersinOriginal
    I = zeros(Int, length(uniquelabels))
    for i = 1:length(uniquelabels)
        I[i] = findfirst(integer_labels .== uniquelabels[i])
    end

    # What is all this aux stuff? Is there a simpler, more readable way to do this?
    Aux = splits[I, :]
    V = repeat(uniquelabels, 1, dim + 1)
    aux = repeat((((dim + 1)^k - 1) ./ dim ) * transpose(collect(0:dim)), size(uniquelabels, 1), 1)

    aux = heaviside0(-abs.(V - aux)) .* repeat(collect(1:n_rows), 1, dim + 1)
    aux = aux[findall(aux .> 0)]

    aux = round.(Int, aux)
    Caux = round.(Int, setdiff(1:n_rows, aux))
    Ipp = [aux; Caux]

    Ip = zeros(Int, size(I))

    # for any i\in Ipp, Ip(i) is the number j\in {1,...,size(Ipp,1)} such that Ipp(j)=i
    Ip[Ipp] = transpose(collect(1:size(I, 1)))

    # Rules for forming the strictly new vertices of the subtriangulation.
    # The i-th row in this array contains what?
    # Indices of the non-repeated vertices that are not part of the original simplex.
    # These are the truly new vertices generated by the splitting.
    rules_new_vertices = Aux[Caux, :]

    allverts = vcat(repeat(collect(1:(dim + 1)), 1, k), rules_new_vertices)

    # Subarrays of dimension E+1 x 1 vertically concatenated. Each subarray contains the
    # indices in the rows of NewVert corresponding to the vertices furnishing the corresponding
    # simplex (analogous to the output of SimplexSplitting, check notes).
    SubTriangulation = round.(Int, Ip[J])
    SubTriangulation = reshape(SubTriangulation, dim + 1, k^dim)

    return allverts, SubTriangulation
end

export simplex_splitting_rules


text:
## Linear Algebra ##

function turing_chol(A::AbstractMatrix, check)
    chol = cholesky(A, check=check)
    (chol.factors, chol.info)
end
function turing_chol_back(A::AbstractMatrix, check)
    C, chol_pullback = rrule(cholesky, A, Val(false), check=check)
    function back(Δ)
        Ȳ = Tangent{typeof(C)}((U=Δ[1]))
        ∂C = chol_pullback(Ȳ)[2]
        (∂C, nothing)
    end
    (C.factors,C.info), back
end

function symm_turing_chol(A::AbstractMatrix, check, uplo)
    chol = cholesky(Symmetric(A, uplo), check=check)
    (chol.factors, chol.info)
end
function symm_turing_chol_back(A::AbstractMatrix, check, uplo)
    C, chol_pullback = rrule(cholesky, Symmetric(A,uplo), Val(false), check=check)
    function back(Δ)
        Ȳ = Tangent{typeof(C)}((U=Δ[1]))
        ∂C = chol_pullback(Ȳ)[2]
        (∂C, nothing, nothing)
    end
    (C.factors, C.info), back
end

# Tracker's implementation of ldiv isn't good. We'll use Zygote's instead.
zygote_ldiv(A::AbstractMatrix, B::AbstractVecOrMat) = A \ B

# fixes `randn` on GPU (https://github.com/TuringLang/DistributionsAD.jl/pull/108)
function adapt_randn(rng::AbstractRNG, x::AbstractArray, dims...)
    return adapt_randn(rng, eltype(x), x, dims...)
end
function adapt_randn(rng::AbstractRNG, ::Type{T}, x::AbstractArray, dims...) where {T}
    return adapt(parameterless_type(x), randn(rng, T, dims...))
end

# required by Adapt >= 3.3.0: https://github.com/SciML/OrdinaryDiffEq.jl/issues/1369
Base.@pure __parameterless_type(T) = Base.typename(T).wrapper
parameterless_type(x) = parameterless_type(typeof(x))
parameterless_type(x::Type) = __parameterless_type(x)

@non_differentiable adapt_randn(::Any...)


text:
# testing poisson time sampler

a, b = 1.1, 0.0
n = 5000
Λ0(a, b, T) = a*T + b*(T)^2/2
P(a, b, T) = 1 - exp(-Λ0(a, b, T))
T = 0.7

for (a, b, pt) in ((1.1, 0.0, NaN), (1.1, 0.3, NaN), (0.0, 0.3, NaN), (1.1, -0.5, NaN),
    (-0.5, 1, P(0, 1, T-0.5)), (-1, -2, 0.0))
    p = mean(poisson_time(a, b, rand()) < T for i in 1:n)
    if isnan(pt)
        pt = P(a, b, T)
    end
    @test abs(p - pt) < 2/sqrt(n)
end

# negative log-density with respect to Lebesgue
# ϕ(x) =  (x - π)^2/2 # not needed
# gradient of ϕ(x)
const σ2 = 1.3
∇ϕ(x) = (x - π/2)/σ2
∇ϕhat(x) = (x - π/2)/σ2 + 0.1(rand()-0.5)

x0, θ0 = 1.01, -1.5
T = 8000.0
out1, _ = ZigZagBoomerang.pdmp(∇ϕhat, x0, θ0, T, 2.5, ZigZag1d())

@testset "ZigZag1d" begin
    @test T/10 < length(out1) < T*10
    est = 1/T*sum((eventposition.(out1)[1:end-1] + eventposition.(out1)[2:end])/2 .* diff(eventtime.(out1)))
    @test abs(est - pi/2) < 2/sqrt(length(out1))
    dt = 0.01
    traj = ZigZagBoomerang.discretize(out1, ZigZag1d(), dt)
    @test abs(-(extrema(diff(traj.t[1:end÷3]))...)) < 1e-10
    est = mean(traj.x)
    @test abs(est - pi/2) < 2/sqrt(length(out1))
    est2 = var(traj.x)
    @test abs(est2 - σ2) < 2/sqrt(length(out1))
end

x0, θ0 = 1.41, +0.1
B = Boomerang1d(1.0)
out2, _ = ZigZagBoomerang.pdmp(∇ϕ, x0, θ0, T, 1.6, B)

@testset "Boomerang1d" begin
    @test T/10 < length(out2) < T*10
    dt = 0.01
    traj = ZigZagBoomerang.discretize(out2, B, dt)
    @test abs(-(extrema(diff(traj.t[1:end÷3]))...)) < 1e-10
    est = mean(traj.x)
    @test abs(est - pi/2) < 5/sqrt(length(out2))
    est2 = var(traj.x)
    @test abs(est2 - σ2) < 11/sqrt(length(out2))
end

B = Boomerang1d(1.1, 1.2, 0.5)
out2, _ = ZigZagBoomerang.pdmp(∇ϕhat, x0, θ0, T, 10.0, B)

@testset "Boomerang1dNoncentredSub" begin
    @test T/10 < length(out2) < T*10
    dt = 0.01
    traj = ZigZagBoomerang.discretize(out2, B, dt)
    @test abs(-(extrema(diff(traj.t[1:end÷3]))...)) < 1e-10
    est = mean(traj.x)
    @test abs(est - pi/2) < 5/sqrt(length(out2))
    est2 = var(traj.x)
    @test abs(est2 - σ2) < 5/sqrt(length(out2))
end


text:
using Colors, CSV, DataFrames, PGFPlotsX
deaths_raws = [normpath("data", "2020-11-7034-4A.csv"),
    normpath("data", "2021-3-7289-4A.csv"),
    normpath("data", "2021-11-7640-4A.csv")]
pop_raw = normpath("data", "BE0101A6_20201118-131030.txt")
pop_form = tempname()
run(pipeline(`iconv -f iso-8859-1 -t utf-8`, stdin=pop_raw,
    stdout=pipeline(`./sepop.awk`, pop_form)))

function load_deaths_frames(raws)
    dfs = []
    for f in raws
        deaths_form = tempname()
        run(pipeline(`./sedor20.awk`, stdin=f, stdout=deaths_form))
        df = CSV.File(deaths_form) |> DataFrame
        dfs = vcat(dfs, df)
    end
    dfs
end

deaths_frames = load_deaths_frames(deaths_raws)

df_dates = ["2020-01-01–2020-06-30", "2020-07-01–2020-12-31", "2021-01-01–2021-06-30"]

ages = DataFrame(age = collect(0:21),
    agest = [0;0;1;collect(5:5:95)],
    age2 =[0;1;fill(2,3);fill(3,6);fill(4,4);fill(5,2);fill(6,5)])
ages2 = combine(groupby(ages, :age2), df -> DataFrame(agest=minimum(df[:, :agest])))

pop = DataFrame()
pop1 = CSV.File(pop_form) |> DataFrame
ages1 = DataFrame(age1 = collect(0:101),
    age = [0;1;fill(2,4);map((x)->fld(x,5)+2, collect(5:94));fill(21,6)])
pop1_age = innerjoin(pop1, ages1, on = :age1)

function popsexadd(sex)
    pop1_sex = pop1_age[pop1_age[:, :sex].==sex, [:n, :age]]
    pop_sex = combine(groupby(pop1_sex, :age), df -> DataFrame(sex = sex, n = sum(df[:, :n])))
end

pop = vcat(popsexadd("Kv"), popsexadd("M"))

function regsub(regexpr, df)
    infile = tempname()
    CSV.write(infile, df)
    outfile = tempname()
    run(pipeline(`awk -v reg=$regexpr -f regsub.awk`, stdin=infile, stdout=outfile))
    CSV.File(outfile) |> DataFrame
end

function df_reg(sex, caexpr, ageaggr, deaths_frame)
    deaths = deaths_frames[deaths_frame]
    if caexpr == "pop"
        df_sub = pop[pop[:, :sex].==sex, [:age, :n]]
        sort!(df_sub, (:age))
    else
        df_sub = regsub(caexpr, deaths[deaths[:, :sex].==sex, [:icd10, :age, :n]])
    end
    df_aggr = combine(groupby(df_sub, :age), df -> DataFrame(n=sum(df[:, :n])))
    df_aggr_age = innerjoin(df_aggr, ages, on = :age)
    if ageaggr
	df_aggr_age2 = combine(groupby(df_aggr_age[:, [:n, :age2]], :age2), df -> DataFrame(n=sum(df[:, :n])))
        return innerjoin(df_aggr_age2, ages2, on = :age2)
    else
        return df_aggr_age
    end
end

function propframe(sex, ca1; ca2 = "A00-Y98", cax = "U07-U07", ageaggr = false, deaths_frame = 1)
    ca1frame = df_reg(sex, ca1, ageaggr, deaths_frame)
    ca2frame = df_reg(sex, ca2, ageaggr, deaths_frame)
    propframe = DataFrame(age = ca1frame[:, 1], agest = ca1frame[:, :agest],
        ca1n = ca1frame[:, :n],
        ca2n = ca2frame[:, :n],
        prop = ca1frame[:, :n] ./ ca2frame[:, :n],
    )
    if cax != ""
        caxframe = df_reg(sex, cax, ageaggr, deaths_frame)
        propframe[!, :caxn] = caxframe[:, :n]
        propframe[!, :xprop] = ca1frame[:, :n] ./ (ca2frame[:, :n] .- caxframe[:, :n])
    end
    return Dict(:sex => sex, :ca1 => ca1, :ca2 => ca2, :cax => cax, :frame => propframe)
end

function cas_sexes_prop_plot(cas, denomca, cax = "", deaths_frame = 1, ymode = "linear")
    if cax == ""
        ylab = "döda/$(denomca[:label])"
    else
        ylab = "döda/($(denomca[:label])−$(cax[:label]))"
    end
    sax = @pgf Axis({"yticklabel style"={"/pgf/number format/use comma"},
        ymode = ymode, legend_pos="outer north east",
        title="Döda Sverige $(df_dates[deaths_frame])",
        xlabel="ålder", ylabel=ylab, xmajorgrids, ymajorgrids, yminorgrids, "minor tick num"=1})
    plotcolors = distinguishable_colors(length(cas)*2 + 1, [RGB(1,1,1)])[2:end]
    for (i, sex) in enumerate(["Kv", "M"])
        for (j, ca) in enumerate(cas)
            ca_sex_prop_plot(sax, ca, denomca, cax, sex, deaths_frame, plotcolors[i*j])
        end
    end
    sax
end

function ca_sex_prop_plot(sax, ca, denomca, cax, sex, deaths_frame, col)
    if cax == ""
        caxpr = ""
        plotcol = :prop
    else
        caxpr = cax[:caexpr]
        plotcol = :xprop
    end
    df = propframe(sex, ca[:caexpr], ca2 = denomca[:caexpr], cax = caxpr,
                   deaths_frame = deaths_frame)[:frame]
    @pgf push!(sax, PlotInc({color = col},
        Table([df[:, :agest][2:end], df[:, plotcol][2:end]])),
        LegendEntry("$(ca[:label]) $(lowercase(sex))"))
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
# Simulate a Poisson point process on a triangle
# Author: H. Paul Keeler, 2019.
# Website: hpaulkeeler.com
# Repository: github.com/hpaulkeeler/posts
# For more details, see the post:
# hpaulkeeler.com/testing-the-julia-language-with-point-process-simulations/

#Note: Need the .+ for adding a scalar to an array
#Also need . for sqrt, exp, cos, sin etc and assinging scalars to arrays
#Don't confuse 1.-V with 1 .-V for array V
#Best to use vectors instead of 1-D matrices eg x=rand(n),  NOT x=rand(n,1).

using Distributions #for random simulations
using Plots #for plotting

#Simulation window parameters -- points A,B,C of a triangle
xA=0;xB=1;xC=1; #x values of three points
yA=0;yB=0;yC=1; #y values of three points;

#Point process parameters
lambda=100; #ntensity (ie mean density) of the Poisson process

#calculate sides of trinagle
a=sqrt((xA-xB)^2+(yA-yB)^2);
b=sqrt((xB-xC)^2+(yB-yC)^2);
c=sqrt((xC-xA)^2+(yC-yA)^2);
s=(a+b+c)/2; #calculate semi-perimeter

#Use Herron's forumula to calculate area
areaTotal=sqrt(s*(s-a)*(s-b)*(s-c)); #area of triangle

#Simulate a Poisson point process
numbPoints=rand(Poisson(areaTotal*lambda)); #Poisson number of points
U=rand(numbPoints);#uniform random variables
V=rand(numbPoints);#uniform random variables

xx=sqrt.(U).*xA .+sqrt.(U).*(1 .-V).*xB .+sqrt.(U).*V.*xC;#x coordinates of points
yy=sqrt.(U).*yA .+sqrt.(U).*(1 .-V).*yB .+sqrt.(U).*V.*yC;#y coordinates of points

#Plotting
plot1=scatter(xx,yy,xlabel ="x",ylabel ="y", leg=false);
display(plot1);

max_stars_repo_name: hpkeeler/posts

max_issues_repo_name: hpkeeler/posts

max_forks_repo_name: hpkeeler/posts


text:
"""
    tags_match_re_im(old_tag, new_tag)

If `old_tag` and `new_tag` tag the real part and imaginary part of the same
quantity, such as `old_tag = "test/struct/re"` and `new_tag = "test/struct/im"`
then returns true. False otherwise
"""
function tags_match_re_im(old_tag, new_tag)
    if old_tag[end-2:end] == "/re"
        old_tag_parts = split(old_tag, "/")
        new_tag_parts = split(new_tag, "/")
        if ( all(old_tag_parts[1:end-1] .== new_tag_parts[1:end-1]) &&
            new_tag_parts[end] == "im" )
            return true
        end
    end
    return false
end

"""
    lookahead_deserialize(old_tag, old_val, summaries, state, type)

Looks ahead at state, and checks that the tag of `state` matches `old_tag`, and
if so, attempts to combine the value of `state` with that of `old_val`.

For example, if `old_tag = "someval/re"` then if `state` contains "someval/im" the
two values will be combined as real and immaginary part.

Returns a possibly modified tag, a possibily modified value, and the number
of states that it has looked ahead.
If the lookahead fails to find matching states, he returns (old_tag, old_val, 0)
"""
function lookahead_deserialize(old_tag, old_val, evs::Summary, state, typ)
    res = old_tag, old_val, state

    if typ == :histo
        res = lookahead_deserialize_histogram_summary(old_tag, old_val, evs, state)
    elseif typ === :image
        res = lookahead_deserialize_image_summary(old_tag, old_val, evs, state)
    elseif typ == :simple_value
        res = lookahead_deserialize_simple_value_summary(old_tag, old_val, evs, state)
    end

    return res
end

# Lookahead for simple values.
# Checks the next entry if it is /re and /im, and combines them
function lookahead_deserialize_simple_value_summary(old_tag, old_val, evs::Summary,
                                                 state)
    # prepare the default output (when no action is taken)
    result = old_tag, old_val, state
    # iterate to the next element
    res = iterate(evs, state + 1)
    res == nothing && return result

    # if the next event is identified, check its type
    (new_tag, summary), i_state = res
    typ = summary_type(summary)

    # if types match, check tags. If they match, return modified structure
    if typ === :simple_value && tags_match_re_im(old_tag, new_tag)
        val_im = summary.simple_value
        result = new_tag[1:end-3], old_val + im*val_im, i_state - 1
    else
        return result
    end
end


text:
using Plots
using Reproduce
using JLD2
using FileIO
using ProgressMeter

previous_GKSwstype = get(ENV, "GKSwstype", "")
ENV["GKSwstype"] = "100"

p = plot()
folder = "experiment_data/RP_0_0x5a78768ca53838bb/"
save_file = "plotting/chunlok/generated_plots/state_vistation_heatmap.gif"

results_file = folder * "results.jld2"
@load results_file results
# println(results)

settings_file = folder * "settings.jld2"
settings = FileIO.load(settings_file)["parsed_args"]

state_visitation_data = results[:oned_tmaze_state_visitation]

num_log = size(state_visitation_data)[3]
log_interval = 100
num_frames = 50
frames_between = floor(Int32, num_log / num_frames)

prog = Progress(num_frames, 1)
anim = @animate for i ∈ 1:num_frames
    next!(prog)
    init_frame = i * frames_between
    last_frame = min((i+1) * frames_between, num_log)
    frame_info = state_visitation_data[:, :, init_frame:last_frame]

    frame_info = dropdims(sum(frame_info, dims=3), dims=3)
    # println(frame_info)
    # display("text/plain", frame_info)
    # asdasd

    # println(sum(frame_info))
    frame_info = frame_info ./ sum(frame_info)

    # println(size(frame_info))

    frame_info = permutedims(frame_info, [2, 1])
    cgradient = cgrad(:rainbow, scale=:exp)
    heatmap!(frame_info, title="Step: $(i * frames_between * log_interval)", c=cgradient, clim=(0, 1), zscale=:log10)
end

gif(anim, save_file, fps = 2)
ProgressMeter.finish!(prog)

# gif(anim, save_file, fps=15)
# savefig(save_file)

ENV["GKSwstype"] = previous_GKSwstype

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using PlayGios
using Test

@testset "PlayGios.jl" begin
    # Write your tests here.
    @test true
end


text:
abstract type RegressionType end

struct OLS
    robust::Bool

    function OLS(robust=true)
        new(robust)
    end
end


text:
# # Linear regression

using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

# Import the dataset.

# ### snippet 5.1

df = DataFrame(CSV.read(joinpath("..", "..", "data", "WaffleDivorce.csv"), delim=';'))
mean_ma = mean(df[:, :MedianAgeMarriage])
df[!, :MedianAgeMarriage_s] = convert(Vector{Float64},
  (df[:, :MedianAgeMarriage]) .- mean_ma)/std(df[:, :MedianAgeMarriage]);

# Model ``y ∼ Normal(y - Xβ, σ)``. Flat prior `β`, half-T for `σ`.

Base.@kwdef mutable struct WaffleDivorce{Ty <: AbstractVector,
  Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
end

# Write a function to return a properly dimensioned transformation.

function make_transformation(model::WaffleDivorce)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

# Instantiate the model with data and inits.

x = hcat(ones(size(df, 1)), df[:, :MedianAgeMarriage_s]);
model = WaffleDivorce(;y=df[:, :Divorce], x=x);

# Make tmodel callable with the parameters *as a single argument*.

function (model::WaffleDivorce)(θ)
    @unpack y, x = model   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(10, 10), x[1]) # alpha
    ll += logpdf(Normal(0, 1), x[2]) # beta
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- x*β)
    ll
end

println()
model((β = [1.0, 2.0], σ = 1.0)) |> display
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display
println()

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display
println()

a3d = Array{Float64, 3}(undef, 1000, 3, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:2, j] = values(posterior[i].β)
    a3d[i, 3, j] = values(posterior[i].σ)
  end
end

pnames = ["β[1]", "β[2]", "σ"]
sections =   Dict(
  :parameters => pnames,
)
chns = create_mcmcchains(a3d, pnames, sections, start=1);
chns = set_names(chns, Dict("β[1]" => "α", "β[2]" => "β"))

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
         Mean        SD       Naive SE       MCSE      ESS
    a  9.6882466 0.22179190 0.0035068378 0.0031243061 1000
   bA -1.0361742 0.21650514 0.0034232469 0.0034433245 1000
sigma  1.5180337 0.15992781 0.0025286807 0.0026279593 1000

Quantiles:
         2.5%      25.0%     50.0%      75.0%       97.5%
    a  9.253141  9.5393175  9.689585  9.84221500 10.11121000
   bA -1.454571 -1.1821025 -1.033065 -0.89366925 -0.61711705
sigma  1.241496  1.4079225  1.504790  1.61630750  1.86642750
";

# Describe chains

describe(chns)

# end of m4.5d.jl


text:
lenn_file   = string(path, "/SampleFiles/ASCII/0215162000.c00")
geocsv_file = string(path, "/SampleFiles/ASCII/geo-tspair.csv")
slist_file  = string(path, "/SampleFiles/ASCII/2m-62.5hz.slist")
lenn_pat    = string(path, "/SampleFiles/ASCII/021516*c00")
geocsv_pat  = string(path, "/SampleFiles/ASCII/geo-tspair.*")
slist_pat   = string(path, "/SampleFiles/ASCII/*.slist")
geoslist_f  = string(path, "/SampleFiles/ASCII/geo-slist.csv")

printstyled("  Lennartz ASCII\n", color=:light_green)
C = verified_read_data("lennartz", lenn_file, vl=true)[1]
@test ≈(C.fs, 62.5)

printstyled("    wildcard support\n", color=:light_green)
S = verified_read_data("lennartz", lenn_pat)

printstyled("  GeoCSV timeseries\n", color=:light_green)
S = verified_read_data("geocsv", geocsv_file, vl=true)
@test S.n == 8
if Sys.iswindows()
  println("    IDs = ", join(S.id, ","))
end
i = findid("CC.JRO..BHZ", S.id)
if i > 0
  @test ≈(S.fs[i], 50.0)
  @test ==(S.loc[i], GeoLoc(lat=46.275269, lon=-122.218262, el=1219.0, inc=180.0))
end

printstyled("  GeoCSV slist\n", color=:light_green)
S = verified_read_data("geocsv.slist", geoslist_f)

printstyled("    wildcard support\n", color=:light_green)
S = verified_read_data("geocsv", geocsv_pat)

printstyled("  slist\n", color=:light_green)
S = verified_read_data("slist", slist_file, vl=true)
nx = length(S.x[1])
@test S.id[1] == "YY.ERTA..EHZ"
@test ≈(S.fs[1], 62.5)
@test isapprox(C.x[1:nx], S.x[1])

printstyled("    empty channel + check_for_gap!\n", color=:light_green)
S.t[1] = Array{Int64}(undef, 0, 2)
S.x[1] = Float32[]
verified_read_data!(S, "slist", slist_file)
@test S.id[1] == "YY.ERTA..EHZ"
@test ≈(S.fs[1], 62.5)
@test C.t[1,2] == S.t[1][1,2]
@test isapprox(C.x[1:nx], S.x[1])

printstyled("    wildcard support\n", color=:light_green)
S = verified_read_data("slist", slist_pat)
@test S.id[1] == "YY.ERTA..EHZ"
@test ≈(S.fs[1], 62.5)
@test isapprox(C.x[1:nx], S.x[1])

printstyled("  channel continuation (Issue 34)\n", color=:light_green)
printstyled("    GeoCSV.tspair\n", color=:light_green)
test_chan_ext(geocsv_file, "geocsv.tspair", "CC.JRO..BHZ", 50.0, 1, 1554777720010000)

printstyled("    GeoCSV.slist\n", color=:light_green)
test_chan_ext(geoslist_f, "geocsv.slist", "IU.ANMO.00.LHZ", 1.0, 3, 1551249000000000)

printstyled("    slist\n", color=:light_green)
test_chan_ext(slist_file, "slist", "YY.ERTA..EHZ", 62.5, 1, 1013790000000000)

printstyled("    lennartz\n", color=:light_green)
test_chan_ext(lenn_file, "lennartz", ".ERTA..c00", 62.5, 1, 1013790000000000)


text:
"""
    AbstractCubedSphere <: AbstractMesh2D

This is an abstract type of cubed-sphere meshes on `SphereDomain`s. A
cubed-sphere mesh has 6 panels, laid out as follows:
```
                                          :   Panel 1   :
                            +-------------+-------------+
                            |     +x1     |     +x1     |
                            |             |             |
                            |    Panel    |    Panel    |
                            |+x3   5   -x3|-x2   6   +x2|
                            |     -x2     |     -x3     |
                            |             |             |
                            |     -x1     |     -x1     |
              +-------------+-------------+-------------+
              |     -x2     |     -x2     |
              |             |             |
              |    Panel    |    Panel    |
              |+x1   3   -x1|+x3   4   -x3|
              |     +x3     |     -x1     |
              |             |             |
              |     +x2     |     +x2     |
+-------------+-------------+-------------+
|     +x3     |     +x3     |
|             |             |
|    Panel    |    Panel    |
|-x2   1   +x2|+x1   2   -x1|
|     +x1     |     +x2     |
|             |             |
|     -x3     |     -x3     |
+-------------+-------------+
:   Panel 6   :
```

This is the same panel ordering used by the S2 Geometry library (though we use 1-based
instead of 0-based numering).

Elements are indexed by a `CartesianIndex{3}` object, where the components are:
- horizontal element index (left to right) within each panel.
- vertical element index (bottom to top) within each panel.
- panel number

Subtypes should have the following fields:
- `domain`: a `SphereDomain`
- `ne`: number of elements across each panel

# External links
- [S2Geometry library](https://s2geometry.io/devguide/s2cell_hierarchy)
- [MIT GCM exch2](https://mitgcm.readthedocs.io/en/latest/phys_pkgs/exch2.html?highlight=cube%20sphere#fig-48tile)
"""
abstract type AbstractCubedSphere <: AbstractMesh2D end

Base.show(io::IO, mesh::AbstractCubedSphere) = print(
    io,
    mesh.ne,
    "×",
    mesh.ne,
    "×",
    6,
    "-element ",
    nameof(typeof(mesh)),
    " of ",
    mesh.domain,
)

domain(mesh::AbstractCubedSphere) = mesh.domain
elements(mesh::AbstractCubedSphere) = CartesianIndices((mesh.ne, mesh.ne, 6))

is_boundary_face(mesh::AbstractCubedSphere, elem, face) = false
boundary_face_name(mesh::AbstractCubedSphere, elem, face) = nothing

function opposing_face(
    mesh::AbstractCubedSphere,
    element::CartesianIndex{3},
    face::Int,
)
    ne = mesh.ne
    x, y, panel = element.I
    if face == 1
        if y > 1
            (CartesianIndex(x, y - 1, panel), 3, true)
        elseif isodd(panel)
            (CartesianIndex(x, ne, mod1(panel - 1, 6)), 3, true)
        else
            (CartesianIndex(ne, ne - x + 1, mod1(panel - 2, 6)), 2, true)
        end
    elseif face == 2
        if x < ne
            (CartesianIndex(x + 1, y, panel), 4, true)
        elseif isodd(panel)
            (CartesianIndex(1, y, mod1(panel + 1, 6)), 4, true)
        else
            (CartesianIndex(ne - y + 1, 1, mod1(panel + 2, 6)), 1, true)
        end
    elseif face == 3
        if y < ne
            (CartesianIndex(x, y + 1, panel), 1, true)
        elseif isodd(panel)
            (CartesianIndex(1, ne - x + 1, mod1(panel + 2, 6)), 4, true)
        else
            (CartesianIndex(x, 1, mod1(panel + 1, 6)), 1, true)
        end
    elseif face == 4
        if x > 1
            (CartesianIndex(x - 1, y, panel), 2, true)
        elseif isodd(panel)
            (CartesianIndex(ne - y + 1, ne, mod1(panel - 2, 6)), 3, true)
        else
            (CartesianIndex(ne, y, mod1(panel - 1, 6)), 2, true)
        end
    end
end

"""
    panel_to_coordinates(panel, coord1::Geometry.Cartesian123Point)

Given a point at `coord1` on panel 1 of a sphere, transform
it to panel `panel`.
"""
function panel_to_coordinates(panel::Integer, coord::Geometry.Cartesian123Point)
    ζ0, ζx, ζy = Geometry.components(coord)
    if panel == 1
        return Geometry.Cartesian123Point(ζ0, ζx, ζy)
    elseif panel == 2
        return Geometry.Cartesian123Point(-ζx, ζ0, ζy)
    elseif panel == 3
        return Geometry.Cartesian123Point(-ζx, -ζy, ζ0)
    elseif panel == 4
        return Geometry.Cartesian123Point(-ζ0, -ζy, -ζx)
    elseif panel == 5
        return Geometry.Cartesian123Point(ζy, -ζ0, -ζx)
    elseif panel == 6
        return Geometry.Cartesian123Point(ζy, ζx, -ζ0)
    end
    error("invalid panel $panel")
end

"""
    panel, coord1 = coordinates_to_panel(coord::Geometry.Cartesian123Point)

Given a point `coord`, return its panel number (`panel`), and its coordinates in panel 1
(`coord1`).
"""
function coordinates_to_panel(coord::Geometry.Cartesian123Point)
    maxdim = argmax(abs.(Geometry.components(coord)))
    if maxdim == 1 && coord.x1 > 0
        return 1, Geometry.Cartesian123Point(coord.x1, coord.x2, coord.x3)
    elseif maxdim == 2 && coord.x2 > 0
        return 2, Geometry.Cartesian123Point(coord.x2, -coord.x1, coord.x3)
    elseif maxdim == 3 && coord.x3 > 0
        return 3, Geometry.Cartesian123Point(coord.x3, -coord.x1, -coord.x2)
    elseif maxdim == 1 && coord.x1 < 0
        return 4, Geometry.Cartesian123Point(-coord.x1, -coord.x3, -coord.x2)
    elseif maxdim == 2 && coord.x2 < 0
        return 5, Geometry.Cartesian123Point(-coord.x2, -coord.x3, coord.x1)
    elseif maxdim == 3 && coord.x3 < 0
        return 6, Geometry.Cartesian123Point(-coord.x3, coord.x2, coord.x1)
    end
    error("invalid coordinates $x")
end

"""
    EquiangularCubedSphere <: AbstractCubedSphere

An equiangular gnomonic mesh proposed by [Ronchi1996](@cite).
Uses the element indexing convention of [`AbstractCubedSphere`](@ref).

# Constructors

    EquiangularCubedSphere(domain::Domains.SphereDomain, ne::Integer)

Constuct an `EquiangularCubedSphere` on `domain` with `ne` elements across
each panel.
"""
struct EquiangularCubedSphere{S <: SphereDomain} <: AbstractCubedSphere
    domain::S
    ne::Int
end

# not yet provided by Base Julia
# https://github.com/JuliaLang/julia/issues/28943
# this appears to be more accurate than tan(pi * x)
tanpi(x) = sinpi(x) / cospi(x)

function coordinates(mesh::EquiangularCubedSphere, elem, (ξ1, ξ2)::NTuple{2})
    radius = mesh.domain.radius
    ne = mesh.ne
    x, y, panel = elem.I
    ξx = (2 * x - ne - 1 + ξ1) / ne
    ξy = (2 * y - ne - 1 + ξ2) / ne
    ux = tanpi(ξx / 4)
    uy = tanpi(ξy / 4)
    ζ0 = radius / hypot(ux, uy, 1)
    ζx = ux * ζ0
    ζy = uy * ζ0
    panel_to_coordinates(panel, Geometry.Cartesian123Point(ζ0, ζx, ζy))
end

function containing_element(
    mesh::EquiangularCubedSphere,
    coord::Geometry.Cartesian123Point,
)
    ne = mesh.ne
    panel, coord1 = coordinates_to_panel(coord)
    ζ0, ζx, ζy = Geometry.components(coord1)
    ux = ζx / ζ0
    uy = ζy / ζ0
    ξx = 4 * atan(ux) / pi
    ξy = 4 * atan(uy) / pi

    x, ξ1 = split_refcoord(ξx, ne)
    y, ξ2 = split_refcoord(ξy, ne)

    return CartesianIndex(x, y, panel), (ξ1, ξ2)
end

"""
    EquidistantCubedSphere <: AbstractCubedSphere

An equidistant gnomonic mesh outlined in [Rancic1996](@cite) and [Nair2005](@cite).
Uses the element indexing convention of [`AbstractCubedSphere`](@ref).

# Constructors

    EquidistantCubedSphere(domain::Domains.SphereDomain, ne::Integer)

Constuct an `EquidistantCubedSphere` on `domain` with `ne` elements across
each panel.
"""
struct EquidistantCubedSphere{S <: SphereDomain} <: AbstractCubedSphere
    domain::S
    ne::Int
end
function coordinates(mesh::EquidistantCubedSphere, elem, (ξ1, ξ2)::NTuple{2})
    radius = mesh.domain.radius
    ne = mesh.ne
    x, y, panel = elem.I
    ξx = (2 * x - ne - 1 + ξ1) / ne
    ξy = (2 * y - ne - 1 + ξ2) / ne
    ζ0 = radius / hypot(ξx, ξy, 1)
    ζx = ξx * ζ0
    ζy = ξy * ζ0
    panel_to_coordinates(panel, Geometry.Cartesian123Point(ζ0, ζx, ζy))
end
function containing_element(
    mesh::EquidistantCubedSphere,
    coord::Geometry.Cartesian123Point,
)
    ne = mesh.ne
    panel, coord1 = coordinates_to_panel(coord)
    ζ0, ζx, ζy = Geometry.components(coord1)
    ξx = ζx / ζ0
    ξy = ζy / ζ0

    x, ξ1 = split_refcoord(ξx, ne)
    y, ξ2 = split_refcoord(ξy, ne)

    return CartesianIndex(x, y, panel), (ξ1, ξ2)
end

"""
    ConformalCubedSphere <: AbstractCubedSphere

A conformal mesh outlined in [Rancic1996](@cite).
Uses the element indexing convention of [`AbstractCubedSphere`](@ref).

# Constructors

    ConformalCubedSphere(domain::Domains.SphereDomain, ne::Integer)

Constuct an `ConformalCubedSphere` on `domain` with `ne` elements across
each panel.
"""
struct ConformalCubedSphere{S <: SphereDomain} <: AbstractCubedSphere
    domain::S
    ne::Int
end
function coordinates(mesh::ConformalCubedSphere, elem, (ξ1, ξ2)::NTuple{2})
    radius = mesh.domain.radius
    ne = mesh.ne
    x, y, panel = elem.I
    ξx = (2 * x - ne - 1 + ξ1) / ne
    ξy = (2 * y - ne - 1 + ξ2) / ne
    ζx, ζy, ζ0 = CubedSphere.conformal_cubed_sphere_mapping(ξx, ξy)
    panel_to_coordinates(
        panel,
        Geometry.Cartesian123Point(radius * ζ0, radius * ζx, radius * ζy),
    )
end
function containing_element(
    mesh::ConformalCubedSphere,
    coord::Geometry.Cartesian123Point,
)
    ne = mesh.ne
    panel, coord1 = coordinates_to_panel(coord)
    ζ0, ζx, ζy = Geometry.components(coord1)
    R = hypot(ζx, ζy, ζ0)
    ξx, ξy = CubedSphere.conformal_cubed_sphere_inverse_mapping(
        abs(ζx) / R,
        abs(ζy) / R,
        ζ0 / R,
    )

    x, ξ1 = split_refcoord(copysign(ξx, ζx), ne)
    y, ξ2 = split_refcoord(copysign(ξy, ζy), ne)

    return CartesianIndex(x, y, panel), (ξ1, ξ2)
end


text:
abstract type Encoding end
struct Binary <: Encoding end
struct JSON <: Encoding end

# types for parsing schemas
abstract type SchemaType end

StructTypes.StructType(::Type{SchemaType}) = StructTypes.AbstractType()
StructTypes.omitempties(::Type{<:SchemaType}) = true
StructTypes.StructType(::Type{<:SchemaType}) = StructTypes.Struct()

StructTypes.subtypes(::Type{SchemaType}) = (
    null=NullType,
    boolean=BooleanType,
    int=IntType,
    long=LongType,
    float=FloatType,
    double=DoubleType,
    bytes=BytesType,
    string=StringType,
    record=RecordType,
    enum=EnumType,
    array=ArrayType,
    map=MapType,
    fixed=FixedType
)

abstract type LogicalType end
StructTypes.StructType(::Type{LogicalType}) = StructTypes.AbstractType()
StructTypes.subtypekey(::Type{LogicalType}) = :logicalType
StructTypes.StructType(::Type{<:LogicalType}) = StructTypes.Struct()

StructTypes.subtypes(::Type{LogicalType}) = NamedTuple{(
    :decimal,
    :uuid,
    :date,
    Symbol("time-millis"),
    Symbol("time-micros"),
    Symbol("timestamp-millis"),
    Symbol("timestamp-micros"),
    Symbol("local-timestamp-millis"),
    Symbol("local-timestamp-micros"),
    :duration
)}((
    DecimalType,
    UUIDType,
    DateType,
    TimeMillisType,
    TimeMicrosType,
    TimestampMillisType,
    TimestampMicrosType,
    LocalTimestampMillisType,
    LocalTimestampMicrosType,
    DurationType
))

const UnionType = Vector{Union{String, SchemaType, LogicalType}}
const Schema = Union{String, LogicalType, SchemaType, UnionType}

StructTypes.StructType(::Type{Schema}) = StructTypes.Struct()

abstract type PrimitiveType <: SchemaType end

for (sym, JT) in (("Null", Missing), ("Boolean", Bool),
                  ("Int", Int32), ("Long", Int64),
                  ("Float", Float32), ("Double", Float64),
                  ("Bytes", Vector{UInt8}), ("String", String))
    T = Symbol(sym, :Type)
    lower = lowercase(sym)
    @eval begin
        struct $T <: PrimitiveType
            type::String
        end
        const $(Symbol(lower)) = $T($lower)
        juliatype(::$T) = $JT
    end
end

==(x::PrimitiveType, y::String) = x.type == y
==(x::String, y::PrimitiveType) = x == y.type

juliatype(s::String) = s == "null" ? Missing :
                       s == "boolean" ? Bool :
                       s == "int" ? Int32 :
                       s == "long" ? Int64 :
                       s == "float" ? Float32 :
                       s == "double" ? Float64 :
                       s == "bytes" ? Vector{UInt8} : String

PrimitiveType(x::String) = schematype(juliatype(x))

function schematype end

schematype(sch::Union{SchemaType, UnionType, LogicalType}) = sch
schematype(::Type{T}) where {T} = schematype(StructTypes.StructType(T), T)

schematype(::StructTypes.NullType, T) = null
schematype(::StructTypes.BoolType, T) = boolean
schematype(S::StructTypes.NumberType, T) = schematype(S, StructTypes.numbertype(T))
schematype(::StructTypes.NumberType, U::Union) = schematype(StructTypes.Struct(), U)
schematype(::StructTypes.NumberType, ::Type{T}) where {T <: Integer} = int
schematype(::StructTypes.NumberType, ::Type{T}) where {T <: Union{Int64, UInt64}} = long
schematype(::StructTypes.NumberType, ::Type{Float32}) = float
schematype(::StructTypes.NumberType, ::Type{Float64}) = double
schematype(::StructTypes.ArrayType, ::Type{<:AbstractVector{UInt8}}) = bytes
schematype(::StructTypes.StringType, T) = string

abstract type NamedType <: SchemaType end

function fullname(x::NamedType, enclosing_namespace="")
    # If the name specified contains a dot, then it is assumed to be a fullname,
    # and any namespace also specified is ignored
    any(==(UInt8('.')), codeunits(x.name)) && return x.name
    # if namespace is empty, use enclosing (if any)
    isempty(x.namespace) && return Base.string(enclosing_namespace, x.name)
    return Base.string(x.namespace, x.name)
end

include("types/binary.jl")
include("types/records.jl")
include("types/enums.jl")
include("types/arrays.jl")
include("types/maps.jl")
include("types/fixed.jl")
include("types/unions.jl")
include("types/logical.jl")
include("types/rows.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ 1f74c496-6726-11eb-3894-119ff752210e
begin
	using SymPy;
	using Test;
	using LinearAlgebra;
	using Plots;
	using Calculus;
	using PlutoUI;

	include("../../lib/KTBC.jl");
	using .KTBC: CreateKTBC;
end

# ╔═╡ 0668a58e-6727-11eb-2b6c-41626284bd0d
md"""
## Init deps
"""

# ╔═╡ 1c3315a0-6774-11eb-0081-f328cca176ed
begin
	import Pkg
	Pkg.add("PlutoUI")
	Pkg.add("SymPy")
	Pkg.add("Test")
	Pkg.add("LinearAlgebra")
	Pkg.add("Plots")
	Pkg.add("Calculus")
end

# ╔═╡ 9020f922-6725-11eb-2bc5-c190c647a78a
md"""
# Task \#1
### Symbol calculations
"""

# ╔═╡ 4acbb302-6726-11eb-1af2-0157c1a07336
begin
	@vars u₁ u₂
	u⁺ = [u₁ u₂]
	u = [
	  u₁
	  u₂
	]

	T₂ = convert(Array{Int}, CreateKTBC(2)[2])

	u⁺T₂u = u⁺ * T₂ * u
end

# ╔═╡ e6a566ba-6726-11eb-0ef2-e148ca198a89
md"""
#### Find the result
"""

# ╔═╡ 1cac6756-6727-11eb-2d9a-3776545b32da
result₁ = [u₁ * (u₁ - u₂) + u₂ * (-u₁ + 2 * u₂)]

# ╔═╡ 263872ca-6727-11eb-1233-c37ef9371211
md"""
### The result is:

`` u^T \cdot T₂ \cdot u = u₁^2 - 2u₁u₂ + u₂^2 ``

### Or you can see it by diff

`` u^T \cdot T₂ \cdot u - (u₁^2 - 2u₁u₂ + u₂^2) = 0 ``

#### If result is equal, difference is equal zero

"""

# ╔═╡ 70f8a80e-6727-11eb-0e9c-833f93a4c9f5
@test u⁺T₂u[1] - result₁[1] == 0

# ╔═╡ 7ded0a64-6727-11eb-1dd7-a353f7a80ba7
md"""
# Task \#3
### Produce Circulant matrix C
"""

# ╔═╡ 965b12bc-6727-11eb-0099-ddd1176df887
A = [
  1 -1 0
  0 1 -1
  -1 0 1
]

# ╔═╡ 9d00019a-6727-11eb-2a6f-3fda51c74b01
C = (A') * A

# ╔═╡ b3db8876-6727-11eb-1e64-1fa9949feee3
z = zeros(Int, 3)

# ╔═╡ c2c35d50-6727-11eb-1f1e-1161a06fa745
md"""
#### A contains dependent columns

```math
\begin{bmatrix}
  a_{11} \\
  a_{21} \\
  a_{31}
 \end{bmatrix} + \begin{bmatrix}
  a_{12} \\
  a_{22} \\
  a_{32}
 \end{bmatrix} + \begin{bmatrix}
  a_{13} \\
  a_{23} \\
  a_{33}
 \end{bmatrix} = \begin{bmatrix}
  0 \\
  0 \\
  0
\end{bmatrix}

```

"""

# ╔═╡ 8b584ee2-6728-11eb-2015-17f4ce49fc83
@test A[1:3] + A[4:6] + A[7:9] == z

# ╔═╡ 93d52f1a-6728-11eb-2733-cd6a61aa2b89
md"""
#### The same for C

```math
 \begin{bmatrix}
  c_{11} \\
  c_{21} \\
  c_{31}
 \end{bmatrix} + \begin{bmatrix}
  c_{12} \\
  c_{22} \\
  c_{32}
 \end{bmatrix} + \begin{bmatrix}
  c_{13} \\
  c_{23} \\
  c_{33}
 \end{bmatrix} = \begin{bmatrix}
  0 \\
  0 \\
  0
 \end{bmatrix}
```
"""

# ╔═╡ cdc18b56-6728-11eb-1b95-c340bfb9d559
@test C[1:3] + C[4:6] + C[7:9] == z

# ╔═╡ e2ccc9c8-6728-11eb-11c0-7fe79dc8a96d
md"""
#### Find the vector u
"""

# ╔═╡ f2cb5eac-6728-11eb-1d1f-fbb306febaef
begin
	v = ones(Int, 3)
	@test A * v == z
end

# ╔═╡ 06eb6980-672a-11eb-0d64-5dc376fb893b
md"""
#### C -> is not Positive definite for Cholesky
"""

# ╔═╡ 1eff6daa-672a-11eb-2df1-b54942784659
@test_throws PosDefException cholesky(C)

# ╔═╡ 2f1b2972-672a-11eb-1df4-11ba022cdd62
md"""
# Task \#11
### Find symmetric matrix for function
`` f(x,y) = 2xy ``
"""

# ╔═╡ 6a32b714-672a-11eb-2e95-d99358a79a68
md"""
#### Partial derivatives table `` f(x, y) = 2xy ``

```math
	\begin{matrix}
		\frac{\mathrm \partial}{\mathrm \partial x}
			\left( f(x, y) \right) = 2y
		&
		\frac{\mathrm \partial}{\mathrm \partial y}
			\left( f(x, y) \right) = 2x
	\\
		\frac{\mathrm \partial}{\mathrm \partial x \partial y}
			\left( f(x, y) \right) = 2
	 	&
		\frac{\mathrm \partial}{\mathrm \partial y \partial x}
			\left( f(x, y) \right) = 2
	\\
		\frac{\mathrm \partial^2}{\mathrm \partial x^2}
			\left( f(x, y) \right) = 0
		&
		\frac{\mathrm \partial^2}{\mathrm \partial y^2}
			\left( f(x, y) \right) = 0
	\\
 \end{matrix}
```

#### Hessian matrix

```math

\begin{bmatrix}
	\frac{\mathrm \partial}{\mathrm \partial x}
		\left( f(x, y) \right)
	&
	\frac{\mathrm \partial}{\mathrm \partial x \partial y}
		\left( f(x, y) \right)
	\\
	\frac{\mathrm \partial}{\mathrm \partial y \partial x}
		\left( f(x, y) \right)
	&
	\frac{\mathrm \partial}{\mathrm \partial y}
		\left( f(x, y) \right)
\end{bmatrix} \cdot \begin{bmatrix}
	x
	\\
	y
\end{bmatrix}

=

\begin{bmatrix}
	2 & 0
	\\
	0 & 2
\end{bmatrix} \cdot \begin{bmatrix}
	x
	\\
	y
\end{bmatrix}

=

\begin{bmatrix}
	2x
	\\
	2y
\end{bmatrix}

```

"""

# ╔═╡ 60f181b4-672a-11eb-16d8-bd6b5b5725cd
f(x, y) = 2 * x * y

# ╔═╡ 561285de-672a-11eb-33f9-3d9265291867
@vars x y a b c

# ╔═╡ 476a3532-6779-11eb-1ff1-9d670dda615f
md"""
#### Let's init symbolic ``S^{symbolic}=S^s``
"""

# ╔═╡ 71a3c0f2-6779-11eb-2fd1-cb4e54b8cfa8
begin
	w⁺ = [x y]
	w = [
	  x
	  y
	]
	Sˢ = [
		a b
		b c
	]
	w⁺ * Sˢ * w
end

# ╔═╡ 01eec5e6-672b-11eb-1bda-b38590341dac
md"""

#### Let's find a, b, c

```math
	\begin{matrix}
		x \left(a x + b y\right) + y \left(b x + c y\right)
		&
		=
		&
		ax^2 + 2byx + cy^2
	\\
		ax^2 = cy^2 = 0
		&
		=>
		&
		a = c = 0
	\\
		2byx = 2yx
		&
		=>
		&
		b = 1
 \end{matrix}

```

"""

# ╔═╡ 4eee9e0a-672a-11eb-2ab1-39f67d4ff891
S = [
  0 1
  1 0
]

# ╔═╡ eacb5054-672a-11eb-0cc7-4b6b9661e804
@test (w⁺ * S * w) - [2*x*y] == [0]

# ╔═╡ 817df440-672c-11eb-16ef-eb398ab32a0c
begin
	vals, vects = eigen(S)
	@test vals == [
		-1.0
		1.0
	]
end

# ╔═╡ c2f65df8-6782-11eb-34eb-a914ea97d400
md"""
### Task \# 24
"""

# ╔═╡ d40dd6e8-6782-11eb-2460-6dec703edafc
begin
	K = CreateKTBC(5)[1]
end

# ╔═╡ 17e154ac-6789-11eb-1936-7176d5492f77
round.(inv(K), digits=2)

# ╔═╡ 17620080-678c-11eb-142c-0be670b9ec02
md"""
``K`` is symmetric as ``K^{-1}``

``K^{-1} = (K^{-1})^{T}``
"""

# ╔═╡ 28d4134e-6789-11eb-1561-dfad67ceb90c
@test round.(inv(K), digits=2) == round.(inv(K)', digits=2)

# ╔═╡ Cell order:
# ╟─0668a58e-6727-11eb-2b6c-41626284bd0d
# ╠═1c3315a0-6774-11eb-0081-f328cca176ed
# ╠═1f74c496-6726-11eb-3894-119ff752210e
# ╟─9020f922-6725-11eb-2bc5-c190c647a78a
# ╠═4acbb302-6726-11eb-1af2-0157c1a07336
# ╟─e6a566ba-6726-11eb-0ef2-e148ca198a89
# ╠═1cac6756-6727-11eb-2d9a-3776545b32da
# ╟─263872ca-6727-11eb-1233-c37ef9371211
# ╠═70f8a80e-6727-11eb-0e9c-833f93a4c9f5
# ╟─7ded0a64-6727-11eb-1dd7-a353f7a80ba7
# ╠═965b12bc-6727-11eb-0099-ddd1176df887
# ╠═9d00019a-6727-11eb-2a6f-3fda51c74b01
# ╠═b3db8876-6727-11eb-1e64-1fa9949feee3
# ╟─c2c35d50-6727-11eb-1f1e-1161a06fa745
# ╠═8b584ee2-6728-11eb-2015-17f4ce49fc83
# ╟─93d52f1a-6728-11eb-2733-cd6a61aa2b89
# ╠═cdc18b56-6728-11eb-1b95-c340bfb9d559
# ╟─e2ccc9c8-6728-11eb-11c0-7fe79dc8a96d
# ╠═f2cb5eac-6728-11eb-1d1f-fbb306febaef
# ╟─06eb6980-672a-11eb-0d64-5dc376fb893b
# ╠═1eff6daa-672a-11eb-2df1-b54942784659
# ╟─2f1b2972-672a-11eb-1df4-11ba022cdd62
# ╟─6a32b714-672a-11eb-2e95-d99358a79a68
# ╠═60f181b4-672a-11eb-16d8-bd6b5b5725cd
# ╠═561285de-672a-11eb-33f9-3d9265291867
# ╟─476a3532-6779-11eb-1ff1-9d670dda615f
# ╠═71a3c0f2-6779-11eb-2fd1-cb4e54b8cfa8
# ╟─01eec5e6-672b-11eb-1bda-b38590341dac
# ╠═4eee9e0a-672a-11eb-2ab1-39f67d4ff891
# ╠═eacb5054-672a-11eb-0cc7-4b6b9661e804
# ╠═817df440-672c-11eb-16ef-eb398ab32a0c
# ╟─c2f65df8-6782-11eb-34eb-a914ea97d400
# ╠═d40dd6e8-6782-11eb-2460-6dec703edafc
# ╠═17e154ac-6789-11eb-1936-7176d5492f77
# ╟─17620080-678c-11eb-142c-0be670b9ec02
# ╠═28d4134e-6789-11eb-1561-dfad67ceb90c


text:
using IncrementalInference
using Test

##

@testset "Test solveCliqueUp! and solveCliqDown!" begin

##

N=8
fg = generateCanonicalFG_lineStep(N;
                                  graphinit=false,
                                  poseEvery=1,
                                  landmarkEvery=N+1,
                                  posePriorsAt=[0],
                                  landmarkPriorsAt=[],
                                  sightDistance=N+1)
#

deleteFactor!.(fg, [Symbol("x$(i)lm0f1") for i=1:(N-1)])

# test the initAll! separately anyway
initAll!(fg)

tree = buildTreeReset!(fg)

# for debuggin use
# ENV["JULIA_DEBUG"] = :csm_2

# solve clique up tests
a,b = solveCliqUp!(fg, tree, 2)
a,b = solveCliqUp!(fg, tree, 2; recordcliq = true)
@test length(a) > 0

# solve clique down tests
a,b = solveCliqDown!(fg, tree, 2)
a,b = solveCliqDown!(fg, tree, 2; recordcliq = true)
@test length(a) > 0

##

end


text:
process = arr ->
    if length(arr) == 1
        arr[1]
    else
        arr
    end


text:
abstract AbstractCell{T}

type RNNCell{T<:AbstractFloat} <: AbstractCell{T}
  Wh::AbstractMatrix{T}
  Wx::AbstractMatrix{T}
  Wo::AbstractMatrix{T}
  ht
end

function compute(c::RNNCell, xt)
    xt = reshape(xt, (1, input_size))
    ht = tanh.(xt * c.Wx .+ c.ht * c.Wh)
    return c.ht * c.Wo
end

type LSTMCell{T<:AbstractFloat} <: AbstractCell
end

function compute(c::LSTMCell, x, h)
end


text:
abstract Status

type Passed <: Status end
type Failed <: Status end
type Skipped <: Status end

type Result{S <: Status}
    codeblock::AbstractString
    exception::Exception
    location::@compat(Tuple{Any, Int})

    Result(codeblock, location) = new(codeblock, ErrorException(""), location)
    Result(codeblock, exception, location) = new(codeblock, exception, location)
end

function writemime{S <: Status}(io::IO, mime::MIME"text/plain", r::Result{S})
    println(io, """
     Entry: $(r.location[1])
     Block: $(r.location[2])
    Status: $(S)
    """)
    indented(io, r.codeblock)
    printexception(io, mime, r)
end

function printexception(io::IO, ::MIME"text/plain", r::Result{Failed})
    println(io, "\n", "Exception:", "\n")
    indented(io, string(r.exception))
end
printexception(io::IO, ::MIME"text/plain", r::Result{Passed}) = println(io)
printexception(io::IO, ::MIME"text/plain", r::Result{Skipped}) = nothing

function indented(io::IO, str::AbstractString, n::Int = 4)
    for line in split(str, "\n")
        println(io, " "^n, line)
    end
end

type Results{S <: Status}
    results::Vector{Result{S}}

    Results() = new(Results{S}[])
end

function writemime(io::IO, mime::MIME"text/plain", rs::Results)
    for r in rs.results
        writemime(io, mime, r)
    end
end

length(rs::Results) = length(rs.results)

type Summary
    modname::Module
    passed::Results{Passed}
    failed::Results{Failed}
    skipped::Results{Skipped}

    Summary(modname) = new(modname, Results{Passed}(), Results{Failed}(), Results{Skipped}())
end

for (T, f, docs) in [(Passed, :passed, "passed"),
                     (Failed, :failed, "failed"),
                     (Skipped, :skipped, "were skipped during")]
    @eval begin
        $(f)(s::Summary) = s.$(f)
        push!(s::Summary, r::Result{$(T)}) = push!($(f)(s).results, r)
    end
end

function writemime(io::IO, mime::MIME"text/plain", s::Summary)
    println(io, "[module summary]\n")
    println(io, """
     *   module: $(s.modname)
     * exported: $(length(names(s.modname)))
     *  entries: $(length(entries(documentation(s.modname))))
    """)

    npassed, nfailed, nskipped = map(length, (passed(s), failed(s), skipped(s)))
    total = npassed + nfailed + nskipped

    println(io, "[doctest summary]\n")
    @printf(io, " * pass: %5d / %d\n", npassed, total)
    @printf(io, " * fail: %5d / %d\n", nfailed, total)
    @printf(io, " * skip: %5d / %d\n", nskipped, total)
end

"""
Run code blocks in the docstrings of the specified module `modname`.
Returns a `Summary` of the results.

Code blocks may be skipped by adding an extra newline at the end of the block.

**Example:**

```julia_skip
doctest(Lexicon)
```
"""
function doctest(modname::Module)
    doc = documentation(modname)
    println("running doctest on $(modname)...")
    summ = Summary(modname)
    for (obj, entry) in entries(doc)
        isa(docs(entry), Docile.Interface.Docs{:md}) || continue # Markdown is the only supported format.
        count = 0
        for block in parsed(docs(entry)).content
            if iscode(block)
                count += 1
                try
                    if !runnable(block) # skip code block with trailing newline
                        push!(summ, Result{Skipped}(block.code, (obj, count)))
                        continue
                    end
                    runcode(block, modname)
                catch err
                    push!(summ, Result{Failed}(block.code, err, (obj, count)))
                    continue
                end
                push!(summ, Result{Passed}(block.code, (obj, count)))
            end
        end
    end
    summ
end

iscode(::Markdown.Code) = true
iscode(::Any)           = false

runnable(code::Markdown.Code) =
    contains(code.language, "julia") &&
    !(contains(code.language, "skip"))

function runcode(block, modname)
    sandbox = Module(:__SANDBOX__)
    eval(sandbox, Expr(:toplevel, Expr(:using, symbol("$(modname)"))))
    i = 1
    while i < length(block.code)
        ex, i = parse(block.code, i)
        eval(sandbox, ex)
    end
end


text:
using DataDeps

register(DataDep(
  "housing",
  "https://archive.ics.uci.edu/ml/datasets/Housing",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data",
  "baadf72995725d76efe787b664e1f083388c79ba21ef9a7990d87f774184735a",
  post_fetch_method=(path -> begin
    # First char of each row is whitespace, remove
    write(path, join(strip.(readlines(path)), "\n"))

    UCIData.process_dataset(path,
        target_index=14,
        feature_indices=1:13,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
function x, cost, iters_used, last_Delta_x = one_d_minimizer(seed, func; tol=1e-5, maxiter=100, start_eta=0.1)
Minimizes a 1-d function using constrained Hessian minimization.
We don't trust the long-range info from the Hessian too much, meaning that there's
a given (adaptive) step size. If Newton's method suggests a step smaller than that step
size, we take it. Otherwise, we only move in the direction of the gradient by the stepsize.
Adaptive step size: Every step that the cost function gets smaller, the step grows by a factor
of 1.1. If a step would have led to a larger cost function, the step is not taken, and
the size falls by a factor of 2.
PARAMETERS:
===========
seed       A float, the starting point for the minimization
func       A one dimensional function, takes a float and returns a float that represents the current cost.
OPTIONAL PARAMETERS:
====================
tol=1e-5   If a step would lead to a change in cost that is smaller in magnitude than tol, stop the minimization.
maxiter=100   Maximum number of iterations for the minimization
start_eta=0.1  The starting value of the step size
RETURNS:
========
x0     the value of x that minimizes f(x)
cost   The value of f(x0)
niters  The number of iterations done
Dparam  The last change to the parameter value
"""
function one_d_minimizer(seed, func; tol=1e-5, maxiter=100, start_eta=0.1)
    eta = start_eta;
    lambdavec = [seed]

    cost, grad, hess = vgh(func, lambdavec)

    i = delta_lambda = 0;  # declare it so we can use it after the if
    for i in [1:maxiter;]
        h_delta = - grad./hess;
        if abs(h_delta[1]) < eta
            new_lambdavec = lambdavec + h_delta
        else
            new_lambdavec = lambdavec - eta*sign.(grad)
        end

        delta_lambda = new_lambdavec[1] - lambdavec[1]
        if abs(delta_lambda) < tol
            break
        end

        new_cost, new_grad, new_hess = vgh(func, new_lambdavec)

        if new_cost .< cost
            lambdavec[1] = new_lambdavec[1]
            cost = new_cost;
            grad = new_grad;
            hess = new_hess;

            eta = eta*1.1
        else
            eta = eta/2
        end

        @printf "%d: cost=%.3f lambda=%.3f\n" i cost lambdavec[1]
    end

    return lambdavec[1], cost, i, delta_lambda
end

# DON'T MODIFY THIS FILE -- the source is in file Constrained Parabolic Minimization.ipynb. Look there for further documentation and examples of running the code.

"""
function x, J, lambdas, costs, niters, Dlambda = constrained_parabolic_minimization(H, G, r; tol=1e-6,
    min_only=true, do_plot=false,
    verbose=false, efactor=3.0, max_efactor_tries=10,
    lambdastepsize=0.003, minimum_tol=1e-24, tol_delta=1e-3, maxiter=200)
Given a Hessian matrix, a gradient vector, and a desired radius from the origin, finds the vector
that minimizes the parabola defined by the Hessian and the gradient, subject to the constraint that the
vector's length equals the desired radius.
# PARAMETERS:
- H      A square symmetric matrix.
- G      A vector, length equal to the size(H,2)
- r      desired radius
# OPTIONAL PARAMETERS:
- tol=1e-6        Numerical tolerance on the computations
- min_only=true   Return only the minimum, or, if false, all xs for all lambdas that match x'*x = r^2
                    (some of thos might be maxima, not minima)
- efactor=3       The initial exploration of lambdas will go from -efactor(max(absolute(eig(H)))) to +efactor(max(absolute(eig(H))))
                If that does not produce a solution within the indicates tolerance, then efactor is multiplied by efactor_growth
                and we try again, up to max_efactor_tries.
- efactor_growth  see efactor
- max_efactor_tries   see efactor
- minimum_tol     Each candidate value of lambda, identified as a minimum in the grid scan and
                proposed as satisfying |x|^2=r^2, is then used to seed a Newton's method one-d minimization with a certain tolerance.
                If the results of that search do not satisfy |x|^2=r^2 within the desired tolerance, then the tolerance
                for the 1-d search is reduced (to make the one-d search more exact), by a factor of tol_delta,
                up to minimum_tol
- tol_delta       see minimum_tol
- max_iter        Maximum number of iterations in the 1-d search (see minimum_tol).
- lambdastepsize=0.003    The step size for initial exploration of lambdas, un units of efactor. It sshould
                probably scale with the smallest difference in the eigenvalues of H; that has not been implemented yet.
- do_plot         If true, produces plot of the initial gridscans of lambda versus (x'*x - r^2)^2
- verbose         If true, produces diagnostic info as it goes.
RETURNS:
========
- x        The vector that minimizes 0.5*x'*H*x + x'*G subject to x'*x = r

- J        0.5*x'*H*x + x'*G at the returned x
- lambda   value of the Lagrange multiplier at which the radius constraint is satisfied
- c        The squared difference between the length of x and r. Should be small, otherwise somthing went wrong!
- niters   The number of iterations used in the one-d minimization that identified the output lambdas.
- Dlambda  The last change in lambda during the one-d minimization. If the one-d minimization did not reach
         its maxiters, then this will be smaller than the one-d minimization's tolerance.
"""
function constrained_parabolic_minimization(H, G, r; tol=1e-6, min_only=true,
    do_plot=false, verbose=false, efactor=3.0, efactor_growth=4, max_efactor_tries=10,
    lambdastepsize=0.003, minimum_tol=1e-24, tol_delta=1e-3, maxiter=200)

    #  --- First a couple of helper functions ----

    """
    function x_of_lambda(lambda)
    Given square matrix H, vector G, and passed scalar lambda, returns the vector x that minimizes

    0.5 x'*H*x + x'*G - lambda *x'*x
    """
    function x_of_lambda(lambda)
        return inv(H - lambda*eye(size(H,1)))*(-G)
    end

    """
    function q(lambda, r)
    Returns the squared difference between r and the norm of x_of_lambda(lambda).
    """
    function q(lambda, r)
        return (r - norm(x_of_lambda(lambda)))^2
    end

    # efactor is the factor that multiplies the biggest eigenvalue of H, to determine the range over which we'll
    # look for a lambda that satisfies the norm(x)==r requirement. If we don't find a solution, we iteratively
    # increase efactor to try to get there, for a maximum of max_efactor_tries
    for m=1:max_efactor_tries
        # First scan lambda to find good candidates for minimizing the parabolic
        # surface under the x'*x = r^2 constraint
        L = eig(H)[1]
        L0 = maximum(abs.(L))
        lambdas = L0*efactor*[-1.0:lambdastepsize:1.0;]
        costs = zeros(size(lambdas))
        for i in [1:length(lambdas);]
            try
                costs[i] = q(lambdas[i], r)
            catch
                costs[i] = Inf
            end
        end

        # if do_plot
        #     figure(2); clf();
        #     plot(lambdas, costs, "b.-")
        #     xlabel("lambda")
        #     ylabel("cost")
        # end
        println("here")

        # Take all candidates where the derivative of costs changes sign
        # from negative to positive (those would be minima),
        # plus the smallest and the largest lambdas tested, as candidates
        g = append!(prepend!(find(diff(sign.(diff(costs))) .> 0.99), [1]), [length(lambdas)])
        if verbose
            @printf("cpm: g (candidate indices) are : ");           print_vector_g(g);        print("\n")
            @printf("cpm: and their corresponding costs are : ");   print_vector(costs[g]);   print("\n");
            @printf("cpm: and their corresponding lambdas are : "); print_vector(lambdas[g]); print("\n");
            @printf("cpm: the minimum cost was : %g\n", minimum(costs[g]))
        end
        # found_it_flag = 0  # A flag for when we've found at least one lambda that satisfies the r constraint
        mytol = tol

        while mytol > minimum_tol
            lambdas_out = zeros(size(g))
            costs_out   = zeros(size(g))
            niters_out  = zeros(size(g))
            Dlambda_out = zeros(size(g))
            for i in [1:length(g);]
                lambdas_out[i], costs_out[i], niters_out[i], Dlambda_out[i] =
                one_d_minimizer(lambdas[g[i]], x -> q(x[1], r), start_eta=1, tol=mytol, maxiter=maxiter)
            end

            # Eliminate any lambdas where x'*x doesn't match our desired value r
            I = find(costs_out .< tol)
            lambdas_out = lambdas_out[I]; costs_out = costs_out[I];
            niters_out  = niters_out[I];  Dlambda_out = Dlambda_out[I]

            if length(I) > 0; break; end

            mytol *= tol_delta
        end
        if verbose
            @printf("%d : After searching for lambdas with efactor=%g, we found these : ", m, efactor)
            print_vector_g(lambdas_out); print("\n")
        end
        if length(lambdas_out) > 0; break; end;
        efactor = efactor*efactor_growth
    end

    # Eliminate any repeated lambdas, to within the specified numerical tolerance.
    I = setdiff(1:length(lambdas_out), find(diff(lambdas_out) .< tol))
    lambdas_out = lambdas_out[I]; costs_out = costs_out[I];
    niters_out  = niters_out[I];  Dlambda_out = Dlambda_out[I]

    # Find the parabolic estimate of the cost function at these points
    J  = zeros(size(lambdas_out))
    xs = zeros(length(G), length(lambdas_out))
    for i in [1:length(J);]
        xs[:,i] = x_of_lambda(lambdas_out[i])
        J[i] = (0.5*xs[:,i]'*H*xs[:,i] + xs[:,i]'*G)[1]
    end

    # Find and return only the x that has the smallest J
    if min_only
        I = indmin(J)
    else
        I = 1:length(J)
    end
    return xs[:,I], J[I], lambdas_out[I], costs_out[I], niters_out[I], Dlambda_out[I]
end


text:
struct Error{T}
    value::T
end

struct EventResult{F}
    f::F
end


text:
include("wavelet.jl")

using Plots; contourf, plot

@doc raw"""
Represents a wavelet transformed. Beside the wavelet coefficients, this object contains additional
information about the wavelet transformed including the scales at which the transform was performed
as well as the wavelet used. This object extends `AbstractArray{Complex{Float64}, 2}` thus allowing
for direct access to the wavelet coefficients. Each voice is stored as a column, thus the number of
rows matches the number of samples of the original time-series.

# Fields
- `wX`: The actual wavelet coefficients, each column represents a single voice.
- `scales`: The scales at which the wavelet transform was performed.
- `wavelet`: The wavelet of the transform.
"""
struct ContinuousWaveletTransformed <: AbstractArray{Complex{Float64}, 2}
    wX::Array{Complex{Float64}, 2}
    scales::Array{Float64, 1}
    wavelet::GenericContinuousWavelet
end

Base.size(A::ContinuousWaveletTransformed) = size(A.wX)
Base.getindex(A::ContinuousWaveletTransformed, I::Vararg{Int, N}) where {N} = getindex(A.wX, I...);
Base.setindex!(A::ContinuousWaveletTransformed, v, I::Vararg{Int, N}) where {N} = setindex!(A.wX, v, I...)

@doc raw"""
    Plots.contourf(A::ContinuousWavelet.ContinuousWaveletTransformed; drawvalid=true, shadevalid=true, kw...)

Plots the modulus (absolute value) of the wavelet transformed using a filled-contour plot.

If the optional keyword argument `drawvalid` is `true` (default), the area of valid wavelet
coefficients is marked by a thick black line. If the optional keyword arguement `shadevalid` is
also `true` (default) the invalid wavelet coefficients are shaded gray. A valid wavelet coefficient
is one, where the accociated wavelet does not overlapp with the edges of the time-series.

Additional keyword arguments are passed to the default implementation of `Plots.contourf()`.
"""
function Plots.contourf(A::ContinuousWaveletTransformed; drawvalid=true, shadevalid=true, kw...)
    # Draw transformed as heatmap
    N,M = size(A)
    cn = contourf(1:N, A.scales, transpose(map(abs, A.wX));
         linewidth=0, linetype=:heatmap, xlims=(1,N), ylims=(minimum(A.scales),maximum(A.scales)),
         kw...);
    # draw valid range
    if drawvalid
        x = [1.0]; y = Float64[A.scales[1]]
        for i in 1:M
            tc = cutoff_time(A.wavelet)*A.scales[i]
            if (tc < N÷2)
                push!(y, A.scales[i])
                push!(x, tc)
            end
        end
        append!(x, reverse(N.-x))
        append!(y, reverse(y))
        if shadevalid
            plot!(x,y; linecolor=:black, linewidth=3, label="",
                  fillrange=A.scales[end], fillalpha=0.5, fillcolor=:black)
        else
            plot!(x,y; linecolor=:black, linewidth=3, label="")
        end
    end
    return cn
end


text:
# It is suggested that test code for MLJ.jl include files be placed in
# a file of the same name under "test/" (and included below) and that
# this test code be wrapped in a module. Any new module name will do -
# eg, `module TestDatasets` for code testing `datasets.jl`.

# using Revise
using MLJBase
using Test

@testset "utilities" begin
  @test include("utilities.jl")
end

@testset "parameters" begin
  @test include("parameters.jl")
end

@testset "scitypes" begin
  @test include("scitypes.jl")
end

@testset "distributions" begin
  @test include("distributions.jl")
end

@testset "data" begin
  @test include("data.jl")
end

@testset "info" begin
  @test include("info.jl")
end

@testset "datasets" begin
  @test include("datasets.jl")
end

@testset "tasks" begin
  @test include("tasks.jl")
end