text:
using Test, Plots, Parameters, LaTeXStrings
# transform to cartesian => polar coordinates for angular integrals
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2858241/

Nx = 200
xm = 10.
x = LinRange(-xm,xm,Nx); y = x

# create some test data
psi = @. 10*exp(-.1*(x^2+y'^2)) + 15*exp(-(x^2+(y-5)'^2)) # + randn()+im*randn()
heatmap(x,y,abs2.(psi),transpose=true)
xlabel!(L"x");ylabel!(L"y")

abstract type Basis end
struct Oscillator <: Basis
    n::Int64
    ω::Float64
end

Oscillator(n::Int64) = Oscillator(n,1.0)
index(b::Oscillator) = 1:b.n
spectrum(b::Oscillator) = (0:b.n-1)*b.ω
qnumbers(b::Oscillator) = 0:b.n-1

function hermite(x,n)
    @assert n >= 1
    h1 = exp(-x^2 / 2) * (1 / π)^(1 / 4)
    n == 1 && (return h1)
    h2 = sqrt(2) * x * h1
    n == 2 && (return h2)
    h3 = 0.0
    for j ∈ 3:n
        h3 = sqrt(2/(j - 1)) * x * h2 - sqrt((j-2)/(j-1)) * h1
        h1 = h2; h2 = h3
    end
    n > 2 && (return h3)
end
hermite(x,n,ω) = hermite(√ω * x,n) * ω ^ ( 1 / 4 )

function (b::Oscillator)(x::Array{Float64,1})
    @unpack ω,n = b
    T = x * ones(1, n) |> zero
    T[:,1] = @. hermite(x,1,ω)
    n >= 2 && (T[:,2] = @. sqrt(2 * ω) * x * T[:,1])
    for j ∈ index(b)[3:end]
        m = qnumbers(b)[j]
        T[:,j] = @. sqrt(2 / m) * (√ω * x) * T[:,j-1] - sqrt((m-1) / m) * T[:,j-2]
    end
    return T
end

(b::Oscillator)(x::LinRange{Float64}) = b(x |> collect)
(b::Oscillator)(x::Float64,n::Int64) = hermite(x,n,b.ω)

function hermite_polar(x,n,ω)
    T = zeros(size(x)...,n)
    T[:,:,1] = @. hermite(x,1,ω)
    n >= 2 && (T[:,:,2] = @. sqrt(2 * ω) * x * T[:,:,1])
    for j in 3:n
        T[:,:,j] = @. sqrt(2 / (j-1)) * (√ω * x) * T[:,:,j-1] - sqrt((j-2) / (j-1)) * T[:,:,j-2]
    end
    return T
end
hermite_polar(x,n) = hermite_polar(x,n,1.)

# test that we can create the basis, or any one of oscillator modes
N = 30
Nv = 1:N
b1 = Oscillator(N,2.)
@time b1.(x,Nv')
plot(x,hermite.(x,Nv'),legend=false,grid=false)
n = 30
b1 = Oscillator(n)
@time b1(x)
plot(x,b1(x),legend=false,grid=false)

# broadcast over any input parameters
plot(x,hermite.(x,2,[1. .5 .2 .1]),legend=false,grid=false)

# filtering cartesian data
function filterH(psi,x,N)
    H = Oscillator(N)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T
    psiF = zero(psi)
        for m = 1:N, n = 1:(N+1-m)
            @. psiF += H[:,m]*F̂[m,n]*H[:,n]'
        end
    return psiF
end

# filter data by projecting onto n oscillator modes
n = 30
@time psiF = filterH(psi,x,n)
heatmap(x,y,abs2.(psiF),transpose=true)
xlabel!(L"x");ylabel!(L"y")

# coversion to polar coords
function slowpolar(psi,x,N)
    H = Oscillator(N)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T

    Nx = length(x)
    θ = LinRange(0,2*pi,2*Nx)
    r = LinRange(0,last(x),Nx/2 |> Int)'

    psiP = zero(θ*r)
    for m = 1:N, n = 1:(N + 1 - m)
        @. psiP += hermite(r*cos(θ),m) * F̂[m,n] * hermite(r*sin(θ),n)
    end
    return psiP
end

function init_polar(x,N)
    θ = LinRange(0,2*pi,2*Nx)
    r = LinRange(0,last(x),Nx/2 |> Int)'
    hx = hermite_polar((@. r*cos(θ)),N)
    hy = hermite_polar((@. r*sin(θ)),N)
    return hx,hy
end

@time hx,hy = init_polar(x,n)

function polar(psi,x,hx,hy)
    N = size(hx)[3]
    H = Oscillator(N)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T

    Nx = length(x)
    θ = LinRange(0,2*pi,2*Nx)
    r = LinRange(0,last(x),Nx/2 |> Int)'

    psiFP = zero(θ*r)
    for m = 1:N, n = 1:(N + 1 - m)
        Hx = @view hx[:,:,m]
        Hy = @view hy[:,:,n]
        @. psiFP += Hx * F̂[m,n] * Hy
    end
    return psiFP
end

# convert to polar coords
r = LinRange(0,last(x),Nx/2 |> Int)'
θ = LinRange(0,2*pi,2*Nx)

@time psiP = slowpolar(psi,x,30)
@time psifp = polar(psi,x,hx,hy)

# plot in polar coordinates
heatmap(θ,r',abs2.(psiP),transpose=true)
xlabel!(L"\theta");ylabel!(L"r")

heatmap(θ,r',abs2.(psifp),transpose=true)
xlabel!(L"\theta");ylabel!(L"r")

dθ = θ[2]-θ[1]
dr = r[2]-r[1]
dx = x[2]-x[1]
psiPth = sum(psiP,dims=1)*dθ
plot(r',psiPth'.*r')
xlabel!(L"r")

# check total integral preserved (< 1%)
sum(psiPth.*r)*dr
sum(psiF)*dx^2

# ===================================
# ==== Now a more realistic example for trapped BEC
using FourierGPE, VortexDistributions
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

# ==== Initialize simulation
L = (25.0,25.0)
N = (256,256)
sim = Sim(L,N)
@unpack_Sim sim

# ==== set simulation parameters
μ = 25.0

# Time dependent potential function (here trivial t dep)
import FourierGPE.V
V(x,y,t) = 0.5*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

# ==== make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim

# ==== evolve
@time sol = runsim(sim)

# speed regression after here?
# ==== ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

# ==== set simulation parameters
γ = 0.0
t = LinRange(ti,tf,Nt)
ϕi = kspace(ψg,sim)
# reltol = 1e-7
# alg = DP5()

# TF parameters
R(w) = sqrt(2*μ/w^2)
R(1)
rv = 3.
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
ξ0 = healinglength.(0.,0.,μ,g)
ξ = healinglength(rv,0.,μ,g)

# vcore = Exact(VortexDistributions.ψi,ξ)
vcore = Exact(ξ)
pv = PointVortex(rv,0.,1)
vi = ScalarVortex(vcore,pv)
psi = Torus(copy(ψg),x,y)
vortex!(psi,vi)
showpsi(x,y,psi.ψ)
ψi .= psi.ψ
ϕi = kspace(ψi,sim)

# ===== compare with Fetter JLTP 2010
ξ = 1/sqrt(μ)
Rtf = R(1)
Ωm = 3*log(Rtf/ξ/sqrt(2))/2/Rtf^2
Ωv = Ωm/(1-rv^2/Rtf^2)

# or a precession period
Tv = 2*pi/Ωv

ti = 0.; tf = Tv
t = LinRange(ti,tf,Nt)

@pack_Sim! sim

# ==== evolve
solv = runsim(sim)

# ===================================
ϕf = solv[end]
ψf = xspace(ϕf,sim)
showpsi(x,y,ψf)

anim = @animate for i in eachindex(t)
    ψ = xspace(solv[i],sim)
    showpsi(x,y,ψ)
end

gif(anim,"./examples/vortexprecession.gif",fps=25)


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

abstract type Instance end

ext: jl


text:
"""
    Transit_Struct
Structure to hold arrays and other quantities for computing transit.
# Members
- `r::T`: radius ratio
- `b::T`: impact parameter
- `jmax::Int64`: maximum number of terms in series expansions of ``I_v`` and ``J_v``

"""

mutable struct Transit_Struct{T<:Real}
	tt0::Array{T,1}
	ttv::Array{T,1}
	jmax::Int64
end

n = Transit_Struct()

max_stars_repo_name: bmlindor/TTVFaster

max_issues_repo_name: bmlindor/TTVFaster

max_forks_repo_name: bmlindor/TTVFaster


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GibbsSeaWater"
version = v"3.5.0"

# Collection of sources required to build GibbsSeaWater
# note 3.05.0-4 is too old to be cross-compiled on FreeBSD
sources = [
    "https://github.com/TEOS-10/GSW-C/archive/d392e91cb63341f543ed1609c4eff613055ab3cb.zip" =>
    "454c46ec00468aeba048fa4f3cee095f927925ab699b1d1b3a25a838cef2d22c",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/GSW-C-*
cc $CFLAGS -fPIC -c -O3 -Wall gsw_oceanographic_toolbox.c gsw_saar.c
cc $LDFLAGS -fPIC -shared -o libgswteos.$dlext gsw_oceanographic_toolbox.o gsw_saar.o -lm
mkdir -p ${libdir}
cp libgswteos.$dlext ${libdir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libgswteos", :libgswteos)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module TestDeprecated

using Test, DataFrames, Random, Logging
import DataFrames: identifier

const ≅ = isequal

old_logger = global_logger(NullLogger())

# old sort(df; cols=...) syntax
df = DataFrame(a=[1, 3, 2], b=[6, 5, 4])
@test sort(df; cols=[:a, :b]) == DataFrame(a=[1, 2, 3], b=[6, 4, 5])
sort!(df; cols=[:a, :b])
@test df == DataFrame(a=[1, 2, 3], b=[6, 4, 5])

df = DataFrame(a=[1, 3, 2], b=[6, 5, 4])
@test sort(df; cols=[:b, :a]) == DataFrame(a=[2, 3, 1], b=[4, 5, 6])
sort!(df; cols=[:b, :a])
@test df == DataFrame(a=[2, 3, 1], b=[4, 5, 6])

@test first.(collect(pairs(DataFrameRow(df, 1, :)))) == [:a, :b]
@test last.(collect(pairs(DataFrameRow(df, 1, :)))) == [df[1, 1], df[1, 2]]

# deprecated combine

df = DataFrame(a=[1, 1, 2, 2, 2], b=1:5)
gd = groupby(df, :a)
@test combine(gd) == combine(identity, gd)

@testset "categorical constructor" begin
    df = DataFrame([Int, String], [:a, :b], [false, true], 3)
    @test !(df[:a] isa CategoricalVector)
    @test df[:b] isa CategoricalVector
    @test_throws DimensionMismatch DataFrame([Int, String], [:a, :b], [true], 3)
end

@testset "DataFrame constructors" begin
    df = DataFrame(Union{Int, Missing}, 10, 3)
    @test size(df, 1) == 10
    @test size(df, 2) == 3
    @test typeof(df[1]) == Vector{Union{Int, Missing}}
    @test typeof(df[2]) == Vector{Union{Int, Missing}}
    @test typeof(df[3]) == Vector{Union{Int, Missing}}
    @test all(ismissing, df[1])
    @test all(ismissing, df[2])
    @test all(ismissing, df[3])
    @test typeof(df[:, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 2]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 3]) == Vector{Union{Int, Missing}}
    @test all(ismissing, df[:, 1])
    @test all(ismissing, df[:, 2])
    @test all(ismissing, df[:, 3])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[1]) == Vector{Union{Int, Missing}}
    @test typeof(df[2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[3]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[1])
    @test all(ismissing, df[2])
    @test all(ismissing, df[3])
    @test typeof(df[:, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[:, 3]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[:, 1])
    @test all(ismissing, df[:, 2])
    @test all(ismissing, df[:, 3])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                [:A, :B, :C], [false, false, true], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[1]) == Vector{Union{Int, Missing}}
    @test typeof(df[2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[3]) <: CategoricalVector{Union{String, Missing}}
    @test all(ismissing, df[1])
    @test all(ismissing, df[2])
    @test all(ismissing, df[3])
    @test typeof(df[:, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[:, 3]) <: CategoricalVector{Union{String, Missing}}
    @test all(ismissing, df[:, 1])
    @test all(ismissing, df[:, 2])
    @test all(ismissing, df[:, 3])
end

df = DataFrame(Union{Int, Missing}, 2, 2)
@test size(df) == (2, 2)
@test eltype.(eachcol(df)) == [Union{Int, Missing}, Union{Int, Missing}]

@test df ≅ DataFrame([Union{Int, Missing}, Union{Float64, Missing}], 2)

@testset "colwise" begin
    Random.seed!(1)
    df = DataFrame(a = repeat(Union{Int, Missing}[1, 3, 2, 4], outer=[2]),
                   b = repeat(Union{Int, Missing}[2, 1], outer=[4]),
                   c = Vector{Union{Float64, Missing}}(randn(8)))

    missingfree = DataFrame([collect(1:10)], [:x1])

    @testset "::Function, ::AbstractDataFrame" begin
        cw = colwise(sum, df)
        answer = [20, 12, -0.4283098098931877]
        @test isa(cw, Vector{Real})
        @test size(cw) == (ncol(df),)
        @test cw == answer

        cw = colwise(sum, missingfree)
        answer = [55]
        @test isa(cw, Array{Int, 1})
        @test size(cw) == (ncol(missingfree),)
        @test cw == answer
    end

    @testset "::Function, ::GroupedDataFrame" begin
        gd = groupby(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise(length, gd) == [[2,2], [2,2]]
    end

    @testset "::Vector, ::AbstractDataFrame" begin
        cw = colwise([sum], df)
        answer = [20 12 -0.4283098098931877]
        @test isa(cw, Array{Real, 2})
        @test size(cw) == (length([sum]),ncol(df))
        @test cw == answer

        cw = colwise([sum, minimum], missingfree)
        answer = reshape([55, 1], (2,1))
        @test isa(cw, Array{Int, 2})
        @test size(cw) == (length([sum, minimum]), ncol(missingfree))
        @test cw == answer

        cw = colwise([Vector{Union{Int, Missing}}], missingfree)
        answer = reshape([Vector{Union{Int, Missing}}(1:10)], (1,1))
        @test isa(cw, Array{Vector{Union{Int, Missing}},2})
        @test size(cw) == (1, ncol(missingfree))
        @test cw == answer

        @test_throws MethodError colwise(["Bob", :Susie], DataFrame(A = 1:10, B = 11:20))
    end

    @testset "::Vector, ::GroupedDataFrame" begin
        gd = groupby(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise([length], gd) == [[2 2], [2 2]]
    end

    @testset "::Tuple, ::AbstractDataFrame" begin
        cw = colwise((sum, length), df)
        answer = Any[20 12 -0.4283098098931877; 8 8 8]
        @test isa(cw, Array{Real, 2})
        @test size(cw) == (length((sum, length)), ncol(df))
        @test cw == answer

        cw = colwise((sum, length), missingfree)
        answer = reshape([55, 10], (2,1))
        @test isa(cw, Array{Int, 2})
        @test size(cw) == (length((sum, length)), ncol(missingfree))
        @test cw == answer

        cw = colwise((CategoricalArray, Vector{Union{Int, Missing}}), missingfree)
        answer = reshape([CategoricalArray(1:10), Vector{Union{Int, Missing}}(1:10)],
                         (2, ncol(missingfree)))
        @test typeof(cw) == Array{AbstractVector,2}
        @test size(cw) == (2, ncol(missingfree))
        @test cw == answer

        @test_throws MethodError colwise(("Bob", :Susie), DataFrame(A = 1:10, B = 11:20))
    end

    @testset "::Tuple, ::GroupedDataFrame" begin
        gd = groupby(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise((length), gd) == [[2,2],[2,2]]
    end
end

@testset "empty!" begin
    df = DataFrame(a=[1, 2], b=[3.0, 4.0])
    @test !isempty(df)

    dfv = view(df, 1:2, 1:2)

    @test empty!(df) === df
    @test isempty(eachcol(df))
    @test isempty(df)
    @test isempty(DataFrame(a=[], b=[]))
end

@testset "deletecols and deletecols!" begin
    df = DataFrame(a=[1,2], b=[3.0, 4.0])
    @test deletecols(df, :a) == DataFrame(b=[3.0, 4.0])
    deletecols!(df, :a)
    @test df == DataFrame(b=[3.0, 4.0])

    df = DataFrame(a=[1,2], b=[3.0, 4.0])
    @test deletecols(df, Not(2)) == DataFrame(b=[3.0, 4.0])
    deletecols!(df, Not(2))
    @test df == DataFrame(b=[3.0, 4.0])

    df = DataFrame(a=[1,2], b=[3.0, 4.0])
    @test deletecols(df, :a, copycols=false)[1] === df.b
    @test deletecols(df, []) == df
    @test deletecols(df, Not([])) == DataFrame()
end

@testset "haskey" begin
    df = DataFrame(x=1:3)
    @test haskey(df, 1)
    @test !haskey(DataFrame(), 1)
    @test !haskey(df, 2)
    @test !haskey(df, 0)
    @test haskey(df, :x)
    @test !haskey(df, :a)
    @test !haskey(df, "a")
end

@testset "df[col] and df[col] for getindex, view, and setindex" begin
    @testset "getindex DataFrame" begin
        df = DataFrame(a=1:3, b=4:6, c=7:9)
        @test df[1] == [1, 2, 3]
        @test df[1] === eachcol(df)[1]
        @test df[1:2] == DataFrame(a=1:3, b=4:6)
        @test df[r"[ab]"] == DataFrame(a=1:3, b=4:6)
        @test df[Not(Not(r"[ab]"))] == DataFrame(a=1:3, b=4:6)
        @test df[Not(3)] == DataFrame(a=1:3, b=4:6)
        @test eachcol(df, false)[1] === df[1]
        @test eachcol(view(df,1:2), false)[1] == eachcol(df, false)[1]
        @test eachcol(df[1:2], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[r"[ab]"], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[Not(Not(r"[ab]"))], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[Not(r"[c]")], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[1:2], false)[1] !== eachcol(df, false)[1]
        @test df[:] == df
        @test df[r""] == df
        @test df[Not(Not(r""))] == df
        @test df[Not(1:0)] == df
        @test df[:] !== df
        @test df[r""] !== df
        @test df[Not(Not(r""))] !== df
        @test df[Not(1:0)] !== df
        @test eachcol(view(df, :), false)[1] == eachcol(df, false)[1]
        @test eachcol(df[:], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[r""], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[Not(1:0)], false)[1] == eachcol(df, false)[1]
        @test eachcol(df[:], false)[1] !== eachcol(df, false)[1]
        @test eachcol(df[r""], false)[1] !== eachcol(df, false)[1]
        @test eachcol(df[Not(1:0)], false)[1] !== eachcol(df, false)[1]
    end
    @testset "getindex df[col] and df[cols]" begin
        x = [1, 2, 3]
        df = DataFrame(x=x, copycols=false)
        @test df[:x] === x
        @test df[[:x]].x !== x
        @test df[:].x !== x
        @test df[r"x"].x !== x
        @test df[r""].x !== x
        @test df[Not(1:0)].x !== x
    end
    @testset "view DataFrame" begin
        df = DataFrame(a=1:3, b=4:6, c=7:9)
        @test view(df, 1) == [1, 2, 3]
        @test view(df, 1) isa SubArray
        @test view(df, 1:2) isa SubDataFrame
        @test view(df, 1:2) == df[1:2]
        @test view(df, r"[ab]") isa SubDataFrame
        @test view(df, r"[ab]") == df[1:2]
        @test view(df, Not(Not(r"[ab]"))) isa SubDataFrame
        @test view(df, Not(Not(r"[ab]"))) == df[1:2]
        @test view(df, :) isa SubDataFrame
        @test view(df, :) == df
        @test parent(view(df, :)) === df
        @test view(df, r"") isa SubDataFrame
        @test view(df, r"") == df
        @test parent(view(df, r"")) === df
        @test view(df, Not(1:0)) isa SubDataFrame
        @test view(df, Not(1:0)) == df
        @test parent(view(df, Not(1:0))) === df
    end
    @testset "getindex SubDataFrame" begin
        df = DataFrame(x=-1:3, a=0:4, b=3:7, c=6:10, d=9:13)
        sdf = view(df, 2:4, 2:4)
        @test sdf[1] == [1, 2, 3]
        @test sdf[1] isa SubArray
        @test sdf[1:2] == DataFrame(a=1:3, b=4:6)
        @test sdf[1:2] isa SubDataFrame
        @test sdf[r"[ab]"] == DataFrame(a=1:3, b=4:6)
        @test sdf[r"[ab]"] isa SubDataFrame
        @test sdf[Not(Not(r"[ab]"))] == DataFrame(a=1:3, b=4:6)
        @test sdf[Not(Not(r"[ab]"))] isa SubDataFrame
        @test sdf[:] == df[2:4, 2:4]
        @test sdf[:] isa SubDataFrame
        @test sdf[r""] == df[2:4, 2:4]
        @test sdf[r""] isa SubDataFrame
        @test sdf[Not(1:0)] == df[2:4, 2:4]
        @test sdf[Not(1:0)] isa SubDataFrame
        @test parent(sdf[:]) === parent(sdf)
        @test parent(sdf[r""]) === parent(sdf)
        @test parent(sdf[Not([])]) === parent(sdf)
    end
    @testset "view SubDataFrame" begin
        df = DataFrame(x=-1:3, a=0:4, b=3:7, c=6:10, d=9:13)
        sdf = view(df, 2:4, 2:4)
        @test view(sdf, 1) == [1, 2, 3]
        @test view(sdf, 1) isa SubArray
        @test view(sdf, 1:2) isa SubDataFrame
        @test view(sdf, 2:3) == df[2:4, 3:4]
        @test view(sdf, r"[ab]") isa SubDataFrame
        @test view(sdf, r"[ab]") == df[2:4, r"[ab]"]
        @test view(sdf, Not(Not(r"[ab]"))) isa SubDataFrame
        @test view(sdf, Not(Not(r"[ab]"))) == df[2:4, Not(Not(r"[ab]"))]
        @test view(sdf, :) isa SubDataFrame
        @test view(sdf, :) == df[2:4, 2:4]
        @test view(sdf, r"") isa SubDataFrame
        @test view(sdf, r"") == df[2:4, 2:4]
        @test view(sdf, Not(1:0)) isa SubDataFrame
        @test view(sdf, Not(1:0)) == df[2:4, 2:4]
        @test parent(view(sdf, :)) == parent(sdf)
        @test parent(view(sdf, r"")) == parent(sdf)
        @test parent(view(sdf, Not(1:0))) == parent(sdf)
    end

    @testset "old setindex! tests" begin
        df = DataFrame(reshape(1:12, 4, :))
        df[1, :] = df[1:1, :]

        df = DataFrame(reshape(1:12, 4, :))

        # Scalar broadcasting assignment of rows
        df[1:2, :] = 1
        df[[true,false,false,true], :] = 3

        # Vector broadcasting assignment of rows
        df[1:2, :] = [2,3]
        df[[true,false,false,true], :] = [2,3]

        # Broadcasting assignment of columns
        df[:, 1] = 1
        df[:x3] = 2

        # assignment of subtables
        df[1, 1:2] = df[2:2, 2:3]
        df[1:2, 1:2] = df[2:3, 2:3]
        df[[true,false,false,true], 2:3] = df[1:2,1:2]

        # scalar broadcasting assignment of subtables
        df[1:2, 1:2] = 3
        df[[true,false,false,true], 2:3] = 3

        # vector broadcasting assignment of subtables
        df[1:2, 1:2] = [3,2]
        df[[true,false,false,true], 2:3] = [2,3]

        # test of 1-row DataFrame assignment
        df = DataFrame([1 2 3])
        df[1, 2:3] = DataFrame([11 12])
        @test df == DataFrame([1 11 12])

        df = DataFrame([1 2 3])
        df[1, [false, true, true]] = DataFrame([11 12])
        @test df == DataFrame([1 11 12])
    end

    @testset "old test/dataframes.jl tests" begin
        df = DataFrame(a = Union{Int, Missing}[2, 3],
                    b = Union{DataFrame, Missing}[DataFrame(c = 1), DataFrame(d = 2)])
        df[1, :b][:e] = 5

        x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])

        #test_group("DataFrame assignment")
        # Insert single column
        x0 = x[Int[], :]
        @test_throws ArgumentError x0[:d] = [1]
        @test_throws ArgumentError x0[:d] = 1:3

        # Insert single value
        x[:d] = 3
        @test x[:d] == [3, 3, 3]
    end

    @testset "test getindex using df[col] and df[cols] syntax" begin
        x = [1]
        y = [1]
        df = DataFrame(x=x, y=y, copycols=false)
        @test df.x === x
        @test df[:y] === y
        @test df[1] === x
        @test df[1:1][1] == x
        @test df[r"x"][1] == x
        @test df[1:1][1] !== x
        @test df[r"x"][1] !== x
        @test df[1:2][:y] == y
        @test df[1:2][:y] !== y
        @test df[r""][:y] == y
        @test df[r""][:y] !== y
        @test df[:][:x] == x
        @test df[:][:x] !== x
        @test df[[:y,:x]][:x] == x
        @test df[[:y,:x]][:x] !== x
    end

    @testset "setindex! special cases" begin
        df = DataFrame(rand(3,2), [:x3, :x3_1])
        @test_throws ArgumentError df[3] = [1, 2]
        @test_throws ArgumentError df[4] = [1, 2, 3]
        df[3] = [1,2,3]
        df[4] = [1,2,3]
        @test names(df) == [:x3, :x3_1, :x3_2, :x4]
        df = DataFrame()
        @test_throws MethodError df[true] = 1
        @test_throws MethodError df[true] = [1,2,3]
        @test_throws MethodError df[1:2, true] = [1,2]
        @test_throws MethodError df[1, true] = 1
        @test_throws ArgumentError df[1, 100] = 1
        @test_throws BoundsError df[1:2, 100] = [1,2]
    end

    @testset "handling of end in indexing" begin
        z = DataFrame(rand(4,5))
        for x in [z, view(z, 1:4, :)]
            y = deepcopy(x)
            @test x[end] == x[5]
            @test x[end:end] == x[5:5]
            @test x[end, :] == x[4, :]
            x[end] = 1:4
            y[5] = 1:4
            @test x == y
            x[4:end] = DataFrame([11:14, 21:24])
            y[4] = [11:14;]
            y[5] = [21:24;]
            @test x == y
        end
    end

    @testset "aliasing in indexing" begin
        # columns should not alias if scalar broadcasted
        df = DataFrame(A=[0], B=[0])
        df[1:end] = 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0

        df = DataFrame(A=[0], B=[0])
        df[:, 1:end] = 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0

        # columns should not alias if vector assigned
        df = DataFrame(A=[0], B=[0])
        x = [0.0]
        df[1:end] = x
        x[1] = 1.0
        @test df[1, :A] === 0.0
        @test df[1, :B] === 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0.0

        df = DataFrame(A=[0], B=[0])
        x = [0.0]
        df[:, 1:end] = x
        x[1] = 1.0
        @test df[1, :A] === 0.0
        @test df[1, :B] === 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0.0
    end
end

@testset "eltypes" begin
    @test eltypes(DataFrame(x=[1], y=["a"])) == [Int, String]
end

@testset "melt" begin
    mdf = DataFrame(id=[missing,1,2,3], a=1:4, b=1:4)
    @test unstack(melt(mdf, :id), :id, :variable, :value)[1:3,:] == sort(mdf)[1:3,:]
    @test unstack(melt(mdf, :id), :id, :variable, :value)[:, 2:3] == sort(mdf)[:, 2:3]
    @test unstack(melt(mdf, Not(Not(:id))), :id, :variable, :value)[1:3,:] == sort(mdf)[1:3,:]
    @test unstack(melt(mdf, Not(Not(:id))), :id, :variable, :value)[:, 2:3] == sort(mdf)[:, 2:3]

    Random.seed!(1234)
    x = DataFrame(rand(100, 50))
    x[!, :id] = [1:99; missing]
    x[!, :id2] = string.("a", x[!, :id])
    x[!, :s] = [i % 2 == 0 ? randstring() : missing for i in 1:100]
    allowmissing!(x, :x1)
    x[1, :x1] = missing
    y = melt(x, [:id, :id2])
    @test y ≅ melt(x, r"id")
    @test y ≅ melt(x, Not(Not(r"id")))

    d1 = DataFrame(a = Array{Union{Int, Missing}}(repeat([1:3;], inner = [4])),
                b = Array{Union{Int, Missing}}(repeat([1:4;], inner = [3])),
                c = Array{Union{Float64, Missing}}(randn(12)),
                d = Array{Union{Float64, Missing}}(randn(12)),
                e = Array{Union{String, Missing}}(map(string, 'a':'l')))
    d1s = stack(d1, [:a, :b])
    d1m = melt(d1, [:c, :d, :e])
    @test d1m == melt(d1, r"[cde]")
    @test d1s == d1m
    d1m = melt(d1[:, [1,3,4]], :a)
    @test names(d1m) == [:variable, :value, :a]
    d1m_named = melt(d1[:, [1,3,4]], :a, variable_name=:letter, value_name=:someval)
    @test names(d1m_named) == [:letter, :someval, :a]
    dx = melt(d1, [], [:a])
    @test dx == melt(d1, r"xxx", r"a")
    @test size(dx) == (12, 2)
    @test names(dx) == [:variable, :value]
    dx = melt(d1, :a, [])
    @test dx == stack(d1, r"xxx", r"a")
    @test size(dx) == (0, 3)
    @test names(dx) == [:variable, :value, :a]
    d1m = melt(d1, [:c, :d, :e], view=true)
    @test d1m == melt(d1, r"[cde]", view=true)
    d1m = melt(d1[:, [1,3,4]], :a, view=true)
    @test names(d1m) == [:variable, :value, :a]
    d1m_named = melt(d1, [:c, :d, :e], variable_name=:letter, value_name=:someval, view=true)
    @test d1m_named == melt(d1, r"[cde]", variable_name=:letter, value_name=:someval, view=true)
    @test names(d1m_named) == [:letter, :someval, :c, :d, :e]
    df1 = melt(DataFrame(rand(10,10)))
    df1[!, :id] = 1:100
    @test size(unstack(df1, :variable, :value)) == (100, 11)
    @test_throws ArgumentError unstack(melt(DataFrame(rand(3,2))), :variable, :value)
end

global_logger(old_logger)

end # module


text:
include("../../src/EDKit.jl")
include("ExactDiagonalization.jl")
using .EDKit
import .OldED as OED
using LinearAlgebra
using Test
#-------------------------------------------------------------------------------------------------------------------------
const X = [0 1; 1 0]
const Y = [0 -1; 1 0] * 1im
const Z = [1 0; 0 -1]
#-------------------------------------------------------------------------------------------------------------------------
@testset "Random Real 4×4 Matrices" begin
    L = 14
    println("Random Real 4×4 Matrices, L = $L")
    mats = [rand(4, 4) |> Hermitian |> Array for i=1:L]
    inds = [mod.([i-1, i], L) .+ 1 for i=1:L]
    print("New :")
    @time op1 = operator(mats, inds, L) |> Array
    print("Old :")
    @time op2 = OED.operation(mats, inds, L) |> Array
end
#-------------------------------------------------------------------------------------------------------------------------
@testset "Random Complex 4×4 Matrices" begin
    L = 14
    println("Random Complex 4×4 Matrices, L = $L")
    mats = [rand(ComplexF64, 4, 4) |> Hermitian |> Array for i=1:L]
    inds = [mod.([i-1, i], L) .+ 1 for i=1:L]
    print("New :")
    @time op1 = operator(mats, inds, L) |> Array
    print("Old :")
    @time op2 = OED.operation(mats, inds, L) |> Array
end
#-------------------------------------------------------------------------------------------------------------------------
@testset "Sparse Matrix IXYZ" begin
    L = 14
    println("Sparse Matrix IXYZ, L = $L")
    mat = kron(I(2), X, Y, Z)
    print("New :")
    @time op1 = trans_inv_operator(mat, 4, L) |> Array
    print("Old :")
    @time op2 = OED.trans_inv_operation(mat, 1:4, L) |> Array
end
#-------------------------------------------------------------------------------------------------------------------------
@testset "Operator Multiplication" begin
    L = 14
    println("Operator Multiplication, L = $L")
    mat = kron(I(2), X, Y, Z)
    v = rand(2^L)
    opt = trans_inv_operator(mat, 4, L)
    print("Operator : ")
    @time v2 = opt * v
    print("Dense    : ")
    optm = Array(opt)
    @time v3 = optm * v
    @test v2 ≈ v3
end


text:
module TestModelsWithError

using Modia
using Modia.Electric

println("MODELS WITH ERRORS")
println()

Modia.ModiaLogging.resetTestStatus()

@model WrongName begin
  R=Resist()
end
checkSimulation(WrongName, 1)

@model WrongParameterName begin
  R=Resistor(res=100.0)
end
# Pantelides loops!!!
# checkSimulation(WrongParameterName, 1)

@model WrongConnect1 begin
  R=Resistor(R=100.0)
  C=Capacitor()
  @equations begin
    connect(R.n, C.pin)
  end
end
checkSimulation(WrongConnect1, 1)

@model WrongConnect2 begin
  R=Resistor(R=100.0)
  C=Capacitor()
  @equations begin
    connect(R.n)
  end
end
checkSimulation(WrongConnect2, 1)

@model WrongEquation1 begin
  R=Resistor(res=100.0)
  y = Var()
  @equations begin
  R.I = y
  end
end
println("Should complain about I not declared!")
checkSimulation(WrongEquation1, 1)

@model WrongEquation2 begin
  R=Resistor(res=100.0)
  @equations begin
  R.i = 0
  end
end
checkSimulation(WrongEquation2, 1)

@model Base begin
  v = Var(start=0.0)
end

@model Extended1 begin
  @extends MyBasic # Forgot()
end
checkSimulation(Extended1, 1)

@model Extended2 begin
  @extends Base # Forgot ()
end
checkSimulation(Extended2, 1)

println("UndefinedVariable test disabled since it stops script. See model code.")
#=
@model UndefinedVariable begin
  @equations begin
    v = 0.0 # ERROR: LoadError: UndefVarError: v not defined
  end
end
checkSimulation(UndefinedVariable, 1)
=#

println("Extended3 test disabled since it stops script. See model code.")
@model Extended3 begin
  @extends Base()
  @equations begin
#    v = 0.0 # ERROR: LoadError: UndefVarError: v not defined
    this.v = 0.0
  end
end
checkSimulation(Extended3, 1)

println("Does not complain about undefined dummy!")
@model Extended4 begin
  @extends Base()
  @inherits v, dummy
end
checkSimulation(Extended4, 1)

@model Extended5 begin
  @extends Base()
  @inherits v, dummy
  w = Var()
end
checkSimulation(Extended5, 1)

Modia.ModiaLogging.printTestStatus()

end


text:
using VLKeggSDK
using DataFrames

# setup an array of ec numbers -
rn_number_array = [
    "rn:R00299"
    "rn:R00771"
    "rn:R04779"
    "rn:R01068"
]

# initialize -
compound_record_array = Array{KEGGCompound,1}()
for rn_number in rn_number_array

    # get the compound records -
    compound_record = get_compound_records_for_reaction(rn_number) |> check
    if (isnothing(compound_record) == false)
        append!(compound_record_array, compound_record)
    end
end

# update -
# unique!(compound_record_array)

unique_compound_array = Array{KEGGCompound,1}()
for compound_object in compound_record_array
    if (in(compound_object, unique_compound_array) == false)
        push!(unique_compound_array, compound_object)
    end
end


text:
TEST_MSG = "test log message"

@testset "Test configure_logging" begin
    # Verify logging to a file.
    logfile = "testlog.txt"
    isfile(logfile) && rm(logfile)

    orig_logger = global_logger()
    try
        logger = configure_logging(
            console_level = Logging.Warn,
            file_level = Logging.Debug,
            filename = logfile,
        )
        with_logger(logger) do
            @info TEST_MSG
        end

        close(logger)
        @test isfile(logfile)
        lines = readlines(logfile)
        @test length(lines) == 2  # two lines per message
        @test occursin(TEST_MSG, lines[1])
    finally
        rm(logfile)
        global_logger(orig_logger)
    end

    # Verify logging with no file.
    try
        logger = configure_logging(
            console_level = Logging.Warn,
            file_level = Logging.Info,
            filename = nothing,
        )
        with_logger(logger) do
            @info TEST_MSG
        end

        close(logger)
        @test !isfile(logfile)
    finally
        global_logger(orig_logger)
    end
end


text:
abstract type AbstractStepper end

export StrangSplitting

struct StrangSplitting <: AbstractStepper end

max_stars_repo_name: JuliaVlasov/VlasovSolvers

max_issues_repo_name: JuliaVlasov/VlasovSolvers

max_forks_repo_name: JuliaVlasov/VlasovSolvers


text:
@testset "SOCtoPSD" begin
    bridged_mock = MOIB.SOCtoPSD{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[√2/2, -1/2, √2/4, -1/2, √2/4, √2/4]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)                            => [[-√2]])
    MOIT.soc1vtest(bridged_mock, config)
    MOIT.soc1ftest(bridged_mock, config)
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone}()))
    test_delete_bridge(bridged_mock, ci, 3, ((MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0),))
end

@testset "RSOCtoPSD" begin
    bridged_mock = MOIB.RSOCtoPSD{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2, 0.5, 1.0],
                          (MOI.SingleVariable,                MOI.EqualTo{Float64})       => [-√2, -1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[√2, -1/2, √2/8, -1/2, √2/8, √2/8]])
    MOIT.rotatedsoc1vtest(bridged_mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[√2, -1/2, √2/8, -1/2, √2/8, √2/8]])
    MOIT.rotatedsoc1ftest(bridged_mock, config)
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}()))
    test_delete_bridge(bridged_mock, ci, 2, ((MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0),))
end


text:
macro queue(expr...)
	esc(genqueue(collect(expr)))
end

function genqueue(e::Vector)
	q  = Expr(:new, :List , nothing  )
	num = length(e)
	for i in num:-1:1
		q = Expr(:new, :List , e[i] , q  )
	end
	return q
end

max_stars_repo_name: frankill/ftest

max_issues_repo_name: frankill/ftest

max_forks_repo_name: frankill/ftest


text:
export
    get_control,
    TVLQR,
    LQR,
    MLQR,
    simulate,
    SE3Tracking,
    HFCA,
    test_ICs

import TrajectoryOptimization.Dynamics: trim_controls, build_state
import TrajectoryOptimization: state_diff, state_diff_jacobian
abstract type AbstractController end
abstract type LQRController <: AbstractController end
abstract type TimeVaryingController <: AbstractController end
abstract type TrackingController <: TimeVaryingController end

function get_k(cntrl::TimeVaryingController, t)
    times = get_times(cntrl)
    searchsortedlast(times,t)
end

""" TVLQR
"""
struct TVLQR{L,TK} <: TimeVaryingController
    model::L
    Z::Traj
    K::Vector{TK}
    t::Vector{Float64}
end

""" Assumes the states of the trajectory Z match the states of the model
"""
function TVLQR(model::AbstractModel, Q, R, Z::Traj)
    A,B = linearize(model, Z)
    K = tvlqr(A,B,Q,R)
    t = [z.t for z in Z]
    TVLQR(model,Z,K,t)
end

function TVLQR(model::RigidBody, Q, R, X::Vector{<:RBState}, U::Vector{<:AbstractVector}, dt)
    # Convert the trajectory to the orientation representation of the model
    Z = Traj(model, X, U, dt)

    # Call as normal
    TVLQR(model, Q, R, Z::Traj)
end

@inline get_times(cntrl::TVLQR) = cntrl.t

function get_control(cntrl::TVLQR, x, t)
    k = get_k(cntrl, t)
    zref = cntrl.Z[k]
    dx = state_diff(cntrl.model, x, state(zref))
    return cntrl.K[k]*dx + control(zref)
end

function tvlqr(A, B, Q, R)
    n,m = size(B[1])
    K = [@SMatrix zeros(m,n) for k = 1:length(A)]
    tvlqr!(K, A, B, Q, R)
end

function tvlqr!(K, A, B, Q, R)
    N = length(A)

    # Solve infinite-horizon at goal state
    # Qf = dare(A[N], B[N], Q, R)
    Qf = Q

    P_ = similar_type(A[N])
    P = copy(Qf)
    for k = N-1:-1:1
        P,K[k] = riccati(P,A[k],B[k],Q,R)
    end
    return K
end

function riccati(P,A,B,Q,R)
    P_ = Q + A'P*A - A'P*B*((R+B'P*B)\(B'P*A))
    K = -(R+B'P*B)\(B'P*A)
    return P_,K
end

""" Discrete LQR
"""
struct LQR{T,N,M,TK} <: LQRController
    K::TK
    xref::SVector{N,T}
    uref::SVector{M,T}
end

function LQR(model::AbstractModel, dt::Real, Q, R,
        xeq=zeros(model)[1], ueq=trim_controls(model))
    # Linearize the model
    A,B = linearize(model, xeq, ueq, dt)

    # Calculate the optimal control gain
    K = calc_LQR_gain(A,B,Q,R)
    LQR(K, xeq, ueq)
end

function get_control(cntrl::LQR, x, t)
    dx = x - cntrl.xref
    return cntrl.K*dx + cntrl.uref
end

""" Multiplicative LQR
"""
struct MLQR{T,N,M,TK} <: LQRController
    model::AbstractModel
    K::TK
    xref::SVector{N,T}
    uref::SVector{M,T}
end

function MLQR(model::AbstractModel, dt::Real, Q::AbstractMatrix, R::AbstractMatrix,
        xeq=zeros(model)[1], ueq=trim_controls(model))
    # Linearize the model, accounting for attitude state
    A,B = linearize(model, xeq, ueq, dt)
    G1 = state_diff_jacobian(model, xeq)
    G2 = state_diff_jacobian(model, discrete_dynamics(RK3, model, xeq, ueq, 0.0, dt))
    A = G2'A*G1
    B = G2'B

    # Calculate the optimal control gain
    @assert size(Q) == size(A)
    K = calc_LQR_gain(A,B,Q,R)
    MLQR(model,K,xeq,ueq)
end

function get_control(cntrl::MLQR, x, t)
    dx = state_diff(cntrl.model, x, cntrl.xref)
    return cntrl.K*dx + cntrl.uref
end

""" SE(3) Quadrotor controller
from [Lee 2010]
"""
struct SE3Tracking{T,N} <: TrackingController
    model::AbstractModel
    kx::T
    kv::T
    kR::T
    kΩ::T
    Xref::Vector{SVector{N,T}}
    Xdref::Vector{SVector{N,T}}
    bref::Vector{SVector{3,T}}
    t::Vector{T}
end

function SE3Tracking(model::AbstractModel, Xref, Xdref, bref, t;
        # kx=59.08, kv=24.3, kR=8.81, kO=1.54)
        kx=2.71, kv=1.01, kR=2.26, kO=0.1)
    SE3Tracking(model, kx, kv, kR, kO, Xref, Xdref, bref, t)
end

@inline get_times(cntrl::SE3Tracking) = cntrl.t

function get_control(c::SE3Tracking, x, t)
    # Get model params
    g = c.model.gravity[3]
    e3 = @SVector [0,0,1.]
    mass = c.model.mass
    J = c.model.J

    # Get time step
    k = get_k(c, t)
    b1d = c.bref[k]
    xd = c.Xref[k]

    # Parse the state
    r,q,v,Ω = Dynamics.parse_state(c.model, x)
    rd,qd,vd,Ωd = Dynamics.parse_state(c.model, xd)
    rdd,qdd,vdd,Ωdd = Dynamics.parse_state(c.model, c.Xdref[k])
    R = rotmat(q)

    # Calculate the linear errors
    ex = r-rd
    ev = v-vd
    xdd = rdd
    Ωdotd = Ωdd

    # Calculate desired attitude
    a = -c.kx*ex - c.kv*ev - mass*g*e3 + mass*xdd
    b3 = normalize(a)
    b2 = normalize(b3 × b1d)
    b1 = b2 × b3
    Rd = @SMatrix [b1[1] b2[1] b3[1];
                   b1[2] b2[2] b3[2];
                   b1[3] b2[3] b3[3]]

    # Calculate the attitude errors
    eR = 0.5*vee(Rd'R - R'Rd)
    eΩ = Ω - R'Rd*Ωd

    # Desired Forces and moments
    f = -a'R*e3
    M = -c.kR*eR - c.kΩ*eΩ + Ω × (J*Ω) - J*(skew(Ω)*R'Rd*Ωd - R'Rd*Ωdotd)

    # Convert to motor thrusts
    C = Dynamics.forceMatrix(c.model)
    u = C\(@SVector [f, M[1], M[2], M[3]])
end

""" Hopf Fibration control
"""
struct HFCA{N,T} <: TrackingController
    model::AbstractModel
    kx::T
    kv::T
    kR::Diagonal{T,SVector{3,T}}
    kΩ::Diagonal{T,SVector{3,T}}
    Xref::Vector{SVector{N,T}}
    Xdref::Vector{SVector{N,T}}
    bref::Vector{SVector{3,T}}
    t::Vector{T}
end

function HFCA(model::AbstractModel, Xref, Xdref, bref, t;
        # kx=59.08, kv=24.3, kR=8.81, kO=1.54)
        kx=11.7, kv=4.0,
        kR=3.0*Diagonal(@SVector ones(3)),
        kO=0.7*Diagonal(@SVector ones(3)))
    HFCA(model, kx, kv, kR, kO, Xref, Xdref, bref, t)
end

get_times(c::HFCA) = c.t

function get_control(c::HFCA, x, t)
    # Get model params
    g = c.model.gravity[3]
    e3 = @SVector [0,0,1.]
    mass = c.model.mass
    J = c.model.J

    # Get time step
    k = get_k(c, t)
    b1d = c.bref[k]
    xd = c.Xref[k]
    xdd = c.Xdref[k]

    # Parse the state
    r,q,v,Ω = Dynamics.parse_state(c.model, x)
    rd,qd,vd,Ωd = Dynamics.parse_state(c.model, xd)
    rdd,qdd,vdd,Ωdd = Dynamics.parse_state(c.model, xdd)
    R = rotmat(q)

    # Calculate the linear errors
    ex = r-rd
    ev = v-vd
    xdd = rdd
    Ωdotd = Ωdd

    # Get the desired heading
    ψ0 = atan(c.bref[1][2], c.bref[1][1])
    ψ = atan(b1d[2], b1d[1]) - ψ0

    # Calculate desired angular velocities and quaternion
    ζ = -c.kx*ex - c.kv*ev + mass*rdd - mass*g*e3  # eq (26)
    ζdot = -c.kx*ev + vdd  # is this right?
    abc = normalize(ζ)  # eq (19)
    abc_dot = (ζ'ζ*I - ζ*ζ')*ζdot/(norm(ζ)^3)  # eq (20)
    ψdot = -(abc[2]*abc_dot[1] - abc[1]*abc_dot[2])/(1+abc[3])
    ω = (@SVector [
        sin(ψ)*abc_dot[1] - cos(ψ)*abc_dot[2] - (abc[1]*sin(ψ) - abc[2]*cos(ψ)),
        cos(ψ)*abc_dot[1] + sin(ψ)*abc_dot[2] - (abc[1]*cos(ψ) + abc[2]*sin(ψ)),
        0
    ]) * (abc_dot[3]/(abc[3] + 1))
    q_abc = @SVector [1+abc[3], -abc[2], abc[1], 0]
    q_ψ = @SVector [cos(ψ/2), 0, 0, sin(ψ/2)]
    qd = Lmult(q_abc)*q_ψ
    Rd = rotmat(UnitQuaternion(qd))

    # Calculate angular errors
    eΩ = Ω - R'Rd*Ωd
    eR = 0.5*vee(Rd'R - R'Rd)

    # Calculate controls
    f = -ζ'R*e3
    M = -c.kR*eR - c.kΩ*eΩ

    # Convert to motor thrusts
    C = Dynamics.forceMatrix(c.model)
    u = C\(@SVector [f, M[1], M[2], M[3]])
end

############################################################################################
############################################################################################

function linearize(model::AbstractModel, xeq, ueq, dt)
    # Linearize the system about the given point
    n,m = size(model)
    z = KnotPoint(xeq, ueq, dt)
    ∇f = zeros(n,n+m+1)
    discrete_jacobian!(RK3, ∇f, model, z)
    ix,iu = z._x, z._u
    A = ∇f[ix,ix]
    B = ∇f[ix,iu]
    return A,B
end

# TODO: redo without static matrices
function linearize(model::AbstractModel, Z::Traj)
    N = length(Z)
    D = [TO.SizedDynamicsExpansion(model) for k = 1:N]
    for k = 1:N
        _linearize!(D[k], model, Z[k])
    end
    A = [SMatrix(d.A) for d in D]
    B = [SMatrix(d.B) for d in D]
    return A,B
end

function _linearize!(D::TO.SizedDynamicsExpansion, model::AbstractModel, z::KnotPoint)
    discrete_jacobian!(RK3, D.∇f, model, z)
	D.tmpA .= D.A_  # avoids allocations later
	D.tmpB .= D.B_
    x2 = discrete_dynamics(RK3, model, z)
    G1 = TO.state_diff_jacobian(model, state(z))
    G2 = TO.state_diff_jacobian(model, x2)
    TO.error_expansion!(D,G1,G2)
    return nothing
end

function calc_LQR_gain(A::AbstractMatrix, B::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix;
        max_iters=100, tol=1e-4)
    P_ = similar_type(A)
    P = P_(A)
    for k = 1:max_iters
        P_ = Q + A'P*A - A'P*B*((R + B'P*B)\(B'P*A))
        err = norm(P-P_)
        P = copy(P_)
        if err < tol
            break
        end
    end
    K = -(R + B'P*B)\(B'P*A)
end

function dare(A,B,Q,R,
        max_iters=100, tol=1e-4)
    P_ = similar_type(A)
    P = P_(A)
    for k = 1:max_iters
        P_ = Q + A'P*A - A'P*B*((R + B'P*B)\(B'P*A))
        err = norm(P-P_)
        P = copy(P_)
        if err < tol
            break
        end
    end
    return P
end

simulate(model::RigidBody, cntrl, x0::RBState, tf; kwargs...) =
    simulate(model, cntrl, build_state(model, x0), tf; kwargs...)

function simulate(model::AbstractModel, cntrl, x0, tf;
        dt=1e-4, w=1e-2)
    N = Int(round(tf/dt)) + 1
    dt = tf/(N-1)
    x = copy(x0)
    t = 0.0

    n,m = size(model)

    # Allocate storage for the state trajectory
    X = [@SVector zeros(n) for k = 1:N]

    for k = 1:N
        # Get control from the controller
        u = get_control(cntrl, x, t)

        # Add disturbances
        u += randn(m)*w

        # Simulate the system forward
        x = discrete_dynamics(RK3, model, x, u, t, dt)
        t += dt  # advance time

        # Store info
        X[k] = x
    end
    return X
end

function test_ICs(model::AbstractModel, cntrl::AbstractController, ICs::Vector{<:RBState};
        tf=10., dt=1e-4, xref=(@SVector zeros(N))) where N
    L = length(ICs)
    Xf = deepcopy(ICs)
    data = Dict{Symbol,Vector}(:Xf=>deepcopy(ICs),
        :max_err=>zeros(L), :avg_err=>zeros(L), :term_err=>zeros(L))
    for i = 1:L
        x0 = Dynamics.build_state(model, ICs[i])
        X = simulate(model, cntrl, x0, tf, dt=dt, w=0.0)
        err = map(X) do x
            dx = Dynamics.state_diff(model, x, x0)
            norm(dx)
        end
        data[:max_err][i] = maximum(err)
        data[:avg_err][i] = mean(err)
        data[:term_err][i] = err[end]
        data[:Xf][i] = RBState(model, X[end])
    end
    data
end


text:
module DominantSetTests

using Test
using WeaklyHard

const n_tests   = 10_000
const k_max     = 200
@testset "2 Random Constraints" begin
    constraint_names = ["RowMissConstraint", "RowHitConstraint", "AnyMissConstraint", "AnyHitConstraint"]

    for _ in 1:n_tests

        k1 = rand(1:k_max)
        k2 = rand(1:k_max)
        x1 = rand(0:k1)
        x2 = rand(0:k2)

        l1_type = rand(constraint_names)
        l2_type = rand(constraint_names)
        if l1_type == "RowHitConstraint"
            x1 = floor(Int, x1/2)
        end
        if l2_type == "RowHitConstraint"
            x2 = floor(Int, x2/2)
        end
        lambda_1 = getfield(WeaklyHard, Symbol(l1_type))(x1, k1)
        lambda_2 = getfield(WeaklyHard, Symbol(l2_type))(x2, k2)

        Lambda = Set{Constraint}([lambda_1, lambda_2])
        Lambda_star = dominant_set(Lambda)

        if is_equivalent(lambda_1, lambda_2)
            @test length(Lambda_star) == 1
        elseif is_dominant(lambda_1, lambda_2)
            @test length(Lambda_star) == 1 && (lambda_1 in Lambda_star)
        elseif !is_dominant(lambda_2, lambda_1)
            @test length(Lambda_star) == 2 && (lambda_1 in Lambda_star) && (lambda_2 in Lambda_star)
        else
            @test length(Lambda_star) == 1 && (lambda_2 in Lambda_star)
        end #if
    end # for
end # testset

@testset "Arbitrary Random Constraints" begin
    constraint_names = ["RowMissConstraint", "RowHitConstraint", "AnyMissConstraint", "AnyHitConstraint"]
    n_constraints   = 30

    for _ in 1:n_tests
        Lambda = Set{Constraint}()
        n = rand(1:n_constraints)
        for _ in 1:n
            k = rand(1:k_max)
            x = rand(0:k)
            lambda_type = rand(constraint_names)
            if lambda_type == "RowHitConstraint"
                x = floor(Int, x/2)
            end
            push!(Lambda, getfield(WeaklyHard, Symbol(lambda_type))(x, k))
        end # for

        Lambda_star = dominant_set(Lambda)
        Lambda_comp = setdiff(Lambda, Lambda_star)

        test_res = true
        for lambda_c in Lambda_comp
            exist_dominant = false
            for lambda_s in Lambda_star
                if is_dominant(lambda_s, lambda_c)
                    exist_dominant = true
                    break
                end
            end
            if !exist_dominant
                @show Lambda
                test_res = false
                break
            end # if
        end # for
        @test test_res
    end # for
end # testset

end # module


text:
import Pkg
Pkg.activate("..")

using CrimsonSkyline
using Distributions: Normal, Gamma, LogNormal
using Logging
using StatsBase: mean, std
using Random: seed!

seed!(2021)

function random_walk(t :: Trace, data :: Array{Float64, 2})
    loc = sample(t, :loc, Normal())
    scale = sample(t, :scale, Gamma(3.0, 5.0))
    shape = size(data)
    for n in 1:shape[1]
        x_n = sample(t, (:ic, n), Normal())
        for t_ix in 1:shape[2]
            x_n = observe(t, (:data, n, t_ix), Normal(loc + x_n, scale), data[n, t_ix])
        end
    end
end

function global_trend(t :: Trace, data :: Array{Float64, 2})
    α = sample(t, :α, Normal())
    β = sample(t, :β, Normal())
    σ = sample(t, :σ, LogNormal())
    shape = size(data)
    for n in 1:shape[1]
        for t_ix in 1:shape[2]
            observe(t, (:data, n, t_ix), Normal(α + β * t_ix, σ), data[n, t_ix])
        end
    end
end

function main()
    n = 10
    t = 20
    t_pred = 10

    @info "Random walk inference"
    loc = 1.0
    scale = 1.5
    noise = randn(n, t + t_pred)
    data = cumsum(loc .+ scale .* noise, dims=2)
    test_data = data[:, t:end]
    data = data[:, 1:t]
    @info "True loc = $loc, true scale = $scale"
    @time results = mh(random_walk; params=(data,), burn=1000, thin=50, num_iterations=50000)
    @info "Estimated loc = ($(mean(results, :loc)) ± $(2.0 * std(results, :loc)))"
    @info "Estimated scale = ($(mean(results, :scale)) ± $(2.0 * std(results, :scale)))"

    @info begin
        """
        Time series model comparison
        Comparing (prior) evidence for random walk model vs global trend
        """
    end
    nsamples = 5000
    rw_results = likelihood_weighting(random_walk, data; nsamples = nsamples)
    gt_results = likelihood_weighting(global_trend, data; nsamples = nsamples)
    log_evidence_rw = log_evidence(rw_results)
    log_evidence_gt = log_evidence(gt_results)
    @info begin
        """
        log p(x | random walk) ≈ $log_evidence_rw
        log p(x | global trend) ≈ $log_evidence_gt
        """
    end
    @info "Comparing model quality using AIC"
    aic_rw = aic(rw_results)
    aic_gt = aic(gt_results)
    @info "AIC(random walk) ≈ $aic_rw, AIC(global trend) ≈ $aic_gt"

end

main()


text:
# sort an array - aplying a function to the sort process
a1 = [1, 2, 0, -2, 7, 3, -11, 12, -43]
res = sort(a1, by=abs) # sort on the basis of the elements absolute value ...
println("$a1  / $res")

# is this array sorted
a2 = [1, 3, 6, 3]
println(issorted(a1))

a3 = [(1, 2), (2, 4), (3, 1), (5, 2)]
println(issorted(a3, by=x -> x[1])) # is this array sorted on the basis of the lambda function provided - get the first element of each tuple

# sorting custom types (structs)
struct Rect
    length::Int
    width::Int
end

a4 = [Rect(12, 22), Rect(10, 18), Rect(7, 12), Rect(9, 30)]

# create a custom sort function
function comparearea(r1::Rect, r2::Rect)
    area1 = r1.length * r1.width
    area2 = r2.length * r2.width

    return area1 < area2
end

sort!(a4, lt=comparearea) # sort 'in place', provide a 'less than' function (lt)
println(a4)

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
module ColoredLLCodes

using InteractiveUtils
import InteractiveUtils: _dump_function, code_llvm, code_native

const IO_ = Union{Base.AbstractPipe, Base.LibuvStream} # avoid overwriting

if VERSION > v"1.5" && :highlighting in names(InteractiveUtils, all=true)
    global highlighting = InteractiveUtils.highlighting
else
    global highlighting = Dict{Symbol, Bool}(
        :llvm => true,
        :native => true,
    )
end

if VERSION > v"1.5" && :llstyle in names(InteractiveUtils, all=true)
    global llstyle = InteractiveUtils.llstyle
else
    global llstyle = Dict{Symbol, Tuple{Bool, Union{Symbol, Int}}}(
        :default     => (false, :normal), # e.g. comma, equal sign, unknown token
        :comment     => (false, :light_black),
        :label       => (false, :light_red),
        :instruction => ( true, :light_cyan),
        :type        => (false, :cyan),
        :number      => (false, :yellow),
        :bracket     => (false, :yellow),
        :variable    => (false, :normal), # e.g. variable, register
        :keyword     => (false, :light_magenta),
        :funcname    => (false, :light_yellow),
    )
end

const num_regex = r"^(?:\$?-?\d+|0x[0-9A-Fa-f]+|-?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?)$"

function printstyled_ll(io::IO, x, s::Symbol, trailing_spaces="")
    printstyled(io, x, bold=llstyle[s][1], color=llstyle[s][2])
    print(io, trailing_spaces)
end

function code_llvm(io::IO_, @nospecialize(f), @nospecialize(types),
                   raw::Bool, dump_module::Bool=false, optimize::Bool=true,
                   debuginfo::Symbol=:default)
    if VERSION >= v"1.1"
        d = _dump_function(f, types, false, false, !raw, dump_module, :att, optimize,
                           debuginfo)
    else
        d = _dump_function(f, types, false, false, !raw, dump_module, :att, optimize)
    end
    if highlighting[:llvm] && get(io, :color, false)
        print_llvm(io, d)
    else
        print(io, d)
    end
end

function print_llvm(io::IO, code::String)
    buf = IOBuffer(code)
    for line in eachline(buf)
        m = match(r"^(\s*)((?:[^;]|;\")*)(.*)$", line)
        m === nothing && continue
        indent, tokens, comment = m.captures
        print(io, indent)
        print_llvm_tokens(io, tokens)
        printstyled_ll(io, comment, :comment)
        println(io)
    end
end

const llvm_types =
    r"^(?:void|half|float|double|x86_\w+|ppc_\w+|label|metadata|type|opaque|token|i\d+)$"
const llvm_cond = r"^(?:[ou]?eq|[ou]?ne|[uso][gl][te]|ord|uno)$" # true|false

function print_llvm_tokens(io, tokens)
    m = match(r"^((?:[^\s:]+:)?)(\s*)(.*)", tokens)
    if m !== nothing
        label, spaces, tokens = m.captures
        printstyled_ll(io, label, :label, spaces)
    end
    m = match(r"^(%[^\s=]+)(\s*)=(\s*)(.*)", tokens)
    if m !== nothing
        result, spaces, spaces2, tokens = m.captures
        printstyled_ll(io, result, :variable, spaces)
        printstyled_ll(io, '=', :default, spaces2)
    end
    m = match(r"^([a-z]\w*)(\s*)(.*)", tokens)
    if m !== nothing
        inst, spaces, tokens = m.captures
        iskeyword = occursin(r"^(?:define|declare|type)$", inst) || occursin("=", tokens)
        printstyled_ll(io, inst, iskeyword ? :keyword : :instruction, spaces)
    end

    print_llvm_operands(io, tokens)
end

function print_llvm_operands(io, tokens)
    while !isempty(tokens)
        tokens = print_llvm_operand(io, tokens)
    end
    return tokens
end

function print_llvm_operand(io, tokens)
    islabel = false
    while !isempty(tokens)
        m = match(r"^,(\s*)(.*)", tokens)
        if m !== nothing
            spaces, tokens = m.captures
            printstyled_ll(io, ',', :default, spaces)
            break
        end
        m = match(r"^(\*+|=)(\s*)(.*)", tokens)
        if m !== nothing
            sym, spaces, tokens = m.captures
            printstyled_ll(io, sym, :default, spaces)
            continue
        end
        m = match(r"^(\"[^\"]*\")(\s*)(.*)", tokens)
        if m !== nothing
            str, spaces, tokens = m.captures
            printstyled_ll(io, str, :variable, spaces)
            continue
        end
        m = match(r"^([({\[<])(\s*)(.*)", tokens)
        if m !== nothing
            bracket, spaces, tokens = m.captures
            printstyled_ll(io, bracket, :bracket, spaces)
            tokens = print_llvm_operands(io, tokens) # enter
            continue
        end
        m = match(r"^([)}\]>])(\s*)(.*)", tokens)
        if m !== nothing
            bracket, spaces, tokens = m.captures
            printstyled_ll(io, bracket, :bracket, spaces)
            break # leave
        end

        m = match(r"^([^\s,*=(){}\[\]<>]+)(\s*)(.*)", tokens)
        m === nothing && break
        token, spaces, tokens = m.captures
        if occursin(llvm_types, token)
            printstyled_ll(io, token, :type)
            islabel = token == "label"
        elseif occursin(llvm_cond, token) # condition code is instruction-level
            printstyled_ll(io, token, :instruction)
        elseif occursin(num_regex, token)
            printstyled_ll(io, token, :number)
        elseif occursin(r"^@.+$", token)
            printstyled_ll(io, token, :funcname)
        elseif occursin(r"^%.+$", token)
            islabel |= occursin(r"^%[^\d].*$", token) & occursin(r"^\]", tokens)
            printstyled_ll(io, token, islabel ? :label : :variable)
            islabel = false
        elseif occursin(r"^[a-z]\w+$", token)
            printstyled_ll(io, token, :keyword)
        else
            printstyled_ll(io, token, :default)
        end
        print(io, spaces)
    end
    return tokens
end

# code_native
# ===========

function code_native(io::IO_, @nospecialize(f), @nospecialize(types=Tuple);
                     syntax::Symbol=:att, debuginfo::Symbol=:default)
    if VERSION >= v"1.1"
        d = _dump_function(f, types, true, false, false, false, syntax, true, debuginfo)
    else
        d = _dump_function(f, types, true, false, false, false, syntax, true)
    end
    if highlighting[:native] && get(io, :color, false)
        print_native(io, d)
    else
        print(io, d)
    end
end

function print_native(io::IO, code::String, arch::Symbol=sys_arch_category())
    archv = Val(arch)
    buf = IOBuffer(code)
    for line in eachline(buf)
        m = match(r"^(\s*)((?:[^;#/]|#\S|;\"|/[^/])*)(.*)$", line)
        m === nothing && continue
        indent, tokens, comment = m.captures
        print(io, indent)
        print_native_tokens(io, tokens, archv)
        printstyled_ll(io, comment, :comment)
        println(io)
    end
end

function sys_arch_category()
    if Sys.ARCH === :x86_64 || Sys.ARCH === :i686
        :x86
    elseif Sys.ARCH === :aarch64 || startswith(string(Sys.ARCH), "arm")
        :arm
    else
        :unsupported
    end
end

print_native_tokens(io, line, ::Val) = print(io, line)

const x86_ptr = r"^(?:(?:[xyz]mm|[dq])?word|byte|ptr|offset)$"
const avx512flags = r"^(?:z|r[nduz]-sae|sae|1to1?\d)$"
const arm_cond = r"^(?:eq|ne|cs|ho|cc|lo|mi|pl|vs|vc|hi|ls|[lg][te]|al|nv)$"
const arm_keywords = r"^(?:lsl|lsr|asr|ror|rrx|!|/[zm])$"

function print_native_tokens(io, tokens, arch::Union{Val{:x86}, Val{:arm}})
    x86 = arch isa Val{:x86}
    m = match(r"^((?:[^\s:]+:|\"[^\"]+\":)?)(\s*)(.*)", tokens)
    if m !== nothing
        label, spaces, tokens = m.captures
        printstyled_ll(io, label, :label, spaces)
    end
    haslabel = false
    m = match(r"^([a-z][\w.]*)(\s*)(.*)", tokens)
    if m !== nothing
        instruction, spaces, tokens = m.captures
        printstyled_ll(io, instruction, :instruction, spaces)
        haslabel = occursin(r"^(?:bl?|bl?\.\w{2,5}|[ct]bn?z)?$", instruction)
    end

    isfuncname = false
    while !isempty(tokens)
        m = match(r"^([,:*])(\s*)(.*)", tokens)
        if m !== nothing
            sym, spaces, tokens = m.captures
            printstyled_ll(io, sym, :default, spaces)
            isfuncname = false
            continue
        end
        m = match(r"^([(){}\[\]])(\s*)(.*)", tokens)
        if m !== nothing
            bracket, spaces, tokens = m.captures
            printstyled_ll(io, bracket, :bracket, spaces)
            continue
        end
        m = match(r"^#([0-9a-fx.-]+)(\s*)(.*)", tokens)
        if !x86 && m !== nothing && occursin(num_regex, m.captures[1])
            num, spaces, tokens = m.captures
            printstyled_ll(io, "#" * num, :number, spaces)
            continue
        end

        m = match(r"^([^\s,:*(){}\[\]][^\s,:*/(){}\[\]]*)(\s*)(.*)", tokens)
        m === nothing && break
        token, spaces, tokens = m.captures
        if occursin(num_regex, token)
            printstyled_ll(io, token, :number)
        elseif x86 && occursin(x86_ptr, token) || occursin(avx512flags, token)
            printstyled_ll(io, token, :keyword)
            isfuncname = token == "offset"
        elseif !x86 && (occursin(arm_keywords, token) || occursin(arm_cond, token))
            printstyled_ll(io, token, :keyword)
        elseif occursin(r"^L.+$", token)
            printstyled_ll(io, token, :label)
        elseif occursin(r"^\$.+$", token)
            printstyled_ll(io, token, :funcname)
        elseif occursin(r"^%?(?:[a-z][\w.]+|\"[^\"]+\")$", token)
            islabel = haslabel & !occursin(',', tokens)
            printstyled_ll(io, token, islabel ? :label : isfuncname ? :funcname : :variable)
            isfuncname = false
        else
            printstyled_ll(io, token, :default)
        end
        print(io, spaces)
    end
end

end # module

max_stars_repo_name: kimikage/ColoredLLCodes

max_forks_repo_name: kimikage/ColoredLLCodes


text:
using PowerAnalyses
using Test: @testset, @test

@testset "t-test" begin
    es = 0.5
    alpha = 0.05
    power = 0.95
    n = 50
    N = [n, n]

    # Values are obtained via `pwr` version 1.3 and `G*Power` version 3.1.9.7.
    @test get_alpha(OneSampleTTest(two_tails); es, power, n) ≈ 0.067 atol=0.001
    @test get_alpha(OneSampleTTest(one_tail); es, power, n) ≈ 0.033 atol=0.001
    @test get_power(OneSampleTTest(two_tails); es, alpha, n) ≈ 0.933 atol=0.001
    @test get_power(OneSampleTTest(one_tail); es, alpha, n) ≈ 0.967 atol=0.001

    # For methods which require `find_zero`, compare the outcome with pwr; G*Power is far off.
    @test get_es(OneSampleTTest(two_tails); alpha, power, n) ≈ 0.520 atol=0.001
    @test get_es(OneSampleTTest(one_tail); alpha, power, n) ≈ 0.471 atol=0.001
    @test get_n(OneSampleTTest(two_tails); alpha, power, es) ≈ 53.941 atol=0.001
    @test get_n(OneSampleTTest(one_tail); alpha, power, es) ≈ 44.679 atol=0.001

    # G*Power gives 0.77; they seem to calculate the noncentrality_parameter differently.
    # pwr.t.test(n=50, d=0.5, sig.level=NULL, power=0.95, type="two.sample", alternative="two.sided")
    @test get_alpha(IndependentSamplesTTest(two_tails); es, power, n) ≈ 0.392 atol=0.01
    # Same as the one sample t-test.
    @test get_alpha(DependentSamplesTTest(two_tails); es, power, n) ≈ 0.067 atol=0.01

    n_predictors = 8
    # Beware that G*Power assumes input to be f² and not f.
    @test get_alpha(DeviationFromZeroMultipleRegression(n_predictors); es=sqrt(es), power, n) ≈ 0.072 atol=0.1
    n_tested_predictors = 6
    @test get_alpha(IncreaseMultipleRegression(n_predictors, n_tested_predictors); es=sqrt(es), power, n) ≈ 0.043 atol=0.1
    n_groups = 2
    @test get_alpha(OneWayANOVA(n_groups); es, power, n) ≈ 0.067 atol=0.001
    df = 10
    @test get_alpha(MultifactorFixedEffectsANOVA(n_groups, df); es, power, n) ≈ 0.461 atol=0.01

    n_variables = 2
    @test get_alpha(ConstantVectorHotellingTsqTest(n_variables); es, power, n) ≈ 0.132 atol=0.001
    @test get_alpha(TwoVectorsHotellingTsqTest(n_variables); es, power, n=N) ≈ 0.511 atol=0.001

    df = 5
    @test get_power(GoodnessOfFitChisqTest(df); es, alpha, n) ≈ 0.787 atol=0.001
    @test get_alpha(GoodnessOfFitChisqTest(df); es, power, n) ≈ 0.253 atol=0.001
    @test get_es(GoodnessOfFitChisqTest(df); alpha, power, n) ≈ 0.629 atol=0.001
    @test get_n(GoodnessOfFitChisqTest(df); alpha, power, es) ≈ 79.12 atol=0.001

    @test get_alpha(ConstantVarianceChisqTest(one_tail); es, power, n) ≈ 0.038 atol=0.02

    # TODO: Add test for PointBiseralTTest
end


text:
module Embed

using Distributions

struct Embedding
    embedding::AbstractArray{Float64}
end

"""
    embed(ts::Vector{Float64}, E::Int, tau::Int)

Embed a time series `ts` in `E` dimensions with embedding lag `tau`.
"""
function embed(ts::Vector{Float64}, E::Int, tau::Int)
  n::Int = length(ts)

  # Initialize
  embedded_ts = zeros(Float64, E, n - ((E - 1) * tau))
  l::Int = n - ((E - 1) * tau)

  for i in 1:E
    start_index = 1 + (i - 1) * tau
    stop_index  = start_index + l - 1
    embedded_ts[i, :] = ts[start_index:stop_index]
  end
  return copy(transpose(embedded_ts))
end

"""
    embed(ts::Vector{Float64}, E::Int, tau::Int)

Embed a time series `ts` in `E` dimensions with embedding lag `tau`.
"""
function embedding(ts::Vector{Float64}, E::Int, tau::Int)
  n::Int = length(ts)

  # Initialize
  embedded_ts = zeros(Float64, E, n - ((E - 1) * tau))
  l::Int = n - ((E - 1) * tau)

  for i in 1:E
    start_index = 1 + (i - 1) * tau
    stop_index  = start_index + l - 1
    embedded_ts[i, :] = ts[start_index:stop_index]
  end
  return Embedding(copy(transpose(embedded_ts)))
end
"""
    invariantset(embedding, tolerance)

Determine whether an embedding forms an invariant set.
"""
function invariantset(embedding)
    lastpoint = embedding[end, :]
    dim = size(embedding, 2)

    # Triangulate the embedding using all points but the last
    triangulation = Simplices.Delaunay.delaunayn(embedding[1:end-1, :])

    points = embedding[1:end-1, :]
    simplex_indices = triangulation

    # Centroids and radii of simplices in the triangulation
    centroids, radii = centroids_radii2(points, simplex_indices)

    lastpoint_matrix = repeat(lastpoint, size(centroids, 1), 1)

    # Find simplices that can contain the last point (not all can)
    dists_lastpoint_and_centroids = sum((lastpoint_matrix - centroids).^2, dims=2)
    distdifferences = radii.^2 - dists_lastpoint_and_centroids

    # Find the row indices of the simplices that possibly contain the last point (meaning that
    # dist(simplex_i, lastpoint) <= radius(simplex), so the circumsphere of the simplex
    # contains the last point.
    valid_simplex_indices = findall((heaviside0(distdifferences) .* collect(1:size(triangulation[2], 1))) .> 0)

    n_validsimplices = size(valid_simplex_indices, 1)

    # Loop over valid simplices and check whether the corresponding simplex actually
    # contains the last point.
    i = 1
    lastpoint_contained = false

    # Continue checking if simplices contain the last point until some simplex contains it.
    # Then the point must necessarily be inside the convex hull of the triangulation formed
    # by those simplices.
    while i <= n_validsimplices && !lastpoint_contained
        # Subsample the valid simplex and compute its orientation
        simplex = embedding[simplex_indices[valid_simplex_indices[i], :], :]
        orientation_simplex = det([ones(dim + 1, 1) simplex])

        beta = 1 # Convex expansion coefficient

        j = 1
        while j <= dim + 1 && beta >= 0
            tmp = copy(simplex)
            tmp[j, :] = copy(lastpoint)
            beta = det([ones(dim + 1, 1) tmp]) * sign(orientation_simplex)
            j = j + 1
        end

        # If the last convex expansion coefficient is positive, the last point is contained
        # in the triangulation (because all the previous coefficients must have been nonnegative)
        if beta >= 0
            lastpoint_contained = true
        end

        i = i + 1
    end
    # If the last point is contained in the triangulation, the set is invariant.
    return lastpoint_contained
end

function gaussian_embedding(npts::Int; covariance::Float64 = 0.4)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        source = rand(dist, npts, 1)
        dest = covariance .* source[1:end] .+ (1.0 - covariance) .* rand(dist, npts, 1)
        # Embedding
        embedd = hcat(source[1:end-1], source[2:end], dest[2:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return Embedding(embedd)
        end
    end
end

function gaussian_embedding(npts::Int, covariance::Float64)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        source = rand(dist, npts, 1)
        dest = covariance .* source[1:end] .+ (1.0 - covariance) .* rand(dist, npts, 1)
        # Embedding
        embedd = hcat(source[1:end-1], source[2:end], dest[2:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return Embedding(embedd)
        end
    end
end

function gaussian_embedding_arr(npts::Int, covariance::Float64 = 0.4)
    dist = Normal()
    # Create an uncorrelated source and target
    invariant_embedding_found = false
    while !invariant_embedding_found
        source = rand(dist, npts, 1)
        dest = covariance .* source[1:end] .+ (1.0 - covariance) .* rand(dist, npts, 1)
        # Embedding
        embedd = hcat(source[1:end-1], source[2:end], dest[2:end])

        if invariantset(embedd)
            invariant_embedding_found = true
            return embedd
        end
    end
end

"""
    embedding_example(n_points::Int, m::Int, tau::Int)

Create an example embedding consisting of `n_points` points in `E` dimension space,
using embedding lag `tau`.
"""
function embedding_example(n_points::Int, E::Int, tau::Int)
  ts = randn(n_points)
  embedding = embed(ts, E, tau)
end

"""
    embedding_example(n_points::Int, m::Int, tau::Int)

Create an example embedding consisting of `n_points` points in `E` dimension space,
using embedding lag `tau`.
"""
function embedding_ex(n_points::Int, E::Int, tau::Int)
  ts = randn(n_points)
  embedding = embedding(ts, E, tau)
end

export embed, embedding_example,  gaussian_embedding, gaussian_embedding_arr

end