text:
##############################################################################
#
# Metadata field getters and setters
#
##############################################################################

import Languages.name

title(d::AbstractDocument) = d.metadata.title
language(d::AbstractDocument) = d.metadata.language
author(d::AbstractDocument) = d.metadata.author
timestamp(d::AbstractDocument) = d.metadata.timestamp

function title!(d::AbstractDocument, nv::AbstractString)
    d.metadata.title = nv
end

function language!(d::AbstractDocument, nv::T) where T <: Language
    d.metadata.language = nv
end

function author!(d::AbstractDocument, nv::AbstractString)
    d.metadata.author = nv
end

function timestamp!(d::AbstractDocument, nv::AbstractString)
    d.metadata.timestamp = nv
end

##############################################################################
#
# Vectorized getters for an entire Corpus
#
##############################################################################

titles(c::Corpus) = map(d -> title(d), documents(c))
languages(c::Corpus) = map(d -> language(d), documents(c))
authors(c::Corpus) = map(d -> author(d), documents(c))
timestamps(c::Corpus) = map(d -> timestamp(d), documents(c))

titles!(c::Corpus, nv::AbstractString) = title!.(documents(c), nv)
languages!(c::Corpus, nv::T) where {T <: Language} = language!.(documents(c), Ref(nv)) #Ref to force scalar broadcast
authors!(c::Corpus, nv::AbstractString) = author!.(documents(c), Ref(nv))
timestamps!(c::Corpus, nv::AbstractString) = timestamp!.(documents(c), Ref(nv))

function titles!(c::Corpus, nvs::Vector{String})
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        title!(d, nvs[i])
    end
end

function languages!(c::Corpus, nvs::Vector{T}) where T <: Language
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        language!(d, nvs[i])
    end
end

function authors!(c::Corpus, nvs::Vector{String})
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        author!(d, nvs[i])
    end
end

function timestamps!(c::Corpus, nvs::Vector{String})
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        timestamp!(d, nvs[i])
    end
end


text:
@testset "Point distances" begin
    p1_2d = SVector(0.0, 0.0)
    p2_2d = SVector(1.0, 1.0)
    @test distance(p1_2d, p2_2d) == sqrt(2.0)

    p1_3d = SVector(0.0, 0.0, 1.0)
    p2_3d = SVector(1.0, 1.0, 2.0)
    @test distance(p1_3d, p2_3d) == sqrt(3.0)

    @test distance(p1_3d, p2_2d) == sqrt(2.0)
    @test distance(p1_2d, p2_3d) == sqrt(2.0)
end

@testset "Catenary distances" begin

    # z = cosh(x) - 1 for x ∈ -1.0 .. 1.0
    cat = Catenary(IdentityTransformation(), -1.0, 1.0, 1.0)

    @test distance(cat, SVector(0.0, 0.0, 0.0)) === 0.0
    @test distance(SVector(0.0, 0.0, 0.0), cat) === 0.0

    @test powerline_distances(cat, SVector(0.0, 0.0, 0.0)) === (0.0, 0.0, 0.0, 0.0)
    @test powerline_distances(SVector(0.0, 0.0, 0.0), cat) === (0.0, 0.0, 0.0, 0.0)

    # Distance to catenary body
    # below
    @test powerline_distances(cat, SVector(0.5, 0.0, 0.0)) === (0.11454038802773013, 0.0, -(cosh(0.5)-1), 0.5)
    @test powerline_distances(cat, SVector(-0.5, 0.0, 0.0)) === (0.11454038802773013, 0.0, -(cosh(0.5)-1), -0.5)
    @test powerline_distances(cat, SVector(0.5, 0.5, 0.0)) === (nextfloat(sqrt(0.5^2 + 0.11454038802773013^2)), 0.5, -(cosh(0.5)-1), 0.5)
    @test powerline_distances(cat, SVector(-0.5, 0.5, 0.0)) === (nextfloat(sqrt(0.5^2 + 0.11454038802773013^2)), 0.5, -(cosh(0.5)-1), -0.5)
    @test powerline_distances(cat, SVector(0.5, -0.5, 0.0)) === (nextfloat(sqrt(0.5^2 + 0.11454038802773013^2)), -0.5, -(cosh(0.5)-1), 0.5)
    @test powerline_distances(cat, SVector(-0.5, -0.5, 0.0)) === (nextfloat(sqrt(0.5^2 + 0.11454038802773013^2)), -0.5, -(cosh(0.5)-1), -0.5)

    # above
    @test powerline_distances(cat, SVector(0.5, 0.0, 0.5)) === (0.31834847276128164, 0.0, -(cosh(0.5)-1.5), 0.5)
    @test powerline_distances(cat, SVector(-0.5, 0.0, 0.5)) === (0.31834847276128164, 0.0, -(cosh(0.5)-1.5), -0.5)
    @test powerline_distances(cat, SVector(0.5, 0.5, 0.5)) === (nextfloat(sqrt(0.5^2 + 0.31834847276128164^2)), 0.5, -(cosh(0.5)-1.5), 0.5)
    @test powerline_distances(cat, SVector(-0.5, 0.5, 0.5)) === (nextfloat(sqrt(0.5^2 + 0.31834847276128164^2)), 0.5, -(cosh(0.5)-1.5), -0.5)
    @test powerline_distances(cat, SVector(0.5, -0.5, 0.5)) === (nextfloat(sqrt(0.5^2 + 0.31834847276128164^2)), -0.5, -(cosh(0.5)-1.5), 0.5)
    @test powerline_distances(cat, SVector(-0.5, -0.5, 0.5)) === (nextfloat(sqrt(0.5^2 + 0.31834847276128164^2)), -0.5, -(cosh(0.5)-1.5), -0.5)

    # distance to catenary end points
    # left
    @test powerline_distances(cat, SVector(-1.5, 0.0, 0.5430806348152437)) === (0.5, 0.0, 0.0, -1.5)

    # right
    @test powerline_distances(cat, SVector(1.5, 0.0, 0.5430806348152437)) === (0.5, 0.0, 0.0, 1.5)
end

ext: jl


text:
struct Sample{A<:Random.AbstractRNG,M<:AbstractModel,S<:AbstractSampler,K} <: Transducers.Transducer
    rng::A
    model::M
    sampler::S
    kwargs::K
end

function Sample(model::AbstractModel, sampler::AbstractSampler; kwargs...)
    return Sample(Random.GLOBAL_RNG, model, sampler; kwargs...)
end

"""
    Sample([rng, ]model, sampler; kwargs...)

Create a transducer that returns samples from the `model` with the Markov chain Monte Carlo
`sampler`.

# Examples

```jldoctest; setup=:(using AbstractMCMC: Sample)
julia> struct MyModel <: AbstractMCMC.AbstractModel end

julia> struct MySampler <: AbstractMCMC.AbstractSampler end

julia> function AbstractMCMC.step(rng, ::MyModel, ::MySampler, state=nothing; kwargs...)
           # all samples are zero
           return 0.0, state
       end

julia> transducer = Sample(MyModel(), MySampler());

julia> collect(transducer(1:10)) == zeros(10)
true
```
"""
function Sample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler;
    kwargs...
)
    return Sample(rng, model, sampler, kwargs)
end

function Transducers.start(rf::Transducers.R_{<:Sample}, result)
    sampler = Transducers.xform(rf)
    return Transducers.wrap(
        rf,
        step(sampler.rng, sampler.model, sampler.sampler; sampler.kwargs...),
        Transducers.start(Transducers.inner(rf), result),
    )
end

function Transducers.next(rf::Transducers.R_{<:Sample}, result, input)
    t = Transducers.xform(rf)
    Transducers.wrapping(rf, result) do (sample, state), iresult
        iresult2 = Transducers.next(Transducers.inner(rf), iresult, sample)
        return step(t.rng, t.model, t.sampler, state; t.kwargs...), iresult2
    end
end

function Transducers.complete(rf::Transducers.R_{Sample}, result)
    _private_state, inner_result = Transducers.unwrap(rf, result)
    return Transducers.complete(Transducers.inner(rf), inner_result)
end


text:
using LinearAlgebra
using Einsum
using ASE
using JuLIP
using JuLIP.Potentials

using JuLIPMaterials.CLE

println("=========================")
println(" Testing Elasticity      ")
println("=========================")

println("Testing Voigt conversions")
σ, ϵ = rand(6), rand(6)
σₘ = cauchy_stress(σ)
ϵₘ = cauchy_strain(ϵ)
@test voigt_stress(σₘ) ≈ σ
@test voigt_strain(ϵₘ) ≈ ϵ

# check invariance of scalar product σ⋅ϵ = sum_ij σ_ij⋅ϵ_ij =
@test  σ' * ϵ ≈ sum([σₘ[i,j] * ϵₘ[i,j] for i=1:3, j=1:3])

C11, C12, C44 = rand(3)
C = voigt_moduli(C11, C12, C44)
c = elastic_moduli(C)

# check invertability of c_ijkl <--> C_ij conversion
@test voigt_moduli(c) ≈ C

# check σ_ij = c_ijkl ϵ_kl <--> σ = C⋅ϵ
ϵ = rand(6)
ϵₘ = cauchy_strain(ϵ)
σ = zeros(3,3)
@einsum σ[i,j] = c[i,j,k,l] * ϵₘ[k,l]
@test cauchy_stress(C * ϵ) ≈ σ
@test voigt_stress(σ) ≈ C * ϵ

println("Testing `elastic_moduli`")

eam = JuLIP.Potentials.EAM(test_pots * "w_eam4.fs")
at = bulk(:W)
set_calculator!(at, eam)
variablecell!(at)
minimise!(at)
E0 = energy(at)
c = elastic_moduli(eam, at)
C = voigt_moduli(c)

# apply small random strain and check CLE prediction holds
ϵ = 1e-4*rand(6)
ϵₘ = cauchy_strain(ϵ)
apply_defm!(at, I(3) + ϵₘ)
σₘ = stress(at)
σ = voigt_stress(σₘ)

# check elastic energy matches potential energy
E = energy(at)
@test E - E0 ≈ 0.5 * ϵ' * σ * volume(at) atol=1e-7

# Hooke's law in Voigt form: σ = C ϵ
@test σₘ ≈ cauchy_stress(C * ϵ) atol=1e-5

# check Hooke's law in tensor form: σ_ij = c_ijkl ϵ_kl
σₜ = zeros(3,3)
@einsum σₜ[i,j] = c[i,j,k,l] * ϵₘ[k,l]
@test σₜ ≈ σₘ atol=1e-3

# check C has cubic symmetry and that C11, C12, C44 elastic constants match expected values for this potential
@show cubic_moduli(C, tol=1e-3)
@test all(isapprox.(cubic_moduli(C, tol=1e-3), [3.2644382942616716, 1.2618231295748539, 1.004134784461818], atol=1e-3))

println("Testing rotation of elastic moduli")

A = rotation_matrix(x=[1,1,1], y=[1,-1,0])
at = bulk(:W)
set_calculator!(at, eam)
c = elastic_moduli(at)
c_rot = rotate_moduli(c, A)
C = voigt_moduli(at)
C_rot = rotate_moduli(C, A)

# now we rotate the Atoms directly and compare
apply_defm!(at, A)
c_direct = elastic_moduli(at)
C_direct = voigt_moduli(at)

@test c_rot ≈ c_direct atol=1e-3
@test C_rot ≈ C_direct atol=1e-3


text:
"""
Comprises everything needed for path following: `Project`, `SystemCore`,
`ContinuationMethod`, and `Visualization`. Sessions can be managed via `create`, `save`,
and `load`.
"""
type Session
	P::Project
	core::SystemCore
	cont::ContinuationMethod
	viz::Visualization
	Session() = new()
end


text:
@testset "meta information for evaluation" begin

doc_body = """
```julia
include("test_include.jl")
```

```julia
@__MODULE__
```

```julia
@__DIR__
```

```julia
@__FILE__
```

```julia
@__LINE__ # broken
```

```julia
read("./test_include.jl", String)
```
"""
doc_dir = normpath(@__DIR__, "..", "mocks")
doc_path = normpath(doc_dir, "test_meta.jmd")
write(doc_path, doc_body)

script_line = ":include_me"
script_body = "$script_line"
script_path = normpath(@__DIR__, "..", "mocks", "test_include.jl")
write(script_path, script_body)

m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
mock = run_doc(WeaveDoc(doc_path); mod = m)
check_output(i, s) = occursin(s, mock.chunks[i].output)

@test check_output(1, script_line)
@test check_output(2, string(m))
@test check_output(3, doc_dir)
@test check_output(4, doc_path)
@test_broken check_output(5, 18)
@test check_output(6, string('"', script_line, '"')) # current working directory

end  # @testset "meta information for evaluation"


text:
using MaximalCouplings
using Distributions
using Test
using Random
using Random: rand, rand!
using LinearAlgebra: norm
import HypothesisTests
import StatsBase

function test_distr(x, d::ContinuousDistribution)
    test = HypothesisTests.ApproximateOneSampleKSTest(x, d)
    pval = HypothesisTests.pvalue(test)
    return pval
end
function test_distr(x, d::DiscreteDistribution)
    theta0 = pdf.(d, minimum(x):maximum(x))
    theta0 ./= sum(theta0)
    contingency = StatsBase.counts(x)
    test = HypothesisTests.ChisqTest(contingency, theta0)
    pval = HypothesisTests.pvalue(test)
    return pval
end

# For two distributions that are the same, the maximal coupling should
# be perfect.
@testset "Equal Distributions" begin
    c_equal = MaximalCoupling(Normal(1,5), Normal(1,5))
    @test prob_couple(c_equal) ≈ 1.0
    @testset "samples" begin
        for _ in 1:1000
            x, y = rand(c_equal)
            @test x == y
        end
    end # testset
end # testset

# Probabilistic test that the probability given by `prob_couple`
# matches the fraction of samples yielded by `rand` that are coupled.
@testset "Probability of Coupling" begin
    coup = MaximalCoupling(Normal(1,2), Gamma(1, 2))
    nsamples = 10^7
    X, Y = rand(coup, nsamples)
    pcouple = prob_couple(coup)
    frac_coupled = mean(X .== Y)
    se = sqrt(pcouple  * (1-pcouple) / nsamples)
    deviation = (pcouple - frac_coupled) / se
    @test abs(deviation) < 5
    @test eltype(coup) == Float64
    @test eltype(coup) == eltype(X)
    @test eltype(coup) == eltype(Y)
    @test eltype(coup) == eltype(rand(coup)[1])
    @test eltype(coup) == eltype(rand(coup)[2])
end # testset

@testset "Multivariate Couplings" begin
    d1=MultivariateNormal([1., 1.], [[3., 1.]'; [1., 4.]'])
    d2=MvLogNormal([0., 0.], [[1., 0.1]'; [0.1, 1.]'])
    coup = MaximalCoupling(d1, d2)
    # Dimensionality checks
    @test length(coup) == 2
    @test length(rand(coup)[1]) == 2
    @test length(rand(coup)[2]) == 2
    @test size(rand(coup, 5)[1]) == size(rand(d1, 5))
    @test size(rand(coup, 5)[2]) == size(rand(d2, 5))
    @test size(rand(coup, 5)[1]) == (2, 5)
    @test size(rand(coup, 5)[2]) == (2, 5)

    X, Y = rand(coup, 10^4)
    equal13 = X[1, :] .== Y[1, :]
    equal24 = X[2, :] .== Y[2, :]
    @test all(equal13 .== equal24)
    # test that some samples are coupled
    @test sum(equal13) > 0
    @test eltype(coup) == Float64
    @test eltype(coup) == eltype(X)
    @test eltype(coup) == eltype(Y)
    @test eltype(coup) == eltype(rand(coup)[1])
    @test eltype(coup) == eltype(rand(coup)[2])

    # KS test on the marginals
    @test test_distr(X[1,:], Normal(mean(d1)[1], √d1.Σ.mat[1,1])) > 1e-5
    @test test_distr(X[2,:], Normal(mean(d1)[2], √d1.Σ.mat[2,2])) > 1e-5
    @test test_distr(Y[1,:], LogNormal(d2.normal.μ[1], √d2.normal.Σ.mat[1,1])) > 1e-5
    @test test_distr(Y[2,:], LogNormal(d2.normal.μ[2], √d2.normal.Σ.mat[2,2])) > 1e-5
end
@testset "Discrete Couplings - finite+finite support " begin
    d1 = Categorical([0.05, 0.05, 0.2, 0.3, 0.3, 0.1])
    d2 = Categorical([0.1, 0.2, 0.3, 0.4])
    coup = MaximalCoupling(d1, d2)
    # Dimensionality checks
    @test length(coup) == 1
    @test length(rand(coup)[1]) == 1
    @test length(rand(coup)[2]) == 1
    @test size(rand(coup, 5)[1]) == size(rand(d1, 5))
    @test size(rand(coup, 5)[2]) == size(rand(d2, 5))
    @test size(rand(coup, 5)[1]) == (5,)
    @test size(rand(coup, 5)[2]) == (5,)
    pcouple = prob_couple(coup)

    nsamples = 10^7
    X, Y = rand(coup, nsamples)
    frac_coupled = mean(X .== Y)
    se = sqrt(pcouple  * (1-pcouple) / nsamples)
    deviation = (pcouple - frac_coupled) / se
    @test abs(deviation) < 5
    @test eltype(coup) == Int64
    @test eltype(coup) == eltype(X)
    @test eltype(coup) == eltype(Y)
    @test eltype(coup) == eltype(rand(coup)[1])
    @test eltype(coup) == eltype(rand(coup)[2])

    @test test_distr(X, d1) > 1e-5
    @test test_distr(Y, d2) > 1e-5
end

@testset "Discrete Couplings - infinite+finite support " begin
    d1 = Poisson(4.2)
    d2 = Categorical([0.1, 0.2, 0.3, 0.4])
    coup = MaximalCoupling(d1, d2)
    # Dimensionality checks
    @test length(coup) == 1
    @test length(rand(coup)[1]) == 1
    @test length(rand(coup)[2]) == 1
    @test size(rand(coup, 5)[1]) == size(rand(d1, 5))
    @test size(rand(coup, 5)[2]) == size(rand(d2, 5))
    @test size(rand(coup, 5)[1]) == (5,)
    @test size(rand(coup, 5)[2]) == (5,)
    pcouple = prob_couple(coup)

    nsamples = 10^7
    X, Y = rand(coup, nsamples)
    frac_coupled = mean(X .== Y)
    se = sqrt(pcouple  * (1-pcouple) / nsamples)
    deviation = (pcouple - frac_coupled) / se
    @test abs(deviation) < 5
    @test eltype(coup) == Int64
    @test eltype(coup) == eltype(X)
    @test eltype(coup) == eltype(Y)
    @test eltype(coup) == eltype(rand(coup)[1])
    @test eltype(coup) == eltype(rand(coup)[2])

    @test test_distr(X, d1) > 1e-5
    @test test_distr(Y, d2) > 1e-5
end

@testset "Discrete Couplings - finite+infinite support" begin
    # just switch the order
    d1 = Categorical([0.1, 0.2, 0.3, 0.4])
    d2 = Poisson(4.2)
    coup = MaximalCoupling(d1, d2)
    # Dimensionality checks
    @test length(coup) == 1
    @test length(rand(coup)[1]) == 1
    @test length(rand(coup)[2]) == 1
    @test size(rand(coup, 5)[1]) == size(rand(d1, 5))
    @test size(rand(coup, 5)[2]) == size(rand(d2, 5))
    @test size(rand(coup, 5)[1]) == (5,)
    @test size(rand(coup, 5)[2]) == (5,)
    pcouple = prob_couple(coup)

    nsamples = 10^7
    X, Y = rand(coup, nsamples)
    frac_coupled = mean(X .== Y)
    se = sqrt(pcouple  * (1-pcouple) / nsamples)
    deviation = (pcouple - frac_coupled) / se
    @test abs(deviation) < 5
    @test eltype(coup) == Int64
    @test eltype(coup) == eltype(X)
    @test eltype(coup) == eltype(Y)
    @test eltype(coup) == eltype(rand(coup)[1])
    @test eltype(coup) == eltype(rand(coup)[2])
end

@testset "Discrete Couplings - infinite support" begin
    # both infinite support
    d1 = Geometric(0.1)
    d2 = Poisson(4.2)
    coup = MaximalCoupling(d1, d2)
    # Dimensionality checks
    @test length(coup) == 1
    @test length(rand(coup)[1]) == 1
    @test length(rand(coup)[2]) == 1
    @test size(rand(coup, 5)[1]) == size(rand(d1, 5))
    @test size(rand(coup, 5)[2]) == size(rand(d2, 5))
    @test size(rand(coup, 5)[1]) == (5,)
    @test size(rand(coup, 5)[2]) == (5,)
    pcouple = prob_couple(coup)

    nsamples = 10^7
    X, Y = rand(coup, nsamples)
    frac_coupled = mean(X .== Y)
    se = sqrt(pcouple  * (1-pcouple) / nsamples)
    deviation = (pcouple - frac_coupled) / se
    @test abs(deviation) < 5
    @test eltype(coup) == Int64
    @test eltype(coup) == eltype(X)
    @test eltype(coup) == eltype(Y)
    @test eltype(coup) == eltype(rand(coup)[1])
    @test eltype(coup) == eltype(rand(coup)[2])
end

@testset "Multivariate Normal" begin
    mvn1=MultivariateNormal([1., 1.], [[3., 1.]'; [1., 4.]'])
    mvn2=MultivariateNormal([1.5, 0.8], [[3., 1.]'; [1., 4.]'])
    @test MaximalCouplings.samecov(mvn1, mvn2)

    coup = MaximalCoupling(mvn1, mvn2)
    # Dimensionality checks
    @test length(coup) == 2
    @test length(rand(coup)[1]) == 2
    @test length(rand(coup)[2]) == 2
    @test size(rand(coup, 5)[1]) == size(rand(mvn1, 5))
    @test size(rand(coup, 5)[2]) == size(rand(mvn2, 5))
    @test size(rand(coup, 5)[1]) == (2, 5)
    @test size(rand(coup, 5)[2]) == (2, 5)

    nsamples = 10^5
    pcouple_BouRabee = mean(isequal(MaximalCouplings._rand_BouRabee(mvn1, mvn2)...) for _ in 1:nsamples)
    pcouple_continus = mean(isequal(MaximalCouplings._rand_continuous(mvn1, mvn2)...) for _ in 1:nsamples)
    pcouple_diff = abs(pcouple_BouRabee-pcouple_continus)

    pcouple_SE = sqrt(pcouple_continus * (1-pcouple_continus) / nsamples)
    @test pcouple_diff < pcouple_SE*10

    # test mean distance is lower with Bou-Rabee by a significant amount
    dist(x1,x2) = norm(x1-x2)
    mdist_BouRabee = mean(dist(MaximalCouplings._rand_BouRabee(mvn1, mvn2)...) for _ in 1:nsamples)
    mdist_continus = mean(dist(MaximalCouplings._rand_continuous(mvn1, mvn2)...) for _ in 1:nsamples)
    vdist_BouRabee = var(dist(MaximalCouplings._rand_BouRabee(mvn1, mvn2)...) for _ in 1:nsamples)
    vdist_continus = var(dist(MaximalCouplings._rand_continuous(mvn1, mvn2)...) for _ in 1:nsamples)
    mdist_SE = sqrt( (vdist_BouRabee + vdist_continus)/2 / nsamples)
    @test mdist_BouRabee < mdist_continus - 4*mdist_SE

    # test samples come from the correct distribution
    nsamples = 10^5
    X, Y = rand(coup, nsamples)
    marginal(d, j) = Normal(d.μ[j], √d.Σ.mat[j,j])
    @test test_distr(X[1,:], marginal(mvn1, 1)) > 1e-5
    @test test_distr(X[2,:], marginal(mvn1, 2)) > 1e-5
    @test test_distr(Y[1,:], marginal(mvn2, 1)) > 1e-5
    @test test_distr(Y[2,:], marginal(mvn2, 2)) > 1e-5
end


text:
using GLMdiagnostics
using Test
using Distributions
using Random
using DataFrames
using StatsBase
using GLM
using Gadfly
using Statistics
using Compose

##############################################
# - Simulate data and linear model
##############################################

# Source simulated linear model to test functions
include("test_simulated_model.jl")

##############################################
# - Run unit tests
##############################################

@testset "GLMdiagnostics.jl" begin

# Each function contains unit tests in its own specific .jl files
# - Source unit tests for each function

# Source tests for plot_qq function
include("test_plot_qq.jl")

# Source tests for plot_fit_res function
include("test_plot_fit_res.jl")

# Source tests for plot_scale_loc function
include("test_plot_scale_loc.jl")

end

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
function get_status_string(status::Status, status_string::Ref{Cstring})
    ccall((:amd_dbgapi_get_status_string, "librocm-dbgapi"), Status,
          (Status, Ref{Cstring}), status, status_string)
end

function get_version(major::Ref{UInt32}, minor::Ref{UInt32}, patch::Ref{UInt32})
    ccall((:amd_dbgapi_get_version, "librocm-dbgapi"), Cvoid,
          (Ref{UInt32}, Ref{UInt32}, Ref{UInt32}), major, minor, patch)
end

function get_build_name()
    ccall((:amd_dbgapi_get_build_name, "librocm-dbgapi"), Cstring, ())
end

function initialize(callbacks::Ref{Callbacks})
    ccall((:amd_dbgapi_initialize, "librocm-dbgapi"), Status,
          (Ref{Callbacks},), callbacks)
end

function finalize()
    ccall((:amd_dbgapi_finalize, "librocm-dbgapi"), Status, ())
end

function architecture_get_info(architecture_id::ArchitectureID, query::ArchitectureInfo,
                               value_size::Csize_t, value::Ref{Cvoid})
    ccall((:amd_dbgapi_architecture_get_info, "librocm-dbgapi"), Status,
          (ArchitectureID, ArchitectureInfo, Csize_t, Ref{Cvoid}),
          architecture_id, query, value_size, value)
end

function get_architecture(elf_amdgpu_machine::UInt32, architecture_id::Ref{ArchitectureID})
    ccall((:amd_dbgapi_get_architecture, "librocm-dbgapi"), Status,
          (UInt32, Ref{ArchitectureID}), elf_amdgpu_machine, architecture_id)
end

function disassemble_instruction(architecture_id::ArchitectureID, address::GlobalAddress,
                                 size::Ref{UInt64}, memory::Ref{Cvoid},
                                 instruction_text::Ref{Cstring},
                                 address_operand_count::Ref{Csize_t},
                                 address_operands::Ref{Ptr{GlobalAddress}})

    ccall((:amd_dbgapi_disassemble_instruction, "librocm-dbgapi"), Status,
          (ArchitectureID, GlobalAddress, Ref{UInt64}, Ref{Cvoid}, Ref{Cstring},
           Ref{Csize_t}, Ref{Ptr{GlobalAddress}}), architecture_id, address, size, memory,
          instruction_text, address_operand_count, address_operands)
end

function process_get_info(process_id::ProcessID, query::ProcessInfo, value_size::Csize_t,
                          value::Ref{Cvoid})

    ccall((:amd_dbgapi_process_get_info, "librocm-dbgapi"), Status,
          (ProcessID, ProcessInfo, Csize_t, Ref{Cvoid}),
          process_id, query, value_size, value)
end

function process_attach(client_process_id::ClientProcessID, process_id::Ref{ProcessID})
    ccall((:amd_dbgapi_process_attach, "librocm-dbgapi"), Status,
          (ClientProcessID, Ref{ProcessID}), client_process_id, process_id)
end

function process_detach(process_id::ProcessID)
    ccall((:amd_dbgapi_process_detach, "librocm-dbgapi"), Status, (ProcessID,), process_id)
end

function process_set_progress(process_id::ProcessID, progress::Progress)
    ccall((:amd_dbgapi_process_set_progress, "librocm-dbgapi"), Status, (ProcessID, Progress), process_id, progress)
end

function code_object_get_info(process_id::ProcessID, code_object_id::CodeObjectID,
                              query::CodeObjectInfo, value_size::Csize_t, value::Ref{Cvoid})
    ccall((:amd_dbgapi_code_object_get_info, "librocm-dbgapi"), Status,
          (ProcessID, CodeObjectID, CodeObjectInfo, Csize_t, Ref{Cvoid}),
          process_id, code_object_id, query, value_size, value)
end

function code_object_list(process_id::ProcessID, code_object_count::Ref{Csize_t},
                          code_objects::Ref{Ptr{CodeObjectID}}, changed::Ref{Changed})
    ccall((:amd_dbgapi_code_object_list, "librocm-dbgapi"), Status,
          (ProcessID, Ref{Csize_t}, Ref{Ptr{CodeObjectID}}, Ref{Changed}),
          process_id, code_object_count, code_objects, changed)
end

function agent_get_info(process_id::ProcessID, agent_id::AgentID, query::AgentInfo,
                        value_size::Csize_t, value::Ref{Cvoid})
    ccall((:amd_dbgapi_agent_get_info, "librocm-dbgapi"), Status,
          (ProcessID, AgentID, AgentInfo, Csize_t, Ref{Cvoid}),
          process_id, agent_id, query, value_size, value)
end

function agent_list(process_id::ProcessID, agent_count::Ref{Csize_t},
                    agents::Ref{Ptr{AgentID}}, changed::Ref{Changed})
    ccall((:amd_dbgapi_agent_list, "librocm-dbgapi"), Status,
          (ProcessID, Ref{Csize_t}, Ref{Ptr{AgentID}}, Ref{Changed}),
          process_id, agent_count, agents, changed)
end

function queue_get_info(process_id::ProcessID, queue_id::QueueID, query::QueueInfo,
                        value_size::Csize_t, value::Ref{Cvoid})
    ccall((:amd_dbgapi_queue_get_info, "librocm-dbgapi"), Status,
          (ProcessID, QueueID, QueueInfo, Csize_t, Ref{Cvoid}),
          process_id, queue_id, query, value_size, value)
end

function queue_list(process_id::ProcessID, queue_count::Ref{Csize_t},
                    queues::Ref{Ptr{QueueID}}, changed::Ref{Changed})

    ccall((:amd_dbgapi_queue_list, "librocm-dbgapi"), Status,
          (ProcessID, Ref{Csize_t}, Ref{Ptr{QueueID}}, Ref{Changed}),
          process_id, queue_count, queues, changed)
end

function queue_packet_list(process_id::ProcessID, queue_id::QueueID,
                           first_packet_id::Ref{QueuePacketID},
                           packets_byte_size::Ref{UInt64},
                           packets_bytes::Ref{Ptr{Cvoid}})

    ccall((:amd_dbgapi_queue_packet_list, "librocm-dbgapi"), Status,
          (ProcessID, QueueID, Ref{QueuePacketID}, Ref{UInt64}, Ref{Ptr{Cvoid}}),
          process_id, queue_id, first_packet_id, packets_byte_size, packets_bytes)
end

function dispatch_get_info(process_id::ProcessID, dispatch_id::DispatchID,
                           query::DispatchInfo, value_size::Csize_t, value::Ref{Cvoid})

    ccall((:amd_dbgapi_dispatch_get_info, "librocm-dbgapi"), Status,
          (ProcessID, DispatchID, DispatchInfo, Csize_t, Ref{Cvoid}),
          process_id, dispatch_id, query, value_size, value)
end

function dispatch_list(process_id::ProcessID, dispatch_count::Ref{Csize_t},
                       dispatches::Ref{Ptr{DispatchID}}, changed::Ref{Changed})

    ccall((:amd_dbgapi_dispatch_list, "librocm-dbgapi"), Status,
          (ProcessID, Ref{Csize_t}, Ref{Ptr{DispatchID}}, Ref{Changed}),
          process_id, dispatch_count, dispatches, changed)
end

function wave_get_info(process_id::ProcessID, wave_id::WaveID, query::WaveInfo,
                       value_size::Csize_t, value::Ref{Cvoid})
    ccall((:amd_dbgapi_wave_get_info, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, WaveInfo, Csize_t, Ref{Cvoid}),
          process_id, wave_id, query, value_size, value)
end

function wave_list(process_id::ProcessID, wave_count::Ref{Csize_t}, waves::Ref{Ptr{WaveID}},
                   changed::Ref{Changed})

    ccall((:amd_dbgapi_wave_list, "librocm-dbgapi"), Status,
          (ProcessID, Ref{Csize_t}, Ref{Ptr{WaveID}}, Ref{Changed}),
          process_id, wave_count, waves, changed)
end

function wave_stop(process_id::ProcessID, wave_id::WaveID)
    ccall((:amd_dbgapi_wave_stop, "librocm-dbgapi"), Status,
          (ProcessID, WaveID), process_id, wave_id)
end

function wave_resume(process_id::ProcessID, wave_id::WaveID, resume_mode::ResumeMode)
    ccall((:amd_dbgapi_wave_resume, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, ResumeMode), process_id, wave_id, resume_mode)
end

function displaced_stepping_start(process_id::ProcessID, wave_id::WaveID,
                                  saved_instruction_bytes::Ref{Cvoid},
                                  displaced_stepping::Ref{DisplacedSteppingID})

    ccall((:amd_dbgapi_displaced_stepping_start, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, Ref{Cvoid}, Ref{DisplacedSteppingID}),
          process_id, wave_id, saved_instruction_bytes, displaced_stepping)
end

function displaced_stepping_complete(process_id::ProcessID, wave_id::WaveID,
                                     displaced_stepping::DisplacedSteppingID)

    ccall((:amd_dbgapi_displaced_stepping_complete, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, DisplacedSteppingID), process_id, wave_id, displaced_stepping)
end

function set_watchpoint(process_id::ProcessID, agent_id::AgentID, address::GlobalAddress,
                        size::UInt64, kind::WatchpointKind,
                        watchpoint_id::Ref{WatchpointID},
                        watchpoint_address::Ref{GlobalAddress},
                        watchpoint_size::Ref{UInt64})

    ccall((:amd_dbgapi_set_watchpoint, "librocm-dbgapi"), Status,
          (ProcessID, AgentID, GlobalAddress, UInt64, WatchpointKind, Ref{WatchpointID},
           Ref{GlobalAddress}, Ref{UInt64}),
          process_id, agent_id, address, size, kind, watchpoint_id, watchpoint_address,
          watchpoint_size)
end

function remove_watchpoint(process_id::ProcessID, agent_id::AgentID,
                           watchpoint_id::WatchpointID)

    ccall((:amd_dbgapi_remove_watchpoint, "librocm-dbgapi"), Status,
          (ProcessID, AgentID, WatchpointID),
          process_id, agent_id, watchpoint_id)
end

function architecture_register_class_get_info(architecture_id::ArchitectureID,
                                              register_class_id::RegisterClassID,
                                              query::RegisterClassInfo,
                                              value_size::Csize_t, value::Ref{Cvoid})

    ccall((:amd_dbgapi_architecture_register_class_get_info, "librocm-dbgapi"), Status,
          (ArchitectureID, RegisterClassID, RegisterClassInfo, Csize_t, Ref{Cvoid}),
          architecture_id, register_class_id, query, value_size, value)
end

function architecture_register_class_list(architecture_id::ArchitectureID,
                                          register_class_count::Ref{Csize_t},
                                          register_classes::Ref{Ptr{RegisterClassID}})

    ccall((:amd_dbgapi_architecture_register_class_list, "librocm-dbgapi"), Status,
          (ArchitectureID, Ref{Csize_t}, Ref{Ptr{RegisterClassID}}),
          architecture_id, register_class_count, register_classes)
end

function architecture_register_get_info(architecture_id::ArchitectureID,
                                        register_id::RegisterID, query::RegisterInfo,
                                        value_size::Csize_t, value::Ref{Cvoid})

    ccall((:amd_dbgapi_architecture_register_get_info, "librocm-dbgapi"), Status,
          (ArchitectureID, RegisterID, RegisterInfo, Csize_t, Ref{Cvoid}),
          architecture_id, register_id, query, value_size, value)
end

function wave_register_get_info(process_id::ProcessID, wave_id::WaveID,
                                register_id::RegisterID, query::RegisterInfo,
                                value_size::Csize_t, value::Ref{Cvoid})

    ccall((:amd_dbgapi_wave_register_get_info, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, RegisterID, RegisterInfo, Csize_t, Ref{Cvoid}),
          process_id, wave_id, register_id, query, value_size, value)
end

function architecture_register_list(architecture_id::ArchitectureID,
                                    register_count::Ref{Csize_t},
                                    registers::Ref{Ptr{RegisterID}})

    ccall((:amd_dbgapi_architecture_register_list, "librocm-dbgapi"), Status,
          (ArchitectureID, Ref{Csize_t}, Ref{Ptr{RegisterID}}),
          architecture_id, register_count, registers)
end

function wave_register_list(process_id::ProcessID, wave_id::WaveID,
                            register_count::Ref{Csize_t},
                            registers::Ref{Ptr{RegisterID}})

    ccall((:amd_dbgapi_wave_register_list, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, Ref{Csize_t}, Ref{Ptr{RegisterID}}),
          process_id, wave_id, register_count, registers)
end

function dwarf_register_to_register(architecture_id::ArchitectureID, dwarf_register::UInt64,
                                    register_id::Ref{RegisterID})

    ccall((:amd_dbgapi_dwarf_register_to_register, "librocm-dbgapi"), Status,
          (ArchitectureID, UInt64, Ref{RegisterID}),
          architecture_id, dwarf_register, register_id)
end

function register_is_in_register_class(architecture_id::ArchitectureID,
                                       register_id::RegisterID,
                                       register_class_id::RegisterClassID,
                                       register_class_state::Ref{RegisterClassState})

    ccall((:amd_dbgapi_register_is_in_register_class, "librocm-dbgapi"), Status,
          (ArchitectureID, RegisterID, RegisterClassID, Ref{RegisterClassState}),
          architecture_id, register_id, register_class_id, register_class_state)
end

function read_register(process_id::ProcessID, wave_id::WaveID, register_id::RegisterID,
                       offset::UInt64, value_size::UInt64, value::Ref{Cvoid})

    ccall((:amd_dbgapi_read_register, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, RegisterID, UInt64, UInt64, Ref{Cvoid}),
          process_id, wave_id, register_id, offset, value_size, value)
end

function write_register(process_id::ProcessID, wave_id::WaveID, register_id::RegisterID,
                        offset::UInt64, value_size::UInt64, value::Ref{Cvoid})

    ccall((:amd_dbgapi_write_register, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, RegisterID, UInt64, UInt64, Ref{Cvoid}),
          process_id, wave_id, register_id, offset, value_size, value)
end

function prefetch_register(process_id::ProcessID, wave_id::WaveID, register_id::RegisterID,
                           register_count::UInt64)

    ccall((:amd_dbgapi_prefetch_register, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, RegisterID, UInt64),
          process_id, wave_id, register_id, register_count)
end

function architecture_address_class_get_info(architecture_id::ArchitectureID,
                                             address_class_id::AddressClassID,
                                             query::AddressClassInfo, value_size::Csize_t,
                                             value::Ref{Cvoid})

    ccall((:amd_dbgapi_architecture_address_class_get_info, "librocm-dbgapi"), Status,
          (ArchitectureID, AddressClassID, AddressClassInfo, Csize_t, Ref{Cvoid}),
          architecture_id, address_class_id, query, value_size, value)
end

function architecture_address_class_list(architecture_id::ArchitectureID,
                                         address_class_count::Ref{Csize_t},
                                         address_classes::Ref{Ptr{AddressClassID}})

    ccall((:amd_dbgapi_architecture_address_class_list, "librocm-dbgapi"), Status,
          (ArchitectureID, Ref{Csize_t}, Ref{Ptr{AddressClassID}}),
          architecture_id, address_class_count, address_classes)
end

function dwarf_address_class_to_address_class(architecture_id::ArchitectureID,
                                              dwarf_address_class::UInt64,
                                              address_class_id::Ref{AddressClassID})

    ccall((:amd_dbgapi_dwarf_address_class_to_address_class, "librocm-dbgapi"), Status,
          (ArchitectureID, UInt64, Ref{AddressClassID}),
          architecture_id, dwarf_address_class, address_class_id)
end

function address_space_get_info(architecture_id::ArchitectureID,
                                address_space_id::AddressSpaceID,
                                query::AddressSpaceInfo, value_size::Csize_t,
                                value::Ref{Cvoid})

    ccall((:amd_dbgapi_address_space_get_info, "librocm-dbgapi"), Status,
          (ArchitectureID, AddressSpaceID, AddressSpaceInfo, Csize_t, Ref{Cvoid}),
          architecture_id, address_space_id, query, value_size, value)
end

function architecture_address_space_list(architecture_id::ArchitectureID,
                                         address_space_count::Ref{Csize_t},
                                         address_spaces::Ref{Ptr{AddressSpaceID}})

    ccall((:amd_dbgapi_architecture_address_space_list, "librocm-dbgapi"), Status,
          (ArchitectureID, Ref{Csize_t}, Ref{Ptr{AddressSpaceID}}),
          architecture_id, address_space_count, address_spaces)
end

function dwarf_address_space_to_address_space(architecture_id::ArchitectureID,
                                              dwarf_address_space::UInt64,
                                              address_space_id::Ref{AddressSpaceID})

    ccall((:amd_dbgapi_dwarf_address_space_to_address_space, "librocm-dbgapi"), Status,
          (ArchitectureID, UInt64, Ref{AddressSpaceID}),
          architecture_id, dwarf_address_space, address_space_id)
end

function address_spaces_may_alias(architecture_id::ArchitectureID,
                                  address_space_id1::AddressSpaceID,
                                  address_space_id2::AddressSpaceID,
                                  address_space_alias::Ref{AddressSpaceAlias})

    ccall((:amd_dbgapi_address_spaces_may_alias, "librocm-dbgapi"), Status,
          (ArchitectureID, AddressSpaceID, AddressSpaceID, Ref{AddressSpaceAlias}),
          architecture_id, address_space_id1, address_space_id2, address_space_alias)
end

function convert_address_space(process_id::ProcessID, wave_id::WaveID, lane_id::LaneID,
                               source_address_space_id::AddressSpaceID,
                               source_segment_address::SegmentAddress,
                               destination_address_space_id::AddressSpaceID,
                               destination_segment_address::Ref{SegmentAddress})

    ccall((:amd_dbgapi_convert_address_space, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, LaneID, AddressSpaceID, SegmentAddress, AddressSpaceID,
           Ref{SegmentAddress}),
          process_id, wave_id, lane_id, source_address_space_id, source_segment_address,
          destination_address_space_id, destination_segment_address)
end

function address_is_in_address_class(process_id::ProcessID, wave_id::WaveID,
                                     lane_id::LaneID, address_space_id::AddressSpaceID,
                                     segment_address::SegmentAddress,
                                     address_class_id::AddressClassID,
                                     address_class_state::Ref{AddressClassState})

    ccall((:amd_dbgapi_address_is_in_address_class, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, LaneID, AddressSpaceID, SegmentAddress, AddressClassID,
           Ref{AddressClassState}),
          process_id, wave_id, lane_id, address_space_id, segment_address, address_class_id,
          address_class_state)
end

function read_memory(process_id::ProcessID, wave_id::WaveID, lane_id::LaneID,
                     address_space_id::AddressSpaceID, segment_address::SegmentAddress,
                     value_size::Ref{UInt64}, value::Ref{Cvoid})

    ccall((:amd_dbgapi_read_memory, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, LaneID, AddressSpaceID, SegmentAddress, Ref{UInt64},
           Ref{Cvoid}),
          process_id, wave_id, lane_id, address_space_id, segment_address, value_size,
          value)
end

function write_memory(process_id::ProcessID, wave_id::WaveID, lane_id::LaneID,
                      address_space_id::AddressSpaceID, segment_address::SegmentAddress,
                      value_size::Ref{UInt64}, value::Ref{Cvoid})

    ccall((:amd_dbgapi_write_memory, "librocm-dbgapi"), Status,
          (ProcessID, WaveID, LaneID, AddressSpaceID, SegmentAddress, Ref{UInt64},
           Ref{Cvoid}),
          process_id, wave_id, lane_id, address_space_id, segment_address, value_size,
          value)
end

function set_memory_precision(process_id::ProcessID, agent_id::AgentID,
                              memory_precision::MemoryPrecision)

    ccall((:amd_dbgapi_set_memory_precision, "librocm-dbgapi"), Status,
          (ProcessID, AgentID, MemoryPrecision), process_id, agent_id, memory_precision)
end

function next_pending_event(process_id::ProcessID, event_id::Ref{EventID},
                            kind::Ref{EventKind})

    ccall((:amd_dbgapi_next_pending_event, "librocm-dbgapi"), Status,
          (ProcessID, Ref{EventID}, Ref{EventKind}), process_id, event_id, kind)
end

function event_get_info(process_id::ProcessID, event_id::EventID, query::EventInfo,
                        value_size::Csize_t, value::Ref{Cvoid})

    ccall((:amd_dbgapi_event_get_info, "librocm-dbgapi"), Status,
          (ProcessID, EventID, EventInfo, Csize_t, Ref{Cvoid}),
          process_id, event_id, query, value_size, value)
end

function event_processed(process_id::ProcessID, event_id::EventID)
    ccall((:amd_dbgapi_event_processed, "librocm-dbgapi"), Status,
          (ProcessID, EventID),
          process_id, event_id)
end

function set_log_level(level::LogLevel)
    ccall((:amd_dbgapi_set_log_level, "librocm-dbgapi"), Cvoid,
          (LogLevel,), level)
end

function report_shared_library(process_id::ProcessID, shared_library_id::SharedLibraryID,
                               shared_library_state::SharedLibraryState)

    ccall((:amd_dbgapi_report_shared_library, "librocm-dbgapi"), Status,
          (ProcessID, SharedLibraryID, SharedLibraryState),
          process_id, shared_library_id, shared_library_state)
end

function report_breakpoint_hit(process_id::ProcessID, breakpoint_id::BreakpointID,
                               client_thread_id::ClientThreadID,
                               breakpoint_action::Ref{BreakpointAction})

    ccall((:amd_dbgapi_report_breakpoint_hit, "librocm-dbgapi"), Status,
          (ProcessID, BreakpointID, ClientThreadID, Ref{BreakpointAction}),
          process_id, breakpoint_id, client_thread_id, breakpoint_action)
end


text:
using CSV
using DataFrames
using Interpolations
using PyCall
using SiriusB
using Statistics

pro = pyimport("proplot")
# rcParams
pro.use_style("ggplot")
pro.rc["image.origin"] = "lower"
pro.rc["image.cmap"] = "inferno"
pro.rc["grid"] = false

contrast_curves = [
    DataFrame(CSV.File(datadir("epoch_2020feb04", "residuals", "2020feb04_contrast-curve_median.csv"))),
    DataFrame(CSV.File(datadir("epoch_2020nov21", "residuals", "2020nov21_contrast-curve_median.csv"))),
    DataFrame(CSV.File(datadir("epoch_2020nov28", "residuals", "2020nov28_contrast-curve_annular-pca-2.csv"))),
]

# convert distances from pixels to AU
distances = [curve.distance .* SiriusB.auscale for curve in contrast_curves]

function contrast_to_absmag(y)
    Δmag = contrast_to_dmag(y)
    appmag = SiriusB.appmag + Δmag
    return appmag - distance_modulus(SiriusB.distance)
end

## plot each epoch as contrast/absolute mag
py"""
import proplot as pro
fig, axs = pro.subplots(width="7.5in", height="4in")

# curves
axs.plot($(distances[1]), $(contrast_curves[1].contrast), c="C0", label="2020-02-04")
axs.plot($(distances[1]), $(contrast_curves[1].contrast_corr), c="C0", ls="--")

axs.plot($(distances[2]), $(contrast_curves[2].contrast), c="C1", label="2020-11-21")
axs.plot($(distances[2]), $(contrast_curves[2].contrast_corr), c="C1", ls="--")

axs.plot($(distances[3]), $(contrast_curves[3].contrast), c="C5", label="2020-11-28")
axs.plot($(distances[3]), $(contrast_curves[3].contrast_corr), c="C5", ls="--")

# formatting
axs.dualx(lambda x: x * $(SiriusB.parallax), label="separation [arcsec]")
axs.format(
    yscale="log",
    xlabel="projected separation [AU]",
    ylabel="contrast ",
    grid=True,
    yformatter="log",
    xlim=(None, 2),
)

# he's mad, he's gona and made a custom legend!
from matplotlib.lines import Line2D
leg1 = axs.legend(ncol=1, loc="ur")
elements = [
    Line2D([0], [0], color="k", alpha=0.6, label="Gaussian"),
    Line2D([0], [0], color="k", alpha=0.6, ls="--", label="Student-t"),
]
axs.legend(handles=elements, queue=True, loc="uc")
axs.add_artist(leg1)

# stability limit
ylims=axs.get_ylim()
axs.vlines(1.5, *ylims, color="k", alpha=0.4, ls="-.")
mid = (ylims[0] + ylims[1]) / 40
# print(mid)
axs.text(1.45, mid, "dynamical stability limit", color="k", alpha=0.5, va="center", ha="left", rotation="vertical")
axs.set_ylim(ylims)

ax2 = axs.alty(reverse=True, label="abs. mag [$\mathrm{L^\prime}$]")
ax2.plot($(distances[2]), $(contrast_to_absmag.(contrast_curves[2].contrast)), alpha=0)

fig.savefig($(figuredir("contrast_curves.pdf")))
pro.close(fig)
"""

for curve in contrast_curves
    curve.absmags = contrast_to_absmag.(curve.contrast)
    curve.masses_225 = map(curve.absmags) do absmag
        table_interpolate(ATMO2020, :Age => 0.225, :MKO_Lp => absmag, :Mass_MJ)
    end
    curve.masses_126 = map(curve.absmags) do absmag
        table_interpolate(ATMO2020, :Age => 0.126, :MKO_Lp => absmag, :Mass_MJ)
    end
    curve.absmags_corr = contrast_to_absmag.(curve.contrast_corr)
    curve.masses_corr_225 = map(curve.absmags_corr) do absmag
        table_interpolate(ATMO2020, :Age => 0.225, :MKO_Lp => absmag, :Mass_MJ)
    end
    curve.masses_corr_126 = map(curve.absmags_corr) do absmag
        table_interpolate(ATMO2020, :Age => 0.126, :MKO_Lp => absmag, :Mass_MJ)
    end
    # these try catch curves are because Sonora grid bottoms out
    curve.masses_225_SonoraSolar = map(curve.absmags) do absmag
        table_interpolate(SonoraSolar, :age_Gyr => 0.225, :Keck_Lp => absmag, :M_Mjup)
    end
    curve.masses_126_SonoraSolar = map(curve.absmags) do absmag
        table_interpolate(SonoraSolar, :age_Gyr => 0.126, :Keck_Lp => absmag, :M_Mjup)
    end
    curve.masses_corr_225_SonoraSolar = map(curve.absmags_corr) do absmag
        table_interpolate(SonoraSolar, :age_Gyr => 0.225, :Keck_Lp => absmag, :M_Mjup)
    end
    curve.masses_corr_126_SonoraSolar = map(curve.absmags_corr) do absmag
        table_interpolate(SonoraSolar, :age_Gyr => 0.126, :Keck_Lp => absmag, :M_Mjup)
    end
    # metalrich
    curve.masses_225_SonoraMetalRich = map(curve.absmags) do absmag
        table_interpolate(SonoraMetalRich, :age_Gyr => 0.225, :Keck_Lp => absmag, :M_Mjup)
    end
    curve.masses_126_SonoraMetalRich = map(curve.absmags) do absmag
        table_interpolate(SonoraMetalRich, :age_Gyr => 0.126, :Keck_Lp => absmag, :M_Mjup)
    end
    curve.masses_corr_225_SonoraMetalRich = map(curve.absmags_corr) do absmag
        table_interpolate(SonoraMetalRich, :age_Gyr => 0.225, :Keck_Lp => absmag, :M_Mjup)
    end
    curve.masses_corr_126_SonoraMetalRich = map(curve.absmags_corr) do absmag
        table_interpolate(SonoraMetalRich, :age_Gyr => 0.126, :Keck_Lp => absmag, :M_Mjup)
    end
end

## plot best epoch as mass limits
py"""
import proplot as pro
fig, axs = pro.subplots(ncols=2, share=2, width="7.5in", height="4in")

# atmo2020 curves
axs[0].plot($(distances[2]), $(contrast_curves[2].masses_225), c="C2", label="225 Myr")
axs[0].plot($(distances[2]), $(contrast_curves[2].masses_corr_225), c="C2", ls="--")
axs[0].fill_between($(distances[2]), $(contrast_curves[2].masses_225), $(contrast_curves[2].masses_corr_225), color="C2", lw=0, alpha=0.1)

axs[0].plot($(distances[2]), $(contrast_curves[2].masses_126), c="C3", label="126 Myr")
axs[0].plot($(distances[2]), $(contrast_curves[2].masses_corr_126), c="C3", ls="--")
axs[0].fill_between($(distances[2]), $(contrast_curves[2].masses_126), $(contrast_curves[2].masses_corr_126), color="C3", lw=0, alpha=0.1)

# curves
axs[1].plot($(distances[2]), $(contrast_curves[2].masses_225_SonoraSolar), c="C4", label="225 Myr")
axs[1].plot($(distances[2]), $(contrast_curves[2].masses_126_SonoraSolar), c="C5", label="126 Myr")
axs[1].plot($(distances[2]), $(contrast_curves[2].masses_225_SonoraMetalRich), c="C4", ls=":")
axs[1].plot($(distances[2]), $(contrast_curves[2].masses_126_SonoraMetalRich), c="C5", ls=":")

# # model mass limit
[ax.hlines([1, 3, 5, 6, 7, 9, 10], *ax.get_xlim(), color="w", alpha=0.5, lw=0.75) for ax in axs]

# formatting
axs.format(
    abc=True,
    yscale=pro.LogScale(base=2),
    xlabel="projected separation [AU]",
    ylabel="companion mass [$M_J$] ",
    grid=True,
    yformatter="auto",
    xlim=(None, 1.5),
)
axs[0].format(title="ATMO2020")
axs[1].format(title="Sonora Bobcat")

# custom legend
axs.legend(ncol=1, loc="ur")
elements = [
    Line2D([0], [0], color="k", alpha=0.6, label="Gaussian"),
    Line2D([0], [0], color="k", alpha=0.6, ls="--", label="Student-t"),
]
axs[0].legend(handles=elements, queue=True, loc="t")
elements = [
    Line2D([0], [0], color="k", alpha=0.6, label="[M/H] = 0.0"),
    Line2D([0], [0], color="k", alpha=0.6, ls=":", label="[M/H] = +0.5"),
]
axs[1].legend(handles=elements, queue=True, loc="t")

fig.savefig($(figuredir("mass_curves.pdf")))
pro.close(fig)
"""

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
using IncDataStats
using Base.Test

grb = GroupBy([:a, :b, :c, :d], [:a, :c])
update!(grb, [1, 2, 3, 4])
update!(grb, [1, 3.2, 3, 7])
update!(grb, [2, 2, 1, 4])
update!(grb, [2, 2, 7, 4])

@test mean(grb, :a) == 1.5
@test mean(grb, :b) == 2.3
@test sum(grb, :c) == 14
@test count(grb, :d) == 4

bins, h = hist(grb, :c)
@test bins == [1.0,3.0,7.0]
@test h == [1,2,1]

@test quantile(grb, :c, .25) == 1
@test quantile(grb, :c, .5) == 3
@test quantile(grb, :c, .75) == 3
@test quantile(grb, :c, .9) == 7

# groupers test

grb = GroupBy([:a, :b, :c, :d], [:a => identity, :c => stepgrouper(3)])
update!(grb, [1, 2, 3, 4])
update!(grb, [1, 3.2, 3, 7])
update!(grb, [2, 2, 1, 4])
update!(grb, [2, 2, 7, 4])


text:
timet = LinRange(0,1,101);
M = length(timet);
N = 30;
lam = 0.001;
center = [.35; .5; .65];
center2 = [4; 3.7; 4];
sd1 = .05;
gam_sd = 2;
num_comp = 5;
f_orig = zeros(M,N*length(center));
omega = 2*pi;
cnt = 1;
for ii = 1:length(center)
    global cnt
    tmp = Normal(center[ii], .075);
    tmp2 = Normal(center2[ii], sd1);
    for jj = 1:N
        f_orig[:, cnt] = rand(tmp2) * pdf.(tmp, timet);
        cnt += 1;
    end
end

q_orig = f_to_srsf(f_orig,collect(timet));
y_orig = zeros(size(q_orig,2));
alpha_t = 0;
b1 = sin.(omega*timet);
b2 = cos.(omega*timet);
B = hcat(b1,b2);
bt = .5 * b1 + 0.9 * b2;
cnt = 1;
for ii = 1:length(center)
    global cnt
    tmp = Normal(center[ii], .075);
    d = Normal(center2[ii], sd1);
    d1 = Normal(0, 0.01);
    tmp2 = f_to_srsf(rand(d) * pdf.(tmp, timet), collect(timet));
    for jj = 1:N
        y_orig[cnt] = alpha_t + trapz(collect(timet), tmp2 .* bt) + rand(d);
        cnt += 1;
    end
end

f = zeros(M, size(f_orig,2));
q = zeros(M, size(f_orig,2));
cnt = 1;
for ii = 1:length(center)
    global cnt
    gam_orig = rgam(M, gam_sd, N);
    for jj = 1:N
        f[:,cnt] = warp_f_gamma(collect(timet), f_orig[:,cnt], gam_orig[:,jj]);
        q[:,cnt] = f_to_srsf(f[:,cnt], collect(timet));
        cnt += 1;
    end
end


text:
# TODO: add impulse response iterators

# features of the input output arrays to plan(w::𝕎)
# ====================================

# features of the grid from w.period
# =================================

function Δpix(w::𝕎{Tf}) where {Tf}
    rTf = real(Tf)
    rTf.(w.period ./ w.sz)
end

function Δfreq(w::𝕎{Tf}) where {Tf}
    rTf = real(Tf)
    rTf.((2π) ./ w.period)
end

function nyq(w::𝕎{Tf}) where {Tf}
    rTf = real(Tf)
    rTf.(π ./ Δpix(w))
end

# Note: this gives the area element of
# only the fourier tranformed coordinates
Ωpix(w::𝕎) = prod(Δr[1] for Δr in zip(Δpix(w), w.region) if Δr[2])

Ωfreq(w::𝕎) = prod(Δr[1] for Δr in zip(Δfreq(w), w.region) if Δr[2])

# 𝕎 scalings
# =================================

function unscale(w::𝕎{Tf,d}) where {Tf,d}
    𝕎{Tf,d}(w.sz, w.region, true, w.period)
end

function Base.real(w::𝕎{Tf,d}) where {Tf,d}
    𝕎{real(Tf),d}(w.sz, w.region, w.scale, w.period)
end

function Base.complex(w::𝕎{Tf,d}) where {Tf,d}
    𝕎{Complex{real(Tf)},d}(w.sz, w.region, w.scale, w.period)
end

"`nv_scale(w::𝕎)->Number` returns the multiplicative scale of the inverse of w::𝕎"
function inv_scale(w::𝕎{Tf}) where {Tf}
    rTf = real(Tf)
    ifft_normalization = FFTW.normalization(
                rTf,
                w.sz,
                tuple(findall(w.region)...)
            )
    return rTf(ifft_normalization / w.scale)
end

function unitary_scale(w::𝕎{Tf}) where {Tf}
    rTf = real(Tf)
    return rTf(prod(1/√i[1] for i in zip(w.sz, w.region) if i[2]))
end

function ordinary_scale(w::𝕎{Tf}) where {Tf}
    rTf = real(Tf)
    rTf(Ωpix(w) / ((2π)^(sum(w.region)/2)))
end

# Test that inv_scale(ordinary_scale(w))== Ωfreq(w) / ((2π)^(sum(w.region)/2))

# pix grid
# =================================

# low level
# -----------

function pix(n::Int, p::Real)
    Δx = p/n
    x  = (0:n-1) * Δx |> collect
    return x
end

# main API
# -----------

function pix(w::𝕎{Tf,d}) where {d,Tf}
    rTf = real(Tf)
    return map(w.sz, w.period) do nᵢ, pᵢ
        rTf.(pix(nᵢ, pᵢ))
    end::NTuple{d,Array{rTf,1}}
end

function fullpix(i::Int, w::𝕎{Tf,d,Tsf,Tp}) where {Tf,d,Tsf,Tp}
    xi     = pix(w)[i]
    rTf    = real(Tf)
    xifull = zeros(rTf, w.sz)
    for Ic ∈ CartesianIndices(xifull)
        xifull[Ic] = getindex(xi, Tuple(Ic)[i])
    end
    return xifull::Array{rTf,d}
end

function fullpix(w::𝕎{Tf,d,Tsi,Tp}) where {Tf,d,Tsi,Tp}
    rTf = real(Tf)
    return map(tuple(1:d...)) do i
        fullpix(i, w)
    end::NTuple{d,Array{rTf,d}}
end

# frequency grid
# =================================

# low level
# -----------

function freq(n::Int, p::Real)
    k    = _fft_output_index_2_freq.(1:n, n, p)
    return k
end

function rfreq(n::Int, p::Tp) where Tp
    freq(n, p)[1:(n÷2+1)]
end

# main API
# -----------

function freq(w::𝕎{Tf,d}) where {d, Tf<:FFTC}
    rTf = real(Tf)
    return map(w.sz, w.period, w.region) do nᵢ, pᵢ, rᵢ
        kᵢ = rᵢ ? freq(nᵢ, pᵢ) : pix(nᵢ, pᵢ)
        rTf.(kᵢ)
    end::NTuple{d,Array{rTf,1}}
end

function freq(w::𝕎{Tf,d}) where {d, Tf<:FFTR}
    ir = findfirst(w.region)
    return map(w.sz, w.period, w.region, tuple(1:d...)) do nᵢ, pᵢ, rᵢ, i
        kᵢ = (i==ir) ? rfreq(nᵢ, pᵢ) :
                  rᵢ ?  freq(nᵢ, pᵢ) : pix(nᵢ, pᵢ)
        Tf.(kᵢ)
    end::NTuple{d,Array{Tf,1}}
end

function fullfreq(i::Int, w::𝕎{Tf,d,Tsf,Tp}) where {Tf,d,Tsf,Tp}
    ki  = freq(w)[i]
    rTf = real(Tf)
    kifull = zeros(rTf, size_out(w))
    for Ic ∈ CartesianIndices(kifull)
        kifull[Ic] = getindex(ki, Tuple(Ic)[i])
    end
    return kifull::Array{rTf,d}
end

function fullfreq(w::𝕎{Tf,d,Tsf,Tp}) where {Tf,d,Tsf,Tp}
    return map(tuple(1:d...)) do i
        fullfreq(i,w)
    end::NTuple{d,Array{real(Tf),d}}
end

function wavenum(w::𝕎{Tf,d,Tsf,Tp}) where {Tf,d,Tsf,Tp}
    k   = freq(w)
    rTf = real(Tf)
    λ   = zeros(rTf, size_out(w))
    for Ic ∈ CartesianIndices(λ)
        λ[Ic] = sqrt(sum(abs2, getindex.(k,Tuple(Ic))))
    end
    return λ::Array{rTf,d}
end

# Internal function
# -----------------------------------
#`k_pre = _fft_output_index_2_freq.(1:n, n, p)` computes the 1-d freq
function _fft_output_index_2_freq(ind, nside, period)
    kpre = (2π / period) * (ind - 1)
    nyq  = (2π / period) * (nside/2)

    # • Both of the following options are valid
    # • Both options return the same value when nside is odd
    # • Using (kpre <= nyq) sets the convention  that the
    #   frequencies fall in (-nyq, nyq] when nside is even
    # • Using (kpre < nyq) sets the convention  that the
    #   frequencies fall in [-nyq, nyq) when nside is even
    #----------------
    # • Here are the two options:
    return  ifelse(kpre <= nyq, kpre, kpre - 2nyq) # option 1
    # return ifelse(kpre < nyq, kpre, kpre - 2nyq)  # option 2
end

# function _get_npd(;nᵢ, pᵢ=nothing, Δxᵢ=nothing)
#     @assert !(isnothing(pᵢ) & isnothing(Δxᵢ)) "either pᵢ or Δxᵢ needs to be specified (note: pᵢ = Δxᵢ .* nᵢ)"
#     d = length(nᵢ)
#     if isnothing(pᵢ)
#         @assert d == length(Δxᵢ) "Δxᵢ and nᵢ need to be tuples of the same length"
#         pᵢ = tuple((prod(xn) for xn in zip(Δxᵢ,nᵢ))...)
#     end
#     @assert d == length(pᵢ) "pᵢ and nᵢ need to be tuples of the same length"
#     nᵢ, pᵢ, d
# end

# #%% Used for constructing the covariance matrix of a subset of frequencies
# function get_rFFTimpulses(::Type{F}) where {T,nᵢ,pᵢ,dnᵢ,F<:rFFTgeneric{T,nᵢ,pᵢ,dnᵢ}}
#     g  = Grid(F)
#     CI = CartesianIndices(Base.OneTo.(g.nki))
#     LI = LinearIndices(Base.OneTo.(g.nki))

#     function _get_dual_k(k,n)
#         dk = n-k+2
#         mod1(dk,n)
#     end

#     function get_dual_ci(ci::CartesianIndex{dnᵢ})
#         return CartesianIndex(map(_get_dual_k, ci.I, g.nxi))
#     end

#     function rFFTimpulses(ci::CartesianIndex{dnᵢ})
#         rimpls = zeros(Complex{T}, g.nki...)
#         cimpls = zeros(Complex{T}, g.nki...)
#         dual_ci = get_dual_ci(ci)
#         if (ci==first(CI)) || (ci==dual_ci)
#             rimpls[ci]  = 1
#         elseif dual_ci ∈ CI
#             rimpls[ci]  = 1/2
#             cimpls[ci]  = im/2
#             rimpls[dual_ci]  =  1/2
#             cimpls[dual_ci]  = -im/2
#         else
#             rimpls[ci]  = 1/2
#             cimpls[ci]  = im/2
#         end
#         return rimpls, cimpls
#     end

#     return rFFTimpulses, CI, LI, get_dual_ci
# end


text:
using Documenter
using Algames  # your package name here

makedocs(
    sitename = "Algames",  # your package name here
    format = Documenter.HTML(prettyurls = false),  # optional
    pages = [
        "Introduction" => "index.md",
        "Basics" => [
            "Creating a New Package" => "new_repo.md",
            "Adding Dependencies" => "deps.md",
            "Adding Tests" => "tests.md",
            "Setting up CI" => "CI.md",
            "Adding Code Coverage" => "codecov.md",
            "Building documentation" => "docs.md",
            "Registering a Package" => "registering.md"
        ],
        "Sample Documentation" => "package_docs.md"
    ]
)

deploydocs(
    repo = "github.com/RoboticExplorationLab/Algames.jl.git",
)


text:
function fib(n)
    if n < 3
        return 1
    end

    return fib(n-1) + fib(n-2)
end

n = parse(Int64, ARGS[1])
println(fib(n))

max_stars_repo_name: juliancoffee/fib

max_issues_repo_name: juliancoffee/fib

max_forks_repo_name: juliancoffee/fib


text:
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges

struct ConstraintDummyBridge <: MOIB.Constraint.AbstractBridge
    id::Int
end

map = MOIB.Constraint.Map()
@testset "Empty" begin
    @test isempty(map)
    @test !MOIB.Constraint.has_bridges(map)
end
x = MOI.VariableIndex(1)
y = MOI.VariableIndex(2)
fx = MOI.SingleVariable(x)
b1 = ConstraintDummyBridge(1)
c1 = MOIB.Constraint.add_key_for_bridge(map, b1, fx, MOI.EqualTo(0.0))
@testset "SingleVariable" begin
    @test c1.value == x.value
    @test haskey(map, c1)
    @test map[c1] == b1
    @test MOIB.Constraint.number_of_type(map, typeof(c1)) == 1
    @test collect(MOIB.Constraint.keys_of_type(map, typeof(c1))) == [c1]

    @test length(map) == 1
    @test !isempty(map)
    @test MOIB.Constraint.has_bridges(map)
end

b2 = ConstraintDummyBridge(2)
vov = MOI.VectorOfVariables([x, y])
c2 = MOIB.Constraint.add_key_for_bridge(map, b2, vov, MOI.SecondOrderCone(2))
@testset "VectorOfVariables" begin
    @test c2.value == -1
    @test haskey(map, c2)
    @test map[c2] == b2
    @test MOIB.Constraint.number_of_type(map, typeof(c2)) == 1
    @test collect(MOIB.Constraint.keys_of_type(map, typeof(c2))) == [c2]

    @test length(map) == 2
    @test !isempty(map)
    @test MOIB.Constraint.has_bridges(map)
end

b3 = ConstraintDummyBridge(3)
c3 = MOIB.Constraint.add_key_for_bridge(map, b3, 1.0fx + 2.0, MOI.EqualTo(0.0))
@testset "ScalarAffineFunction" begin
    @test haskey(map, c3)
    @test map[c3] == b3
    @test MOIB.Constraint.number_of_type(map, typeof(c3)) == 1
    @test collect(MOIB.Constraint.keys_of_type(map, typeof(c3))) == [c3]

    @test length(map) == 3
    @test !isempty(map)
    @test MOIB.Constraint.has_bridges(map)

    bridges = collect(values(map))
    @test sort([b.id for b in bridges]) == 1:3
    elements = sort(collect(map), by = el -> el.second.id)
    @test elements[1].first == c1
    @test elements[1].second == b1
    @test elements[2].first == c2
    @test elements[2].second == b2
    @test elements[3].first == c3
    @test elements[3].second == b3

    @test MOIB.Constraint.list_of_key_types(map) == Set([
        (MOI.VectorOfVariables, MOI.SecondOrderCone),
        (MOI.SingleVariable, MOI.EqualTo{Float64}),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
    ])
    @test MOIB.Constraint.variable_constraints(map, x) == [c1]
    @test isempty(MOIB.Constraint.variable_constraints(map, y))
end

@testset "Delete" begin
    delete!(map, c1)
    @test length(map) == 2
    @test sprint(MOIB.print_num_bridges, map) == "\nwith 2 constraint bridges"
    @test !isempty(map)
    @test MOIB.Constraint.has_bridges(map)

    @test MOIB.Constraint.list_of_key_types(map) == Set([
        (MOI.VectorOfVariables, MOI.SecondOrderCone),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
    ])
    @test isempty(MOIB.Constraint.variable_constraints(map, x))
    @test isempty(MOIB.Constraint.variable_constraints(map, y))
    @test MOIB.Constraint.number_of_type(map, typeof(c1)) == 0
    @test MOIB.Constraint.number_of_type(map, typeof(c2)) == 1
    @test MOIB.Constraint.number_of_type(map, typeof(c3)) == 1

    delete!(map, c2)
    @test length(map) == 1
    @test sprint(MOIB.print_num_bridges, map) == "\nwith 1 constraint bridge"
    @test !isempty(map)
    @test MOIB.Constraint.has_bridges(map)

    @test MOIB.Constraint.list_of_key_types(map) ==
          Set([(MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})])
    @test MOIB.Constraint.number_of_type(map, typeof(c1)) == 0
    @test MOIB.Constraint.number_of_type(map, typeof(c2)) == 0
    @test MOIB.Constraint.number_of_type(map, typeof(c3)) == 1

    delete!(map, c3)
    @test length(map) == 0
    @test isempty(map)
    @test MOIB.Constraint.has_bridges(map)

    @test isempty(MOIB.Constraint.list_of_key_types(map))
    @test MOIB.Constraint.number_of_type(map, typeof(c1)) == 0
    @test MOIB.Constraint.number_of_type(map, typeof(c2)) == 0
    @test MOIB.Constraint.number_of_type(map, typeof(c3)) == 0

    @test isempty(collect(map))
end

@testset "EmptyMap" begin
    map = MOIB.Constraint.EmptyMap()
    empty!(map)
    @test isempty(map)
    @test isempty(keys(map))
    @test isempty(values(map))
    @test !MOIB.Constraint.has_bridges(map)
    @test iszero(MOIB.Constraint.number_of_type(map, typeof(c1)))
    @test iszero(MOIB.Constraint.number_of_type(map, typeof(c2)))
    @test iszero(MOIB.Constraint.number_of_type(map, typeof(c3)))
    @test sprint(MOIB.print_num_bridges, map) == ""
end


text:
function generate_Cn(n)
    symels = []
    axis = [0 0 1]'
    cn_r = Cn(axis, n)
    for i = 1:n-1
        a, b = reduce(n, i)
        push!(symels, Symel("C_$a^$b", cn_r^i)) # Cns
    end
    return symels
end

function generate_Sn(n)
    generate_Sn(n, false)
end

function generate_Sn(n, S2n)
    symels = []
    axis = [0 0 1]'
    σh = σ(axis)
    cn_r = Cn(axis, n)
    if S2n # Generating improper rotations for S2n PG
        for i = 1:n-1
            if i % 2 == 0
                continue
            else
                a, b = reduce(n, i)
                if a == 2
                    continue
                else
                    push!(symels, Symel("S_$a^$b", cn_r^i * σh))
                end
            end
        end
        return symels
    end
    for i = 1:n-1
        a, b = reduce(n, i)
        if b % 2 == 0
            b += n
        end
        if a == 2
            continue
        else
            push!(symels, Symel("S_$a^$b", cn_r^i * σh)) # Sns
        end
    end
    return symels
end

function generate_σv(n)
    symels = []
    x_axis = [1 0 0]' # Orient C2 and σv along x-axis
    rot_mat = Cn([0 0 1]', n)
    for i = 1:n
        axis = (rot_mat ^ i) * x_axis
        push!(symels, Symel("sigmav_$(i % n + 1)", σ(axis)))
    end
    return symels
end

function generate_σd(n)
    symels = []
    x_axis = [1 0 0]' # Orient C2 and σv along x-axis
    z_axis = [0 0 1]'
    rot_mat = Cn(z_axis, 2*n)
    base_axis = Cn(z_axis, 4*n)*x_axis # Rotate x-axis by Cn/2 to produce an axis for σd's
    for i = 1:n
        axis = (rot_mat ^ i) * base_axis
        push!(symels, Symel("sigmad_$(i % n + 1)", σ(axis))) # Take moduli so that 360° rotated axis is 0
    end
    return symels
end

function generate_C2(n)
    if n % 2 == 0
        nn = 2*n
    else
        nn = n
    end
    symels = []
    x_axis = [1 0 0]' # Orient C2 and σv along x-axis
    rot_mat = Cn([0 0 1]', n)
    for i = 1:n
        axis = (rot_mat ^ i) * x_axis
        push!(symels, Symel("C_2($(i % n + 1))", Cn(axis, 2)))
    end
    return symels
end

function generate_T()
    """
        Assume a tetrahedron contained in a cube, then we can easily generate
        the vectors for the rotation elements.
    """
    # Generate C3's
    symels = []
    C3_1v = normalize!([1.0 1.0 1.0])
    C3_2v = normalize!([-1.0 1.0 -1.0])
    C3_3v = normalize!([-1.0 -1.0 1.0])
    C3_4v = normalize!([1.0 -1.0 -1.0])
    C3list = (C3_1v, C3_2v, C3_3v, C3_4v)
    namelist = ("α", "β", "γ", "δ")
    for i = 1:4
        C3 = Cn(C3list[i], 3)
        C3s = C3^2
        push!(symels, Symel("C_3($(namelist[i]))",C3))
        push!(symels, Symel("C_3^2($(namelist[i]))",C3s))
    end
    # Generate C2's
    C2_x = [1.0 0.0 0.0]
    C2_y = [0.0 1.0 0.0]
    C2_z = [0.0 0.0 1.0]
    C2list = (C2_x, C2_y, C2_z)
    namelist = ("x", "y", "z")
    for i = 1:3
        C2 = Cn(C2list[i], 2)
        push!(symels, Symel("C_2($(namelist[i]))", C2))
    end
    return symels
end

function generate_Td()
    symels = generate_T()
    # σd's
    σd_1v = normalize!([1.0 1.0 0.0])
    σd_2v = normalize!([1.0 -1.0 0.0])
    σd_3v = normalize!([1.0 0.0 1.0])
    σd_4v = normalize!([1.0 0.0 -1.0])
    σd_5v = normalize!([0.0 1.0 1.0])
    σd_6v = normalize!([0.0 1.0 -1.0])
    σs = (σd_1v,σd_2v,σd_3v,σd_4v,σd_5v,σd_6v)
    namelist = ("xyp","xym","xzp","xzm","yzp","yzm")
    for i = 1:6
        σd = σ(σs[i])
        push!(symels, Symel("σd($(namelist[i]))", σd))
    end
    # S4's
    S4_1v = [1.0 0.0 0.0]
    S4_2v = [0.0 1.0 0.0]
    S4_3v = [0.0 0.0 1.0]
    S4vlist = (S4_1v, S4_2v, S4_3v)
    namelist = ("x","y","z")
    for i = 1:3
        S4 = Sn(S4vlist[i], 4)
        S43 = S4^3
        push!(symels, Symel("S_4($(namelist[i]))", S4))
        push!(symels, Symel("S_4^3($(namelist[i]))", S43))
    end
    return symels
end

function generate_Th()
    symels = generate_T()
    # i
    push!(symels, Symel("i", i()))
    # S6
    S6_1v = normalize!([1.0 1.0 1.0])
    S6_2v = normalize!([-1.0 1.0 -1.0])
    S6_3v = normalize!([-1.0 -1.0 1.0])
    S6_4v = normalize!([1.0 -1.0 -1.0])
    S6list = (S6_1v, S6_2v, S6_3v, S6_4v)
    namelist = ("α", "β", "γ", "δ")
    for i = 1:4
        S6 = Sn(S6list[i], 6)
        S65 = S6^5
        push!(symels, Symel("S_6($(namelist[i]))", S6))
        push!(symels, Symel("S_6^5($(namelist[i]))", S65))
    end
    # 3σh
    σh_xv = [1.0 0.0 0.0]
    σh_yv = [0.0 1.0 0.0]
    σh_zv = [0.0 0.0 1.0]
    σlist = (σh_xv,σh_yv,σh_zv)
    namelist = ("x", "y", "z")
    for i = 1:3
        σh = σ(σlist[i])
        push!(symels, Symel("σh($(namelist[i]))", σh))
    end
    return symels
end

function generate_O()
    symels = []
    # C4
    C4_xv = [1.0 0.0 0.0]
    C4_yv = [0.0 1.0 0.0]
    C4_zv = [0.0 0.0 1.0]
    C4list = (C4_xv, C4_yv, C4_zv)
    namelist = ("x", "y", "z")
    for i = 1:3
        C4 = Cn(C4list[i], 4)
        C42 = C4^2
        C43 = C4^3
        push!(symels, Symel("C_4($(namelist[i]))", C4))
        push!(symels, Symel("C_2($(namelist[i]))", C42))
        push!(symels, Symel("C_4^3($(namelist[i]))", C43))
    end
    # C3
    C3_1v = normalize!([1.0 1.0 1.0])
    C3_2v = normalize!([1.0 -1.0 1.0])
    C3_3v = normalize!([1.0 1.0 -1.0])
    C3_4v = normalize!([1.0 -1.0 -1.0])
    C3list = (C3_1v, C3_2v, C3_3v, C3_4v)
    namelist = ("α", "β", "γ", "δ")
    for i = 1:4
        C3 = Cn(C3list[i], 3)
        C32 = C3^2
        push!(symels, Symel("C_3($(namelist[i]))", C3))
        push!(symels, Symel("C_3^2($(namelist[i]))", C32))
    end
    # C2
    C2_1v = normalize!([1.0 0.0 1.0])
    C2_2v = normalize!([1.0 0.0 -1.0])
    C2_3v = normalize!([1.0 1.0 0.0])
    C2_4v = normalize!([1.0 -1.0 0.0])
    C2_5v = normalize!([0.0 -1.0 1.0])
    C2_6v = normalize!([0.0 -1.0 -1.0])
    C2list = (C2_1v, C2_2v, C2_3v, C2_4v, C2_5v, C2_6v)
    namelist = ("xzp", "xzm", "xyp", "xym", "yzm", "yzp")
    for i = 1:6
        C2 = Cn(C2list[i],2)
        push!(symels, Symel("C_2($(namelist[i]))", C2))
    end
    return symels
end

function generate_Oh()
    symels = generate_O()
    symels = vcat(symels, Symel("i",i()))
    # S4 and σh
    S4_xv = [1.0 0.0 0.0]
    S4_yv = [0.0 1.0 0.0]
    S4_zv = [0.0 0.0 1.0]
    S4list = (S4_xv, S4_yv, S4_zv)
    namelist = ("x", "y", "z")
    for i = 1:3
        S4 = Sn(S4list[i], 4)
        σh = σ(S4list[i])
        S43 = S4^3
        push!(symels, Symel("S_4($(namelist[i]))", S4))
        push!(symels, Symel("σh($(namelist[i]))", σh))
        push!(symels, Symel("S_4^3($(namelist[i]))", S43))
    end
    # S6
    S6_1v = normalize!([1.0 1.0 1.0])
    S6_2v = normalize!([1.0 -1.0 1.0])
    S6_3v = normalize!([1.0 1.0 -1.0])
    S6_4v = normalize!([1.0 -1.0 -1.0])
    S6list = (S6_1v, S6_2v, S6_3v, S6_4v)
    namelist = ("α", "β", "γ", "δ")
    for i = 1:4
        S6 = Sn(S6list[i], 6)
        S65 = S6^5
        push!(symels, Symel("S_6($(namelist[i]))", S6))
        push!(symels, Symel("S_6^5($(namelist[i]))", S65))
    end
    # C2
    σd_1v = normalize!([1.0 0.0 1.0])
    σd_2v = normalize!([1.0 0.0 -1.0])
    σd_3v = normalize!([1.0 1.0 0.0])
    σd_4v = normalize!([1.0 -1.0 0.0])
    σd_5v = normalize!([0.0 -1.0 1.0])
    σd_6v = normalize!([0.0 -1.0 -1.0])
    σdlist = (σd_1v, σd_2v, σd_3v,σd_4v,σd_5v,σd_6v)
    namelist = ("xzp", "xzm", "xyp", "xym", "yzm", "yzp")
    for i = 1:6
        σd = σ(σdlist[i])
        push!(symels, Symel("σd($(namelist[i]))", σd))
    end
    return symels
end

function generate_I()
    symels = []
    faces, vertices, edgecenters = generate_I_vectors()

    # C5 (face vectors)
    for i = 1:6
        C5 = Cn(faces[i],5)
        C52 = C5^2
        C53 = C5^3
        C54 = C5^4
        push!(symels, Symel("C_5($i)", C5))
        push!(symels, Symel("C_5^2($i)", C52))
        push!(symels, Symel("C_5^3($i)", C53))
        push!(symels, Symel("C_5^4($i)", C54))
    end

    # C3 (vertex vectors)
    for i = 1:10
        C3 = Cn(vertices[i],3)
        C32 = C3^2
        push!(symels, Symel("C_3($i)", C3))
        push!(symels, Symel("C_3^2($i)", C32))
    end

    # C2 (edge vectors)
    for i = 1:15
        C2 = Cn(edgecenters[i],2)
        push!(symels, Symel("C_2($i)", C2))
    end

    return symels
end

function generate_Ih()
    symels = generate_I()
    faces, vertices, edgecenters = generate_I_vectors()
    push!(symels, Symel("i", i()))
    # S10 (face vectors)
    for i = 1:6
        S10 = Sn(faces[i],10)
        S103 = S10^3
        S107 = S10^7
        S109 = S10^9
        push!(symels, Symel("S_10($i)", S10))
        push!(symels, Symel("S_10^3($i)", S103))
        push!(symels, Symel("S_10^7($i)", S107))
        push!(symels, Symel("S_10^9($i)", S109))
    end

    # S6 (vertex vectors)
    for i = 1:10
        S6 = Sn(vertices[i],6)
        S65 = S6^5
        push!(symels, Symel("S_6($i)", S6))
        push!(symels, Symel("S_6^5($i)", S65))
    end

    # σ (edge vectors)
    for i = 1:15
        σi = σ(edgecenters[i])
        push!(symels, Symel("σ($i)", σi))
    end

    return symels
end

function generate_I_vectors()
    ϕ = (1+(5^0.5))/2
    ϕi = 1/ϕ
    faces = [[1.0 ϕ 0.0],[1.0 -ϕ 0.0],[-1.0 ϕ 0.0],[-1.0 -ϕ 0.0],
             [0.0 1.0 ϕ],[0.0 1.0 -ϕ],[0.0 -1.0 ϕ],[0.0 -1.0 -ϕ],
             [ϕ 0.0 1.0],[-ϕ 0.0 1.0],[ϕ 0.0 -1.0],[-ϕ 0.0 -1.0]]
    vertices = [[1.0 1.0 1.0],[1.0 1.0 -1.0],[1.0 -1.0 1.0],[-1.0 1.0 1.0],
                [1.0 -1.0 -1.0],[-1.0 1.0 -1.0],[-1.0 -1.0 1.0],[-1.0 -1.0 -1.0],
                [0.0 ϕ ϕi],[0.0 ϕ -ϕi],[0.0 -ϕ ϕi],[0.0 -ϕ -ϕi],
                [ϕi 0.0 ϕ],[-ϕi 0.0 ϕ],[ϕi 0.0 -ϕ],[-ϕi 0.0 -ϕ],
                [ϕ ϕi 0.0],[ϕ -ϕi 0.0],[-ϕ ϕi 0.0],[-ϕ -ϕi 0.0]]
    l = length(vertices)
    edglen = 2*ϕi
    edgecenters = []
    for i = 1:l
        for j = i+1:l
            if isapprox(distance(vertices[i], vertices[j]), edglen)
                v = normalize(vertices[i]+vertices[j])
                for k in edgecenters
                    if isapprox(abs(dot(k,v)), 1.0)
                        @goto here
                    end
                end
                push!(edgecenters, v)
                @label here
            end
        end
    end
    for i in faces
        normalize!(i)
    end
    face_vectors = []
    for i in faces
        v = normalize(i)
        for j in face_vectors
            if isapprox(abs(dot(v,j)), 1.0)
                @goto here2
            end
        end
        push!(face_vectors, v)
        @label here2
    end
    vertex_vectors = []
    for i in vertices
        v = normalize(i)
        for j in vertex_vectors
            if isapprox(abs(dot(v,j)), 1.0)
                @goto here3
            end
        end
        push!(vertex_vectors, v)
        @label here3
    end
    return (face_vectors, vertex_vectors, edgecenters)
end

function distance(a,b)
    return (sum((a-b).^2))^0.5
end

function reduce(n, i)
    g = gcd(n, i)
    return div(n, g), div(i, g)
end

function gcd(A, B)
    # A quick implementation of the Euclid algorithm for finding the greatest common divisor
    a = max(A,B)
    b = min(A,B)
    if a == 0
        return b
    elseif b == 0
        return a
    else
        r = a % b
        gcd(b, r)
    end
end


text:
function GP(y::REALVEC,
            x::REALMAT,
            z::REALMAT,
            K::Kernel,
            s2n::Real)

    @assert s2n > 0 "Output noise, s2n, must be positive."
    T, D = size(x)

    # Kernel
    kbb = kernel(x, K)    # Training
    kbt = kernel(x, z, K) # Training vs test
    # Cholesky
    L = cholesky(kbb + I*s2n).L
    # Predictive mean
    α = L' \ (L \ y)
    μ = kbt' * α
    # Predictive variance
    v = L \ kbt
    Σ = kernel(z, z, kbb) - v' * v
    # Log-likelihood
    LL = -dot(y, α)/2 - sum(diag(L)) - T/2 * log(2π)

    return (predictions = μ, covariance = Σ, loglikelihood = LL)
end

function GP(y::REALVEC,
            x::REALVEC,
            z::REALVEC,
            K::Kernel,
            s2n::Real = 1.0)
    GP(y, reshape(x, :, 1), reshape(z, :, 1), K, s2n)
end


text:
# Easing equation based on Robert Penner's work:
# http://robertpenner.com/easing/
# author
#    Aurelien Ribon | http://www.aurelienribon.com/ (Original java code)
#    Xavier Guzman (dart port)
#    Cleveland (julia port)

abstract type AbstractTweenEquation end

const EASE_IN = 0
const EASE_OUT = 1
const EASE_INOUT = 2

include("linear.jl")
include("expo.jl")
include("quad.jl")

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
IMGUI_DEMO_REPEATER_STATES =
    Dict("basic_repeater" => AppDemo.Repeaters.State("repeater"))


text:
@testset "TD_OPs" begin
#test transform-domain operators (linear operators)

# 2D discrete gradient
  n1=9
  n2=6
  h1=0.99
  h2=1.123
  D2D = get_discrete_Grad(n1,n2,h1,h2,"TV")
  D2x = get_discrete_Grad(n1,n2,h1,h2,"D_x")
  D2z = get_discrete_Grad(n1,n2,h1,h2,"D_z")

  x=zeros(n1,n2) #test on a 'cross' image
  x[:,3].=1.0
  x[4,:].=1.0

  a1 = D2x*vec(x); a1=reshape(a1,n1-1,n2)
  a2 = D2z*vec(x); a2=reshape(a2,n1,n2-1)
  a3 = D2D*vec(x); a3a=a3[1:(n2-1)*n1]; a3b=a3[1+(n2-1)*n1:end];
  a3a = reshape(a3a,n1,n2-1)
  a3b = reshape(a3b,n1-1,n2)

  #this test depends on the values of h1 and h2, as well as the type of derivative
  @test a1==diff(x, dims=1)./h1
  @test a2==diff(x, dims=2)./h2

  #some more general tests:
  @test count(!iszero, a1[:,3])==0
  for i in [1 2 4 5 6]
    @test a1[:,1]==a1[:,i]
  end

  #some more general tests:
  @test count(!iszero, a2[4,:])==0
  for i in [1 2 3 5 6 7 8 9]
   @test a2[1,:]==a2[i,:]
  end

  @test a3a==a2
  @test a3b==a1

  # 3D discrete gradient
    n1=4
    n2=6
    n3=5
    h1=0.99
    h2=1.123
    h3=1.0
    D3D = get_discrete_Grad(n1,n2,n3,h1,h2,h3,"TV")
    D3x = get_discrete_Grad(n1,n2,n3,h1,h2,h3,"D_x")
    D3y = get_discrete_Grad(n1,n2,n3,h1,h2,h3,"D_y")
    D3z = get_discrete_Grad(n1,n2,n3,h1,h2,h3,"D_z")

    x=zeros(n1,n2,n3) #test on a 'cross' image
    x[2,:,:].=1.0
    x[:,4,:].=1.0
    x[:,:,3].=1.0

    a1 = D3x*vec(x); a1=reshape(a1,n1-1,n2,n3)
    a2 = D3y*vec(x); a2=reshape(a2,n1,n2-1,n3)
    a3 = D3z*vec(x); a3=reshape(a3,n1,n2,n3-1)
    for i=1:n2
      @test a1[:,i,:]==diff(x[:,i,:], dims=1)./h1
    end
    for i=1:n3
      @test a1[:,:,i]==diff(x[:,:,i], dims=1)./h1
    end

    for i=1:n1
      @test a2[i,:,:]==diff(x[i,:,:], dims=1)./h2
    end
    for i=1:n3
      @test a2[:,:,i]==diff(x[:,:,i], dims=2)./h2
    end

    for i=1:n1
      @test a3[i,:,:]==diff(x[i,:,:], dims=2)./h3
    end
    for i=1:n2
      @test a3[:,i,:]==diff(x[:,i,:], dims=2)./h3
    end

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module Rauschen
    include("AdaptiveFilter.jl")
    # Empty
end


text:
include("football_form.jl")
using ..FootballForm

leagues = [
    "Barclays Premier League",
    "French Ligue 1",
    "German Bundesliga",
    "Italy Serie A",
    "Spanish Primera Division",
    "UEFA Champions League",
    "UEFA Europa League"
]

for league ∈ leagues
    games = load_league(league) |> load_games
    results = run_form_analysis(games)
    log_all_plots(results, games, league)
end

run(`git add -A`)
run(`git commit -a -m "Finished Run"`)
run(`git push`)

max_stars_repo_name: shalinkpatel/football-form

max_issues_repo_name: shalinkpatel/football-form

max_forks_repo_name: shalinkpatel/football-form


text:
using MortalityTables, Plots

tables = MortalityTables.tables()
cso_2001 = tables["2001 CSO Super Preferred Select and Ultimate - Male Nonsmoker, ANB"]
cso_2017 = tables["2017 Loaded CSO Preferred Structure Nonsmoker Super Preferred Male ANB"]

issue_ages = 18:80
durations = 1:40

a = [cso_2001.select[issue_age][issue_age + duration - 1] for issue_age in issue_ages,duration in durations]
b = [cso_2017.select[issue_age][issue_age + duration - 1] for issue_age in issue_ages,duration in durations]

# compute the relative rates with the element-wise division ("brodcasting" in Julia)
function rel_diff(a, b, issue_age,duration)
        att_age = issue_age + duration - 1
        return a[issue_age][att_age] / b[issue_age][att_age]
end

diff = [rel_diff(cso_2017.select,cso_2001.select,ia,dur) for ia in issue_ages, dur in durations]
x = contour(durations,
        issue_ages,
        diff,
        xlabel="duration",ylabel="issue ages",
        title="Relative difference between 2017 and 2001 CSO \n M PFN",
        fill=true
        )

savefig(x,"contour.png")


text:
export measure

"""
Computes the outcome probabilities for a quantum measurement. The conditional
probabilities are determined by the Born rule, ``P(i|j) = \\text{Tr}[\\Pi_i \\rho_j]``,
where ``\\Pi_j`` is a POVM element and ``\\rho_j`` is a density matrix.

Measurement of a single `Ket` or  `DensityMatrix`:

    measure(
        Π :: POVM,
        ρ :: State
    ) :: Conditionals

    measure(
        Π :: POVM,
        ψ :: Ket
    ) :: Conditionals

Measurement of an ensemble of `Ket` or `DensityMatrix` types:

    measure(
        Π :: POVM,
        ρ_states :: Vector{<:State}
    ) :: Conditionals

    measure(
        Π :: POVM,
        ψ_kets :: Vector{<:Ket}
    ) :: Conditionals
"""
function measure(
    Π :: POVM,
    ρ :: State
) :: Probabilities
    Probabilities( real.(map(Π_el -> tr(Π_el * ρ), Π)) )
end
function measure(
    Π :: POVM,
    states :: Vector{<:State}
) :: Conditionals
    Conditionals( real.(tr.(transpose(states * transpose(Π))) ))
end
function measure(
    Π :: POVM,
    ψ :: Ket
) :: Probabilities
    Probabilities( real.(map( Π_el -> ψ' * Π_el * ψ, Π)) )
end
function measure(
    Π :: POVM,
    ψ_kets :: Vector{<:Ket}
) :: Conditionals
    measure(Π, pure_state.(ψ_kets))
end


text:
# mpiexec.exe -n 4 julia .\mpisum.jl
using ADCME

mpi_init()
r = mpi_rank()
a = constant(ones(10) * r)
a = mpi_sendrecv(a, 0, 2)

# Equivalently, we can use
# if r==2
#     global a
#     a = mpi_send(a, 0)
# end
# if r==0
#     global a
#     a = mpi_recv(a,2)
# end

L = sum(a^2)
g = gradients(L, a)

sess = Session(); init(sess)
v, G = run(sess, [a,g])

# processor 0 should have the same `v` as processor 2
@info r, v

G = run(sess, g)

# gradients on processor 0 should be the same as processor 2
# because `a` is received from 2
@info r, G

mpi_finalize()


text:
export readModelAndGenerateMeshMref
# 1.75 and 2.9 are the velocities that are suitable for the SEG salt model

function readModelAndGenerateMeshMref(readModelFolder::String,modelFilename::String,dim::Int64,pad::Int64,domain::Vector{Float64},newSize::Vector=[],velBottom::Float64=1.75,velHigh::Float64=2.9,
	doTranspose=true)
########################## m,mref are in Velocity here. ###################################

if dim==2
	# SEGmodel2Deasy.dat
	m = readdlm(string(readModelFolder,"/",modelFilename));
	m = m*1e-3;
	if doTranspose
		m = Matrix(m');
	else
		m = Matrix(m);
	end
	mref = getSimilarLinearModel(m,velBottom,velHigh);
else
	# 3D SEG slowness model
	file = matopen(string(readModelFolder,"/",modelFilename)); DICT = read(file); close(file);
	m = DICT["VELs"];
	m = m*1e-3;
	mref = getSimilarLinearModel(m,velBottom,velHigh);
end

sea = abs.(m[:] .- minimum(m)) .< 7e-2;
mref[sea] = m[sea];
if newSize!=[]
	m    = expandModelNearest(m,   collect(size(m)),newSize);
	mref = expandModelNearest(mref,collect(size(mref)),newSize);
end

Minv = getRegularMesh(domain,collect(size(m)));

(mPadded,MinvPadded) = addAbsorbingLayer(m,Minv,pad);
(mrefPadded,MinvPadded) = addAbsorbingLayer(mref,Minv,pad);

N = prod(MinvPadded.n);
boundsLow  = minimum(mPadded);
boundsHigh = maximum(mPadded);

boundsLow  = ones(N)*boundsLow;
boundsLow = convert(Array{Float32},boundsLow);
boundsHigh = ones(N)*boundsHigh;
boundsHigh = convert(Array{Float32},boundsHigh);

return (mPadded,MinvPadded,mrefPadded,boundsHigh,boundsLow);
end


text:
module PlayGround

    import Base: product
    using StaticArrays

    abstract type Grid{d} end

    struct UCGrid{d} <: Grid{d}
        a::SVector{d, Float64}
        b::SVector{d, Float64}
        n::SVector{d, Int64}
        nodes:: Vector{SVector{d, Float64}}
    end

    ndims(g::Grid{d}) where d = d

    import Base
    Base.iterate(ucg::UCGrid{d}, args...) where d = iterate(ucg.nodes, args...)
    Base.length(ucg::UCGrid{d}) where d = length(ucg.nodes)
    Base.eltype(ucg::UCGrid{d}) where d = SVector{d, Float64}

    function UCGrid{d}(a::SVector{d, Float64},
                         b::SVector{d, Float64},
                         n::SVector{d, Int64}
    ) where d
        intervals = tuple( [range(a[i], b[i]; length=n[i]) for i=1:d]... )
        # print( [product(intervals)...])
        # nn = [(e) for e in product(intervals...)]
        # return nn
        cnodes = [SVector(e...) for e in product(intervals...)]
        # TODO: can we avoid the cnodes[:] allocation
        return UCGrid{d}(a,b,n,cnodes[:])
    end

    function UCGrid(
        a::AbstractVector{Float64},
        b::AbstractVector{Float64},
        n::AbstractVector{Int64})

        d = length(n)
        @assert (length(a)==d) & (length(b)==d)

        aa = SVector(a...)
        bb = SVector(b...)
        nn = SVector(n...)

        return UCGrid{d}(aa, bb, nn)
    end

    # Concatenate SVectors
end

import Main.PlayGround

grid_1 = Main.PlayGround.UCGrid([0.0, 0.0], [1.0, 1.0], [2, 2])
grid_2 = Main.PlayGround.UCGrid([0.0, 0.0], [1.0, 1.0], [2, 2])

pg = Main.PlayGround.PGrid(grid_1, grid_2)

show(pg)

for k in Main.PlayGround.iter(pg)
    println(k)
end

for k in Main.PlayGround.enumerate_product(pg)
    println(k)
end


text:
@inline function f_inline(x)
a=x*5
b=a+3
c=a-4
d=b/c
end
g_inline(x) = f_inline(2*x)

max_stars_repo_path: Module 2/Chapter04/f_inline.jl

max_issues_repo_path: Module 2/Chapter04/f_inline.jl

max_forks_repo_path: Module 2/Chapter04/f_inline.jl


text:
using DataFrames

@testset "calibrate" begin
    @testset "basic" begin
        @system SCalibrate(Controller) begin
            a => 0 ~ preserve(parameter)
            b(a) ~ accumulate
        end
        n = 10
        t, a, b = 10.0u"hr", 20, 200
        A = (0.0, 100.0)
        obs = DataFrame(time=[t], b=[b])
        p = calibrate(SCalibrate, obs, stop=n, target=:b, parameters=("SCalibrate.a" => A))
        @test p[:SCalibrate][:a] == a
        r = simulate(SCalibrate, stop=n, config=p)
        @test r[r.time .== t, :][1, :b] == b
    end

    @testset "unit" begin
        @system SCalibrateUnit(Controller) begin
            a => 0 ~ preserve(parameter, u"m/hr")
            b(a) ~ accumulate(u"m")
        end
        n = 10
        t, a, b = 10.0u"hr", 20u"m/hr", 200u"m"
        #FIXME: parameter range units are just ignored
        A = [0.0, 100.0]u"m/hr"
        obs = DataFrame(time=[t], b=[b])
        p = calibrate(SCalibrateUnit, obs, stop=n, target=:b, parameters=("SCalibrateUnit.a" => A))
        @test p[:SCalibrateUnit][:a] == Cropbox.deunitfy(a)
        r = simulate(SCalibrateUnit, stop=n, config=p)
        @test r[r.time .== t, :][1, :b] == b
    end

    @testset "config" begin
        @system SCalibrateConfig(Controller) begin
            a => 0 ~ preserve(parameter)
            w => 1 ~ preserve(parameter)
            b(a, w) => w*a ~ accumulate
        end
        n = 10
        t, a, b = 10.0u"hr", 20, 200
        w1, w2 = 1, 2
        A = (0.0, 100.0)
        obs = DataFrame(time=[t], b=[b])
        params = :SCalibrateConfig => :a => A
        p1 = calibrate(SCalibrateConfig, obs, stop=n, target=:b, config=(:SCalibrateConfig => :w => w1), parameters=params)
        @test p1[:SCalibrateConfig][:a] == a/w1
        p2 = calibrate(SCalibrateConfig, obs, stop=n, target=:b, config=(:SCalibrateConfig => :w => w2), parameters=params)
        @test p2[:SCalibrateConfig][:a] == a/w2
    end

    @testset "configs as index" begin
        @system SCalibrateConfigsIndex(Controller) begin
            a => 0 ~ preserve(parameter)
            w => 1 ~ preserve(parameter)
            b(a, w) => w*a ~ accumulate
        end
        n = 10
        t, w, b = [10.0u"hr", 10.0u"hr"], [1, 2], [100, 200]
        A = (0.0, 100.0)
        obs = DataFrame(; time=t, w, b)
        configs = [
            :SCalibrateConfigsIndex => :w => 1,
            :SCalibrateConfigsIndex => :w => 2,
        ]
        index = [:time => "context.clock.time", :w]
        target = :b
        params = :SCalibrateConfigsIndex => :a => A
        p = calibrate(SCalibrateConfigsIndex, obs, configs; stop=n, index, target, parameters=params)
        @test p[:SCalibrateConfigsIndex][:a] == 10
    end
end


text:
@testset "Encoding: One-Hot" begin
    x = [0,1,2,3,4]
    b = [1,0,0,0,0, 0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1]

    @test collect(binarize(x)) == b
    @test debinarize(b) == x

    d = domain(x)
    @test debinarize(collect(binarize(x, d)), d) == x
end

@testset "Encoding: Domain Wall" begin
    x = [0,1,2,3,4,5]

    d = domain([0,1,2,3,4])
    y = [0,1,2,3,4]

    a = [0,0,0,0,0, 1,0,0,0,0, 1,1,0,0,0, 1,1,1,0,0, 1,1,1,1,0, 1,1,1,1,1]

    @test collect(binarize(x; binarization = :domain_wall)) == a
    @test debinarize(a; binarization = :domain_wall) == x

    @test debinarize(collect(binarize(y, d; binarization = :domain_wall)), d; binarization = :domain_wall) == y
end


text:
include("led.jl")
include("sensors.jl")


text:
using DReal
using Base.Test

# Ackley's: bounds -5 <= x,y <= 5, min: f(0,0) = 0
begin
println("Trying Ackleys")
ackleys(x::Array) = -20*exp(-0.2*sqrt(0.5(x[1]*x[1] + x[2]*x[2]))) - exp(0.5(cos(2*π*x[1])+cos(2*π*x[2]))) + e + 20
x = Var(Float64,"x2",-5.12,5.12)
y = Var(Float64,"y2",-5.12,5.12)
f = Var(Float64,"f2",-10.,10.)
add!(f == ackleys([x,y]))
cost, optimal_model = minimize(f,x,y; lb=-10.,ub = 10.)

reset_ctx!(default_global_context)

# Commented out by Soonho Kong, 2016/06/09
# The following optimization query returns UNSAT, which is wrong.
# Until I fix the problem, we don't run this test
#
# # # Beales function
# println("Trying Beales")
# beale(x::Array) =  (1.5 - x[1] + x[1]x[2])^2 +
#                 (2.25 - x[1] + x[1]x[2]^2)^2 +
#                 (2.625 - x[1] + x[1]x[2]^3)^2
# x = Var(Float64,"x3",-5.12,5.12)
# y = Var(Float64,"y3",-5.12,5.12)
# f = Var(Float64,"f3",-10.,10.)
# add!(f == beale([x,y]))
# cost, optimal_model = minimize(f,x,y; lb=-10.,ub = 10.)
# @show cost, optimal_model
#
# reset_global_ctx!()

# # rastrigin
println("Trying rastrigin")

rastrigin(x::Array) = 10 * length(x) + sum([xi*xi - 10*cos(2pi*xi) for xi in x])
x = Var(Float64,"x4",-5.12,5.12)
y = Var(Float64,"y4",-5.12,5.12)
f = Var(Float64,"f4",-10.,10.)
add!(f == rastrigin([x,y]))
cost, optimal_model = minimize(f,x,y; lb=-10.,ub = 10.)
@show cost, optimal_model

reset_global_ctx!()

end


text:
module OnlineStats

using Statistics, LinearAlgebra, Dates

import OnlineStatsBase: value, name, _fit!, _merge!, bessel, pdf, probs, smooth, smooth!,
    smooth_syr!, nvars, Weight, eachrow, eachcol, TwoThings
import StatsBase: fit!, nobs, autocov, autocor, confint, skewness, kurtosis, entropy, midpoints,
    fweights, sample, coef, predict, Histogram, ecdf, transform, log1p

import AbstractTrees: AbstractTrees

using OrderedCollections: OrderedDict

using OnlineStatsBase, RecipesBase

@static if VERSION < v"1.1.0"
    isnothing(x) = x === nothing
end

export
# Statistics
    mean, var, std, cov, cor,
# functions
    fit!, nobs, value, autocov, autocor, predict, confint, probs, skewness, kurtosis,
    classify, coef, ecdf, eachrow, eachcol,
# weights
    EqualWeight, ExponentialWeight, LearningRate, LearningRate2, HarmonicWeight,
    McclainWeight, Bounded, Scaled,
# algorithms
    ADAGRAD, ADAM, ADAMAX, ADADELTA, MSPI, OMAS, OMAP, RMSPROP, SGD,
# stats
    Ash,
    AutoCov,
    Bootstrap,
    CallFun, Counter, CountMap, CountMissing, CovMatrix, CCIPCA,
    Diff,
    Extrema, ExpandingHist,
    FitBeta, FitCauchy, FitGamma, FitLogNormal, FitNormal, FitMultinomial, FitMvNormal,
    FastNode, FastTree, FastForest, FTSeries,
    Group, GroupBy,
    HeatMap, Hist, HyperLogLog,
    IndexedPartition,
    KHist, KHist2D, KMeans, KahanSum, KahanMean, KahanVariance, KIndexedPartition,
    Lag, LinReg, LinRegBuilder,
    Mean, Moments, ModelSchema, Mosaic, MovingTimeWindow, MovingWindow,
    NBClassifier,
    OrderStats,
    Part, Partition, PlotNN, ProbMap, P2Quantile,
    Quantile,
    ReservoirSample,
    Series, SGDStat, StatHistory, StatLag, StatLearn, Sum,
    Variance,
# other
    OnlineStat, BiasVec

#-----------------------------------------------------------------------------# utils
const Tup = Union{Tuple, NamedTuple}
const VectorOb = Union{AbstractVector, Tup}
const XY{T,S} = Union{Tuple{T,S}, Pair{T,S}, NamedTuple{names,Tuple{T,S}}} where {names,T<:AbstractVector{<:Number},S<:Number}

const ϵ = 1e-7  # avoid dividing by 0 in some cases

neighbors(x) = ((x[i], x[i+1]) for i in eachindex(x)[1:end-1])

#-----------------------------------------------------------------------# BiasVec
"""
    BiasVec(x)

Lightweight wrapper of a vector which adds a bias/intercept term at the end.

# Example

    BiasVec(rand(5))
"""
struct BiasVec{T, A <: VectorOb} <: AbstractVector{T}
    x::A
    bias::T
end
BiasVec(x::AbstractVector{T}) where {T} = BiasVec(x, one(T))
BiasVec(x::Tup) = BiasVec(x, one(typeof(first(x))))

Base.length(v::BiasVec) = length(v.x) + 1
Base.size(v::BiasVec) = (length(v), )
Base.getindex(v::BiasVec, i::Int) = i > length(v.x) ? v.bias : v.x[i]
Base.IndexStyle(::Type{<:BiasVec}) = IndexLinear()

#-----------------------------------------------------------------------------# stats
include("algorithms.jl")

include("stats/histograms.jl")
include("stats/stats.jl")
include("stats/distributions.jl")
include("stats/nbclassifier.jl")
include("stats/fasttree.jl")
include("stats/linreg.jl")
include("stats/kahan.jl")
include("stats/pca.jl")
include("stats/statlearn.jl")
#-----------------------------------------------------------------------------# viz
include("viz/khist.jl")
include("viz/khist2d.jl")
include("viz/partition.jl")
include("viz/mosaicplot.jl")
include("viz/recipes.jl")
include("viz/heatmap.jl")
include("viz/hexlattice.jl")
include("viz/ash.jl")
end


text:
@testset "verify perron frobenius " begin
    A = [1 2;3 4]
    ρ = bound_perron_frobenius_eigenvalue(A)
    @test  (5+sqrt(big(5)))/2 ≤ ρ
end

@testset "verified eigenvalues" begin
    n = 5 # matrix size

    # symmetric case
    ev = sort(randn(n))
    D = Diagonal(ev)
    Q, _ = qr(rand(n, n))
    A = Symmetric(IA.Interval.(Q) * D * IA.Interval.(Q)')

    evals, evecs, cert = verify_eigen(A)
    @test all(cert)
    @test all(ev .∈ evals)

    # real eigenvalues case
    P = rand(n, n)
    Pinv, _ = epsilon_inflation(P, Diagonal(ones(n)))
    A = IA.Interval.(P) * D * Pinv

    evals, evecs, cert = verify_eigen(A)
    @test all(cert)
    @test all(ev .∈ evals)

    # test complex eigenvalues
    ev = sort(rand(Complex{Float64}, n), by = x -> (real(x), imag(x)))
    A = IA.Interval.(P) * Matrix(Diagonal(ev)) * Pinv

    evals, evecs, cert = verify_eigen(A)
    @test all(cert)
    @test all(ev .∈ evals)
end


text:
function create(c::Client, ::Type{Ban}, g::AbstractGuild, u::User; kwargs...)
    return create_guild_ban(c, g.id, u.id; kwargs...)
end

function retrieve(c::Client, ::Type{Ban}, g::AbstractGuild, u::User)
    return get_guild_ban(c, g.id, u.id)
end
retrieve(c::Client, ::Type{Ban}, g::AbstractGuild) = get_guild_bans(c, g.id)

function delete(c::Client, b::Ban, g::AbstractGuild)
    return remove_guild_ban(c, g.id, b.user.id)
end