text:
"""
    MushroomTools
Module containing many functions helpful in simulating (perfect) mushroom
billiards, see [`billiard_mushroom`](@ref).
Contains stuff like initializing efficiently regular or chaotic particles
and functions that return the corresponding chaotic or regular phase-space volumes
or portions.
The functions [`V_3D_tot`](@ref) and [`V_3D_reg`](@results) use equations derived
in ref. [1].

Made by Lukas Hupe.

## References
[1] A. Barnett & T. Betcke, [Chaos **17**, 043125 (2007)](https://doi.org/10.1063/1.2816946).
"""
module MushroomTools

using DynamicalBilliards
using DynamicalBilliards: SV, cossin
using LinearAlgebra

#="""
    is_regular(p::Particle, o::Semicircle, w)
checks whether the particle `p` is on an integrable orbit in a mushroom billiard
with cap `o` and stem width `w`
"""=#
function is_regular(p::Particle{T}, l, w, r) where {T}
    δ = SV{T}(0.0, l) - p.pos
    vperp = [-p.vel[2],p.vel[1]]
    if norm(δ) > r
        #debug && println("not inside circle: pos = $(p.pos)")
        return false
    elseif p.pos[2] < l
        #debug && println("wrong half of circle")
        return false
    elseif abs(dot(vperp, δ)) < w/2
        #debug && println("not outside critical radius")
        return false
    else
        #debug && println("Perfectly happy with $(p.pos) and $(p.vel)")
        return true
    end
end

#Formulae from quantum mushroom paper
"""
    V_3D_tot(l,w,r)
Return the total phasespace volume (3D) of a [`billiard_mushroom`](@ref)
parameterized by `(l,w,r)`.
"""
V_3D_tot(l,w,r) = 2π*(l*w + (1/2)*π*r^2)
V_3D_reg(l,w,r) = 2π*r^2*(acos(w/(2r)) - w/(2r)*sqrt(1 - (w^2)/(4r^2)))
V_3D_cha(l,w,r) = V_3D_tot(l,w,r) - V_3D_reg(l,w,r)

"""
    g_r_3D(l, w, r)
Return the regular phasespace portion of the full (3D) phase-space of a
[`billiard_mushroom`](@ref) with stem length `l`, stem width `w` and cap radious `r`.
This result is known analytically, see [`MushroomTools`](@ref) for references.
"""
g_r_3D(l,w,r) = V_3D_reg(l,w,r)/V_3D_tot(l,w,r)

"""
    g_c_3D(l, w, r)
Return the chaotic phasespace portion of the full (3D) phase-space of a
[`billiard_mushroom`](@ref) with stem length `l`, stem width `w` and cap radious `r`.
This result is known analytically, see [`MushroomTools`](@ref) for references.
"""
g_c_3D(l,w,r) = 1 - g_r_3D(l,w,r)

"""
    V_2D_tot(l,w,r)
Return the total boundary map volume (2D) of a [`billiard_mushroom`](@ref)
parameterized by `(l,w,r)`.
"""
V_2D_tot(l, w, r) = 2(π*r + 2r + 2l)
V_2D_reg(l, w, r) = 2π*r*(1 - w/2r) + 2sqrt(4r*r - w*w) -2w*acos(w/2r)
V_2D_cha(l,w,r) = V_2D_tot(l,w,r) - V_2D_reg(l,w,r)

"""
    g_r_2D(l, w, r)
Return the regular phasespace portion of the boundary map (2D) of a
[`billiard_mushroom`](@ref) with stem length `l`, stem width `w` and cap radious `r`.
This result is known analytically, see [`MushroomTools`](@ref) for references.
"""
g_r_2D(l,w,r) = V_2D_reg(l,w,r)/V_2D_tot(l,w,r)

"""
    g_c_2D(l, w, r)
Return the chaotic phasespace portion of the boundary map (2D) of a
[`billiard_mushroom`](@ref) with stem length `l`, stem width `w` and cap radious `r`.
This result is known analytically, see [`MushroomTools`](@ref) for references.
"""
g_c_2D(l,w,r) = 1 - g_r_2D(l,w,r)

#applying Kac's Lemma to the stem
τ_r_2D_kac(l,w,r) = (V_2D_tot(l,w,r) - V_2D_reg(l,w,r)) / (2w + 4l)

#="""
    insidemushroom(pos::StaticVector, l, w, r)
Returns `true` if pos is within the mushroom parameterised by `l`,  `w` and `r`
"""=#
function insidemushroom(pos::SV{T}, l::T, w::T, r::T) where {T <: AbstractFloat}
    if pos[2] > 0
        if pos[2] <= l
            if abs(pos[1]) <= w/2
                return true
            end
        elseif pos[2] <= l + r
            if pos[1]^2 + (pos[2] - l)^2 <= r^2
                return true
            end
        end
    end
    return false
end

"""
    randin_mushroom(l, w, r [, ω])
Generate a random particle within the [`billiard_mushroom`](@ref)
parameterised by `l`, `w` and `r`.
If `ω` is given the particle is magnetic instead.

This function is much more efficient than [`randominside`](@ref).
"""
function randin_mushroom(l::T = 1.0, w::T = 0.2, r::T = 1.0) where {T <: AbstractFloat}
    pos, vel = _randin_mushroom(l, w, r)
    return Particle(pos, vel, zero(SV{T}))
end
function randin_mushroom(l::T, w::T, r::T, ω::T) where {T <: AbstractFloat}
    pos, vel = _randin_mushroom(l, w, r)
    return MagneticParticle(pos, vel, zero(SV{T}), T(ω))
end
function _randin_mushroom(l::T = 1.0, w::T = 0.2, r::T = 1.0) where {T <: AbstractFloat}
    x = (-r, r, 2r)
    y = (0, l + r, l + r)
    pos = SV{T}(rand(T)*x[3] + x[1], rand(T)*y[3] + y[1])
    while ! insidemushroom(pos, l, w, r)
        pos = SV{T}(rand(T)*x[3] + x[1], rand(T)*y[3] + y[1])
    end

    φ = T(2π * rand(T))
    vel = SV{T}(cossin(φ)...)

    return pos, vel
end

"""
    randomchaotic(l, w, r)
Generate a chaotic particle, i.e. not trapped in the cap.
"""
function randomchaotic(l, w, r)
    p = randin_mushroom(l, w, r)
    while is_regular(p, l, w, r)
        p = randin_mushroom(l, w, r)
    end
    return p
end

"""
    randomregular(l, w, r)
Generate a regular particle (i.e. trapped in the cap).
"""
function randomregular(l, w, r)
    p = randin_mushroom(l, w, r)
    while !is_regular(p, l, w, r)
        p = randin_mushroom(l, w, r)
    end
    return p
end

end # module


text:
struct WeightingPotential{T, N, S} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S}
end

@inline size(wp::WeightingPotential{T, N, S}) where {T, N, S} = size(wp.data)
@inline length(wp::WeightingPotential{T, N, S}) where {T, N, S} = length(wp.data)
@inline getindex(wp::WeightingPotential{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(wp.data, I...)
@inline getindex(wp::WeightingPotential{T, N, S}, i::Int) where {T, N, S} = getindex(wp.data, i)
@inline getindex(wp::WeightingPotential{T, N, S}, s::Symbol) where {T, N, S} = getindex(wp.grid, s)

function WeightingPotential(fss::PotentialSimulationSetup{T, 3, :cylindrical}; kwargs...)::WeightingPotential{T, 3, :cylindrical} where {T <: SSDFloat}
    return get_2π_potential(WeightingPotential{T, 3, :cylindrical}(fss.potential, fss.grid); kwargs...)
end
function WeightingPotential(fss::PotentialSimulationSetup{T, 3, :cartesian})::WeightingPotential{T, 3, :cartesian} where {T <: SSDFloat}
    return WeightingPotential{T, 3, :cartesian}(fss.potential, fss.grid)
end

function WeightingPotential(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    WeightingPotential{T, N, S}( nt.values, grid)
end
Base.convert(T::Type{WeightingPotential}, x::NamedTuple) = T(x)

function NamedTuple(ep::WeightingPotential{T, 3}) where {T}
    return (
        grid = NamedTuple(ep.grid),
        values = ep.data,
    )
end
Base.convert(T::Type{NamedTuple}, x::WeightingPotential) = T(x)

"""
    PointTypes(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::PointTypes{T, 3, :cylindrical}

Extracts the electric potential from `setup` and extrapolate it to an 2π grid.

For 2D grids (r and z) the user has to set the keyword `n_points_in_φ::Int`, e.g.: `n_points_in_φ = 36`.
"""
function PointTypes(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::PointTypes{T, 3, :cylindrical} where {T}
    return get_2π_potential(PointTypes{T, 3, :cylindrical}(setup.pointtypes, setup.grid); kwargs...)
end

ext: jl


text:
using Test
using NeXLSpectrum
using NeXLCore
using Statistics
using Distributions

@testset "Spectrum" begin
    les = LinearEnergyScale(-495.0, 5.0)

    @testset "Energy Scale" begin
        @test NeXLSpectrum.energy(200, les) == -495.0 + 5.0 * (200 - 1)
        @test NeXLSpectrum.energy(1, les) == -495.0 + 5.0 * (1 - 1)
        @test NeXLSpectrum.energy(500, les) == -495.0 + 5.0 * (500 - 1)

        @test channel(NeXLSpectrum.energy(200, les) + 1.0, les) == 200
        @test channel(NeXLSpectrum.energy(200, les) - 1.0, les) == 199
        @test channel(NeXLSpectrum.energy(200, les), les) == 200

        @test channel(600.0, les) == 220
        @test channel(1600.0, les) == round((1600.0 + 495.0) / 5.0) + 1
    end

    mnk = MnKaResolution(130.0)

    @testset "Resolution" begin
        @test resolution(5898.7, mnk) == 130.0
        @test isapprox(resolution(500.0, mnk), 60.0, atol = 1.0)
    end

    @testset "Spectrum" begin
        s1 = Spectrum(
            LinearEnergyScale(0.0, 10.0),
            collect(1:2048),
            Dict( :BeamEnergy => 20.0e3, :LiveTime => 72.0, :ProbeCurrent => 1.2 ),
        )

        @test s1[:LiveTime] == 72.0
        @test s1[100] == 100
        @test s1[1000] == 1000

        @test s1[100:10:150] == [100, 110, 120, 130, 140, 150]
        s1[100] = 200
        s1[101:103] = [300, 400, 500]

        @test s1[99] == 99
        @test s1[100] == 200
        @test s1[101] == 300
        @test s1[102] == 400
        @test s1[103] == 500
        @test s1[104] == 104

        s1[:RealTime] = 90.0
        @test s1[:RealTime] == 90.0

        s1n = dosenormalize(s1)
        @test isapprox(dose(s1n), 60.0, rtol = 1.0e-10)
        @test isapprox(s1[100]/dose(s1), s1n[100]/dose(s1n),rtol = 1.0e-8)
        @test isapprox(s1[200]/dose(s1), s1n[200]/dose(s1n),rtol = 1.0e-8)
    end

    @testset "SimpleEDS" begin
        eds = simpleEDS(4096, 5.0, -495.0, 125.0)
        @test channelcount(eds) == 4096
        @test isapprox(resolution(10000.0, eds), 160.0, atol = 1.0)
        @test NeXLSpectrum.energy(2000, eds) == NeXLSpectrum.energy(2000, les)
        @test channel(NeXLSpectrum.energy(2000, eds), eds) == 2000
    end
    @testset "BasicEDS" begin
        det3 = BasicEDS(
            4095,
            0.0,
            10.0,
            126.0,
            10,
            Dict(
                Shell(1) => n"Be",
                Shell(2) => n"Sc",
                Shell(3) => n"Cs",
                Shell(4) => n"Pu",
            ),
        )
        @test all(ch -> channel(energy(ch, det3), det3) == ch, 1:4096)
        @test resolution(energy(n"Mn K-L3"), det3) == 126.0
        @test isapprox(resolution(energy(n"C K-L2"), det3), 45.867, atol = 0.001)
        @test isvisible(n"C K-L2", det3)
        @test !isvisible(n"Ca L3-M3", det3)
        @test length(isvisible(characteristic(n"Sc", ltransitions), det3)) >= 7
        @test !isvisible(n"Xe M3-N5", det3)
        @test isvisible(n"Cs M3-N5", det3)
        @test length(isvisible(characteristic(n"U", ntransitions), det3)) == 0
        @test length(isvisible(characteristic(n"Pu", ntransitions), det3)) == 0 # Database doesn't contain any...
        @test channel(n"Mn K-L3", det3) == channel(energy(n"Mn K-L3"), det3)
    end
    @testset "readEMSA" begin
        sp = loadspectrum(joinpath(@__DIR__, "K412 spectra", "Al2O3 std.msa"))
        @test sp[:BeamEnergy] == 20.0e3
        @test sp[:ProbeCurrent] == 1.10989
        @test sp[:Elevation] == deg2rad(35.0)
        @test sp[:LiveTime] == 1172.19288
        @test sp[:RealTime] == 1491.4828
        @test sp[:Owner] == "Unknown"
        @test sp[:Name] == "Al2O3 std"
        comp = sp[:Composition]

        @test isapprox(comp[n"O"], 0.470749, atol = 1.0e-6)
        @test isapprox(comp[n"Al"], 0.529251, atol = 1.0e-6)
        @test comp[n"Zr"] == 0.0
        @test density(comp) == 4.00

        chs = channel(1335.0, sp):channel(1688.0, sp)
        @test integrate(sp, chs) == 33143490

        low, high = 120:133, 172:180
        @test isapprox(
            NeXLUncertainties.value(integrate(sp, low, chs, high)),
            3.22027e7,
            rtol = 3.0e-3,
        )

        ss = [subsample(sp, 0.1) for _ = 1:10]

        @test isapprox(
            sp[:LiveTime],
            sum(ss[i][:LiveTime] for i in eachindex(ss)),
            rtol = 1.0e-9,
        )
        @test isapprox(
            sp[:RealTime],
            sum(ss[i][:RealTime] for i in eachindex(ss)),
            rtol = 1.0e-9,
        )

        ch = channel(1450.0, sp)
        @test isapprox(
            sum(ss[i][ch] for i in eachindex(ss)),
            sp[ch],
            atol = 4.0 * sqrt(sp[ch]),
        ) # can fail occasionally...

    end
    @testset "SLR" begin
        x = [ 5.7, 2.5, 3.8, 2.8, 4.6, 1.6, 1.6, 1.4 ]
        y = [ 18, 19, 20, 19, 17, 32, 29, 30 ]
        slr = NeXLSpectrum.simple_linear_regression(x, y)
        @test isapprox(slr[1],-3.181, atol=1.0e-3)
        @test isapprox(slr[2], 32.542, atol=1.0e-3)
        @test isapprox(slr[3], -0.8072, atol=1.0e-4)
        @test isapprox(slr[4], -3.3497, atol=1.0e-4)
    end
    @testset "Duane-Hunt" begin
        fns = ( "Al2O3 std", "CaF2 std", "Fe std", "MgO std", "SiO2 std", ("III-E K412[$i][4]" for i in 0:4)...)
        specs = map(fn->loadspectrum(joinpath(@__DIR__,"K412 spectra",fn*".msa")),fns)
        @test all(isapprox.(NeXLSpectrum.duane_hunt.(specs), 20.0e3, atol=100.0))
    end

    @testset "sigma" begin
        specs = loadspectrum.(joinpath(@__DIR__,"K412 spectra","III-E K412[$i][4].msa") for i in 0:4)
        ss = [ sigma(specs[i], specs, 1:2000) for i in eachindex(specs) ]
        @test isapprox(mean(ss[1]), -0.05240393674216127, atol=1.0e-8)
        @test isapprox(mean(ss[2]), -0.02711953072380802, atol=1.0e-8)
        @test isapprox(mean(ss[3]), 0.00900794633309199, atol=1.0e-8)
        @test isapprox(mean(ss[4]), 0.1068253967760764, atol=1.0e-8)
        @test isapprox(mean(ss[5]), 0.11115156896689808, atol=1.0e-8)
        @test isapprox(std(ss[1]), 1.00731886928230457, atol=1.0e-8)
        @test isapprox(std(ss[2]), 1.0422450376085455, atol=1.0e-8)
        @test isapprox(std(ss[3]), 1.0294978144777747, atol=1.0e-8)
        @test isapprox(std(ss[4]), 0.9931443781679111, atol=1.0e-8)
        @test isapprox(std(ss[5]), 0.9858796014395503, atol=1.0e-8)

        n=fit(Normal, ss[1])
        @test isapprox(n.μ, -0.05240393674216127, atol = 1.0e-8)
        @test isapprox(n.σ, 1.007067008078397, atol = 1.0e-8)

        s1 = convert(Spectrum{Float64}, specs[1])
        @test s1 isa Spectrum{Float64}
        @test eltype(s1) == Float64
        ns1 = dosenormalize(s1, 10.0)
        @test isapprox(dose(ns1), 10.0, atol=1.0e-6)
        @test isapprox(ns1[100], s1[100]*10.0/dose(s1), rtol = 1.0e-6)
        @test isapprox(ns1[20], s1[20]*10.0/dose(s1), rtol = 1.0e-6)

        ss1 = similar(s1, Int)
        @test ss1 isa Spectrum{Int}
        @test eltype(ss1) == Int
        @test length(ss1) == length(s1)
        @test ss1.energy === s1.energy
        @test isequal(ss1.properties, s1.properties)
    end

    @testset "uv" begin
        spec = loadspectrum(joinpath(@__DIR__,"K412 spectra","III-E K412[0][4].msa"))
        u = uv(spec)
        @test eltype(u) == UncertainValue
        @test length(u) == length(spec)
        @test all(value(u[i]) == spec[i] for i in eachindex(spec))
        @test all(value(u[i])==0.0 || σ(u[i]) == sqrt(spec[i]) for i in eachindex(spec))
        @test all(value(u[i])!=0.0 || σ(u[i]) == 1.0 for i in eachindex(spec))
    end
end


text:
#
# InfiniteMPO
#

mutable struct InfiniteMPO <: AbstractInfiniteMPS
  data::CelledVector{ITensor}
  llim::Int #RealInfinity
  rlim::Int #RealInfinity
  reverse::Bool
end

## struct InfiniteSumLocalOps <: AbstractInfiniteMPS
##   data::CelledVector{ITensor}
## end
## InfiniteSumLocalOps(N::Int) = InfiniteSumLocalOps(Vector{ITensor}(undef, N))
## InfiniteSumLocalOps(data::Vector{ITensor}) = InfiniteSumLocalOps(CelledVector(data))
## getindex(l::InfiniteSumLocalOps, n::Integer) = ITensors.data(l)[n]

# TODO? Instead of having a big quasi empty ITensor, store only the non zero blocks

translator(mpo::InfiniteMPO) = mpo.data.translator


text:
function del_met!(metnet::MetNet, met::IDER_TYPE)
    meti = metindex(metnet, met)
    metnet.mets[meti] = EMPTY_SPOT
    z = zero(eltype(metnet.S))
    metnet.S[meti, :] .= z
    metnet.b[meti] = z

    return metnet
end

function del_met!(metnet::MetNet, mets::Vector)
    for met in mets
        del_met!(metnet, met)
    end
    return metnet
end


text:
@testset "DAR Model" begin
    # Univariate
    sunspot_year = CSV.File(StateSpaceModels.SUNSPOTS_YEAR) |> DataFrame

    @test has_fit_methods(DAR)

    model = DAR(sunspot_year.value, 9)
    fit!(model)
    # Runned on Python statsmodels
    @test loglike(model) ≈ -1175.9129 atol = 1e-5 rtol = 1e-5

    # forecasting
    @test_throws ErrorException forecast(model, 10)
end


text:
using GCMAES
using Zygote
using Distributed
using Test

@everywhere rastrigin(x) = 10length(x) + sum(x.^2 .- 10 .* cos.(2π .* x))

D = 200
x0 = fill(0.3, D)
σ0 = 0.2
lo = fill(-5.12, D)
hi = fill(5.12, D)

GCMAES.minimize(rastrigin, x0, σ0, lo, hi, maxiter = 5, autodiff = false)
GCMAES.minimize(rastrigin, x0, σ0, lo, hi, maxiter = 5, autodiff = true)
GCMAES.minimize(rastrigin, x0, σ0, lo, hi, maxiter = 5, constr = BoxConstraint(lo, hi, 2))
GCMAES.minimize(rastrigin, x0, σ0, lo, hi, maxiter = 5, constr = NormConstraint(2))
GCMAES.minimize(rastrigin, x0, σ0, lo, hi, maxiter = 5, trans = BoxLinQuadTransform(lo, hi, scaling = false))
GCMAES.minimize(rastrigin, x0, σ0, lo, hi, maxiter = 5, trans = BoxLinQuadTransform(lo, hi, scaling = true))

rm("CMAES.bson", force = true)


text:
#=
Data:  Yℓ = Cθℓ + ℘ (dCθωℓ_dω) * ω + nℓ, where Cθℓ is LCDM and dCθωℓ_dω are the linearized exotic model
       and ℘ projects out the LCDM fluxtuations. ω ∼ N(0, σ2 I)

Model: Yℓ = Cθℓ + X * β + nℓ the columns of X are from the top LDA modes
	from S_P/(κ S_ + N), S_P =  dCθωℓ_dω * dCθωℓ_dωᵀ. the prior on β is non-informative and looks like
	Put g-prior on β ∼ N(0, g * (Xᵀ(N⁻)X)⁻) with π(g) ∝ ?, which is independent of π(θ).
	(.... maybe β ∼ N(0, g * (Xᵀ(Cθℓ + N)⁻X))⁻) is better instead).
	Use Metropolis Hasitings to sample from the posterior on (β, θ).

Notes:
	We use a subjective prior on ω for constructing the LDA design matrix X,
	but use a non-informative g-prior the posterior on β.
	I think this solves the problem that a non-informative prior for ω
	will not generate meaningful LDA modes (since a non-informative prior on ω make Cθωℓ behave like white noise)

	kill -9 -1 ...this kills all procs from a given user

	run with  include("scripts/MetropolisForLDA.jl")

ToDo:
	- eigs(Symmetric(...)) in place of custom_real_eigs
	- start a new git repo, move source code to src, clean up src, make a ipython notebook

# The order of  LCDM parameters in each file is omega_b*h^2, omega_c*h^2, theta, A_s, n_s, tau.

=#

using LDAforCMB
using Optim, Distributions, PyCall, PyPlot
@pyimport seaborn as sns
@pyimport pandas as pd

# --- set the seed
if false # < --- set to false if you want to use a specific seed
	seedstart = rand(UInt64)
else
	seedstart = 0x6d800cba0725d111
end
srand(seedstart)

const mcmc_runs     = 5_000
const mcmc_thin     = 5
const mcmc_burn     = 500

const exotic_on    = true
const exotic_scale = 0.4
const n_basis      = 5 	       # <----- number of columns of X in Bayesian chain
const kappa        = 1_000.0   # <---- factor to force LDA modes to be approx orthogonal to S_L

const lmax         = 2400
const indexrange   = collect(10:lmax)            # <----  the files start at ell = 0
const ell          = collect(0:5000)[indexrange]
const n            = length(indexrange)
const fwhm         = 20 			              # <---- radius of ell smoothing of dCl_dPa
const lcdm_fid_par = [0.0825,  0.125,   3.218,  0.010413, 0.97,  0.022]
const lcdm_names   = [:tau,  :omch2,  :logA, :theta  , :ns,  :ombh2]
const cl_fid       = trans_to_pico(lcdm_fid_par)[indexrange]

# ----- make noise_var
function var_of_cls(ell, cls)
	Nell = 0.00015048 * exp( ell.^2 * 7.69112e-7 ) .* (ell.^2) / 2 / pi
	fsky = 0.3
	2 * ((cls .+ Nell).^2 ./ (2*ell .+ 1)) / fsky # the` broadcasts Nell
end

# get the cov matrix for lcdm
const lcdm_cov = let
		covpath = joinpath(LDAforCMB.pathtosrc, "chain_data/proposal.covmat")
		propcov = map(Float64, readdlm(covpath))
		propcov_names =
			covpath |>
			open |>
			readline |>
			split |>
			x->convert(Array{ASCIIString},x) |>
			x->x[2:end]

		lcdmi = [1, 9, 12, 14, 15, 16]
		propcov[lcdmi,lcdmi]  # <-- gives the corresponding cov mat
	end
const lcdm_chol = chol(lcdm_cov, Val{:L})

# -----  make S_L, N (LCDM signal cov and noise cov matrix)
const S_L, N, dCl_dtheta =  let
	dx = diagm(1e-6 * ones(6))
	dCl_dtheta = Array(Float64, length(ell), 6)
	for k=1:6
		dCl_dtheta[:,k] =  (trans_to_pico(lcdm_fid_par + dx[:,k]) - trans_to_pico(lcdm_fid_par))[indexrange] /  dx[k,k]
	end
	S_L  = dCl_dtheta  * dCl_dtheta.'
	N    = diagm(var_of_cls(ell,cl_fid))

	S_L, N, dCl_dtheta
end

# -----  make S_P, X (exotic cov matrix, LDA exotic design)
const S_P, S_Pforsim, X, Xfull, eigfull = let
	responses_Pa       = readcsv(joinpath(LDAforCMB.pathtosrc, "covariances/responses.csv"))[indexrange, 30:(end-60)]
	responses_transfer = transpose(readcsv(joinpath(LDAforCMB.pathtosrc, "covariances/responses_transfer.csv")))[indexrange, :]
	lcdm_fid           = readcsv(joinpath(LDAforCMB.pathtosrc, "covariances/lcdm_fid.csv"))[indexrange]

	dCl_dT   = responses_transfer .- lcdm_fid
	dCl_dPa  = responses_Pa .- lcdm_fid
	smoothmat = zeros(n, n)
	for k in -fwhm:fwhm
		smoothmat += diagm(ones(n-abs(k)),k)
	end
	smoothmat = smoothmat ./ sum(smoothmat, 2)
	dCl_dPa   = smoothmat * dCl_dPa

	## !!!!! here I take out the LCDM looking modes
	## dCl_dPa -=  dCl_dtheta * ((dCl_dtheta.' * dCl_dtheta) \ (dCl_dtheta.' * dCl_dPa))

	# ---- generate the LDA vectors

	### option:    random raw modes
	# XX   =  dCl_dPa[:, rand(1:size(dCl_dPa,2), n_ω)]]
	### option:    random raw modes with LCDM projected out
	# XX  =  dCl_dPa[:,rand(1:size(dCl_dPa,2), n_ω)]
	# XX =  XX - dCl_dtheta * ((dCl_dtheta.' * dCl_dtheta) \ (dCl_dtheta.' * XX))
	### option:    LDA modes (more than used for estimation)
	# _, XX = custom_real_eigs((N + kappa * S_L) \ dCl_dPa * dCl_dPa.' , n_ω)
	### option:     PCA then project out dCl_dtheta
	#_, XX = custom_real_eigs(dCl_dPa * dCl_dPa.', n_ω)
	#XX =  XX - dCl_dtheta * ((dCl_dtheta.' * dCl_dtheta) \ (dCl_dtheta.' * XX))
	XX   =  dCl_dPa
	S_P  =  XX * XX.'
	eigfull, Xfull = custom_real_eigs((N + kappa * S_L) \ S_P, n_basis + 50)

	# ---- This is used for the simulations (and is slightly smoother)
	# ---- The hope is that the LDA will automatically look for smoother detectable modes even when
	# ---- generated from a non-smooth prior.
	#
	# ---- This also has the benefit of generating simulations which have a large detectable smooth fluctuation
	# ---- without large visible spikes at low ell.
	nm = size(dCl_dPa, 2)
	Σ = Float64[exp(-abs(i-j)^1.9/25) for i=1:nm, j=1:nm]
	#plot(chol(Σ, Val{:L})*randn(111))
	XXforsim = XX * chol(Σ, Val{:L})  # correlates the primordial top-hat modes
	S_Pforsim = XXforsim * XXforsim.'

	S_P, S_Pforsim, Xfull[:,1:n_basis], Xfull, eigfull
end

# ---- set the prior and a proposal distribution
# const β_prior_RV    = MvTDist(20.0, zeros(n_basis), 10*inv(X.' * ((S_L + N) \ X)))
const β_prior_RV    = MultivariateNormal(zeros(n_basis), abs2(25)*inv(X.' * (N \ X)))
const lcdm_prior_RV = MultivariateNormal(lcdm_fid_par, abs2(1)*lcdm_cov)

# ------  choose a proposal distribution
const β_proposal_RV    = MultivariateNormal(zeros(n_basis), abs2(0.5) * inv(X.' * (N \ X)) )
const lcdm_proposal_RV = MultivariateNormal(zeros(lcdm_fid_par), abs2(0.5) * lcdm_cov)
prop_β(β_curr)	       = rand(β_proposal_RV) + β_curr
prop_lcdm(lcdm_curr)   = rand(lcdm_proposal_RV) + lcdm_curr

# ------ loglikelihood
function loglike(lcdm_pars, β, cl_obs)
	clTT = try trans_to_pico(lcdm_pars)[indexrange] catch zeros(n) end
	clTT += exotic_on ? (X * β) : zeros(n)
	noise_var_new =  var_of_cls(ell, clTT)
	logl = 0.0
	for k in 1:length(indexrange)
		logl += -0.5 * ((cl_obs[k] - clTT[k])^2) / noise_var_new[k]
	end
	logl += logpdf(lcdm_prior_RV, lcdm_pars)
	logl += exotic_on ? logpdf(β_prior_RV, β) : zero(Float64)
	logl
end

# find mle now make an mle for LCDM parameters which uses LDA to avoid exotic directions
let
	# topXfullmodes is the number of modes used for the max_scan statistic
	const topXfullmodes = findfirst(cumsum(eigfull) .> (.99 * sum(eigfull)) )
	#const topXfullmodes = findfirst(eigfull .< 0.1maximum(eigfull))
	const v = nullspace(Xfull[:,1:topXfullmodes].')
	# ------ loglikelihood
	function loglike_maxscan(lcdm_pars, cl_obs)
		clTT          = try trans_to_pico(lcdm_pars)[indexrange] catch zeros(n) end
		noise_var_new =  var_of_cls(ell, clTT)
		vcldiff       = v.' * (clTT - cl_obs)
		Nnew          = diagm(var_of_cls(ell,clTT))
		vNv           = v.' * (Nnew * v)
		L             = chol(vNv, Val{:L})
		logl          = -0.5 * sumabs2(L \ vcldiff)
		logl          += logpdf(lcdm_prior_RV, lcdm_pars)
		logl
	end
	global findmle
	function findmle(cl_obs)
		minimize_me(x) = -loglike_maxscan(lcdm_chol * x, cl_obs)
		initalx  = lcdm_chol \ lcdm_fid_par
		optout   = Optim.optimize(minimize_me, initalx)
		x_mle    = optout.minimum
		lcdm_mle = lcdm_chol * x_mle
		return lcdm_mle
	end
end

function mcmc_chain(cl_obs, numruns::Int64, thin::Int64, burn::Int64, lcdm_start = lcdm_fid_par, β_start = zeros(n_basis))
	lcdmmat      = Array(Float64, 6, numruns)
	βmat         = Array(Float64, n_basis, numruns)
	lcdmmat[:,1] = copy(lcdm_start)
	βmat[:,1] 	 = copy(β_start)
	curr_ll		 = loglike(lcdm_start, β_start, cl_obs)
	for j = 2:numruns
		lcdm_prop = prop_lcdm(lcdmmat[:, j-1])
		β_prop    = prop_β(βmat[:, j-1])
		prop_ll   = loglike(lcdm_prop, β_prop, cl_obs)
		prob_accept = min(1.0, exp(prop_ll - curr_ll))
		if j % 250 == 1
			@show prob_accept
		end
		if rand() < prob_accept
			lcdmmat[:,j] = copy(lcdm_prop)
			βmat[:,j]	 = copy(β_prop)
			curr_ll		 = prop_ll
		else
			lcdmmat[:,j] = copy(lcdmmat[:, j-1])
			βmat[:,j]	 = copy(βmat[:, j-1])
		end
	end
	return lcdmmat[:,burn:thin:end], βmat[:,burn:thin:end]
end

# ------ make data --- nonlinear
function gendata(scalefactor)
	cl_lcdm_pars = rand(MultivariateNormal(lcdm_fid_par, abs2(0.25)*lcdm_cov))
	cl_lcdm      = trans_to_pico(cl_lcdm_pars)[indexrange]
	dp           = samp_cov(S_Pforsim)  # <---- exotic fluctuation
	cl_exotic    = cl_lcdm
	cl_exotic   += exotic_on ? (dp .* scalefactor) : zeros(cl_exotic)
	cl_obs    	 = cl_exotic + √(var_of_cls(ell, cl_exotic)) .* randn(length(ell))
	cl_obs, cl_exotic, cl_lcdm, cl_lcdm_pars
end

cl_obs, cl_exotic, cl_lcdm, cl_lcdm_pars = gendata(exotic_scale)
#=
semilogy(cl_exotic)
semilogy(cl_lcdm)

semilogy(abs(cl_exotic-cl_lcdm))
semilogy(abs(cl_obs -cl_exotic))

plot(cl_exotic-cl_lcdm)
=#

#  ----- run the mcmc
lcdm_mle = findmle(cl_obs)
clmle = trans_to_pico(lcdm_mle)[indexrange]
@time lcdmmat, βmat = mcmc_chain(cl_obs, mcmc_runs, mcmc_thin, mcmc_burn, lcdm_mle)

#=  ----- make some plots of the output

# ----- trace plots for the lcdm parameters
figure(figsize = (14,10))
for cntr = 1:6
	subplot(3,2,cntr)
	plot(lcdmmat[cntr,:].', label="$(lcdm_names[cntr])")
	plot(zeros(lcdmmat[cntr,:].')+cl_lcdm_pars[cntr], color="r", label="true $(lcdm_names[cntr])")
	plot(zeros(lcdmmat[cntr,:].')+lcdm_fid_par[cntr], color="k", label="fid $(lcdm_names[cntr])")

	legend()
end

# -----  trace plots for the βs
figure(figsize = (14,10))
for cntr = 1:min(n_basis,4)
	subplot(2,2,cntr)
	plot(βmat[cntr,:].')
	plot(zeros(βmat[cntr,:].'), color="k", label="fid")
	legend()
end

# -----  pairplots for the lcdm parameters
lcdmdf = pd.DataFrame(lcdmmat.', columns=lcdm_names)
g = sns.PairGrid(lcdmdf)
g[:map_diag](sns.distplot, kde=false, rug=false)
g[:map_lower](sns.kdeplot, n_levels = 3)
g[:map_upper](plt[:scatter], alpha = 0.2)
centroids = pd.DataFrame(cl_lcdm_pars.', columns=lcdm_names)
g[:data] = centroids
g[:map_offdiag](plt[:scatter], s=200, marker=".", color="red")
centroids = pd.DataFrame(lcdm_mle.', columns=lcdm_names)
g[:data] = centroids
g[:map_offdiag](plt[:scatter], s=50, marker="*", color="k")

# -----  pairplots for the βs
βdf = pd.DataFrame(βmat[1:4,:].')
g = sns.PairGrid(βdf)
g[:map_diag](sns.distplot, kde=false, rug=false)
g[:map_lower](sns.kdeplot, n_levels = 3)
g[:map_upper](plt[:scatter], alpha = 0.2)
centroids = pd.DataFrame(zeros(1,n_basis))
g[:data] = centroids
g[:map_offdiag](plt[:scatter], s=200, marker=".", color="red")
# centroids = pd.DataFrame(β_mle.')
# g[:data] = centroids
# g[:map_offdiag](plt[:scatter], s=50, marker="*", color="k")

#  ----- best fit
lcdm_fit = vec(mean(lcdmmat,2))
β_fit = vec(mean(βmat,2))
cl_fit = trans_to_pico(lcdm_fit)[indexrange] +  X * β_fit
noise_level = sqrt(var_of_cls(ell, cl_fit))
plot((cl_exotic-cl_lcdm), label = "exotic fluctuation")
plot((cl_exotic-cl_fit),  label = "best fit")
plot(noise_level,  label = "noise level")
#plot(cl_exotic,  label = "true model")
legend()
=#

# ------ poterior PTE for β

figure()
const βΣ = cov(βmat.')
const βm = vec(mean(βmat,2))
β_functional(β) = (β.-βm).' * (βΣ \  (β.-βm))
post_β_functional = mapslices(β_functional, βmat, 1)
pte = sum(post_β_functional .>= β_functional(zeros(βm)))/size(βmat,2)
sns.distplot(post_β_functional, label = "posterior distance from mean")
sns.rugplot(β_functional(zeros(βm)), color = "red", linewidth = 2, label = "distance of 0 from mean, pte = $(round(pte,4))")
legend()

# ---- max scan test:

function sim_maxscanChi2(Xfull, clmle, n_sims)
	sims = zeros(n_sims)
	sims_maxat = zeros(n_sims)
	for k = 1:n_sims
		noisesim = √(var_of_cls(ell, clmle)) .* randn(length(ell))
		sims[k], sims_maxat[k] = maxscanChi2(Xfull, clmle + noisesim, clmle)
	end
	return sims, sims_maxat
end
function sim_redularChi2(Xfull,  clmle, n_sims)
	sims = zeros(n_sims)
	for k = 1:n_sims
		noisesim = √(var_of_cls(ell, clmle)) .* randn(length(ell))
		sims[k] = redularChi2(Xfull, clmle + noisesim, clmle)
	end
	return sims
end

function maxscanChi2(Xfull, cl_obs, clmle)
	Xclmle = Xfull.' * clmle
	Xcl_obs = Xfull.' * cl_obs
	Nmle = diagm(var_of_cls(ell, clmle))
	XNX = Xfull.' * (Nmle * Xfull)
	teststat, maxat = 0.0, 0
	for modes = 1:size(XNX,2)
		tmp = reducdedChi2(Xclmle[1:modes], Xcl_obs[1:modes], XNX[1:modes, 1:modes])
		if tmp >= teststat
			teststat = tmp
			maxat = modes
		end
	end
	return teststat, maxat
end

function redularChi2(Xfull, cl_obs, clmle)
	Xclmle = Xfull.' * clmle
	Xcl_obs = Xfull.' * cl_obs
	Nmle = diagm(var_of_cls(ell, clmle))
	XNX = Xfull.' * (Nmle * Xfull)
	return reducdedChi2(Xclmle, Xcl_obs, XNX)
end

function reducdedChi2(Xclmle, Xcl_obs, XNX)
	chival =  ((Xcl_obs - Xclmle).' * (XNX \ (Xcl_obs - Xclmle)))[1]
	# ---- option 1
	chival -= size(XNX,1)
	chival /= √(2*size(XNX,1))
	# ---- option 2
	# chival = chi2_to_normal(chival, size(XNX, 1))
	return chival
end

function chi2_to_normal(x, df)
	cdfx = logcdf(Chisq(df), x)
	invcdfcdfx = invlogcdf(Normal(0, 1), cdfx)
	return invcdfcdfx
end

number_of_modes_maxing = findfirst(cumsum(eigfull) .> (0.99 * sum(eigfull)) )
# number_of_modes_maxing = findfirst(eigfull .< 0.01maximum(eigfull))
# number_of_modes_maxing = 50
obs_test_stat, obs_maxat = maxscanChi2(Xfull[:,1:number_of_modes_maxing], cl_obs, clmle)
sim_test_stat, sim_maxat = sim_maxscanChi2(Xfull[:,1:number_of_modes_maxing], clmle, 5000)
pte = sum(sim_test_stat .>= obs_test_stat)/length(sim_test_stat)

figure()
sns.distplot(sim_test_stat, kde = false, label = "null model")
sns.rugplot([obs_test_stat], color = "red", linewidth = 2, label = "obs, pte = $(round(pte,4))")
legend()

#=
figure()
sns.distplot(sim_maxat, kde = false, label = "where max occurs in null",kde=false)
sns.rugplot([obs_maxat], color = "red", linewidth = 2, label = "obs max")
legend()

=#

obs_test_statr = redularChi2(Xfull[:,1:number_of_modes_maxing], cl_obs, clmle)
sim_test_statr = sim_redularChi2(Xfull[:,1:number_of_modes_maxing], clmle, 5000)
pte = sum(sim_test_statr .>= obs_test_statr) / length(sim_test_statr)

figure()
sns.distplot(sim_test_statr, kde = false, label = "null model")
sns.rugplot([obs_test_statr], color = "red", linewidth = 2, label = "obs, pte = $(round(pte,4))")
legend()

max_stars_repo_name: EthanAnderes/LDAforCMB

max_issues_repo_name: EthanAnderes/LDAforCMB

max_forks_repo_name: EthanAnderes/LDAforCMB


text:
using Test
using CompatDevTools

function simulate_keystrokes(keys...)
    keydict =  Dict(:up => "\e[A",
                    :down => "\e[B",
                    :enter => "\r")
    for key in keys
        if isa(key, Symbol)
            write(stdin.buffer, keydict[key])
        elseif isa(key, Char)
            write(stdin.buffer, "$key")
        else
            write(stdin.buffer, key)
        end
    end
end

function generate_env(filename::S, pkg_versions::Tuple{S,S,S}) where {S <: String}
    mkpath(dirname(filename))
    open(filename, "w") do io
        println(io, "[deps]")
        println(io, "SomePkgA = \"commitSHAA\"")
        println(io, "SomePkgB = \"commitSHAB\"")
        println(io, "SomePkgC = \"commitSHAC\"")
        println(io, "")
        println(io, "[compat]")
        println(io, "SomePkgA = \"$(pkg_versions[1])\"")
        println(io, "SomePkgB = \"$(pkg_versions[2])\"")
        println(io, "SomePkgC = \"$(pkg_versions[3])\"")
    end
end

@testset "CompatDevTools Unit tests - changes" begin
    mktempdir() do path
        tomlPkg = joinpath(path, "Pkg.jl", "Project.toml")
        tomlA   = joinpath(path, "Pkg.jl", "envA", "Project.toml")
        tomlB   = joinpath(path, "Pkg.jl", "envB", "Project.toml")
        generate_env(tomlPkg, ("1.2.3", "4.5.6", "7.8.9"))
        generate_env(tomlA,   ("1.2.4", "4.6.6", "8.8.9"))
        generate_env(tomlB,   ("1.3.3", "4.6.6", "7.8.9"))

        simulate_keystrokes(:down, :down, :enter, 'd')
        simulate_keystrokes(:down, :down, :enter, 'd')
        simulate_keystrokes(:down, :down, :enter, 'd')
        CompatDevTools.synchronize_compats(path)

        contents = join(readlines(tomlPkg; keep=true))
        @test contents == "[deps]\nSomePkgA = \"commitSHAA\"\nSomePkgB = \"commitSHAB\"\nSomePkgC = \"commitSHAC\"\n\n[compat]\nSomePkgA = \"1.3.3\"\nSomePkgB = \"4.6.6\"\nSomePkgC = \"7.8.9\"\n"
    end
end

@testset "CompatDevTools Unit tests - compat kickstart" begin
    @test CompatDevTools.compat_kick_start(dirname(@__DIR__)) == "[compat]\nOrderedCollections = \"1.4.1\"\njulia = \"1\""
end


text:
"""
`Rotate270 <: ImageTransformation`

Description
============

Rotates the image upwards 270 degrees, which can also be described
as rotating the image downwards 90 degrees. This is a special case
rotation, because it can be performed very efficiently by simply
rearranging the existing pixels. However, it is generally not the
case that the output image will have the same size as the input
image, which is something to be aware of.

If created using the parameter `chance`, the operation will be
lifted into a `Either`. See the documentation of
`Either` for more information.

Usage
======

    Rotate270()

    Rotate270(chance)

Arguments
==========

- **`chance`** : The probability of the operation to be applied to
the image if the operation is used in `transform`. Must be within
the set [0, 1].

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Images
    using Augmentor

    # Create a test image with 4 pixel.
    img = grayim(UInt8[200 150; 50 1])

    # Apply the image operation to the image
    img_new = transform(Rotate270(), img)

    # Check manually if the result is as expected
    @assert img_new == grayim(UInt8[150 1; 200 50])

see also
=========

`ImageTransformation`, `Either`, `transform`
"""
immutable Rotate270 <: ImageTransformation
end

Rotate270(chance) = ProbableOperation(Rotate270(); chance = chance)

Base.showcompact(io::IO, op::Rotate270) = print(io, "Rotate 270 degrees.")
multiplier(::Rotate270) = 1

function transform{T<:AbstractImage}(op::Rotate270, img::T)
    result = copyproperties(img, rotate_expand(img, -π/2))::T
    _log_operation!(result, op)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_reader

"""
    vtk_reader(;kwargs...)
    vtk_reader(children::Any;kwargs...)
    vtk_reader(children_maker::Function;kwargs...)

A Reader component.
Reader is exposing a reader to a downstream filter
It takes the following set of properties:
  - vtkClass: vtk.js reader class name
  - url: string
  - parseAsText: string
  - parseAsArrayBuffer: base64String
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `port` (Real; optional): downstream connection port
- `vtkClass` (String; optional): vtkClass name
- `url` (String; optional): set of url to fetch data from
- `parseAsText` (String; optional): set text data to process
- `parseAsArrayBuffer` (String; optional): set binary data to process from base64 string
- `renderOnUpdate` (Bool; optional): Automatically render on data loaded
- `resetCameraOnUpdate` (Bool; optional): Automatically reset camera on data loaded
"""
function vtk_reader(; kwargs...)
        available_props = Symbol[:children, :id, :port, :vtkClass, :url, :parseAsText, :parseAsArrayBuffer, :renderOnUpdate, :resetCameraOnUpdate]
        wild_props = Symbol[]
        return Component("vtk_reader", "Reader", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_reader(children::Any; kwargs...) = vtk_reader(;kwargs..., children = children)
vtk_reader(children_maker::Function; kwargs...) = vtk_reader(children_maker(); kwargs...)

max_stars_repo_name: ericbdaniels/dash-vtk

max_issues_repo_name: ericbdaniels/dash-vtk

max_forks_repo_name: ericbdaniels/dash-vtk


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

mutable struct IoK8sApiStorageV1CSINodeDriver <: SwaggerModel
    allocatable::Any # spec type: Union{ Nothing, IoK8sApiStorageV1VolumeNodeResources } # spec name: allocatable
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    nodeID::Any # spec type: Union{ Nothing, String } # spec name: nodeID
    topologyKeys::Any # spec type: Union{ Nothing, Vector{String} } # spec name: topologyKeys

    function IoK8sApiStorageV1CSINodeDriver(;allocatable=nothing, name=nothing, nodeID=nothing, topologyKeys=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("allocatable"), allocatable)
        setfield!(o, Symbol("allocatable"), allocatable)
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("nodeID"), nodeID)
        setfield!(o, Symbol("nodeID"), nodeID)
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("topologyKeys"), topologyKeys)
        setfield!(o, Symbol("topologyKeys"), topologyKeys)
        o
    end
end # type IoK8sApiStorageV1CSINodeDriver

const _property_map_IoK8sApiStorageV1CSINodeDriver = Dict{Symbol,Symbol}(Symbol("allocatable")=>Symbol("allocatable"), Symbol("name")=>Symbol("name"), Symbol("nodeID")=>Symbol("nodeID"), Symbol("topologyKeys")=>Symbol("topologyKeys"))
const _property_types_IoK8sApiStorageV1CSINodeDriver = Dict{Symbol,String}(Symbol("allocatable")=>"IoK8sApiStorageV1VolumeNodeResources", Symbol("name")=>"String", Symbol("nodeID")=>"String", Symbol("topologyKeys")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiStorageV1CSINodeDriver }) = collect(keys(_property_map_IoK8sApiStorageV1CSINodeDriver))
Swagger.property_type(::Type{ IoK8sApiStorageV1CSINodeDriver }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1CSINodeDriver[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1CSINodeDriver }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1CSINodeDriver[property_name]

function check_required(o::IoK8sApiStorageV1CSINodeDriver)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("nodeID")) === nothing) && (return false)
    true
end

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


text:
module StandardMap
export standard_map
using ..ExampleBase: LEDemo, DiscreteExample

@inline function phase_dynamics!(u_next, u, k, t)
    u_next[2] = (u[2] + k * sin(u[1])) % 2π
    u_next[1] = (u[1] + u_next[2]) % 2π
end

@inline @views function tangent_dynamics!(u_next, u, k, t)
    phase_dynamics!(u_next[:, 1], u[:, 1], k, t)
    u_next[2, 2:end] .= u[2, 2:end] .+ k * cos(u[1, 1]) .* u[1, 2:end]
    u_next[1, 2:end] .= u[1, 2:end] .+ u_next[2, 2:end]
end

"""
Return a [`LEDemo`](@ref) for the Chirikov standard map.

* B. V. Chirikov, Physics Reports 52, 263-379 (1979)
* <http://sprott.physics.wisc.edu/chaos/comchaos.htm>
* <https://en.wikipedia.org/wiki/Standard_map>
* <http://www.scholarpedia.org/article/Chirikov_standard_map>
"""
function standard_map(;
        u0=[2.68156, 2.31167],
        t_renorm=10,
        t_attr=1000000,
        atol=0, rtol=0.01,
        terminator_options = [:atol => atol, :rtol => rtol,
                              :max_tail_corr => rtol,
                              :first_check => 10000],
        kwargs...)
    # TODO: Improve the accuracy. Check the paper.  It looks like
    # `t_attr=1000000` is required to see some kind of convergence.
    k = 1
    LEDemo(DiscreteExample(
        "Chirikov standard map",
        phase_dynamics!, u0, t_renorm, k,
        tangent_dynamics!,
        t_attr,
        [0.10497, -0.10497],   # known_exponents
        atol, rtol,
        terminator_options,
    ); kwargs...)
end

end


text:
#hello

#okay, this is a very quick test of a basic blockchain in julia
# for educational purpoess probably

#first we set up our types

abstract type CryptoPrimitive end

abstract type Hash <: CryptoPrimitive end
abstract type MAC <: CryptoPrimitive end
abstract type Signature <: CryptoPrimitive end

abstract type AbstractBlockChain end
abstract type AbstractBlock end

abstract type AbstractGenesisBlock <: AbstractBlock end

abstract type AbstractBlockData end
abstract type AbstractMessage <:AbstractBlockData end

primitive type 256Hash <: Hash 256 end
primitive type 512Hash <: Hash 512 end

primitive type 256MAC <: MAC 256 end
primitive type 512MAC <: MAC 512 end

primitive type 256Signature <: Signature end
primitive type 512Signature <: Signature end

#now for our more complex types

const genesis_message = "This is how the world begins: not with a whimper, but with a bang!"

type Blockchain <: AbstractBlockChain
	blocks:: Array{AbstractBlock}
	timeCreated:: Date
	Description:: AbstractString
	id:: Hash
	sig:: Signature
	MAC:: MAC
end

type Block <: AbstractBlock
	index:: Integer
	hash:: Hash
	previousHash:: Hash
	timestamp:: Date
	data:: Array{Message}
	difficulty:: Integer
	nonce:: Integer

	function Block(index, hash, previousHash, timestamp, data, difficulty, nonce)
		return new(index, hash, previousHash, timestamp, data, difficulty,nonce)
	end

end

type GenesisBlock <: AbstractGenesisBlock
	index:: Integer
	hash:: Hash
	previousHash:: Hash
	timestamp:: Date
	data::Array{Message}
	difficulty:: Integer
	nonce:: Integer

	function GenesisBlock()
		const genesis_message = getGenesisMessage()
		index = 0
		data = Message(0, genesis_message)
		difficulty = 0
		nonce = 0
		previousHash = Hash(0)
		timestamp = 1 # pass
		hash=  caluclate_hash(index, data[0].data, difficulty, nonce, previousHash, timestamp)
		return new(index, hash, previousHash, timestamp, data, difficulty, nonce)
	end
end

type Message <: AbstractMessage
	id:: AbstractString
	data:: AbstractString
	MAC:: MAC

	function Message(id:: Integer, data::AbstactString)
		mac = calculate_MAC(id, data)
		return new(id, data, mac)
	end
	function Message(id:: Integer, data::AbstractString, encrypt_func:: function) {
		encrypt_data = encrypt_func(data)
		mac = calculate_MAC(id, encrypt_data)
		return new(id, encrypt_data, mac)
	end

	function Message(data)
		id = calculate_id()
		mac = calculate_MAC(id, data)
		return new(id, data, mac)
	end
	function Message(data, encrypt_func:: function)
		id = calculate_id()
		encrypt_data = encrypt_func(data)
		mac = calculate_MAC(id, encrypt_data)
		return new(id, encrypt_data, mac)
	end

end

function calculate_hash(inputs)
	# so originally we're going to d something here
	# but I don't know how
	#either use a crypto implementation or roll our own
	# so now we're just going to return
	return inputs
	end

function calculate_MAC(id, data)
	#no idea how to implement this either
	return MAC()
end

function getGenesisMessage()
	return genesis_message
end
# okay, now we do operations on the block chain and functoins thereby

function getLatestBlock(blockchain:: AbstractBlockChain)
	return blockchain.blocks[-1]
end

function getBlocks(blockchain:: AbstractBlockChain)
	return blockchain.blocks
end

function generateNextBlock(blockchain:: AbstractBlockChain, blockData:: Array{Message})
	prevBlock = getLatestBlock(blockchain)
	nextIndex = prevBlock.index +1
	timestamp = getDateTime()
	nextHash = calculate_hash((nextIndex, prevBlock.hash,timestamp, blockData))
	difficulty = getDifficulty()
	nonce = getNonce()
	newBlock = Block(index, hash, previousHash, timestamp, blockData, difficulty,nonce)
	return newBlock
end

function getDateTime()
	Error("not implemented yet")
end

function getDifficulty()
	Error("not implemented yet")
end

function getNonce()
	Error("not implemented yet")
end

function isValidBlock(block:: AbstractBlock, prevBlock::AbstractBlock)
	if prevBlock.index +1 != block.index
		Error('New block index does not match')
		return false
	end
	if block.previousHash != prevBlock.hash
		Error('Block previous hash does not match')
		return false
	end
	if calculateHash(block) != block.hash
		Error('Block hash invalid')
		return false
	end
	if !isValidTimestamp(block.timestamp, prevBlock.timestamp)
		Error('Block timestamp is invalid')
		return false
	end
	return true
end

function isValidGenesis(proposedGenesis:: AbstractGenesisBlock)
	if proposedGenesis.index != 0
		Error('Genesis does not start at 0')
		return false
	end
	if proposedGenesis.hash != calculate_hash(getGenesisMessage())
		Error('Genesis block hash invalid')
		return false
	end
	if proposedGenesis.previousHash != Hash(0)
		Error('Genesis block previous hash invalid')
		return false
	end
	#we might need to validate timestamp, but not sure how
	return true
end

function validateBlockChainBlocks(blocks:: Array{AbstractBlock})
	len = lengh(blocks)
	if !isValidGenesis(blocks[0])
		return false
	end
	for i in 1:len
		if !isValidBlock(blocks[i], blocks[i-1])
			return false
		end
	end
	return true
end

function validateBlockChainProperties(blockchain:: AbstractBlockChain)
	#need to implement properly
	return true
end

function validateBlockChain(blockchain:: AbstractBlockChain)
	return validateBlockChainProperties(blockchain) && validateBlockChainBlocks(blockchain.blocks)
end

function addBlockToChain(blockchain:: AbstractBlockChain, block:: AbstractBlock)
	#do we assume chain is valid here - we might as well
	if !isValidBlock(block, getLatestBlock(blockchain))
		Error('Tried to add invalid block')
		return false
	end
	push!(blockchain.blocks, block)
	# we may need to recalibrate the MAC or something of the chain
	# if so do that here!
	return true
end

function calculateBlockchainSignature(blockchain<: AbstractBlockChain)
	#we're not sure how to implement that AbstractMessage
	# so instead we're just filling this with a placeholder
	return blockchain.Signature
end

# get the cumulative difficulty of a blockchain to see if it works!
function getCumulativeDifficulty = function(blockchain:: AbstractBlockChain)
	blocks = blockchain.blocks
	if length(blocks) <=0
		Error('Blockchain has no length!')
	end
	total = 0
	for block in blocks
		hash = block.hash
		diff = block.difficulty
		if hashMatchesDifficulty(hash, diff)
			total += 2**diff
		end
	else
		Error('not matching hash/difficulty found')
		return -1
	end
	return total
end

function getDifficulty(blockchain:: AbstractBlockChain, difficulty_adjustment_interval, block_generation_interval)
	latestBlock = getLatestBlock(blockchain)
	if latestBlock.index % difficulty_adjustment_interval ==0 && latestBlock.index>=0
		return getAdjustedDifficulty(latestBlock, blockchain,difficulty_adjustment_interval,block_generation_interval)
	end
	return latestBlock.difficulty
end

function getAdjustedDifficulty(latestBlock:: AbstractBlock, blockchain::AbstractBlockChain, difficulty_adjustment_interval, block_generation_interval)
	const blockTimeAdjustmentFactor = 2
	blocks = blockchain.blocks
	prevAdjustmentBlock = blocks[length(blocks)-difficulty_adjustment_interval]
	timeExpected = block_generation_interval* difficulty_adjustment_interval
	if isValidTimestamp(latestBlock) && isValidTimestamp(prevAdjustmentBlock)
		timeTaken = latestBlock.timestamp - prevAdjustmentBlock.timestamp
		if timeTaken < timeExpected /blockTimeAdjustmentFactor
			return prevAdjustmentBlock.difficulty +1
		end
		if timeTaken > timeExpected * blockTimeAdjustmentFactor
			return prevAdjustmentBlock.difficulty -1
		end
		return prevAdjustmentBlock.difficulty
	Error('One of the block timestamps is not valid')
end

max_stars_repo_name: Bmillidgework/Misc-Maths

max_issues_repo_name: Bmillidgework/Misc-Maths

max_forks_repo_name: Bmillidgework/Misc-Maths


text:
using Distributions
using Turing
using Test

# Define model
@model ad_test2(xs) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt.(s))
  xs[1] ~ Normal(m, sqrt.(s))
  xs[2] ~ Normal(m, sqrt.(s))
  s, m
end

# Run HMC with chunk_size=1
chain = sample(ad_test2([1.5 2.0]), HMC(300, 0.1, 1))


text:
module Day08

using LinearAlgebra: dot

export day08

const TESTINPUT = """be cfbegad cbdgef fgaecd cgeb fdcge agebfd fecdb fabcd edb | fdgacbe cefdb cefbgd gcbe
edbfga begcd cbg gc gcadebf fbgde acbgfd abcde gfcbed gfec | fcgedb cgb dgebacf gc
fgaebd cg bdaec gdafb agbcfd gdcbef bgcad gfac gcb cdgabef | cg cg fdcagb cbg
fbegcd cbd adcefb dageb afcb bc aefdc ecdab fgdeca fcdbega | efabcd cedba gadfec cb
aecbfdg fbg gf bafeg dbefa fcge gcbea fcaegb dgceab fcbdga | gecf egdcabf bgf bfgea
fgeab ca afcebg bdacfeg cfaedg gcfdb baec bfadeg bafgc acf | gebdcfa ecba ca fadegcb
dbcfg fgd bdegcaf fgec aegbdf ecdfab fbedc dacgb gdcebf gf | cefg dcbef fcge gbcadfe
bdfegc cbegaf gecbf dfcage bdacg ed bedf ced adcbefg gebcd | ed bcgafe cdgba cbgef
egadfb cdbfeg cegd fecab cgb gbdefca cg fgcdab egfdb bfceg | gbdfcae bgc cg cgb
gcafb gcf dcaebfg ecagb gf abcdeg gaef cafbge fdbac fegbdc | fgae cfgab fg bagce"""

function parseLine(line)
    combinations, output = split.(split(line, "|"))
    keys = _solve(Set.(combinations))
    return getindex.(Ref(keys),Set.(output))
end

function _solve(combinations)
    find(pred, list) = list[findfirst(pred, list)]
    haslength(L) = x -> length(x) == L
    subsetof(s1) = s2 -> issubset(s2, s1)
    contains(s1) = s2 -> issubset(s1, s2)
    unequalto(s1) = s2 -> s1 !== s2 # can check identity in this case
    and(preds...) = x -> all(pred(x) for pred in preds)

    one = find(haslength(2), combinations)
    seven = find(haslength(3), combinations)
    four = find(haslength(4), combinations)
    eight = find(haslength(7), combinations)
    three = find(and(haslength(5), contains(one)), combinations)
    nine = find(and(haslength(6), contains(three)), combinations)
    five = find(and(haslength(5), subsetof(nine), unequalto(three)), combinations)
    two = find(and(haslength(5), unequalto(five), unequalto(three)), combinations)
    six = find(and(haslength(6), contains(five), unequalto(nine)), combinations)
    zero = find(and(haslength(6), unequalto(six), unequalto(nine)), combinations)
    return Dict([zero=>0, one=>1, two=>2, three=>3, four=>4, five=>5, six=>6, seven=>7, eight=>8, nine=>9])
end

function readInput(inputfile="inputs/input08.txt")
    linedata = [parseLine(line) for line in eachline(inputfile)]
    return linedata # transpose
end

function day08(inputdata = readInput("inputs/input08.txt"))
    return (day08Part1(inputdata), day08Part2(inputdata))
end

function day08Part1(outputs)
    # count 1,4,7,8
    count = 0
    for line in outputs
        for out in line
            if out==1 || out==4 || out==7 || out==8
                count += 1
            end
        end
    end
    return count
end

function day08Part2(outputs)
    magnitudes = collect(10 .^ (3:-1:0))
    count = 0
    for line in outputs
        count += dot(line, magnitudes)
    end
    return count
end

end #module


text:
export spine, components_complement, ronkin_polynomial, hypersurface, amoeba_approximation

function Spine2D(polynomial::AbstractTropicalPolynomial, ccs::Vector{<:ComponentComplement}, aa::Bitmap2D)
    Spine2D(polynomial, TropicalCurve(polynomial), ccs, aa)
end

Spine2D(f::MP.AbstractPolynomial; kwargs...) = spine(f; kwargs...)

"""
    spine(f::MP.AbstractPolynomial; options...)

Compute the spine of the amoeba ``\\mathcal{A}(f)``. This algorithm computes first an approximation of
``\\mathcal{A}(f)`` and from this the spine. Returns a Spine2D.

## Example
```julia
@polyvar x y
# use all the defaults
spine(x^2 + y^2 + 1)

# Maybe we think that the compleement of the amoeba has some very small components
# and we want to use an explicit domain
spine(x^2 + y^2 + 1, domain=(-5, 5, -5, 5), minimal_component_size=0.001)
```

Optional arguments:
* `minimal_component_size=0.01`: A guarantee that in each component of the complement of ``\\mathcal{A}(f)`` fits a ball with this diameter.
If this does not hold the algorithm can return a wrong result.
* `domain`: A tuple in the form `(xmin, xmax, ymin, ymax)` which defines a section Ω
for which the amoeba ``\\mathcal{A}(f)`` is computed. This domain has to be such that the intersection
Ω ∩ ``\\mathcal{A}(f)`` still captures the correct topology of ``\\mathcal{A}(f)``.
* `grid`: Based on `minimal_component_size` and `domain` a grid can be computed automatically.
This can also be overwritten with this option.
* `membership_options=[MembershipTestOptions](@ref)()`: Options for the membership test.
* `nsamples=1024`: For the computation we numerically evaluate intergrals. This is the number
of sample points used in the computation.

    spine(f::MP.AbstractPolynomial, A::Bitmap2D; nsamples=1024)

Compute the spine based on the given amoeba approximation `A`.
"""
function spine(f::MP.AbstractPolynomial;
    domain=amoeba_carcase_domain_heuristic(f),
    minimal_component_size=0.01,
    grid = spine_grid(domain..., minimal_component_size),
    membership_options=MembershipTestOptions(),
    kwargs...)
    A = amoeba(f, alg=Greedy(), grid=grid, membership_options=membership_options)
    spine(f, A; kwargs...)
end

function spine(f::MP.AbstractPolynomial, A::Bitmap2D; nsamples=1024)
    p = SP.Polynomial((1.0+0.0im)*f)
    cc, _ = components_complement(A, AmoebaFiber2D(f), p, nsamples=nsamples)

    ronkincoeffs = ronkincoefficients(p, cc)
    ronkinpoly = ronkinpolynomial(f, ronkincoeffs, cc)
    Spine2D(ronkinpoly, cc, A)
end

function spine_grid(f, δ; kwargs...)
    spine_grid(amoeba_carcase_domain_heuristic(f; kwargs...)..., δ)
end
function spine_grid(xmin, xmax, ymin, ymax, δ)
    rx = ceil(Int, (xmax-xmin) / √(0.5δ))
    ry = ceil(Int, (ymax-ymin) / √(0.5δ))
    Grid2D(xmin, xmax, ymin, ymax, rx, ry)
end

"""
    components_complement(S::Spine2D)

Get the ComponentComplement of the spine `S`.
"""
components_complement(S::Spine2D) = S.components_complement

"""
    ronkin_polynomial(S::Spine2D)

Get the tropical polynomial which defines this hypersurface.
"""
ronkin_polynomial(S::Spine2D) = S.ronkin_polynomial

"""
    hypersurface(S::Spine2D)

Get the tropical curve which describes the spine.
"""
hypersurface(S::Spine2D) = S.curve

"""
    amoeba_approximation(S::Spine2D)

Get the approximated amoeba used to compute the spine.
"""
amoeba_approximation(S::Spine2D) = S.amoeba_approximation

function Base.show(io::IO, mime::MIME"text/html", S::Spine2D)
    if inIJulia() && plotsdefined()
        show(io, mime, Main.Plots.plot(S))
    else
        show(io, S)
    end
end
function Base.show(io::IO, S::Spine2D)
    nbounded = length(filter(isbounded, components_complement(S)))
    println(io, typeof(S), " with:")
    println(io, " * ", length(vertices(S.curve)), " vertices")
    println(io, " * ", length(halfrays(S.curve)), " half-rays")
    print(io, " * ", nbounded, " bounded cells")
end

@recipe function plot(S::Spine2D)
    xmin, xmax, ymin, ymax = limits(S.amoeba_approximation.grid)
    xlims --> (xmin, xmax)
    ylims --> (ymin, ymax)
    @series begin
        S.curve
    end
end


text:
# -log(density)
using LinearAlgebra

function normal_(x::Vector{Float64}, μ::Vector{Float64}, Σ::Matrix{Float64})
    # TODO allow Σ a parameter
    d = x - μ
    return 0.5 * (d' * (Σ \ d))
end

function normal_(x::Real, μ::Real, σ::Real)
    # TODO allow σ a parameter
    d = x - μ
    return 0.5 * (d' * d) / (σ * σ)
end

# TODO figure out how to incorporate reparameterizations? (is that what we should call these)
function t_(nu::Float64, s::Float64)
    return function(x)
        n = 0.5 * x[1] * x[1]
        ig = 0.5 * nu * x[2] + 0.5 * nu * s * s * exp(-x[2])
        return n + ig
    end
end

function normal_precision_(x::Vector{Float64}, μ::Vector{Float64}, Ω::Matrix{Float64})
    d = x - μ
    return 0.5 * (d' * Ω * d)
end

function normal_precision_(x::Real, μ::Real, ω::Real)
    d = x - μ
    return 0.5 * d * d * ω
end

## below needs lots of reviewing
# logistic
function logistic(x)
    return 1 / (1 + exp(-x))
end

function genlogistic(N = 1001, K = 3, a = Normal(0, 10), b = Normal(0, 2.5), x = Normal(0, 10))
    α = rand(a)
    β = rand(b, K)
    θ = [α; β]
    X = [ones(N) rand(x, (N, K))]
    μ = X * θ
    p = logistic.(μ)
    y = convert(Vector{Int}, rand.(Bernoulli.(p)))
    return y, X, θ
end

# df = DataFrame(y=Y, x2 = X[:, 2], x3 = X[:, 3], x4 = X[:, 4])
# fit = glm(@formula(y ~ x2 + x3 + x4), df, Bernoulli())

function logisticmodel(y, X, θ)
    μ = X * θ
    β = θ[2:end-1]
    return -y' * μ + sum(log1p.(exp.(μ))) + 0.5 + θ[1] * θ[1] / 100 + 0.5 * β' * β / 10
end
# logisticreg(θ) = logisticmodel(y, X, θ)

function hierlogisticmodel(y, X, θ)
    β = θ[1:(end - 1)]
    K = length(β)
    μ = X * β
    likelihood =
    σ = θ[end]
    κ = exp(σ)
    prior = -0.5 * β' * β / κ ^ 2 - K * σ - 0.1 * κ + σ
    return likelihood + prior
end

heirlogisticreg(θ) = hierlogisticmodel(θ, y, X)

# linear
function genlinear(N = 1001, K = 3, a = Normal(0, 10), b = Normal(0, 2.5), x = Normal(0, 10))
    α = rand(a)
    β = rand(b, K)
    θ = [α; β]
    X = [ones(N) rand(x, (N, K))]
    μ = X * θ
    y = rand.(Normal.(μ, 1))
    return y, X, θ
end

# df = DataFrame(y=Y, x2 = X[:, 2], x3 = X[:, 3], x4 = X[:, 4])
# fit = glm(@formula(y ~ x2 + x3 + x4), df, Normal())

function linearmodel(y, X, θ)
    e = y - X * θ
    return 0.5 * e' * e + 0.5 * θ' * θ / 25
end
# linearreg(θ) = linearmodel(y, X, θ)

function hierlinearmodel(y, X, θ)
    β = θ[1:(end - 1)]
    K = length(β)
    σ = θ[end]
    e = y - X * β
    return 0.5 * e' * e + K*σ + 0.5 * θ' * θ / exp(2 * σ) + 0.5 * exp(2 * σ) / 25 - σ
end

max_stars_repo_name: roualdes/yastan

max_issues_repo_name: roualdes/yastan

max_forks_repo_name: roualdes/yastan


text:
# Define an abstract type to represent color difference metrics
@compat abstract type DifferenceMetric end

# TODO?: make the DifferenMetrics parametric, to preserve type-stability

# CIE Delta E 2000 recommendation
immutable DE_2000 <: DifferenceMetric
    kl::Float64
    kc::Float64
    kh::Float64
    DE_2000(kl,kc,kh) = new(kl,kc,kh)
    DE_2000() = new(1,1,1)
end

# CIE Delta E 94 recommendation
immutable DE_94 <: DifferenceMetric
    kl::Float64
    kc::Float64
    kh::Float64
    DE_94(kl,kc,kh) = new(kl,kc,kh)
    DE_94() = new(1,1,1)
end

# McDonald "JP Coates Thread Company" formulation
immutable DE_JPC79 <: DifferenceMetric

end

# CMC recommendation
immutable DE_CMC <: DifferenceMetric
    kl::Float64
    kc::Float64
    DE_CMC(kl,kc) = new(kl,kc)
    DE_CMC() = new(1,1)
end

# BFD recommendation
immutable DE_BFD <: DifferenceMetric
    wp::XYZ{Float64}
    kl::Float64
    kc::Float64
    DE_BFD(wp,kl,kc) = new(wp,kl,kc)
    DE_BFD() = new(WP_DEFAULT,1,1)
    DE_BFD(kl, kc) = new(WP_DEFAULT,kl, kc)
end

# The original CIE Delta E equation (Euclidian)
immutable DE_AB <: DifferenceMetric

end

# DIN99 color difference (Euclidian)
immutable DE_DIN99 <: DifferenceMetric

end

# DIN99d color difference (Euclidian)
immutable DE_DIN99d <: DifferenceMetric

end

# DIN99o color difference (Euclidian)
immutable DE_DIN99o <: DifferenceMetric

end

# Compute the mean of two hue angles
function mean_hue(h1, h2)
    if abs(h2 - h1) > 180
        if h1 + h2 < 360
            mh = (h1 + h2 + 360) / 2
        else
            mh = (h1 + h2 - 360) / 2
        end
    else
        mh = (h1 + h2) / 2
    end

    mh
end

# Color difference metrics
# ------------------------

# Evaluate the CIEDE2000 color difference formula, implemented according to:
#   Klaus Witt, CIE Color Difference Metrics, Colorimetry: Understanding the CIE
#   System. 2007
#
# Args:
#   a, b: Any two colors.
#
# Returns:
#   The CIEDE2000 color difference metric evaluated between a and b.
#

pow7(x) = (y = x*x*x; y*y*x)
const twentyfive7 = 25^7

# Delta E 2000
function colordiff(ai::Color, bi::Color, m::DE_2000)
    # Ensure that the input values are in L*a*b* space
    a_Lab = convert(Lab, ai)
    b_Lab = convert(Lab, bi)

    # Calculate some necessary factors from the L*a*b* values
    ac, bc = sqrt(a_Lab.a^2 + a_Lab.b^2), sqrt(b_Lab.a^2 + b_Lab.b^2)
    mc = (ac + bc)/2
    g = (1 - sqrt(pow7(mc) / (pow7(mc) + twentyfive7))) / 2
    a_Lab = Lab(a_Lab.l, a_Lab.a * (1 + g), a_Lab.b)
    b_Lab = Lab(b_Lab.l, b_Lab.a * (1 + g), b_Lab.b)

    # Convert to L*C*h, where the remainder of the calculations are performed
    a = convert(LCHab, a_Lab)
    b = convert(LCHab, b_Lab)

    # Calculate the delta values for each channel
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H*
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Calculate mean L* and C* values
    ml, mc = (a.l + b.l) / 2, (a.c + b.c) / 2

    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # lightness weight
    mls = (ml - 50)^2
    sl = 1.0 + 0.015 * mls / sqrt(20 + mls)

    # chroma weight
    sc = 1 + 0.045mc

    # hue weight
    t = 1 - 0.17 * cosd(mh - 30) +
            0.24 * cosd(2mh) +
            0.32 * cosd(3mh + 6) -
            0.20 * cosd(4mh - 63)
    sh = 1 + 0.015 * mc * t

    # rotation term
    dtheta = 30 * exp(-((mh - 275)/25)^2)
    cr = 2 * sqrt(pow7(mc) / (pow7(mc) + twentyfive7))
    tr = -sind(2*dtheta) * cr

    # Final calculation
    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2 +
         tr * (dc/(m.kc*sc)) * (dh/(m.kh*sh)))
end

# Delta E94
function colordiff(ai::Color, bi::Color, m::DE_94)

    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate the delta values for each channel
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H*
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Calculate geometric mean of chroma
    mc = sqrt(a.c*b.c)

    # Lightness, hue, chroma correction terms
    sl = 1
    sc = 1+0.045*mc
    sh = 1+0.015*mc

    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2)

end

# Metric form of jpc79 color difference equation (mostly obsolete)
function colordiff(ai::Color, bi::Color, m::DE_JPC79)

    # Convert directly into LCh
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate deltas in each direction
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H* from C*'s and h
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    #Calculate mean lightness
    ml = (a.l + b.l)/2
    mc = (a.c + b.c)/2
    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # L* adjustment term
    sl = 0.08195*ml/(1+0.01765*ml)

    # C* adjustment term
    sc = 0.638+(0.0638*mc/(1+0.0131*mc))

    # H* adjustment term
    if (mc < 0.38)
        sh = sc
    elseif (mh >= 164 && mh <= 345)
        sh = sc*(0.56+abs(0.2*cosd(mh+168)))
    else
        sh = sc*(0.38+abs(0.4*cosd(mh+35)))
    end

    # Calculate the final difference
    sqrt((dl/sl)^2 + (dc/sc)^2 + (dh/sh)^2)

end

# Metric form of the cmc color difference
function colordiff(ai::Color, bi::Color, m::DE_CMC)

    # Convert directly into LCh
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate deltas in each direction
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H* from C*'s and h
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Find the mean value of the inputs to use as the "standard"
    ml, mc = (a.l + b.l)/2, (a.c + b.c)/2

    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # L* adjustment term
    if (a.l <= 16)
        sl = 0.511
    else
        sl = 0.040975*ml/(1+0.01765*ml)
    end

    # C* adjustment term
    sc = 0.0638*mc/(1+0.0131*mc)+0.638

    f = sqrt((mc^4)/(mc^4 + 1900))

    if (mh >= 164 && mh < 345)
        t = 0.56 + abs(0.2*cosd(mh+168))
    else
        t = 0.36 + abs(0.4*cosd(mh+35))
    end

    # H* adjustment term
    sh = sc*(t*f+1-f)

    # Calculate the final difference
    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/sh)^2)

end

# The BFD color difference equation
function colordiff(ai::Color, bi::Color, m::DE_BFD)

    # We have to start back in XYZ because BFD uses a different L equation
    a_XYZ = convert(XYZ, ai, m.wp)
    b_XYZ = convert(XYZ, bi, m.wp)

    la = 54.6*log10(a_XYZ.y+1.5)-9.6
    lb = 54.6*log10(b_XYZ.y+1.5)-9.6

    # Convert into LCh with the proper white point
    a_Lab = convert(Lab, a_XYZ, m.wp)
    b_Lab = convert(Lab, b_XYZ, m.wp)
    a1 = convert(LCHab, a_Lab)
    b1 = convert(LCHab, b_Lab)

    # Substitute in the different L values into the L*C*h values
    a = LCHab(la, a1.c, a1.h)
    b = LCHab(lb, b1.c, b1.h)

    # Calculate deltas in each direction
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end

    # Calculate H* from C*'s and h
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Find the mean value of the inputs to use as the "standard"
    ml, mc = (a.l + b.l)/2, (a.c + b.c)/2

    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # Correction terms for a variety of nonlinearities in CIELAB.
    g = sqrt(mc^4/(mc^4 + 14000))

    t = 0.627 + 0.055*cosd(mh - 245) - 0.04*cosd(2*mh - 136) + 0.07*cosd(3*mh - 32) + 0.049*cosd(4*mh + 114) - 0.015*cosd(5*mh + 103)

    rc = sqrt(mc^6/(mc^6 + 7e7))

    rh = -0.26cosd(mh-308) - 0.379cosd(2*mh-160) - 0.636*cosd(3*mh - 254) + 0.226cosd(4*mh + 140) - 0.194*cosd(5*mh + 280)

    dcc = 0.035*mc/(1+0.00365*mc) + 0.521
    dhh = dcc*(g*t+1-g)
    rt = rc*rh

    # Final calculation
    sqrt((dl/m.kl)^2 + (dc/(m.kc*dcc))^2 + (dh/dhh)^2 + rt*((dc*dh)/(dcc*dhh)))
end

# Delta E*ab (the original)
function colordiff(ai::Color, bi::Color, m::DE_AB)

    # Convert directly into L*a*b*
    a = convert(Lab, ai)
    b = convert(Lab, bi)

    dl, da, db = (b.l - a.l), (b.a - a.a), (b.b - a.b)

    sqrt(dl^2 + da^2 + db^2)

end

# Evaluate the DIN99 color difference formula, implemented according to the
# DIN 6176 specification.
#
# Args:
#   a, b: Any two colors.
#
# Returns:
#   The DIN99 color difference metric evaluated between a and b.
function colordiff(ai::Color, bi::Color, m::DE_DIN99)

    a = convert(DIN99, ai)
    b = convert(DIN99, bi)

    sqrt((a.l - b.l)^2 + (a.a - b.a)^2 + (a.b - b.b)^2)

end

# A color difference formula for the DIN99d uniform color space
function colordiff(ai::Color, bi::Color, m::DE_DIN99d)

    a = convert(DIN99d, ai)
    b = convert(DIN99d, bi)

    sqrt((a.l - b.l)^2 + (a.a - b.a)^2 + (a.b - b.b)^2)

end

# The DIN99o color difference metric evaluated between colors a and b.
function colordiff(ai::Color, bi::Color, m::DE_DIN99o)

    a = convert(DIN99o, ai)
    b = convert(DIN99o, bi)

    sqrt((a.l - b.l)^2 + (a.a - b.a)^2 + (a.b - b.b)^2)

end

# Default to Delta E 2000
"""
    colordiff(a, b)
    colordiff(a, b, metric)

Compute an approximate measure of the perceptual difference between
colors `a` and `b`.  Optionally, a `metric` may be supplied, chosen
among `DE_2000` (the default), `DE_94`, `DE_JPC79`, `DE_CMC`,
`DE_BFD`, `DE_AB`, `DE_DIN99`, `DE_DIN99d`, `DE_DIN99o`.
"""
colordiff(ai::Color, bi::Color) = colordiff(ai::Color, bi::Color, DE_2000())


text:
include("actorcritic.jl")
include("reinforce.jl")
include("sarsa.jl")
include("q.jl")
include("fdpolicysearch.jl")
include("nactd.jl")
include("ppo.jl")

export ActorCritic, REINFORCE, NACTD
export Sarsa, ExpectedSarsa, QLearning
export PPO
export AbstractFiniteDifferencePolicySearch, BasicRandomSearch, AugmentedRandomSearch

function act!(agent::TA, env::TE, rng::AbstractRNG, learn::Bool) where {TA<:AbstractPolicyAgent,TE<:AbstractEnvironment} # do not learn if learn flag is set
    logp = get_action!(agent.π, env.state, rng)
    r = step!(env, agent.π.action, rng)
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# dependencies:
# - none

module datastructs

include("misc.jl")
include("helpers.jl")
include("error_handling.jl")
include("stack.jl")
include("queue.jl")
include("linked_list.jl")

include("tree_search_helper.jl")
include("binary_tree.jl")
include("red_black_tree.jl")
include("arbitrary_tree.jl")
include("graph.jl")

export Stack, MatrixStack
export Zero, One

end # datastructs

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module MosekMathProgTest

import Mosek
using MathProgBase
using MathProgBase.SolverInterface

include(joinpath(Pkg.dir("MathProgBase"),"test","linprog.jl"))
include(joinpath(Pkg.dir("MathProgBase"),"test","quadprog.jl"))
include(joinpath(Pkg.dir("MathProgBase"),"test","mixintprog.jl"))
include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
include(joinpath(Pkg.dir("MathProgBase"),"test","linproginterface.jl"))
include(joinpath(Pkg.dir("MathProgBase"),"test","nlp.jl"))

solver = Mosek.MosekSolver(QUIET=true)

linprogtest(solver)
mixintprogtest(solver)
quadprogtest(solver)
convexnlptest(solver)

#socptest(MosekSolver())
coniclineartest(solver,duals=true)
conicSOCtest(solver,duals=true)
conicSOCRotatedtest(solver,duals=true)
conicSDPtest(solver,duals=true)

linprogsolvertest(solver)

# linprogsolvertestextra(solver)
end


text:
module RenderManager

export Renderer

using CoreExtended

using ..GraphicsManager
using ..StorageManager
using ..ModelManager
using ..CameraManager
using ..TransformManager
using ..ShaderManager
using ..TextureManager
using ..GameObjectManager

const Transform=TransformManager.Typ
const Model=ModelManager.Typ
const Camera=CameraManager.Typ
const Texture=TextureManager.Typ
const ShaderProgram=ShaderManager.Typ

const DEFAULT_FUNCTION = ()->nothing

""" TODO """
OnPreDraw = DEFAULT_FUNCTION

""" TODO """
OnPostDraw = DEFAULT_FUNCTION

""" TODO """
OnPreRender = DEFAULT_FUNCTION

""" TODO """
OnPostRender = DEFAULT_FUNCTION

""" TODO """
OnRender = DEFAULT_FUNCTION

""" TODO """
function setOnDraw(f1::Union{Void,Function}, f2::Union{Void,Function})
	global OnPreDraw = f1 != nothing ? f1 : DEFAULT_FUNCTION
	global OnPostDraw = f2 != nothing ? f2 : DEFAULT_FUNCTION
end

""" TODO """
function setOnRender(f1::Union{Void,Function}, f2::Union{Void,Function}, f3::Union{Void,Function})
	global OnPreRender = f1 != nothing ? f1 : DEFAULT_FUNCTION
	global OnPostRender = f2 != nothing ? f2 : DEFAULT_FUNCTION
	global OnRender = f3 != nothing ? f3 : DEFAULT_FUNCTION
end

""" TODO """
type RenderGroup
	enabled::Bool
	transform::Transform
	model::Model
	shaderProgram::ShaderProgram
	textures::Dict{Symbol,Texture}
	modes::Dict{Symbol,Any}
end

""" TODO """
type RenderGroups
	enabled::Bool
	camera::Camera
	groups::SortedDict{Symbol,RenderGroup}
end

""" TODO """
type Renderer
	id::Symbol
	#gameObject::GameObject

	order::Symbol
	enabled::Bool
	transform::Transform
	model::Model
	camera::Camera
	shaderProgram::ShaderProgram
	textures::Dict{Symbol,Texture}
	modes::Dict{Symbol,Any}
	setModes::Function

	function Renderer(id=:_)
		#obj=GameObjectManager.create(id)
		this=new(id,:_,true,nothing,nothing,nothing,nothing,Dict(),Dict(),()->nothing)
		#GameObjectManager.setComponent(obj,this)
		this
	end
end

""" TODO """
function link(this::Renderer)
	ModelManager.linkTo(this.model)
	ShaderManager.linkTo(this.shaderProgram)
	if this.camera != nothing CameraManager.update(this.camera, this.transform) end
end

#createTransform(this::Renderer) = this.transform = TransformManager.create(this.id)
#createModel(this::Renderer) = this.model = ModelManager.create(this.id)
#createCamera(this::Renderer) = this.camera = CameraManager.create(this.id)

""" TODO """
function register(this::Renderer)
	this.model = ModelManager.getSelected()
	this.camera = CameraManager.getSelected()
	this.shaderProgram = ShaderManager.getSelected()
	#RenderManager.createShaderProgram(:FG, "shaders/default.glsl")
end

""" TODO """
function addTexture(this::Renderer, path::String)
	id=Symbol(path)
	texture = TextureManager.create(id)
	this.textures[id] = texture
	TextureManager.load(texture, path)
	texture
end

""" TODO """
function update()
	ShaderManager.watchShaderReload(RenderManager.reload)
end

""" TODO """
function init(this::Renderer)
	setMode(this, :POLYGONMODE,[:FRONT_AND_BACK,:FILL])
	this.transform=TransformManager.create(this.id)
end

""" TODO """
function render()
	OnPreRender()
	OnRender()
	OnPostRender()
end

""" TODO """
function render(this::Renderer)
	if !this.enabled return end
	setSelected(this)
	link() # shaderprogram, model, texture, ...
	getShaderProperties(this) 	# set propberties
	#this.setModes()
	OnPreDraw()
	draw(this)
	OnPostDraw()
end

""" TODO """
function resets()
	ModelManager.reset()
	ShaderManager.reset()
end

""" TODO """
function unlinks()
	ModelManager.unlink()
	ShaderManager.unlink()
end

# preset for manager: includes code here
presetManager(Renderer)

useSortedList = false
renderSortedList(use::Bool) = (global useSortedList = use)

""" TODO """
renderDefault() = for (_,this) in list; render(this); end

""" TODO """
renderAll = DEFAULT_FUNCTION

""" TODO """
function start()
	reset()
	RenderManager.reload()

	# shader reload
	for (k,this) in list
		if this.shaderProgram != nothing ShaderManager.reload(this.shaderProgram) end
	end
end

""" TODO """
function reload()
	if !useSortedList
		renderAll = stabilize(renderDefault)
	else
		#all=Expr(:block,:())

		sortedList=SortedDict{Symbol, Array{Renderer,1}}(Forward)

		order=0
		for (k,this) in list
			order+=1
			ex=:()
			for (k,v) in this.modes
				ex=:($(ex.args...); $(GRAPHICSDRIVER().getMode(k,v).args...);)
			end

			this.setModes=()->nothing #eval(:(function(); $(ex.args...); end;))

			if this.order == :_ this.order=Symbol(order) end
			if !haskey(sortedList,this.order) sortedList[this.order]=Array{Renderer,1}(); end
			push!(sortedList[this.order], this)
		end

		global renderAll = stabilize(eval(:(function();
			for (_,a) in $sortedList; for this in a; render(this); end; end;
		end;)))
		println("sortedList generated.")
	end

	global OnRender = OnRender != DEFAULT_FUNCTION ? OnRender : renderAll
end

""" TODO """
function copy(l::Renderer,r::Renderer)
	l.model = r.model
	l.camera = r.camera
	l.transform = r.transform
	l.shaderProgram = r.shaderProgram
end

""" TODO """
function unlink(typ::DataType)
	#obj = typ == Camera ? EMPTY_CAMERA : (typ == Model ? EMPTY_MODEL : (typ == Model ? nothing : nothing))
	if obj != nothing link(obj) end
end

""" TODO """
function link(ent::Transform)
	this = getSelected()
	this.transform = ent
end

""" TODO """
function link(camera::Camera)
	this = getSelected()
	this.camera = camera
end

""" TODO """
function link(model::Model)
	this = getSelected()
	this.model = model
end

""" TODO """
function link(program::ShaderProgram)
	this = getSelected()
	this.shaderProgram = program
end

""" TODO """
is(id::Symbol) = RenderManager.getSelected().id == id

""" TODO """
isTransform(id::Symbol) = RenderManager.getSelected().transform.id == id

""" TODO """
isModel(id::Symbol) = RenderManager.getSelected().model.id == id

""" TODO """
isCamera(id::Symbol) = RenderManager.getSelected().camera.id == id

""" TODO """
getMVP() = RenderManager.getSelected().camera.mvpMat

""" TODO """
getData(k::Symbol) = ModelManager.getData(getSelected().model,k)

""" TODO """
function init()
	ShaderManager.setListenerOnShaderPropertyUpdate(:OnShaderPropertyUpdateAttributes, (
		"iVertex" => () ->  getData(:VERTEX)
		,"iUV" => () -> (d=getData(:UV); d != nothing ? d : getData(:VERTEX))
		,"iNormal" => () -> (d=getData(:NORMAL); d != nothing ? d : getData(:VERTEX))
		)
	)
end

""" TODO """
function getShaderProperties(this::Renderer)
	ShaderManager.getShaderProperties(this.shaderProgram, ModelManager.getMainBuffer(this.model).data)
end

""" TODO """
setRenderMode(this::Renderer, mode::Symbol) = GRAPHICSDRIVER().update(ModelManager.getGroup(this.model,:VERTEX).data, mode)

""" TODO """
setMode(this::Renderer, mode::Symbol, value::Any) = this.modes[mode]=value

""" TODO """
setRenderMode(mode::Symbol) = setRenderMode(getSelected(), mode)

""" TODO """
setMode(mode::Symbol, value::Any) = GRAPHICSDRIVER().setMode(mode,value)

""" TODO """
function draw(this::Renderer)
	model=this.model
	buffer=ModelManager.getMainBuffer(model)
	if model.enabled && ModelManager.isLinked(model) && buffer != nothing && buffer.data != nothing && (data=StorageManager.getData(buffer, 1)) != nothing
		ShaderManager.render(this.shaderProgram, buffer.data, data)
	end
end

end #RenderManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using Tokens
using Test

function testtoken()

    t = BToken(T_TEXT,"äöü")
    s = convert(SubString{String},t)
    @test s == t
end


text:
facts("Dipoles") do

    dips = Dipole[]

    context("Create") do

        dip1 = Dipole("Talairach", 1, 2, 1, 0, 0, 0, 1, 1, 1)
        dip2 = Dipole("Talairach", 1, 2, 3, 0, 0, 0, 2, 2, 2)

        dips = push!(dips, dip1)
        dips = push!(dips, dip2)

        @fact length(dips) --> 2
    end

    context("Show") do

        show(dips[1])
        show(dips)
    end

    context("Mean") do

        b = mean(dips)
        @fact b.x --> 1.0 * SIUnits.ShortUnits.m
        @fact b.y --> 2.0 * SIUnits.ShortUnits.m
        @fact b.z --> 2.0 * SIUnits.ShortUnits.m
    end

    context("Std") do

        b = std(dips)
        @fact b.x --> 0.0 * SIUnits.ShortUnits.m
        @fact b.y --> 0.0 * SIUnits.ShortUnits.m
        @fact b.z --> std([1, 3]) * SIUnits.ShortUnits.m
    end
end


text:
# Matern 3/2 ARD covariance function

@doc """
# Description
Constructor for the ARD Matern 3/2 kernel (covariance)

k(x,x') = σ²(1+√3*d/L^2)exp(-√3*d/L^2), where d = |x-x'| and L = diag(ℓ₁,ℓ₂,...)
# Arguments:
* `ll::Vector{Float64}`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type Mat32Ard <: StationaryARD
    ℓ2::Vector{Float64}     # Log of Length scale
    σ2::Float64            # Log of Signal std
    Mat32Ard(ll::Vector{Float64}, lσ::Float64) = new(exp(2.0*ll), exp(2.0*lσ))
end

function set_params!(mat::Mat32Ard, hyp::Vector{Float64})
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat32 kernel only has $(num_params(mat)) parameters"))
    d=length(mat.ℓ2)
    mat.ℓ2 = exp(2.0*hyp[1:d])
    mat.σ2 = exp(2.0*hyp[d+1])
end

get_params(mat::Mat32Ard) = [log(mat.ℓ2)/2.0; log(mat.σ2)/2.0]
get_param_names(mat::Mat32Ard) = [get_param_names(mat.ℓ2, :ll); :lσ]
num_params(mat::Mat32Ard) = length(mat.ℓ2) + 1

metric(mat::Mat32Ard) = WeightedEuclidean(1.0./(mat.ℓ2))
cov(mat::Mat32Ard, r::Float64) = mat.σ2*(1+sqrt(3)*r)*exp(-sqrt(3)*r)

function grad_kern(mat::Mat32Ard, x::Vector{Float64}, y::Vector{Float64})
    r = distance(mat, x, y)
    exp_r = exp(-sqrt(3)*r)
    wdiff = ((x-y).^2)./mat.ℓ2

    g1 = 3*mat.σ2*wdiff*exp_r
    g2 = 2.0*mat.σ2*(1+sqrt(3)*r)*exp_r

    return [g1; g2]
end

function grad_stack!(stack::AbstractArray, mat::Mat32Ard, X::Matrix{Float64}, data::StationaryARDData)
    d = size(X,1)
    R = distance(mat, X, data)
    exp_R = exp(-sqrt(3)*R)
    broadcast!(*, view(stack, :, :, 1:d), mat.σ2*3, data.dist_stack, reshape(1.0./mat.ℓ2, (1,1,d)), exp_R)
    stack[:,:, d+1] = 2.0 * mat.σ2 * (1 + sqrt(3) * R) .* exp_R
    return stack
end


text:
@testset "Overrides" begin
    g=MetaDiGraph(star_digraph(4), 8)

    set_indexing_prop!(g,:nodelabel)
    foreach(v->set_indexing_prop!(g,v,:nodelabel,"Vertex $v"),vertices(g))
    foreach(e->set_prop!(g,e,:edgelabel,"Edge ($(e.src),$(e.dst))") ,edges(g))

    h=reverse(g)

    for v in 1:3
        @test get_prop(g, v, :nodelabel) == get_prop(h, v, :nodelabel)
    end

    for sd in 2:4
        @test get_prop(g, 1, sd, :edgelabel) == get_prop(h, sd, 1, :edgelabel)
    end

    @test weightfield(h) == weightfield(g) == :weight
    @test defaultweight(h) == defaultweight(g) == 8

end


text:
function derivative!(df::AbstractArray{<:Number}, f, x::Union{Number,AbstractArray{<:Number}}, fx::AbstractArray{<:Number}, integrator, grad_config)
    alg = unwrap_alg(integrator, true)
    tmp = length(x) # We calculate derivtive for all elements in gradient
    if alg_autodiff(alg)
        ForwardDiff.derivative!(df, f, fx, x, grad_config)
        integrator.destats.nf += 1
    else
        DiffEqDiffTools.finite_difference_gradient!(df, f, x, grad_config, dir = diffdir(integrator))
        fdtype = alg.diff_type
        if fdtype == Val{:forward} || fdtype == Val{:central}
            tmp *= 2
            if eltype(df)<:Complex
              tmp *= 2
            end
        end
        integrator.destats.nf += tmp
    end
    nothing
end

function derivative(f, x::Union{Number,AbstractArray{<:Number}},
                    integrator)
    local d
    tmp = length(x) # We calculate derivtive for all elements in gradient
    alg = unwrap_alg(integrator, true)
    if alg_autodiff(alg)
      integrator.destats.nf += 1
      d = ForwardDiff.derivative(f, x)
    else
      d = DiffEqDiffTools.finite_difference_derivative(f, x, alg.diff_type, dir = diffdir(integrator))
      if alg.diff_type == Val{:central} || alg.diff_type == Val{:forward}
          tmp *= 2
      end
      integrator.destats.nf += tmp
      d
    end
end

function jacobian(f, x, integrator)
    alg = unwrap_alg(integrator, true)
    local tmp
    if alg_autodiff(alg)
      if DiffEqBase.has_colorvec(integrator.f)
        J,tmp = jacobian_autodiff(f, x, integrator)
      else
        J,tmp = jacobian_autodiff(f, x)
      end
    else
      if DiffEqBase.has_colorvec(integrator.f)
        J,tmp = jacobian_finitediff(f, x, alg.diff_type, integrator, integrator.f.colorvec)
      else
        J,tmp = jacobian_finitediff(f, x, alg.diff_type, integrator)
      end
    end
    integrator.destats.nf += tmp
    J
end

jacobian_autodiff(f, x) = (ForwardDiff.derivative(f,x),1)
jacobian_autodiff(f, x::AbstractArray) = (ForwardDiff.jacobian(f, x),1)
function jacobian_autodiff(f, x::AbstractArray, integrator)
  colorvec=integrator.f.colorvec
  jac=integrator.f.jac_prototype
  J=jac isa SparseMatrixCSC ? similar(jac) : fill(0.,size(jac))
  (forwarddiff_color_jacobian!(J,f,x,colorvec=colorvec,sparsity=jac),1)
end
#jacobian_autodiff(f, x::AbstractArray, colorvec) = (ForwardDiff.jacobian(f, x, colorvec = colorvec),1)

function _nfcount(N,diff_type)
  if diff_type==Val{:complex}
    tmp = N
  elseif diff_type==Val{:forward}
    tmp = N + 1
  else
    tmp = 2N
  end
  tmp
end

jacobian_finitediff(f, x, diff_type, integrator) =
    (DiffEqDiffTools.finite_difference_derivative(f, x, diff_type, eltype(x), dir = diffdir(integrator)),2)
jacobian_finitediff(f, x::AbstractArray, diff_type, integrator) =
    (DiffEqDiffTools.finite_difference_jacobian(f, x, diff_type, eltype(x), Val{false},
      dir = diffdir(integrator)),_nfcount(length(x),diff_type))
jacobian_finitediff(f, x::AbstractArray, diff_type, integrator, colorvec) =
    (DiffEqDiffTools.finite_difference_jacobian(f, x, diff_type, eltype(x), Val{false},
      dir = diffdir(integrator), colorvec = colorvec, sparsity = integrator.f.jac_prototype),_nfcount(maximum(colorvec),diff_type))

jacobian_finitediff_forward!(J,f,x,jac_config,forwardcache,integrator)=(DiffEqDiffTools.finite_difference_jacobian!(J,f,x,jac_config,forwardcache,dir=diffdir(integrator));length(x))
jacobian_finitediff_forward!(J,f,x,jac_config,forwardcache,integrator,colorvec)=
  (DiffEqDiffTools.finite_difference_jacobian!(J,f,x,jac_config,forwardcache,
    dir=diffdir(integrator),colorvec=colorvec,sparsity=integrator.f.jac_prototype);maximum(colorvec))
jacobian_finitediff!(J,f,x,jac_config,integrator)=(DiffEqDiffTools.finite_difference_jacobian!(J,f,x,jac_config,dir=diffdir(integrator));2*length(x))
jacobian_finitediff!(J,f,x,jac_config,integrator,colorvec)=
  (DiffEqDiffTools.finite_difference_jacobian!(J,f,x,jac_config,
    dir=diffdir(integrator),colorvec=colorvec,sparsity=integrator.f.jac_prototype);2*maximum(colorvec))
jacobian_autodiff!(J,f,x,jac_config)=forwarddiff_color_jacobian!(J,f,x,jac_config)#J::SparseMatrixCSC
jacobian_autodiff!(J,f,fx,x,jac_config)=ForwardDiff.jacobian!(J,f,fx,x,jac_config)

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, fx::AbstractArray{<:Number}, integrator::DiffEqBase.DEIntegrator, jac_config)
    alg = unwrap_alg(integrator, true)
    if alg_autodiff(alg)
      if DiffEqBase.has_colorvec(integrator.f)
        jacobian_autodiff!(J,f,x,jac_config)
      else
        jacobian_autodiff!(J,f,fx,x,jac_config)
      end
      integrator.destats.nf += 1
    else
      isforward = alg.diff_type === Val{:forward}
      if isforward
        forwardcache = get_tmp_cache(integrator, alg, unwrap_cache(integrator, true))[2]
        f(forwardcache, x)
        integrator.destats.nf += 1
        if DiffEqBase.has_colorvec(integrator.f)
          tmp=jacobian_finitediff_forward!(J, f, x, jac_config, forwardcache, integrator, integrator.f.colorvec)
        else
          tmp=jacobian_finitediff_forward!(J, f, x, jac_config, forwardcache, integrator)
        end
      else # not forward difference
        if DiffEqBase.has_colorvec(integrator.f)
          tmp=jacobian_finitediff!(J, f, x, jac_config, integrator, integrator.f.colorvec)
        else
          tmp=jacobian_finitediff!(J, f, x, jac_config, integrator)
        end
      end
      integrator.destats.nf += tmp
    end
    nothing
end

function DiffEqBase.build_jac_config(alg::Union{OrdinaryDiffEqAlgorithm,DAEAlgorithm},f,uf,du1,uprev,u,tmp,du2)
  if !DiffEqBase.has_jac(f)
    if alg_autodiff(alg)
      if DiffEqBase.has_colorvec(f)
        jac_config = ForwardColorJacCache(uf,uprev,colorvec=f.colorvec,sparsity=f.jac_prototype)
      else
        jac_config = ForwardDiff.JacobianConfig(uf,du1,uprev,ForwardDiff.Chunk{determine_chunksize(u,alg)}())
      end
    else
      if alg.diff_type != Val{:complex}
        jac_config = DiffEqDiffTools.JacobianCache(tmp,du1,du2,alg.diff_type)
      else
        jac_config = DiffEqDiffTools.JacobianCache(Complex{eltype(tmp)}.(tmp),Complex{eltype(du1)}.(du1),nothing,alg.diff_type,eltype(u))
      end
    end
  else
    jac_config = nothing
  end
  jac_config
end

get_chunksize(jac_config::ForwardDiff.JacobianConfig{T,V,N,D}) where {T,V,N,D} = N

function DiffEqBase.resize_jac_config!(jac_config::ForwardDiff.JacobianConfig, i)
  for j in eachindex(jac_config.duals)
    resize!(jac_config.duals[j], i)
  end
  jac_config
end

function DiffEqBase.resize_jac_config!(jac_config::DiffEqDiffTools.JacobianCache, i)
  resize!(jac_config, i)
  jac_config
end

function resize_grad_config!(grad_config::ForwardDiff.DerivativeConfig, i)
  resize!(grad_config.duals, i)
  grad_config
end

function resize_grad_config!(grad_config::DiffEqDiffTools.GradientCache, i)
  @unpack fx, c1, c2 = grad_config
  fx !== nothing && resize!(fx, i)
  c1 !== nothing && resize!(c1, i)
  c2 !== nothing && resize!(c2, i)
  grad_config
end

function build_grad_config(alg,f,tf,du1,t)
  if !DiffEqBase.has_tgrad(f)
    if alg_autodiff(alg)
      grad_config = ForwardDiff.DerivativeConfig(tf,du1,t)
    else
      grad_config = DiffEqDiffTools.GradientCache(du1,t,alg.diff_type)
    end
  else
    grad_config = nothing
  end
  grad_config
end


text:
module Sog

export sog
"""
    sog(iVector=[], orders = eachindex(iVector))

The function `sog()` takes a `Vector` of any type `iVector` as its input with any possible length, and will return a `Bool Vector` with the same length as input. Values
 `true` or `false` in the output depends on whether one element in the `Vector` is the start of a group or not. The sog() will return
 'nothing' if the length of input is 0. As any type can be applied in the input, be careful! To check whether two elements `a` and `b` are treated
as the same, please call `isequal(a, b)` to see the result.

When taking `Vector` of `Vectors`, the length of output should the same as the length of every `Vector`, and what `sog()` do is to compare
one element with the previous one on the same position for each `Vector`.

The additional argument `orders` can work for a `Vector` of `Vectors` if you are sure that the change of order can get a faster speed.
This is not perfect for users who are not familar with types of data so just leave it there is fine. This part might be removed to let
the program does the thing in the future.

Parallel computation is used so, by setting your environment variable "JULIA_NUM_THREADS" to more than one, you can get a faster speed with this package.

Examples:

    sog([1,2,2,2,1,3,3,1,1]) returns: [true, true, false, false, true, true, false, true, false]
    sog() returns: nothing
    sog(["a", "a", "a", 'a']) returns: [true, false, false, true]
    sog([[1, 1, 1, 2, 2], [1.5, 1.5, 20.0, 3.0, 3.0]]) returns: [true, false, true, true, false]
    sog([[1, 1, 1, 2, 2], [1.5, 1.5, 20.0, 3.0, 3.0]], [2, 1]) returns: [true, false, true, true, false]

"""
function sog(iVector=[], orders = eachindex(iVector))
    if iVector == []
        return nothing
    end

    len = length(iVector)

    if typeof(iVector[1]) <: Vector{}           #For Vector of Vectors
        lenv = length(iVector[1])
        if lenv == 0
            return nothing
        end

        if lenv == 1
            return [true]
        end

        if lenv >= 2
            if length(eachindex(iVector)) != length(orders)
                println("Please specify the order of columns correctly! Number of columns is the length of order.")
                return nothing
            end
            oVector = ini0(iVector[orders[1]])
            for i in orders[2:end]
                oVector = zore_or_one(oVector, iVector[i])
            end
            return oVector
        end

    elseif len == 1
        return [true]

    else
        oVector = ini0(iVector)
    end

    return oVector
end

function ini0(x1)
    len = length(x1)
    re = zeros(Bool, len)
    re[1] = 1
    Threads.@threads for j in 2:len
        @inbounds re[j] = !(x1[j]===x1[j-1])
    end
    return re
end

function zore_or_one(re, xi)
    Threads.@threads for j in 2:length(re)
        @inbounds re[j] = re[j]==1 ? 1 : !(xi[j]===xi[j-1])
    end
    return re
end

end