text:
"""
    ConvertNotFound

Conversion is not defined via [`convert_to_option`](@ref). One should let the conversion fallback
to `Base.convert` when see this.
"""
struct ConvertNotFound end

"""
    convert_to_option(::Type{OptionType}, ::Type{ValueType}, x) where {OptionType, ValueType}

Convert `x` to type `ValueType` for option type `OptionType`. This is similar to `Base.convert`,
but will not error if the conversion is not overloaded, and will return a [`ConvertNotFound`](@ref)
object, but one can use this to avoid type piracy and define contextual conversion based on option
types.

# Example

One may have different string syntax for a `Symbol`, e.g

```julia
@option struct StringIsSymbol
    name::Symbol
end

@option struct SymbolNeedsColon
    name::Symbol
end

Configurations.convert_to_option(::Type{StringIsSymbol}, ::Type{Symbol}, x::String) = Symbol(x)

function Configurations.convert_to_option(::Type{SymbolNeedsColon}, ::Type{Symbol}, x::String)
    if startswith(x, ':')
        Symbol(x[2:end])
    else
        error("expect a Symbol, got String")
    end
end
```

then if we run it, we will have different behaviour by context.

```julia
julia> from_dict(StringIsSymbol, d)
StringIsSymbol(:ccc)

julia> from_dict(SymbolNeedsColon, d)
ERROR: expect a Symbol, got String
```
"""
convert_to_option(::Type, ::Type{T}, x) where T = convert_to_option(T, x)
convert_to_option(::Type{T}, x) where T = ConvertNotFound()


text:
# Interpolated Factored Green Function (IFGF) Method interface

wavenumber(p::AbstractPDE)    = abstractmethod(p)
wavenumber(::Laplace)         = nothing
wavenumber(p::Helmholtz)      = p.k
wavenumber(::Elastostatic)    = nothing
wavenumber(p::Maxwell)        = p.k
wavenumber(k::AbstractKernel) = wavenumber(pde(k))

function IFGF.IFGFOperator(iop::IntegralOperator;
                           p=(3,5,5),
                           nmax=100,
                           _splitter=IFGF.DyadicSplitter,
                           _profile=false)
    K = kernel(iop)
    Ypts = iop |> source_surface |> dofs
    Xpts = iop |> target_surface |> dofs
    splitter = _splitter(;nmax)
    p_func = (node) -> p
    k = wavenumber(K)
    ds_func = IFGF.cone_domain_size_func(k)
    ifgf = IFGF.IFGFOperator(K,Ypts,Xpts;splitter,p_func,ds_func,_profile)
    # wrap ifgf and qweights together as a DiscreteOp
    # (kernels do not include qweights)
    w = iop |> source_surface |> qweights |> collect |> Diagonal
    op = DiscreteOp(ifgf)*DiscreteOp(w)
    return op
end

function IFGFCompressor(;p=(3,5,5),
                        nmax=100,
                        _splitter=IFGF.DyadicSplitter,
                        _profile=false)
    return (iop::IntegralOperator) -> IFGFOperator(iop;p,nmax,_splitter,_profile)
end

# AbstractKernel
IFGF.centered_factor(K::AbstractKernel,x,y::IFGF.SourceTree) = abstractmethod(K)

# Laplace
function IFGF.centered_factor(K::SingleLayerKernel{T,<:Laplace},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    return K(x,yc)
end
function IFGF.centered_factor(::DoubleLayerKernel{T,<:Laplace},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    # TODO: pick a better centered_factor
    r = coords(x)-yc
    d = norm(r)
    return 1/d
end

# Helmholtz
function IFGF.centered_factor(K::SingleLayerKernel{T,<:Helmholtz{3}},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    # return exp(im*k*d)/d
    return K(x,yc)
end
function IFGF.centered_factor(K::DoubleLayerKernel{T,<:Helmholtz{3}},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    k = wavenumber(K)
    r = coords(x)-yc
    d = norm(r)
    return exp(im*k*d)/d*(-im*k+1/d)
end

# Elastostatic
function IFGF.centered_factor(::SingleLayerKernel{T,<:Elastostatic{3}},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    r = coords(x)-yc
    d = norm(r)
    return 1/d
end
function IFGF.centered_factor(::DoubleLayerKernel{T,<:Elastostatic{3}},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    # TODO: pick a better centered_factor
    r = coords(x)-yc
    d = norm(r)
    return 1/d
end

# Maxwell
function IFGF.centered_factor(K::SingleLayerKernel{T,<:Maxwell},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    # TODO: pick a better centered_factor
    k = wavenumber(K)
    r = coords(x)-yc
    d = norm(r)
    g = exp(im*k*d)/(4π*d)
    # gp  = im*k*g - g/d
    # gpp = im*k*gp - gp/d + g/d^2
    # RRT = rvec*transpose(rvec) # rvec ⊗ rvecᵗ
    return g
end
function IFGF.centered_factor(K::DoubleLayerKernel{T,<:Maxwell},x,y::IFGF.SourceTree) where T
    yc = IFGF.center(y)
    # TODO: pick a better centered_factor
    k = wavenumber(K)
    r = coords(x)-yc
    d = norm(r)
    g = exp(im*k*d)/(4π*d)
    return g
end

ext: jl

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using Test

using CompScienceMeshes
using BEAST
using StaticArrays
using LinearAlgebra

c = 3e8
μ = 4*π*1e-7
ε = 1/(μ*c^2)
f = 5e7
λ = c/f
k = 2*π/λ
ω = k*c
η = sqrt(μ/ε)

a = 1
Γ_orig = CompScienceMeshes.meshcuboid(a,a,a,0.1)
Γ = translate(Γ_orig,SVector(-a/2,-a/2,-a/2))

Φ, Θ = [0.0], range(0,stop=π,length=100)
pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for ϕ in Φ for θ in Θ]

# This is an electric dipole
# The pre-factor (1/ε) is used to resemble
# (9.18) in Jackson's Classical Electrodynamics
E = (1/ε) * dipolemw3d(location=SVector(0.4,0.2,0),
                       orientation=1e-9.*SVector(0.5,0.5,0),
                       wavenumber=k)

n = BEAST.NormalVector()

𝒆 = (n × E) × n
H = (-1/(im*μ*ω))*curl(E)
𝒉 = (n × H) × n

𝓣 = Maxwell3D.singlelayer(wavenumber=k)
𝓝 = BEAST.NCross()
𝓚 = Maxwell3D.doublelayer(wavenumber=k)

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

T = Matrix(assemble(𝓣,X,X))
e = Vector(assemble(𝒆,X))
j_EFIE = T\e

nf_E_EFIE = potential(MWSingleLayerField3D(wavenumber=k), pts, j_EFIE, X)
nf_H_EFIE = potential(BEAST.MWDoubleLayerField3D(wavenumber=k), pts, j_EFIE, X) ./ η
ff_E_EFIE = potential(MWFarField3D(wavenumber=k), pts, j_EFIE, X)

@test norm(nf_E_EFIE - E.(pts))/norm(E.(pts)) ≈ 0 atol=0.01
@test norm(nf_H_EFIE - H.(pts))/norm(H.(pts)) ≈ 0 atol=0.01
@test norm(ff_E_EFIE - E.(pts, isfarfield=true))/norm(E.(pts, isfarfield=true)) ≈ 0 atol=0.001

K_bc = Matrix(assemble(𝓚,Y,X))
G_nxbc_rt = Matrix(assemble(𝓝,Y,X))
h_bc = η*Vector(assemble(𝒉,Y))
M_bc = -0.5*G_nxbc_rt + K_bc
j_BCMFIE = M_bc\h_bc

nf_E_BCMFIE = potential(MWSingleLayerField3D(wavenumber=k), pts, j_BCMFIE, X)
nf_H_BCMFIE = potential(BEAST.MWDoubleLayerField3D(wavenumber=k), pts, j_BCMFIE, X) ./ η
ff_E_BCMFIE = potential(MWFarField3D(wavenumber=k), pts, j_BCMFIE, X)

@test norm(nf_E_BCMFIE - E.(pts))/norm(E.(pts)) ≈ 0 atol=0.01
@test norm(nf_H_BCMFIE - H.(pts))/norm(H.(pts)) ≈ 0 atol=0.01
@test norm(ff_E_BCMFIE - E.(pts, isfarfield=true))/norm(E.(pts, isfarfield=true)) ≈ 0 atol=0.01

H = dipolemw3d(location=SVector(0.0,0.0,0.3),
               orientation=1e-9.*SVector(0.5,0.5,0),
               wavenumber=k)

𝒉 = (n × H) × n
E = (1/(im*ε*ω))*curl(H)
𝒆 = (n × E) × n

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

T = Matrix(assemble(𝓣,X,X))
e = Vector(assemble(𝒆,X))
j_EFIE = T\e

nf_E_EFIE = potential(MWSingleLayerField3D(wavenumber=k), pts, j_EFIE, X)
nf_H_EFIE = potential(BEAST.MWDoubleLayerField3D(wavenumber=k), pts, j_EFIE, X) ./ η
ff_E_EFIE = potential(MWFarField3D(wavenumber=k), pts, j_EFIE, X)

@test norm(nf_E_EFIE - E.(pts))/norm(E.(pts)) ≈ 0 atol=0.01
@test norm(nf_H_EFIE - H.(pts))/norm(H.(pts)) ≈ 0 atol=0.01
@test norm(ff_E_EFIE - E.(pts, isfarfield=true))/norm(E.(pts, isfarfield=true)) ≈ 0 atol=0.01

K_bc = Matrix(assemble(𝓚,Y,X))
G_nxbc_rt = Matrix(assemble(𝓝,Y,X))
h_bc = η*Vector(assemble(𝒉,Y))
M_bc = -0.5*G_nxbc_rt + K_bc
j_BCMFIE = M_bc\h_bc

nf_E_BCMFIE = potential(MWSingleLayerField3D(wavenumber=k), pts, j_BCMFIE, X)
nf_H_BCMFIE = potential(BEAST.MWDoubleLayerField3D(wavenumber=k), pts, j_BCMFIE, X) ./ η
ff_E_BCMFIE = potential(MWFarField3D(wavenumber=k), pts, j_BCMFIE, X)

@test norm(j_BCMFIE - j_EFIE)/norm(j_EFIE) ≈ 0 atol=0.02
@test norm(nf_E_BCMFIE - E.(pts))/norm(E.(pts)) ≈ 0 atol=0.01
@test norm(nf_H_BCMFIE - H.(pts))/norm(H.(pts)) ≈ 0 atol=0.01
@test norm(ff_E_BCMFIE - E.(pts, isfarfield=true))/norm(E.(pts, isfarfield=true)) ≈ 0 atol=0.01


text:
using Markowitz
using LinearAlgebra
import Plots

assets = [ "Bonds - US Government"
           "Bonds - US Corporate"
           "Bonds - International"
           "Bonds - High Yield"
           "Bonds - Bank Loans"
           "Bonds - Emerging USD Debt"
           "Bonds - Emerging Local Debt"
           "Alternative - Emerging FX"
           "Alternative - Commodities"
           "Alternative - REITs"
           "Stocks - US Large"
           "Stocks - US Small"
           "Stocks - International"
           "Stocks - Emerging" ]

V = [ 133.0  39.0  36.0  -17.0  -28.0   31.0   -5.0   -6.0  -34.0  -10.0  -46.0  -68.0  -52.0  -63.0
       39.0  16.0  17.0   10.0    0.0   22.0   13.0    7.0    1.0   19.0    3.0    0.0    5.0    9.0
       36.0  17.0  38.0   19.0    3.0   34.0   45.0   29.0   37.0   48.0   20.0   21.0   42.0   50.0
      -17.0  10.0  19.0   98.0   69.0   65.0   73.0   43.0   68.0  151.0   99.0  130.0  121.0  165.0
      -28.0   0.0   3.0   69.0   84.0   34.0   42.0   25.0   55.0   99.0   64.0   93.0   87.0  119.0
       31.0  22.0  34.0   65.0   34.0   93.0   83.0   48.0   57.0  115.0   76.0   85.0   99.0  145.0
       -5.0  13.0  45.0   73.0   42.0   83.0  154.0   87.0  111.0  168.0  113.0  143.0  165.0  225.0
       -6.0   7.0  29.0   43.0   25.0   48.0   87.0   57.0   75.0   95.0   71.0   89.0  104.0  142.0
      -34.0   1.0  37.0   68.0   55.0   57.0  111.0   75.0  286.0  117.0  101.0  125.0  164.0  239.0
      -10.0  19.0  48.0  151.0   99.0  115.0  168.0   95.0  117.0  491.0  231.0  327.0  259.0  322.0
      -46.0   3.0  20.0   99.0   64.0   76.0  113.0   71.0  101.0  231.0  214.0  256.0  217.0  269.0
      -68.0   0.0  21.0  130.0   93.0   85.0  143.0   89.0  125.0  327.0  256.0  380.0  265.0  342.0
      -52.0   5.0  42.0  121.0   87.0   99.0  165.0  104.0  164.0  259.0  217.0  265.0  297.0  359.0
      -63.0   9.0  50.0  165.0  119.0  145.0  225.0  142.0  239.0  322.0  269.0  342.0  359.0  556.0 ]

E = [ 0.1 0.7 0.8 2.3 2.2 1.9 5.6 5.6 2.2 1.3 0.7 -0.1 4.1 7.2 ]

class = vec([ :FI :FI :FI :FI :FI :FI :FI :ALT :ALT :ALT :EQ :EQ :EQ :EQ ])
colors = map(c-> c==:FI ? "blue" : (c==:ALT ? "green" : "red"),class)
alphas = 0.5 .+ 0.5*[ (1:7)/7 ; (1:3)/3 ; (1:4)/4 ]

function plot_frontier()
    m = f.problem
    unused = vec(sum(f.weights, dims=1).==0)
    xrange = [0, max(maximum(sqrt.(diag(m.V))),maximum(f.vol))]*1.1
    yrange = [minimum(m.E),max(maximum(m.E),maximum(f.ret))]+(maximum(m.E)-minimum(m.E))*[-1,1]*0.1
    dense = smooth(f)
    p1 = Plots.plot(dense[:,1], dense[:,2], legend=false, label="", color="orange", width=2)
    for i = 1:length(m.E)
        Plots.scatter!([sqrt.(diag(m.V)[i])], [m.E[i]], label=m.names[i],
                       markersize=4-1.5*unused[i], markercolor=colors[i], markeralpha=alphas[i],
                       markerstrokecolor=colors[i], markerstrokealpha=alphas[i])
    end
    Plots.scatter!(f.vol, f.ret, legend=false, label="", color="orange", markerstrokecolor="orange", markersize=1)
    Plots.xlabel!("Volatility")
    Plots.ylabel!("Return")
    Plots.xlims!(Tuple(xrange))
    Plots.ylims!(Tuple(yrange))
    weights = f.weights[:,.!unused]
    indices = [1]
    for i in 2:size(weights,1)
        w0 = weights[i-1,:]
        w1 = weights[i,:]
        if sum((sign.(w0).*sign.(w1)).==-1)>0
            for j in 1:length(w0)
                if sign.(w0[j]) != sign.(w1[j])
                    indices = [indices ; i-1+abs(w0[j]/(w1[j]-w0[j]))]
                end
            end
        end
        indices = [indices ; i]
    end
    extweights = zeros(0,size(weights,2))
    extret = []
    for idx in sort(indices)
        if rem(idx,1) == 0
            extweights = [ extweights ; weights[[Int(idx)],:] ]
            extret = [ extret ; f.ret[Int(idx)] ]
        else
            w0 = weights[[Int(floor(idx))],:]
            w1 = weights[[Int(ceil(idx))],:]
            r0 = f.ret[[Int(floor(idx))]]
            r1 = f.ret[[Int(ceil(idx))]]
            extweights = [ extweights ; w0+rem(idx,1)*(w1-w0) ]
            extret = [ extret ; r0+rem(idx,1)*(r1-r0) ]
        end
    end
    wrange = [ minimum(sum((extweights.<0) .* (extweights), dims=2)), maximum(sum((extweights.>0) .* (extweights), dims=2)) ]
    wrange = [ wrange[1], wrange[1] + 2.5*(wrange[2]-wrange[1])]
    ylong = [ zeros(size(extweights,1),1) cumsum((extweights.>0) .* (extweights), dims=2) ]
    yshort = [ zeros(size(extweights,1),1) cumsum((extweights.<0) .* (extweights), dims=2) ]
    p2 = Plots.plot(legend=:right, background_color_legend="#ffffff00", legendfont=Plots.font(5), grid=nothing)
    Plots.xticks!([0,0.5,1])
    for i = 1:(size(ylong,2)-1)
        Plots.plot!(Plots.Shape([ vec(ylong[:,i]); reverse(vec(ylong[:,i+1])) ], [ extret; reverse(extret)]),
                    color=colors[.!unused][i], alpha=alphas[.!unused][i], label=m.names[.!unused][i])
    end
    for i = 1:(size(yshort,2)-1)
        Plots.plot!(Plots.Shape([ vec(yshort[:,i]); reverse(vec(yshort[:,i+1])) ], [ extret; reverse(extret)]),
                    color=colors[.!unused][i],alpha=alphas[.!unused][i], label="")
    end
    Plots.xlims!(Tuple(wrange))
    Plots.ylims!(Tuple(yrange))
    Plots.plot(p1, p2, size=(800, 400))
end

m = markowitz(E, V, names=assets)
unit_sum(m) # total weight = 100%
f=frontier(m)
plot_frontier()
optimal(f) # volatility, return and weights for the minimum variance portofolio
optimal(f,4) # volatility, return and weights for the optimal portofolio with return = 4

m=markowitz(E, V, names=assets, lower=0.05, upper=0.15) # min 5%, max 15% per position
unit_sum(m)
f=frontier(m)
plot_frontier()

m=markowitz(E, V, names=assets, # asset bounds by class: stocks -10/30, bonds 0/20, alt. 0/10
            lower = -0.1 * (class .== :EQ),
            upper = 0.3 * (class .== :EQ) + 0.2 * (class .== :FI) + 0.1 * (class .== :ALT))
unit_sum(m)
add_constraint(m, 1 * (class .== :EQ), '>', 0.3) # net equity exposure between 30% and 60%
add_constraint(m, 1 * (class .== :EQ), '<', 0.6)
add_constraint(m, [1 1 0 0 0 0 0 0 0 0 0 0 0 0], '=', 0.25) # US govt + Investment Grade = 25%
f=frontier(m)
plot_frontier()


text:
using Documenter, Weave

makedocs( modules = Weave, sitename="Weave.jl",
    pages = ["index.md", "getting_started.md", "usage.md",
    "publish.md", "chunk_options.md", "notebooks.md",
    "function_index.md"]
)

include("make_examples.jl")

deploydocs(
    repo = "github.com/mpastell/Weave.jl.git",
    target = "build"
)


text:
"""
Constraint handler enforcing tree topologies.
"""
mutable struct TreeTopoHdlr <: SCIP.AbstractConstraintHandler
    scip::SCIP.Optimizer
end

"""
Constraint data, referencing the topology and binary edges variables.
"""
mutable struct TreeTopoCons <: SCIP.AbstractConstraint{TreeTopoHdlr}
    topo::Topology
    y::Array{JuMP.VariableRef}
end

function SCIP.check(ch::TreeTopoHdlr,
                    constraints::Array{Ptr{SCIP.SCIP_CONS}},
                    sol::Ptr{SCIP.SCIP_SOL},
                    checkintegrality::Bool,
                    checklprows::Bool,
                    printreason::Bool,
                    completely::Bool)
    @assert length(constraints) == 1
    cons::TreeTopoCons = SCIP.user_constraint(constraints[1])
    ysol = SCIP.sol_values(ch.scip, cons.y, sol)

    # no need to find out more details
    return (is_tree(topology_from_candsol(cons.topo, ysol))
            ? SCIP.SCIP_FEASIBLE
            : SCIP.SCIP_INFEASIBLE)
end

"""
Enforce tree topology.

Case 1: Topology is disconnected --> add no-good on all edges.
Case 2: Topology contains cycle --> add no-good on cycle edges.
"""
function enforce_tree_topo(ch::TreeTopoHdlr, cons::TreeTopoCons)
    ysol = SCIP.sol_values(ch.scip, cons.y)

    # Tree topology is fine, nothing to do.
    if is_tree(topology_from_candsol(cons.topo, ysol))
        return SCIP.SCIP_FEASIBLE
    end

    # Find cycle in topology.
    cycle = find_cycle(topology_from_candsol(cons.topo, ysol, true))

    # Case distinction
    if length(cycle) == 0
        # 1) There is no cycle. So, topology must be disconnected (forest).
        # TODO: Refactor no-good specific code for reuse.
        # TODO: Can we do better than no-good by forcing cut edges only?
        active = ysol .> 0.5
        coefs = 2.0 * active .- 1.0
        ci = @cb_constraint(ch.scip, coefs ⋅ cons.y ≤ sum(active) - 1)
    else
        # 2) There is a cycle. Let's forbid it (with anti-parallel arcs).
        arcidx = arcindex(cons.topo)
        antidx = antiparallelindex(cons.topo)
        fwd = [arcidx[e] for e in cycle]
        bwd = [antidx[e] for e in fwd]
        arcs = vcat(fwd, bwd)
        ci = @cb_constraint(
            ch.scip, sum(cons.y[a] for a in arcs) ≤ length(cycle) - 1)
    end

    return SCIP.SCIP_CONSADDED
end

function SCIP.enforce_lp_sol(ch::TreeTopoHdlr,
                             constraints::Array{Ptr{SCIP.SCIP_CONS}},
                             nusefulconss::Cint,
                             solinfeasible::Bool)
    @assert length(constraints) == 1
    cons::TreeTopoCons = SCIP.user_constraint(constraints[1])
    return enforce_tree_topo(ch, cons)
end

function SCIP.enforce_pseudo_sol(ch::TreeTopoHdlr,
                                 constraints::Array{Ptr{SCIP.SCIP_CONS}},
                                 nusefulconss::Cint,
                                 solinfeasible::Bool,
                                 objinfeasible::Bool)
    @assert length(constraints) == 1
    cons::TreeTopoCons = SCIP.user_constraint(constraints[1])
    return enforce_tree_topo(ch, cons)
end

function SCIP.lock(ch::TreeTopoHdlr,
                   constraint::Ptr{SCIP.SCIP_CONS},
                   locktype::SCIP.SCIP_LOCKTYPE,
                   nlockspos::Cint,
                   nlocksneg::Cint)
    cons::TreeTopoCons = SCIP.user_constraint(constraint)
    for y in cons.y
        # TODO: understand why lock is called during SCIPfree, after the
        # constraint should have been deleted already. Does it mean we should
        # implement CONSTRANS?
        var_::Ptr{SCIP.SCIP_VAR} = SCIP.var(ch.scip, JuMP.index(y))
        var_ != C_NULL || continue  # avoid segfault!

        SCIP.@SC SCIP.SCIPaddVarLocksType(
            ch.scip, var_, locktype, nlockspos + nlocksneg, nlockspos + nlocksneg)
    end
end


text:
#using GRUtils
#using Plots
#pyplot()
function locustest()
    lambdas = collect(360.0:10.0:830.0)
    col=cmf(CMF31)
    #hold(false)
    #fig = gcf()
    x=zeros(length(lambdas))
    y=zeros(length(lambdas))
    z=zeros(length(lambdas))
    @inbounds for w in 0.1:10.0:400.1
        @inbounds for i in 1:length(lambdas)
            ill=led_spec(360.0,830.0,1.0,lambdas[i],w)
            color = convert(xyY,ill*col)
            x[i] = color.x
            y[i] = color.y
            z[i] = color.Y
            #color = convert(DIN99o,ill*col)
            #x[i] = color.a
            #y[i] = color.b
            #z[i] = color.l

        end
        #plot3!(fig,x,y,z,":b")
        #wireframe!(fig,x, y, z)
        #hold(true)
    end
    #aspectratio(1.0)
    #viewpoint(30, -20)
    #fig
    #videofile("rotationtestYxy.mp4") do
    #    for rot = -90:1:90
    #        viewpoint(rot, 45)
    #        draw(fig)
    #    end
    #    for rot = 89:-1:-89
    #        viewpoint(rot,45)
    #        draw(fig)
    #    end
    #end
end


text:
mutable struct Image
    uri::Union{Nothing,String}
    mimeType::Union{Nothing,String}
    bufferView::Union{Nothing,Int}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Image()
        obj = new()
        obj.uri = nothing
        obj.mimeType = nothing
        obj.bufferView = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Image, sym::Symbol, x)
    if sym === :mimeType && x !== nothing
        x == "image/jpeg" || x == "image/png" ||
            throw(ArgumentError("""the image's MIME type should be "image/jpeg" or "image/png" """))
    elseif sym === :bufferView && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the bufferView that contains the image should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Image}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Image}) = (:uri, :mimeType, :bufferView, :name, :extensions, :extras)


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the contract_left(B, C) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using LinearAlgebra: ⋅, det, diag, qr
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

@testset "contract_left(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "contract_left(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "contract_left(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::Blade, C::Blade)" begin
    # --- Preparations

    # Dimension of embedding space
    test_dim = 10
    test_grade_1 = 3
    test_grade_2 = 5

    # Test vectors
    test_vector_1 = rand(test_dim)
    test_vector_2 = rand(test_dim)

    # Test bases
    test_basis_1 = rand(test_dim, test_grade_1)
    test_basis_2 = rand(test_dim, test_grade_2)

    # --- Test cases

    # grade(B) == grade(C) == 1
    B = Blade(test_vector_1)
    C = Blade(test_vector_2)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Scalar((basis(B) ⋅ basis(C)) * volume(B) * volume(C))
    @test (B < C) == B_contract_left_C

    # grade(B) == 1, grade(C) > 1
    B = Blade(test_vector_1)

    for test_grade in 5:8
        C = Blade(test_basis_2)

        F = qr(test_basis_2)
        Q = Matrix(F.Q)
        projection = Blade(Q * transpose(Q) * test_vector_1)
        expected_volume_C = prod(diag(F.R))
        expected_result = mod(grade(C), 4) < 2 ?
            expected_volume_C * dual(projection, C) :
           -expected_volume_C * dual(projection, C)

        B_contract_left_C = contract_left(B, C)
        @test B_contract_left_C ≈ expected_result
        @test (B < C) == B_contract_left_C
    end

    # grade(B) > 1, grade(C) == 1
    B = Blade(test_basis_1)
    C = Blade(test_vector_2)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # grade(B) == grade(C) > 1
    for test_grade in 5:8
        B = Blade(rand(test_dim, test_grade))
        C = Blade(rand(test_dim, test_grade))

        B_contract_left_C = contract_left(B, C)
        @test B_contract_left_C ≈ volume(B) * volume(C) *
                              det(transpose(basis(B)) * basis(C))

        @test (B < C) == B_contract_left_C
    end

    # grade(B) < grade(C), grade(B) > 1, grade(C) > 1,
    B = Blade(test_basis_1)

    for test_grade in 5:8
        test_basis_C = rand(test_dim, test_grade)
        C = Blade(test_basis_C)

        F = qr(test_basis_C)
        Q = Matrix(F.Q)
        projection = Blade(Q * transpose(Q) * test_basis_1)
        expected_volume_C = prod(diag(F.R))
        expected_result = mod(grade(C), 4) < 2 ?
            expected_volume_C * dual(projection, C) :
           -expected_volume_C * dual(projection, C)

        B_contract_left_C = contract_left(B, C)
        @test B_contract_left_C ≈ expected_result
        @test (B < C) == B_contract_left_C
    end

    # grade(B) > grade(C), grade(B) > 1, grade(C) > 1,
    B = Blade(test_basis_2)
    C = Blade(test_basis_1)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # dim(B) != dim(C)
    B = Blade(rand(test_dim, 3))
    C = Blade(rand(test_dim + 1, 4))
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

@testset "contract_left(B::Blade, C::Pseudoscalar)" begin
    # ---Preparations

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Tests

    # dim(B) == dim(C)
    test_grade = 3
    for test_dim in 10:13
        B = Blade(rand(test_dim, test_grade))
        C = Pseudoscalar(test_dim, test_value)

        B_contract_left_C = contract_left(B, C)

        expected_result = mod(test_dim, 4) < 2 ?
            test_value * dual(B) :
           -test_value * dual(B)

        @test B_contract_left_C ≈ expected_result
        @test (B < C) == B_contract_left_C
    end

    # dim(B) != dim(C)
    test_dim = 8
    B = Blade(rand(test_dim, test_grade))
    C = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

@testset "contract_left(B::Blade, C::Scalar)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    test_vector = rand(test_dim)

    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(B) == 1
    B = Blade(test_vector)
    C = Scalar(test_value)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # grade(B) > 1
    B = Blade(test_basis)
    C = Scalar(test_value)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Blade, C::One)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3
    test_vector = rand(test_dim)
    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(B) == 1
    B = Blade(test_vector)
    C = One()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # grade(B) > 1
    B = Blade(test_basis)
    C = One()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Blade, C::Zero)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3
    test_vector = rand(test_dim)
    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(B) == 1
    B = Blade(test_vector)
    C = Zero()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # grade(B) > 1
    B = Blade(test_basis)
    C = Zero()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Blade, C::Real)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    test_vector = rand(test_dim)

    test_basis = rand(test_dim, test_grade)

    # --- Tests
    #     C::Real, B::Blade

    # grade(B) == 1
    B = Blade(test_vector)
    C = test_value

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # grade(B) > 1
    B = Blade(test_basis)
    C = test_value

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Blade, C::Vector)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3

    test_vector_1 = rand(test_dim)
    test_vector_2 = rand(test_dim)

    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(B) == 1, dim(B) == length(C)
    B = Blade(test_vector_1)
    C = test_vector_2

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ test_vector_1 ⋅ test_vector_2
    @test (B < C) == B_contract_left_C

    # grade(B) == 1, dim(B) != length(C)
    B = Blade(rand(test_dim + 1))
    C = test_vector_2
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C

    # grade(B) > 1, dim(B) == length(C)
    C = test_vector_2
    for test_grade in 5:8
        B = Blade(rand(test_dim, test_grade))

        @test iszero(contract_left(B, C))
        @test iszero(B < C)
    end

    # grade(B) > 1, dim(B) != length(C)
    B = Blade(rand(test_dim + 1, 3))
    C = test_vector_2
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

# ------ B::Pseudoscalar

@testset "contract_left(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::Pseudoscalar, C::Blade)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 5

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    C = Blade(rand(test_dim, test_grade))

    # --- Tests

    # dim(B) == dim(C)
    B = Pseudoscalar(test_dim, test_value)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # dim(B) != dim(C)
    B = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch (B < C)
end

@testset "contract_left(B::Pseudoscalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    test_value_2 = get_random_value(1)  # add 1 to keep value away from 0

    # --- Tests

    # dim(B) == dim(C)
    for test_dim = 5:8
        B = Pseudoscalar(test_dim, test_value_1)
        C = Pseudoscalar(test_dim, test_value_2)

        B_contract_left_C = contract_left(B, C)

        expected_result = mod(test_dim, 4) < 2 ?
            test_value_1 * test_value_2 :
           -test_value_1 * test_value_2

        @test B_contract_left_C isa Scalar
        @test B_contract_left_C == expected_result
        @test (B < C) == B_contract_left_C
    end

    # dim(B) != dim(C)
    test_dim = 5
    B = Pseudoscalar(test_dim, test_value_1)
    C = Pseudoscalar(test_dim + 1, test_value_2)
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

@testset "contract_left(B::Pseudoscalar, C::Scalar)" begin
    test_dim = 11
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Pseudoscalar, C::One)" begin
    test_dim = 12
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = One()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Pseudoscalar, C::Zero)" begin
    test_dim = 12
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = Zero()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Pseudoscalar, C::Real)" begin
    test_dim = 12
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Pseudoscalar, C::Vector)" begin
    # --- Preparations

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Tests

    # dim(v) == dim(B)
    for test_dim in 5:8
        B = Pseudoscalar(test_dim, test_value)
        C = Vector(rand(test_dim))

        @test iszero(contract_left(B, C))
        @test iszero(B < C)
    end

    # dim(B) != dim(C)
    test_dim = 10
    B = Pseudoscalar(test_dim + 1, test_value)
    C = rand(test_dim)

    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

# ------ B::Scalar

@testset "contract_left(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::Scalar, C::Blade)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    test_vector = rand(test_dim)

    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(C) == 1
    B = Scalar(test_value)
    C = Blade(test_vector)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Blade(basis(C), volume=test_value * volume(C))
    @test (B < C) == B_contract_left_C

    # grade(C) > 1
    B = Scalar(test_value)
    C = Blade(test_basis)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Blade(basis(C), volume=test_value * volume(C))
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Scalar, C::Pseudoscalar)" begin
    test_dim = 10

    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Scalar(test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C == Pseudoscalar(test_dim, test_value_1 * test_value_2)
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Scalar, C::Scalar)" begin
    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    test_value_2 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value_2)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C isa Scalar
    @test B_contract_left_C == test_value_1 * test_value_2
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Scalar, C::One)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = One()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === B
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Scalar, C::Zero)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = Zero()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === C
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Scalar, C::Real)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Scalar(test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C isa AbstractScalar
    @test B_contract_left_C == test_value_1 * test_value_2
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Scalar, C::Vector)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Scalar(test_value)

    C = rand(5)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Blade(C, volume=norm(C) * test_value)
    @test (B < C) == B_contract_left_C
end

# ------ B::One

@testset "contract_left(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::One, C::Blade)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3
    test_vector = rand(test_dim)
    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(B) == 1
    B = One()
    C = Blade(test_vector)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C == C
    @test (B < C) == B_contract_left_C

    # grade(B) > 1
    B = One()
    C = Blade(test_basis)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C == C
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::One, C::Pseudoscalar)" begin
    B = One()

    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C == Pseudoscalar(test_dim, test_value)
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::One, C::Scalar)" begin
    B = One()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === C
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::One, C::One)" begin
    B = One()
    C = One()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === B
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::One, C::Zero)" begin
    B = One()
    C = Zero()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === C
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::One, C::Real)" begin
    B = One()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C isa AbstractScalar
    @test B_contract_left_C == test_value
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::One, C::Vector)" begin
    B = One()
    C = rand(5)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Blade(C, volume=norm(C))
    @test (B < C) == B_contract_left_C
end

# ------ B::Zero

@testset "contract_left(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::Zero, C::Blade)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3
    test_vector = rand(test_dim)
    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(B) == 1
    B = Zero()
    C = Blade(test_vector)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)

    # grade(B) > 1
    B = Zero()
    C = Blade(test_basis)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Zero, C::Pseudoscalar)" begin
    B = Zero()

    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Zero, C::Scalar)" begin
    B = Zero()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === B
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Zero, C::One)" begin
    B = Zero()
    C = One()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === B
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Zero, C::Zero)" begin
    B = Zero()
    C = Zero()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C === B
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Zero, C::Real)" begin
    B = Zero()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Zero, C::Vector)" begin
    B = Zero()
    C = rand(5)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

# ------ B::Real

@testset "contract_left(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::Real, C::Blade)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3

    # Test values
    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    test_vector = rand(test_dim)

    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(C) == 1
    B = test_value
    C = Blade(test_vector)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Blade(basis(C), volume=test_value * volume(C))
    @test (B < C) == B_contract_left_C

    # grade(C) > 1
    B = test_value
    C = Blade(test_basis)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ Blade(basis(C), volume=test_value * volume(C))
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Real, C::Pseudoscalar)" begin
    test_dim = 12
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = test_value_1

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C == Pseudoscalar(test_dim, test_value_1 * test_value_2)
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Real, C::Scalar)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = test_value_1

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C isa AbstractScalar
    @test B_contract_left_C == test_value_1 * test_value_2
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Real, C::One)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    C = One()

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C isa AbstractScalar
    @test B_contract_left_C == test_value
    @test (B < C) == B_contract_left_C
end

@testset "contract_left(B::Real, C::Zero)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    C = Zero()

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

# ------ B::Vector

@testset "contract_left(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "contract_left(B::Vector, C::Blade)" begin
    # --- Preparations

    test_dim = 10
    test_grade = 3

    test_vector_1 = rand(test_dim)
    test_vector_2 = rand(test_dim)

    test_basis = rand(test_dim, test_grade)

    # --- Tests

    # grade(C) == 1, dim(B) == length(C)
    B = test_vector_1
    C = Blade(test_vector_2)

    B_contract_left_C = contract_left(B, C)
    @test B_contract_left_C ≈ test_vector_1 ⋅ test_vector_2
    @test (B < C) == B_contract_left_C

    # grade(C) == 1, dim(B) != length(C)
    C = Blade(rand(test_dim + 1))
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C

    # grade(C) > 1, dim(B) == length(C)
    B = test_vector_1
    for test_grade in 5:8
        C = Blade(rand(test_dim, test_grade))

        B_contract_left_C = contract_left(B, C)
        @test B_contract_left_C ≈ contract_left(Blade(B), C)
        @test (B < C) == B_contract_left_C
    end

    # grade(C) > 1, dim(B) != length(C)
    B = test_vector_1
    C = Blade(rand(test_dim + 1, 3))
    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

@testset "contract_left(B::Vector, C::Pseudoscalar)" begin
    # --- Preparations

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Tests

    # dim(B) == dim(C)
    for test_dim in 5:8
        B = Vector(rand(test_dim))
        C = Pseudoscalar(test_dim, test_value)

        B_contract_C = contract_left(B, C)
        @test B_contract_C ≈ contract_left(Blade(B), C)
        @test (B< C) == B_contract_C
    end

    # dim(B) != dim(C)
    test_dim = 10
    B = Pseudoscalar(test_dim + 1, test_value)
    C = rand(test_dim)

    @test_throws DimensionMismatch contract_left(B, C)
    @test_throws DimensionMismatch B < C
end

@testset "contract_left(B::Vector, C::Scalar)" begin
    B = rand(5)

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Scalar(test_value)

    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Vector, C::One)" begin
    B = rand(5)
    C = One()
    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end

@testset "contract_left(B::Vector, C::Zero)" begin
    B = rand(5)
    C = Zero()
    @test iszero(contract_left(B, C))
    @test iszero(B < C)
end


text:
using Test

include("./Tools.jl")
@testset "Test Tools" begin
#evalpoly testing
@testset "Test evalpoly function" begin
coefs = [1.0,7.0,9.0]
x = 3.0
@test evalpoly(x,coefs) ≈ sum([x^(i-1)*coefs[i] for i = 1:length(coefs)])
end;

@testset "Test get index first greater function" begin
#getIndexFirstGreater testing
a = [-4.0,1.0,3.0,7.0]
@test getBoundingIndsSorted(2.0,a) == (2,3)
@test getBoundingIndsSorted(-5.0,a) == (1,-1)
@test getBoundingIndsSorted(8.0,a) == (length(a),-1)
end;
end;


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

using ParallelProcessingTools
using Test
using Base.Threads

@testset "getlocalvalue" begin
    @testset "ThreadLocal" begin
        tl = @inferred ThreadLocal{Float32}(undef)
        @test typeof(tl) <: ThreadLocal{Float32}
        @test length(tl.value) == nthreads()

        @test (@inferred getallvalues(ThreadLocal{Int}(threadid))) == 1:nthreads()

        tmp = 2.5
        tl = @inferred ThreadLocal(tmp)
        tmp = 0.0
        @test tl[] ≈ 2.5

        tl[] = 1.0
        @test tl[] ≈ 1.0

        tmpF = () -> 3.0

        @test get(tl) == tl[]
        @test get(tl, tmpF) ≈ 1.0
        @test get!(tmpF, tl) ≈ 1.0
        @test get!(tl, 2.0) ≈ 1.0

        tmpF = () -> "test"
        tl = @inferred ThreadLocal{String}(undef)

        @test get(tmpF, tl) == "test"
        @test get(tl, "default") == "default"

        @test get!(tmpF, tl) == "test"
        @test get(tl) == "test"

        tl = @inferred ThreadLocal{String}(undef)
        @test get!(tl, "default") == "default"
        @test get(tl) == "default"

        @test getlocalvalue(3) == 3
        @test getlocalvalue(tl) == "default"
    end
end


text:
#!/usr/bin/env julia
using Base.Test
using RdRand
using RNGTest

RNGTest.smallcrush(rdrand)


text:
import Base:lastindex, getindex
export
constant,
Variable,
cell,
set_shape,
get_dtype,
get_variable,
placeholder,
variable_scope,
gradients,
gradient_magnitude,
hessian,
constant_initializer,
glorot_normal_initializer,
random_normal_initializer,
random_uniform_initializer,
truncated_normal_initializer,
uniform_unit_scaling_initializer,
variance_scaling_initializer,
sym,
spd,
tensor,
convert_to_tensor,
hessian_vector,
TensorArray,
gradient_checkpointing,
zeros_like,
ones_like,
gradients_colocate

"""
    constant(value; kwargs...)

Constructs a non-trainable tensor from `value`.
"""
function constant(value; kwargs...)
    if isa(value, PyObject)
        return value
    end
    if isa(value, Char)
        value = string(value)
    end
    kwargs = jlargs(kwargs)
    if !(:dtype in keys(kwargs))
        kwargs[:dtype] = DTYPE[eltype(value)]
    end
    tf.constant(value; kwargs...)
end

"""
    Variable(initial_value;kwargs...)

Constructs a trainable tensor from `value`.
"""
function Variable(initial_value;kwargs...)
    kwargs = jlargs(kwargs)
    if !(:dtype in keys(kwargs))
        kwargs[:dtype] = DTYPE[eltype(initial_value)]
    end
    tf.Variable(initial_value; kwargs...)
end

"""
    cell(arr::Array, args...;kwargs...)

Construct a cell tensor.
# Example
```julia-REPL
julia> r = cell([[1.],[2.,3.]])
julia> run(sess, r[1])
1-element Array{Float32,1}:
 1.0
julia> run(sess, r[2])
2-element Array{Float32,1}:
 2.0
 3.0
```
"""
function cell(arr::Array, args...;kwargs...)
    kwargs = jlargs(kwargs)
    if !(:dtype in keys(kwargs))
        kwargs[:dtype] = DTYPE[eltype(arr[1])]
    end
    tf.ragged.constant(arr, args...;kwargs...)
end

"""
    copy(o::PyObject)

Creates a tensor that has the same value that is currently stored in a variable.

!!! note
    The output is a graph node that will have that value when evaluated. Any time you evaluate it, it will grab the current value of `o`.
"""
function Base.:copy(o::PyObject)
    return tf.identity(o)
end

function get_variable(name::String; kwargs...)
    kwargs = jlargs(kwargs)
    tf.compat.v1.get_variable(name;kwargs...)
end

"""
    get_variable(o::Union{PyObject, Bool, Array{<:Number}};
        name::Union{String, Missing} = missing,
        scope::String = "")

Creates a new variable with initial value `o`. If `name` exists, `get_variable` returns the variable instead of create a new one.
"""
function get_variable(o::Union{PyObject, Bool, Number, Array{<:Number}};
    name::Union{String, Missing} = missing,
    scope::String = "")
    local v
    o = constant(o)
    if ismissing(name)
        name = "unnamed_"*randstring(10)
    end
    variable_scope(scope) do
        v = tf.compat.v1.get_variable(name = name, initializer=o, dtype=DTYPE[get_dtype(o)])
    end
    return v
end

"""
    get_variable(dtype::Type;
    shape::Union{Array{<:Integer}, NTuple{N, <:Integer}},
    name::Union{Missing,String} = missing
    scope::String = "")

Creates a new variable with initial value `o`. If `name` exists, `get_variable` returns the variable instead of create a new one.
"""
function get_variable(dtype::Type;
     shape::Union{Array{<:Integer}, NTuple{N, <:Integer}},
     name::Union{Missing,String} = missing,
     scope::String = "") where N
    local v
    dtype = DTYPE[dtype]
    if ismissing(name)
        name = "unnamed_"*randstring(10)
    end
    variable_scope(scope) do
        v = tf.compat.v1.get_variable(name = name, shape=shape, dtype = dtype)
    end
    return v
end

get_variable(dtype::Type,
    shape::Union{Array{<:Integer}, NTuple{N, <:Integer}},
    name::Union{Missing,String} = missing,
    scope::String = "") where N = get_variable(dtype, shape=shape, name=name, scope = scope)

"""
    placeholder(dtype::Type; kwargs...)

Creates a placeholder of the type `dtype`.
# Example
```julia
a = placeholder(Float64, shape=[20,10])
b = placeholder(Float64, shape=[]) # a scalar
c = placeholder(Float64, shape=[nothing]) # a vector
```
"""
function placeholder(dtype::Type; kwargs...)
    dtype = DTYPE[dtype]
    kwargs = Dict{Any,Any}(kwargs)
    if !(:shape in keys(kwargs))
        kwargs[:shape] = []
    end
    tf.compat.v1.placeholder(dtype;kwargs...)
end

"""
    placeholder(o::Union{Number, Array, PyObject}; kwargs...)

Creates a placeholder of the same type and size as `o`. `o` is the default value.
"""
function placeholder(o::Union{Number, Array, PyObject}; kwargs...)
    o = convert_to_tensor(o; kwargs...)
    tf.compat.v1.placeholder_with_default(o, shape=size(o))
end

function variable_scope(f, name_or_scope; reuse=AUTO_REUSE, kwargs...)
    @pywith tf.variable_scope(name_or_scope;reuse = reuse, kwargs...) begin
        f()
    end
end

function set_shape(o::PyObject, shape)
    o.set_shape(shape)
    o
end

function get_dtype(o::PyObject)
    for (k,v) in DTYPE
        if occursin(string(v),string(o.dtype))
            return k
        end
    end
    return nothing
end

Base.:eltype(o::PyObject) = get_dtype(o)

@deprecate get_shape size
get_shape(o::PyObject, i::Union{Int64,Nothing}) = size(o,i)

# compatible with Julia size
function PyCall.:size(o::PyObject, i::Union{Int64, Nothing}=nothing)
    d = o.shape.dims
    if d===nothing
        return nothing
    end
    if length(d)==0
        if i!=nothing
            return 1
        else
            return ()
        end
    else
        s = o.get_shape().as_list()
        s = [isnothing(x) ? x : Int64(x) for x in s]
        if i==nothing
            return Tuple(s)
        elseif 0<i<=length(s)
            return s[i]
        else
            return 1
        end
    end
end

function PyCall.:length(o::PyObject)
    # If `o.shape.dims` is invalid, it is not a TensorFlow object.
    if !hasproperty(o, :graph)
        return PyCall.@pycheckz ccall((PyCall.@pysym :PySequence_Size), Int, (PyCall.PyPtr,), o)
    end
    if any(isnothing.(size(o)))
        return nothing
    else
        return prod(size(o))
    end
end

function Base.:ndims(o::PyObject)
    return length(size(o))
end

@doc raw"""
    gradients(ys::PyObject, xs::PyObject; kwargs...)

Computes the gradients of `ys` w.r.t `xs`.

- If `ys` is a scalar, `gradients` returns the gradients with the same shape as `xs`.
- If `ys` is a vector, `gradients` returns the Jacobian $\frac{\partial y}{\partial x}$

!!! note
    The second usage is not suggested since `ADCME` adopts reverse mode automatic differentiation.
    Although in the case `ys` is a vector and `xs` is a scalar, `gradients` cleverly uses forward mode automatic differentiation,
    it requires that the second order gradients are implemented for relevant operators.
"""
function gradients(ys::PyObject, xs::PyObject; kwargs...)
    s1 = size(ys)
    s2 = size(xs)
    kwargs = jlargs(kwargs)
    if isnothing(s1) && isnothing(s2)
        error("s1, s2 should be rank 0, 1, 2")
    end
    if length(s1)==0
        ret = tf.gradients(ys, xs; kwargs...)
        if ret != nothing
            return ret[1]
        else
            return nothing
        end
    elseif length(s1)==1 && length(s2)==0
        return gradients10(ys, xs; kwargs...)
    elseif length(s1)==1 && length(s2)==1
        return gradients11(ys, xs)
    elseif length(s1)==2 && length(s2)==0
        grad = gradients(vec(ys), xs; kwargs...)
        reshape(grad, s1...)
    else
        error("gradients: Invalid argument")
    end
end

function gradients(ys::PyObject, xs::Array{T}; kwargs...) where T <: Union{Any, PyObject}
    zs = Array{PyObject}(undef, length(xs))
    for i = 1:length(zs)
        zs[i] = gradients(ys, xs[i]; kwargs...)
    end
    zs
end

function gradients(ys::Array{T}, xs::PyObject; kwargs...) where T <: Union{Any, PyObject}
    zs = Array{PyObject}(undef, length(ys))
    for i = 1:length(ys)
        zs[i] = gradients(ys[i], xs; kwargs...)
    end
    zs
end

function gradients10(ys::PyObject, xs::PyObject; kwargs...)
    kwargs = jlargs(kwargs)
    try
        u = Variable(rand(length(ys)), trainable=false)
        g = tf.gradients(ys, xs, grad_ys=u; kwargs...)
        g==nothing && (return nothing)
        r = tf.gradients(g[1], u; unconnected_gradients="zero", kwargs...)[1]
    catch
        gradients_v(ys, xs; kwargs...)
    end
end

function gradients_v(ys::PyObject, xs::PyObject;kwargs...)
    kwargs = jlargs(kwargs)
    if length(size(ys))!=1
        error("ys should be a n dimensional vector function")
    end
    if length(size(xs))!=0
        error("xs should be a scalar")
    end
    n = length(ys)
    function condition(i, ta)
        i<=n
    end
    function body(i,ta)
        ta = write(ta, i, tf.gradients(ys[i], xs, unconnected_gradients="zero", kwargs...)[1])
        i+1, ta
    end
    ta = TensorArray(n)
    i = constant(1, dtype=Int32)
    _, ta = while_loop(condition, body, [i,ta], parallel_iterations=10)
    out = stack(ta)
end

# https://stackoverflow.com/questions/50244270/computing-jacobian-matrix-in-tensorflow
function gradients11(ys::PyObject, xs::PyObject; kwargs...)
    n = size(ys,1)
    function condition(i, ta)
        i <= n
    end
    function body(i, ta)
        ta = write(ta, i, tf.convert_to_tensor(gradients(ys[i], xs; unconnected_gradients="zero")))
        i+1, ta
    end
    ta = TensorArray(n)
    i = constant(1, dtype=Int32)
    _, ta = while_loop(condition, body, [i,ta],parallel_iterations=10)
    stack(ta)
end

"""
    gradients_colocate(loss::PyObject, xs::Union{PyObject, Array{PyObject}}, args...;use_locking::Bool = true, kwargs...)

Computes the gradients of a **scalar** loss function `loss` with respect to `xs`. The gradients are colocated with respect to the forward pass.
This function is usually in distributed computing.
"""
function gradients_colocate(loss::PyObject, xs::Union{PyObject, Array{PyObject}}, args...;use_locking::Bool = true, kwargs...)
    flag = false
    if isa(xs, PyObject)
        xs = [xs]
        flag = true
    end
    opt = tf.train.Optimizer(use_locking, "default_optimizer_"*randstring(8))
    grads_and_vars = opt.compute_gradients(
        loss, var_list=xs,
        colocate_gradients_with_ops=true)
    grads = [x[1] for x in grads_and_vars]
    flag && (grads = grads[1])
    return grads
end

function hessian_vector(f, xs, v; kwargs...)
    kwargs = jlargs(kwargs)
    gradients_impl = pyimport("tensorflow.python.ops.gradients_impl")
    gradients_impl._hessian_vector_product(f, [xs], [v]; kwargs...)[1]
end

"""
`hessian` computes the hessian of a scalar function f with respect to vector inputs xs
"""
function hessian(ys::PyObject, xs::PyObject; kwargs...)
    kwargs = jlargs(kwargs)
    s1 = size(ys)
    s2 = size(xs)
    if s1==nothing || s2 == nothing || (length(s1)!=0 && length(s2)!=1)
        error("Invalid input arguments")
    end
    h = tf.gradients(ys, xs; kwargs...)
    if h==nothing
        return nothing
    else
        h = h[1]
    end
    vs = Array{PyObject}(undef, s2[1])
    for i = 1:s2[1]
        # verbose && (@info "_hessian... $i/$(s2[1])")
        vs[i] = gradients(get(h,i-1), xs)
        if isnothing(vs[i])
            vs[i] = constant(zeros(get_dtype(ys), s2[1]))
        end
    end
    stack(vs, dims=1)
end

# initializers
constant_initializer(args...;kwargs...) = tf.constant_initializer(args...;kwargs...)
glorot_normal_initializer(args...;kwargs...) = tf.glorot_normal_initializer(args...;kwargs...)
glorot_uniform_initializer(args...;kwargs...) = tf.glorot_uniform_initializer(args...;kwargs...)
random_normal_initializer(args...;kwargs...) = tf.random_normal_initializer(args...;kwargs...)
random_uniform_initializer(args...;kwargs...) = tf.random_uniform_initializer(args...;kwargs...)
truncated_normal_initializer(args...;kwargs...) = tf.truncated_normal_initializer(args...;kwargs...)
uniform_unit_scaling_initializer(args...;kwargs...) = tf.uniform_unit_scaling_initializer(args...;kwargs...)
variance_scaling_initializer(args...;kwargs...) = tf.uniform_unit_scaling_initializer(args...;kwargs...)

#--------------------------------------------------------------------------------------------------------
# Indexing

function PyCall.:lastindex(o::PyObject)
    return size(o,1)
end

function PyCall.:lastindex(o::PyObject, i::Int64)
    return size(o,i)
end

# rank 1 tensor
function getindex(o::PyObject, r::Union{Colon, Array{Bool,1}, BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}})
    if typeof(r)==Colon
        return vec(o)
    elseif typeof(r)==Array{Bool,1} || typeof(r)==BitArray{1}
        return getindex(o, findall(r))
    elseif typeof(r)==UnitRange{Int64} || typeof(r)==StepRange{Int64, Int64}
        return getindex(o, collect(r))
    elseif typeof(r)==Array{Int64,1}
        return tf.gather(o, r.-1)
    end
end

function getindex(o::PyObject, i::PyObject)
    return tf.gather(o, i-1)
end

# rank 2 tensor
# drawback: only access the 1st dimension of a tensor
function getindex(o::PyObject, i1::Union{Int64, Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}}, c::Colon)
    if typeof(i1)==Colon
        return o
    elseif typeof(i1)==Array{Bool,1}|| typeof(i1)==BitArray{1}
        return getindex(o, findall(i1), c)
    elseif typeof(i1)==UnitRange{Int64} || typeof(i1)==StepRange{Int64, Int64}
        i1 = (i1|>collect)
    end
    tf.gather(o, i1.-1)
end

# less efficient
function _to_range_array(o::PyObject,
    i::Union{Int64, Colon, Array{Bool,1},BitArray{1},  UnitRange{Int64}, Array{Int64,1}, StepRange{Int64, Int64}},
    d::Int64)
    if typeof(i)==Int64
        return [i]
    elseif typeof(i)==Colon
        if isnothing(size(o,d))
            error(ArgumentError("Dimension $d of `o` is `nothing`. You need to set a concrete dimension, e.g., using `set_shape`"))
        end
        return collect(1:size(o,d))
    elseif typeof(i)<:StepRange || typeof(i)<:UnitRange
        return collect(i)
    elseif typeof(i)==Array{Bool,1}|| typeof(i)==BitArray{1}
        return findall(i)
    else
        return i
    end
end

function getindex(o::PyObject, i1::Union{Int64, Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}},
    i2::Union{Int64, Array{Bool,1},Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}})
    sdim1 = typeof(i1)==Int64
    sdim2 = typeof(i2)==Int64
    if typeof(i1)==Int64 && typeof(i2)==Int64
        return squeeze(tf.strided_slice(o, (i1[1]-1,i2[1]-1),(i1[1],i2[1]),(1,1)), dims=(1,2))
    end
    i1 = _to_range_array(o, i1, 1)
    i2 = _to_range_array(o, i2, 2)
    temp = Base.Iterators.product(i2,i1)|>collect
    indices = [vec([x[2] for x in temp]) vec([x[1] for x in temp])] .- 1
    p = tf.gather_nd(o, indices)
    p = tf.reshape(p', (length(i1), length(i2)))

    if sdim1
        return squeeze(p, dims=1)
    elseif sdim2
        return squeeze(p, dims=2)
    else
        return p
    end
end

function getindex(o::PyObject, i1, i2, i3)
    sdims = Int64[]
    if isa(i1, Int64)
         push!(sdims, 1)
         i1 = i1:i1
    end
    if isa(i2, Int64)
        push!(sdims, 2)
        i2 = i2:i2
   end
    if isa(i3, Int64)
        push!(sdims, 3)
        i3 = i3:i3
    end
    res = getindex(o, i1, i2, i3)
    squeeze(res, dims=sdims)
end

function getindex(o::PyObject,
    i1::Union{Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}},
    i2::Union{Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}},
    i3::Union{Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}})
    if length(size(o))!=3
        error(DimensionMismatch("input tensor is $(length(size(o))) dimensional, but expected 3"))
    end
    i1 = _to_range_array(o, i1, 1)
    i2 = _to_range_array(o, i2, 2)
    i3 = _to_range_array(o, i3, 3)
    idx = Int64[]
    for i = 1:length(i1)
        for j = 1:length(i2)
            for k = 1:length(i3)
                ii = i3[k] + (i2[j]-1)*size(o,3) + (i1[i]-1)*size(o,2)*size(o,3)
                push!(idx, ii)
            end
        end
    end
    p = reshape(o, (-1,))[idx]
    reshape(p, (length(i1), length(i2), length(i3)))
end

function Base.:getindex(o::PyObject, i::PyObject, j::Union{Int64, Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}})
    flag = false
    if isa(j, Colon)
        return o[i]
    end
    if length(size(i))!=0 || !(get_dtype(i)<:Integer)
        error(ArgumentError("Only integer `i` is supported"))
    end
    if isnothing(size(o,2))
        error(ArgumentError("Dimension 2 of `o` is `nothing`. You need to set a concrete dimension, e.g., using `set_shape`"))
    end
    if isa(j, Integer)
        flag = true
    end
    j = _to_range_array(o, j, 2)
    idx = (i-1)*size(o,2) + j
    ret = tf.reshape(o, (-1,))[idx]
    if flag
        return get(ret, 0)
    else
        return ret
    end
end

function Base.:getindex(o::PyObject, i::Union{Int64, Colon, Array{Bool,1},BitArray{1}, Array{Int64,1},UnitRange{Int64}, StepRange{Int64, Int64}}, j::PyObject)
    flag = false
    if length(size(j))!=0 || !(get_dtype(j)<:Integer)
        error(ArgumentError("Only integer `j` is supported"))
    end
    if isnothing(size(o,2))
        error(ArgumentError("Dimension 2 of `o` is `nothing`. You need to set a concrete dimension, e.g., using `set_shape`"))
    end
    if isa(i, Integer)
        flag = true
    end
    i = _to_range_array(o, i, 1)
    idx = (i .- 1)*size(o,2) + j
    ret = tf.reshape(o, (-1,))[idx]
    if flag
        return get(ret, 0)
    else
        return ret
    end
end

#-------------------------------------------------------------------------------------------------------

# https://stackoverflow.com/questions/46718356/tensorflow-symmetric-matrix
function sym(o::Union{Array{<:Real}, PyObject})
    convert_to_tensor(1/2 * (o + o'))
end

function spd(o::Union{Array{<:Real}, PyObject})
    if length(size(o))!=2 || size(o,1)!=size(o,2)
        error("Input `o` must be a square matrix")
    end
    convert_to_tensor(o * o')
end

"""
    tensor(v::Array{T,2}; dtype=Float64, sparse=false) where T
"""
function tensor(v::Array{T,1}; dtype=Float64, sparse=false) where T
    local ret
    N = length(v)
    ret = Variable(zeros(dtype, N), trainable=false)
    for i = 1:N
        if sparse && isa(v[i], Number) && v[i]≈0
            continue
        end
        if isa(v[i], Number)
            v[i] = dtype(v[i])
        end
        ret = scatter_add(ret, i, v[i])
    end
    ret
end

"""
    tensor(v::Array{T,2}; dtype=Float64, sparse=false) where T

Convert a generic array `v` to a tensor. For example,
```julia
v = [0.0 constant(1.0) 2.0
    constant(2.0) 0.0 1.0]
u = tensor(v)
```
`u` will be a ``2\\times 3`` tensor.
!!! note
    This function is expensive. Use with caution.
"""
function tensor(v::Array{T,2}; dtype=Float64, sparse=false) where T
    local ret
    M, N = size(v)
    ret = Variable(zeros(dtype, M, N), trainable=false)
    for i = 1:M
        for j = 1:N
            if sparse && isa(v[i,j], Number) && v[i,j]≈0
                continue
            end
            if isa(v[i,j], Number)
                v[i,j] = dtype(v[i,j])
            end
            ret = scatter_add(ret, i, j, v[i, j])
        end
    end
    ret
end

"""
    TensorArray(size_::Int64=0, args...;kwargs...)

Constructs a tensor array for [`while_loop`](@ref).
"""
function TensorArray(size_::Int64=0, args...;kwargs...)
    kwargs = jlargs(kwargs)
    if !(haskey(kwargs, :dtype))
        kwargs[:dtype] = tf.float64
    end
    if !haskey(kwargs, :dynamic_size)
        kwargs[:dynamic_size] = false
    end
    if !haskey(kwargs, :clear_after_read)
        kwargs[:clear_after_read] = false
    end
    kwargs[:size] = size_
    tf.TensorArray(args...;kwargs...)
end

"""
    read(ta::PyObject, i::Union{PyObject,Integer})

Reads data from [`TensorArray`](@ref) at index `i`.
"""
function Base.:read(ta::PyObject, i::Union{PyObject,Integer})
    ta.read(i-1)
end

"""
    write(ta::PyObject, i::Union{PyObject,Integer}, obj)

Writes data `obj` to [`TensorArray`](@ref) at index `i`.
"""
function Base.:write(ta::PyObject, i::Union{PyObject,Integer}, obj::PyObject)
    ta.write(i-1, obj)
end

Base.:write(ta::PyObject, i::Union{PyObject,Integer}, obj::Union{Array{<:Real}, Real}) = write(ta, i, constant(obj))

"""
    convert_to_tensor(o::Union{PyObject, Number, Array{T}, Missing, Nothing}; dtype::Union{Type, Missing}=missing) where T<:Number
    convert_to_tensor(os::Array, dtypes::Array)

Converts the input `o` to tensor. If `o` is already a tensor and `dtype` (if provided) is the same as that of `o`, the operator does nothing.
Otherwise, `convert_to_tensor` converts the numerical array to a constant tensor or casts the data type.
`convert_to_tensor` also accepts multiple tensors.

# Example
```julia
convert_to_tensor([1.0, constant(rand(2)), rand(10)], [Float32, Float64, Float32])
```
"""
function convert_to_tensor(o::Union{PyObject, Number, Array{T}, Missing, Nothing};
    dtype::Union{Type, Missing}=missing) where T<:Number
    if ismissing(o) || isnothing(o)
        return o
    end
    if isa(o, PyObject)
        if ismissing(dtype) || dtype==eltype(o)
            return o
        else
            return cast(o, dtype)
        end
    else
        if !ismissing(dtype)
            return constant(o, dtype=dtype)
        else
            return constant(o)
        end
    end
end

function convert_to_tensor(os::Array, dtypes::Array)
    [convert_to_tensor(o, dtype=d) for (o, d) in zip(os, dtypes)]
end

"""
    gradient_checkpointing(type::String="speed")

Uses checkpointing scheme for gradients.
- 'speed':  checkpoint all outputs of convolutions and matmuls. these ops are usually the most expensive,
    so checkpointing them maximizes the running speed
    (this is a good option if nonlinearities, concats, batchnorms, etc are taking up a lot of memory)
- 'memory': try to minimize the memory usage
    (currently using a very simple strategy that identifies a number of bottleneck tensors in the graph to checkpoint)
- 'collection': look for a tensorflow collection named 'checkpoints', which holds the tensors to checkpoint
"""
function gradient_checkpointing(type::String="speed")
    pyfile = "$(ADCME.LIBDIR)/memory_saving_gradients.py"
    if !isfile(pyfile)
        @info "Downloading memory_saving_gradients.py..."
        download("https://raw.githubusercontent.com/cybertronai/gradient-checkpointing/master/memory_saving_gradients.py",
                pyfile)
    end
    py"""exec(open($pyfile).read())"""
    gradients_speed = py"gradients_speed"
    gradients_memory = py"gradients_memory"
    gradients_collection = py"gradients_collection"
    if type=="speed"
        tf.__dict__["gradients"] = gradients_speed
    elseif type=="memory"
        tf.__dict__["gradients"] = gradients_memory
    elseif type=="collection"
        tf.__dict__["gradients"] = gradients_collection
    else
        error("ADCME: $type not defined")
    end
    @info "Loaded: $type"
end

@doc raw"""
    gradient_magnitude(l::PyObject, o::Union{Array, PyObject})

Returns the gradient sum

$$\sqrt{\sum_{i=1}^n \|\frac{\partial l}{\partial o_i}\|^2}$$

This function is useful for debugging the training
"""
function gradient_magnitude(l::PyObject, o::Union{Array, PyObject})
    if isa(o, PyObject)
        o = [o]
    end
    tf.global_norm(gradients(l, o))
end

"""
    zeros_like(o::Union{PyObject,Real, Array{<:Real}}, args...; kwargs...)

Returns a all-zero tensor, which has the same size as `o`.

# Example
```julia
a = rand(100,10)
b = zeros_like(a)
@assert run(sess, b)≈zeros(100,10)
```
"""
function zeros_like(o::Union{PyObject,Real, Array{<:Real}}, args...; kwargs...)
    kwargs = jlargs(kwargs)
    tf.zeros_like(o, args...;kwargs...)
end

"""
    ones_like(o::Union{PyObject,Real, Array{<:Real}}, args...; kwargs...)

Returns a all-one tensor, which has the same size as `o`.

# Example
```julia
a = rand(100,10)
b = ones_like(a)
@assert run(sess, b)≈ones(100,10)
```
"""
function ones_like(o::Union{PyObject,Real, Array{<:Real}}, args...; kwargs...)
    kwargs = jlargs(kwargs)
    tf.ones_like(o, args...;kwargs...)
end


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

mutable struct IoK8sApiNodeV1alpha1Scheduling <: SwaggerModel
    nodeSelector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: nodeSelector
    tolerations::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Toleration} } # spec name: tolerations

    function IoK8sApiNodeV1alpha1Scheduling(;nodeSelector=nothing, tolerations=nothing)
        o = new()
        validate_property(IoK8sApiNodeV1alpha1Scheduling, Symbol("nodeSelector"), nodeSelector)
        setfield!(o, Symbol("nodeSelector"), nodeSelector)
        validate_property(IoK8sApiNodeV1alpha1Scheduling, Symbol("tolerations"), tolerations)
        setfield!(o, Symbol("tolerations"), tolerations)
        o
    end
end # type IoK8sApiNodeV1alpha1Scheduling

const _property_map_IoK8sApiNodeV1alpha1Scheduling = Dict{Symbol,Symbol}(Symbol("nodeSelector")=>Symbol("nodeSelector"), Symbol("tolerations")=>Symbol("tolerations"))
const _property_types_IoK8sApiNodeV1alpha1Scheduling = Dict{Symbol,String}(Symbol("nodeSelector")=>"Dict{String, String}", Symbol("tolerations")=>"Vector{IoK8sApiCoreV1Toleration}")
Base.propertynames(::Type{ IoK8sApiNodeV1alpha1Scheduling }) = collect(keys(_property_map_IoK8sApiNodeV1alpha1Scheduling))
Swagger.property_type(::Type{ IoK8sApiNodeV1alpha1Scheduling }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNodeV1alpha1Scheduling[name]))}
Swagger.field_name(::Type{ IoK8sApiNodeV1alpha1Scheduling }, property_name::Symbol) =  _property_map_IoK8sApiNodeV1alpha1Scheduling[property_name]

function check_required(o::IoK8sApiNodeV1alpha1Scheduling)
    true
end

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


text:
immutable Cauchy <: ContinuousUnivariateDistribution
    μ::Float64
    β::Float64

    function Cauchy(μ::Real, β::Real)
        β > zero(β) || error("Cauchy: scale must be positive")
        new(float64(μ), float64(β))
    end

    Cauchy(μ::Real) = new(float64(μ), 1.0)
    Cauchy() = new(0.0, 1.0)
end

@distr_support Cauchy -Inf Inf

#### Parameters

location(d::Cauchy) = d.μ
scale(d::Cauchy) = d.β

params(d::Cauchy) = (d.μ, d.β)

#### Statistics

mean(d::Cauchy) = NaN
median(d::Cauchy) = location(d)
mode(d::Cauchy) = location(d)

var(d::Cauchy) = NaN
skewness(d::Cauchy) = NaN
kurtosis(d::Cauchy) = NaN

entropy(d::Cauchy) = log(scale(d)) + log4π

#### Functions

zval(d::Cauchy, x::Float64) = (x - d.μ) / d.β
xval(d::Cauchy, z::Float64) = d.μ + z * d.β

pdf(d::Cauchy, x::Float64) = 1.0 / (π * scale(d) * (1 + zval(d, x)^2))
logpdf(d::Cauchy, x::Float64) = - (logπ + log(scale(d)) + log1psq(zval(d, x)))

function cdf(d::Cauchy, x::Float64)
    μ, β = params(d)
    invπ * atan2(x - μ, β) + 0.5
end

function ccdf(d::Cauchy, x::Float64)
    μ, β = params(d)
    invπ * atan2(μ - x, β) + 0.5
end

function quantile(d::Cauchy, p::Float64)
    μ, β = params(d)
    μ + β * tan(π * (p - 0.5))
end

function cquantile(d::Cauchy, p::Float64)
    μ, β = params(d)
    μ + β * tan(π * (0.5 - p))
end

mgf(d::Cauchy, t::Real) = t == zero(t) ? 1.0 : NaN

cf(d::Cauchy, t::Real) = exp(im * t * d.location - d.scale * abs(t))

#### Fitting

# Note: this is not a Maximum Likelihood estimator
function fit{T <: Real}(::Type{Cauchy}, x::Array{T})
    l, u = iqr(x)
    Cauchy(median(x), (u - l) / 2.0)
end


text:
using CGMethod1D
using Documenter, Pkg

push!(LOAD_PATH,"../src/")

makedocs(
         sitename = "CGMethod1D.jl",
         modules  = [CGMethod1D],
         pages=[
                "Home" => "index.md"
               ])

deploydocs(;
    repo="github.com/cfgarciar/CGMethod1D.jl.git",
    target = "build",
    push_preview = true,
)


text:
module BSE_k_ISDF

# package code goes here
include("bse_type.jl")
include("bse_helper.jl")
include("isdf.jl")
include("1d_interface.jl")
include("exciting_interface.jl")
include("bse.jl")
include("lanczos.jl")

end # module


text:
# Cells that are commented out ( that is, starting with `#` ) are supposed to be run once, if the cells underneath are not running correctly.

## using Pkg
## pkg"add PyCall RCall"
#----------------------------------------------------------------------------

# This is loading the `PyCall` and `RCall` packages.

using PyCall,RCall;
#----------------------------------------------------------------------------

# We can inspect the names exported by a Package with the `names` function. We use the `show` function to view the full output instead of the truncated one.

show(names(RCall))
#----------------------------------------------------------------------------

# The documentation string of a function can be viewed by prefixing a `?`

?@rput
#----------------------------------------------------------------------------

# ## Using Python from Julia

# Python code can be written via the python string macro ( `@py_str` ) as shown below.

py"""
def fib(n):
    if n<2:
        return 1
    x,y = 1,1
    for i in range(n-2):
        x,y = y,x+y
    return y
"""
#----------------------------------------------------------------------------

# Python libraries can be imported in julia via the `@pyimport` macro.

@pyimport numpy
#----------------------------------------------------------------------------

##Pkg.add("PyPlot");
#----------------------------------------------------------------------------

using PyPlot ## that's using matplotlib from julia
#----------------------------------------------------------------------------

# The multiline string version ( `py""" ... """` ) does return `nothing` whereas the inline string version  ( `py" ... "` ) returns the last expression.

py"""
import numpy as np
X = np.arange(1,30)
Y = np.array([fib(x) for x in X])
ax.scatter(X, Y)
ax.set(xlabel='n', ylabel='fib(n)',
       title='The Fibonacci sequence grows fast !')
"""
X, Y = py"X", py"Y"
function fibplot( x, y )
    scatter( x, y )
    xlabel("n")
    ylabel("fib(n)")
    title("The Fibonacci sequence grows fast !")
end
fibplot( X, Y );
#----------------------------------------------------------------------------

# ## Using R from Julia
# Here we are passing the arrays that we created with python and plotted in Julia to R

@rput X Y;
#----------------------------------------------------------------------------

R"""
my_summary = summary(lm(log(Y)~X))
val <- my_summary$coefficients

plot(X, log(Y))
abline(my_summary)
""";
#----------------------------------------------------------------------------

# ## Writing, compiling and executing C from Julia

# And now we write and compile the same function in C and wrap it in Julia

## Pkg.add("Libdl")
using Libdl
C_code = """
int cfib(int n)
{
    int res=0;
    if (n <= 2){
        res = 1;
    } else {
        res = cfib(n-1)+cfib(n-2);
    }
    return res;
}
"""

const Clib = tempname()   ## make a temporary file

## compile to a shared library by piping C_code to gcc
## (works only if you have gcc installed):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

## define a Julia function that calls the C function:
c_fib(X::Int) = ccall(("cfib", Clib), Int, (Int,), X)
#----------------------------------------------------------------------------

x = 1:30
y = [c_fib(ξ) for ξ in x] ## ξ is \xi<tab>
fibplot( x, y );
#----------------------------------------------------------------------------

max_stars_repo_name: BeastyBlacksmith/julia-users-workshop

max_issues_repo_name: BeastyBlacksmith/julia-users-workshop

max_forks_repo_name: BeastyBlacksmith/julia-users-workshop


text:
module UnitConversions

export change_energy_unit,
       energy_from,
       energy_to

"""
    energy_to(e, unit)

Changes energy from hartrees to given unit

Suported units include:
  * cm^-1/cm⁻¹/cm-1
  * eV
  * hartree
  * kcal/mol
  * kJ/mol
  * K

If unit is not recognized defaults to hartree.

Source is [Wikipedia](https://en.wikipedia.org/wiki/Hartree)
"""
function energy_to(e, unit)
    if unit in ["cm-1", "cm^-1", "cm⁻¹"]
        return e*219474.6313702
    elseif unit in ["ev", "eV"]
        return e*27.21138602
    elseif unit in ["kcal/mol"]
        return e*627.509474
    elseif unit in ["kJ/mol", "kj/mol"]
        return e*2625.499639
    elseif unit in ["K"]
        return e*315775.13
    else
        return e
    end
end

"""
    energy_from(e, unit)

Changes given energy unit to hartree

Suported units include:
  * cm^-1/cm⁻¹/cm-1
  * eV
  * hartree
  * kcal/mol
  * kJ/mol
  * K

If unit is not recognized defaults to hartree.

Source is [Wikipedia](https://en.wikipedia.org/wiki/Hartree)
"""
function energy_from(e, unit)
    if unit in ["cm-1", "cm^-1", "cm⁻¹"]
        return e/219474.6313702
    elseif unit in ["ev", "eV"]
        return e/27.21138602
    elseif unit in ["kcal/mol"]
        return e/627.509474
    elseif unit in ["kJ/mol", "kj/mol"]
        return e/2625.499639
    elseif unit in ["K"]
        return e/315775.13
    else
        return e
    end
end

"""
    change_energy(e, from, to)

Changes energy unit

Suported units include:
  * cm^-1/cm⁻¹/cm-1
  * eV
  * hartree
  * kcal/mol
  * kJ/mol
  * K

If unit is not recognized defaults to hartree.

Source is [Wikipedia](https://en.wikipedia.org/wiki/Hartree)
"""
function change_energy_unit(e, from, to)
    tmp = energy_from(e, from)
    return energy_to(tmp,to)
end

end #module


text:
using Test
using DataAndPlotting

filename ="array.h5"
arr = init()

#will not overwrite existing file, therefore remove file first
# should access file in current directory which is test -
# therfore file would be C:/Julia/DataAndPlotting.jl/test/array.h5
rm("array.h5", force=true)

output(arr, "array.h5")
# for windows try - output(arr, "/temp/array.h5")

arr2 = input("array.h5")
@test isequal(arr2, arr)

plotname = "wave.png"		#gets put in current directory "C:/Julia/DataAndPlotting.jl/test"
makeplot(arr, plotname)


text:
using TemporalGPs:
    GaussMarkovModel,
    dim,
    LGSSM,
    Gaussian,
    StorageType,
    ScalarStorage,
    is_of_storage_type,
    storage_type,
    SmallOutputLGC,
    ScalarOutputLGC,
    LargeOutputLGC,
    BottleneckLGC

# Generation of positive semi-definite matrices.

function random_vector(rng::AbstractRNG, N::Int, ::ArrayStorage{T}) where {T<:Real}
    return randn(rng, T, N)
end

function random_vector(rng::AbstractRNG, N::Int, ::SArrayStorage{T}) where {T<:Real}
    return SVector{N}(randn(rng, T, N))
end

function random_matrix(rng::AbstractRNG, M::Int, N::Int, ::ArrayStorage{T}) where {T<:Real}
    return randn(rng, T, M, N)
end

function random_matrix(rng::AbstractRNG, M::Int, N::Int, ::SArrayStorage{T}) where {T<:Real}
    return SMatrix{M, N}(randn(rng, T, M, N))
end

function random_nice_psd_matrix(
    rng::AbstractRNG, N::Integer, ::Val{:dense}, ::ArrayStorage{T},
) where {T}

    # Generate random positive definite matrix.
    S = Symmetric(kernelmatrix(Matern12Kernel(), 5 .* randn(rng, T, N)) + T(1e-3) * I)

    # Centre (make eigenvals N(0, 2^2)) and bound the eigenvalues between 0 and 1.
    λ, Γ = eigen(S)
    m_λ = N > 1 ? mean(λ) : mean(λ) + T(0.1)
    σ_λ = N > 1 ? std(λ) : T(1.0)
    λ .= 2 .* (λ .- (m_λ + 0.1)) ./ σ_λ
    @. λ = T(1 / (1 + exp(-λ)) * 0.9 + 0.1)
    return collect(Symmetric(Γ * Diagonal(λ) * Γ'))
end

function random_nice_psd_matrix(
    rng::AbstractRNG, N::Integer, ::Val{:dense}, ::SArrayStorage{T},
) where {T}
    return SMatrix{N, N, T}(random_nice_psd_matrix(rng, N, Val(:dense), ArrayStorage(T)))
end

function random_nice_psd_matrix(rng::AbstractRNG, N::Integer, storage::StorageType)
    return random_nice_psd_matrix(rng, N, Val(:dense), storage)
end

function random_nice_psd_matrix(
    rng::AbstractRNG, N::Integer, ::Val{:diag}, ::ArrayStorage{T},
) where {T}
    return Diagonal(rand(T, N) .+ T(0.1))
end

function random_nice_psd_matrix(
    rng::AbstractRNG, N::Integer, ::Val{:diag}, ::SArrayStorage{T},
) where {T}
    return Diagonal(rand(SVector{N, T}) .+ T(0.1))
end

function random_nice_psd_matrix(rng::AbstractRNG, ::Integer, ::ScalarStorage{T}) where {T}
    return rand(rng, T) + convert(T, 0.1)
end

#
# Generation of Gaussians.
#

function random_gaussian(rng::AbstractRNG, dim::Int, s::StorageType)
    return Gaussian(random_vector(rng, dim, s), random_nice_psd_matrix(rng, dim, s))
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, d::T) where {T<:Gaussian}
    return Composite{T}(
        m=rand_tangent(rng, d.m),
        P=random_nice_psd_matrix(rng, length(d.m), storage_type(d)),
    )
end

# Generation of SmallOutputLGC.

function random_small_output_lgc(
    rng::AbstractRNG, Dlat::Int, Dobs::Int, Q_type::Val, s::StorageType,
)
    return SmallOutputLGC(
        random_matrix(rng, Dobs, Dlat, s),
        random_vector(rng, Dobs, s),
        random_nice_psd_matrix(rng, Dobs, Q_type, s),
    )
end

function random_scalar_output_lgc(rng::AbstractRNG, Dlat::Int, s::StorageType)
    return ScalarOutputLGC(
        random_vector(rng, Dlat, s)',
        randn(rng, eltype(s)),
        rand(rng, eltype(s)) + 0.1,
    )
end

function lgc_from_scalar_output_lgc(lgc::ScalarOutputLGC)
    return SmallOutputLGC(
        collect(lgc.A), [lgc.a], reshape([lgc.Q], 1, 1),
    )
end

function random_large_output_lgc(
    rng::AbstractRNG, Dlat::Int, Dobs::Int, Q_type::Val, s::StorageType,
)
    return LargeOutputLGC(
        random_matrix(rng, Dobs, Dlat, s),
        random_vector(rng, Dobs, s),
        random_nice_psd_matrix(rng, Dobs, Q_type, s),
    )
end

function random_bottleneck_lgc(
    rng::AbstractRNG, Dlat::Int, Dmid::Int, Dobs::Int, Q_type::Val, s::StorageType,
)
    return BottleneckLGC(
        random_matrix(rng, Dmid, Dlat, s),
        random_vector(rng, Dmid, s),
        random_large_output_lgc(rng, Dmid, Dobs, Q_type, s),
    )
end

function small_output_lgc_from_bottleneck(model::BottleneckLGC)
    return SmallOutputLGC(
        model.fan_out.A * model.H,
        model.fan_out.a + model.fan_out.A * model.h,
        model.fan_out.Q,
    )
end

# Generation of GaussMarkovModels.

function random_tv_gmm(rng::AbstractRNG, ordering, Dlat::Int, N::Int, s::StorageType)

    As = map(_ -> -random_nice_psd_matrix(rng, Dlat, s), 1:N)
    as = map(_ -> random_vector(rng, Dlat, s), 1:N)
    x0 = random_gaussian(rng, Dlat, s)

    # For some reason this operation seems to be _incredibly_ inaccurate. My guess is that
    # the numerics are horrible for some reason, but I'm not sure why. Hence we add a pretty
    # large constant to the diagonal to ensure that all Qs are positive definite.
    Qs = map(n -> x0.P - Symmetric(As[n] * x0.P * As[n]') + eltype(s)(1e-1) * I, 1:N)

    return GaussMarkovModel(ordering, As, as, Qs, x0)
end

function random_ti_gmm(rng::AbstractRNG, ordering, Dlat::Int, N::Int, s::StorageType)

    As = Fill(-random_nice_psd_matrix(rng, Dlat, s), N)
    as = Fill(random_vector(rng, Dlat, s), N)
    x0 = random_gaussian(rng, Dlat, s)

    # For some reason this operation seems to be _incredibly_ inaccurate. My guess is that
    # the numerics are horrible for some reason, but I'm not sure why. Hence we add a pretty
    # large constant to the diagonal to ensure that all Qs are positive definite.
    Qs = Fill(x0.P - As[1] * x0.P * As[1]' + eltype(s)(1e-1) * I, N)
    return GaussMarkovModel(ordering, As, as, Qs, x0)
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, gmm::T) where {T<:GaussMarkovModel}
    return Composite{T}(
        ordering = nothing,
        As = rand_tangent(rng, gmm.As),
        as = rand_tangent(rng, gmm.as),
        Qs = gmm_Qs_tangent(rng, gmm.Qs, storage_type(gmm)),
        x0 = rand_tangent(rng, gmm.x0),
    )
end

function gmm_Qs_tangent(
    rng::AbstractRNG, Qs::T, storage_type::StorageType,
) where {T<:Vector{<:AbstractMatrix}}
    return map(Q -> random_nice_psd_matrix(rng, size(Q, 1), storage_type), Qs)
end

function gmm_Qs_tangent(
    rng::AbstractRNG, Qs::T, storage_type::StorageType,
) where {T<:Fill{<:AbstractMatrix}}
    Δ = random_nice_psd_matrix(rng, size(first(Qs), 1), storage_type)
    return Composite{T}(value=Δ)
end

function gmm_Qs_tangent(
    rng::AbstractRNG, Qs::T, storage_type::StorageType,
) where {T<:Vector{<:Real}}
    return map(Q -> convert(eltype(storage_type), rand(rng) + 0.1), Qs)
end

function gmm_Qs_tangent(
    rng::AbstractRNG, Qs::T, storage_type::StorageType,
) where {T<:Fill{<:Real}}
    return Composite{T}(value=convert(eltype(storage_type), rand(rng) + 0.1))
end

# Generation of LGSSMs.

function random_lgssm(
    rng::AbstractRNG,
    ordering::Union{Forward, Reverse},
    ::Val{:time_varying},
    emission_type::Type{<:Union{SmallOutputLGC, LargeOutputLGC}},
    Dlat::Int,
    Dobs::Int,
    N::Int,
    Q_type::Val=Val{:dense},
    storage::StorageType=ArrayStorage(Float64),
)
    transitions = random_tv_gmm(rng, ordering, Dlat, N, storage)
    Hs = map(_ -> random_matrix(rng, Dobs, Dlat, storage), 1:N)
    hs = map(_ -> random_vector(rng, Dobs, storage), 1:N)
    Σs = map(_ -> random_nice_psd_matrix(rng, Dobs, Q_type, storage), 1:N)
    T = emission_type{eltype(Hs), eltype(hs), eltype(Σs)}
    emissions = StructArray{T}((Hs, hs, Σs))
    return LGSSM(transitions, emissions)
end

function random_lgssm(
    rng::AbstractRNG,
    ordering::Union{Forward, Reverse},
    ::Val{:time_invariant},
    emission_type::Type{<:Union{SmallOutputLGC, LargeOutputLGC}},
    Dlat::Int,
    Dobs::Int,
    N::Int,
    Q_type::Val=Val{:dense},
    storage::StorageType=ArrayStorage(Float64),
)
    transitions = random_ti_gmm(rng, ordering, Dlat, N, storage)
    Hs = Fill(random_matrix(rng, Dobs, Dlat, storage), N)
    hs = Fill(random_vector(rng, Dobs, storage), N)
    Σs = Fill(random_nice_psd_matrix(rng, Dobs, Q_type, storage), N)
    T = emission_type{eltype(Hs), eltype(hs), eltype(Σs)}
    emissions = StructArray{T}((Hs, hs, Σs))
    return LGSSM(transitions, emissions)
end

function random_lgssm(
    rng::AbstractRNG,
    ordering::Union{Forward, Reverse},
    ::Val{:time_varying},
    ::Type{ScalarOutputLGC},
    Dlat::Int,
    Dobs::Int,
    N::Int,
    storage::StorageType,
)
    transitions = random_tv_gmm(rng, ordering, Dlat, N, storage)
    Hs = map(_ -> random_vector(rng, Dlat, storage)', 1:N)
    hs = map(_ -> randn(rng, eltype(storage)), 1:N)
    Σs = map(_ -> convert(eltype(storage), rand(rng) + 0.1), 1:N)
    T = ScalarOutputLGC{eltype(Hs), eltype(hs), eltype(Σs)}
    emissions = StructArray{T}((Hs, hs, Σs))
    return LGSSM(transitions, emissions)
end

function random_lgssm(
    rng::AbstractRNG,
    ordering::Union{Forward, Reverse},
    ::Val{:time_invariant},
    ::Type{ScalarOutputLGC},
    Dlat::Int,
    Dobs::Int,
    N::Int,
    storage::StorageType,
)
    transitions = random_ti_gmm(rng, ordering, Dlat, N, storage)
    Hs = Fill(random_vector(rng, Dlat, storage)', N)
    hs = Fill(randn(rng, eltype(storage)), N)
    Σs = Fill(convert(eltype(storage), rand(rng) + 0.1), N)
    T = ScalarOutputLGC{eltype(Hs), eltype(hs), eltype(Σs)}
    emissions = StructArray{T}((Hs, hs, Σs))
    return LGSSM(transitions, emissions)
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, ssm::T) where {T<:LGSSM}
    Hs = ssm.emissions.A
    hs = ssm.emissions.a
    Σs = ssm.emissions.Q
    return Composite{T}(
        transitions = rand_tangent(rng, ssm.transitions),
        emissions = Composite{typeof(ssm.emissions)}(components=(
            A=rand_tangent(rng, Hs),
            a=rand_tangent(rng, hs),
            Q=gmm_Qs_tangent(rng, Σs, storage_type(ssm)),
        )),
    )
end

# function random_tv_scalar_lgssm(rng::AbstractRNG, Dlat::Int, N::Int, storage)
#     return ScalarLGSSM(random_tv_lgssm(rng, Dlat, 1, N, storage))
# end

# function random_ti_scalar_lgssm(rng::AbstractRNG, Dlat::Int, N::Int, storage)
#     return ScalarLGSSM(random_ti_lgssm(rng, Dlat, 1, N, storage))
# end

# function random_tv_posterior_lgssm(rng::AbstractRNG, Dlat::Int, Dobs::Int, N::Int, storage)
#     lgssm = random_tv_lgssm(rng, Dlat, Dobs, N, storage)
#     y = rand(rng, lgssm)
#     Σs = map(_ -> random_nice_psd_matrix(rng, Dobs, storage), eachindex(y))
#     return posterior(lgssm, y, Σs)
# end

# function random_ti_posterior_lgssm(rng::AbstractRNG, Dlat::Int, Dobs::Int, N::Int, storage)
#     lgssm = random_ti_lgssm(rng, Dlat, Dobs, N, storage)
#     y = rand(rng, lgssm)
#     Σs = Fill(random_nice_psd_matrix(rng, Dobs, storage), length(lgssm))
#     return posterior(lgssm, y, Σs)
# end

#
# Validation of internal consistency.
#

function validate_dims(model::Union{SmallOutputLGC, LargeOutputLGC})
    @test size(model.A) == (dim_out(model), dim_in(model))
    @test size(model.a) == (dim_out(model), )
    @test size(model.Q) == (dim_out(model), dim_out(model))
    return nothing
end

function validate_dims(model::ScalarOutputLGC)
    @test size(model.A) == (1, dim_in(model))
    @test size(model.a) == ()
    @test size(model.Q) == ()
    return nothing
end

function validate_dims(model::BottleneckLGC)
    validate_dims(model.fan_out)
    @test dim_in(model.fan_out) == size(model.H, 1)
    @test dim_in(model.fan_out) == length(model.h)
end

function validate_dims(gmm::GaussMarkovModel)

    # Check all vectors are the correct length.
    @test length(gmm.As) == length(gmm)
    @test length(gmm.as) == length(gmm)
    @test length(gmm.Qs) == length(gmm)

    # Check sizes of each element of the struct are correct.
    @test all(map(n -> dim_out(gmm[n]) == dim(gmm), eachindex(gmm)))
    @test all(map(n -> dim_in(gmm[n]) == dim(gmm), eachindex(gmm)))

    @test size(gmm.x0.m) == (dim(gmm),)
    @test size(gmm.x0.P) == (dim(gmm), dim(gmm))
    return nothing
end

function validate_dims(model::LGSSM)
    validate_dims(model.transitions)

    @test length(model) == length(model.transitions)
    @test length(model) == length(model.emissions)

    map(validate_dims, model.emissions)

    @test all(n -> dim(model.transitions) == dim_in(model.emissions[n]), eachindex(model))

    return nothing
end

# function __verify_model_properties(model, Dlat, Dobs, N, storage_type)
#     @test is_of_storage_type(model, storage_type)
#     @test length(model) == N
#     @test dim_obs(model) == Dobs
#     @test dim_latent(model) == Dlat
#     validate_dims(model)
#     return nothing
# end

# function __verify_model_properties(model, Dlat, N, storage_type)
#     return __verify_model_properties(model, Dlat, 1, N, storage_type)
# end


text:
# Advent of Code, Day 11
#
# Simulate energy growth of a field of "octopusses", with each
# exploding when they reach a level, infecting those around them.
#
# AK, 11 Dec 2021

# Parameters
f = "test_input.txt"
f = "input.txt"
iterations = 1000  # For part 2, need to go beyond 100

# Read data into list of rows, convert to matrix
d = split(strip(read(f, String)), "\n")
l = map(r -> [x - '0' for x in r], d)
m = Matrix(transpose(reshape(vcat(l...), (:, length(l)))))
println("Starting matrix:")
display(m)
nrows, ncols = size(m)

# Run the simulation
function simulate(m)

    # Start the simulation
    totFlashes = 0
    for t in 1:iterations

        # At the beginning of each step, add 1 to each cell
        m = m .+ 1

        # Initialize a matrix of flags indicating whether a cell
        # has already flashed this iteration
        notFlashed = ones(Int64, nrows, ncols)

        # Any cells with a value > 9 flash, incrementing the cells
        # around them, continue until no more flashes
        flashesFound = true
        while flashesFound
            flashesFound = false
            for r in 1:nrows
                for c in 1:ncols
                    if m[r,c] > 9 && notFlashed[r,c] == 1
                        flash(m, r, c)
                        notFlashed[r,c] = 0
                        flashesFound = true
                        totFlashes += 1
                    end
                end
            end
        end

        # Finally, set all the cells that have flashed to zero
        m = m .* notFlashed

        # Show totals and matrix this step
        println("\n\nAfter step $t, Total flashes = $totFlashes:")
        display(m)

        # For part 2, stop when all cells have flashed simulataneously
        if sum(notFlashed) == 0
            println("\nAll cells flashed at t = $t")
            break
        end
    end
end

# Flash one cell in the matrix: this increments all diagonal cells by one
function flash(m, r, c)
    for ri in r-1:r+1
        for ci in c-1:c+1
            if ri >= 1 && ri <= nrows && ci >= 1 && ci <= ncols
                m[ri,ci] += 1
            end
        end
    end
end

# Run the simulation
simulate(m)


text:
"""
    information_ratio(asset_returns, benchmark_returns; multiplier=1.0)

This function calculates the Information Ratio as the active return divided by the tracking error. The calculation equals William F. Sharpe's revision of the original version of the Sharpe Ratio (see function `sharpe_ratio`).

# Formula

    IR = E[asset_returns - benchmark_returns] / std(asset_returns - benchmark_returns) * sqrt(multiplier)

# Arguments
- `asset_returns`:      Vector of asset returns.
- `benchmark_returns`:  Vector or scalar value of benchmark returns having same same frequency (e.g. daily) as the provided returns.
- `multiplier`:         Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.

# Sources
- Sharpe, William F. (1994). "The Sharpe Ratio". The Journal of Portfolio Management.
"""
function information_ratio(asset_returns, benchmark_returns; multiplier=1.0)
    mean(asset_returns .- benchmark_returns) / std(asset_returns .- benchmark_returns) * sqrt(multiplier)
end


text:
function combine(obj::AbstractSingleFrequencyCombination, arg0::ObservationData, arg1::ObservationData)
    return jcall(obj, "combine", CombinedObservationData, (ObservationData, ObservationData), arg0, arg1)
end

function combine(obj::AbstractSingleFrequencyCombination, arg0::ObservationDataSet)
    return jcall(obj, "combine", CombinedObservationDataSet, (ObservationDataSet,), arg0)
end

function get_name(obj::AbstractSingleFrequencyCombination)
    return jcall(obj, "getName", JString, ())
end