text:
# # DYCOMS:
#
## Introduction
#
# This driver defines the initial condition for the
# Dynamics and Chemistry of Marine Stratocumulus (DYCOMS)
# LES test described in
#
# [1] B. Stevens et al. (2005) Evaluation of Large-Eddy Simulations via Observations of Nocturnal Marine Stratocumulus, Mon. Wea. Rev. 133:1443-1462
#
#
# Below is a program interspersed with comments.
#md # The full program, without comments, can be found in the next
#md # [section](@ref ex_001_periodic_advection-plain-program).
#
# ## Commented Program

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

# ### Preliminaries
# Load modules that are used in the CliMA project.
# These are general modules not necessarily specific
# to CliMA
using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using DelimitedFiles
using Dierckx

# Load modules specific to CliMA project
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Vtk
#md nothing # hide

# The prognostic equations are conservations laws solved with respect to
# the dynamics and moisture quantities:
# ```math
#  ρ,\;(ρu),\;(ρv),\;(ρw),;(ρe_tot),\;(ρq_tot)
# ```
#md nothing # hide

# Load the MoistThermodynamics and PlanetParameters modules:
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters: R_d, cp_d, grav, cv_d, MSLP, T_0
#md nothing # hide

# Define the ids for each dynamics and moist states:
const _nstate = 6
const _ρ, _U, _V, _W, _E, _QT = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E, QTid = _QT)
const statenames    = ("ρ", "U", "V", "W", "E", "Qtot")
const auxstatenames = ("ax","ay","az","Qliq")

const _nviscstates = 6
const _τ11, _τ22, _τ33, _τ12, _τ13, _τ23 = 1:_nviscstates

const _ngradstates = 3
const _states_for_gradient_transform = (_ρ, _U, _V, _W)
#md nothing # hide

if !@isdefined integration_testing
    const integration_testing =
        parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
    using Random
end

const γ_exact = 7 // 5
const μ_exact = 10

# Domain:
const xmin =    0
const zmin =    0
const ymin =    0
const xmax = 1500 #domain length
const zmax =  150 #domain depth
const ymax = 1500 #domain height (vertical)
const xc   = (xmax + xmin) / 2
const zc   = (zmax + zmin) / 2
const yc   = (ymax + ymin) / 2

# NOTICE: Gravity in the current setting acts along y
# This will be changed to z in the next revision
#md nothing # hide

# -------------------------------------------------------------------------
# Preflux calculation: This function computes parameters required for the
# DG RHS (but not explicitly solved for as a prognostic variable)
# In the case of the dycoms test: the saturation
# adjusted temperature and pressure are such examples. Since we define
# the equation and its arguments here the user is afforded a lot of freedom
# around its behaviour.
# The preflux function interacts with the following
# Modules: NumericalFluxes.jl
# functions: wavespeed, cns_flux!, bcstate!
# -------------------------------------------------------------------------
@inline function preflux(Q,VF, aux, _...)
    γ::eltype(Q) = γ_exact
    gravity::eltype(Q) = grav
    R_gas::eltype(Q) = R_d
    @inbounds ρ, U, V, W, E, QT = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E], Q[_QT]
    ρinv = 1 / ρ

    xvert = aux[_a_y]

    u, v, w = ρinv * U, ρinv * V, ρinv * W
    e_int = (E - (U^2 + V^2+ W^2)/(2*ρ) - ρ * gravity * xvert) / ρ
    qt = QT / ρ
    # Establish the current thermodynamic state using the prognostic variables
    TS           = PhaseEquil(e_int, qt, ρ)
    T            = air_temperature(TS)
    P            = air_pressure(TS) # Test with dry atmosphere
    q_phase_part = PhasePartition(TS)

    (P, u, v, w, ρinv)
    # Preflux returns pressure, 3 velocity components, and 1/ρ
end

# -------------------------------------------------------------------------
# max eigenvalue
@inline function wavespeed(n, Q, aux, t, P, u, v, w, ρinv)
    γ::eltype(Q) = γ_exact
    @inbounds abs(n[1] * u + n[2] * v + n[3] * w) + sqrt(ρinv * γ * P)
end

# -------------------------------------------------------------------------
# ### read sounding
#md #
#md # The sounding file contains the following quantities along a 1D column.
#md # It needs to have the following structure:
#md #
#md # z[m]   theta[K]  q[g/kg]   u[m/s]   v[m/s]   p[Pa]
#md # ...      ...       ...      ...      ...      ...
#md #
#md #
# -------------------------------------------------------------------------
function read_sounding()
    #read in the original squal sounding
    fsounding  = open(joinpath(@__DIR__, "../soundings/sounding_DYCOMS_TEST1.dat"))
    sounding = readdlm(fsounding)
    close(fsounding)
    (nzmax, ncols) = size(sounding)
    if nzmax == 0
        error("SOUNDING ERROR: The Sounding file is empty!")
    end
    return (sounding, nzmax, ncols)
end

# -------------------------------------------------------------------------
# ### Physical Flux (Required)
#md # Here, we define the physical flux function, i.e. the conservative form
#md # of the equations of motion for the prognostic variables ρ, U, V, W, E, QT
#md # $\frac{\partial Q}{\partial t} + \nabla \cdot \boldsymbol{F} = \boldsymbol {S}$
#md # $\boldsymbol{F}$ contains both the viscous and inviscid flux components
#md # and $\boldsymbol{S}$ contains source terms.
#md # Note that the preflux calculation is splatted at the end of the function call
#md # to cns_flux!
# -------------------------------------------------------------------------
cns_flux!(F, Q, VF, aux, t) = cns_flux!(F, Q, VF, aux, t, preflux(Q,VF, aux)...)
@inline function cns_flux!(F, Q, VF, aux, t, P, u, v, w, ρinv)
    @inbounds begin
        ρ, U, V, W, E, QT = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E], Q[_QT]
        # Inviscid contributions
        F[1, _ρ], F[2, _ρ], F[3, _ρ] = U          , V          , W
        F[1, _U], F[2, _U], F[3, _U] = u * U  + P , v * U      , w * U
        F[1, _V], F[2, _V], F[3, _V] = u * V      , v * V + P  , w * V
        F[1, _W], F[2, _W], F[3, _W] = u * W      , v * W      , w * W + P
        F[1, _E], F[2, _E], F[3, _E] = u * (E + P), v * (E + P), w * (E + P)
        F[1, _QT], F[2, _QT], F[3, _QT] = u * QT  , v * QT     , w * QT
        # Stress tensor
        τ11, τ22, τ33 = VF[_τ11], VF[_τ22], VF[_τ33]
        τ12 = τ21 = VF[_τ12]
        τ13 = τ31 = VF[_τ13]
        τ23 = τ32 = VF[_τ23]
        # Viscous contributions
        F[1, _U] -= τ11; F[2, _U] -= τ12; F[3, _U] -= τ13
        F[1, _V] -= τ21; F[2, _V] -= τ22; F[3, _V] -= τ23
        F[1, _W] -= τ31; F[2, _W] -= τ32; F[3, _W] -= τ33
        # Energy dissipation
        F[1, _E] -= u * τ11 + v * τ12 + w * τ13
        F[2, _E] -= u * τ21 + v * τ22 + w * τ23
        F[3, _E] -= u * τ31 + v * τ32 + w * τ33
    end
end

# -------------------------------------------------------------------------
#md # Here we define a function to extract the velocity components from the
#md # prognostic equations (i.e. the momentum and density variables). This
#md # function is not required in general, but provides useful functionality
#md # in some cases.
# -------------------------------------------------------------------------
# Compute the velocity from the state
@inline function velocities!(vel, Q, _...)
    @inbounds begin
        # ordering should match states_for_gradient_transform
        ρ, U, V, W = Q[_ρ], Q[_U], Q[_V], Q[_W]
        ρinv = 1 / ρ
        vel[1], vel[2], vel[3] = ρinv * U, ρinv * V, ρinv * W
    end
end

# -------------------------------------------------------------------------
#md ### Auxiliary Function (Not required)
#md # In this example the auxiliary function is used to store the spatial
#md # coordinates. This may also be used to store variables for which gradients
#md # are needed, but are not available through the prognostic variable
#md # calculations. (An example of this will follow - in the Smagorinsky model,
#md # where a local Richardson number via potential temperature gradient is required)
# -------------------------------------------------------------------------
const _nauxstate = 4
const _a_x, _a_y, _a_z, _a_ql = 1:_nauxstate
@inline function auxiliary_state_initialization!(aux, x, y, z)
    @inbounds begin
        aux[_a_x]  = x
        aux[_a_y]  = y
        aux[_a_z]  = z
    end
end

# -------------------------------------------------------------------------
#md ### Viscous fluxes.
#md # The viscous flux function compute_stresses computes the components of
#md # the velocity gradient tensor, and the corresponding strain rates to
#md # populate the viscous flux array VF. SijSij is calculated in addition
#md # to facilitate implementation of the constant coefficient Smagorinsky model
#md # (pending)
@inline function compute_stresses!(VF, grad_vel, _...)
    μ::eltype(VF) = μ_exact
    @inbounds begin
        dudx, dudy, dudz = grad_vel[1, 1], grad_vel[2, 1], grad_vel[3, 1]
        dvdx, dvdy, dvdz = grad_vel[1, 2], grad_vel[2, 2], grad_vel[3, 2]
        dwdx, dwdy, dwdz = grad_vel[1, 3], grad_vel[2, 3], grad_vel[3, 3]
        # strains
        ϵ11 = dudx
        ϵ22 = dvdy
        ϵ33 = dwdz
        ϵ12 = (dudy + dvdx) / 2
        ϵ13 = (dudz + dwdx) / 2
        ϵ23 = (dvdz + dwdy) / 2
        # --------------------------------------------
        # SMAGORINSKY COEFFICIENT COMPONENTS
        # --------------------------------------------
        SijSij = (ϵ11 + ϵ22 + ϵ33
                  + 2.0 * ϵ12
                  + 2.0 * ϵ13
                  + 2.0 * ϵ23)
        # mod strain rate ϵ ---------------------------
        # deviatoric stresses
        VF[_τ11] = 2μ * (ϵ11 - (ϵ11 + ϵ22 + ϵ33) / 3)
        VF[_τ22] = 2μ * (ϵ22 - (ϵ11 + ϵ22 + ϵ33) / 3)
        VF[_τ33] = 2μ * (ϵ33 - (ϵ11 + ϵ22 + ϵ33) / 3)
        VF[_τ12] = 2μ * ϵ12
        VF[_τ13] = 2μ * ϵ13
        VF[_τ23] = 2μ * ϵ23

    end
end
# -------------------------------------------------------------------------

# -------------------------------------------------------------------------
# generic bc for 2d , 3d
@inline function bcstate!(QP, VFP, auxP, nM, QM, VFM, auxM, bctype, t, PM, uM, vM, wM, ρinvM)
    @inbounds begin
        x, y, z = auxM[_a_x], auxM[_a_y], auxM[_a_z]
        ρM, UM, VM, WM, EM, QTM = QM[_ρ], QM[_U], QM[_V], QM[_W], QM[_E], QM[_QT]
        UnM = nM[1] * UM + nM[2] * VM + nM[3] * WM
        QP[_U] = UM - 2 * nM[1] * UnM
        QP[_V] = VM - 2 * nM[2] * UnM
        QP[_W] = WM - 2 * nM[3] * UnM
        QP[_ρ] = ρM
        QP[_E] = EM
        QP[_QT] = QTM
        VFP .= VFM
        # To calculate PP, uP, vP, wP, ρinvP we use the preflux function
        nothing
        #preflux(QP, auxP, t)
        # Required return from this function is either nothing or preflux with plus state as arguments
    end
end
# -------------------------------------------------------------------------

@inline stresses_boundary_penalty!(VF, _...) = VF.=0

@inline function stresses_penalty!(VF, nM, velM, QM, aM, velP, QP, aP, t)
    @inbounds begin
        n_Δvel = similar(VF, Size(3, 3))
        for j = 1:3, i = 1:3
            n_Δvel[i, j] = nM[i] * (velP[j] - velM[j]) / 2
        end
        compute_stresses!(VF, n_Δvel)
    end
end
# -----------------------------------------------------------------
# --------------DEFINE SOURCES HERE -------------------------------
# -----------------------------------------------------------------
"""
The function source! collects all the individual source terms
associated with a given problem. We do not define sources here,
rather we only call those source terms which are necessary based
on the governing equations.
by terms defined elsewhere
"""
@inline function source!(S,Q,aux,t)
    # Initialise the final block source term
    S .= 0

    # Typically these sources are imported from modules
    @inbounds begin
        source_geopot!(S, Q, aux, t)
    end
end

# Sponge: classical Rayleigh type absorbing layers:
@inline function sponge_rectangular(S, Q, aux)

    U, V, W = Q[_U], Q[_V], Q[_W]
    x, y, z = aux[_a_x], aux[_a_y], aux[_a_z]

    xmin = brickrange[1][1]
    xmax = brickrange[1][end]
    ymin = brickrange[2][1]
    ymax = brickrange[2][end]
    zmin = brickrange[3][1]
    zmax = brickrange[3][end]

    # Define Sponge Boundaries
    xc       = (xmax + xmin)/2
    yc       = (ymax + ymin)/2

    zsponge  = 0.85 * zmax
    xsponger = xmax - 0.15*abs(xmax - xc)
    xspongel = xmin + 0.15*abs(xmin - xc)
    ysponger = ymax - 0.15*abs(ymax - yc)
    yspongel = ymin + 0.15*abs(ymin - yc)

    csxl, csxr  = 0.0, 0.0
    csyl, csyr  = 0.0, 0.0
    ctop        = 0.0

    csx         = 0.0
    csy         = 0.0
    ct          = 1.0

    #x left and right
    #xsl
    if (x <= xspongel)
        csxl = csx * sinpi(1/2 * (x - xspongel)/(xmin - xspongel))^4
    end
    #xsr
    if (x >= xsponger)
        csxr = csx * sinpi(1/2 * (x - xsponger)/(xmax - xsponger))^4
    end
    #y left and right
    #ysl
    if (y <= yspongel)
        csyl = csy * sinpi(1/2 * (y - yspongel)/(ymin - yspongel))^4
    end
    #ysr
    if (y >= ysponger)
        csyr = csy * sinpi(1/2 * (y - ysponger)/(ymay - ysponger))^4
    end

    #Vertical sponge:
    if (z >= zsponge)
        ctop = ct * sinpi(1/2 * (z - zsponge)/(zmax - zsponge))^4
    end

    beta  = 1.0 - (1.0 - ctop) * (1.0 - csxl)*(1.0 - csxr) * (1.0 - csyl)*(1.0 - csyr)
    beta  = min(beta, 1.0)
    alpha = 1.0 - beta

    @inbounds begin
        S[_U] -= beta * U
        S[_V] -= beta * V
        S[_W] -= beta * W
    end

end
#---END SPONGE

"""
Geopotential source term. Gravity forcing applied to the vertical
momentum equation
"""
@inline function source_geopot!(S,Q,aux,t)
    gravity::eltype(Q) = grav
    @inbounds begin
        ρ, U, V, W, E  = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]
        S[_V] += - ρ * gravity
    end
end

"""
Large scale subsidence common to several atmospheric observational
campaigns. In the absence of a GCM to drive the flow we may need to
specify a large scale forcing function.
"""
@inline function source_ls_subsidence!(S,Q,aux,t)
    @inbounds begin
        nothing
    end
end

# ------------------------------------------------------------------
# -------------END DEF SOURCES-------------------------------------#

# initial condition
"""
User-specified. Required.
This function specifies the initial conditions
for the dycoms driver.
"""
function dycoms!(dim, Q, t, x, y, z, _...)
    DFloat         = eltype(Q)
    p0::DFloat      = MSLP
    gravity::DFloat = grav

    # ----------------------------------------------------
    # GET DATA FROM INTERPOLATED ARRAY ONTO VECTORS
    # This driver accepts data in 6 column format
    # ----------------------------------------------------
    (sounding, _, ncols) = read_sounding()

    # WARNING: Not all sounding data is formatted/scaled
    # the same. Care required in assigning array values
    # height theta qv    u     v     pressure
    zinit, tinit, qinit, uinit, vinit, pinit  = sounding[:, 1],
    sounding[:, 2],
    sounding[:, 3],
    sounding[:, 4],
    sounding[:, 5],
    sounding[:, 6]
    #------------------------------------------------------
    # GET SPLINE FUNCTION
    #------------------------------------------------------
    spl_tinit    = Spline1D(zinit, tinit; k=1)
    spl_qinit    = Spline1D(zinit, qinit; k=1)
    spl_uinit    = Spline1D(zinit, uinit; k=1)
    spl_vinit    = Spline1D(zinit, vinit; k=1)
    spl_pinit    = Spline1D(zinit, pinit; k=1)
    # --------------------------------------------------
    # INITIALISE ARRAYS FOR INTERPOLATED VALUES
    # --------------------------------------------------
    xvert          = y

    datat          = spl_tinit(xvert)
    dataq          = spl_qinit(xvert)
    datau          = spl_uinit(xvert)
    datav          = spl_vinit(xvert)
    datap          = spl_pinit(xvert)
    dataq          = dataq * 1.0e-3

    randnum   = rand(1)[1] / 100

    θ_liq = datat
    q_tot = dataq + randnum * dataq
    P     = datap
    T     = air_temperature_from_liquid_ice_pottemp(θ_liq, P, PhasePartition(q_tot))
    ρ     = air_density(T, P)

    # energy definitions
    u, v, w     = 0*datau, 0*datav, 0.0 #geostrophic. TO BE BUILT PROPERLY if Coriolis is considered
    U           = ρ * u
    V           = ρ * v
    W           = ρ * w
    e_kin       = (u^2 + v^2 + w^2) / 2
    e_pot       = gravity * xvert
    e_int       = internal_energy(T, PhasePartition(q_tot))
    E           = ρ * total_energy(e_kin, e_pot, T, PhasePartition(q_tot))

    #Get q_liq and q_ice
    TS           = PhaseEquil(e_int, q_tot, ρ)
    q_phase_part = PhasePartition(TS)

    @inbounds Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E], Q[_QT]= ρ, U, V, W, E, ρ * q_tot

end

function run(mpicomm, dim, Ne, N, timeend, DFloat, dt)

    ArrayType = Array
    # CuArray option (TODO merge new master)

    brickrange = (range(DFloat(xmin), length=Ne[1]+1, DFloat(xmax)),
                  range(DFloat(ymin), length=Ne[2]+1, DFloat(ymax)),
                  range(DFloat(zmin), length=Ne[3]+1, DFloat(zmax)))

    # User defined periodicity in the topl assignment
    # brickrange defines the domain extents
    topl = BrickTopology(mpicomm, brickrange, periodicity=(false,false,false))

    grid = DiscontinuousSpectralElementGrid(topl,
                                            FloatType = DFloat,
                                            DeviceArray = ArrayType,
                                            polynomialorder = N)

    numflux!(x...)   = NumericalFluxes.rusanov!(x..., cns_flux!, wavespeed, preflux)
    numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(x..., cns_flux!, bcstate!, wavespeed, preflux)

    # spacedisc = data needed for evaluating the right-hand side function
    spacedisc = DGBalanceLaw(grid = grid,
                             length_state_vector = _nstate,
                             flux! = cns_flux!,
                             numerical_flux! = numflux!,
                             numerical_boundary_flux! = numbcflux!,
                             number_gradient_states = _ngradstates,
                             states_for_gradient_transform =
                             _states_for_gradient_transform,
                             number_viscous_states = _nviscstates,
                             gradient_transform! = velocities!,
                             viscous_transform! = compute_stresses!,
                             viscous_penalty! = stresses_penalty!,
                             viscous_boundary_penalty! = stresses_boundary_penalty!,
                             auxiliary_state_length = _nauxstate,
                             auxiliary_state_initialization! =
                             auxiliary_state_initialization!,
                             source! = source!)

    # This is a actual state/function that lives on the grid
    initialcondition(Q, x...) = dycoms!(Val(dim), Q, DFloat(0), x...)
    Q = MPIStateArray(spacedisc, initialcondition)

    lsrk = LowStorageRungeKutta(spacedisc, Q; dt = dt, t0 = 0)

    eng0 = norm(Q)
    @info @sprintf """Starting norm(Q₀) = %.16e""" eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(5, mpicomm) do (s=false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf("""Update
                         simtime = %.16e
                         runtime = %s
                         norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                           Dates.format(convert(Dates.DateTime,
                                                Dates.now()-starttime[]),
                                                Dates.dateformat"HH:MM:SS"),energy)
        end
    end

    step = [0]
    mkpath("vtk")
    cbvtk = GenericCallbacks.EveryXSimulationSteps(10) do (init=false)
        outprefix = @sprintf("vtk/cns_%dD_mpirank%04d_step%04d", dim,
                             MPI.Comm_rank(mpicomm), step[1])
        @debug "doing VTK output" outprefix
        writevtk(outprefix, Q, spacedisc, statenames, spacedisc.auxstate, auxstatenames)

        step[1] += 1
        nothing
    end

    # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, ))
    solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

    # Print some end of the simulation information
    engf = norm(Q)
    if integration_testing
        Qe = MPIStateArray(spacedisc,
                           (Q, x...) -> initialcondition!(Val(dim), Q,
                                                          DFloat(timeend), x...))
        engfe = norm(Qe)
        errf = euclidean_distance(Q, Qe)
        @info @sprintf """Finished
        norm(Q)                 = %.16e
        norm(Q) / norm(Q₀)      = %.16e
        norm(Q) - norm(Q₀)      = %.16e
        norm(Q - Qe)            = %.16e
        norm(Q - Qe) / norm(Qe) = %.16e
        """ engf engf/eng0 engf-eng0 errf errf / engfe
    else
        @info @sprintf """Finished
        norm(Q)            = %.16e
        norm(Q) / norm(Q₀) = %.16e
        norm(Q) - norm(Q₀) = %.16e""" engf engf/eng0 engf-eng0
    end
    integration_testing ? errf : (engf / eng0)
end

using Test
let
    MPI.Initialized() || MPI.Init()
    Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())
    mpicomm = MPI.COMM_WORLD
    if MPI.Comm_rank(mpicomm) == 0
        ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
        loglevel = ll == "DEBUG" ? Logging.Debug :
            ll == "WARN"  ? Logging.Warn  :
            ll == "ERROR" ? Logging.Error : Logging.Info
        global_logger(ConsoleLogger(stderr, loglevel))
    else
        global_logger(NullLogger())
    end
    # User defined number of elements
    # User defined timestep estimate
    # User defined simulation end time
    # User defined polynomial order
    numelem = (10,10, 1)
    dt = 1e-2
    timeend = 10
    polynomialorder = 5
    for DFloat in (Float64,) #Float32)
        for dim = 3:3
            engf_eng0 = run(mpicomm, dim, numelem[1:dim], polynomialorder, timeend,
                            DFloat, dt)
        end
    end
end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
module DiffEqBase

using Requires, IterativeSolvers, RecursiveFactorization, ArrayInterface

using StaticArrays # data arrays

using LinearAlgebra, Printf

using DocStringExtensions

using FunctionWrappers: FunctionWrapper

using MuladdMacro, Parameters

using NonlinearSolve

using Reexport

using Statistics

using FastBroadcast: @.., True, False

using PreallocationTools
import PreallocationTools: get_tmp

import Distributions
import ChainRulesCore
import LabelledArrays
import RecursiveArrayTools

import ChainRulesCore: NoTangent
import ZygoteRules

import DEDataArrays: DEDataArray, DEDataVector, DEDataMatrix, copy_fields!
export DEDataArray, DEDataVector, DEDataMatrix

using Setfield

using ForwardDiff

@reexport using SciMLBase

using SciMLBase: @def, DEIntegrator, DEProblem, AbstractDiffEqOperator,
                 AbstractDiffEqLinearOperator, AbstractDiffEqInterpolation,
                 DECallback, AbstractDEOptions, DECache, AbstractContinuousCallback,
                 AbstractDiscreteCallback, AbstractLinearProblem, AbstractNonlinearProblem,
                 AbstractOptimizationProblem, AbstractQuadratureProblem,
                 AbstractSteadyStateProblem, AbstractJumpProblem,
                 AbstractNoiseProblem, AbstractEnsembleProblem, AbstractDynamicalODEProblem,
                 DEAlgorithm, StandardODEProblem,
                 AbstractSensitivityAlgorithm, AbstractODEAlgorithm,
                 AbstractSDEAlgorithm, AbstractDDEAlgorithm, AbstractDAEAlgorithm,
                 AbstractSDDEAlgorithm, AbstractRODEAlgorithm, DAEInitializationAlgorithm,
                 AbstractSteadyStateAlgorithm, AbstractQuadratureAlgorithm,
                 AbstractODEProblem, AbstractDiscreteProblem,
                 AbstractSDEProblem, AbstractRODEProblem, AbstractDDEProblem,
                 AbstractDAEProblem, AbstractSDDEProblem, AbstractBVProblem,
                 AbstractTimeseriesSolution, AbstractNoTimeSolution, numargs,
                 AbstractODEFunction, AbstractSDEFunction, AbstractRODEFunction,
                 AbstractDDEFunction, AbstractSDDEFunction, AbstractDAEFunction,
                 AbstractNonlinearFunction, AbstractEnsembleSolution,
                 AbstractODESolution, AbstractRODESolution, AbstractDAESolution,
                 AbstractDDESolution,
                 EnsembleAlgorithm, EnsembleSolution, EnsembleSummary,
                 QuadratureSolution, NonlinearSolution,
                 TimeGradientWrapper, TimeDerivativeWrapper, UDerivativeWrapper,
                 UJacobianWrapper, ParamJacobianWrapper, JacobianWrapper,
                 check_error!, has_jac, has_tgrad, has_Wfact, has_Wfact_t, has_paramjac,
                 AbstractODEIntegrator, AbstractSDEIntegrator, AbstractRODEIntegrator,
                 AbstractDDEIntegrator, AbstractSDDEIntegrator,
                 AbstractDAEIntegrator, unwrap_cache, has_reinit, reinit!,
                 postamble!, last_step_failed, islinear, has_destats,
                 initialize_dae!, build_solution, solution_new_retcode,
                 solution_new_tslocation, plot_indices,
                 NullParameters, isinplace, AbstractADType, AbstractDiscretization,
                 DISCRETE_OUTOFPLACE_DEFAULT, DISCRETE_INPLACE_DEFAULT,
                 has_analytic, calculate_solution_errors!, AbstractNoiseProcess,
                 has_colorvec, parameterless_type, undefined_exports,
                 is_diagonal_noise, AbstractDiffEqFunction, sensitivity_solution,
                 interp_summary, AbstractHistoryFunction, LinearInterpolation,
                 ConstantInterpolation, HermiteInterpolation, SensitivityInterpolation,
                 NoAD, @add_kwonly,
                 calculate_ensemble_errors, DEFAULT_UPDATE_FUNC, isconstant,
                 DEFAULT_REDUCTION, isautodifferentiable,
                 isadaptive, isdiscrete, has_syms, AbstractAnalyticalSolution,
                 RECOMPILE_BY_DEFAULT

import SciMLBase: solve, init, solve!, __init, __solve, update_coefficients!,
                  update_coefficients, isadaptive, wrapfun_oop, wrapfun_iip,
                  unwrap_fw, promote_tspan, set_u!, set_t!, set_ut!

SciMLBase.isfunctionwrapper(x::FunctionWrapper) = true

"""
$(TYPEDEF)
"""
abstract type Tableau end

"""
$(TYPEDEF)
"""
abstract type ODERKTableau <: Tableau end

"""
$(TYPEDEF)
"""
abstract type DECostFunction end

include("utils.jl")
include("fastpow.jl")
include("destats.jl")
include("calculate_residuals.jl")
include("tableaus.jl")

include("nlsolve/type.jl")
include("nlsolve/newton.jl")
include("nlsolve/functional.jl")
include("nlsolve/utils.jl")
include("callbacks.jl")
include("linear_nonlinear.jl")
include("common_defaults.jl")
include("solve.jl")
include("internal_euler.jl")
include("init.jl")
include("forwarddiff.jl")
include("chainrules.jl")

include("precompile.jl")

"""
$(TYPEDEF)
"""
abstract type AbstractParameterizedFunction{iip} <: AbstractODEFunction{iip} end

"""
$(TYPEDEF)
"""
struct ConvergenceSetup{P,C}
    probs::P
    convergence_axis::C
end

export initialize!, finalize!

export LinSolveFactorize, LinSolveGPUFactorize, DefaultLinSolve, DEFAULT_LINSOLVE,
       LinSolveGMRES, LinSolveCG, LinSolveBiCGStabl, LinSolveChebyshev,
       LinSolveMINRES, LinSolveIterativeSolvers

export NLNewton, NLFunctional, NLAnderson

export SensitivityADPassThrough

end # module

ext: jl


text:
using BinDeps

@BinDeps.setup

liquiddsp = library_dependency("liquiddsp", aliases = ["libliquid.dylib", "libliquid.so", "libliquid.dll"])

url     = "https://github.com/jgaeddert/liquid-dsp/archive/master.tar.gz"
depsdir = BinDeps.depsdir(liquiddsp)
srcdir  = joinpath(depsdir, "src", "liquid-dsp-master")
prefix  = joinpath(depsdir, "usr")

@unix_only  libfilename = "libliquid.so"
@osx_only   libfilename = "libliquid.dylib"

provides(Sources, URI(url), liquiddsp, unpacked_dir="liquid-dsp-master")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(liquiddsp)
        @build_steps begin
            ChangeDirectory(srcdir)
            FileRule( joinpath(prefix, "lib", libfilename),
            @build_steps begin
                `sh bootstrap.sh`
                `sh configure --prefix=$prefix --enable-fftoverride`
                `make`
                `make install`
            end)
        end
    end), liquiddsp, os = :Unix)

@BinDeps.install @compat Dict(:liquiddsp => :liquiddsp)

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
"Which discrete dictionary matches the coefficients of the given dictionary?"
discrete_set(dict::Dictionary) = discrete_set(dict, containertype(dict))
# We dispatch on the container type:
# - it is a vector
discrete_set(dict::Dictionary, ::Type{Vector{T}}) where {T} =
    DiscreteVectorDictionary{T}(length(dict))
# - it is an array
discrete_set(dict::Dictionary, ::Type{Array{T,N}}) where {T,N} =
    DiscreteArrayDictionary{T}(size(dict))
# - it is something else that we don't know about here: return a vector set
discrete_set(dict::Dictionary, ::Type{F}) where {F} =
    DiscreteVectorDictionary{coefficienttype(dict)}(length(dict))

"""
The synthesis operator maps a discrete set of coefficients to a function
in the span of a dictionary.
"""
struct SynthesisOperator <: AbstractOperator
    dictionary  ::  Dictionary
    measure     ::  Union{Measure,Nothing}
end

SynthesisOperator(dict::Dictionary) = hasmeasure(dict) ?
    SynthesisOperator(dict, measure(dict)) :
    SynthesisOperator(dict, nothing)

dictionary(op::SynthesisOperator) = op.dictionary
measure(op::SynthesisOperator) = op.measure

src(op::SynthesisOperator) = discrete_set(dictionary(op))
src_space(op::SynthesisOperator) = Span(src(op))
dest_space(op::SynthesisOperator) = Span(dictionary(op))

(op::SynthesisOperator)(args...) = apply(op, args...)

# We attempt to convert the given coefficients to the container type of the dictionary
apply(op::SynthesisOperator, coef; options...) = _apply(op, coef, dictionary(op))
_apply(op::SynthesisOperator, coef, dict::Dictionary) = Expansion(dict, convert(containertype(dict), coef))

apply(op::SynthesisOperator, expansion::Expansion{D}; options...) where {D <: DiscreteDictionary} =
    apply(op, coefficients(expansion); options...)
apply(op::SynthesisOperator, coef::Expansion; opts...) =
    error("A synthesis operator applies only to coefficients or expansions in discrete sets.")

show(io::IO, mime::MIME"text/plain", op::SynthesisOperator) = composite_show(io, mime, op)
Display.object_parentheses(op::SynthesisOperator) = false
Display.stencil_parentheses(op::SynthesisOperator) = false
Display.displaystencil(op::SynthesisOperator) = _stencil(op, dictionary(op), measure(op))
_stencil(op::SynthesisOperator, dict, ::Nothing) = ["SynthesisOperator(", dict, ")"]
_stencil(op::SynthesisOperator, dict, measure::Measure) = ["SynthesisOperator(", dict, ", ", measure, ")"]

"""
The desynthesis operator does the inverse of the synthesis operator: it maps
an expansion in a dictionary to an expansion in the corresponding discrete set.
"""
struct DesynthesisOperator <: AbstractOperator
    dictionary  ::  Dictionary
end

dictionary(op::DesynthesisOperator) = op.dictionary

src_space(op::DesynthesisOperator) = Span(dictionary(op))
dest(op::DesynthesisOperator) = discrete_set(dictionary(op))
dest_space(op::DesynthesisOperator) = Span(dest(op))

# We accept any expansion and return its coefficients
apply(op::DesynthesisOperator, expansion::Expansion; opts...) = coefficients(expansion)

(*)(op::DesynthesisOperator, expansion::Expansion) = apply(op, expansion)

inv(op::SynthesisOperator) = DesynthesisOperator(dictionary(op))
inv(op::DesynthesisOperator) = SynthesisOperator(dictionary(op))

"Object that can be indexed with a dictionary to construct a synthesis operator."
struct SynthesisOperatorGenerator
end

getindex(op::SynthesisOperatorGenerator, Φ::Dictionary) = SynthesisOperator(Φ)

# the symbol is \bscrT
export 𝓣
"Generate a synthesis operator of a dictionary by indexing 𝓣 with that dictionary."
const 𝓣 = SynthesisOperatorGenerator()


text:
module PardisoMKL

using Libdl

export new_pardiso_handle
export new_iparm
export new_iparm_64
export pardiso_data_type
export pardisoinit!
export pardiso!
export pardiso_64!
export pardiso_getdiag!

include("load_mkl.jl")

deps_jl = joinpath(@__DIR__, "..", "deps", "deps.jl")

if !isfile(deps_jl)
  s = """
  Package PardisoBindings not installed properly.
  """
  error(s)
end

include(deps_jl)

const pardisoinit_sym = Ref{Ptr}()
const pardiso_sym = Ref{Ptr}()
const pardiso_64_sym = Ref{Ptr}()
#const pardiso_getenv_sym = Ref{Ptr}()
#const pardiso_setenv_sym = Ref{Ptr}()
const pardiso_getdiag_sym = Ref{Ptr}()
#const pardiso_export_sym = Ref{Ptr}()
#const pardiso_handle_store_sym = Ref{Ptr}()
#const pardiso_handle_restore_sym = Ref{Ptr}()
#const pardiso_handle_delete_sym = Ref{Ptr}()
const MKL_PARDISO_LOADED = Ref(false)

function __init__()

  libmkl = load_mkl_gcc(mkllibdir,gcclibdir)

  pardisoinit_sym[] = Libdl.dlsym(libmkl,:pardisoinit)
  pardiso_sym[] = Libdl.dlsym(libmkl,:pardiso )
  pardiso_64_sym[] = Libdl.dlsym(libmkl,:pardiso_64 )
  #pardiso_getenv_sym[] = Libdl.dlsym(libmkl,:pardiso_getenv)
  #pardiso_setenv_sym[] = Libdl.dlsym(libmkl,:pardiso_setenv)
  pardiso_getdiag_sym[] = Libdl.dlsym(libmkl,:pardiso_getdiag)
  #pardiso_export_sym[] = Libdl.dlsym(libmkl,:pardiso_export)
  #pardiso_handle_store_sym[] = Libdl.dlsym(libmkl,:pardiso_handle_store)
  #pardiso_handle_restore_sym[] = Libdl.dlsym(libmkl,:pardiso_handle_restore)
  #pardiso_handle_delete_sym[] = Libdl.dlsym(libmkl,:pardiso_handle_delete)

  MKL_PARDISO_LOADED[] = true

end

include("constructors.jl")

include("bindings.jl")

end # module


text:
# these need to be at the top level:

# A dummy RegressionModel type
struct DummyMod <: RegressionModel
    beta::Vector{Float64}
    x::Matrix
    y::Vector
end

## dumb fit method: just copy the x and y input over
StatsBase.fit(::Type{DummyMod}, x::Matrix, y::Vector) =
    DummyMod(collect(1:size(x, 2)), x, y)
StatsBase.response(mod::DummyMod) = mod.y
## dumb coeftable: just prints the "beta" values
StatsBase.coeftable(mod::DummyMod) =
    CoefTable(reshape(mod.beta, (size(mod.beta,1), 1)),
              ["'beta' value"],
              ["" for n in 1:size(mod.x,2)],
              0)
# dumb predict: return values predicted by "beta" and dummy confidence bounds
function StatsBase.predict(mod::DummyMod;
                           interval::Union{Nothing,Symbol}=nothing)
    pred = mod.x * mod.beta
    if interval === nothing
        return pred
    elseif interval === :prediction
        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
    else
        throw(ArgumentError("value not allowed for interval"))
    end
end
function StatsBase.predict(mod::DummyMod, newX::Matrix;
                           interval::Union{Nothing,Symbol}=nothing)
    pred = newX * mod.beta
    if interval === nothing
        return pred
    elseif interval === :prediction
        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
    else
        throw(ArgumentError("value not allowed for interval"))
    end
end

# A dummy RegressionModel type that does not support intercept
struct DummyModNoIntercept <: RegressionModel
    beta::Vector{Float64}
    x::Matrix
    y::Vector
end

StatsModels.drop_intercept(::Type{DummyModNoIntercept}) = true

## dumb fit method: just copy the x and y input over
StatsBase.fit(::Type{DummyModNoIntercept}, x::Matrix, y::Vector) =
    DummyModNoIntercept(collect(1:size(x, 2)), x, y)
StatsBase.response(mod::DummyModNoIntercept) = mod.y
## dumb coeftable: just prints the "beta" values
StatsBase.coeftable(mod::DummyModNoIntercept) =
    CoefTable(reshape(mod.beta, (size(mod.beta,1), 1)),
              ["'beta' value"],
              ["" for n in 1:size(mod.x,2)],
              0)
# dumb predict: return values predicted by "beta" and dummy confidence bounds
function StatsBase.predict(mod::DummyModNoIntercept;
                           interval::Union{Nothing,Symbol}=nothing)
    pred = mod.x * mod.beta
    if interval === nothing
        return pred
    elseif interval === :prediction
        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
    else
        throw(ArgumentError("value not allowed for interval"))
    end
end
function StatsBase.predict(mod::DummyModNoIntercept, newX::Matrix;
                           interval::Union{Nothing,Symbol}=nothing)
    pred = newX * mod.beta
    if interval === nothing
        return pred
    elseif interval === :prediction
        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
    else
        throw(ArgumentError("value not allowed for interval"))
    end
end

## Another dummy model type to test fall-through show method
struct DummyModTwo <: RegressionModel
    msg::String
end

StatsBase.fit(::Type{DummyModTwo}, ::Matrix, ::Vector) = DummyModTwo("hello!")
Base.show(io::IO, m::DummyModTwo) = println(io, m.msg)

@testset "stat model types" begin

    ## Test fitting
    d = DataFrame(y = 1:4,
                  x1 = allowmissing(5:8),
                  x2 = 9:12,
                  x3 = 13:16,
                  x4 = 17:20)

    d.x1p = categorical(d.x1)

    f = @formula(y ~ x1 * x2)
    m = fit(DummyMod, f, d)
    @test response(m) == Array(d.y)

    ## coefnames delegated to model frame by default
    @test coefnames(m) == coefnames(ModelFrame(f, d)) == ["(Intercept)", "x1", "x2", "x1 & x2"]

    ## test prediction method
    ## vanilla
    @test predict(m) == [ ones(size(d,1)) Array(d.x1) Array(d.x2) Array(d.x1).*Array(d.x2) ] * collect(1:4)

    ## new data from matrix
    mm = ModelMatrix(ModelFrame(f, d))
    p = predict(m, mm.m)
    @test p == mm.m * collect(1:4)
    p2 = predict(m, mm.m, interval=:prediction)
    @test p2 isa NamedTuple
    @test p2.prediction == p
    @test p2.lower == p .- 1
    @test p2.upper == p .+ 1

    ## new data from DataFrame (via ModelMatrix)
    @test predict(m, d) == p
    p3 = predict(m, d, interval=:prediction)
    @test p3 isa DataFrame
    @test p3.prediction == p
    @test p3.lower == p .- 1
    @test p3.upper == p .+ 1

    d2 = deepcopy(d)
    d2[3, :x1] = missing
    @test length(predict(m, d2)) == 4

    ## test copying of names from Terms to CoefTable
    ct = coeftable(m)
    @test ct.rownms == ["(Intercept)", "x1", "x2", "x1 & x2"]

    ## show with coeftable defined
    io = IOBuffer()
    show(io, m)

    ## with categorical variables
    f2 = @formula(y ~ x1p)
    m2 = fit(DummyMod, f2, d)

    @test coeftable(m2).rownms == ["(Intercept)", "x1p: 6", "x1p: 7", "x1p: 8"]

    ## predict w/ new data missing levels
    @test predict(m2, d[2:4, :]) == predict(m2)[2:4]

    ## predict w/ new data with _extra_ levels (throws an error)
    d3 = deepcopy(d)
    d3[1, :x1] = 0
    d3.x1p = categorical(d3.x1)
    # TODO: check for level mismatch earlier...this throws a KeyError when it
    # goes to do the lookup in the contrasts matrix from the previously
    # generated categorical term.
    @test_throws KeyError predict(m2, d3)
    # @test_throws ArgumentError predict(m2, d3)

    ## predict with dataframe that doesn't have the dependent variable
    d4 = deepcopy(d)
    select!(d4, Not(:y))
    @test predict(m, d4) == predict(m, d)

    ## attempting to fit with d4 should fail since it doesn't have :y
    @test_throws ErrorException fit(DummyMod, f, d4)

    ## fit with contrasts specified
    d.x2p = categorical(d.x2)
    f3 = @formula(y ~ x1p + x2p)
    m3 = fit(DummyMod, f3, d)
    fit(DummyMod, f3, d, contrasts = Dict(:x1p => EffectsCoding()))
    fit(DummyMod, f3, d, contrasts = Dict(:x1p => EffectsCoding(),
                                          :x2p => DummyCoding()))
    @test_throws Exception fit(DummyMod, f3, d, contrasts = Dict(:x1p => EffectsCoding(),
                                                                 :x2p => 1))

    f1 = @formula(y ~ 1 + x1 * x2)
    f2 = @formula(y ~ 0 + x1 * x2)
    f3 = @formula(y ~ x1 * x2)
    @test_throws ArgumentError m1 = fit(DummyModNoIntercept, f1, d)
    m2 = fit(DummyModNoIntercept, f2, d)
    m3 = fit(DummyModNoIntercept, f3, d)
    ct2 = coeftable(m2)
    ct3 = coeftable(m3)
    @test ct3.rownms == ct2.rownms == ["x1", "x2", "x1 & x2"]
    @test predict(m2, d[2:4, :]) == predict(m2)[2:4]
    @test predict(m3, d[2:4, :]) == predict(m3)[2:4]

    f1 = @formula(y ~ 1 + x1p)
    f2 = @formula(y ~ 0 + x1p)
    f3 = @formula(y ~ x1p)
    @test_throws ArgumentError m1 = fit(DummyModNoIntercept, f1, d)
    m2 = fit(DummyModNoIntercept, f2, d)
    m3 = fit(DummyModNoIntercept, f3, d)
    ct2 = coeftable(m2)
    ct3 = coeftable(m3)
    @test ct2.rownms == ct3.rownms == ["x1p: 6", "x1p: 7", "x1p: 8"]
    m4 = fit(DummyModNoIntercept, f3, d, contrasts = Dict(:x1p => EffectsCoding()))
    @test predict(m2, d[2:4, :]) == predict(m2)[2:4]
    @test predict(m3, d[2:4, :]) == predict(m3)[2:4]
    @test predict(m4, d[2:4, :]) == predict(m4)[2:4]

    m2 = fit(DummyModTwo, f, d)
    show(io, m2)

end


text:
# This file will load some things in the REPL for you automatically
# Run it using `julia -L repl_init.jl`

using Pkg

Pkg.activate(".")

using PhysicsEngine

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
import ScikitLearnBase: fit!

export fit!, fit_line_search!

function fit!(model::MatFacModel, A::AbstractMatrix; method="adagrad", kwargs...)

    if method == "adagrad"
        fit_adagrad!(model, A; kwargs...)
    elseif method == "line_search"
        fit_line_search!(model, A; kwargs...)
    end

end

function compute_Z!(Z, X, Y, C, beta, mu, theta)
    Z .= transpose(X)*Y
    Z .+= C*beta
    Z .+= mu
    Z .+= transpose(theta)
end

function adagrad_update!(value, grad, G, lr)

   # complete value update
   G .+= (grad.^2.0)

   value .-= grad .* (lr ./ sqrt.(G))

   return value, G

end

function compute_factor_importances(X, Y)
    X_norms = sqrt.(sum(X.^2, dims=2))
    Y_norms = sqrt.(sum(Y.^2, dims=2))

    importances = X_norms .* Y_norms

    return reshape(importances, (size(importances,1),))
end

function normalize_factors!(X,Y; target_size=nothing)

    N = size(Y, 2)

    if target_size == nothing
        target_size = N/10.0 # This is kind of arbitrary
    end

    Y_norms = sqrt.(sum(Y.^2, dims=2))

    multipliers = target_size ./ Y_norms
    Y .*= multipliers
    X ./= multipliers
end

function fit_adagrad!(model::MatFacModel, A::AbstractMatrix;
                      instance_covariates::Union{Nothing,AbstractMatrix}=nothing,
                      inst_reg_weight::Real=1.0, feat_reg_weight::Real=1.0,
                      a_0_tau::Real=1.0, b_0_tau::Real=1e-3,
                      max_iter::Integer=1000,
                      lr::Real=0.01, eps::Real=1e-8,
                      lr_halflife=Inf,
                      abs_tol::Real=1e-3, rel_tol::Real=1e-7,
                      loss_iter::Integer=10,
                      normalize_Y::Bool=true)

    # Setup
    iter = 0
    cur_loss = Inf
    new_loss = Inf

    M = size(A,1)
    N = size(A,2)
    K = size(model.X, 1)

    @assert size(model.X,1) == size(model.Y,1)

    lr = Float32(lr)
    lr_shrink = Float32(1.0)
    if lr_halflife != Inf
        lr_shrink = Float32(0.5^(1.0/lr_halflife))
    end

    if instance_covariates == nothing
        instance_covariates = zeros(M,1)
        instance_covariate_coeff = zeros(1,N)
    else
        @assert size(instance_covariates,1) == M
        instance_covariate_coeff = randn(size(instance_covariates,2), N)
    end

    covariate_dim = size(instance_covariates, 2)

    # Move data and model to the GPU.
    # Float32 is sufficiently precise.
    X_d = CuArray{Float32}(model.X)
    Y_d = CuArray{Float32}(model.Y)
    beta_d = CuArray{Float32}(instance_covariate_coeff)
    beta_reg_d = CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}(feat_reg_weight .* model.instance_covariate_coeff_reg)

    mu_d = CuArray{Float32}(model.instance_offset)
    mu_reg_d = CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}(inst_reg_weight .* model.instance_offset_reg)
    theta_d = CuArray{Float32}(model.feature_offset)
    theta_reg_d = CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}(feat_reg_weight .* model.feature_offset_reg)
    tau_d = CuArray{Float32}(model.feature_precision)

    A_d = CuArray{Float32}(A)
    C_d = CuArray{Float32}(instance_covariates)
    Z_d = CUDA.zeros(Float32, M, N)

    # Some bookkeeping for missing values.
    obs_mask = (!isnan).(A_d)
    missing_mask = (isnan.(A_d) .* Float32(0.5))
    # Convert NaNs in the data --> zeros so CuArray arithmetic works
    mask_func(x) = isnan(x) ? Float32(0.5) : Float32(x)
    map!(mask_func, A_d, A_d)

    # Bookkeeping for the different noise models/loss functions
    ql_idx = CuVector{Int64}(findall(typeof.(model.losses) .== QuadLoss))
    ll_idx = CuVector{Int64}(findall(typeof.(model.losses) .== LogisticLoss))
    pl_idx = CuVector{Int64}(findall(typeof.(model.losses) .== PoissonLoss))

    Z_ql_view = view(Z_d, :, ql_idx)
    Z_ll_view = view(Z_d, :, ll_idx)
    Z_pl_view = view(Z_d, :, pl_idx)
    A_ql_view = view(A_d, :, ql_idx)
    A_ll_view = view(A_d, :, ll_idx)
    A_pl_view = view(A_d, :, pl_idx)

    tau_ql_view = view(tau_d, ql_idx)

    # Convert regularizers to CuSparse matrices
    inst_reg_mats_d = [CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}(mat .* inst_reg_weight) for mat in model.instance_reg_mats]
    feat_reg_mats_d = [CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}(mat .* feat_reg_weight) for mat in model.feature_reg_mats]

    # Arrays for holding gradients and summed gradients
    grad_X = CUDA.zeros(Float32, (K, M))
    G_X = CuArray{Float32}(fill(eps, K, M))

    grad_mu = CUDA.zeros(Float32, M)
    G_mu = CuArray{Float32}(fill(eps, M))

    grad_Y = CUDA.zeros(Float32, (K, N))
    G_Y = CuArray{Float32}(fill(eps, K, N))

    grad_beta = CUDA.zeros(Float32, (covariate_dim, N))
    G_beta = CuArray{Float32}(fill(eps, covariate_dim, N))

    grad_theta = CUDA.zeros(Float32, N)
    G_theta = CuArray{Float32}(fill(eps, N))

    while iter < max_iter

        ################################
        # Update Row quantities (X, mu)

        compute_Z!(Z_d, X_d, Y_d, C_d, beta_d, mu_d, theta_d)
        Z_d .*= obs_mask
        Z_d .+= missing_mask

        # update X (NOTE: this stores Z-A in Z_d)
        compute_grad_X!(grad_X, X_d, Y_d, Z_d,
                        Z_ql_view, Z_ll_view, Z_pl_view,
                        A_ql_view, A_ll_view, A_pl_view,
                        inst_reg_mats_d)
        X_d, G_X = adagrad_update!(X_d, grad_X, G_X, lr)

        # Update mu (Reuse the Z-A stored in Z_d)
        compute_grad_mu!(grad_mu, mu_d, Z_d, tau_d, mu_reg_d)
        mu_d, G_mu = adagrad_update!(mu_d, grad_mu, G_mu, lr)

        ############################
        # Update Column quantities (Y, beta, theta, tau)

        compute_Z!(Z_d, X_d, Y_d, C_d, beta_d, mu_d, theta_d)
        Z_d .*= obs_mask
        Z_d .+= missing_mask

        # Update Y
        compute_grad_Y!(grad_Y, X_d, Y_d, Z_d,
                        Z_ql_view, Z_ll_view, Z_pl_view,
                        A_ql_view, A_ll_view, A_pl_view,
                        feat_reg_mats_d)
        Y_d, G_Y = adagrad_update!(Y_d, grad_Y, G_Y, lr)

        # Update beta
        compute_grad_beta!(grad_beta, beta_d, C_d, Z_d, tau_d, beta_reg_d)
        beta_d, G_beta = adagrad_update!(beta_d, grad_beta, G_beta, lr)

        # Update theta
        compute_grad_theta!(grad_theta, theta_d, Z_d, tau_d, theta_reg_d)
        theta_d, G_theta = adagrad_update!(theta_d, grad_theta, G_theta, lr)

        # Update tau (closed-form exact update)
        update_tau!(tau_ql_view, Z_ql_view, a_0_tau, b_0_tau)

        ############################
        # Every so often, compute the loss

        iter += 1
        print_str = "Iteration: $iter"

        if (iter % loss_iter == 0)

            compute_Z!(Z_d, X_d, Y_d, C_d, beta_d, mu_d, theta_d)
            Z_d .*= obs_mask
            Z_d .+= missing_mask

            new_loss = compute_loss!(X_d, Y_d, Z_ql_view, Z_ll_view, Z_pl_view,
                                               A_ql_view, A_ll_view, A_pl_view,
                                               beta_d, beta_reg_d,
                                               mu_d, mu_reg_d, theta_d, theta_reg_d,
                                               tau_ql_view, a_0_tau, b_0_tau,
                                               inst_reg_mats_d, feat_reg_mats_d;
                                               verbose=true)
            print_str = string(print_str, "\tLoss: $new_loss")
            println(print_str)

            if abs(new_loss - cur_loss) < abs_tol
                println(string("Absolute change <", abs_tol, ". Terminating."))
                break
            end
            if abs((new_loss - cur_loss)/cur_loss) < rel_tol
                println(string("Relative change <", rel_tol, ". Terminating."))
                break
            end
            cur_loss = new_loss
        end

        # Decrease the learning rate
        lr *= lr_shrink

    end # while

    importances = compute_factor_importances(X_d, Y_d)

    if normalize_Y
        normalize_factors!(X_d, Y_d)
    end

    # Move model parameters back to CPU
    model.X = Array{Float32}(X_d)
    model.Y = Array{Float32}(Y_d)
    model.instance_covariate_coeff = Array{Float32}(beta_d)
    model.instance_offset = Array{Float32}(mu_d)
    model.feature_offset = Array{Float32}(theta_d)
    model.feature_precision = Array{Float32}(tau_d)
    model.factor_importances = Array{Float32}(importances)

    return model
end

# Some precompilation
precompile(compute_Z!, (CuArray{Float32,2}, CuArray{Float32,2},
                        CuArray{Float32,2}, CuArray{Float32,2},
                        CuArray{Float32,2},
                        CuArray{Float32,1}, CuArray{Float32,1}))

for t in (CuArray{Float32,2}, CuArray{Float32,1})
    precompile(adagrad_update!, (t, t, t, Float32))
end

precompile(normalize_factors!, (CuArray{Float32}, CuArray{Float32}))
precompile(compute_factor_importances, (CuArray{Float32}, CuArray{Float32}))


text:
#### Pauli

mutable struct Pauli{T, B}
    x::B
    z::B
    phase::T
end

Base.length(p::Pauli) = length(p.x)

"""
    symplectic_to_char(x::Integer, z::Integer)

Return a character in the set `IXYZ` given by the symplectic representation `x, z`,
where `x` and `z` are typically `Bool`s.
"""
function symplectic_to_char(x::Integer, z::Integer)
    (x & z == 1) && return 'Y'
    x == 1 && return 'X'
    z == 1 && return 'Z'
    return 'I'
end

function Base.show(io::IO, p::Pauli)
    phase_strings = ("", "-i", "-", "i")
    print(io, phase_strings[p.phase + 1])
    for (x, z) in zip(p.x, p.z)
        print(io, symplectic_to_char(x, z))
    end
end

Pauli(n::Integer) = Pauli(Vector{Bool}(undef, n), Vector{Bool}(undef, n), 0)
#Pauli(n::Integer) = Pauli(BitVector(undef, n), BitVector(undef, n), 0)

function Base.setindex!(p::Pauli, (x, z), i::Integer)
    p.x[i] = x
    p.z[i] = z
    return (x, z)
end

function symplectic(c::AbstractChar)
    c == 'I' && return (false, false)
    c == 'X' && return (true, false)
    c == 'Y' && return (true, true)
    c == 'Z' && return (false, true)
    throw(ArgumentError("Expecting one of IXYZ"))
end

function symplectic(s::AbstractString)
    p = Pauli(length(s))
    for (i, c) in enumerate(s)
        p[i] = symplectic(c)
    end
    return p
end

function Pauli(s::AbstractString)
    phase = 0
    offset = 0
    if s[1] == '+'
        if s[2] == 'i'
            phase = 3
            offset = 2
        else
            phase = 0
            offset = 1
        end
    elseif s[1] == '-'
        if s[2] == 'i'
            phase = 1
            offset = 2
        else
            phase = 2
            offset = 1
        end
    elseif s[1] == 'i'
        phase = 3
        offset = 1
    end
    if offset == 0 # There is an overhead for SubString :(
        p = symplectic(s)
    else
        p = symplectic(SubString(s, offset + 1, length(s)))
    end
    p.phase = phase
    return p
end


text:
module TestGLM

# using Revise
using Test
using MLJBase

import MLJModels
import GLM # MLJModels.GLM_ now available for loading
using MLJModels.GLM_

task = load_boston()
X, y = X_and_y(task)

train, test = partition(eachindex(y), 0.7)

atom_ols = OLSRegressor()

fitresult, cache, report = MLJBase.fit(atom_ols, 1, MLJBase.selectrows(X, train), y[train])
p = predict(atom_ols, fitresult, MLJBase.selectrows(X, test))

typeof(fitresult) isa  MLJBase.fitresult_type(atom_ols)

# ols = machine(atom_ols, MLJBase.selectrows(X, train), y[train])
# fit!(ols)

# p = predict(ols, MLJBase.selectrows(X, test))

# hand made regression to compare

Xa = MLJBase.matrix(X) # convert(Matrix{Float64}, X)
Xa1 = hcat(Xa, ones(size(Xa, 1)))
coefs = Xa1[train, :] \ y[train]

p2 = Xa1[test, :] * coefs

@test p ≈ p2

info(atom_ols)

end # module
true

max_stars_repo_name: ablaom/MLJModels

max_issues_repo_name: ablaom/MLJModels

max_forks_repo_name: ablaom/MLJModels


text:
abstract type AbstractMsgPackFormat end

#####
##### `int` family
#####

# `fixint`

struct IntFixPositiveFormat <: AbstractMsgPackFormat
    byte::UInt8
end

struct IntFixNegativeFormat <: AbstractMsgPackFormat
    byte::Int8
end

Base.@pure magic_byte_min(::Type{IntFixPositiveFormat}) = 0x00
Base.@pure magic_byte_max(::Type{IntFixPositiveFormat}) = 0x7f
Base.@pure magic_byte_min(::Type{IntFixNegativeFormat}) = 0xe0
Base.@pure magic_byte_max(::Type{IntFixNegativeFormat}) = 0xff

# `uint`

struct UInt8Format <: AbstractMsgPackFormat end
struct UInt16Format <: AbstractMsgPackFormat end
struct UInt32Format <: AbstractMsgPackFormat end
struct UInt64Format <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{UInt8Format}) = 0xcc
Base.@pure magic_byte(::Type{UInt16Format}) = 0xcd
Base.@pure magic_byte(::Type{UInt32Format}) = 0xce
Base.@pure magic_byte(::Type{UInt64Format}) = 0xcf

# `int`

struct Int8Format <: AbstractMsgPackFormat end
struct Int16Format <: AbstractMsgPackFormat end
struct Int32Format <: AbstractMsgPackFormat end
struct Int64Format <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{Int64Format}) = 0xd3
Base.@pure magic_byte(::Type{Int32Format}) = 0xd2
Base.@pure magic_byte(::Type{Int16Format}) = 0xd1
Base.@pure magic_byte(::Type{Int8Format}) = 0xd0

#####
##### `nil` family
#####

struct NilFormat <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{NilFormat}) = 0xc0

#####
##### `bool` family
#####

struct FalseFormat <: AbstractMsgPackFormat end
struct TrueFormat <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{FalseFormat}) = 0xc2
Base.@pure magic_byte(::Type{TrueFormat}) = 0xc3

#####
##### `float` family
#####

struct Float32Format <: AbstractMsgPackFormat end
struct Float64Format <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{Float32Format}) = 0xca
Base.@pure magic_byte(::Type{Float64Format}) = 0xcb

#####
##### `str` family
#####

struct StrFixFormat <: AbstractMsgPackFormat
    byte::UInt8
end

struct Str8Format <: AbstractMsgPackFormat end
struct Str16Format <: AbstractMsgPackFormat end
struct Str32Format <: AbstractMsgPackFormat end

Base.@pure magic_byte_min(::Type{StrFixFormat}) = 0xa0
Base.@pure magic_byte_max(::Type{StrFixFormat}) = 0xbf
Base.@pure magic_byte(::Type{Str8Format}) = 0xd9
Base.@pure magic_byte(::Type{Str16Format}) = 0xda
Base.@pure magic_byte(::Type{Str32Format}) = 0xdb

#####
##### `bin` family
#####

struct Bin8Format <: AbstractMsgPackFormat end
struct Bin16Format <: AbstractMsgPackFormat end
struct Bin32Format <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{Bin8Format}) = 0xc4
Base.@pure magic_byte(::Type{Bin16Format}) = 0xc5
Base.@pure magic_byte(::Type{Bin32Format}) = 0xc6

#####
##### `array` family
#####

struct ArrayFixFormat <: AbstractMsgPackFormat
    byte::UInt8
end

struct Array16Format <: AbstractMsgPackFormat end
struct Array32Format <: AbstractMsgPackFormat end

Base.@pure magic_byte_min(::Type{ArrayFixFormat}) = 0x90
Base.@pure magic_byte_max(::Type{ArrayFixFormat}) = 0x9f
Base.@pure magic_byte(::Type{Array16Format}) = 0xdc
Base.@pure magic_byte(::Type{Array32Format}) = 0xdd

#####
##### `map` family
#####

struct MapFixFormat <: AbstractMsgPackFormat
    byte::UInt8
end

struct Map16Format <: AbstractMsgPackFormat end
struct Map32Format <: AbstractMsgPackFormat end

Base.@pure magic_byte_min(::Type{MapFixFormat}) = 0x80
Base.@pure magic_byte_max(::Type{MapFixFormat}) = 0x8f
Base.@pure magic_byte(::Type{Map16Format}) = 0xde
Base.@pure magic_byte(::Type{Map32Format}) = 0xdf

#####
##### `ext` family
#####

struct ExtFix1Format <: AbstractMsgPackFormat end
struct ExtFix2Format <: AbstractMsgPackFormat end
struct ExtFix4Format <: AbstractMsgPackFormat end
struct ExtFix8Format <: AbstractMsgPackFormat end
struct ExtFix16Format <: AbstractMsgPackFormat end
struct Ext8Format <: AbstractMsgPackFormat end
struct Ext16Format <: AbstractMsgPackFormat end
struct Ext32Format <: AbstractMsgPackFormat end

Base.@pure magic_byte(::Type{ExtFix1Format}) = 0xd4
Base.@pure magic_byte(::Type{ExtFix2Format}) = 0xd5
Base.@pure magic_byte(::Type{ExtFix4Format}) = 0xd6
Base.@pure magic_byte(::Type{ExtFix8Format}) = 0xd7
Base.@pure magic_byte(::Type{ExtFix16Format}) = 0xd8
Base.@pure magic_byte(::Type{Ext8Format}) = 0xc7
Base.@pure magic_byte(::Type{Ext16Format}) = 0xc8
Base.@pure magic_byte(::Type{Ext32Format}) = 0xc9


text:
@testset "FixedHomotopy" begin
    PolyImpl.@polyvar x y a
    H = StraightLineHomotopy([y^2-x], [x^2+3.0y])

    H1 = FixedHomotopy(H, 1.0)
    H0 = FixedHomotopy(H, 0.0)
    HH = StraightLineHomotopy(H1, H0)

    @test typeof(H1)<:FixedHomotopy{Float64, typeof(H), Float64}

    @test evaluate(H1, [2.0, 1.0]) == evaluate(H, [2.0, 1.0], 1.0)
    @test H1([2.0, 1.0]) == H([2.0, 1.0], 1.0)
    @test H0([2.0, 1.0]) == H([2.0, 1.0], 0.0)
    @test HH([2.0, 1.0], 1.0) == H([2.0, 1.0], 1.0)
    @test HH([2.0, 1.0], 0.0) == H([2.0, 1.0], 0.0)

    J_H = differentiate(H)
    J_H1 = differentiate(H1)
    J_H0 = differentiate(H0)

    @test J_H0([1.0, 1.0]) == J_H([1.0, 1.0], 0.0)
    @test J_H1([1.0, 1.0]) == J_H([1.0, 1.0], 1.0)

    @test degrees(H1) == [2]

    H = FixedHomotopy(StraightLineHomotopy([y^2-x, 3x + y], [x^2+3y, 2y-3x]), 0.34)

    @test homogenized(H) == false
    K = homogenize(H)
    @test ishomogenous(H) == false
    @test homogenized(K)
    @test ishomogenous(K)
    @test nvariables(K) == 3
end


text:
#gtk_editdlg.jl: Edit dialog
#-------------------------------------------------------------------------------

#==Types
===============================================================================#

mutable struct ExploreEditDlg <: AbstractDialog
	data::FieldViewEditData
	explore::ExploreWnd #Reference
	sel::ExploreSelection #Store desired selection (copied so it doesn't change)
	wnd::_Gtk.Window
	ent_fieldinfo::_Gtk.Entry
	frame_fieldview::_Gtk.Frame
	frame_shortdescr::_Gtk.Frame
	tv_shortdescr::_Gtk.TextView
	frame_descr::_Gtk.Frame
	tv_descr::_Gtk.TextView
	btn_ok::_Gtk.Button
	btn_apply::_Gtk.Button
end

#==Helper functions
===============================================================================#
function _apply(dlg::ExploreEditDlg)
	tb = dlg.tv_shortdescr.buffer[_Gtk.TextBuffer]
		dlg.data.shortdescr = get_gtk_property(tb, "text", String)
	tb = dlg.tv_descr.buffer[_Gtk.TextBuffer]
		dlg.data.descr = get_gtk_property(tb, "text", String)
	try
		_write(dlg.data, dlg.explore.db, dlg.sel)
	catch e
		@warn(e)
		return false
	end
	return true
end

#==Callback functions
===============================================================================#
@guarded function cb_wnddestroyed(w::Ptr{Gtk.GObject}, dlg::ExploreEditDlg)
	dlg.explore.editdlg = NoDialog()
	refresh!(dlg.explore)
	return #Known value
end
@guarded function cb_ok_clicked(w::Ptr{Gtk.GObject}, dlg::ExploreEditDlg)
	if _apply(dlg)
		close(dlg) #Only close if succeeds
	end
	return #Known value
end
@guarded function cb_cancel_clicked(w::Ptr{Gtk.GObject}, dlg::ExploreEditDlg)
	close(dlg)
	return #Known value
end
@guarded function cb_apply_clicked(w::Ptr{Gtk.GObject}, dlg::ExploreEditDlg)
	_apply(dlg)
	refresh!(dlg.explore)
	return #Known value
end

#==Constructors
===============================================================================#
function ExploreEditDlg(explore::ExploreWnd)
	border_width = 10
	_spacing = 10

	vbox = _Gtk.Box(true, _spacing) #Main vbox for all widgets
		set_gtk_property!(vbox, "border-width", border_width)

	#Textbox explaining which field we are changing:
	ent_fieldinfo = _Gtk.Entry()
		push!(vbox, ent_fieldinfo)
		set_gtk_property!(ent_fieldinfo, "editable", false)

	frame_fieldview = _Gtk.Frame()
		push!(vbox, frame_fieldview)
	vbox_fieldview = _Gtk.Box(true, 0)
		push!(frame_fieldview, vbox_fieldview)

	#Short description
	frame_shortdescr = _Gtk.Frame()
		push!(vbox_fieldview, frame_shortdescr)
		set_gtk_property!(frame_shortdescr, "label", "Description (courte):")
		set_gtk_property!(frame_shortdescr, "border-width", border_width)
	vbox_pad = _Gtk.Box(true, 0)
		push!(frame_shortdescr, vbox_pad)
		set_gtk_property!(vbox_pad, "border-width", border_width)
	tv_shortdescr = _Gtk.TextView()
		push!(vbox_pad, tv_shortdescr)
		set_gtk_property!(tv_shortdescr, "editable", true)
		set_gtk_property!(tv_shortdescr, "wrap-mode", GtkWrapMode.WORD_CHAR)

	#Long description
	frame_descr = _Gtk.Frame()
		push!(vbox_fieldview, frame_descr)
		set_gtk_property!(frame_descr, "label", "Description:")
		set_gtk_property!(frame_descr, "border-width", border_width)
	vbox_pad = _Gtk.Box(true, 0)
		push!(frame_descr, vbox_pad)
		set_gtk_property!(vbox_pad, "border-width", border_width)
	tv_descr = _Gtk.TextView()
		push!(vbox_pad, tv_descr)
		set_gtk_property!(tv_descr, "editable", true)
		set_gtk_property!(tv_descr, "wrap-mode", GtkWrapMode.WORD_CHAR)
#		GAccessor.size_request(tv_descr, -1, 100)
		set_gtk_property!(tv_descr, "vexpand", true)
		set_gtk_property!(tv_descr, "vexpand-set", true)

	btnbox = _Gtk.Box(false, _spacing)
		push!(vbox, btnbox)
	filler = _Gtk.Box(true, 0)
		push!(btnbox, filler)
		set_gtk_property!(filler, "hexpand", true)
		set_gtk_property!(filler, "hexpand-set", true)
	btn_ok = _Gtk.Button("Ok")
		push!(btnbox, btn_ok)
	btn_apply = _Gtk.Button("Appliquer")
		push!(btnbox, btn_apply)
	btn_cancel = _Gtk.Button("Annuler")
		push!(btnbox, btn_cancel)
#	gtk_container_child_set_property #Wrapped in GAccessor.property:
#	GAccessor.property(btnbox, btn_ok, "pack-type", GtkPackType.END) #Crashes; expects ptr for value
#	GAccessor.property(btnbox, btn_cancel, "pack-type", GtkPackType.END)

	wnd = Gtk.Window(vbox, "", 400, 300)#, true)
		set_gtk_property!(wnd, "title", "ExploreCMEO")

	dlg = ExploreEditDlg(FieldViewEditData(),
		explore, ExploreSelection(), wnd,
		ent_fieldinfo, frame_fieldview,
		frame_shortdescr, tv_shortdescr, frame_descr, tv_descr, btn_ok, btn_apply
	)

	signal_connect(cb_wnddestroyed, wnd, "destroy", Nothing, (), false, dlg)
	signal_connect(cb_ok_clicked, btn_ok, "clicked", Nothing, (), false, dlg)
	signal_connect(cb_cancel_clicked, btn_cancel, "clicked", Nothing, (), false, dlg)
	signal_connect(cb_apply_clicked, btn_apply, "clicked", Nothing, (), false, dlg)

	return dlg
end

#==_show() algorithm
===============================================================================#
function _show(dlg::ExploreEditDlg)
	fieldview = _activefieldview(dlg.explore)
	dlg.sel = deepcopy(dlg.explore.sel) #Get new selection
	success = true
	try
		dlg.data = FieldViewEditData(fieldview) #Set appropriate type
		_read(dlg.data, dlg.explore.db, dlg.sel)
	catch e
		#rethrow(e) #for debug purposes
		success = false
	end
	fv_id = getlabelledid(fieldview, dlg.sel)

	#Update ent_fieldinfo:
	_subject = strip(dlg.sel.subject)
	labelfieldinfo = ""
	if isempty_subject(_subject)
		success = false
		labelfieldinfo = MSG_NOSUBJECT
	elseif !success
		labelfieldinfo = string(fv_id, ": ", MSG_CANNOTEDIT)
	else
		labelfieldinfo = string(
			"Sujet: ", _subject,
			", Niveau: ", getid_grade(dlg.sel.grade_idx),
		)
	end

	set_gtk_property!(dlg.ent_fieldinfo, "text", labelfieldinfo)
	set_gtk_property!(dlg.frame_fieldview, "label", fv_id)
	set_gtk_property!(dlg.frame_shortdescr, "visible", false)
	set_gtk_property!(dlg.frame_shortdescr, "no-show-all", !dlg.data.has_shortdescr)
	tb = dlg.tv_shortdescr.buffer[_Gtk.TextBuffer]
		set_gtk_property!(tb, "text", dlg.data.shortdescr)
	tb = dlg.tv_descr.buffer[_Gtk.TextBuffer]
		set_gtk_property!(tb, "text", dlg.data.descr)
		#tb.text[String] = "some text" #Could use this syntax - not sure if supported, though.
	set_gtk_property!(dlg.btn_ok, "sensitive", success) #Conditionaly disable
	set_gtk_property!(dlg.btn_apply, "sensitive", success) #Conditionaly disable

	Gtk.showall(dlg.wnd)
	refresh!(dlg.explore) #Just in case out of sync (might have created data)
	nothing
end

#==_show() wrapper (Create dialog if not registered with .editdlg)
===============================================================================#
_show_editdlg(explore::ExploreWnd, editdlg::ExploreEditDlg) = _show(editdlg)
function _show_editdlg(explore::ExploreWnd, editdlg) #Have not specified edit dialog yet
	explore.editdlg = ExploreEditDlg(explore)
	_show(explore.editdlg)
end
show_editdlg(explore::ExploreWnd) = _show_editdlg(explore, explore.editdlg)

#Last line


text:
export SparseMatrixCOO

"""
    SparseMatrixCOO(is::Vector, js::Vector, vs::Vector, m::Int, n::Int) -> SparseMatrixCOO
    SparseMatrixCOO{Tv, Ti}(is::Vector{Ti}, js::Vector{Ti}, vs::Vector{Tv}, m::Int, n::Int) -> SparseMatrixCOO

A sparse matrix in COOrdinate format.

Also known as the ‘ijv’ or ‘triplet’ format.

# Notes

COO matrices should not be used in arithmetic operations like addition, subtraction, multiplication, division, and matrix power.

### Advantages of the COO format
* facilitates fast conversion among sparse formats
* permits duplicate entries (see example)
* very fast conversion to and from CSR/CSC formats (CSR is not implemented)

### Disadvantages of the COO format
does not directly support:
* arithmetic operations
* slicing

### Intended Usage
* COO is a fast format for constructing sparse matrices
* Once a matrix has been constructed, convert to CSR or CSC format for fast arithmetic and matrix vector operations
* By default when converting to CSR or CSC format, duplicate (i,j) entries will be summed together. This facilitates efficient construction of finite element matrices and the like. (see example)
"""
mutable struct SparseMatrixCOO{Tv,Ti} <: AbstractSparseMatrix{Tv,Ti}
    is::Vector{Ti}
    js::Vector{Ti}
    vs::Vector{Tv}
    m::Int
    n::Int

    function SparseMatrixCOO{Tv,Ti}(
        is::Vector{Ti},
        js::Vector{Ti},
        vs::Vector{Tv},
        m::Int,
        n::Int,
    ) where {Ti,Tv}
        length(is) == length(js) == length(vs) ||
        throw(ArgumentError("Input row, col, data should be equal size."))
        new{Tv,Ti}(is, js, vs, m, n)
    end
end

SparseMatrixCOO(is::Vector{Ti}, js::Vector{Ti}, vs::Vector{Tv}, m, n) where {Ti,Tv} =
    SparseMatrixCOO{Tv,Ti}(is, js, vs, m, n)

copy(coo::SparseMatrixCOO{Tv,Ti}) where {Tv,Ti} =
    SparseMatrixCOO{Tv,Ti}(copy(coo.is), copy(coo.js), copy(coo.vs), coo.m, coo.n)
function copyto!(A::SparseMatrixCOO{Tv,Ti}, B::SparseMatrixCOO{Tv,Ti}) where {Tv,Ti}
    size(A) == size(B) && nnz(A) == nnz(B) ||
    throw(MethodError("size/nnz of two coo matrices do not match!"))
    copyto!(A.is, B.is)
    copyto!(A.js, B.js)
    copyto!(A.vs, B.vs)
    A
end

"""
    allocated_coo(::Type, M::Int, N::Int, nnz::Int) -> SparseMatrixCOO

Construct a preallocated `SparseMatrixCOO` instance.
"""
function allocated_coo(::Type{T}, M::Int, N::Int, nnz::Int) where {T}
    is = Vector{Int}(undef, nnz)
    js = Vector{Int}(undef, nnz)
    vs = Vector{T}(undef, nnz)
    SparseMatrixCOO(is, js, vs, M, N)
end

function getindex(coo::SparseMatrixCOO{Tv,Ti}, i::Ti, j::Ti) where {Tv,Ti}
    res = zero(Tv)
    for k = 1:nnz(coo)
        if coo.is[k] == i && coo.js[k] == j
            res += coo.vs[k]
        end
    end
    res
end

size(coo::SparseMatrixCOO) = (coo.m, coo.n)
size(coo::SparseMatrixCOO, axis::Int) =
    axis == 1 ? coo.m : (axis == 2 ? coo.n : throw(MethodError("invalid axis parameter")))

# SparseArrays: SparseMatrixCSC, nnz, nonzeros, dropzeros!, findnz
nnz(coo::SparseMatrixCOO) = coo.is |> length
nonzeros(coo::SparseMatrixCOO) = coo.vs

function dropzeros!(coo::SparseMatrixCOO{Tv,Ti}; trim::Bool = false) where {Tv,Ti}
    mask = abs.(coo.vs) .> 1e-15
    SparseMatrixCOO{Tv,Ti}(coo.is[mask], coo.js[mask], coo.vs[mask], coo.m, coo.n)
end

findnz(coo::SparseMatrixCOO) = (coo.is, coo.js, coo.vs)
isdense(::SparseMatrixCOO) = false


text:
#####
##### simple problems for testing
#####

@testset "linear problem" begin
    for _ in 1:100
        n = rand(2:10)
        J = randn(n, n)
        x0 = randn(n)
        b = J * x0
        for l in (Dogleg(), )
            @show l
            result = trust_region_solver(x -> (residual = J * x .- b, Jacobian = J),
                                         x0 .* 1000, local_method = l)
            @test result.x ≈ x0 atol = √eps() * n
            @test norm(result.fx.residual, 2) ≈ 0 atol = √eps()
            @test norm(result.fx.residual, 2) == result.residual_norm
            @test result.fx.Jacobian == J
            @test result.converged
            display(result)
        end
    end
end

@testset "infeasible region" begin
    @testset "bounded away from solution" begin
        # the solution x = 0 is infeasible, but do we get close?
        function f(x)
            (residual = x[1] ≥ 1 ? x : x .+ NaN, Jacobian = Diagonal(ones(length(x))))
        end
        result = trust_region_solver(f, [3.0])
        @test !result.converged
        @test result.x ≈ [1.0]
    end

    @testset "jump over narrow infeasible region" begin
        history = Bool[]
        function f3(x)
            is_feasible = abs(x[1] - 2) ≥ 0.1 # first step takes us here
            push!(history, is_feasible)
            (residual = is_feasible ? x .^ 3  : x .+ NaN, Jacobian = Diagonal(@. 3 * abs2(x)))
        end
        result = trust_region_solver(f3, [3.0]; Δ = 5)
        @test result.converged
        @test result.x ≈ [0.0] atol = 0.03
        @test any(!, history)   # check that infeasible region was visited
    end
end

#####
##### nonlinear test problems
#####

TEST_FUNCTIONS = [F_NWp281(),
                  Rosenbrock(),
                  PowellSingular(),
                  PowellBadlyScaled(),
                  HelicalValley()]

@testset "solver tests" begin
    for f in TEST_FUNCTIONS
        for local_method in (Dogleg(), GeneralizedEigenSolver())
            @info "solver test" f local_method
            if false            # condition on broken tests here
                @warn "skipping because broken"
                continue
            end
            res = trust_region_solver(ForwardDiff_wrapper(f, domain_dimension(f)),
                                      starting_point(f); local_method = local_method)
            @test res.converged
            @test res.x ≈ root(f) atol = 1e-4 * domain_dimension(f)
        end
    end
end


text:
## replace base functions with libdevice alternatives

jaxfunc(f) = f
jaxfunc(::Type{T}) where {T} = (x...) -> T(x...) # broadcasting type ctors isn't GPU compatible

const libdevice_1 = Dict(
#complex
:abs=>:abs,
:angle=>:angle,
:conj=>:conj,
:real=>:real, :imag=>:imag,
#trig
:cos=>:cos, :cosh=>:cosh,
:acos=>:arccos, :acosh=>:arccosh,
:sin=>:sin, :sinh=>:sinh,
:asin=>:arcsin, :asinh=>:arcsinh,
:tan=>:tan, :tanh=>:tanh,
:atan=>:arctan, :atanh=>:arctanh,
#
:ceil=>:ceil, :floor=>:floor,
#math
:exp=>:exp, :expm1=>:expm1, :exp2=>:exp2,
:log=>:log, :log1p=>:log1p, :log10=>:log10, :log2=>:log2,
:sqrt=>:sqrt, :inv=>:reciprocal,
:sign=>:sign)

const libdevice_2 = Dict(
#comparison
:(==)=>:equal, :!= => :not_equal,
:> => :greater, :>= => :greater_equal,
:< => :less, :<= => :less_equal,
)

#const libdevice_3 = Dict(
#
#)

const libdevice = merge(libdevice_1,
                        libdevice_2,
                        #libdevice_3
                        )

for (f_jl, f_np) in libdevice
    isdefined(Base, f_jl) || continue
    @eval jaxfunc(::typeof(Base.$f_jl)) = numpy.$f_np
    #@eval Base.$f(x::AbstractJaxArray) = $f(x)
end

jaxliteral_pow(::typeof(^), x::AbstractJaxArray, ::Val{0}) = one(x)
jaxliteral_pow(::typeof(^), x::AbstractJaxArray, ::Val{1}) = x
jaxliteral_pow(::typeof(^), x::AbstractJaxArray, ::Val{2}) = x * x
jaxliteral_pow(::typeof(^), x::AbstractJaxArray, ::Val{3}) = x * x * x
jaxliteral_pow(::typeof(^), x::AbstractJaxArray, ::Val{p}) where p = numpy.power(x, Int32(p))

# ??
jaxliteral_pow(::Base.RefValue{typeof(^)}, x::AbstractJaxArray, ::Base.RefValue{Val{p}}) where p = numpy.power(x, Int32(p))

jaxfunc(::typeof(Base.literal_pow)) = jaxliteral_pow
jaxfunc(::typeof(Base.:(^))) = power

# test
using MacroTools

const _jaxfuncs = [collect(keys(libdevice)); :^; copy(overridenbfuncs)]
jaxfuncs() = (global _jaxfuncs; _jaxfuncs)

_jaxint(x::Int) = Int32(x)#JaxArray{Int32,0}(Jax.int32(x), tuple())
_jaxint(x::Expr) = x.head == :call && x.args[1] == :Int32 && x.args[2] isa Int ?
    Int32(x.args[2]) : x
_jaxint(x) = x

function _jaxpowliteral(x::Expr)
    if x.head == :call && x.args[1] == :(Jax.jaxfunc(^)) && x.args[3] isa Int32
        num = x.args[3]
        if 0 <= num <= 3
            sym = gensym(:x)
            new_x = Expr(:block, :($sym = $(x.args[2])))

            if iszero(num)
                #push!(new_x.args, :(one($sym)))
                push!(new_x.args, :(one(eltype($sym))))
            elseif isone(num)
                push!(new_x.args, sym)
            else
                unroll = Expr(:call, :(Jax.jaxfunc(*)))
                for x = one(num):num
                    push!(unroll.args, sym)
                end
                push!(new_x.args, unroll)
            end

            x = new_x
        end
    end
    x
end
_jaxpowliteral(x) = x

function replace_device(ex)
    global _jaxfuncs
    MacroTools.postwalk(ex) do x
        x = _jaxpowliteral(x)
        x = x in _jaxfuncs ? :(Jax.jaxfunc($x)) : x
        x = _jaxint(x)
        x
    end
end

macro jaxfunc(ex)
    global _jaxfuncs
    def = MacroTools.splitdef(ex)
    f = def[:name]
    def[:name] = Symbol(:jax, f)
    def[:body] = replace_device(def[:body])
    push!(_jaxfuncs, f)
    quote
        $(esc(MacroTools.combinedef(def)))
        Jax.jaxfunc(::typeof($(esc(f)))) = $(esc(def[:name]))
    end
end

for (op, fun) in libdevice_2
    @eval begin
        Base.$op(a::AbstractJaxArray, b::AbstractJaxArray) = numpy.all(numpy.$fun(a.o, b.o))
    end
end


text:
##=================####==============##

# OD Agents Interactions

##=================####==============##

"""
This function looks at the infected neighbors of each agent to compute wether or not
she will get infected. If both agents adopt a cooperative behavior then the probability
of infection gets reduced by both terms. As the SIS and SIR models are basically the same
excep for reinfections then only a flag for SIR is needed.
"""
function SI_coop!(agent, agents;inf_prob=0.1, rec_prob=0.3, coop_red=0.7, R=false)
    infected_coop = [ag.state for ag in agents[agent.coopf] if ag.state == "I"] |> length
    infected_noncoop = [ag.state for ag in agents[agent.non_coopf] if ag.state == "I"] |> length
    if agent.state == "S"
        if agent.at_home
            inf_prob_co = inf_prob * (1-coop_red)^2
            inf_prob_no = inf_prob * (1-coop_red)
        else
            inf_prob_co = inf_prob * (1-coop_red)
            inf_prob_no = inf_prob
        end
        oddsc = [sample([true,false],Weights([inf_prob_co,1-inf_prob_co])) for i in 1:infected_coop] |> sum
        oddsn = [sample([true,false],Weights([inf_prob_no,1-inf_prob_no])) for i in 1:infected_noncoop] |> sum
        odds = oddsc + oddsn
        if odds > 0
            agent.new_state = "I"
        else
            agent.new_state = "S"
        end
    elseif agent.state == "I"
        if !R
            if rand() <= rec_prob
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        else
            if rand() <= rec_prob
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

##=================####==============##
function SI_attitude_old!(agent, agents;inf_prob=0.1, rec_prob=0.3, R=false)
    #  sus = findall(x -> x.state == "S",agents)
    if agent.state == "S"
        infec = [ag.coop_effect for ag in agents[agent.contacts_t] if ag.state == "I"]
        the_odds = @. inf_prob * (1-agent.coop_effect) * (1-infec)
        odds = the_odds |> f -> map(x -> sample([true,false], Weights([x,1-x])),f) |> sum
        dis = truncated(Normal(7,5),1,Inf)
        if odds > 0
            agent.new_state = "I"
            agent.days = Int(ceil(rand(dis)))
        else
            agent.new_state = "S"
        end
    elseif agent.state == "I"
        if !R
            if rand() <= rec_prob
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        else
            if rand() <= rec_prob
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

##=================####==============##

function SI_attitude!(agent, agents;inf_prob=0.3, rec_prob=0.03, R=false)
    #  sus = findall(x -> x.state == "S",agents)
    if agent.state == "S"
        infec = [ag.coop_effect for ag in agents[agent.contacts_t] if ag.state == "I"]
        the_odds = @. inf_prob * (1-agent.coop_effect) * (1-infec)
        odds = the_odds |> f -> map(x -> sample([true,false], Weights([x,1-x])),f) |> sum
        if odds > 0
            agent.new_state = "I"
        else
            agent.new_state = "S"
        end
    elseif agent.state == "I"
        if !R
            if agent.counter >= agent.recovery_t
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        else
            if agent.counter >= agent.recovery_t
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

##=================####==============##

function next_state!(agents;fun=SI_coop!,kwargs...)
    for agent in agents
        fun(agent,agents;kwargs...)
    end
end

##=================####==============##

function SI_next!(agents; inf_prob=0.1, rec_prob=0.3)
    for agent in agents
        if agent.state == "S"
            infected = [ag.state for ag in agents[agent.contacts_t] if ag.state == "I"] |> length
            odds = [sample([true,false],Weights([inf_prob,1-inf_prob])) for i in 1:infected] |> sum
            if odds > 0
                agent.new_state = "I"
            else
                agent.new_state = "S"
            end
        elseif agent.state == "I"
            if rand() <= rec_prob
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        end
    end
end

##=================####==============##

function SIR_next!(agents; inf_prob=0.1, rec_prob=0.3)
    for agent in agents
        if agent.state == "S"
            infected = [ag.state for ag in agents[agent.contacts_t] if ag.state == "I"] |> length
            odds = [sample([true,false],Weights([inf_prob,1-inf_prob])) for i in 1:infected] |> sum
            if odds > 0
                agent.new_state = "I"
            else
                agent.new_state = "S"
            end
        elseif agent.state == "I"
            if rand() <= rec_prob
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

##=================####==============##


text:
# syntax: proto2
# COV_EXCL_START
using ProtoBuf
import ProtoBuf.meta

mutable struct TokenProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TokenProto(; kwargs...)
        obj = new(meta(TokenProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TokenProto
const __meta_TokenProto = Ref{ProtoMeta}()
function meta(::Type{TokenProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TokenProto)
            __meta_TokenProto[] = target = ProtoMeta(TokenProto)
            req = Symbol[:identifier,:password,:kind,:service]
            allflds = Pair{Symbol,Union{Type,String}}[:identifier => Array{UInt8,1}, :password => Array{UInt8,1}, :kind => AbstractString, :service => AbstractString]
            meta(target, TokenProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TokenProto[]
    end
end
function Base.getproperty(obj::TokenProto, name::Symbol)
    if name === :identifier
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :password
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :kind
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :service
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct GetDelegationTokenRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GetDelegationTokenRequestProto(; kwargs...)
        obj = new(meta(GetDelegationTokenRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct GetDelegationTokenRequestProto
const __meta_GetDelegationTokenRequestProto = Ref{ProtoMeta}()
function meta(::Type{GetDelegationTokenRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GetDelegationTokenRequestProto)
            __meta_GetDelegationTokenRequestProto[] = target = ProtoMeta(GetDelegationTokenRequestProto)
            req = Symbol[:renewer]
            allflds = Pair{Symbol,Union{Type,String}}[:renewer => AbstractString]
            meta(target, GetDelegationTokenRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GetDelegationTokenRequestProto[]
    end
end
function Base.getproperty(obj::GetDelegationTokenRequestProto, name::Symbol)
    if name === :renewer
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct GetDelegationTokenResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GetDelegationTokenResponseProto(; kwargs...)
        obj = new(meta(GetDelegationTokenResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct GetDelegationTokenResponseProto
const __meta_GetDelegationTokenResponseProto = Ref{ProtoMeta}()
function meta(::Type{GetDelegationTokenResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GetDelegationTokenResponseProto)
            __meta_GetDelegationTokenResponseProto[] = target = ProtoMeta(GetDelegationTokenResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[:token => TokenProto]
            meta(target, GetDelegationTokenResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GetDelegationTokenResponseProto[]
    end
end
function Base.getproperty(obj::GetDelegationTokenResponseProto, name::Symbol)
    if name === :token
        return (obj.__protobuf_jl_internal_values[name])::TokenProto
    else
        getfield(obj, name)
    end
end

mutable struct RenewDelegationTokenRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RenewDelegationTokenRequestProto(; kwargs...)
        obj = new(meta(RenewDelegationTokenRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RenewDelegationTokenRequestProto
const __meta_RenewDelegationTokenRequestProto = Ref{ProtoMeta}()
function meta(::Type{RenewDelegationTokenRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RenewDelegationTokenRequestProto)
            __meta_RenewDelegationTokenRequestProto[] = target = ProtoMeta(RenewDelegationTokenRequestProto)
            req = Symbol[:token]
            allflds = Pair{Symbol,Union{Type,String}}[:token => TokenProto]
            meta(target, RenewDelegationTokenRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RenewDelegationTokenRequestProto[]
    end
end
function Base.getproperty(obj::RenewDelegationTokenRequestProto, name::Symbol)
    if name === :token
        return (obj.__protobuf_jl_internal_values[name])::TokenProto
    else
        getfield(obj, name)
    end
end

mutable struct RenewDelegationTokenResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RenewDelegationTokenResponseProto(; kwargs...)
        obj = new(meta(RenewDelegationTokenResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RenewDelegationTokenResponseProto
const __meta_RenewDelegationTokenResponseProto = Ref{ProtoMeta}()
function meta(::Type{RenewDelegationTokenResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RenewDelegationTokenResponseProto)
            __meta_RenewDelegationTokenResponseProto[] = target = ProtoMeta(RenewDelegationTokenResponseProto)
            req = Symbol[:newExpiryTime]
            allflds = Pair{Symbol,Union{Type,String}}[:newExpiryTime => UInt64]
            meta(target, RenewDelegationTokenResponseProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RenewDelegationTokenResponseProto[]
    end
end
function Base.getproperty(obj::RenewDelegationTokenResponseProto, name::Symbol)
    if name === :newExpiryTime
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    else
        getfield(obj, name)
    end
end

mutable struct CancelDelegationTokenRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CancelDelegationTokenRequestProto(; kwargs...)
        obj = new(meta(CancelDelegationTokenRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CancelDelegationTokenRequestProto
const __meta_CancelDelegationTokenRequestProto = Ref{ProtoMeta}()
function meta(::Type{CancelDelegationTokenRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CancelDelegationTokenRequestProto)
            __meta_CancelDelegationTokenRequestProto[] = target = ProtoMeta(CancelDelegationTokenRequestProto)
            req = Symbol[:token]
            allflds = Pair{Symbol,Union{Type,String}}[:token => TokenProto]
            meta(target, CancelDelegationTokenRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CancelDelegationTokenRequestProto[]
    end
end
function Base.getproperty(obj::CancelDelegationTokenRequestProto, name::Symbol)
    if name === :token
        return (obj.__protobuf_jl_internal_values[name])::TokenProto
    else
        getfield(obj, name)
    end
end

mutable struct CancelDelegationTokenResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CancelDelegationTokenResponseProto(; kwargs...)
        obj = new(meta(CancelDelegationTokenResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CancelDelegationTokenResponseProto
const __meta_CancelDelegationTokenResponseProto = Ref{ProtoMeta}()
function meta(::Type{CancelDelegationTokenResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CancelDelegationTokenResponseProto)
            __meta_CancelDelegationTokenResponseProto[] = target = ProtoMeta(CancelDelegationTokenResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, CancelDelegationTokenResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CancelDelegationTokenResponseProto[]
    end
end

export TokenProto, GetDelegationTokenRequestProto, GetDelegationTokenResponseProto, RenewDelegationTokenRequestProto, RenewDelegationTokenResponseProto, CancelDelegationTokenRequestProto, CancelDelegationTokenResponseProto
# COV_EXCL_STOP