text:
Base.@kwdef struct LearningMethodRegistry
    methods::Dict{Any,Any} = Dict()
end

"""
    registerlearningmethod!(registry, methodfn, blocktypes)

Register a learning method constructor `methodfn` as compatible with
`blocktypes` in `registry::LearningMethodRegistry`.

`blocks` should be the least specific set of types that a `methodfn`
can handle. `methodfn` needs to have a method that takes concrete block
instances as the only non-keyword argument, i.e. `methodfn(blocks; kwargs...)`.
"""
function registerlearningmethod!(reg::LearningMethodRegistry, f, blocktypes)
	reg.methods[f] = typify(blocktypes)
	return reg
end

"""
	findlearningmethods(blocktypes)
	findlearningmethods(registry, blocktypes)

Find learning methods compatible with block types `TBlocks` in
`registry::LearningMethodRegistry`.

#### Examples

```julia
julia> findlearningmethods((Image, Label))
[ImageClassificationSingle,]

julia> findlearningmethods((Image, Any))
[ImageClassificationSingle, ImageClassificationMulti, ImageSegmentation, ImageKeypointRegression, ...]
```
"""
function findlearningmethods(reg::LearningMethodRegistry, blocktypes=Any)
	return [methodfn for (methodfn, methodblocks) in reg.methods if blocktypesmatch(methodblocks, blocktypes)]
end

function blocktypesmatch(
        BSupported::Type{<:AbstractBlock},
        BWanted::Type{<:AbstractBlock})
    BWanted <: BSupported
end
function blocktypesmatch(B1::Type{<:Tuple}, B2::Type{<:Tuple})
    all(blocktypesmatch(b1, b2) for (b1, b2) in zip(B1.types, B2.types))
end

blocktypesmatch(BSupported::Type, _::Type{Any}) = true
blocktypesmatch(bsupported, bwanted) = blocktypesmatch(typify(bsupported), bwanted)
blocktypesmatch(BSupported::Type, bwanted) = blocktypesmatch(BSupported, typify(bwanted))

@testset "`blocktypesmatch`" begin
    @test blocktypesmatch(FastAI.Image, FastAI.Image{2})
    @test blocktypesmatch(FastAI.Image, FastAI.Image)
    @test !blocktypesmatch(FastAI.Image{2}, FastAI.Image)
    @test blocktypesmatch(Tuple{FastAI.Image}, Tuple{FastAI.Image})
    @test !blocktypesmatch(Tuple{FastAI.Image{2}}, Tuple{FastAI.Image})
    @test blocktypesmatch(Tuple{FastAI.Image{2}}, Any)
    @test blocktypesmatch(FastAI.Image{2}(), FastAI.Image{2})
    @test blocktypesmatch(FastAI.Image, FastAI.Image{2}())
    @test blocktypesmatch((FastAI.Image, FastAI.Label), (FastAI.Image{2}(), FastAI.Label(1:10)))
    @test blocktypesmatch((FastAI.Image, FastAI.Label), (FastAI.Image{2}(), Any))
end

const FASTAI_METHOD_REGISTRY = LearningMethodRegistry()


text:
function accelerated_gradient_line_search(fcn, x, step_direction, f, search_options)
    return 1.0
end

mutable struct NesterovGradientDescentMethod <: DescentMethod
    x
    f

    g
    v
    v_prev
    step_direction

    alpha # learning rate
    beta  # momentum decay
    max_v # maximum speed
end
function NesterovGradientDescentMethod(;alpha = 0.8, beta = 0.61, max_v = 10)
    empty_array = zeros(Float64, 0)
    empty_scalar = zero(Float64)
    NesterovGradientDescentMethod(
        empty_array, empty_scalar, empty_array,
        empty_array, empty_array, empty_array,
        alpha,
        beta,
        max_v)
end

function init!(m::NesterovGradientDescentMethod, x0)
    z        = zero(eltype(x0))
    za       = zeros(eltype(x0), size(x0))

    m.x      = x0
    m.f      = z

    m.g      = copy(za)
    m.v      = copy(za)
    m.v_prev = copy(za)
    m.step_direction = copy(za)

    # m.alpha   = 0.8
    # m.beta   = 0.61
end

function step_direction_fcn_nesterov!(m::NesterovGradientDescentMethod, fcn, gradient_fcn)
    # accelerated steepest gradient
    # nesterov momentum

    m.g                 = gradient_fcn( fcn, m.x + m.beta * m.v )
    v_                  = m.beta * m.v_prev - m.alpha * m.g
    max_v_norm          = min(m.max_v, norm(v_))
    normalize!(v_)
    m.v                 = max_v_norm * v_
    m.step_direction    = m.v
    m.x                 = m.x + m.v
end

function update_fcn_nesterov!(m::NesterovGradientDescentMethod, x1, f1)
    # m.x      = x1
    m.f      = f1
    m.v_prev = m.v
end

function nesterov_momentum_gradient_descent(
    fcn,
    x0;
    tol_options::ToleranceOptions       = ToleranceOptions(),
    search_options::LineSearchOptions   = LineSearchOptions(),
    gradient_fcn                        = numder.gradient_fw,
    line_search_fcn                     = accelerated_gradient_line_search,
    log_path                            = false,
    method_parameters                   = nothing)

    if isnothing(method_parameters)
        m = NesterovGradientDescentMethod();
    else
        m = method_parameters;
    end
    init_fcn            = init!
    step_direction_fcn  = step_direction_fcn_nesterov!
    update_fcn          = update_fcn_nesterov!

    stat = gradient_descent_template(
        m,
        init_fcn,
        step_direction_fcn,
        update_fcn,
        fcn,
        x0;
        tol_options     = tol_options,
        search_options  = search_options,
        gradient_fcn    = gradient_fcn,
        line_search_fcn = line_search_fcn,
        log_path        = log_path)

    return stat

end

#=
# something is wrong with the implementation

mutable struct HyperNesterovGradientDescentMethod <: DescentMethod
    x
    f

    g
    g_prev
    v
    v_prev
    step_direction

    alpha # learning rate
    beta  # momentum decay
    max_v # maximum speed

    mu    # learning rate of the learning rate

end
function HyperNesterovGradientDescentMethod(;alpha = 0.8, beta = 0.61, max_v = 10, mu = 0.7)
    empty_array = zeros(Float64, 0)
    empty_scalar = zero(Float64)
    HyperNesterovGradientDescentMethod(
        empty_array, empty_scalar, empty_array,
        empty_array, empty_array, empty_array, empty_array,
        alpha,
        beta,
        max_v,
        mu)
end

function hyper_nesterov_momentum_gradient_descent(
    fcn,
    x0;
    tol_options::ToleranceOptions       = ToleranceOptions(),
    search_options::LineSearchOptions   = LineSearchOptions(),
    gradient_fcn                        = numder.gradient_fw,
    line_search_fcn                     = accelerated_gradient_line_search,
    log_path                            = false,
    method_parameters                   = nothing)

    if isnothing(method_parameters)
        m = HyperNesterovGradientDescentMethod();
    else
        m = method_parameters;
    end
    init_fcn            = init!
    step_direction_fcn  = step_direction_fcn_nesterov!
    update_fcn          = update_fcn_nesterov!

    stat = gradient_descent_template(
        m,
        init_fcn,
        step_direction_fcn,
        update_fcn,
        fcn,
        x0;
        tol_options     = tol_options,
        search_options  = search_options,
        gradient_fcn    = gradient_fcn,
        line_search_fcn = line_search_fcn,
        log_path        = log_path)

    return stat

end

function init!(m::HyperNesterovGradientDescentMethod, x0)
    z        = zero(eltype(x0))
    za       = zeros(eltype(x0), size(x0))

    m.x      = x0
    m.f      = z

    m.g      = copy(za)
    m.g_prev      = copy(za)

    m.v      = copy(za)
    m.v_prev = copy(za)
    m.step_direction = copy(za)

    # m.alpha   = 0.8
    # m.beta   = 0.61
end

function step_direction_fcn_nesterov!(m::HyperNesterovGradientDescentMethod, fcn, gradient_fcn)
    # accelerated steepest gradient
    # hyper nesterov momentum

    # m.g                 = gradient_fcn( fcn, m.x + m.beta * m.v )
    m.g                 = gradient_fcn( fcn, m.x )

    tmp = m.g .* (m.g_prev + m.beta * m.v)
    stmp = sum(tmp)
    m.alpha = m.alpha + m.mu * stmp

    v_                  = m.beta * m.v_prev - m.alpha * m.g
    max_v_norm          = min(m.max_v, norm(v_))
    normalize!(v_)
    m.v                 = max_v_norm * v_

    m.step_direction    = m.v
    # m.x                 = m.x + m.v
    m.x                 = m.x - m.alpha * (m.g + m.beta * m.v)
end

function update_fcn_nesterov!(m::HyperNesterovGradientDescentMethod, x1, f1)
    m.f      = f1
    m.v_prev = m.v
    m.g_prev = m.g
end

=#

ext: jl

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Model: Anand, McFarland, Kihm, Wong (Optimization of Aerosol Penetration through Transport Lines, 1992)
# Based on: straightTubePenetrationAnand1992.m

dir = "./src"
push!(LOAD_PATH,dir)

module tubeStraight

import physicalConstants: k, g, p, MM

struct Parameters
 # CarrierFluidProfile
    T :: Float64        # Carrier fluid temperature (T) SI units: K
    Ro :: Float64       # Carrier fluid density (ro) SI units: kg m-3
    Q :: Float64        # Carrier fluid flow rate (Q) SI units: m3 s-1
    mu :: Float64       # Carrier fluid dynamic viscosity (mu) SI units: N s m-2
end

## Input
print("Parameter details \n")
print("Enter T, Ro, Q, mu: ")
params = split(readline(), ",")

T = params[1]
Ro = params[2]
Q = params[3]
mu = params[4]

# print(T)
# print(Ro)
# print(Q)
# print(mu)
# print("\n")

mutable struct Properties
 # Droplet
    ro_d :: Float64                         # Droplet density (ro_d) SI units: kg m-3
    mu_d :: Float64                         # Droplet dynamic viscosity (mu_d) SI units: N s m-2
    d_d :: Float64                          # Droplet diameter (dd) SI units: m
    md :: Float64                           # Droplet mass (md) SI units: kg

 # Tube
    d :: Float64                            # Tube inner diameter (d) SI units: m
    phi :: Float64                          # Inclination angle in radians (phi) SI units: rad
    L :: Float64                            # Distance from the tube inlet (L) SI units: m
end

## Input
print("Droplet Properties \n")
print("Enter ro_d, mu_d, d_d: ")
propsDroplet = split(readline(), ",")

ro_d = propsDroplet[1]
mu_d = propsDroplet[2]
d_d = propsDroplet[3]

# print(ro_d)
# print(mu_d)
# print(d_d)
# print("\n")

print("Tube Properties \n")
print("Enter d, phi, L: ")
propsTube = split(readline(), ",")

function meanFlow(Q,  d)
        U=4*Q/(pi*d^2)                       # Fluid flow mean velocity (U) SI units: m s-1
end

function reynold(U,  d, ro, mu)
        Re=(U*d*ro)/mu                       # Reynolds number (Re)
end

function fff(Re)
        f=0.316/(4*Re^0.25)                  # Fanning friction factor, Blasius equation for turbulent pipe flow (f) dimensionless
end

# The Blasius equation is not applicable for the following range of Reynolds number:
function blasius(Re)
    if (Re .< 3000 || Re .> 100000)
        extrapolation=extrapolation+1
end

function frictionV(U, f)
        Vf = U*sqrt(f/2)                      # Friction velocity (Vf) SI units: m s-1
end

function CunninghamSlipCorrection(mu, T, d_d)
         lambda = (mu/p)*sqrt((pi*k*T)/(2*MM))                   # Model for the mean free path of the carrier fluid (lambda)
         C = 1 + (lambda/dd)*(2.34+1.05*exp(-.39*d_d/lambda))    #Cunningham slip correction (C)
end

function relaxationTime(C, ro_d, d_d,  mu)
        tau = (C*ro_d*d_d^2)/(18*mu)          # Droplet relaxation time (tau)
end

function Rplus(d_d, Vf, ro, mu)
        Rplus = (d_d*Vf*ro)/(2*mu)            # Parameter Rplus (Rplus)
end

function paramS(tau, Vf, d_d)
        ParamS = 0.9*tau*Vf + d_d/2            # Parameter S (ParamS)
end

function Splus(ParamS, Vf, ro, mu)
        Splus = (ParamS*Vf*ro)/(mu)            # Parameter Splus (Splus)
end

# Dimensionless depositional velocity (Vs)
if (Splus >= 0 && Splus <= 10)
    Vs = 0.05*Splus
else
    Vs = 0.5 + 0.0125*(Splus - 10)
    extrapolation=extrapolation+1
end

# Dimensionless depositional velocity (Vr)
if (Rplus >= 0 && Rplus <= 10)
    Vr = 0.05*Rplus
else
    Vr = 0.5 + 0.0125*(Splus - 10)
    extrapolation=extrapolation+1
end

function diffVelTur(Vf, Vs, Vr)
        Vt = (Vf*(Vs+Vr))/(4)                  # Turbulent diffusional velocity (Vt)
end

function gSettlingVel(tau, g, phi)
        Vg = tau*g*cos(phi)                    # Gravitational settling velocity (Vg)
end

function depositionV(Vt, Vb, K)
        Vd = Vt + Vb*K                         # Sum of Vt and Vb is defined as Vd (Vd)
end

function critAngle(Vd, Vg)
        TETAc = asin(Vd/Vg)                    # Critical angle (TETAc)
end

function effectiveDepositionV(Vd, TETAc, Vg)
        Ve=Vd*TETAc/pi + Vd/2 + Vg*cos(TETAc)/pi        # Ve is the effective depositional velocity, given by the vector sum of velocities due to turbulent diffusion (Vt),
end

# If Vd>Vg, there is no need to use TETAc, since the integral that originated the Ve equation can be evaluated from 0 to 2pi. Also, arcsine(x) for x>1 is not defined in the real domain, therefore, TETAc would be imaginary.
if (Vd>Vg)
    Ve=Vd
end

function penetration(d, Ve, L, Q)
        P = -exp(-(pi*d*Ve*L)/(Q))             # P = penetration = C/C0 = ratio between concentration at distance x and initial concentration
end

    d = parse(Float64, propsTube[1])
phi = parse(Float64, propsTube[2])
L = parse(Float64, propsTube[3])

print(typeof(d))

# print(d)
# print(phi)
# print(L)
# print("\n")

# function dropletMass(ro_d, d_d)
#     md = ro_d*(4/3)*(pi*(d_d/2)^3)     # Droplet mass (md) SI units: kg
# end

# md = ro_d*(4/3)*(pi*(d_d/2)^3)     # Droplet mass (md) SI units: kg

# dropletMass(ro_d, d_d)

# print(md)

# function velocityBr(k, T, md)
#     Vb = sqrt((k*T)/(2*pi*md))          # Brownian diffusion velocity (Vb) SI units: m s-1
# end

# test = velocityBr(k, T, md)

# print(test)

function velocityTu(Q, d)
end
#Debug test
# println(k)

end

end


text:
function parse_add_node!(st::Structure,vals::Vector{String})
    id=vals[1]
    x,y,z=parse.(Float64,vals[2:4])
    add_node!(st,id,x,y,z)
end

function parse_add_quad!(st::Structure,vals::Vector{String})
    id=vals[1]
    n1,n2,n3,n4=vals[2:5]
    mat=vals[6]
    t=parse(Float64,vals[7])
    add_quad!(st,id,n1,n2,n3,n4,mat,t,elm_type="TMGQ")
end

function parse_add_tria!(st::Structure,vals::Vector{String})
    id=vals[1]
    n1,n2,n3=vals[2:4]
    mat=vals[5]
    t=parse(Float64,vals[6])
    add_tria!(st,id,n1,n2,n3,mat,t,elm_type="TMGT")
end

function parse_line!(st::Structure,lcset::LoadCaseSet,line::String)
    ss=split(line,",")
    key=lowercase(strip(ss[1]))
    vals=lowercase.(strip.(ss[2:end]))
    if key=="mat"
    elseif key=="sec"
    elseif key=="node"
        parse_add_node!(st,vals)
    elseif key=="beam"
    elseif key=="quad"
        parse_add_quad!(st,vals)
    elseif key=="tria"
        parse_add_tria!(st,vals)
    elseif key=="lc"
    elseif key=="beamload"
    elseif key=="quadload"
    end
end

function parse_block!(st,lcset,block)
    if length(block)==1
        parse_line!(st,lcset,block[1])
    end
end


text:
@testset "graphs" begin
    in_channel = 3
    out_channel = 5
    N = 6

    adjs = Dict(
        :simple => [0. 1. 1. 0. 0. 0.;
                    1. 0. 1. 0. 1. 0.;
                    1. 1. 0. 1. 0. 1.;
                    0. 0. 1. 0. 0. 0.;
                    0. 1. 0. 0. 0. 0.;
                    0. 0. 1. 0. 0. 0.],
        :weight => [0. 2. 2. 0. 0. 0.;
                    2. 0. 1. 0. 2. 0.;
                    2. 1. 0. 5. 0. 2.;
                    0. 0. 5. 0. 0. 0.;
                    0. 2. 0. 0. 0. 0.;
                    0. 0. 2. 0. 0. 0.],
    )

    degs = Dict(
        :simple => [2. 0. 0. 0. 0. 0.;
                    0. 3. 0. 0. 0. 0.;
                    0. 0. 4. 0. 0. 0.;
                    0. 0. 0. 1. 0. 0.;
                    0. 0. 0. 0. 1. 0.;
                    0. 0. 0. 0. 0. 1.],
        :weight => [4. 0. 0. 0. 0. 0.;
                    0. 5. 0. 0. 0. 0.;
                    0. 0. 10. 0. 0. 0.;
                    0. 0. 0. 5. 0. 0.;
                    0. 0. 0. 0. 2. 0.;
                    0. 0. 0. 0. 0. 2.]
    )

    laps = Dict(
        :simple => [2. -1. -1. 0. 0. 0.;
                    -1. 3. -1. 0. -1. 0.;
                    -1. -1. 4. -1. 0. -1.;
                    0. 0. -1. 1. 0. 0.;
                    0. -1. 0. 0. 1. 0.;
                    0. 0. -1. 0. 0. 1.],
        :weight => [4. -2. -2. 0. 0. 0.;
                    -2. 5. -1. 0. -2. 0.;
                    -2. -1. 10. -5. 0. -2.;
                    0. 0. -5. 5. 0. 0.;
                    0. -2. 0. 0. 2. 0.;
                    0. 0. -2. 0. 0. 2.],
    )

    norm_laps = Dict(
        :simple => [1. -1/sqrt(2*3) -1/sqrt(2*4) 0. 0. 0.;
                    -1/sqrt(2*3) 1. -1/sqrt(3*4) 0. -1/sqrt(3) 0.;
                    -1/sqrt(2*4) -1/sqrt(3*4) 1. -1/2 0. -1/2;
                    0. 0. -1/2 1. 0. 0.;
                    0. -1/sqrt(3) 0. 0. 1. 0.;
                    0. 0. -1/2 0. 0. 1.],
        :weight => [1. -2/sqrt(4*5) -2/sqrt(4*10) 0. 0. 0.;
                    -2/sqrt(4*5) 1. -1/sqrt(5*10) 0. -2/sqrt(2*5) 0.;
                    -2/sqrt(4*10) -1/sqrt(5*10) 1. -5/sqrt(5*10) 0. -2/sqrt(2*10);
                    0. 0. -5/sqrt(5*10) 1. 0. 0.;
                    0. -2/sqrt(2*5) 0. 0. 1. 0.;
                    0. 0. -2/sqrt(2*10) 0. 0. 1.]
    )

    ugs = Dict(
        :simple => SimpleGraph(N),
        :weight => SimpleWeightedGraph(N),
    )

    dgs = Dict(
        :simple => SimpleDiGraph(N),
        :weight => SimpleWeightedDiGraph(N),
    )

    add_edge!(ugs[:simple], 1, 2); add_edge!(ugs[:simple], 1, 3);
    add_edge!(ugs[:simple], 2, 3); add_edge!(ugs[:simple], 3, 4)
    add_edge!(ugs[:simple], 2, 5); add_edge!(ugs[:simple], 3, 6)

    add_edge!(dgs[:simple], 1, 3); add_edge!(dgs[:simple], 2, 3)
    add_edge!(dgs[:simple], 1, 6); add_edge!(dgs[:simple], 2, 5)
    add_edge!(dgs[:simple], 3, 4); add_edge!(dgs[:simple], 3, 5)

    add_edge!(ugs[:weight], 1, 2, 2); add_edge!(ugs[:weight], 1, 3, 2)
    add_edge!(ugs[:weight], 2, 3, 1); add_edge!(ugs[:weight], 3, 4, 5)
    add_edge!(ugs[:weight], 2, 5, 2); add_edge!(ugs[:weight], 3, 6, 2)

    add_edge!(dgs[:weight], 1, 3, 2); add_edge!(dgs[:weight], 2, 3, 2)
    add_edge!(dgs[:weight], 1, 6, 1); add_edge!(dgs[:weight], 2, 5, -2)
    add_edge!(dgs[:weight], 3, 4, -2); add_edge!(dgs[:weight], 3, 5, -1)

    @testset "simplegraphs" begin
        g = ugs[:simple]
        for T in [Int8, Int16, Int32, Int64, Int128, Float16, Float32, Float64]
            @test GraphLaplacians.degrees(g, T, dir=:out) == T.(diag(degs[:simple]))
            @test GraphLaplacians.degree_matrix(g, T, dir=:out) == T.(degs[:simple])
            @test GraphLaplacians.degree_matrix(g, T, dir=:out) == GraphLaplacians.degree_matrix(adjs[:simple], T, dir=:in)
            @test GraphLaplacians.degree_matrix(g, T, dir=:out) == GraphLaplacians.degree_matrix(adjs[:simple], T, dir=:both)
            @test Graphs.laplacian_matrix(g, T) == T.(laps[:simple])
        end
        for T in [Float16, Float32, Float64]
            @test GraphLaplacians.normalized_laplacian(g, T) ≈ T.(norm_laps[:simple])
        end
    end

    @testset "weightedgraphs" begin
        g = ugs[:weight]
        for T in [Int8, Int16, Int32, Int64, Int128, Float16, Float32, Float64]
            @test GraphLaplacians.degrees(g, T, dir=:out) == T.(diag(degs[:weight]))
            @test GraphLaplacians.degree_matrix(g, T, dir=:out) == T.(degs[:weight])
            @test GraphLaplacians.degree_matrix(g, T, dir=:out) == GraphLaplacians.degree_matrix(adjs[:weight], T, dir=:in)
            @test GraphLaplacians.degree_matrix(g, T, dir=:out) == GraphLaplacians.degree_matrix(adjs[:weight], T, dir=:both)
            @test Graphs.laplacian_matrix(g, T) == T.(laps[:weight])
        end
        for T in [Float16, Float32, Float64]
            @test GraphLaplacians.normalized_laplacian(g, T) ≈ T.(norm_laps[:weight])
        end
    end
end


text:
module SpectralMethodFFT

include("core.jl")
include("exceptions.jl")
include("supplimentals.jl")
include("specific_dimensions.jl")

export ConfigFFT, XFunc, KFunc

export x_Xgen, xy_Xgen, xy_Ygen, xyz_Xgen, xyz_Ygen, xyz_Zgen,
       k_Kgen, kl_Kgen, kl_Lgen, klm_Kgen, klm_Lgen, klm_Mgen

export pass_K!, highpass_K!, K_highpass_K, lopass_K!, K_lowpass_K,
       k_highpass_k, kl_highpass_kl, klm_highpass_klm,
       k_lowpass_k, kl_lowpass_kl, klm_lowpass_klm

export K_prod_32_K_K, K_dot_32_Ks_Ks, K3_cross_32_K3_K3,
       ∗, ⋅, ×

export K_prod_23_K_K, K_dot_23_Ks_Ks, K3_cross_23_K3_K3,
       ⊛, ⊙, ⊗

export K_X, k_x, kl_xy, klm_xyz, X_K, x_k, xy_kl, xyz_klm

export K_∂Xaxis_K, k_∂x_k, kl_∂x_kl, kl_∂y_kl,
       klm_∂x_klm, klm_∂y_klm, klm_∂z_klm,
       X_∂Xaxis_X, x_∂x_x, xy_∂x_xy, xy_∂y_xy,
       xyz_∂x_xyz, xyz_∂y_xyz, xyz_∂z_xyz

export K_laplacian_K, K_Δ_K,
       k_laplacian_k, kl_laplacian_kl, klm_laplacian_klm,
       k_Δ_k, kl_Δ_kl, klm_Δ_klm,
       X_laplacian_X, X_Δ_X,
       x_laplacian_x, xy_laplacian_xy, xyz_laplacian_xyz,
       x_Δ_x, xy_Δ_xy, xyz_Δ_xyz

export K_laplainv_K, K_Δ⁻¹_K,
       k_laplainv_k, kl_laplainv_kl, klm_laplainv_klm,
       k_Δ⁻¹_k, kl_Δ⁻¹_kl, klm_Δ⁻¹_klm,
       X_laplainv_X, X_Δ⁻¹_X,
       x_laplainv_x, xy_laplainv_xy, xyz_laplainv_xyz,
       x_Δ⁻¹_x, xy_Δ⁻¹_xy, xyz_Δ⁻¹_xyz

export kl2_grad_kl, kl_rot_kl2, kl_div_kl2,
       klm3_grad_klm, klm3_rot_klm3, klm_div_klm3

export integ_X, integ_x, integ_xy, integ_xyz, ∫,
       norm_X, norm_x, norm_xy, norm_xyz,
       l2inpr_X_X, l2inpr_x_x, l2inpr_xy_xy, l2inpr_xyz_xyz

export xlens_xranges

end

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
using DataFrames
using Dates
using Documenter
using FeatureTransforms

makedocs(;
    modules=[FeatureTransforms],
    authors="Invenia Technical Computing Corporation",
    repo="https://github.com/invenia/FeatureTransforms.jl/blob/{commit}{path}#L{line}",
    sitename="FeatureTransforms.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://invenia.github.io/FeatureTransforms.jl",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Guide to Transforms" => "transforms.md",
        "Transform Interface" => "transform_interface.md",
        "Examples" => "examples.md",
        "TestUtils" => "test_utils.md",
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
)

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


text:
# The glaciers for ESA:

println("Time to load BITEModel.jl")
@time using BITEModel
const BM=BITEModel
using Parameters
import PyPlot
using Base.Test

"""
    get_glaciers(region)

Return the glacier numbers to be run.  Sorted with
large glaciers first.

Return:
- glnrs1 -- with radar
- glnrs2 -- without radar
- glnrs_all -- both combined.
"""
function get_glaciers(region)

    glnrs, lon, lat, subregion, area  = BM.read_attrib_file(region)
    # purge some
    inds = 1:length(glnrs)
    # if region==7
    #     inds = 1:length(glnrs)
    #     # inds = (.!( (lon.>18.5) .& (lat.>79.7) ) .& # rm Westfonna
    #     #         .!( (lon.>20) .& (lat.>79.2) )) # rm Austfonna
    # elseif region==11
    #     # Alps only
    #     inds = subregion.==1  # only Alps: 3892 glaciers
    if region==14
        inds = subregion.==2  # only Karakoram: 13757 glaciers
    end
    # else
    #     inds = 1:length(glnrs)
    # end;
    glnrs_all = glnrs[inds]

    # now get radar vs non-radar
    glnrs1 = BM.get_glaciers_with_radar(region)
    # purge glnrs1 which are not in glnrs_all
    glnrs1 = [n for n in glnrs1 if (n in glnrs_all)]
    glnrs2 = [n for n in glnrs_all if !(n in glnrs1)]

    # sort with big glaciers first
    area = BM.read_attrib_file(region, glnrs1)[5]
    glnrs1 = glnrs1[reverse(sortperm(area))]

    area = BM.read_attrib_file(region, glnrs2)[5]
    glnrs2 = glnrs2[reverse(sortperm(area))]

    area = BM.read_attrib_file(region, glnrs_all)[5]
    glnrs_all = glnrs_all[reverse(sortperm(area))]

    return glnrs1, glnrs2, glnrs_all
end

function fit_it(nr, region,
                sigma_of_model, fit_vars, fit_target, runtyp, posterior_2D_or_1D::Val{I},
                expect=nothing, stdev=nothing;
                store_thetas=false,
                update_cache=true, prior_only_resamples=300,
                retsol=false,
                dir_output=error("Provide dir_output as kwarg"),
                save_geotiff=true,
                save_sol=false,
                flroot="",
                fit_sigma=false,
                debugging=true,
                plkws...
                ) where I
    (gid, gl,gb,pp,pm,pn,pl,th0d,
     theta0, logposterior, logposterior1d, logprior, logprior_debug,
     loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
     pmcmc_defaults, fit_target) = BM.init(nr, region, fit_vars, fit_target,
                                           sigma_of_model, runtyp;
                                           fit_sigma=fit_sigma,
                                           update_cache=update_cache,
                                           plkws...)

    # Update prior (if needed) to take values of expect and stdev into account.
    # This assumes that the previous posteriors had independent normal distributions.
    if expect!=nothing
        @assert length(expect)==length(theta0.th0)
        logprior_ = logprior
        # update priors to take expect and stdev into account
        logprior = theta -> logprior_(theta) +
            sum( - (theta[i]-expect[i])^2/(2*stdev[i]^2) for i=1:length(expect))
        # update the logposterior and logposterior1d to take the new logprior into accout:
        logposterior, logposterior1d = make_logposterior(loglikelihood, loglikelihood1d, logprior)
    end
    if I==2 && !(runtyp=="test" || runtyp=="testmid")
        merge!(pmcmc_defaults, Dict(:niter=>2*10^4,
                                    :nthin =>10,
                                    :nburnin=>5*10^3,
                                    :nchains=>50))
    end

    tosample = (logposterior, logposterior1d, logprior)[I]
    pmcmc = BM.MCMCNum(;pmcmc_defaults...)

        # some tests:
    if debugging
        fwdm_fn(theta0.th0);
        print("    Time to run forward model:")
        @time fwdsol = fwdm_fn(theta0.th0);
        print("    Posterior value:")
        @show logposterior(theta0.th0)[1]
        print("    Time to calculate posterior:")
        @time logposterior(theta0.th0)[1]
        print("    Starting MCMC:")
    end

    if logprior(theta0.th0)==-Inf
        logprior_debug(theta0.th0)
        error("Zero density of prior of $gl at $(theta0.th0)")
    end

    sol = BM.mcmc(tosample, theta0, pmcmc;
                  verbose=debugging,
                  use_progress_meter=progress_bar,
                  misc=Dict(:logposterior=>logposterior,
                            :logposterior1d=>logposterior1d,
                            :logprior=>logprior,
                            :logprior_debug=>logprior_debug,
                            :loglikelihood=>loglikelihood,
                            :loglikelihood1d=>loglikelihood1d,
                            :fwdm_fn=>fwdm_fn,
                            :fwdm1d_fn=>fwdm1d_fn,
                            :fit_target=>fit_target,
                            :runtyp=>runtyp,
                            :posterior_or_prior=>I,
                            :expect=>expect,
                            :stdev=>stdev,
                            # pp, pm, pn, etc are part of FWDSol
                            ))

    # if we sampled the logprior or logposterior1d, then run the forward model for
    # prior_only_resamples to get ice thickness samples
    if I==2
        tosample = logposterior
        sol.misc[:fit_target] = BM.FitTarget.prior
        sol = BM.retrace_samples(tosample, sol, pmcmc, nsamples=prior_only_resamples, misc=sol.misc)
    end

    # postproc

    debugging && println("MCMC of experiment $exp of $(gl.gid) done in $(signif(sol.solverstats[:wall_time]/60,4)) minutes")
    debugging && println("   sampling $(["log-posterior", "log-prior"][I])")
    #debugging && I==1 && println("   Rhat=$(mean(sol.solverstats[:Rhat])), $(maximum(sol.solverstats[:Rhat]))")
    debugging && I==1 && println("   eff_sample_size=$(sol.solverstats[:eff_sample_size]), $(maximum(sol.solverstats[:eff_sample_sizes]))")
    debugging && println("\n")

    if dir_output!="" && save_geotiff
        BM.make_geotiff(sol, dir_output*"/h-$flroot-$(BM.getrgi(gid)).tif")
    end
    if dir_output!="" && save_sol
        open(dir_output*"/sol-$flroot-$(BM.getrgi(gid)).jls", "w") do io
            serialize(io, sol)
        end
    end
    if retsol
        return BM.condense_MCMCSol(sol, store_thetas=store_thetas), sol
    else
        return BM.condense_MCMCSol(sol, store_thetas=store_thetas)
    end
end

var2filename(gid, sigma_h_model, sigma_iv_model) =
    dir_output*"$(BM.getrgi(gid))_sh=$(sigma_h_model)_siv=$(sigma_iv_model)"

import OnlineStats
function calculate_prior_update(out1::Vector, drop_sigma)
    out1_ = filter(o-> o isa Dict, out1)
    # if there is no out1 data, return nothing:
    length(out1_)==0 && return nothing, nothing, nothing, nothing

    # find out length without sigma
    if drop_sigma
        d = Dict(out1_[1][:theta0][:names]...)
        # remove all sigmas
        haskey(d, :sigma_iv_model) && delete!(d, :sigma_iv_model)
        haskey(d, :sigma_h_model) && delete!(d, :sigma_h_model)
        # and see what's the last index
        len = maximum(map(maximum, values(d)))
    else
        len = length(out1_[1][:thetas_expect])
        @assert all([length(o[:thetas_expect]) for o in out1_].==len) "Not all o[:thetas_expect] the same length"
    end

    # calculate expectation.
    expect_mean = mean([o[:thetas_expect][1:len] for o in out1_])
    expect_median = median(hcat([o[:thetas_expect][1:len] for o in out1_]...),2)[:] # might be more robust because of outliers
    expect = expect_median

    # use OnlineStats to get the total variance (and quantiles (todo))
    v = OnlineStats.Group([OnlineStats.Variance() for i=1:len])
    q = OnlineStats.Group([OnlineStats.Quantile() for i=1:len])
    for o in out1_
        ns = o[:solverstats][:nsamples]
        v1 = OnlineStats.Group([OnlineStats.Variance(o[:thetas_std][i]^2, o[:thetas_expect][i], OnlineStats.EqualWeight, ns)
                                for i=1:len])
        merge!(v, v1)
    end
    stdev = sqrt.(map(x->OnlineStats.value(x), v))

    return expect, expect_median, expect_mean, stdev
end


text:
tests = [
    "edgelist",
    "adjlist",
    "bellman_test",
    "inclist",
    "graph",
    "gmatrix",
    "bfs",
    "dfs",
    "conn_comp",
    "dijkstra",
    "a_star_spath",
    "mst",
    "floyd",
    "dot",
    "cliques",
    "random" ]

for t in tests
    tp = joinpath(Pkg.dir("Graphs"),"test","$(t).jl")
    println("running $(tp) ...")
    include(tp)
end


text:
"""
Write the wrapper to `destfile`.
"""
function Base.write(vw::VulkanWrapper, destfile, docfile)
    exprs = vcat(vw.handles, vw.structs, vw.funcs, vw.enums)
    ordered_exprs = sort_expressions(exprs)
    structs = filter(is_category(:struct), ordered_exprs)
    funcs = filter(is_category(:function), ordered_exprs)
    enums = filter(is_category(:enum), ordered_exprs)

    open(destfile, "w+") do io
        print_block(io, enums)
        print_block(io, structs)
        print_block(io, funcs)
        print_block(io, setdiff(ordered_exprs, vcat(structs, funcs, enums)))

        write_exports(io, exprs)
    end

    open(docfile, "w+") do io
        print_block(io, vw.docs)
    end
end

function sort_expressions(exprs)
    exprs_order = resolve_dependencies(name.(exprs), exprs)
    ordered_exprs = exprs[exprs_order]
    check_dependencies(ordered_exprs)

    ordered_exprs
end

is_category(cat) = x -> category(x) == cat

function print_block(io::IO, exs)
    print.(Ref(io), block.(exs))
    println(io)
end

function block(ex::Expr)
    str = @match category(ex) begin
        :doc => string('\"'^3, ex.args[3], '\"'^3, '\n', ex.args[4])
        :enum => string(ex)
        _ => string(prettify(ex))
    end
    str * spacing(ex)
end

spacing(ex::Expr) = spacing(category(ex))

spacing(cat::Symbol) = @match cat begin
    :struct => '\n'^2
    :function => '\n'^2
    :const => '\n'
    :enum => '\n'^2
    :doc => '\n'^2
end

function write_exports(io::IO, decls)
    println(io)

    ignored_symbols = vcat(:(Base.convert), :Base)

    decl_symbols = sort(filter(!in(ignored_symbols), unique(Iterators.flatten(exported_names.(decls)))))

    exports = :(export $([decl_symbols; getproperty.(spec_all_semantic_enums, :name)]...))

    println(io, string(exports))
end

function exported_names(decl)
    @match category(decl) begin
        :enum => begin
            p = deconstruct(decl)
            [p[:name]; name.(p[:values])]
        end
        _ => [name(decl)]
    end
end


text:
## context

## A context stores objects where named values are looked up.
mutable struct Context
    view ## of what? a Dict, Module, CompositeKind, DataFrame.parent.
    parent ## a context or nothing
    _cache::Dict
end
function Context(view, parent=nothing)
    d = Dict()
    d["."] = isa(view, AnIndex) ? view.value : view
    Context(view, parent, d)
end

## used by renderTokens
function ctx_push(ctx::Context, view)
    Context(view, ctx) ## add context as a parent
end

# we have some rules here
# * Each part of a dotted name should resolve only against its parent.
# * Any falsey value prior to the last part of the name should yield ''.
# * The first part of a dotted name should resolve as any other name.
function lookup_dotted(ctx::Context, dotted)
    for key in split(dotted, ".")
        nctx = lookup(Context(ctx.view), key)
        falsy(nctx) && return nothing
        ctx = Context(nctx, ctx)
    end
    ctx.view
end

## Lookup value by key in the context
function lookup(ctx::Context, key)

    if haskey(ctx._cache, key)
        value = ctx._cache[key]
    else
        context = ctx
        value = nothing
        while value == nothing && context != nothing
            ## does name have a .?
            if occursin(r"\.", key)

                value = lookup_dotted(context, key)
                value != nothing && break

                ## do something with "."
                ## we use .[ind] to refer to value in parent of given index;
                m = match(r"^\.\[(.*)\]$", key)
                m == nothing && break
#                m == nothing && error("Not implemented. Can use Composite Kinds in the view.")

                idx = m.captures[1]
                vals = context.parent.view

                # this has limited support for indices: "end", or a number, but no
                # arithmetic, such as `end-1`.
                if isa(vals, Vector) # supports getindex(v, i)?
                   if idx == "end"
                       value = AnIndex(-1, vals[end])
                   else
                       ind = Base.parse(Int, idx)
                       value = AnIndex(ind, vals[ind])
                   end
                    break
                end
            else
                ## strip leading, trailing whitespace in key
                value = lookup_in_view(context.view, stripWhitespace(key))
            end
            context = context.parent

        end

        ## cache
        ctx._cache[key] = value
    end

    return(value)
end

## Lookup value in an object by key
## This of course varies based on the view.
## we special case dataframes here, so that we don't have to assume package is loaded
function lookup_in_view(view, key)
    if is_dataframe(view)
        if occursin(r":", key)  key = key[2:end] end
        key = Symbol(key)
        out = nothing
        if haskey(view, key)
            out = view[1, key] ## first element only
        end
        return out
    else
        _lookup_in_view(view, key)
    end
end

function _lookup_in_view(view::AbstractDict, key)
    ## is it a symbol?
    if occursin(r"^:", key)
        key = Symbol(key[2:end])
    end

    get(view, key, nothing)

end

function _lookup_in_view(view::NamedTuple, key)
    ## is it a symbol?
    if occursin(r"^:", key)
        key = Symbol(key[2:end])
    end
     if haskey(view, key)
        getindex(view, key)
    else
        nothing
    end
end

function _lookup_in_view(view::Module, key)

    hasmatch = false
    re = Regex("^$key\$")
    for i in names(view, all=true)
        if occursin(re, string(i))
            hasmatch = true
            break
        end
    end

    out = nothing
    if hasmatch
        out = getfield(view, Symbol(key))  ## view.key
    end
    out

end

## Default is likely not great, but we use CompositeKind
function _lookup_in_view(view, key)
    nms = fieldnames(typeof(view))
    re = Regex(key)
    has_match = false
    for i in nms
        if occursin(Regex(key), string(i))
            has_match=true
            break
        end
    end

    out = nothing
    if has_match
        out = getfield(view, Symbol(key))  ## view.key
    end

    out
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Rotations from an Earth-Fixed, Earth-Centered (ECEF) reference frame to an
#   Earth-Fixed Inertial (ECI) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export rECEFtoECI

"""
    rECEFtoECI([T,] ECEF, ECI, JD_UTC::Number [, eop_data])

Compute the rotation from an Earth-Centered, Earth-Fixed (`ECEF`) reference
frame to an Earth-Centered Inertial (`ECI`) reference frame at the Julian Day
[UTC] `JD_UTC`. The rotation description that will be used is given by `T`,
which can be `DCM` or `Quaternion`. The ECEF frame is selected by the input
`ECEF` and the `ECI` frame is selected by the input `ECI`. The possible values
are listed below. The model used to compute the rotation is specified by the
selection of the origin and destination frames. Currently, there are two models
supported: IAU-76/FK5 and IAU-2006 with 2010 conventions (CIO approach only).

# Rotation description

The rotations that aligns the ECEF with ECI can be described by Direction Cosine
Matrices or Quaternions. This is selected by the parameter `T`. The possible
values are:

* `DCM`: The rotation will be described by a Direction Cosine Matrix.
* `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported yet.**

# ECEF Frame

The ECEF frame is selected by the parameter `ECEF`. The possible values are:

* `ITRF()`: ECEF will be selected as the International Terrestrial Reference
            Frame (ITRF).
* `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
           frame.
* `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
            System (TIRS).

# ECI Frame

The ECI frame is selected by the parameter `ECI`. The possible values are:

* `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
            reference frame.
* `TOD()`: ECI will be selected as the True of Date (TOD).
* `MOD()`: ECI will be selected as the Mean of Date (MOD).
* `J2000()`: ECI will be selected as the J2000 reference frame.
* `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
            (GCRF).
* `CIRS()`: ECEF will be selected as the Celestial Intermediate Reference System
            (CIRS).

# EOP Data

The conversion between the frames depends on EOP Data (see `get_iers_eop` and
`read_iers_eop`). If IAU-76/FK5 model is used, then the type of `eop_data` must
be `EOPData_IAU1980`. Otherwise, if IAU-2006/2010 model is used, then the type
of `eop_data` must be `EOPData_IAU2000A`. The following table shows the
requirements for EOP data given the selected frames.

|   Model       |  ECEF  |   ECI   |    EOP Data     |
|:--------------|:-------|:--------|:---------------=|
| IAU-76/FK5    | `ITRF` | `GCRF`  | EOP IAU1980     |
| IAU-76/FK5    | `ITRF` | `J2000` | EOP IAU1980     |
| IAU-76/FK5    | `ITRF` | `MOD`   | EOP IAU1980     |
| IAU-76/FK5    | `ITRF` | `TOD`   | EOP IAU1980     |
| IAU-76/FK5    | `ITRF` | `TEME`  | EOP IAU1980     |
| IAU-76/FK5    | `PEF`  | `GCRF`  | EOP IAU1980     |
| IAU-76/FK5    | `PEF`  | `J2000` | Not required¹   |
| IAU-76/FK5    | `PEF`  | `MOD`   | Not required¹   |
| IAU-76/FK5    | `PEF`  | `TOD`   | Not required¹   |
| IAU-76/FK5    | `PEF`  | `TEME`  | Not required¹   |
| IAU-2006/2010 | `ITRF` | `CIRS`  | EOP IAU2000A    |
| IAU-2006/2010 | `ITRF` | `GCRF`  | EOP IAU2000A    |
| IAU-2006/2010 | `TIRS` | `CIRS`  | Not required¹   |
| IAU-2006/2010 | `TIRS` | `GCRF`  | Not required¹ ² |

`¹`: In this case, the Julian Time UTC will be assumed equal to Julian Time UT1
to compute the Greenwich Mean Sidereal Time. This is an approximation, but
should be sufficiently accurate for some applications. Notice that, if EOP Data
is provided, the Julian Day UT1 will be accurately computed.

`²`: In this case, the terms that account for the free-core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the ECEF reference
frame into alignment with the ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(:IAU1980);

julia> rECEFtoECI(DCM, ITRF(), GCRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33492e-5
 -0.000797313   0.00106478   0.999999

julia> rECEFtoECI(ITRF(), GCRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33492e-5
 -0.000797313   0.00106478   0.999999

julia> rECEFtoECI(PEF(), J2000(), DatetoJD(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619271      0.785176    -0.00133066
 -0.785177     -0.619272     3.45854e-5
 -0.000796885   0.00106622   0.999999

julia> rECEFtoECI(PEF(), J2000(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267      0.78518     -0.00133066
 -0.78518      -0.619267     3.45854e-5
 -0.000796879   0.00106623   0.999999

julia> rECEFtoECI(Quaternion, ITRF(), GCRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.4363098936462618 - 0.0005909969666939257.i + 0.00030510511316206974.j + 0.8997962182293519.k

julia> eop_IAU2000A = get_iers_eop(:IAU2000A);

julia> rECEFtoECI(ITRF(), GCRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33502e-5
 -0.000797312   0.00106478   0.999999

julia> rECEFtoECI(TIRS(), GCRF(), DatetoJD(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619271      0.785176    -0.00133066
 -0.785177     -0.619272     3.45884e-5
 -0.000796885   0.00106623   0.999999

julia> rECEFtoECI(Quaternion, ITRF(), GCRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.4363098936309669 - 0.000590996988144556.i + 0.0003051056555230158.j + 0.8997962182365703.k
```
"""
@inline rECEFtoECI(T_ECEF::T_ECEFs,
                   T_ECI::T_ECIs,
                   JD_UTC::Number,
                   eop_data::EOPData_IAU1980) =
    rECEFtoECI(DCM, T_ECEF, T_ECI, JD_UTC, eop_data)

@inline rECEFtoECI(T_ECEF::T_ECEFs_IAU_2006,
                   T_ECI::T_ECIs_IAU_2006,
                   JD_UTC::Number,
                   eop_data::EOPData_IAU2000A) =
    rECEFtoECI(DCM, T_ECEF, T_ECI, JD_UTC, eop_data)

# Specializations for those cases that EOP Data is not needed.
@inline rECEFtoECI(T_ECEF::Val{:PEF},
                   T_ECI::Union{Val{:J2000}, Val{:MOD}, Val{:TOD}, Val{:TEME}},
                   JD_UTC::Number) =
    rECEFtoECI(DCM, T_ECEF, T_ECI, JD_UTC)

@inline rECEFtoECI(T_ECEF::Val{:TIRS}, T_ECI::Union{Val{:CIRS}, Val{:GCRF}},
                   JD_UTC::Number) =
    rECEFtoECI(DCM, T_ECEF, T_ECI, JD_UTC)

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                 ITRF => GCRF
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:GCRF}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Compute the rotation.
    return rITRFtoGCRF_fk5(T, JD_UT1, JD_TT, x_p, y_p, δΔϵ_1980, δΔψ_1980)
end

#                                ITRF => J2000
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:J2000}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad

    # Compute the rotation.
    return rITRFtoGCRF_fk5(T, JD_UT1, JD_TT, x_p, y_p, 0, 0)
end

#                                 ITRF => MOD
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:MOD}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_PEF_ITRF = rITRFtoPEF_fk5(T, x_p, y_p)
    r_MOD_PEF  = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

    return compose_rotation(r_PEF_ITRF, r_MOD_PEF)
end

#                                 ITRF => TOD
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:TOD}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_PEF_ITRF = rITRFtoPEF_fk5(T, x_p, y_p)
    r_TOD_PEF  = rPEFtoTOD_fk5(T, JD_UT1, JD_TT, δΔψ_1980)

    return compose_rotation(r_PEF_ITRF, r_TOD_PEF)
end

#                                 ITRF => TEME
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:TEME}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_PEF_ITRF = rITRFtoPEF_fk5(T, x_p, y_p)
    r_TEME_PEF = rPEFtoTEME(T, JD_UT1)

    return compose_rotation(r_PEF_ITRF, r_TEME_PEF)
end

#                                 PEF => GCRF
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:GCRF}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_MOD_PEF  = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
    r_GCRF_MOD = rMODtoGCRF_fk5(T, JD_TT)

    return compose_rotation(r_MOD_PEF, r_GCRF_MOD)
end

#                                PEF => J2000
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:J2000}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute the rotation.
    r_MOD_PEF  = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_GCRF_MOD = rMODtoGCRF_fk5(T, JD_TT)

    return compose_rotation(r_MOD_PEF, r_GCRF_MOD)
end

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:J2000}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute the rotation.
    r_MOD_PEF  = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_GCRF_MOD = rMODtoGCRF_fk5(T, JD_TT)

    return compose_rotation(r_MOD_PEF, r_GCRF_MOD)
end

#                                 PEF => MOD
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:MOD}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Compute the rotation.
    return rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:MOD}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute the rotation.
    return rPEFtoMOD_fk5(T, JD_UT1, JD_TT, 0, 0)
end

#                                 PEF => TOD
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:TOD}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    x_p      = eop_data.x(JD_UTC)*arcsec2rad
    y_p      = eop_data.y(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Compute the rotation.
    return rPEFtoTOD_fk5(T, JD_UT1, JD_TT, δΔψ_1980)
end

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:TOD}, JD_UTC::Number)

    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute the rotation.
    return rPEFtoTOD_fk5(T, JD_UT1, JD_TT, 0)
end

#                                 PEF => TEME
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:TEME}, JD_UTC::Number,
                    eop_data::EOPData_IAU1980)

    # Get the time in UT1.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)

    # Compute the rotation.
    return rPEFtoTEME(T, JD_UT1)
end

function rECEFtoECI(T::T_ROT, ::Val{:PEF}, ::Val{:TEME}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC

    # Compute the rotation.
    return rPEFtoTEME(T, JD_UT1)
end

################################################################################
#                                IAU-2006/2010
################################################################################

#                                 ITRF => CIRS
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:CIRS}, JD_UTC::Number,
                    eop_data::EOPData_IAU2000A)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC)*arcsec2rad
    y_p = eop_data.y(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_TIRS_ITRF = rITRFtoTIRS_iau2006(T, JD_TT, x_p, y_p)
    r_CIRS_TIRS = rTIRStoCIRS_iau2006(T, JD_UT1)

    return compose_rotation(r_TIRS_ITRF, r_CIRS_TIRS)
end

#                                 ITRF => GCRF
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:ITRF}, ::Val{:GCRF}, JD_UTC::Number,
                    eop_data::EOPData_IAU2000A)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC)*arcsec2rad
    y_p = eop_data.y(JD_UTC)*arcsec2rad
    dX  = eop_data.dX(JD_UTC)*arcsec2rad
    dY  = eop_data.dY(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_TIRS_ITRF = rITRFtoTIRS_iau2006(T, JD_TT, x_p, y_p)
    r_CIRS_TIRS = rTIRStoCIRS_iau2006(T, JD_UT1)
    r_GCRF_CIRS = rCIRStoGCRF_iau2006(T, JD_TT, dX, dY)

    return compose_rotation(r_TIRS_ITRF, r_CIRS_TIRS, r_GCRF_CIRS)
end

#                                 TIRS => CIRS
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:TIRS}, ::Val{:CIRS}, JD_UTC::Number,
                    eop_data::EOPData_IAU2000A)

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)

    # Compute the rotation.
    return rTIRStoCIRS_iau2006(T, JD_UT1)
end

function rECEFtoECI(T::T_ROT, ::Val{:TIRS}, ::Val{:CIRS}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC

    # Compute the rotation.
    return rTIRStoCIRS_iau2006(T, JD_UT1)
end

#                                 TIRS => GCRF
# ==============================================================================

function rECEFtoECI(T::T_ROT, ::Val{:TIRS}, ::Val{:GCRF}, JD_UTC::Number,
                    eop_data::EOPData_IAU2000A)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    dX  = eop_data.dX(JD_UTC)*arcsec2rad
    dY  = eop_data.dY(JD_UTC)*arcsec2rad

    # Compute the rotation.
    r_CIRS_TIRS = rTIRStoCIRS_iau2006(T, JD_UT1)
    r_GCRF_CIRS = rCIRStoGCRF_iau2006(T, JD_TT, dX, dY)

    return compose_rotation(r_CIRS_TIRS, r_GCRF_CIRS)
end

function rECEFtoECI(T::T_ROT, ::Val{:TIRS}, ::Val{:GCRF}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC

    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute the rotation.
    r_CIRS_TIRS = rTIRStoCIRS_iau2006(T, JD_UT1)
    r_GCRF_CIRS = rCIRStoGCRF_iau2006(T, JD_TT)

    return compose_rotation(r_CIRS_TIRS, r_GCRF_CIRS)
end


text:
import ContinuumArrays: apply, MulQuasiMatrix
import FastTransforms: chebyshevpoints
import OffsetArrays: OffsetArray

@testset "$(rpad("Chebyshev Basis Tests",80))" begin

    t = ChebyshevT(2)
    u = ChebyshevU(2)
    d = Derivative(axes(t,1))

    @test t isa Chebyshev
    @test u isa Chebyshev

    @test ChebyshevTDerivative <: ChebyshevDerivative
    @test ChebyshevUDerivative <: ChebyshevDerivative

    @test apply(*,d,t) isa ChebyshevTDerivative
    @test apply(*,d,t) isa ChebyshevDerivative
    @test d*t isa ChebyshevTDerivative
    @test d*t isa ChebyshevDerivative
    @test t' isa ChebyshevTDerivative
    @test t' isa ChebyshevDerivative

    @test apply(*,d,u) isa ChebyshevUDerivative
    @test apply(*,d,u) isa ChebyshevDerivative
    @test d*u isa ChebyshevUDerivative
    @test d*u isa ChebyshevDerivative
    @test u' isa ChebyshevUDerivative
    @test u' isa ChebyshevDerivative

    @test grid(t) == t.x
    @test basis(t) == t.b
    @test nodes(t) == t.x
    @test nbasis(t) == 2
    @test nnodes(t) == 2
    @test eachindex(t) == 0:1
    @test order(t)  == 2
    @test degree(t) == 1

    @test grid(u) == u.x
    @test basis(u) == u.b
    @test nodes(u) == u.x
    @test nbasis(u) == 2
    @test nnodes(u) == 2
    @test eachindex(u) == 0:1
    @test order(u)  == 2
    @test degree(u) == 1

    t0 = Chebyshev(Float64, 2, Val(1))
    t1 = ChebyshevT(Float64, 2)
    t2 = ChebyshevT(Float32, 2)
    t3 = ChebyshevT(Float64, 3)
    t4 = Chebyshev(2, Val(1))

    @test hash(t) == hash(t0)
    @test hash(t) == hash(t1)
    @test hash(t) != hash(t2)
    @test hash(t) != hash(t3)
    @test hash(t) == hash(t4)

    @test t == t0
    @test t == t1
    @test t ≈  t2
    @test t ≠  t3
    @test t == t4

    @test  isequal(t, t0)
    @test  isequal(t, t1)
    @test !isequal(t, t2)
    @test !isequal(t, t3)
    @test  isequal(t, t4)

    u0 = Chebyshev(Float64, 2, Val(2))
    u1 = ChebyshevU(Float64, 2)
    u2 = ChebyshevU(Integer, 2)
    u3 = ChebyshevU(Float64, 3)
    u4 = Chebyshev(2, Val(2))

    @test hash(u) == hash(u0)
    @test hash(u) == hash(u1)
    @test hash(u) == hash(u2)
    @test hash(u) != hash(u3)
    @test hash(u) == hash(u4)

    @test u == u0
    @test u == u1
    @test u ≈  u2
    @test u ≠  u3
    @test u == u4

    @test  isequal(u, u0)
    @test  isequal(u, u1)
    @test !isequal(u, u2)
    @test !isequal(u, u3)
    @test  isequal(u, u4)

    @test hash(t) != hash(u)

    y = rand(5)

    z1 = parent([ t[y[i], j] for i in eachindex(y), j in eachindex(t)])
    z2 =   hcat([ t[y[i], :] for i in eachindex(y)]...)'
    z3 =   hcat([ t[y,    j] for j in eachindex(t)]...)
    z4 = parent(  t[y,    :] )

    @test z1 == z2 == z3 == z4

    z1 = parent([ u[y[i], j] for i in eachindex(y), j in eachindex(u)] )
    z2 =   hcat([ u[y[i], :] for i in eachindex(y)]...)'
    z3 =   hcat([ u[y,    j] for j in eachindex(u)]...)
    z4 = parent(  u[y,    :] )

    @test z1 == z2 == z3 == z4

    z1 = parent([ (d*t)[y[i], j] for i in eachindex(y), j in eachindex(t)])
    z2 =   hcat([ (d*t)[y[i], :] for i in eachindex(y)]...)'
    z3 =   hcat([ (d*t)[y,    j] for j in eachindex(t)]...)
    z4 = parent(  (d*t)[y,    :] )

    @test z1 == z2 == z3 == z4

    z1 = parent([ (d*u)[y[i], j] for i in eachindex(y), j in eachindex(u)])
    z2 =   hcat([ (d*u)[y[i], :] for i in eachindex(y)]...)'
    z3 =   hcat([ (d*u)[y,    j] for j in eachindex(u)]...)
    z4 = parent(  (d*u)[y,    :] )

    @test z1 == z2 == z3 == z4

end


text:
struct DirectedGraph
    nodes::Set{Int}
    incoming_edges::Vector{Set{Int}}
    outgoing_edges::Vector{Set{Int}}
end

DirectedGraph() = DirectedGraph(Set{Int}(), Set{Int}[], Set{Int}[])

nodes(graph::DirectedGraph) = graph.nodes

incoming_edges(graph::DirectedGraph, node::Int) = graph.incoming_edges[node]

outgoing_edges(graph::DirectedGraph, node::Int) = graph.outgoing_edges[node]

Base.copy(graph::DirectedGraph) = DirectedGraph(copy())

function create_node!(graph)
    push!(graph.incoming_edges, Set{Int}())
    push!(graph.outgoing_edges, Set{Int}())

    node = length(graph.incoming_edges)
    push!(graph.nodes, node)

    return node
end

function add_edge!(graph::DirectedGraph, edge::Pair{Int,Int})
    from, to = edge

    push!(graph.incoming_edges[to], from)
    push!(graph.outgoing_edges[from], to)

    return graph
end

function remove_edge!(graph::DirectedGraph, edge::Pair{Int,Int})
    from, to = edge

    delete!(graph.incoming_edges[to], from)
    delete!(graph.outgoing_edges[from], to)

    return graph
end

function remove_node!(graph::DirectedGraph, node::Int)
    for e in incoming_edges(graph, node)
        remove_edge!(graph, e => node)
    end

    for e in outgoing_edges(graph, node)
        remove_edge!(graph, node => e)
    end

    delete!(graph.nodes, node)

    return graph
end

function remove_nodes!(graph::DirectedGraph, nodes)
    for node in nodes
        remove_node!(graph, node)
    end

    return graph
end

# returns ancestors including the nodes themselves
function ancestors(graph::DirectedGraph, nodes::Set{Int})
    ancestors = copy(nodes)

    queue = Int[ancestors...]

    while !isempty(queue)
        node = pop!(queue)

        for ancestor in incoming_edges(graph, node)
            if !(ancestor in ancestors)
                push!(ancestors, ancestor)
                push!(queue, ancestor)
            end
        end
    end

    return ancestors
end

function descendants(graph::DirectedGraph, nodes::Set{Int})
    descendants = Set{Int}()

    queue = Int[nodes...]

    while !isempty(queue)
        node = pop!(queue)

        for descendant in outgoing_edges(graph, node)
            if !(descendant in descendants)
                push!(descendants, descendant)
                push!(queue, descendant)
            end
        end
    end

    return descendants
end

# https://en.wikipedia.org/wiki/Topological_sorting#Kahn's_algorithm
function topological_sort!(graph)
    # TODO: this is valid thanks to the way we currently manipulate the graph
    return sort(collect(nodes(graph)))

#=
    L = Int[]

    S = filter(node -> isempty(incoming_edges(graph, node)), nodes(graph))

    while !isempty(S)
        n = pop!(S)
        push!(L, n)

        for m in collect(outgoing_edges(graph, n))
            remove_edge!(graph, n => m)

            if isempty(incoming_edges(graph, m))
                push!(S, m)
            end
        end
    end

    #if !isempty(edges(graph))
    #    error("graph has at least one cycle")
    #end

    return L
=#
end

topological_sort(graph) = topological_sort!(deepcopy(graph))


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

module Base64
import Base: read, write, close, eof, empty!
export Base64EncodePipe, Base64DecodePipe, base64encode, base64decode

# Base64EncodePipe is a pipe-like IO object, which converts into base64 data sent
# to a stream. (You must close the pipe to complete the encode, separate from
# closing the target stream).  We also have a function base64encode(f,
# args...) which works like sprint except that it produces
# base64-encoded data, along with base64encode(args...)  which is equivalent
# to base64encode(write, args...), to return base64 strings.
# A Base64DecodePipe object can be used to decode base64-encoded data read from a stream
# , while function base64decode is useful for decoding strings
#############################################################################

"""
    Base64EncodePipe(ostream)

Returns a new write-only I/O stream, which converts any bytes written to it into
base64-encoded ASCII bytes written to `ostream`.
Calling [`close`](@ref) on the `Base64EncodePipe` stream
is necessary to complete the encoding (but does not close `ostream`).

# Examples
```jldoctest
julia> io = IOBuffer();

julia> iob64_encode = Base64EncodePipe(io);

julia> write(iob64_encode, "Hello!")
6

julia> close(iob64_encode);

julia> str = String(take!(io))
"SGVsbG8h"

julia> String(base64decode(str))
"Hello!"
```
"""
mutable struct Base64EncodePipe <: IO
    io::IO
    # writing works in groups of 3, so we need to cache last two bytes written
    b0::UInt8
    b1::UInt8
    nb::UInt8 # number of bytes in cache: 0, 1, or 2

    function Base64EncodePipe(io::IO)
        b = new(io,0,0,0)
        finalizer(b, close)
        return b
    end
end

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

# Based on code by Stefan Karpinski from https://github.com/hackerschool/WebSockets.jl (distributed under the same MIT license as Julia)

const b64chars = ['A':'Z';'a':'z';'0':'9';'+';'/']

const base64_pad = UInt8('=')

function b64(x::UInt8, y::UInt8, z::UInt8)
    n = Int(x)<<16 | Int(y)<<8 | Int(z)
    b64chars[(n >> 18) + 1],
    b64chars[(n >> 12) & 0b111111 + 1],
    b64chars[(n >> 6) & 0b111111 + 1],
    b64chars[(n ) & 0b111111 + 1]
end

function b64(x::UInt8, y::UInt8)
    a, b, c = b64(x, y, 0x0)
    a, b, c, base64_pad
end

function b64(x::UInt8)
    a, b = b64(x, 0x0, 0x0)
    a, b, base64_pad, base64_pad
end

const sentinel = typemax(UInt8)
const revb64chars = fill(sentinel, 256)
# Fill revb64chars
for (val, ch) in enumerate(b64chars)
    revb64chars[UInt8(ch)] = UInt8(val - 1)
end

# Decode a block of at least 2 and at most 4 bytes, received in encvec
# Returns the first decoded byte and stores up to two more in cache
function b64decode!(encvec::Vector{UInt8}, cache::Vector{UInt8})
    if length(encvec) < 2
        throw(ArgumentError("incorrect base64 format, block must be at least 2 and at most 4 bytes"))
    end
    @inbounds u = revb64chars[encvec[1]]
    @inbounds v = revb64chars[encvec[2]]
    empty!(cache)
    res = (u << 2) | (v >> 4)
    if length(encvec) > 2
        @inbounds w = revb64chars[encvec[3]]
        push!(cache, (v << 4) | (w >> 2))
    end
    if length(encvec) > 3
        @inbounds z = revb64chars[encvec[4]]
        push!(cache, (w << 6) | z)
    end
    res
end

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

function unsafe_write(b::Base64EncodePipe, x::Ptr{UInt8}, n::UInt)
    s = 1 # starting index
    # finish any cached data to write:
    if b.nb == 1
        if n >= 2
            write(b.io, b64(b.b0, unsafe_load(x, 1), unsafe_load(x, 2))...)
            s = 3
        elseif n == 1
            b.b1 = unsafe_load(x, 1)
            b.nb = 2
            return
        else
            return
        end
    elseif b.nb == 2
        if n >= 1
            write(b.io, b64(b.b0, b.b1, unsafe_load(x, 1))...)
            s = 2
        else
            return
        end
    end
    # write all groups of three bytes:
    while s + 2 <= n
        write(b.io, b64(unsafe_load(x, s), unsafe_load(x, s + 1), unsafe_load(x, s + 2))...)
        s += 3
    end
    # cache any leftover bytes:
    if s + 1 == n
        b.b0 = unsafe_load(x, s)
        b.b1 = unsafe_load(x, s + 1)
        b.nb = 2
    elseif s == n
        b.b0 = unsafe_load(x, s)
        b.nb = 1
    else
        b.nb = 0
    end
    n
end

function write(b::Base64EncodePipe, x::UInt8)
    if b.nb == 0
        b.b0 = x
        b.nb = 1
    elseif b.nb == 1
        b.b1 = x
        b.nb = 2
    else
        write(b.io, b64(b.b0,b.b1,x)...)
        b.nb = 0
    end
    1
end

function close(b::Base64EncodePipe)
    if b.nb > 0
        # write leftover bytes + padding
        if b.nb == 1
            write(b.io, b64(b.b0)...)
        else # b.nb == 2
            write(b.io, b64(b.b0, b.b1)...)
        end
        b.nb = 0
    end
    nothing
end

# like sprint, but returns base64 string
"""
    base64encode(writefunc, args...)
    base64encode(args...)

Given a [`write`](@ref)-like function `writefunc`, which takes an I/O stream as its first argument,
`base64encode(writefunc, args...)` calls `writefunc` to write `args...` to a base64-encoded
string, and returns the string. `base64encode(args...)` is equivalent to `base64encode(write, args...)`:
it converts its arguments into bytes using the standard [`write`](@ref) functions and returns the
base64-encoded string.

See also [`base64decode`](@ref).
"""
function base64encode(f::Function, args...)
    s = IOBuffer()
    b = Base64EncodePipe(s)
    f(b, args...)
    close(b)
    String(take!(s))
end
base64encode(x...) = base64encode(write, x...)

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

"""
    Base64DecodePipe(istream)

Returns a new read-only I/O stream, which decodes base64-encoded data read from `istream`.

# Examples
```jldoctest
julia> io = IOBuffer();

julia> iob64_decode = Base64DecodePipe(io);

julia> write(io, "SGVsbG8h")
8

julia> seekstart(io);

julia> String(read(iob64_decode))
"Hello!"
```
"""
mutable struct Base64DecodePipe <: IO
    io::IO
    # reading works in blocks of 4 characters that are decoded into 3 bytes and 2 of them cached
    cache::Vector{UInt8}
    encvec::Vector{UInt8}

    function Base64DecodePipe(io::IO)
        b = new(io,[],[])
        finalizer(b, close)
        return b
    end
end

function read(b::Base64DecodePipe, t::Type{UInt8})
    if !isempty(b.cache)
        return shift!(b.cache)
    else
        empty!(b.encvec)
        while !eof(b.io) && length(b.encvec) < 4
            c::UInt8 = read(b.io, t)
            @inbounds if revb64chars[c] != sentinel
                push!(b.encvec, c)
            end
        end
        return b64decode!(b.encvec,b.cache)
    end
end

eof(b::Base64DecodePipe) = isempty(b.cache) && eof(b.io)
close(b::Base64DecodePipe) = nothing

# Decodes a base64-encoded string

"""
    base64decode(string)

Decodes the base64-encoded `string` and returns a `Vector{UInt8}` of the decoded bytes.

See also [`base64encode`](@ref)

# Examples
```jldoctest
julia> b = base64decode("SGVsbG8h")
6-element Array{UInt8,1}:
 0x48
 0x65
 0x6c
 0x6c
 0x6f
 0x21

julia> String(b)
"Hello!"
```
"""
function base64decode(s)
    b = IOBuffer(s)
    try
        return read(Base64DecodePipe(b))
    finally
        close(b)
    end
end

end # module

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
"""
A rule contains all the information required to run a rule in a cellular
simulation, given an initial array. Rules can be chained together sequentially.

The output of the rule for an Rule is allways written to the current cell in the grid.
"""
abstract type Rule end

show(io::IO, rule::R) where R <: Rule = begin
    indent = get(io, :indent, "")
    printstyled(io, indent, Base.nameof(typeof(rule)); color=:red)
    if nfields(rule) > 0
        printstyled(io, " :\n"; color=:red)
        for fn in fieldnames(R)
            if fieldtype(R, fn) <: Union{Number,Symbol,String}
                println(io, indent, "    ", fn, " = ", repr(getfield(rule, fn)))
            else
                # Avoid prining arrays etc. Just show the type.
                println(io, indent, "    ", fn, " = ", fieldtype(R, fn))
            end
        end
    end
end

"""
PartialRule is for rules that manually write to whichever cells of the grid
that they choose, instead of updating every cell with their output.

Updates to the destination array (`dest(data)`) must be performed manually, while
the source array can be accessed with `source(data)`.

The dest array is copied from the source prior to running the `applyrule!` method.
"""
abstract type PartialRule <: Rule end

"""
A Rule That only accesses a neighborhood, defined by its radius distance from the current cell.

For each cell a buffer will be populated containing the neighborhood cells, accessible with
`buffer(data)`. This allows memory optimisations and the use of BLAS routines on the neighborhood.
It also means that and no bounds checking is required.

`NeighborhoodRule` must read only from the state variable and the
neighborhood_buffer array, and never manually write to the `dest(data)` array.
Its return value is allways written to the central cell.

Custom Neighborhood rules must return their radius with a `radius()` method.
"""
abstract type NeighborhoodRule{R} <: Rule end

neighborhood(rule::NeighborhoodRule) = rule.neighborhood

"""
A Rule that only writes to its neighborhood, defined by its radius distance from the current point.
TODO: should this exist?

Custom PartialNeighborhood rules must return their radius with a `radius()` method.
"""
abstract type PartialNeighborhoodRule{R} <: PartialRule end

neighborhood(rule::PartialNeighborhoodRule) = rule.neighborhood

"""
A Rule that only writes and accesses a single cell: its return value is the new
value of the cell. This limitation can be useful for performance optimisations.

Accessing the `data.source` and `data.dest` arrays directly is not guaranteed to have
correct results, and should not be done.
"""
abstract type CellRule <: Rule end


text:
# This file is required for Julia packages.
module Blog

import TikzPictures

using AlgebraOfGraphics
using Base64: base64encode
using CairoMakie

# Don't export methods to avoid confusing readers of the code.
# export

const PKGDIR = string(pkgdir(Blog))::String

"""
    blue

Blue from ColorBrewer's "Blues".
"""
blue = "#3182bd"

"""
    causal_preamble

Preamble to be used for causal graphs.
Source: https://dkumor.com/posts/technical/2018/08/15/causal-tikz/.
"""
causal_preamble = raw"""
\usepackage{tikz}
\usetikzlibrary{shapes,decorations,arrows,calc,arrows.meta,fit,positioning}
\tikzset{
    -Latex,auto,node distance =2 cm and 2 cm,semithick,
    state/.style ={ellipse, draw, minimum width = 0.7 cm},
    point/.style = {circle, draw, fill, inner sep=0.04cm, node contents={}},
    box/.style ={rectangle, draw, fill, node contents={}},
    bidirected/.style={Latex-Latex,dashed},
    el/.style = {inner sep=2pt, align=left, sloped}
}
"""

"""
  print_graph(name::String, width::String, tex::String, out_path::String)

Print a causal graph with a name without extension `name`, width `width` (usually in `px`) and
the Tikz code `text`.
The `out_path` is expected to be the output of `@OUTPUT`; use `print_graph_partial` as a
convinience helper.
"""
function print_graph(name::String, width::String, text::String, out_path::String)
  tp = TikzPictures.TikzPicture(text, options="scale=1.0", preamble=causal_preamble)
  path = joinpath(out_path, "$name.svg")
  TikzPictures.save(TikzPictures.SVG(path), tp)
  println("""
  ~~~

  ~~~
  """)
end

"""
    print_graph_partial(out_path)

Partial function application on `print_graph` to fix `out_path`.
"""
function print_graph_partial(out_path)
    return function(name::String, width::String, tex::String)
        print_graph(name, width, tex, out_path)
    end
end

function aog_svg(dir::String, name::String, fg;
        literate=false)
    file = "$name.svg"
    path = joinpath(dir, file)
    AlgebraOfGraphics.save(joinpath(dir, file), fg)
    if !literate
        println("\\fig{./code/$file}")
    end
end

const OG_IMAGE_DIR = joinpath(PKGDIR, "__site", "assets", "og-image")

"""
    aog_og_image(fg, filename::String)

Write `fg` to "__site/assets/og-image/\$(filename).png".
Writing directly to "/__site/assets" and not "_assets" to avoid having to run two Franklin passes.
This is useful for social media images (`og:image`).
"""
function aog_og_image(fg, filename::String)
    dir = OG_IMAGE_DIR
    mkpath(dir)
    image_path = joinpath(dir, "$(filename).png")
    AlgebraOfGraphics.save(image_path, fg)
    return image_path
end

function makie_svg(dir::String, name::String, scene;
        literate=false)
    file = "$name.svg"
    path = joinpath(dir, file)
    Makie.FileIO.save(joinpath(dir, file), scene)
    if !literate
        println("\\fig{./code/$file}")
    end
end

function makie_og_image(scene, filename::String)
    dir = OG_IMAGE_DIR
    mkpath(dir)
    image_path = joinpath(dir, "$(filename).png")
    Makie.FileIO.save(image_path, scene)
    return image_path
end

end # module

max_forks_repo_name: rikhuijzer/franklin-blog


text:
module TikhonovEnsembleKalmanSampling

using Distributed
using Statistics
using LinearAlgebra
using Optim

include("sampling.jl")

export teks_update, whteks_update
export teks, whteks

end # module


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: March 2022
# Copyright: Georgia Institute of Technology, 2022

module ConditionalNFs4Imaging

using DrWatson
using Flux
using JLD2
using JSON
using HDF5
using ArgParse
using Random
using DataFrames
using LinearAlgebra
using Distributions
using Statistics
using JUDI
using JUDI.TimeModeling
using ProgressMeter
using PyPlot
using Seaborn
using InvertibleNetworks
using CUDA

import DrWatson: _wsave
import Random: rand
import Base.getindex
import Distributions: logpdf, gradlogpdf
import InvertibleNetworks:
    forward, inverse, ResidualBlock, CouplingLayerBasic, ActivationFunction, conv1x1_grad_v

# Utilities.
include("./utils/load_experiment.jl")
include("./utils/data_loader.jl")
include("./utils/savefig.jl")
include("./utils/logpdf.jl")
include("./utils/config.jl")

# Network.
include("./network/conditional_network_hint.jl")
include("./network/invertible_layer_basic.jl")
include("./network/invertible_layer_conv1x1.jl")
include("./network/put_param.jl")

# Objective functions.
include("./objectives/objectives.jl")

# Forward modeling operator.
include("./born_modeling/create_operator.jl")

end


text:
struct DiscreteMDP
    # TODO: Use sparse matrices?
    T::Array{Float64, 3} # T(s,a,s′)
    R::Array{Float64, 2} # R(s,a) = ∑_s' R(s,a,s')*T(s,a,s′)
    γ::Float64
end

n_states(mdp::DiscreteMDP) = size(mdp.T, 1)
n_actions(mdp::DiscreteMDP) = size(mdp.T, 2)
discount(mdp::DiscreteMDP) = mdp.γ
ordered_states(mdp::DiscreteMDP) = collect(1:n_states(mdp))
ordered_actions(mdp::DiscreteMDP) = collect(1:n_actions(mdp))
state_index(mdp::DiscreteMDP, s::Int) = s

function transition(mdp::DiscreteMDP, s::Int, a::Int)
    return Categorical(mdp.T[s,a,:])
end
function generate_s(mdp::DiscreteMDP, s::Int, a::Int)
    s′ = 1
    r = rand() - mdp.T[s,a,s′]
    while r > 0.0 && s′ < size(mdp.T, 3)
        s′ += 1
        r -= mdp.T[s,a,s′]
    end
    return s′
end
reward(mdp::DiscreteMDP, s::Int, a::Int) = mdp.R[s,a]

function MDP(mdp::DiscreteMDP; γ::Float64=discount(mdp))
    return MDP(
        γ,
        ordered_states(mdp),
        ordered_actions(mdp),
        (s,a, s′=nothing) -> begin
            S′ = transition(mdp, s, a)
            if s′ == nothing
                return S′
            end
            return pdf(S′, s′)
        end,
        (s,a) -> reward(mdp, s, a),
        (s, a)->begin
            s′ = rand(transition(mdp,s,a))
            r = reward(mdp, s, a)
            return (s′, r)
        end
    )
end


text:
export cross_validate, cv_by_iter, regularization_path, get_train_and_test, precision_at_k

# the loss function evaluates the objective minus the regularization
# it is the default error metric
loss_fn(args...; kwargs...) = objective(args...; include_regularization=false, kwargs...)

# to use with error_metric when we have domains in the namespace, call as:
# cross_validate(glrm, error_fn = error_metric(glrm,domains,glrm.X,glrm.Y))
function cross_validate(glrm::AbstractGLRM;
                        nfolds=5, params=Params(),
                        verbose=true, use_folds=nfolds,
                        error_fn=loss_fn,
                        init=nothing,
                        do_obs_check = false)
    if verbose println("flattening observations") end
#    obs = flattenarray(map(ijs->map(j->(ijs[1],j),ijs[2]),zip(1:length(glrm.observed_features),glrm.observed_features)))
    obs = flatten_observations(glrm.observed_features)
    if verbose println("computing CV folds") end
    folds = getfolds(obs, nfolds, size(glrm.A)..., do_check = do_obs_check)
    train_glrms = Array(typeof(glrm), nfolds)
    test_glrms = Array(typeof(glrm), nfolds)
    train_error = Array(Float64, nfolds)
    test_error = Array(Float64, nfolds)
    for ifold=1:use_folds
        if verbose println("\nforming train and test GLRM for fold $ifold") end
        train_observed_features, train_observed_examples, test_observed_features, test_observed_examples = folds[ifold]
	    ntrain = sum(map(length, train_observed_features))
    	ntest = sum(map(length, test_observed_features))
	    if verbose println("training model on $ntrain samples and testing on $ntest") end
        # form glrm on training dataset
        train_glrms[ifold] = copy_estimate(glrm)
        train_glrms[ifold].observed_examples = train_observed_examples
        train_glrms[ifold].observed_features = train_observed_features
        # form glrm on testing dataset
        test_glrms[ifold] = copy_estimate(glrm)
        test_glrms[ifold].observed_examples = test_observed_examples
        test_glrms[ifold].observed_features = test_observed_features
        # evaluate train and test error
        if verbose println("fitting train GLRM for fold $ifold") end
        if init != nothing
            init(train_glrms[ifold])
        end
        fit!(train_glrms[ifold], params, verbose=verbose)
        if verbose println("computing train and test error for fold $ifold:") end
        train_error[ifold] = error_fn(train_glrms[ifold],
            parameter_estimate(train_glrms[ifold])...) / ntrain
        if verbose println("\ttrain error: $(train_error[ifold])") end
        test_error[ifold] = error_fn(test_glrms[ifold],
            parameter_estimate(train_glrms[ifold])...) / ntest
        if verbose println("\ttest error:  $(test_error[ifold])") end
    end
    return train_error, test_error, train_glrms, test_glrms
end

@compat function getfolds(obs::Array{Tuple{Int,Int},1}, nfolds, m, n; ntrials = 5, do_check = true)
    # partition elements of obs into nfolds groups
    groups = Array(Int, size(obs))
    rand!(groups, 1:nfolds)  # fill an array with random 1 through N
    # create the training and testing observations for each fold
    folds = Array(Tuple, nfolds)
    for itrial = 1:ntrials
        enough_observations = 0
        for ifold=1:nfolds
            train = obs[filter(i->groups[i]!=ifold, 1:length(obs))] # all the obs that didn't get the ifold label
            train_observed_features, train_observed_examples = sort_observations(train,m,n)
            if !do_check ||
                (check_enough_observations(train_observed_features) && check_enough_observations(train_observed_examples))
                    enough_observations += 1
            else
                warn("Not enough data to cross validate; one of the cross validation folds has no observations in one row or column. Trying again...")
                break
            end
            test = obs[filter(i->groups[i]==ifold, 1:length(obs))] # all the obs that did
            test_observed_features, test_observed_examples = sort_observations(test,m,n,check_empty=false)
            folds[ifold] = (train_observed_features, train_observed_examples,
                        test_observed_features,  test_observed_examples)
        end
        if enough_observations == nfolds
            return folds
        end
    end
    error("Not enough data to cross validate automatically.")
end

function check_enough_observations(observed_examples_or_features)
    all(map(length, observed_examples_or_features) .> 0)
end

function get_train_and_test(obs, m, n, holdout_proportion=.1)

    # generate random uniform number for each observation
    groups = Array(Float64, size(obs))
    rand!(groups)

    # create the training and testing observations
    # observation is in test set if random number < holdout_proportion
    train = obs[filter(i->(groups[i]>=holdout_proportion), 1:length(obs))]
    train_observed_features, train_observed_examples = sort_observations(train,m,n)
    test = obs[filter(i->(groups[i]<holdout_proportion), 1:length(obs))]
    test_observed_features, test_observed_examples = sort_observations(test,m,n,check_empty=false)
    return (train_observed_features, train_observed_examples,
                test_observed_features,  test_observed_examples)
end

function flatten_observations(observed_features::ObsArray)
    obs = @compat Array(Tuple{Int,Int}, 0)
    for (i, features_in_example_i) in enumerate(observed_features)
        for j in features_in_example_i
            push!(obs, (i,j))
        end
    end
    return obs
end

function flatten(x, y)
    state = start(x)
    if state==false
        push!(y, x)
    else
        while !done(x, state)
          (item, state) = next(x, state)
          flatten(item, y)
        end
    end
    y
end
flatten{T}(x::Array{T})=flatten(x,Array(T, 0))

function flattenarray(x, y)
    if typeof(x)<:Array
        for xi in x
          flattenarray(xi, y)
        end
    else
        push!(y, x)
    end
    y
end
flattenarray{T}(x::Array{T})=flattenarray(x,Array(T, 0))

function cv_by_iter(glrm::AbstractGLRM, holdout_proportion=.1,
                    params=Params(100,max_iter=1,abs_tol=.01,min_stepsize=.01),
                    ch = ConvergenceHistory("cv_by_iter");
                    verbose=true)
    # obs = flattenarray(map(ijs->map(j->(ijs[1],j),ijs[2]),zip(1:length(glrm.observed_features),glrm.observed_features)))
    obs = flatten_observations(glrm.observed_features)

    train_observed_features, train_observed_examples, test_observed_features, test_observed_examples =
        get_train_and_test(obs, size(glrm.A)..., holdout_proportion)

    # form glrm on training dataset
    train_glrm = copy_estimate(glrm)
    train_glrm.observed_examples = train_observed_examples
    train_glrm.observed_features = train_observed_features

    # form glrm on testing dataset
    test_glrm = copy_estimate(glrm)
    test_glrm.observed_examples = test_observed_examples
    test_glrm.observed_features = test_observed_features

    ntrain = sum(map(length, train_glrm.observed_features))
    ntest = sum(map(length, test_glrm.observed_features))

    niters = params.maxiters
    params.maxiters = 1
    train_error = Array(Float64, niters)
    test_error = Array(Float64, niters)
    if verbose
        @printf("%12s%12s%12s\n", "train error", "test error", "time")
        t0 = time()
    end
    for iter=1:niters
        # evaluate train and test error
        fit!(train_glrm, params, ch=ch, verbose=false)
        train_error[iter] = ch.objective[end] # objective(train_glrm, parameter_estimate(train_glrm)..., include_regularization=false)/ntrain
        test_error[iter] = objective(test_glrm, parameter_estimate(train_glrm)..., include_regularization=false)/ntest
        if verbose
            @printf("%12.4e%12.4e%12.4e\n", train_error[iter], test_error[iter], time() - t0)
        end
    end
    return train_error, test_error
end

function regularization_path(glrm::AbstractGLRM; params=Params(), reg_params=logspace(2,-2,5),
                                         holdout_proportion=.1, verbose=true,
                                         ch::ConvergenceHistory=ConvergenceHistory("reg_path"))
    if verbose println("flattening observations") end
    # obs = flattenarray(map(ijs->map(j->(ijs[1],j),ijs[2]),zip(1:length(glrm.observed_features),glrm.observed_features)))
    obs = flatten_observations(glrm.observed_features)

    if verbose println("splitting train and test sets") end
    train_observed_features, train_observed_examples, test_observed_features, test_observed_examples =
        get_train_and_test(obs, size(glrm.A)..., holdout_proportion)

    if verbose println("forming train and test GLRMs") end
    # form glrm on training dataset
    train_glrm = copy_estimate(glrm)
    train_glrm.observed_examples = train_observed_examples
    train_glrm.observed_features = train_observed_features

    # form glrm on testing dataset
    test_glrm = copy_estimate(glrm)
    test_glrm.observed_examples = test_observed_examples
    test_glrm.observed_features = test_observed_features

    return regularization_path(train_glrm, test_glrm; params=params, reg_params=reg_params,
                                         verbose=verbose,
                                         ch=ch)
end

# For each value of the regularization parameter,
# compute the training error, ie, average error (sum over (i,j) in train_glrm.obs of L_j(A_ij, x_i y_j))
# and the test error, ie, average error (sum over (i,j) in test_glrm.obs of L_j(A_ij, x_i y_j))
function regularization_path(train_glrm::AbstractGLRM, test_glrm::AbstractGLRM; params=Params(), reg_params=logspace(2,-2,5),
                                         verbose=true,
                                         ch::ConvergenceHistory=ConvergenceHistory("reg_path"))
    train_error = Array(Float64, length(reg_params))
    test_error = Array(Float64, length(reg_params))
    ntrain = sum(map(length, train_glrm.observed_features))
    ntest = sum(map(length, test_glrm.observed_features))
    if verbose println("training model on $ntrain samples and testing on $ntest") end
    @show params
    train_time = Array(Float64, length(reg_params))
    for iparam=1:length(reg_params)
        reg_param = reg_params[iparam]
        # evaluate train and test error
        if verbose println("fitting train GLRM for reg_param $reg_param") end
        scale_regularizer!(train_glrm, reg_param)
        # no need to restart glrm X and Y even if they went to zero at the higher regularization
        # b/c fit! does that automatically
        fit!(train_glrm, params, ch=ch, verbose=verbose)
        train_time[iparam] = ch.times[end]
        if verbose println("computing mean train and test error for reg_param $reg_param:") end
        train_error[iparam] = objective(train_glrm, parameter_estimate(train_glrm)..., include_regularization=false) / ntrain
        if verbose println("\ttrain error: $(train_error[iparam])") end
        test_error[iparam] = objective(test_glrm, parameter_estimate(train_glrm)..., include_regularization=false) / ntest
        if verbose println("\ttest error:  $(test_error[iparam])") end
    end
    return train_error, test_error, train_time, reg_params
end

function precision_at_k(train_glrm::GLRM, test_observed_features; params=Params(), reg_params=logspace(2,-2,5),
                        holdout_proportion=.1, verbose=true,
                        ch::ConvergenceHistory=ConvergenceHistory("reg_path"), kprec=10)
    m,n = size(train_glrm.A)
    ntrain = sum(map(length, train_glrm.observed_features))
    ntest = sum(map(length, test_observed_features))
    train_observed_features = train_glrm.observed_features
    train_error = Array(Float64, length(reg_params))
    test_error = Array(Float64, length(reg_params))
    prec_at_k = Array(Float64, length(reg_params))
    @compat solution = Array(Tuple{Float64,Float64}, length(reg_params))
    train_time = Array(Float64, length(reg_params))
    test_glrm = GLRM(train_glrm.A, train_glrm.losses, train_glrm.rx, train_glrm.ry, train_glrm.k,
                     X=copy(train_glrm.X), Y=copy(train_glrm.Y),
                     observed_features = test_observed_features)
    for iparam=1:length(reg_params)
        reg_param = reg_params[iparam]
        # evaluate train error
        if verbose println("fitting train GLRM for reg_param $reg_param") end
        scale!(train_glrm.rx, reg_param)
        scale!(train_glrm.ry, reg_param)
        train_glrm.X, train_glrm.Y = randn(train_glrm.k,m), randn(train_glrm.k,n) # this bypasses the error checking in GLRM(). Risky.
        X, Y, ch = fit!(train_glrm, params, ch=ch, verbose=verbose)
        train_time[iparam] = ch.times[end]
        if verbose println("computing train error and precision at k for reg_param $reg_param:") end
        train_error[iparam] = objective(train_glrm, X, Y, include_regularization=false) / ntrain
        if verbose println("\ttrain error: $(train_error[iparam])") end
        test_error[iparam] = objective(test_glrm, X, Y, include_regularization=false) / ntrain
        if verbose println("\ttest error: $(test_error[iparam])") end
        # precision at k
        XY = X'*Y
        q = sort(XY[:],rev=true)[ntrain] # the ntest+ntrain largest value in the model XY
        true_pos = 0; false_pos = 0
        kfound = 0
        for i=1:m
            if kfound >= kprec
                break
            end
            for j=1:n
                if kfound >= kprec
                    break
                end
                if XY[i,j] >= q
                    # i predict 1 and (i,j) was in my test set and i observed 1
                    if j in test_observed_features[i]
                        true_pos += 1
                        kfound += 1
                    # i predict 1 and i did not observe a 1 (in either my test *or* train set)
                    elseif !(j in train_observed_features[i])
                        false_pos += 1
                        kfound += 1
                    end
                end
            end
        end
        prec_at_k[iparam] = true_pos / (true_pos + false_pos)
        if verbose println("\prec_at_k:  $(prec_at_k[iparam])") end
        solution[iparam] = (sum(X)+sum(Y), sum(abs(X))+sum(abs(Y)))
        if verbose println("\tsum of solution, one norm of solution:  $(solution[iparam])") end
    end
    return train_error, test_error, prec_at_k, train_time, reg_params, solution
end


text:
module Geometry

using ..RecursiveApply
import LinearAlgebra

import StaticArrays: SVector

export ⊗

include("coordinates.jl")
include("axistensors.jl")
include("localgeometry.jl")
include("conversions.jl")
include("globalgeometry.jl")

end # module


text:
###
### COIN-OR Clp API Wrapper
###

__precompile__()

module Clp

include("ClpCInterface.jl")
include("ClpSolverInterface.jl")

using Clp.ClpMathProgSolverInterface
export ClpSolver

end # module


text:
export Add

struct Add
    blocks
    activation
    function Add(blocks...; activation=nothing)
        new(blocks, activation)
    end
end

function (r::Add)(x)
    y = r.blocks[1](x)
    for i in 2:length(r.blocks)
        y = y .+ r.blocks[i](x)
    end
    r.activation === nothing ? y : r.activation(y)
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames

@testset "prob circuit structure learn tests" begin
    function test_vtree_variable_partition(v, partition)
        if v isa PlainVtreeLeafNode
            @test variable(v) == partition
        else
            test_vtree_variable_partition(v.left, partition[1])
            test_vtree_variable_partition(v.right, partition[2])
        end
    end
    data = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
                1 1 1 1 1 1 1 1 1 1;
                0 0 0 0 0 0 0 0 0 0;
                0 1 1 0 1 0 0 1 0 1]))
    p1 = [[[[1,7],[9,10]],[4,6]], [[2,8],[3,5]]]
    p2 = [[[6,4],[2,[10,8]]],[[5,3],[1,[9,7]]]]
    p3 = [1, [[[3,[[2,[8,10]],5]],[4,6]],[7, 9]]]
    @test_throws ErrorException learn_vtree(data; alg=nothing)
    for (alg, vars) in zip([:bottomup, :topdown, :clt], [p1, p2, p3])
        v = learn_vtree(data; alg=alg)
        test_vtree_variable_partition(v, vars)
        @test num_variables(v) == 10
        @test num_edges(v) == 18
        @test num_nodes(v) == 19
    end
end


text:
using Pkg
Pkg.instantiate()
Pkg.develop(PackageSpec(path=pwd()))
Pkg.develop(PackageSpec(url="https://github.com/sunoru/Random123.jl"))
Pkg.build("Random123")
Pkg.build("RandomNumbers")


text:
###
### COIN-OR Clp API Wrapper
###

VERSION < v"0.7.0-beta2.199" && __precompile__()
module Clp

using Compat
# This 'using' is required to suppress a warning about Clp not having Libdl in its
# dependencies (Libdl is used by BinaryProvider), e.g.: bicycle1885/CodecZlib.jl#26.
using Libdl

include("ClpCInterface.jl")
include("ClpSolverInterface.jl")
include("MOIWrapper.jl")

using Clp.ClpMathProgSolverInterface
export ClpSolver

end # module