text:
using OrdinaryDiffEq, DiffEqProblemLibrary,Base.Test, DiffEqBase

bools = Vector{Bool}(0)
prob = prob_ode_linear

sol =solve(prob,Euler(),dt=1//2^(2),dense=true)

interpd_1d = sol(0:1//2^(4):1)

sol2 =solve(prob,Euler(),dt=1//2^(4),dense=true)

sol3 =solve(prob,Euler(),dt=1//2^(5),dense=true)

prob = prob_ode_2Dlinear
sol =solve(prob,Euler(),dt=1//2^(2),dense=true)

interpd = sol(0:1//2^(4):1)

interpd_idxs = sol(0:1//2^(4):1,idxs=1:2:5)

@test minimum([interpd_idxs[i] == interpd[i][1:2:5] for i in 1:length(interpd)])

interpd_single = sol(0:1//2^(4):1,idxs=1)

@test typeof(interpd_single) <: Vector{Float64}

@test typeof(sol(0.5,idxs=1)) <: Float64

A = rand(4,2)
sol(A,0.777)
A == sol(0.777)
sol2 =solve(prob,Euler(),dt=1//2^(4),dense=true)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < .2

sol(interpd,0:1//2^(4):1)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < .2

sol =solve(prob,Euler(),dt=1//2^(2),dense=false)

@test_throws BoundsError sol(0.5)

prob = prob_ode_linear

sol =solve(prob,Midpoint(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Midpoint(),dt=1//2^(4),dense=true)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 2e-2

prob = prob_ode_2Dlinear

sol =solve(prob,Midpoint(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Midpoint(),dt=1//2^(4),dense=true)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 2.3e-2

prob = prob_ode_linear

sol =solve(prob,RK4(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,RK4(),dt=1//2^(4),dense=true)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 5e-5

prob = prob_ode_2Dlinear

sol =solve(prob,RK4(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,RK4(),dt=1//2^(4),dense=true)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-4

prob = prob_ode_linear

sol =solve(prob,DP5(),dense=true)

sol2 =solve(prob,DP5(),dt=1//2^(4),dense=true,adaptive=false)

sol(interpd_1d,0:1//2^(4):1)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 5e-6

prob = prob_ode_2Dlinear

sol =solve(prob,DP5(),dense=true)

sol2 =solve(prob,DP5(),dt=1//2^(4),dense=true,adaptive=false)

sol(interpd,0:1//2^(4):1)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-5

prob = prob_ode_linear

sol =solve(prob,BS3(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,BS3(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 1e-3

prob = prob_ode_2Dlinear

sol =solve(prob,BS3(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,BS3(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-3

prob = prob_ode_linear

sol =solve(prob,Tsit5(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Tsit5(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 1e-5

prob = prob_ode_2Dlinear

sol =solve(prob,Tsit5(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Tsit5(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-5

prob = prob_ode_linear

sol =solve(prob,TanYam7(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,TanYam7(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 1e-3

prob = prob_ode_2Dlinear

sol =solve(prob,TanYam7(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,TanYam7(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-3

prob = prob_ode_linear

sol =solve(prob,TsitPap8(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,TsitPap8(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 2e-2

prob = prob_ode_2Dlinear

sol =solve(prob,TsitPap8(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,TsitPap8(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-2

prob = prob_ode_linear

sol =solve(prob,Feagin10(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Feagin10(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 1e-3

prob = prob_ode_2Dlinear

sol =solve(prob,Feagin10(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Feagin10(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-3

### Vern6()
const linear_bigα = parse(BigFloat,"1.01")
f = (t,u) -> (linear_bigα*u)
prob_ode_bigfloatlinear = ODEProblem(f,parse(BigFloat,"0.5"),(0.0,1.0))
prob = prob_ode_bigfloatlinear

sol =solve(prob,Vern6(),dt=1//2^(2),dense=true)

interpd_1d_big = sol(0:1//2^(7):1)

sol2 =solve(prob,Vern6(),dt=1//2^(7),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d_big)) < 1e-7

#plot(sol2.t,interpd)
#plot!(sol.t,sol[:])
#scatter!(sol.t,sol[:])

prob_ode_bigfloatveclinear = ODEProblem(f,[parse(BigFloat,"0.5")],(0.0,1.0))
prob = prob_ode_bigfloatveclinear

sol =solve(prob,Vern6(),dt=1//2^(2),dense=true)

interpd_big = sol(0:1//2^(4):1)

sol2 =solve(prob,Vern6(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_big)) < 1e-7

### BS5()

prob = prob_ode_linear

sol =solve(prob,BS5(),dt=1//2^(1),dense=true,adaptive=false)

interpd_1d_long = sol(0:1//2^(7):1)

sol2 =solve(prob,BS5(),dt=1//2^(7),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d_long)) < 2e-7

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,BS5(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,BS5(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 2e-7

### Vern7()

prob = prob_ode_linear

sol =solve(prob,Vern7(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Vern7(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 3e-9

prob = prob_ode_2Dlinear

sol =solve(prob,Vern7(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Vern7(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 5e-9

### Vern8()

prob = prob_ode_linear

sol =solve(prob,Vern8(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Vern8(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 1e-7

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,Vern8(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Vern8(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-7

### Vern9()

prob = prob_ode_linear

sol =solve(prob,Vern9(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Vern9(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 1e-9

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,Vern9(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Vern9(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 2e-9

### Rosenbrock23()

prob = prob_ode_linear

sol =solve(prob,Rosenbrock23(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Rosenbrock23(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 3e-3

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,Rosenbrock23(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Rosenbrock23(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 6e-3

### Rosenbrock32()

prob = prob_ode_linear

sol =solve(prob,Rosenbrock32(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Rosenbrock32(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 3e-4

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,Rosenbrock32(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Rosenbrock32(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 6e-4

### ImplicitEuler()

prob = prob_ode_linear

sol =solve(prob,ImplicitEuler(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,ImplicitEuler(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 0.2

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,ImplicitEuler(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,ImplicitEuler(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < .354

### Trapezoid()

prob = prob_ode_linear

sol =solve(prob,Trapezoid(),dt=1//2^(2),dense=true)

sol(interpd_1d,0:1//2^(4):1)

sol2 =solve(prob,Trapezoid(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d)) < 7e-3

# plot(sol2.t,interpd)

prob = prob_ode_2Dlinear

sol =solve(prob,Trapezoid(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,Trapezoid(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1.4e-2

### DP8()

prob = prob_ode_linear

sol =solve(prob,DP8(),dt=1//2^(2),dense=true)

sol(interpd_1d_long,0:1//2^(7):1)

sol2 =solve(prob,DP8(),dt=1//2^(7),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d_long)) < 2e-7

prob = prob_ode_2Dlinear

sol =solve(prob,DP8(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,DP8(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 2.01e-7

### ExplicitRK()

prob = prob_ode_linear

sol =solve(prob,ExplicitRK(),dt=1//2^(2),dense=true)

sol(interpd_1d_long,0:1//2^(7):1)

sol2 =solve(prob,ExplicitRK(),dt=1//2^(7),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d_long)) < 6e-5

prob = prob_ode_2Dlinear

sol =solve(prob,ExplicitRK(),dt=1//2^(2),dense=true)

sol(interpd,0:1//2^(4):1)

sol2 =solve(prob,ExplicitRK(),dt=1//2^(4),dense=true,adaptive=false)

@test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd)) < 1e-3


text:
module CategoricalDistributions

import Distributions
import ScientificTypesBase: Finite, Multiclass, OrderedFactor
using OrderedCollections
using CategoricalArrays
import Missings
using Random
using UnicodePlots

const Dist = Distributions
const MAX_NUM_LEVELS_TO_SHOW_BARS = 12

import Distributions: pdf, logpdf, support, mode

include("utilities.jl")
include("types.jl")
include("methods.jl")
include("arrays.jl")
include("arithmetic.jl")

export UnivariateFinite, UnivariateFiniteArray, UnivariateFiniteVector

# re-eport from Distributions:
export pdf, logpdf, support, mode

# re-export from ScientificTypesBase:
export Multiclass, OrderedFactor

end

ext: jl


text:
const PyJuliaIOValue_Type__ref = Ref(PyNULL)
PyJuliaIOValue_Type() = begin
    ptr = PyJuliaIOValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.IOValue",
                base = base,
                iter = pyjlio_iter,
                iternext = pyjlio_next,
                methods = [
                    (name = "close", flags = Py_METH_NOARGS, meth = pyjlio_close),
                    (name = "fileno", flags = Py_METH_NOARGS, meth = pyjlio_fileno),
                    (name = "flush", flags = Py_METH_NOARGS, meth = pyjlio_flush),
                    (name = "isatty", flags = Py_METH_NOARGS, meth = pyjlio_isatty),
                    (name = "readable", flags = Py_METH_NOARGS, meth = pyjlio_readable),
                    (name = "writable", flags = Py_METH_NOARGS, meth = pyjlio_writable),
                    (name = "tell", flags = Py_METH_NOARGS, meth = pyjlio_tell),
                    (name = "seek", flags = Py_METH_VARARGS, meth = pyjlio_seek),
                    (name = "writelines", flags = Py_METH_O, meth = pyjlio_writelines),
                    (name = "readlines", flags = Py_METH_VARARGS, meth = pyjlio_readlines),
                    (name = "truncate", flags = Py_METH_VARARGS, meth = pyjlio_truncate),
                    (name = "seekable", flags = Py_METH_NOARGS, meth = pyjlio_seekable),
                    (name = "__enter__", flags = Py_METH_NOARGS, meth = pyjlio_enter),
                    (name = "__exit__", flags = Py_METH_VARARGS, meth = pyjlio_exit),
                    # extras
                    (name = "torawio", flags = Py_METH_NOARGS, meth = pyjlio_torawio),
                    (
                        name = "tobufferedio",
                        flags = Py_METH_NOARGS,
                        meth = pyjlio_tobufferedio,
                    ),
                    (name = "totextio", flags = Py_METH_NOARGS, meth = pyjlio_totextio),
                ],
                getset = [(name = "closed", get = pyjlio_closed)],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyIOBase_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaIOValue_Type__ref[] = ptr
    end
    ptr
end

const PyJuliaRawIOValue_Type__ref = Ref(PyNULL)
PyJuliaRawIOValue_Type() = begin
    ptr = PyJuliaRawIOValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaIOValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.RawIOValue",
                base = base,
                methods = [
                    (name = "write", flags = Py_METH_O, meth = pyjlio_write_bytes),
                    (name = "readall", flags = Py_METH_VARARGS, meth = pyjlio_read_bytes),
                    (name = "read", flags = Py_METH_VARARGS, meth = pyjlio_read_bytes),
                    (
                        name = "readline",
                        flags = Py_METH_VARARGS,
                        meth = pyjlio_readline_bytes,
                    ),
                    (name = "readinto", flags = Py_METH_O, meth = pyjlio_readinto),
                    # extras
                    (name = "torawio", flags = Py_METH_NOARGS, meth = pyjlio_identity),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyRawIOBase_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaRawIOValue_Type__ref[] = ptr
    end
    ptr
end

const PyJuliaBufferedIOValue_Type__ref = Ref(PyNULL)
PyJuliaBufferedIOValue_Type() = begin
    ptr = PyJuliaBufferedIOValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaIOValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.BufferedIOValue",
                base = base,
                methods = [
                    (name = "detach", flags = Py_METH_NOARGS, meth = pyjlio_detach),
                    (name = "write", flags = Py_METH_O, meth = pyjlio_write_bytes),
                    (name = "read", flags = Py_METH_VARARGS, meth = pyjlio_read_bytes),
                    (name = "read1", flags = Py_METH_VARARGS, meth = pyjlio_read_bytes),
                    (
                        name = "readline",
                        flags = Py_METH_VARARGS,
                        meth = pyjlio_readline_bytes,
                    ),
                    (name = "readinto", flags = Py_METH_O, meth = pyjlio_readinto),
                    (name = "readinto1", flags = Py_METH_O, meth = pyjlio_readinto),
                    # extras
                    (name = "tobufferedio", flags = Py_METH_NOARGS, meth = pyjlio_identity),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyBufferedIOBase_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaBufferedIOValue_Type__ref[] = ptr
    end
    ptr
end

const PyJuliaTextIOValue_Type__ref = Ref(PyNULL)
PyJuliaTextIOValue_Type() = begin
    ptr = PyJuliaTextIOValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaIOValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.TextIOValue",
                base = base,
                methods = [
                    (name = "detach", flags = Py_METH_NOARGS, meth = pyjlio_detach),
                    (name = "write", flags = Py_METH_O, meth = pyjlio_write_str),
                    (name = "read", flags = Py_METH_VARARGS, meth = pyjlio_read_str),
                    (
                        name = "readline",
                        flags = Py_METH_VARARGS,
                        meth = pyjlio_readline_str,
                    ),
                    # extras
                    (name = "totextio", flags = Py_METH_NOARGS, meth = pyjlio_identity),
                ],
                getset = [
                    (name = "encoding", get = pyjlio_encoding),
                    (name = "errors", get = pyjlio_errors),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyTextIOBase_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaTextIOValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaIOValue_New(x::IO) = PyJuliaValue_New(PyJuliaIOValue_Type(), x)
PyJuliaRawIOValue_New(x::IO) = PyJuliaValue_New(PyJuliaRawIOValue_Type(), x)
PyJuliaBufferedIOValue_New(x::IO) = PyJuliaValue_New(PyJuliaBufferedIOValue_Type(), x)
PyJuliaTextIOValue_New(x::IO) = PyJuliaValue_New(PyJuliaTextIOValue_Type(), x)
PyJuliaValue_From(x::IO) = PyJuliaBufferedIOValue_New(x)

pyjlio_torawio(xo::PyPtr, ::PyPtr) = PyJuliaRawIOValue_New(PyJuliaValue_GetValue(xo)::IO)
pyjlio_tobufferedio(xo::PyPtr, ::PyPtr) =
    PyJuliaBufferedIOValue_New(PyJuliaValue_GetValue(xo)::IO)
pyjlio_totextio(xo::PyPtr, ::PyPtr) = PyJuliaTextIOValue_New(PyJuliaValue_GetValue(xo)::IO)
pyjlio_identity(xo::PyPtr, ::PyPtr) = (Py_IncRef(xo); xo)

pyjlio_iter(xo::PyPtr) = (Py_IncRef(xo); xo)

pyjlio_next(xo::PyPtr) = begin
    rl = PyObject_GetAttrString(xo, "readline")
    isnull(rl) && return PyNULL
    line = PyObject_CallNice(rl)
    Py_DecRef(rl)
    isnull(line) && return PyNULL
    if PyObject_Length(line) > 0
        line
    else
        Py_DecRef(line)
        PyNULL
    end
end

pyjlio_enter(xo::PyPtr, ::PyPtr) = (Py_IncRef(xo); xo)

pyjlio_exit(xo::PyPtr, args::PyPtr) = begin
    cl = PyObject_GetAttrString(xo, "close")
    isnull(cl) && return PyNULL
    v = PyObject_CallNice(cl)
    Py_DecRef(cl)
    Py_DecRef(v)
    isnull(v) && return PyNULL
    PyNone_New()
end

pyjlio_close(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        close(x)
        PyNone_New()
    catch err
        if err isa MethodError && err.f === close
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_closed(xo::PyPtr, ::Ptr{Cvoid}) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(!isopen(x))
    catch err
        if err isa MethodError && err.f === isopen
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_fileno(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(fd(x))
    catch err
        if err isa MethodError && err.f === fd
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_flush(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        flush(x)
        PyNone_New()
    catch err
        if err isa MethodError && err.f === flush
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_isatty(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(x isa Base.TTY)
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlio_readable(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(isreadable(x))
    catch err
        if err isa MethodError && err.f === isreadable
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_writable(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(iswritable(x))
    catch err
        if err isa MethodError && err.f === iswritable
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_tell(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(position(x))
    catch err
        if err isa MethodError && err.f === position
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_seek(xo::PyPtr, args::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    ism1(PyArg_CheckNumArgsBetween("seek", args, 1, 2)) && return PyNULL
    ism1(PyArg_GetArg(Int, "seek", args, 0)) && return PyNULL
    n = takeresult(Int)
    ism1(PyArg_GetArg(Int, "seek", args, 1, 0)) && return PyNULL
    w = takeresult(Int)
    try
        if w == 0
            seek(x, n)
        elseif w == 1
            seek(x, position(x) + n)
        elseif w == 2
            seekend(x)
            seek(x, position(x) + n)
        else
            PyErr_SetString(PyExc_ValueError(), "'whence' argument must be 0, 1 or 2 (got $w)")
            return PyNULL
        end
        PyObject_From(position(x))
    catch err
        if err isa MethodError && (err.f === seek || err.f === position || err.f === seekend)
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_writelines(xo::PyPtr, lines::PyPtr) = begin
    wr = PyObject_GetAttrString(xo, "write")
    isnull(wr) && return PyNULL
    r = PyIterable_Map(lines) do line
        r = PyObject_CallNice(wr, PyObjectRef(line))
        isnull(r) ? -1 : 1
    end
    r == -1 ? PyNULL : PyNone_New()
end

pyjlio_readlines(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("readlines", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Int, "readlines", args, 0, -1)) && return PyNULL
    limit = takeresult(Int)
    rl = PyObject_GetAttrString(xo, "readline")
    isnull(rl) && return PyNULL
    lines = PyList_New(0)
    isnull(lines) && (Py_DecRef(rl); return PyNULL)
    nread = 0
    while limit < 0 || nread < limit
        line = PyObject_CallNice(rl)
        isnull(line) && (Py_DecRef(lines); Py_DecRef(rl); return PyNULL)
        len = PyObject_Length(line)
        ism1(len) && (Py_DecRef(line); Py_DecRef(lines); Py_DecRef(rl); return PyNULL)
        len == 0 && (Py_DecRef(line); break)
        err = PyList_Append(lines, line)
        Py_DecRef(line)
        ism1(err) && (Py_DecRef(lines); Py_DecRef(rl); return PyNULL)
        nread += len
    end
    Py_DecRef(rl)
    lines
end

pyjlio_truncate(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("truncate", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Union{Int,Nothing}, "truncate", args, 0, nothing)) &&
        return PyNULL
    n = takeresult(Union{Int,Nothing})
    x = PyJuliaValue_GetValue(xo)::IO
    try
        m = n === nothing ? position(x) : n
        truncate(x, m)
        PyObject_From(m)
    catch err
        if err isa MethodError && (err.f === truncate || err.f === position)
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_seekable(xo::PyPtr, ::PyPtr) = begin
    T = typeof(PyJuliaValue_GetValue(xo)::IO)
    PyObject_From(hasmethod(seek, Tuple{T,Int}) && hasmethod(position, Tuple{T}))
end

pyjlio_encoding(::PyPtr, ::Ptr{Cvoid}) = PyUnicode_From("UTF-8")
pyjlio_errors(::PyPtr, ::Ptr{Cvoid}) = PyUnicode_From("strict")
pyjlio_detach(o::PyPtr, ::PyPtr) = (
    PyErr_SetString(PyExc_ValueError(), "cannot detach '$(PyType_Name(Py_Type(o)))'"); PyNULL
)

pyjlio_write_bytes(xo::PyPtr, bo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    b = Ref(Py_buffer())
    ism1(PyObject_GetBuffer(bo, b, PyBUF_SIMPLE)) && return PyNULL
    s = unsafe_wrap(Vector{UInt8}, Ptr{UInt8}(b[].buf), b[].len)
    try
        write(x, s)
        PyObject_From(length(s))
    catch err
        if err isa MethodError && err.f === write
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    finally
        PyBuffer_Release(b)
    end
end

pyjlio_write_str(xo::PyPtr, so::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::IO
    ism1(PyObject_Convert(so, String)) && return PyNULL
    s = takeresult(String)
    sep = ""
    havesep = false
    try
        i = firstindex(s)
        while true
            j = findnext('\n', s, i)
            if j === nothing
                write(x, SubString(s, i))
                break
            else
                write(x, SubString(s, i, prevind(s, j)))
                if !havesep
                    m = Py_OSModule()
                    isnull(m) && return PyNULL
                    sepo = PyObject_GetAttrString(m, "linesep")
                    isnull(sepo) && return PyNULL
                    r = PyObject_TryConvert(sepo, String)
                    Py_DecRef(sepo)
                    r == -1 && return PyNULL
                    r == 0 && (
                        PyErr_SetString(PyExc_TypeError(), "os.linesep must be a str"); return PyNULL
                    )
                    sep = takeresult(String)
                    havesep = true
                end
                write(x, sep)
                i = nextind(s, j)
            end
        end
        PyObject_From(length(s))
    catch err
        if err isa MethodError && err.f === write
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlio_readinto(xo::PyPtr, bo::PyPtr) = begin
    b = Ref(Py_buffer())
    ism1(PyObject_GetBuffer(bo, b, PyBUF_WRITABLE)) && return PyNULL
    x = PyJuliaValue_GetValue(xo)::IO
    a = unsafe_wrap(Vector{UInt8}, Ptr{UInt8}(b[].buf), b[].len)
    try
        PyObject_From(readbytes!(x, a, length(a)))
    catch err
        if err isa MethodError && err.f === readbytes!
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    finally
        PyBuffer_Release(b)
    end
end

readpybytes(io::IO, limit = nothing) = begin
    buf = UInt8[]
    len = 0
    while (limit === nothing || len < limit) && !eof(io)
        c = read(io, UInt8)
        push!(buf, c)
    end
    buf
end

pyjlio_read_bytes(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("read", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Union{Int,Nothing}, "read", args, 0, nothing)) && return PyNULL
    limit = takeresult(Union{Int,Nothing})
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyBytes_From(readpybytes(x, (limit === nothing || limit < 0) ? nothing : limit))
    catch err
        if err isa MethodError && (err.f === eof || err.f === read || err.f === peek)
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

readpystr(io::IO, limit = nothing) = begin
    buf = IOBuffer()
    len = 0
    while (limit === nothing || len < limit) && !eof(io)
        c = read(io, Char)
        if c == '\n'
            write(buf, '\n')
        elseif c == '\r'
            !eof(io) && peek(io) == 0x0A && read(io, UInt8)
            write(buf, '\n')
        else
            write(buf, c)
            len += 1
        end
    end
    String(take!(buf))
end

pyjlio_read_str(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("read", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Union{Int,Nothing}, "read", args, 0, nothing)) && return PyNULL
    limit = takeresult(Union{Int,Nothing})
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(readpystr(x, (limit === nothing || limit < 0) ? nothing : limit))
    catch err
        if err isa MethodError && (err.f === eof || err.f === read || err.f === peek)
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

readpybyteline(io::IO, limit = nothing) = begin
    buf = UInt8[]
    len = 0
    while (limit === nothing || len < limit) && !eof(io)
        c = read(io, UInt8)
        push!(buf, c)
        c == 0x0A && break
    end
    buf
end

pyjlio_readline_bytes(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("readline", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Union{Int,Nothing}, "readline", args, 0, nothing)) &&
        return PyNULL
    limit = takeresult(Union{Int,Nothing})
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyBytes_From(readpybyteline(x, (limit === nothing || limit < 0) ? nothing : limit))
    catch err
        if err isa MethodError && (err.f === eof || err.f === read || err.f === peek)
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

readpyline(io::IO, limit = nothing) = begin
    buf = IOBuffer()
    len = 0
    while (limit === nothing || len < limit) && !eof(io)
        c = read(io, Char)
        if c == '\n'
            write(buf, '\n')
            break
        elseif c == '\r'
            !eof(io) && peek(io) == 0x0A && read(io, UInt8)
            write(buf, '\n')
            break
        else
            write(buf, c)
            len += 1
        end
    end
    String(take!(buf))
end

pyjlio_readline_str(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("readline", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Union{Int,Nothing}, "readline", args, 0, nothing)) &&
        return PyNULL
    limit = takeresult(Union{Int,Nothing})
    x = PyJuliaValue_GetValue(xo)::IO
    try
        PyObject_From(readpyline(x, (limit === nothing || limit < 0) ? nothing : limit))
    catch err
        if err isa MethodError && (err.f === eof || err.f === read || err.f === peek)
            PyErr_SetStringFromJuliaError(PyExc_ValueError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end


text:
# This file was generated, do not modify it. # hide
ypred = predict(svc, X)
misclassification_rate(ypred, y)


text:
"""

The exported names are
$(EXPORTS)
"""
module ReferenceFEs

using Test
using DocStringExtensions
using LinearAlgebra
using Combinatorics

using InplaceArrays.Helpers
using InplaceArrays.TensorValues
using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Polynomials

import InplaceArrays.Arrays: kernel_cache
import InplaceArrays.Arrays: apply_kernel!
import InplaceArrays.Arrays: kernel_return_type
import InplaceArrays.Fields: evaluate
import InplaceArrays.Polynomials: MonomialBasis

import Base: ==

export Polytope
export ExtrusionPolytope
export get_faces
export get_dimrange
export get_vertex_coordinates
export get_facet_normals
export get_facet_orientations
export get_edge_tangents
export get_vertex_permutations
export get_vtkid
export get_vtknodes
export num_dims
export num_vertices
export num_faces
export num_facets
export num_edges
export get_facedims
export get_offsets
export get_offset
export test_polytope
export VERTEX
export SEGMENT
export TRI
export QUAD
export TET
export HEX
export WEDGE
export PYRAMID
export HEX_AXIS
export TET_AXIS
export INVALID_PERM

export Dof
export evaluate_dof!
export evaluate_dof
export dof_cache
export dof_return_type
export test_dof
export evaluate_dof_array

export ReferenceFE
export GenericRefFE
export get_polytope
export get_prebasis
export get_dofs
export get_face_own_dofids
export get_own_dofs_permutations
export get_shapefuns
export compute_shapefuns
export test_reference_fe
export num_dofs

export LagrangianRefFE
export LagrangianDofBasis
export compute_monomial_basis
export compute_own_nodes
export compute_face_orders
export compute_nodes
export compute_own_nodes_permutations
export compute_lagrangian_reffaces
export get_node_coordinates
export get_dof_to_node
export get_dof_to_comp
export get_node_and_comp_to_dof

export SerendipityRefFE

include("Polytopes.jl")

include("ExtrusionPolytopes.jl")

include("Dofs.jl")

include("MockDofs.jl")

include("LagrangianDofBases.jl")

include("ReferenceFEInterfaces.jl")

include("LagrangianRefFEs.jl")

include("SerendipityRefFEs.jl")

end # module


text:
function gss(arr::Vector{<:Number})
    smax = hmax = tmax = 0
    for head in eachindex(arr), tail in head:length(arr)
        s = sum(arr[head:tail])
        if s > smax
            smax = s
            hmax, tmax = head, tail
        end
    end
    return arr[hmax:tmax]
end

arr = [-1, -2, 3, 5, 6, -2, -1, 4, -4, 2, -1]
subseq = gss(arr)
s = sum(subseq)

println("Greatest subsequential sum of $arr:\n → $subseq with sum $s")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ExplainableAI
using Suppressor

# Flux layers
unknown_function(x) = x
@test check_model(:LRP, Chain(Dense(2, 2, relu)))
@test_throws ArgumentError("Unknown or unsupported activation functions found in model") check_model(
    :LRP, Chain(Dense(2, 2, softmax)); verbose=false
)
@test_throws ArgumentError("Unknown layers found in model") check_model(
    :LRP, Chain(unknown_function); verbose=false
)
@test_throws ArgumentError("Unknown layers found in model") @suppress check_model(
    :LRP,
    Chain(
        unknown_function,
        Chain(unknown_function),
        Parallel(+, unknown_function, unknown_function),
    );
    verbose=false,
)

# Custom layers
## Test using a simple wrapper
struct MyLayer{T}
    x::T
end
TestLayer = MyLayer(Dense(2, 2, relu))
@test_throws ArgumentError("Unknown layers found in model") check_model(
    :LRP, Chain(TestLayer); verbose=false
)
@test_throws ArgumentError("Unknown layers found in model") LRPZero(
    Chain(TestLayer); verbose=false
)
@test_nowarn LRPZero(Chain(TestLayer); skip_checks=true)
## Test should pass after registering the layer
LRP_CONFIG.supports_layer(::MyLayer) = true
@test check_model(:LRP, Chain(TestLayer); verbose=false) == true
@test_nowarn LRPZero(Chain(TestLayer))
## ...repeat for layers that are functions
@test_throws ArgumentError("Unknown layers found in model") check_model(
    :LRP, Chain(unknown_function); verbose=false
)
LRP_CONFIG.supports_layer(::typeof(unknown_function)) = true
@test check_model(:LRP, Chain(unknown_function); verbose=false) == true
## ...repeat for activation functions
@test_throws ArgumentError("Unknown or unsupported activation functions found in model") check_model(
    :LRP, Chain(Dense(2, 2, unknown_function)); verbose=false
)
LRP_CONFIG.supports_activation(::typeof(unknown_function)) = true
@test check_model(:LRP, Chain(Dense(2, 2, unknown_function)); verbose=false) == true


text:
"""
    Elliptic{T <: Number} <: Construct{T}

An elliptic Cayley-Dickson construct as a pair.
Elliptic constructs use the elliptic multiplication operation.
Examples of elliptic Cayley-Dickson constructs are binions (complex numbers), quaternions,
octonions, and sedenions.
"""
struct Elliptic{T <: Number} <: Construct{T}
    l::T
    r::T
end

"""
    (==)(x::CayleyDickson.Elliptic, y::CayleyDickson.Elliptic)
    (==)(x::CayleyDickson.Parabolic, y::CayleyDickson.Parabolic)
    (==)(x::CayleyDickson.Hyperbolic, y::CayleyDickson.Hyperbolic)

Return true if the left components of `x` and `y` are equal,
and also the right components of `x` and `y` are equal.
"""
function (==)(x::Elliptic, y::Elliptic)
    x.l == y.l && x.r == y.r
end

"""
    real(z::CayleyDickson.Elliptic)
    real(z::CayleyDickson.Parabolic)
    real(z::CayleyDickson.Hyperbolic)
For Cayley-Dickson constructs, return the real part of `z`.
This is equivalent to the real part of the left component of `z`.
"""
function real(z::Elliptic)
    real(z.l)
end

"""
    quadrance(z::CayleyDickson.Elliptic)

For elliptic Cayley-Dickson constructs, return the sum of the squared absolute value of the
left and right components of `z`. Note that this is always non-negative.

For the `Binion`, `Quaternion`, and `Octonion` types, `quadrance(z)` satisfies the composition
property:
```math
    quadrance(x * y) == quadrance(x) * quadrance(y)
```
The composition property does not hold true for the `Sedenion` type.
"""
function quadrance(z::Elliptic)
    quadrance(z.l) + quadrance(z.r)
end

"""
    iszerodivisor(z::CayleyDickson.Elliptic)

For `Binion`, `Quaternion`, and `Octonion` types, return `true` if `z` is zero.
"""
function iszerodivisor(z::Elliptic)
    iszero(z)
end

"""
    zero(::Type{CayleyDickson.Elliptic{T}}) where T <: Number
    zero(::Type{CayleyDickson.Parabolic{T}}) where T <: Number
    zero(::Type{CayleyDickson.Hyperbolic{T}}) where T <: Number

For Cayley-Dickson construct types, return a Cayley-Dickson
construct consisting of both left and right components being zero of type `T`.
"""
function zero(::Type{Elliptic{T}}) where T <: Number
    Elliptic{T}(zero(T), zero(T))
end

"""
    zero(::CayleyDickson.Elliptic{T}) where T <: Number
    zero(::CayleyDickson.Parabolic{T}) where T <: Number
    zero(::CayleyDickson.Hyperbolic{T}) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson
construct consisting of both left and right components being zero of type `T`.
"""
function zero(::Elliptic{T}) where T <: Number
    zero(Elliptic{T})
end

"""
    one(::Type{CayleyDickson.Elliptic{T}}) where T <: Number
    one(::Type{CayleyDickson.Parabolic{T}}) where T <: Number
    one(::Type{CayleyDickson.Hyperbolic{T}}) where T <: Number

For Cayley-Dickson construct types, return a Cayley-Dickson
construct where the left component is one of type `T`, and the right component
is zero of type `T`.
"""
function one(::Type{Elliptic{T}}) where T <: Number
    Elliptic{T}(one(T), zero(T))
end

"""
    one(::CayleyDickson.Elliptic{T}) where T <: Number
    one(::CayleyDickson.Parabolic{T}) where T <: Number
    one(::CayleyDickson.Hyperbolic{T}) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson
construct where the left component is one of type `T`, and the right component
is zero of type `T`.
"""
function one(::Elliptic{T}) where T <: Number
    one(Elliptic{T})
end

"""
    (-)(z::CayleyDickson.Elliptic)
    (-)(z::CayleyDickson.Parabolic)
    (-)(z::CayleyDickson.Hyperbolic)

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
minus the left component of `z`, and the right component is minus the right component of
`z`.
"""
function (-)(z::Elliptic)
    Elliptic(-(z.l), -(z.r))
end

"""
    conj(z::CayleyDickson.Elliptic)
    conj(z::CayleyDickson.Parabolic)
    conj(z::CayleyDickson.Hyperbolic)

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the conjugate of the left component of `z`, and the right component is minus the right component of `z`.
"""
function conj(z::Elliptic)
    Elliptic(conj(z.l), -(z.r))
end

"""
    (+)(x::Elliptic{T}, y::Elliptic{T}) where T <: Number
    (+)(x::Parabolic{T}, y::Parabolic{T}) where T <: Number
    (+)(x::Hyperbolic{T}, y::Hyperbolic{T}) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the sum of the left components of `x` and `y`, and the right component is the sum of the
right components of `x` and `y`.
"""
function (+)(x::Elliptic{T}, y::Elliptic{T}) where T <: Number
    Elliptic{T}(x.l + y.l, x.r + y.r)
end

"""
    (+)(a::T, z::Elliptic{T}) where T <: Number
    (+)(z::Elliptic{T}, a::T) where T <: Number
    (+)(a::T, z::Parabolic{T}) where T <: Number
    (+)(z::Parabolic{T}, a::T) where T <: Number
    (+)(a::T, z::Hyperbolic{T}) where T <: Number
    (+)(z::Hyperbolic{T}, a::T) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the sum of `a` and the left component of `z`, and the right component is the right component
of `z`. Note that `a` is of the same type as the components of `z`.
"""
function (+)(a::T, z::Elliptic{T}) where T <: Number
    Elliptic{T}(a + z.l, z.r)
end

function (+)(z::Elliptic{T}, a::T) where T <: Number
    Elliptic{T}(z.l + a, z.r)
end

"""
    (+)(a::Real, z::Elliptic)
    (+)(z::Elliptic, a::Real)
    (+)(a::Real, z::Parabolic)
    (+)(z::Parabolic, a::Real)
    (+)(a::Real, z::Hyperbolic)
    (+)(z::Hyperbolic, a::Real)

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the sum of `a` and the left component of `z`, and the right component is the right component
of `z` converted to the same type as the left component of the result.
Note that `a` is not necessarily of the same type as the real components of `z`.
"""
function (+)(a::Real, z::Elliptic)
    l = a + z.l
    Elliptic(l, oftype(l, z.r))
end

function (+)(z::Elliptic, a::Real)
    l = z.l + a
    Elliptic(l, oftype(l, z.r))
end

"""
    (-)(x::Elliptic{T}, y::Elliptic{T}) where T <: Number
    (-)(x::Parabolic{T}, y::Parabolic{T}) where T <: Number
    (-)(x::Hyperbolic{T}, y::Hyperbolic{T}) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the differemce of the left components of `x` and `y`, and the right component is the
differemce of the right components of `x` and `y`.
"""
function (-)(x::Elliptic{T}, y::Elliptic{T}) where T <: Number
    Elliptic{T}(x.l - y.l, x.r - y.r)
end

"""
    (-)(a::T, z::Elliptic{T}) where T <: Number
    (-)(a::T, z::Parabolic{T}) where T <: Number
    (-)(a::T, z::Hyperbolic{T}) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the difference of `a` and the left component of `z`, and the right component is minus the
right component of `z`. Note that `a` is of the same type as the components of `z`.
"""
function (-)(a::T, z::Elliptic{T}) where T <: Number
    Elliptic{T}(a - z.l, -(z.r))
end

"""
    (-)(z::Elliptic{T}, a::T) where T <: Number
    (-)(z::Parabolic{T}, a::T) where T <: Number
    (-)(z::Hyperbolic{T}, a::T) where T <: Number

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the difference of the left component of `z` and `a`, and the right component is the right
component of `z`. Note that `a` is of the same type as the components of `z`.
"""
function (-)(z::Elliptic{T}, a::T) where T <: Number
    Elliptic{T}(z.l - a, z.r)
end

"""
    (-)(a::Real, z::Elliptic)
    (-)(a::Real, z::Parabolic)
    (-)(a::Real, z::Hyperbolic)

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the difference of `a` and the left component of `z`, and the right component is minus the
right component of `z` converted to the same type as the left component of the result. Note
that `a` and the components of `z` are not necessarily of the same type.
"""
function (-)(a::Real, z::Elliptic)
    l = a - z.l
    Elliptic(l, -(oftype(l, z.r)))
end

"""
    (-)(z::Elliptic, a::Real)
    (-)(z::Parabolic, a::Real)
    (-)(z::Hyperbolic, a::Real)

For Cayley-Dickson constructs, return a Cayley-Dickson construct where the left component is
the difference of the left component of `z` and `a`, and the right component is the right
component of `z` converted to the same type as the left component of the result. Note
that `a` and the components of `z` are not necessarily of the same type.
"""
function (-)(z::Elliptic, a::Real)
    l = z.l - a
    Elliptic(l, oftype(l, z.r))
end

"""
    (*)(x::Elliptic{T}, y::Elliptic{T})

For elliptic Cayley-Dickson constructs, the elliptic multiplication operation.

The elliptic multiplication of two `Binion` types is commutative and associative.

The elliptic multiplication of two `Quaternion` types is non-commutative but associative.

The elliptic multiplication of two `Octonion` types is non-commutative, non-associative,
but alternative.

The elliptic multiplication of two `Sedenion` types is non-commutative, non-associative,
non-alternative but flexible.
"""
function (*)(x::Elliptic{T}, y::Elliptic{T}) where T <: Number
    Elliptic{T}(
        (x.l * y.l) - (conj(y.r) * x.r),
        (y.r * x.l) + (x.r * conj(y.l))
    )
end

"""
    (*)(a::T, z::Elliptic{T}) where T <: Number
    (*)(z::Elliptic{T}, a::T) where T <: Number
    (*)(a::T, z::Parabolic{T}) where T <: Number
    (*)(z::Parabolic{T}, a::T) where T <: Number
    (*)(a::T, z::Hyperbolic{T}) where T <: Number
    (*)(z::Hyperbolic{T}, a::T) where T <: Number
"""
function (*)(a::T, z::Elliptic{T}) where T <: Number
    Elliptic{T}(a * z.l, z.r * a)
end

function (*)(z::Elliptic{T}, a::T) where T <: Number
    Elliptic{T}(z.l * a, z.r * conj(a))
end

"""
    (*)(a::Real, z::Elliptic)
    (*)(z::Elliptic, a::Real)
    (*)(a::Real, z::Parabolic)
    (*)(z::Parabolic, a::Real)
    (*)(a::Real, z::Hyperbolic)
    (*)(z::Hyperbolic, a::Real)
"""
function (*)(a::Real, z::Elliptic)
    Elliptic(a * z.l, a * z.r)
end

function (*)(z::Elliptic, a::Real)
    Elliptic(z.l * a, z.r * a)
end

"""
    (/)(z::Elliptic, a::Real)
    (/)(z::Parabolic, a::Real)
    (/)(z::Hyperbolic, a::Real)
"""
function (/)(z::Elliptic, a::Real)
    if iszero(a)
        error(ZeroDenominator)
    end

    Elliptic(z.l / a, z.r / a)
end

"""
    (\\)(a::Real, z::Elliptic)
    (\\)(a::Real, z::Parabolic)
    (\\)(a::Real, z::Hyperbolic)
"""
function (\)(a::Real, z::Elliptic)
    if iszero(a)
        error(ZeroDenominator)
    end

    Elliptic(a \ z.l, a \ z.r)
end

"""
    widen(::Type{Elliptic{T}}) where T
    widen(::Type{Parabolic{T}}) where T
    widen(::Type{Hyperbolic{T}}) where T

For Cayley-Dickson construct types, return a wider Cayley-Dickson construct type.
"""
function widen(::Type{Elliptic{T}}) where T
    Elliptic{widen(T)}
end

"""
    float(::Type{Elliptic{T}}) where T <: AbstractFloat
    float(::Type{Parabolic{T}}) where T <: AbstractFloat
    float(::Type{Hyperbolic{T}}) where T <: AbstractFloat

For Cayley-Dickson construct of floating-point types, this is equivalent to the identity
operation.
"""
function float(::Type{Elliptic{T}}) where T <: AbstractFloat
    Elliptic{T}
end

"""
    float(::Type{Elliptic{T}}) where T
    float(::Type{Parabolic{T}}) where T
    float(::Type{Hyperbolic{T}}) where T

For Cayley-Dickson construct types, return an appropriate floating-point type to represent a
given Cayley-Dickson construct type.
"""
function float(::Type{Elliptic{T}}) where T
    Elliptic{float(T)}
end

"""
    cloak(a::Elliptic)

For elliptic Cayley-Dickson constructs, equivalent to the minus operation.
"""
function cloak(a::Elliptic)
    -(a)
end

"""
    dagger(a::Elliptic)

For elliptic Cayley-Dickson constructs, equivalent to the identity operation.
"""
function dagger(a::Elliptic)
    a
end

"""
    star(a::Elliptic)

For elliptic Cayley-Dickson constructs, the identity operation.
"""
function star(a::Elliptic)
    a
end


text:
function NOD(a,b)
    m = a
    n = b
    while n != 0
        m, n = n, m % n
    end
    return m
end


text:
using Base.Iterators
using Zygote

"""
    xlogx(x::Number)

Compute `x * log(x)`, returning zero if `x` is zero. (Copied from StatsFuns.jl package.)
"""
function xlogx(x::Number)
    result = x * log(x)
    ifelse(iszero(x), zero(result), result)
end

"""
Compute the von Neumann entropy of a density matrix `ρ`.
"""
neumann_entropy(ρ::AbstractMatrix) = -sum(xlogx.(real(eigvals(Matrix(ρ)))))

"""
    binary_entropy(p::Number)

Compute the binary entropy of `p` (natural logarithm).
"""
binary_entropy(p::Number) = -(xlogx(p) + xlogx(1 - p))

eig_entropy(θlist) = -sum(xlogx.(sum(θlist .* eigvs)))

"""
Compute `-log(Z)`, with `Z` the partition function.
"""
thermal_logZ(βH::AbstractMatrix) = -log(tr(exp(-Matrix(βH))))

"""
Use representation of density matrix in terms of Pauli matrices; cos(θ) is Bloch vector coefficient of Pauli-Z for each qubit
"""
latent_density(θlist) = DensityMatrix(kron([[1, 0, 0, cos(θ)] for θ in θlist]...), length(θlist))

"""
Trace distance to target density matrix `σ`.
"""
trdistance(θlist, σ, cgc) = 0.5 * opnorm(real.(matrix(apply(latent_density(θlist), cgc))) - σ, 1)

"""
Calculate fidelity as defined in arxiv:1910.02071
"""
function fidelity(θlist, σtherm, cgc)
    ρsqrt = sqrt(real.(matrix(apply(latent_density(θlist), cgc))))
    return real(tr(sqrt(ρsqrt * σtherm * ρsqrt))^2)
end

function trans_mat(L::Int)
    A = [1 0; 0 0; 0 0; 0 1]
    Avec = fill(A, L)
    kron(Avec...)
end

function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2-y1)/δ
        if eltype(x) <: Complex
            # derivative with respect to imaginary part
            x[i] = tmp - im*δ/2
            y1 = f(xs...)
            x[i] = tmp + im*δ/2
            y2 = f(xs...)
            x[i] = tmp
            Δ[i] += im*(y2-y1)/δ
        end
    end
    return grads
end

sigmoid(z::Real) = 1.0 / (1.0 + exp(-z))

max_stars_repo_name: cmendl/density-matrix-bloch-qcircuit

max_issues_repo_name: cmendl/density-matrix-bloch-qcircuit

max_forks_repo_name: cmendl/density-matrix-bloch-qcircuit


text:
using Plots
gr(xlims=[-pi,3*pi],ylims=[-0.4,1.2],border=:none)

Np_F = 4
h = 2*pi/Np_F

Sinc(x) = @. sin(pi*x/h)/(2*pi/h*tan(x/2))
r_plot = -pi:pi/200:3*pi
plot(r_plot,Sinc(r_plot),lw=3,legend=nothing)
plot!(r_plot,Sinc(r_plot.-h),lw=3,legend=nothing)
plot!(r_plot,Sinc(r_plot.-2*h),lw=3,legend=nothing)
plot!(r_plot,Sinc(r_plot.-3*h),lw=3,legend=nothing)

y = ones(size(r_plot))
l = 0.3
u = 0.3
plot!([-pi,0],[2,2],fill=(0,:grey),fillalpha=0.3)
plot!([-pi,0],-[2,2],fill=(0,:grey),fillalpha=0.3)
plot!([2*pi,3*pi],[2,2],fill=(0,:grey),fillalpha=0.3)
plot!([2*pi,3*pi],-[2,2],fill=(0,:grey),fillalpha=0.3)
scatter!([0;pi/2;pi;3*pi/2;2*pi],[0;0;0;0;0],markerstrokewidth=0,markersize=5,markercolor=:black)

savefig("sinc.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
x = [1.0, -2.0, 3.0, -4.0, 5.0, -6.0]
f = Conjugate(IndZero()) # = IndFree
grad_f_x, f_x = ProximalOperators.gradient(f, x)
@test iszero(grad_f_x)
@test iszero(f_x)
g = Conjugate(IndFree()) # = IndZero
prox_g_x, g_y = prox(g, x)
@test iszero(prox_g_x)
@test iszero(g_y)
l = Conjugate(SqrNormL2())
grad_l_x, l_x = ProximalOperators.gradient(l, x)
@test isequal(grad_l_x, x)
lam = 2
l = Conjugate(SqrNormL2(lam))
grad_l_x, l_x = ProximalOperators.gradient(l, x)
@test isequal(grad_l_x, x/lam)


text:
import StableRNGs: StableRNG

# create scattered data
npts = 50
xmin = 0.0
xspan = 3/2*π
x = collect(range(xmin, xmin+xspan, length=npts))
rng = StableRNG(655)
x = x + xspan/npts*(rand(rng,npts) .- 0.5)
# select a subset randomly
idx = unique(rand(rng, collect(eachindex(x)), 20))
t = x[unique(idx)]
npts = length(t)
ut = sin.(t)
stdev = 1e-1*maximum(ut)
u = ut + stdev*randn(rng, npts)
# data must be ordered if t̂ is not provided
idx = sortperm(t)
tₒ = t[idx]
uₒ = u[idx]

tolerance = 1e-3

@testset "Direct smoothing" begin
    # fixed with default λ = 1.0
    A = RegularizationSmooth(uₒ,tₒ; alg=:fixed)
    ans = [0.6456173647252937 0.663974701324226 0.7631218523665086 0.778654700697601 0.7489958320589535 0.7319087707475104 0.6807082599508811 0.6372557895089508 0.5832859790765743 0.5021274805916013 0.3065928203396211 0.1353332321156384 -0.3260000640060584 -0.6557906092739154 -0.9204882447932498]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # non-default d and λ
    A = RegularizationSmooth(uₒ,tₒ, 4; λ=1e-2, alg=:fixed)
    ans = [0.19865190868740357 0.2885349151737291 0.6756699442978945 0.9165887141895426 0.9936113717653254 1.0042825002191034 0.9768118192829827 0.9184595331808411 0.8214983284892922 0.6538356458824783 0.28295521578898 0.018060767871253963 -0.5301723647977373 -0.8349855890541111 -1.1085048455468356]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # GCV (default) to determine λ
    A = RegularizationSmooth(uₒ,tₒ)
    @test isapprox(A.λ, 0.12788440382063268, rtol=tolerance)
    ans = [0.21974931848164914 0.2973284508009968 0.6908546278415386 0.9300465474303226 0.9741453042418977 0.9767572556868123 0.9432951659303452 0.8889834700087442 0.804842790047182 0.6603217445567791 0.30341652659101737 0.05924456463634589 -0.5239939779242144 -0.8421768233191822 -1.107517099580091]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # L-curve to determine λ
    A = RegularizationSmooth(uₒ,tₒ; alg=:L_curve)
    @test isapprox(A.λ, 0.9536286111306728, rtol=tolerance)
    ans = [0.6261657429321232, 0.6470204841904836, 0.7599270022828396, 0.7835725197598925, 0.7567105872094757, 0.7404815750685363, 0.6906841961987067, 0.647628105931872, 0.5937273796308717, 0.512087780658067, 0.3136272387739983, 0.1392761732695201, -0.3312498167413961, -0.6673268474631847, -0.9370342562716745]
    @test isapprox(A.û, ans, rtol=tolerance)
end

@testset "Smoothing with weights" begin
    # midpoint rule integration
    A = RegularizationSmooth(uₒ,tₒ, nothing, :midpoint)
    @test isapprox(A.λ, 0.10787235405005478, rtol=tolerance)
    ans = [0.3068904607028622, 0.3637388879266782, 0.6654462500501238, 0.9056440536733456, 0.9738150157541853, 0.9821315604309402, 0.9502526946446999, 0.8953643918063283, 0.8024431779821514, 0.6415812230114304, 0.2834706832220367, 0.05281575111822609, -0.5333542714497277, -0.8406745098604134, -1.0983391396173634]
    @test isapprox(A.û, ans, rtol=tolerance)
    # arbitrary weights for wls (and fixed λ, GCV not working well for some of these)
    A = RegularizationSmooth(uₒ,tₒ, nothing, collect(1:npts); λ=1e-1, alg=:fixed)
    ans = [0.24640196218427968, 0.3212059975226125, 0.6557626475144205, 0.9222911426465459, 0.9913331910731215, 1.0072241662103494, 0.9757899817730779, 0.935880516370941, 0.8381074902073471, 0.6475589703422522, 0.2094170714475404, 0.09102085384961625, -0.5640882848240228, -0.810519277110118, -1.1159124134900906]
    @test isapprox(A.û, ans, rtol=tolerance)
    # arbitrary weights for wls and wr
    nhalf = Int(floor(npts/2))
    wls = vcat(ones(nhalf),10*ones(npts-nhalf))
    wr = collect(1:npts-2)
    A = RegularizationSmooth(uₒ,tₒ, nothing, wls, wr; λ=1e-1, alg=:fixed)
    ans = [0.21878709713242372, 0.3118480645325099, 0.7669822464946172, 1.0232343854914931, 1.0526513115274412, 1.0469579284244412, 0.9962426294084775, 0.9254407155702626, 0.8204764044515936, 0.6514510142804217, 0.27796896299068763, 0.04756024028728636, -0.5301034620974782, -0.8408107101140526, -1.1058428573417736]
    @test isapprox(A.û, ans, rtol=tolerance)
end

@testset "Smoothing with t̂ provided" begin
    N̂ = 20
    t̂ = collect(range(xmin, xmin+xspan, length=N̂))
    # with t̂, no weights
    A = RegularizationSmooth(u,t, t̂)
    @test isapprox(A.λ, 0.138273889585313, rtol=tolerance)
    ans = [0.21626377852882872 0.39235926952322575 0.5573848799950002 0.7072474496656729 0.8361906119247042 0.9313473799797176 0.9809844353757837 0.9750833208625507 0.9096038940899813 0.7816929736202427 0.6052694276527628 0.4015903497629387 0.1913719025253403 -0.01979786871512895 -0.23400354001942947 -0.44481229967011127 -0.6457913359497256 -0.8405146928672158 -1.0367229293434395 -1.2334090099343238]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # t̂ and wls
    A = RegularizationSmooth(u,t, t̂, collect(1:npts))
    @test isapprox(A.λ, 0.26746430253489195, rtol=tolerance)
    ans = [0.3118247878815087 0.44275860852897864 0.5705834985506882 0.6979119448253899 0.8234189540704866 0.9289458273102476 0.9970803470992273 1.0071205506077525 0.9443157518324818 0.7954860908242515 0.5847385548859145 0.34813493129868633 0.1237494751337505 -0.0823517516424196 -0.28265170846635246 -0.4760833187699964 -0.6615795059853024 -0.844779821396189 -1.0341162283806349 -1.225270266213379]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # t̂, wls, and wr
    nhalf = Int(floor(npts/2))
    wls = vcat(ones(nhalf),10*ones(npts-nhalf))
    wr = collect(1:N̂-2)
    A = RegularizationSmooth(u,t, t̂, wls, wr)
    @test isapprox(A.λ, 0.04555080890920959, rtol=tolerance)
    ans = [0.2799800686914433, 0.4627548444527547, 0.5611922868318674, 0.6647761469309206, 0.7910803348948329, 0.9096001134420562, 1.0067644979677808, 1.0541868144785513, 0.9889720466386331, 0.8088479651943575, 0.5677592185997403, 0.31309698432269184, 0.08587106716465115, -0.11476265128730469, -0.30749376694236485, -0.4942769809562725, -0.676806367664006, -0.8587832527770329, -1.0443430843364814, -1.2309001260104093]
    @test isapprox(A.û, ans, rtol=tolerance)
end


text:
export LinearVApproximator, LinearQApproximator

using LinearAlgebra: dot

"""
    LinearQApproximator(weights::Vector{Float64}, feature_func::F, actions::Vector{Int}) -> LinearQApproximator{F}

Use weighted sum to represent the estimation given a state and an action.

# Fields

- `weights::Vector{Float64}`: the weight of each feature.
- `feature_func::Function`: decide how to generate a feature vector of `length(weights)` given a state and an action as parameters.
- `actions::Vector{Int}`: all possible actions.

See also [`LinearVApproximator`](@ref).
"""
Base.@kwdef struct LinearQApproximator{F} <: AbstractQApproximator
    weights::Vector{Float64}
    feature_func::F
    actions::Vector{Int}
end

(Q::LinearQApproximator)(s, a::Int) = dot(Q.weights, Q.feature_func(s, a))

(Q::LinearQApproximator)(s) = [Q(s, a) for a in Q.actions]

function update!(Q::LinearQApproximator, correction::Pair)
    (s, a), e = correction
    xs = Q.feature_func(s, a)
    Q.weights .+= xs .* e
end


text:
function rot_n_stride!(traj::ContactTraj, stride::SizedArray)
    rotate!(traj)
    mpc_stride!(traj, stride)
    return nothing
end

function rotate!(traj::ContactTraj)
    H = traj.H
	# TODO: remove copy
    q1 = copy(traj.q[1])
    u1 = copy(traj.u[1])
    w1 = copy(traj.w[1])
    γ1 = copy(traj.γ[1])
    b1 = copy(traj.b[1])
    z1 = copy(traj.z[1])
    θ1 = copy(traj.θ[1])

    for t = 1:H+1
        traj.q[t] .= traj.q[t+1]
    end

    for t = 1:H-1
        traj.u[t] .= traj.u[t+1]
        traj.w[t] .= traj.w[t+1]
        traj.γ[t] .= traj.γ[t+1]
        traj.b[t] .= traj.b[t+1]
        traj.z[t] .= traj.z[t+1]
        traj.θ[t] .= traj.θ[t+1]
    end

    traj.q[end] .= q1
    traj.u[end] .= u1
    traj.w[end] .= w1
    traj.γ[end] .= γ1
    traj.b[end] .= b1
    traj.z[end] .= z1
    traj.θ[end] .= θ1

    return nothing
end

"""
    Update the last two cofiguaraton to be equal to the fisrt two up to an constant offset.
"""
function mpc_stride!(traj::ContactTraj, stride::SizedArray)
    H = traj.H

    for t = H+1:H+2
        traj.q[t] .= deepcopy(traj.q[t-H] + stride)
        update_z!(traj, t - 2)
        update_θ!(traj, t - 2)
        update_θ!(traj, t - 3)
    end

    return nothing
end

function get_stride(model::ContactModel, traj::ContactTraj) #TODO: dispatch over environment / model
    stride = zeros(SizedVector{model.dim.q})
    stride[1] = traj.q[end-1][1] - traj.q[1][1]
    return stride
end

# TODO: make more efficient / allocation free
function update_altitude!(alt, s::Simulation,  traj, t, N_sample;
	threshold = 1.0, verbose = false)

	# idx1 = max(0, t - N_sample) + 1
	idx1 = max(0, t-1 - N_sample) + 1 #@@@ simplify

	for i = 1:s.model.dim.c
		γ_max = 0.0
		idx_max = 0

		# for j = idx1:t
		for j = idx1:t-1
			if traj.γ[j][i] > γ_max
				γ_max =  traj.γ[j][i]
				idx_max = j
			end
		end
		# @show ([γ[1] for γ in traj.γ[max(1, t-3): t]])
		if γ_max > threshold
			alt[i] = ϕ_func(s.model, s.env, traj.q[idx_max+2])[i] #TODO check this this is correct
			# alt[i] = ϕ_func(model, traj.q[idx_max+2])[i]
			verbose && println(" ")
			verbose && println("point $i in contact")
			verbose && println("sim_step : $idx_max")
			verbose && println("alt      : $(ϕ_func(s.model, s.env, traj.q[idx_max+2])[i])")
			verbose && println("force    : $(traj.γ[idx_max][i])")
		end
	end
end

function live_plotting(model::ContactModel, ref_traj::ContactTraj,
		sim_traj::ContactTraj, newton::Newton, q0::AbstractVector{T},
		q1::AbstractVector{T}, t::Int) where {T}
	nq = model.dim.q
	nu = model.dim.u

	ql = 1
	qu = nq
	ul = 1
	uu = nu
	plt = plot(layout=grid(2,1,heights=[0.7, 0.3], figsize=[(1000, 1000),(400,400)]), legend=:bottomright, xlims=(0,20))
	plot!(plt[1,1], hcat([Vector(x[ql:qu]) for x in newton.traj.q]...)', color=:blue, linewidth=1.0, label=nothing)
	plot!(plt[1,1], hcat([Vector(x[ql:qu]) for x in ref_traj.q]...)', linestyle=:dot, color=:red, linewidth=3.0, label=nothing)

	scatter!((2-1/N_sample)ones(nq), sim_traj.q[t+0], markersize=8.0, color=:lightgreen, label="simul. q0 q1")
	scatter!(2*ones(nq), sim_traj.q[t+1], markersize=8.0, color=:lightgreen, label=nothing)
	# scatter!(plt[1,1], (2+1/N_sample)ones(nq), sim_traj.q[t+2], markersize=8.0, color=:lightgreen)

	# @show norm(sim_traj.q[t+1] - q1, Inf)
	scatter!(plt[1,1], 1*ones(nq), q0, markersize=6.0, color=:blue, label="newton q0 q1")
	scatter!(plt[1,1], 2*ones(nq), q1, markersize=6.0, color=:blue, label=nothing)

	scatter!(plt[1,1], 1*ones(nq), ref_traj.q[1], markersize=4.0, color=:red, label="refer. q0 q1 q2")
	scatter!(plt[1,1], 2*ones(nq), ref_traj.q[2], markersize=4.0, color=:red, label=nothing)
	scatter!(plt[1,1], 3*ones(nq), ref_traj.q[3], markersize=4.0, color=:red, label=nothing)

	plot!(plt[2,1], hcat([Vector(x[ul:uu]) for x in newton.traj.u]...)', color=:blue, linewidth=1.0, label="u newton")
	plot!(plt[2,1], hcat([Vector(x[ul:uu]) for x in ref_traj.u[1:end]]...)', linestyle=:dot, color=:red, linewidth=3.0, label="u refer.")
	display(plt)
end


text:
#=
# 2. Basic flow with a stationary body
In this notebook we will simulate the flow past a stationary body.
=#

using ViscousFlow
#-
using Plots

#=
### The basic steps
From the previous notebook, we add one additional step:
* **Specify the problem**: Set the Reynolds number and free stream
* **Discretize**: Set up a solution domain, grid cell size, time step size
* **Set up bodies**: *Create the body or bodies and specify their motions, if any*
* **Construct the system structure**: Create the operators that will be used to perform the simulation
* **Initialize**: Set the initial flow field and initialize the integrator
* **Solve**: Solve the flow field
* **Examine**: Examine the results
We will go through all of these here. For the examples we will carry out in this notebook,
the first three steps need only be carried out once.
=#

# ### Problem specification
# Set the Reynolds number and free stream
Re = 200 # Reynolds number
U = 1.0 # Free stream velocity
U∞ = (U,0.0);

#=
### Discretize
We set the grid Re to 4.0 here to get a quicker solution, though it is generally
better to make this smaller (it defaults to 2.0).
=#
xlim = (-1.0,3.0)
ylim = (-1.5,1.5)
Δx, Δt = setstepsizes(Re,gridRe=4.0)

#=
### Set up bodies
Here, we will set up a rectangle of half-height 0.5 and half-width 0.25
at 45 degrees angle of attack
=#
body = Rectangle(0.5,0.25,1.5Δx)

#=
We place the body at a desired location and orientation with the `RigidTransform`
function. This function creates an operator `T` that acts in-place on the body:
after the operation is applied, `body` is transformed to the correct location/orientation.
=#
cent = (0.0,0.0) # center of body
α = 45π/180 # angle
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration

# Let's plot it just to make sure
plot(body,xlim=xlim,ylim=ylim)

#=
### Construct the system structure
This step is like the previous notebook, but now we also provide the body and
the freestream:
=#
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,body,freestream = U∞)

#=
### Initialize
Now, we initialize with zero vorticity. Note that we do this by calling
`newstate` with no argument except for `sys` itself.
=#
u0 = newstate(sys)

#=
and now create the integrator, with a long enough time span to hold the whole
solution history:
=#
tspan = (0.0,20.0)
integrator = init(u0,tspan,sys)

#=
### Solve
Now we are ready to solve the problem. Let's advance the solution to $t = 1$.
=#
@time step!(integrator,1.0)

#=
### Examine
Let's look at the flow field at the end of this interval
=#
plot(
plot(vorticity(integrator),sys,title="Vorticity",clim=(-10,10),levels=range(-10,10,length=30), color = :RdBu,ylim=ylim),
plot(streamfunction(integrator),sys,title="Streamlines",ylim=ylim,color = :Black),
    size=(700,300)
    )

#=
Now let's make a movie, like we did last time.
=#

sol = integrator.sol;
@gif for (u,t) in zip(sol.u,sol.t)
    plot(vorticity(u,sys,t),sys,clim=(-10,10),levels=range(-10,10,length=30), color = :RdBu)
end every 5

#=
#### Compute the force history
To do this, we supply the solution history `sol`, the system `sys`, and the index
of the body (1).
=#
fx, fy = force(sol,sys,1);

#=
Plot the histories. Note that we are actually plotting the drag and lift
coefficient histories here:
$$ C_D = \dfrac{F_x}{\frac{1}{2}\rho U_\infty^2 L}, \quad C_L = \dfrac{F_y}{\frac{1}{2}\rho U_\infty^2 L} $$
Since the quantities in this simulation are already scaled by $\rho$, $U_\infty$, and $L$
(because $\rho$ has been scaled out of the equations, and the free stream speed is
set to 1 and the height of the shape to 1), then we obtain these coefficients by
simply dividing by 1/2, or equivalently, by multiplying by 2:
=#
plot(
plot(sol.t,2*fx,xlim=(0,Inf),ylim=(0,6),xlabel="Convective time",ylabel="\$C_D\$",legend=:false),
plot(sol.t,2*fy,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=:false),
    size=(800,350)
)

# The mean drag and lift coefficients are
meanCD = GridUtilities.mean(2*fx)
#-
meanCL = GridUtilities.mean(2*fy)


text:
export integratePosterior, wrapperIntegratePosterior, integratePosteriorParallel, iphelper

"""
function integratePosterior(param,rowfuns;K,wquad)

Compute posterior mean of a separable function

Input:
	1. param   - StrainReconParam
	2. rowfuns - functions phi_1,...,phi_m in an array, where phi_i is a function of the i'th block of the binary matrix and the frequency vector w (the target function whose mean is computed is the sum of these functions)

Keyword arguments:
	3. K	 - number of Monte Carlo points in Omega_W
	4. wquad - number of starting points w in Omega_W used to determine K for method overloading

Output:
	1. postmean - value of the expectation integral

Examples:
	Compute the posterior mean of the frequency vector w:
	rowfuns = [(r,w)->w;repmat([(r,w)->zeros(n,1)],m-1)]

	Compute the posterior mean of the binary matrix M:
	rowfuns = Array{Function,1}()
	for i = 1:m
		z1 = zeros((i-1)*(p-1),n)
		z2 = zeros((m-i)*(p-1),n)
		push!(rowfuns, (r,w)->[z1;r;z2])
	end

    Compute the expected misfit 0.5*dot(wres,wres), wres = gamma^0.5 * (Mw-d)
	rowfuns = Array{Function,1}()
	Rd = diag(param.R)
    for i = 1:m
		iblock = (i-1)*(p-1)+1 : i*(p-1)
        push!(rowfuns, (r,w)->[0.5*norm(Rd[iblock] .* (r*w-d[iblock]))^2])
	end
"""

function integratePosterior(param::StrainReconParam,rowfuns::Vector{Function};K::Int64=10000,wquad=randomw(param.n,K))
	K = size(wquad,2)

	n = param.n
	p = param.p
	d = param.d
	m = param.m
	Rd = diag(param.R)
	lambda = param.lambda

	# output size of rowfuns (assuming each has the same dimensions)
	funsize   = size(rowfuns[1](ones(p-1,n),ones(n,1)))
	funlength = prod(funsize)

	# Monte Carlo points, corresponding integrand values and normalization constants
	integrands = zeros(K,funlength)
	normconsts = zeros(K,1)
	logconsts = zeros(K,1)

	# generate all possible rows of M
	Mrows = getB(param)

	# loop through Monte Carlo points
	for idx1 = 1:K
		w = wquad[:,idx1]
		logposteriorvals = zeros(m,p^n)
		funvals = zeros(m,p^n,funlength)

		# loop through block indices
		for idx2 = 1:m
			blockidx2 = (idx2-1)*(p-1)+1 : idx2*(p-1)
			di = d[blockidx2]
			Wdi = Rd[blockidx2]

			# loop through all possible blocks
			for idx3 = 1:p^n
				blockidx3 = (idx3-1)*(p-1)+1 : idx3*(p-1)
				r = Mrows[blockidx3[1]:blockidx3[end],:]
				logposteriorvals[idx2,idx3] = -1/2 * norm(Wdi.*(r*w-di))^2 - lambda * sum(r)
				funvals[idx2,idx3,:] = rowfuns[idx2](r,w)[:]

			end

		end

		maxlogs = maximum(logposteriorvals,2)
		posteriorvals = exp.(logposteriorvals .- maxlogs)
		funposteriorvals = posteriorvals .* funvals

		a = sum(posteriorvals,2) # size = [m,1]
		b = squeeze(sum(funposteriorvals,2),2) # size = [m,funlength]

		c = zeros(m,funlength)
		for idx2 = 1:m
			c[idx2,:] = prod(a[[1:idx2-1;idx2+1:m]]) * b[idx2,:]
		end

		integrands[idx1,:] = sum(c,1)
		normconsts[idx1] = prod(a)
		logconsts[idx1] = sum(maxlogs)

	end

	maxmaxlog = maximum(logconsts)
	logscaling = logconsts - maxmaxlog
	integrands = integrands .* exp.(logscaling)
	normconsts = normconsts .* exp.(logscaling)

	postmean = sum(integrands,1) / sum(normconsts)
	postmean = reshape(postmean,funsize)

	return postmean

end

"""
	wrapperIntegratePosterior(param;K,wquad)

    Wrapper function for integratePosterior() to compute expected (M,w), misfit, and variance of (M,w)

    Input:
		1. param - StrainReconParam

    Optional arguments:
        2. K 	 - number of Monte Carlo points in Omega_W
		3. wquad - number of starting points w in Omega_W used to determine K for method overloading

    Output:
        1. Mmean  - E(M)
        2. wmean  - E(w)
        3. misfit - E(misfit)
		4. Mvar   - Var(M), i.e., E((M-E(M)).^2)
		5. wvar   - Var(w), i.e., E((w-E(w)).^2)
"""

function wrapperIntegratePosterior(param;K::Int64=10000,wquad=randomw(param.n,K))
	K = size(wquad,2)

	d = param.d
	m = param.m
	n = param.n
	p = param.p
	Rd = diag(param.R)

    # E(w)
    rowfuns = [(r,w)->w;repmat([(r,w)->zeros(n,1)],m-1)]
    wmean = integratePosteriorParallel(param,rowfuns,K=K,wquad=wquad)

	# Var(w)
	rowfuns = [(r,w)->(w-wmean).^2;repmat([(r,w)->zeros(n,1)],m-1)]
	wvar = integratePosteriorParallel(param,rowfuns,K=K,wquad=wquad)

    # E(M)
	mrowfuns = Array{Function,1}()
	for i = 1:m
        z1 = zeros((i-1)*(p-1),n)
        z2 = zeros((m-i)*(p-1),n)
        push!(mrowfuns, (r,w)->[z1;r;z2])
	end
    Mmean = integratePosteriorParallel(param,mrowfuns,K=K,wquad=wquad)

	# Var(M)
	mrowfuns = Array{Function,1}()
	for i = 1:m
        z1 = zeros((i-1)*(p-1),n)
        z2 = zeros((m-i)*(p-1),n)
        push!(mrowfuns, (r,w)->[z1;(r-Mmean[(i-1)*(p-1)+1 : i*(p-1), :]).^2;z2])
	end
    Mvar = integratePosteriorParallel(param,mrowfuns,K=K,wquad=wquad)

    # E(misfit)
    prowfuns = Array{Function,1}()
    for i = 1:m
		iblock = (i-1)*(p-1)+1 : i*(p-1)
        push!(prowfuns, (r,w)->[0.5*norm(Rd[iblock] .* (r*w-d[iblock]))^2])
	end
    misfit = integratePosteriorParallel(param,prowfuns,K=K,wquad=wquad)

    return Mmean, wmean, misfit, Mvar, wvar
end

"""
	function integratePosteriorParallel(param,rowfuns;K,wquad)

    Parallelized version of integratePosterior()

	Input:
		1. param   - StrainReconParam
		2. rowfuns - functions phi_1,...,phi_m in an array, where phi_i is a function of the i'th block of the binary matrix and the frequency vector w (the target function whose mean is computed is the sum of these functions)

	Keyword arguments:
		3. K	 - number of Monte Carlo points in Omega_W
		4. wquad - number of starting points w in Omega_W used to determine K for method overloading

    Output:
        1. postmean - value of expected variable
"""

function integratePosteriorParallel(param::StrainReconParam,rowfuns::Vector{Function};K::Int64=10000,wquad=randomw(param.n,K))
	K = size(wquad,2)

	n = param.n
	p = param.p
	d = param.d
	m = param.m
	Rd = diag(param.R)
	lambda = param.lambda

	# output size of rowfuns (assuming each has the same dimensions)
	funsize   = size(rowfuns[1](ones(p-1,n),ones(n,1)))
	funlength = prod(funsize)

	# Monte Carlo points, corresponding integrand values and normalization constants
	integrands = zeros(K,funlength)
	normconsts = zeros(K,1)
	logconsts = zeros(K,1)

	# generate all possible rows of M
	Mrows = getB(param)

    # update integrand values and normalization constants from MC points
    updateRes = (a,c,maxlogs,idx1)->(integrands[idx1,:] = sum(c,1);normconsts[idx1] = prod(a);logconsts[idx1] = sum(maxlogs))

	# loop through Monte Carlo points
    i=1
    nextidx = () -> (idx1=i; i+=1; idx1)
    @sync begin
        for pwork = workers()
            @async begin
                while true
                    idx1 = nextidx()
                    if isempty(idx1) || idx1 > K
                        break
                    end
                    a,c,maxlogs = remotecall_fetch(iphelper,pwork,m,n,p,d,Rd,rowfuns,Mrows,wquad[:,idx1],lambda,funlength)
                    updateRes(a,c,maxlogs,idx1)
                end
            end
        end
    end

	maxmaxlog = maximum(logconsts)
	logscaling = logconsts - maxmaxlog
	integrands = integrands .* exp.(logscaling)
	normconsts = normconsts .* exp.(logscaling)

	postmean = sum(integrands,1) / sum(normconsts)
	postmean = reshape(postmean,funsize)

	return postmean
end

"""
	function iphelper(m,n,p,d,Rd,rowfuns,Mrows,w,lambda,funlength)

    Helper function for integratePosteriorParallel

    Input:
        1. m 		  - number of SNP sites
        2. n 		  - number of strains
		3. p 		  - number of categories
        4. d 		  - data vector
		5. Rd 		  - diagonal of the inverse square root of the covariance matrix
        6. rowfuns	  - functions phi_1,...,phi_m in an array, where phi_i is a function of the i'th row of the binary matrix and the frequency vector w (the target function whose mean is computed is the sum of these functions)
        7. Mrows 	  - all possible rows of binary matrix
        8. w 		  - initial starting point of weight vector
        9. lambda 	  - sparsity parameter
        10. funlength - output size of rowfuns

	Output:
		1. a	   - sum of posterior values
		2. c	   - product of a and sum of normalized posterior values
		3. maxlogs - maximum of log posterior values to normalize values

"""
function iphelper(m,n,p,d,Rd,rowfuns,Mrows,w,lambda,funlength)

    logposteriorvals = zeros(m,p^n)
    funvals = zeros(m,p^n,funlength)

	# loop through block indices
	for idx2 = 1:m
		blockidx2 = (idx2-1)*(p-1)+1 : idx2*(p-1)
		di = d[blockidx2]
		Wdi = Rd[blockidx2]

		# loop through all possible blocks
		for idx3 = 1:p^n
			blockidx3 = (idx3-1)*(p-1)+1 : idx3*(p-1)
			r = Mrows[blockidx3[1]:blockidx3[end],:]
			logposteriorvals[idx2,idx3] = -1/2 * norm(Wdi.*(r*w-di))^2 - lambda * sum(r)
			funvals[idx2,idx3,:] = rowfuns[idx2](r,w)[:]

		end

	end

    maxlogs = maximum(logposteriorvals,2)
    posteriorvals = exp.(logposteriorvals .- maxlogs)
    funposteriorvals = posteriorvals .* funvals

    a = sum(posteriorvals,2) # size = [m,1]
    b = squeeze(sum(funposteriorvals,2),2) # size = [m,funlength]

    c = zeros(m,funlength)
    for idx2 = 1:m
        c[idx2,:] = prod(a[[1:idx2-1;idx2+1:m]]) * b[idx2,:]
    end

    return a,c,maxlogs
end


text:
using FactCheck
using PETSc
facts("  ---Checking Petsc data types---") do

  @fact PetscScalar => Complex64
  @fact PetscReal => Float32
  @fact PetscInt => Int64

end

FactCheck.exitstatus()


text:
module PairwiseMatrix

using Merlin
using HDF5
using MLDatasets

include("data.jl")
include("model.jl")
include("train.jl")

end

max_stars_repo_name: hshindo/PairwiseMatrix

max_issues_repo_name: hshindo/PairwiseMatrix

max_forks_repo_name: hshindo/PairwiseMatrix


text:
module LoLStatusV3 # LOLTools

using ..LOLTools: AbstractDTO, lol_api_server, http_action, call_api, nothing_in_event
using HTTP

struct ShardStatus <: AbstractDTO
    name
    platform_tag
    hostname
    services
    slug
    locales
    ShardStatus(name, platform_tag, hostname, services, slug, locales) = new(name, platform_tag, hostname, services, slug, locales)
end

"""
         shard_data(api_key::String,
                    platform::String ;
                    endpoint::HTTP.URI = lol_api_server(platform),
                    action::Function = http_action,
                    event::Function = nothing_in_event)::ShardStatus
"""
function shard_data(api_key::String,
                    platform::String ;
                    endpoint::HTTP.URI = lol_api_server(platform),
                    action::Function = http_action,
                    event::Function = nothing_in_event)::ShardStatus
    call_api(ShardStatus, api_key, action, endpoint, "/lol/status/v3/shard-data", event, shard_data)
end

end # module LOLTools.LoLStatusV3