text:
"""
    NormSpectralBridge{T}

The `NormSpectralCone` is representable with a PSD constraint, since
``t \\ge \\sigma_1(X)`` if and only if ``[tI X^\\top; X tI] \\succ 0``.
"""
struct NormSpectralBridge{T, F, G} <: AbstractBridge
    row_dim::Int # row dimension of X
    column_dim::Int # column dimension of X
    psd_index::CI{F, MOI.PositiveSemidefiniteConeTriangle}
end
function bridge_constraint(::Type{NormSpectralBridge{T, F, G}}, model::MOI.ModelLike, f::G, s::MOI.NormSpectralCone) where {T, F, G}
    f_scalars = MOIU.eachscalar(f)
    t = f_scalars[1]
    row_dim = s.row_dim
    column_dim = s.column_dim
    side_dim = row_dim + column_dim
    psd_set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
    psd_func = MOIU.zero_with_output_dimension(F, MOI.dimension(psd_set))
    for i in 1:side_dim
        MOIU.operate_output_index!(+, T, trimap(i, i), psd_func, t)
    end
    X_idx = 2
    for j in 1:column_dim, i in (column_dim + 1):side_dim
        MOIU.operate_output_index!(+, T, trimap(i, j), psd_func, f_scalars[X_idx])
        X_idx += 1
    end
    psd_index = MOI.add_constraint(model, psd_func, psd_set)
    return NormSpectralBridge{T, F, G}(row_dim, column_dim, psd_index)
end

MOI.supports_constraint(::Type{NormSpectralBridge{T}}, ::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormSpectralCone}) where T = true
MOIB.added_constrained_variable_types(::Type{<:NormSpectralBridge}) = Tuple{DataType}[]
MOIB.added_constraint_types(::Type{NormSpectralBridge{T, F, G}}) where {T, F, G} = [(F, MOI.PositiveSemidefiniteConeTriangle)]
function concrete_bridge_type(::Type{<:NormSpectralBridge{T}}, G::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormSpectralCone}) where T
    F = MOIU.promote_operation(vcat, T, MOIU.scalar_type(G), T)
    return NormSpectralBridge{T, F, G}
end

# Attributes, Bridge acting as a model
MOI.get(bridge::NormSpectralBridge{T, F, G}, ::MOI.NumberOfConstraints{F, MOI.PositiveSemidefiniteConeTriangle}) where {T, F, G} = 1
MOI.get(bridge::NormSpectralBridge{T, F, G}, ::MOI.ListOfConstraintIndices{F, MOI.PositiveSemidefiniteConeTriangle}) where {T, F, G} = [bridge.psd_index]

# References
MOI.delete(model::MOI.ModelLike, bridge::NormSpectralBridge) = MOI.delete(model, bridge.psd_index)

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintFunction, bridge::NormSpectralBridge{T, F, G}) where {T, F, G}
    psd_func = MOIU.eachscalar(MOI.get(model, MOI.ConstraintFunction(), bridge.psd_index))
    t = psd_func[1]
    side_dim = bridge.row_dim + bridge.column_dim
    X = psd_func[[trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim + 1):side_dim]]
    return MOIU.convert_approx(G, MOIU.operate(vcat, T, t, X))
end
MOI.get(model::MOI.ModelLike, ::MOI.ConstraintSet, bridge::NormSpectralBridge) = MOI.NormSpectralCone(bridge.row_dim, bridge.column_dim)
MOI.supports(::MOI.ModelLike, ::MOI.ConstraintPrimalStart, ::Type{<:NormSpectralBridge}) = true
function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintPrimalStart}, bridge::NormSpectralBridge)
    primal = MOI.get(model, attr, bridge.psd_index)
    t = primal[1]
    side_dim = bridge.row_dim + bridge.column_dim
    X = primal[[trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim + 1):side_dim]]
    return vcat(t, X)
end
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintPrimalStart, bridge::NormSpectralBridge{T}, value) where T
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    primal = zeros(T, div(side_dim * (side_dim + 1), 2))
    for i in 1:side_dim
        primal[trimap(i, i)] = value[1]
    end
    X_idx = 2
    for j in 1:column_dim, i in (column_dim + 1):side_dim
        primal[trimap(i, j)] = value[X_idx]
        X_idx += 1
    end
    MOI.set(model, MOI.ConstraintPrimalStart(), bridge.psd_index, primal)
    return
end
# Given [U X'; X V] is dual on PSD constraint, the dual on NormSpectralCone
# constraint is (tr(U) + tr(V), 2X) in NormNuclearCone.
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintDual, bridge::NormSpectralBridge)
    dual = MOI.get(model, MOI.ConstraintDual(), bridge.psd_index)
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    t = sum(dual[trimap(i, i)] for i in 1:side_dim)
    X = 2 * dual[[trimap(i, j) for j in 1:column_dim for i in (column_dim + 1):side_dim]]
    return vcat(t, X)
end

"""
    NormNuclearBridge{T}

The `NormNuclearCone` is representable with an SDP constraint and extra variables,
since ``t \\ge \\sum_i \\sigma_i (X)`` if and only if there exists symmetric
matrices ``U, V`` such that ``[U X^\\top; X V] \\succ 0`` and ``t \\ge (tr(U) + tr(V)) / 2``.
"""
struct NormNuclearBridge{T, F, G, H} <: AbstractBridge
    row_dim::Int # row dimension of X
    column_dim::Int # column dimension of X
    U::Vector{MOI.VariableIndex}
    V::Vector{MOI.VariableIndex}
    ge_index::CI{F, MOI.GreaterThan{T}}
    psd_index::CI{G, MOI.PositiveSemidefiniteConeTriangle}
end
function bridge_constraint(::Type{NormNuclearBridge{T, F, G, H}}, model::MOI.ModelLike, f::H, s::MOI.NormNuclearCone) where {T, F, G, H}
    f_scalars = MOIU.eachscalar(f)
    row_dim = s.row_dim
    column_dim = s.column_dim
    side_dim = row_dim + column_dim
    U_dim = div(column_dim * (column_dim + 1), 2)
    V_dim = div(row_dim * (row_dim + 1), 2)
    U = MOI.add_variables(model, U_dim)
    V = MOI.add_variables(model, V_dim)
    diag_vars = vcat([U[trimap(i, i)] for i in 1:column_dim], [V[trimap(i, i)] for i in 1:row_dim])
    ge_index = MOIU.normalize_and_add_constraint(model, MOIU.operate(-, T, f_scalars[1], MOIU.operate!(/, T, MOIU.operate(sum, T, diag_vars), T(2))), MOI.GreaterThan(zero(T)), allow_modify_function=true)
    psd_set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
    psd_func = MOI.VectorOfVariables(U)
    nuc_dim = 1 + row_dim * column_dim
    for i in 1:row_dim
        row_i = (1 + i):row_dim:nuc_dim
        psd_func = MOIU.operate(vcat, T, psd_func, f_scalars[row_i])
        psd_func = MOIU.operate(vcat, T, psd_func, MOI.VectorOfVariables(V[[trimap(i, j) for j in 1:i]]))
    end
    psd_index = MOI.add_constraint(model, psd_func, psd_set)
    return NormNuclearBridge{T, F, G, H}(row_dim, column_dim, U, V, ge_index, psd_index)
end

MOI.supports_constraint(::Type{NormNuclearBridge{T}}, ::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormNuclearCone}) where T = true
MOIB.added_constrained_variable_types(::Type{<:NormNuclearBridge}) = Tuple{DataType}[]
MOIB.added_constraint_types(::Type{NormNuclearBridge{T, F, G, H}}) where {T, F, G, H} = [(F, MOI.GreaterThan{T}), (G, MOI.PositiveSemidefiniteConeTriangle)]
function concrete_bridge_type(::Type{<:NormNuclearBridge{T}}, H::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormNuclearCone}) where T
    S = MOIU.scalar_type(H)
    F = MOIU.promote_operation(-, T, S, MOIU.promote_operation(/, T, MOIU.promote_operation(+, T, T, T), T))
    G = MOIU.promote_operation(vcat, T, MOI.VectorOfVariables, H)
    return NormNuclearBridge{T, F, G, H}
end

# Attributes, Bridge acting as a model
MOI.get(bridge::NormNuclearBridge, ::MOI.NumberOfVariables) = length(bridge.U) + length(bridge.V)
MOI.get(bridge::NormNuclearBridge, ::MOI.ListOfVariableIndices) = vcat(bridge.U, bridge.V)
MOI.get(bridge::NormNuclearBridge{T, F, G, H}, ::MOI.NumberOfConstraints{F, MOI.GreaterThan{T}}) where {T, F, G, H} = 1
MOI.get(bridge::NormNuclearBridge{T, F, G, H}, ::MOI.NumberOfConstraints{G, MOI.PositiveSemidefiniteConeTriangle}) where {T, F, G, H} = 1
MOI.get(bridge::NormNuclearBridge{T, F, G, H}, ::MOI.ListOfConstraintIndices{F, MOI.GreaterThan{T}}) where {T, F, G, H} = [bridge.ge_index]
MOI.get(bridge::NormNuclearBridge{T, F, G, H}, ::MOI.ListOfConstraintIndices{G, MOI.PositiveSemidefiniteConeTriangle}) where {T, F, G, H} = [bridge.psd_index]

# References
function MOI.delete(model::MOI.ModelLike, bridge::NormNuclearBridge)
    MOI.delete(model, bridge.ge_index)
    MOI.delete(model, bridge.psd_index)
    MOI.delete(model, bridge.U)
    MOI.delete(model, bridge.V)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintFunction, bridge::NormNuclearBridge{T, F, G, H}) where {T, F, G, H}
    ge_func = MOI.get(model, MOI.ConstraintFunction(), bridge.ge_index)
    psd_func = MOIU.eachscalar(MOI.get(model, MOI.ConstraintFunction(), bridge.psd_index))
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    t = MOIU.operate(+, T, ge_func, MOIU.operate(/, T, MOIU.operate(+, T, [psd_func[trimap(i, i)] for i in 1:side_dim]...), T(2)))
    t = MOIU.remove_variable(MOIU.remove_variable(t, bridge.U), bridge.V)
    X = psd_func[[trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim + 1):side_dim]]
    return MOIU.convert_approx(H, MOIU.operate(vcat, T, t, X))
end
MOI.get(model::MOI.ModelLike, ::MOI.ConstraintSet, bridge::NormNuclearBridge) = MOI.NormNuclearCone(bridge.row_dim, bridge.column_dim)
MOI.supports(::MOI.ModelLike, ::MOI.ConstraintDualStart, ::Type{<:NormNuclearBridge}) = true
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintPrimal, bridge::NormNuclearBridge)
    ge_primal = MOI.get(model, MOI.ConstraintPrimal(), bridge.ge_index)
    psd_primal = MOI.get(model, MOI.ConstraintPrimal(), bridge.psd_index)
    side_dim = bridge.row_dim + bridge.column_dim
    t = ge_primal + sum(psd_primal[trimap(i, i)] for i in 1:side_dim) / 2
    X = psd_primal[[trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim + 1):side_dim]]
    return vcat(t, X)
end
# Given t is dual on GreaterThan constraint and [U X'; X V] is dual on PSD constraint,
# the dual on NormNuclearCone constraint is (t, 2X) in NormNuclearCone.
function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintDual, MOI.ConstraintDualStart}, bridge::NormNuclearBridge)
    t = MOI.get(model, attr, bridge.ge_index)
    psd_dual = MOI.get(model, attr, bridge.psd_index)
    side_dim = bridge.row_dim + bridge.column_dim
    X = 2 * psd_dual[[trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim + 1):side_dim]]
    return vcat(t, X)
end
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintDualStart, bridge::NormNuclearBridge{T}, value) where T
    MOI.set(model, MOI.ConstraintDualStart(), bridge.ge_index, value[1])
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    dual = zeros(T, div(side_dim * (side_dim + 1), 2))
    for i in 1:side_dim
        dual[trimap(i, i)] = value[1]
    end
    X_idx = 2
    for j in 1:column_dim, i in (column_dim + 1):side_dim
        dual[trimap(i, j)] = value[X_idx] / 2
        X_idx += 1
    end
    MOI.set(model, MOI.ConstraintDualStart(), bridge.psd_index, dual)
    return
end


text:
#
# features.jl --
#
# Definitions of features for Andor cameras.
#
#-------------------------------------------------------------------------------
#
# This file is part of "AndorCameras.jl" released under the MIT license.
#
# Copyright (C) 2017-2021, Éric Thiébaut.
#

"""
    using AndorCameras.Features

makes all Andor cameras feature definitions available.

"""
module Features

import
    ..AT,
    ..AbstractFeature,
    ..BooleanFeature,
    ..BooleanOrEnumeratedFeature,
    ..CommandFeature,
    ..EnumeratedFeature,
    ..FloatingPointFeature,
    ..IntegerFeature,
    ..StringFeature,
    ..isavailable,
    ..isimplemented

# Export the following methods for those who `using AndorCameras.Features`
# (feature constants will be exported while being defined below):
export
    isavailable,
    isimplemented

for (sym, T) in ((:AccumulateCount, IntegerFeature),
                 (:AOIBinning, EnumeratedFeature),
                 (:AOIHBin, IntegerFeature),
                 (:AOIHeight, IntegerFeature),
                 (:AOILayout, EnumeratedFeature),
                 (:AOILeft, IntegerFeature),
                 (:AOIStride, IntegerFeature),
                 (:AOITop, IntegerFeature),
                 (:AOIVBin, IntegerFeature),
                 (:AOIWidth, IntegerFeature),
                 (:AcquisitionStart, CommandFeature),
                 (:AcquisitionStop,  CommandFeature),
                 (:AlternatingReadoutDirection, BooleanFeature),
                 (:AuxiliaryOutSource, EnumeratedFeature),
                 (:AuxiliaryOutSourceTwo, EnumeratedFeature),
                 (:BackoffTemperatureOffset, FloatingPointFeature),
                 (:Baseline, IntegerFeature),
                 (:BaselineLevel, IntegerFeature), # FIXME: deprecated?
                 (:BitDepth, EnumeratedFeature),
                 (:BufferOverflowEvent, IntegerFeature),
                 (:BytesPerPixel, FloatingPointFeature),
                 (:CameraAcquiring, BooleanFeature),
                 (:CameraDump, CommandFeature),
                 (:CameraFamily, StringFeature),
                 (:CameraMemory, IntegerFeature),
                 (:CameraModel, StringFeature),
                 (:CameraName, StringFeature),
                 (:CameraPresent, BooleanFeature),
                 (:ColourFilter, EnumeratedFeature),
                 (:ControllerID, StringFeature),
                 (:CoolerPower, FloatingPointFeature),
                 (:CycleMode, EnumeratedFeature),
                 (:DDGIOCEnable, BooleanFeature),
                 (:DDGIOCNumberOfPulses, IntegerFeature),
                 (:DDGIOCPeriod, IntegerFeature),
                 (:DDGInsertionDelay, EnumeratedFeature),
                 (:DDGOpticalWidthEnable, BooleanFeature),
                 (:DDGOutputDelay, IntegerFeature),
                 (:DDGOutputEnable, BooleanFeature),
                 (:DDGOutputPolarity, EnumeratedFeature),
                 (:DDGOutputSelector, EnumeratedFeature),
                 (:DDGOutputStepEnable, BooleanFeature),
                 (:DDGOutputWidth, IntegerFeature),
                 (:DDGQueue, CommandFeature),
                 (:DDGStepCount, IntegerFeature),
                 (:DDGStepDelayCoefficientA, FloatingPointFeature),
                 (:DDGStepDelayCoefficientB, FloatingPointFeature),
                 (:DDGStepDelayMode, EnumeratedFeature),
                 (:DDGStepEnabled, BooleanFeature),
                 (:DDGStepUploadModeValues, CommandFeature),
                 (:DDGStepUploadProgress, IntegerFeature),
                 (:DDGStepUploadRequired, BooleanFeature),
                 (:DDGStepWidthCoefficientA, FloatingPointFeature),
                 (:DDGStepWidthCoefficientB, FloatingPointFeature),
                 (:DDGStepWidthMode, EnumeratedFeature),
                 (:DDR2Type, StringFeature),
                 (:DeviceCount, IntegerFeature),
                 (:DeviceVideoIndex, IntegerFeature),
                 (:DisableShutter, BooleanFeature),
                 (:DriverVersion, StringFeature),
                 (:ElectronicShutteringMode, EnumeratedFeature),
                 (:EventEnable, BooleanFeature),
                 (:EventSelector, EnumeratedFeature),
                 (:EventsMissedEvent, IntegerFeature),
                 (:ExposedPixelHeight, IntegerFeature),
                 (:ExposureEndEvent, IntegerFeature),
                 (:ExposureStartEvent, IntegerFeature),
                 (:ExposureTime, FloatingPointFeature),
                 (:ExternalIOReadout, BooleanFeature),
                 (:ExternalTriggerDelay, FloatingPointFeature),
                 (:FanSpeed, EnumeratedFeature),
                 (:FastAOIFrameRateEnable, BooleanFeature),
                 (:FirmwareVersion, StringFeature),
                 (:ForceShutterOpen, BooleanFeature),
                 (:FrameCount, IntegerFeature),
                 (:FrameInterval, FloatingPointFeature),
                 (:FrameIntervalTiming, BooleanFeature),
                 (:FrameRate, FloatingPointFeature),
                 (:FullAOIControl, BooleanFeature),
                 (:GateMode, EnumeratedFeature),
                 (:HeatSinkTemperature, FloatingPointFeature),
                 (:IOControl, EnumeratedFeature),
                 (:IODirection, BooleanOrEnumeratedFeature),
                 (:IOInvert, BooleanFeature),
                 (:IOSelector, EnumeratedFeature),
                 (:IOState, BooleanFeature),
                 (:IRPreFlashEnable, BooleanFeature),
                 (:ImageSizeBytes, IntegerFeature),
                 (:InputVoltage, FloatingPointFeature),
                 (:InterfaceType, StringFeature),
                 (:KeepCleanEnable, BooleanFeature),
                 (:KeepCleanPostExposureEnable, BooleanFeature),
                 (:LUTIndex, IntegerFeature),
                 (:LUTValue, IntegerFeature),
                 (:LineScanSpeed, FloatingPointFeature),
                 (:MCPGain, IntegerFeature),
                 (:MCPIntelligate, BooleanFeature),
                 (:MCPVoltage, IntegerFeature),
                 (:MaxInterfaceTransferRate, FloatingPointFeature),
                 (:MetadataEnable, BooleanFeature),
                 (:MetadataFrame, BooleanFeature),
                 (:MetadataTimestamp, BooleanFeature),
                 (:MicrocodeVersion, StringFeature),
                 (:MultitrackBinned, BooleanFeature),
                 (:MultitrackCount, IntegerFeature),
                 (:MultitrackEnd, IntegerFeature),
                 (:MultitrackSelector, IntegerFeature),
                 (:MultitrackStart, IntegerFeature),
                 (:Overlap, BooleanFeature),
                 (:PixelCorrection, EnumeratedFeature),
                 (:PixelEncoding, EnumeratedFeature),
                 (:PixelHeight, FloatingPointFeature),
                 (:PixelReadoutRate, EnumeratedFeature),
                 (:PixelWidth, FloatingPointFeature),
                 (:PortSelector, IntegerFeature),
                 (:PreAmpGain, EnumeratedFeature),
                 (:PreAmpGainChannel, EnumeratedFeature),
                 (:PreAmpGainControl, EnumeratedFeature),
                 (:PreAmpGainSelector, EnumeratedFeature),
                 (:PreAmpGainValue, IntegerFeature),
                 (:PreAmpOffsetValue, IntegerFeature),
                 (:PreTriggerEnabled, BooleanFeature),
                 (:ReadoutTime, FloatingPointFeature),
                 (:RollingShutterGlobalClear, BooleanFeature),
                 (:RowNExposureEndEvent, IntegerFeature),
                 (:RowNExposureStartEvent, IntegerFeature),
                 (:RowReadTime, FloatingPointFeature),
                 (:ScanSpeedControlEnable, BooleanFeature),
                 (:SensorCooling, BooleanFeature),
                 (:SensorHeight, IntegerFeature),
                 (:SensorModel, StringFeature),
                 (:SensorReadoutMode, EnumeratedFeature),
                 (:SensorTemperature, FloatingPointFeature),
                 (:SensorType, EnumeratedFeature),
                 (:SensorWidth, IntegerFeature),
                 (:SerialNumber, StringFeature),
                 (:ShutterAmpControl, BooleanFeature),
                 (:ShutterMode, EnumeratedFeature),
                 (:ShutterOutputMode, EnumeratedFeature),
                 (:ShutterState, BooleanFeature),
                 (:ShutterStrobePeriod, FloatingPointFeature),
                 (:ShutterStrobePosition, FloatingPointFeature),
                 (:ShutterTransferTime, FloatingPointFeature),
                 (:SimplePreAmpGainControl, EnumeratedFeature),
                 (:SoftwareTrigger, CommandFeature),
                 (:SoftwareVersion, StringFeature),
                 (:SpuriousNoiseFilter, BooleanFeature),
                 (:StaticBlemishCorrection, BooleanFeature),
                 (:SynchronousTriggering, BooleanFeature),
                 (:TargetSensorTemperature, FloatingPointFeature),
                 (:TemperatureControl, EnumeratedFeature),
                 (:TemperatureStatus, EnumeratedFeature),
                 (:TimestampClock, IntegerFeature),
                 (:TimestampClockFrequency, IntegerFeature),
                 (:TimestampClockReset, CommandFeature),
                 (:TransmitFrames, BooleanFeature),
                 (:TriggerMode, EnumeratedFeature),
                 (:UsbDeviceId, IntegerFeature),
                 (:UsbProductId, IntegerFeature),
                 (:VerticallyCenterAOI, BooleanFeature), # FIXME: typo?
                 (:VerticallyCentreAOI, BooleanFeature))
    @eval begin
        const $sym = $(T(sym))
        export $sym
    end
end

end # module Features

ext: jl


text:
module CrossEntropys

using ExprRules
using ..ProbabilisticExprRules

using ExprOptimization: ExprOptAlgorithm, ExprOptResult, BoundedPriorityQueue, enqueue!
import ExprOptimization: optimize

export CrossEntropy

abstract type InitializationMethod end
abstract type TrackingMethod end

"""
    CrossEntropy

Cross Entropy method.
# Arguments
- `pop_size::Int`: population size
- `iterations::Int`: number of iterations
- `max_depth::Int`: maximum depth of derivation tree
- `top_k::Int`: top k elite samples used in selection
- `p_init::Float64`: initial value when fitting MLE
- `init_method::InitializationMethod`: Initialization method
- `track_method::TrackingMethod`: additional tracking, e.g., track top k exprs (default: no additional tracking)
"""
struct CrossEntropy <: ExprOptAlgorithm
    pop_size::Int
    iterations::Int
    max_depth::Int
    top_k::Int
    p_init::Float64
    init_method::InitializationMethod
    track_method::TrackingMethod

    function CrossEntropy(
        pop_size::Int,                                  #population size
        iterations::Int,                                #number of iterations
        max_depth::Int,                                 #maximum depth of derivation tree
        top_k::Int,                                     #top k elite samples used in selection
        p_init::Float64=0.0,                            #initial value when fitting MLE
        init_method::InitializationMethod=RandomInit(); #initialization method
        track_method::TrackingMethod=NoTracking())   #tracking method

        new(pop_size, iterations, max_depth, top_k, p_init, init_method, track_method)
    end
end

"""
    RandomInit

Uniformly random initialization method.
"""
struct RandomInit <: InitializationMethod end

"""
    NoTracking

No additional tracking of expressions.
"""
struct NoTracking <: TrackingMethod end

"""
    TopKTracking

Track the top k expressions.
"""
struct TopKTracking <: TrackingMethod
    k::Int
    q::BoundedPriorityQueue{RuleNode,Float64}

    function TopKTracking(k::Int)
        q = BoundedPriorityQueue{RuleNode,Float64}(k,Base.Order.Reverse) #lower is better
        obj = new(k, q)
        obj
    end
end

"""
    optimize(p::CrossEntropy, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)

Expression tree optimization using the cross-entropy method with parameters p, grammar 'grammar', and start symbol typ, and loss function 'loss'.  Loss function has the form: los::Float64=loss(node::RuleNode, grammar::Grammar)

See: Rubinstein, "Optimization of Computer Simulation Models with Rare Events", European Journal of Operations Research, 99, 89-112, 1197
"""
function optimize(p::CrossEntropy, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)
    cross_entropy(p, grammar, typ, loss; kwargs...)
end

"""
    cross_entropy(p::CrossEntropy, grammar::Grammar, typ::Symbol)

Expression tree optimization using cross-entropy method with parameters p, grammar 'grammar', and start symbol typ, and loss function 'loss'.  Loss function has the form: los::Float64=loss(node::RuleNode, grammar::Grammar)

See: Rubinstein, "Optimization of Computer Simulation Models with Rare Events", European Journal of Operations Research, 99, 89-112, 1197
"""
function cross_entropy(p::CrossEntropy, grammar::Grammar, typ::Symbol, loss::Function;
    verbose::Bool=false)
    iseval(grammar) && error("Cross-entropy does not support _() functions in the grammar")

    dmap = mindepth_map(grammar)
    losses = Vector{Float64}(undef,p.pop_size)
    pcfg = ProbabilisticGrammar(grammar)
    pop = initialize(p.init_method, p.pop_size, pcfg, typ, dmap, p.max_depth)
    best_tree, best_loss = evaluate!(p, loss, grammar, pop, losses, RuleNode(0), Inf)
    for iter = 1:p.iterations
        verbose && println("iterations: $iter of $(p.iterations)")
        for i in eachindex(pop)
            pop[i] = rand(RuleNode, pcfg, typ, dmap, p.max_depth)
        end
        fit_mle!(pcfg, pop[1:p.top_k], p.p_init)
        best_tree, best_loss = evaluate!(p, loss, grammar, pop, losses, best_tree, best_loss)
    end
    alg_result = Dict{Symbol,Any}()
    _add_result!(alg_result, p.track_method)
    ExprOptResult(best_tree, best_loss, get_executable(best_tree, grammar), alg_result)
end

"""
    _add_result!(d::Dict{Symbol,Any}, t::NoTracking)

Add tracking results to alg_result.  No op for NoTracking.
"""
_add_result!(d::Dict{Symbol,Any}, t::NoTracking) = nothing
"""
    _add_result!(d::Dict{Symbol,Any}, t::TopKTracking)

Add tracking results to alg_result.
"""
function _add_result!(d::Dict{Symbol,Any}, t::TopKTracking)
    d[:top_k] = collect(t.q)
    d
end

"""
    initialize(::RandomInit, pop_size::Int, grammar::Grammar, typ::Symbol, dmap::AbstractVector{Int},
        max_depth::Int)

Random population initialization.
"""
function initialize(::RandomInit, pop_size::Int, pcfg::ProbabilisticGrammar, typ::Symbol,
            dmap::AbstractVector{Int}, max_depth::Int)
    [rand(RuleNode, pcfg, typ, dmap, max_depth) for i = 1:pop_size]
end

"""
    evaluate!(p::CrossEntropy, loss::Function, grammar::Grammar, pop::Vector{RuleNode},
        losses::Vector{Float64}, best_tree::RuleNode, best_loss::Float64)

Evaluate the loss function for population and sort.  Update the globally best tree, if needed.
"""
function evaluate!(p::CrossEntropy, loss::Function, grammar::Grammar, pop::Vector{RuleNode},
    losses::Vector{Float64}, best_tree::RuleNode, best_loss::Float64)

    Threads.@threads for i in 1:length(losses)
        losses[i] = loss(pop[i], grammar)
    end

    perm = sortperm(losses)
    pop[:], losses[:] = pop[perm], losses[perm]
    if losses[1] < best_loss
        best_tree, best_loss = pop[1], losses[1]
    end
    _update_tracker!(p.track_method, pop, losses)
    (best_tree, best_loss)
end

"""
    _update_tracker!(t::NoTracking, pop::Vector{RuleNode}, losses::Vector{Union{Float64,Missing}})

Update the tracker.  No op for NoTracking.
"""
function _update_tracker!(t::NoTracking, pop::Vector{RuleNode}, losses::Vector{Float64})
    nothing
end
"""
    _update_tracker!(t::TopKTracking, pop::Vector{RuleNode}, losses::Vector{Union{Float64,Missing}})

Update the tracker.  Track top k expressions.
"""
function _update_tracker!(t::TopKTracking, pop::Vector{RuleNode}, losses::Vector{Float64})
    n = 0
    for i = 1:length(pop)
        r = enqueue!(t.q, pop[i], losses[i])
        r >= 0 && (n += 1) #no clash, increment counter
        n >= t.k && break
    end
end

end #module


text:
using VortexDistributions, Test

# @testset "Multivortex creation and detection" begin
@test found_near(30, 10)
# end


text:
## ---------------------------------------------------------------------
# Base intake info
# The Base model will have a medium (expressed as open intake fluxes) that resamble the cultivation
# at xi = 1 using the set up in Rath 2017 exp A. So the lb of the intakes will be directly the (negative)
# concentration in 42_MAX_UB standard medium (see Cossio's paper). Also, a few intakes, not justified in
# the standard medium will be add based in the intakes of the original model FBA analysis.
function load_base_intake_info(; inf_medium = true)

    base_intake_info = Dict()
    mets_map = load_mets_map()
    met_readable_ids = load_met_readable_ids()
    exch_met_map = load_exch_met_map()

    # From 42_MAX_UB standard medium
    for rath_met in RathData.ALL_METS

        # base_model id
        model_met = mets_map[rath_met]
        exch_rxn = exch_met_map[model_met]

        # 42_MAX_UB standard medium
        # we take the openest version of the intakes for building the
        # base model
        conc = maximum(RathData.cval(rath_met, RathData.STSTS, 0.0))
        lb = -ABS_MAX_BOUND # intake bound
        base_intake_info[exch_rxn] = Dict("c" => conc, "lb" => lb)
    end

    # From ham's medium see Hams_medium.jl
    for (Ham_id, conc) in load_Ham_medium(;inf_medium = inf_medium)
        model_met = met_readable_ids[Ham_id]
        exch_rxn = exch_met_map[model_met]
        haskey(base_intake_info, exch_rxn) && continue # Not overwrite 42_MAX_UB standard medium

        lb = -ABS_MAX_BOUND # intake bound
        base_intake_info[exch_rxn] = Dict("c" => conc, "lb" => lb)
    end

    return base_intake_info

end

"""
    returns a copy of the base_intake_info but with the
    feed medium concentration of a given Rath steady state
"""
function stst_base_intake_info(stst)
    intake_info = load_base_intake_info()
    mets_map = load_mets_map()
    exch_met_map = load_exch_met_map()

    # The feed medium of each steady state only vary
    # in the composition of this mets (see Rath2017)
    for rath_met in ["GLC", "GLN", "GAL"]
        model_met = mets_map[rath_met]
        model_exch = exch_met_map[model_met]
        intake_info[model_exch]["c"] = RathData.cval(rath_met, stst)
    end
    return intake_info
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#=
Emitters are defined by Pixels and SpatialLayouts. An emitter has a spectrum, and an optical power distribution over the hemisphere.
These are intrinsic physical properties of the emitter.
=#

module Emitters
export generate, visual_size, apply
export right, up, forward
export Spectrum, Directions, Origins, AngularPower, Sources
export pointemitter, collimatedemitter

using ...OpticSim, ..Geometry
using LinearAlgebra
import Unitful.Length
using Unitful.DefaultSymbols

# define some utility functions
right(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,1], t[2,1], t[3,1]))
up(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,2], t[2,2], t[3,2]))
forward(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,3], t[2,3], t[3,3]))
OpticSim.origin(t::OpticSim.Geometry.Transform{<:Real}) = Vec3(t[1,4], t[2,4], t[3,4])

# defining name placeholders to override in nested modules

"""
    generate(???)

[TODO]
"""
function generate end

"""
    visual_size(???)

[TODO]
"""
function visual_size end

"""
    apply(???)

[TODO] Returns ray power
"""
function apply end

include("Spectrum.jl")
include("Directions.jl")
include("Origins.jl")
include("AngularPower.jl")
include("Sources.jl")
include("StandardEmitters.jl")

end # module Emitters

export Emitters


text:
struct TVec{
    K,V,D<:Storage{K,V},A<:AbstractVector{D},S<:StochasticStyle,
    NR,ID
} <: AbstractDVec{K,V}
    segments::A
    style::S
end

function TVec{K,V}(
    ; style=default_style(V), num_segments=4,
    _num_ranks=nothing,
    _rank_id=nothing,
) where {K,V}
    comm = MPI.COMM_WORLD

    # Setting the following only supported for debugging.
    num_ranks = isnothing(_num_ranks) ? MPI.Comm_size(comm) : _num_ranks
    rank_id = isnothing(_rank_id) ? MPI.Comm_rank(comm) : _rank_id
    segments = [Storage{K,V}() for _ in 1:num_segments * Threads.nthreads()]

    TVec{K,V,eltype(segments),typeof(segments),typeof(style),num_ranks,rank_id}(
        segments, style
    )
end
function TVec(ps::Vararg{Pair{K,V}}; kwargs...) where {K,V}
    t = TVec{K,V}(; kwargs...)
    for (k, v) in ps
        t[k] = v
    end
    return t
end
function TVec(pairs; kwargs...)
    K = typeof(first(pairs)[1])
    V = typeof(first(pairs)[2])
    t = TVec{K,V}(; kwargs...)
    for (k, v) in pairs
        t[k] = v
    end
    return t
end

# Properties
rank_id(t::TVec{<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,ID}) where {ID} = ID
num_ranks(t::TVec{<:Any,<:Any,<:Any,<:Any,<:Any,NR}) where {NR} = NR
is_distributed(t::TVec) = num_ranks(t) > 1
num_segments(t::TVec) = length(t.segments)
Rimu.StochasticStyle(s::TVec) = s.style

function Base.length(s::TVec)
    if is_distributed(s)
        MPI.Allreduce(sum(length, s.segments), +, MPI.COMM_WORLD)
    else
        sum(length, s.segments)
    end
end

function Base.copyto!(dst::TVec, src::TVec)
    if num_segments(dst) == num_segments(src)
        Folds.foreach(enumerate(src.segments)) do (i, s_s)
            copy!(dst.segments[i], s_s)
        end
        return dst
    else
        # Mismatched number of segments - copy by moving each element over.
        return invoke(copyto!, Tuple{AbstractDVec, TVec}, dst, src)
    end
end
function Base.copy!(dst::TVec, src::TVec)
    return copyto!(dst, src)
end
function Base.copy(dst::TVec, src::TVec)
    return copy!(empty(dst), src)
end

"""
     target_segment(t::TVec, hash::UInt64)

Determine the target segment from key hash. For MPI distributed vectors, this may return
numbers that are out of range.
"""
function target_segment(t, hash::UInt64)
    total_segments = num_segments(t) * num_ranks(t)
    return mod1(hash, total_segments) % Int - rank_id(t) * num_segments(t)
end

function Base.getindex(t::TVec{K}, key::K) where {K}
    h = hash(key)
    target = target_segment(t, h)
    if 1 ≤ target ≤ num_segments(t)
        getindex(t.segments[target_segment(t, h)], key, h)
    else
        error("attempted to `getindex` from different MPI rank")
    end
end

function Base.setindex!(t::TVec{K,V}, val, key::K) where {K,V}
    h = hash(key)
    target = target_segment(t, h)
    if 1 ≤ target ≤ num_segments(t)
        return setindex!(t.segments[target_segment(t, h)], V(val), key, h)
    else
        return V(val)
    end
end
function Rimu.deposit!(t::TVec{K,V}, key::K, val, parent) where {K,V}
    h = hash(key)
    target = target_segment(t, h)
    if 1 ≤ seg ≤ num_segments(t)
        deposit!(t.segments[target_segment(t, h)], key, V(val), h, parent)
    end
    return nothing
end

###
### empty(!)
###
function Base.empty(s::TVec{K,V}; style=s.style) where {K,V}
    return TVec{K,V}(; style, num_segments=num_segments(s) ÷ Threads.nthreads())
end
function Base.empty(s::TVec{K}, ::Type{V}; style=s.style) where {K,V}
    return TVec{K,V}(; style, num_segments=num_segments(s) ÷ Threads.nthreads())
end
function Base.empty(s::TVec, ::Type{K}, ::Type{V}; style=s.style) where {K,V}
    return TVec{K,V}(; style, num_segments=num_segments(s) ÷ Threads.nthreads())
end

function Base.empty!(s::TVec)
    Folds.foreach(empty!, s.segments)
    return s
end
function Base.sizehint!(s::TVec, n)
    n_per_segment = cld(n, length(s.segments))
    Folds.foreach(d -> sizehint!(d, n_per_segment), s.segments)
    return s
end

###
### Iterators and parallel reduction
###
struct TVecIterator{F,V,T,D}
    fun::F
    segments::V

    function TVecIterator(fun::F, ::Type{T}, tv::TVec) where {F,T}
        return new{F,typeof(tv.segments),T,is_distributed(tv)}(fun, tv.segments)
    end
end
Base.values(s::TVec) = TVecIterator(values, valtype(s), s)
Base.keys(s::TVec) = TVecIterator(keys, keytype(s), s)
Base.pairs(s::TVec) = TVecIterator(pairs, eltype(s), s)

Base.eltype(s::Type{<:TVecIterator{<:Any,<:Any,T}}) where {T} = T
Base.length(s::TVecIterator) = sum(length, s.segments)
is_distributed(s::TVecIterator{<:Any,<:Any,<:Any,D}) where {D} = D

function Base.iterate(s::TVecIterator, i::Int=1)
    if is_distributed(s)
        @warn "Vector is distributed, iterating over local entries only." maxlog=1
    end
    if i > length(s.segments)
        return nothing
    end
    it = iterate(s.fun(s.segments[i]))
    if isnothing(it)
        return iterate(s, i + 1)
    else
        return it[1], (i, it[2])
    end
end
function Base.iterate(s::TVecIterator, (i, st))
    it = iterate(s.fun(s.segments[i]), st)
    if isnothing(it)
        return iterate(s, i + 1)
    else
        return it[1], (i, it[2])
    end
end

function Base.mapreduce(f, op, s::TVecIterator; kwargs...)
    result = Folds.mapreduce(op, s.segments; kwargs...) do v
        mapreduce(f, op, s.fun(v); kwargs...)
    end
    # The compiler will know whether it's distributed or not and elide this statement.
    if is_distributed(s)
        return MPI.Allreduce(result, op, MPI.COMM_WORLD)
    else
        return result
    end
end
function Base.map!(f, s::TVecIterator{typeof(values)})
    Folds.foreach(s.segments) do v
        map!(f, values(v))
    end
    return s
end
function Base.map!(f, dst::AbstractDVec, src::TVecIterator{typeof(values)})
    @assert length(dst.segments) == length(src.segments)
    Folds.foreach(zip(dst.segments, src.segments)) do (dst_s, src_s)
        map!(f, dst_s, values(src_s))
    end
    return dst
end

function Base.:*(α::Number, dv::TVec)
    T = promote_type(typeof(α), valtype(dv))
    if T == valtype(dv)
        if iszero(α)
            result = similar(dv)
        else
            result = copy(dv)
            map!(x -> α * x, result, values(dv))
        end
    else
        result = similar(dv, T)
        if !iszero(α)
            map!(x -> α * x, result, values(dv))
        end
    end
    return result
end

function LinearAlgebra.rmul!(dv::TVec, α::Number)
    map!(x -> x * α, values(dv))
    return dv
end
function LinearAlgebra.mul!(dst::TVec, src::TVec, α::Number)
    return map!(x -> α * x, dst, values(src))
end

function LinearAlgebra.axpby!(α, v::TVec, β::Number, w::TVec)
    rmul!(w, β)
    axpy!(α, v, w)
end
function LinearAlgebra.axpy!(α, v::TVec, w::TVec)
    Folds.foreach(zip(w.segments, v.segments)) do (w_s, v_s)
        add!(w_s, v_s, α)
    end
    return w
end
function LinearAlgebra.dot(v::TVec, w::TVec)
    # TODO Check for match or fall back to default implementation
    T = promote_type(valtype(v), valtype(w))
    result = Folds.sum(zip(v.segments, w.segments)) do (v_s, w_s)
        dot(v_s, w_s)
    end::T
    if is_distributed(v)
        return MPI.Allreduce(result, +, MPI.COMM_WORLD)
    else
        return result
    end
end
function LinearAlgebra.dot(v::AbstractDVec, w::TVec)
    result = invoke(dot, Tuple{AbstractDVec,AbstractDVec}, v, w)
    if is_distributed(w)
        return MPI.Allreduce(result, +, MPI.COMM_WORLD)
    else
        return result
    end
end
LinearAlgebra.dot(w::TVec, v::AbstractDVec) = dot(v, w)

function Base.real(v::TVec)
    dst = similar(v, real(valtype(v)))
    map!(real, dst, values(v))
end
function Base.imag(v::TVec)
    dst = similar(v, real(valtype(v)))
    map!(imag, dst, values(v))
end

function Base.:+(v::AbstractDVec, w::AbstractDVec)
    result = similar(v, promote_type(valtype(v), valtype(w)))
    copy!(result, v)
    add!(result, w)
    return result
end
function Base.:-(v::AbstractDVec, w::AbstractDVec)
    result = similar(v, promote_type(valtype(v), valtype(w)))
    copy!(result, v)
    return axpy!(-one(valtype(result)), w, result)
end
function LinearAlgebra.normalize!(v::AbstractDVec, p=2)
    n = norm(v, p)
    rmul!(v, inv(n))
    return v
end
function LinearAlgebra.normalize(v::AbstractDVec, p=2)
    res = copy(v)
    return normalize!(res, p)
end


text:
"Functions for parsing Crystals / SymOps from text / .cif files"

"Strips trailing uncertainty values from a String, then parses as a Float"
function _parse_cif_float(str::String) :: Float64
    i = findfirst('(', str)
    str = isnothing(i) ? str : str[1:i-1]
    return parse(Float64, str)
end

function parse_number_or_fraction(s)
    # Parse a number or fraction
    cnt = count(==('/'), s)
    if cnt == 0
        return parse(Float64, s)
    elseif cnt == 1
        n, d = parse.(Float64, split(s, '/'))
        return n/d
    else
        error("Cannot parse '$s' as a number.")
    end
end

function format_number_or_fraction(x)
    # check for fraction with denominator 2, 3, 4, 6, or 12
    c = 12
    xc′ = round(Int, x*c)
    if abs(x*c - xc′) < 1e-12
        frac = xc′ // c
        if frac.den == 1
            return @sprintf "%+d" frac.num
        else
            return @sprintf "%+d/%d" frac.num frac.den
        end
    else
        return @sprintf "%+.4g" x
    end
end

function _parse_op(str::AbstractString) :: SymOp
    D = 3
    R = zeros(D, D)
    T = zeros(D)

    @assert length(split(str, ",")) == D

    # The substring s describes the i'th row of the output
    for (i, s) in enumerate(split(str, ","))
        # Remove all spaces
        s = replace(s, " " => "")
        # In weird circumstances, double negatives may be generated. Get rid of them.
        s = replace(s, "--" => "+")
        # Trick to split on either + or - symbols
        s = replace(s, "-" => "+-")

        # Each term t describes the numerical value for one element of R or T
        for t in split(s, '+', keepempty=false)
            j = findfirst(last(t), "xyz")
            if isnothing(j)
                T[i] += parse_number_or_fraction(t)
            else
                coeff = t[begin:end-1]
                R[i, j] += begin
                    if coeff == ""
                        1
                    elseif coeff == "-"
                        -1
                    else
                        parse_number_or_fraction(coeff)
                    end
                end
            end
        end
    end

    # Wrap translations
    T = mod.(T, 1)

    return SymOp(R, T)
end

"""
    Crystal(filename::AbstractString; symprec=1e-5)

Reads the crystal from a `.cif` file located at the path `filename`.
"""
function Crystal(filename::AbstractString; symprec=nothing)
    cif = Cif(Path(filename))
    # For now, assumes there is only one data collection per .cif
    cif = cif[first(keys(cif))]

    a = _parse_cif_float(cif["_cell_length_a"][1])
    b = _parse_cif_float(cif["_cell_length_b"][1])
    c = _parse_cif_float(cif["_cell_length_c"][1])
    α = _parse_cif_float(cif["_cell_angle_alpha"][1])
    β = _parse_cif_float(cif["_cell_angle_beta"][1])
    γ = _parse_cif_float(cif["_cell_angle_gamma"][1])
    lat_vecs = lattice_vectors(a, b, c, α, β, γ)

    geo_table = get_loop(cif, "_atom_site_fract_x")
    xs = _parse_cif_float.(geo_table[:, "_atom_site_fract_x"])
    ys = _parse_cif_float.(geo_table[:, "_atom_site_fract_y"])
    zs = _parse_cif_float.(geo_table[:, "_atom_site_fract_z"])
    unique_atoms = Vec3.(zip(xs, ys, zs))

    sitetypes = String.(geo_table[:, :_atom_site_type_symbol])

    multiplicities = nothing
    if "_atom_site_symmetry_multiplicity" in names(geo_table)
        multiplicities = parse.(Int, geo_table[:, "_atom_site_symmetry_multiplicity"])
    end

    wyckoffs = nothing
    if "_atom_site_wyckoff_symbol" in names(geo_table)
        wyckoffs = String.(geo_table[:, "_atom_site_wyckoff_symbol"])
    end

    # Try to infer symprec from coordinate strings
    # TODO: Use uncertainty information if available from .cif
    if isnothing(symprec)
        strs = vcat(geo_table[:, "_atom_site_fract_x"], geo_table[:, "_atom_site_fract_y"], geo_table[:, "_atom_site_fract_z"])
        elems = vcat(xs, ys, zs)
        # guess fractional errors by assuming each elem is a fraction with simple denominator (2, 3, or 4)
        errs = map(elems) do x
            c = 12
            return abs(rem(x*c, 1, RoundNearest)) / c
        end
        (err, i) = findmax(errs)
        if err < 1e-12
            println("Precision parameter is unspecified, but all coordinates seem to be simple fractions.")
            println("Setting symprec=1e-12.")
            symprec = 1e-12
        elseif 1e-12 < err < 1e-4
            @printf "Precision parameter is unspecified, but coordinate string '%s' seems to have error %.1e.\n" strs[i] err
            symprec = 15err
            @printf "Setting symprec=%.1e.\n" symprec
        else
            println("Error: Cannot infer precision. Please provide an explicit `symprec` parameter to load '$filename'")
            return Nothing
        end
    end

    symmetries = nothing
    for sym_header in ("_space_group_symop_operation_xyz", "_symmetry_equiv_pos_as_xyz")
        if sym_header in keys(cif)
            sym_table = get_loop(cif, sym_header)
            symmetries = _parse_op.(sym_table[:, sym_header])
        end
    end

    groupnum = nothing
    for group_header in ("_space_group_it_number", "_symmetry_int_tables_number")
        if group_header in keys(cif)
            groupnum = parse(Int, cif[group_header][1])
        end
    end

    hall_symbol = nothing
    hall_header = "_space_group_name_hall"
    if hall_header in keys(cif)
        hall_symbol = cif[hall_header][1]
    end

    hm_symbol = nothing
    hm_header = "_space_group_name_h-m_alt"
    if hm_header in keys(cif)
        hm_symbol = cif[hm_header][1]
    end

    spacegroup = if !isnothing(hm_symbol)
        if !isnothing(groupnum)
            "HM symbol '$hm_symbol' ($groupnum)"
        else
            "HM symbol '$hm_symbol'"
        end
    else
        if !isnothing(groupnum)
            "Spacegroup number $groupnum"
        else
            ""
        end
    end

    if !isnothing(symmetries)
        # Use explicitly provided symmetries
        return crystal_from_symops(lat_vecs, unique_atoms, sitetypes, symmetries, spacegroup; symprec)
    elseif !isnothing(hall_symbol)
        # Use symmetries for Hall symbol
        return Crystal(lat_vecs, unique_atoms, sitetypes, hall_symbol; symprec)
    elseif !isnothing(groupnum)
        # Use symmetries for international group number
        return Crystal(lat_vecs, unique_atoms, sitetypes, groupnum; symprec)
    else
        # Infer the symmetries automatically
        return Crystal(lat_vecs, unique_atoms, sitetypes; symprec)
    end
end


text:
function ueberpruefeParameter(M,b,nx1,nx2,nx3,u,v,b3)

f = [2.35; 9.14; -0.31];
u1 = u[1];
u2 = u[2];
u3 = u[3];
#dritte Spalte der Matrix
v1 = v[1];
v2 = v[2];
v3 = v[3];

A = zeros(3,3);
#erste Spalte der Matrix
A[1,1] = -nx1;
A[2,1] = -nx2;
A[3,1] = -nx3;
#zweite Spalte der Matrix
A[1,2] = u1;
A[2,2] = u2;
A[3,2] = u3;
#dritte Spalte der Matrix
A[1,3] = v1;
A[2,3] = v2;
A[3,3] = v3;

SuS = norm(M*f);
LuL = norm(A*f);

if SuS-LuL <= 0.1 && norm(b)-norm(b3) <= 0.1
    println("Ihr habt das korrekte Gleichungssystem aufgestellt. Die Lösung eures Gleichungssystems lautet:")
    Lvec = M\b3;
    r = Lvec[1];
    s = Lvec[2];
    t = Lvec[3];
    print(round(r;digits=4),"\n");
    print(round(s;digits=4),"\n");
    print(round(t;digits=4));
elseif norm(b)-norm(b3) <= 0.1
    println("Überprüft erneut eure Matrix.")
    r = [];
    s = [];
    t = [];
elseif SuS-LuL <= 0.1
    println("Überprüft erneut euren Vektor b.")
    r = [];
    s = [];
    t = [];
elseif b[1]-b3[1] < 0.1 && b[2]-b3[2] < 0.1 && b[3]-b3[3] < 0.1
    println("Überprüft erneut euren Vektor b.")
    r = [];
    s = [];
    t = [];
elseif round(A[1,1];digits=2) == round(-M[1,1];digits=2) || round(A[1,1];digits=2) == round(-M[1,2];digits=2) || round(A[1,1];digits=2) == round(-M[3,1];digits=2)
    println("Überprüft noch einmal eure Umformungen. Ihr müsst alle Vektoren, die mit Parametern multipliziert werden auf die eine und den einzelnen auf die andere Seite stellen.")
    r = [];
    s = [];
    t = [];
else
    println("Überprüft erneut eure Eingabe und, ob ihr alle Umformungen korrekt gemacht habt.")
    r = [];
    s = [];
    t = [];
end
return r,s,t;
end


text:
###############################################################################
## General Utilities
###############################################################################

logit(X) = log.(X ./ (1 .- X))
logit(X::AbstractNode) = node(x->logit(x), X)

expit(X) = 1 ./ (1 .+ exp.(-X))
expit(X::AbstractNode) = node(x->expit(x), X)

"""

Let's default to no warnings for now.
"""
MLJBase.check(model::TMLEstimator, args... ; full=false) = true

Base.merge(ndt₁::AbstractNode, ndt₂::AbstractNode) =
    node((ndt₁, ndt₂) -> merge(ndt₁, ndt₂), ndt₁, ndt₂)

"""

Adapts the type of the treatment variable passed to the G learner
"""
adapt(T::NamedTuple{<:Any, NTuple{1, Z}}) where Z = T[1]
adapt(T) = T
adapt(T::AbstractNode) = node(adapt, T)

function log_over_threshold(covariate::AbstractNode, threshold)
    node(cov -> findall(x -> x >= 1/threshold, cov), covariate)
end

###############################################################################
## Interactions Generation
###############################################################################

"""
    interaction_combinations(query::NamedTuple{names,})
Returns a generator over the different combinations of interactions that
can be built from a query.
"""
function interaction_combinations(query::NamedTuple{names,}) where names
    return (NamedTuple{names}(query) for query in Iterators.product(query...))
end

"""
    indicator_fns(query::NamedTuple{names,})

Implements the (-1)^{n-j} formula representing the cross-value of
indicator functions,  where:
    - n is the order of interaction considered
    - j is the number of treatment variables different from the "case" value
"""
function indicator_fns(query::NamedTuple{nms,}) where nms
    case = NamedTuple{nms}([v[1] for v in query])
    interactionorder = length(query)
    return Dict(q => (-1)^(interactionorder - sum(q[key] == case[key] for key in nms))
                for q in interaction_combinations(query))
end

###############################################################################
## Offset
###############################################################################
expected_value(ŷ, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Finite}}) = pdf.(ŷ, levels(first(ŷ))[2])
expected_value(ŷ::AbstractNode, t::Type{<:Probabilistic}, s::Type{<:AbstractArray{<:Finite}}) =
    node(ŷ->expected_value(ŷ, t, s), ŷ)

expected_value(ŷ, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:MLJBase.Continuous}}) = mean.(ŷ)
expected_value(ŷ::AbstractNode, t::Type{<:Probabilistic}, s::Type{<:AbstractArray{<:MLJBase.Continuous}}) =
    node(ŷ->expected_value(ŷ, t, s), ŷ)

expected_value(ŷ, ::Type{<:Deterministic}, ::Type{<:AbstractArray{<:MLJBase.Continuous}}) = ŷ
expected_value(ŷ::AbstractNode, t::Type{<:Deterministic}, s::Type{<:AbstractArray{<:MLJBase.Continuous}}) =
    node(ŷ->expected_value(ŷ, t, s), ŷ)

maybelogit(x, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Finite}}) = logit(x)
maybelogit(x, _, _) = x

function compute_offset(mach::Machine, X)
    ŷ = MLJBase.predict(mach, X)
    expectation = expected_value(ŷ, typeof(mach.model), target_scitype(mach.model))
    return maybelogit(expectation, typeof(mach.model), target_scitype(mach.model))
end

###############################################################################
## Covariate
###############################################################################

function indicator_values(indicators, T)
    covariate = zeros(nrows(T))
    for (i, row) in enumerate(Tables.namedtupleiterator(T))
        if haskey(indicators, row)
            covariate[i] = indicators[row]
        end
    end
    covariate
end
indicator_values(indicators, T::AbstractNode) =
    node(t -> indicator_values(indicators, t), T)

plateau_likelihood(likelihood, threshold) = max.(threshold, likelihood)
plateau_likelihood(likelihood::AbstractNode, threshold) =
    node(l -> plateau_likelihood(l, threshold), likelihood)

elemwise_divide(x, y) = x ./ y
elemwise_divide(x::AbstractNode, y::AbstractNode) =
    node((x, y) -> elemwise_divide(x,y), x, y)

"""
For each data point, computes: (-1)^(interaction-oder - j)
Where j is the number of treatments different from the reference in the query.
"""
function compute_covariate(Gmach::Machine, W, T, indicators; threshold=0.005)
    # Compute the indicator value
    indic_vals = indicator_values(indicators, T)

    # Compute density and truncate
    likelihood = density(Gmach, W, T)

    likelihood = plateau_likelihood(likelihood, threshold)

    return elemwise_divide(indic_vals, likelihood)
end

###############################################################################
## Fluctuation
###############################################################################

fluctuation_input(covariate, offset) = (covariate=covariate, offset=offset)
fluctuation_input(covariate::AbstractNode, offset::AbstractNode) =
    node((c, o) -> fluctuation_input(c, o), covariate, offset)

function counterfactualTreatment(vals, T)
    names = keys(vals)
    n = nrows(T)
    NamedTuple{names}(
            [categorical(repeat([vals[name]], n), levels=levels(Tables.getcolumn(T, name)))
                            for name in names])
end

function compute_fluctuation(Fmach::Machine,
                             Q̅mach::Machine,
                             Gmach::Machine,
                             indicators,
                             W,
                             T,
                             X;
                             threshold=0.005)
    offset = compute_offset(Q̅mach, X)
    covariate = compute_covariate(Gmach, W, T, indicators;
                                    threshold=threshold)
    Xfluct = fluctuation_input(covariate, offset)
    return predict_mean(Fmach, Xfluct)
end

###############################################################################
## Report Generation
###############################################################################

function estimation_report(Fmach::Machine,
    Q̅mach::Machine,
    Gmach::Machine,
    Hmach::Machine,
    W::AbstractNode,
    T::AbstractNode,
    observed_fluct::AbstractNode,
    ys::AbstractNode,
    covariate::AbstractNode,
    indicators,
    threshold,
    query)

    node((w, t, o, y, c) -> estimation_report(Fmach, Q̅mach, Gmach, Hmach, w, t, o, y, c, indicators, threshold, query),
                                W, T, observed_fluct, ys, covariate)
end

"""

For a given query, identified by `indicators`, reports the different quantities of
interest. An important intermediate quantity is obtained by aggregation of
E[Y|T, W] evaluated at the various counterfactual values of the treatment.
For instance, if the order of Interaction is 2 with binary variables, this is computed as:

AggregatedCounterfactual = Fluctuation(t₁=1, t₂=1, W=w) - Fluctuation(t₁=1, t₂=0, W=w)
                - Fluctuation(t₁=0, t₂=1, W=w) + Fluctuation(t₁=0, t₂=0, W=w)
"""
function estimation_report(Fmach::Machine,
                            Q̅mach::Machine,
                            Gmach::Machine,
                            Hmach::Machine,
                            W,
                            T,
                            observed_fluct,
                            ys,
                            covariate,
                            indicators,
                            threshold,
                            query)

    tmle_ct_agg = zeros(nrows(T))
    initial_ct_agg = zeros(nrows(T))
    for (vals, sign) in indicators
        counterfactualT = counterfactualTreatment(vals, T)
        Thot = transform(Hmach, counterfactualT)
        X = merge(Thot, W)

        initial_expectation = expected_value(MLJBase.predict(Q̅mach, X), typeof(Q̅mach.model), target_scitype(Q̅mach.model))
        initial_ct_agg += sign*initial_expectation

        tmle_ct_agg += sign*compute_fluctuation(Fmach,
                    Q̅mach,
                    Gmach,
                    indicators,
                    W,
                    counterfactualT,
                    X;
                    threshold=threshold)
    end

    initial_estimate = mean(initial_ct_agg)
    tmle_estimate = mean(tmle_ct_agg)
    inf_curve = influencecurve(covariate, ys, observed_fluct, tmle_ct_agg, tmle_estimate)

    return QueryReport(query, inf_curve, tmle_estimate, initial_estimate)
end

max_stars_repo_name: olivierlabayle/TMLE

max_issues_repo_name: olivierlabayle/TMLE

max_forks_repo_name: olivierlabayle/TMLE


text:
include("finitedifference.jl")
include("tdist_mle.jl")
using Plots

function flatten_reshape(A::Array{T} where T<:Array{G} where G)
    B = reshape(A, length(A), 1)
    C = zeros(length(B), length(B[1, 1]))
    for i = 1:size(C, 1)
        for j = 1:size(C, 2)
            C[i, j] = B[i][j]
        end
    end
    return C
end

btg1 = btg(trainingData([1 2; 3 4; 5 6.0], reshape([1.0 ; 1.0; 1.0], 3, 1), [2, 3, 4.0]), [1 10.0; 2 12.0], [-0.5 0.5]);

f = x ->  (w = Σθ(btg1, x); reshape(w, length(w), 1));
g = x -> (w = partial_Σθ_theta(btg1, x); flatten_reshape(w));

(r1, r2, plt1, pol1) = checkDerivative(f, g, [1.1, 2.1], nothing);

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


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

JLLWrappers.@generate_main_file_header("JSBSim")
JLLWrappers.@generate_main_file("JSBSim", UUID("1f2de530-32af-5b65-a726-9f0f4dbc15ee"))
end  # module JSBSim_jll


text:
using Pkg; Pkg.activate(".")
using Revise
using ArrayLisp
# %%
ArrayLisp.@dump_eval function add(x, y)
    println("x is $x and y is $y")

    # Functions return the value of their last statement
    x + y
end

sexpr"(function (add x y) (println x y) (+ x y))"

# %%
ArrayLisp.@dump_eval(x = 3)
sexpr"(= x 1)"

# %%
ArrayLisp.@dump_eval(
    begin
        b = [4, 5, 6] # => 3-element Array{Int64,1}: [4, 5, 6]
        b = [4; 5; 6] # => 3-element Array{Int64,1}: [4, 5, 6]
        b[1]    # => 4
        b[end]  # => 6
    end
)
sexpr"""
(= b (vect 4 5 6))
(= b (vcat 4 5 6))
[b 1]
[b end]
"""

# %%
x = 1
ArrayLisp.@dump_eval(@show(x))
sexpr"(@show x)"

# %%
ArrayLisp.@dump_eval(b = Int8[4, 5, 6])
sexpr"(= b [Any 4 5 6])"

# %%
ArrayLisp.@dump_eval(begin
    struct A{T}
        x::T
    end
    A(x) = A{Int}(x)
    x = A{Float64}(3.0)
    y = A(4)
end)
sexpr"({Array Int} undef 2 3)"

# %%
ArrayLisp.@dump_eval(using Dates, Pkg)
sexpr"(using Dates ..Pkg)"

# %%
ArrayLisp.@dump_eval(import Dates, Pkg)
sexpr"(import Dates ..Pkg)"

# %%
ArrayLisp.@dump_eval(x = Dates.Date(2019, 6, 12))
sexpr"(import Dates ..Pkg)"

# %%
ArrayLisp.@dump_eval(begin
    struct A{T}
        x::T
    end
    A(x) = A{Int}(x)
    y = A{A{Int}}(A(1))
    y.x.x
end)
sexpr"(= y ({A {A Int}} (A 2))) (. (. y x) x)"

# %%
ArrayLisp.@dump_eval function f(x; y=3, z)
    print("$(x) = ")
    println(y)
end
sexpr"(function (f x (y 3) (z)) (println x y z))"

# %%
ArrayLisp.@dump_eval(f(1, y=2, z=3))
sexpr"(f 1 (kw y 2) (kw z 3))"

# %%

ArrayLisp.@dump_eval(
    for i in 1:4
        println(i)
    end)

sexpr"(for ((i (: 1 4))) (println i))"

ArrayLisp.@dump_eval(
    for i in 1:4, j = 5:8
        println(i, j)
    end)

sexpr"(for ((i (: 1 4)) (j (: 5 8))) (println i j))"

# %%
using ThinkJulia
ArrayLisp.@dump_eval(@svg begin
    t = Turtle()
    for i in 1:4
        forward(t, 100)
        turn(t, -90)
    end
end)

# %%
sexpr"
(@svg (begin
  (= t (Turtle))
  (for ((i (: 1 4)))
    (forward t 100)
    (turn t -90))))
"

# %%
ArrayLisp.@dump_eval(
begin
    """Some documentation"""
    function fab()
        println("Hello")
    end
end
)

sexpr"""
(function
  """Some documentation"""
  (fab)
  (println "Hello"))
"""

# %%
sexpr"""
(println "hello")
;(println "hello")
"""

# %%
ArrayLisp.@dump_eval(
if 1 > 2
    println(1)
end
)

sexpr"""
(if (> 1 2) (println 1))
(if (< 1 2) (println 2))
"""

# %%
ArrayLisp.@dump_eval(
if 1 > 2
    println(1)
else
    println(2)
end
)

sexpr"""
(if (> 1 2)
    (println 1)
    (println 2))
"""

# %%
ArrayLisp.@dump_eval(
if 1 > 2
    println(1)
elseif 2 > 3
    println(3)
elseif 2 < 3
    println(4)
else
    println(2)
end
)

sexpr"""
(if (> 1 2) (println 1)
 elseif (> 2 3) (println 3)
 elseif (< 2 3) (println 4)
 else (println 2))
"""


text:
using LinearAlgebra;
using BivMatFun;
using SpecialFunctions;
using Printf;
using MAT;
using BenchmarkTools;
using DelimitedFiles;

include("experiments_common.jl");

function run_test(use_grcar = false)

    success = true;

    a = 1/2;
    f = (x,y,i,j) -> (x+y).^(-a-i-j) * gamma(a+i+j) ./ gamma(a) * (-1)^(i+j);

    ntests = 5;

    # n, errD, errT, nblocksDA, nblocksDB, nblocksTA, nblocksTB, digitsD, digitsT, errest, timeD, timeT
    data = zeros(2*ntests, 12);

    for i = 1 : ntests
        n = 32 * i;

        vars = matread("$(data_directory)/experiment1_$(n)_random.mat")
        A = vars["A"]; B = vars["B"]; C = vars["C"]; Y = vars["Y"]; errest1 = vars["errest"]

        @printf("n = %d, A,B,C random\n", n);

        XD, infod = fun2m(f, A, B, C, method = BivMatFun.Diag);
        if ! is_ci_test()
            tdiag = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.Diag);
            tdiag = median(tdiag.times) / 1e9
        else
            tdiag = 0.0
        end

        @printf("  Diag: err = %e, nblocks = %d %d, max digits = %d, time = %f, errest = %e\n", norm(XD - Y) / norm(Y), infod.nblocksA, infod.nblocksB, infod.digits, tdiag, errest1);

        XT, infot = fun2m(f, A, B, C, method = BivMatFun.Taylor);

        if ! is_ci_test()
            tt = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.Taylor);
            tt = median(tt.times) / 1e9;
        else
            tt = 0.0
        end

        @printf("  Taylor: err = %e, nblocks = %d %d, max deg = %d, time = %f, errest = %e\n", norm(XT - Y) / norm(Y), infot.nblocksA, infot.nblocksB, infot.digits, tt, errest1);

        success &= norm(XD - Y) / norm(Y) < errest1 * 100;
        success &= norm(XT - Y) / norm(Y) < errest1 * 100;

        # GRCAR
        vars = matread("$(data_directory)/experiment1_$(n)_grcar.mat")
        A = vars["A"]; B = vars["B"]; C = vars["C"]; Y2 = vars["Y"]; errest2 = vars["errest"]

        @printf("n = %d, A grcar, B,C random\n", n);

        XD2, infod2 = fun2m(f, A, B, C, method = BivMatFun.Diag);

        if ! is_ci_test()
            tdiag2 = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.Diag);
            tdiag2 = median(tdiag2.times) / 1e9;
        else
            tdiag2 = 0.0
        end

        @printf("  Diag: err = %e, nblocks = %d %d, max digits = %d, time = %f, errest = %e\n", Float64(norm(XD2 - Y2) / norm(Y2)), infod2.nblocksA, infod2.nblocksB, infod2.digits, tdiag2, errest2);

        BivMatFun.warnings(false);
        XT2, infot2 = fun2m(f, A, B, C, method = BivMatFun.Taylor);
        if ! is_ci_test()
            tt2 = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.Taylor);
            tt2 = median(tt2.times) / 1e9;
        else
            tt2 = 0.0
        end

        BivMatFun.warnings(true);

        @printf("  Taylor: err = %e, nblocks = %d %d, max deg = %d, time = %f, errest = %e\n", Float64(norm(XT2 - Y2) / norm(Y2)), infot2.nblocksA, infot2.nblocksB, infot2.digits, tt2, errest2);

        # For the GRCAR matrix, which involves a-posteriori re-blocking due to a
        # highly non-normal eigenvector matrix, we expect Taylor to fail, hence
        # we do not check it.
        success &= norm(XD2 - Y2) / norm(Y) < errest2 * 100;

        #
        # % n, errD, errT, nblocksDA, nblocksDB, nblocksTA, nblocksTB, digitsD, digitsT, errest, timeD, timeT
        data[2*i-1, 1] = n;
        data[2*i-1, 2] = norm(XD - Y) / norm(Y);
        data[2*i-1, 3] = norm(XT - Y) / norm(Y);
        data[2*i-1, 4] = infod.nblocksA;
        data[2*i-1, 5] = infod.nblocksB;
        data[2*i-1, 6] = infot.nblocksA;
        data[2*i-1, 7] = infot.nblocksB;
        data[2*i-1, 8] = infod.digits;
        data[2*i-1, 9] = infot.digits;
        data[2*i-1, 10] = errest1;
        data[2*i-1, 11] = tdiag;
        data[2*i-1, 12] = tt;
        #
        data[2*i, 1] = n;
        data[2*i, 2] = norm(XD2 - Y2) / norm(Y2);
        data[2*i, 3] = norm(XT2 - Y2) / norm(Y2);
        data[2*i, 4] = infod2.nblocksA;
        data[2*i, 5] = infod2.nblocksB;
        data[2*i, 6] = infot2.nblocksA;
        data[2*i, 7] = infot2.nblocksB;
        data[2*i, 8] = infod2.digits;
        data[2*i, 9] = infot2.digits;
        data[2*i, 10] = errest2;
        data[2*i, 11] = tdiag2;
        data[2*i, 12] = tt2;

        writedlm("exp1.dat", data, '\t');
    end

    return success;

end

run_test(true);

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
"""
    easter_rata(y::Dates.Year) → Int

Returns Easter date as a *[Rata Die](https://en.wikipedia.org/wiki/Rata_Die)* number.

Based on *[Algo R](http://www.linuxtopia.org/online_books/programming_books/python_programming/python_ch38.html)*.
"""
function easter_rata(y::Dates.Year)

    # Algo R only works after 1582
    if y.value < 1582
        # Are you using this? Send me a postcard!
        error("Year cannot be less than 1582. Provided: $(y.value).")
    end

    # Century
    c = div(y.value , 100) + 1

    # Shifted Epact
    local se::Int = mod(14 + 11*(mod(y.value, 19)) - div(3*c, 4) + div(5+8*c, 25), 30)

    # Adjust Epact
    if (se == 0) || ((se == 1) && ( 10 < mod(y.value, 19) ))
       se += 1
    end

    # Paschal Moon
    p = Dates.Date(y.value, 4, 19).instant.periods.value - se

    # Easter: locate the Sunday after the Paschal Moon
    return p + 7 - mod(p, 7)
end

"""
    easter_date(y::Dates.Year) → Dates.Date

Returns result of `easter_rata` as a `Dates.Date` instance.
"""
@inline easter_date(y::Dates.Year) = Dates.Date(Dates.rata2datetime(easter_rata(y)))

"""
    findweekday(weekday_target::Integer, yy::Integer, mm::Integer, occurrence::Integer, ascending::Bool) → Date

Given a year `yy` and month `mm`, finds a date where a choosen weekday occurs.

`weekday_target` values are declared in module `Base.Dates`:
`Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7`.

If `ascending` is true, searches from the beginning of the month. If false, searches from the end of the month.

If `occurrence` is `2` and `weekday_target` is `Monday`, searches the 2nd Monday of the given month, and so on.
"""
function findweekday(weekday_target::Integer, yy::Integer, mm::Integer, occurrence::Integer, ascending::Bool) :: Dates.Date

    dt = Dates.Date(yy, mm, 1)
    local dt_dayofweek::Integer
    local offset::Integer

    @assert occurrence > 0 "occurrence must be > 0. Provided $(occurrence)."

    if ascending
        dt_dayofweek = Dates.dayofweek(dt)
        offset = mod(weekday_target + 7 - dt_dayofweek, 7)
    else
        dt = Dates.lastdayofmonth(dt)
        dt_dayofweek = Dates.dayofweek(dt)
        offset = mod(dt_dayofweek + 7 - weekday_target, 7)
    end

    if occurrence > 1
        offset += 7 * (occurrence - 1)
    end

    if ascending
        return dt + Dates.Day(offset)
    else
        return dt - Dates.Day(offset)
    end
end

"""

   adjustweekendholidayPost(dt, [adjust_saturdays]) → Date

In the UK and Canada, if a holiday falls on Saturday or Sunday, it's observed on the next business day.
This function will adjust to the next Monday.

`adjust_saturdays` kwarg defaults to `true`.
"""
function adjustweekendholidayPost(dt::Dates.Date; adjust_saturdays::Bool = true) :: Dates.Date

    if adjust_saturdays && (Dates.dayofweek(dt) == Dates.Saturday)
        return dt + Dates.Day(2)
    end

    if Dates.dayofweek(dt) == Dates.Sunday
        return dt + Dates.Day(1)
    end

    return dt
end

"""
    adjustweekendholidayUS(dt) → Date

In the United States, if a holiday falls on Saturday, it's observed on the preceding Friday.
If it falls on Sunday, it's observed on the next Monday.
"""
function adjustweekendholidayUS(dt::Dates.Date) :: Dates.Date

    if Dates.dayofweek(dt) == Dates.Saturday
        return dt - Dates.Day(1)
    end

    if Dates.dayofweek(dt) == Dates.Sunday
        return dt + Dates.Day(1)
    end

    return dt
end


text:
include("module.jl")

using Logging
Logging.global_logger(Logging.ConsoleLogger(Logging.Info))

F5_DETECTED_USELESS_NF = 0
F5_USELESS_NF   = 0
F5_NF        = 0
F5_BASIS_SIZE      = 0

function zerof5()
    global F5_DETECTED_USELESS_NF
    global F5_USELESS_NF
    global F5_NF
    global F5_BASIS_SIZE
    F5_DETECTED_USELESS_NF = 0
    F5_USELESS_NF        = 0
    F5_NF      = 0
    F5_BASIS_SIZE   = 0
end

function katsuran(n; ground=QQ)
    R, x = PolynomialRing(ground, ["x$i" for i in 0:n])

    return [
        (sum(x[abs(l)+1]*x[abs(m-l)+1] for l=-n:n if abs(m-l)<=n) -
        x[m+1] for m=0:n-1)...,
        x[1] + 2sum(x[i+1] for i=1:n) - 1
    ]
end

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

function homogenize(F)
    base = parent(F[1])
    s = copy(symbols(base))
    push!(s, :h)
    R, xsu = PolynomialRing(base_ring(base), s, ordering=ordering(base))
    xs = xsu[1:end-1]
    u = xsu[end]

    homF = Vector{eltype(F)}(undef, 0)
    for f in F
        tdeg = total_degree(leading_monomial(f))
        f = evaluate(f, xs)
        homf = zero(R)
        for t in terms(f)
            homf += t * u^(tdeg - total_degree(t))
        end
        push!(homF, homf)
    end

    homF
end

function dehomogenize(homF)
    R = parent(homF[1])
    xs = collect(gens(R))
    xs[end] = one(R)

    F = Vector{eltype(homF)}(undef, 0)
    for homf in homF
        f = evaluate(homf, xs)
        push!(F, f)
    end

    F
end

#-----------------------------------------------------------------------------
# Regular reductions

# try to reduce f with g in a regular way and return result
function regular_reduction_step(f::ModuleElement, g::ModuleElement)
    evalf = f.ev
    evalg = g.ev

    # @info "reducing $(f) with $(g)" evalf evalg

    for t in terms(evalf)
        # if divides
        success, u = divides(t, leading_term(evalg))
        if success
            u = divexact(t, leading_term(evalg))
            @info "$(leading_term(evalg)) | $t"

            # if reduction is regular (and does not change signature)
            # not really correct!
            if f > u*g
                newf = f - u*g
                # @info "regular! $(f) - $(u)*$(g) --> $newf"
                return true, newf
            end
        end
    end
    return false, f
end

# // same as above, but in terms of set G //
# try to reduce f with G in a regular way and return result
function regular_reduction_step(f::ModuleElement, G)
    evalf = f.ev
    for g in G
        success, newf = regular_reduction_step(f, g)
        if success
            return true, newf
        end
    end
    return false, f
end

# regular normal form of f w.r.t. G
function regular_normal_form(f::ModuleElement, G)
    # @info "computing reg normalform of $f w.r.t $G.."

    success = true
    newf = copy(f)
    while success
        success, newf = regular_reduction_step(newf, G)
        # @info "reduction $success"
    end
    newf
end

#-----------------------------------------------------------------------------
# Singular reductions

# if f is singurarly top reducible by G
function issingularlytopreducible(f::ModuleElement, G)
    leadevalf = leading_monomial(f.ev)
    for g in G
        evalg = g.ev
        success, u = divides(leadevalf, leading_monomial(evalg))
        if success
            if (g*u).sgn == f.sgn
                return true
            end
        end
    end
    return false
end

# if f is singurarly top reducible by G
function syzygy_criterion(f::ModuleElement, syzygies)
    for syz in syzygies
        success, u = divides(f.sgn.monom, syz.sgn.monom)
        if success && f.sgn.index == syz.sgn.index
            # @warn "DISCARDED $f by $syz"
            global F5_DETECTED_USELESS_NF
            F5_DETECTED_USELESS_NF += 1
            return true
        end
    end
    return false
end

#-----------------------------------------------------------------------------
# S-polynomial

# returns u, v,  such that
# u*a - v*b is spoly(a, b)
function mults(a::ModuleElement, b::ModuleElement)
    u = lcm(leading_monomial(a.ev), leading_monomial(b.ev))
    u = divexact(u, leading_term(a.ev))

    v = lcm(leading_monomial(a.ev), leading_monomial(b.ev))
    v = divexact(v, leading_term(b.ev))

    u, v
end

# S-polynomial of a and b
function spoly(a::ModuleElement, b::ModuleElement)
    u, v = mults(a, b)
    u*a - v*b
end

#-----------------------------------------------------------------------------
# Groebner basis things

function construct_module(F)
    [
        ModuleElement(f, Signature(i, one(f)))
        for (i, f) in enumerate(F)
    ]
end

function signature_groebner_basis(F)
    sort!(F, by=leading_monomial, lt=f5_total_degree_lead_cmp)

    F = map(f -> map_coefficients(c -> c // leading_coefficient(f), f), F)

    G = construct_module(F)

    syzygies = []
    for i in 1:length(G)
        for j in i+1:length(G)
            s1 = (G[i]*leading_monomial(G[j].ev)).sgn
            s2 = (G[j]*leading_monomial(G[i].ev)).sgn
            push!(syzygies, ModuleElement( zero(F[1]), max(s1, s2) ) )
        end
    end

    P = []
    for i in 1:length(G)
        for j in i+1:length(G)
            if syzygy_criterion(spoly(G[i], G[j]), syzygies)
                continue
            end

            push!(P, spoly(G[i], G[j]))
        end
    end

    # @info "generated initial G and P:"
    # @info "F = $F"
    # @info "G = $G"
    # @info "P = $P"
    # @info "syz = $syzygies"

    while !isempty(P)
        sort!(P, rev=false)
        f = popfirst!(P)
        @info "selected $f"

        if syzygy_criterion(f, syzygies)
            continue
        end

        # global F5_NF
        # F5_NF += 1
        fNF = regular_normal_form(f, G)
        global F5_NF
        F5_NF += 1

        if issyzygy(fNF)
            # @warn "Reduction to zero!"
            # global F5_REDUCED
            # F5_REDUCED += 1
            push!(syzygies, fNF)

            global F5_USELESS_NF
            F5_USELESS_NF += 1
        elseif !issingularlytopreducible(fNF, G)
            # update P
            for fj in G

                if syzygy_criterion(spoly(fNF, fj), syzygies)
                    continue
                end

                push!(P, spoly(fNF, fj))

                # @info "$fNF  $fj"
                # @info "SPOLY $(last(P))"
            end

            # update G
            fNFnormalzed = fNF * inv(leading_coefficient(fNF.ev))
            push!(G, fNFnormalzed)
        else
            global F5_USELESS_NF
            F5_USELESS_NF += 1
        end

        # @info "updated G and P"
        # @info "G = $G"
        # @info "P = $P"
        # @info "syz = $syzygies"
    end

    # global F5_SIZE
    # F5_SIZE = length(G) - length(F)

    # global F5_REDUCED
    # global F5_DISCARDED
    global F5_BASIS_SIZE
    F5_BASIS_SIZE = length(G)

    G
end

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

using Logging
Logging.global_logger(ConsoleLogger(Logging.Warn))

R, (x,y,z) = PolynomialRing(QQ, ["x","y", "z"], ordering=:degrevlex)

F = [x*z + 1, y*z + 1]

# evaluation
G = signature_groebner_basis(F)
Gev = [g.ev for g in G]
println("############")
println(G)
println(Gev)


text:
## SINGULARITY DISTRIBUTIONS
#============================================#

## Helper functions

function polar_transform(x, z, x1, x2)
    s1, s2 = x - x1, x - x2
    r1, r2 = √(s1^2 + z^2), √(s2^2 + z^2)
    θ1, θ2 = atan(z, s1), atan(z, s2)

    r1, r2, θ1, θ2
end

# "Fundamental" terms
term_1(xp, x, z, d, r1, r2, θ1, θ2) = d - z * (θ2 - θ1) + (x - xp) * log(r2 / r1)
term_2(xp, x, z, r1, r2, θ1, θ2) = (x - xp) * (θ2 - θ1) + z * log(r2 / r1)

boundary_term(xp, x, d) = (1/2π, 1/2 * (x - xp) / d)

# VORTEX SINGULARITIES
#============================================#

## Linear-strength singularities

# Streamfunctions
function linear_vortex_stream_minus(γ, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1

    # Singularity checks
    ϵ = 1e-10
    if r1 < ϵ; logr1 = 0.; else logr1 = log(r1) end
    if r2 < ϵ; logr2 = 0.; else logr2 = log(r2) end

    ψ_m = γ / 2π * (z * (θ2 - θ1) - d + x * logr1 - (x - d) * logr2)
end

function linear_vortex_stream_plus(γ, x, z, x1, x2)
    r1, r2, _, _= polar_transform(x, z, x1, x2)
    d   = x2 - x1
    ψ_m = linear_vortex_stream_minus(1., x, z, x1, x2)

    ϵ = 1e-10
    if r1 < ϵ; logr1 = 0.; else logr1 = log(r1) end
    if r2 < ϵ; logr2 = 0.; else logr2 = log(r2) end

    ψ_p = γ / d * (x * ψ_m + 1/4π * (r2^2 * (logr2 - 1/2) - r1^2 * (logr1 - 1/2)))
end

# Velocities
function linear_vortex_velocity_a(γ1, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1

    # ua = term_2(x2, x, z, r1, r2, θ1, θ2)
    # wa = term_1(x2, x, z, d, r1, r2, θ1, θ2)

    dlogr = log(r1 / r2)
    dθ    = (θ2 - θ1)

    u =     dθ / 2π + (z * dlogr     - x * dθ) / (2π * d)
    w = -dlogr / 2π + (x * dlogr - d + z * dθ) / (2π * d)

    γ1 * SVector(u, w)
end

function linear_vortex_velocity_b(γ2, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1

    # ub = term_2(x1, x, z, r1, r2, θ1, θ2)
    # wb = term_1(x1, x, z, d, r1, r2, θ1, θ2)

    dlogr = log(r1 / r2)
    dθ    = (θ2 - θ1)

    u = -(z * dlogr     - x * dθ) / (2π * d)
    w = -(x * dlogr - d + z * dθ) / (2π * d)

    γ2 * SVector(u, w)
end

# SOURCE SINGULARITIES
#============================================#

## Constant-strength singularities

# Streamfunctions
function constant_source_stream(σ, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1

    # Singularity checks
    ϵ = 1e-10;
    if r1 < ϵ; logr1 = 0.; θ1 = π; θ2 = π; else logr1 = log(r1); end
    if r2 < ϵ; logr2 = 0.; θ1 = 0.; θ2 = 0.; else logr2 = log(r2); end

    ψ = σ / 2π * (x * (θ1 - θ2) + d * θ2 + z * (logr1 - logr2))

    # Branch cuts
    dψ = x2
    if (θ1 + θ2) > π
        ψ -= 0.25 * dψ
    else
        ψ += 0.75 * dψ
    end
end

# Velocity
function constant_source_velocity(σ, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    @. σ / 2π * SVector(-log(r2 / r1), θ2 - θ1)
end

## Linear-strength singularities

# Streamfunctions
function linear_source_stream_minus(σ, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1

    # Branch cuts
    if θ1 < 0.; θ1 += 2π; end
    if θ2 < 0.; θ2 += 2π; end

    # Singularity checks
    ϵ = 1e-10;
    if r1 < ϵ; logr1 = 0.; θ1 = π; θ2 = π; else logr1 = log(r1); end
    if r2 < ϵ; logr2 = 0.; θ1 = 0.; θ2 = 0.; else logr2 = log(r2); end

    ψ = σ / 2π * (x * (θ1 - θ2) + d * θ2 + z * (logr1 - logr2))
end

function linear_source_stream_plus(σ, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1
    ψ_m = linear_source_stream_minus(1., x, z, x1, x2)

    # Branch cuts
    if θ1 < 0.; θ1 += 2π; end
    if θ2 < 0.; θ2 += 2π; end

    # Singularity checks
    ϵ = 1e-10;
    if r1 < ϵ; logr1 = 0.; θ1 = π; θ2 = π; else logr1 = log(r1); end
    if r2 < ϵ; logr2 = 0.; θ1 = 0.; θ2 = 0.; else logr2 = log(r2); end

    ψ_p = σ / d * (x * ψ_m + 1/4π * (r2^2 * θ2 - r1^2 * θ1 - z * d))
end

# Velocities
function linear_source_velocity_a(σ1, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1
    # if abs(z) <= 1e-12
        # ua, wa = # boundary_term(x2, x, d)
        # σ1 .* SVector(ua, -wa)
    # else
        # ua = term_1(x2, x, z, d, r1, r2, θ1, θ2)
        # wa = term_2(x2, x, z, r1, r2, θ1, θ2)

        # @. σ1 / (2π * (x2 - x1)) * SVector(ua, -wa)
    # end

    dlogr = log(r1 / r2)
    dθ    = (θ2 - θ1)

    u = dlogr / 2π - (x * dlogr - d + z * dθ) / (2π * d)
    w =    dθ / 2π + (z * dlogr     - x * dθ) / (2π * d)

    σ1 * SVector(u, w)
end

function linear_source_velocity_b(σ2, x, z, x1, x2)
    r1, r2, θ1, θ2 = polar_transform(x, z, x1, x2)
    d = x2 - x1
    # if abs(z) <= 1e-12
    #     ub, wb = # boundary_term(x1, x, d)
    #     σ2 .* SVector(-ub, wb)
    # else
    #     ub = term_1(x1, x, z, d, r1, r2, θ1, θ2)
    #     wb = term_2(x1, x, z, r1, r2, θ1, θ2)

    #     @. σ2 / (2π * (x2 - x1)) * SVector(-ub, wb)
    # end

    dlogr = log(r1 / r2)
    dθ    = (θ2 - θ1)

    u =  (x * dlogr - d + z * dθ) / (2π * d)
    w = -(z * dlogr     - x * dθ) / (2π * d)

    σ2 * SVector(u, w)
end


text:
export page, row, cell

function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                                channel::String = Genie.config.webchannels_default_route , head_content::String = "", kwargs...)
  Stipple.Layout.layout(Genie.Renderer.Html.div(id=elemid, args...; kwargs...), partial=partial, title=title, class=class,
                        style=style, head_content=head_content, channel=channel)
end

function row(args...; kwargs...)
  kwargs = NamedTuple(Dict(kwargs...), :class, "row")
  Genie.Renderer.Html.div(args...; kwargs...)
end

function cell(args...; size::Int=0, kwargs...)
  kwargs = NamedTuple(Dict(kwargs...), :class, "col col-12 col-sm$(size > 0 ? "-$size" : "")")
  Genie.Renderer.Html.div(args...; kwargs...)
end


text:
"""
    TestState

                          +->Passed+--+
                          |           |
None+-->Queued+-->Running+-->Failed+--+
           ^              |           |
           |              +->Errored+-+
           |                          |
           +--------------------------+
"""
@enum TestState Passed Failed Errored Queued Running Initial
@with_kw mutable struct TestResult
    name::String
    state::TestState = Initial
    output::String = ""
end
TestResult(name::String) = TestResult(name = name)

function runtests(pkgfolder)::TestResult
    prevpkg = Pkg.project().path

    # Activate test subpackage
    testfolder = joinpath(pkgfolder, "test")
    @suppress Pkg.activate(testfolder)
    @suppress Pkg.instantiate()

    # Run tests for package and capture output
    state = Passed
    output = @capture_out try
        @suppress_err include(joinpath(testfolder, "runtests.jl"))
    catch e
        if (e isa LoadError) && (occursin("did not pass", string(e.error)))
            state = Failed
        else
            state = Errored
        end
        println(e)
    end

    # Restore active project
    @suppress Pkg.activate(prevpkg)

    return TestResult(
        splitpath(pkgfolder)[end],
        state,
        output,
    )
end


text:
"""
    BuoyancyTracer <: AbstractBuoyancyModel{Nothing}

Type indicating that the tracer `b` represents buoyancy.
"""
struct BuoyancyTracer <: AbstractBuoyancyModel{Nothing} end

const BuoyancyTracerModel = Buoyancy{<:BuoyancyTracer}

required_tracers(::BuoyancyTracer) = (:b,)

@inline buoyancy_perturbation(i, j, k, grid, ::BuoyancyTracer, C) = @inbounds C.b[i, j, k]

@inline ∂x_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂xᶠᵃᵃ(i, j, k, grid, C.b)
@inline ∂y_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂yᵃᶠᵃ(i, j, k, grid, C.b)
@inline ∂z_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂zᵃᵃᶠ(i, j, k, grid, C.b)


text:
function nfactory()
    o1 = construct_organ(params=Params(y_EC_ECT=0.8, y_EN_ENT=0.8),
               vars=Vars(assimilation=NitrogenVars()));;
    p1 = o1.params
    o2 = construct_organ(params=Params(y_EC_ECT=0.8, y_EN_ENT=0.8))
    p2 = o2.params;
    u1 = [9.0,8.0,7.0,6.0,5.0,4.0]u"mol"
    u2 = [9.0,8.0,7.0,6.0,5.0,4.0]u"mol"
    u1 = LVector{eltype(u1),typeof(u1),DynamicEnergyBudgets.STATE}(u1)
    u2 = LVector{eltype(u2),typeof(u2),DynamicEnergyBudgets.STATE}(u2)
    u2 .*= 2.7
    du1 = fill(0.0u"mol*hr^-1", 6)
    du2 = fill(0.0u"mol*hr^-1", 6)
    du1 = LVector{eltype(du1),typeof(du1),DynamicEnergyBudgets.STATE}(du1)
    du2 = LVector{eltype(du2),typeof(du2),DynamicEnergyBudgets.STATE}(du2)
    set_var!(o1.vars, :scale, scaling(o1.params.scaling, u1[:V]))
    set_var!(o2.vars, :scale, scaling(o2.params.scaling, u2[:V]))
    f = NAssim();

    o1, p1, u1, du1, o2, p2, u2, du2, f
end

function cfactory()
    o1 = construct_organ(params=Params(y_EC_ECT = 0.8, y_EN_ENT = 0.8),
               vars=Vars(assimilation=CarbonVars()));;
    p1 = o1.params;
    o2 = construct_organ(params=Params(y_EC_ECT = 0.8, y_EN_ENT = 0.8));
    p2 = o2.params;
    u1 = [9.0,8.0,7.0,6.0,5.0,4.0]u"mol"
    u2 = [9.0,8.0,7.0,6.0,5.0,4.0]u"mol"
    u1 = LVector{eltype(u1),typeof(u1),DynamicEnergyBudgets.STATE}(u1)
    u2 = LVector{eltype(u2),typeof(u2),DynamicEnergyBudgets.STATE}(u2)
    u2 .*= 1.9
    du1 = fill(0.0u"mol*hr^-1", 6)
    du2 = fill(0.0u"mol*hr^-1", 6)
    du1 = LVector{eltype(du1),typeof(du1),DynamicEnergyBudgets.STATE}(du1)
    du2 = LVector{eltype(du2),typeof(du2),DynamicEnergyBudgets.STATE}(du2)
    set_var!(o1.vars, :scale, scaling(o1.params.scaling, u1[:V]))
    set_var!(o2.vars, :scale, scaling(o2.params.scaling, u2[:V]))
    f = KooijmanSLAPhotosynthesis();

    o1, p1, u1, du1, o2, p2, u2, du2, f
end

@testset "N assimilation" begin
    global o1, p1, u1, du1, o2, p2, u2, du2, f = nfactory();
    uptake_nitrogen(f, o1, u1)

    @testset "N assimilation depends on soil N concentration sublinealy" begin
        global a = uptake_nitrogen(f, o1, u1)
        # increase soil N
        o1.vars.assimilation.X_NO *= 2
        global b = uptake_nitrogen(f, o1, u1)
        # More N, more uptake
        @test a < b
        # But double N has less than double uptake
        @test 2a > b
    end

    @testset "N assimilation depends on root surface scaling" begin
        global a = uptake_nitrogen(f, o1, u1)
        o1.vars.scale /= 2
        global b = uptake_nitrogen(f, o1, u1)
        @test a > b
    end

    # @testset "N assimilation depends inversely on shoot surface scaling, but maybe shouldn't" begin
        # Shoot scaling inversley affects water uptake, which affects N concentration.
        # The larger the shoot is, the smaller its surface scaling coefficient.
        # This means that except at the very start of growth, less N is assimilated by
        # the roots the larger the shoot gets. Absolutely, not realtively. That doesn't make sense.

        # a = uptake_nitrogen(f, o1, u1)
        # o2.vars.scale /= 2
        # b = uptake_nitrogen(f, o1, u1)
        # @test a < b
    # end

    @testset "N assimilation depends on structure linearly (ignoring scaling)" begin
        global a = uptake_nitrogen(f, o1, u1)
        u1[:V] *= 2
        global b = uptake_nitrogen(f, o1, u1)
        # Lower scaling, less uptake
        @test 2a == b
    end

    @testset "N assimilation flux is added correctly" begin
        @test uptake_nitrogen(f, o1, u1) > 0.0u"μmol/s"
        assimilation!(f, o1, u1)
        sum_flux!(du1, o1, 0)
        sum_flux!(du2, o2, 0)
        global c1 = sum(du1)
        global c2 = sum(du2)
        # @test n1 == o1.J[:N,:ass]
        # @test n1 == uptake_nitrogen(f, o1, u1)
    end

    @testset "N assimilation flux is merged correctly" begin

        set_rej!(o1, o2) = begin
            c_trans = 1.9oneunit(o1.J1[:1,:1])
            o1.J[:C,:tra] = o2.J[:C,:tra] = c_trans
            o1.J[:C,:rej] = o2.J[:C,:rej] = -c_trans
        end

        global o1a, p1a, u1a, du1a, o2a, p2a, u2a, du2a, fa = nfactory();
        set_rej!(o1a, o2a)
        sum_flux!(du1a, o1a, 0); sum_flux!(du2a, o2a, 0)

        global o1b, p1b, u1b, du1b, o2b, p2b, u2b, du2b, fa = nfactory();
        set_rej!(o1b, o2b)
        sum_flux!(du1b, o1b, 0); sum_flux!(du2b, o2b, 0)

        @test o1a.J[:C,:tra] == o1b.J[:C,:tra]
        @test du1a[:C] == du1b[:C]
        @test du1a[:N] == du1b[:N]
        @test du1a[:E] == du1b[:E]
        @test sum(du2a) == sum(du2b)

        global uptake_n = uptake_nitrogen(f, o1b, u1b)
        assimilation!(f, o1b, u1b)
        sum_flux!(du1b, o1b, 0); sum_flux!(du2b, o2b, 0)

        @test du1a[:C] > du1b[:C]
        @test du1a[:N] < du1b[:N]
        @test du1a[:E] < du1b[:E]

        # Actual C and N calcs
        global c1a, c2a, c1b, c2b = sum.((du1a, du2a, du1b, du2b))
        global n1a, n2a, n1b, n2b = dot.((du1a, du2a, du1b, du2b), (n_ratios,))

        global c_loss = o1b.J1[:C,:los] + o2b.J1[:C,:los]
        global n_loss = o1b.J1[:N,:los] + o2b.J1[:N,:los]

        @test c1a > c1b
        @test n1a < n1b

        @testset "N assimilation should balance" begin
            @test -c_loss != zero(c_loss)
            @test c1b + c2b ≈ -c_loss
            @test_broken upreferred(n1b + n2b) ≈ upreferred(uptake_n) - upreferred(n_loss)
        end
    end
end

@testset "C assimilation" begin
    global o1, p1, u1, du1, o2, p2, u2, du2, f = cfactory();
    photosynthesis(f, o1, u1)

    @testset "photosynthesis depends on air C concentration sub-linealy" begin
        global a = photosynthesis(f, o1, u1);
        # increase soil C
        o1.vars.assimilation.X_C /= 2
        global b = photosynthesis(f, o1, u1)
        # Less C, less uptake
        @test a > b
        # But half C does not half uptake
        @test a < 2b
    end

    @testset "photosynthesis depends inversely on air O concentration, sub-linealy" begin
        global a = photosynthesis(f, o1, u1);
        # increase soil C
        o1.vars.assimilation.X_O /= 2
        global b = photosynthesis(f, o1, u1)
        # Less C, less uptake
        @test a < b
        # But half C does not half uptake
        @test 2a > b
    end

    @testset "photosynthesis depends on Light concentration sub-linealy" begin
        global a = photosynthesis(f, o1, u1);
        # increase soil C
        o1.vars.assimilation.J_L_F /= 2
        global b = photosynthesis(f, o1, u1)
        # Less C, less uptake
        @test a > b
        # But half C does not half uptake
        @test a < 2b
    end

    @testset "photosynthesis depends on shoot scaling" begin
        global a = photosynthesis(f, o1, u1)
        o1.vars.scale /= 10
        global b = photosynthesis(f, o1, u1)
        @test a > b
    end

    # @testset "C assimilation does not depend on root scaling, but maybe it should" begin
    #     # Most stomatatal conducance models incorporate root water uptake.
    #     # Less water means closed stomata, and less CO2.
    #     # That is not considered in this model.

    #     a = photosynthesis(f, o1, u1)
    #     o2.vars.scale /= 2
    #     b = photosynthesis(f, o1, u1)
    #     @test a == b
    # end

    @testset "photosynthesis depends on structure linearly (ignoring scaling)" begin
        global a = photosynthesis(f, o1, u1)
        u1[:V] *= 2
        global b = photosynthesis(f, o1, u1)
        # Lower scaling, less uptake
        @test 2a == b
    end

    @testset "C assimilation flux is added correctly" begin
        @test photosynthesis(f, o1, u1) > 0.0u"μmol/s"
        assimilation!(f, o1, u1)
        sum_flux!(du1, o1, 0)
        sum_flux!(du2, o1, 0)
        global c1 = sum(du1)
        global c2 = sum(du2)
        @test c1 == o1.J[:C,:ass]
        @test c1 == photosynthesis(f, o1, u1)
    end

    @testset "C assimilation flux is merged correctly" begin

        global o1a, p1a, u1a, du1a, o2a, p2a, u2a, du2a, fa = cfactory();
        global n_trans = 1.9oneunit(o1a.J1[:1,:1])
        o1a.J[:N,:tra] = o2a.J[:N,:tra] = o1b.J[:N,:tra] = o2b.J[:N,:tra] = n_trans
        o1a.J[:N,:rej] = o2a.J[:N,:rej] = o1b.J[:N,:rej] = o2b.J[:N,:rej] = -n_trans
        sum_flux!(du1a, o1a, 0); sum_flux!(du2a, o2a, 0)

        global o1b, p1b, u1b, du1b, o2b, p2b, u2b, du2b, fa = cfactory();
        sum_flux!(du1b, o1b, 0); sum_flux!(du2b, o2b, 0)

        @test du2a[:C] == du2b[:C]
        @test du2a[:N] == du2b[:N]
        @test du2a[:E] == du2b[:E]

        # assimilation only on o1b
        global uptake_c = photosynthesis(f, o1b, u1b)
        assimilation!(KooijmanSLAPhotosynthesis(), o1b, u1b)
        upreferred(o1b.J1[:N,:los])
        o1b.J
        sum_flux!(du1b, o1b, 0); sum_flux!(du2b, o2b, 0)
        du1a
        du1b

        @test du1a[:C] < du1b[:C]
        @test du1a[:N] > du1b[:N]
        @test du1a[:E] < du1b[:E]

        global c1a, c2a, c1b, c2b = sum.((du1a, du2a, du1b, du2b))
        global n1a, n2a, n1b, n2b = dot.((du1a, du2a, du1b, du2b), (n_ratios,))
        global c_loss = o1b.J1[:C,:los] + o2b.J1[:C,:los]
        global n_loss = o1b.J1[:N,:los] + o2b.J1[:N,:los]

        @testset "C assimilation should have gained total cmols" begin
            @test c1a < c1b
            @test c2a == c2b
        end

        @testset "C assimilation should balance" begin
            @test upreferred(c1b + c2b) ≈ upreferred(-c_loss + uptake_c)
            @test upreferred(n1b + n2b) ≈ upreferred(-n_loss)
        end

    end
end
nothing