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

@doc raw"""Represents a Persistent Disk resource in AWS.  An AWS EBS disk must exist before mounting to a container. The disk must also be in the same AWS zone as the kubelet. An AWS EBS disk can only be mounted as read/write once. AWS EBS volumes support ownership management and SELinux relabeling.

    IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource(;
        fsType=nothing,
        partition=nothing,
        readOnly=nothing,
        volumeID=nothing,
    )

    - fsType::String : Filesystem type of the volume that you want to mount. Tip: Ensure that the filesystem type is supported by the host operating system. Examples: \"ext4\", \"xfs\", \"ntfs\". Implicitly inferred to be \"ext4\" if unspecified. More info: https://kubernetes.io/docs/concepts/storage/volumes#awselasticblockstore
    - partition::Int32 : The partition in the volume that you want to mount. If omitted, the default is to mount by volume name. Examples: For volume /dev/sda1, you specify the partition as \"1\". Similarly, the volume partition for /dev/sda is \"0\" (or you can leave the property empty).
    - readOnly::Bool : Specify \"true\" to force and set the ReadOnly property in VolumeMounts to \"true\". If omitted, the default is \"false\". More info: https://kubernetes.io/docs/concepts/storage/volumes#awselasticblockstore
    - volumeID::String : Unique ID of the persistent disk resource in AWS (Amazon EBS volume). More info: https://kubernetes.io/docs/concepts/storage/volumes#awselasticblockstore
"""
mutable struct IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    partition::Any # spec type: Union{ Nothing, Int32 } # spec name: partition
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    volumeID::Any # spec type: Union{ Nothing, String } # spec name: volumeID

    function IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource(;fsType=nothing, partition=nothing, readOnly=nothing, volumeID=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("partition"), partition)
        setfield!(o, Symbol("partition"), partition)
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("volumeID"), volumeID)
        setfield!(o, Symbol("volumeID"), volumeID)
        o
    end
end # type IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource

const _property_map_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("partition")=>Symbol("partition"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("volumeID")=>Symbol("volumeID"))
const _property_types_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("partition")=>"Int32", Symbol("readOnly")=>"Bool", Symbol("volumeID")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource[property_name]

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

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


text:
int main() {
    boolean b=false;
    while(b){}
    return 0;
}

ext: jl

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
function img_to_data(img::AbstractArray{T,N}) where T<:Colorant where N
    AT = accum_type(T)
    aimg = AT.(img)
    pimg = parent(aimg)
    ch = channelview(pimg)
    data = reshape(ch, :, *(size(pimg)...))
end

kmeans(img::AbstractArray{T,N}, args...; kwargs...) where {T<:Colorant,N} =
    kmeans(img_to_data(img), args...; kwargs...)

fuzzy_cmeans(img::AbstractArray{T,N}, args...; kwargs...) where {T<:Colorant,N} =
    fuzzy_cmeans(img_to_data(img), args...; kwargs...)


text:
module DifferentialDynamicProgramming
using Requires
const DEBUG = false # Set this flag to true in order to print debug messages
# package code goes here

export QPTrace, boxQP, demoQP, Trace, iLQG, demo_linear

include("boxQP.jl")
include("iLQG.jl")
include("demo_linear.jl")
@require Control include("system_pendcart.jl")

function debug(x)
    DEBUG && print_with_color(:blue, string(x),"\n")
end

end # module


text:
using StochasticGroundMotionSimulation
using Plots
using QuadGK
using BenchmarkTools
using Test
# using ForwardDiff
using Polynomials
using PolynomialRoots

fi = exp10.(range(-2.0, stop=2.0, length=2001))

m = 6.0
r = 100.0
fas = FASParams(100.0, [1.0, 50.0, Inf], [1.0, 0.5], 200.0, 0.4, 0.039)
sdof = Oscillator(1.0)

Afi = map(f -> fourier_spectral_ordinate(f, m, r, fas), fi)
Hfi = transfer(fi, sdof)

fac = 1.5

plot(fi, Afi.^2, lab="Af^2", box=:true, grid=:on, xscale=:log10)
plot!(fi, (Hfi .* Afi).^2, lab="Hf^2Af^2")
vline!([sdof.f_n, sdof.f_n/fac, sdof.f_n*fac])
xlims!(sdof.f_n/ceil(fac),sdof.f_n*ceil(fac))

fun(f) = squared_transfer(f, sdof)

fac = 1.5

fj = collect(range(sdof.f_n/fac, stop=sdof.f_n*fac, length=101))
H2j = map(f -> fun(f), fj)

@time igk = quadgk(fun, sdof.f_n/fac, sdof.f_n*fac)[1]
@time isr = simpsons_rule(fj, H2j)
@time itr = trapezoidal_rule(fj, H2j)

@test igk ≈ isr rtol=1e-5
@test igk ≈ itr rtol=1e-5

fun(f) = squared_transfer(f, sdof) * ( fourier_spectral_ordinate(f, m, r, fas)^2 )

fj = collect(range(sdof.f_n/fac, stop=sdof.f_n*fac, length=101))
H2j = map(f -> fun(f), fj)

@time igk = quadgk(fun, sdof.f_n/fac, sdof.f_n*fac)[1]
@time isr = simpsons_rule(fj, H2j)
@time itr = trapezoidal_rule(fj, H2j)

@test igk ≈ isr rtol=1e-5
@test igk ≈ itr rtol=1e-4

plot(fi, Afi.^2, lab="Af^2", box=:true, grid=:on)
plot!(fi, (Hfi .* Afi).^2, lab="Hf^2Af^2")
# vline!([sdof.f_n, sdof.f_n/fac, sdof.f_n*fac])
xlims!(1e-3, sdof.f_n/fac)

fj = collect(range(1e-3, stop=sdof.f_n/fac, length=101))
H2j = map(f -> fun(f), fj)

@time igk = quadgk(fun, 1e-3, sdof.f_n/fac)[1]
@time isr = simpsons_rule(fj, H2j)
@time itr = trapezoidal_rule(fj, H2j)

@test igk ≈ isr rtol=1e-5
@test igk ≈ itr rtol=1e-3

fj = collect(range(sdof.f_n*fac, stop=1e1, length=101))
H2j = map(f -> fun(f), fj)

@time igk = quadgk(fun, sdof.f_n*fac, 1e1)[1]
@time isr = simpsons_rule(fj, H2j)
@time itr = trapezoidal_rule(fj, H2j)

@test igk ≈ isr rtol=1e-3
@test igk ≈ itr rtol=1e-3

fun(f) = (2*π*f)^2 * squared_transfer(f, sdof) * ( fourier_spectral_ordinate(f, m, r, fas)^2 )

fj = collect(range(sdof.f_n/fac, stop=sdof.f_n*fac, length=101))
H2j = map(f -> fun(f), fj)
fk = collect(range(sdof.f_n/fac, stop=sdof.f_n*fac, length=71))
H2k = map(f -> fun(f), fk)

@time igk = quadgk(fun, sdof.f_n/fac, sdof.f_n*fac)[1]
@time isr = simpsons_rule(fj, H2j)
@time itr = trapezoidal_rule(fk, H2k)

@test igk ≈ isr rtol=1e-5
@test igk ≈ itr rtol=1e-4

function fourier_derivative(f, m, r, fas)
    βQ = π * r / ( fas.Q0 * fas.cQ )
    βκ = π * fas.κ0
    fc, ~, ~ = corner_frequency(m, fas)

    return -βκ*f^4 - βQ*(1.0-fas.η)*f^(4-fas.η) + 2*(1.0-fc^2)*f^3 - βκ*fc^2*f^2 - βQ*(1.0-fas.η)*fc^2*f^(2-fas.η) + 2*fc^2*f
end

fas = FASParams(100.0, [1.0, 50.0, Inf], [1.0, 0.5], 200.0, 0.0, 0.039)

fi = exp10.(range(-2, stop=2, length=101))
Afi = fourier_spectrum(fi, m, r, fas)
dAfi = map(f -> fourier_derivative(f, m, r, fas), fi)

# FAS(f) = fourier_spectral_ordinate(f, m, r, fas)
# dFAS(f) = ForwardDiff.derivative(FAS, f)
#
# gAfi = map(f -> dFAS(f), fi)

βQ = π * r / ( fas.Q0 * fas.cQ )
βκ = π * fas.κ0
fc, ~, ~ = corner_frequency(m, fas)

dAfp = Polynomial([ 0, 2*fc^2, -fc^2*(βQ+βκ), 2*(1.0-fc^2), -(βQ+βκ)])
ri = PolynomialRoots.roots([ 0, 2*fc^2, -fc^2*(βQ+βκ), 2*(1.0-fc^2), -(βQ+βκ)])

real.(ri)
imag.(ri)

plot(fi, Afi, lab="FAS")
vline!([abs(ri[1])])

plot(fi, dAfi, lab="dFAS")
ylims!(-1e2,1e2)


text:
using DiffEqJump, DiffEqBase
using Test
using StableRNGs
rng = StableRNG(12345)

rate = (u, p, t) -> u[1]
affect! = function (integrator)
    integrator.u[1] -= 1
    integrator.u[2] += 1
end
jump = ConstantRateJump(rate, affect!)

prob = DiscreteProblem([0.0, 0.0], (0.0, 10.0))
jump_prob = JumpProblem(prob, Direct(), jump; rng=rng)

sol = solve(jump_prob, SSAStepper())

@test sol.t == [0.0, 10.0]
@test sol.u == [[0.0, 0.0], [0.0, 0.0]]

condition(u,t,integrator) = t == 5
function fuel_affect!(integrator)
  integrator.u[1] += 100
  reset_aggregated_jumps!(integrator)
end
cb = DiscreteCallback(condition, fuel_affect!, save_positions=(false, true))

sol = solve(jump_prob, SSAStepper(), callback=cb, tstops=[5])

@test sol.t[1:2] == [0.0, 5.0] # no jumps between t=0 and t=5
@test sol(5 + 1e-10) == [100, 0] # state just after fueling before any decays can happen

# test that callback initializer/finalizer is called and add_tstop! works as expected
random_tstops = rand(rng,100) .* 10 # 100 random Float64 between 0.0 and 10.0

function fuel_init!(cb,u,t,integrator)
  for tstop in random_tstops
    add_tstop!(integrator, tstop)
  end
  @test issorted(integrator.tstops)
end
finalizer_called = 0
fuel_finalize(cb, u, t, integrator) = global finalizer_called += 1

cb2 = DiscreteCallback(condition, fuel_affect!, initialize=fuel_init!, finalize=fuel_finalize)
sol = solve(jump_prob, SSAStepper(), callback=cb2)
for tstop in random_tstops
  @test tstop ∈ sol.t
end
@test finalizer_called == 1

# test for updating MassActionJump parameters
rs = [[1 => 1],[2=>1]]
ns = [[1 => -1, 2 => 1],[1=>1,2=>-1]]
p  = [1.0,0.0]
maj = MassActionJump(rs, ns; param_idxs=[1,2])
u₀ = [100,0]
tspan = (0.0,2000.0)
dprob = DiscreteProblem(u₀,tspan,p)
jprob = JumpProblem(dprob,Direct(),maj, save_positions=(false,false), rng=rng)
pcondit(u,t,integrator) = t==1000.0
function paffect!(integrator)
  integrator.p[1] = 0.0
  integrator.p[2] = 1.0
  reset_aggregated_jumps!(integrator)
end
sol = solve(jprob, SSAStepper(), tstops=[1000.0], callback=DiscreteCallback(pcondit,paffect!))
@test all(p .== [0.0,1.0])
@test sol[1,end] == 100

p   .= [1.0,0.0]
maj1 = MassActionJump([1 => 1],[1 => -1, 2 => 1]; param_idxs=1)
maj2 = MassActionJump([2 => 1],[1 => 1, 2 => -1]; param_idxs=2)
jprob = JumpProblem(dprob, Direct(), maj1, maj2, save_positions=(false,false), rng=rng)
sol = solve(jprob, SSAStepper(), tstops=[1000.0], callback=DiscreteCallback(pcondit,paffect!))
@test all(p .== [0.0,1.0])
@test sol[1,end] == 100

p2    = [1.0, 0.0, 0.0]
maj3  = MassActionJump([1 => 1],[1 => -1, 2 => 1]; param_idxs=3)
dprob = DiscreteProblem(u₀,tspan,p2)
jprob = JumpProblem(dprob, Direct(), maj1, maj2, maj3, save_positions=(false,false), rng=rng)
sol = solve(jprob, SSAStepper(), tstops=[1000.0], callback=DiscreteCallback(pcondit,paffect!))
@test all(p2 .== [0.0,1.0,0.0])
@test sol[1,end] == 100

p2    .= [1.0, 0.0, 0.0]
jprob = JumpProblem(dprob, Direct(), JumpSet(; massaction_jumps=[maj1, maj2, maj3]), save_positions=(false,false), rng=rng)
sol = solve(jprob, SSAStepper(), tstops=[1000.0], callback=DiscreteCallback(pcondit,paffect!))
@test all(p2 .== [0.0,1.0,0.0])
@test sol[1,end] == 100

p   .= [1.0,0.0]
dprob = DiscreteProblem(u₀,tspan,p)
maj4 = MassActionJump([[1 => 1],[2 => 1]],[[1 => -1, 2 => 1],[1 => 1, 2 => -1]]; param_idxs=[1,2])
jprob = JumpProblem(dprob, Direct(), maj4, save_positions=(false,false), rng=rng)
sol = solve(jprob, SSAStepper(), tstops=[1000.0], callback=DiscreteCallback(pcondit,paffect!))
@test all(p .== [0.0,1.0])
@test sol[1,end] == 100

# test scale_rates kwarg
p .= [1.0]
dprob = DiscreteProblem(u₀,tspan,p)
maj5 = MassActionJump([[1 => 2]],[[1 => -1, 2 => 1]]; param_idxs=[1])
jprob = JumpProblem(dprob, Direct(), maj5, save_positions=(false,false), rng=rng)
@test all(jprob.massaction_jump.scaled_rates .== [.5])
jprob = JumpProblem(dprob, Direct(), maj5, save_positions=(false,false), rng=rng, scale_rates=false)
@test all(jprob.massaction_jump.scaled_rates .== [1.0])


text:
using DataFrames
using CSV
using Plots
#%%
data_rk = open("rk_error.dat")
data_rk4 = open("rk4_error.dat")
data_eu = open("eu_error.dat")

df_rk=CSV.read(data_rk, header=false, delim=" ")
df_rk4=CSV.read(data_rk4, header=false, delim=" ")
df_eu=CSV.read(data_eu, header=false, delim=" ")
#%%
x_=4:2048
eu=x_.^-1
rk=x_.^-2
rk4=x_.^-4
#%%
plot(df_rk[1],df_rk[2], label="Runge-Kutta", color="black", dpi=600)
plot!(x_,rk, label="y=x^-2", color="black", line=(:dash))

plot!(df_rk4[1],df_rk4[2], label="RK 4th order", color="red")
plot!(x_,rk4, label="y=x^-4", color="red" , line=(:dash))

plot!(df_eu[1],df_eu[2], label="Euler", color="blue")
plot!(x_,eu, label="y=x^-1",color="blue", line=(:dash))

xaxis!(:log10, xlims=(4,2048))
yaxis!(:log10, ylims=(1e-15,Inf))
title!("error as a function of n_steps")
xlabel!("n_steps")
ylabel!("error")
plot!(legend=:bottom)
savefig("errors.png")

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
println("------------------------------------")
println("Sectors")
println("------------------------------------")
ti = time()
@timedtestset "Sector properties of $(TensorLabXD.type_repr(I))" for I in sectorlist
    Istr = TensorLabXD.type_repr(I)
    @testset "Sector $Istr: Basic properties" begin
        s = (randsector(I), randsector(I), randsector(I))
        @test eval(Meta.parse(sprint(show, I))) == I
        @test eval(Meta.parse(TensorLabXD.type_repr(I))) == I
        @test eval(Meta.parse(sprint(show,s[1]))) == s[1]
        @test @constinferred(hash(s[1])) == hash(deepcopy(s[1]))
        @test @constinferred(one(s[1])) == @constinferred(one(I))
        @constinferred dual(s[1])
        @constinferred dim(s[1])
        @constinferred frobeniusschur(s[1])
        @constinferred Nsymbol(s...)
        @constinferred Rsymbol(s...)
        @constinferred Bsymbol(s...)
        @constinferred Fsymbol(s..., s...)
        it = @constinferred s[1] ⊗ s[2]
        @constinferred ⊗(s..., s...)
    end
    @testset "Sector $Istr: Value iterator" begin
        @test eltype(values(I)) == I
        sprev = one(I)
        for (i, s) in enumerate(values(I))
            @test !isless(s, sprev) # confirm compatibility with sort order
            if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
                @test_throws ArgumentError values(I)[i]
                @test_throws ArgumentError TensorLabXD.findindex(values(I), s)
            elseif hasmethod(Base.getindex, Tuple{typeof(values(I)), Int})
                @test s == @constinferred (values(I)[i])
                @test TensorLabXD.findindex(values(I), s) == i
            end
            sprev = s
            i >= 10 && break
        end
        @test one(I) == first(values(I))
        if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
            @test_throws ArgumentError TensorLabXD.findindex(values(I), one(I))
        elseif hasmethod(Base.getindex, Tuple{typeof(values(I)), Int})
            @test (@constinferred TensorLabXD.findindex(values(I), one(I))) == 1
            for s in smallset(I)
                @test (@constinferred values(I)[TensorLabXD.findindex(values(I), s)]) == s
            end
        end
    end
    if hasfusiontensor(I)
        @testset "Sector $I: fusion tensor and F-move and R-move" begin
            for a in smallset(I), b in smallset(I)
                for c in ⊗(a,b)
                    X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4))
                    X2 = fusiontensor(b, a, c)
                    l =dim(a)*dim(b)*dim(c)
                    R = LinearAlgebra.transpose(Rsymbol(a, b, c))
                    sz = (l, convert(Int, Nsymbol(a, b, c)))
                    @test reshape(X1, sz) ≈ reshape(X2, sz) * R
                end
            end
            for a in smallset(I), b in smallset(I), c in smallset(I)
                for e in ⊗(a, b), f in ⊗(b, c)
                    for d in intersect(⊗(e, c), ⊗(a, f))
                        X1 = fusiontensor(a, b, e)
                        X2 = fusiontensor(e, c, d)
                        Y1 = fusiontensor(b, c, f)
                        Y2 = fusiontensor(a, f, d)
                        @tensor f1[-1,-2,-3,-4] := conj(Y2[a,f,d,-4])*conj(Y1[b,c,f,-3])*
                                                    X1[a,b,e,-1] * X2[e,c,d,-2]
                        if FusionStyle(I) isa MultiplicityFreeFusion
                            f2 = fill(Fsymbol(a,b,c,d,e,f)*dim(d), (1,1,1,1))
                        else
                            f2 = Fsymbol(a,b,c,d,e,f)*dim(d)
                        end
                        @test isapprox(f1, f2; atol = 1e-12, rtol = 1e-12)
                    end
                end
            end
        end
    end
    @testset "Sector $Istr: Unitarity of F-move" begin
        for a in smallset(I), b in smallset(I), c in smallset(I)
            for d in ⊗(a,b,c)
                es = collect(intersect(⊗(a,b), map(dual, ⊗(c,dual(d)))))
                fs = collect(intersect(⊗(b,c), map(dual, ⊗(dual(d),a))))
                if FusionStyle(I) isa MultiplicityFreeFusion
                    @test length(es) == length(fs)
                    F = [Fsymbol(a,b,c,d,e,f) for e in es, f in fs]
                else
                    Fblocks = Vector{Any}()
                    for e in es
                        for f in fs
                            Fs = Fsymbol(a,b,c,d,e,f)
                            push!(Fblocks, reshape(Fs, (size(Fs, 1)*size(Fs, 2), size(Fs, 3)*size(Fs, 4))))
                        end
                    end
                    F = hvcat(length(fs), Fblocks...)
                end
                @test isapprox(F'*F, one(F); atol = 1e-12, rtol = 1e-12)
            end
        end
    end
    @testset "Sector $Istr: Pentagon equation" begin
        for a in smallset(I), b in smallset(I), c in smallset(I), d in smallset(I)
            for f in ⊗(a,b), h in ⊗(c,d)
                for g in ⊗(f,c), i in ⊗(b,h)
                    for e in intersect(⊗(g,d), ⊗(a,i))
                        if FusionStyle(I) isa MultiplicityFreeFusion
                            p1 = Fsymbol(f,c,d,e,g,h) * Fsymbol(a,b,h,e,f,i)
                            p2 = zero(p1)
                            for j in ⊗(b,c)
                                p2 += Fsymbol(a,b,c,g,f,j) *
                                        Fsymbol(a,j,d,e,g,i) *
                                        Fsymbol(b,c,d,i,j,h)
                            end
                            @test isapprox(p1, p2; atol = 1e-12, rtol = 1e-12)
                        else
                            @tensor p1[λ,μ,ν,κ,ρ,σ] := Fsymbol(f,c,d,e,g,h)[λ,μ,ν,τ]*
                                                        Fsymbol(a,b,h,e,f,i)[κ,τ,ρ,σ]
                            p2 = zero(p1)
                            for j in ⊗(b,c)
                                @tensor p2[λ,μ,ν,κ,ρ,σ]  += Fsymbol(a,b,c,g,f,j)[κ,λ,α,β]*
                                                            Fsymbol(a,j,d,e,g,i)[β,μ,τ,σ]*
                                                            Fsymbol(b,c,d,i,j,h)[α,τ,ν,ρ]
                            end
                            @test isapprox(p1, p2; atol = 1e-12, rtol = 1e-12)
                        end
                    end
                end
            end
        end
    end
    @testset "Sector $Istr: Hexagon equation" begin
        for a in smallset(I), b in smallset(I), c in smallset(I)
            for e in ⊗(c,a), g in ⊗(c,b)
                for d in intersect(⊗(e,b), ⊗(a,g))
                    if FusionStyle(I) isa MultiplicityFreeFusion
                        p1 = Rsymbol(c,a,e)*Fsymbol(a,c,b,d,e,g)*Rsymbol(b,c,g)
                        p2 = zero(p1)
                        for f in ⊗(a,b)
                            p2 += Fsymbol(c,a,b,d,e,f)*Rsymbol(c,f,d)*Fsymbol(a,b,c,d,f,g)
                        end
                        @test isapprox(p1, p2; atol = 1e-12, rtol = 1e-12)
                    else
                        @tensor p1[α,β,μ,ν] := Rsymbol(c,a,e)[α,λ]*
                                                Fsymbol(a,c,b,d,e,g)[λ,β,γ,ν]*
                                                Rsymbol(b,c,g)[γ,μ]
                        p2 = zero(p1)
                        for f in ⊗(a,b)
                            @tensor p2[α,β,μ,ν] += Fsymbol(c,a,b,d,e,f)[α,β,δ,σ]*
                                                    Rsymbol(c,f,d)[σ,ψ]*
                                                    Fsymbol(a,b,c,d,f,g)[δ,ψ,μ,ν]
                        end
                        @test isapprox(p1, p2; atol = 1e-12, rtol = 1e-12)
                    end
                end
            end
        end
    end
end
tf = time()
printstyled("Finished sector tests in ",
            string(round(tf-ti; sigdigits=3)),
            " seconds."; bold = true, color = Base.info_color())
println()


text:
function central_diff(u)
    du1 = diff(u; dims=1) / 2
    a = [du1[[1], :]; du1]
    a .+= [du1; du1[[end], :]]

    du2 = diff(u; dims=2) / 2
    b = [du2[[1], :]; du2]
    b .+= [du2; du2[[end], :]]

    return a, b
end

"""
    normalize(img)
Normalize an image based on the mean and the standard deviation

# TODO: check this
"""
normalize(img) = (img .- minimum(img)) ./ (maximum(img) .- minimum(img))


text:
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

using AbstractGPsMakie

DocMeta.setdocmeta!(
    AbstractGPsMakie, :DocTestSetup, :(using AbstractGPsMakie); recursive=true
)

makedocs(;
    modules=[AbstractGPsMakie],
    authors="David Widmann",
    repo="https://github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl/blob/{commit}{path}#{line}",
    sitename="AbstractGPsMakie.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliagaussianprocesses.github.io/AbstractGPsMakie.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md", "api.md"],
    #strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl",
    push_preview=true,
    devbranch="main",
)


text:
# Alex Hoffman
#cd("/Users/hoffm/.julia/dev")
using astrofunk
using Test

#variables for testing
Earth_a, Mars_a, Sun_mu = 149600000, 227920000, 1.3271*10^11
pos = [-0.270; -0.420; 0] #nondimensional
vel = [0.300; -1.00; 0] #nondimensional
μ₃ = 0.012150664267189

#test loop
@testset "astrofunk.jl" begin
    # Write your tests here.
    @test DCM(pi,'x',false) ≈ [1 0 0; 0 -1 0; 0 0 -1]
    @test circHohmann(Earth_a, Mars_a, Sun_mu) ≈ [2.944004579970198, 2.648345397475325, 44729392.19214477]
    @test JC([pos pos], [vel vel], μ₃) ≈ permutedims([3.186469735193683; 3.186469735193683])
    @test CR3BP_EOM([pos;vel],μ₃,'x') ≈ -0.135541250504792
    @test U_star(pos,μ₃,"yy") ≈ 10.724399149271521
end


text:
module Optimizers
  #collection of optimizers.

  import ..MLAlgorithm

  include("gradientdescent.jl")
  include("particleswarm.jl")
  include("backpropagation.jl")

end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
type MBIE<:FiniteHorizonAgent
    policy::Matrix{UInt64}
    ns::Array{UInt64, 4}
    phat::Array{Float64, 4}
    n::Array{UInt64, 3}
    rewards::Array{Float64, 3}
    Vopt::Matrix{Float64}
    δ::Float64
    ϵ::Float64
    rewards_known::Bool
    maxobs::Float64
    maxRet::Float64
    maxR::Float64
    keep_updating::Bool
end

function MBIE(S, A, H, rewards, δ, ϵ, keep_updating=false)
    policy = zeros(UInt64, S, H)
    ns = zeros(UInt64, S, S, A, H)
    phat = zeros(S, S, A, H)
    phat[1, :, :, :] = 1.
    n = zeros(UInt64, S, A, H)
    maxRet = sum(maximum(reshape(rewards, S*A, H), 1))
    # use asymptotic form here, may be off by a constant factor
    maxobs = 4 * S / ϵ^2 * H^4 + 1 / ϵ^2 * H^4 * log(S * A * H / ϵ / δ)

    m = MBIE(policy, ns, phat, n, rewards, zeros(S, H+1), δ, ϵ, true, maxobs, maxRet, maximum(rewards[:]), keep_updating)
    rand!(m.policy, 1:A) # Initialize with uniformly random policy
    m
end
function observe!(m::MBIE, s, a, t, r, sn)
    if m.keep_updating || m.n[s,a,t] < m.maxobs
        m.ns[sn, s, a, t] += 1
        m.n[s, a, t] +=  1
        if !m.rewards_known
            m.rewards[s, a, t] += r
        end
    end
end
maxV(m::MBIE) = m.maxRet
maxR(m::MBIE) = m.maxR
function update_policy!(m::MBIE)
    S = nS(m)
    H = horizon(m)
    A = nA(m)
    δ = m.δ
    ϵ = m.ϵ
    Vmax = maxV(m)
    Rmax = maxR(m)
    Q = zeros(A)

    # precomputations for confidence intervals
    m_val = m.maxobs
    lR = log(4 * S * A * m_val / δ) / 2
    lC = if S < 30
        log(2^S - 2) + log(2 * S * A * m_val / δ)
    else
        S*log(2) + log(2 * S * A * m_val / δ)
    end
    m.phat .= m.ns ./ reshape(m.n, 1, S, A, H)
    for t=H:-1:1
        V = m.Vopt[:, t+1]
        ind = sortperm(V)
        curmaxV = min(V[ind[end]], Vmax)

        for s ∈ 1:S
            for a ∈ 1:A
                slack = √(2 * lC / m.n[s, a, t])
                mass_left = 1.
                EV = 0.
                for i ∈ S:-1:1
                    ii = ind[i]
                    ph = m.phat[ii, s, a, t]
                    pmax = min(mass_left, ph + slack)
                    slack -= pmax - ph
                    mass_left -= pmax
                    EV += pmax * V[ii]
                    if mass_left < 3 * eps(typeof(mass_left))
                        break
                    end
                end

                r = m.rewards[s, a, t]
                if !m.rewards_known
                    min(m.rewards[s, a, t] / n[s, a, t] + Rmax * √(lR / m.n[s, a, t]), Rmax)
                end

                Q[a] = r + EV
            end
            bestA = indmax(Q)
            m.policy[s, t] = bestA
            m.Vopt[s, t] = Q[bestA]
        end
    end
end


text:
function longest_common_prefix(prefixa, prefixb)
    maxplength = min(sizeof(prefixa), sizeof(prefixb))
    maxplength == 0 && return ""
    idx = findfirst(i -> (prefixa[i] != prefixb[i]), 1:maxplength)
    idx = idx === nothing ? maxplength : idx - 1
    prefixa[1:idx]
end

function skip_to_nl(str, idxend)
    while (idxend < sizeof(str)) && str[idxend] != '\n'
        idxend = nextind(str, idxend)
    end
    idxend > sizeof(str) ? prevind(str, idxend) : idxend
end

tostr(buf::IOBuffer) = _unescape_string(String(take!(buf)))

"""
parse_string_or_cmd(ps)

When trying to make an `INSTANCE` from a string token we must check for
interpolating opoerators.
"""
function parse_string_or_cmd(ps::ParseState, prefixed=false)
    sfullspan = ps.nt.startbyte - ps.t.startbyte
    sspan = 1 + ps.t.endbyte - ps.t.startbyte

    istrip = (kindof(ps.t) === Tokens.TRIPLE_STRING) || (kindof(ps.t) === Tokens.TRIPLE_CMD)
    iscmd = kindof(ps.t) === Tokens.CMD || kindof(ps.t) === Tokens.TRIPLE_CMD

    lcp = nothing
    exprs_to_adjust = []
    function adjust_lcp(expr::EXPR, last=false)
        if isliteral(expr)
            push!(exprs_to_adjust, expr)
            str = valof(expr)
            (isempty(str) || (lcp !== nothing && isempty(lcp))) && return
            (last && str[end] == '\n') && return (lcp = "")
            idxstart, idxend = 2, 1
            prevpos = idxend
            while nextind(str, idxend) - 1 < sizeof(str) && (lcp === nothing || !isempty(lcp))
                idxend = skip_to_nl(str, idxend)
                idxstart = nextind(str, idxend)
                prevpos1 = idxend
                while nextind(str, idxend) - 1 < sizeof(str)
                    c = str[nextind(str, idxend)]
                    if c == ' ' || c == '\t'
                        idxend += 1
                    elseif c == '\n'
                        # All whitespace lines in the middle are ignored
                        idxend += 1
                        idxstart = idxend + 1
                    else
                        prefix = str[idxstart:idxend]
                        lcp = lcp === nothing ? prefix : longest_common_prefix(lcp, prefix)
                        break
                    end
                    if idxend <= prevpos1
                        throw(CSTInfiniteLoop("Infinite loop in adjust_lcp"))
                    else
                        prevpos1 = idxend
                    end
                end
                if idxend < prevpos
                    throw(CSTInfiniteLoop("Infinite loop in adjust_lcp"))
                else
                    prevpos = idxend
                end
            end
            if idxstart != nextind(str, idxend)
                prefix = str[idxstart:idxend]
                lcp = lcp === nothing ? prefix : longest_common_prefix(lcp, prefix)
            end
        end
    end

    isinterpolated = false
    erroredonlast = false

    t_str = val(ps.t, ps)
    if istrip && length(t_str) == 6
        if iscmd
            return wrapwithcmdmacro(EXPR(:TRIPLESTRING , sfullspan, sspan, ""))
        else
            return EXPR(:TRIPLESTRING, sfullspan, sspan, "")
        end
    elseif length(t_str) == 2
        if iscmd
            return wrapwithcmdmacro(EXPR(:STRING , sfullspan, sspan, ""))
        else
            return EXPR(:STRING, sfullspan, sspan, "")
        end
    elseif prefixed != false
        _val = istrip ? t_str[4:prevind(t_str, sizeof(t_str), 3)] : t_str[2:prevind(t_str, sizeof(t_str))]
        if iscmd
            _val = replace(_val, "\\\\" => "\\")
            _val = replace(_val, "\\`" => "`")
        else
            _val = unescape_prefixed(_val)
        end
        expr = EXPR(istrip ? :TRIPLESTRING : :STRING, sfullspan, sspan, _val)
        if istrip
            adjust_lcp(expr)
            ret = EXPR(:string, EXPR[expr], nothing, sfullspan, sspan)
        else
            return iscmd ? wrapwithcmdmacro(expr) : expr
        end
    else
        ret = EXPR(:string, EXPR[], EXPR[], sfullspan, sspan)
        input = IOBuffer(t_str)
        startbytes = istrip ? 3 : 1
        seek(input, startbytes)
        b = IOBuffer()
        prevpos = position(input)
        while !eof(input)
            c = read(input, Char)
            if c == '\\'
                write(b, c)
                write(b, read(input, Char))
            elseif c == '$'
                isinterpolated = true
                lspan = position(b)
                str = tostr(b)
                ex = EXPR(:STRING, lspan + startbytes, lspan + startbytes, str)
                if position(input) == (istrip ? 3 : 1) + 1
                    # Need to add empty :STRING at start to account for \"
                    pushtotrivia!(ret, ex)
                elseif !isempty(str)
                    push!(ret, ex)
                end
                istrip && adjust_lcp(ex)
                startbytes = 0
                op = EXPR(:OPERATOR, 1, 1, "\$")
                if peekchar(input) == '('
                    skip(input, 1) # skip past '('
                    lpfullspan = -position(input)
                    if iswhitespace(peekchar(input)) || peekchar(input) === '#'
                        read_ws_comment(input, readchar(input))
                    end
                    lparen = EXPR(:LPAREN, lpfullspan + position(input) + 1, 1)
                    rparen = EXPR(:RPAREN, 1, 1)

                    prev_input_size = input.size
                    input.size = input.size - (istrip ? 3 : 1)
                    # We're reusing a portion of the string from `ps` so we need to make sure `ps1` knows where the end of the string is.
                    ps1 = ParseState(input)

                    if kindof(ps1.nt) === Tokens.RPAREN
                        push!(ret, EXPR(:ERRORTOKEN, EXPR[], nothing))
                        pushtotrivia!(ret, op)
                        pushtotrivia!(ret, lparen)
                        pushtotrivia!(ret, rparen)
                        seek(input, ps1.nt.startbyte + 1)
                    else
                        interp_val = @closer ps1 :paren parse_expression(ps1, true)
                        push!(ret, interp_val)
                        pushtotrivia!(ret, op)
                        pushtotrivia!(ret, lparen)
                        if kindof(ps1.nt) === Tokens.RPAREN
                            # Need to check the parenthese were actually closed.
                            pushtotrivia!(ret, rparen)
                            seek(input, ps1.nt.startbyte + 1)
                        else
                            pushtotrivia!(ret, EXPR(:RPAREN, 0, 0))
                            seek(input, ps1.nt.startbyte) # We don't skip ahead one as there wasn't a closing paren
                        end
                    end
                    # Compared to flisp/JuliaParser, we have an extra lookahead token,
                    # so we need to back up one here
                    input.size = prev_input_size
                elseif Tokenize.Lexers.iswhitespace(peekchar(input)) || peekchar(input) === '#'
                    pushtotrivia!(ret, op)
                    push!(ret, mErrorToken(ps, StringInterpolationWithTrailingWhitespace))
                elseif sspan == position(input) + (istrip ? 3 : 1)
                    # Error. We've hit the end of the string
                    pushtotrivia!(ret, op)
                    push!(ret, mErrorToken(ps, StringInterpolationWithTrailingWhitespace))
                else
                    pos = position(input)
                    ps1 = ParseState(input)
                    next(ps1)
                    if kindof(ps1.t) === Tokens.WHITESPACE
                        error("Unexpected whitespace after \$ in String")
                    else
                        t = INSTANCE(ps1)
                    end
                    # Attribute trailing whitespace to the string
                    t = adjustspan(t)
                    push!(ret, t)
                    pushtotrivia!(ret, op)
                    seek(input, pos + t.fullspan)
                end
            else
                write(b, c)
            end
            prevpos = loop_check(input, prevpos)
        end

        # handle last String section
        lspan = position(b)
        if erroredonlast
            ex = EXPR(istrip ? :TRIPLESTRING : :STRING, (istrip ? 3 : 1) + (ps.nt.startbyte - ps.t.endbyte - 1), istrip ? 3 : 1, "")
            pushtotrivia!(ret, ex)
        elseif b.size == 0
            ex = mErrorToken(ps, Unknown)
            push!(ret, ex)
        else
            str = tostr(b)
            if istrip
                str = str[1:prevind(str, lastindex(str), 3)]
                # only mark non-interpolated triple strings
                ex = EXPR(length(ret) == 0 ? :TRIPLESTRING : :STRING, lspan + ps.nt.startbyte - ps.t.endbyte - 1 + startbytes, lspan + startbytes, str)
                adjust_lcp(ex, true)
            else
                str = str[1:prevind(str, lastindex(str))]
                ex = EXPR(:STRING, lspan + ps.nt.startbyte - ps.t.endbyte - 1 + startbytes, lspan + startbytes, str)
            end
            if isempty(str)
                pushtotrivia!(ret, ex)
            else
                push!(ret, ex)
            end
        end
    end

    single_string_T = (:STRING, :TRIPLESTRING, literalmap(kindof(ps.t)))
    if istrip
        if lcp !== nothing && !isempty(lcp)
            for expr in exprs_to_adjust
                for (i, a) in enumerate(ret.args)
                    if expr == a
                        ret.args[i].val = replace(valof(expr), "\n$lcp" => "\n")
                        break
                    end
                end
            end
        end
        # Drop leading newline
        if !isempty(ret.args) && isliteral(ret.args[1]) && headof(ret.args[1]) in single_string_T &&
                !isempty(valof(ret.args[1])) && valof(ret.args[1])[1] == '\n'
            ret.args[1] = dropleadlingnewline(ret.args[1])
        end
    end

    if (length(ret.args) == 1 && isliteral(ret.args[1]) && headof(ret.args[1]) in single_string_T) && !isinterpolated
        ret = ret.args[1]
    end
    update_span!(ret)

    return iscmd ? wrapwithcmdmacro(ret) : ret
end

function adjustspan(x::EXPR)
    x.fullspan = x.span
    return x
end

dropleadlingnewline(x::EXPR) = EXPR(headof(x), x.fullspan, x.span, valof(x)[2:end])

wrapwithcmdmacro(x) =EXPR(:macrocall, EXPR[EXPR(:globalrefcmd, 0, 0), EXPR(:NOTHING, 0, 0), x])

"""
    parse_prefixed_string_cmd(ps::ParseState, ret::EXPR)

Parse prefixed strings and commands such as `pre"text"`.
"""
function parse_prefixed_string_cmd(ps::ParseState, ret::EXPR)
    arg = parse_string_or_cmd(next(ps), ret)

    if ret.head === :IDENTIFIER && valof(ret) == "var" && isstringliteral(arg) && VERSION > v"1.3.0-"
        return EXPR(:NONSTDIDENTIFIER, EXPR[ret, arg], nothing)
    elseif headof(arg) === :macrocall && headof(arg.args[1]) === :globalrefcmd
        mname = EXPR(:IDENTIFIER, ret.fullspan, ret.span, string("@", valof(ret), "_cmd")) # NOTE: sizeof(valof(mname)) != mname.span
        return EXPR(:macrocall, EXPR[mname, EXPR(:NOTHING, 0, 0), arg.args[3]], nothing)
    elseif is_getfield(ret)
        if headof(ret.args[2]) === :quote || headof(ret.args[2]) === :quotenode
            str_type = valof(ret.args[2].args[1]) isa String ? valof(ret.args[2].args[1]) : "" # to handle some malformed case
            ret.args[2].args[1] = setparent!(EXPR(:IDENTIFIER, ret.args[2].args[1].fullspan, ret.args[2].args[1].span, string("@", str_type, "_str")), ret.args[2])
        else
            str_type = valof(ret.args[2]) isa String ? valof(ret.args[2]) : "" # to handle some malformed case
            ret.args[2] = EXPR(:IDENTIFIER, ret.args[2].fullspan, ret.args[2].span, string("@", str_type, "_str"))
        end

        return EXPR(:macrocall, EXPR[ret, EXPR(:NOTHING, 0, 0), arg], nothing)
    else
        return EXPR(:macrocall, EXPR[EXPR(:IDENTIFIER, ret.fullspan, ret.span, string("@", valof(ret), "_str")), EXPR(:NOTHING, 0, 0), arg], nothing)
    end
end

function unescape_prefixed(str)
    edits = UnitRange{Int}[]
    start = -1
    for (i, c) in enumerate(str)
        if start == -1 && c === '\\'
            start = i
        end
        if start > -1
            if c === '\\'
            elseif c === '\"'
                push!(edits, start:i)
                start = -1
            else
                start = -1
            end
        end
    end

    if !isempty(edits) || start > -1
        str1 = deepcopy(str)
        if start > -1
            # slashes preceding closing '"'
            n = div(length(start:length(str)), 2) - 1
            str1 = string(str1[1:prevind(str1, start)], repeat("\\", n + 1))
        end

        for e in reverse(edits)
            n = div(length(e), 2) - 1
            str1 = string(str1[1:prevind(str1, first(e))], string(repeat("\\", n), "\""), str1[nextind(str1, last(e)):lastindex(str1)])

        end
        return str1
    end
    return str
end


text:
"""
    PythonDisplay()

Like TextDisplay() but prints to Python's stdout.
"""
struct PythonDisplay <: AbstractDisplay end

function Base.display(d::PythonDisplay, m::MIME, @nospecialize(x))
    istextmime(m) || throw(MethodError(display, (d, m, x)))
    buf = IOBuffer()
    io = IOContext(buf, :limit=>true)
    try
        show(io, m, x)
    catch
        throw(MethodError(display, (d, m, x)))
    end
    data = String(take!(buf))
    pyprint(data)
    return
end

function Base.display(d::PythonDisplay, @nospecialize(x))
    display(d, MIME("text/plain"), x)
end

"""
    IPythonDisplay()

For displaying multimedia with IPython's display mechanism.
"""
struct IPythonDisplay <: AbstractDisplay end

function Base.display(d::IPythonDisplay, m::MIME, @nospecialize(x))
    ipy = pyimport("IPython")
    buf = IOBuffer()
    io = IOContext(buf, :limit=>true)
    dict = pydict()
    try
        show(io, m, x)
    catch
        throw(MethodError(display, (d, m, x)))
    end
    data = take!(buf)
    dict[string(m)] = istextmime(m) ? pystr_fromUTF8(data) : pybytes(data)
    ipy.display.display(dict, raw=true)
    return
end

function Base.display(d::IPythonDisplay, @nospecialize(x))
    ipy = pyimport("IPython")
    if ispy(x)
        ipy.display.display(x)
        return
    end
    buf = IOBuffer()
    io = IOContext(buf, :limit=>true)
    dict = pydict()
    for m in Utils.mimes_for(x)
        try
            show(io, MIME(m), x)
        catch
            continue
        end
        data = take!(buf)
        dict[m] = istextmime(m) ? pystr_fromUTF8(data) : pybytes(data)
    end
    length(dict) == 0 && throw(MethodError(display, (d, x)))
    ipy.display.display(dict, raw=true)
    return
end

function init_ipython()
    # EXPERIMENTAL: IPython integration
    if C.CTX.is_embedded && CONFIG.auto_ipython_display
        is_ipython = ("IPython" in pysysmodule.modules) && !pyisnone(pysysmodule.modules["IPython"].get_ipython())
        if is_ipython
            pushdisplay(PythonDisplay())
            pushdisplay(IPythonDisplay())
        end
    end
end


text:
# Definition of StrBootTest "class" for holding intermediate results, with associated utilities and get functions

"Auxilliary weight types: `rademacher`, `mammen`, `webb`, `normal`, `gamma`"
@enum AuxWtType rademacher mammen webb normal gamma

"p value types: `symmetric`, `equaltail`, `lower`, `upper`"
@enum PType symmetric equaltail lower upper

"Multiple hypothesis adjustment types: `nomadj`, `bonferroni`, `sidak`"
@enum MAdjType nomadj bonferroni sidak

"Bootstrap distribution statistics optionally returned"
@enum DistStatType nodist t numer

struct StrClust{T<:Real}
	N::Int; multiplier::T; even::Bool
	order::Vector{Int64}
	info::Vector{UnitRange{Int64}}
end

struct StrFE{T<:Real}
	is::SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}
  wt::Vector{T}
end

mutable struct StrEstimator{T<:AbstractFloat}
  isDGP::Bool; LIML::Bool; Fuller::T; κ::T
  R₁perp::Matrix{T}; Rpar::Matrix{T}

  kZ::Int64
  y₁::Vector{T}; ü₁::Vector{T}; u⃛₁::Vector{T}; β̈::Vector{T}; γ̈::Vector{T}; β̂₀::Vector{T}; invXXXy₁par::Vector{T}
  Yendog::Vector{Bool}
  invZperpZperp::Symmetric{T,Matrix{T}}; XZ::Matrix{T}; PXZ::Matrix{T}; YPXY::Symmetric{T,Matrix{T}}; R₁invR₁R₁::Matrix{T}
	restricted::Bool; RperpX::Matrix{T}; RperpXperp::Matrix{T}; RRpar::Matrix{T}; RparY::Matrix{T}; RR₁invR₁R₁::Matrix{T}
	∂β̂∂r::Matrix{T}; YY::Symmetric{T,Matrix{T}}; AR::Matrix{T}; XAR::Matrix{T}; R₁invR₁R₁Y::Matrix{T}; invXXXZ::Matrix{T}; Ü₂::Matrix{T}; XinvXX::Matrix{T}; Rt₁::Vector{T}
	invXX::Symmetric{T,Matrix{T}}; Y₂::Matrix{T}; X₂::Matrix{T}; invH::Symmetric{T,Matrix{T}}
	y₁par::Vector{T}; Xy₁par::Vector{T}
	A::Symmetric{T,Matrix{T}}; Z::Matrix{T}; Zperp::Matrix{T}; X₁::Matrix{T}
	FillingT₀::Matrix{Matrix{T}}
	WXAR::Matrix{T}; S⋂PXYZperp::Vector{Matrix{T}}; S⋂YX::Vector{Matrix{T}}; CT_XAR::Vector{Matrix{T}}; CT_FE⋂PY::Vector{Matrix{T}}

  # IV/GMM only
  ZZ::Symmetric{T,Matrix{T}}; XY₂::Matrix{T}; XX::Symmetric{T,Matrix{T}}; H_2SLS::Symmetric{T,Matrix{T}}; V::Matrix{T}; ZY₂::Matrix{T}; X₂Y₂::Matrix{T}; X₁Y₂::Matrix{T}; ZR₁ZR₁::Symmetric{T,Matrix{T}}; X₂ZR₁::Matrix{T}; ZR₁Y₂::Matrix{T}; X₁ZR₁::Matrix{T}
  ZZR₁::Matrix{T}; X₂y₁::Vector{T}; X₁y₁::Vector{T}; Zy₁::Vector{T}; ZXinvXXXZ::Matrix{T}; H_2SLSmZZ::Symmetric{T,Matrix{T}}
  ZXinvXXXy₁par::Vector{T}; t₁Y::Vector{T}
  Y₂y₁::Vector{T}; twoR₁Zy₁::Vector{T}
  y₁y₁::T; y₁pary₁par::T
  X₂y₁par::Vector{T}; X₁y₁par::Vector{T}; Zy₁par::Vector{T}
  Y₂y₁par::Vector{T}
  Rperp::Matrix{T}; ZR₁::Matrix{T}
  kX::Int64
	S⋂y₁X₁::Matrix{T}; S⋂y₁X₂::Matrix{T}; invZperpZperpZperpy₁::Vector{T}; invZperpZperpZperpX::Matrix{T}; invZperpZperpZperpZ::Matrix{T}; invZperpZperpZperpZR₁::Matrix{T}; invZperpZperpZperpY₂::Matrix{T}; S⋂XZperp::Array{T,3}; S⋂Xy₁::Array{T,3}; S⋂ZperpZperp::Array{T,3}; ZperpX::Matrix{T}; S⋂Zperpy₁::Array{T,3}; S⋂X_ZR₁::Array{T,3}
	Zperp_ZR₁::Matrix{T}; S⋂Zperp_ZR₁::Array{T,3}; FillingT₀₀::Array{T,3}; ∂FillingT₀∂r::Array{T,3}; S⋂y₁X₀::Matrix{T}; ∂S⋂y₁X∂r::Array{T,3}
Π̂::Matrix{T}; _X₁::Matrix{T}; _Z::Matrix{T}; ZperpZ::Matrix{T}; ZperpY₂::Matrix{T}; _ZR₁::Matrix{T}

  StrEstimator{T}(isDGP, LIML, Fuller, κ) where T<:AbstractFloat = new(isDGP, LIML, Fuller, κ, Matrix{T}(undef,0,0))
end

mutable struct StrBootTest{T<:AbstractFloat}
  R::Matrix{T}; r::Vector{T}; R₁::Matrix{T}; r₁::Vector{T}
  y₁::Vector{T}; X₁::Matrix{T}; Y₂::Matrix{T}; X₂::Matrix{T}
  wt::Vector{T}; fweights::Bool
  LIML::Bool; Fuller::T; κ::T; ARubin::Bool
  B::Int64; auxtwtype::AuxWtType; rng::AbstractRNG; maxmatsize::Float16
  ptype::PType; null::Bool; bootstrapt::Bool
	ID::Matrix{Int64}; nbootclustvar::Int8; nerrclustvar::Int8; issorted::Bool; small::Bool
  FEID::Vector{Int64}; FEdfadj::Bool
  level::T; rtol::T
  madjtype::MAdjType; NH₀::Int16
  ML::Bool; β̈::Vector{T}; A::Symmetric{T,Matrix{T}}; sc::Matrix{T}
  willplot::Bool; gridmin::Vector{T}; gridmax::Vector{T}; gridpoints::Vector{Float32}

  q::Int16; twotailed::Bool; scorebs::Bool; robust::Bool

  WRE::Bool; initialized::Bool; NFE::Int64; FEboot::Bool; granular::Bool; NErrClustCombs::Int16; subcluster::Int8; BFeas::Int64; interpolating::Bool
  v_sd::T; notplotted::Bool
  confpeak::Vector{T}
  ID✻::Vector{Int64}; ID✻_✻⋂::Vector{Int64}
  anchor::Vector{T}; poles::Vector{T}; numer::Matrix{T}
  CI::Matrix{T}
  peak::NamedTuple{(:X, :p), Tuple{Vector{T}, T}}

	Nobs::Int64; NClustVar::Int8; kX₁::Int64; kX₂::Int64; kY₂::Int64; WREnonARubin::Bool; boottest!::Function
	coldotplus!::Function; colquadformminus!::Function; matmulplus!::Function; panelsum!::Function

  sqrt::Bool; _Nobs::T; kZ::Int64; sumwt::T; haswt::Bool; REst::Bool; multiplier::T; smallsample::T
		dof::Int64; dof_r::T; p::T; BootClust::Int8
		purerobust::Bool; N✻::Int64; Nw::Int64; enumerate::Bool; interpolable::Bool; interpolate_u::Bool; kX::Int64
  _FEID::Vector{Int64}; AR::Matrix{T}; v::Matrix{T}; u✻::Matrix{T}; CT_WE::Matrix{T}
  info✻::Vector{UnitRange{Int64}}; info⋂_✻⋂::Vector{UnitRange{Int64}}
  JN⋂N✻::Matrix{T}; statDenom::Matrix{T}; uXAR::Matrix{T}; SuwtXA::Matrix{T}; numer₀::Matrix{T}; β̂dev::Matrix{T}; δdenom_b::Matrix{T}
	_J⋂::Matrix{T}; YY✻_b::Matrix{T}; YPXY✻_b::Matrix{T}; numerw::Matrix{T}; Zyg::Vector{Matrix{T}}; numer_b::Vector{T}; dist::Matrix{T}

	distCDR::Matrix{T}; plotX::Vector{Vector{T}}; plotY::Vector{T}; ClustShare::Vector{T}; WeightGrp::Vector{UnitRange{Int64}}
  numersum::Vector{T}; ü₀::Vector{T}; invFEwt::Vector{T}
	β̂s::Matrix{T}; As::Matrix{T}
	info✻⋂::Vector{UnitRange{Int64}}; info⋂::Vector{UnitRange{Int64}}; ID✻⋂::Matrix{T}
	ü::Vector{T}
	DGP::StrEstimator{T}; Repl::StrEstimator{T}; M::StrEstimator{T}
	clust::Vector{StrClust{T}}
	denom::Matrix{Matrix{T}}; Kcd::Matrix{Matrix{T}}; Jcd::Matrix{Matrix{T}}; denom₀::Matrix{Matrix{T}}; Jcd₀::Matrix{Matrix{T}}; S✻⋂u₁XinvXX::Matrix{Matrix{T}}; S✻UU::Matrix{Vector{T}}; CTUX::Matrix{Matrix{T}}
	∂u∂r::Vector{Matrix{T}}; ∂numer∂r::Vector{Matrix{T}}; IDCT⋂✻::Vector{Vector{Int64}}; infoCT⋂✻::Vector{Vector{UnitRange{Int64}}}; S✻UX::Vector{Matrix{T}}; S✻UXinvXX::Vector{Matrix{T}}; S✻UZperpinvZperpZperp::Vector{Matrix{T}}; S✻UY::Matrix{Vector{T}}; S✻UMZperp::Vector{Matrix{T}}; S✻UPX::Vector{Matrix{T}}; S✻UZperp::Vector{Matrix{T}}; CTFEU::Vector{Matrix{T}}
  ∂denom∂r::Vector{Matrix{Matrix{T}}}; ∂Jcd∂r::Vector{Matrix{Matrix{T}}}
  ∂²denom∂r²::Matrix{Matrix{Matrix{T}}}
	FEs::Vector{StrFE{T}}
  T1L::Vector{Matrix{T}}; T1R::Vector{Matrix{T}}
	crosstab⋂✻ind::Vector{Int64}; crosstabBootind::Vector{Int64}
  seed::UInt64

	S✻u₁X0::Array{T,3}; ∂S✻u₁X∂β̂::Array{T,3}; ∂S✻u₁X∂γ̈::Array{T,3}; ∂S✻u₁X∂γ̈Π̂::Array{T,3}; ∂S✻u₁X∂r::Array{T,3}; S✻U₂X0::Array{T,3}; ∂S✻U₂X∂Π̂::Array{T,3}
	S✻Zperp_DGPy₁::Array{T,3}; S✻Zperp_DGPZ::Array{T,3}; S✻Zperp_DGPY₂::Array{T,3}; S✻Zperp_DGPX::Array{T,3}; S✻Zperp_DGPZR₁::Array{T,3}
	S✻⋂u₁XinvXX0::Array{T,3}; ∂S✻⋂u₁XinvXX∂β̂::Array{T,3}; ∂S✻⋂u₁XinvXX∂γ̈::Array{T,3}; ∂S✻⋂u₁XinvXX∂γ̈Π̂::Array{T,3}; ∂S✻⋂u₁XinvXX∂r::Array{T,3}; S✻⋂U₂XinvXX0::Array{T,3}; ∂S✻⋂U₂XinvXX0∂Π̂::Array{T,3}
	S✻y₁u₁0::Array{T,3}; ∂S✻y₁u₁∂β̂::Array{T,3}; ∂S✻y₁u₁∂γ̈::Array{T,3}; ∂S✻y₁u₁∂γ̈Π̂::Array{T,3}; ∂S✻y₁u₁∂r::Array{T,3}; S✻y₁U₂0::Array{T,3}; ∂S✻y₁U₂∂r::Array{T,3}; ∂S✻y₁U₂∂β̂::Array{T,3}; ∂S✻y₁U₂∂γ̈Π̂::Array{T,3}; ∂S✻y₁U₂∂γ̈::Array{T,3}; ∂S✻y₁U₂∂Π̂::Array{T,3}
	Reply₁_DGPy₁::Array{T,3}; Reply₁_DGPZR₁::Array{T,3}; Reply₁_DGPZ::Array{T,3}; Reply₁_DGPY₂::Array{T,3}; Reply₁_DGPX::Array{T,3}
	ReplZ_DGPy₁::Array{T,3}; ReplZ_DGPZR₁::Array{T,3}; ReplZ_DGPZ::Array{T,3}; ReplZ_DGPY₂::Array{T,3}; ReplZ_DGPX::Array{T,3}

	_ID✻⋂::Vector{Int}
	StrBootTest{T}(R, r, R₁, r₁, y₁, X₁, Y₂, X₂, wt, fweights, LIML,
	               Fuller, κ, ARubin, B, auxtwtype, rng, maxmatsize, ptype, null, scorebs, bootstrapt, ID, nbootclustvar, nerrclustvar, issorted, robust, small, FEID, FEdfadj, level, rtol, madjtype, NH₀, ML,
								 β̈, A, sc, willplot, gridmin, gridmax, gridpoints, turbo) where T<:Real =

		begin
			kX₂ = ncols(X₂)
			scorebs = scorebs || iszero(B) || ML
			WREnonARubin = !(iszero(kX₂) || scorebs) && !ARubin

			new(R, r, R₁, r₁, y₁, X₁, Y₂, X₂, wt, fweights, LIML || !iszero(Fuller),
					Fuller, κ, ARubin, B, auxtwtype, rng, maxmatsize, ptype, null, bootstrapt, ID, nbootclustvar, nerrclustvar, issorted, small, FEID, FEdfadj, level, rtol, madjtype, NH₀, ML,
					β̈, A, sc, willplot, gridmin, gridmax, gridpoints,
				nrows(R), ptype==symmetric || ptype==equaltail, scorebs, robust || nerrclustvar>0,
				false, false, 0, false, false, 0, 0, 0, false,
				one(T), true,
				[zero(T)],
				Vector{Int64}(undef,0), Vector{Int64}(undef,0),
				Vector{T}(undef,0), Vector{T}(undef,0), Matrix{T}(undef,0,0),
				Matrix{T}(undef,0,0),
				(X = Vector{T}(undef,0), p = T(NaN)),
				nrows(X₁), ncols(ID), ncols(X₁), kX₂, ncols(Y₂), WREnonARubin, WREnonARubin ? boottestWRE! : boottestOLSARubin!,
				turbo ? coldotplus_turbo! : coldotplus_nonturbo!, turbo ? colquadformminus_turbo! : colquadformminus_nonturbo!, turbo ? matmulplus_turbo! : matmulplus_nonturbo!, turbo ? panelsum_turbo! : panelsum_nonturbo!)
		end
end

# cross-tab sum of a column vector w.r.t. given panel info and fixed-effect var
# one row per FE, one col per other grouping
function crosstabFE(o::StrBootTest{T}, v::AbstractVector{T}, info::Vector{UnitRange{Int64}}) where T
  dest = zeros(T, o.NFE, nrows(info))
  if length(info)>0
		@inbounds for i ∈ 1:nrows(info)
	    FEIDi = @view o._FEID[info[i]]
	    vi    = @view       v[info[i]]
	    @simd for j in eachindex(vi, FEIDi)
	  	  dest[FEIDi[j],i] += vi[j]
	    end
	  end
  else  # "robust" case, no clustering
	  @inbounds @simd for i in eachindex(v,o._FEID)
	    dest[o._FEID[i],i] = v[i]
	  end
  end
  dest
end
# same, transposed
function crosstabFEt(o::StrBootTest{T}, v::AbstractVector{T}, info::Vector{UnitRange{Int64}}) where T
  dest = zeros(T, nrows(info), o.NFE)
  if length(info)>0
		@inbounds for i ∈ 1:nrows(info)
	    FEIDi = @view o._FEID[info[i]]
	    vi    = @view       v[info[i]]
	    @simd for j ∈ eachindex(vi, FEIDi)
	  	  dest[i,FEIDi[j]] += vi[j]
	    end
	  end
  else  # "robust" case, no clustering
	  @inbounds @simd for i ∈ eachindex(v,o._FEID)
	    dest[i,o._FEID[i]] = v[i]
	  end
  end
  dest
end
# same, but handles multiple columns in v
# *2nd* dimension of resulting 3-D array corresponds to cols of v;
# this facilitates reshape() to 2-D array in which results for each col of v are stacked vertically
function crosstabFEt(o::StrBootTest{T}, v::AbstractMatrix{T}, info::Vector{UnitRange{Int64}}) where T
  dest = zeros(T, nrows(info), ncols(v), o.NFE)
  if length(info)>0
	  @inbounds for i ∈ 1:nrows(info)
	    FEIDi = @view o._FEID[info[i]]
	    vi    = @view       v[info[i],:]
	    @simd for j ∈ 1:length(FEIDi)
	  	  dest[i,:,FEIDi[j]] += @view vi[j,:]
	    end
	  end
  else  # "robust" case, no clustering
	  @inbounds @simd for i ∈ 1:length(FEIDi)
	    dest[i,:,o._FEID[i]] = @view v[i,:]
	  end
  end
  dest
end

# partial any fixed effects out of a data matrix
function partialFE!(o::StrBootTest, In::AbstractArray)
  if length(In)>0
	  for f ∈ o.FEs
	    tmp = @view In[f.is,:]
	    tmp .-= f.wt'tmp
    end
  end
	nothing
end
function partialFE(o::StrBootTest, In::AbstractArray)
  Out = similar(In)
  if length(In)>0
	  for f ∈ o.FEs
	    tmp = @view In[f.is,:]
	    Out[f.is,:] .= tmp .- f.wt'tmp
	  end
  end
  Out
end

macro storeWtGrpResults!(dest, content)  # poor hygiene in referencing caller's o and w
  if dest == :(o.dist)
		return quote
			if isone($(esc(:o)).Nw)
				$(esc(dest)) .= $(esc(content))
			else
				$(esc(dest))[$(esc(:o)).WeightGrp[$(esc(:w))]] .= $(esc(content))
			end
		end
  else
	  return quote
			local _content = $(esc(content))
	    if isone($(esc(:o)).Nw)
	  	  $(esc(dest)) = _content
	    else
	  	  $(esc(dest))[:,$(esc(:o)).WeightGrp[$(esc(:w))]] = _content
	    end
	  end
  end
end

macro clustAccum!(X, c, Y)  # efficiently add a cluster combination-specific term, factoring in the needed multiplier and sign
  return quote
		local _Y = $(esc(Y))
	  if isone($(esc(c)))
	    if isone($(esc(:o)).clust[1].multiplier)
	  	  $(esc(X)) = $(esc(:o)).clust[1].even ? _Y : -_Y
	    else
	  	  $(esc(X)) = _Y * ($(esc(:o)).clust[1].even ? $(esc(:o)).clust[1].multiplier : -$(esc(:o)).clust[1].multiplier)
	    end
	  elseif $(esc(:o)).clust[$(esc(c))].even
	    if isone($(esc(:o)).clust[$(esc(c))].multiplier)
	  	  $(esc(X)) .+= _Y
	    else
	  	  $(esc(X)) .+= _Y .* $(esc(:o)).clust[$(esc(c))].multiplier
	    end
	  elseif isone($(esc(:o)).clust[$(esc(c))].multiplier)
	    $(esc(X)) .-= _Y
	  else
	    $(esc(X)) .-= _Y .* $(esc(:o)).clust[$(esc(c))].multiplier
	  end
  end
end

function getdist(o::StrBootTest, diststat::DistStatType=nodist)
  if diststat == numer
	  _numer = isone(o.v_sd) ? o.numer : o.numer / o.v_sd
	  o.distCDR = (@view _numer[:,2:end]) .+ o.r
	  sort!(o.distCDR)
  elseif nrows(o.distCDR)==0  # return test stats
    if length(o.dist) > 1
	    o.distCDR = (@view o.dist[2:end])' * o.multiplier
	    sort!(o.distCDR, dims=1)
	  else
	    o.distCDR = zeros(0,1)
	  end
  end
  o.distCDR
end

function sumgreater(x, v)
  dest = zero(Int64)
  @inbounds @simd for i in v
	  x > i && (dest += 1)
  end
  dest
end
function sumless(x, v)
  dest = zero(Int64)
  @inbounds @simd for i in v
	  x < i && (dest += 1)
  end
  dest
end
function sumlessabs(x, v)
  dest = zero(Int64)
  @inbounds @simd for i in v
	  x < abs(i) && (dest += 1)
  end
  dest
end

# store p value in o.p. Return optionally-multiple-hypothesis-adjusted p value. Robust to missing bootstrapped values interpreted as +infinity.
function getp(o::StrBootTest{T}) where T
  o.boottest!(o)
  tmp = o.dist[1]
  isnan(tmp) && return tmp
  if o.B>0
  	if o.sqrt && o.ptype ≠ upper
  	  if o.ptype==symmetric
  	    n = sumlessabs(abs(tmp), o.dist)   # symmetric p value; do so as not to count missing entries in *dist
  	  elseif o.ptype==equaltail
  	    n = 2min(sumgreater(tmp, o.dist) , sumless(tmp, o.dist))
  	  else
  		  n = sumgreater(tmp,  o.dist)  # lower-tailed p value
      end
  	else
  	  n = sumless(tmp, o.dist)  # upper-tailed p value or p value based on squared stats
    end
    o.p = n / o.BFeas |> T
  else
		tmp *= o.multiplier
		o.p = ccdf(o.small ? FDist{T}(T(o.dof), o.dof_r) : Chisq{T}(T(o.dof)), o.sqrt ? tmp^2 : tmp)
		if o.sqrt && !o.twotailed
			o.p /= 2
			(o.ptype==upper) == (tmp<0) && (o.p = 1 - o.p)
		end
  end

	if o.madjtype==bonferroni min(one(T), o.NH₀ * o.p)
  elseif o.madjtype==sidak  one(T) - (one(T) - o.p) ^ o.NH₀
  else o.p
  end
end

getb(o::StrBootTest) = isone(o.v_sd) ? o.numer[:,1] : o.numer[:,1] / o.v_sd  # numerator for full-sample test stat
getV(o::StrBootTest) = o.statDenom / ((isone(o.v_sd) ? o.smallsample : o.v_sd^2 * o.smallsample) * (o.sqrt ? o.multiplier^2 : o.multiplier) * o.dof)  # denominator for full-sample test stat
getv(o::StrBootTest) = @views isone(o.v_sd) ? o.v[:,2:end] : o.v[:,2:end] / o.v_sd  # wild weights
@inline getstat(o::StrBootTest) = o.multiplier * o.dist[1]


text:
include("../utils.jl")

using Manifolds: TFVector, CoTFVector

@testset "Circle" begin
    M = Circle()
    @testset "Real Circle Basics" begin
        @test repr(M) == "Circle(ℝ)"
        @test representation_size(M) == ()
        @test manifold_dimension(M) == 1
        @test !is_point(M, 9.0)
        @test !is_point(M, zeros(3, 3))
        @test_throws DomainError is_point(M, 9.0, true)
        @test_throws DomainError is_point(M, zeros(3, 3), true)
        @test !is_vector(M, 9.0, 0.0)
        @test !is_vector(M, zeros(3, 3), zeros(3, 3))
        @test_throws DomainError is_vector(M, 9.0, 0.0, true)
        @test_throws DomainError is_vector(M, zeros(3, 3), zeros(3, 3), true)
        @test_throws DomainError is_vector(M, 0.0, zeros(3, 3), true)
        @test is_vector(M, 0.0, 0.0)
        @test get_coordinates(M, Ref(0.0), Ref(2.0), DefaultOrthonormalBasis())[] ≈ 2.0
        @test get_coordinates(
            M,
            Ref(0.0),
            Ref(2.0),
            DiagonalizingOrthonormalBasis(Ref(1.0)),
        )[] ≈ 2.0
        @test get_coordinates(
            M,
            Ref(0.0),
            Ref(-2.0),
            DiagonalizingOrthonormalBasis(Ref(1.0)),
        )[] ≈ -2.0
        @test get_coordinates(
            M,
            Ref(0.0),
            Ref(2.0),
            DiagonalizingOrthonormalBasis(Ref(-1.0)),
        )[] ≈ -2.0
        @test get_coordinates(
            M,
            Ref(0.0),
            Ref(-2.0),
            DiagonalizingOrthonormalBasis(Ref(-1.0)),
        )[] ≈ 2.0
        y = [0.0]
        get_coordinates!(M, y, Ref(0.0), Ref(2.0), DiagonalizingOrthonormalBasis(Ref(1.0)))
        @test y ≈ [2.0]
        @test get_vector(M, Ref(0.0), [2.0], DefaultOrthonormalBasis())[] ≈ 2.0
        @test get_vector(M, [0.0], [2.0], DefaultOrthonormalBasis())[] ≈ 2.0
        @test get_vector(M, Ref(0.0), [2.0], DiagonalizingOrthonormalBasis(Ref(1.0)))[] ≈
              2.0
        @test get_vector(M, Ref(0.0), [-2.0], DiagonalizingOrthonormalBasis(Ref(1.0)))[] ≈
              -2.0
        @test get_vector(M, Ref(0.0), [2.0], DiagonalizingOrthonormalBasis(Ref(-1.0)))[] ≈
              -2.0
        @test get_vector(M, Ref(0.0), [-2.0], DiagonalizingOrthonormalBasis(Ref(-1.0)))[] ≈
              2.0
        @test number_of_coordinates(M, DiagonalizingOrthonormalBasis(Ref(-1.0))) == 1
        @test number_of_coordinates(M, DefaultOrthonormalBasis()) == 1
        rrcv = Manifolds.RieszRepresenterCotangentVector(M, 0.0, 1.0)
        @test flat(M, 0.0, 1.0) == rrcv
        @test sharp(M, 0.0, rrcv) == 1.0
        B_cot = Manifolds.dual_basis(M, 0.0, DefaultOrthonormalBasis())
        @test get_coordinates(M, 0.0, rrcv, B_cot) ≈ @SVector [1.0]
        @test get_vector(M, 0.0, 1.0, B_cot) isa Manifolds.RieszRepresenterCotangentVector
        a = fill(NaN)
        get_coordinates!(M, a, 0.0, rrcv, B_cot)
        @test a[] ≈ 1.0
        rrcv2 = Manifolds.RieszRepresenterCotangentVector(M, fill(0.0), fill(-10.0))
        get_vector!(M, rrcv2, 0.0, 1.0, B_cot)
        @test isapprox(rrcv.X, rrcv2.X[])
        @test vector_transport_to(M, 0.0, 1.0, 1.0, ParallelTransport()) == 1.0
        @test retract(M, 0.0, 1.0) == exp(M, 0.0, 1.0)
        @test injectivity_radius(M) ≈ π
        @test injectivity_radius(M, Ref(-2.0)) ≈ π
        @test injectivity_radius(M, Ref(-2.0), ExponentialRetraction()) ≈ π
        @test injectivity_radius(M, ExponentialRetraction()) ≈ π
        @test mean(M, [-π / 2, 0.0, π]) ≈ -π / 2
        @test mean(M, [-π / 2, 0.0, π], [1.0, 1.0, 1.0]) == -π / 2
        z = project(M, 1.5 * π)
        z2 = fill(0.0)
        project!(M, z2, 1.5 * π)
        @test z2[1] == z
        @test project(M, z) == z
        @test project(M, 1.0, 2.0) == 2.0
    end
    TEST_STATIC_SIZED && @testset "Real Circle and static sized arrays" begin
        X = @MArray fill(0.0)
        p = @SArray fill(0.0)
        log!(M, X, p, @SArray fill(π / 4))
        @test norm(M, p, X) ≈ π / 4
        @test is_vector(M, p, X)
        @test is_vector(M, [], X)
        @test project(M, 1.0) == 1.0
        p = @MArray fill(0.0)
        project!(M, p, p)
        @test p == @MArray fill(0.0)
        p .+= 2 * π
        project!(M, p, p)
        @test p == @MArray fill(0.0)
        @test project(M, 0.0, 1.0) == 1.0
    end
    types = [Float64]
    TEST_FLOAT32 && push!(types, Float32)

    basis_types = (DefaultOrthonormalBasis(),)
    basis_types_real = (
        DefaultOrthonormalBasis(),
        DiagonalizingOrthonormalBasis(Ref(-1.0)),
        DiagonalizingOrthonormalBasis(Ref(1.0)),
    )
    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [-π / 4, 0.0, π / 4])
            test_manifold(
                M,
                pts,
                test_vector_spaces=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                test_vee_hat=false,
                is_mutating=false,
                test_rand_point=true,
                test_rand_tvector=true,
            )
            ptsS = map(p -> (@SArray fill(p)), pts)
            test_manifold(
                M,
                ptsS,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                test_vee_hat=true,
                basis_types_vecs=basis_types_real,
                basis_types_to_from=basis_types_real,
                test_rand_point=true,
                test_rand_tvector=true,
            )
        end
    end
    @testset "Mutating Rand for real Circle" begin
        p = fill(NaN)
        X = fill(NaN)
        rand!(M, p)
        @test is_point(M, p)
        rand!(M, X; vector_at=p)
        @test is_vector(M, p, X)

        rng = MersenneTwister()
        rand!(rng, M, p)
        @test is_point(M, p)
        rand!(rng, M, X; vector_at=p)
        @test is_vector(M, p, X)
    end
    @testset "Test sym_rem" begin
        p = 4.0 # not a point
        p = sym_rem(p) # modulo to a point
        @test is_point(M, p)
    end
    Mc = Circle(ℂ)
    @testset "Complex Circle Basics" begin
        @test repr(Mc) == "Circle(ℂ)"
        @test representation_size(Mc) == ()
        @test manifold_dimension(Mc) == 1
        @test is_vector(Mc, 1im, 0.0)
        @test is_point(Mc, 1im)
        @test !is_point(Mc, 1 + 1im)
        @test_throws DomainError is_point(Mc, 1 + 1im, true)
        @test !is_vector(Mc, 1 + 1im, 0.0)
        @test_throws DomainError is_vector(Mc, 1 + 1im, 0.0, true)
        @test !is_vector(Mc, 1im, 2im)
        @test_throws DomainError is_vector(Mc, 1im, 2im, true)
        rrcv = Manifolds.RieszRepresenterCotangentVector(Mc, 0.0 + 0.0im, 1.0im)
        @test flat(Mc, 0.0 + 0.0im, 1.0im) == rrcv
        @test sharp(Mc, 0.0 + 0.0im, rrcv) == 1.0im
        @test norm(Mc, 1.0, log(Mc, 1.0, -1.0)) ≈ π
        @test is_vector(Mc, 1.0, log(Mc, 1.0, -1.0))
        X = @MArray fill(0.0 + 0.0im)
        p = @SArray fill(1.0 + 0.0im)
        log!(Mc, X, p, @SArray fill(-1.0 + 0.0im))
        @test norm(Mc, (@SArray fill(1.0)), X) ≈ π
        @test is_vector(Mc, p, X)
        @test project(Mc, 1.0) == 1.0
        @test project(Mc, 1 / sqrt(2.0) + 1 / sqrt(2.0) * im) ≈
              1 / sqrt(2.0) + 1 / sqrt(2.0) * im
        p = @MArray fill(1.0 + 0.0im)
        project!(Mc, p, p)
        @test p == @MArray fill(1.0 + 0.0im)
        p .*= 2
        project!(Mc, p, p)
        @test p == @MArray fill(1.0 + 0.0im)

        angles = map(pp -> exp(pp * im), [-π / 2, 0.0, π])
        @test mean(Mc, angles) ≈ exp(-π * im / 2)
        @test mean(Mc, angles, [1.0, 1.0, 1.0]) ≈ exp(-π * im / 2)
        @test_throws ErrorException mean(Mc, [-1.0 + 0im, 1.0 + 0im])
        @test_throws ErrorException mean(Mc, [-1.0 + 0im, 1.0 + 0im], [1.0, 1.0])
    end
    types = [Complex{Float64}]
    TEST_FLOAT32 && push!(types, Complex{Float32})

    for T in types
        @testset "Type $T" begin
            a = 1 / sqrt(2.0)
            pts = convert.(Ref(T), [a - a * im, 1 + 0im, a + a * im])
            test_manifold(
                Mc,
                pts,
                test_vector_spaces=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                is_mutating=false,
                test_vee_hat=false,
                exp_log_atol_multiplier=2.0,
                is_tangent_atol_multiplier=2.0,
            )
            ptsS = map(p -> (@SArray fill(p)), pts)
            test_manifold(
                Mc,
                ptsS,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                test_vee_hat=true,
                exp_log_atol_multiplier=2.0,
                is_tangent_atol_multiplier=2.0,
                basis_types_vecs=basis_types,
                basis_types_to_from=basis_types,
            )
        end
    end
end

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
const I2 = SMatrix{2,2,Float64}(1, 0, 0, 1)

function check_points(topo::Symbol, p0::AbstractVector{MVector{2,Float64}})
    if topo == :O || topo == :C
        @assert p0[1] === p0[end] "the first and last node must be the same object for `:O` and `:C` type grids"
    end
    np = length(p0)
    for i in 2:np
        @assert !isapprox(p0[i],p0[i - 1],atol=1e-10)  "duplicated nodes $(i-1) and $(i) (distance < 1e-10)"
    end
    return nothing
end

function gen_steps(ss::StepSize, height::Float64, mg_level::Int)::Vector{Float64}
    i::Int = 0
    h::Float64 = 0.0
    s = Vector{Float64}()
    while h < height
        i += 1
        push!(s, ss(i))
        h += s[end]
    end
    if mg_level > 0
        mg_step = 2^mg_level
        ns1 = ceil(Int, length(s)//mg_step) * mg_step
        dnp = ns1 - length(s)
        append!(s, fill(s[end], dnp))
    end
    return s
end

@generated function diag_add(a::A, val::T) where {N, T, A<:MMatrix{N,N,T}}
    quote
        $(map(i->:(@inbounds a[$i,$i] += val), 1:N)...)
        :(return a)
    end
end

@generated function diag_sub(a::A, val::T) where {N, T, A<:MMatrix{N,N,T}}
    quote
        $(map(i->:(@inbounds a[$i,$i] -= val), 1:N)...)
        :(return a)
    end
end

function create_wake(
    w::WakeInfo, surf::V; mg_level::Int=0
) where {V<:AbstractVector{MVector{2,Float64}}}
    ds0 = norm(surf[1] - surf[2])
    ds1 = norm(surf[end] - surf[end - 1])
    @assert isapprox(ds0, ds1, rtol=0.5)
    ds = sqrt(ds0 * ds1)
    ss = ExpStepSize(ds, w.ratio)
    if w.dsmax > 0
        ss = MinStepSize(ss, w.dsmax)
    end
    s = gen_steps(ss, w.len, mg_level)
    dir = SVector(cosd(w.aoa), sind(w.aoa))
    ns = length(s)
    points = Vector{MVector{2,Float64}}(undef, ns)
    points[1] = surf[end] + s[1]*dir
    for i in 2:ns
        points[i] = points[i-1] + s[i]*dir
    end

    return points
end

join_wake(surf, wake) = vcat(reverse(wake), surf, wake)


text:
export isnorm, norm_equation

#TODO: verbose printing

function norm_1_generators(A::Vector{NfOrdIdl})
  @assert all(isprime, A)
  @assert all(x->x.gen_one == A[1].gen_one, A)

  f = matrix(FlintZZ, 1, length(A), [degree(x) for x = A])
  k = nullspace(f)[2]

  id = [FacElem(A, [k[i,j] for i=1:length(A)]) for j=1:ncols(k)]
  return id
end

@doc Markdown.doc"""
    norm_equation(K::AnticNumerField, a) -> nf_elem

For $a$ an integer or rational, try to find $T \in K$ s.th.
$N(T) = a$. Raises an error if unsuccessful.
"""
function norm_equation(K::AnticNumberField, a)
  fl, s = isnorm(K, a)
  if fl
    return evaluate(s)
  end
  error("no solution")
end

@doc Markdown.doc"""
    isnorm(K::AnticNumberField, a) -> Bool, nf_elem

For $a$ an integer or rational, try to find $T \in K$ s.th. $N(T) = a$
holds. If successful, return true and $T$, otherwise false and some element.
The element will be returned in factored form.
"""
function isnorm(K::AnticNumberField, a::Integer)
  return isnorm(K, fmpz(a))
end
function isnorm(K::AnticNumberField, a::fmpq)
  fl, s = isnorm(K, numerator(a)*denominator(a)^(degree(K)-1))
  return fl, s * FacElem(Dict(K(denominator(a)) => fmpz(-1)))
end
function isnorm(K::AnticNumberField, a::Rational)
  return isnorm(K, fmpq(a))
end

@doc Markdown.doc"""
    isnorm(K::AnticNumberField, a::fmpz; extra::Vector{fmpz}) -> Bool, nf_elem

For a fmpz $a$, try to find $T \in K$ s.th. $N(T) = a$
holds. If successful, return true and $T$, otherwise false and some element.
In \testtt{extra} one can pass in additional prime numbers that
are allowed to occur in the solution. This will then be supplemented.
The element will be returned in factored form.
"""
function isnorm(K::AnticNumberField, a::fmpz; extra::Vector{fmpz}=fmpz[])
  L = lll(maximal_order(K))
  C, mC = narrow_class_group(L)
#  println("narrow group is : $C")
  S = union(Set(keys(factor(a).fac)), Set(keys(factor(discriminant(L)).fac)))
  S = union(S, Set(extra))

  g = Set(elem_type(C)[])
  for p = S
    P = prime_ideals_over(L, typeof(p)[p])
    P1 = norm_1_generators(P)
    for x = P1
      push!(g, sum([e*preimage(mC, k) for (k, e) = x.fac]))
    end
  end

#  println("found $g")
  sC, msC = sub(C, collect(g))

  #think: might tbe better using quo, test for zero and a function
  #       that adds relations to a quo.
  if order(sC) != order(C)
    p = 2
    stable = 0
    while true
      while p in S
        p = next_prime(p)
      end
      change = false
      P1 =  norm_1_generators(prime_ideals_over(L, [p]))
      for x = P1
        y = sum([e*preimage(mC, k) for (k, e) = x.fac]) # move to log function(?)
        if !haspreimage(msC, y)[1]
          change = true
          push!(g, y)
          push!(S, p)
          sC, msC = sub(C, collect(g))
          if order(sC) == order(C)
#            println("groups are the same")
            break
          end
        end
      end
      p = next_prime(p)
      if change
        stable = 0
      else
        stable += 1
        if stable > 5
          break
        end
      end
    end
  else
#    println("groups are the same")
  end

  SP = prime_ideals_over(L, [x for x = S])

  U, mU = sunit_group_fac_elem(SP)
  u, mu = sunit_group_fac_elem(collect(S))

  h = elem_type(u)[]
  for i=1:ngens(U)
    x = mU(U[i])
    y = norm(x)
    push!(h, preimage(mu, y))
    #@assert norm(evaluate(x)) == evaluate(image(mu, h[end]))
  end
  s, ms = sub(u, h)
  mp = GrpAbFinGenMap(U, u, vcat([x.coeff for x=h]))

  fl, p = haspreimage(mp, preimage(mu, a))
  if fl
    return true, mU(p)
  else
    return false, FacElem(one(K))
  end
end


text:
# # Compute the solitary wave of the Whitham equation with prescribed velocity
export PlotSolitaryWaveKdV,PlotSolitaryWaveWhitham
using WaterWaves1D,FFTW,Plots;

#---- KdV
"""
	PlotSolitaryWaveKdV()

A proof of concept: numerically compute the solitary wave of the KdV equation with prescribed velocity `c=5`
(augmenting progressively the velocity to generate initial guesses for the iterative scheme)
and then computes the difference with the exact solution.
"""
function PlotSolitaryWaveKdV(c = 5)
	param = ( μ  = 1,
			ϵ  = 1,
        	N  = 2^10,
            L  = 10,
			)
	mesh = Mesh(param)

	function sol(x,α)
		2*α*sech.(sqrt(3*2*α)/2*x).^2
	end

	u₀, = SolitaryWaveWhitham( merge(param,(c=2.5,)); guess = sol(mesh.x,1.5)+.1*exp.(-(mesh.x).^2), iterative = false, KdV = true, α=1,tol = 1e-8)
	u₁, = SolitaryWaveWhitham( merge(param,(c=2.6,)); guess = sol(mesh.x,1.6)+.1*exp.(-(mesh.x).^2), iterative = false, KdV = true, α=1,tol = 1e-8)
	u₂, = SolitaryWaveWhitham( merge(param,(c=2.7,)); guess = sol(mesh.x,1.7)+.5*exp.(-(mesh.x).^2), iterative = false, KdV = true, α=1,tol = 1e-8)

	# barycentric Lagrange interpolation (4.2) in https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
	function P(X,dc,u₀,u₁,u₂)
		(u₀/(2*(X/dc+2))-u₁/(X/dc+1)+u₂/(2*X/dc))/(1/(2*(X/dc+2))-1/(X/dc+1)+1/(2*X/dc))
	end

	function solve!(u₀,u₁,u₂,c,dc)
		tu₂ = copy(u₂)
		u₂ .= SolitaryWaveWhitham( merge(param,(c=c,)); guess = P.(dc,dc,u₀,u₁,u₂), iterative = false, α=1, KdV = true, verbose = true, max_iter = 5, tol = 1e-10)[1]
		u₀ .= u₁
		u₁ .= tu₂
		end

		for cs = range(2.8; step = 0.1, stop = c)
			solve!(u₀,u₁,u₂,cs,0.1)
		end

		ucomp = u₂
		uexact = sol(mesh.x,c-1)
		error = maximum(abs.(ucomp-uexact))
		@info("Final error : $error")
		plt = plot(layout=(2,2))

		plot!(plt[1,1], mesh.x, [ucomp uexact];
			title=string("c=$c"),
			label=["computed" "exact"])

		plot!(plt[2,1], fftshift(mesh.k),
			[log10.(abs.(fftshift(fft(ucomp)))) log10.(abs.(fftshift(fft(uexact))))];
			title="frequency",
			label=["computed" "exact"])

		plot!(plt[1,2], mesh.x, ucomp-uexact;
			title=string("c=$c"),
			label="difference")

		plot!(plt[2,2], fftshift(mesh.k),
			log10.(abs.(fftshift(fft(ucomp-uexact))));
			title="frequency",
			label="difference")
			display(plt)
		return error
	end

#----
"""
	`PlotSolitaryWaveWhitham(c)`

Compute the solitary wave of the Whitham equation with prescribed velocity.

`c` is the velocity, and should be more between `1` and `1.2290408`.
"""
function PlotSolitaryWaveWhitham(c=1.1)
 	function solKdV(x,α)
   	 	2*α*sech.(sqrt(3*2*α)/2*x).^2
 	end
 	if c<=1 || c>1.2290408
	 	error("The velocity should be between `1` and `1.2290408`")
 	elseif c<1.15
		param = ( μ  = 1,
 				ϵ  = 1,
         		N  = 2^10,
            	L  = 10/sqrt(c-1) )
 		mesh = Mesh(param)
		u, = SolitaryWaveWhitham( merge(param,(c=c,)); iterative = false)

 	else
		param = ( μ  = 1,
			ϵ  = 1,
        	N  = 2^10,
            L  = 15,
			)
		mesh = Mesh(param)

		# barycentric Lagrange inerpolation (4.2) in https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
		function P(X,dc,u₀,u₁,u₂)
			(u₀/(2*(X/dc+2))-u₁/(X/dc+1)+u₂/(2*X/dc))/(1/(2*(X/dc+2))-1/(X/dc+1)+1/(2*X/dc))
		end

		function solve2!(u₀,u₁,u₂,c,dc)
			tu₂ = copy(u₂)
			u₂ .= SolitaryWaveWhitham(merge(param,(c=c,)); guess = P.(dc,dc,u₀,u₁,u₂), iterative = false, verbose = true, max_iter = 5)[1]
			u₀ .= u₁
			u₁ .= tu₂
		end

		function solve2!(u₀,u₁,u₂,c,old_dc,new_dc)
			tu₀ = copy(u₀)
			tu₁ = copy(u₁)
			tu₂ = copy(u₂)
			u₀ .= tu₂
			u₁ .= SolitaryWaveWhitham( merge(param,(c=c+1*new_dc,)); guess = P.(1*new_dc,old_dc,tu₀,tu₁,tu₂), iterative = false, verbose = true)[1]
			u₂ .= SolitaryWaveWhitham( merge(param,(c=c+2*new_dc,)); guess = P.(2*new_dc,old_dc,tu₀,tu₁,tu₂), iterative = false, verbose = true)[1]
		end

		u₀, = SolitaryWaveWhitham( merge(param,(c=1.10,)); iterative = false)
		u₁, = SolitaryWaveWhitham( merge(param,(c=1.11,)); iterative = false)
		u₂, = SolitaryWaveWhitham( merge(param,(c=1.12,)); iterative = false)

		if c < 1.22
			for cs = range(1.13; step = 0.01, stop = c)
				solve2!(u₀,u₁,u₂,cs,0.01)
			end
		else
			for cs = range(1.13; step = 0.01, stop = 1.22)
				print(string("c = ",cs,"\n"))
				solve2!(u₀,u₁,u₂,cs,0.01)
			end

			solve2!(u₀,u₁,u₂,1.221,0.01,1e-4)

			if c < 1.229
				for cs = range(1.221+3e-4; step = 1e-4, stop = c)
					solve2!(u₀,u₁,u₂,cs,1e-4)
				end
			else
				c=1.2290408
				for cs = range(1.221+3e-4; step = 1e-4, stop = 1.229)
					solve2!(u₀,u₁,u₂,cs,1e-4)
				end

				solve2!(u₀,u₁,u₂,1.229,1e-4,1e-6)

				for cs = range(1.229+3e-6; step = 1e-6, stop = 1.22904)
					solve2!(u₀,u₁,u₂,cs,1e-6)
				end

				solve2!(u₀,u₁,u₂,1.22904,1e-6,1e-8)

				for cs = range(1.22904+3e-8; step = 1e-8, stop = 1.2290407)
					solve2!(u₀,u₁,u₂,cs,1e-8)
				end
				for cs = range(1.2290407+1e-8; step = 1e-8, stop = 1.2290408)
					solve2!(u₀,u₁,u₂,cs,1e-8)
				end
			end
		end
		u = u₂
 	end
 	plt = plot(layout=(1,2))
 	plot!(plt[1,1], mesh.x, [u solKdV(mesh.x,c-1)];
  		title=string("c=",c),
  		label=["Whitham" "KdV"])

 	plot!(plt[1,2], fftshift(mesh.k),
  		log10.(abs.(fftshift(fft(u))));
  		title="frequency",
  		label="Whitham")

  	return 	maximum(abs.(u-solKdV(mesh.x,c-1)))
end


text:
using DataFrames
using Gadfly
using Lumberjack
using Docile

"""
Abstract:  parses a bismark alignment report for the alignment statistics
Arguments: bismark_report_path - string to the bismark report path
Return:    Dictionary{String,Number}
              keys are string description of the alignment statistic
"""
function parse_bismark_report(bismark_report_path::String)
    bismark_line_itr = eachline(open(bismark_report_path))
    tab_lines = filter(x->ismatch(r"\t",x), collect(bismark_line_itr) ) # if there is a tab keep it
    rgx=r"([A-Za-z\'\s]+):\t([0-9\.]+)"
    report_matches = map(x-> match(rgx,x), tab_lines )
    report_dict=Dict()
    for report_match in report_matches
       statistic_desc  = report_match.captures[1]
       statistic_value = report_match.captures[2]
       report_dict[ statistic_desc   ]= parsefloat( statistic_value)
    end
    return report_dict
end

"""
Abstract:  parses a set of bismark report files
Arguments:
    bismark_report_filenames: Vector of strings that are path names to bismark report files
Returns:
    Dictionary; key: statistic name, value: Vectors of Numbers; Each Vector holds the values for each sample in sampleinfo
"""
function parse_bismark_reports{T<:String}(bismark_report_paths::Vector{T} )
  master_report_dict=Dict()
  for bismark_report_path in bismark_report_paths
      report_dict = parse_bismark_report(bismark_report_path)
      Lumberjack.info("parsing: $bismark_report_path")
      for (key,value) in report_dict
          #assign to master dict, and push on number value
          if !haskey(master_report_dict,key)
             master_report_dict[key] = Any[]
          end
          push!( master_report_dict[key], value)
      end
  end
  return master_report_dict
end

"""
Abstract: appends bismark report statistics to sampleinfo
Arguments:
    sampleinfo   - DataFrame of samples information
    report_dict  - Dictionary{String, Vector{Number}}; Key: statistic name, Value: statistic values for each sample in sampleinfo
"""
function append_report_info_to_sampleinfo!(sampleinfo::DataFrame,report_dict::Dict)
    nrow_sampleinfo = nrow(sampleinfo)
    for (key, values) in report_dict
       num_values = length(values)
       if nrow_sampleinfo == num_values
           sampleinfo[symbol(key) ]=values
       else
           Lumberjack.error("Cannot add $key to metadata, number of values ($num_values) does not equal metadata rows ($nrow_sampleinfo) ")
       end
    end
end

# bismark_stats are the bismark outputs such as mapping efficiency
# datasource_groupby_field: e.g. strain,sex
# files - the filenames to the bismark report file
# for each of the bismark stats, create a box plot based on each of the
#  datasource_groupby_fields
"""
Abstract: Creates basic gadfly box plots of bismark statistics over various sample attributes
e.g. sex, strain, or whatever other condition is described in the sample information.
Arguments::
    sampleinfo                - DataFrame of sample information and bismark statistics
    bismark_stats             - Vector of names of the bismark statistic
    sampleinfo_groupby_fields - Vector of fields that bismark_stats should be plotted over
"""
function plot_bismark_summary{T<:String}(sampleinfo::DataFrame, bismark_stats::Vector{T},
                                         sampleinfo_groupby_fields::Vector{T})

    plots = Any[]
    theme = Theme(
         major_label_color=color("black"),
         minor_label_color=color("black"),
         major_label_font_size=18pt,
         minor_label_font_size=14pt
    )
    for statistic in bismark_stats
        for groupby_field in sampleinfo_groupby_fields
            p = Gadfly.plot( sampleinfo, x=groupby_field, y=statistic,Geom.boxplot, theme )
            push!(plots,p)
        end
    end
    return plots
end

"""
Abstract: bismark summary report.
Arguments:
    reportdir                 - path
    bismark_report_paths      - Vector of paths to bismark report files
    sampleinfo_path           - path to datasource
"""
function bismark_summary_report{T<:String}(reportdir::String, bismark_report_paths::Vector{T},
                                        sampleinfo_path::String, sampleinfo_groupby_fields::Vector{T})
    report_dict = parse_bismark_reports(bismark_report_paths)
    println("REPORT_DICT: $report_dict")
    sampleinfo = readtable(sampleinfo_path,separator='\t',header=true)
    append_report_info_to_sampleinfo!(sampleinfo,report_dict)
    return sampleinfo
end

"""
Abstract: Produce bismark summary report.
Arguments:
    reportdir                 - path
    bismark_report_paths      - Vector of paths to bismark report files
    sampleinfo_path           - path to datasource
"""
function make_bismark_summary_report{T<:String}(reportdir::String, bismark_report_paths::Vector{T},
                                        sampleinfo_path::String, sampleinfo_groupby_fields::Vector{T};
                                        overwrite=false,
                                        bismark_stats=["Sequence pairs analysed in total","Mapping efficiency"],
                                        reportfilename="datasource.tsv",
                                        plotfilename="plots.png")

    sampleinfo = bismark_summary_report(reportdir,bismark_report_paths,sampleinfo_path,sampleinfo_groupby_fields)
    plots=plot_bismark_summary(sampleinfo, bismark_stats,sampleinfo_groupby_fields)

    # if directory does not exist create it
    if(!isdir( reportdir))
       mkdir(reportdir)
    else
      Lumberjack.warn("$reportdir exists, overwrite is set to $overwrite")
      if overwrite
         error("Overwriting report directory disallowed.")
      end
    end

    # write sampleinfo+bismark_stats to directory
    datatable_path = joinpath(reportdir,reportfilename)
    writetable(datatable_path,sampleinfo) # will automatically use tabs as separator if .tsv
    println("DATATABLE SAVED TO $datatable_path")
    # draw plots to directory
    plotpath = joinpath(reportdir,plotfilename)
    Gadfly.draw( PNG(plotpath,24cm, 12cm)  , hstack(plots...) )
    println("PLOTS SAVED TO $plotpath")
end

# shift to Misc
# For each row in the datasource extract information into a dataframe
function get_bismark_report_filenames_from_plite_pipeline(pipeline_path,run_number,num_jobs,step_name="bismark",report_name="file1_1_val_1.fq_bismark_bt2_PE_report.txt")
    files = map( x->joinpath(pipeline_path, "output","run$run_number", "job$x",step_name, report_name), 0:(num_jobs-1) )
    return files
end


text:
# the IncrementalInference API

# DFG SpecialDefinitions
export AbstractDFG,
  InMemDFGType,
  getSolverParams,
  LightDFG,
  findShortestPathDijkstra, isPathFactorsHomogeneous,
  getSolvedCount, isSolved, setSolvedCount!,
  listSupersolves, listSolveKeys,
  deepcopySolvekeys!, deepcopySupersolve!,
  diagm,
  listDataEntries,
  FolderStore,
  addBlobStore!,
  addData!,
  getData,
  DFGVariable,
  DFGVariableSummary,
  DFGFactor,
  DFGFactorSummary,
  deleteVariableSolverData!
  # listDataBlobs  # ERROR: LightDFG{} doesn't override 'listDataBlobs'.

# Inference types
export FunctorInferenceType, PackedInferenceType
export AbstractPrior, AbstractRelative, AbstractRelativeRoots, AbstractRelativeMinimize

# not sure if this is necessary
export convert, *

export
  CSMHistory,
  # getTreeCliqsSolverHistories,

  AbstractBayesTree,
  BayesTreeNodeData,
  PackedBayesTreeNodeData,

  # state machine methods
  StateMachine,
  exitStateMachine,
  print,
  getGraphFromHistory,
  getCliqSubgraphFromHistory,
  sandboxStateMachineStep,
  # draw and animate state machine
  getStateLabel,
  histStateMachineTransitions,
  histGraphStateMachineTransitions,
  drawStateTransitionStep,
  drawStateMachineHistory,
  animateStateMachineHistoryByTime,
  animateStateMachineHistoryByTimeCompound,
  animateCliqStateMachines,
  makeCsmMovie,
  areSiblingsRemaingNeedDownOnly,

  # general types for softtyping of variable nodes
  BeliefArray,
  InferenceVariable,
  ContinuousScalar,
  SamplableBelief,
  PackedSamplableBelief,
  Prior,
  PackedPrior,
  MsgPrior,
  PackedMsgPrior,
  PartialPrior,
  PackedPartialPrior,

  ls2,
  # lsRear,
  # from DFG
  ls,
  lsf,
  listVariables,
  listFactors,
  exists,
  sortDFG,
  getLabel,
  getVariables,
  getVariableOrder,
  getPPE,
  getPPEDict,
  getVariablePPE,
  isVariable,
  isFactor,
  getFactorType,
  getSofttype,
  getVariableType,
  getLogPath,
  joinLogPath,
  lsfPriors,
  isPrior,
  lsTypes,
  lsfTypes,
  findClosestTimestamp,
  printVariable,
  printFactor,
  getTimestamp,
  deepcopyGraph,
  deepcopyGraph!,
  copyGraph!,
  getSolverData,

  # using either dictionary or cloudgraphs
  FactorMetadata,
  FunctionNodeData,
  PackedFunctionNodeData, # moved to DFG
  normalfromstring,
  categoricalfromstring,
  # extractdistribution,

  SolverParams,
  getSolvable,
  setSolvable!,
  addVariable!,
  deleteVariable!,
  addFactor!,
  deleteFactor!,
  addMsgFactors!,
  deleteMsgFactors!,
  factorCanInitFromOtherVars,
  doautoinit!,
  initManual!,
  initVariableManual!,
  resetInitialValues!,
  resetInitValues!,
  # asyncTreeInferUp!,
  # initInferTreeUp!,
  solveCliqWithStateMachine!,
  resetData!,
  resetTreeCliquesForUpSolve!,
  resetFactorGraphNewTree!,
  setVariableInitialized!,
  setVariableInferDim!,
  resetVariable!,
  getFactor,
  getFactorDim,
  getVariableDim,
  getVariable,
  getCliqueData,
  setCliqueData!,
  getManifold,  # new Manifolds.jl based operations
  getVal,
  getBW,
  setVal!,
  getNumPts,
  getBWVal,
  setBW!,
  setValKDE!,
  buildCliqSubgraph,

  #
  isPartial,
  isInitialized,
  isTreeSolved,
  isUpInferenceComplete,
  isCliqInitialized,
  isCliqUpSolved,
  areCliqVariablesAllInitialized,
  ensureSolvable!,
  initAll!,
  cycleInitByVarOrder!,
  BayesTree, MetaBayesTree,
  TreeBelief,
  LikelihoodMessage,
  initfg,
  buildSubgraph,
  buildCliqSubgraph!,
  transferUpdateSubGraph!,
  getEliminationOrder,
  buildBayesNet!,
  buildTree!,
  buildTreeReset!,
  buildCliquePotentials,

  getCliqDepth,
  getTreeAllFrontalSyms,
  getTreeCliqUpMsgsAll,
  childCliqs,
  getChildren,
  parentCliq,
  getParent,
  getCliqSiblings,
  getNumCliqs,
  getBelief,
  CliqStateMachineContainer,

  solveCliqUp!,
  solveCliqDown!,
  fifoFreeze!,

  #functors need
  getSample,
  sampleFactor!,
  sampleFactor,

  #Visualization
  drawGraph,
  drawGraphCliq,
  drawCliqSubgraphUpMocking,
  drawTree,
  drawTreeAsyncLoop,

  # Bayes (Junction) Tree
  evalFactor,
  approxConvBelief,
  approxConv,

  # more debugging tools
  localProduct,
  treeProductUp,
  approxCliqMarginalUp!,
  dontMarginalizeVariablesAll!,
  unfreezeVariablesAll!,
  resetVariableAllInitializations!,
  isMarginalized,
  setMarginalized!,
  isMultihypo,
  getMultihypoDistribution,
  getHypothesesVectors,

  # weiged sampling
  AliasingScalarSampler,
  rand!,
  rand,
  fastnorm,

  # new wrapper (experimental)
  CommonConvWrapper,

  getCliqVarInitOrderUp,
  getCliqNumAssocFactorsPerVar,

  # introduced for approximate convolution operations
  setThreadModel!,
  SingleThreaded,
  MultiThreaded,

  # user functions
  predictbelief,
  propagateBelief,
  getCliqMat,
  getCliqAssocMat,
  getCliqMsgMat,
  getCliqFrontalVarIds,
  getFrontals,
  getCliqSeparatorVarIds,
  getCliqAllVarIds,
  getCliqVarIdsAll,
  getCliqAllVarSyms,
  getCliqVarIdsPriors,
  getCliqVarSingletons,
  getCliqAllFactIds,
  getCliqFactorIdsAll,
  getCliqFactors,
  areCliqVariablesAllMarginalized,

  # generic marginal used during elimitation game
  GenericMarginal,
  PackedGenericMarginal,

  # factor graph operating system utils (fgos)
  saveTree,
  loadTree,

  # Temp placeholder for evaluating string types to real types
  saveDFG,
  loadDFG!,  loadDFG,
  rebuildFactorMetadata!,

  getCliqVarSolveOrderUp,
  getFactorsAmongVariablesOnly,
  setfreeze!,

  #internal dev functions for recycling cliques on tree
  attemptTreeSimilarClique,

  # some utils
  compare,
  compareAllSpecial,
  getMeasurements,
  findFactorsBetweenFrom,
  addDownVariableFactors!,
  getDimension,
  getPointType,
  getPointIdentity,
  setVariableRefence!,
  reshapeVec2Mat

export  buildCliqSubgraph_StateMachine


text:
# Vertex and Edge types; useful for wrapping an existing type with the edge interface.
# Note that DirectedGraph does not require using these types; just implement the edge interface.
for typename in (:Edge, :Vertex)
    getid = Symbol(lowercase(string(typename)) * "_id")
    setid = Symbol("set_" * lowercase(string(typename)) * "_id!")
    @eval begin
        mutable struct $typename{T}
            data::T
            id::Int64
        end
        $typename(data) = $typename(data, -1)
        $getid(x::$typename) = x.id
        $setid(x::$typename, index::Int64) = (x.id = index)
    end
end


text:
"""
This module implements a value iteration solver that uses the interface defined in POMDPs.jl
"""
module DiscreteValueIteration

using Printf
using POMDPs
using POMDPModelTools
using POMDPPolicies
using SparseArrays
import POMDPLinter: @POMDP_require, @req, @subreq, @warn_requirements

import POMDPs: Solver, solve, Policy, action, value

export
    ValueIterationPolicy,
    ValueIterationSolver,
    SparseValueIterationSolver,
    solve,
    action,
    value,
    locals

include("common.jl")
include("vanilla.jl")
include("sparse.jl")

end # module


text:
using PQPolygonSDK
using Dates

# First build an options contract model -
contract_model_options = Dict{String,Any}()
contract_model_options["underlying"] = "MU"
contract_model_options["expiration"] = Date(2022, 09, 16)
contract_model_options["strike"] = 55.0
put_contract_model = model(PolygonPutOptionContractModel, contract_model_options)

# what is the ticker symbol for this option?
options_ticker = ticker(put_contract_model)

# ok, now that we have the options ticker, let's do a hit to see some data -
# build a user model -
options = Dict{String,Any}()
options["email"] = "jvarner@paliquant.com"
options["apikey"] = "abc1234" # do _not_ check in a real API key

# build the user model -
user_model = model(PQPolygonSDKUserModel, options)

# build the endpoint model -
endpoint_options = Dict{String,Any}()
endpoint_options["ticker"] = options_ticker
endpoint_model = model(PolygonOptionsLastTradeEndpointModel, user_model, endpoint_options)

# build the url -
base = "https://api.polygon.io"
my_url_string = url(base, endpoint_model)

# make the call -
(header, data) = api(PolygonOptionsLastTradeEndpointModel, my_url_string)


text:
using JuliaFormatter

format_text("""

module X

struct Foo
    x::Int
end

function foo(x::Foo)
    return nothing
end

bar() = foo(Foo(1))
bar(f::Foo) = foo(f)

end

""")


text:
using Documenter#, HydrationOfCement

makedocs(modules=[HydrationOfCement],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/thonytruand/HydrationOfCement.git",
    julia  = "1.6.2",
    osname = "windows")


text:
# -----------------------------------------------------------------------
# Total Radiative Forcing From Methane (based on FUND model equations).
# -----------------------------------------------------------------------

@defcomp rf_ch4_total_fund begin

    CH₄_0           = Parameter()             # Initial (pre-industrial) methane concentration (ppb).
    ϕ               = Parameter()             # Scaling term to represent indirect radiative forcing from methane (25% increase from tropospheric ozone, 15% increase from stratospheric water vapor).
    CH₄             = Parameter(index=[time]) # Atmospheric methane concetration for current period (ppb).
    rf_ch4_direct   = Parameter(index=[time]) # Direct forcing from atmospheric methane concentrations (Wm⁻²).

    rf_ch4_indirect = Variable(index=[time])  # Indirect radiative forcing from atmospheric methane concentrations (Wm⁻²).
    rf_ch4_total    = Variable(index=[time])  # Total radiative forcing (direct and indirect effects) from methane (Wm⁻²).

    function run_timestep(p, v, d, t)

        if is_first(t)
            # Set initial radiative forcing to 0.
            v.rf_ch4_indirect[t] = 0.0
            v.rf_ch4_total[t] = 0.0
        else
            # Indirect CH₄ forcing from FUND (following FUND equations, this is 40% of direct forcing without accounting for CH₄-N₂O absorption overlap).
            v.rf_ch4_indirect[t] = 0.036 * p.ϕ * (sqrt(p.CH₄[t]) - sqrt(p.CH₄_0))

            # Calculate total CH₄ radiative forcing.
            v.rf_ch4_total[t] = p.rf_ch4_direct[t] + v.rf_ch4_indirect[t]
        end
    end
end