text:
# MultipleBranches.jl
#  Extending the Feynman theory to multiple phonon branches

"""
frohlichPartial((f, ϵ_mode); ϵ_o,ϵ_s,meff)

Calculate a (partial) dielectric electron-phonon coupling element.
ϵ_o - optical dielectric
ϵ_s - total static dielectric contribution
ϵ_mode - this mode's contribution to dielectric
"""
function frohlichPartial((f, ϵ_mode); ϵ_o,ϵ_s,meff)
    ω=f*1E12*2π
    α=1/(4*π*ϵ_0) * ϵ_mode/(ϵ_o*ϵ_s) * (q^2/ħ) * sqrt(meff*me/(2*ω*ħ))
    return α
end

# deprecated signature wrapped via multiple dispatch
frohlichPartial(ϵ_o,ϵ_s,ϵ_mode,f,meff) = frohlichPartial((f, ϵ_mode), ϵ_o=ϵ_o,ϵ_s=ϵ_s,meff=meff)

rows(M::Matrix) = map(x->reshape(getindex(M, x, :), :, size(M)[2]), 1:size(M)[1])
"""
IRtoDielectric(IR,V)

From absolute value of IR activities of phonon modes, generate a per-mode
contribution to the low-frequencydielectric constant.
"""
function IRtoDielectric(IR,volume)
    ϵ=0.0 #* q^2 * THz^-2 * amu^-1 * m^-3
    for r in rows(IR)
        f,S=r # frequency in THz; activity in e^2 amu^-1
        f=f * 1E12 #* THz
        ω=2π * f
        S=S * q^2 / amu
        ϵ_mode = S / ω^2 / volume
        ϵ_mode /= 3 # take isotropic average = divide by 3
        ϵ+=ϵ_mode
        println("Mode f= $f S= $S ϵ_mode = $(ϵ_mode/ϵ_0)")
    end
    println("Raw ionic dielectric contribution: $ϵ absolute $(ϵ/ϵ_0) relative")
    return ϵ/ϵ_0
end

"""
IRtoalpha(IR,volume)

Calculates contribution to dielectric constant for each polar phonon mode, and
thereby the Frohlich alpha contribution for this mode.
"""
function IRtoalpha(IR; volume, ϵ_o,ϵ_s,meff)
    ϵ=0.0 #* q^2 * THz^-2 * amu^-1 * m^-3
    α_sum=0.0
    for r in rows(IR)
        f,S=r # frequency in THz; activity in e^2 amu^-1
        f=f * 1E12 #* THz
        ω=2π * f
        S=S * q^2 / amu
        ϵ_mode = S / ω^2 / volume
        ϵ_mode /= 3 # take isotropic average = divide by 3
        ϵ_mode /= ϵ_0 # reduced dielectric constant
        ϵ+=ϵ_mode
        #println("Mode f= $f S= $S ϵ_mode = $(upreferred(ϵ_mode/u"ϵ0"))")
        α=frohlichPartial(ϵ_o,ϵ_s,ϵ_mode,f/1E12,meff)
        α_sum+=α
        if (α>0.1)
            println("Notable Mode f= $f    α_partial=$α")
        end
    end
    println("Sum alpha: $(α_sum)")
    return f, α_sum
end

"""
DieletricFromIRmode(IRmode)

Calculate dielectric from an individual mode.
"""
function DielectricFromIRmode(IRmode; volume)
    f,S=IRmode # Assumes Frequency in THz; activity in e^2 amu^-1
    f*=1E12 # convert to Hz from THz
    ω=2π * f
    S=S * q^2 / amu
    ϵ_mode = S / ω^2 / volume
    ϵ_mode /= 3 # take isotropic average = divide by 3
    ϵ_mode /= ϵ_0 # reduced dielectric
    return ϵ_mode
end

function Hellwarth1999mobilityRHS( (α, (v,w) ,f) ,effectivemass, T)
    mb=effectivemass*MassElectron
    ω=f*1E12*2π
    βred=ħ*ω/(kB*T)

    R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)
    b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)
    a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
    k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
    K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

    #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
    RHS=α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
    μ=RHS^-1 * (q)/(ω*mb)

    return 1/μ
end


text:
using Test

@test 2 + 2 == 4

ext: jl


text:
"""
# References
[1] T. Bian and Z.-P. Jiang, “Value Iteration, Adaptive Dynamic Programming, and Optimal Control of Nonlinear Systems,” in 2016 IEEE 55th Conference on Decision and Control (CDC), Las Vegas, NV, USA, Dec. 2016, pp. 3375–3380. doi: 10.1109/CDC.2016.7798777.
"""
struct TwoDimensionalNonlinearPolynomialSystem <: AbstractEnv
end

function State(env::TwoDimensionalNonlinearPolynomialSystem)
    return function (x1::Real=-2.9, x2::Real=-2.9)
        ComponentArray(x1=x1, x2=x2)
    end
end

function Dynamics!(env::TwoDimensionalNonlinearPolynomialSystem)
    @Loggable function dynamics!(dx, x, p, t; u)
        @assert length(u) == 1
        @log state = x
        @log input = u
        _u = u[1]  # Real or Array
        @unpack x1, x2 = x
        dx.x1 = -(1/2)*x1^3 - x1 - 2*x2
        dx.x2 = (1/8)*x2^3 - x2 + (1/2)*_u^3
        nothing
    end
end

function OptimalControl(env::TwoDimensionalNonlinearPolynomialSystem)
    return function (x::ComponentArray, p, t)
        @unpack x1, x2 = x
        -x2
    end
end

function OptimalValue(env::TwoDimensionalNonlinearPolynomialSystem)
    return function (x::ComponentArray)
        @unpack x1, x2 = x
        x1^2 + x2^2
    end
end

function RunningCost(env::TwoDimensionalNonlinearPolynomialSystem)
    return function (x::ComponentArray, _u)
        @assert length(_u) == 1
        u = _u[1]  # Real or Array
        @unpack x1, x2 = x
        x1^4 + 2*(x1+x2)^2 + (3/4)*u^4
    end
end


text:
export localSympFromQuad, getDecoupleSequence, getInterferenceBasedSequence, getSequence

# From (source quadrautre, target quadratre) to Local Symplectic Matrix transforming source to Ω target
↑(m, i) = 2*m - (i % 2)
↓(m, i) = 2*m -1 + (i  % 2)
getPattern(v::Vector)=(n->n>tolerance).( (diag∘sqrt∘Diagonal∘(x->x'*x)∘(v->reshape(v, 2, :)))(v) )
function Lkmij(u, v, m, i, j)
    if norm(v[2*m-1:2*m])<tolerance
        return i==j ? 1 : 0
    else
        return (-1)^(j+1)*v[↑(m, i)]*u[↓(m, j)]/( v[2*m]^2 + v[2*m-1]^2 ) + (-1)^i*v[↓(m, i)]*u[↑(m, j)]/( u[2*m]^2 + u[2*m-1]^2 )
    end
end
function localSympFromQuad(u::Vector, v::Vector; lambda::Real = 1)::Symp
    if length(u) % 2 !=0 || length(u) != length(v) || getPattern(u) != getPattern(v)
        return error("input not valid")
    end
    n = length(u) ÷ 2
    if n == 1
        return Symp([ Lkmij(lambda * u, v, 1, i, j) for i in 1:2, j in 1:2 ])
    else
        return ⊕([ Symp([ Lkmij(lambda * u, v, m, i, j) for i in 1:2, j in 1:2 ]) for m in 1:n ]...)
    end
end
# Get Local Operations L3 L2 L1 for the Decoupling Sequence S4 L3 S3 L2 S2 L1 S1
function getDecoupleSequence(S4::Symp, S3::Symp, S2::Symp, S1::Symp, m::Int=1)::Vector{Symp}
    # if !(all(x->isGeneric(x), [S1, S2, S3, S4]) ) return error("not generic") end
    L1 = localSympFromQuad(S1.S[:,2*m-1], S2.S[2*m-1,:])
    L3 = localSympFromQuad(S3.S[:,2*m-1], S4.S[2*m-1,:])
    T1 = S2 * L1 * S1
    T2 = S4 * L3 * S3
    u = T1.S[:,2*m]
    v = T2.S[2*m,:] + sum( (k->T1.S[k,2*m-1]*T1.S[k,2*m]-T2.S[2*m-1,k]*T2.S[2*m,k]).(1:length(u))  ) * T2.S[2*m-1, :]
    L2 = localSympFromQuad(u, v; lambda = lambda)
    return [S4, L3, S3, L2, S2, L1, S1]
end

# Get local Operations L16, L15, ..., L1 for the interfernce-basded sequence LR S L15 S ... S L1 S, with a given 4×4 target Symplectic Matrix S⊙
function getSequence(Ss::Vector{Symp{T}}, modeToDecouple::Int)::Vector{Symp} where T
    if length(Ss) == 4
        return getDecoupleSequence(Ss..., modeToDecouple)
    else
        step = length(Ss) ÷ 4
        Rs = [ getSequence(Ss[i:i+step-1], modeToDecouple-1 ) for i in 1:step:length(Ss) ]

        Ls = getDecoupleSequence([ *(RLst...) for RLst in Rs]..., modeToDecouple)
        return vcat([ (i%2==0) ? [ Ls[i] ] : Rs[(i+1) ÷ 2] for i in 1:7 ]...)
    end
end

function getInterferenceBasedSequence(S::Symp, ST::Symp)::Vector{Symp}
    if !(isGeneric(S))  return error("not generic") end
    Ss = [S for i in 1:4^nModes(ST)]
    Ss[end] =  Ss[end] * (ST ⊕ Id(nModes(S) - nModes(ST)))^-1
    Ss = getSequence(Ss, nModes(ST))
    LR = Symp( ( *(Ss...).S[1:2*nModes(ST), 1:2*nModes(ST)]) )^-1  ⊕ Id(nModes(S) - nModes(ST))
    Ss[end] =  Ss[end] * (ST ⊕ Id(nModes(S) - nModes(ST)))
    return vcat([ LR ] , Ss)
end

# Interference-Based Sequence with Automatic Local Randomization
function getSequenceWithRandomization(Ss::Vector{Symp{T}}, modeToDecouple::Int)::Vector{Symp} where T
    if length(Ss) == 4
        return getDecoupleSequence(Ss..., modeToDecouple)
    else
        step = length(Ss) ÷ 4
        Rs = [ getSequence(Ss[i:i+step-1], modeToDecouple-1 ) for i in 1:step:length(Ss) ]

        Ls = getDecoupleSequence([ *(RLst...) for RLst in Rs]..., modeToDecouple)
        return vcat([ (i%2==0) ? [ Ls[i] ] : Rs[(i+1) ÷ 2] for i in 1:7 ]...)
    end
end
function getIQSWithRandomization(S::Symp, ST::Symp)::Vector{Symp}
    if !(isGeneric(S))  return error("not generic") end
    Ss = [S for i in 1:4^nModes(ST)]
    Ss[end] =  Ss[end] * (ST ⊕ Id(nModes(S) - nModes(ST)))^-1
    Ss = getSequenceWithRandomization(Ss, nModes(ST))
    LR = Symp( ( *(Ss...).S[1:2*nModes(ST), 1:2*nModes(ST)]) )^-1  ⊕ Id(nModes(S) - nModes(ST))
    Ss[end] =  Ss[end] * (ST ⊕ Id(nModes(S) - nModes(ST)))
    return vcat([ LR ] , Ss)
end


text:
## School groups

# remove ++++++

# n_school(s :: SchoolGroups3) = UInt8(3);
# s_groups(s :: SchoolGroups3) = [:HSG, :SC, :CG];

# # Index from Symbol. E.g., :HSG => 1
# s_idx(s :: SchoolGroups3, edLevel :: Symbol) =
#     findfirst(s_groups(s) .== edLevel);

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


text:
# indicator of the (primal) exponential cone
# the dual exponential cone is obtained through calculus rules

"""
  IndExpPrimal()

Returns the indicator function of the primal exponential cone, that is
`cl{(r,s,t) : s > 0, s⋅exp(r/s) ⩽ t}`.
"""

immutable IndExpPrimal <: IndicatorConvexCone end

"""
  IndExpDual()

Returns the indicator function of the dual exponential cone, that is
`cl{(u,v,w) : u < 0, -u⋅exp(v/u) ⩽ w⋅exp(1)}`.
"""

IndExpDual() = Precompose(Conjugate(IndExpPrimal()), -1.0)

EXP_PRIMAL_CALL_TOL = 1e-6
EXP_POLAR_CALL_TOL = 1e-3
EXP_PROJ_TOL = 1e-15
EXP_PROJ_MAXIT = 100

function (f::IndExpPrimal){R <: Real}(x::AbstractArray{R,1})
  if (x[2] > 0.0 && x[2]*exp(x[1]/x[2]) <= x[3]+EXP_PRIMAL_CALL_TOL) ||
     (x[1] <= EXP_PRIMAL_CALL_TOL && abs(x[2]) <= EXP_PRIMAL_CALL_TOL && x[3] >= -EXP_PRIMAL_CALL_TOL)
    return 0.0
  end
  return +Inf
end

function (f::Conjugate{IndExpPrimal}){R <: Real}(x::AbstractArray{R,1})
  if (x[1] > 0.0 && x[1]*exp(x[2]/x[1]) <= -exp(1)*x[3]+EXP_POLAR_CALL_TOL) ||
     (abs(x[1]) <= EXP_POLAR_CALL_TOL && x[2] <= EXP_POLAR_CALL_TOL && x[3] <= EXP_POLAR_CALL_TOL)
    return 0.0
  end
  return +Inf
end

# Projection onto the cone is performed as in SCS (https://github.com/cvxgrp/scs).
# See the following copyright and permission notices.

# The MIT License (MIT)
#
# Copyright (c) 2012 Brendan O'Donoghue (bodonoghue85@gmail.com)
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function prox!{R <: Real}(f::IndExpPrimal, x::AbstractArray{R,1}, y::AbstractArray{R,1}, gamma::Real=1.0)
  r = x[1]
  s = x[2]
  t = x[3]
  if (s*exp(r/s) <= t && s > 0) || (r <= 0 && s == 0 && t >= 0)
    # x in the cone
    y[:] = x
  elseif (-r < 0 && r*exp(s/r) <= -exp(1)*t) || (-r == 0 && -s >= 0 && -t >= 0)
    # -x in the dual cone (x in the polar cone)
    y[:] = 0.0
  elseif r < 0 && s < 0
    # analytical solution
    y[1] = x[1]
    y[2] = max(x[2], 0.0)
    y[3] = max(x[3], 0.0)
  else
    v = x
    ub, lb = getRhoUb(x)
    for iter = 1:EXP_PROJ_MAXIT
      rho = (ub + lb)/2
      g, v = calcGrad(x,rho)
      if g > 0
        lb = rho
      else
        ub = rho
      end
      if ub - lb <= EXP_PROJ_TOL
        break
      end
    end
    y[:] = v
  end
  return 0.0
end

function getRhoUb(v)
  lb = 0
  rho = 2.0^(-3)
  g, z = calcGrad(v, rho)
  while g > 0
    lb = rho
    rho = rho*2
    g, z = calcGrad(v, rho)
  end
  ub = rho
  return ub, lb
end

function calcGrad(v, rho)
  x = solve_with_rho(v, rho)
  if x[2] == 0.0
    g = x[1]
  else
    g = x[1] + x[2]*log(x[2]/x[3])
  end
  return g, x
end

function solve_with_rho(v, rho)
  x = zeros(3)
  x[3] = newton_exp_onz(rho, v[2], v[3])
  x[2] = (1/rho)*(x[3] - v[3])*x[3]
  x[1] = v[1] - rho
  return x
end

function newton_exp_onz(rho, y_hat, z_hat)
  t = max(-z_hat,EXP_PROJ_TOL)
  for iter=1:EXP_PROJ_MAXIT
    f = (1.0/rho^2)*t*(t + z_hat) - y_hat/rho + log(t/rho) + 1.0
    fp = (1.0/rho^2)*(2.0*t + z_hat) + 1.0/t
    t = t - f/fp
    if t <= -z_hat
      t = -z_hat
      break
    elseif t <= 0
      t = 0
      break
    elseif abs(f) <= EXP_PROJ_TOL
      break
    end
  end
  z = t + z_hat
  return z
end

fun_name(f::IndExpPrimal) = "indicator exponential cone (primal)"
fun_dom(f::IndExpPrimal) = "AbstractArray{Real}"
fun_expr(f::IndExpPrimal) = "x ↦ 0 if x ∈ cl{(r,s,t) : s > 0, s*exp(r/s) ⩽ t}, +∞ otherwise"
fun_params(f::IndExpPrimal) = "none"

fun_name{R <: Real}(f::Precompose{Conjugate{IndExpPrimal}, R}) = "indicator of the exponential cone (dual)"
fun_expr{R <: Real}(f::Precompose{Conjugate{IndExpPrimal}, R}) = "x ↦ 0 if x ∈ cl{(u,v,w) : u < 0, -u*exp(v/u) ⩽ w*exp(1)}, +∞ otherwise"
fun_params{R <: Real}(f::Precompose{Conjugate{IndExpPrimal}, R}) = "none"

prox_naive{R <: Real}(f::IndExpPrimal, x::AbstractArray{R}, gamma::Real=1.0) =
  prox(f, x, gamma) # we don't have a much simpler way to do this yet

prox_naive{R <: Real}(f::Precompose{Conjugate{IndExpPrimal}}, x::AbstractArray{R}, gamma::Real=1.0) =
  prox(f, x, gamma) # we don't have a much simpler way to do this yet


text:
using ProjAssistant
@quickactivate

# ------------------------------------------------------------------
@time begin

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014
    const Fd = ChF.FolsomData # experimental data

    using Plots
    import GR
    !isinteractive() && GR.inline("png")
end

## -------------------------------------------------------------------
# D vs X
let
    p = plot(; title = "Kayser",
        xlabel = string("X (", Fd.unit(:X), ")"),
        ylabel = string("D (", Fd.unit(:D), ")")
    )
    plot!(p, Fd.val(:D), Fd.val(:X);
        label = "", ls = :dash, lw = 3, color = :black,
        alpha = 0.6
    )
    scatter!(p, Fd.val(:D), Fd.val(:X);
        label = "", m = 8, color = :black
    )
    sfig(Fd, p,
        @fileid, "X_vs_D", ".png"
    )
end

## -------------------------------------------------------------------
# BALANCE
let
    ps = Plots.Plot[]
    for met in [:GLC]
        p = plot(; title = string("Balance: ", met),
            xlabel = "feed",
            ylabel = "exch + drain"
        )

        exps = 1:4

        feed = Fd.cval.(met, exps) .* Fd.val.(:D, exps)
        exch = Fd.uval.(met, exps) .* Fd.val.(:X, exps) .|> abs
        drain = zeros(length(exch)) # This experiments hase sg == 0

        scatter!(p, feed, exch .+ drain;
            label = "", m = 8, color = :black
        )
        vals = [feed; exch .+ drain] |> sort
        plot!(p, vals, vals;
            label = "", ls = :dash, lw = 3, alpha = 0.6, color = :black
        )
        push!(ps, p)
    end
    sfig(Fd, ps,
        @fileid, "Balances", ".png"
    )
end


text:
exp = Expr(:call, *, 3, 4) #> :((*)(3,4))
eval(exp) #> 12

# exp = Expr(:call, *, 3, a) #> ERROR: a not defined
exp = Expr(:call, *, 3, :a) #> :((*)(3,a))
# eval(exp) #> ERROR: a not defined
a = 4 #> 4
eval(exp) #> 12

# b #> ERROR: b not defined
exp = :(b = 1)
eval(exp) #> 1
b #> 1

exp1 = :(a + b) #> :(a + b)
exp2 = :($a + b) #> :(4 + b)
eval(exp1) #> 5
eval(exp2) #> 5

max_stars_repo_path: Module 1/Chapter07/Old/eval.jl

max_issues_repo_path: Module 1/Chapter07/Old/eval.jl

max_forks_repo_path: Module 1/Chapter07/Old/eval.jl


text:
using DiffEqGPU, CUDA, OrdinaryDiffEq, Test

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

CUDA.allowscalar(false)
u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = (10.0f0,28.0f0,8/3f0)
prob = ODEProblem(lorenz,u0,tspan,p)
const pre_p = [rand(Float32,3) for i in 1:10]
prob_func = (prob,i,repeat) -> remake(prob,p=pre_p[i].*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func)

@info "Explicit Methods"

#Performance check with nvvp
# CUDAnative.CUDAdrv.@profile
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)
@test length(filter(x -> x.u != sol.u[1].u, sol.u)) != 0 # 0 element array
@time sol = solve(monteprob,ROCK4(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)
@time sol2 = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,
                                                 batch_size=5,saveat=1.0f0)

@test length(filter(x -> x.u != sol.u[1].u, sol.u)) != 0 # 0 element array
@test length(filter(x -> x.u != sol2.u[6].u, sol.u)) != 0 # 0 element array
@test all(all(sol[i].prob.p .== pre_p[i].*p) for i in 1:10)
@test all(all(sol2[i].prob.p .== pre_p[i].*p) for i in 1:10)

@time solve(monteprob,Tsit5(),EnsembleCPUArray(),trajectories=10,saveat=1.0f0)
@time solve(monteprob,Tsit5(),EnsembleThreads(), trajectories=10,saveat=1.0f0)
@time solve(monteprob,Tsit5(),EnsembleSerial(),  trajectories=10,saveat=1.0f0)

#=
solve(monteprob,TRBDF2(),EnsembleCPUArray(),dt=0.1,trajectories=2,saveat=1.0f0)
solve(monteprob,TRBDF2(),EnsembleGPUArray(),dt=0.1,trajectories=2,saveat=1.0f0)
@test_broken solve(monteprob,TRBDF2(),EnsembleGPUArray(),dt=0.1,trajectories=2,saveat=1.0f0)
=#

@info "Implicit Methods"

function lorenz_jac(J,u,p,t)
    σ = p[1]
    ρ = p[2]
    β = p[3]
    x = u[1]
    y = u[2]
    z = u[3]
    J[1,1] = -σ
    J[2,1] = ρ - z
    J[3,1] = y
    J[1,2] = σ
    J[2,2] = -1
    J[3,2] = x
    J[1,3] = 0
    J[2,3] = -x
    J[3,3] = -β
end

function lorenz_tgrad(J,u,p,t)
    nothing
end

func = ODEFunction(lorenz,jac=lorenz_jac,tgrad=lorenz_tgrad)
prob_jac = ODEProblem(func,u0,tspan,p)
monteprob_jac = EnsembleProblem(prob_jac, prob_func = prob_func)

@time solve(monteprob_jac,Rodas5(),EnsembleCPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)
@time solve(monteprob_jac,Rodas5(),EnsembleGPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(),EnsembleCPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(),EnsembleGPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)

@info "Callbacks"

condition = function (u,t,integrator)
    @inbounds u[1] > 5
end

affect! = function (integrator)
    @inbounds integrator.u[1] = -4
end

# test discrete
discrete_callback = DiscreteCallback(condition,affect!,save_positions=(false,false))
callback_prob = ODEProblem(lorenz,u0,tspan,p,callback=discrete_callback)
callback_monteprob = EnsembleProblem(callback_prob, prob_func = prob_func)
@time solve(callback_monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)

c_condition = function (u,t,integrator)
    @inbounds u[1] - 3
end

c_affect! = function (integrator)
    @inbounds integrator.u[1] += 20
end

# test continuous
continuous_callback = ContinuousCallback(c_condition,c_affect!,save_positions=(false,false))
callback_prob = ODEProblem(lorenz,u0,tspan,p,callback=continuous_callback)
callback_monteprob = EnsembleProblem(callback_prob, prob_func = prob_func)
solve(callback_monteprob,Tsit5(),EnsembleGPUArray(),trajectories=2,saveat=1.0f0)

# test callback set
callback_set = CallbackSet(discrete_callback, continuous_callback)
callback_prob = ODEProblem(lorenz,u0,tspan,p,callback=callback_set)
callback_monteprob = EnsembleProblem(callback_prob, prob_func = prob_func)
solve(callback_monteprob,Tsit5(),EnsembleGPUArray(),trajectories=2,saveat=1.0f0)

# test merge
callback_prob = ODEProblem(lorenz,u0,tspan,p,callback=discrete_callback)
callback_monteprob = EnsembleProblem(callback_prob, prob_func = prob_func)
solve(callback_monteprob,Tsit5(),EnsembleGPUArray(),trajectories=2,saveat=1.0f0, callback=continuous_callback)

@info "ROBER"

#=
using OrdinaryDiffEq, LinearAlgebra, ParameterizedFunctions
LinearAlgebra.BLAS.set_num_threads(1)
rober = @ode_def begin
  dy₁ = -k₁*y₁+k₃*y₂*y₃
  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
=#

function rober_f(internal_var___du,internal_var___u,internal_var___p,t)
    @inbounds begin
        internal_var___du[1] = -(internal_var___p[1]) * internal_var___u[1] + internal_var___p[3] * internal_var___u[2] * internal_var___u[3]
        internal_var___du[2] = (internal_var___p[1] * internal_var___u[1] - internal_var___p[2] * internal_var___u[2] ^ 2) - internal_var___p[3] * internal_var___u[2] * internal_var___u[3]
        internal_var___du[3] = internal_var___p[2] * internal_var___u[2] ^ 2
    end
    nothing
end

function rober_jac(internal_var___J,internal_var___u,internal_var___p,t)
    @inbounds begin
        internal_var___J[1, 1] = -(internal_var___p[1])
        internal_var___J[1, 2] = internal_var___p[3] * internal_var___u[3]
        internal_var___J[1, 3] = internal_var___p[3] * internal_var___u[2]
        internal_var___J[2, 1] = internal_var___p[1] * 1
        internal_var___J[2, 2] = -2 * internal_var___p[2] * internal_var___u[2] - internal_var___p[3] * internal_var___u[3]
        internal_var___J[2, 3] = -(internal_var___p[3]) * internal_var___u[2]
        internal_var___J[3, 1] = 0 * 1
        internal_var___J[3, 2] = 2 * internal_var___p[2] * internal_var___u[2]
        internal_var___J[3, 3] = 0 * 1
    end
    nothing
end

function rober_tgrad(J,u,p,t)
 nothing
end

rober_prob = ODEProblem(ODEFunction(rober_f,jac=rober_jac,tgrad=rober_tgrad),
                        Float32[1.0,0.0,0.0],(0.0f0,1f5),(0.04f0,3f7,1f4))
sol = solve(rober_prob,Rodas5(),abstol=1f-8,reltol=1f-8)
sol = solve(rober_prob,TRBDF2(),abstol=1f-4,reltol=1f-1)
rober_monteprob = EnsembleProblem(rober_prob, prob_func = prob_func)

@time sol = solve(rober_monteprob,Rodas5(),
                  EnsembleGPUArray(),trajectories=10,saveat=1.0f0,abstol=1f-8,
                  reltol=1f-8)
@time sol = solve(rober_monteprob,TRBDF2(),
                  EnsembleGPUArray(),trajectories=10,saveat=1.0f0,abstol=1f-4,
                  reltol=1f-1)
@time sol = solve(rober_monteprob,TRBDF2(),EnsembleThreads(),trajectories=10,
                  abstol=1e-4,reltol=1e-1,saveat=1.0f0)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Ncurses_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Ncurses")
JLLWrappers.@generate_main_file("Ncurses", UUID("68e3532b-a499-55ff-9963-d1c0c0748b3a"))
end  # module Ncurses_jll


text:
# vcf.jl
# ======
#
# A submodule for reading, writing, and working with VCF formatted files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE

module VCF

import Automa
import Automa.RegExp: @re_str
import BioCore: BioCore, isfilled
import BioCore.Exceptions: missingerror
import BufferedStreams

include("record.jl")
include("metainfo.jl")
include("header.jl")
include("reader.jl")
include("writer.jl")

end


text:
# Currently there is an inefficiency in Julia v0.5 and possibly v0.6 with ntuple
# that will hopefully be fixed in the next version. To avoid
# slowdown, we implement a tuple constructor independent of Julia's
# metaprogramming.
#
# The method below is used for building the multi-dimensional
# derivative matrix:
@generated function make_cartesian_index(i::Int, arr1::CartesianIndex{1}, arr2::CartesianIndex{D}) where D
    levs = [:($q == i ? arr1[1] : arr2[$q]) for q=1:D]
    quote
        $(Expr(:meta, :inline))
        CartesianIndex{D}($(levs...))
    end
end

# Get the vector index for a 1D (level, cell, mode)
# Here l starts at 1
@inline function get_index_1D(k::Int, l::Int, c::Int, m::Int)
    return k * (1<<((l-2)%UInt) + (c-1)) + m
end

# Threshold method for full matrices, producing a sparse one
function threshold(mat::AbstractArray{T, 2}, atol = eps(T)) where T <: Real
    for i in eachindex(mat)
        val = mat[i]
        if abs(val) < atol
            mat[i] = 0
        end
    end
    return sparse(mat)
end

# Threshold method for sparse matrices
function threshold(mat::SparseMatrixCSC{T, Int}, atol = eps(T)) where T <: Real
    (I, J, V) = findnz(mat)
    for (index, val) in enumerate(V)
        if abs(val) < atol
            V[index] = zero(T)
        end
    end
    return dropzeros!(sparse(I, J, V))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Test skipchars for IOStreams
mktemp() do path, file
    function append_to_file(str)
        mark(file)
        print(file, str)
        flush(file)
        reset(file)
    end

    # test it doesn't error on eof
    @test eof(skipchars(file, isspace))

    # test it correctly skips
    append_to_file("    ")
    @test eof(skipchars(file, isspace))

    # test it correctly detects comment lines
    append_to_file("#    \n   ")
    @test eof(skipchars(file, isspace, linecomment='#'))

    # test it stops at the appropriate time
    append_to_file("   not a space")
    @test !eof(skipchars(file, isspace))
    @test read(file, Char) == 'n'

    # test it correctly ignores the contents of comment lines
    append_to_file("  #not a space \n   not a space")
    @test !eof(skipchars(file, isspace, linecomment='#'))
    @test read(file, Char) == 'n'

    # test it correctly handles unicode
    for (byte,char) in zip(1:4, ('@','߷','࿊','𐋺'))
        append_to_file("abcdef$char")
        @test Base.codelen(char) == byte
        @test !eof(skipchars(file, isalpha))
        @test read(file, Char) == char
    end
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
using CUFSM

prop = [100 29500.00 29500.00 0.30 0.30 11346.15]

node = [1 5.0000 1.0000 1 1 1 1 -38.889
2 5.0000 0.0000 1 1 1 1 -50.000
3 2.5000 0.0000 1 1 1 1 -50.000
4 0.0000 0.0000 1 1 1 1 -50.000
5 0.0000 3.0000 1 1 1 1 -16.667
6 0.0000 6.0000 1 1 1 1 16.667
7 0.0000 9.0000 1 1 1 1 50.000
8 2.5000 9.0000 1 1 1 1 50.000
9 5.0000 9.0000 1 1 1 1 50.000
10 5.0000 8.0000 1 1 1 1 38.889]

elem = [1 1 2 0.100000 100
2 2 3 0.100000 100
3 3 4 0.100000 100
4 4 5 0.100000 100
5 5 6 0.100000 100
6 6 7 0.100000 100
7 7 8 0.100000 100
8 8 9 0.100000 100
9 9 10 0.100000 100]

lengths = [1.00
2.00
3.00
4.00
5.00
6.00
7.00
8.00
9.00
10.00
11.00
12.00
13.00
14.00
15.00
20.00
30.00
40.00
50.00
60.00
70.00
80.00
90.00
100.00
200.00
300.00
400.00
500.00
600.00
700.00
800.00
900.00
1000.00]

springs = []

constraints = []

neigs = 20

curve, shapes = CUFSM.strip(prop, node, elem, lengths, springs, constraints, neigs)

max_stars_repo_name: runtosolve/CUFSM

max_issues_repo_name: runtosolve/CUFSM


text:
export div_update!

"""
    div_update!

Inplace solve (left or right matrix divide). Usage:

    div_update!(B,α,/,A)    # B = (αB)/A

    div_update!(B,/,A)      # B = B/A

    div_update!(B,α,\\,A)    # B = A\\(αB)

    div_update!(B,\\,A)      # B = A\\B

Matrix B is updated inplace. α is scalar. A must be triangular or simpler.
"""
function div_update!
end

for (div, side, fun) in ( ( :/, 'R', :(rdiv!(B,A)) ),
                          ( :\, 'L', :(ldiv!(A,B)) ) )
    @eval begin
        function div_update!(B::BlasMatrix{T}, α::Number, ::typeof($div), A::BlasTriangular{T}) where T
            try
                α = convert(T,α)
                BLAS.trsm!($side, A.uplo, A.trans, A.diag, α, A.blasnode, B)
            catch err
                ip_error(err)
            end
        end

        function div_update!(B::BlasArray, ::typeof($div), A)
            try
                $fun
            catch err
                ip_error(err)
            end
        end

    end
end

div_update!(B,α,f::Function,A) = ip_error("solve method not available for this combination of types .")
div_update!(B,f::Function,A) = ip_error("solve method not available for this combination of types .")


text:
import SynapseClient: Wiki

@testset "Wiki" begin
    # """Test the construction and accessors of Wiki objects."""

    #Wiki contstuctor only takes certain values
    @test_pythrows ValueError Wiki(title="foo")

    #Construct a wiki and test uri's
    wiki = Wiki(title ="foobar2", markdown="bar", owner=Dict("id"=>"5"))
end


text:
ex=:x^3
@test simplify(ex,Pow)==Pow(:x,3)
ex=:x*3*:x*:x
@test simplify(ex,Pow)==Expression(Any[3,Pow(:x,3)])
ex=:x*:y*3*:y*:x
@test simplify(ex,Pow)==Expression(Any[3,Pow(Expression(Any[:x,:y]),2)])
@test simplify(:x^3*:y^2*:x*3*:y,Pow)==Expression(Any[3,:x,Pow(Expression(Any[:x,:y]),3)]) #You are right dear Equations, that is simpler than the previous suggestion of Expression({3,Pow(Expression({:x,:x,:y}),2),:y})
@test simplify(:x^5,Pow)==Pow(:x,5)


text:
"""
Test function of a parabola
"""
function parabola(t, x; active::Bool=false)
    if active == true
        return [3*t^2]
    else
        return [0.0]
    end
end

"""
Simple test function for spherical earth dynamics
"""
function point_earth(epc::Epoch, x::Array{<:Real, 1})
    # Restrict inputs to position only. Considered in body frame
    r  = x[1:3]
    v  = x[4:6]

    # Acceleration
    a = -GM_EARTH * r/norm(r)^3

    return vcat(v, a)
end

# Test basic integration of parabola
let
    rk4 = RK4(parabola, active=false)
    @test istep(rk4, 0.0, 1.0, [0.0])[1] == 0.0

    rk4 = RK4(parabola, active=true)
    @test istep(rk4, 0.0, 1.0, [0.0])[1] == 1.0

    t = 0.0
    for i in 1:10
        t = istep(rk4, 0, i, [0.0])[1]
    end

    @test t == 1000
end

# Test stability of integration of conservative motion (point-mass orbit)
let
    rk4 = RK4(point_earth)

    epc0 = Epoch(2019, 1, 1)
    oe0  = [R_EARTH + 500e3, 0.01, 90.0, 0, 0, 0]
    eci0 = sOSCtoCART(oe0, use_degrees=true)

    epcf = epc0 + orbit_period(oe0[1])

    # Manually propagate orbit
    epc = deepcopy(epc0)
    eci = copy(eci0)
    while epc < epcf
        dt  = min(1.0, epcf-epc) # Set step-size
        eci = istep(rk4, epc, dt, eci)
        epc += dt
    end

    # Test that initial state equals intial propagated state
    @test isapprox(eci[1, 1], eci0[1], atol=1e-5)
    @test isapprox(eci[2, 1], eci0[2], atol=1e-5)
    @test isapprox(eci[3, 1], eci0[3], atol=1e-5)
    @test isapprox(eci[4, 1], eci0[4], atol=1e-5)
    @test isapprox(eci[5, 1], eci0[5], atol=1e-5)
    @test isapprox(eci[6, 1], eci0[6], atol=1e-5)

    oe = sCARTtoOSC(eci[:, 1], use_degrees=true)
    @test isapprox(oe[1], oe0[1], atol=1e-5)
    @test isapprox(oe[2], oe0[2], atol=1e-5)
    @test isapprox(oe[3], oe0[3], atol=1e-5)
    @test isapprox(oe[4], oe0[4], atol=1e-5)
    @test isapprox(oe[5], oe0[5], atol=1e-5)
    @test isapprox(oe[6], oe0[6], atol=1e-5)
end

let
    # Initialize integrator
    rk4 = RK4(point_earth)

    epc = Epoch(2019, 1, 1)
    eci = sOSCtoCART([R_EARTH + 500e3, 0.01, 90.0, 0, 0, 0], use_degrees=true)
    phi = diagm(0 => ones(Float64, length(eci)))

    # No step variational equations STM should be identity
    xu, phiu = istep(rk4, epc, 0.0, eci, phi)

    for i in 1:length(eci)
        @test phiu[i,i] == 1.0
    end

    # Small step STM should be nearly identity but non-zero
    xu, phiu = istep(rk4, epc, 1.0, eci, phi)

    for i in 1:length(eci)
        @test phiu[i, i] != 0.0
        @test isapprox(phiu[i, i], 1.0, atol=1.0e-5)
    end

    # Actually performing a perturbed step should equal the STM output
    pert = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    # Perturbed position from integrating pertturbed initial conditions
    ecip   = istep(rk4, epc, 1.0, eci + pert)

    # Pertrubed position from variational equation approximation
    ecistm = istep(rk4, epc, 1.0, eci) + phiu*pert

    # Assert equal
    for i in 1:6
        @test isapprox(ecip[i], ecistm[i], atol=1.0e-9)
    end
end


text:
using Test
using ScottishTaxBenefitModel
using .ModelHousehold
using .Results
using .Intermediate: MTIntermediate, make_intermediate
using .ScottishBenefits
using .STBIncomes
using .ExampleHelpers

sys = get_system( scotland=true )

@testset "Scottish Child Payment" begin

    for (hht,hh) in get_all_examples()
        bus = get_benefit_units( hh )
        intermed = make_intermediate(
            hh,
            sys.hours_limits,
            sys.age_limits,
            sys.child_limits )
        hres = init_household_result( hh )
        hhead = get_head( hh )
        hpid = hhead.pid
        for buno in eachindex(bus)
            bures = hres.bus[buno]
            bu = bus[buno]
            bint = intermed.buint[buno]
            ncs = num_children( bu )
            bhead = get_head( bus[buno ])
            bspouse = get_spouse( bus[buno ])
            if ncs > 0
                println( "on hhld $(hht)")
                # no qualifying benefit
                ages = fill( 6, ncs ) # 1 year too old
                set_childrens_ages!( hh, ages... )
                calc_scottish_child_payment!(
                    bures,
                    bu,
                    bint,
                    sys.scottish_child_payment
                )
                aggregate!( bures )
                @test ! has_any( bures, SCOTTISH_CHILD_PAYMENT )
                ti = rand( sys.scottish_child_payment.qualifying_benefits )
                bures.pers[bhead.pid].income[ti] = 1.0
                calc_scottish_child_payment!(
                    bures,
                    bu,
                    bint,
                    sys.scottish_child_payment
                )
                aggregate!( bures )
                # still should be zero since all over age
                @test ! has_any( bures, SCOTTISH_CHILD_PAYMENT )
                # all qualify on age
                ages = fill( 3, ncs )
                set_childrens_ages!( hh, ages... )
                calc_scottish_child_payment!(
                    bures,
                    bu,
                    bint,
                    sys.scottish_child_payment
                )
                aggregate!( bures )
                @test bures.income[SCOTTISH_CHILD_PAYMENT] == ncs*10
                # 1 qualifying
                set_childrens_ages!( hh, 3,8,9,10,12 )
                calc_scottish_child_payment!(
                    bures,
                    bu,
                    bint,
                    sys.scottish_child_payment
                )
                aggregate!( bures )
                @test bures.income[SCOTTISH_CHILD_PAYMENT] == 10
            end # with children
        end
    end

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Gtk

function flush_keyboard()
    win = GtkWindow("", 1, 1)
    keyget(w, event) = Int32(0)
    signal_connect(keyget, win, "key-press-event")
    visible(win, false)
    sleep(0.25)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# analyzer analyze proof space in core

##### STATIC Analysis
## Pairing

struct LIDPAIR
  plid
  nlid
end

function pairing(plist, nlist)
  pairs=[]
  for ps in plist
  for ns in nlist
    push!(pairs, LIDPAIR(ps, ns))
  end
  end
  pairs
end

function pairmaker(core)
 pairs = []
 for psym in core.allpsym
  pos = oppositof(:-, psym, core)
  neg = oppositof(:+, psym, core)

  append!(pairs, pairing(pos, neg))

 end
 pairs
end

function resolvepair(pn, core)
 plid = pn.plid
 nlid = pn.nlid
 varp = varsof(cidof(plid, core), core)
 varn = varsof(cidof(nlid, core), core)
 vars = vcat(varp, varn)
 plit = literalof(plid, core)
 nlit = literalof(nlid, core)
 mgu = unify(vars, plit.body.args[2], nlit.body.args[2])
 return pn, vars, mgu

end

function pairmakefromcore(path::String)
 core = readcore(path)
 pairs = pairmaker(core)
 core, map(pn->resolvepair(pn,core), pairs)
end

function pairmakefromcore(core::CORE)
 pairs = pairmaker(core)
 core, map(pn->resolvepair(pn,core), pairs)
end

function printvm(vs, ts)
  cflag=false
  for ix in 1:length(vs)
    v = vs[ix]
    t = ts[ix]
    if v==t; continue end
    cflag && print(",")
    cflag=true
    print("$v/$t")
  end
  println()
end

function printlform2(lfm)
  print("$(lfm.lid):$(lfm.body)")
end

function printpair(pair, core)
 if cidof(pair.plid,core) == cidof(pair.nlid,core)
   print("*")
 else
   print(" ")
 end
 print("<")
 printlform2(literalof(pair.plid, core))
 print("; ")
 printlform2(literalof(pair.nlid, core))
 println(">")
end

function printpvms(pvms, core)
  for pvm in pvms
    println("---")
    printpair(pvm[1], core)
    printvm(pvm[2], pvm[3])
  end
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using Distributions, Plots, CSV, Pandas

# Sample x_{k+1} given x_k
# Parameters:
# - x_k: The kth X sample. One of -1, 0, or 1
# Outputs:
# - x_{k+1}: The (k+1)th X sample
function sample_x(xk)
    # Define q parameter
    q = 0.3

    # Probability that x_{k+1} = 1
    p = q*(xk == 1) + (1-q)*(xk == 0)

    return rand() < p
end

# Obtains n samples of X based on HMM
# Inputs
# - n: An integer representing the number of samples to return
# Outputs
# A vector of length n where the kth element is the kth sample of X
function get_x_samples(n)
    # Store x samples
    samples = zeros(n)

    # Randomly set x1
    samples[1] = rand(0:1)

    # Get samples x_2 through x_n
    for i=2:n
        samples[i] = sample_x(samples[i-1])
    end

    return samples
end

# Obtains n samples of Y based on HMM
# Inputs
# - x_samples: A vector of length n representing the samples of X
# Outputs
# A vector of length n where the kth element is the kth sample of Y
function get_y_samples(x_samples)
    # Define normal standard deviation
    σ = 1.5

    # Store y samples
    n = length(x_samples)
    y_samples = zeros(n)

    for i=1:n
        𝒩 = Normal(x_samples[i], σ)
        y_samples[i] = rand(𝒩)
    end

    return y_samples
end

# Obtian and plot samples
n = 1500
x_samples = get_x_samples(n)
y_samples = get_y_samples(x_samples)

plot(1:n,x_samples, seriestype=:scatter, bg = RGB(247/255, 236/255, 226/255), color = RGB(0,191/255,255/255), label = "x", alpha = 0.2)
plot!(1:n,y_samples, seriestype=:scatter, color = RGB(191/255,1,0), label =
"y", alpha = 0.2)

cd("/Users/elvis/Documents/DSI/courses/content/stochastic-approximations/images/")

# Note X and Y variables swapped here from what was presented in example
df = DataFrame(Dict(:X=>y_samples))

CSV.write("/Users/elvis/Documents/DSI/courses/content/bayesian-inference-and-graphical-models/code/hmm_observations.csv", df)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using Distributions: UnivariateDistribution, ContinuousUnivariateDistribution, DiscreteUnivariateDistribution
using LinearAlgebra: Cholesky, cholesky

abstract type AbstractCorrelation end
struct Pearson <: AbstractCorrelation end
struct Spearman <: AbstractCorrelation end
struct Kendall <: AbstractCorrelation end

struct CorMat{T<:Real, S<:AbstractMatrix, C<:Union{AbstractCorrelation, Nothing}} <: AbstractMatrix{T}
    mat::S
    chol::Cholesky{T,S}
    type::C
end
CorMat(mat::Matrix) = CorMat(mat, cholesky(mat), nothing)

r = [1.0 0.5; 0.5 1.0]
CorMat(r, cholesky(r), Pearson)
CorMat(r)

struct MvDist{T<:Real, S<:AbstractMatrix, C<:AbstractCorrelation}
    margins::Vector{<:UnivariateDistribution}
    target_cor::CorMat{T, S, C}
    adjust_cor::CorMat{T, S, Pearson}

    function MvDist{T,S,C}(
        m::Vector{<:UnivariateDistribution},
        t::CorMat{T,S,C},
        a::CorMat{T,S,Pearson}
    ) where {T,S,C}
        return new{T,S,C}(m, t, a)
    end
end


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

using EuclidianNormalizingFlows
using Test

using LinearAlgebra
using InverseFunctions, ChangesOfVariables
using ForwardDiff, Zygote

using EuclidianNormalizingFlows: householder_trafo, chained_householder_trafo, HouseholderTrafo

@testset "householder_trafo" begin
    v = rand(5)
    x = rand(5)
    X = rand(5, 3)

    householder_matrix(v::AbstractVector{<:Real}) = I - 2 * (v*v') / (v'*v)

    @test @inferred(householder_trafo(v, x)) ≈ householder_matrix(v) * x
    @test @inferred(householder_trafo(v, householder_trafo(v, x))) ≈ x

    @test @inferred(householder_trafo(v, X)) ≈ hcat(householder_trafo.(Ref(v), eachcol(X))...)
    @test @inferred(householder_trafo(v, X)) ≈ householder_matrix(v) * X
    @test @inferred(householder_trafo(v, householder_trafo(v, X))) ≈ X

    @test Zygote.pullback(householder_trafo, v, x)[1] == householder_trafo(v, x)
    @test Zygote.jacobian(householder_trafo, v, x)[1] ≈ ForwardDiff.jacobian(v -> householder_matrix(v) * x, v)
    @test Zygote.jacobian(householder_trafo, v, x)[2] ≈ ForwardDiff.jacobian(x -> householder_matrix(v) * x, x)

    @test Zygote.pullback(householder_trafo, v, X)[1] == householder_trafo(v, X)
    @test Zygote.jacobian(householder_trafo, v, X)[1] ≈ ForwardDiff.jacobian(v -> householder_matrix(v) * X, v)
    @test Zygote.jacobian(householder_trafo, v, X)[2] ≈ ForwardDiff.jacobian(X -> householder_matrix(v) * X, X)

    V = rand(5, 3)

    chained_householder_matrix(v::AbstractMatrix{<:Real}) = *(reverse(householder_matrix.(eachcol(V)))...)

    @test chained_householder_trafo(V, x) ≈ *(reverse(householder_matrix.(eachcol(V)))...) * x
    @test chained_householder_trafo(reverse(V, dims = 2), chained_householder_trafo(V, x)) ≈ x

    @test @inferred(chained_householder_trafo(V, X)) ≈ hcat(chained_householder_trafo.(Ref(V), eachcol(X))...)

    # # Yields incorrect results:
    # ForwardDiff.jacobian(V -> chained_householder_matrix(V) * x, V)
    # Zygote.jacobian(V -> chained_householder_matrix(V) * x, V)

    @test Zygote.pullback(chained_householder_trafo, V, x)[1] == chained_householder_trafo(V, x)
    @test Zygote.jacobian(chained_householder_trafo, V, x)[1] ≈ ForwardDiff.jacobian(V -> chained_householder_trafo(V, x), V)
    @test Zygote.jacobian(chained_householder_trafo, V, x)[2] ≈ ForwardDiff.jacobian(x -> chained_householder_trafo(V, x), x)

    @test Zygote.pullback(chained_householder_trafo, V, X)[1] == chained_householder_trafo(V, X)
    @test Zygote.jacobian(chained_householder_trafo, V, X)[1] ≈ ForwardDiff.jacobian(V -> chained_householder_trafo(V, X), V)
    @test Zygote.jacobian(chained_householder_trafo, V, X)[2] ≈ ForwardDiff.jacobian(X -> chained_householder_trafo(V, X), X)

    @test HouseholderTrafo(V) == HouseholderTrafo(V)
    @test isequal(HouseholderTrafo(V), HouseholderTrafo(V))
    @test hash(HouseholderTrafo(V)) == hash(HouseholderTrafo(V))

    for arg in (x, X)
        InverseFunctions.test_inverse(HouseholderTrafo(V), arg)
    end

    @test @inferred(with_logabsdet_jacobian(HouseholderTrafo(V), x)) == (chained_householder_trafo(V, x), 0)
    @test @inferred(with_logabsdet_jacobian(HouseholderTrafo(V), X)) == (chained_householder_trafo(V, X), fill(0, size(X, 2))')
end


text:
#######################################################################
# A demo of a clustered category chart in XlsxWriter.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_clustered.xlsx")
	ws = add_worksheet!(wb)
	bold = add_format!(wb, Dict("bold" => 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Types" "Sub Type" "Value 1" "Value 2" "Value 3"]

	data = [
		["Type 1" "Sub Type A" 5000      8000      6000];
		[""       "Sub Type B" 2000      3000      4000];
		[""       "Sub Type C" 250       1000      2000];
		["Type 2" "Sub Type D" 6000      6000      6500];
		[""       "Sub Type E" 500       300        200];
	]

	write_row!(ws, "A1", headings, fmt=bold)

	for row_num in 1:size(data,1)
		write_row!(ws, row_num, 0, data[row_num,:])
	end
	# Create a new chart object. In this case an embedded chart.
	chart = add_chart!(wb, Dict("type" => "column"))

	# Configure the series. Note, that the categories are 2D ranges (from column A
	# to column B). This creates the clusters. The series are shown as formula
	# strings for clarity but you can also use the list syntax. See the docs.
	add_series!(chart, Dict(
		"name" => "=Sheet1!\$C\$1",
		"categories" => "=Sheet1!\$A\$2:\$B\$6",
		"values" => "=Sheet1!\$C\$2:\$C\$6",
	))

	add_series!(chart, Dict(
		"name" => "=Sheet1!\$D\$1",
		"categories" => "=Sheet1!\$A\$2:\$B\$6",
		"values" => "=Sheet1!\$D\$2:\$D\$6",
	))

	add_series!(chart, Dict(
		"name" => "=Sheet1!\$E\$1",
		"categories" => "=Sheet1!\$A\$2:\$B\$6",
		"values" => "=Sheet1!\$E\$2:\$E\$6",
	))

	# Set the Excel chart style.
	set_style!(chart, 37)

	# Turn off the legend.
	set_legend!(chart, Dict("position" => "none"))

	# Insert the chart into the worksheet.
	insert_chart!(ws, "G3", chart)

	close(wb)

	isfile("chart_clustered.xlsx")
end

test()


text:
"""
```julia
vvar(A; dims=:, mean=nothing, corrected=true)
```
Compute the variance of all elements in `A`, optionally over dimensions specified by `dims`.
As `Statistics.var`, but vectorized and (optionally) multithreaded.

A precomputed `mean` may optionally be provided, which results in a somewhat faster
calculation. If `corrected` is `true`, then _Bessel's correction_ is applied, such
that the sum is divided by `n-1` rather than `n`.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> vvar(A, dims=1)
1×2 Matrix{Float64}:
 2.0  2.0

julia> vvar(A, dims=2)
2×1 Matrix{Float64}:
 0.5
 0.5
```
"""
function vvar(A; dims=:, mean=nothing, corrected=true, multithreaded=:auto)
    if (multithreaded===:auto && length(A) > 4095) || multithreaded===true
        _vtvar(mean, corrected, A, dims)
    else
        _vvar(mean, corrected, A, dims)
    end
end
export vvar

## Singlethreaded implementation
# If dims is an integer, wrap it in a tuple
_vvar(μ, corrected::Bool, A, dims::Int) = _vvar(μ, corrected, A, (dims,))

# If the mean isn't known, compute it
_vvar(::Nothing, corrected::Bool, A, dims::Tuple) = _vvar!(_vmean(A, dims), corrected, A, dims)
function _vvar(::Nothing, corrected::Bool, A, ::Colon)
    # Reduce all the dims!
    n = length(A)
    Tₒ = Base.promote_op(/, eltype(A), Int)
    Σ = zero(Tₒ)
    @turbo for i ∈ eachindex(A)
            Σ += A[i]
    end
    μ = Σ / n
    σ² = zero(typeof(μ))
    @turbo for i ∈ eachindex(A)
            δ = A[i] - μ
            σ² += δ * δ
    end
    return σ² / (n-corrected)
end

# If the mean is known, pass it on in the appropriate form
_vvar(μ, corrected::Bool, A, dims::Tuple) = _vvar!(collect(μ), corrected, A, dims)
_vvar(μ::Array, corrected::Bool, A, dims::Tuple) = _vvar!(copy(μ), corrected, A, dims)
_vvar(μ::Number, corrected::Bool, A, dims::Tuple) = _vvar!([μ], corrected, A, dims)
function _vvar(μ::Number, corrected::Bool, A, ::Colon)
    # Reduce all the dims!
    n = length(A)
    σ² = zero(typeof(μ))
    @turbo for i ∈ eachindex(A)
        δ = A[i] - μ
        σ² += δ * δ
    end
    return σ² / (n-corrected)
end

# Chris Elrod metaprogramming magic:
# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_var_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the var over.
  Bᵥ = Expr(:call, :view, :B)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(A[])
  Bind = :(Bᵥ[])
  inds = Vector{Symbol}(undef, N)
  len = Expr(:call, :*)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    push!(Aind.args, ind)
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Bᵥ.args, :(firstindex(B,$n)))
      push!(len.args, :(size(A, $n)))
    else
      push!(nonreduct_inds, n)
      push!(Bᵥ.args, :)
      push!(Bind.args, ind)
    end
  end
  # Secondly, build up our set of loops
  if !isempty(nonreduct_inds)
    firstn = first(nonreduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = indices((A,B),$firstn)), block)
    if length(nonreduct_inds) > 1
      for n ∈ @view(nonreduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = indices((A,B),$n)), newblock))
        block = newblock
      end
    end
    rblock = block
    # Push more things here if you want them at the beginning of the reduction loop
    push!(rblock.args, :(μ = $Bind))
    push!(rblock.args, :(σ² = zero(eltype(Bᵥ))))
    # Build the reduction loop
    for n ∈ reduct_inds
      newblock = Expr(:block)
      push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
      block = newblock
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(δ = $Aind - μ))
    push!(block.args, :(σ² += δ * δ))
    # Push more things here if you want them at the end of the reduction loop
    push!(rblock.args, :($Bind = σ² * invdenom))
    # Put it all together
    return quote
      invdenom = inv(($len) - corrected)
      Bᵥ = $Bᵥ
      @turbo $loops
      return B
    end
  else
    firstn = first(reduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = axes(A,$firstn)), block)
    if length(reduct_inds) > 1
      for n ∈ @view(reduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
        block = newblock
      end
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(δ = $Aind - μ))
    push!(block.args, :(σ² += δ * δ))
    # Put it all together
    return quote
      invdenom = inv(($len) - corrected)
      Bᵥ = $Bᵥ
      μ = Bᵥ[]
      σ² = zero(eltype(Bᵥ))
      @turbo $loops
      Bᵥ[] = σ² * invdenom
      return B
    end
  end
end

# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_var_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return _vvar!(B, corrected, A, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      # Else, if dimm ∉ 1:N, drop it from list and continue
      tc = copy(t)
      for r ∈ m+1:M
        push!(tc.args, :(dims[$r]))
      end
      push!(qold.args, Expr(:block, :(return _vvar!(B, corrected, A, $tc))))
      return q
    end
  end
  staticdim_var_quote(static_dims, N)
end

# Efficient @generated in-place var
@generated function _vvar!(B::AbstractArray{Tₒ,N}, corrected::Bool, A::AbstractArray{T,N}, dims::D) where {Tₒ,T,N,M,D<:Tuple{Vararg{Integer,M}}}
  branches_var_quote(N, M, D)
end
@generated function _vvar!(B::AbstractArray{Tₒ,N}, corrected::Bool, A::AbstractArray{T,N}, dims::Tuple{}) where {Tₒ,T,N}
  :(fill!(B, Tₒ(NaN)); return B)
end

## Multithreaded implementation
# If dims is an integer, wrap it in a tuple
_vtvar(μ, corrected::Bool, A, dims::Int) = _vtvar(μ, corrected, A, (dims,))

# If the mean isn't known, compute it
_vtvar(::Nothing, corrected::Bool, A, dims::Tuple) = _vtvar!(_vtmean(A, dims), corrected, A, dims)
function _vtvar(::Nothing, corrected::Bool, A, ::Colon)
    # Reduce all the dims!
    n = length(A)
    Tₒ = Base.promote_op(/, eltype(A), Int)
    Σ = zero(Tₒ)
    @tturbo for i ∈ eachindex(A)
            Σ += A[i]
    end
    μ = Σ / n
    σ² = zero(typeof(μ))
    @tturbo for i ∈ eachindex(A)
            δ = A[i] - μ
            σ² += δ * δ
    end
    return σ² / (n-corrected)
end

# If the mean is known, pass it on in the appropriate form
_vtvar(μ, corrected::Bool, A, dims::Tuple) = _vtvar!(collect(μ), corrected, A, dims)
_vtvar(μ::Array, corrected::Bool, A, dims::Tuple) = _vtvar!(copy(μ), corrected, A, dims)
_vtvar(μ::Number, corrected::Bool, A, dims::Tuple) = _vtvar!([μ], corrected, A, dims)
function _vtvar(μ::Number, corrected::Bool, A, ::Colon)
    # Reduce all the dims!
    n = length(A)
    σ² = zero(typeof(μ))
    @tturbo for i ∈ eachindex(A)
        δ = A[i] - μ
        σ² += δ * δ
    end
    return σ² / (n-corrected)
end

# Chris Elrod metaprogramming magic:
# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_tvar_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the var over.
  Bᵥ = Expr(:call, :view, :B)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(A[])
  Bind = :(Bᵥ[])
  inds = Vector{Symbol}(undef, N)
  len = Expr(:call, :*)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    push!(Aind.args, ind)
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Bᵥ.args, :(firstindex(B,$n)))
      push!(len.args, :(size(A, $n)))
    else
      push!(nonreduct_inds, n)
      push!(Bᵥ.args, :)
      push!(Bind.args, ind)
    end
  end
  # Secondly, build up our set of loops
  if !isempty(nonreduct_inds)
    firstn = first(nonreduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = indices((A,B),$firstn)), block)
    if length(nonreduct_inds) > 1
      for n ∈ @view(nonreduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = indices((A,B),$n)), newblock))
        block = newblock
      end
    end
    rblock = block
    # Push more things here if you want them at the beginning of the reduction loop
    push!(rblock.args, :(μ = $Bind))
    push!(rblock.args, :(σ² = zero(eltype(Bᵥ))))
    # Build the reduction loop
    for n ∈ reduct_inds
      newblock = Expr(:block)
      push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
      block = newblock
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(δ = $Aind - μ))
    push!(block.args, :(σ² += δ * δ))
    # Push more things here if you want them at the end of the reduction loop
    push!(rblock.args, :($Bind = σ² * invdenom))
    # Put it all together
    return quote
      invdenom = inv(($len) - corrected)
      Bᵥ = $Bᵥ
      @tturbo $loops
      return B
    end
  else
    firstn = first(reduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = axes(A,$firstn)), block)
    if length(reduct_inds) > 1
      for n ∈ @view(reduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
        block = newblock
      end
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(δ = $Aind - μ))
    push!(block.args, :(σ² += δ * δ))
    # Put it all together
    return quote
      invdenom = inv(($len) - corrected)
      Bᵥ = $Bᵥ
      μ = Bᵥ[]
      σ² = zero(eltype(Bᵥ))
      @tturbo $loops
      Bᵥ[] = σ² * invdenom
      return B
    end
  end
end

# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_tvar_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return _vtvar!(B, corrected, A, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      # Else, if dimm ∉ 1:N, drop it from list and continue
      tc = copy(t)
      for r ∈ m+1:M
        push!(tc.args, :(dims[$r]))
      end
      push!(qold.args, Expr(:block, :(return _vtvar!(B, corrected, A, $tc))))
      return q
    end
  end
  staticdim_tvar_quote(static_dims, N)
end

# Efficient @generated in-place var
@generated function _vtvar!(B::AbstractArray{Tₒ,N}, corrected::Bool, A::AbstractArray{T,N}, dims::D) where {Tₒ,T,N,M,D<:Tuple{Vararg{Integer,M}}}
  branches_tvar_quote(N, M, D)
end
@generated function _vtvar!(B::AbstractArray{Tₒ,N}, corrected::Bool, A::AbstractArray{T,N}, dims::Tuple{}) where {Tₒ,T,N}
  :(fill!(B, Tₒ(NaN)); return B)
end


text:
## Pressure Field Demonstration

## Preamble
using Plots

include("../AcousticPropagation.jl")

## Lloyds Mirror

## Convergence Zone
cVec = [1520, 1500, 1515, 1495, 1545.]
zVec = [0., 300., 1200., 2e3, 5000.]

cMat = vcat([0 0 250e3], hcat(zVec, cVec, cVec))

ati = AcousticPropagation.Boundary(0)
bty = AcousticPropagation.Boundary(5e3)
ocn = AcousticPropagation.Medium(cMat, 250e3)
src = AcousticPropagation.Entity(0, 0)
θ_crit = acos(ocn.c(0, 0)/ocn.c(0, 5e3))
# θ₀ = range(0.01θ_crit, θ_crit, length = 10)
θ₀ = θ_crit

ray = AcousticPropagation.Ray(θ₀, src, ocn, bty, ati)

pt = plot(yaxis = :flip)
plot!(ray.Sol, vars = (1, 2), label = "")
display(pt)

##
s = range(0, ray.S, length = 1000)
plot(s, ray.A₀)


text:
# This routine was made by Leandro Martinez https://github.com/leandromartinez98

using DelimitedFiles
include("./align.jl")

#a = readdlm("./1.xyz")
#b = readdlm("./2.xyz")

function rmsd(x,y)
	# Align x to y (same size, correspondence given by sequence)

	z = align(x,y)

	# Compute rmsd before and after alignment

	rmsd_before = 0.
	rmsd_after = 0.
	n = size(x,1)
	for i in size(x,1)
#	  global rmsd_before, rmsd_after, n
	  rmsd_before = rmsd_before + (x[i,1] - y[i,1])^2 + (x[i,2] - y[i,2])^2 + (x[i,3] - y[i,3])^2
	  rmsd_after = rmsd_after + (z[i,1] - y[i,1])^2 + (z[i,2] - y[i,2])^2 + (z[i,3] - y[i,3])^2
	end
	rmsd_before = sqrt(rmsd_before/n)
	rmsd_after = sqrt(rmsd_after/n)
#	println(" rmsd before = ", rmsd_before)
#	println(" rmsd after = ", rmsd_after)
	return [rmsd_before,rmsd_after]
end


text:
using StingerGraphs
using BenchmarkTools
using JLD

const seed = 10^4

function generate!(edges, numvertices::Integer, numedges::Integer, rng=srand(seed))
    for i in 1:numedges
        src, dst = rand(rng, 0:numvertices-1), rand(rng, 0:numvertices-1)
        while src==dst
            dst = rand(rng, 0:numvertices-1)
        end
        edges[:,i] = [0, src, dst, 0, 0]
    end
end

function setup_stinger(edges, numedges)
    s = Stinger()
    insert_edges!(s, edges, numedges)
    return s
end

function insert_remove(s, edges, numedges)
    insert_edges!(s, edges, numedges)
    remove_edges!(s, edges, numedges)
end

function bench(numvertices=10000, numedges=1000000, filename="add_remove_bench.jld")

    edges = zeros(Int64, 5, numedges)
    generate!(edges, numvertices, numedges)
    removals = unique(edges, 2)
    numremovals = size(removals, 2)
    info("$(size(removals,2)) unique edges being inserted and removed")

    insert_bench = @benchmarkable insert_edges!(s, $edges, $numedges) seconds=10 setup=(s=Stinger())
    remove_bench = @benchmarkable remove_edges!(s, $removals, $numremovals) seconds=10 setup=(s=setup_stinger($edges, $numedges))
    insert_remove_bench = @benchmarkable insert_remove(s, $removals, $numremovals) seconds=10 setup=(s=Stinger())

    info("Running insert benchmark")
    insert_trial = run(insert_bench)
    info("Running removals benchmark")
    remove_trial = run(remove_bench)
    info("Running insert and remove benchmark")
    insert_remove_trial = run(insert_remove_bench)

    jldopen(filename, "w") do f
        write(f, "insert_trial", insert_trial)
        write(f, "remove_trial", remove_trial)
        write(f, "insert_remove_trial", insert_remove_trial)
    end

    return insert_trial, remove_trial, insert_remove_trial
end


text:
using CausalInference
using Graphs
using Test

@testset "chickering" begin
    global g = randdag(200)
    @test collect(CausalInference.ordered_edges(g)) == map(Edge,  CausalInference.chickering_order(g))
end

# http://www.multimedia-computing.de/mediawiki/images/5/55/SS08_BN-Lec2-BasicProbTheory_3.pdf
E1a = [1=>2, 1=>3, 2=>4, 3=>4]
E1b = [2=>1, 1=>3, 2=>4, 3=>4]
E1c = [1=>2, 3=>1, 2=>4, 3=>4]
C1 = [1=>2, 2=>1, 1=>3, 3=>1, 2=>4, 3=>4]

# Corrected from
E2a = [1=>2, 2=>3, 3=>5, 1=>4, 4=>5]
C2 = [1=>2, 2=>1, 2=>3, 3=>2, 3=>5, 1=>4, 4=>1, 4=>5]

E3 = [6=>4, 4=>3, 3=>2, 3=>1, 2=>1, 3=>5]
C3 = copy(E3)
append!(C3, reverse.(E3))

h, S = skel_oracle(digraph(E1a))
@test vpairs(h) == E1a
@test length(S) == 2
@test S[Edge(2, 3)] == [1]
@test S[Edge(1, 4)] == [2, 3]

h = pc_oracle(digraph(E1a))
@test sort(C1) == sort(vpairs(h))
@test digraph(C1) == cpdag(digraph(E1a))
@test digraph(C1) == cpdag(digraph(E1b))
@test digraph(C1) == cpdag(digraph(E1c))

for (C, Es) in [(C1, [E1a, E1b, E1c]),
                (C2, [E2a]),
                (C3, [E3]),
    ]
    for E in Es
        global h = pc_oracle(digraph(E))
        @test sort(C) == sort(vpairs(h))
        @test digraph(C) == cpdag(digraph(E))
    end
end

for n in 0:10

    alpha = 0.1
    @testset "randdag($n)" begin
        for k = 1:1000
            global g = randdag(n, alpha)
            h1 = pc_oracle(g)
            h2 = cpdag(g)
            h1 == h2  || println(vpairs(g))
            @test vpairs(h1) ⊆ vpairs(h2)
            @test vpairs(h2) ⊆ vpairs(h1)

            #@test h1 == h2
        end
    end
end

@testset "randdag(50)" begin # up to level 10 or so
    n = 50
    global g = randdag(n, 0.07)
    h1 = pc_oracle(g)
    h2 = cpdag(g)
    h1 == h2  || println(vpairs(g))
    @test vpairs(h1) ⊆ vpairs(h2)
    @test vpairs(h2) ⊆ vpairs(h1)
end

@testset "orient after v structre" begin
    global g = digraph([1=>4, 2=>1, 3=>1])
    @test sort(vpairs(vskel(g))) ==  [1=>4, 2=>1, 3=>1, 4=>1]

    h1 = pc_oracle(g)
    h2 = cpdag(g)
    @test h1 == h2
end

@testset "check that edge is still undirected" begin
    global g = digraph([1=>4, 2=>1, 3=>1, 3=>4])
    @test sort(vpairs(vskel(g))) == [1=>4, 2=>1, 3=>1, 3=>4, 4=>1, 4=>3]

    h1 = pc_oracle(g)
    h2 = cpdag(g)
    @test h1 == h2
end


text:
#####
##### generate the example for the README
#####
##### NOTE run with its own current directory, use `make`

using BrowseTables, Tables
# make example table, but any table that supports Tables.jl will work
table = Tables.columntable(collect(i == 5 ? (a = missing, b = "string", c = nothing) :
                                   (a = i, b = Float64(i), c = 'a'-1+i) for i in 1:10))
file = "readme_example.html"
write_html_table(file, table)


text:
using ADCME
using AdFem

include("../chip_unstructured_solver.jl")
include("../chip_unstructured_geometry.jl")

k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

function k_exact(x, y)
    1 + 1 / (1+x^2)
end

k_chip = eval_f_on_dof_pts(k_exact, mesh)[chip_fem_idx]

nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

NT = 15    # number of iterations for Newton's method

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)

B = constant(compute_interaction_matrix(mesh))
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
      bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

S0 = constant(zeros(nelem+3*ndof))
S = solve_navier_stokes(S0, NT, k_chip)

sess = Session(); init(sess)
output = run(sess, S)

matwrite("data.mat",
    Dict(
        "V"=>output[end, :]
    ))

# TODO: fix plot
u_out, v_out, p_out, T_out = output[NT+1,1:nnode], output[NT+1,ndof+1:ndof+nnode],
                             output[NT+1,2*ndof+1:2*ndof+nelem],output[NT+1,2*ndof+nelem+1:2*ndof+nelem+nnode]

figure(figsize=(10,10))
subplot(221)
title("u velocity")
visualize_scalar_on_fem_points(u_out .* u_std, mesh);#gca().invert_yaxis()
subplot(222)
title("v velocity")
visualize_scalar_on_fem_points(v_out .* u_std, mesh);#gca().invert_yaxis()
subplot(223)
visualize_scalar_on_fvm_points(p_out .* p_std, mesh);#gca().invert_yaxis()
title("pressure")
subplot(224)
title("temperature")
visualize_scalar_on_fem_points(T_out.* T_infty .+ T_infty, mesh);#gca().invert_yaxis()
tight_layout()
savefig("forward_solution_unstructured.pdf")

print("Solution range:",
    "\n [u velocity] \t min:", minimum(u_out .* u_std), ",\t max:", maximum(u_out .* u_std),
    "\n [v velocity] \t min:", minimum(v_out .* u_std), ",\t max:", maximum(v_out .* u_std),
    "\n [pressure]   \t min:", minimum(p_out .* p_std), ",\t max:", maximum(p_out .* p_std),
    "\n [temperature]\t min:", minimum(T_out.* T_infty .+ T_infty), ",\t\t\t max:", maximum(T_out.* T_infty .+ T_infty))


text:
# Core traits and accessor functions

"""
`alpha(p)` extracts the alpha component of a color. For a color
without an alpha channel, it will always return 1.
"""
alpha(c::TransparentColor) = c.alpha
alpha(c::Color)   = oneunit(eltype(c))
alpha(c::RGB24)   = N0f8(1)
alpha(c::ARGB32)  = N0f8((c.color & 0xff000000)>>24, 0)
alpha(c::AGray32) = N0f8((c.color & 0xff000000)>>24, 0)

"`red(c)` returns the red component of an `AbstractRGB` opaque or transparent color."
red(c::AbstractRGB   ) = c.r
red(c::TransparentRGB) = c.r
red(c::RGB24)  = N0f8((c.color & 0x00ff0000)>>16, 0)
red(c::ARGB32) = N0f8((c.color & 0x00ff0000)>>16, 0)

"`green(c)` returns the green component of an `AbstractRGB` opaque or transparent color."
green(c::AbstractRGB   ) = c.g
green(c::TransparentRGB) = c.g
green(c::RGB24)  = N0f8((c.color & 0x0000ff00)>>8, 0)
green(c::ARGB32) = N0f8((c.color & 0x0000ff00)>>8, 0)

"`blue(c)` returns the blue component of an `AbstractRGB` opaque or transparent color."
blue(c::AbstractRGB   ) = c.b
blue(c::TransparentRGB) = c.b
blue(c::RGB24)  = N0f8(c.color & 0x000000ff, 0)
blue(c::ARGB32) = N0f8(c.color & 0x000000ff, 0)

"`gray(c)` returns the gray component of a grayscale opaque or transparent color."
gray(c::Gray)    = c.val
gray(c::TransparentGray) = c.val
gray(c::Gray24)  = N0f8(c.color & 0x000000ff, 0)
gray(c::AGray32) = N0f8(c.color & 0x000000ff, 0)
gray(x::Number)  = x

"""
`chroma(c)` returns the chroma of a `Lab`, `Luv` or their variants color.
!!! note
    The other color types (e.g. `RGB`, `HSV` or `YCbCr`) are not supported
    because their definitions of *chroma* are not clear. Colorfulness, chroma
    and saturation are defined as distinct aspects by the CIE.
"""
chroma(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{Lab,DIN99,DIN99o,DIN99d}} = sqrt(c.a^2 + c.b^2)
chroma(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Luv} = sqrt(c.u^2 + c.v^2)
chroma(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{LCHab,LCHuv}} = c.c

"""
`hue(c)` returns the hue in degrees. This function does not guarantee that the
return value is in [0, 360].
"""
hue(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{HSV,HSL,HSI,LCHab,LCHuv}} = c.h
hue(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{Lab,DIN99,DIN99o,DIN99d}} =
    (h = atand(c.b, c.a); h < 0 ? h + 360 : h)
hue(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Luv} =
    (h = atand(c.v, c.u); h < 0 ? h + 360 : h)

# fallbacks for compN
_comp(::Val{N}, c::Colorant) where N = getfield(c, N)
_comp(::Val{N}, c::AlphaColor) where N = getfield(c, N + 1)
_comp(::Val{N}, c::AlphaColorN{N}) where N = alpha(c)

@noinline function _comp_error(c::ColorantN{N}, n::Int) where N
    io = IOBuffer()
    print(io, "attempt to access ", N, "-component color ")
    print(IOContext(io, :compact=>true), typeof(c), " with `comp", n, "`")
    throw(ArgumentError(String(take!(io))))
end

"""
`comp1(c)` extracts the first component you'd pass to the constructor
of the corresponding object.  For most color types without an alpha
channel, this is just the first field, but for types like `BGR` that
reverse the internal storage order this provides the value that you'd
use to reconstruct the color.

Specifically, for any `Color{T,3}`,

    c == typeof(c)(comp1(c), comp2(c), comp3(c))

returns true.
"""
comp1(c::Colorant) = _comp(Val(1), c)
comp1(c::Union{AbstractRGB, TransparentRGB}) = red(c)
comp1(c::Union{AbstractGray, TransparentGray}) = gray(c)

"`comp2(c)` extracts the second constructor argument (see `comp1`)."
comp2(c::Colorant) = _comp(Val(2), c)
comp2(c::ColorantN{1}) = _comp_error(c, 2)
comp2(c::Union{AbstractRGB, TransparentRGB}) = green(c)

"`comp3(c)` extracts the third constructor argument (see `comp1`)."
comp3(c::Colorant) = _comp(Val(3), c)
comp3(c::Union{ColorantN{1}, ColorantN{2}}) = _comp_error(c, 3)
comp3(c::Union{AbstractRGB, TransparentRGB}) = blue(c)

"`comp4(c)` extracts the fourth constructor argument (see `comp1`)."
comp4(c::Colorant) = _comp(Val(4), c)
comp4(c::Union{ColorantN{1}, ColorantN{2}, ColorantN{3}}) = _comp_error(c, 4)

"`comp5(c)` extracts the fifth constructor argument (see `comp1`)."
comp5(c::Colorant) = _comp(Val(5), c)
comp5(c::Union{ColorantN{1}, ColorantN{2}, ColorantN{3}, ColorantN{4}}) = _comp_error(c, 5)

"`color(c)` extracts the opaque color component from a Colorant (e.g., omits the alpha channel, if present)."
color(c::Color) = c
color(c::TransparentColorN{2,C}) where {C} = C(comp1(c))
color(c::TransparentColorN{3,C}) where {C} = C(comp1(c), comp2(c))
color(c::TransparentColorN{4,C}) where {C} = C(comp1(c), comp2(c), comp3(c))
color(c::TransparentColorN{5,C}) where {C} = C(comp1(c), comp2(c), comp3(c), comp4(c))

# Some of these traits exploit a nice trick: for subtypes, walk up the
# type hierarchy until we get to a stage where we can define the
# function in general

# recurse up the type hierarchy until you get to Colorant{T,N} for
# specific T,N.
to_top(::Type{Colorant{T,N}}) where {T,N} = Colorant{T,N}
@pure to_top(::Type{C}) where {C<:Colorant} = to_top(supertype(C))

to_top(c::Colorant) = to_top(typeof(c))

# eltype(RGB{Float32}) -> Float32
eltype(::Type{C}) where {C<:Colorant{T}} where {T} = T

eltype(c::Colorant) = eltype(typeof(c))

@pure function _parameter_upper_bound(t::UnionAll, idx)
    Base.rewrap_unionall((Base.unwrap_unionall(t)::DataType).parameters[idx], t)
end

# eltypes_supported(Colorant{T<:X}) -> X
function eltypes_supported(::Type{C}) where {C<:Colorant}
    Cb = base_colorant_type(C)
    isconcretetype(C) && C === Cb && return eltype(C)
    _eltypes_supported(Cb)
end
@pure function _eltypes_supported(::Type{C}) where {C<:Colorant}
    nameof(C) === :Colorant || return _eltypes_supported(supertype(C))
    _parameter_upper_bound(C, 1)
end

eltypes_supported(c::Colorant) = eltypes_supported(typeof(c))

"""
    issupported(C::Type, T::Type)::Bool

Returns `true` if `T` is a valid numeric eltype for `C<:Colorant`.
"""
issupported(::Type{C}, ::Type{T}) where {C<:Colorant,T} = T <: eltypes_supported(C)

"""
    CT = color_type(C::Type)
    CT = color_type(c::Colorant)

Return the type of the Color object, discarding any alpha channel.
For example,

    color_type(RGB)          === RGB
    color_type(RGB{Float32}) === RGB{Float32}
    color_type(ARGB{N0f8})   === RGB{N0f8}
    color_type(RGB(1,0,0))   === RGB{N0f8}

`color_type`, and related functions like [`base_color_type`](@ref),
[`base_colorant_type`](@ref), and [`ccolor`](@ref) are useful for manipulating types
when writing generic code.
"""
color_type(::Type{C}) where {C<:Color} = C
color_type(::Type{C}) where {C<:AlphaColor} = color_type(supertype(C))
color_type(::Type{C}) where {C<:ColorAlpha} = color_type(supertype(C))
function color_type(::Type{TC}) where {TC <: TransparentColor}
    _color_type(::Type{TC}) where {C, TC <: TransparentColor{C}} = C
    if TC isa UnionAll && nameof(TC) === :TransparentColor
        C1 = _parameter_upper_bound(TC, 1)
        color_type(C1) # simplify
    else
        _color_type(TC)
    end
end
color_type(::Type{Colorant{T,N}}) where {T,N} = Color{T,N}
color_type(::Type{ColorantN{N}}) where {N} = ColorN{N}
color_type(::Type{Colorant{T}}) where {T} = Color{T}
color_type(::Type{Colorant}) = Color
color_type(::Type{X}) where {X<:Number} = typeof(Gray(zero(X)))

color_type(c::Union{Colorant, Number}) = color_type(typeof(c))

"""
    Cbase = base_color_type(C::Type)
    Cbase = base_color_type(c::Colorant)

Return the Color type without alpha channel or a specified numeric element type.
`base_color_type` is similar to [`color_type`](@ref), except it "strips off" the
element type. It always returns a parametric type, which makes it convenient
for switching the numeric element type.

For example, compare

    base_color_type(RGB{N0f8})  === RGB
    base_color_type(RGBA{N0f8}) === RGB

vs

    color_type(RGB{N0f8})       === RGB{N0f8}

Switching element types can be done as follows:

    c64 = base_color_type(c){Float64}(color(c))

`base_color_type` discards any alpha channel. See [`base_colorant_type`](@ref)
if you want to preserve the alpha channel.
"""
base_color_type(::Type{C}) where {C<:Colorant} = base_colorant_type(color_type(C))
base_color_type(::Type{<:Number}) = Gray

base_color_type(x::Union{Colorant,Number}) = base_color_type(typeof(x))

## base_colorant_type implementation

base_colorant_type(::Type{C}) where {C<:Colorant} = isabstracttype(C) ? abstract_basetype(C) : basetype(C)
base_colorant_type(::Type{<:Number}) = Gray

@pure basetype(@nospecialize(C)) = Base.typename(C).wrapper

abstract_basetype(::Type{<:AbstractRGB}) = AbstractRGB
abstract_basetype(::Type{<:ColorN{N}}) where N = ColorN{N}
abstract_basetype(::Type{<:Color}) = Color
function abstract_basetype(::Type{TC}) where {TC<:TransparentColor}
    P = TC <: AlphaColor ? AlphaColor : TC <: ColorAlpha ? ColorAlpha : TransparentColor
    _abstract_transparent_basetype(TC, P, base_color_type(TC))
end
function _abstract_transparent_basetype(::Type{TC}, P, Cb) where {N, TC <: TransparentColorN{N}}
    Cb === Color && return P{C,T,N} where {T,C<:ColorN{N-1,T}}
    isabstracttype(Cb) || return P{Cb{T},T,N} where T
    return P{C,T,N} where {T,C<:Cb{T}}
end
function _abstract_transparent_basetype(::Type{TC}, P, Cb) where {TC <: TransparentColor}
    TC === P && return P
    isabstracttype(Cb) || return P{Cb{T},T} where T
    return P{C,T} where {T,C<:Cb{T}}
end
# This fallback dispatches to a separate function to control method ordering.
# Otherwise the generic ColorantN methods supersede the AlphaColor/ColorAlpha/Transparent methods
abstract_basetype(::Type{C}) where C = abstract_colorant_basetype(C)

abstract_colorant_basetype(::Type{<:ColorantN{N}}) where N = ColorantN{N}
abstract_colorant_basetype(::Type{<:Colorant})             = Colorant

"""
    Cbase = base_colorant_type(C::Type)
    Cbase = base_colorant_type(c::Colorant)

Return the Colorant type without specified numeric element type.

For example, compare

    base_colorant_type(ARGB{N0f8}) === ARGB

vs

    base_color_type(ARGB{N0f8})    === RGB
    color_type(ARGB{N0f8})         === RGB{N0f8}

When possible, `base_colorant_type` returns a parametric `UnionAll`, so
that `Cbase{T}` can be used to specify a colorant with element type `T`.
However, `base_colorant_type` is not guaranteed to return a parametric result,
for example

```jldoctest; setup = :(using ColorTypes)
julia> base_colorant_type(RGB24)
RGB24
```

In generic code, you can check whether `Cbase` is a `DataType` or `UnionAll`;
`Cbase{T}` will be valid only if it's a `UnionAll`.

# Usage example

A safe and easy way to switch the element type of a colorant value `c` to be `Float64` is

    Cbase = base_colorant_type(parametric_colorant(typeof(c)))
    c64 = Cbase{Float64}(c)

"""
base_colorant_type(c::Union{Colorant, Number}) = base_colorant_type(typeof(c))

"""
    Cp = parametric_colorant(C::Type)

Return a parametric colorant type.

# Examples

```jldoctest; setup = :(using ColorTypes, ColorTypes.FixedPointNumbers)
julia> parametric_colorant(RGB)
RGB

julia> parametric_colorant(RGB{Float32})
RGB{Float32}

julia> parametric_colorant(BGR)
BGR

julia> parametric_colorant(RGB24) == RGB{N0f8}
true
```
"""
parametric_colorant(::Type{C}) where C<:Colorant = C
parametric_colorant(::Type{RGB24}) = RGB{N0f8}
parametric_colorant(::Type{Gray24}) = Gray{N0f8}
parametric_colorant(::Type{ARGB32}) = ARGB{N0f8}
parametric_colorant(::Type{AGray32}) = AGray{N0f8}

"""
    CF = floattype(C::Type)

Promote storage data type of colorant type `C` to `AbstractFloat` while keeping the
`base_colorant_type` the same.

!!! info

    Non-parametric colorants will be promoted to corresponding parametric
    colorants. For example, `floattype(RGB24) == RGB{Float32}`.
"""
floattype(::Type{T}) where T <: Colorant =
    base_colorant_type(T){floattype(eltype(T))} # 1 parameter
# 0 parameter
floattype(::Type{RGB24})   = RGB{floattype(N0f8)}
floattype(::Type{Gray24})  = Gray{floattype(N0f8)}
floattype(::Type{ARGB32})  = ARGB{floattype(N0f8)}
floattype(::Type{AGray32}) = AGray{floattype(N0f8)}

@pure pureintersect(::Type{C1}, ::Type{C2}) where {C1,C2} = typeintersect(C1, C2)

"""
    Calpha, Cbase, T = colorsplit(C)

Split a color type `C` into three components: `T` is the numeric type,
`Cbase` is the `base_color_type(C)`, and `Calpha` is one of `Nothing`
(when `C` has no transparency), `AlphaColor`, or `ColorAlpha`.

If `C` is an abstract type, `Calpha` may be `TransparentColor` or `Any`.
"""
function colorsplit(::Type{C}) where C<:Colorant
    Calpha = C <: AlphaColor ? AlphaColor :
             C <: ColorAlpha ? ColorAlpha :
             C <: Color ? Nothing :
             C <: TransparentColor ? TransparentColor : Any
    Cbase = C <: Union{Color,TransparentColor} ? base_color_type(C) : Color
    return Calpha, Cbase, eltype(C)
end

"""
    C = ccolor(Cdest, Csrc)

`ccolor` ("concrete color") supports independent selection of colorspace and
numeric element type. `ccolor` chooses the numeric element type from
`Cdest` if available, but if not specified gets it from `Csrc`.

```jldoctest; setup = :(using ColorTypes, ColorTypes.FixedPointNumbers)
julia> ccolor(RGB, Gray{Float32})
RGB{Float32}

julia> ccolor(RGB, Gray{N0f8}) == RGB{N0f8}
true

julia> ccolor(RGB{Float32}, Gray{N0f8})
RGB{Float32}
```

Some colorspaces don't support `FixedPoint` numeric element types;
in such cases the `eltype_default` for `Cdest` is chosen:

```jldoctest; setup = :(using ColorTypes, ColorTypes.FixedPointNumbers)
julia> ccolor(HSV, RGB{N0f8})
HSV{Float32}
```

`Cdest` can be an abstract type, in which case `Csrc` must be in that
abstract color space.

```jldoctest; setup = :(using ColorTypes, ColorTypes.FixedPointNumbers)
julia> ccolor(Color{Float32}, RGB{N0f8})
RGB{Float32}

julia> ccolor(AbstractRGB{Float32}, BGR{N0f8})
BGR{Float32}

julia> ccolor(AbstractRGB{Float32}, HSV{Float32})
ERROR: in ccolor, empty intersection between AbstractRGB and HSV
[...]
```
`ccolor` will throw an error when no unambiguous concrete type can be determined.
In the previous case, `AbstractRGB{Float32}` encompases `RGB{Float32}` and `BGR{Float32}`.

`ccolor` is most useful in defining other methods.
For example, to allow users to write `convert(RGB, c)` without having to specify the
numeric element type of `RGB`, define

```
convert(::Type{C}, p::Colorant) where C<:Colorant = cnvt(ccolor(C,typeof(p)), p)
```

where `cnvt` is the function that performs explicit conversion.
"""
function ccolor(::Type{Cdest}, ::Type{Csrc}) where {Cdest<:Colorant, Csrc<:Union{Number,Colorant}}
    Cdestalpha, Cdestbase, Tdest = colorsplit(Cdest)
    if Csrc <: Number
        Cdestbase <: Union{AbstractGray, AbstractRGB} || throw(ColorTypeResolutionError(:ccolor, "no automatic conversion from", Csrc, Cdestbase))
        Csrcalpha, Csrcbase, Tsrc = Any, Color, Csrc
    else
        Csrcalpha, Csrcbase, Tsrc = colorsplit(Csrc)
    end
    # Step 1: pick the base color type
    if isabstracttype(Cdestbase)
        C = pureintersect(Cdestbase, Csrcbase)
        C === Union{} && throw(ColorTypeResolutionError(:ccolor, "empty intersection between", Cdestbase, Csrcbase))
        isabstracttype(C) && throw(ColorTypeResolutionError(:ccolor, "abstract intersection between", Cdestbase, Csrcbase))
    else
        C = Cdestbase
    end
    # Step 2: combine with the alpha representation
    if Cdestalpha === Nothing
    else
        Calpha = (Cdestalpha === TransparentColor || Cdestalpha === Any) ? Csrcalpha : Cdestalpha
        Cdestalpha === TransparentColor && Calpha === Nothing && throw(ColorTypeResolutionError(:ccolor, "ambiguous alpha storage between", Cdest, Csrc))
        if Calpha !== Nothing
            Calpha === TransparentColor && throw(ColorTypeResolutionError(:ccolor, "ambiguous alpha storage between", Cdest, Csrc))
            C = Calpha <: AlphaColor ? alphacolor(C) :
                Calpha <: ColorAlpha ? coloralpha(C) : error("unexpected Calpha ", Calpha)
        end
    end
    # Step 3: assign the eltype
    if !isa(C, UnionAll)
        eltype(C) <: Tdest && return C   # RGB24 etc.
        error("nonparametric type ", C, " has ambiguous destination ", Cdest)
    end
    isabstracttype(Tdest) || !issupported(C, Tdest) || return C{Tdest}
    T = pureintersect(Tdest, Tsrc)
    T === Union{} && throw(ColorTypeResolutionError(:ccolor, "empty intersection between", Tdest, Tsrc))
    T === Any && return C
    issupported(C, T) && return C{T}
    Tdef = eltype_default(C)
    T<:Integer && return C{Tdef}
    Tf = floattype(Tsrc)
    issupported(C, Tf) && return C{Tf}
    return C{Tdef}
end

### Equality
function ==(c1::AbstractRGB, c2::AbstractRGB)
    red(c1) == red(c2) && green(c1) == green(c2) && blue(c1) == blue(c2)
end
==(c1::HSV, c2::HSV) = c1.h == c2.h && c1.s == c2.s && c1.v == c2.v
==(c1::HSI, c2::HSI) = c1.h == c2.h && c1.s == c2.s && c1.i == c2.i
==(c1::HSL, c2::HSL) = c1.h == c2.h && c1.s == c2.s && c1.l == c2.l
==(c1::XYZ, c2::XYZ) = c1.x == c2.x && c1.y == c2.y && c1.z == c2.z
==(c1::xyY, c2::xyY) = c1.x == c2.x && c1.y == c2.y && c1.Y == c2.Y
==(c1::Lab, c2::Lab) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::Luv, c2::Luv) = c1.l == c2.l && c1.u == c2.u && c1.v == c2.v
==(c1::LCHab, c2::LCHab) = c1.l == c2.l && c1.c == c2.c && c1.h == c2.h
==(c1::LCHuv, c2::LCHuv) = c1.l == c2.l && c1.c == c2.c && c1.h == c2.h
==(c1::DIN99, c2::DIN99) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::DIN99d, c2::DIN99d) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::DIN99o, c2::DIN99o) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::LMS, c2::LMS) = c1.l == c2.l && c1.m == c2.m && c1.s == c2.s
==(c1::YIQ, c2::YIQ) = c1.y == c2.y && c1.i == c2.i && c1.q == c2.q
==(c1::YCbCr, c2::YCbCr) = c1.y == c2.y && c1.cb == c2.cb && c1.cr == c2.cr

==(x::AbstractGray, y::AbstractGray) = gray(x) == gray(y)
==(x::Number, y::AbstractGray) = x == gray(y)
==(x::AbstractGray, y::Number) = ==(y, x)

function ==(x::TransparentColor, y::TransparentColor)
    color(x) == color(y) && alpha(x) == alpha(y)
end

struct BoolTuple end
@inline BoolTuple(args::Bool...) = (args...,)

function _isapprox(a::Colorant, b::Colorant; kwargs...)
    componentapprox(x, y) = isapprox(x, y; kwargs...)
    all(ColorTypes._mapc(BoolTuple, componentapprox, a, b))
end
isapprox(a::C, b::C; kwargs...) where {C<:Colorant} =
    _isapprox(a, b; kwargs...)
isapprox(a::Colorant, b::Colorant; kwargs...) =
    _isapprox(base_colorant_type(a), base_colorant_type(b), a, b; kwargs...)
_isapprox(::Type{C}, ::Type{C}, a, b; kwargs...) where {C<:Colorant} =
    _isapprox(a, b; kwargs...)
_isapprox(::Type{<:AbstractRGB}, ::Type{<:AbstractRGB}, a, b; kwargs...) =
    _isapprox(RGB(a), RGB(b); kwargs...)
_isapprox(::Type{<:AbstractGray}, ::Type{<:AbstractGray}, a, b; kwargs...) =
    isapprox(gray(a), gray(b); kwargs...)
_isapprox(TA::Type, TB::Type, a, b; kwargs...) = false
isapprox(x::Number, y::AbstractGray; kwargs...) =
    isapprox(x, gray(y); kwargs...)
isapprox(x::AbstractGray, y::Number; kwargs...) =
    isapprox(y, x; kwargs...)

zero(::Type{C}) where {C<:Gray} = C(0)
oneunit(::Type{C}) where {C<:Gray} = C(1)

function Base.one(::Type{C}) where {C<:Gray}
    Base.depwarn("one($C) will soon switch to returning 1; you might need to switch to `oneunit`", :one)
    C(1)
end

Base.broadcastable(x::Colorant) = Ref(x)

Base.isless(a::AbstractGray, b::AbstractGray) =
    isless(gray(a), gray(b))


text:
function sortcolumns(a)

n = size(a,2);
b = sortslices([a; reshape(collect(1:n),1,n)], dims=2);
ind = Int.(b[end,:]);
b = b[1:end-1,:];

return b,ind

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using ParameterSchedulers
using Test

using InfiniteArrays: OneToInf

@testset "Decay" begin
    include("decay.jl")
end
@testset "Cyclic" begin
    include("cyclic.jl")
end
@testset "Complex" begin
    using ParameterSchedulers: Stateful, next!, reset!
    include("complex.jl")
end


text:
function pin_piece!(cur_piece::AbstractPiece)
  cur_player = cur_piece.owner

  cur_piece.i = cur_player.grid.rows
  cur_piece.j = cur_player.grid.cols / 2

  isodd(cur_piece.width) && ( cur_piece.i += 1 )
  ( cur_piece.width == 2 ) && ( cur_piece.i += 1 )

  return
end


text:
include("Tests/Tests.jl")
include("utilities.jl")

using Test, JuMP

include("square.jl")
include("invariant.jl")
include("controlled_invariant.jl")
include("switched_system.jl")


text:
using LearnBase
using Base.Test

# Test if types are exported properly
@test Cost <: Any
@test Penalty <: Cost
@test Loss <: Cost
@test UnsupervisedLoss <: Loss
@test SupervisedLoss <: Loss
@test MarginLoss <: SupervisedLoss
@test DistanceLoss <: SupervisedLoss

@test Minimizable <: Any
@test Transformation <: Any
@test StochasticTransformation <: Transformation

# Test if functions are exported properly
@test typeof(transform) <: Function
@test typeof(transform!) <: Function
@test typeof(getobs) <: Function
@test typeof(getobs!) <: Function
@test typeof(learn) <: Function
@test typeof(learn!) <: Function
@test typeof(update) <: Function
@test typeof(update!) <: Function

@test typeof(scaled) <: Function

@test typeof(value) <: Function
@test typeof(value!) <: Function
@test typeof(meanvalue) <: Function
@test typeof(sumvalue) <: Function
@test typeof(meanderiv) <: Function
@test typeof(sumderiv) <: Function
@test typeof(deriv) <: Function
@test typeof(deriv2) <: Function
@test typeof(deriv!) <: Function
@test typeof(value_deriv) <: Function
@test typeof(grad) <: Function
@test typeof(grad!) <: Function
@test typeof(addgrad!) <: Function
@test typeof(value_grad) <: Function
@test typeof(value_grad!) <: Function
@test typeof(prox) <: Function
@test typeof(prox!) <: Function

@test typeof(isminimizable) <: Function
@test typeof(isdifferentiable) <: Function
@test typeof(istwicedifferentiable) <: Function
@test typeof(isconvex) <: Function
@test typeof(isstrictlyconvex) <: Function
@test typeof(isstronglyconvex) <: Function
@test typeof(isnemitski) <: Function
@test typeof(isunivfishercons) <: Function
@test typeof(isfishercons) <: Function
@test typeof(islipschitzcont) <: Function
@test typeof(islocallylipschitzcont) <: Function
@test typeof(islipschitzcont_deriv) <: Function
@test typeof(isclipable) <: Function
@test typeof(ismarginbased) <: Function
@test typeof(isclasscalibrated) <: Function
@test typeof(isdistancebased) <: Function

@test typeof(issymmetric) <: Function

@test typeof(getobs) <: Function
@test typeof(getobs!) <: Function
@test typeof(datasubset) <: Function

@test typeof(targets) <: Function
@test typeof(LearnBase.gettarget) <: Function
@test typeof(LearnBase.gettargets) <: Function

@test DataView <: AbstractVector
@test DataView <: AbstractDataIterator
@test DataView{Int} <: AbstractVector{Int}
@test DataView{Int,Vector{Int}} <: AbstractDataIterator{Int,Vector{Int}}
@test AbstractObsView <: DataView
@test AbstractObsView <: AbstractObsIterator
@test AbstractObsView{Int,Vector{Int}} <: DataView{Int,Vector{Int}}
@test AbstractObsView{Int,Vector{Int}} <: AbstractObsIterator{Int,Vector{Int}}
@test AbstractBatchView <: DataView
@test AbstractBatchView <: AbstractBatchIterator
@test AbstractBatchView{Int,Vector{Int}} <: DataView{Int,Vector{Int}}
@test AbstractBatchView{Int,Vector{Int}} <: AbstractBatchIterator{Int,Vector{Int}}

@test DataIterator <: AbstractDataIterator
@test DataIterator{Int,Vector{Int}} <: AbstractDataIterator{Int,Vector{Int}}
@test ObsIterator <: DataIterator
@test ObsIterator <: AbstractObsIterator
@test ObsIterator{Int,Vector{Int}} <: DataIterator{Int,Vector{Int}}
@test ObsIterator{Int,Vector{Int}} <: AbstractObsIterator{Int,Vector{Int}}
@test BatchIterator <: DataIterator
@test BatchIterator <: AbstractBatchIterator
@test BatchIterator{Int,Vector{Int}} <: DataIterator{Int,Vector{Int}}
@test BatchIterator{Int,Vector{Int}} <: AbstractBatchIterator{Int,Vector{Int}}

@test typeof(fit) <: Function
@test typeof(fit!) <: Function
@test typeof(nobs) <: Function

@test ObsDim.Constant <: LearnBase.ObsDimension
@test ObsDim.First <: LearnBase.ObsDimension
@test ObsDim.Last <: LearnBase.ObsDimension
@test ObsDim.Undefined <: LearnBase.ObsDimension
@test typeof(ObsDim.Constant(2)) <: ObsDim.Constant{2}

# Test that LearnBase reuses StatsBase functions
using StatsBase
@test StatsBase.fit  == LearnBase.fit
@test StatsBase.fit! == LearnBase.fit!
@test StatsBase.nobs == LearnBase.nobs

# Test superset fallbacks
immutable MyStronglyConvexType end
LearnBase.isstronglyconvex(::MyStronglyConvexType) = true
LearnBase.islipschitzcont(::MyStronglyConvexType) = true
@test isstronglyconvex(MyStronglyConvexType())
@test isstrictlyconvex(MyStronglyConvexType())
@test isconvex(MyStronglyConvexType())
@test islipschitzcont(MyStronglyConvexType())
@test islocallylipschitzcont(MyStronglyConvexType())

# IntervalSet
let s = IntervalSet(-1,1)
    @test typeof(s) == IntervalSet{Int}
    @test typeof(s) <: AbstractSet
    for x in (-1,0,0.5,1,1.0)
        @test x in s
    end
    for x in (-1-1e-10, 1+1e-10, -Inf, Inf, 2, NaN)
        @test !(x in s)
    end
    for i=1:10
        x = rand(s)
        @test typeof(x) == Float64
        @test x in s
    end
    xs = rand(s, 10)
    @test typeof(xs) == Vector{Float64}
    for x in xs
        @test typeof(x) == Float64
        @test x in s
    end
    @test LearnBase.randtype(s) == Float64
    # @show s LearnBase.randtype(s)
end
let s = IntervalSet(-1,1.0)
    @test typeof(s) == IntervalSet{Float64}
    @test typeof(s) <: AbstractSet
    @test 1 in s
    # @show s LearnBase.randtype(s)
    @test length(s) == 1
end

# IntervalSet{Vector}
let s = IntervalSet([-1.,0.], [1.,1.])
    @test typeof(s) == IntervalSet{Vector{Float64}}
    @test typeof(s) <: AbstractSet
    @test LearnBase.randtype(s) == Vector{Float64}
    @test typeof(rand(s)) == Vector{Float64}
    @test rand(s) in s
    @test [-1, 0] in s
    @test !([-1.5,0] in s)
    @test !([0,2] in s)
    @test length(s) == 2
end

# DiscreteSet
let s = DiscreteSet([-1,1])
    @test typeof(s) == DiscreteSet{Vector{Int}}
    @test typeof(s) <: AbstractSet
    for x in (-1, 1, -1.0, 1.0)
        @test x in s
    end
    for x in (0, Inf, -Inf, NaN)
        @test !(x in s)
    end
    for i=1:10
        x = rand(s)
        @test typeof(x) == Int
        @test x in s
    end
    xs = rand(s, 10)
    @test typeof(xs) == Vector{Int}
    for x in xs
        @test typeof(x) == Int
        @test x in s
    end
    @test LearnBase.randtype(s) == Int
    @test length(s) == 2
    @test s[1] == -1
    # @show s LearnBase.randtype(s)
end
let s = DiscreteSet([-1,1.0])
    @test typeof(s) == DiscreteSet{Vector{Float64}}
    @test typeof(s) <: AbstractSet
    @test typeof(rand(s)) == Float64
    @test typeof(rand(s, 2)) == Vector{Float64}
    # @show s LearnBase.randtype(s)
end

# TupleSet
let s = TupleSet(IntervalSet(0,1), DiscreteSet([0,1]))
    @test typeof(s) == TupleSet{Tuple{IntervalSet{Int}, DiscreteSet{Vector{Int}}}}
    @test typeof(s) <: AbstractSet
    for x in ([0,0], [0.0,0.0], [0.5,1.0])
        @test x in s
    end
    for x in ([0,0.5], [-1,0])
        @test !(x in s)
    end
    @test typeof(rand(s)) == Vector{Float64}
    @test typeof(rand(s, 2)) == Vector{Vector{Float64}}
    @test typeof(rand(s, Tuple)) == Tuple{Float64,Int}
    @test typeof(rand(s, Tuple, 2)) == Vector{Tuple{Float64,Int}}
    @test LearnBase.randtype(s) == Vector{Float64}
    # @show s LearnBase.randtype(s)

    tot = 0
    for (i,x) in enumerate(s)
        @test x == s.sets[i]
        tot += length(x)
    end
    @test length(s) == tot
end

# arrays of sets
let s = [IntervalSet(0,1), DiscreteSet([0,1])]
    @test typeof(s) == Vector{AbstractSet}
    for x in ([0,0], [0.0,0.0], [0.5,1.0])
        @test x in s
    end
    for x in ([0,0.5], [-1,0])
        @test !(x in s)
    end
    @test typeof(rand(s)) == Vector{Float64}
    @test typeof(rand(s, 2)) == Vector{Vector{Float64}}
    # @show s LearnBase.randtype(s)
end

@testset "obsdim typetree and Constructor" begin
    @test_throws MethodError ObsDim.Constant(2.0)

    @test typeof(ObsDim.First()) <: LearnBase.ObsDimension
    @test typeof(ObsDim.First()) <: ObsDim.First
    @test typeof(ObsDim.First()) <: ObsDim.Constant{1}

    @test typeof(ObsDim.Last()) <: LearnBase.ObsDimension
    @test typeof(ObsDim.Last()) <: ObsDim.Last

    @test typeof(ObsDim.Constant(2)) <: LearnBase.ObsDimension
    @test typeof(ObsDim.Constant(2)) <: ObsDim.Constant{2}
end

@testset "obs_dim helper constructor" begin
    @test_throws ArgumentError convert(LearnBase.ObsDimension, "test")
    @test_throws ArgumentError convert(LearnBase.ObsDimension, 1.0)

    @test @inferred(convert(LearnBase.ObsDimension, ObsDim.First())) === ObsDim.First()
    @test @inferred(convert(LearnBase.ObsDimension, ObsDim.First())) === ObsDim.Constant(1)
    @test @inferred(convert(LearnBase.ObsDimension, ObsDim.Last()))  === ObsDim.Last()
    @test @inferred(convert(LearnBase.ObsDimension, ObsDim.Constant(2))) === ObsDim.Constant(2)

    @test_throws ErrorException @inferred convert(LearnBase.ObsDimension, 1)
    @test_throws ErrorException @inferred convert(LearnBase.ObsDimension, 6)
    @test convert(LearnBase.ObsDimension, 1) === ObsDim.First()
    @test convert(LearnBase.ObsDimension, 2) === ObsDim.Constant(2)
    @test convert(LearnBase.ObsDimension, 6) === ObsDim.Constant(6)
    @test_throws ErrorException @inferred convert(LearnBase.ObsDimension, :first)
    @test_throws ErrorException @inferred convert(LearnBase.ObsDimension, "first")
    @test convert(LearnBase.ObsDimension, (:first,:last))  === (ObsDim.First(),ObsDim.Last())
    @test convert(LearnBase.ObsDimension, :first)  === ObsDim.First()
    @test convert(LearnBase.ObsDimension, :begin)  === ObsDim.First()
    @test convert(LearnBase.ObsDimension, "first") === ObsDim.First()
    @test convert(LearnBase.ObsDimension, "BEGIN") === ObsDim.First()
    @test convert(LearnBase.ObsDimension, :end)   === ObsDim.Last()
    @test convert(LearnBase.ObsDimension, :last)  === ObsDim.Last()
    @test convert(LearnBase.ObsDimension, "End")  === ObsDim.Last()
    @test convert(LearnBase.ObsDimension, "LAST") === ObsDim.Last()
    @test convert(LearnBase.ObsDimension, :nothing) === ObsDim.Undefined()
    @test convert(LearnBase.ObsDimension, :none) === ObsDim.Undefined()
    @test convert(LearnBase.ObsDimension, :na) === ObsDim.Undefined()
    @test convert(LearnBase.ObsDimension, :null) === ObsDim.Undefined()
    @test convert(LearnBase.ObsDimension, :undefined) === ObsDim.Undefined()
    @test convert(LearnBase.ObsDimension, nothing) === ObsDim.Undefined()

    # TODO: remove test after deprecation phase
    @test LearnBase.obs_dim(:last) === ObsDim.Last()
end

immutable SomeType end
@testset "obsdim default values" begin
    @testset "Arrays, SubArrays, and Sparse Arrays" begin
        @test @inferred(LearnBase.default_obsdim(rand(10))) === ObsDim.Last()
        @test @inferred(LearnBase.default_obsdim(view(rand(10),:))) === ObsDim.Last()
        @test @inferred(LearnBase.default_obsdim(rand(10,5))) === ObsDim.Last()
        @test @inferred(LearnBase.default_obsdim(view(rand(10,5),:,:))) === ObsDim.Last()
        @test @inferred(LearnBase.default_obsdim(sprand(10,0.5))) === ObsDim.Last()
        @test @inferred(LearnBase.default_obsdim(sprand(10,5,0.5))) === ObsDim.Last()
    end

    @testset "Types with no specified default" begin
        @test @inferred(LearnBase.default_obsdim(SomeType())) === ObsDim.Undefined()
    end

    @testset "Tuples" begin
        @test @inferred(LearnBase.default_obsdim((SomeType(),SomeType()))) === (ObsDim.Undefined(), ObsDim.Undefined())
        @test @inferred(LearnBase.default_obsdim((SomeType(),rand(2,2)))) === (ObsDim.Undefined(), ObsDim.Last())
        @test @inferred(LearnBase.default_obsdim((rand(10),SomeType()))) === (ObsDim.Last(), ObsDim.Undefined())
        @test @inferred(LearnBase.default_obsdim((rand(10),rand(2,2)))) === (ObsDim.Last(), ObsDim.Last())
    end
end