text:
###################
## general GARCH ##
###################

## type GARCH{T} <: AbstrUnivarModel
##     mu::Float64
##     kappa::Float64
##     alphas::Array{Float64, 1}
##     betas::Array{Float64, 1}
##     distr::T
## end

## type GARCH <: AbstrUnivarModel
## - how to deal with different lag lengths?!
## end

#############
## display ##
#############

## normal innovations
##-------------------

import Base.Multimedia.display
function display(gMod::GARCH_1_1{Normal})
    nDigits = 3
    typNam = string(typeof(gMod))
    println("\n$typNam :: Univariate GARCH(1,1), εₜ~N(0,1)")
    println("    μ: $(round(gMod.μ, nDigits))")
    println("    ĸ: $(round(gMod.κ, nDigits))")
    println("    α: $(round(gMod.α, nDigits))")
    println("    β: $(round(gMod.β, nDigits))")
end

## t innovations
##--------------

function display(gMod::GARCH_1_1{TDist})
    nDigits = 3
    typNam = string(typeof(gMod))
    nuVal = round(dof(gMod.distr), nDigits)
    println("\n$typNam :: Univariate GARCH(1,1), εₜ~t($nuVal)")
    println("    μ: $(round(gMod.μ, nDigits))")
    println("    ĸ: $(round(gMod.κ, nDigits))")
    println("    α: $(round(gMod.α, nDigits))")
    println("    β: $(round(gMod.β, nDigits))")
end

## GARCH_1_1_Fit
##--------------

function display(gModFit::GARCH_1_1_Fit)
    nDigits = 3
    typNam = string(typeof(gModFit))
    println("\n$typNam with model:")
    display(gModFit.model)
    println("    data: $(typeof(gModFit.data))")
    println("  sigmas: $(typeof(gModFit.sigmas))")
    println("    nllh: $(gModFit.nllh)")
end

#################
## description ##
#################

## normal innovations
##-------------------

function description(gMod::GARCH_1_1{Normal})
    nDigits = 3
    println("\nUnivariate GARCH(1,1) model, εₜ~N(0,1): ")
    println("    Xₜ = $(round(gMod.μ, nDigits)) + σₜ*εₜ, εₜ~N(0,1)")
    println("    σₜ² = $(round(gMod.κ, nDigits)) + $(round(gMod.α, nDigits))*(Xₜ₋₁ - μ)² + $(round(gMod.β, nDigits))*σₜ₋₁²")
end

## t innovations
##--------------

function description(gMod::GARCH_1_1{TDist})
    nDigits = 3
    nuVal = round(dof(gMod.distr), nDigits)
    println("\nUnivariate GARCH(1,1), εₜ~t($nuVal): ")
    println("    Xₜ = $(round(gMod.μ, nDigits)) + σₜ*εₜ, εₜ~N(0,1)")
    println("    σₜ² = $(round(gMod.κ, nDigits)) + $(round(gMod.α, nDigits))*(Xₜ₋₁ - μ)² + $(round(gMod.β, nDigits))*σₜ₋₁²")
end

##############
## simulate ##
##############

function simulate(gMod::GARCH_1_1, nObs::Int,
                  initSigma::Float64 = NaN)

    simVals = Array(Float64, nObs)
    innovs = simulateInnovations(gMod.distr, nObs)

    if isnan(initSigma)
        ## take long term variance as starting value
        initSigma = sqrt(gMod.κ/(1-gMod.α-gMod.β))
    end

    currSigma = initSigma
    simVals[1] = innovs[1]*initSigma
    for ii=2:nObs
        currSigma = sqrt(gMod.κ + gMod.α*simVals[ii-1]^2 +
                          gMod.β*currSigma^2)
        simVals[ii] = currSigma*innovs[ii]
    end
    return simVals
end

function simulateInnovations(d::Normal, nObs::Int)
    return rand(d, nObs)
end

function simulateInnovations(d::TDist, nObs::Int)
    nu = dof(d)
    tVariance = nu/(nu-2)
    return rand(d, nObs)/sqrt(tVariance)
end

###############
## getSigmas ##
###############

@doc doc"""
Retrieve hidden sigma series from given dataset with given parameters
and given initial sigma.
"""->
function getSigmas(params::Array{Float64, 1}, data::Array{Float64, 1})
    μ, κ, α, β, σ0 = params

    nObs = length(data)

    ## normalize data
    centered = data - μ

    ## get sigma series
    sigmas = Array(Float64, nObs)
    sigmas[1] = σ0
    for ii=2:nObs
        sigmas[ii] = sqrt(κ + α*centered[ii-1]^2 +
                          β*sigmas[ii-1]^2)
    end
    return sigmas
end

@doc doc"""
Retrieve hidden sigma series from given dataset and given GARCH model.
Initial sigma either can be manually specified or it will be set
automatically.
"""->
function getSigmas(mod::GARCH_1_1, data::Array{Float64, 1},
                   initSigma::Float64 = NaN)
    μ, κ, α, β = mod.μ, mod.κ, mod.α, mod.β
    ## automatically get initial sigma if not specified
    if isnan(initSigma)
        ## alternative automatic setting:
        ## initSigma = std(data)
        initSigma = sqrt(κ/(1 - α - β))
    end
    params = [μ, κ, α, β, initSigma]
    return getSigmas(params, data)
end

##########
## nllh ##
##########

@doc doc"""
Negative log-likelihood of GARCH(1,1) with normally distributed
innovations.
"""->
function garch_norm_nllh(params::Array{Float64, 1},
                         data::Array{Float64, 1})
    ## get sigmas
    sigmas = getSigmas(params, data)

    ## normalize data
    centered = data - params[1]

    return sum(0.5*log(sigmas.^2*2*pi) + 0.5*(centered.^2./sigmas.^2))
end

@doc doc"""
Negative log-likelihood of GARCH(1,1) with t distributed innovations.
"""->
function garch_t_nllh(params::Array{Float64, 1},
                      data::Array{Float64, 1})
    μ, κ, α, β, ν, σ0 = params
    if isnan(ν)
        ## warn("\nν is NaN")
        return NaN
    end
    nObs = length(data)

    ## get sigmas
    sigmas = getSigmas([μ, κ, α, β, σ0], data)

    ## normalize data
    centered = data - μ

    ## get standard deviation of unscaled t innovations
    stdT = sqrt(ν/(ν-2))

    ## get scaling factor
    a = 1/stdT

    llhs = Float64[log(pdf(TLSDist(ν, 0.0, a*sigmas[ii]), centered[ii])) for ii=1:nObs]

    ## println(ν)
    return -sum(llhs)
end

##############
## estimate ##
##############

@doc doc"""
Estimate GARCH coefficients with normally distributed innovations. By
default, the estimation includes μ and the initial volatility σ0 but
it will not be returned.
"""->
function estimate(dt::Type{GARCH_1_1{Normal}}, data::Array{Float64, 1})
    μ, κ, α, β, σ0 = garchFit_norm(data)[1]
    return GARCH_1_1(μ, κ, α, β, Normal())
end

@doc doc"""
Estimate GARCH coefficients with t distributed innovations. By
default, the estimation includes μ and the initial volatility σ0 but
it will not be returned.
"""->
function estimate(dt::Type{GARCH_1_1{TDist}}, data::Array{Float64, 1})
    μ, κ, α, β, ν, σ0 = garchFit_t(data)[1]
    return GARCH_1_1(μ, κ, α, β, TDist(ν))
end

@doc doc"""
Estimate GARCH coefficients with normally distributed innovations. By
default, the estimation includes μ and the initial volatility σ0. σ0
also will be returned.
"""->
function garchFit_norm(data::Array{Float64, 1})
    opt = Opt(:LN_COBYLA, 5)

    function objFun(x::Vector, grad::Vector)
        if length(grad) > 0
            ## no partial derivative given
        end
        return garch_norm_nllh(x, data)
    end

    min_objective!(opt, objFun)
    xtol_rel!(opt, 1e-4)

    ## inequality constraints
    ineqConstraint(params, g) = params[3] + params[4] - 0.9999
    inequality_constraint!(opt, ineqConstraint, 1e-4)

    ## μ, κ, α, β, σ0
    lower_bounds!(opt, [-Inf, 0.00001, 0, 0, 0])
    upper_bounds!(opt, [Inf, Inf, Inf, Inf, Inf])

    ## come up with initial values
    initSigma = std(data)
    initVal = [mean(data), initSigma^2, 0.2, 0.8, initSigma]
    (minf, minx, ret) = optimize(opt, initVal)

    return (minx, minf)
end

@doc doc"""
Estimate GARCH coefficients with t distributed innovations. By
default, the estimation includes μ and the initial volatility σ0. σ0
also will be returned.
"""->
function garchFit_t(data::Array{Float64, 1})
    opt = Opt(:LN_COBYLA, 6)

    function objFun(x::Vector, grad::Vector)
        if length(grad) > 0
            ## no partial derivative given
        end
        return garch_t_nllh(x, data)
    end

    min_objective!(opt, objFun)
    xtol_rel!(opt, 1e-4)

    ## inequality constraints
    ineqConstraint(params, g) = params[3] + params[4] - 0.9999
    inequality_constraint!(opt, ineqConstraint, 1e-4)

    ## μ, κ, α, β, ν, σ0
    lower_bounds!(opt, [-Inf, 0.00001, 0, 0, 2.8, 0])
    upper_bounds!(opt, [Inf, Inf, Inf, Inf, Inf, Inf])

    ## come up with initial values
    initSigma = std(data)
    initVal = [mean(data), initSigma^2, 0.2, 0.8, 8.0, initSigma]
    (minf, minx, ret) = optimize(opt, initVal)

    return (minx, minf)
end

#########
## fit ##
#########

function fit(dt::Type{GARCH_1_1{Normal}}, data::Array{Float64, 1})
    optOutput = garchFit_norm(data)
    μ, κ, α, β, σ0 = optOutput[1]
    nllh = optOutput[2]

    ## retrieve sigma series
    sigmas = getSigmas([μ, κ, α, β, σ0], data)

    ## construct garch model
    garchMod = GARCH_1_1(μ, κ, α, β, Normal())

    return GARCH_1_1_Fit(garchMod, data, sigmas, nllh)
end

function fit(dt::Type{GARCH_1_1{TDist}}, data::Array{Float64, 1})
    optOutput = garchFit_t(data)
    μ, κ, α, β, ν, σ0 = optOutput[1]
    nllh = optOutput[2]

    ## retrieve sigma series
    sigmas = getSigmas([μ, κ, α, β, σ0], data)

    ## construct garch model
    garchMod = GARCH_1_1(μ, κ, α, β, TDist(ν))

    return GARCH_1_1_Fit(garchMod, data, sigmas, nllh)
end

## fit with Timenum / Timematr
##----------------------------

function fit(dt::Type{GARCH_1_1{Normal}}, tn::AbstractTimematr)
    nObs, nAss = size(tn)

    if nAss != 1
        error("Original data must be univariate for univariate models.")
    end

    ## extract data
    data = asArr(tn[1], Float64)[:]

    gFit = fit(dt, data)

    ## determine name for sigma column
    colName = symbol(string(names(tn.vals)[1], "_sigmas"))

    sigmasTm = Timematr(gFit.sigmas, [colName], idx(tn))

    return GARCH_1_1_Fit(gFit.model, tn, sigmasTm, gFit.nllh)
end

function fit(dt::Type{GARCH_1_1{TDist}}, tn::AbstractTimematr)
    nObs, nAss = size(tn)

    if nAss != 1
        error("Original data must be univariate for univariate models.")
    end

    ## extract data
    data = asArr(tn[1], Float64)[:]

    gFit = fit(dt, data)

    ## determine name for sigma column
    colName = symbol(string(names(tn.vals)[1], "_sigmas"))

    sigmasTm = Timematr(gFit.sigmas, [colName], idx(tn))

    return GARCH_1_1_Fit(gFit.model, tn, sigmasTm, gFit.nllh)
end

function fit(dt::Type{GARCH_1_1{Normal}}, tn::AbstractTimenum)
    nObs, nAss = size(tn)

    if nAss != 1
        error("Original data must be univariate for univariate models.")
    end

    ## extract data
    data = asArr(tn[1], Float64, NaN)[:]

    gFit = fit(dt, data[!isnan(data), 1])

    ## fill sigmas
    noNAInds = !isna(tn.vals[1])
    sigmas = DataArray(Float64, nObs)
    sigmas[noNAInds] = gFit.sigmas

    sigmasDf = DataFrame()
    sigmasDf[1] = sigmas

    ## determine name for sigma column
    colName = symbol(string(names(tn.vals)[1], "_sigmas"))
    rename!(sigmasDf, :x1, colName)

    sigmasTn = Timenum(sigmasDf, idx(tn))

    return GARCH_1_1_Fit(gFit.model, tn, sigmasTn, gFit.nllh)
end

function fit(dt::Type{GARCH_1_1{TDist}}, tn::AbstractTimenum)
    nObs, nAss = size(tn)

    if nAss != 1
        error("Original data must be univariate for univariate models.")
    end

    ## extract data
    data = asArr(tn[1], Float64, NaN)[:]

    gFit = fit(dt, data[!isnan(data), 1])

    ## fill sigmas
    noNAInds = !isna(tn.vals[1])
    sigmas = DataArray(Float64, nObs)
    sigmas[noNAInds] = gFit.sigmas

    sigmasDf = DataFrame()
    sigmasDf[1] = sigmas

    ## determine name for sigma column
    colName = symbol(string(names(tn.vals)[1], "_sigmas"))
    rename!(sigmasDf, :x1, colName)

    sigmasTn = Timenum(sigmasDf, idx(tn))

    return GARCH_1_1_Fit(gFit.model, tn, sigmasTn, gFit.nllh)
end


text:
struct BikeAttributes
    make :: String
    model :: String
    weight :: Float16
    drag_coefficient :: Float16
end

ext: jl


text:
# Edgeworth expansion approximations of sums, means and z statistics of iid variables.
# Quantiles are computed via the Cornish-Fisher expansion
# TODO: make the expansion of arbitrary order, using cumulants and Hermite polynomials

# Edgeworth approximation of the Z statistic
# EdgeworthSum and EdgeworthMean are both defined in terms of this
abstract EdgeworthAbstract <: ContinuousUnivariateDistribution

skewness(d::EdgeworthAbstract) = skewness(d.dist) / sqrt(d.n)
kurtosis(d::EdgeworthAbstract) = kurtosis(d.dist) / d.n

immutable EdgeworthZ{D<:UnivariateDistribution} <: EdgeworthAbstract
    dist::D
    n::Float64

    function EdgeworthZ{T<:UnivariateDistribution}(d::T, n::Real)
        @check_args(EdgeworthZ, n > zero(n))
        new(d, n)
    end
end
EdgeworthZ(d::UnivariateDistribution,n::Real) = EdgeworthZ{typeof(d)}(d,n)

mean(d::EdgeworthZ) = 0.0
var(d::EdgeworthZ) = 1.0

function pdf(d::EdgeworthZ, x::Float64)
    s = skewness(d)
    k = kurtosis(d)
    x2 = x*x
    pdf(Normal(0,1),x)*(1 + s*x*(x2-3.0)/6.0 +
                        k*(x2*(x2-6.0)+3.0)/24.0
                        + s*s*(x2*(x2*(x2-15.0)+45.0)-15.0)/72.0)
end
function cdf(d::EdgeworthZ, x::Float64)
    s = skewness(d)
    k = kurtosis(d)
    x2 = x*x
    cdf(Normal(0,1),x) -
    pdf(Normal(0,1),x)*(
                        s*(x2-1.0)/6.0 +
                        k*x*(x2-3.0)/24.0
                        + s*s*x*(x2*(x2-10.0)+15.0)/72.0)
end
function ccdf(d::EdgeworthZ, x::Float64)
    s = skewness(d)
    k = kurtosis(d)
    x2 = x*x
    ccdf(Normal(0,1),x) +
    pdf(Normal(0,1),x)*(
                        s*(x2-1.0)/6.0 +
                        k*x*(x2-3.0)/24.0
                        + s*s*x*(x2*(x2-10.0)+15.0)/72.0)
end

# Cornish-Fisher expansion.
function quantile(d::EdgeworthZ, p::Float64)
    s = skewness(d)
    k = kurtosis(d)
    z = quantile(Normal(0,1),p)
    z2 = z*z
    z + s*(z2-1)/6.0 + k*z*(z2-3)/24.0 - s*s/36.0*z*(2.0*z2-5.0)
end

function cquantile(d::EdgeworthZ, p::Float64)
    s = skewness(d)
    k = kurtosis(d)
    z = cquantile(Normal(0,1),p)
    z2 = z*z
    z + s*(z2-1)/6.0 + k*z*(z2-3)/24.0 - s*s/36.0*z*(2.0*z2-5.0)
end

# Edgeworth approximation of the sum
immutable EdgeworthSum{D<:UnivariateDistribution} <: EdgeworthAbstract
    dist::D
    n::Float64
    function EdgeworthSum{T<:UnivariateDistribution}(d::T, n::Real)
        @check_args(EdgeworthSum, n > zero(n))
        new(d, n)
    end
end
EdgeworthSum(d::UnivariateDistribution, n::Real) = EdgeworthSum{typeof(d)}(d,n)

mean(d::EdgeworthSum) = d.n*mean(d.dist)
var(d::EdgeworthSum) = d.n*var(d.dist)

# Edgeworth approximation of the mean
immutable EdgeworthMean{D<:UnivariateDistribution} <: EdgeworthAbstract
    dist::D
    n::Float64
    function EdgeworthMean{T<:UnivariateDistribution}(d::T, n::Real)
        # although n would usually be an integer, no methods are require this
        n > zero(n) ||
            error("n must be positive")
        new(d, Float64(n))
    end
end
EdgeworthMean(d::UnivariateDistribution,n::Real) = EdgeworthMean{typeof(d)}(d,n)

mean(d::EdgeworthMean) = mean(d.dist)
var(d::EdgeworthMean) = var(d.dist) / d.n

function pdf(d::EdgeworthAbstract, x::Float64)
    m, s = mean(d), std(d)
    pdf(EdgeworthZ(d.dist,d.n),(x-m)/s)/s
end

cdf(d::EdgeworthAbstract, x::Float64) = cdf(EdgeworthZ(d.dist,d.n), (x-mean(d))/std(d))

ccdf(d::EdgeworthAbstract, x::Float64) = ccdf(EdgeworthZ(d.dist,d.n), (x-mean(d))/std(d))

quantile(d::EdgeworthAbstract, p::Float64) = mean(d) + std(d)*quantile(EdgeworthZ(d.dist,d.n), p)
cquantile(d::EdgeworthAbstract, p::Float64) = mean(d) + std(d)*cquantile(EdgeworthZ(d.dist,d.n), p)


text:
using Test, CSV, DataFrames, StanQuap

ProjDir = @__DIR__

df = CSV.read(joinpath(ProjDir,  "Howell1.csv"), DataFrame)
df = filter(row -> row[:age] >= 18, df);

stan4_1 = "
// Inferring the mean and std
data {
  int N;
  real<lower=0> h[N];
}
parameters {
  real<lower=0.1> sigma;
  real<lower=100,upper=180> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 20);
  sigma ~ exponential(1);

  // Observed heights
  h ~ normal(mu, sigma);
}
";

data = (N = size(df, 1), h = df.height)
init = (mu = 160.0, sigma = 10.0)

qm, sm, om = stan_quap("s4.1s", stan4_1; data, init)
#qm, sm, om = stan_quap("s4.1s", stan4_1; data)

println()
qm |> display
println()

println()
om.optim |> display
println()


text:
##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test01.jl"))
omib_sys_file = System(PowerModelsData(omib_file_dir), runchecks = false)

##################################################
############### SOLVE PROBLEM ####################
##################################################
#Define Fault: Change of YBus
Ybus_change = NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault,
) #New YBus

@testset "Prog meter enabling" begin
    @test !PSID._prog_meter_enabled()
end

@testset "Make Simulation" begin
    path1 = (joinpath(pwd(), "test-Base-1"))
    !isdir(path1) && mkdir(path1)
    try
        #Define Simulation Problem
        sim = Simulation(
            ResidualModel,
            omib_sys, #system
            path1,
            (0.0, 30.0), #time span
            Ybus_change;
            system_to_file = true,
        )
        @test sim.status == PSID.BUILT
        # Test accessor functions
        dyn_wrapper = PowerSimulationsDynamics.get_dynamic_wrapper(sim, "generator-102-1")
        @test isa(dyn_wrapper, PowerSimulationsDynamics.DynamicWrapper)
        dic_init_conds = read_initial_conditions(sim)
        @test isa(dic_init_conds, Dict)
        dic_control_refs = get_setpoints(sim)
        @test isa(dic_control_refs, Dict)

        o_system = System(joinpath(path1, "input_system.json"))
        for b in get_components(Bus, o_system)
            b_sys = get_component(Bus, omib_sys, get_name(b))
            b_file = get_component(Bus, omib_sys_file, get_name(b))
            @test get_angle(b) == get_angle(b_sys)
            @test get_angle(b) == get_angle(b_file)
        end

    finally
        @info("removing test files")
        rm(path1, force = true, recursive = true)
    end

    path2 = (joinpath(pwd(), "test-Base-2"))
    !isdir(path2) && mkdir(path2)
    try
        #Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            omib_sys, #system
            path2,
            (0.0, 30.0), #time span
            Ybus_change;
            system_to_file = true,
        )
        @test sim.status == PSID.BUILT
        m_system = System(joinpath(path2, "initialized_system.json"))
        for b in get_components(Bus, m_system)
            b_sys = get_component(Bus, omib_sys, get_name(b))
            b_file = get_component(Bus, omib_sys_file, get_name(b))
            @test get_angle(b) == get_angle(b_sys)
            if get_bustype(b) == PSY.BusTypes.REF
                @test get_angle(b) == get_angle(b_file)
            else
                @test get_angle(b) != get_angle(b_file)
            end
        end

    finally
        @info("removing test files")
        rm(path2, force = true, recursive = true)
    end
end

@testset "Hybrid Line Indexing" begin
    ## Create threebus system with more dyn lines ##
    three_bus_file_dir = joinpath(TEST_FILES_DIR, "data_tests/ThreeBusInverter.raw")
    threebus_sys_dyns = System(three_bus_file_dir, runchecks = false)
    add_source_to_ref(threebus_sys_dyns)

    dyn_branch12 =
        DynamicBranch(get_component(Branch, threebus_sys_dyns, "BUS 1-BUS 2-i_2"))
    dyn_branch23 =
        DynamicBranch(get_component(Branch, threebus_sys_dyns, "BUS 2-BUS 3-i_3"))

    add_component!(threebus_sys_dyns, dyn_branch12)
    add_component!(threebus_sys_dyns, dyn_branch23)

    # Attach dyn devices
    for g in get_components(Generator, threebus_sys_dyns)
        if get_number(get_bus(g)) == 102
            case_gen = dyn_gen_second_order(g)
            add_component!(threebus_sys_dyns, case_gen, g)
        elseif get_number(get_bus(g)) == 103
            case_inv = inv_case78(g)
            add_component!(threebus_sys_dyns, case_inv, g)
        end
    end

    # Tests for all Dynamic Lines
    sim = Simulation(ResidualModel, threebus_sys_dyns, mktempdir(), (0.0, 10.0))
    @test sim.status == PSID.BUILT
    sim_inputs = sim.inputs
    DAE_vector = PSID.get_DAE_vector(sim_inputs)
    @test all(DAE_vector)
    @test all(LinearAlgebra.diag(sim_inputs.mass_matrix) .> 0)
    total_shunts = PSID.get_total_shunts(sim_inputs)
    # Total shunts matrix follows same pattern as the rectangular Ybus
    for v in LinearAlgebra.diag(total_shunts[1:3, 4:end])
        @test v > 0
    end

    for v in LinearAlgebra.diag(total_shunts[4:end, 1:3])
        @test v < 0
    end

    for entry in LinearAlgebra.diag(sim_inputs.mass_matrix)
        @test entry > 0
    end
    voltage_buses_ix = PSID.get_voltage_buses_ix(sim_inputs)
    @test length(voltage_buses_ix) == 3
    @test isempty(PSID.get_current_buses_ix(sim_inputs))

    # Tests for dynamic lines with b = 0
    set_b!(dyn_branch12, (from = 0.0, to = 0.0))
    set_b!(dyn_branch23, (from = 0.0, to = 0.0))
    sim = Simulation(ResidualModel, threebus_sys_dyns, pwd(), (0.0, 10.0))
    @test sim.status == PSID.BUILT
    sim_inputs = sim.inputs
    DAE_vector = PSID.get_DAE_vector(sim_inputs)
    @test sum(.!DAE_vector) == 6
    for (ix, entry) in enumerate(DAE_vector)
        if !entry
            @test LinearAlgebra.diag(sim_inputs.mass_matrix)[ix] == 0
        elseif entry
            @test LinearAlgebra.diag(sim_inputs.mass_matrix)[ix] == 1
        else
            @test false
        end
    end

    total_shunts = PSID.get_total_shunts(sim_inputs)
    @test sum(LinearAlgebra.diag(total_shunts)) == 0.0
    voltage_buses_ix = PSID.get_voltage_buses_ix(sim_inputs)
    @test isempty(voltage_buses_ix)
end

@testset "Test initial conditions externally" begin
    #Create system with 2 GENROU generators
    threebus_file_dir =
        joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROU/ThreeBusMulti_LessLoad.raw")
    dyr_file_dir =
        joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROU/ThreeBus_GENROU_SEXS.dyr")
    sys = System(threebus_file_dir, dyr_file_dir)
    x0_test = zeros(23)
    x0_test[1:6] .= 1.0
    #Initialize System from given initial condition
    sim = Simulation(
        ResidualModel,
        sys,
        pwd(),
        (0.0, 20.0),
        BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1");
        initialize_simulation = false,
        initial_conditions = x0_test,
    )
    @test LinearAlgebra.norm(sim.x0_init - x0_test) <= 1e-6
    #Initialize System normally
    sim_normal = Simulation(
        ResidualModel,
        sys,
        pwd(),
        (0.0, 20.0),
        BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"),
    )
    #Save states without generator at bus 2
    x0 = sim_normal.x0_init
    x0_no_gen = vcat(x0[1:6], x0[13:end])
    #Make Generator 2 unavailable and transform bus into PQ bus
    gen = PSY.get_component(ThermalStandard, sys, "generator-102-1")
    PSY.set_available!(gen, false)
    b = PSY.get_component(Bus, sys, "BUS 2")
    PSY.set_bustype!(b, PSY.BusTypes.PQ)
    #Create Simulation without Gen 2 starting from steady-state with Gen 2
    sim_trip_gen = Simulation(
        ResidualModel,
        sys,
        pwd(),
        (0.0, 20.0);
        initialize_simulation = false,
        initial_conditions = x0_no_gen,
    )
    @test LinearAlgebra.norm(sim_trip_gen.x0_init - x0_no_gen) <= 1e-6
    #Create Simulation without Gen 2 at steady state
    sim_normal_no_gen = Simulation(
        ResidualModel,
        sys,
        pwd(),
        (0.0, 20.0),
        BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"),
    )
    @test length(sim_normal_no_gen.x0_init) == 17
    #Ignore Initial Conditions without passing initialize_simulation = false
    sim_ignore_init =
        Simulation(ResidualModel, sys, pwd(), (0.0, 20.0); initial_conditions = x0_no_gen)
    @test LinearAlgebra.norm(sim_ignore_init.x0_init - sim_normal_no_gen.x0_init) <= 1e-6
    #Pass wrong vector size
    x0_wrong = zeros(20)
    @test_logs (:error, "Build failed") match_mode = :any Simulation(
        ResidualModel,
        sys,
        pwd(),
        (0.0, 20.0);
        initialize_simulation = false,
        initial_conditions = x0_wrong,
    )
    #Flat start initialization
    sim_flat =
        Simulation(ResidualModel, sys, pwd(), (0.0, 20.0); initialize_simulation = false)
    x0_flat = zeros(17)
    x0_flat[1:3] .= 1.0
    @test LinearAlgebra.norm(sim_flat.x0_init - x0_flat) <= 1e-6
end

@testset "Test Network Kirchoff Calculation" begin
    voltages = [
        1.05
        1.019861574381897
        0.0
        -0.016803841801155062
    ]
    V_r = voltages[1:2]
    V_i = voltages[3:end]
    ybus_ = PSY.Ybus(omib_sys).data
    I_balance_ybus = -1 * ybus_ * (V_r + V_i .* 1im)
    inputs = PSID.SimulationInputs(ResidualModel, omib_sys, ConstantFrequency)
    I_balance_sim = zeros(4)
    PSID.network_model(inputs, I_balance_sim, voltages)
    for i in 1:2
        @test isapprox(real(I_balance_ybus[i]), I_balance_sim[i])
        @test isapprox(imag(I_balance_ybus[i]), I_balance_sim[i + 2])
    end
end

@testset "Test Network Modification Callback Affects" begin
    three_bus_file_dir = joinpath(TEST_FILES_DIR, "data_tests/ThreeBusInverter.raw")
    threebus_sys = System(three_bus_file_dir, runchecks = false)
    add_source_to_ref(threebus_sys)
    # Attach dyn devices
    for g in get_components(Generator, threebus_sys)
        if get_number(get_bus(g)) == 102
            case_gen = dyn_gen_second_order(g)
            add_component!(threebus_sys, case_gen, g)
        elseif get_number(get_bus(g)) == 103
            case_inv = inv_case78(g)
            add_component!(threebus_sys, case_inv, g)
        end
    end

    ybus_original = PSY.Ybus(threebus_sys)

    inputs = PSID.SimulationInputs(ResidualModel, threebus_sys, ConstantFrequency)

    for i in 1:3, j in 1:3
        complex_ybus = ybus_original.data[i, j]
        @test inputs.ybus_rectangular[i, j] == real(complex_ybus)
        @test inputs.ybus_rectangular[i + 3, j + 3] == real(complex_ybus)
        @test inputs.ybus_rectangular[i + 3, j] == -imag(complex_ybus)
        @test inputs.ybus_rectangular[i, j + 3] == imag(complex_ybus)
    end

    br = get_component(Line, threebus_sys, "BUS 1-BUS 3-i_1")
    PSID.ybus_update!(inputs, br, -1.0)

    remove_component!(threebus_sys, br)
    ybus_line_trip = PSY.Ybus(threebus_sys)

    # Use is approx because the inversion of complex might be different in
    # floating point than the inversion of a single float
    for i in 1:3, j in 1:3
        complex_ybus = ybus_line_trip.data[i, j]
        @test isapprox(inputs.ybus_rectangular[i, j], real(complex_ybus), atol = 1e-10)
        @test isapprox(
            inputs.ybus_rectangular[i + 3, j + 3],
            real(complex_ybus),
            atol = 1e-10,
        )
        @test isapprox(inputs.ybus_rectangular[i + 3, j], -imag(complex_ybus), atol = 1e-10)
        @test isapprox(inputs.ybus_rectangular[i, j + 3], imag(complex_ybus), atol = 1e-10)
    end

    threebus_sys = System(three_bus_file_dir, runchecks = false)
    ybus_original = PSY.Ybus(threebus_sys)
    cb1 = NetworkSwitch(1.0, ybus_original)

    @test all(
        cb1.ybus_rectangular .== NetworkSwitch(1.0, ybus_original.data).ybus_rectangular,
    )

    for i in 1:3, j in 1:3
        complex_ybus = ybus_original.data[i, j]
        @test cb1.ybus_rectangular[i, j] == real(complex_ybus)
        @test cb1.ybus_rectangular[i + 3, j + 3] == real(complex_ybus)
        @test cb1.ybus_rectangular[i + 3, j] == -imag(complex_ybus)
        @test cb1.ybus_rectangular[i, j + 3] == imag(complex_ybus)
    end
end

@testset "Test Generation perturbations callback affects" begin
    three_bus_file_dir = joinpath(TEST_FILES_DIR, "data_tests/ThreeBusInverter.raw")
    threebus_sys = System(three_bus_file_dir, runchecks = false)
    add_source_to_ref(threebus_sys)
    # Attach dyn devices
    for g in get_components(Generator, threebus_sys)
        if get_number(get_bus(g)) == 102
            case_gen = dyn_gen_second_order(g)
            add_component!(threebus_sys, case_gen, g)
        elseif get_number(get_bus(g)) == 103
            case_inv = inv_case78(g)
            add_component!(threebus_sys, case_inv, g)
        end
    end

    mach = get_component(DynamicGenerator, threebus_sys, "generator-102-1")
    inv = get_component(DynamicInverter, threebus_sys, "generator-103-1")

    cref = ControlReferenceChange(1.0, mach, :P_ref, 10.0)
    ωref = ControlReferenceChange(1.0, inv, :ω_ref, 0.9)

    inputs = PSID.SimulationInputs(ResidualModel, threebus_sys, ConstantFrequency)
    integrator_for_test = MockIntegrator(inputs)

    cref_affect_f = PSID.get_affect(inputs, threebus_sys, cref)
    ωref_affect_f = PSID.get_affect(inputs, threebus_sys, ωref)

    cref_affect_f(integrator_for_test)
    ωref_affect_f(integrator_for_test)

    @test PSID.get_P_ref(inputs.dynamic_injectors[1]) == 10.0
    @test PSID.get_ω_ref(inputs.dynamic_injectors[2]) == 0.9

    inputs = PSID.SimulationInputs(ResidualModel, threebus_sys, ConstantFrequency)
    integrator_for_test = MockIntegrator(inputs)

    mach_trip = PSID.GeneratorTrip(1.0, mach)
    inv_trip = PSID.GeneratorTrip(1.0, inv)

    mtrip_affect_f = PSID.get_affect(inputs, threebus_sys, mach_trip)
    itrip_affect_f = PSID.get_affect(inputs, threebus_sys, inv_trip)

    mtrip_affect_f(integrator_for_test)
    itrip_affect_f(integrator_for_test)

    @test PSID.get_connection_status(inputs.dynamic_injectors[1]) == 0.0
    @test PSID.get_connection_status(inputs.dynamic_injectors[2]) == 0.0
end

@testset "Test Load perturbations callback affects" begin
    three_bus_file_dir = joinpath(TEST_FILES_DIR, "data_tests/ThreeBusInverter.raw")
    threebus_sys = System(three_bus_file_dir, runchecks = false)
    add_source_to_ref(threebus_sys)
    # Attach dyn devices
    for g in get_components(Generator, threebus_sys)
        if get_number(get_bus(g)) == 102
            case_gen = dyn_gen_second_order(g)
            add_component!(threebus_sys, case_gen, g)
        elseif get_number(get_bus(g)) == 103
            case_inv = inv_case78(g)
            add_component!(threebus_sys, case_inv, g)
        end
    end

    load_1 = get_component(ElectricLoad, threebus_sys, "load1011")
    load_2 = get_component(ElectricLoad, threebus_sys, "load1021")

    load_val = LoadChange(1.0, load_1, :P_ref, 10.0)
    load_trip = LoadTrip(1.0, load_2)

    inputs = PSID.SimulationInputs(ResidualModel, threebus_sys, ConstantFrequency)
    integrator_for_test = MockIntegrator(inputs)

    lref_affect_f = PSID.get_affect(inputs, threebus_sys, load_val)
    ltrip_affect_f = PSID.get_affect(inputs, threebus_sys, load_trip)

    lref_affect_f(integrator_for_test)
    @test PSID.get_P_ref(inputs.static_loads[1]) == 10.0
    ltrip_affect_f(integrator_for_test)
    @test PSID.get_connection_status(inputs.static_loads[3]) == 0.0
end

@testset "Global Index" begin
    three_bus_file_dir = joinpath(TEST_FILES_DIR, "data_tests/ThreeBusInverter.raw")
    threebus_sys_dyns = System(three_bus_file_dir, runchecks = false)
    add_source_to_ref(threebus_sys_dyns)

    dyn_branch12 =
        DynamicBranch(get_component(Branch, threebus_sys_dyns, "BUS 1-BUS 2-i_2"))
    dyn_branch23 =
        DynamicBranch(get_component(Branch, threebus_sys_dyns, "BUS 2-BUS 3-i_3"))

    add_component!(threebus_sys_dyns, dyn_branch12)
    add_component!(threebus_sys_dyns, dyn_branch23)

    # Attach dyn devices
    for g in get_components(Generator, threebus_sys_dyns)
        if get_number(get_bus(g)) == 102
            case_gen = dyn_gen_second_order(g)
            add_component!(threebus_sys_dyns, case_gen, g)
        elseif get_number(get_bus(g)) == 103
            case_inv = inv_case78(g)
            add_component!(threebus_sys_dyns, case_inv, g)
        end
    end

    # Tests for all Dynamic Lines
    sim = Simulation(ResidualModel, threebus_sys_dyns, mktempdir(), (0.0, 10.0))
    global_index = PSID.make_global_state_map(sim.inputs)
    @test_throws ErrorException PSID.get_state_from_ix(global_index, 40)
    @test ("V_2", :R) == PSID.get_state_from_ix(global_index, 2)
    @test ("BUS 2-BUS 3-i_3", :Il_I) == PSID.get_state_from_ix(global_index, 10)
    @test ("generator-102-1", :Vr1) == PSID.get_state_from_ix(global_index, 16)
    @test ("generator-103-1", :ε_pll) == PSID.get_state_from_ix(global_index, 30)
end

@testset "Frequency Reference" begin
    sys = System(
        joinpath(TEST_FILES_DIR, "data_tests/240busWECC_2018_PSS32_fixed_shunts.raw"),
        joinpath(TEST_FILES_DIR, "data_tests/240busWECC_2018_PSS.dyr"),
        bus_name_formatter = x -> string(strip(x["name"])) * "-" * string(x["index"]),
    )
    sim = Simulation(
        ResidualModel,
        sys, #system
        mktempdir(),
        (0.0, 20.0), #time span
        # Not initialized to speed up the test
        initialize_simulation = false,
    )

    @test PSID.get_global_vars_update_pointers(sim.inputs)[1] != 0

    sim = Simulation(
        ResidualModel,
        sys, #system
        mktempdir(),
        (0.0, 20.0),
        # Not initialized to speed up the test
        initialize_simulation = false,
        frequency_reference = ConstantFrequency, #time span
    )

    @test PSID.get_global_vars_update_pointers(sim.inputs)[1] == 0

    ref_bus = get_component(Bus, sys, "TESLA-3933")
    devs_in_ref = get_components(StaticInjection, sys, x -> get_bus(x) == ref_bus)
    for d in devs_in_ref
        remove_component!(sys, get_dynamic_injector(d))
        remove_component!(sys, d)
    end

    @test_logs (:error, "Build failed") match_mode = :any Simulation(
        ResidualModel,
        sys, #system
        mktempdir(),
        (0.0, 20.0),
        # Not initialized to speed up the test
        initialize_simulation = false,
        frequency_reference = ConstantFrequency, #time span
    )
end

@testset "Line and Branch transformation" begin
    #Define Simulation Problem
    omib_sys_copy = deepcopy(omib_sys)
    sim = Simulation(
        ResidualModel,
        omib_sys_copy, #system
        mktempdir(),
        (0.0, 30.0); #time span
        all_lines_dynamic = true,
    )
    @test sim.status == PSID.BUILT
    inputs = PSID.get_simulation_inputs(sim)
    @test !isempty(PSID.get_dynamic_branches(inputs))
    @test sum(PSID.get_ybus(inputs)) == 0.0

    sim = Simulation(
        ResidualModel,
        omib_sys_copy, #system
        mktempdir(),
        (0.0, 30.0); #time span
        all_branches_dynamic = true,
    )
    @test sim.status == PSID.BUILT
    inputs = PSID.get_simulation_inputs(sim)
    @test !isempty(PSID.get_dynamic_branches(inputs))
    @test sum(PSID.get_ybus(inputs)) == 0.0

    sim = Simulation!(
        ResidualModel,
        omib_sys_copy, #system
        mktempdir(),
        (0.0, 30.0); #time span
        all_lines_dynamic = true,
    )
    @test sim.status == PSID.BUILT
    inputs = PSID.get_simulation_inputs(sim)
    @test !isempty(PSID.get_dynamic_branches(inputs))
    @test sum(PSID.get_ybus(inputs)) == 0.0
    @test !isempty(PSY.get_components(PSY.DynamicBranch, omib_sys_copy))
    @test isempty(PSY.get_components(PSY.Line, omib_sys_copy))
end


text:
@testset "DirectCodec" begin
    dc = DirectCodec(1.0)
    @test encode(dc, 3) == 3
    @test encode(dc, [-2, 0]) == [-2, 0]
    @test decode(dc, 5) == 5
    @test decode(dc, [1, 2]) == [1, 2]
end


text:
using InterpretableValidation
using Distributions
using ExprRules
using Random
using POMDPModels
using POMDPs
using POMDPSimulators
using POMDPPolicies

mdp = SimpleGridWorld(size = (9,9), rewards = Dict(GWPos(5,5) => 1, GWPos(6,5) => 0), tprob = .7, discount=1)

mvts = MvTimeseriesDistribution(:x => IID(collect(1:1.), Categorical(4)))
s_fn(s) = SymbolTable(:sx => s[1], :sy => s[2])
x_fn(v) = actions(mdp)[v[:x][1]]

# Demonstrate that a certain choice of policy obtains high reward
expr = Meta.parse("[(sx .== 5) .& (sy .<= 5), x .== 1] && [(sx .<= 4) .& (sy .<= 5), x .== 4] && [(sx .>= 6) .& (sy .<= 5), x .== 3]")
p = get_policy(expr, s_fn, x_fn, mvts)
mean([simulate(RolloutSimulator(), mdp, RandomPolicy(mdp)) for i=1:1000])
mean([simulate(RolloutSimulator(), mdp, FunctionPolicy(p)) for i=1:1000])

# Setup the grammar
x_samp_dist = default_comparison_distribution(mvts)
x_comp = default_comparisons(mvts)
s_samp_dist = Dict{Symbol, Distribution}(:sx => Categorical(9), :sy => Categorical(9))
all_syms = [Symbol(".<="), Symbol(".=="), Symbol(".>=")]
s_comp = Dict(:sx => all_syms, :sy => all_syms)
g = create_policy_grammar(1, x_samp_dist, x_comp, s_samp_dist, s_comp)

# Setup the loss function
lf = policy_loss_fn(mdp, s_fn, x_fn, mvts)

# optimize
optimize(()->nothing, mvts, Npop = 1000, Niter = 30, loss = lf, grammar = g, max_depth = 5)


text:
import Symbolics: get_variables; export get_variables
import Base.isequal; export isequal

# pretty-printing
display(var::HarmonicVariable) = display(var.name)
display(var::Vector{HarmonicVariable}) = display.(getfield.(var, Symbol("name")))

function _coordinate_transform(new_var, ω, t, type)
    coords = Dict([
        "u" =>  new_var * cos(ω*t),
        "v" => new_var * sin(ω*t),
        "a" => new_var])
    return coords[type]
end

function _create_harmonic_variable(nat_var::Num, ω::Num, t::Num, type::String; new_symbol::String)
    new_var = declare_variable(new_symbol, t) # this holds the internal symbol
    name = type * "_{" * var_name(nat_var) * "," * Base.replace(string(ω),"*"=>"") * "}"
    rule = _coordinate_transform(new_var, ω, t, type) # contribution of this harmonic variable to the natural variable
    hvar = HarmonicVariable(new_var, name, type, ω, nat_var)
    rule, hvar
end

###
# Functions for variable substutions and manipulation of HarmonicVariable
###

# when HV is used for substitute, substitute its symbol
substitute_all(eq::Union{Num, Equation}, rules::Dict{HarmonicVariable}) = substitute(eq, Dict(zip(getfield.(keys(rules), :symbol), values(rules))))

function substitute_all(var::HarmonicVariable, rules)
    sym, freq = var.symbol, var.ω
    HarmonicVariable(substitute_all(sym, rules), var.name, var.type, substitute_all(freq, rules), var.natural_variable)
end

substitute_all(vars::Vector{HarmonicVariable}, rules) = [substitute_all(var, rules) for var in vars]

"Returns the symbols of a `HarmonicVariable`."
get_variables(vars::Vector{Num}) = unique(flatten([Num.(get_variables(x)) for x in vars]))

get_variables(var::HarmonicVariable) = Num.(get_variables(var.symbol))

isequal(v1::HarmonicVariable, v2::HarmonicVariable) = isequal(v1.symbol, v2.symbol)


text:
using AccurateArithmetic
using Test

using AccurateArithmetic.Summation: accumulate
using AccurateArithmetic.Summation: sumAcc, compSumAcc, dotAcc, compDotAcc
using AccurateArithmetic.Test: generate_sum, generate_dot
using LinearAlgebra

@testset "AccurateArithmetic" begin
    @testset "Tests" begin
        @testset "generate_sum" begin
            @testset "vector length" begin
                for n in 100:110
                    let (x, s, c_) = generate_sum(n, 1e10)
                        @test length(x) == n
                    end
                end
            end

            @testset "condition number" begin
                for c in (1e10, 1e20)
                    let (x, s, c_) = generate_sum(101, c)
                        @test c/100 < c_ < 1000c
                    end
                end
            end
        end

        @testset "generate_dot" begin
            @testset "vector length" begin
                for n in 100:110
                    let (x, y, s, c_) = generate_dot(n, 1e10)
                        @test length(x) == n
                        @test length(y) == n
                    end
                end
            end

            @testset "condition number" begin
                for c in (1e10, 1e20)
                    let (x, y, s, c_) = generate_dot(101, c)
                        @test c/100 < c_ < 1000c
                    end
                end
            end
        end
    end

    @testset "summation" begin
        @testset "naive" begin
            for N in 100:110
                x = rand(N)
                ref = sum(x)
                @test ref ≈ sum_naive(x)

                acc = sumAcc
                @test ref ≈ accumulate((x,), acc, Val(:scalar), Val(2))
                @test ref ≈ accumulate((x,), acc, Val(:mask),   Val(2))
            end
        end

        @testset "compensated" begin
            for N in 100:110
                x, ref, _ = generate_sum(N, 1e10)
                @test ref == sum_oro(x)
                @test ref == sum_kbn(x)

                acc = compSumAcc(two_sum)
                @test ref == accumulate((x,), acc, Val(:scalar), Val(2))
                @test ref == accumulate((x,), acc, Val(:mask),   Val(2))
            end
        end
    end

    @testset "dot product" begin
        @testset "naive" begin
            for N in 100:110
                x = rand(N)
                y = rand(N)
                ref = dot(x, y)
                @test ref ≈ dot_naive(x, y)

                acc = dotAcc
                @test ref ≈ accumulate((x,y), acc, Val(:scalar), Val(2))
                @test ref ≈ accumulate((x,y), acc, Val(:mask),   Val(2))
            end
        end

        @testset "compensated" begin
            for N in 100:110
                x, y, ref, _ = generate_dot(N, 1e10)
                @test ref == dot_oro(x, y)

                acc = compDotAcc
                @test ref == accumulate((x,y), acc, Val(:scalar), Val(2))
                @test ref == accumulate((x,y), acc, Val(:mask),   Val(2))
            end
        end
    end
end

using BenchmarkTools

acc = compSumAcc(two_sum)

BenchmarkTools.DEFAULT_PARAMETERS.evals = 1000
println("\nsize 10_000")
print("  sum_oro"); @btime sum_oro($(rand(10_000)))
print("  sum_kbn"); @btime sum_kbn($(rand(10_000)))

BenchmarkTools.DEFAULT_PARAMETERS.evals = 10
println("\nsize 1_000_000")
print("  sum_oro"); @btime sum_oro($(rand(1_000_000)))
print("  sum_kbn"); @btime sum_kbn($(rand(1_000_000)))

println("\nsize 100_000_000")
x = rand(100_000_000)
print("  sum_oro"); @btime sum_oro($x)
print("  sum_kbn"); @btime sum_naive($x)
print("  sum    "); @btime sum($x)

println()
y = rand(100_000_000)
print("  dot_oro"); @btime dot_oro($x, $y)
print("  dot_kbn"); @btime dot_naive($x, $y)
print("  dot    "); @btime dot($x, $y)


text:
struct BSONConversionError <: Exception
    src_type::UInt8
    src_subtype::Union{Nothing, UInt8}
    dst_type::DataType
end

BSONConversionError(src_type::UInt8, dst_type::DataType) = BSONConversionError(src_type, nothing, dst_type)

struct BSONValidationError <: Exception
    msg::String
end


text:
import Base: (*), (==)

# ---------------------
# Concatenation
# ---------------------

*(a::Org) = a
*(a::Org, b::Org) = Org(vcat(a.contents, b.contents))
*(orgs::Org...) = Org(getproperty.(orgs, :content) |> Iterators.flatten |> collect)

# Section

*(hs::Heading...) = Org(hs)
*(a::Section, b::Section) = Section(vcat(a.contents, b.contents))
function *(s::Section, o::OrgComponent)
    sc = deepcopy(s)
    if o isa Object
        if sc.section.contents[end] isa Paragraph
            sc.section.contents[end].contents =
                vcat(sc.section.contents[end].contents, o)
        else
            sc.section.contents = vcat(sc.section.contents, Paragraph(o))
        end
    else
        sc.section.contents = vcat(sc.section.contents, o)
    end
    sc
end
function *(h::Heading, o::OrgComponent)
    hc = copy(h)
    if isnothing(hc.section)
        hc.section = Section([]) * o
    else
        hc.section = h.section * o
    end
    hc
end
function *(h::Heading, s::Section)
    hc = copy(h)
    hc.section = vcat(hc.section.contents, s.contents)
    hc
end

# Greater Element

*(components::Vector{Element}...) =
    Section([components])

# Element

*(a::Paragraph, b::Paragraph) = Paragraph([a.contents, b.contents])

*(a::Comment, b::Comment) = Comment(vcat(a.contents, b.contents))
*(a::FixedWidth, b::FixedWidth) = Comment(vcat(a.contents, b.contents))

*(rows::TableRow...) = Table(rows, [])

# Object

*(objects::Object...) = Paragraph(objects)
*(cells::TableCell...) = TableRow(cells)

function *(a::TextPlain{SubString}, b::TextPlain{SubString})
    if a.text.string === b.text.string &&
        a.text.offset + a.text.ncodeunits == b.text.offset
        TextPlain(@inbounds SubString(a.text.string, 1 + a.text.offset,
                                      b.text.offset + lastindex(b.text)))
    else
        TextPlain(a.text * b.text)
    end
end

*(a::TextPlain, b::TextPlain) = TextPlain(a.text * b.text)

# ---------------------
# Equality
# ---------------------

function ==(a::T, b::T) where {T <: OrgComponent}
    for f in fieldnames(T)
        if getproperty(a, f) != getproperty(b, f)
            return false
        end
    end
    true
end

==(a::Org, b::Org) = a.contents == b.contents
==(a::TextPlain, b::TextPlain) = a.text == b.text

## conversion

# ---------------------
# Terminality
# ---------------------

terminal(::Element) = false
terminal(::LesserElement) = true
terminal(::Clock) = false
terminal(::Paragraph) = false
terminal(::TableRow) = false

terminal(::Object) = true
terminal(::Citation) = false
terminal(::CitationReference) = false
# terminal(::Link) = false
# terminal(::RadioTarget) = false
terminal(::TableCell) = false
terminal(::TextMarkup{Vector{Object}}) = false

# ---------------------
# Iteration
# ---------------------

import Base: iterate, length

length(org::Org) = length(org.contents)
iterate(org::Org) =
    if length(org.contents) > 0
        (org.contents[1], 2)
    end
iterate(org::Org, index::Integer) =
    if index <= length(org.contents)
        (org.contents[index], index + 1)
    end

# Section

length(h::Heading) = length(h.title) + if isnothing(h.section) 0 else 1 end
iterate(h::Heading) = if length(h) > 0 (h.title[1], 2) end
iterate(h::Heading, index::Integer) =
    if index <= length(h.title)
        (h.title[index], index + 1)
    elseif index == 1 + length(h.title) && !isnothing(h.section)
        (h.section, index + 1)
    end

length(s::Section) = length(s.contents)
iterate(s::Section) = if length(s) > 0 (s.contents[1], 2) end
iterate(s::Section, index::Integer) =
    if index <= length(s.contents)
        (s.contents[index], index + 1)
    end

length(a::AffiliatedKeywordsWrapper) = 1 + length(a.keywords)
iterate(a::AffiliatedKeywordsWrapper) = (a.element, 1)
iterate(a::AffiliatedKeywordsWrapper, index::Integer) =
    if index <= length(a.keywords)
        (a.keywords[index], index + 1)
    end

# Greater Element

length(g::GreaterBlock) = length(g.contents)
iterate(g::GreaterBlock) = if length(g) > 0 (g.contents[1], 2) end
iterate(g::GreaterBlock, index::Integer) =
    if index <= length(g.contents)
        (g.contents[index], index + 1)
    end

length(d::Drawer) = length(d.contents)
iterate(d::Drawer) = if length(d) > 0 (d.contents[1], 2) end
iterate(d::Drawer, index::Integer) =
    if index <= length(d.contents)
        (d.contents[index], index + 1)
    end

length(d::DynamicBlock) = length(d.contents)
iterate(d::DynamicBlock) = if length(d) > 0 (d.contents[1], 2) end
iterate(d::DynamicBlock, index::Integer) =
    if index <= length(d.contents)
        (d.contents[index], index + 1)
    end

length(f::FootnoteDefinition) = length(f.definition)
iterate(f::FootnoteDefinition) =
    if length(f) > 0 (f.definition[1], 2) end
iterate(f::FootnoteDefinition, index::Integer) =
    if index <= length(f.definition)
        (f.definition[index], index + 1)
    end

length(l::List) = length(l.items)
iterate(l::List) = if length(l) > 0 (l.items[1], 2) end
iterate(l::List, index::Integer) =
    if index <= length(l.items)
        (l.items[index], index + 1)
    end

length(i::Item) = length(i.contents)
iterate(i::Item) = if length(i) > 0 (i.contents[1], 2) end
iterate(i::Item, index::Integer) =
    if index <= length(i.contents)
        (i.contents[index], index + 1)
    end

length(p::PropertyDrawer) = length(p.contents)
iterate(p::PropertyDrawer) = if length(p) > 0 (p.contents[1], 2) end
iterate(p::PropertyDrawer, index::Integer) =
    if index <= length(p.contents)
        (p.contents[index], index + 1)
    end

length(t::Table) = length(t.rows)
iterate(t::Table) = if length(t) > 0 (t.rows[1], 2) end
iterate(t::Table, index::Integer) =
    if index <= length(t.rows)
        (t.rows[index], index + 1)
    end

# Lesser Element

length(c::Clock) = 1
iterate(c::Clock) = (c.timestamp, 0)
iterate(::Clock, ::Integer) = nothing

length(p::Paragraph) = length(p.contents)
iterate(p::Paragraph) = if length(p) > 0 (p.contents[1], 2) end
iterate(p::Paragraph, index::Integer) =
    if index <= length(p.contents)
        (p.contents[index], index + 1)
    end

length(r::TableRow) = length(r.cells)
iterate(r::TableRow) = if length(r) > 0 (r.cells[1], 2) end
iterate(r::TableRow, index::Integer) =
    if index <= length(r.cells)
        (r.cells[index], index + 1)
    end

# Object

length(k::CitationReference) =
    length(k.prefix) + length(k.suffix)
iterate(k::CitationReference) =
    if length(k.prefix) > 0
         (k.prefix[1], 2)
     elseif length(k.suffix) > 0
         (k.suffix[1], 2)
     end
iterate(k::CitationReference, index::Integer) =
    if index <= length(k.prefix)
        (k.prefix[index], index + 1)
    elseif index <= length(k.prefix) + length(k.suffix)
        (k.suffix[index-length(k.prefix)], index + 1)
    end

length(c::Citation) =
    length(c.globalprefix) + length(c.citerefs) + length(c.globalsuffix)
iterate(c::Citation) =
    (if length(c.globalprefix) > 0
         c.globalprefix[1]
     else
         c.citerefs[1]
     end, 2)
iterate(c::Citation, index::Integer) =
    if index <= length(c.globalprefix)
        (c.globalprefix[index], index + 1)
    elseif index <= length(c.globalprefix) + length(c.citerefs)
        (c.citerefs[index-length(c.globalprefix)], index + 1)
    elseif index <= length(c.globalprefix) + length(c.citerefs) + length(c.globalsuffix)
        (c.globalsuffix[index-length(c.globalprefix)-length(c.citerefs)], index + 1)
    end

length(f::FootnoteReference{<:Any, Vector{Object}}) = length(f.definition)
iterate(f::FootnoteReference{<:Any, Vector{Object}}) =
    if length(f) > 0 (f.definition[1], 2) end
iterate(f::FootnoteReference{<:Any, Vector{Object}}, index::Integer) =
    if index <= length(f.definition)
        (f.definition[index], index + 1)
    end

length(c::TableCell) = length(c.contents)
iterate(c::TableCell) = if length(c) > 0 (c.contents[1], 2) end
iterate(c::TableCell, index::Integer) =
    if index <= length(c.contents)
        (c.contents[index], index + 1)
    end

length(t::TextMarkup{Vector{Object}}) = length(t.contents)
iterate(t::TextMarkup{Vector{Object}}) =
    if length(t) > 0 (t.contents[1], 2) end
iterate(t::TextMarkup{Vector{Object}}, index::Integer) =
    if index <= length(t.contents)
        (t.contents[index], index + 1)
    end

# More iterating

struct OrgIterator
    o::Org
end

Base.IteratorSize(::OrgIterator) = Base.SizeUnknown()
iterate(it::OrgIterator) =
    if length(it.o.contents) > 0
        el, state = iterate(it.o)
        (el, Vector{Tuple}([(it.o, state), (el, nothing)]))
    end
iterate(it::OrgIterator, stack::Vector) =
    if length(stack) > 0
        next = if isnothing(stack[end][2])
            iterate(stack[end][1])
        else
            iterate(stack[end][1], stack[end][2])
        end
        if isnothing(next)
            pop!(stack)
            iterate(it, stack)
        else
            el, state = next
            stack[end] = (stack[end][1], state)
            if ! terminal(el)
                push!(stack, (el, nothing))
            end
            (el, stack)
        end
    end

struct OrgElementIterator
    o::Org
end

Base.IteratorSize(::OrgElementIterator) = Base.SizeUnknown()
iterate(it::OrgElementIterator) =
    if length(it.o.contents) > 0
        el, state = iterate(it.o)
        (el, Vector{Tuple}([(it.o, state), (el, nothing)]))
    end
iterate(it::OrgElementIterator, stack::Vector) =
    if length(stack) > 0
        next = if isnothing(stack[end][2])
            iterate(stack[end][1])
        else
            iterate(stack[end][1], stack[end][2])
        end
        if isnothing(next)
            pop!(stack)
            iterate(it, stack)
        elseif next[1] isa Object
            el, state = next
            stack[end] = (stack[end][1], state)
            iterate(it, stack)
        else
            el, state = next
            stack[end] = (stack[end][1], state)
            if ! terminal(el)
                push!(stack, (el, nothing))
            end
            (el, stack)
        end
    end

# ---------------------
# Accessors
# ---------------------

import Base.getindex

getindex(o::Org, i::Integer) = o.contents[i]

function getindex(props::PropertyDrawer, name::AbstractString)
    additive = if endswith(name, '+')
        name = name[1:end-1]
        true
    else
        false
    end
    matches = filter(n -> n.name == name && n.additive == additive, props.contents)
    if length(matches) == 1
        matches[1].value
    else
        nothing
    end
end

function getindex(h::Heading, prop::AbstractString)
    if isnothing(h.properties)
        nothing
    else
        getindex(h.properties, prop)
    end
end


text:
import Eglob: eglob, eglobt

using Base.Test

include("inc_helpers.jl")

function run_eglob_tests(globber::Function)
    # Basic glob test against generated files.
    withtempfiles("foo/x.c", "bar/y.c", "foo/bar/baz.c") do tempdir, files
        cd(tempdir) do
            c_files = Set([abspath(f) for f in files if endswith(f, ".c")])
            @test Set(abspath.(globber("**/*.c"))) == c_files

            foo_c_files = Set([f for f in files if ismatch(r"^foo/.*\.c$", f)])
            @test Set(globber("f*/**/*.c")) == foo_c_files
        end
    end

    # Basic glob test against existing files.
    expected = Set(abspath.(shellglob("**/*.jl")))
    @test Set(abspath.(globber("**/*.jl"))) == expected

    # Glob test that is guaranteed not to match.
    withtempfiles("quux/foo", "quux/baz/bar") do tempdir, files
        @test isempty(globber("quux/**/*.c"))
    end

    # Case-sensitivity check.
    files = ["foo/first.cc", "bar/second.c", "bar/third.cc",
             "foo/bar/fourth.c", "foo/bar/baz/FIFTH.C", "bar/sixth.cc"]
    withtempfiles(files...) do tempdir,  files
      cd(tempdir) do
          expected = Set([abspath(f) for f in files if ismatch(r"\.[Cc]$", f)])
          @test Set(abspath.(globber("**/*.[Cc]"))) == expected

          expected = Set([abspath(f) for f in files if ismatch(r"^f.*\.c$", f)])
          @test Set(abspath.(globber("f*/**/*.c"))) == expected
      end
    end

    # Check with complicated patterns.
    files = ["aaa/kjhasd.txt", "aaa/bbb/ioua.txt", "abc/ccc/ddd/foo.c",
             "xyz/ddd/eee/bar.java", "abc/ddd/eee/baz.c", "aaaa/quux.java",
             "abcc/ccc/ddd/foo.c", "bbb/ccccccc/ddddddd/e/hello.java",
             "bbb/iswadf/abs.txt", "ghsadfkjhasdf/kd/8wequ/wog.txt",
             "ba/lkjjd/yuio/poiu/fred",
             "argh/this/is/a/silly/path/to/a//file.txt"]
    withtempfiles(files...) do tempdir, files
        cd(tempdir) do
            @test Set(globber("a??/*/[cd]*/**/*.c")) == Set([
                "abc/ccc/ddd/foo.c"
            ])
            @test Set(globber("[ab]*/**/*.java")) == Set([
                "aaaa/quux.java",
                "bbb/ccccccc/ddddddd/e/hello.java"
            ])
            @test Set(globber("[!a]*/**/*.txt")) == Set([
                "bbb/iswadf/abs.txt",
                "ghsadfkjhasdf/kd/8wequ/wog.txt"
            ])
            @test Set(globber("[ab][!ga]*/**/f*")) == Set([
                "abcc/ccc/ddd/foo.c",
                "abc/ccc/ddd/foo.c",
                "argh/this/is/a/silly/path/to/a/file.txt"
            ])
            @test Set(globber("?[!b]*/**/f*")) == Set([
                "argh/this/is/a/silly/path/to/a/file.txt",
                "ba/lkjjd/yuio/poiu/fred"
            ])
            @test Set(globber("[agx][!z]*/**/g*")) == Set{String}([])
        end
    end
end

@testset "eglob" begin
    run_eglob_tests(eglob)
end

@testset "eglobt" begin
    run_eglob_tests(eglob)
end


text:
# type Gate{T<:AbstractArray}
#     Wx::T
#     Wh::T
#     b::T
# end

# function lstm_cell(input::Gate, forget::Gate, output::Gate, candidate::Gate, ct, ht)
#     ft = sigmoid(xt * forget.Wx + ht * forget.Wh .+ forget.b)
#     it = sigmoid(xt * input.Wx + ht * input.Wh .+ input.b)
#     ot = sigmoid(xt * output.Wx + ht * output.Wh .+ output.b)
#     candidate_memory = tanh(xt * candidate.Wx + ht * candidate.Wh .+ candidate.b)
#     ct = ft * ct + it * candidate_memory
#     ht = ot * tanh(ct)
#     return (ht, ct, ht)
# end


text:
export antiplane_viscosity, mantle_viscosity
@doc raw"""
    antiplane_viscosity(ε::Union{PyObject, Array{Float64}}, σ::Union{PyObject, Array{Float64}},
    μ::Union{PyObject, Float64}, η::Union{PyObject, Float64}, Δt::Float64)

Calculates the stress at time $t_{n+1}$ given the strain at $t_{n+1}$ and stress at $t_{n}$. The governing equation is
```math
\dot\sigma + \frac{\mu}{\eta}\sigma = 2\mu \dot\epsilon
```
The discretization form is
```math
\sigma^{n+1} = \frac{1}{\frac{1}{\Delta t}+\frac{\mu}{\eta}}(2\mu\dot\epsilon^{n+1} + \frac{\sigma^n}{\Delta t})
```
"""
function antiplane_viscosity(dε::Union{PyObject, Array{Float64}}, σ::Union{PyObject, Array{Float64}},
        μ::Union{PyObject, Float64, Array{Float64}}, η::Union{PyObject, Float64, Array{Float64}}, Δt::Union{PyObject,Float64})
    dε, σ, μ, η = convert_to_tensor([dε, σ, μ,η], [Float64, Float64, Float64, Float64])
    if length(size(μ))==0
        η/(η/Δt + μ)*(2μ*dε+σ/Δt)
    elseif length(size(μ))==1
        repeat(η/(η/Δt + μ),1,2).*(2repeat(μ, 1, 2).*dε+σ/Δt)
    end
end

@doc raw"""
    mantle_viscosity(u::Union{Array{Float64}, PyObject},
        T::Union{Array{Float64}, PyObject}, m::Int64, n::Int64, h::Float64;
        σ_yield::Union{Float64, PyObject} = 300e6,
        ω::Union{Float64, PyObject},
        η_min::Union{Float64, PyObject} = 1e18,
        η_max::Union{Float64, PyObject} = 1e23,
        E::Union{Float64, PyObject} = 9.0,
        C::Union{Float64, PyObject} = 1000., N::Union{Float64, PyObject} = 2.)

```math
\eta = \eta_{\min} + \min\left( \frac{\sigma_{\text{yield}}}{2\sqrt{\epsilon_{II}}}, \omega\min(\eta_{\max}, \eta) \right)
```
with
```math
\epsilon_{II} = \frac{1}{2} \epsilon(u)\qquad \eta = C e^{E(0.5-T)} (\epsilon_{II})^{(1-n)/2n}
```

Here $\epsilon_{II}$ is the second invariant of the strain rate tensor, $C > 0$ is a viscosity pre-factor, $E > 0$ is the non-dimensional activation energy,
$n > 0$ is the nonlinear exponent, $η_\min$, $η_\max$ act as minimum and maximum bounds for the effective viscosity, and $σ_{\text{yield}} > 0$ is the yield
stress. $w\in (0, 1]$ is the weakening factor, which is used to incorporate phenomenological aspects that
cannot be represented in a purely viscous flow model, such as processes which govern mega-thrust faults along the subduction interface, or
partial melting near a mid-ocean ridge.

The viscosity of the mantle is governed by the high-temperature creep of silicates, for which laboratory experiments show that the creep
strength is temperature-, pressure-, compositional- and stress-dependent.

The output is a length $4mn$ vector.

!!! info
    See __Towards adjoint-based inversion of time-dependent mantle convection with nonlinear viscosity__ for details.
"""
function mantle_viscosity(u::Union{Array{Float64}, PyObject},
     T::Union{Array{Float64}, PyObject}, m::Int64, n::Int64, h::Float64;
     σ_yield::Union{Float64, PyObject} = 300e6,
     ω::Union{Float64, PyObject} = 1.0,
     η_min::Union{Float64, PyObject} = 1e-2,
     η_max::Union{Float64, PyObject} = 1e3,
     E::Union{Float64, PyObject} = 9.0,
     C::Union{Float64, PyObject} = 1000., N::Union{Float64, PyObject} = 2.)
    u, T, σ_yield, ω, η_min, η_max, E, C, N = convert_to_tensor(
        [u, T, σ_yield, ω, η_min, η_max, E, C, N], [Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64]
    )
    ε = eval_strain_on_gauss_pts(u, m, n, h)
    εII = 0.5 * sum(ε .* ε, dims = 2)
    T = reshape(repeat(T, 1, 4), (-1,))
    η = C * exp(E*(0.5 - T))*εII^((1-N)/2N)
    η_min + min(σ_yield/2/sqrt(εII), ω*min(η_max, η))
end


text:
module JModels

include("interface.jl")

end # module


text:
# test breaking a component and examine the deterministic costs
addprocs(30);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

N = 5;

caseList = [13,33,123];
Δt = 0.25;
iterMax = 20;
T = 24;

for ci in 1:length(caseList)
    dataList = Dict();
    τ = Int64(1/6*T);

    for ω in keys(pDistr.ωDistrn)
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T);
            if length(ω) == 1
                remotecall_fetch(breakComponent, j, fData, ω, "g");
            else
                remotecall_fetch(breakComponent, j, fData, ω, "l");
            end
        end

        # break the component and then calculate the deterministic cost
        detSol,detObj = detBuild(Δt, T, fData, bData, dData);
        dataList[ω] = [detSol,detOb];
    end

    save("breakResults_$(ci).jld","data",dataList);
end


text:
#!/usr/bin/env julia

#
# Task sink
# Binds PULL socket to tcp://localhost:5558
# Collects results from workers via that socket
#

using ZMQ

context = Context()

# Socket to receive messages on
receiver = Socket(context, PULL)
ZMQ.bind(receiver, "tcp://*:5558")

# Wait for start of batch
s = ZMQ.recv(receiver)

# Start our tic toc clock
tic()

# Process 100 confirmations
total_msec = 0
for task_nbr in [1:100]
    s = ZMQ.recv(receiver)
    if task_nbr % 10 == 0
        write(STDOUT, ":")
    else
        write(STDOUT, ".")
    end
    flush(STDOUT)
end

# Calculate and report duration of batch
tend = toq()
println("\nTotal elapsed time: $(tend * 1000) msec")

# Making a clean exit.
ZMQ.close(receiver)
ZMQ.close(context)

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

module TestBase

using Test
using RelationalDatasets
import Base.Filesystem: expanduser, joinpath

@testset "Base Functions" begin
    @test RelationalDatasets.get_data_home() == expanduser(joinpath("~", "relational_datasets"))
end

end # module


text:
function topoint2d(ift::ImageFeatures.Feature)
    return Point(ift.keypoint.I[2],ift.keypoint.I[1])
end

function drawfeaturecircle2d!(im, ift::Feature; radius=5, boundary=6)
    h,w = size(im)
    pt = topoint2d(ift)
    boundr = radius < boundary ? boundary : radius
    checkboun = (0 <= pt.x-radius-boundr && pt.x+radius+boundr <= w && 0 < pt.y-radius-boundr && pt.y+radius+boundr <= h)
    # rad = T ? radius : 0
    if checkboun
      draw!(im, CirclePointRadius(pt,radius))
    end
end

"""
    drawfeatureX!(image,feature [,crosslength = 2])

Draw a + on a feature.

# Examples
```julia-repl
julia> map(ft->drawfeatureX!(image, ft, length=5),features)
```
"""
function drawfeatureX!(im, pt::Keypoint; length = 2)
    h,w = size(im)
    boundr = 2
    checkboun = (length < pt[2] < w-length) && (length < pt[1] < h-length)
    # rad = T ? radius : 0
    if checkboun
      draw!(im, LineSegment( pt- CartesianIndex(0,length), pt + CartesianIndex(0,length)))
      draw!(im, LineSegment( pt- CartesianIndex(length,0), pt + CartesianIndex(length,0)))
    end
end

drawfeatureX!(im, ift::Feature; length = 2) = drawfeatureX!(im, Keypoint(ift), length = length)

"""
    drawfeatureLine!(image,feature1, feature2)

Draw a line between 2 features.
"""
function drawfeatureLine!(im, ift1::Feature, ift2::Feature)
    h,w = size(im)
    pt1 = Keypoint(ift1)
    pt2 = Keypoint(ift2)

    checkboun = (0 < pt1[2] <= w) && (0 < pt1[1] <= h) && (0 < pt2[2] <= w) && (0 < pt2[1] <= h)
    # rad = T ? radius : 0
    if checkboun
      draw!(im, LineSegment( pt1, pt2 ))
    end
end

function getharriscorners(iml;fracmax=0.5, k=0.04)
  harris_response = harris(iml,k=k)
  mth = fracmax*Base.maximum(harris_response)
  map(i -> i > mth, harris_response);
end
# using BenchmarkTools
# @btime mth = 0.9*maximum(harris_response)
# @btime mth = 100*Base.median(harris_response)

function getapproxbestharris(iml,nfeatures=100;k=0.04, stepguess=0.4)
    harris_response = harris(iml, k=k);
    # mth = 100*Base.mean(harris_response)
    mth = maximum(harris_response)
    # mnth = minimum(harris_response)
    # get top 100 features
    nfea = 0
    targnfea = nfeatures
    ftsl = nothing
    while nfea < targnfea
        mth *= stepguess
        resp = map(i -> i > mth, harris_response);
        ftsl = Features(resp)
        nfea = length(ftsl)
    end
    return ftsl
end

"""
    getApproxBestHarrisInWindow(image,[n=100, windowSize = 9, k=0.04, stepguess=0.4])

Return the n best Harris features in a window
"""
function getApproxBestHarrisInWindow(iml;nfeatures=100, window = 9, k=0.04, stepguess=0.9,threshold = 0.0)
    harris_response = harris(iml)

    maxima = mapwindow(maximum,harris_response,(window,window))
    window_max = (harris_response .== maxima).*harris_response

    mth = maximum(window_max)

    nfea = 0
    targnfea = nfeatures
    ftsl = nothing
    while nfea < targnfea
        mth *= stepguess
        resp = map(i -> i > mth, window_max);
        ftsl = Features(resp)
        nfea = length(ftsl)
        if mth < threshold break end
    end
    return ftsl[1: min(nfeatures, nfea)]
end

"""
    getApproxBestShiTomasi(image,[n=100, windowSize = 9, k=0.04, stepguess=0.4, threshold = 1e-4])

Return the n aproxamate best Shi Tomasi features in a window
"""
function getApproxBestShiTomasi(iml; nfeatures=100, window = 9, k=0.04, stepguess=0.7, threshold = 1e-4)
    shi_tomasi_response = shi_tomasi(iml)

    # maxima = mapwindow(maximum,shi_tomasi_response,(window,window))
    #the mapwindow extrema function is faster, but returns a tuple, the map slows it down again
    maxima = map(tu -> tu[2], mapwindow(extrema,shi_tomasi_response,(window,window)))

    shi_tomasi_response .*= (shi_tomasi_response .== maxima)# .& (shi_tomasi_response .> threshold)

    mth = maximum(shi_tomasi_response)

    nfea = 0
    targnfea = nfeatures
    ftsl = nothing
    resp = zeros(Bool,size(iml))
    while nfea < targnfea
        mth *= stepguess
        resp .= shi_tomasi_response .> mth  #map(i -> i > mth, shi_tomasi_response);
        # ftsl = Keypoints(resp)
        nfea = length(find(resp))
        if nfea >= targnfea
            # ftsl = Keypoints(resp)
            return Keypoints(resp)[1:nfeatures]
            break
        end
        if mth < threshold
            # ftsl = Keypoints(resp)
            return Keypoints(resp)[1:nfea]
            break
        end
    end
    # return ftsl[1: min(nfeatures, nfea)]
end

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
struct CPolygon
  vertexlist::Ptr{Cint}
  numberofvertices::Cint
end

struct CFacet{T}
  polygonlist::Ptr{CPolygon}
  numberofpolygons::Cint
  holelist::Ptr{T}
  numberofholes::Cint
end

struct CPPTetGenIO{T}
    firstnumber::Cint # 0 or 1, default 0.
    mesh_dim::Cint # must be 3.

    pointlist::Ptr{T}
    pointattributelist::Ptr{T}
    pointmtrlist::Ptr{T}

    pointmarkerlist::Ptr{Cint}
    numberofpoints::Cint
    numberofpointattributes::Cint
    numberofpointmtrs::Cint

    tetrahedronlist::Ptr{Cint}
    tetrahedronattributelist::Ptr{T}
    tetrahedronvolumelist::Ptr{T}
    neighborlist::Ptr{Cint}
    numberoftetrahedra::Cint
    numberofcorners::Cint
    numberoftetrahedronattributes::Cint

    facetlist::Ptr{CFacet{T}}
    facetmarkerlist::Ptr{Cint}
    numberoffacets::Cint

    holelist::Ptr{T}
    numberofholes::Cint

    regionlist::Ptr{T}
    numberofregions::Cint

    facetconstraintlist::Ptr{T}
    numberoffacetconstraints::Cint

    segmentconstraintlist::Ptr{T}
    numberofsegmentconstraints::Cint

    trifacelist::Ptr{Cint}
    trifacemarkerlist::Ptr{Cint}
    numberoftrifaces::Cint

    edgelist::Ptr{Cint}
    edgemarkerlist::Ptr{Cint}
    numberofedges::Cint
end

"""
   Error struct for TetGen
"""
struct TetGenError <: Exception
    rc::Cint
end

"""
   Show TetGen error, messages have been lifted
   from TetGen
"""
function Base.show(io::IO, e::TetGenError)
    if e.rc==1
        println(io,"TetGen error $(e.rc): out of memory.");
    elseif e.rc==2
        println(io,"TetGen error $(e.rc): internal error.")
    elseif e.rc==3
        println(io,"TetGen error $(e.rc): a self-intersection was detected. Hint: use -d option to detect all self-intersections.");
    elseif e.rc==4
        println(io,"TetGen error $(e.rc): a very small input feature size was detected. Hint: use -T option to set a smaller tolerance.")
    elseif e.rc==5
        println(io,"TetGen error $(e.rc): two very close input facets were detected. Hint: use -Y option to avoid adding Steiner points in boundary.\n");
    elseif e.rc==10
        println(io,"TetGen error $(e.rc): an input error was detected.\n");
    else
        println(io,"TetGen error $(e.rc): unknown error.\n");
    end
end

"""
$(SIGNATURES)

Tetrahedralization with error handling
"""
function tetrahedralize(input::CPPTetGenIO{Float64}, command::String)
    rc=Cint[0]
    output = ccall((:tetrahedralize2_f64, libtet), CPPTetGenIO{Float64}, (CPPTetGenIO{Float64}, Cstring,  Ptr{Cint}), input, command, rc)
    if rc[1]!=0
        throw(TetGenError(rc[1]))
    end
    output
end

"""
   Trivial Julia tetunsuitable function
"""
my_jl_tetunsuitable=(pa,pb,pc,pd)->0

"""
   Tetunsuitable function called from C wrapper
"""
function jl_wrap_tetunsuitable(pa::Ptr{Float64}, pb::Ptr{Float64}, pc::Ptr{Float64}, pd::Ptr{Float64})
    pax=Base.unsafe_wrap(Array,pa,(3,),own=false)
    pbx=Base.unsafe_wrap(Array,pb,(3,),own=false)
    pcx=Base.unsafe_wrap(Array,pc,(3,),own=false)
    pdx=Base.unsafe_wrap(Array,pd,(3,),own=false)
    Cint(my_jl_tetunsuitable(pax,pbx,pcx,pdx))
end

"""
   Set tetunsuitable function called from C wrapper.
   Setting this function is valid only for one subsequent
   call to tetrahedralize
"""
function tetunsuitable(unsuitable::Function;check_signature=true)
    if check_signature
        unsuitable(rand(3),rand(3),rand(3),rand(3))
    end
    global my_jl_tetunsuitable
    my_jl_tetunsuitable=unsuitable
    c_wrap_tetunsuitable=@cfunction(jl_wrap_tetunsuitable, Cint, (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))
    ccall((:tetunsuitable_callback,libtet),Cvoid,(Ptr{Cvoid},),c_wrap_tetunsuitable)
end


text:
@testset "Package Base.in" begin
    p = CompatHelper.Package("Foobar", UUID(1))
    ce = CompatHelper.DepInfo(p)

    @testset "Exists" begin
        s = Set([ce])

        @test p in s
    end

    @testset "DNE" begin
        p2 = CompatHelper.Package("BizBaz", UUID(0))
        @test !(p2 in Set([ce]))
    end
end


text:
using Compat  # for startswith

function test_file_string(s)
    if !startswith("test_", s)
        s = string("test_", s)
    end

    if !endswith(s, ".jl")
        s = string(s, ".jl")
    end
    return Pkg.dir("QuantEcon", "test", s)
end

if length(ARGS) > 0
    tests = map(test_file_string, ARGS)
else
    tests = map(test_file_string, ["arma", "compute_fp", "discrete_rv", "ecdf",
                                   "estspec", "kalman", "lae", "lqcontrol",
                                   "lqnash", "lss", "markov_approx",
                                   "matrix_eqn", "mc_tools", "quad", "quadsum",
                                   "random_mc", "robustlq"])
end

# n = min(8, CPU_CORES, length(tests))
# n > 1 && addprocs(n)

# @everywhere using FactCheck
# @everywhere include(Pkg.dir("QuantEcon", "test", "util.jl"))
using FactCheck
include(Pkg.dir("QuantEcon", "test", "util.jl"))

for test_file in tests
    println("running for $test_file")
    include(test_file)
end

exitstatus()


text:
using Bridge, StaticArrays, Bridge.Models
using Test, LinearAlgebra

iterations = 5000
rho = 0.05 #
independent = false # true independent proposals
t = 1.0
T = 1.75
n = 5001 # total imputed length
m = 10 # number of segments
M = div(n-1,m)
skiplast = 0
skippoints = 2
dt = (T-t)/(n-1)
tt = t:dt:T

si = 3.
P = Bridge.Models.Lorenz(ℝ{3}(10, 20, 8/3), ℝ{3}(si,si,si))

x0 = Models.x0(P)

W = sample(tt, Wiener{ℝ{3}}())
X = SamplePath(tt, zeros(ℝ{3}, length(tt)))
Bridge.solve!(Euler(), X, x0, W, P)

Y = Bridge.solve!(Bridge.R3(), copy(X), x0, P)

# Observation scheme and subsample
_pairs(collection) = Base.Generator(=>, keys(collection), values(collection))

L = I
Σ = SDiagonal(1., 1., 1.)
lΣ = cholupper(Σ)'
RV = ℝ{3}

V = SamplePath(collect(_pairs(X))[1:M:end])
#Vo = copy(V)
map!(y -> L*y + lΣ*randn(RV), V.yy, V.yy)

XX = Vector{typeof(X)}(undef, m)
XXmean = Vector{typeof(X)}(undef, m)
XXᵒ = Vector{typeof(X)}(undef, m)
WW = Vector{typeof(W)}(undef, m)
WWᵒ = Vector{typeof(W)}(undef, m)

Xtrue = copy(X)

# Create linear noise approximations

forward = false # direction of the linear noise approximation

TPt = Any
TPᵒ = Any

𝕃 = typeof(Bridge.outer(zero(x0)))

Phi = Vector{Any}(undef, m)
Pt = Vector{TPt}(undef, m)
Pᵒ = Vector{TPᵒ}(undef, m)
H♢ = Bridge.outer(zero(x0))
#v = Xtrue.yy[end]
v = ( L' * inv(Σ) * L)\(L' * inv(Σ) *  V.yy[end])
H♢ = one(Bridge.outer(zero(x0)))*inv( L' * inv(Σ) * L)

for i in m:-1:1
    XX[i] = SamplePath(X.tt[1 + (i-1)*M:1 + i*M], X.yy[1 + (i-1)*M:1 + i*M])
    WW[i] = SamplePath(W.tt[1 + (i-1)*M:1 + i*M], W.yy[1 + (i-1)*M:1 + i*M])
    # short-cut, take v later
    Pt[i] =  Bridge.LinearAppr(Y[1 + (i-1)*M:1 + i*M], P)
    Phi[i] = Bridge.solvei(Bridge.R3(), (i,C,iP) -> Bridge.Bi(iP[1],iP[2])*C, XX[i].tt, one(𝕃), Pt[i])
    #Bridge.LinearNoiseAppr(XX[i].tt, P, XX[i].yy[end], Bridge.a(XX[i].tt[end], XX[i].yy[end], P), forward)
    #Pt[i] = Bridge.LinearNoiseAppr(XX[i].tt, P, v, Bridge.a(XX[i].tt[end], v, P), forward)
    Pᵒ[i] = Bridge.GuidedBridge(XX[i].tt, P, Pt[i], v, H♢)
    H♢, v = Bridge.gpupdate(Pᵒ[i], L, Σ, V.yy[i])
end

v0 = V.yy[1]
#v0 = v

y = x0
for i in 1:m
    sample!(WW[i], Wiener{ℝ{3}}())
    y = Bridge.solve!(Euler(), XX[i], y, WW[i], Pᵒ[i])
end
XXmean = [zero(XX[i]) for i in 1:m]

μ0 = ℝ{3}(0,0,0)
Σ0 = SDiagonal(40., 40., 40.)
π0 = Bridge.Gaussian(μ0 + Σ0*L'*inv(L*Σ0*L' + Σ)*(v0 - L*μ0), Σ0 - Σ0*L'*inv(L*Σ0*L' + Σ)*L*Σ0)

H♢, v = Bridge.gpupdate(Pᵒ[1], L, Σ0, V.yy[1])
π0 = Bridge.Gaussian(v, H♢)

if false
    using GLVisualize, Colors, GeometryTypes
    window = glscreen()
    sca = 1/15
    phi = Float32(0.2*2pi)

    rotation = @SMatrix Float32[ sca*cos(phi)  -sca*sin(phi)    0.0       0.0
    sca*sin(phi)   sca*cos(phi)    0.0       0.0
    0.0            0.0             sca*1.0  -1.25
    0.0            0.0             0.0       1.0]

    X3d = visualize(
        Point3f0.(X.yy), :lines,
        color = RGBA{Float32}(0.04, 0.15, 0.44, 0.6),
        model = rotation,
    )

    XX3d = visualize(
        Point3f0.(vcat(XX...).yy), :lines,
        color = RGBA{Float32}(0.04, 0.44, 0.12, 0.6),
        model = rotation,
    )

    Y3d = visualize(
        Point3f0.(Y.yy), :lines,
        color = RGBA{Float32}(0.54, 0.15, 0.14, 0.6),
        model = rotation,
    )
    circle = Sphere(Point2f0(0), 1f0)
    V3d = visualize(
        (circle, Point3f0.(V.yy)),
        scale = fill(sca*0.5Point2f0(1,1), length(V)),
        color = RGBA{Float32}(0.7, 0.3, 0., 0.8),
        model = rotation,
    )

    _view(X3d, window, camera = :perspective)
    _view(XX3d, window, camera = :perspective)
    _view(Y3d, window, camera = :perspective)
    _view(V3d, window, camera = :perspective)
    renderloop(window)
end
h = sqrt(eps())
@test norm(SMatrix{3,3}(Iterators.flatten(h\(Bridge.b(0, x0 + h*eye(SMatrix{3,3,Float64})[i,:], P) - Bridge.b(0, x0, P)) for i in 1:3)...) - Bridge.bderiv(0, x0, P)) < 100h


text:
"""
    write_smac2_par([...])

Writes a P-Smac2 parameter file. Not all relevant parameters implemented yet!
"""
function write_smac2_par(x, y, z,
                         euler_angle_0, euler_angle_1, euler_angle_2,
                         xy_size, z_depth, xy_pix::Integer,
                         input_file, output_file, path,
                         effect_module::Integer=0, effect_flag::Integer=0,
                         ν_obs::Real=1.44e9, cosmology::Integer=0,
                         CR_pmin::Real=10.0, CR_pmax::Real=1.e7)

    open(path * "smac2.par", "w") do f
        write(f,
        """%% PSmac Parameter File %%

        Cosmology $cosmology         % set to 0 for non-cosmological sims.  Cosmo.h = 1

        %% Image Properties, comoving %%

        Center_X $x
        Center_Y $y
        Center_Z $z

        Use_Barycenter 0    % 1 - set to, 2 - set relative to barycenter

        XYSize $xy_size        % set these two ==0 to project whole box
        ZDepth $z_depth
        XYPix $xy_pix

        %% Projection, see description at EOF %%
        Euler_Angle_0 $euler_angle_0    % [deg] along x : (0,90,90)
        Euler_Angle_1 $euler_angle_1    % [deg] along y : (0,90,0)
        Euler_Angle_2 $euler_angle_2     % [deg] along z : (0,0,0)

        %% HEALPIX %%

        NSide 32
        Rmin 1000
        Rmax 1000000

        %% I/O Options %%

        Input_File $input_file
        N_IOTasks 1

        NoClobber 0
        Output_File $output_file

        %% Effect Selection & Options %%

        Effect_Module $effect_module
        Effect_Flag $effect_flag

        E_min 5e4             % Energy Range [eV]: Xray, Gamma, Synchro
        E_max 1e10

        Freq $ν_obs          % Frequency [Hz]: Sz,Synchro

        a_cr 2.375          % CR spectral index
        X_cr 0.01           % CR Normalisation rel. thermal
        IntrRM 0            % Toggle Intrinsic RM
        PitchAngInt 1       % Toggle Pitch angle integration

        t_turb 1e7          % Reacceleration Timescale Cassano 05
        eta_t 0.26          % Fraction turb. Energy in magnetosonic waves
        turb_scale 100      % scale to extrapolate turbulent velocity to

        CR_Emin $CR_pmin       % Minimum energy for BP_REAL_CRs, as in OpenGadget3 parameterfile
        CR_Emax $CR_pmax       % Maximum energy for BP_REAL_CRs, as in OpenGadget3 parameterfile

        %Gadget
        UnitLength_in_cm 			3.085678e21        %  1.0 kpc
        UnitMass_in_g 				1.989e43           %  1.0e10 solar masses
        UnitVelocity_in_cm_per_s 	1e5                %  1 km/sec
        """
        )
    end
end

"""
    read_smac2_info(filename::String)

Returns the image info in a `mappingParameters` struct.
"""
function read_smac2_info(filename::String)

    f = FITS(filename)

    center = [read_key(f[1], 23)[2],
              read_key(f[1], 24)[2],
              read_key(f[1], 25)[2] ]

    xy_size = read_key(f[1], 27)[2]
    z_size  = read_key(f[1], 28)[2]
    Npixels = read_key(f[1], 29)[2]

    close(f)

    return mappingParameters(center=center,
									x_size=xy_size,
									y_size=xy_size,
									z_size=z_size,
									Npixels=Npixels)

end

"""
    read_smac2_image(filename::String)

Returns the image of a Smac2 FITS file.
"""
function read_smac2_image(filename::String, num_image::Integer=1)

    f = FITS(filename)

    image = read(f[1])[:,:,num_image]

    close(f)

    return image

end


text:
# extension of NFFT to parallel workers

using NFFT

function pnfft_combi_worker(myz::Vector{Int}, data::Array{Complex64,3},
    K::Array{Float32,3}, pstart::Vector{Int},
    pend::Vector{Int}, kernwidth=3, oversamp=1.5f0; verbose=false)
  nchan, nro, npe = size(data)
  nz = length(pstart)
  n = int(nro / 2)
  fApprox = Array(Complex64, n, n, nz)
  sliceImgs = Array(Complex64, nchan, n, n)
  for z in myz   # for each "slice", recon 2D image
    if verbose
      println("slice $z/$nz, $nchan channels")
    end

    # initialize NFFT
    Kslice = K[:,:,pstart[z]:pend[z]]
    nfft_plan = NFFTPlan(Kslice[:,:], (n, n), kernwidth, oversamp)

    # SDC
    mynpe = pend[z] - pstart[z] + 1
    weights = radialsdc(nro, mynpe; golden=true, offset=pstart[z])
    #weights = sdc(nfft_plan) # Pipe & Menon SDC, not necessary for golden angle radial

    # apply NFFT plan to each channel
    for c in 1:size(data,1)
      # pull out the relevant profiles
      #fhat = vec(data[c,:,pstart[z]:pend[z]])
      mydata = data[c,:,pstart[z]:pend[z]]
      # radial phase correction
      fhat, slope, intercepts = unshift(mydata; golden=true)
      # density compensate
      fhat = fhat[:] .* weights[:]
      # transform and save the resulting image
      sliceImgs[c,:,:] = nfft_adjoint(nfft_plan, fhat)
    end
    # optimal channel combination
    fApprox[:,:,z] = walshcombine(sliceImgs)
  end
  fApprox
end

function pnfft_combi(data::Array{Complex64,3}, K::Array{Float32,3},
                     pstart::Vector{Int},
                     pend::Vector{Int}, kernwidth=3, oversamp=1.5f0;
                     verbose=false)

  # set up parallel stuff
  nz = length(pstart)
  nchan, nro, npe = size(data)
  n = int(nro / 2)
  nw = nworkers()
  zs = [1:nz]
  zperw = int(nz/nw)
  myz = Any[zs[(w - 1)*zperw + 1 : min(w*zperw, nz)] for w in  1:nw]
  println("processing $nz slices on $nw workers ($zperw slices each)")
  if nw == 1
    fApprox = pnfft_combi_worker(myz[1], data, K, pstart, pend,
                          kernwidth, oversamp; verbose=verbose)
  else
    reflist = Any[]
    wids = workers()
    for w in 1:nw
      println("spawning job on worker $(wids[w])")
      r = @spawnat wids[w] pnfft_combi_worker(myz[w], data, K,
             pstart, pend, kernwidth, oversamp; verbose=verbose)
      push!(reflist, r)
    end
    fApprox = Array(Complex64, n, n, nz)
    for w in 1:nw
      r = fetch(reflist[w])
      println("fetching zs $(myz[w])")
      fApprox[:,:,myz[w]] = r[:,:,myz[w]]
    end
  end
  fApprox
end

max_stars_repo_path: code/Julia Recon Scripts/pnfft.jl

max_issues_repo_path: code/Julia Recon Scripts/pnfft.jl

max_forks_repo_path: code/Julia Recon Scripts/pnfft.jl


text:
module ClimateMaps

using PyCall
using Reexport
@reexport using ClimateTypes

end # module


text:
include("autotuning_utils.jl")
function update_hyperparameters!(model::Union{GP,VGP})
    update_hyperparameters!.(model.f,get_Z(model))
    model.inference.HyperParametersUpdated = true
end

function update_hyperparameters!(model::SVGP)
    update_hyperparameters!.(model.f,[model.inference.xview],∇E_μ(model.likelihood,model.inference.vi_opt[1],get_y(model)),∇E_Σ(model.likelihood,model.inference.vi_opt[1],get_y(model)),model.inference,model.inference.vi_opt)
    model.inference.HyperParametersUpdated = true
end

function update_hyperparameters!(model::MOSVGP)
    update_hyperparameters!.(model.f,[model.inference.xview],∇E_μ(model),∇E_Σ(model),model.inference,model.inference.vi_opt)
    model.inference.HyperParametersUpdated = true
end

## Update all hyperparameters for the full batch GP models ##
function update_hyperparameters!(gp::Union{_GP{T},_VGP{T}},X::AbstractMatrix) where {T}
    if !isnothing(gp.opt_ρ) || !isnothing(gp.opt_σ) || !isnothing(get_opt(gp.μ₀))
        f_l,f_v,f_μ₀ = hyperparameter_gradient_function(gp,X)
        if !isnothing(gp.opt_ρ)
            Jnn = kernelderivative(gp.kernel,X)
            grads_ρ = compute_hyperparameter_gradient(gp.kernel,f_l,Jnn)
        end
        if !isnothing(gp.opt_σ)
            grads_σ_k = f_v(gp.kernel,gp.σ_k)
        end
        if !isnothing(get_opt(gp.μ₀))
            grads_μ₀ = f_μ₀()
        end

        if !isnothing(gp.opt_ρ)
            apply_gradients_lengthscale!(gp.opt_ρ,gp.kernel,grads_ρ) #Send the derivative of the matrix to the specific gradient of the model
        end
        if !isnothing(gp.opt_σ)
            apply_gradients_variance!(gp,grads_σ_k) #Send the derivative of the matrix to the specific gradient of the model
        end
        if !isnothing(get_opt(gp.μ₀))
            apply_gradients_mean_prior!(gp.μ₀,grads_μ₀,X)
        end
    end
end

## Update all hyperparameters for the sparse variational GP models ##
function update_hyperparameters!(gp::_SVGP{T},X,∇E_μ::AbstractVector{T},∇E_Σ::AbstractVector{T},i::Inference,opt::AbstractOptimizer) where {T}
    if !isnothing(gp.opt_ρ) || !isnothing(gp.opt_σ) || !isnothing(get_opt(gp.μ₀))
        f_ρ,f_σ_k,f_μ₀ = hyperparameter_gradient_function(gp)
        if !isnothing(gp.opt_ρ)
            Jmm = kernelderivative(gp.kernel,gp.Z.Z)
            Jnm = kernelderivative(gp.kernel,X,gp.Z.Z)
            Jnn = kerneldiagderivative(gp.kernel,X)
            grads_l = compute_hyperparameter_gradient(gp.kernel,f_ρ,Jmm,Jnm,Jnn,∇E_μ,∇E_Σ,i,opt)
        end
        if !isnothing(gp.opt_σ)
            grads_σ_k = f_σ_k(gp.kernel,gp.σ_k,∇E_Σ,i,opt)
        end
        if !isnothing(get_opt(gp.μ₀))
            grads_μ₀ = f_μ₀()
        end
        if !isnothing(gp.Z.opt)
            Z_gradients = inducingpoints_gradient(gp,X,∇E_μ,∇E_Σ,i,opt) #Compute the gradient given the inducing points location
            gp.Z.Z .+= GradDescent.update(gp.Z.opt,Z_gradients) #Apply the gradients on the location
        end
        if !isnothing(gp.opt_ρ)
            apply_gradients_lengthscale!(gp.opt_ρ,gp.kernel,grads_l)
        end
        if !isnothing(gp.opt_σ)
            apply_gradients_variance!(gp,grads_σ_k)
        end
        if !isnothing(get_opt(gp.μ₀))
            apply_gradients_mean_prior!(gp.μ₀,grads_μ₀,gp.Z.Z)
        end
    end
end

## Return the derivative of the KL divergence between the posterior and the GP prior ##
function hyperparameter_KL_gradient(J::AbstractMatrix{T},A::AbstractMatrix{T}) where {T<:Real}
    return 0.5*opt_trace(J,A)
end

function hyperparameter_gradient_function(gp::_GP{T},X::AbstractMatrix) where {T}
    A = (inv(gp.K)+ (-gp.μ*transpose(gp.μ)))
    return (function(Jnn)
                return -hyperparameter_KL_gradient(Jnn,A)
            end,
            function(kernel,σ_k)
                return -one(T)/σ_k*hyperparameter_KL_gradient(gp.K.mat,A)
            end,
            function()
                return -gp.K\(gp.μ₀(X)-gp.y)
            end)
end

## Return functions computing gradients of the ELBO given the kernel hyperparameters for a non-sparse model ##
function hyperparameter_gradient_function(gp::_VGP{T},X::AbstractMatrix) where {T<:Real}
    A = (I-gp.K\(gp.Σ+(gp.µ-gp.μ₀(X))*transpose(gp.μ-gp.μ₀(X))))/gp.K.mat
    return (function(Jnn)
                return -hyperparameter_KL_gradient(Jnn,A)
            end,
            function(kernel,σ_k)
                return -one(T)/σ_k*hyperparameter_KL_gradient(gp.K.mat,A)
            end,
            function()
                return -gp.K\(gp.μ₀(X)-gp.μ)
            end)
end

## Return functions computing gradients of the ELBO given the kernel hyperparameters for a non-sparse latent GP ##
function hyperparameter_gradient_function(gp::_SVGP{T}) where {T<:Real}
    A = (Diagonal{T}(I,gp.dim).-gp.K\(gp.Σ.+(gp.µ-gp.μ₀(gp.Z.Z))*transpose(gp.μ-gp.μ₀(gp.Z.Z))))/gp.K.mat
    ι = similar(gp.κ) #Empty container to save data allocation
    κΣ = gp.κ*gp.Σ
    return (function(Jmm,Jnm,Jnn,∇E_μ,∇E_Σ,i,opt)
                return (hyperparameter_expec_gradient(gp,∇E_μ,∇E_Σ,i,opt,ι,κΣ,Jmm,Jnm,Jnn)-hyperparameter_KL_gradient(Jmm,A))
            end,
            function(kernel::Kernel,σ_k::Real,∇E_Σ,i,opt)
                return one(T)/σ_k*(
                        - i.ρ*dot(∇E_Σ,gp.K̃)
                        - hyperparameter_KL_gradient(gp.K.mat,A))
            end,
            function()
                return -gp.K\(gp.μ₀(gp.Z.Z)-gp.μ)
            end)
end

function hyperparameter_gradient_function(model::VStP{T},X::AbstractMatrix) where {T<:Real}
    A = ([Diagonal{T}(I,model.nFeatures)].-invK(model).*(model.Σ.+(model.µ.-model.μ₀(X)).*transpose.(model.μ.-model.μ₀(X)))).*invK(model)
    return (function(Jnn,index)
                return -hyperparameter_KL_gradient(Jnn,A[index])
            end,
            function(kernel,index)
                return -1.0/getvariance(kernel)*hyperparameter_KL_gradient(model.Knn[index].*model.χ[index],A[index])
            end,
            function(index)
                return -invK(model,index)*(model.μ₀(X)-model.μ[index])
            end)
end

## Gradient with respect to hyperparameter for analytical VI ##
function hyperparameter_expec_gradient(gp::_SVGP{T},∇E_μ::AbstractVector{T},∇E_Σ::AbstractVector{T},i::AnalyticVI,opt::AVIOptimizer,ι::AbstractMatrix{T},κΣ::AbstractMatrix{T},Jmm::AbstractMatrix{T},Jnm::AbstractMatrix{T},Jnn::AbstractVector{T}) where {T<:Real}
    ι .= (Jnm-gp.κ*Jmm)/gp.K.mat
    Jnn .-= opt_diag(ι,gp.Knm) + opt_diag(gp.κ,Jnm)
    dμ = dot(∇E_μ,ι*gp.μ)
    dΣ = -dot(∇E_Σ,Jnn)
    dΣ += -dot(∇E_Σ,2.0*(opt_diag(ι,κΣ)))
    dΣ += -dot(∇E_Σ,2.0*(ι*gp.μ).*(gp.κ*gp.μ))
    return i.ρ*(dμ+dΣ)
end

## Gradient with respect to hyperparameters for numerical VI ##
function hyperparameter_expec_gradient(gp::_SVGP{T},∇E_μ::AbstractVector{T},∇E_Σ::AbstractVector{T},i::NumericalVI,opt::NVIOptimizer,ι::AbstractMatrix{T},κΣ::AbstractMatrix{T},Jmm::AbstractMatrix{T},Jnm::AbstractMatrix{T},Jnn::AbstractVector{T}) where {T<:Real}
    ι .= (Jnm-gp.κ*Jmm)/gp.K.mat
    Jnn .-= opt_diag(ι,gp.Knm) + opt_diag(gp.κ,Jnm)
    dμ = dot(∇E_μ,ι*gp.μ)
    dΣ = dot(∇E_Σ,Jnn+2.0*opt_diag(ι,κΣ))
    return i.ρ*(dμ+dΣ)
end

## Return a function computing the gradient of the ELBO given the inducing point locations ##
function inducingpoints_gradient(gp::_SVGP{T},X,∇E_μ::AbstractVector{T},∇E_Σ::AbstractVector{T},i::Inference,opt::AbstractOptimizer) where {T<:Real}
    gradient_inducing_points = similar(gp.Z.Z)
    A = (I-gp.K\(gp.Σ+gp.µ*transpose(gp.μ)))/gp.K
    #preallocation
    ι = similar(gp.κ)
    Jmm,Jnm = indpoint_derivative(gp.kernel,gp.Z),indpoint_derivative(gp.kernel,X,gp.Z)
    κΣ = gp.κ*gp.Σ
    for j in 1:gp.dim #Iterate over the points
        for k in 1:size(gp.Z,2) #iterate over the dimensions
            @views ι = (Jnm[:,:,j,k]-gp.κ*Jmm[:,:,j,k])/gp.K
            @views gradient_inducing_points[j,k] = hyperparameter_expec_gradient(gp,∇E_μ,∇E_Σ,i,opt,ι,κΣ,Jmm[:,:,j,k],Jnm[:,:,j,k],zero(gp.K̃))-hyperparameter_KL_gradient(Jmm[:,:,j,k],A)
        end
    end
    return gradient_inducing_points
end


text:
"""
# module FiniteDifference

# Examples

```jldoctest
julia>
```
"""
module FiniteDifference

using BandedMatrices

export forward_difference_matrix,
    backward_difference_matrix,
    cyclic_difference_matrix,
    centered_difference_matrix

"""
    forward_difference_matrix(T::Type{<:Number}, shape::NTuple{N,<:Number})

Create a forward difference matrix.

The default type of returned value `A` is a `BandedMatrix`.
If you want a regular `Matrix`, run `Matrix(A)`.
"""
function forward_difference_matrix(T::Type{<:Number}, shape::NTuple{N,<:Number}) where {N}
    A = BandedMatrix(-Ones{T}(shape), (0, 1))
    A[band(1)] .= 1
    return A
end

"""
    backward_difference_matrix(T::Type{<:Number}, shape::NTuple{N,<:Number})

Create a backward difference matrix.

The default type of returned value `A` is a `BandedMatrix`.
If you want a regular `Matrix`, run `Matrix(A)`.
"""
function backward_difference_matrix(T::Type{<:Number}, shape::NTuple{N,<:Number}) where {N}  # I am unsure about this construction.
    A = BandedMatrix(Ones{T}(shape), (1, 0))
    A[band(-1)] .= -1
    return A
end

function cyclic_difference_matrix(T::Type{<:Number}, shape::NTuple{N,<:Number}) where {N}
    A = BandedMatrix(Ones{T}(shape), (1, 0))
    A[band(-1)] .= -1
    B = Matrix(A)
    B[1, end] = -1
    return B
end

function centered_difference_matrix(T::Type{<:Number}, shape::NTuple{N,<:Number}) where {N}
    A = BandedMatrix(Zeros{T}(shape), (1, 1))
    A[band(-1)] .= -1
    A[band(1)] .= 1
    return A
end

end


text:
export CubicSpline

struct CubicSpline end

"""

    compute_interpolants( n, f)

Compute interpolation coefficients

This function is a julia version of a Fortran code written by
Edwin Chacon Golcher (Institute of Physics of the Czech Academy of Sciences)

"""
function compute_interpolants( n::Int, f::Array{Float64})

    @assert (n > 27)

    num_terms = 27

    coeffs = zeros(Float64,n+3)
    d      = zeros(Float64,n)

    a   = sqrt((2.0+sqrt(3.0))/6.0)
    r_a = 1.0/a
    b   = sqrt((2.0-sqrt(3.0))/6.0)
    b_a = b/a

    d1 = f[1]
    coeff_tmp = 1.0
    for i in 0:num_terms-1
        coeff_tmp *= (-b_a)
        d1 += coeff_tmp*f[n-1-i]
    end

    d[1] = d1*r_a
    for i in 2:n-1
        d[i] = r_a*(f[i] - b*d[i-1])
    end

    d1        = d[n-1]
    coeff_tmp = 1.0
    for i in 1:num_terms
        coeff_tmp *= (-b_a)
        d1 += coeff_tmp*d[i]
    end

    coeffs[n] = d1*r_a

    for i in n-2:-1:1
        coeffs[i+1] = r_a*(d[i] - b*coeffs[i+2])
    end

    coeffs[1]   = coeffs[n]
    coeffs[n+1] = coeffs[2]
    coeffs[n+2] = coeffs[3]
    coeffs[n+3] = coeffs[4]
    coeffs
end

export interpolate

"""

    interpolate( coeffs, n, f)

Compute interpolatted value at `x` from coefficients get
from `compute_interpolants`.

This function is a julia version of a Fortran code written by

Edwin Chacon Golcher (Institute of Physics of the Czech Academy of Sciences)

"""
function interpolate( coeffs :: Array{Float64,1},
                      n1     :: Int,
                      x1min   :: Float64,
                      x1max   :: Float64,
                      x      :: Float64 )

      rh        = (n1-1) / (x1max - x1min)
      t0        = (x-x1min)*rh
      cell      = floor(t0)
      delta1        = t0 - cell
      cdelta1       = 1.0 - delta1
      icell     = trunc(Int, cell+1)
      cim1      = coeffs[icell]
      ci        = coeffs[icell+1]
      cip1      = coeffs[icell+2]
      cip2      = coeffs[icell+3]
      t1        = 3.0*ci
      t3        = 3.0*cip1
      t2        = cdelta1*(cdelta1*(cdelta1*(cim1 - t1) + t1) + t1) + ci
      t4        =  delta1*( delta1*( delta1*(cip2 - t3) + t3) + t3) + cip1

      (1.0/6.0) * (t2 + t4)

end

"""
     advection!( f, mesh1, v,  dt)

Semi-lagrangian advection function of 2D distribution function represented
by array `f`. The advection operates along `axis` (=1 is most efficient)
with speed `v` during `dt`.

It uses cubic splines interpolation.

"""
function advection!( f::Array{Float64,2},
		     mesh::UniformMesh,
                     v::Vector{Float64},
		     dt::Float64,
		     interp::CubicSpline,
		     axis::Int64 )

    if (axis == 1)

        @assert ( mesh.length == size(f)[1] )
        @assert ( size(v)[1]  == size(f)[2] )
        lx = mesh.stop - mesh.start
        @simd for j in 1:size(v)[1]
            coeffs = compute_interpolants(mesh.length, f[:,j])
            @inbounds for i in 1:mesh.length
                x_new = mesh.points[i] - dt * v[j]
                x_new = mesh.start + mod(x_new - mesh.start, lx)
                f[i,j] = interpolate(coeffs, mesh.length, mesh.start, mesh.stop, x_new)
            end
        end

    else

        @assert ( mesh.length == size(f)[2] )
        @assert ( size(v)[1]  == size(f)[1] )

        lv = mesh.stop - mesh.start
        @simd for i in 1:size(v)[1]
            coeffs = compute_interpolants(mesh.length, f[i,:])
            @inbounds for j in 1:mesh.length
                v_new = mesh.points[j] - dt * v[i]
                v_new = mesh.start + mod(v_new - mesh.start, lv)
                f[i,j] = interpolate(coeffs, mesh.length, mesh.start, mesh.stop, v_new)
            end
        end

    end

end