text:
# # Advection-diffusion on a Sphere
# ## Setup
# - Passive advection of an anomaly using non-divergent horizontal flow on the sphere (currently 300km/s but switching to Williamson stream function)
# - Homogeneous tracer concentration of the outer shell couples to an inhomogeneous concentration of the inner shell

# ## Equations
# The RHS evaluation follows:
# ```math
# \frac{\partial \theta}{\partial t} = - \nabla \cdot (\vec{u} \theta + \kappa \nabla \theta) \qquad
# ```
# where
#  -  $\theta$ is the tracer (e.g. potential temperature)
#  -  $\vec{u}(u,0,0)$ is a constant zonal velocity
#  -  $\kappa$ is the diffusivity tensor.
# **Boundary Conditions**
#  - External boundary conditions (Neumann):
# ```math
# \vec{n} \cdot \kappa \nabla \theta = 0 \qquad \text{at } z = 4\text{km}, -4\text{km}
# ```
# where $\vec{n}$ is the upward pointing unit normal vector at the boundary.
#  - Primary coupled boundary condition (forced by secondary domain):
# ```math
# \vec{n} \cdot \kappa_a \nabla \theta =f_a= \lambda (\theta - \theta_{secondary}) \qquad \text{at atmos } z = 0
# ```
# where $\lambda$ is a constant but will be replaced by a roughness coefficient factor in more complex test cases.
# $f_a$ and $\kappa_a$ are, respectively, the flux and diffusivity of the atmosphere component.
#  - Secondary coupled boundary condition (accumulated flux by primary domain or prescribed flux):
# If accumulated,
# ```math
# \vec{n} \cdot \kappa_o\nabla \theta = f_o \qquad \text{at ocean } z = 0 \text{ (or pseudo subsurface)}
# and
# ```
# ```math
# f_o = f_{accum} = \int_{\tau_0}^{\tau} f_a dt
# ```
# where $τ-τ_0$ gives the coupling cycle length.
# $f_o$ and $\kappa_o$ are, respectively, the flux and diffusivity of the atmosphere component, with
# $f_{accum}$ being the accumulated flux.

# ```
# ATMOS
# ------------- θ @z=0 [lowest atmos level]
#   A
#  / \          f = n⋅κ_a ∇θ = λ(θ−θ_secondary)
#   |
# ------------- θ_secondary @z=0 [highest ocean level]
# OCEAN
# ```

# # Import packages
using ClimateMachine, MPI
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods

using ClimateMachine.ODESolvers
using ClimateMachine.ODESolvers: solve!

using ClimateMachine.Atmos: SphericalOrientation, latitude, longitude, AtmosFilterPerturbations

using CLIMAParameters
using CLIMAParameters.Planet: MSLP, R_d, day, grav, Omega, planet_radius

using ClimateMachine.Mesh.Filters

using ClimateMachine.BalanceLaws:
    BalanceLaw,
    AbstractStateType,
    Prognostic,
    Auxiliary,
    Gradient,
    GradientFlux,
    GradientLaplacian,
    Hyperdiffusive,
    UpwardIntegrals,
    DownwardIntegrals,
    vars_state,
    number_states

using ClimaCoupler

using Unitful
using Dates: DateTime
using Statistics
using StaticArrays

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

ClimateMachine.init()
FT = Float64

## Shared functions
include("domains.jl")
include("abstractions.jl")
include("callbacks.jl")

## Main balance law and its components
include("CplMainBL.jl")

# # Set simulation parameters
nstepsA = 10
nstepsO = 5

##  Background atmos and ocean diffusivities
const κᵃʰ = FT(1e4) * 0.0
const κᵃᶻ = FT(1e-1)
const κᵒʰ = FT(1e3) * 0.0
const κᵒᶻ = FT(1e-4)
const τ_airsea = FT(60 * 86400)
const L_airsea = FT(500)
const λ_airsea = FT(L_airsea / τ_airsea)
coupling_lambda() = (λ_airsea)
const u_max = FT(1e-5) # max. advective velocity in radians

# # Set up coupled model
# Define component models and initialize the coupler

function main(::Type{FT}) where {FT}

    ## Domain - Spherical Shells
    ## confusing name - better might be to use something like DeepSphericalShellDomain directly?
    ΩO = DeepSphericalShellDomain(radius = FT(planet_radius(param_set)) - FT(4e3), height = FT(4e3))
    ΩA = DeepSphericalShellDomain(radius = FT(planet_radius(param_set)), height = FT(4e3))

    ## Grid
    nelem = (; horizontal = 8, vertical = 4)
    polynomialorder = (; horizontal = 5, vertical = 5)
    overintegrationorder = (; horizontal = 1, vertical = 1)

    gridA = DiscontinuousSpectralElementGrid(ΩA, nelem, polynomialorder)
    gridO = DiscontinuousSpectralElementGrid(ΩO, nelem, polynomialorder)

    ## Numerics-specific options
    numerics = (
        NFfirstorder = CentralNumericalFluxFirstOrder(),
        NFsecondorder = PenaltyNumFluxDiffusive(),
        overint_params = (overintegrationorder, polynomialorder),
    )

    ## Timestepping
    Δt_ = calculate_dt(gridA, wavespeed = u_max * (ΩA.radius), diffusivity = maximum([κᵃʰ, κᵃᶻ]))

    t_time, end_time = (0, 20Δt_)

    ## Collect spatial info, timestepping, balance law and DGmodel for the two components
    boundary_mask(xc, yc, zc) = @. (xc^2 + yc^2 + zc^2)^0.5 ≈ planet_radius(param_set)
    #+
    # **1.** Atmos component

    ## Prop atmos functions to override defaults
    atmos_structure(λ, ϕ, r) = FT(30)#30.0 + 10.0 * cos(ϕ) * sin(5λ)
    atmos_θⁱⁿⁱᵗ(npt, el, x, y, z) = atmos_structure(lon(x, y, z), lat(x, y, z), rad(x, y, z))                # Set atmosphere initial state function
    ##atmos_θ_shadowflux(θᵃ, θᵒ, npt, el, xc, yc, zc) = FT(0.0)
    atmos_θ_shadowflux(θᵃ, θᵒ, npt, el, xc, yc, zc) = is_surface(xc, yc, zc) ? (1.0 / τ_airsea) * (θᵃ - θᵒ) : 0.0 # Set atmosphere shadow boundary flux function
    atmos_calc_kappa_diff(_...) = κᵃʰ, κᵃʰ, κᵃᶻ               # Set atmos diffusion coeffs
    atmos_source_θ(θᵃ, npt, el, xc, yc, zc, θᵒ) = FT(0.0)     # Set atmos source!
    atmos_get_penalty_tau(_...) = FT(3.0 * 0.0)               # Set penalty term tau (for debugging)

    ## Set atmos advective velocity (constant in time) and convert to Cartesian
    uˡᵒⁿ(λ, ϕ, r) = u_max * r * cos(ϕ)
    atmos_uⁱⁿⁱᵗ(npt, el, x, y, z) =
        (0 * r̂(x, y, z) + 0 * ϕ̂(x, y, z) + uˡᵒⁿ(lon(x, y, z), lat(x, y, z), rad(x, y, z)) * λ̂(x, y, z))

    ## Collect atmos props
    bl_propA = prop_defaults()
    bl_propA = (;
        bl_propA...,
        init_theta = atmos_θⁱⁿⁱᵗ,
        theta_shadow_boundary_flux = atmos_θ_shadowflux,
        calc_kappa_diff = atmos_calc_kappa_diff,
        source_theta = atmos_source_θ,
        get_penalty_tau = atmos_get_penalty_tau,
        coupling_lambda = coupling_lambda,
        init_u = atmos_uⁱⁿⁱᵗ,
    )
    ## Setup atmos component model
    mA = CplModel(;
        grid = gridA,
        equations = CplMainBL(bl_propA, (CoupledPrimaryBoundary(), ExteriorBoundary()), param_set),
        boundary_z = boundary_mask,
        nsteps = nstepsA,
        dt = Δt_ / nstepsA,
        timestepper = LSRK54CarpenterKennedy,
        numerics...,
    )
    #+
    # **2.** Ocean component
    ## Prop ocean functions to override defaults
    tropical_heating_1(λ, ϕ, r) = 30.0 + 10.0 * cos(ϕ) * sin(5λ)
    tropical_heating_2(λ, ϕ, r) = 30.0 + 10.0 * cos(ϕ) + 1 * sin(5λ) * cos(ϕ)
    tropical_heating(λ, ϕ, r) = tropical_heating_1(λ, ϕ, r)
    ocean_θⁱⁿⁱᵗ(npt, el, x, y, z) = tropical_heating(lon(x, y, z), lat(x, y, z), rad(x, y, z))                    # Set ocean initial state function
    ocean_calc_kappa_diff(_...) = κᵒʰ, κᵒʰ, κᵒᶻ               # Set ocean diffusion coeffs
    ocean_source_θ(θᵃ, npt, el, xc, yc, zc, θᵒ) = FT(0.0)     # Set ocean source!
    ocean_get_penalty_tau(_...) = FT(0.15 * 0.0)               # Set penalty term tau (for debugging)
    ocean_uⁱⁿⁱᵗ(xc, yc, zc, npt, el) = SVector(FT(0.0), FT(0.0), FT(0.0)) # Set ocean advective velocity

    ## Collect ocean props
    bl_propO = prop_defaults()
    bl_propO = (;
        bl_propO...,
        init_theta = ocean_θⁱⁿⁱᵗ,
        calc_kappa_diff = ocean_calc_kappa_diff,
        source_theta = ocean_source_θ,
        get_penalty_tau = ocean_get_penalty_tau,
        coupling_lambda = coupling_lambda,
        init_u = ocean_uⁱⁿⁱᵗ,
    )
    ## Setup ocean component model
    mO = CplModel(;
        grid = gridO,
        equations = CplMainBL(bl_propO, (ExteriorBoundary(), CoupledSecondaryBoundary()), param_set),
        boundary_z = boundary_mask,
        nsteps = nstepsO,
        dt = Δt_ / nstepsO,
        timestepper = LSRK54CarpenterKennedy,
        numerics...,
    )
    #+
    # Create the coupler object for holding import/export fields and performs mappings
    # and instantiate the coupled timestepper:
    coupler = CplState()
    coupler_register!(coupler, :Ocean_SST, deepcopy(mO.state.θ[mO.boundary]), mO.grid, DateTime(0), u"°C")
    coupler_register!(
        coupler,
        :Atmos_MeanAirSeaθFlux,
        deepcopy(mA.state.F_accum[mA.boundary]),
        mA.grid,
        DateTime(0),
        u"°C",
    )

    # Instantiate a coupled timestepper that steps forward the components and
    # implements mapings between components export bondary states and
    # other components imports.
    compA = (pre_step = preatmos, component_model = mA, post_step = postatmos)
    compO = (pre_step = preocean, component_model = mO, post_step = postocean)
    component_list = (atmosphere = compA, ocean = compO)
    cpl_solver = CplSolver(component_list = component_list, coupler = coupler, coupling_dt = Δt_, t0 = 0.0)

    #+
    # Creat callbacks and bundle the simulation:
    ## For now applying callbacks only to atmos.
    callbacks = (
        ExponentialFiltering(),
        VTKOutput(
            (iteration = string(5Δt_) * "ssecs", overdir = "output", overwrite = true, number_sample_points = 0)...,
        ),
    )

    simulation = (;
        coupled_odesolver = cpl_solver,
        odesolver = cpl_solver.component_list.atmosphere.component_model.stepper,
        state = cpl_solver.component_list.atmosphere.component_model.state,
        dgmodel = cpl_solver.component_list.atmosphere.component_model.discretization,
        callbacks = callbacks,
        simtime = (t_time, end_time),
        name = "Coupler_UnitTest_ctrl",
    )

    return simulation
end

function run(cpl_solver, numberofsteps, cbvector)
    ## Run the model
    solve!(nothing, cpl_solver; numberofsteps = numberofsteps, callbacks = cbvector)
end

# # Run simulation

simulation = main(Float64);
nsteps = Int(simulation.simtime[2] / simulation.coupled_odesolver.dt)
cbvector = create_callbacks(simulation, simulation.odesolver)
println("Initialized. Running...")
@time run(simulation.coupled_odesolver, nsteps, cbvector)

max_stars_repo_name: CliMA/ClimaCoupler


text:
module Density

# This implements different methods to compute density
#
# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using ..PhaseDiagrams
#using ..MaterialParameters
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
import Base.show

abstract type AbstractDensity <: AbstractMaterialParam end

export  ComputeDensity,         # calculation routines
        ComputeDensity!,        # in place calculation
        ConstantDensity,        # constant
        PT_Density              # P & T dependent density

# Constant Density -------------------------------------------------------
"""
    ConstantDensity(ρ=2900kg/m^3)

Set a constant density:
```math
    \\rho  = cst
```
where ``\\rho`` is the density [``kg/m^3``].
"""
@with_kw_noshow mutable struct ConstantDensity <: AbstractDensity
    equation::LaTeXString   =   L"\rho = cst"
    ρ::GeoUnit              =   2900kg/m^3                # density
end

# Calculation routines
function ComputeDensity(P,T, s::ConstantDensity)
    @unpack ρ   = s
    if length(T)>1
        return Value(ρ).*ones(size(T))
    else
        return ρ*1.0
    end
end

function ComputeDensity!(rho::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::ConstantDensity) where N
    @unpack ρ   = s

    rho .= NumValue(ρ)

    return nothing
end

# Print info
function show(io::IO, g::ConstantDensity)
    print(io, "Constant density: ρ=$(g.ρ.val)")
end
#-------------------------------------------------------------------------

# Pressure & Temperature dependent density -------------------------------
"""
    PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-9/Pa, T0=0C, P=0MPa)

Set a pressure and temperature-dependent density:
```math
    \\rho  = \\rho_0 (1.0 - \\alpha (T-T_0) + \\beta  (P-P_0) )
```
where ``\\rho_0`` is the density [``kg/m^3``] at reference temperature ``T_0`` and pressure ``P_0``,
``\\alpha`` is the temperature dependence of density and ``\\beta`` the pressure dependence.

"""
@with_kw_noshow mutable struct PT_Density <: AbstractDensity
    equation::LaTeXString   =   L"\rho = \rho_0(1.0-\alpha (T-T_0) + \beta (P-P_0)"
    ρ0::GeoUnit             =   2900kg/m^3                  # density
    α::GeoUnit              =   3e-5/K                      # T-dependence of density
    β::GeoUnit              =   1e-9/Pa                     # P-dependence of density
    T0::GeoUnit             =   0C                          # Reference temperature
    P0::GeoUnit             =   0MPa                        # Reference pressure
end

# Calculation routine
function ComputeDensity(P,T, s::PT_Density)
    @unpack ρ0,α,β,P0, T0   = s

    ρ = ρ0*(1.0 - α*(T - T0) + β*(P - P0) )

    return ρ
end

function ComputeDensity!(ρ, P,T, s::PT_Density)
    @unpack ρ0,α,β,P0, T0   = s

    ρ[:] = ρ0*(1.0 - α*(T-T0) + β*(P-P0) )

    return ρ
end

function ComputeDensity!(ρ::AbstractArray{<:AbstractFloat},P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PT_Density)
    @unpack ρ0,α,β,P0, T0   = s
    ρ0 = NumValue(ρ0)
    α  = NumValue(α)
    β  = NumValue(β)
    P0 = NumValue(P0)
    T0 = NumValue(T0)

    ρ  .= ρ0*(1.0 .- α*( T .- T0) + β*(P .- P0) )

    return nothing
end

# Print info
function show(io::IO, g::PT_Density)
    print(io, "P/T-dependent density: ρ0=$(g.ρ0.val), α=$(g.α.val), β=$(g.β.val), T0=$(g.T0.val), P0=$(g.P0.val)")
end
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
# Phase diagrams
"""
    ComputeDensity(P,T, s::PhaseDiagram_LookupTable)

Interpolates density as a function of `T,P`
"""
function ComputeDensity(P,T, s::PhaseDiagram_LookupTable)
    return s.Rho.(T,P)
end

"""
    ComputeDensity!(rho::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PhaseDiagram_LookupTable)

In-place computation of density as a function of `T,P`, in case we are using a lookup table.
"""
function ComputeDensity!(rho::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, s::PhaseDiagram_LookupTable)
    rho[:] = s.Rho.(T,P)
    return nothing
end

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

"""
    ComputeDensity!(rho::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `Phase` of every point is specified as an Integer in the `Phases` array.

# Example
```julia
julia> MatParam    =   Array{MaterialParams, 1}(undef, 2);
julia> MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=1,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = PerpleX_LaMEM_Diagram("test_data/Peridotite.in"));
julia> MatParam[2] =   SetMaterialParams(Name="Crust", Phase=2,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = ConstantDensity(ρ=2900kg/m^3));
julia> Phases = ones(Int64,400,400);
julia> Phases[:,20:end] .= 2
julia> rho     = zeros(size(Phases))
julia> T       =  ones(size(Phases))
julia> P       =  ones(size(Phases))*10
julia> ComputeDensity!(rho, Phases, P,T, MatParam)
julia> rho
400×400 Matrix{Float64}:
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46  …  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46  …  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
    ⋮                                            ⋮     ⋱                     ⋮
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46  …  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 3334.46  3334.46  3334.46  3334.46  3334.46  3334.46     2900.0  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
```
The routine is made to minimize allocations:
```julia
julia> julia> @time ComputeDensity!(rho, Phases, P,T, MatParam)
0.003121 seconds (49 allocations: 3.769 MiB)
```
"""
function ComputeDensity!(rho::AbstractArray{<:AbstractFloat, N}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where N

    for i = 1:length(MatParam)

        if !isnothing(MatParam[i].Density)
            # Create views into arrays (so we don't have to allocate)
            ind = Phases .== i;
            rho_local   =   view(rho, ind )
            P_local     =   view(P  , ind )
            T_local     =   view(T  , ind )

            ComputeDensity!(rho_local, P_local, T_local, MatParam[i].Density[1] )
        end

    end

end

"""
    ComputeDensity!(rho::AbstractArray{<:AbstractFloat,N}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)

"""
function ComputeDensity!(rho::AbstractArray{<:AbstractFloat, N}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {N,M}

    if M!=(N+1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    rho .= 0.0;
    for i = 1:length(MatParam)

        rho_local   = zeros(size(rho))
        Fraction    = selectdim(PhaseRatios,M,i);
        if (maximum(Fraction)>0.0) & (!isnothing(MatParam[i].Density))

            ComputeDensity!(rho_local, P, T, MatParam[i].Density[1] )

            rho .= rho .+ rho_local.*Fraction
        end

    end

end

end

ext: jl


text:
module MutateModule

import ..CoreModule: Node, copy_node, Options, Dataset, RecordType
import ..EquationUtilsModule: compute_complexity, count_constants, count_depth
import ..LossFunctionsModule: score_func, score_func_batch
import ..CheckConstraintsModule: check_constraints
import ..PopMemberModule: PopMember
import ..MutationFunctionsModule:
    gen_random_tree_fixed_size,
    mutate_constant,
    mutate_operator,
    append_random_op,
    prepend_random_op,
    insert_random_op,
    delete_random_op,
    crossover_trees
import ..SimplifyEquationModule: simplify_tree, combine_operators
import ..RecorderModule: @recorder

# Go through one simulated options.annealing mutation cycle
#  exp(-delta/T) defines probability of accepting a change
function next_generation(
    dataset::Dataset{T},
    baseline::T,
    member::PopMember,
    temperature::T,
    curmaxsize::Int,
    frequencyComplexity::AbstractVector{T},
    options::Options;
    tmp_recorder::RecordType,
)::Tuple{PopMember,Bool,Float64} where {T<:Real}
    prev = member.tree
    parent_ref = member.ref
    tree = prev
    mutation_accepted = false
    num_evals = 0.0

    #TODO - reconsider this
    if options.batching
        beforeScore, beforeLoss = score_func_batch(dataset, baseline, prev, options)
        num_evals += (options.batchSize / dataset.n)
    else
        beforeScore = member.score
        beforeLoss = member.loss
    end

    nfeatures = dataset.nfeatures

    mutationChoice = rand()
    #More constants => more likely to do constant mutation
    weightAdjustmentMutateConstant = min(8, count_constants(prev)) / 8.0
    cur_weights = copy(options.mutationWeights) .* 1.0
    cur_weights[1] *= weightAdjustmentMutateConstant
    n = compute_complexity(prev, options)
    depth = count_depth(prev)

    # If equation too big, don't add new operators
    if n >= curmaxsize || depth >= options.maxdepth
        cur_weights[3] = 0.0
        cur_weights[4] = 0.0
    end
    cur_weights /= sum(cur_weights)
    cweights = cumsum(cur_weights)

    successful_mutation = false
    #TODO: Currently we dont take this \/ into account
    is_success_always_possible = true
    attempts = 0
    max_attempts = 10

    #############################################
    # Mutations
    #############################################
    while (!successful_mutation) && attempts < max_attempts
        tree = copy_node(prev)
        successful_mutation = true
        if mutationChoice < cweights[1]
            tree = mutate_constant(tree, temperature, options)
            @recorder tmp_recorder["type"] = "constant"

            is_success_always_possible = true
            # Mutating a constant shouldn't invalidate an already-valid function

        elseif mutationChoice < cweights[2]
            tree = mutate_operator(tree, options)

            @recorder tmp_recorder["type"] = "operator"

            is_success_always_possible = true
            # Can always mutate to the same operator

        elseif mutationChoice < cweights[3]
            if rand() < 0.5
                tree = append_random_op(tree, options, nfeatures)
                @recorder tmp_recorder["type"] = "append_op"
            else
                tree = prepend_random_op(tree, options, nfeatures)
                @recorder tmp_recorder["type"] = "prepend_op"
            end
            is_success_always_possible = false
            # Can potentially have a situation without success
        elseif mutationChoice < cweights[4]
            tree = insert_random_op(tree, options, nfeatures)
            @recorder tmp_recorder["type"] = "insert_op"
            is_success_always_possible = false
        elseif mutationChoice < cweights[5]
            tree = delete_random_op(tree, options, nfeatures)
            @recorder tmp_recorder["type"] = "delete_op"
            is_success_always_possible = true
        elseif mutationChoice < cweights[6]
            tree = simplify_tree(tree, options) # Sometimes we simplify tree
            tree = combine_operators(tree, options) # See if repeated constants at outer levels
            @recorder tmp_recorder["type"] = "partial_simplify"
            mutation_accepted = true
            return (
                PopMember(
                    tree,
                    beforeScore,
                    beforeLoss;
                    parent=parent_ref,
                    deterministic=options.deterministic,
                ),
                mutation_accepted,
                num_evals,
            )

            is_success_always_possible = true
            # Simplification shouldn't hurt complexity; unless some non-symmetric constraint
            # to commutative operator...

        elseif mutationChoice < cweights[7]
            # Sometimes we generate a new tree completely tree
            # We select a random size, though the generated tree
            # may have fewer nodes than we request.
            tree_size_to_generate = rand(1:curmaxsize)
            tree = gen_random_tree_fixed_size(tree_size_to_generate, options, nfeatures)
            @recorder tmp_recorder["type"] = "regenerate"

            is_success_always_possible = true
        else # no mutation applied
            @recorder begin
                tmp_recorder["type"] = "identity"
                tmp_recorder["result"] = "accept"
                tmp_recorder["reason"] = "identity"
            end
            mutation_accepted = true
            return (
                PopMember(
                    tree,
                    beforeScore,
                    beforeLoss;
                    parent=parent_ref,
                    deterministic=options.deterministic,
                ),
                mutation_accepted,
                num_evals,
            )
        end

        successful_mutation =
            successful_mutation && check_constraints(tree, options, curmaxsize)

        attempts += 1
    end
    #############################################

    if !successful_mutation
        @recorder begin
            tmp_recorder["result"] = "reject"
            tmp_recorder["reason"] = "failed_constraint_check"
        end
        mutation_accepted = false
        return (
            PopMember(
                copy_node(prev),
                beforeScore,
                beforeLoss;
                parent=parent_ref,
                deterministic=options.deterministic,
            ),
            mutation_accepted,
            num_evals,
        )
    end

    if options.batching
        afterScore, afterLoss = score_func_batch(dataset, baseline, tree, options)
        num_evals += (options.batchSize / dataset.n)
    else
        afterScore, afterLoss = score_func(dataset, baseline, tree, options)
        num_evals += 1
    end

    if isnan(afterScore)
        @recorder begin
            tmp_recorder["result"] = "reject"
            tmp_recorder["reason"] = "nan_loss"
        end
        mutation_accepted = false
        return (
            PopMember(
                copy_node(prev),
                beforeScore,
                beforeLoss;
                parent=parent_ref,
                deterministic=options.deterministic,
            ),
            mutation_accepted,
            num_evals,
        )
    end

    probChange = 1.0
    if options.annealing
        delta = afterScore - beforeScore
        probChange *= exp(-delta / (temperature * options.alpha))
    end
    if options.useFrequency
        oldSize = compute_complexity(prev, options)
        newSize = compute_complexity(tree, options)
        old_frequency = (oldSize <= options.maxsize) ? frequencyComplexity[oldSize] : 1e-6
        new_frequency = (newSize <= options.maxsize) ? frequencyComplexity[newSize] : 1e-6
        probChange *= old_frequency / new_frequency
    end

    if probChange < rand()
        @recorder begin
            tmp_recorder["result"] = "reject"
            tmp_recorder["reason"] = "annealing_or_frequency"
        end
        mutation_accepted = false
        return (
            PopMember(
                copy_node(prev),
                beforeScore,
                beforeLoss;
                parent=parent_ref,
                deterministic=options.deterministic,
            ),
            mutation_accepted,
            num_evals,
        )
    else
        @recorder begin
            tmp_recorder["result"] = "accept"
            tmp_recorder["reason"] = "pass"
        end
        mutation_accepted = true
        return (
            PopMember(
                tree,
                afterScore,
                afterLoss;
                parent=parent_ref,
                deterministic=options.deterministic,
            ),
            mutation_accepted,
            num_evals,
        )
    end
end

"""Generate a generation via crossover of two members."""
function crossover_generation(
    member1::PopMember,
    member2::PopMember,
    dataset::Dataset{T},
    baseline::T,
    curmaxsize::Int,
    options::Options,
)::Tuple{PopMember,PopMember,Bool,Float64} where {T<:Real}
    tree1 = member1.tree
    tree2 = member2.tree
    crossover_accepted = false

    # We breed these until constraints are no longer violated:
    child_tree1, child_tree2 = crossover_trees(tree1, tree2)
    num_tries = 1
    max_tries = 10
    num_evals = 0.0
    while true
        # Both trees satisfy constraints
        if check_constraints(child_tree1, options, curmaxsize) &&
            check_constraints(child_tree2, options, curmaxsize)
            break
        end
        if num_tries > max_tries
            crossover_accepted = false
            return member1, member2, crossover_accepted, num_evals  # Fail.
        end
        child_tree1, child_tree2 = crossover_trees(tree1, tree2)
        num_tries += 1
    end
    if options.batching
        afterScore1, afterLoss1 = score_func_batch(dataset, baseline, child_tree1, options)
        afterScore2, afterLoss2 = score_func_batch(dataset, baseline, child_tree2, options)
        num_evals += 2 * (options.batchSize / dataset.n)
    else
        afterScore1, afterLoss1 = score_func(dataset, baseline, child_tree1, options)
        afterScore2, afterLoss2 = score_func(dataset, baseline, child_tree2, options)
        num_evals += options.batchSize / dataset.n
    end

    baby1 = PopMember(
        child_tree1,
        afterScore1,
        afterLoss1;
        parent=member1.ref,
        deterministic=options.deterministic,
    )
    baby2 = PopMember(
        child_tree2,
        afterScore2,
        afterLoss2;
        parent=member2.ref,
        deterministic=options.deterministic,
    )

    crossover_accepted = true
    return baby1, baby2, crossover_accepted, num_evals
end

end


text:
using WoodburyMatrices
using LinearAlgebra, SparseArrays, Test
using Random: seed!

include("woodbury.jl")
include("symwoodbury.jl")


text:
# System tests.
module TestSystem

using Orchestra.Types
using Orchestra.System
using Orchestra.Transformers
importall Orchestra.Util

include("fixture_learners.jl")
using .FixtureLearners
nfcp = NumericFeatureClassification()
fcp = FeatureClassification()

function all_concrete_subtypes(a_type::Type)
  a_subtypes = Type[]
  for a_subtype in subtypes(a_type)
    if isleaftype(a_subtype)
      push!(a_subtypes, a_subtype)
    else
      append!(a_subtypes, all_concrete_subtypes(a_subtype))
    end
  end
  return a_subtypes
end

concrete_learner_types = setdiff(
  all_concrete_subtypes(Learner),
  all_concrete_subtypes(TestLearner)
)

using FactCheck

facts("Orchestra system") do
  context("All learners train and predict on fixture data.") do
    for concrete_learner_type in concrete_learner_types
      learner = concrete_learner_type()
      fit_and_transform!(learner, nfcp)
    end

    @fact 1 => 1
  end

  context("All learners train and predict on iris dataset.") do
    # Get data
    dataset = readcsv(joinpath(dirname(@__FILE__), "iris.csv"))
    instances = dataset[:,1:(end-1)]
    labels = dataset[:, end]
    (train_ind, test_ind) = holdout(size(instances, 1), 0.3)
    train_instances = instances[train_ind, :]
    test_instances = instances[test_ind, :]
    train_labels = labels[train_ind]
    test_labels = labels[test_ind]

    # Test all learners
    for concrete_learner_type in concrete_learner_types
      learner = concrete_learner_type()
      fit!(learner, train_instances, train_labels)
      transform!(learner, test_instances)
    end

    @fact 1 => 1
  end

  context("Ensemble with learners from different libraries work.") do
    learners = Learner[]
    push!(learners, RandomForest())
    push!(learners, StackEnsemble())
    if LIB_SKL_AVAILABLE
      push!(learners, SKLLearner())
    end
    if LIB_CRT_AVAILABLE
      push!(learners, CRTLearner())
    end
    ensemble = VoteEnsemble({:learners => learners})
    predictions = fit_and_transform!(ensemble, nfcp)

    @fact 1 => 1
  end

  context("Pipeline works with fixture data.") do
    transformers = [
      OneHotEncoder(),
      Imputer(),
      StandardScaler(),
      BestLearner()
    ]
    pipeline = Pipeline({:transformers => transformers})
    predictions = fit_and_transform!(pipeline, fcp)

    @fact 1 => 1
  end
end

end # module


text:
function solve_with_loop(input)
	pc = 1
	visited = falses(length(input))
	ret = 0
	while !visited[pc]
		visited[pc] = true
		instruction, value = split(input[pc])
		if instruction == "acc"
			ret += parse(Int64, value)
			pc += 1
		elseif instruction == "jmp"
			pc += parse(Int64, value)
		elseif instruction == "nop"
			pc += 1
		end
	end
	ret
end

function solve_fix_loop(input)
	pc = saved_pc = 1
	visited = saved_visited = falses(length(input))
	ret = saved_acc = 0
	fix = fixable = true
	while pc <= length(input)
		if !visited[pc]
			visited[pc] = true
			instruction, value = split(input[pc])
			if instruction == "acc"
				ret += parse(Int64, value)
				pc += 1
			elseif instruction == "jmp"
				if fixable
					if fix
						fixable = false
						saved_pc = pc
						saved_visited = copy(visited)
						saved_acc = ret
						pc += 1
					else
						pc += parse(Int64, value)
						fix = true
					end
				else
					pc += parse(Int64, value)
				end
			elseif instruction == "nop"
				if fixable
					if fix
						fixable = false
						saved_pc = pc
						saved_visited = copy(visited)
						saved_acc = ret
						pc += parse(Int64, value)
					else
						fix = true
						pc += 1
					end
				else
					pc += 1
				end
			end
		else
			fixable = true
			pc = saved_pc
			visited = copy(saved_visited)
			visited[pc] = false
			fix = false
			ret = saved_acc
		end
	end
	ret
end

@assert solve_with_loop(readlines("input/test/test_day8.txt")) == 5 ["test case 1 failed"]
@assert solve_fix_loop(readlines("input/test/test_day8.txt")) == 8 ["test case 2 failed"]

println("First star: " * string(solve_with_loop(readlines("input/day8.txt"))))
println("Second star: " * string(solve_fix_loop(readlines("input/day8.txt"))))


text:
import ImageMagick
using HTTP
using Test
using FileIO

if !isdir("data")
  @info "download data.zip"
  HTTP.open("GET", "http://media.tuhh.de/ibi/MPIReco/data.zip") do http
    open("data.zip", "w") do file
        write(file, http)
    end
  end
  @info "extracting data.zip"
  run(`unzip -oq data.zip`)
  rm("data.zip")
end

mkpath("./img/")

function exportImage(filename, I::AbstractMatrix)
  Iabs = abs.(I)
  Icolored = colorview(Gray, Iabs./maximum(Iabs))
  save(filename, Icolored )
end

include("Cartesian.jl")
include("MotionCompensation.jl")
include("Reconstruction.jl")
include("MultiPatch.jl")
include("MultiGradient.jl")
include("Sparse.jl")
include("SMCenter.jl")
include("LoadSaveMDF.jl")
include("SMRecovery.jl")


text:
using SPDEPrecisionMatrices
using Random
using SparseArrays
using TriangleMesh
using LinearAlgebra
using SpecialFunctions
using Test

@testset "Meshes" begin
    Random.seed!(0)
    n = 200
    points = rand(2, n) .* [100, 50]
    mesh = refine(create_mesh(collect(points')))

    A = observation_matrix(mesh, points)
    @test size(A) == (n, mesh.n_point)
    @test nnz(A) == n * 3
    @test all(sum(A, dims=2) .≈ 1)
end

@testset "Precision matrices" begin
    Random.seed!(0)
    n = 200
    points = rand(2, n) .* [100, 50]
    mesh = refine(create_mesh(collect(points')))

    r = 15 # decorrelation range
    σ = 5 # marginal variance
    ν = 2  # smoothness parameter
    d = size(mesh.point, 1)

    for ν in 1:3
        println("ν = $ν")
        κ = calculate_κ(ν, r)
        τ = calculate_τ(ν, d, κ, σ)
        @test κ ≈ sqrt(8ν) / r
        @test τ ≈ sqrt(gamma(ν) / (gamma(ν + d/2) * (4π)^(d/2))) / (σ * κ^ν)

        Q0 = unscaled_precision_matrix(mesh, κ, ν)
        Q = precision_matrix(mesh, r, σ, ν)
        @test all((Q .- τ^2 * Q0) .< 1e-6)
        @test isposdef(Q0 + 1e-3*I)
        @test isposdef(Q + 1e-3*I)
    end
end


text:
# * Carriers

carrier_types = Dict{Symbol,Any}()

max_frequency(carrier::AbstractCarrier) = frequency(carrier)

# ** Fixed carrier
@doc raw"""
    FixedCarrier <: LinearCarrier

The carrier is fixed in the sense that the instantaneous frequency is
constant throughout the pulse, i.e. no chirp/dispersion:

```math
\Im\{\exp[\im(\omega t + \phi)]\} = \sin(\omega t + \phi).
```
"""
struct FixedCarrier{Λ,Tt,Ω,Φ} <: LinearCarrier
    λ::Λ
    T::Tt
    ω::Ω
    ϕ::Φ # Carrier–envelope phase, in radians
end

(carrier::FixedCarrier)(t) = sin(carrier.ω*t + carrier.ϕ)

carrier_types[:fixed] = FixedCarrier

wavelength(carrier::FixedCarrier) = carrier.λ
period(carrier::FixedCarrier) = carrier.T

frequency(carrier::FixedCarrier) = 1/carrier.T
wavenumber(carrier::FixedCarrier) = 1/carrier.λ
fundamental(carrier::FixedCarrier) = carrier.ω
photon_energy(carrier::FixedCarrier) = carrier.ω
phase(carrier::FixedCarrier) = carrier.ϕ

phase_shift(c::FixedCarrier, δϕ) =
    FixedCarrier(c.λ, c.T, c.ω, c.ϕ+δϕ)

function FixedCarrier(field_params::Dict{Symbol,Any})
    @unpack λ, T, ω = field_params
    ϕ = get(field_params, :ϕ, 0)
    FixedCarrier(λ, T, austrip(ω), ϕ)
end

function show(io::IO, carrier::FixedCarrier)
    printfmt(io, "Fixed carrier @ λ = {1:s} (T = {2:s}, ω = {3:.4f} Ha = {4:s})",
             si_round(u"m"(carrier.λ)),
             si_round(u"s"(carrier.T)),
             carrier.ω, au2si_round(carrier.ω, u"eV"))
    !iszero(carrier.ϕ) &&
        printfmt(io, "; CEP = {1:0.2f}π", carrier.ϕ/π)
end

# ** Transverse carriers

# *** Linear

@doc raw"""
    LinearTransverseCarrier

The carrier is identical to [`FixedCarrier`](@ref), but explicitly in
3d, with the polarization taken along the ``z`` axis:

```math
\vec{C}(t) = \bmat{0\\0\\\sin(\omega t + \phi)}.
```
"""
struct LinearTransverseCarrier{Carrier<:LinearCarrier} <: TransverseCarrier
    carrier::Carrier
end

LinearTransverseCarrier(field_params::Dict{Symbol,Any}) =
    LinearTransverseCarrier(FixedCarrier(field_params))

function show(io::IO, carrier::LinearTransverseCarrier)
    write(io, "LinearTransverseCarrier: ")
    show(io, carrier.carrier)
end

(carrier::LinearTransverseCarrier)(t) = SVector(0, 0, carrier.carrier(t))

carrier_types[:linear] = LinearTransverseCarrier

for fun in [:wavelength, :period, :frequency, :wavenumber, :fundamental, :photon_energy, :phase]
    @eval $fun(carrier::LinearTransverseCarrier) = $fun(carrier.carrier)
end

phase_shift(c::LinearTransverseCarrier, δϕ) =
    LinearTransverseCarrier(phase_shift(c.carrier, δϕ))

# *** Elliptical

@doc raw"""
    EllipticalCarrier

An elliptically polarized field (of which circularly and linearly
polarized are special cases) is given in the canonical coordinate system as

```math
\vec{C}(t) = \frac{1}{\sqrt{1+\xi^2}}
\bmat{\xi\cos\theta\\0\\\sin\theta},
\quad
\theta \defd \omega t + \phi,
```
i.e. with the principal axis of the ellipse along ``z``.
"""
struct EllipticalCarrier{Λ,Tt,Ω,Φ,Ξ} <: TransverseCarrier
    λ::Λ
    T::Tt
    ω::Ω
    ϕ::Φ # Carrier–envelope phase, in radians
    ξ::Ξ # Amount of ellipticity, minor/major axis ratio
end

function EllipticalCarrier(field_params::Dict{Symbol,Any})
    @unpack λ, T, ω = field_params
    ϕ = get(field_params, :ϕ, 0)
    ξ = get(field_params, :ξ, 0)
    EllipticalCarrier(λ, T, austrip(ω), ϕ, ξ)
end

function show(io::IO, carrier::EllipticalCarrier)
    ξ = carrier.ξ
    pol = if iszero(ξ)
        ""
    elseif ξ == 1
        " (RCP)"
    elseif ξ == -1
        " (LCP)"
    elseif ξ > 0
        " (right)"
    else
        " (left)"
    end
    printfmt(io, "Elliptical carrier with ξ = {1:.2f}{2:s} @ λ = {3:s} (T = {4:s}, ω = {5:.4f} Ha = {6:s})",
             ξ, pol,
             si_round(u"m"(carrier.λ)),
             si_round(u"s"(carrier.T)),
             carrier.ω, au2si_round(carrier.ω, u"eV"))
    !iszero(carrier.ϕ) &&
        printfmt(io, "; CEP = {1:0.2f}π", carrier.ϕ/π)
end

function (carrier::EllipticalCarrier)(t)
    ξ = carrier.ξ
    ϕ = carrier.ω*t + carrier.ϕ
    s,c = sincos(ϕ)
    SVector(ξ*c, 0, s)/√(1 + ξ^2)
end

carrier_types[:elliptical] = EllipticalCarrier

wavelength(carrier::EllipticalCarrier) = carrier.λ
period(carrier::EllipticalCarrier) = carrier.T

frequency(carrier::EllipticalCarrier) = 1/carrier.T
wavenumber(carrier::EllipticalCarrier) = 1/carrier.λ
fundamental(carrier::EllipticalCarrier) = carrier.ω
photon_energy(carrier::EllipticalCarrier) = carrier.ω
phase(carrier::EllipticalCarrier) = carrier.ϕ

phase_shift(c::EllipticalCarrier, δϕ) =
    EllipticalCarrier(c.λ, c.T, c.ω, c.ϕ+δϕ, c.ξ)


text:
abstract type AbstractNonLinearLeastSquares end

mutable struct NonLinearLeastSquares{BFVT, BFJVT} <: AbstractNonLinearLeastSquares
    # Basis functions where f[i] ∀ i ∈ [1, k] should be a function
    # which takes the state and independant
    # variable as arguments (i.e f[i](x,tₖ) where tₖ is the time at
    # the kth epoch. For now, must be in vector form even if there
    # is only one nonlinear basis function.
    f::BFVT

    # User provided jacobian of f w.r.t. x where h[i,j] ∀ i ∈ [1, k], j ∈ [1, n]
    # computes the partial derivative of f[i] w.r.t. x[j] and takes the state
    # and independant variabels as arguments(i.e. h[i,j](x, tₖ)). For now
    # must be in matrix form even if there is only one nonlinear basis function.
    h::BFJVT

    # Jacobian matrix
    H::Matrix{Float64}

    # Estimate vector
    xhat::Vector{Float64}

    # Measurements
    y::Vector{Float64}

    # Measurement Weights
    W::Symmetric{Float64, Matrix{Float64}}

    # Independant variables
    t::Vector{Float64}

    # Tolerance
    ϵ::Float64

    # Max iterations
    iMax::Int
end

mutable struct ADNonLinearLeastSquares{BFVT} <: AbstractNonLinearLeastSquares
    # Basis functions where f[i] ∀ i ∈ [0, n] should be a function
    # which takes the state and independant
    # variable as arguments (i.e f[i](x,tₖ) where tₖ is the time at
    # the kth epoch.
    f::BFVT

    # Jacobian matrix
    H::Matrix{Float64}

    # Estimate vector
    xhat::Vector{Float64}

    # Measurements
    y::Vector{Float64}

    # Measurement Weights
    W::Symmetric{Float64, Matrix{Float64}}

    # Independant variables
    t::Vector{Float64}

    # Tolerance
    ϵ::Float64

    # Max iterations
    iMax::Int
end

function NonLinearLeastSquares(f::AbstractVector, h::AbstractMatrix, xhat::AbstractVector;
                               ϵ = 1e-6, iMax = 100)
    k = length(f)
    n = length(xhat)
    kh, nh = size(h)
    if k != kh
        throw(ArgumentError("Provided jacobian does not have the correct number of rows."))
    end
    if n != nh
        throw(ArgumentError("Provided jacobian does not have the correct number of columns."))
    end

    NonLinearLeastSquares(f,h,Matrix{Float64}(undef, (0,0)), xhat, Vector{Float64}(undef, 0),
        Symmetric(Matrix{Float64}(undef, (0,0))), Vector{Float64}(undef, 0), ϵ, iMax)
end

function NonLinearLeastSquares(f::AbstractVector, xhat::AbstractVector; ϵ = 1e-6, iMax = 100)
    ADNonLinearLeastSquares(f, Matrix{Float64}(undef, (0,0)), xhat, Vector{Float64}(undef, 0),
        Symmetric(Matrix{Float64}(undef, (0,0))), Vector{Float64}(undef, 0), ϵ, iMax)
end

function FeedMeasurementBatch!(nlls::AbstractNonLinearLeastSquares, ys::AbstractVector, W::AbstractMatrix, ts::AbstractVector)
    m = length(ys)
    n = length(nlls.xhat)
    if m / length(nlls.f) != length(ts)
        throw(ArgumentError("Measurement vector and indepandant variable vector lengths are not appropriate."))
    end
    r, c = size(W)
    if r != c
        throw(ArgumentError("Weight matrix must be square!"))
    end

    # Set measurements and indepandant variable vectors
    nlls.y = ys
    nlls.W = Symmetric(W)
    nlls.t = ts

    # Size H matrix
    nlls.H = zeros(m, n)

    return nothing
end

function FeedMeasurementBatch!(nlls::AbstractNonLinearLeastSquares, ys::AbstractVector, ts::AbstractVector)
    m = length(ys)
    W = zeros(m,m)
    for i in 1:m; W[i,i] = 1.0; end
    FeedMeasurementBatch!(nlls, ys, W, ts)
end

function ComputeEstimate!(nlls::AbstractNonLinearLeastSquares)

    # Get requirements
    conv = false
    n  = length(nlls.xhat)
    ml = length(nlls.t)
    m  = length(nlls.y)
    p  = length(nlls.f)
    iter = 0

    # Allocate memory
    Jco  = Inf
    Δyc  = zeros(m)
    Δx   = zeros(n)
    fxc  = zeros(m)
    tv1  = zeros(m)
    tv2  = zeros(n)
    tm1  = zeros(m,n)
    tm2  = zeros(n,n)
    while conv == false && iter < nlls.iMax
        # Increment iteration counter
        iter += 1

        # Fill H
        fillH!(nlls)

        # Fill expected measurements
        @inbounds for i in 1:ml
            for j in 1:p
                fxc[(i-1)*p + j] = nlls.f[j](nlls.xhat, nlls.t[i])
            end
        end

        # Compute residual and local cost function
        Δyc .= nlls.y .- fxc
        mul!(tv1, nlls.W, Δyc)
        Jc = dot(Δyc, tv1)

        # Compute change in local cost function and print
        δJc = abs(Jc - Jco)
        Jco = Jc
        println("Iteration #: " * string(iter) *
                " δJ: " * string(δJc))

        # Check for convergence
        if δJc < nlls.ϵ / norm(nlls.W)
            conv = true
        else
            # Compute correction
            mul!(tv2, transpose(nlls.H), tv1)
            mul!(tm1, nlls.W, nlls.H)
            mul!(tm2, transpose(nlls.H), tm1)
            fm = factorize(tm2)
            ldiv!(Δx, fm, tv2)

            # Update
            nlls.xhat .+= Δx
        end
    end

    return conv
end

function fillH!(nlls::NonLinearLeastSquares)
    n  = length(nlls.xhat)
    ml = length(nlls.t)
    p  = length(nlls.f)

    # Fill
    @inbounds for i in 1:ml
        for j in 1:p
            for k in 1:n
                nlls.H[(i - 1)*p + j,k] = nlls.h[j,k](nlls.xhat, nlls.t[i])
            end
        end
    end

    return nothing
end

function fillH!(nlls::ADNonLinearLeastSquares)
    n  = length(nlls.xhat)
    ml = length(nlls.t)
    p  = length(nlls.f)

    # Fill
    @inbounds for i in 1:ml
        for j in 1:p
            grad = gradient(x -> nlls.f[j](x, nlls.t[i]), nlls.xhat)[1]
            for k in 1:n
                nlls.H[(i - 1)*p + j,k] = grad[k]
            end
        end
    end

    return nothing
end

function GetMeasurementEstimate(nlls::AbstractNonLinearLeastSquares, t::AbstractFloat)
    return [nlls.f[i](nlls.xhat, t) for i in 1:length(nlls.f)]
end

function GetEstimate(nlls::AbstractNonLinearLeastSquares)
    return nlls.xhat
end


text:
using Defaults, Test
using Defaults: get_default

@default_kw struct TestStruct
    foo::Int    | 1
    bar::Symbol | :bar
end

@test get_default(TestStruct, Val{:foo}) == 1
@test get_default(TestStruct, Val{:bar}) == :bar

t = TestStruct()
@test t.foo == 1
@test t.bar == :bar

t = TestStruct(foo = 2, bar=:foo)
@test t.foo == 2
@test t.bar == :foo

@default_kw struct TestStruct
    foo::Int    | 3
    bar::Symbol | :foobar
end

t = TestStruct()
@test t.foo == 3
@test t.bar == :foobar


text:
abstract type Sheeps{L} end

function setNSheep!( sheep::Sheeps, i, j, n )
    oldSheep = _setNSheep!( sheep, i, j, n )
    sheep.current_sheep[] += n - oldSheep
    return nothing
end # function

function kickSheep!( sheep::Sheeps )
    while !isempty(sheep.diffusion_candidates)
        coords = popfirst!( sheep.diffusion_candidates )
        n_sheep = getNSheep( sheep, coords... )
        if n_sheep <= sheep.capacity
            continue
        end # if
        diffuseSheep!( sheep, coords... )
    end # while
end # function

function diffuseSheep!( sheep::Sheeps{L}, i, j ) where L
    #TODO: improve this (more generic, iterators)
    n_sheep = getNSheep( sheep, i, j )
    if n_sheep == 0
        return nothing
    end # if
    setNSheep!( sheep, i, j, n_sheep - 1 )
    r = rand(1:4)
    if r == 1
        i += 1
        if i > L
            i = 1
        end # if
    elseif r == 2
        if i == 1
            i = L
        end # if
        i -= 1
    elseif r == 3
        j += 1
        if j > L
            j = 1
        end # if
    elseif r == 4
        if j == 1
            j = L
        end # if
        j -= 1
    end # if

    new_sheep = getNSheep( sheep, i, j )
    setNSheep!( sheep, i, j, new_sheep + 1 )
    return nothing
end # function


text:
#=

  Labeled dice and Building blocks problems in Julia ConstraintSolver.jl

  Labeled dice
  --------------
  From Jim Orlin "Colored letters, labeled dice: a logic puzzle"
  http://jimorlin.wordpress.com/2009/02/17/colored-letters-labeled-dice-a-logic-puzzle/
  """
  My daughter Jenn bough a puzzle book, and showed me a cute puzzle.  There
  are 13 words as follows:  BUOY, CAVE, CELT, FLUB, FORK, HEMP, JUDY,
  JUNK, LIMN, QUIP, SWAG, VISA, WISH.

  There are 24 different letters that appear in the 13 words.  The question
  is:  can one assign the 24 letters to 4 different cubes so that the
  four letters of each word appears on different cubes.  (There is one
  letter from each word on each cube.)  It might be fun for you to try
  it.  I'll give a small hint at the end of this post. The puzzle was
  created by Humphrey Dudley.
  """

  Building blocks
  ---------------
  From http://brownbuffalo.sourceforge.net/BuildingBlocksClues.html
  """
  Each of four alphabet blocks has a single letter of the alphabet on each
  of its six sides. In all, the four blocks contain every letter but
  Q and Z. By arranging the blocks in various ways, you can spell all of
  the words listed below. Can you figure out how the letters are arranged
  on the four blocks?

  BAKE ONYX ECHO OVAL
  GIRD SMUG JUMP TORN
  LUCK VINY LUSH WRAP
  """

  Note: This is a somewhat generalized model for solving both
        Building blocks and Labeled Dice problems.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function block_problem(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 3, # default, 10
                                                            "activity.max_confidence_deviation" => 3, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    # a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z = 1:26
    println("problem ", problem[:name])
    num_cubes = problem[:num_cubes]
    num_sides = problem[:num_sides]
    letters = split(problem[:letters],"")
    num_chars = length(letters)

    # letters = problem[:letters]
    words = problem[:words]
    num_words = length(words)

    d = Dict(zip(letters,1:length(letters)))
    d_rev = [k=>v for (v,k) in d] # reversed dict

    cube_len = num_cubes*num_sides
    @variable(model, 1 <= cube[1:cube_len] <= num_cubes, Int)

    # each letters in a word must be on a different die
    for i in 1:num_words
      word = [d[w] for w in split(words[i,:][1],"")]
      this_die = @variable(model, [1:length(word)],CS.Integers(1:num_cubes))
      for i in 1:length(word)
          my_element(model,word[i],cube,this_die[i])
      end
      @constraint(model, this_die in CS.AllDifferent())
    end

    # there must be exactly NumSides (6) letters of each die (num_cubes)
    global_cardinality_count(model, cube, [num_sides for _ in 1:num_cubes])

    # simple symmetry breaking: place first letter on cube 1
    @constraint(model, cube[1] == 1)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("\nsolution #$sol")
                cube_val = convert.(Integer,JuMP.value.(cube; result=sol))
                println("cube:$cube_val")
                println("words:")
                for word in words
                    print(word,": ")
                    for w in split(word,"")
                        dd = d[w]
                        print("$w: $(cube_val[dd]) ")
                    end
                    println()
                end
                println("\ndice:")
                dice = [ [] for _ in 1:num_cubes]
                for (c,i) in collect(zip(letters,cube_val))
                    push!(dice[i],c)
                end
                join.(sort.(unique.(dice))," ").|>println
            end
        end
    else
        println("status:$status")
    end

    return status
end

function block_problems(problem)
    all_block_problems = Dict(
      :labeled_dice =>  Dict(
           :name => "labeled_dice",
           # number of cubes
           :num_cubes => 4,
           # number of sides of of a cube
           :num_sides => 6,
           # the letters to use
           :letters => "abcdefghijklmnopqrstuvwy",
           # the words to place
           :words => ["buoy", "cave", "celt", "flub",
                      "fork", "hemp", "judy","junk",
                      "limn", "quip", "swag", "visa",
                      "wish"]
    ),

    :building_blocks => Dict(
           :name => "building_blocks",
           :num_cubes => 4,
           :num_sides => 6,
           # Note: 'f' is not in any word but it needed to get 24 chars
           :letters => "abcdefghijklmnoprstuvwxy",
           :words => ["bake", "onyx", "echo", "oval",
                      "gird", "smug", "jump", "torn",
                      "luck", "viny", "lush", "wrap"]
        )

   )

   return all_block_problems[problem]

end

@time block_problem(block_problems(:labeled_dice),true,true)
println()
@time block_problem(block_problems(:building_blocks),true,true)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Plots, OrdinaryDiffEq, DiffEqPhysics, RecursiveArrayTools
gr()

G = 2.95912208286e-4
M = [1.00000597682, 0.000954786104043, 0.000285583733151, 0.0000437273164546, 0.0000517759138449, 1/1.3e8]
planets = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"]

pos_x = [0.0,-3.5023653,9.0755314,8.3101420,11.4707666,-15.5387357]
pos_y = [0.0,-3.8169847,-3.0458353,-16.2901086,-25.7294829,-25.2225594]
pos_z = [0.0,-1.5507963,-1.6483708,-7.2521278,-10.8169456,-3.1902382]
pos = ArrayPartition(pos_x,pos_y,pos_z)

vel_x = [0.0,0.00565429,0.00168318,0.00354178,0.00288930,0.00276725]
vel_y = [0.0,-0.00412490,0.00483525,0.00137102,0.00114527,-0.00170702]
vel_z = [0.0,-0.00190589,0.00192462,0.00055029,0.00039677,-0.00136504]
vel = ArrayPartition(vel_x,vel_y,vel_z)

tspan = (0.,200_000)

const ∑ = sum
const N = 6
potential(p, t, x, y, z, M) = -G*∑(i->∑(j->(M[i]*M[j])/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2), 1:i-1), 2:N)

nprob = NBodyProblem(potential, M, pos, vel, tspan)
sol = solve(nprob,Yoshida6(), dt=100);

orbitplot(sol,body_names=planets)

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
abstract type Dataset end

struct Lee08 <: Dataset end

function raw_table(::Lee08)
    _path = joinpath(dirname(@__FILE__), "..", "data", "lee08.csv")
    matrix, header = readdlm(_path, ',', Float64; header=true)
    columntable( (voteshare=matrix[:,1], margin = matrix[:,2]))
end

function RDData(dataset::Lee08)
    lee08 = raw_table(dataset)
    ZsR = RunningVariable(lee08.margin ./ 100; cutoff = 0.0, treated = :≥)
    Ys = lee08.voteshare ./ 100
    RDData(Ys, ZsR)
end

#oreopoulos_matrix, orepoulos_header = readdlm(oreopoulos_path, ',', Float64; header=true)
# bla =


text:
using HDF5,JLD,KUparser
ftstr = length(ARGS) >= 1 ? "$(ARGS[1])" : "zn11orig"
ptstr = length(ARGS) >= 2 ? "$(ARGS[2])" : "ArcEager13"
@show (ftstr, ptstr)
ft = eval(parse("Flist.$ftstr"))
pt = eval(parse(ptstr))
@date dev=load("acl11.dev.jld4") # corpus, deprel, postag
@date tst=load("acl11.tst.jld4")
@date trn=load("acl11.trn.jld4")
@assert trn["deprel"] == dev["deprel"] == tst["deprel"]
@assert trn["postag"] == dev["postag"] == tst["postag"]
nd = length(dev["deprel"])
@date (pdev,xdev,ydev) = oparse(pt, dev["corpus"], nd, ft)
@date (ptst,xtst,ytst) = oparse(pt, tst["corpus"], nd, ft)
@date (ptrn,xtrn,ytrn) = oparse(pt, trn["corpus"], nd, ft)
if isa(xtrn, SparseMatrixCSC)
    nrows = max(maximum(xtrn.rowval), maximum(xtst.rowval), maximum(xdev.rowval))
    xtrn.m = xtst.m = xdev.m = nrows
end
@save "$ftstr.jld" xtrn ytrn xdev ydev xtst ytst


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

# # Natural frequency analysis of a formula frame

# The model is a 3d formula frame shown in picture. Geometry is from Formula
# Studen Oulu. Material and cross-section values used in this example
# aren't the same ones used in the real frame. Boundary conditions may also
# differ from the analysis done by Formula Student Oulu.

# ![](f_frame/model.png)

using JuliaFEM
using FEMBase.Test
using Logging
Logging.configure(level=INFO)
add_elements! = JuliaFEM.add_elements!

# Creating geometry.

L1 = 270; L2 = 150; L3 = 400; H = 530 # [mm]

X = Dict(1 => [ 0.0,  0.0,  0.0],
         2 => [L1+L2, 0.0,  0.0],
         3 => [L1,    H,    0.0],
         4 => [L1+L2+L3, H, 0.0])

# Creating five beam elements

beam1 = Element(Seg2, [1,2])
beam2 = Element(Seg2, [1,3])
beam3 = Element(Seg2, [2,3])
beam4 = Element(Seg2, [3,5])
beam5 = Element(Seg2, [2,5])

beam_elements = [beam1, beam2, beam3, beam4, beam5]
info("Number of elements: ", length(beam_elements))

# Defining material values for beam elements. In this case structural streel
# values are used.
update!(beam_elements, "geometry", X)
update!(beam_elements, "youngs modulus", 210.0e3)
update!(beam_elements, "shear modulus", 80.77e3)
update!(beam_elements, "density", 7.800e-9)

# Defining cross-section values. We have to define cross-section area,
# moment of inertias in both local directions and polar moment of inertia.
# In this case values for 25mm x 2.5 circular tube (pipe) are used.

# To calculate inertias see:
# https://calcresource.com/moment-of-inertia-ctube.html
# For more information about moment of inertias:
# https://en.wikipedia.org/wiki/List_of_second_moments_of_area

update!(beam_elements, "cross-section area", 176.715)
update!(beam_elements, "torsional moment of inertia 1", 11320.778)
update!(beam_elements, "torsional moment of inertia 2", 11320.778)
update!(beam_elements, "polar moment of inertia", 22641.556)

# The direction of beam is defined in the same way as in ABAQUS.
# That is, we have a tangent direction and one normal direction.
# The third direction is then cross product of tangent and normal.
# This is a little shortcut to define all the beam orientations with one for
# loop. Because all the elements are pipes, we only need to define a vector
# which is orthogonal to tangent of the beam. This for loop finds the right
# vector for every beam element.

update!(beam_elements, "normal", [0.0, 0.0, -1.0])

# Create boundary conditions: fix all degrees of freedom for nodes in
# a set FIXED. Here we first create elements of type `Poi1` for each
# node j in set FIXED, update geometry field and then create new fields
# `fixed displacmeent 1`, `fixed displacement 2`, and so on, where the
# displacement / rotation is prescribed.

# Create a problem, containing beam elements and boundary conditions:

frame = Problem(Beam, "3d frame", 6)
add_elements!(frame, beam_elements)
add_elements!(frame, bc_elements)

# Perform modal analysis

step = Analysis(Modal)
step.properties.nev = 10
xdmf = Xdmf(joinpath(datadir, "f_frame_results"); overwrite=true)
add_results_writer!(step, xdmf)
add_problems!(step, [frame])
run!(step)
close(xdmf.hdf)

# Each `Analysis` can have properties, e.g. time, maximum number of iterations,
# convergence tolerance and so on. Eigenvalues of calculation are stored as a
# properties of analysis:

eigvals = step.properties.eigvals
freqs = sqrt(eigvals)/(2*pi)

println("Five lowest natural frequencies [Hz]:")
println(round.(freqs[1:5], 2))

# [![mode5](3d_frame/natfreq.png)](https://www.youtube.com/watch?v=GzktCqeASmo)


text:
abstract type AbstractSampler
end

struct GibbsSampler <: AbstractSampler
   nsteps::Int
end

struct TemperatureDrivenSampler <: AbstractSampler
   nsteps::Int
   autocorcoeff::Float64
   betasd::Float64
end

function TemperatureDrivenSampler(;
      nsteps::Int = 5, autocorcoeff::Float64 = 0.9,
      betasd::Float64 = 0.9)

   TemperatureDrivenSampler(nsteps, autocorcoeff, betasd)
end

struct NoSampler <: AbstractSampler
end

function sample!(particles::Particles, sampler::NoSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0)
   # do nothing
end

function sample!(vmodel::M, hmodel::M, sampler::NoSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0
      ) where {M <:AbstractArray{Float64}}
   # do nothing
end

function sample!(vmodel::M, hmodel::M, sampler::AbstractSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0
      ) where {M <:AbstractArray{Float64}}

   particles = Particles(2)
   particles[1] = vmodel
   particles[2] = hmodel
   sample!(particles, sampler, rbm, upfactor, downfactor)
   nothing
end

function sample!(vmodel::M, hmodel::M, sampler::GibbsSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0
      ) where {M <:AbstractArray{Float64}}

   for i = 1:sampler.nsteps
      samplevisible!(vmodel, rbm, hmodel, downfactor)
      samplehidden!(hmodel, rbm, vmodel, upfactor)
   end
   nothing
end

function sample!(particles::Particles, sampler::GibbsSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0)

   gibbssample!(particles, rbm, sampler.nsteps)
end

function sample!(particles::Particles, sampler::TemperatureDrivenSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0)
   # TODO respect upfactor and downfactor
   gibbssample_gamma!(particles, rbm, sampler.nsteps;
         autocorcoeff = sampler.autocorcoeff, betasd = sampler.betasd)
end


text:
module LinearSolvers

using LinearAlgebra
using Printf
using SparseArrays

import Base: show

using CUDA
using KernelAbstractions
import CUDA.CUBLAS
import CUDA.CUSOLVER
import CUDA.CUSPARSE
import IterativeSolvers
import Krylov
import LightGraphs
import Metis

import ..ExaPF: xnorm, csclsvqr!

const KA = KernelAbstractions

export bicgstab, list_solvers
export DirectSolver, BICGSTAB, EigenBICGSTAB, KrylovBICGSTAB
export get_transpose

@enum(
    SolveStatus,
    Unsolved,
    MaxIterations,
    NotANumber,
    Converged,
    Diverged,
)

include("preconditioners.jl")
include("bicgstab.jl")
include("bicgstab_eigen.jl")

abstract type AbstractLinearSolver end
abstract type AbstractIterativeLinearSolver <: AbstractLinearSolver end

"""
    list_solvers(::KernelAbstractions.Device)

List linear solvers available on current device. Currently,
only `CPU()` and `CUDADevice()` are supported.

"""
function list_solvers end

get_transpose(::AbstractLinearSolver, M::AbstractMatrix) = transpose(M)

"""
    ldiv!(solver, x, A, y)
    ldiv!(solver, x, y)

* `solver::AbstractLinearSolver`: linear solver to solve the system
* `x::AbstractVector`: Solution
* `A::AbstractMatrix`: Input matrix
* `y::AbstractVector`: RHS

Solve the linear system ``A x = y`` using the algorithm
specified in `solver`. If `A` is not specified, the function
will used directly the factorization stored inplace.

"""
function ldiv! end

"""
    rdiv!(solver, x, A, y)
    rdiv!(solver, x, y)

* `solver::AbstractLinearSolver`: linear solver to solve the system
* `x::AbstractVector`: Solution
* `A::AbstractMatrix`: Input matrix
* `y::AbstractVector`: RHS

Solve the linear system ``A^⊤ x = y`` using the algorithm
specified in `solver`. If `A` is not specified, the function
will used directly the factorization stored inplace.

"""
function rdiv! end

"""
    DirectSolver <: AbstractLinearSolver

Solve linear system ``A x = y`` with direct linear algebra.

* On the CPU, `DirectSolver` uses UMFPACK to solve the linear system
* On CUDA GPU, `DirectSolver` redirects the resolution to the method `CUSOLVER.csrlsvqr`

"""
struct DirectSolver{Fac<:Union{Nothing, LinearAlgebra.Factorization}} <: AbstractLinearSolver
    factorization::Fac
end

DirectSolver(J::SparseMatrixCSC) = DirectSolver(lu(J))
DirectSolver() = DirectSolver(nothing)
DirectSolver(precond::AbstractPreconditioner) = DirectSolver(nothing)

# Reuse factorization in update
function ldiv!(s::DirectSolver{<:LinearAlgebra.Factorization}, y::AbstractVector, J::AbstractMatrix, x::AbstractVector)
    lu!(s.factorization, J) # Update factorization inplace
    LinearAlgebra.ldiv!(y, s.factorization, x) # Forward-backward solve
    return 0
end
# Solve system Ax = y
function ldiv!(s::DirectSolver{<:LinearAlgebra.Factorization}, y::AbstractArray, x::AbstractArray)
    LinearAlgebra.ldiv!(y, s.factorization, x) # Forward-backward solve
    return 0
end
# Solve system A'x = y
function rdiv!(s::DirectSolver{<:LinearAlgebra.Factorization}, y::Array, x::Array)
    LinearAlgebra.ldiv!(y, s.factorization', x) # Forward-backward solve
    return 0
end

function ldiv!(::DirectSolver{Nothing}, y::Vector, J::AbstractMatrix, x::Vector)
    y .= J \ x
    return 0
end

function ldiv!(::DirectSolver{Nothing},
    y::CUDA.CuVector, J::CUSPARSE.CuSparseMatrixCSR, x::CUDA.CuVector,
)
    CUSOLVER.csrlsvqr!(J, x, y, 1e-8, one(Cint), 'O')
    return 0
end
function ldiv!(::DirectSolver{Nothing},
    y::CUDA.CuVector, J::CUSPARSE.CuSparseMatrixCSC, x::CUDA.CuVector,
)
    csclsvqr!(J, x, y, 1e-8, one(Cint), 'O')
    return 0
end
get_transpose(::DirectSolver, M::CUSPARSE.CuSparseMatrixCSR) = CUSPARSE.CuSparseMatrixCSC(M)

function update!(solver::AbstractIterativeLinearSolver, J)
    update(J, solver.precond)
end

"""
    BICGSTAB <: AbstractIterativeLinearSolver
    BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)

Custom BICGSTAB implementation to solve iteratively the linear system
``A  x = y``.
"""
struct BICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    maxiter::Int
    tol::Float64
    verbose::Bool
end
BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false) = BICGSTAB(precond, maxiter, tol, verbose)
function ldiv!(solver::BICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P

    y[:], n_iters, status = bicgstab(J, x, P, y; maxiter=solver.maxiter,
                                        verbose=solver.verbose, tol=solver.tol)
    if status != Converged
        @warn("BICGSTAB failed to converge. Final status is $(status)")
    end
    return n_iters
end

"""
    EigenBICGSTAB <: AbstractIterativeLinearSolver
    EigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)

Julia's port of Eigen's BICGSTAB to solve iteratively the linear system
``A x = y``.
"""
struct EigenBICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    maxiter::Int
    tol::Float64
    verbose::Bool
end
EigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false) = EigenBICGSTAB(precond, maxiter, tol, verbose)
function ldiv!(solver::EigenBICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P

    y[:], n_iters, status = bicgstab_eigen(J, x, P, y; maxiter=solver.maxiter,
                                            verbose=solver.verbose, tol=solver.tol)
    if status != Converged
        error("EigenBICGSTAB failed to converge. Final status is $(status)")
    end

    return n_iters
end

struct RefBICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    verbose::Bool
end
RefBICGSTAB(precond; verbose=true) = RefBICGSTAB(precond, verbose)
function ldiv!(solver::RefBICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    y[:], history = IterativeSolvers.bicgstabl(P*J, P*x, log=solver.verbose)
    return history.iters
end

struct RefGMRES <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    restart::Int
    verbose::Bool
end
RefGMRES(precond; restart=4, verbose=true) = RefGMRES(precond, restart, verbose)
function ldiv!(solver::RefGMRES,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    y[:], history = IterativeSolvers.gmres(P*J, P*x, restart=solver.restart, log=solver.verbose)
    return history.iters
end

struct DQGMRES <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    memory::Int
    verbose::Bool
end
DQGMRES(precond; memory=4, verbose=false) = DQGMRES(precond, memory, verbose)
function ldiv!(solver::DQGMRES,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    (y[:], status) = Krylov.dqgmres(J, x, N=P, memory=solver.memory)
    return length(status.residuals)
end

"""
    KrylovBICGSTAB <: AbstractIterativeLinearSolver
    KrylovBICGSTAB(precond; verbose=0, rtol=1e-10, atol=1e-10)

Wrap `Krylov.jl`'s BICGSTAB algorithm to solve iteratively the linear system
``A x = y``.
"""
struct KrylovBICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    verbose::Int
    atol::Float64
    rtol::Float64
end
KrylovBICGSTAB(precond; verbose=0, rtol=1e-10, atol=1e-10) = KrylovBICGSTAB(precond, verbose, atol, rtol)
function ldiv!(solver::KrylovBICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    (y[:], status) = Krylov.bicgstab(J, x, N=P,
                                     atol=solver.atol,
                                     rtol=solver.rtol,
                                     verbose=solver.verbose)
    return length(status.residuals)
end

"""
    list_solvers(::KA.CPU)

List all linear solvers available solving the power flow on the CPU.
"""
list_solvers(::KA.CPU) = [RefBICGSTAB, RefGMRES, DQGMRES, BICGSTAB, EigenBICGSTAB, DirectSolver, KrylovBICGSTAB]

"""
    list_solvers(::KA.CUDADevice)

List all linear solvers available solving the power flow on an NVIDIA GPU.
"""
list_solvers(::KA.CUDADevice) = [BICGSTAB, DQGMRES, EigenBICGSTAB, DirectSolver, KrylovBICGSTAB]

end


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

JLLWrappers.@generate_main_file_header("samtools")
JLLWrappers.@generate_main_file("samtools", UUID("cafa7ba0-c216-50a9-a84c-e0b2e0e60c48"))
end  # module samtools_jll


text:
"""
Implementation of `Paxos.Transports.Common.Transport` for TCP-based communication
between participants.
"""
module TCP

export tcp

import Sockets

using Serialization

using ..Common
using ...Utils

struct TCPTransport <: Transport
end

struct TCPConnection <: Connection
    socket::Sockets.TCPSocket
end

function tcp()
    TCPTransport()
end

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

"""
Create a connection (if one is possible) that enables sending and receiving
messages to the recipient. The provided handler should take 1 argument,
a `Connection` for sending or receiving messages with the recipient.
"""
function Common.connectTo(transport::TCPTransport, recipient)
    host, port = recipient
    @debug "Connecting to $recipient over TCP"
    TCPConnection(Sockets.connect(host, port))
end

"""
Listen for incoming connections and invoke the indicated function when they appear.
Returns a `Listener`
"""
function Common.listenOn(handler::Function, transport::TCPTransport, address)
    @debug "Beginning to listen for TCP connections on $address"
    host, port = address
    @sync try
        server = Sockets.listen(port)
        @debug "Waiting for TCP connections"
        while true
            client = Sockets.accept(server)
            @async handler(TCPConnection(client))
        end
        @debug "Finished waiting for TCP connections"
    finally
        close(server)
       @debug "Finished listening for TCP connections on $address"
    end
end

"""
Send a message on the indicated connection, using the specified Transport
"""
function Common.sendTo(connection::TCPConnection, message)
    serialize(connection.socket, message)
end

"""
Receive a message over the indicated connection
"""
function Common.receiveFrom(connection::TCPConnection)
    try
        deserialize(connection.socket)
    catch ex
        if isa(ex, EOFError)
            @debug "Connection closed"
        else
            rethrow()
        end
    end
end

end


text:
function pick_integer_bytes(W, preferred)
    # SIMD quadword integer support requires AVX512DQ
    if !AVX512DQ
        preferred = min(4, preferred)
    end
    max(1,min(preferred, prevpow2(REGISTER_SIZE ÷ W)))
end
function integer_of_bytes(bytes)
    if bytes == 8
        Int64
    elseif bytes == 4
        Int32
    elseif bytes == 2
        Int16
    elseif bytes == 1
        Int8
    else
        throw("$bytes is an invalid number of bytes for integers.")
    end
end
function pick_integer(W, pref)
    integer_of_bytes(pick_integer_bytes(W, pref))
end
pick_integer(::Val{W}) where {W} = pick_integer(W, sizeof(Int))

@generated function vrange(::Val{W}, ::Type{T}, ::Val{O}, ::Val{F}) where {W,T,O,F}
    if T <: Integer
        _T2 = pick_integer(W, sizeof(T))
        T2 = T <: Signed ? _T2 : unsigned(_T2)
    else
        T2 = T
    end
    t = Expr(:tuple)
    foreach(w -> push!(t.args, Expr(:call, :(Core.VecElement), T2(F*w + O))), 0:W-1)
    Expr(:block, Expr(:meta, :inline), Expr(:call, :Vec, t))
end

"""
  vrange(::Val{W}, i::I, ::Val{O}, ::Val{F})

W - Vector width
i::I - dynamic offset
O - static offset
F - static multiplicative factor
"""
@generated function vrangeincr(::Val{W}, i::I, ::Val{O}, ::Val{F}) where {W,I<:Integer,O,F}
    isone(W) && return Expr(:block, Expr(:meta,:inline), :(vadd(i, $(O % I))))
    bytes = pick_integer_bytes(W, sizeof(I))
    bits = 8bytes
    jtypesym = Symbol(:Int, bits)
    iexpr = bytes == sizeof(I) ? :i : Expr(:call, :%, :i, jtypesym)
    typ = "i$(bits)"
    vtyp = vtype(W, typ)
    rangevec = join(("$typ $(F*w + O)" for w ∈ 0:W-1), ", ")
    instrs = """
        %ie = insertelement $vtyp undef, $typ %0, i32 0
        %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
        %res = add nsw $vtyp %v, <$rangevec>
        ret $vtyp %res
    """
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($instrs, _Vec{$W,$jtypesym}, Tuple{$jtypesym}, $iexpr))
    end
end
@generated function vrangeincr(::Val{W}, i::T, ::Val{O}, ::Val{F}) where {W,T<:FloatingTypes,O,F}
    isone(W) && return Expr(:block, Expr(:meta,:inline), :(Base.FastMath.add_fast(i, $(T(O)))))
    typ = LLVM_TYPES[T]
    vtyp = vtype(W, typ)
    rangevec = join(("$typ $(F*w+O).0" for w ∈ 0:W-1), ", ")
    instrs = """
        %ie = insertelement $vtyp undef, $typ %0, i32 0
        %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
        %res = fadd fast $vtyp %v, <$rangevec>
        ret $vtyp %res
    """
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($instrs, _Vec{$W,$T}, Tuple{$T}, i))
    end
end
@generated function vrangemul(::Val{W}, i::I, ::Val{O}, ::Val{F}) where {W,I<:Integer,O,F}
    isone(W) && return Expr(:block, Expr(:meta,:inline), :(vmul(i, $(O % I))))
    bytes = pick_integer_bytes(W, sizeof(T))
    bits = 8bytes
    jtypesym = Symbol(:Int, bits)
    iexpr = bytes == sizeof(I) ? :i : Expr(:call, :%, :i, jtypesym)
    typ = "i$(bits)"
    vtyp = vtype(W, typ)
    rangevec = join(("$typ $(F*w+O)" for w ∈ 0:W-1), ", ")
    instrs = """
        %ie = insertelement $vtyp undef, $typ %0, i32 0
        %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
        %res = mul nsw $vtyp %v, <$rangevec>
        ret $vtyp %res
    """
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall(instrs, _Vec{$W,$jtypesym}, Tuple{$jtypesym}, $iexpr))
    end
end
@generated function vrangemul(::Val{W}, i::T, ::Val{O}, ::Val{F}) where {W,T<:FloatingTypes,O,F}
    isone(W) && return Expr(:block, Expr(:meta,:inline), :(Base.FastMath.mul_fast(i, $(T(O)))))
    typ = LLVM_TYPES[T]
    vtyp = vtype(W, typ)
    rangevec = join(("$typ $(F*w+O).0" for w ∈ 0:W-1), ", ")
    instrs = """
        %ie = insertelement $vtyp undef, $typ %0, i32 0
        %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
        %res = fmul fast $vtyp %v, <$rangevec>
        ret $vtyp %res
    """
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall(instrs, _Vec{$W,$T}, Tuple{$T}, i))
    end
end

@inline Vec(i::MM{W,X}) where {W,X} = vrangeincr(Val{W}(), data(i), Val{0}(), Val{X}())
@inline Vec(i::MM{W,X,StaticInt{N}}) where {W,X,N} = vrange(Val{W}(), Int, Val{N}(), Val{X}())
@inline Vec(i::MM{1}) = data(i)
@inline Vec(i::MM{1,<:Any,StaticInt{N}}) where {N} = N
@inline Base.convert(::Type{Vec{W,T}}, i::MM{W,X}) where {W,X,T} = vrangeincr(Val{W}(), T(data(i)), Val{0}(), Val{X}())
@inline Base.convert(::Type{T}, i::MM{W,X}) where {W,X,T<:NativeTypes} = vrangeincr(Val{W}(), T(data(i)), Val{0}(), Val{X}())

# Addition
# @inline Base.:(+)(i::MM{W}, j::MM{W}) where {W} = vadd(vrange(i), vrange(j))
# @inline Base.:(+)(i::MM{W,X}, j::MM{W,Y}) where {W} = vrange(vrangeincr(Val{W}(), data(i) + data(j), Val{0}(), StaticInt{X}() + StaticInt{Y}()))
@inline Base.:(+)(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} = MM{W}(vadd(data(i), data(j)), StaticInt{X}() + StaticInt{Y}())
@inline Base.:(+)(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vadd(Vec(i), j)
@inline Base.:(+)(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vadd(i, Vec(j))

# @inline vadd(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vadd(i.i, j))
# @inline vadd(j::Integer, i::MM{W,X}) where {W,X} = MM{W,X}(vadd(i.i, j))
@inline vadd(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} = MM{W}(vadd(data(i), data(j)), StaticInt{X}() + StaticInt{Y}())
@inline vadd(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vadd(Vec(i), j)
@inline vadd(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vadd(i, Vec(j))
# Subtraction
@inline Base.:(-)(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} = MM{W}(vsub(data(i), data(j)), StaticInt{X}() - StaticInt{Y}())
# @inline Base.:(-)(i::MM{W,X}, j::MM{W,X}) where {W,X} = vsub(data(i), data(j))
@inline Base.:(-)(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vsub(Vec(i), j)
@inline Base.:(-)(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vsub(i, Vec(j))
@inline vsub(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} = MM{W}(vsub(data(i), data(j)), StaticInt{X}() + StaticInt{Y}())
@inline vsub(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vsub(Vec(i), j)
@inline vsub(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vsub(i, Vec(j))
# Multiplication
@inline Base.:(*)(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vmul(Vec(i), j)
@inline Base.:(*)(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vmul(i, Vec(j))
@inline Base.:(*)(i::MM{W}, j::MM{W}) where {W} = vmul(Vec(i), Vec(j))
@inline vmul(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vmul(Vec(i), j)
@inline vmul(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vmul(i, Vec(j))
@inline vmul(i::MM{W}, j::MM{W}) where {W} = vmul(Vec(i), Vec(j))
@inline vmul(i::MM, j::Integer) = vmul(Vec(i), j)
@inline vmul(j::Integer, i::MM) = vmul(j, Vec(i))

# Multiplication without promotion
@inline vmul_no_promote(a, b) = vmul(a, b)
@inline vmul_no_promote(a::MM{W}, b) where {W} = MM{W}(vmul(a.i, b))
@inline vmul_no_promote(a, b::MM{W}) where {W} = MM{W}(vmul(a, b.i))
@inline vmul_no_promote(a::MM{W}, b::MM{W}) where {W} = vmul(a, b) # must promote
vmul_no_promote(a::MM, b::MM) = throw("Dimension mismatch.")

# Division
@generated function floattype(::Val{W}) where {W}
    (REGISTER_SIZE ÷ W) ≥ 8 ? :Float64 : :Float32
end
@inline Base.float(i::MM{W,X}) where {W,X} = Vec(MM{W,X}(floattype(Val{W}())(i.i)))
@inline Base.:(/)(i::MM, j::T) where {T<:Real} = float(i) / j
@inline Base.:(/)(j::T, i::MM) where {T<:Real} = j / float(i)
@inline Base.:(/)(i::MM, j::MM) = float(i) / float(j)
@inline Base.inv(i::MM) = inv(float(i))
@inline Base.:(/)(vu::VecUnroll, m::MM) = vu * inv(m)
@inline Base.:(/)(m::MM, vu::VecUnroll) = Vec(m) / vu

# @inline Base.:(<<)(i::MM, j::IntegerTypes) = Vec(i) << j
# @inline Base.:(>>)(i::MM, j::IntegerTypes) = Vec(i) >> j
# @inline Base.:(>>>)(i::MM, j::IntegerTypes) = Vec(i) >>> j
# @inline Base.:(<<)(i::MM, j::Vec) = Vec(i) << j
# @inline Base.:(>>)(i::MM, j::Vec) = Vec(i) >> j
# @inline Base.:(>>>)(i::MM, j::Vec) = Vec(i) >>> j

@inline Base.:(<<)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:IntegerTypes} = MM{W}(i.i << j, StaticInt{X}() << j)
@inline Base.:(>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:IntegerTypes} = MM{W}(i.i >> j, StaticInt{X}() >> j)
@inline Base.:(>>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:IntegerTypes} = MM{W}(i.i >>> j, StaticInt{X}() >>> j)

@inline Base.:(<<)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:StaticInt} = MM{W}(i.i << j, StaticInt{X}() << j)
@inline Base.:(>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:StaticInt} = MM{W}(i.i >> j, StaticInt{X}() >> j)
@inline Base.:(>>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:StaticInt} = MM{W}(i.i >>> j, StaticInt{X}() >>> j)

for (f,op) ∈ [
    (:scalar_less, :(<)), (:scalar_greater,:(>)), (:scalar_greaterequal,:(≥)), (:scalar_lessequal,:(≤)), (:scalar_equal,:(==)), (:scalar_notequal,:(!=))
]
    @eval @inline $f(i::MM, j::Real) = $op(data(i), j)
    @eval @inline $f(i::Real, j::MM) = $op(i, data(j))
    @eval @inline $f(i::MM, ::StaticInt{j}) where {j} = $op(data(i), j)
    @eval @inline $f(::StaticInt{i}, j::MM) where {i} = $op(i, data(j))
    @eval @inline $f(i::MM, j::MM) = $op(data(i), data(j))
    @eval @inline $f(i, j) = $op(i, j)
end

for f ∈ [:(<<), :(>>>), :(>>)]
    @eval begin
        @inline Base.$f(i::MM{W,X,T}, v::SignedHW) where {W,X,T<:SignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T}, v::SignedHW) where {W,X,T<:UnsignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T}, v::UnsignedHW) where {W,X,T<:SignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T}, v::UnsignedHW) where {W,X,T<:UnsignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T}, v::IntegerTypes) where {W,X,T<:StaticInt} = $f(Vec(i), v)

        @inline Base.$f(v::SignedHW, i::MM{W,X,T}) where {W,X,T<:SignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::UnsignedHW, i::MM{W,X,T}) where {W,X,T<:SignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::SignedHW, i::MM{W,X,T}) where {W,X,T<:UnsignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::UnsignedHW, i::MM{W,X,T}) where {W,X,T<:UnsignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::IntegerTypes, i::MM{W,X,T}) where {W,X,T<:StaticInt} = $f(v, Vec(i))

        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:SignedHW,T2<:SignedHW} = $f(Vec(i), Vec(j))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:UnsignedHW,T2<:SignedHW} = $f(Vec(i), Vec(j))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:SignedHW,T2<:UnsignedHW} = $f(Vec(i), Vec(j))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:UnsignedHW,T2<:UnsignedHW} = $f(Vec(i), Vec(j))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:StaticInt,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:IntegerTypes,T2<:StaticInt} = $f(Vec(i), Vec(j))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:StaticInt,T2<:StaticInt} = $f(Vec(i), Vec(j))

        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:SignedHW,T2<:SignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:UnsignedHW,T2<:SignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:SignedHW,T2<:UnsignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:UnsignedHW,T2<:UnsignedHW} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:StaticInt,T2<:IntegerTypes} = $f(Vec(i), v)

        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:SignedHW,T2<:SignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:UnsignedHW,T2<:SignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:SignedHW,T2<:UnsignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:UnsignedHW,T2<:UnsignedHW} = $f(v, Vec(i))
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:IntegerTypes,T2<:StaticInt} = $f(v, Vec(i))
    end
end
for f ∈ [ :(÷), :(&), :(|), :(⊻), :(%), :(<), :(>), :(≥), :(≤), :(==), :(!=), :min, :max, :copysign]
    @eval begin
        @inline Base.$f(i::MM{W,X,T}, v::IntegerTypes) where {W,X,T<:IntegerTypes} = $f(Vec(i), v)
        @inline Base.$f(v::IntegerTypes, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:IntegerTypes,T2<:IntegerTypes} = $f(Vec(i), v)
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:IntegerTypes,T2<:IntegerTypes} = $f(v, Vec(i))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:IntegerTypes,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
    end
end
for f ∈ [:(<), :(>), :(≥), :(≤), :(==), :(!=), :min, :max, :copysign]

    @eval begin
        # left floating
        @inline Base.$f(i::MM{W,X,T}, v::IntegerTypes) where {W,X,T<:FloatingTypes} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:FloatingTypes,T2<:IntegerTypes} = $f(Vec(i), v)
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:FloatingTypes,T2<:IntegerTypes} = $f(v, Vec(i))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:FloatingTypes,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
        # right floating
        @inline Base.$f(i::MM{W,X,T}, v::FloatingTypes) where {W,X,T<:IntegerTypes} = $f(Vec(i), v)
        @inline Base.$f(v::IntegerTypes, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:IntegerTypes,T2<:FloatingTypes} = $f(Vec(i), v)
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:IntegerTypes,T2<:FloatingTypes} = $f(v, Vec(i))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:IntegerTypes,T2<:FloatingTypes} = $f(Vec(i), Vec(j))
        # both floating
        @inline Base.$f(i::MM{W,X,T}, v::FloatingTypes) where {W,X,T<:FloatingTypes} = $f(Vec(i), v)
        @inline Base.$f(i::MM{W,X,T1}, v::AbstractSIMDVector{W,T2}) where {W,X,T1<:FloatingTypes,T2<:FloatingTypes} = $f(Vec(i), v)
        @inline Base.$f(v::AbstractSIMDVector{W,T1}, i::MM{W,X,T2}) where {W,X,T1<:FloatingTypes,T2<:FloatingTypes} = $f(v, Vec(i))
        @inline Base.$f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:FloatingTypes,T2<:FloatingTypes} = $f(Vec(i), Vec(j))
    end
    if f === :copysign
        @eval begin
            @inline Base.$f(v::Float32, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
            @inline Base.$f(v::Float32, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
            @inline Base.$f(v::Float64, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
            @inline Base.$f(v::Float64, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
        end
    else
        @eval begin
            @inline Base.$f(v::FloatingTypes, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
            @inline Base.$f(v::FloatingTypes, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
        end
    end
end

# @inline vadd(::MM{W,Zero}, v::AbstractSIMDVector{W,T}) where {W,T} = vadd(vrange(Val{W}(), T, Val{0}(), Val{1}()), v)
# @inline vadd(v::AbstractSIMDVector{W,T}, ::MM{W,Zero}) where {W,T} = vadd(vrange(Val{W}(), T, Val{0}(), Val{1}()), v)
@inline vadd(i::MM{W,Zero}, j::MM{W,Zero}) where {W} = vrange(Val{W}(), Int, Val{0}(), Val{2}())
# @inline vadd(a::MM, ::Zero) = a
# @inline vadd(::Zero, a::MM) = a
# @inline Base.:(+)(a::MM, ::Zero) = a
# @inline Base.:(+)(::Zero, a::MM) = a
# # @inline vmul(::MM{W,Zero}, i) where {W} = svrangemul(Val{W}(), i, Val{0}())
# @inline vmul(i, ::MM{W,Zero}) where {W} = svrangemul(Val{W}(), i, Val{0}())


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1APIGroup <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    preferredVersion::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery } # spec name: preferredVersion
    serverAddressByClientCIDRs::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1ServerAddressByClientCIDR} } # spec name: serverAddressByClientCIDRs
    versions::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery} } # spec name: versions

    function IoK8sApimachineryPkgApisMetaV1APIGroup(;apiVersion=nothing, kind=nothing, name=nothing, preferredVersion=nothing, serverAddressByClientCIDRs=nothing, versions=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroup, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroup, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroup, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroup, Symbol("preferredVersion"), preferredVersion)
        setfield!(o, Symbol("preferredVersion"), preferredVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroup, Symbol("serverAddressByClientCIDRs"), serverAddressByClientCIDRs)
        setfield!(o, Symbol("serverAddressByClientCIDRs"), serverAddressByClientCIDRs)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroup, Symbol("versions"), versions)
        setfield!(o, Symbol("versions"), versions)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1APIGroup

const _property_map_IoK8sApimachineryPkgApisMetaV1APIGroup = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("preferredVersion")=>Symbol("preferredVersion"), Symbol("serverAddressByClientCIDRs")=>Symbol("serverAddressByClientCIDRs"), Symbol("versions")=>Symbol("versions"))
const _property_types_IoK8sApimachineryPkgApisMetaV1APIGroup = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("preferredVersion")=>"IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery", Symbol("serverAddressByClientCIDRs")=>"Vector{IoK8sApimachineryPkgApisMetaV1ServerAddressByClientCIDR}", Symbol("versions")=>"Vector{IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery}")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroup }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1APIGroup))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroup }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1APIGroup[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroup }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1APIGroup[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1APIGroup)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("versions")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroup }, name::Symbol, val)
end


text:
#springboard file that contains gnum implementation parts.

include("gnum-constants.jl")
include("gnum-utilities.jl")

include("gnum-addition.jl")
include("gnum-subtraction.jl")
include("gnum-multiplication.jl")
include("gnum-division.jl")

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Blink
using Test

# IMPORTANT: Window(...) cannot appear inside of a @testset for as-of-yet
# unknown reasons.
w = Window(Blink.@d(:show => false, :width=>150, :height=>100), async=false);
@testset "size Tests" begin
    @test size(w) == [150,100]

    size(w, 200,200)
    @test size(w) == [200,200]
end

# @testset "async" begin
#     # Test that async Window() creation is faster than synchronous creation.
#     # (Repeat the test a few times, just to be sure it's consistent.)
#     for _ in 1:5
#         (@timed Window(Blink.@d(:show => false), async=true))[2] <
#          (@timed Window(Blink.@d(:show => false), async=false))[2]
#     end
# # end