text:
#=
Macros for the interface.
Many of these will eventually be removed.
Try to use the regular function interface in finch_interface.jl
=#
export @generateFor, @domain, @mesh, @solver, @stepper, @setSteps, @functionSpace, @trialFunction, @matrixFree,
        @testFunction, @nodes, @order, @boundary, @addBoundaryID, @referencePoint, @variable, @coefficient,
        @parameter, @testSymbol, @initial, @preStepFunction, @postStepFunction, @callbackFunction,
        @timeInterval, @weakForm, @fluxAndSource, @LHS, @RHS, @exportCode, @importCode,
        @customOperator, @customOperatorFile,
        @outputMesh, @useLog, @finalize

###############################################################################
# These will be kept.

#= Assume something like this
myf = @callbackFunction(
    function f(a,b,c)
        body...
    end
)
=#
macro callbackFunction(f)
    fname = string(f.args[1].args[1]);
    fargs = [string(f.args[1].args[i]) for i=2:length(f.args[1].args)];
    fbody = string(f.args[2]);
    return esc(:(callbackFunction($f, name=$fname, args=$fargs, body=$fbody)));
end

macro preStepFunction(f)
    return esc(quote
        function PRE_STEP_FUNCTION()
            $f
        end
        preStepFunction(PRE_STEP_FUNCTION);
    end)
end

macro postStepFunction(f)
    return esc(quote
        function POST_STEP_FUNCTION()
            $f
        end
        postStepFunction(POST_STEP_FUNCTION);
    end)
end

###############################################################################
# These three could be kept for a minor convenience, but may be removed eventually.

macro variable(var) return esc(:(@variable($var, SCALAR, NODAL))); end
macro variable(var, type) return esc(:(@variable($var, $type, NODAL))); end
macro variable(var, type, location)
    varsym = string(var);
    return esc(quote
        $var = variable($varsym, $type, $location);
    end)
end

macro coefficient(c, val) return esc(:(@coefficient($c, SCALAR, $val))); end
macro coefficient(c, type, val)
    csym = string(c);
    return esc(quote
        coefficient($csym, $val, $type);
    end)
end

macro parameter(p,val) return esc(:(@parameter($p, SCALAR, $val))); end
macro parameter(p, type, val)
    psym = string(p);
    return esc(quote
        parameter($psym, $val, $type);
    end)
end

#################################################################################
# These below will eventually be removed

# Initializes code generation module
macro generateFor(lang) return esc(:(@generateFor($lang, Finch.project_name, ""))); end
macro generateFor(lang, filename) return esc(:(@generateFor($lang, $filename, ""))); end
macro generateFor(lang, filename, header)
    return esc(quote
        generateFor($lang, filename=$filename, header=$header);
    end)
end

# Optionally create a log file
macro useLog() return :(@uselog(Finch.project_name)) end
macro useLog(name) return esc(:(@useLog($name, Finch.output_dir))) end
macro useLog(name, dir)
    return esc(quote
        useLog($name, dir=$dir);
    end)
end

# Sets dimension, domain shape, discretization type
macro domain(dims) return :(@domain($dims, SQUARE, UNIFORM_GRID)) end
macro domain(dims, geometry, mesh)
    return esc(quote
        domain($dims, shape=$geometry, grid=$mesh);
    end)
end

# Set the solver type: CG, etc.
macro solver(type)
    return esc(quote
        solverType($type);
    end)
end

# Set basis function type and order
macro functionSpace(space) return esc(:(@functionSpace($space, 1, 1))) end
macro functionSpace(space, N) return esc(:(@functionSpace($space, $N, $N))) end
macro functionSpace(space, Nmin, Nmax)
    return esc(quote
        functionSpace(space=$space, orderMin=$Nmin, orderMax=$Nmax);
    end)
end

macro trialFunction(space, N) return esc(:(@trialFunction($space, $N, $N))) end
macro trialFunction(space, Nmin, Nmax)
    return esc(quote
        trialSpace(space=$space, orderMin=$Nmin, orderMax=Nmax);
    end)
end

macro testFunction(space, N) return esc(:(@testFunction($space, $N, $N))) end
macro testFunction(space, Nmin, Nmax)
    return esc(quote
        testSpace(space=$space, orderMin=$Nmin, orderMax=Nmax);
    end)
end

# Sets elemental node locations: Lobatto, etc.
macro nodes(nodetype)
    return esc(quote
        nodeType($nodetype);
    end)
end

# Sets time stepper type
macro stepper(s) return esc(:(@stepper($s, 0))) end
macro stepper(s, cfl)
    return esc(quote
        timeStepper($s, cfl=$cfl);
    end)
end

# Sets time stepper dt and Nsteps to specified values
macro setSteps(dt, steps)
    return esc(quote
        setSteps($dt, $steps);
    end)
end

# Selects matrix free
macro matrixFree() return esc(:(@matrixFree(1000, 1e-6))) end
macro matrixFree(max, tol)
    return esc(quote
        matrixFree(maxiters=$max, tol=$tol);
    end)
end

# Adds a custom operator
macro customOperator(s, handle)
    symb = string(s);
    return esc(quote
        customOperator($symb, $handle);
    end)
end

# Adds a set of custom operators defined in a file
macro customOperatorFile(file)
    return esc(quote
        customOperatorFile($file);
    end)
end

# ============= end config , begin prob ==============

# This one-argument version should only be used for importing from a file
macro mesh(m)
    return esc(quote
        mesh($m);
    end)
end
# More than one argument is for generating a simple mesh.
macro mesh(m,N) return esc(quote @mesh($m,$N,1,[0,1]); end) end
macro mesh(m,N,bids) return esc(quote @mesh($m,$N,$bids,[0,1]); end) end
macro mesh(m, N, bids, interval)
    return esc(quote
        mesh($m, elsperdim=$N, bids=$bids, interval=$interval);
    end)
end

# Write the mesh to a MSH file
macro outputMesh(file) return esc(:(@outputMesh($file, MSH_V2))); end
macro outputMesh(file, format)
    return esc(quote
        exportMesh($file, $format);
    end)
end

macro addBoundaryID(bid, expression)
    return esc(quote
        trueOnBdry(x, y=0, z=0) = $expression; # points with x,y,z on this bdry segment evaluate true here
        addBoundaryID($bid, trueOnBdry);
    end)
end

macro boundary(var, bid, bc_type) return esc(:(@boundary($var,$bid,$bc_type, 0))); end
macro boundary(var, bid, bc_type, bc_exp)
    return esc(quote
        boundary($var, $bid, $bc_type, $bc_exp);
    end)
end

macro referencePoint(var, pos, val)
    return esc(quote
        referencePoint($var, $pos, $val);
    end)
end

macro timeInterval(t)
    return esc(quote
        timeInterval($t);
    end)
end

macro initial(var, ic)
    return esc(quote
        initial($var, $ic)
    end)
end

macro testSymbol(var) return esc(:(@testSymbol($var, SCALAR))); end
macro testSymbol(var, type)
    varsym = string(var);
    return esc(quote
        testSymbol($varsym, $type);
    end)
end

macro weakForm(var, ex)
    return esc(quote
        weakForm($var, $ex);
    end)
end

macro fluxAndSource(var, fex, sex)
    return esc(quote
        fluxAndSource($var, $fex, $sex);
    end)
end

# Import and export the code layer functions to allow manual changes
macro exportCode(file)
    return esc(quote
        exportCode($file);
    end)
end

macro importCode(file)
    return esc(quote
        importCode($file);
    end)
end

macro finalize()
    return esc(:(Finch.finalize()));
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


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

@doc raw"""Represents a Rados Block Device mount that lasts the lifetime of a pod. RBD volumes support ownership management and SELinux relabeling.

    IoK8sApiCoreV1RBDVolumeSource(;
        fsType=nothing,
        image=nothing,
        keyring=nothing,
        monitors=nothing,
        pool=nothing,
        readOnly=nothing,
        secretRef=nothing,
        user=nothing,
    )

    - fsType::String : Filesystem type of the volume that you want to mount. Tip: Ensure that the filesystem type is supported by the host operating system. Examples: \"ext4\", \"xfs\", \"ntfs\". Implicitly inferred to be \"ext4\" if unspecified. More info: https://kubernetes.io/docs/concepts/storage/volumes#rbd
    - image::String : The rados image name. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
    - keyring::String : Keyring is the path to key ring for RBDUser. Default is /etc/ceph/keyring. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
    - monitors::Vector{String} : A collection of Ceph monitors. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
    - pool::String : The rados pool name. Default is rbd. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
    - readOnly::Bool : ReadOnly here will force the ReadOnly setting in VolumeMounts. Defaults to false. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
    - secretRef::IoK8sApiCoreV1LocalObjectReference : SecretRef is name of the authentication secret for RBDUser. If provided overrides keyring. Default is nil. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
    - user::String : The rados user name. Default is admin. More info: https://examples.k8s.io/volumes/rbd/README.md#how-to-use-it
"""
mutable struct IoK8sApiCoreV1RBDVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    image::Any # spec type: Union{ Nothing, String } # spec name: image
    keyring::Any # spec type: Union{ Nothing, String } # spec name: keyring
    monitors::Any # spec type: Union{ Nothing, Vector{String} } # spec name: monitors
    pool::Any # spec type: Union{ Nothing, String } # spec name: pool
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1LocalObjectReference } # spec name: secretRef
    user::Any # spec type: Union{ Nothing, String } # spec name: user

    function IoK8sApiCoreV1RBDVolumeSource(;fsType=nothing, image=nothing, keyring=nothing, monitors=nothing, pool=nothing, readOnly=nothing, secretRef=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("image"), image)
        setfield!(o, Symbol("image"), image)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("keyring"), keyring)
        setfield!(o, Symbol("keyring"), keyring)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("monitors"), monitors)
        setfield!(o, Symbol("monitors"), monitors)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("pool"), pool)
        setfield!(o, Symbol("pool"), pool)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1RBDVolumeSource, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiCoreV1RBDVolumeSource

const _property_map_IoK8sApiCoreV1RBDVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("image")=>Symbol("image"), Symbol("keyring")=>Symbol("keyring"), Symbol("monitors")=>Symbol("monitors"), Symbol("pool")=>Symbol("pool"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiCoreV1RBDVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("image")=>"String", Symbol("keyring")=>"String", Symbol("monitors")=>"Vector{String}", Symbol("pool")=>"String", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1LocalObjectReference", Symbol("user")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1RBDVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1RBDVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1RBDVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1RBDVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1RBDVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1RBDVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1RBDVolumeSource)
    (getproperty(o, Symbol("image")) === nothing) && (return false)
    (getproperty(o, Symbol("monitors")) === nothing) && (return false)
    true
end

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

ext: jl


text:
function elimina_zeros(v)
    p = filter!(x->x!=0,v);
    return p;
end


text:
#= modelcode = "quantum_ising_2" so on...
SpinInfo(modelcode, factorization, externalfield, temperature)
isZeroTemperature(spininfo)
isClassical(spininfo)
getModelname(spininfo)
getStates(spininfo)
isSymmetricFactorization(spininfo)
getTemperature(spininfo)
getExternalfield(spininfo)
getEnvParameters(spininfo)
setTemperature!(spininfo, temperature)
setExternalfield!(spininfo, externalfield)
setEnvParameters!(spininfo, temperature, externalfield)
testSpinInfo()
=#

"""
	SpinInfo{T}(modelcode::AbstractString, factorization::AbstractString, externalfield::T, temperature::T)

Store informations about the spin model.

# Arguments

* `modelcode::AbstractString`: `[quantum,classical]_[modelname]_[numberOfState]`
* `factorization::AbstractString`: `sym` or `asym` (default: `asym`)
* `temperature::T`: (default: one(T))
* `externalfield::T`: (default: one(T))

# Examples

```julia
modelcode = "quantum_ising"
factorization = "sym"
externalfield = 2.0
temperature = 1.5
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode, factorization, externalfield, temperature)
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode, factorization, externalfield)
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode, externalfield, temperature)
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode, externalfield)
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode, factorization)
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode)
```

```jldoctest spininfo
modelcode = "quantum_ising"
factorization = "sym"
externalfield = 2.0
temperature = 1.5
spininfo = HoTRG_lemon.SpinModule.SpinInfo(modelcode, factorization, externalfield, temperature)

# output

HoTRG_lemon.SpinModule.SpinInfo{Float64}("quantum","ising",2,"sym",1.5,2.0)
```

"""
type SpinInfo{T}
	quantumOrClassical::AbstractString
	modelname::AbstractString
	numberOfState::Int
	factorization::AbstractString

	temperature::T
	externalfield::T

	function SpinInfo{T}(modelcode::AbstractString, factorization::AbstractString, externalfield::T, temperature::T)
		this = new{T}()
		this.quantumOrClassical, this.modelname, this.numberOfState = getInfoFromModelcode(modelcode)
		this.factorization = checkFactorization(factorization)
		this.externalfield = externalfield
		this.temperature = temperature
		return this
	end

end

SpinInfo{T}(modelcode::AbstractString, factorization::AbstractString, externalfield::T, temperature::T) =
	SpinInfo{T}(modelcode, factorization, externalfield, temperature)
SpinInfo{T}(modelcode::AbstractString, factorization::AbstractString, externalfield::T) =
	SpinInfo{T}(modelcode, factorization, externalfield, one(T))
SpinInfo{T}(modelcode::AbstractString, externalfield::T, temperature::T) =
	SpinInfo{T}(modelcode, "asym", externalfield, temperature)
SpinInfo{T}(modelcode::AbstractString, externalfield::T) =
	SpinInfo{T}(modelcode, "asym", externalfield, one(T))
SpinInfo(modelcode::AbstractString, factorization::AbstractString) =
	SpinInfo{Float64}(modelcode, factorization, one(Float64), one(Float64))
SpinInfo(modelcode::AbstractString) =
	SpinInfo{Float64}(modelcode, "asym", one(Float64), one(Float64))

#---
# checking for initialize SpinInfo type

function getInfoFromModelcode(modelcode::AbstractString)
	model = split(modelcode,"_")
	quantumOrClassical = checkQuantumOrClassical(model[1])
	modelname = lowercase(model[2])
	numberOfState = try
			parse(Int, model[3])
		catch
			2
		end
	return quantumOrClassical, modelname, numberOfState
end

function checkFactorization(factorization::AbstractString)
	factorization = lowercase(factorization)
	if !(factorization in ("sym", "asym"))
		println("warning: factorization should be either 'sym' or 'asym'. now factorization set to be 'asym'.")
		factorization = "asym"
	end
	return factorization
end

function checkQuantumOrClassical(QOC::AbstractString)
	QOC = lowercase(QOC)
	if !(QOC in ("quantum","classical"))
		error("should be quantum or classical")
	end
	return QOC
end

#---
# functions for SpinInfo type

"""
	getModelname(spininfo)

Return the model name.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.getModelname(spininfo)
"ising"
```
"""
function getModelname(spininfo::SpinInfo)
	return spininfo.modelname
end

"""
	isClassical(spininfo)

Return true when the spin is classical, false for quantum.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.isClassical(spininfo)
false
```
"""
function isClassical(spininfo::SpinInfo)
	if spininfo.quantumOrClassical == "classical"
		return true
	else
		return false
	end
end

"""
	getStates(spininfo)

Return the number of states.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.getStates(spininfo)
2
```
"""
function getStates(spininfo::SpinInfo)
	return spininfo.numberOfState
end

"""
	isSymmetricFactorization(spininfo)

Return `true` for symmetric factorization

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.isSymmetricFactorization(spininfo)
true
```
"""
function isSymmetricFactorization(spininfo::SpinInfo)
	if spininfo.factorization == "sym"
		return true
	else
		return false
	end
end

"""
	getTemperature(spininfo)

Return the temperature of the spin system.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.getTemperature(spininfo)
1.5
```
"""
function getTemperature(spininfo::SpinInfo)
	return spininfo.temperature
end

function isZeroTemperatuer{T}(temperature::T)
	if temperature == zero(T)
		return true
	else
		return false
	end
end

"""
	isZeroTemperature(spininfo::SpinInfo)

Return true when  the system is at the Zero Temperature.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.isZeroTemperature(spininfo)
false
```
"""
function isZeroTemperature{T}(spininfo::SpinInfo{T})
	return isZeroTemperatuer(getTemperature(spininfo))
end

"""
	getExternalfield(spininfo)

Return external field.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.getExternalfield(spininfo)
2.0
```
"""
function getExternalfield(spininfo::SpinInfo)
	return spininfo.externalfield
end

"""
	getEnvParameters(spininfo)

Return temperature and external field.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.getEnvParameters(spininfo)
(1.5,2.0)
```
"""
function getEnvParameters(spininfo::SpinInfo)
	return spininfo.temperature, spininfo.externalfield
end

"""
	setTemperature!(spininfo, temperature)

Set temperature for spininfo to the input temperature.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.setTemperature!(spininfo, 3.0);

julia> HoTRG_lemon.SpinModule.getTemperature(spininfo)
3.0
```
"""
function setTemperature!{T}(spininfo::SpinInfo{T}, temperature::T)
	spininfo.temperature = temperature
end

"""
	setExternalfield!(spininfo, externalfield)

Set spininfo's external field  into the input.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.setExternalfield!(spininfo, 4.0);

julia> HoTRG_lemon.SpinModule.getExternalfield(spininfo)
4.0
```
"""
function setExternalfield!{T}(spininfo::SpinInfo{T}, externalfield::T)
	spininfo.externalfield = externalfield
end

"""
	setEnvParameters!(spininfo, temperature, externalfield)

Set spininfo's temperature and external field.

```jldoctest spininfo
julia> HoTRG_lemon.SpinModule.setEnvParameters!(spininfo, 1.2, 3.4);

julia> HoTRG_lemon.SpinModule.getEnvParameters(spininfo)
(1.2,3.4)
```
"""
function setEnvParameters!{T}(spininfo::SpinInfo{T}, temperature::T, externalfield::T)
	spininfo.temperature = temperature
	spininfo.externalfield = externalfield
end

#---
# test to make a example SpinInfo instatce

"""
	testSpinInfo()

Test SpinInfo for quantum ising model.

`testSpininfo()`
"""
function testSpinInfo()
	modelcode = "quantum_ising"
	factorization = "sym"
	externalfield = 2.0
	temperature = 1.5
	spininfo = SpinInfo(modelcode, factorization, externalfield, temperature)
	return spininfo
end


text:
module OptionsModule

import Distributed: nworkers
import LossFunctions: L2DistLoss
import Zygote: gradient
#TODO - eventually move some of these
# into the SR call itself, rather than
# passing huge options at once.
import ..OperatorsModule:
    plus,
    pow,
    mult,
    sub,
    div,
    log_abs,
    log10_abs,
    log2_abs,
    log1p_abs,
    sqrt_abs,
    acosh_abs,
    atanh_clip
import ..EquationModule: Node, string_tree
import ..OptionsStructModule: Options, ComplexityMapping

"""
         build_constraints(una_constraints, bin_constraints,
                           unary_operators, binary_operators)

Build constraints on operator-level complexity from a user-passed dict.
"""
function build_constraints(
    una_constraints, bin_constraints, unary_operators, binary_operators, nuna, nbin
)::Tuple{Array{Int,1},Array{Tuple{Int,Int},1}}
    # Expect format ((*)=>(-1, 3)), etc.
    # TODO: Need to disable simplification if (*, -, +, /) are constrained?
    #  Or, just quit simplification is constraints violated.

    is_bin_constraints_already_done = typeof(bin_constraints) <: Array{Tuple{Int,Int},1}
    is_una_constraints_already_done = typeof(una_constraints) <: Array{Int,1}

    if typeof(bin_constraints) <: Array && !is_bin_constraints_already_done
        bin_constraints = Dict(bin_constraints)
    end
    if typeof(una_constraints) <: Array && !is_una_constraints_already_done
        una_constraints = Dict(una_constraints)
    end

    if una_constraints === nothing
        una_constraints = [-1 for i in 1:nuna]
    elseif !is_una_constraints_already_done
        una_constraints::Dict
        _una_constraints = Int[]
        for (i, op) in enumerate(unary_operators)
            did_user_declare_constraints = haskey(una_constraints, op)
            if did_user_declare_constraints
                constraint::Int = una_constraints[op]
                push!(_una_constraints, constraint)
            else
                push!(_una_constraints, -1)
            end
        end
        una_constraints = _una_constraints
    end
    if bin_constraints === nothing
        bin_constraints = [(-1, -1) for i in 1:nbin]
    elseif !is_bin_constraints_already_done
        bin_constraints::Dict
        _bin_constraints = Tuple{Int,Int}[]
        for (i, op) in enumerate(binary_operators)
            did_user_declare_constraints = haskey(bin_constraints, op)
            if did_user_declare_constraints
                constraint::Tuple{Int,Int} = bin_constraints[op]
                push!(_bin_constraints, constraint)
            else
                push!(_bin_constraints, (-1, -1))
            end
        end
        bin_constraints = _bin_constraints
    end

    return una_constraints, bin_constraints
end

function binopmap(op)
    if op == plus
        return +
    elseif op == mult
        return *
    elseif op == sub
        return -
    elseif op == div
        return /
    elseif op == ^
        return pow
    end
    return op
end

function unaopmap(op)
    if op == log
        return log_abs
    elseif op == log10
        return log10_abs
    elseif op == log2
        return log2_abs
    elseif op == log1p
        return log1p_abs
    elseif op == sqrt
        return sqrt_abs
    elseif op == acosh
        return acosh_abs
    elseif op == atanh
        return atanh_clip
    end
    return op
end

"""
    Options(;kws...)

Construct options for `EquationSearch` and other functions.
The current arguments have been tuned using the median values from
https://github.com/MilesCranmer/PySR/discussions/115.

# Arguments
- `binary_operators`: Tuple of binary
    operators to use. Each operator should be defined for two input scalars,
    and one output scalar. All operators need to be defined over the entire
    real line (excluding infinity - these are stopped before they are input).
    Thus, `log` should be replaced with `log_abs`, etc.
    For speed, define it so it takes two reals
    of the same type as input, and outputs the same type. For the SymbolicUtils
    simplification backend, you will need to define a generic method of the
    operator so it takes arbitrary types.
- `unary_operators`: Same, but for
    unary operators (one input scalar, gives an output scalar).
- `constraints`: Array of pairs specifying size constraints
    for each operator. The constraints for a binary operator should be a 2-tuple
    (e.g., `(-1, -1)`) and the constraints for a unary operator should be an `Int`.
    A size constraint is a limit to the size of the subtree
    in each argument of an operator. e.g., `[(^)=>(-1, 3)]` means that the
    `^` operator can have arbitrary size (`-1`) in its left argument,
    but a maximum size of `3` in its right argument. Default is
    no constraints.
- `batching`: Whether to evolve based on small mini-batches of data,
    rather than the entire dataset.
- `batchSize`: What batch size to use if using batching.
- `loss`: What loss function to use. Can be one of
    the following losses, or any other loss of type
    `SupervisedLoss`. You can also pass a function that takes
    a scalar target (left argument), and scalar predicted (right
    argument), and returns a scalar. This will be averaged
    over the predicted data. If weights are supplied, your
    function should take a third argument for the weight scalar.
    Included losses:
        Regression:
            - `LPDistLoss{P}()`,
            - `L1DistLoss()`,
            - `L2DistLoss()` (mean square),
            - `LogitDistLoss()`,
            - `HuberLoss(d)`,
            - `L1EpsilonInsLoss(ϵ)`,
            - `L2EpsilonInsLoss(ϵ)`,
            - `PeriodicLoss(c)`,
            - `QuantileLoss(τ)`,
        Classification:
            - `ZeroOneLoss()`,
            - `PerceptronLoss()`,
            - `L1HingeLoss()`,
            - `SmoothedL1HingeLoss(γ)`,
            - `ModifiedHuberLoss()`,
            - `L2MarginLoss()`,
            - `ExpLoss()`,
            - `SigmoidLoss()`,
            - `DWDMarginLoss(q)`.
- `npopulations`: How many populations of equations to use. By default
    this is set equal to the number of cores
- `npop`: How many equations in each population.
- `ncyclesperiteration`: How many generations to consider per iteration.
- `ns`: Number of equations in each subsample during regularized evolution.
- `topn`: Number of equations to return to the host process, and to
    consider for the hall of fame.
- `complexity_of_operators`: What complexity should be assigned to each operator,
    and the occurrence of a constant or variable. By default, this is 1
    for all operators. Can be a real number as well, in which case
    the complexity of an expression will be rounded to the nearest integer.
    Input this in the form of, e.g., [(^) => 3, sin => 2].
- `complexity_of_constants`: What complexity should be assigned to use of a constant.
    By default, this is 1.
- `complexity_of_variables`: What complexity should be assigned to each variable.
    By default, this is 1.
- `alpha`: The probability of accepting an equation mutation
    during regularized evolution is given by exp(-delta_loss/(alpha * T)),
    where T goes from 1 to 0. Thus, alpha=infinite is the same as no annealing.
- `maxsize`: Maximum size of equations during the search.
- `maxdepth`: Maximum depth of equations during the search, by default
    this is set equal to the maxsize.
- `parsimony`: A multiplicative factor for how much complexity is
    punished.
- `useFrequency`: Whether to use a parsimony that adapts to the
    relative proportion of equations at each complexity; this will
    ensure that there are a balanced number of equations considered
    for every complexity.
- `useFrequencyInTournament`: Whether to use the adaptive parsimony described
    above inside the score, rather than just at the mutation accept/reject stage.
- `fast_cycle`: Whether to thread over subsamples of equations during
    regularized evolution. Slightly improves performance, but is a different
    algorithm.
- `migration`: Whether to migrate equations between processes.
- `hofMigration`: Whether to migrate equations from the hall of fame
    to processes.
- `fractionReplaced`: What fraction of each population to replace with
    migrated equations at the end of each cycle.
- `fractionReplacedHof`: What fraction to replace with hall of fame
    equations at the end of each cycle.
- `shouldOptimizeConstants`: Whether to use NelderMead optimization
    to periodically optimize constants in equations.
- `optimizer_nrestarts`: How many different random starting positions to consider
    when using NelderMead optimization.
- `hofFile`: What file to store equations to, as a backup.
- `perturbationFactor`: When mutating a constant, either
    multiply or divide by (1+perturbationFactor)^(rand()+1).
- `probNegate`: Probability of negating a constant in the equation
    when mutating it.
- `mutationWeights`: Relative probabilities of the mutations, in the order: MutateConstant, MutateOperator, AddNode, InsertNode, DeleteNode, Simplify, Randomize, DoNothing.
- `annealing`: Whether to use simulated annealing.
- `warmupMaxsize`: Whether to slowly increase the max size from 5 up to
    `maxsize`. If nonzero, specifies how many cycles (populations*iterations)
    before increasing by 1.
- `verbosity`: Whether to print debugging statements or
    not.
- `bin_constraints`: See `constraints`. This is the same, but specified for binary
    operators only (for example, if you have an operator that is both a binary
    and unary operator).
- `una_constraints`: Likewise, for unary operators.
- `seed`: What random seed to use. `nothing` uses no seed.
- `progress`: Whether to use a progress bar output (`verbosity` will
    have no effect).
- `probPickFirst`: Expressions in subsample are chosen based on, for
    p=probPickFirst: p, p*(1-p), p*(1-p)^2, and so on.
- `earlyStopCondition`: Float - whether to stop early if the mean loss gets below this value.
    Function - a function taking (loss, complexity) as arguments and returning true or false.
- `timeout_in_seconds`: Float64 - the time in seconds after which to exit (as an alternative to the number of iterations).
- `max_evals`: Int (or Nothing) - the maximum number of evaluations of expressions to perform.
- `skip_mutation_failures`: Whether to simply skip over mutations that fail or are rejected, rather than to replace the mutated
    expression with the original expression and proceed normally.
- `enable_autodiff`: Whether to enable automatic differentiation functionality. This is turned off by default.
    If turned on, this will be turned off if one of the operators does not have well-defined gradients.
- `nested_constraints`: Specifies how many times a combination of operators can be nested. For example,
    `[sin => [cos => 0], cos => [cos => 2]]` specifies that `cos` may never appear within a `sin`,
    but `sin` can be nested with itself an unlimited number of times. The second term specifies that `cos`
    can be nested up to 2 times within a `cos`, so that `cos(cos(cos(x)))` is allowed (as well as any combination
    of `+` or `-` within it), but `cos(cos(cos(cos(x))))` is not allowed. When an operator is not specified,
    it is assumed that it can be nested an unlimited number of times. This requires that there is no operator
    which is used both in the unary operators and the binary operators (e.g., `-` could be both subtract, and negation).
    For binary operators, both arguments are treated the same way, and the max of each argument is constrained.
- `deterministic`: Use a global counter for the birth time, rather than calls to `time()`. This gives
    perfect resolution, and is therefore deterministic. However, it is not thread safe, and must be used
    in serial mode.
"""
function Options(;
    binary_operators::NTuple{nbin,Any}=(+, -, /, *),
    unary_operators::NTuple{nuna,Any}=(),
    constraints=nothing,
    loss=L2DistLoss(),
    ns=12, #1 sampled from every ns per mutation
    topn=12, #samples to return per population
    complexity_of_operators=nothing,
    complexity_of_constants::Union{Nothing,Real}=nothing,
    complexity_of_variables::Union{Nothing,Real}=nothing,
    parsimony=0.0032f0,
    alpha=0.100000f0,
    maxsize=20,
    maxdepth=nothing,
    fast_cycle=false,
    migration=true,
    hofMigration=true,
    fractionReplacedHof=0.035f0,
    shouldOptimizeConstants=true,
    hofFile=nothing,
    npopulations=15,
    perturbationFactor=0.076f0,
    annealing=false,
    batching=false,
    batchSize=50,
    mutationWeights=[0.048, 0.47, 0.79, 5.1, 1.7, 0.0020, 0.00023, 0.21],
    crossoverProbability=0.066f0,
    warmupMaxsizeBy=0.0f0,
    useFrequency=true,
    useFrequencyInTournament=true,
    npop=33,
    ncyclesperiteration=550,
    fractionReplaced=0.00036f0,
    verbosity=convert(Int, 1e9),
    probNegate=0.01f0,
    seed=nothing,
    bin_constraints=nothing,
    una_constraints=nothing,
    progress=true,
    terminal_width=nothing,
    warmupMaxsize=nothing,
    optimizer_algorithm="BFGS",
    optimizer_nrestarts=2,
    optimize_probability=0.14f0,
    optimizer_iterations=8,
    recorder=nothing,
    recorder_file="pysr_recorder.json",
    probPickFirst=0.86f0,
    earlyStopCondition::Union{Function,AbstractFloat,Nothing}=nothing,
    stateReturn::Bool=false,
    timeout_in_seconds=nothing,
    max_evals=nothing,
    skip_mutation_failures::Bool=true,
    enable_autodiff::Bool=false,
    nested_constraints=nothing,
    deterministic=false,
    # for thermo constraints
    penalties=nothing,
) where {nuna,nbin}
    if warmupMaxsize !== nothing
        error(
            "warmupMaxsize is deprecated. Please use warmupMaxsizeBy, and give the time at which the warmup will end as a fraction of the total search cycles.",
        )
    end

    if hofFile === nothing
        hofFile = "hall_of_fame.csv" #TODO - put in date/time string here
    end

    @assert maxsize > 3
    @assert warmupMaxsizeBy >= 0.0f0

    # Make sure nested_constraints contains functions within our operator set:
    if nested_constraints !== nothing
        # Check that intersection of binary operators and unary operators is empty:
        for op in binary_operators
            if op ∈ unary_operators
                error(
                    "Operator $(op) is both a binary and unary operator. " *
                    "You can't use nested constraints.",
                )
            end
        end

        # Convert to dict:
        if !(typeof(nested_constraints) <: Dict)
            # Convert to dict:
            nested_constraints = Dict(
                [cons[1] => Dict(cons[2]...) for cons in nested_constraints]...
            )
        end
        for (op, nested_constraint) in nested_constraints
            if !(op ∈ binary_operators || op ∈ unary_operators)
                error("Operator $(op) is not in the operator set.")
            end
            for (nested_op, max_nesting) in nested_constraint
                if !(nested_op ∈ binary_operators || nested_op ∈ unary_operators)
                    error("Operator $(nested_op) is not in the operator set.")
                end
                @assert nested_op ∈ binary_operators || nested_op ∈ unary_operators
                @assert max_nesting >= -1 && typeof(max_nesting) <: Int
            end
        end

        # Lastly, we clean it up into a dict of (degree,op_idx) => max_nesting.
        new_nested_constraints = []
        # Dict()
        for (op, nested_constraint) in nested_constraints
            (degree, idx) = if op ∈ binary_operators
                2, findfirst(isequal(op), binary_operators)
            else
                1, findfirst(isequal(op), unary_operators)
            end
            new_max_nesting_dict = []
            # Dict()
            for (nested_op, max_nesting) in nested_constraint
                (nested_degree, nested_idx) = if nested_op ∈ binary_operators
                    2, findfirst(isequal(nested_op), binary_operators)
                else
                    1, findfirst(isequal(nested_op), unary_operators)
                end
                # new_max_nesting_dict[(nested_degree, nested_idx)] = max_nesting
                push!(new_max_nesting_dict, (nested_degree, nested_idx, max_nesting))
            end
            # new_nested_constraints[(degree, idx)] = new_max_nesting_dict
            push!(new_nested_constraints, (degree, idx, new_max_nesting_dict))
        end
        nested_constraints = new_nested_constraints
    end

    if typeof(constraints) <: Tuple
        constraints = collect(constraints)
    end
    if constraints !== nothing
        @assert bin_constraints === nothing
        @assert una_constraints === nothing
        # TODO: This is redundant with the checks in EquationSearch
        for op in binary_operators
            @assert !(op in unary_operators)
        end
        for op in unary_operators
            @assert !(op in binary_operators)
        end
        bin_constraints = constraints
        una_constraints = constraints
    end

    una_constraints, bin_constraints = build_constraints(
        una_constraints, bin_constraints, unary_operators, binary_operators, nuna, nbin
    )

    # Define the complexities of everything.
    use_complexity_mapping = (
        complexity_of_constants !== nothing ||
        complexity_of_variables !== nothing ||
        complexity_of_operators !== nothing
    )
    if use_complexity_mapping
        if complexity_of_operators === nothing
            complexity_of_operators = Dict()
        else
            # Convert to dict:
            complexity_of_operators = Dict(complexity_of_operators)
        end

        # Get consistent type:
        promoted_type = promote_type(
            (complexity_of_variables !== nothing) ? typeof(complexity_of_variables) : Int,
            (complexity_of_constants !== nothing) ? typeof(complexity_of_constants) : Int,
            (x -> typeof(x)).(values(complexity_of_operators))...,
        )

        # If not in dict, then just set it to 1.
        binop_complexities = promoted_type[
            (haskey(complexity_of_operators, op) ? complexity_of_operators[op] : 1) #
            for op in binary_operators
        ]
        unaop_complexities = promoted_type[
            (haskey(complexity_of_operators, op) ? complexity_of_operators[op] : 1) #
            for op in unary_operators
        ]

        variable_complexity = (
            (complexity_of_variables !== nothing) ? complexity_of_variables : 1
        )
        constant_complexity = (
            (complexity_of_constants !== nothing) ? complexity_of_constants : 1
        )

        complexity_mapping = ComplexityMapping(;
            binop_complexities=binop_complexities,
            unaop_complexities=unaop_complexities,
            variable_complexity=variable_complexity,
            constant_complexity=constant_complexity,
        )
    else
        complexity_mapping = ComplexityMapping(false)
    end
    # Finish defining complexities

    if maxdepth === nothing
        maxdepth = maxsize
    end

    if npopulations === nothing
        npopulations = nworkers()
    end

    binary_operators = map(binopmap, binary_operators)
    unary_operators = map(unaopmap, unary_operators)

    if enable_autodiff
        diff_binary_operators = Any[]
        diff_unary_operators = Any[]

        test_inputs = map(x -> convert(Float32, x), LinRange(-100, 100, 99))
        # Create grid over [-100, 100]^2:
        test_inputs_xy = reduce(
            hcat, reduce(hcat, ([[[x, y] for x in test_inputs] for y in test_inputs]))
        )
        for op in binary_operators
            diff_op(x, y) = gradient(op, x, y)

            test_output = diff_op.(test_inputs_xy[1, :], test_inputs_xy[2, :])
            gradient_exists = all((x) -> x !== nothing, Iterators.flatten(test_output))
            if gradient_exists
                push!(diff_binary_operators, diff_op)
            else
                if verbosity > 0
                    @warn "Automatic differentiation has been turned off, since operator $(op) does not have well-defined gradients."
                end
                enable_autodiff = false
                break
            end
        end

        for op in unary_operators
            diff_op(x) = gradient(op, x)[1]
            test_output = diff_op.(test_inputs)
            gradient_exists = all((x) -> x !== nothing, test_output)
            if gradient_exists
                push!(diff_unary_operators, diff_op)
            else
                if verbosity > 0
                    @warn "Automatic differentiation has been turned off, since operator $(op) does not have well-defined gradients."
                end
                enable_autodiff = false
                break
            end
        end
        diff_binary_operators = Tuple(diff_binary_operators)
        diff_unary_operators = Tuple(diff_unary_operators)
    end

    if !enable_autodiff
        diff_binary_operators = nothing
        diff_unary_operators = nothing
    end

    mutationWeights = map((x,) -> convert(Float64, x), mutationWeights)
    if length(mutationWeights) != 8
        error("Not the right number of mutation probabilities given")
    end

    for (op, f) in enumerate(map(Symbol, binary_operators))
        _f = if f == Symbol(pow)
            Symbol(^)
        else
            f
        end
        if !isdefined(Base, _f)
            continue
        end
        @eval begin
            Base.$_f(l::Node, r::Node)::Node =
                if (l.constant && r.constant)
                    Node($f(l.val, r.val)::AbstractFloat)
                else
                    Node($op, l, r)
                end
            Base.$_f(l::Node, r::AbstractFloat)::Node =
                l.constant ? Node($f(l.val, r)::AbstractFloat) : Node($op, l, r)
            Base.$_f(l::AbstractFloat, r::Node)::Node =
                r.constant ? Node($f(l, r.val)::AbstractFloat) : Node($op, l, r)
        end
    end

    # Redefine Base operations:
    for (op, f) in enumerate(map(Symbol, unary_operators))
        if !isdefined(Base, f)
            continue
        end
        @eval begin
            Base.$f(l::Node)::Node =
                l.constant ? Node($f(l.val)::AbstractFloat) : Node($op, l)
        end
    end

    if progress
        verbosity = 0
    end

    if recorder === nothing
        recorder = haskey(ENV, "PYSR_RECORDER") && (ENV["PYSR_RECORDER"] == "1")
    end

    if typeof(earlyStopCondition) <: AbstractFloat
        # Need to make explicit copy here for this to work:
        stopping_point = Float64(earlyStopCondition)
        earlyStopCondition = (loss, complexity) -> loss < stopping_point
    end

    options = Options{
        typeof(binary_operators),
        typeof(unary_operators),
        typeof(diff_binary_operators),
        typeof(diff_unary_operators),
        typeof(loss),
        eltype(complexity_mapping),
    }(
        binary_operators,
        unary_operators,
        diff_binary_operators,
        diff_unary_operators,
        bin_constraints,
        una_constraints,
        complexity_mapping,
        ns,
        parsimony,
        alpha,
        maxsize,
        maxdepth,
        fast_cycle,
        migration,
        hofMigration,
        fractionReplacedHof,
        shouldOptimizeConstants,
        hofFile,
        npopulations,
        perturbationFactor,
        annealing,
        batching,
        batchSize,
        mutationWeights,
        crossoverProbability,
        warmupMaxsizeBy,
        useFrequency,
        useFrequencyInTournament,
        npop,
        ncyclesperiteration,
        fractionReplaced,
        topn,
        verbosity,
        probNegate,
        nuna,
        nbin,
        seed,
        loss,
        progress,
        terminal_width,
        optimizer_algorithm,
        optimize_probability,
        optimizer_nrestarts,
        optimizer_iterations,
        recorder,
        recorder_file,
        probPickFirst,
        earlyStopCondition,
        stateReturn,
        timeout_in_seconds,
        max_evals,
        skip_mutation_failures,
        enable_autodiff,
        nested_constraints,
        deterministic,
        penalties,
    )

    @eval begin
        Base.print(io::IO, tree::Node) = print(io, string_tree(tree, $options))
        Base.show(io::IO, tree::Node) = print(io, string_tree(tree, $options))
    end

    return options
end

end


text:
using InfinitesimalGenerators

# all these examples are generators with only real eigenvalues so I should choose better ones

##  Ornstein–Uhlenbeck
κx = 0.1
σ = 0.02
x = range(- 10 * sqrt(σ^2 /(2 * κx)), stop = 10 * sqrt(σ^2 /(2 * κx)), length = 1000)
μx = -κx .* x
σx = σ .* ones(length(x))

## stationnary distribution
@time g = stationary_distribution(DiffusionProcess(x, μx, σx))
#   0.002711 seconds (248 allocations: 713.797 KiB)

## Feynman-Kac
ψ = x.^2
t = range(0, stop = 1000, step = 1/10)
@time u = feynman_kac(generator(DiffusionProcess(x, μx, σx)); t = t, ψ = ψ)[:, end]
#   0.019786 seconds (3.09 k allocations: 31.120 MiB, 14.88% gc time)
g'u ≈ g'ψ

## test left and right eigenvector
κx = 0.1
σ = 0.02
x = range(- 3 * sqrt(σ^2 /(2 * κx)), stop = 3 * sqrt(σ^2 /(2 * κx)), length = 500)
μx = -κx .* x
σx = σ .* ones(length(x))
μM = -0.01 .+ x
σM = 0.1 .* ones(length(x))
ρ = 1.0
M = AdditiveFunctionalDiffusion(DiffusionProcess(x, μx, σx), μM, σM; ρ = ρ)
@time ζ = tail_index(M)
#  0.22s
@time l, η, r = cgf(M, eigenvector = :both)(ζ)
#  0.06s


text:
#####################
# REFERENCE DOMAINS #
#####################

"""
````
    reference_domain(EG::Type{<:AbstractElementGeometry}, T::Type{<:Real} = Float64; scale = [1,1,1], shift = [0,0,0]) -> ExtendableGrid{T,Int32}
````

Generates an ExtendableGrid{T,Int32} for the reference domain of the specified Element Geometry. With scale and shift the coordinates can be manipulated.
"""
function reference_domain(EG::Type{<:AbstractElementGeometry}, T::Type{<:Real} = Float64; scale = [1,1,1], shift = [0,0,0])
    xgrid=ExtendableGrid{T,Int32}()
    xCoordinates=Array{T,2}(refcoords_for_geometry(EG))
    for j = 1 : size(xCoordinates,1)
        xCoordinates[j,:] .+= shift[j]
        xCoordinates[j,:] .*= scale[j]
    end
    xgrid[Coordinates] = xCoordinates
    xCellNodes = zeros(Int32,num_nodes(EG),1)
    xCellNodes[:] = 1:num_nodes(EG)
    xgrid[CellNodes] = xCellNodes
    xgrid[CellGeometries] = VectorOfConstants(EG,1);
    xgrid[CellRegions]=ones(Int32,1)
    xgrid[BFaceRegions]=Array{Int32,1}(1:num_faces(EG))
    xgrid[BFaceNodes]=Array{Int32,2}(local_cellfacenodes(EG))
    xgrid[BFaceGeometries]=VectorOfConstants(facetype_of_cellface(EG, 1), num_faces(EG))
    if dim_element(EG) == 0
        xgrid[CoordinateSystem]=Cartesian0D
    elseif dim_element(EG) == 1
        xgrid[CoordinateSystem]=Cartesian1D
    elseif dim_element(EG) == 2
        xgrid[CoordinateSystem]=Cartesian2D
    elseif dim_element(EG) == 3
        xgrid[CoordinateSystem]=Cartesian3D
    end
    return xgrid
end

"""
````
    Triangle(coords) -> ExtendableGrid{T,Int32}
````

Generates a single triangle with the given coordinates, that should be a 2 x 3 array
with the coordinates of the three vertices, e.g. coords = [0.0 0.0; 1.0 0.0; 0.0 1.0]'.
"""
function grid_triangle(coords::AbstractArray{T,2}) where {T}
    xgrid = reference_domain(Triangle2D)
    xgrid[Coordinates]=coords
    return xgrid
end

# unit cube as one cell with six boundary regions (bottom, front, right, back, left, top)
function grid_unitcube(EG::Type{<:Hexahedron3D}; scale = [1,1,1], shift = [0,0,0])
    return reference_domain(EG; scale = scale, shift = shift)
end

# unit cube as six tets with six boundary regions (bottom, front, right, back, left, top)
function grid_unitcube(::Type{Tetrahedron3D}; scale = [1,1,1], shift = [0,0,0])
    xgrid=ExtendableGrid{Float64,Int32}()
    xCoordinates=Array{Float64,2}([0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1]')
    for j = 1 : size(xCoordinates,1)
        xCoordinates[j,:] .+= shift[j]
        xCoordinates[j,:] .*= scale[j]
    end
    xgrid[Coordinates] = xCoordinates

    xCellNodes=Array{Int32,2}([1 2 3 7; 1 3 4 7; 1 5 6 7; 1 8 5 7; 1 6 2 7;1 4 8 7]')
    xgrid[CellNodes] = xCellNodes
    xgrid[CellGeometries] = VectorOfConstants(Tetrahedron3D,6);
    ncells = num_sources(xCellNodes)
    xgrid[CellRegions]=ones(Int32,ncells)
    xgrid[BFaceRegions]=Array{Int32,1}([1,1,2,2,3,3,4,4,5,5,6,6])
    xBFaceNodes=Array{Int32,2}([1 3 2; 1 4 3; 1 2 6;1 6 5;2 3 7;2 7 6;3 4 7;7 4 8;8 4 1;1 5 8; 5 6 7; 5 7 8]')
    xgrid[BFaceNodes]=xBFaceNodes
    nbfaces = num_sources(xBFaceNodes)
    xgrid[BFaceGeometries]=VectorOfConstants(Triangle2D,nbfaces)
    xgrid[CoordinateSystem]=Cartesian3D
    return xgrid
end

# unit square as one cell with four boundary regions (bottom, right, top, left)
function grid_unitsquare(EG::Type{<:Quadrilateral2D}; scale = [1,1], shift = [0,0])
    return reference_domain(EG; scale = scale, shift = shift)
end

# unit square as two triangles with four boundary regions (bottom, right, top, left)
function grid_unitsquare(::Type{<:Triangle2D}; scale = [1,1], shift = [0,0])
    xgrid=ExtendableGrid{Float64,Int32}()
    xCoordinates=Array{Float64,2}([0 0; 1 0; 1 1; 0 1; 0.5 0.5]')
    for j = 1 : size(xCoordinates,1)
        xCoordinates[j,:] .+= shift[j]
        xCoordinates[j,:] .*= scale[j]
    end
    xgrid[Coordinates] = xCoordinates
    xgrid[CellNodes]=Array{Int32,2}([1 2 5; 2 3 5; 3 4 5; 4 1 5]')
    xgrid[CellGeometries]=VectorOfConstants(Triangle2D,4)
    xgrid[CellRegions]=ones(Int32,4)
    xgrid[BFaceRegions]=Array{Int32,1}([1,2,3,4])
    xgrid[BFaceNodes]=Array{Int32,2}([1 2; 2 3; 3 4; 4 1]')
    xgrid[BFaceGeometries]=VectorOfConstants(Edge1D,4)
    xgrid[CoordinateSystem]=Cartesian2D
    return xgrid
end

# unit square as two triangles with four boundary regions (bottom, right, top, left)
function grid_lshape(::Type{<:Triangle2D}; scale = [1,1], shift = [0,0])
    xgrid=ExtendableGrid{Float64,Int32}()
    xCoordinates=Array{Float64,2}([0 0; 1 0; 1 1; 0 1; -1 1; -1 0; -1 -1; 0 -1]')
    for j = 1 : size(xCoordinates,1)
        xCoordinates[j,:] .+= shift[j]
        xCoordinates[j,:] .*= scale[j]
    end
    xgrid[Coordinates] = xCoordinates
    xgrid[CellNodes]=Array{Int32,2}([4 2 3; 2 4 1; 6 4 5; 4 6 1; 6 8 1; 8 6 7]')
    xgrid[CellGeometries]=VectorOfConstants(Triangle2D,6)
    xgrid[CellRegions]=Array{Int32,1}([1,1,1,1,1,1])
    xgrid[BFaceRegions]=Array{Int32,1}([1,2,3,4,5,6,7,8])
    xgrid[BFaceNodes]=Array{Int32,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1]')
    xgrid[BFaceGeometries]=VectorOfConstants(Edge1D,8)
    xgrid[CoordinateSystem]=Cartesian2D
    return xgrid
end

# unit suqare as mixed triangles and squares with four boundary regions (bottom, right, top, left)
function grid_unitsquare_mixedgeometries()

    xgrid=ExtendableGrid{Float64,Int32}()
    xgrid[Coordinates]=Array{Float64,2}([0 0; 4//10 0; 1 0; 0 6//10; 4//10 6//10; 1 6//10;0 1; 4//10 1; 1 1]')
    xCellNodes=VariableTargetAdjacency(Int32)
    xCellGeometries=[Triangle2D, Triangle2D, Parallelogram2D, Parallelogram2D, Triangle2D, Triangle2D];

    append!(xCellNodes,[1,5,4])
    append!(xCellNodes,[1,2,5])
    append!(xCellNodes,[2,3,6,5])
    append!(xCellNodes,[4,5,8,7])
    append!(xCellNodes,[5,6,9])
    append!(xCellNodes,[8,5,9])

    xgrid[CellNodes] = xCellNodes
    xgrid[CellGeometries] = xCellGeometries
    ncells = num_sources(xCellNodes)
    xgrid[CellRegions]=ones(Int32,ncells)
    xgrid[BFaceRegions]=Array{Int32,1}([1,1,2,2,3,3,4,4])
    xBFaceNodes=Array{Int32,2}([1 2; 2 3; 3 6; 6 9; 9 8; 8 7; 7 4; 4 1]')
    xgrid[BFaceNodes]=xBFaceNodes
    nbfaces = num_sources(xBFaceNodes)
    xgrid[BFaceGeometries]=VectorOfConstants(Edge1D,nbfaces)
    xgrid[CoordinateSystem]=Cartesian2D

    return xgrid
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using LinearAlgebra
using Plots
using DirectTrajectoryOptimization

# ##
eval_hess_lag = true

# ## horizon
T = 11

# ## acrobot
num_state = 2
num_action = 1
num_parameter = 0

function pendulum(x, u, w)
    mass = 1.0
    length_com = 0.5
    gravity = 9.81
    damping = 0.1

    [
        x[2],
        (u[1] / ((mass * length_com * length_com))
            - gravity * sin(x[1]) / length_com
            - damping * x[2] / (mass * length_com * length_com))
    ]
end

function midpoint_implicit(y, x, u, w)
    h = 0.05 # timestep
    y - (x + h * pendulum(0.5 * (x + y), u, w))
end

# ## model
dt = Dynamics(
    midpoint_implicit,
    num_state,
    num_state,
    num_action,
    num_parameter=num_parameter,
    evaluate_hessian=eval_hess_lag)
dynamics = [dt for t = 1:T-1]

# ## initialization
x1 = [0.0; 0.0]
xT = [π; 0.0]

# ## objective
ot = (x, u, w) -> 0.1 * dot(x[1:2], x[1:2]) + 0.1 * dot(u, u)
oT = (x, u, w) -> 0.1 * dot(x[1:2], x[1:2])
ct = Cost(ot, num_state, num_action,
    num_parameter=num_parameter,
    evaluate_hessian=eval_hess_lag)
cT = Cost(oT, num_state, 0,
    num_parameter=num_parameter,
    evaluate_hessian=eval_hess_lag)
objective = [[ct for t = 1:T-1]..., cT]

# ## constraints
bnd1 = Bound(num_state, num_action)
bndt = Bound(num_state, num_action)
bndT = Bound(num_state, 0)
bounds = [bnd1, [bndt for t = 2:T-1]..., bndT]

# ## initial
con1 = Constraint((x, u, w) -> x - x1, num_state, num_action,
    evaluate_hessian=eval_hess_lag)
conT = Constraint((x, u, w) -> x - xT, num_state, num_action,
    evaluate_hessian=eval_hess_lag)
constraints = [con1, [Constraint() for t = 2:T-1]..., conT]

# ## problem
solver = Solver(dynamics, objective, constraints, bounds,
    evaluate_hessian=eval_hess_lag,
    options=Options{Float64}())

# ## initialize
x_interpolation = linear_interpolation(x1, xT, T)
u_guess = [1.0 * randn(num_action) for t = 1:T-1]

initialize_states!(solver, x_interpolation)
initialize_controls!(solver, u_guess)

# ## solve
@time solve!(solver)

# ## solution
x_sol, u_sol = get_trajectory(solver)

@show x_sol[1]
@show x_sol[T]

# ## state
plot(hcat(x_sol...)')

# ## control
plot(hcat(u_sol[1:end-1]..., u_sol[end-1])',
    linetype=:steppost)

# DirectTrajectoryOptimization.jl environment

# ## tested w/ pendulum example

# ## assemble full system

# dimensions
nz = solver.nlp.num_variables
ny = solver.nlp.num_constraint
nj = solver.nlp.num_jacobian
nh = solver.nlp.num_hessian_lagrangian

# variables
z = rand(nz)
y = rand(ny)

# data
obj_grad = zeros(nz)
con = zeros(ny)
con_jac = zeros(nj)
hess_lag = zeros(nh)

# objective
MOI.eval_objective(solver.nlp, z)
MOI.eval_objective_gradient(solver.nlp, obj_grad, z)

# constraints
MOI.eval_constraint(solver.nlp, con, z)
MOI.eval_constraint_jacobian(solver.nlp, con_jac, z)
solver.nlp.hessian_lagrangian && MOI.eval_hessian_lagrangian(solver.nlp, hess_lag, z, 1.0, y)

"""
KKT system

H = [
        ∇²L C'
        C   0
    ]

h = [
        ∇L
        c
    ]
"""

H = spzeros(nz + ny, nz + ny)
h = zeros(nz + ny)

# ∇L
for i = 1:nz
    # objective gradient
    h[i] += obj_grad[i]

    # C' y
    cy = 0.0
    for j = 1:ny
        # cy += C[j, i] * y[j]
        if (j, i) in solver.nlp.jacobian_sparsity
            cy += C[j, i] * y[j]
        end
    end
    h[i] += cy
end

# c
for j = 1:ny
    h[nz + j] += con[j]
end

#TODO: only fill in half?
# ∇²L
for (i, idx) in enumerate(solver.nlp.hessian_lagrangian_sparsity)
    H[idx...] = hess_lag[i]
end

for (i, idx) in enumerate(solver.nlp.jacobian_sparsity)
    # C
    H[nz + idx[1], idx[2]] = con_jac[i]
    # C'
    H[idx[2], nz + idx[1]] = con_jac[i]
end

# regularization
primal_reg = 1.0e-5
for i = 1:nz
    H[i, i] += primal_reg
end

dual_reg = 1.0e-5
for j = 1:ny
    H[nz + j, nz + j] -= dual_reg
end

H_dense = Array(H)
eigen(H_dense)
cond(H_dense)
rank(H_dense)

using QDLDL
H
F = qdldl(H)
sol = zeros(nz + ny)
sol = copy(h)
QDLDL.solve!(F, sol)
norm(sol - (H \ h), Inf)


text:
# api2.jl

include("../domain/domain.jl")

# test with default values
include("./myfunctions2.jl")


text:
abstract type AbstractParticleBelief{T} end

sampletype(::Type{B}) where B<:AbstractParticleBelief{T} where T = T
Random.gentype(::Type{B}) where B<:AbstractParticleBelief{T} where T = T

########################
### Belief interface ###
########################

# also rand(), pdf(), and mode() from POMDPs.jl are part of the belief interface.
"""
    n_particles(b::AbstractParticleBelief)

Return the number of particles.
"""
function n_particles end

"""
    particles(b::AbstractParticleBelief)

Return an iterator over the particles.
"""
function particles end

"""
    weights(b::AbstractParticleBelief)

Return an iterator over the weights.
"""
function weights end

"""
    weighted_particles(b::AbstractParticleBelief)

Return an iterator over particle-weight pairs.
"""
function weighted_particles end

"""
    weight_sum(b::AbstractParticleBelief)

Return the sum of the weights of the particle collection.
"""
function weight_sum end

"""
    weight(b::AbstractParticleBelief, i)

Return the weight for particle i.
"""
function weight end

"""
    particle(b::AbstractParticleBelief, i)

Return particle i.
"""
function particle end

"""
    probdict(b::AbstractParticleBelief)

Return a dictionary mapping states to probabilities.

The probability is the normalized sum of the weights for all matching particles.

For ParticleCollection and WeightedParticleBelief, the result is cached for efficiency so the calculation is only performed the first time this is called. There is a default implementation for all AbstractParticleBeliefs, but it is inefficient (it creates a new dictionary every time). New AbstractParticleBelief implementations should provide an efficient implementation.
"""
function probdict end

#############################
### Concrete Belief types ###
#############################

## Unweighted ##

"""
    ParticleCollection{S}

Unweighted particle belief consisting of equally important particles of type `S`.
"""
mutable struct ParticleCollection{T} <: AbstractParticleBelief{T}
    particles::Vector{T}
    _probs::Union{Nothing, Dict{T,Float64}} # a cache for the probabilities

    ParticleCollection{T}() where {T} = new(T[], nothing)
    ParticleCollection{T}(particles) where {T} = new(particles, Dict{T,Float64}())
    ParticleCollection{T}(particles, _probs) where {T} = new(particles, _probs)
end
ParticleCollection(p::AbstractVector{T}) where T = ParticleCollection{T}(p, nothing)

n_particles(b::ParticleCollection) = length(b.particles)
particles(b::ParticleCollection) = b.particles
weighted_particles(b::ParticleCollection) = (p=>1.0/n_particles(b) for p in particles(b))
weight_sum(::ParticleCollection) = 1.0
weight(b::ParticleCollection, i::Int) = 1.0/length(b.particles)
particle(b::ParticleCollection, i::Int) = b.particles[i]
weights(b::ParticleCollection) = fill(1/n_particles(b),n_particles(b))
rand(rng::AbstractRNG, b::ParticleCollection) = b.particles[rand(rng, 1:length(b.particles))]

mean(b::ParticleCollection) = sum(b.particles)/length(b.particles)
iterator(b::ParticleCollection) = particles(b)
function cov(b::ParticleCollection)
    diff = hcat(b.particles...) .- mean(b)
    (diff * diff') / length(b.particles)
end
support(b::ParticleCollection) = unique(particles(b))

## Weighted ##

"""
    WeightedParticleBelief{S}

Weighted particle belief consisting of particles of type `S` and their associated weights.
"""
mutable struct WeightedParticleBelief{T} <: AbstractParticleBelief{T}
    particles::Vector{T}
    weights::Vector{Float64}
    weight_sum::Float64
    _probs::Union{Nothing, Dict{T,Float64}} # this is not used now, but may be later
end
WeightedParticleBelief(particles::AbstractVector{T}, weights::AbstractVector, weight_sum=sum(weights)) where {T} = WeightedParticleBelief{T}(particles, weights, weight_sum, nothing)

n_particles(b::WeightedParticleBelief) = length(b.particles)
particles(p::WeightedParticleBelief) = p.particles
weighted_particles(b::WeightedParticleBelief) = (b.particles[i]=>b.weights[i] for i in 1:length(b.particles))
weight_sum(b::WeightedParticleBelief) = b.weight_sum
weight(b::WeightedParticleBelief, i::Int) = b.weights[i]
particle(b::WeightedParticleBelief, i::Int) = b.particles[i]
weights(b::WeightedParticleBelief) = b.weights

function Random.rand(rng::AbstractRNG, b::WeightedParticleBelief)
    t = rand(rng) * weight_sum(b)
    i = 1
    cw = b.weights[1]
    while cw < t && i < length(b.weights)
        i += 1
        @inbounds cw += b.weights[i]
    end
    return particles(b)[i]
end

mean(b::WeightedParticleBelief) = mean(b.particles, weights(b.weights))
function cov(b::WeightedParticleBelief)
    diff = hcat(b.particles...) .- mean(b)
    (diff .* b.weights') * diff' / weight_sum(b)
end

### Shared implementations ###
function probdict(b::AbstractParticleBelief{S}) where {S}
    probs = Dict{S, Float64}()
    for (i,p) in enumerate(particles(b))
        if haskey(probs, p)
            probs[p] += weight(b, i)/weight_sum(b)
        else
            probs[p] = weight(b, i)/weight_sum(b)
        end
    end
    return probs
end

function probdict(b::Union{WeightedParticleBelief{S}, ParticleCollection{S}}) where {S}
    if b._probs == nothing
        # update the cache
        probs = Dict{S, Float64}()
        for (i,p) in enumerate(particles(b))
            if haskey(probs, p)
                probs[p] += weight(b, i)/weight_sum(b)
            else
                probs[p] = weight(b, i)/weight_sum(b)
            end
        end
        b._probs = probs
    end
    return b._probs
end

pdf(b::AbstractParticleBelief{S}, s::S) where {S} = get(probdict(b), s, 0.0)
var(b::AbstractParticleBelief) = (d=diag(cov(b)); length(d)==1 ? d[1] : d)
std(b::AbstractParticleBelief) = sqrt.(var(b))

mode(b::AbstractParticleBelief) = argmax(probdict(b)) # don't know if this is the most efficient way
support(b::AbstractParticleBelief) = keys(probdict(b))


text:
using Test

@testset "GenerateJuliaPrograms" begin
  include("GenerateJuliaPrograms.jl")
end

@testset "ParseJuliaPrograms" begin
  include("ParseJuliaPrograms.jl")
end

@testset "RelationalPrograms" begin
  include("RelationalPrograms.jl")
end

@testset "DiagrammaticPrograms" begin
  include("DiagrammaticPrograms.jl")
end


text:
struct NamedSparseArray{S<:Tuple, V<:AbstractVector}
    keys::S
    values::V
    function NamedSparseArray(args::AbstractVector...)
        keys = Base.front(args)
        values = last(args)
        return new{typeof(keys), typeof(values)}(keys, values)
    end
end

# Turn named sparse array into a dense one, while remembering the labels.
function dense(kv::NamedSparseArray)
    keys, values = kv.keys, kv.values
    labels = map(collect∘uniquesorted, keys)
    indices = map(Base.OneTo∘length, labels)
    converter = map((k, v) -> Dict(zip(k, v)), labels, indices)
    d = zeros(eltype(values), indices)
    for (k, v) in zip(StructArray(keys), values)
        I = map(getindex, converter, k)
        d[I...] = v
    end
    return labels, d
end

function _frequency(args...)
    s = StructArray(args)
    gp = GroupPerm(fast_sortable(s))
    perm = sortperm(gp)
    itr = (s[perm[first(range)]] => length(range) for range in gp)
    keys, values = fieldarrays(StructArray(itr, unwrap = t -> t <: Tuple))
    namedarray =  NamedSparseArray(fieldarrays(keys)..., values)
    labels, values = dense(namedarray)
    plottypes = [:BarPlot, :Heatmap, :Volume]
    return mapping(labels..., values) * visual(plottypes[length(labels)])
end

"""
    frequency(data...)

Compute a frequency table of the arguments.
"""
const frequency = Analysis(_frequency)


text:
# Workaround for GR warnings
ENV["GKSwstype"] = "100"

using Documenter, Meshes
using DocumenterTools: Themes

istravis = "TRAVIS" ∈ keys(ENV)

Themes.compile(joinpath(@__DIR__,"src/assets/meshes-light.scss"), joinpath(@__DIR__,"src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__,"src/assets/meshes-dark.scss"), joinpath(@__DIR__,"src/assets/themes/documenter-dark.css"))

makedocs(
  format = Documenter.HTML(
    assets = ["assets/favicon.ico", asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css)],
    prettyurls = istravis,
    mathengine = KaTeX(Dict(
      :macros => Dict(
        "\\x" => "\\boldsymbol{x}",
        "\\z" => "\\boldsymbol{z}",
        "\\l" => "\\boldsymbol{\\lambda}",
        "\\c" => "\\boldsymbol{c}",
        "\\C" => "\\boldsymbol{C}",
        "\\g" => "\\boldsymbol{g}",
        "\\G" => "\\boldsymbol{G}",
        "\\f" => "\\boldsymbol{f}",
        "\\F" => "\\boldsymbol{F}",
        "\\R" => "\\mathbb{R}",
        "\\1" => "\\mathbb{1}"
      )
    ))
  ),
  sitename = "Meshes.jl",
  authors = "Júlio Hoffimann and contributors",
  pages = [
    "Home" => "index.md",
    "Reference guide" => [
      "Points" => "points.md",
      "Vectors" => "vectors.md",
      "Angles" => "angles.md",
      "Geometries" => [
        "geometries/primitives.md",
        "geometries/polytopes.md"
      ],
      "Meshes" => "meshes.md",
      "Neighborhoods" => "neighborhoods.md",
      "Algorithms" => [
        "algorithms/sampling.md",
        "algorithms/partitioning.md",
        "algorithms/discretization.md",
        "algorithms/boundbox.md",
        "algorithms/intersection.md",
        "algorithms/neighborsearch.md"
      ]
  ],
    "Contributing" => [
      "contributing/guidelines.md",
      "contributing/devguide.md"
    ],
    "About" => [
      "License" => "about/license.md"
    ],
    "Index" => "links.md"
    ]
)

deploydocs(repo="github.com/JuliaGeometry/Meshes.jl.git")


text:
using SerialPorts
using Test

@show list_serialports()

SerialPorts.in_dialout()

SerialPorts.check_serial_access()


text:
using IntArrays

function bench_getindex(array, n)
    s = zero(eltype(array))
    for i in 1:endof(array)
        s += array[i]
    end
    best = Inf
    for _ in 1:n
        best = min(best, @elapsed for i in 1:endof(array)
            s += array[i]
        end)
    end
    return best / length(array)
end

function bench_setindex(array, n)
    x = rand(0x00:0x01)
    for i in 1:endof(array)
        array[i] = x
    end
    best = Inf
    for _ in 1:n
        best = min(best, @elapsed for i in 1:endof(array)
            array[i] = x
        end)
    end
    return best / length(array)
end

function bench_fill0(array, n)
    fill!(array, 0)
    best = Inf
    for _ in 1:n
        best = min(best, @elapsed fill!(array, 0))
    end
    return best
end

function bench_fill1(array, n)
    fill!(array, 0)
    best = Inf
    for _ in 1:n
        best = min(best, @elapsed fill!(array, 1))
    end
    return best
end

function bench_copy(array, n)
    copy(array)
    best = Inf
    for _ in 1:n
        best = min(best, @elapsed copy(array))
    end
    return best
end

function bench_copysort(array, n)
    tmp = similar(array)
    sort!(copy!(tmp, array))
    best = Inf
    for _ in 1:n
        best = min(best, @elapsed sort!(copy!(tmp, array)))
    end
    return best
end

let
    srand(12345)
    baseline = !isempty(ARGS) && shift!(ARGS) == "--baseline"
    Ts = [UInt8, UInt16, UInt32, UInt64]
    size = 100_000
    small = 5
    large = 10
    benchmarks = [
        (:bench_getindex, large),
        (:bench_setindex, large),
        (:bench_fill0, small),
        (:bench_fill1, small),
        (:bench_copy, small),
        (:bench_copysort, small),
    ]
    if baseline
        columns = ["type", "eltype"]
        println(join(vcat(columns, [name for (name, _) in benchmarks]), '\t'))
        for T in Ts
            array = Vector{T}(rand(typemin(T):typemax(T), size))
            print(typeof(array), '\t', eltype(array))
            for (benchmark, n) in benchmarks
                gc()
                print('\t', eval(benchmark)(array, n))
            end
            println()
        end
    else
        columns = ["type", "eltype", "w"]
        println(join(vcat(columns, [name for (name, _) in benchmarks]), '\t'))
        for T in Ts, w in 1:sizeof(T)*8
            array = IntVector{w,T}(rand(typemin(T):typemax(T), size))
            print(typeof(array), '\t', eltype(array), '\t', w)
            for (benchmark, n) in benchmarks
                gc()
                print('\t', eval(benchmark)(array, n))
            end
            println()
        end
    end
end


text:
module Day4

export get_inputs, get_solution1, get_solution2

using UnPack: @unpack

## Input getting
maybe_parse(T, x) = x
maybe_parse(T, str::AbstractString) = parse(T, str)

fill_missing(pair_array) = fill_missing(; pair_array...)
function fill_missing(; byr=missing,
                        iyr=missing,
                        eyr=missing,
                        hgt=missing,
                        hcl=missing,
                        ecl=missing,
                        pid=missing,
                        cid=missing)

    byr, iyr, eyr = maybe_parse.(Int, (byr, iyr, eyr))
    ecl = ismissing(ecl) ? missing : Symbol(ecl)

    return (; byr, iyr, eyr, hgt, hcl, ecl, pid, cid)
end

parse_passport(str) = split(str) .|>
    Base.Fix2(split, ":") .|>
    (x -> Symbol(x[1]) => x[2]) |>
    fill_missing

read_input(f_name) = read(joinpath(@__DIR__, f_name), String) |>
    Base.Fix2(split, r"\n\n|\r\n\r\n") .|>
    parse_passport

get_inputs() = (;
    test_input1 = read_input("test_input1.txt"),
    test_input2 = read_input("test_input2.txt"),
    test_output1 = 2,
    test_output2 = 4,
    data = read_input("input.txt"),
)

## Solution functions
# Part 1
is_valid1(passport) = all(.!ismissing.(getproperty.(Ref(passport), (:byr, :iyr, :eyr, :hgt, :hcl, :ecl, :pid))))

get_solution1(data) = count(is_valid1, data)

# Part 2
function is_valid2(passport)
    @unpack byr, iyr, eyr, hgt, hcl, ecl, pid = passport

    return (
        is_valid1(passport) &&
        ndigits(byr)==4 && 1920≤byr≤2002 &&
        ndigits(iyr)==4 && 2010≤iyr≤2020 &&
        ndigits(eyr)==4 && 2020≤eyr≤2030 &&
        (
            (contains(hgt, "cm") && 150 ≤ parse(Int, hgt[1:end-2]) ≤ 193) ||
            (contains(hgt, "in") &&  59 ≤ parse(Int, hgt[1:end-2]) ≤  76)
        ) &&
        (
            hcl[1]=='#' &&
            length(hcl)==7 &&
            all(('0'≤x≤'9' || 'a'≤x≤'f') for x in hcl[2:end])
        ) &&
        ecl in (:amb, :blu, :brn, :gry, :grn, :hzl, :oth) &&
        length(pid)==9
    )
end

get_solution2(data) = count(is_valid2, data)

end


text:
module Cracker

using ChainRules: rrule, unthunk
import ChainRules
using LinearAlgebra

export track, untrack

include("types/record.jl")
include("types/number.jl")
include("types/array.jl")

is_tracked(x::Tuple) = any(is_tracked, x)
untrack(x::Tuple) = untrack.(x)
function track(A::Tuple, record::Record=leaf(A))
    return track.(A, Ref(record))
end

const TrackedType = Union{TrackedNumber, TrackedArrayType}
function Base.show(io::IO, x::TrackedType)
    print(io, "track(")
    show(io, x.value)
    print(io, ")")
end

include("trace.jl")
include("rrule.jl")

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2021, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#
import Modia3D.Frames

"""
    joint = Revolute(;obj1, obj2, path="", axis=3, phi=0, w=0, canCollide=false)

Return a Revolute `joint` that rotates `obj1::`[`Object3D`](@ref) into
`obj2::`[`Object3D`](@ref) along the axis `axis` of `obj1` (`axis = 1,2,3,-1,-2,-3`).
The initial start angle is `phi` and the initial angular velocity
is `w`. Negative `axis` values describe axes in negative axis directions, i.e.
the signs of the values of `phi` and `w` are reversed. If `canCollide=false`, no
collision detection will occur between `obj1` and `obj2` (and `Object3D`s that are
directly or indirectly rigidly fixed to `obj1` or `obj2`).

It is currently not supported that a `Revolute` joint *closes a kinematic loop*.
"""
mutable struct Revolute <: Modia3D.AbstractJoint
    path::String

    obj1::Object3D
    obj2::Object3D

    posAxis::Int      # = 1,2,3
    posMovement::Bool # = true, if movement along posAxis, otherwise in negative posAxis
    ndof::Int
    canCollide::Bool  # = false, if no collision between obj1 and obj2

    phi::Float64
    w::Float64
    a::Float64
    tau::Float64
    residue::Float64

    function Revolute(; obj1::Object3D,
                        obj2::Object3D,
                        path::String="",
                        axis::Int = 3,
                        phi::Real = 0.0,
                        w::Real   = 0.0,
                        canCollide::Bool = false)

        (parent,obj,cutJoint) = attach(obj1, obj2)
        if cutJoint
            error("\nError from Revolute joint connecting ", Modia3D.fullName(obj1), " with ", Modia3D.fullName(obj2), ":\n",
                "    This joint is a cut-joint which is currently not supported.!")
        end

        if !(1 <= abs(axis) <= 3)
            error("\nError from Revolute joint connecting ", Modia3D.fullName(obj1), " with ", Modia3D.fullName(obj2), ":\n",
                "    axis = $axis, but must be 1,2,3 or -1,-2-3.!")
        end

        if !(typeof(phi) <: AbstractFloat) || !(typeof(w) <: AbstractFloat)
            error("\nError from Revolute joint connecting ", Modia3D.fullName(obj1), " with ", Modia3D.fullName(obj2), ":\n",
                "    phi and/or w are not <: AbstractFloat!")
        end

        axis = obj === obj2 ? axis : -axis
        phi = Modia3D.convertAndStripUnit(Float64, u"rad"  , phi)
        w   = Modia3D.convertAndStripUnit(Float64, u"rad/s", w)

        posAxis     = abs(axis)
        posMovement = axis > 0

        obj.joint = new(path, parent, obj, posAxis, posMovement, 1, canCollide, phi, w, 0.0, 0.0, 0.0)
        obj.jointKind  = RevoluteKind
        obj.jointIndex = 0
        obj.ndof       = 1
        obj.canCollide = canCollide
        obj.r_rel      = Modia3D.ZeroVector3D
        obj.R_rel      = Frames.rotAxis(posAxis, posMovement, phi)

        parent.hasChildJoint = true
        return obj.joint
    end
end


text:
module SplineInterpolations

using ..BSplines
using ..Collocation
using ..Splines

using BandedMatrices
using LinearAlgebra

import ..BSplines:
    order, knots, basis

import ..Splines:
    coefficients, integral

import Interpolations:
    interpolate, interpolate!

export
    SplineInterpolation, spline, interpolate, interpolate!

"""
    SplineInterpolation

Spline interpolation.

This is the type returned by [`interpolate`](@ref).

A `SplineInterpolation` `I` can be evaluated at any point `x` using the `I(x)`
syntax.

It can also be updated with new data on the same data points using
[`interpolate!`](@ref).

---

    SplineInterpolation(undef, B::AbstractBSplineBasis, x::AbstractVector, [T = eltype(x)])

Initialise a `SplineInterpolation` from B-spline basis and a set of
interpolation (or collocation) points `x`.

Note that the length of `x` must be equal to the number of B-splines.

Use [`interpolate!`](@ref) to actually interpolate data known on the `x`
locations.
"""
struct SplineInterpolation{
        S <: Spline,
        F <: Factorization,
        Points <: AbstractVector,
    } <: SplineWrapper{S}

    spline :: S
    C :: F  # factorisation of collocation matrix
    x :: Points

    function SplineInterpolation(
            B::AbstractBSplineBasis, C::Factorization, x::AbstractVector,
        )
        N = length(B)
        size(C) == (N, N) ||
            throw(DimensionMismatch("collocation matrix has wrong dimensions"))
        length(x) == N ||
            throw(DimensionMismatch("wrong number of collocation points"))
        T = eltype(C)
        s = Spline(undef, B, T)  # uninitialised spline
        new{typeof(s), typeof(C), typeof(x)}(s, C, x)
    end
end

# Construct SplineInterpolation from basis and collocation points.
function SplineInterpolation(
        ::UndefInitializer, B, x::AbstractVector, ::Type{T},
    ) where {T}
    # Here we construct the collocation matrix and its LU factorisation.
    N = length(B)
    if length(x) != N
        throw(DimensionMismatch(
            "incompatible lengths of B-spline basis and collocation points"))
    end
    C = collocation_matrix(B, x, CollocationMatrix{T})
    SplineInterpolation(B, lu!(C), x)
end

interpolation_points(S::SplineInterpolation) = S.x

SplineInterpolation(init, B, x::AbstractVector) =
    SplineInterpolation(init, B, x, eltype(x))

function Base.show(io::IO, I::SplineInterpolation)
    println(io, nameof(typeof(I)), " containing the ", spline(I))
    let io = IOContext(io, :compact => true, :limit => true)
        println(io, " interpolation points: ", interpolation_points(I))
    end
    nothing
end

"""
    interpolate!(I::SplineInterpolation, y::AbstractVector)

Update spline interpolation with new data.

This function allows to reuse a [`SplineInterpolation`](@ref) returned by a
previous call to [`interpolate`](@ref), using new data on the same locations
`x`.

See [`interpolate`](@ref) for details.
"""
function interpolate!(I::SplineInterpolation, y::AbstractVector)
    s = spline(I)
    if length(y) != length(s)
        throw(DimensionMismatch("input data has incorrect length"))
    end
    ldiv!(coefficients(s), I.C, y)  # determine B-spline coefficients
    I
end

"""
    interpolate(x, y, BSplineOrder(k))

Interpolate values `y` at locations `x` using B-splines of order `k`.

Grid points `x` must be real-valued and are assumed to be in increasing order.

Returns a [`SplineInterpolation`](@ref) which can be evaluated at any intermediate
point.

See also [`interpolate!`](@ref).

# Examples

```jldoctest
julia> xs = -1:0.1:1;

julia> ys = cospi.(xs);

julia> itp = interpolate(xs, ys, BSplineOrder(4))
SplineInterpolation containing the 21-element Spline{Float64}:
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3  …  0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-1.0, -1.00111, -0.8975, -0.597515, -0.314147, 1.3265e-6, 0.314142, 0.597534, 0.822435, 0.96683  …  0.96683, 0.822435, 0.597534, 0.314142, 1.3265e-6, -0.314147, -0.597515, -0.8975, -1.00111, -1.0]
 interpolation points: -1.0:0.1:1.0

julia> itp(-1)
-1.0

julia> itp(0)
1.0

julia> itp(0.42)
0.2486897676885842
```
"""
function interpolate(x::AbstractVector, y::AbstractVector, k::BSplineOrder)
    t = make_knots(x, order(k))
    B = BSplineBasis(k, t; augment = Val(false))  # it's already augmented!
    T = float(eltype(y))
    itp = SplineInterpolation(undef, B, x, T)
    interpolate!(itp, y)
end

# Define B-spline knots from collocation points and B-spline order.
# Note that the choice is not unique.
function make_knots(x::AbstractVector{Tx}, k) where {Tx <: Real}
    Base.require_one_based_indexing(x)
    T = float(Tx)  # just in case Tx is Integer...
    N = length(x)
    t = Array{T}(undef, N + k)  # knots

    # First and last knots have multiplicity `k`.
    t[1:k] .= x[1]
    t[(N + 1):(N + k)] .= x[end]

    # Use initial guess for optimal set of knot locations (de Boor 2001, p. 193).
    # One could get the actual optimal using Newton iterations...
    kinv = one(T) / (k - 1)
    @inbounds for i = 1:(N - k)
        ti = zero(T)
        for j = 1:(k - 1)
            ti += x[i + j]
        end
        t[k + i] = ti * kinv
    end

    t
end

end


text:
###############
###############
###############

T = 100                                 # Number of periods in the data
P = 2                                   # Number of lags in the factor equation
L = 1                                   # Number of AR lags in the observation equation
N = 10                                  # Total number of series
K = 1                                   # Total number of factors (only the global factor)

H = [ones(N) ident(N) zeros(N)]
R = zeros(N, N)
A = zeros(N, 2)
Z = zeros(2,2)

F = zeros(P + N, P + N)
F[1,1] = 0.5
F[1,P+N] = 0.1
F[P+N,1] = 1.0
F[2:(P+N-1), 2:(P+N-1)] = 0.25 * ident(N)

Q = zeros(P + N, P + N)
Q[1,1] = 1.0
for i in 1:N
    Q[1+i,1+i] = 1.0
end

μ = zeros(P+N)

num_obs = T

# Gather all SS parameters
ssmodel = SSModel(H, A, F, μ, R, Q, Z)

# Simulate common components model in state space form
data_y, data_z, data_β = simulateSSModel(num_obs, ssmodel)

priors = DFMStruct(P, L, 1000, 50)             # Set model priors

###############
###############
###############

results = OW1LevelEstimator(data_y, priors)

###############
###############
###############

medians = Any[]
quant33 = Any[]
quant66 = Any[]
stds = Any[]

j = 1
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, j, :]))
    push!(quant33, quantile(results.F[i, :], 0.05))
    push!(quant66, quantile(results.F[i, :], 0.95))
    push!(medians, median(results.F[i, :]))
end

plot(results.means.F[:, j])
plot(medians)
plot!(quant33)
plot!(quant66)


text:
const USE_GPU  = false  # Use GPU? If this is set false, then the CUDA packages do not need to be installed! :)
const GPU_ID   = 0
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    CUDA.device!(GPU_ID) # select GPU
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Test, Plots, Printf, Statistics, LinearAlgebra

################################
@parallel function def_err!(errD::Data.Array, D::Data.Array)

    @all(errD) = @all(D)

    return
end

@parallel function chk_err!(errD::Data.Array, D::Data.Array)

    @all(errD) = @all(D) - @all(errD)

    return
end

@parallel function compute_ϕ!(ϕ::Data.Array, rdh::Data.Number, rr::Data.Number, Zb::Data.Array, D::Data.Array, H::Data.Array)

    @all(ϕ) = @all(Zb) + rdh*@all(D) + rr*@all(H)

    return
end

@parallel_indices (ix,iy) function compute_flux!(Ux::Data.Array, Uy::Data.Array, qDx::Data.Array, qDy::Data.Array, k::Data.Number, _dx::Data.Number, _dy::Data.Number, ϕ::Data.Array, D::Data.Array)

    if (ix<=size(Ux,1)  && iy<=size(Ux,2))    Ux[ix,iy] = -k*_dx*(ϕ[ix+1,iy]-ϕ[ix,iy]) end
    if (ix<=size(Uy,1)  && iy<=size(Uy,2))    Uy[ix,iy] = -k*_dy*(ϕ[ix,iy+1]-ϕ[ix,iy]) end
    if (ix<=size(qDx,1) && iy<=size(qDx,2))  qDx[ix,iy] = Ux[ix,iy]*(Ux[ix,iy]>0.0)*D[ix,iy] + Ux[ix,iy]*(Ux[ix,iy]<0.0)*D[ix+1,iy  ] end
    if (ix<=size(qDy,1) && iy<=size(qDy,2))  qDy[ix,iy] = Uy[ix,iy]*(Uy[ix,iy]>0.0)*D[ix,iy] + Uy[ix,iy]*(Uy[ix,iy]<0.0)*D[ix  ,iy+1] end

    return
end

@parallel function compute_∂D!(∂D::Data.Array, _dx::Data.Number, _dy::Data.Number, M::Data.Number, qDx::Data.Array, qDy::Data.Array)

    @inn(∂D) = -(_dx*@d_xi(qDx) + _dy*@d_yi(qDy)) + 1.0/M

    return
end

@parallel function update_D!(∂Ddτ::Data.Array, D::Data.Array, damp::Data.Number, dτ::Data.Number, ∂D::Data.Array)

    @all(∂Ddτ) =  damp*@all(∂Ddτ) + @all(∂D)
    @all(D)    =  @all(D) + dτ*@all(∂Ddτ)

    return
end

@parallel_indices (ix,iy) function set_BCx!(D::Data.Array, val_BC::Data.Number)

    if (ix==1         && iy<=size(D,2)) D[ix,iy] = val_BC end
    if (ix==size(D,1) && iy<=size(D,2)) D[ix,iy] = val_BC end

    return
end

@parallel_indices (ix,iy) function set_BCy!(D::Data.Array, val_BC::Data.Number)
    # fixed height
    if (ix<=size(D,1) && iy==1        ) D[ix,iy] = val_BC end
    if (ix<=size(D,1) && iy==size(D,2)) D[ix,iy] = val_BC end
    # no flux
    # if (ix<=size(D,1) && iy==1        ) D[ix,iy] = D[ix,iy+1] end
    # if (ix<=size(D,1) && iy==size(D,2)) D[ix,iy] = D[ix,iy-1] end
    return
end

# physics scales
const ρ̂i  = 910
const ρ̂w  = 1000
const ĝ   = 9.81
const ϕ̂   = 1e7
const x̂   = 1e6
const M̂   = 1e-6
const D̂   = 1
# physics dependent
const t̂   = D̂/M̂
const Ĥ   = ϕ̂/ρ̂w/ĝ
const Ψ̂   = ϕ̂/x̂
const û   = x̂/t̂
const k̂   = û/Ψ̂
const rr  = ρ̂i/ρ̂w
const rdh = D̂/Ĥ    # fixed ice thikness routing
const s2d = 3600*24

@views function flow_routing2D()
    # nondim
    Lx, Ly = 10e3/x̂, 10e3/x̂
    k      = 0.1
    M      = M̂*s2d
    # numerics
    res    = 2
    nx     = res*128 # fastest if multiple of 16 (as no overlength here)
    ny     = res*128 # fastest if multiple of 16 (as no overlength here)
    nt     = 1
    itrMax = 1e8
    nout   = 1000
    nmax   = 100
    ε      = 1e-8     # aboslute tolerance
    damp   = 0.94  #0.9 for nx=100
    dτ_sc  = res*30.0 #60 for nx,ny=256 #30 for nx,ny=128 #24 for nx,ny=96 #12.0 for nx,ny=64
    dx, dy = Lx/nx, Ly/ny
    _dx, _dy = 1.0/dx, 1.0/dy
    xc     = LinRange(dx/2, Lx-dx/2, nx)
    yc     = LinRange(dy/2, Ly-dy/2, ny)
    # initial
    Zb     = (Lx.-xc)./Lx./100.0 .- 30.0./Ĥ.*exp.(-(xc.-Lx/2.0).^2/(Lx/8.0)^2 .-(yc'.-Ly/2.0).^2/(Ly/8.0)^2)
    H      = 100.0./Ĥ .- xc.*3.0 .+ 0*yc'
    Zb     = Data.Array(Zb)
    H      = Data.Array(H)
    # H      = (Lx-xc)/Lx/100 + 30/H_s*exp(-(xc-Lx/2).^2/(Lx/8)^2) + (xc + 50)/H_s;
    D      = 1.0.*@ones(nx,ny)
    # D      = max(Zb)-Zb + 10;
    ∂D     = @zeros(nx  ,ny  )
    ϕ      = @zeros(nx  ,ny  )
    Ux     = @zeros(nx-1,ny  )
    Uy     = @zeros(nx  ,ny-1)
    qDx    = @zeros(nx-1,ny  )
    qDy    = @zeros(nx  ,ny-1)
    ∂Ddτ   = @zeros(nx  ,ny  )
    errD   = @zeros(nx  ,ny  )
    err1=[]; err2=[]
    # action
    iter = 1; err = 2*ε;  max_U = 1.0
    while err > ε && iter < itrMax
    	@parallel def_err!(errD, D)
    	# flow routing physics
    	@parallel compute_ϕ!(ϕ, rdh, rr, Zb, D, H)
    	@parallel compute_flux!(Ux, Uy, qDx, qDy, k, _dx, _dy, ϕ, D)
    	@parallel compute_∂D!(∂D, _dx, _dy, M, qDx, qDy)
    	if mod(iter,nmax)==0 || iter==1
    		max_U = max(maximum(abs.(Ux)), maximum(abs.(Uy)))
    	end
    	dτ = min(dx,dy)/max_U/dτ_sc
    	@parallel update_D!(∂Ddτ, D, damp, dτ, ∂D)
    	@parallel set_BCx!(D, 0.0)
    	@parallel set_BCy!(D, 0.0)
    	# Check errs
    	if mod(iter,nout)==0
    		@parallel chk_err!(errD, D)
    		ermb  = 1.0./M*dx*dy*(nx-2)*(ny-2) - sum(abs.(qDx[[1, end],2:end-1])) - sum(abs.(qDy[2:end-1,[1, end]]))
    		push!(err1, maximum(abs.(errD[:]))); push!(err2, abs(ermb))
    		err   = err1[end]
    		@printf("iter=%d  errD=%1.3e, errMB=%1.3e \n", iter, err1[end], err2[end])
    	end
    	iter+=1
    end
    # ploting
    xcp = xc*x̂/1e3; ycp = yc*x̂/1e3
    p1 = heatmap(xcp, ycp, Zb'*Ĥ, aspect_ratio=1, xlims=(xcp[1], xcp[end]), ylims=(ycp[1], ycp[end]), c=:inferno, title="Zb")
    # p1 = heatmap(transpose(qDx), aspect_ratio=1)
    p2 = heatmap(xcp, ycp, D'*D̂ + Zb'*Ĥ, aspect_ratio=1, xlims=(xcp[1], xcp[end]), ylims=(ycp[1], ycp[end]), c=:inferno, title="Zb+D")
    p3 = plot(xcp[2:end-1], D[2:end-1,Int(round(ny/2))]*D̂, ylabel="D [m]", yscale=:log10, linewidth=2, framestyle=:box, legend=false)
    l  = @layout [a b; c]
    # display(plot(p1, p2, p3, layout = l))
    savefig(plot(p1, p2, p3, layout = l), joinpath(@__DIR__, "output/o$nx.png"))
    return
end

@time flow_routing2D()


text:
@testset "sizing" begin
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 0, 0, .1))
    test_ref("references/spy/default_0x0.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 10, 10, .15))
    @test p isa Plot
    test_ref("references/spy/default_10x10.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 10, 15, .15))
    test_ref("references/spy/default_10x15.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 15, 10, .15))
    test_ref("references/spy/default_15x10.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 2000, 200, .0001))
    test_ref("references/spy/default_2000x200.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 200, 2000, .0001))
    test_ref("references/spy/default_200x2000.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 2000, 2000, .1))
    test_ref("references/spy/default_overdrawn.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 200, 200, .001))
    test_ref("references/spy/default_200x200_normal.txt", @show_col(p, :displaysize=>T_SZ))
    test_ref("references/spy/default_200x200_normal_nocolor.txt", @show_nocol(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(Matrix(stable_sprand(RNG, 200, 200, .001)))
    test_ref("references/spy/default_200x200_normal.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 200, 200, .001), width=10)
    test_ref("references/spy/default_200x200_normal_small.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 200, 200, .001), height=5)
    test_ref("references/spy/default_200x200_normal_small.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 200, 200, .001), height=5, width=20)
    test_ref("references/spy/default_200x200_normal_misshaped.txt", @show_col(p, :displaysize=>T_SZ))
end

@testset "parameters" begin
    seed!(RNG, 1337)
    p = @inferred spy(stable_sprand(RNG, 200, 200, .001), color=:green)
    test_ref("references/spy/parameters_200x200_green.txt", @show_col(p, :displaysize=>T_SZ))
    test_ref("references/spy/parameters_200x200_green_nocolor.txt", @show_nocol(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = spy(stable_sprand(RNG, 200, 200, .001), title="Custom Title", canvas=DotCanvas, border=:ascii)
    test_ref("references/spy/parameters_200x200_dotcanvas.txt", @show_col(p, :displaysize=>T_SZ))
end


text:
@testset "fastmod" begin
    a = convert(Vector{UInt32},sample(1:typemax(UInt32),10_000))
    a2 = copy(a)
    b = convert(Vector{UInt32},sample(1:typemax(UInt32),1_000))

    @testset for m in b
            modm = Mod(m)
            for i in eachindex(a)
                c1 = mod(a[i],m)
                c2 = fastmod(a[i],modm)
                @test c1==c2
            end
    end
end

@testset "fastmod!" begin
    a = convert(Vector{UInt32},sample(1:typemax(UInt32),10_000))
    a2 = copy(a)
    b = convert(Vector{UInt32},sample(1:typemax(UInt32),1_000))

    @testset for m in b
        modm = Mod(m)
        fastmod!(a,modm)
        inplace_mod!(a2,m)
        @test all(a.==a2)
    end
end

@testset "divtest" begin
    a = convert(Vector{UInt32},sample(1:typemax(UInt32),10_000))
    a2 = copy(a)
    b = convert(Vector{UInt32},sample(1:typemax(UInt32),1_000))
    c1 = similar(a)
    c2 = similar(a)

    @testset for m in b
        modm = Mod(m)
        for i in eachindex(a)
            c1[i] = (mod(a[i],m) == 0)
            c2[i] = (divtest(a[i],modm))
        end
        @test all(c1.==c2)
    end
end