text:
using Tensors
using Parameters
using NLsolve
using ForwardDiff
using Profile

abstract type AbstractMaterial end
abstract type AbstractMaterialState end

@generated function Base.:+(state::T, dstate::T) where {T <: AbstractMaterialState}
   expr = [:(state.$p+ dstate.$p) for p in fieldnames(T)]
   return :(T($(expr...)))
end

@with_kw mutable struct ChabocheDriverState <: AbstractMaterialState
    time :: Float64 = zero(Float64)
    strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
end

@with_kw struct ChabocheParameterState <: AbstractMaterialState
    E :: Float64 = 0.0
    nu :: Float64 = 0.0
    R0 :: Float64 = 0.0
    Kn :: Float64 = 0.0
    nn :: Float64 = 0.0
    C1 :: Float64 = 0.0
    D1 :: Float64 = 0.0
    C2 :: Float64 = 0.0
    D2 :: Float64 = 0.0
    Q :: Float64 = 0.0
    b :: Float64 = 0.0
end

@with_kw struct ChabocheVariableState <: AbstractMaterialState
    stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    X1 :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    X2 :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    plastic_strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    cumeq :: Float64 = zero(Float64)
    R :: Float64 = zero(Float64)
    jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64})
end

@with_kw mutable struct Chaboche <: AbstractMaterial
    drivers :: ChabocheDriverState = ChabocheDriverState()
    ddrivers :: ChabocheDriverState = ChabocheDriverState()
    variables :: ChabocheVariableState = ChabocheVariableState()
    variables_new :: ChabocheVariableState = ChabocheVariableState()
    parameters :: ChabocheParameterState = ChabocheParameterState()
    dparameters :: ChabocheParameterState = ChabocheParameterState()
end

@with_kw mutable struct IdealPlasticDriverState <: AbstractMaterialState
    time :: Float64 = zero(Float64)
    strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
end

@with_kw struct IdealPlasticParameterState <: AbstractMaterialState
    youngs_modulus :: Float64 = zero(Float64)
    poissons_ratio :: Float64 = zero(Float64)
    yield_stress :: Float64 = zero(Float64)
end

@with_kw struct IdealPlasticVariableState <: AbstractMaterialState
    stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    plastic_strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    cumeq :: Float64 = zero(Float64)
end

function update!(material::M) where {M <: AbstractMaterial}
    material.drivers += material.ddrivers
    material.parameters += material.dparameters
    material.variables = material.variables_new
    # material.ddrivers = typeof(material.ddrivers)()
    # material.dparameters = typeof(material.dparameters)()
    # material.variables_new = typeof(material.variables_new)()
    reset!(material)
end

function reset!(material::M) where {M <: AbstractMaterial}
    material.ddrivers = typeof(material.ddrivers)()
    material.dparameters = typeof(material.dparameters)()
    material.variables_new = typeof(material.variables_new)()
end

# @with_kw mutable struct Material{M <: AbstractMaterial}
#     drivers :: DriverState{M} = DriverState{M}()
#     ddrivers :: DriverState{M} = DriverState{M}()
#     variables :: VariableState{M} = VariableState{M}()
#     dvariables :: VariableState{M} = VariableState{M}()
#     parameters :: ParameterState{M} = ParameterState{M}()
#     dparameters :: ParameterState{M} = ParameterState{M}()
# end

@with_kw mutable struct IdealPlastic <: AbstractMaterial
    drivers :: IdealPlasticDriverState = IdealPlasticDriverState()
    ddrivers :: IdealPlasticDriverState = IdealPlasticDriverState()
    variables :: IdealPlasticVariableState = IdealPlasticVariableState()
    variables_new :: IdealPlasticVariableState = IdealPlasticVariableState()
    parameters :: IdealPlasticParameterState = IdealPlasticParameterState()
    dparameters :: IdealPlasticParameterState = IdealPlasticParameterState()
end

mat = Chaboche()
mat2 = IdealPlastic()

function isotropic_elasticity_tensor(lambda, mu)
    delta(i,j) = i==j ? 1.0 : 0.0
    g(i,j,k,l) = lambda*delta(i,j)*delta(k,l) + mu*(delta(i,k)*delta(j,l)+delta(i,l)*delta(j,k))
    jacobian = SymmetricTensor{4, 3, Float64}(g)
    return jacobian
end

function integrate_material!(material::Chaboche)
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = p
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R, jacobian = v

    jacobian = isotropic_elasticity_tensor(lambda, mu)

    stress += dcontract(jacobian, dstrain)
    seff = stress - X1 - X2
    seff_dev = dev(seff)
    f = sqrt(1.5)*norm(seff_dev) - (R0 + R)
    if f > 0.0
        g! = create_nonlinear_system_of_equations(material)
        x0 = [tovoigt(stress); R; tovoigt(X1); tovoigt(X2)]
        F = similar(x0)
        res = nlsolve(g!, x0; autodiff = :forward)
        x = res.zero
        res.f_converged || error("Nonlinear system of equations did not converge!")

        stress = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[1:6])
        R = x[7]
        X1 = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[8:13])
        X2 = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[14:19])
        seff = stress - X1 - X2
        seff_dev = dev(seff)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R)
        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = 1.5*seff_dev/norm(seff_dev)

        plastic_strain += dp*n
        cumeq += dp
        # Compute Jacobian
        function residuals(x)
            F = similar(x)
            g!(F, x)
            return F
        end
        drdx = ForwardDiff.jacobian(residuals, x)
        drde = zeros((length(x),6))
        drde[1:6, 1:6] = -tovoigt(jacobian)
        jacobian = fromvoigt(SymmetricTensor{4,3}, (drdx\drde)[1:6, 1:6])
    end
    variables_new = ChabocheVariableState(stress = stress,
                                          X1 = X1,
                                          X2 = X2,
                                          R = R,
                                          plastic_strain = plastic_strain,
                                          cumeq = cumeq,
                                          jacobian = jacobian)
    material.variables_new = variables_new
end

function create_nonlinear_system_of_equations(material::Chaboche)
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = p
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R = v

    function g!(F, x::Vector{T}) where {T} # System of non-linear equations
        jacobian = isotropic_elasticity_tensor(lambda, mu)
        stress_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[1:6])
        R_ = x[7]
        X1_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[8:13])
        X2_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[14:19])

        seff = stress_ - X1_ - X2_
        seff_dev = dev(seff)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R_)

        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = 1.5*seff_dev/norm(seff_dev)
        dstrain_plastic = dp*n
        tovoigt!(view(F, 1:6), stress - stress_ + dcontract(jacobian, dstrain - dstrain_plastic))
        F[7] = R - R_ + b*(Q-R_)*dp
        if isapprox(C1, 0.0)
            tovoigt!(view(F,8:13),X1 - X1_)
        else
            tovoigt!(view(F,8:13), X1 - X1_ + 2.0/3.0*C1*dp*(n - 1.5*D1/C1*X1_))
        end
        if isapprox(C2, 0.0)
            tovoigt!(view(F,14:19), X2 - X2_)
        else
            tovoigt!(view(F, 14:19), X2 - X2_ + 2.0/3.0*C2*dp*(n - 1.5*D2/C2*X2_))
        end
    end
    return g!
end

function simple_integration_test()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)

    dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
    ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
    chabmat = Chaboche(parameters = parameters, ddrivers = ddrivers)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
end

using DelimitedFiles, Test
path = joinpath(@__DIR__, "one_elem_disp_chaboche", "unitelement_results.rpt")
data = readdlm(path, Float64; skipstart=4)
ts = data[:,1]
s11_ = data[:,2]
s12_ = data[:,3]
s13_ = data[:,4]
s22_ = data[:,5]
s23_ = data[:,6]
s33_ = data[:,7]
e11_ = data[:,8]
e12_ = data[:,9]
e13_ = data[:,10]
e22_ = data[:,11]
e23_ = data[:,12]
e33_ = data[:,13]

strains = [[e11_[i], e22_[i], e33_[i], e23_[i], e13_[i], e12_[i]] for i in 1:length(ts)]

function test_chaboche()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)
    chabmat = Chaboche(parameters = parameters)
    s33s = [chabmat.variables.stress[3,3]]
    for i=2:length(ts)
        dtime = ts[i]-ts[i-1]
        dstrain = fromvoigt(SymmetricTensor{2,3,Float64}, strains[i]-strains[i-1]; offdiagscale=2.0)
        chabmat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain)
        integrate_material!(chabmat)
        update!(chabmat)
        push!(s33s, chabmat.variables.stress[3,3])
    end
    # @test isapprox(s33s, s33_; rtol=0.01)
end

test_chaboche()
Profile.clear_malloc_data()
test_chaboche()
# using BenchmarkTools
# @btime test_chaboche()

function simple_integration_test_fd_tangent()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)

    dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
    ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
    chabmat = Chaboche(parameters = parameters, ddrivers = ddrivers)

    function get_stress(dstrain)
        chabmat.ddrivers.strain = dstrain
        integrate_material!(chabmat)
        return chabmat.variables_new.stress
    end
    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D_mat = $(tovoigt(chabmat.variables_new.jacobian))"
    @info "D = $(tovoigt(D))"

    chabmat.variables_new = typeof(chabmat.variables_new)()
    chabmat.ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    function get_stress(dstrain)
        chabmat.ddrivers.strain = dstrain
        integrate_material!(chabmat)
        return chabmat.variables_new.stress
    end
    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D = $(tovoigt(D))"
    chabmat.variables_new = typeof(chabmat.variables_new)()
    chabmat.ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
    function get_stress(dstrain)
        chabmat.ddrivers.strain = dstrain
        integrate_material!(chabmat)
        return chabmat.variables_new.stress
    end
    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D = $(tovoigt(D))"
    chabmat.variables_new = typeof(chabmat.variables_new)()
    chabmat.ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
    function get_stress(dstrain)
        chabmat.ddrivers.strain = dstrain
        integrate_material!(chabmat)
        return chabmat.variables_new.stress
    end
    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D = $(tovoigt(D))"
end

# simple_integration_test_fd_tangent()

function simple_integration_test_fd_tangent2()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)

    dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
    ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
    chabmat = Chaboche(parameters = parameters, ddrivers = ddrivers)
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    g! = create_nonlinear_system_of_equations(chabmat)
    function residuals(x)
        F = similar(x)
        g!(F, x)
        return F
    end

    x0 = [tovoigt(chabmat.variables_new.stress); chabmat.variables_new.R; tovoigt(chabmat.variables_new.X1); tovoigt(chabmat.variables_new.X2)]
    drdx = ForwardDiff.jacobian(residuals, x0)
    @info "size(drdx) = $(size(drdx))"
    @info "drdx = $drdx"
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = parameters
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))
    jacobian = isotropic_elasticity_tensor(lambda, mu)
    drde = zeros((19,6))
    drde[1:6, 1:6] = -tovoigt(jacobian)
    @info "drde = $drde"
    @info "size(drde) = $(size(drde))"
    jacobian2 = (drdx\drde)[1:6, 1:6]
    @info "jacobian = $(tovoigt(jacobian))"
    @info "jacobian2 = $jacobian2"
    jacobian3 = (drdx[1:6, 1:6] +  drdx[1:6,7:end]*(drdx[7:end,7:end]\-drdx[7:end, 1:6]))\drde[1:6, 1:6]
    @info "jacobian3 = $jacobian3"
    @info "jacobian4 = $(tovoigt(chabmat.variables_new.jacobian))"
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    g! = create_nonlinear_system_of_equations(chabmat)
    function residuals(x)
        F = similar(x)
        g!(F, x)
        return F
    end

    x0 = [tovoigt(chabmat.variables_new.stress); chabmat.variables_new.R; tovoigt(chabmat.variables_new.X1); tovoigt(chabmat.variables_new.X2)]
    drdx = ForwardDiff.jacobian(residuals, x0)
    @info "size(drdx) = $(size(drdx))"
    @info "drdx = $drdx"
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = parameters
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))
    jacobian = isotropic_elasticity_tensor(lambda, mu)
    drde = zeros((19,6))
    drde[1:6, 1:6] = -tovoigt(jacobian)
    @info "drde = $drde"
    @info "size(drde) = $(size(drde))"
    jacobian2 = (drdx\drde)[1:6, 1:6]
    @info "jacobian = $(tovoigt(jacobian))"
    @info "jacobian2 = $jacobian2"
    jacobian3 = (drdx[1:6, 1:6] +  drdx[1:6,7:end]*(drdx[7:end,7:end]\-drdx[7:end, 1:6]))\drde[1:6, 1:6]
    @info "jacobian3 = $jacobian3"
    @info "jacobian4 = $(tovoigt(chabmat.variables_new.jacobian))"
    update!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
end
#simple_integration_test_fd_tangent2()


text:
using DashBootstrapComponents

items = [
    dbc_dropdownmenuitem("Action"),
    dbc_dropdownmenuitem("Another action"),
    dbc_dropdownmenuitem("Something else here"),
    dbc_dropdownmenuitem(divider=true),
    dbc_dropdownmenuitem("Something else here after the divider"),
];

dropdown = dbc_row([
    dbc_col(
        dbc_dropdownmenu(
            label="Left-aligned menu (default)",
            children=items,
            right=false,
        ),
    ),
    dbc_col(dbc_dropdownmenu(label="Right-aligned menu", children=items, right=true)),
]);

ext: jl

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
for P in (:Year, :Month, :Week, :Day, :Hour, :Minute, :Second, :Millisecond)
    @eval begin
        @inline fastpart(p::$P) = Nanosecond(0)
        @inline slowpart(p::$P) = p
        @inline fastslow(p::$P) = Nanosecond(0), p
    end
end

for P in (:Microsecond, :Nanosecond)
    @eval begin
        @inline fastpart(p::$P) = p
        @inline slowpart(p::$P) = Millisecond(0)
        @inline fastslow(p::$P) = p, Millisecond(0)
    end
end

fastpart(cperiod::CompoundPeriod) = Nanosecond(cperiod) + Microsecond(cperiod)
slowpart(cperiod::CompoundPeriod) = cperiod - fastpart(cperiod)

@inline function fastslow(cperiod::CompoundPeriod)
    fast_part = fastpart(cperiod)
    slow_part = cperiod - fast_part
    return fast_part, slow_part
end

@inline function fastslow(x)
    tm = Time(x)
    fast_time = fastpart(tm)
    slow_time = tm - fast_time
    return fast_time, slow_time
end

function CompoundPeriod(x::TimeDate)
    tim, dat = timedate(x)
    fast_time = fastpart(tim)
    slow_time = slowpart(tim)
    tm = CompoundPeriod(slow_time) + fast_time
    dt = CompoundPeriod(dat)
    return dt + tm
end

function CompoundPeriod(x::TimeDateZone)
    tim, dat = timedate(x)
    fast_time = fastpart(tim)
    slow_time = slowpart(tim)
    tm = CompoundPeriod(slow_time) + fast_time
    dt = CompoundPeriod(dat)
    return dt + tm
end

function (+)(x::TimeDate, y::Period)
    tim, dat = timedate(x)
    fasttime, slowtime = fastpart(tim), slowpart(tim)
    fastperiod, slowperiod = fastslow(y)
    fasttime, lessfast = fastslow(fasttime + fastperiod)
    slowtime = slowtime + lessfast
    datetime = dat + slowtime
    timedate = TimeDate(datetime)
    timedate = timedate + fasttime
    return timedate
end

function (-)(x::TimeDate, y::Period)
    return x + (-y)
end

(+)(y::Period, x::TimeDate) = x + y

function (+)(x::TimeDate, z::CompoundPeriod)
    result = x
    z = canonical(z)
    for period in z
        result += period
    end
    return result
end

function (-)(x::TimeDate, z::CompoundPeriod)
    result = x
    z = canonical(z)
    for period in z
        result -= period
    end
    return result
end

(+)(y::CompoundPeriod, x::TimeDate) = x + y

function (+)(x::TimeDateZone, y::Period)
    fast_part, slow_part = fastslow(y)
    zdt = ZonedDateTime(x)
    fast_time = at_time(x) - at_time(zdt)
    fast_time += fast_part
    zdt = zdt + slow_part
    tim, dat = timedate(zdt)
    tim = tim + fast_time
    inzone = in_zone(zdt)
    atzone = at_zone(zdt)
    TimeDateZone(tim, dat, inzone, atzone)
end

function (-)(x::TimeDateZone, y::Period)
    return x + (-y)
end

(+)(y::Period, x::TimeDateZone) = x + y

function (+)(x::TimeDateZone, y::Period, idx::Int)
    fast_part, slow_part = fastslow(y)
    zdt = ZonedDateTime(x)
    zdt = ZonedDateTime(DateTime(x), in_zone(zdt), idx)
    fast_time = at_time(x) - at_time(zdt)
    fast_time += fast_part
    zdt = zdt + slow_part
    tim, dat = timedate(zdt)
    tim = tim + fast_time
    inzone = in_zone(zdt)
    atzone = at_zone(zdt)
    TimeDateZone(tim, dat, inzone, atzone)
end

function (-)(x::TimeDateZone, y::Period, idx::Int)
    return (+)(x, (-y), idx)
end

(+)(y::Period, x::TimeDateZone, idx::Int) = (+)(x, y, idx)

function (+)(x::TimeDateZone, y::CompoundPeriod)
    result = x
    for period in y
        result += period
    end
    return result
end

function (-)(x::TimeDateZone, y::CompoundPeriod)
    result = x
    for period in y
        result -= period
    end
    return result
end

(+)(y::CompoundPeriod, x::TimeDateZone) = x + y

for P in (:Hour, :Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)
   @eval begin
       (+)(dt::Date, period::$P) = TimeDate(dt) + period
       (-)(dt::Date, period::$P) = TimeDate(dt) - period
   end
end


text:
# This file contains functions for the Radial DEA model
"""
    AbstractRadialDEAModel
An abstract type representing a radial DEA model.
"""
abstract type AbstractRadialDEAModel <: AbstractTechnicalDEAModel end

"""
    RadialDEAModel
An data structure representing a radial DEA model.
"""
struct RadialDEAModel <: AbstractRadialDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
end

"""
    dea(X, Y)
Compute the radial model using data envelopment analysis for inputs X and outputs Y.

# Optional Arguments
- `orient=:Input`: chooses the radially oriented input mode. For the radially oriented output model choose `:Output`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> dea(X, Y)
Radial DEA Model
DMUs = 11; Inputs = 2; Outputs = 1
Orientation = Input; Returns to Scale = CRS
──────────────────────────────────────────────────
    efficiency       slackX1      slackX2  slackY1
──────────────────────────────────────────────────
1     1.0        0.0          0.0              0.0
2     0.62229   -4.41868e-15  0.0              0.0
3     0.819856   0.0          8.17926e-15      0.0
4     1.0       -8.03397e-16  0.0              0.0
5     0.310371   1.80764e-15  0.0              0.0
6     0.555556   4.44444      0.0              0.0
7     1.0        0.0          0.0              0.0
8     0.757669   1.60679e-15  0.0              0.0
9     0.820106   1.64021      0.0              0.0
10    0.490566   9.68683e-15  0.0              0.0
11    1.0        0.0          4.0              0.0
──────────────────────────────────────────────────
```
"""
function dea(X::Matrix, Y::Matrix; orient::Symbol = :Input, rts::Symbol = :CRS, slack = true, Xref::Matrix = X, Yref::Matrix = Y)::RadialDEAModel
    # Check parameters
    nx, m = size(X)
    ny, s = size(Y)

    nrefx, mref = size(Xref)
    nrefy, sref = size(Yref)

    if nx != ny
        error("number of observations is different in inputs and outputs")
    end
    if nrefx != nrefy
        error("number of observations is different in inputs reference set and ouputs reference set")
    end
    if m != mref
        error("number of inputs in evaluation set and reference set is different")
    end
    if s != sref
        error("number of outputs in evaluation set and reference set is different")
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        deamodel = Model(with_optimizer(GLPK.Optimizer))
        @variable(deamodel, eff)
        @variable(deamodel, lambda[1:nref] >= 0)

        if orient == :Input
            # Input orientation
            @objective(deamodel, Min, eff)

            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= eff * x0[j])
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= y0[j])
        elseif orient == :Output
            # Output orientation
            @objective(deamodel, Max, eff)

            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= x0[j])
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= eff * y0[j])
        else
            error("Invalid orientation $orient. Orientation should be :Input or :Output")
        end

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            error("Invalid returns to scale $rts. Returns to scale should be :CRS or :VRS")
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        lambdaeff[i,:] = JuMP.value.(lambda)

    end

    # Compute slacks
    if slack == true

        # Get first-stage efficient X and Y
        if orient == :Input
            Xeff = X .* effi
            Yeff = Y
        elseif orient == :Output
            Xeff = X
            Yeff = Y .* effi
        end

        # Use additive model with radial efficient X and Y to get slacks
        radialSlacks = deaadd(Xeff, Yeff, :Ones, rts = rts, Xref = Xref, Yref = Yref)
        slackX = slacks(radialSlacks, :X)
        slackY = slacks(radialSlacks, :Y)
    else
        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    return RadialDEAModel(n, m, s, orient, rts, effi, slackX, slackY, lambdaeff)

end

function dea(X::Vector, Y::Matrix; orient::Symbol = :Input, rts::Symbol = :CRS, slack = true, Xref::Vector = X, Yref::Matrix = Y)::RadialDEAModel
    X = X[:,:]
    Xref = Xref[:,:]
    return dea(X, Y, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref)
end

function dea(X::Matrix, Y::Vector; orient::Symbol = :Input, rts::Symbol = :CRS, slack = true, Xref::Matrix = X, Yref::Vector = Y)::RadialDEAModel
    Y = Y[:,:]
    Yref = Yref[:,:]
    return dea(X, Y, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref)
end

function dea(X::Vector, Y::Vector; orient::Symbol = :Input, rts::Symbol = :CRS, slack = true, Xref::Vector = X, Yref::Vector = Y)::RadialDEAModel
    X = X[:,:]
    Xref = Xref[:,:]
    Y = Y[:,:]
    Yref = Yref[:,:]
    return dea(X, Y, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref)
end

function Base.show(io::IO, x::RadialDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Radial DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if hasslacks == true
            show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ; ["slackY$i" for i in 1:s ]], ["$i" for i in 1:n]))
        else
            show(io, CoefTable(hcat(eff), ["efficiency"], ["$i" for i in 1:n]))
        end
    end

end


text:
@with_kw mutable struct ConstantParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    γ::T = 1.
end

struct Constant <: AbstractStep
    params::ConstantParameters{Float64}

    function (::Type{Constant})(; kwargs...)
        return new(ConstantParameters(; kwargs...))
    end
end

function stepsize(constant::Constant,klocal::Integer,kglobal::Integer,fval::Real,x::AbstractVector,g::AbstractVector)
    return constant.params.γ
end


text:
######################
# interface for spaces
######################

"""
    dimensions(s::Any)

Returns the number of dimensions in space `s`.
"""
function dimensions end

"""
    states(problem::POMDP)
    states(problem::MDP)

Returns the complete state space of a POMDP.
"""
function states end

"""
    states{S,A,O}(problem::POMDP{S,A,O}, state::S)
    states{S,A}(problem::MDP{S,A}, state::S)

Returns a subset of the state space reachable from `state`.
"""
states(problem::Union{POMDP,MDP}, s) = states(problem)
@impl_dep states(::P,::S) where {P<:Union{POMDP,MDP},S} states(::P)

"""
    actions(problem::POMDP)
    actions(problem::MDP)

Returns the entire action space of a POMDP.
"""
function actions end

"""
    actions{S,A,O}(problem::POMDP{S,A,O}, state::S)
    actions{S,A}(problem::MDP{S,A}, state::S)

Return the action space accessible from the given state.
"""
actions(problem::Union{MDP,POMDP}, state) = actions(problem)
@impl_dep actions(::P,::S) where {P<:Union{POMDP,MDP},S} actions(::P)

"""
    actions(problem::POMDP, belief)

Return the action space accessible from the states with nonzero belief.
"""
actions(problem::POMDP, belief) = actions(problem)

"""
    observations(problem::POMDP)

Return the entire observation space.
"""
function observations end

"""
    observations{S,A,O}(problem::POMDP{S,A,O}, state::S)

Return the observation space accessible from the given state and returns it.
"""
observations(problem::POMDP, state) = observations(problem)
@impl_dep observations(::P,::S) where {P<:POMDP,S} observations(::P)


text:
function snow17(forcings::DataFrame, elev::Float64, par::AbstractArray; prcpCol::Symbol=:precip, tavgCol::Symbol=:tavg, dateCol::Symbol=:datetime, initstate::AbstractArray=[0.,0.,0.,0.])
    # Define constants
    dtt = 24  # time interval of temperature data
    dtp = 24  # time interval of prcipitation data

    # SET PARAMETERS
    SCF    = par[1]
    PXTEMP = par[2]
    MFMAX  = par[3]
    MFMIN  = par[4]
    UADJ   = par[5]
    MBASE  = par[6]
    TIPM   = par[7]
    PLWHC  = par[8]
    NMF    = par[9]
    DAYGM  = par[10]

    dts = forcings[!,dateCol]
    tavg = forcings[!,tavgCol]
    prcp = forcings[!,prcpCol]
    len = length(tavg)

    doy = map(dayofyear, dts)

    meltNrain = Array{Float64}(undef,len)

    # Set initial states
    W_i, ATI, W_q, Deficit = initstate

    # LOOP THROUGH EACH PERIOD
    for i in 1:len

        # Set current temperature and precipitation
        Ta = tavg[i]  # Air temperature at this time step [deg C]
        Pr = prcp[i]  # Precipitation at this time step [mm]

        # FORM OF PRECIPITATION
        if Ta <= PXTEMP
            # Air temperature is cold enough for snow to occur
            SNOW = Pr
            RAIN = 0.0
        else
            # Air temperature is warm enough for rain
            SNOW = 0.0
            RAIN =  Pr
        end

        # ACCUMULATION OF THE SNOW COVER
        Pn = SNOW * SCF  # Water equivalent of new snowfall [mm]
        W_i = W_i + Pn  # Water equivalent of the ice portion of the snow cover [mm]
        E = 0  # Excess liquid water in the snow cover

        # ENERGY EXCHANGE AT SNOW/AIR SURFACE DURING NON-MELT PERIODS

        # Seasonal variation in the non-rain melt factor (Assume a year has 365 days)
        N_Mar21 = doy[i] - 80

        Sv = (0.5 * sin((N_Mar21 * 2 * pi)/365)) + 0.5  # Seasonal variation
        Av = 1  # Seasonal variation adjustment, Av<-1.0 when lat < 54N
        Mf = dtt/6 * ((Sv * Av * (MFMAX - MFMIN)) + MFMIN)  # Seasonally varying non-rain melt factor

        # New snow temperature and heat deficit from new snow
        if (Ta < 0)
            T_snow_new = Ta
        else
            T_snow_new = 0
        end

        # Change in the heat deficit due to new snowfall [mm], 80 cal/g: latent heat of fusion, 0.5 cal/g/C:
        # specific heat of ice
        delta_HD_snow = -(T_snow_new * Pn)/(80/0.5)

        # Heat Exchange due to a temperature gradient change in heat deficit due to a temperature gradient
        # [mm]
        delta_HD_T = NMF * dtp/6 * Mf/MFMAX * (ATI - T_snow_new)

        # Update ATI[Antecedent Temperature Index]
        if (Pn > 1.5 * dtp)
            ATI = T_snow_new  #Antecedent temperature index
        else
            TIPM_dtt = 1 - ((1 - TIPM)^(dtt/6))
            ATI = ATI + TIPM_dtt * (Ta - ATI)
        end

        ATI = min(ATI, 0)

        # SNOW MELT
        T_rain = max(Ta, 0)  # Temperature of rain (deg C), Ta or 0C, whichever greater

        if (RAIN > 0.25 * dtp)
            # Rain-on-Snow Melt
            stefan = 6.12 * (10^(-10))  # Stefan-Boltzman constant (mm/K/hr)
            e_sat = 2.7489 * (10^8) * exp((-4278.63/(Ta + 242.792)))  # Saturated vapor pressure at Ta (mb)
            P_atm = 33.86 * (29.9 - (0.335 * (elev/100)) + (0.00022 * ((elev/100)^2.4)))  # Atmospheric pressure (mb) where elevation is in HUNDREDS of meters (this is incorrectly stated in the manual)
            term1 = stefan * dtp * (((Ta + 273)^4) - (273^4))
            term2 = 0.0125 * RAIN * T_rain
            term3 = 8.5 * UADJ * (dtp/6) * ((0.9 * e_sat - 6.11) + (0.00057 * P_atm * Ta))
            Melt = term1 + term2 + term3
            Melt = max(Melt, 0)

        elseif ((RAIN <= 0.25 * dtp) && (Ta > MBASE))
            # Non-Rain Melt
            Melt = (Mf * (Ta - MBASE) * (dtp/dtt)) + (0.0125 * RAIN * T_rain)
            Melt = max(Melt, 0)

        else
            Melt = 0
        end

        # Ripeness of the snow cover W_i : water equivalent of the ice portion of the snow cover W_q :
        # liquide water held by the snow W_qx: liquid water storage capacity Qw : Amount of available water
        # due to melt and rain

        Deficit = max(Deficit + delta_HD_snow + delta_HD_T, 0)  # Deficit <- heat deficit [mm]

        if (Deficit > (0.33 * W_i))
            # limits of heat deficit
            Deficit = 0.33 * W_i
        end

        if (Melt < W_i)
            W_i = W_i - Melt
            Qw = Melt + RAIN
            W_qx = PLWHC * W_i

            if ((Qw + W_q) > (Deficit + Deficit * PLWHC + W_qx))
                # THEN the snow is RIPE

                E = Qw + W_q - W_qx - Deficit - (Deficit * PLWHC)  # Excess liquid water [mm]
                W_i = W_i + Deficit  # W_i increases because water refreezes as heat deficit is decreased
                W_q = W_qx + PLWHC * Deficit  # fills liquid water capacity
                Deficit = 0

            elseif ((Qw + W_q) >= Deficit)

                # & [[Qw + W_q] <= [[Deficit*[1+PLWHC]] + W_qx]] THEN the snow is NOT yet ripe, but ice is being
                # melted

                E = 0
                W_i = W_i + Deficit  # W_i increases because water refreezes as heat deficit is decreased
                W_q = W_q + Qw - Deficit
                Deficit = 0

            elseif ((Qw + W_q) < Deficit)
                # elseif [[Qw + W_q] < Deficit]

                # THEN the snow is NOT yet ripe
                E = 0
                W_i = W_i + Qw + W_q  # W_i increases because water refreezes as heat deficit is decreased
                Deficit = Deficit - Qw - W_q
            end

        else
            Melt = W_i + W_q  # Melt >= W_i
            W_i = 0
            W_q = 0
            Qw = Melt + RAIN
            E = Qw
            # SWE = 0

        end

        if (Deficit == 0)
            ATI = 0
        end

        # Constant daily amount of melt which takes place at the snow-soil interface whenever there is a
        # snow cover
        if (W_i > DAYGM)

            gmwlos = (DAYGM/W_i) * W_q
            gmslos = DAYGM
            gmro = gmwlos + gmslos
            W_i = W_i - gmslos
            W_q = W_q - gmwlos

            E = E + gmro
            SWE = W_i + W_q

        else
          gmro = W_i + W_q
          W_i = 0
          W_q = 0
          E = E + gmro
          SWE = 0
        end

        meltNrain[i] = E
    end

    return meltNrain
end


text:
using OrdinaryDiffEq, Test

f = function (out,du,u,p,t)
  out[1] = - p[1]*u[1]              + p[3]*u[2]*u[3] - du[1]
  out[2] = + p[1]*u[1] - p[2]*u[2]^2 - p[3]*u[2]*u[3] - du[2]
  out[3] = u[1] + u[2] + u[3] - p[4]
end
u₀ = [1.0, 0, 0]
du₀ = [0.0, 0.0, 0.0]
p = [0.04, 3.0e7, 1.0e4, 1.0]
tspan = (0.0,100.0)
differential_vars = [true,true,false]
prob = DAEProblem(f,du₀,u₀,tspan,p,differential_vars=differential_vars)
condition(u,t,integrator) = t in [50.0]
affect!(integrator) = integrator.p[4] = 2.0
cb = DiscreteCallback(condition,affect!)

#=
# Regression test
using Sundials
sol = solve(prob, IDA(), callback=cb, tstops=[50.0],abstol=1e-14,reltol=1e-14)
=#

p = [0.04, 3.0e7, 1.0e4, 1.0]
prob = DAEProblem(f,du₀,u₀,tspan,p,differential_vars=differential_vars)
sol = solve(prob, DFBDF(), callback=cb, tstops=[50.0],abstol=1e-12,reltol=1e-12)
@test sol.t[end] == 100.0
@test sol[end][1] ≈ 0.686300529575259 atol=1e-7
@test sol[end][2] ≈ 2.0797982209353813e-6 atol=1e-7
@test sol[end][3] ≈ 1.31369739062652 atol=1e-7


text:
module Coad

using CPUTime
using Printf

"""
    AbstractCoadModel

A super type for simulating collision/coalescence problems.
"""
abstract type AbstractCoadModel end

"""
    AbstractSizeDist{FT}

A super type for encapsulating size distribution calculations.
"""
abstract type AbstractSizeDist{FT} end

include("util.jl")
include("size_dists.jl")
include("time.jl")
include("coad1d.jl")

export
  Coad1D, set!, step!, ExponentialDist, nc,

  mass_from_r, r_from_mass

end # module


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

# {{{ FIXME: remove this after we've figure out how to pass through to kernel
const _nvgeo = 14
const _ξx, _ηx, _ζx, _ξy, _ηy, _ζy, _ξz, _ηz, _ζz, _MJ, _MJI,
       _x, _y, _z = 1:_nvgeo

const _nsgeo = 5
const _nx, _ny, _nz, _sMJ, _vMJI = 1:_nsgeo
# }}}

"""
    volumerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
               ::Val{nauxstate}, flux!, source!, rhs, Q, Qvisc, auxstate,
               vgeo, t, D, elems) where {dim, N, nstate, nviscstate,

Computational kernel: Evaluate the volume integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function volumerhs!(::Val{dim}, ::Val{N},
                    ::Val{nstate}, ::Val{nviscstate},
                    ::Val{nauxstate},
                    flux!, source!,
                    rhs, Q, Qvisc, auxstate, vgeo, t,
                    D, elems) where {dim, N, nstate, nviscstate,
                                     nauxstate}
  DFloat = eltype(Q)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_F = @shmem DFloat (3, Nq, Nq, Nqk, nstate)
  s_D = @shmem DFloat (Nq, Nq)
  l_rhs = @scratch DFloat (nstate, Nq, Nq, Nqk) 3

  source! !== nothing && (l_S = MArray{Tuple{nstate}, DFloat}(undef))
  l_Q = MArray{Tuple{nstate}, DFloat}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_aux = MArray{Tuple{nauxstate}, DFloat}(undef)
  l_F = MArray{Tuple{3, nstate}, DFloat}(undef)

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_D[i, j] = D[i, j]
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          MJ = vgeo[ijk, _MJ, e]
          ξx, ξy, ξz = vgeo[ijk,_ξx,e], vgeo[ijk,_ξy,e], vgeo[ijk,_ξz,e]
          ηx, ηy, ηz = vgeo[ijk,_ηx,e], vgeo[ijk,_ηy,e], vgeo[ijk,_ηz,e]
          ζx, ζy, ζz = vgeo[ijk,_ζx,e], vgeo[ijk,_ζy,e], vgeo[ijk,_ζz,e]

          @unroll for s = 1:nstate
            l_rhs[s, i, j, k] = rhs[ijk, s, e]
          end

          @unroll for s = 1:nstate
            l_Q[s] = Q[ijk, s, e]
          end

          @unroll for s = 1:nviscstate
            l_Qvisc[s] = Qvisc[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s] = auxstate[ijk, s, e]
          end

          flux!(l_F, l_Q, l_Qvisc, l_aux, t)

          @unroll for s = 1:nstate
            s_F[1,i,j,k,s] = MJ * (ξx*l_F[1,s] + ξy*l_F[2,s] + ξz*l_F[3,s])
            s_F[2,i,j,k,s] = MJ * (ηx*l_F[1,s] + ηy*l_F[2,s] + ηz*l_F[3,s])
            s_F[3,i,j,k,s] = MJ * (ζx*l_F[1,s] + ζy*l_F[2,s] + ζz*l_F[3,s])
          end

          if source! !== nothing
            source!(l_S, l_Q, l_aux, t)

            @unroll for s = 1:nstate
              l_rhs[s, i, j, k] += l_S[s]
            end
          end
        end
      end
    end
    @synchronize

    @unroll for s = 1:nstate
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            ijk = i + Nq * ((j-1) + Nq * (k-1))
            MJI = vgeo[ijk, _MJI, e]
            for n = 1:Nq
              Dni = s_D[n, i]
              Dnj = s_D[n, j]
              Nqk > 1 && (Dnk = s_D[n, k])
              # ξ-grid lines
              l_rhs[s, i, j, k] += MJI * Dni * s_F[1, n, j, k, s]

              # η-grid lines
              l_rhs[s, i, j, k] += MJI * Dnj * s_F[2, i, n, k, s]

              # ζ-grid lines
              Nqk > 1 && (l_rhs[s, i, j, k] += MJI * Dnk * s_F[3, i, j, n, s])
            end
          end
        end
      end
    end
    @unroll for s = 1:nstate
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            ijk = i + Nq * ((j-1) + Nq * (k-1))
            rhs[ijk, s, e] = l_rhs[s, i, j, k]
          end
        end
      end
    end
    @synchronize
  end
  nothing
end

"""
    facerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
             ::Val{nauxstate}, numerical_flux!,
             numerical_boundary_flux!, rhs, Q, Qvisc, auxstate,
             vgeo, sgeo, t, vmapM, vmapP, elemtobndy,
             elems) where {dim, N, nstate, nviscstate, nauxstate}

Computational kernel: Evaluate the surface integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function facerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
                  ::Val{nauxstate}, numerical_flux!, numerical_boundary_flux!,
                  rhs, Q, Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                  elemtobndy, elems) where {dim, N, nstate, nviscstate,
                                            nauxstate}

  DFloat = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  l_QM = MArray{Tuple{nstate}, DFloat}(undef)
  l_QviscM = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_auxM = MArray{Tuple{nauxstate}, DFloat}(undef)

  l_QP = MArray{Tuple{nstate}, DFloat}(undef)
  l_QviscP = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_auxP = MArray{Tuple{nauxstate}, DFloat}(undef)

  l_F = MArray{Tuple{nstate}, DFloat}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    for f = 1:nface
      @loop for n in (1:Nfp; threadIdx().x)
        nM = (sgeo[_nx, n, f, e], sgeo[_ny, n, f, e], sgeo[_nz, n, f, e])
        sMJ, vMJI = sgeo[_sMJ, n, f, e], sgeo[_vMJI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        @unroll for s = 1:nstate
          l_QM[s] = Q[vidM, s, eM]
        end

        @unroll for s = 1:nviscstate
          l_QviscM[s] = Qvisc[vidM, s, eM]
        end

        @unroll for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        # Load plus side data
        @unroll for s = 1:nstate
          l_QP[s] = Q[vidP, s, eP]
        end

        @unroll for s = 1:nviscstate
          l_QviscP[s] = Qvisc[vidP, s, eP]
        end

        @unroll for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        bctype =
            numerical_boundary_flux! === nothing ? 0 : elemtobndy[f, e]
        if bctype == 0
          numerical_flux!(l_F, nM, l_QM, l_QviscM, l_auxM, l_QP, l_QviscP,
                          l_auxP, t)
        else
          numerical_boundary_flux!(l_F, nM, l_QM, l_QviscM, l_auxM, l_QP,
                                   l_QviscP, l_auxP, bctype, t)
        end

        #Update RHS
        @unroll for s = 1:nstate
          # FIXME: Should we pretch these?
          rhs[vidM, s, eM] -= vMJI * sMJ * l_F[s]
        end
      end
      # Need to wait after even faces to avoid race conditions
      f % 2 == 0 && @synchronize
    end
  end
  nothing
end

function volumeviscterms!(::Val{dim}, ::Val{N}, ::Val{nstate},
                          ::Val{states_grad}, ::Val{ngradstate},
                          ::Val{nviscstate}, ::Val{nauxstate},
                          viscous_transform!, gradient_transform!, Q,
                          Qvisc, auxstate, vgeo, t, D,
                          elems) where {dim, N, states_grad, ngradstate,
                                        nviscstate, nstate, nauxstate}
  DFloat = eltype(Q)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  ngradtransformstate = length(states_grad)

  s_G = @shmem DFloat (Nq, Nq, Nqk, ngradstate)
  s_D = @shmem DFloat (Nq, Nq)

  l_Q = @scratch DFloat (ngradtransformstate, Nq, Nq, Nqk) 3
  l_aux = @scratch DFloat (nauxstate, Nq, Nq, Nqk) 3
  l_G = MArray{Tuple{ngradstate}, DFloat}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_gradG = MArray{Tuple{3, ngradstate}, DFloat}(undef)

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_D[i, j] = D[i, j]
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:ngradtransformstate
            l_Q[s, i, j, k] = Q[ijk, states_grad[s], e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s, i, j, k] = auxstate[ijk, s, e]
          end

          gradient_transform!(l_G, l_Q[:, i, j, k], l_aux[:, i, j, k], t)
          @unroll for s = 1:ngradstate
            s_G[i, j, k, s] = l_G[s]
          end
        end
      end
    end
    @synchronize

    # Compute gradient of each state
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          ξx, ξy, ξz = vgeo[ijk, _ξx, e], vgeo[ijk, _ξy, e], vgeo[ijk, _ξz, e]
          ηx, ηy, ηz = vgeo[ijk, _ηx, e], vgeo[ijk, _ηy, e], vgeo[ijk, _ηz, e]
          ζx, ζy, ζz = vgeo[ijk, _ζx, e], vgeo[ijk, _ζy, e], vgeo[ijk, _ζz, e]

          @unroll for s = 1:ngradstate
            Gξ = Gη = Gζ = zero(DFloat)
            @unroll for n = 1:Nq
              Din = s_D[i, n]
              Djn = s_D[j, n]
              Nqk > 1 && (Dkn = s_D[k, n])

              Gξ += Din * s_G[n, j, k, s]
              Gη += Djn * s_G[i, n, k, s]
              Nqk > 1 && (Gζ += Dkn * s_G[i, j, n, s])
            end
            l_gradG[1, s] = ξx * Gξ + ηx * Gη + ζx * Gζ
            l_gradG[2, s] = ξy * Gξ + ηy * Gη + ζy * Gζ
            l_gradG[3, s] = ξz * Gξ + ηz * Gη + ζz * Gζ
          end

          viscous_transform!(l_Qvisc, l_gradG, l_Q[:, i, j, k],
                             l_aux[:, i, j, k], t)

          @unroll for s = 1:nviscstate
            Qvisc[ijk, s, e] = l_Qvisc[s]
          end
        end
      end
    end
    @synchronize
  end
end

function faceviscterms!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{states_grad},
                        ::Val{ngradstate}, ::Val{nviscstate},
                        ::Val{nauxstate}, viscous_penalty!,
                        viscous_boundary_penalty!, gradient_transform!,
                        Q, Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                        elemtobndy, elems) where {dim, N, states_grad,
                                                  ngradstate, nviscstate,
                                                  nstate, nauxstate}
  DFloat = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  ngradtransformstate = length(states_grad)

  l_QM = MArray{Tuple{ngradtransformstate}, DFloat}(undef)
  l_auxM = MArray{Tuple{nauxstate}, DFloat}(undef)
  l_GM = MArray{Tuple{ngradstate}, DFloat}(undef)

  l_QP = MArray{Tuple{ngradtransformstate}, DFloat}(undef)
  l_auxP = MArray{Tuple{nauxstate}, DFloat}(undef)
  l_GP = MArray{Tuple{ngradstate}, DFloat}(undef)

  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    for f = 1:nface
      @loop for n in (1:Nfp; threadIdx().x)
        nM = (sgeo[_nx, n, f, e], sgeo[_ny, n, f, e], sgeo[_nz, n, f, e])
        sMJ, vMJI = sgeo[_sMJ, n, f, e], sgeo[_vMJI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        @unroll for s = 1:ngradtransformstate
          l_QM[s] = Q[vidM, states_grad[s], eM]
        end

        @unroll for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        gradient_transform!(l_GM, l_QM, l_auxM, t)

        # Load plus side data
        @unroll for s = 1:ngradtransformstate
          l_QP[s] = Q[vidP, states_grad[s], eP]
        end

        @unroll for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        gradient_transform!(l_GP, l_QP, l_auxP, t)

        bctype =
            viscous_boundary_penalty! === nothing ? 0 : elemtobndy[f, e]
        if bctype == 0
          viscous_penalty!(l_Qvisc, nM, l_GM, l_QM, l_auxM, l_GP,
                                  l_QP, l_auxP, t)
        else
          viscous_boundary_penalty!(l_Qvisc, nM, l_GM, l_QM, l_auxM,
                                           l_GP, l_QP, l_auxP, bctype, t)
        end

        @unroll for s = 1:nviscstate
          Qvisc[vidM, s, eM] += vMJI * sMJ * l_Qvisc[s]
        end
      end
      # Need to wait after even faces to avoid race conditions
      f % 2 == 0 && @synchronize
    end
  end
  nothing
end

"""
    initauxstate!(::Val{dim}, ::Val{N}, ::Val{nauxstate}, auxstatefun!,
                  auxstate, vgeo, elems) where {dim, N, nauxstate}

Computational kernel: Initialize the auxiliary state

See [`DGBalanceLaw`](@ref) for usage.
"""
function initauxstate!(::Val{dim}, ::Val{N}, ::Val{nauxstate}, auxstatefun!,
                       auxstate, vgeo, elems) where {dim, N, nauxstate}

  DFloat = eltype(auxstate)

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Np = Nq * Nq * Nqk

  l_aux = MArray{Tuple{nauxstate}, DFloat}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      x, y, z = vgeo[n, _x, e], vgeo[n, _y, e], vgeo[n, _z, e]
      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end

      auxstatefun!(l_aux, x, y, z)

      @unroll for s = 1:nauxstate
        auxstate[n, s, e] = l_aux[s]
      end
    end
  end
end

"""
    elem_grad_field!(::Val{dim}, ::Val{N}, ::Val{nstate}, Q, vgeo, D, elems, s,
                     sx, sy, sz) where {dim, N, nstate}

Computational kernel: Compute the element gradient of state `s` of `Q` and store
it in `sx`, `sy`, and `sz` of `Q`.

!!! warning

    This does not compute a DG gradient, but only over the element. If ``Q_s``
    is discontinuous you may want to consider another approach.

"""
function elem_grad_field!(::Val{dim}, ::Val{N}, ::Val{nstate}, Q, vgeo,
                          D, elems, s, sx, sy, sz) where {dim, N, nstate}

  DFloat = eltype(vgeo)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_f = @shmem DFloat (Nq, Nq, Nqk)
  s_D = @shmem DFloat (Nq, Nq)

  l_fξ = @scratch DFloat (Nq, Nq, Nqk) 3
  l_fη = @scratch DFloat (Nq, Nq, Nqk) 3
  l_fζ = @scratch DFloat (Nq, Nq, Nqk) 3

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_D[i, j] = D[i, j]
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          s_f[i, j, k] = Q[ijk, s, e]
        end
      end
    end
    @synchronize

    # reference gradient
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          l_fξ[i, j, k] = 0
          l_fη[i, j, k] = 0
          l_fζ[i, j, k] = 0
          @unroll for n = 1:Nq
            Din = s_D[i, n]
            Djn = s_D[j, n]
            Nqk > 1 && (Dkn = s_D[k, n])

            # ξ-grid lines
            l_fξ[i, j, k] += Din * s_f[n, j, k]

            # η-grid lines
            l_fη[i, j, k] += Djn * s_f[i, n, k]

            # ζ-grid lines
            Nqk > 1 && (l_fζ[i, j, k] += Dkn * s_f[i, j, n])
          end
        end
      end
    end

    # Physical gradient
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))

          ξx, ξy, ξz = vgeo[ijk, _ξx, e], vgeo[ijk, _ξy, e], vgeo[ijk, _ξz, e]
          ηx, ηy, ηz = vgeo[ijk, _ηx, e], vgeo[ijk, _ηy, e], vgeo[ijk, _ηz, e]
          ζx, ζy, ζz = vgeo[ijk, _ζx, e], vgeo[ijk, _ζy, e], vgeo[ijk, _ζz, e]

          Q[ijk, sx, e] = ξx * l_fξ[ijk] + ηx * l_fη[ijk] + ζx * l_fζ[ijk]
          Q[ijk, sy, e] = ξy * l_fξ[ijk] + ηy * l_fη[ijk] + ζy * l_fζ[ijk]
          Q[ijk, sz, e] = ξz * l_fξ[ijk] + ηz * l_fη[ijk] + ζz * l_fζ[ijk]
        end
      end
    end
    @synchronize
  end
end

"""
    knl_dof_iteration!(::Val{dim}, ::Val{N}, ::Val{nRstate}, ::Val{nstate},
                       ::Val{nviscstate}, ::Val{nauxstate}, dof_fun!, R, Q,
                       QV, auxstate, elems) where {dim, N, nRstate, nstate,
                                                   nviscstate, nauxstate}

Computational kernel: fill postprocessing array

See [`DGBalanceLaw`](@ref) for usage.
"""
function knl_dof_iteration!(::Val{dim}, ::Val{N}, ::Val{nRstate}, ::Val{nstate},
                            ::Val{nviscstate}, ::Val{nauxstate}, dof_fun!, R, Q,
                            QV, auxstate, elems) where {dim, N, nRstate, nstate,
                                                        nviscstate, nauxstate}
  DFloat = eltype(R)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  Np = Nq * Nq * Nqk

  nelem = size(auxstate)[end]

  l_R = MArray{Tuple{nRstate}, DFloat}(undef)
  l_Q = MArray{Tuple{nstate}, DFloat}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_aux = MArray{Tuple{nauxstate}, DFloat}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      @unroll for s = 1:nRstate
        l_R[s] = R[n, s, e]
      end

      @unroll for s = 1:nstate
        l_Q[s] = Q[n, s, e]
      end

      @unroll for s = 1:nviscstate
        l_Qvisc[s] = QV[n, s, e]
      end

      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end

      dof_fun!(l_R, l_Q, l_Qvisc, l_aux)

      @unroll for s = 1:nRstate
        R[n, s, e] = l_R[s]
      end
    end
  end
end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
"""
    mel = MELink{Tx,Ty}(dummylabel)

Create a linked-list storing labels and x/y pairs that represent the "minimum edge" of a
set of values, where the stored x and y values are monotonic in both x and y. Insertion
of a new x/y pair expunges any entries with smaller x but larger y. Add new elements with
`insert!(mel, x, label=>y)`.

The first item in the list is a dummy item with `x=0`, so that the head of the list is constant.
"""
mutable struct MELink{Tl,Tw,Tf}
    l::Tl   # label
    w::Tw   # x-coordinate (box width)
    f::Tf   # y-coordinate (function value)
    next::MELink{Tl,Tw,Tf}

    function MELink{Tl,Tw,Tf}(l, w, f) where {Tw, Tf, Tl}
        mel = new{Tl,Tw,Tf}(l, w, f)
        mel.next = mel
        return mel
    end
    function MELink{Tl,Tw,Tf}(l, w, f, next) where {Tw, Tf, Tl}
        mel = new{Tl,Tw,Tf}(l, w, f, next)
        return mel
    end
end

MELink{Tw,Tf}(dummylabel::Tl) where {Tl,Tw,Tf} = MELink{Tl,Tw,Tf}(dummylabel, typemin(Tw), typemin(Tf))

Base.IteratorSize(::Type{<:MELink}) = Base.SizeUnknown()

# Mutable length-3 vector without the storage overhead of Array
mutable struct MVector3{T} <: AbstractVector{T}
    x1::T
    x2::T
    x3::T
end
MVector3{T}(a::AbstractVector) where T = convert(MVector3{T}, a)

Base.IndexStyle(::Type{<:MVector3}) = IndexLinear()
Base.size(v::MVector3) = (3,)
Base.axes1(v::MVector3) = Base.OneTo(3)
function Base.getindex(v::MVector3, i::Int)
    @boundscheck ((1 <= i) & (i <= 3)) || Base.throw_boundserror(v, i)
    ifelse(i == 1, v.x1, ifelse(i == 2, v.x2, v.x3))
end
function Base.setindex!(v::MVector3, val, i::Int)
    @boundscheck ((1 <= i) & (i <= 3)) || Base.throw_boundserror(v, i)
    if i == 1
        v.x1 = val
    elseif i == 2
        v.x2 = val
    else
        v.x3 = val
    end
    v
end
Base.convert(::Type{MVector3{T}}, a::MVector3{T}) where T = a
function Base.convert(::Type{MVector3{T}}, a::AbstractVector) where T
    @boundscheck axes(a) == (Base.OneTo(3),) || throw(DimensionMismatch("vector must have length 3"))
    @inbounds ret = MVector3{T}(a[1], a[2], a[3])
    return ret
end
Base.convert(::Type{MVector3}, a::AbstractVector{T}) where T =
    convert(MVector3{T}, a)

Base.similar(::MVector3{T}, ::Type{S}, dims::Tuple{Vararg{Int}}) where {T,S} =
    Array{S}(undef, dims)

boxeltype(::Type{<:Integer}) = Float64
boxeltype(::Type{T}) where T = T

mutable struct Box{T,N}
    parent::Box{T,N}
    parent_cindex::Int    # of its parent's children, which one is this?
    splitdim::Int         # the dimension along which this box has been split, or 0 if this is a leaf node
    qnconverged::Bool     # true if the box was the minimum-value box in a Quasi-Newton step and the improvement was within convergence criterion
    minmax::Tuple{T,T}    # the outer edges not corresponding to one of the parent's xvalues (splits that occur between dots in Fig 2)
    xvalues::MVector3{T}  # the values of x_splitdim at which f is evaluated
    fvalues::MVector3{T}  # the corresponding values of f
    children::MVector3{Box{T,N}}

    function default!(box::Box{T,N}) where {T,N}
        box.qnconverged = false
        box.minmax = (typemax(T), typemin(T))
        box.xvalues = MVector3{T}(typemax(T), zero(T), typemin(T))
        box.fvalues = MVector3{T}(zero(T), zero(T), zero(T))
        box.children = MVector3(box, box, box)
    end
    function Box{T,N}() where {T,N}
        # Create the root box
        box = new{T,N}()
        box.parent = box
        box.parent_cindex = 0
        box.splitdim = 0
        default!(box)
        return box
    end
    function Box{T,N}(parent::Box, parent_cindex::Integer) where {T,N}
        # Create a new child and store it in the parent
        box = new{T,N}(parent, parent_cindex, 0)
        parent.children[parent_cindex] = box
        default!(box)
        box
    end
end

Box(parent::Box{T,N}, parent_cindex::Integer) where {T,N} = Box{T,N}(parent, parent_cindex)

isroot(box::Box) = box.parent == box
isleaf(box::Box) = box.splitdim == 0
Base.parent(box::Box) = box.parent
Base.ndims(box::Box{T,N}) where {T,N} = N
Base.IteratorSize(::Type{<:Box}) = Base.SizeUnknown()

abstract type WrappedFunction <: Function end

# A function that keeps track of the number of evaluations
mutable struct CountedFunction{F} <: WrappedFunction
    f::F
    evals::Int

    CountedFunction{F}(f) where F = new{F}(f, 0)
end
CountedFunction(f::Function) = CountedFunction{typeof(f)}(f)

function (f::CountedFunction{F})(x::AbstractVector) where F
    f.evals += 1
    return f.f(x)
end

numevals(f::CountedFunction) = f.evals

# A function that keeps track of the number of evaluations
mutable struct LoggedFunction{F} <: WrappedFunction
    f::F
    values::Vector{Float64}

    LoggedFunction{F}(f) where F = new{F}(f, Float64[])
end
LoggedFunction(f::Function) = LoggedFunction{typeof(f)}(f)

function (f::LoggedFunction{F})(x::AbstractVector) where F
    val = f.f(x)
    push!(f.values, val)
    return val
end

numevals(f::LoggedFunction) = length(f.values)

# Quadratic-model Incremental Gaussian Elimination
struct QmIGE{T,N}
    coefs::PermutedDimsArray{T,2,(2, 1),(2, 1),Array{T,2}} # to make row-major operations fast
    rhs::Vector{T}
    dimpiv::Vector{Int}       # order in which dimensions were added
    rowzero::Vector{Bool}     # true if a row in coefs is all-zeros
    rowbox::Vector{Box{T,N}}  # the box that set this row
    ndims::typeof(Ref(1))     # number of dimensions added so far
    nzrows::typeof(Ref(1))
    rowtmp::Vector{T}
end

function QmIGE{T,N}() where {T,N}
    m = ((N+1)*(N+2))÷2 - 1  # the constant is handled separately
    coefs = PermutedDimsArray(zeros(T, m, m), (2, 1))
    rhs = zeros(T, m)
    rowtmp = zeros(T, m)
    dimpiv = zeros(Int, N)
    rowzero = fill(true, m)
    rowbox = Vector{Box{T,N}}(undef, m)
    QmIGE{T,N}(coefs, rhs, dimpiv, rowzero, rowbox, Ref(0), Ref(0), rowtmp)
end


text:
struct CrayonWrapper
    c::Crayon
    v::Vector{Union{CrayonWrapper,String}}
end

function (c::Crayon)(args::Union{CrayonWrapper,AbstractString}...)
    typefix(cw::CrayonWrapper) = cw
    typefix(str) = String(str)

    CrayonWrapper(c, typefix.(collect(args)))
end

Base.show(io::IO, cw::CrayonWrapper) = _show(io, cw, CrayonStack(incremental = true))

_show(io::IO, str::String, stack::CrayonStack) = print(io, str)

function _show(io::IO, cw::CrayonWrapper, stack::CrayonStack)
    print(io, push!(stack, cw.c))
    for obj in cw.v
        _show(io, obj, stack)
    end
    length(stack.crayons) > 1 && print(io, pop!(stack))
    return
end

Base.:*(c::Crayon, cw::CrayonWrapper) = CrayonWrapper(c * cw.c, cw.v)
Base.:*(cw::CrayonWrapper, c::Crayon) = CrayonWrapper(cw.c * c, cw.v)


text:
using DataFrames, Parameters, CSV, Distributed

#install via:
# using Pkg; Pkg.add("https://github.com/pszufe/TBBSim.jl.git")
using TBBSim

addprocs(2) #change to the desired number of processes
@everywhere using DataFrames, Parameters, CSV, Distributed
@everywhere using TBBSim

p = MarketParams(nA=10,home_ln_avg=6.5,home_ln_std=0.10,nB=100)
pp = (nH=[70,80,90,100,110,120,130,140], buyer_value_p=[0.9,0.95,1.0,1.05,1.1,1.15,1.2,1.25], buyer_value_std=[0.05,0.1,0.15,0.2,0.25,0.30])
vals = vec(collect(Base.Iterators.product(pp.nH,pp.buyer_value_p,pp.buyer_value_std)))
sweep = [MarketParams(p,Dict(:nH => val[1], :buyer_value_p=>val[2], :buyer_value_std=>val[3] )) for val in vals]

big_res = @distributed (append!) for i in 1:length(sweep)
    println("############ sweep=$i ")
    res = runsims(sweep[i],nm=30,nN=10)
    tablerow = res_agg(res,sweep[i];sim_params_id=i)
    CSV.write("TAB_2$i", tablerow,delim='\t')
    tablerow
end

println("All OK")
# write the aggregated result to a file
CSV.write("tbb_sim_res_2.txt", big_res,delim='\t')
# write the aggregated result to the standard output
CSV.write(stdout, big_res,delim='\t')


text:
include("symbolic_repr.jl")

# Fillers list
# For learning
α = 0.1

#---
# For buddha
bu = symbolic_repr.filler("bu",true, false,α+.2) # model higher frequency of bu
bhu = symbolic_repr.filler("bhu",true, true,α)
d = symbolic_repr.filler("d",true, false,α)
dh = symbolic_repr.filler("dh", true, true,α)
t = symbolic_repr.filler("t", false, false,α)
dhi = symbolic_repr.filler("dhi", false, true,α)
di = symbolic_repr.filler("di", false, false,α)
ti = symbolic_repr.filler("ti", false, false,α)
dha = symbolic_repr.filler("dha", true, true,α)
ta = symbolic_repr.filler("ta", false, false,α)

#---
# bodhami
bo = symbolic_repr.filler("bo",true, false,α)
bho = symbolic_repr.filler("bho",true, true,α)
da = symbolic_repr.filler("da",true, false,α)
mi = symbolic_repr.filler("mi",true,false,1) # discrete
# Define possible fillers for possible roles

#---
# bhotsyati
syati = symbolic_repr.filler("syati",false,false,1) # discrete
#---

# Define possible fillers for possible roles
initial = [bo, bho]
final = [d, dh,t,dhi,ti,di]
suffix = [syati]
should_win = "bhotsyati"

# Constraints
MAX = 1
DEP = -1
VOICE = -1
LICENSE = -1

candidates = symbolic_repr.build_candidates(initial, final, suffix)

"""Update constraints weights"""
function update_constraints(winner::String, optimal::String, violations::Vector)
    global MAX
    global DEP
    global VOICE
    global LICENSE
    #Update Constraints
    ## Find indices
    optViol_idx = findall(x -> x["candidate"] == optimal, violations)
    winViol_idx = findall(x -> x["candidate"] == winner, violations)
    ## Find dictionaries entries
    optViol = violations[optViol_idx[1]]
    winViol = violations[winViol_idx[1]]
    # Improve constraints
    optViol["MaxIO"] > winViol["MaxIO"] ? MAX -= .5 : MAX += .5
    optViol["Dep"] > winViol["Dep"] ? DEP += .5 : DEP -= .5
    optViol["Voice"] > winViol["Voice"] ? VOICE += .5 : VOICE -= .5
    optViol["License"] > winViol["License"] ? LICENSE += .5 : LICENSE -= .5
end

# Learning algorithm
error = true # control var
while error
# Calculate harmony
global violations = []
for c in candidates
    # Initialize violations
    cand_viols = Dict("candidate" => c.repr,
                      "MaxIO" => 0.0,
                      "Dep" => 0.0,
                      "Voice" => 0.0,
                      "License" => 0.0)

    #maxIO
    maxio = MAX * (c.root_initial.act + c.root_final.act + c.suffix.act)
    cand_viols["MaxIO"] = maxio

    dep = DEP * (3 - (c.root_initial.act + c.root_final.act + c.suffix.act))
    cand_viols["Dep"] = dep

    harm = maxio + dep

    # Add voice violations
    if c.root_final.voiced == !c.suffix.voiced
        cand_viols["Voice"] += VOICE
        harm += VOICE
    end # voice violations

    # Add License(lar) violations
    if c.root_final.aspirated && !c.suffix.voiced
        if c.root_final ∉ [dhi, ti, di]
            cand_viols["License"] += LICENSE
            #display("$(c) has License violations")
            harm += LICENSE
        end
    end

    # Add epenthesis violations
    if c.root_final in [dhi, ti, di]
        cand_viols["Dep"] += DEP
        harm += DEP
    end # end epenthesis check

    # Add candidate and harmony to tableau/candidates
    c.harmony = harm
    push!(violations, cand_viols)
end # candidates loop

# Find winner
winner, maxh = symbolic_repr.maxHarm(candidates)

# Define optimal candidate
optimal_index = findall(x -> x.repr == should_win, candidates)
optimal = candidates[optimal_index][1]

# case 1 : more than one winner:
indices_winners = findall(x -> x.harmony == maxh, candidates)

if length(indices_winners) > 1
    # increase act optimal
    symbolic_repr.increase_activations(optimal)

    # remove optimal candidate
    indices_winners = filter(x -> x ≠ optimal_index, indices_winners)

    for c in candidates[indices_winners]
        symbolic_repr.decrease_activations(c)
    end # decrease suboptimal

    # Update constraints weights
    update_constraints(winner, optimal.repr,violations)

# Case 2 : unique winner but not the right one
elseif winner != optimal.repr
    win_idx = findall(x -> x.repr == winner, candidates)
    symbolic_repr.decrease_activations(candidates[win_idx[1]])
    symbolic_repr.increase_activations(optimal)

    #Update Constraints
    update_constraints(winner, optimal.repr,violations)

else
    error = false
    display("The winner is : $(winner)")
    break
end # check condition
end # while loop

#display(candidates)
open("results_$(should_win).txt", "w") do f
for i in candidates
    println(f,i)
end
end

# Violations
open("violations_$(should_win).txt", "w") do f
for cand in violations
    println(f,"--------------------------------------------")
    for (key, value) in cand
    println(f, "$(key) = $(value)")
    end
    println(f,"--------------------------------------------")
end
end

# Violations
open("Constraints_weights_$(should_win).txt", "w") do f
println(f, "MAX : $(MAX)")
println(f, "DEP : $(DEP)")
println(f, "License : $(LICENSE)")
println(f, "Voice Contrast : $(VOICE)")
end

max_stars_repo_name: amasotti/grascom

max_issues_repo_name: amasotti/grascom

max_forks_repo_name: amasotti/grascom


text:
using GNSSTools
using Statistics
using ProgressMeter
using PyPlot
pygui(true)

CN0s = [25, 30, 35, 40]  # dB⋅Hz
h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
# h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16]./10e6^2   # OCXO

t_length = 5.  # seconds
wait_time = 1.  # seconds (the time to wait before calculating STD of phase)
f_s = 5e6  # Hz
sig_freq = L1_freq  # Hz
f_d = 0.  # Hz
f_if = 0.  # Hz
fd_rates = 0:-5:-350  # Hz/s
integration_times = range(1e-3, 40e-3, step=1e-3)
iterations = 1
acquisition_T = 1e-3  # seconds
fine_acq_T = 10e-3  # seconds
fine_acq_method = :fft
prn = 1

sigtype = define_l1ca_code_type(;sig_freq=sig_freq)
signal = definesignal(sigtype, f_s, t_length; prn=prn, f_d=f_d, f_if=f_if)

results = Array{Float64}(undef, length(CN0s), length(fd_rates), length(integration_times))

p = Progress(iterations*length(CN0s)*length(fd_rates)*length(integration_times), 1, "Processing...")
for i in 1:length(CN0s)
    for j in 1:length(fd_rates)
        for k in 1:length(integration_times)
            Δϕσ = 0.
            for iteration in 1:iterations
                CN0 = CN0s[i]
                fd_rate = fd_rates[j]
                T = integration_times[k]
                ϕ₀ = rand(0:0.0001:2π)
                code_start_idx = rand(1:signal.sample_num)
                definesignal!(signal; fd_rate=fd_rate, CN0=CN0,
                              code_start_idx=code_start_idx, phi=ϕ₀,
                              new_thermal_noise=true, new_phase_noise=true)
                generatesignal!(signal)
                # println([i,j,k])
                acqresults, trackresults = process(signal, sigtype, prn,
                                                   acquisition_T=acquisition_T,
                                                   fine_acq_T=fine_acq_T,
                                                   fine_acq_method=fine_acq_method,
                                                   tracking_T=T, show_plot=false,
                                                   h₋₂=h_parms[1], h₀=h_parms[3])
                dphi_start_idx = floor(Int, wait_time/T)
                Δϕσ += std(trackresults.dphi_meas[dphi_start_idx:end])*180/π
                next!(p)
            end
            results[i,j,k] = Δϕσ/iterations
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Gama"
version = v"2.16.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/gama/gama-$(version.major).$(version.minor).tar.gz","7ced801d99ea46f085d7e025283783df2822bd6df35c2b85d743b69ca2f6096b")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gama-*

if [[ "${target}" == x86_64-linux-musl ]] || [[ "${target}" == x86_64-unknown-freebsd* ]]; then
    export CPPFLAGS="-I${includedir}"
fi

./configure \
--prefix=${prefix} \
--build=${MACHTYPE} \
--host=${target}

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental = true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("gama-local-xml2txt", :gama_local_xml2txt),
    ExecutableProduct("gama-local-xml2sql", :gama_local_xml2sql),
    ExecutableProduct("gama-local-yaml2gkf", :gama_local_yaml2gkf),
    ExecutableProduct("gama-local", :gama_local),
    ExecutableProduct("gama-local-gkf2yaml", :gama_local_gkf2yaml),
    ExecutableProduct("gama-g3", :gama_g3)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Expat_jll"; compat="2.2.10")
    Dependency("yaml_cpp_jll")
    Dependency("SQLite_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
#needs a c++ 14 compliant compiler
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil