text:
# The SLEEF.jl package is licensed under the MIT "Expat" License:

# > Copyright (c) 2016: Mustafa Mohamad and other contributors:
# >
# > https://github.com/musm/SLEEF.jl/graphs/contributors
# >
# > Permission is hereby granted, free of charge, to any person obtaining a copy
# > of this software and associated documentation files (the "Software"), to deal
# > in the Software without restriction, including without limitation the rights
# > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# > copies of the Software, and to permit persons to whom the Software is
# > furnished to do so, subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be included in all
# > copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.
# >

# SLEEF.jl includes ported code from the following project

#     - [SLEEF](https://github.com/shibatch/SLEEF) [public domain] Author Naoki Shibata

using Base.Math: IEEEFloat
for (op, f, ff) ∈ [
  ("fadd", :add_ieee, :(+)),
  ("fsub", :sub_ieee, :(-)),
  ("fmul", :mul_ieee, :(*)),
  ("fdiv", :fdiv_ieee, :(/)),
  ("frem", :rem_ieee, :(%)),
]
  @eval begin
    @generated $f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} =
      VectorizationBase.binary_op($op, W, T)
    @inline $f(s1::T, s2::T) where {T<:Union{Float32,Float64}} = $ff(s1, s2)
    @inline $f(args::Vararg{Any,K}) where {K} = $f(promote(args...)...)
    @inline $f(a::VecUnroll, b::VecUnroll) = VecUnroll(
      VectorizationBase.fmap($f, VectorizationBase.data(a), VectorizationBase.data(b)),
    )
  end
end
@inline add_ieee(a, b, c) = add_ieee(add_ieee(a, b), c)
@inline add_ieee(a, b, c, d::Vararg{Any,K}) where {K} =
  add_ieee(add_ieee(a, b), add_ieee(c, d...))
function sub_ieee!(ex)
  ex isa Expr || return
  if ex.head === :call
    _f = ex.args[1]
    if _f isa Symbol
      f::Symbol = _f
      if f === :(+)
        ex.args[1] = :(VectorizationBase.add_ieee)
      elseif f === :(-)
        ex.args[1] = :(VectorizationBase.sub_ieee)
      elseif f === :(*)
        ex.args[1] = :(VectorizationBase.mul_ieee)
      elseif f === :(/)
        ex.args[1] = :(VectorizationBase.fdiv_ieee)
      elseif f === :(%)
        ex.args[1] = :(VectorizationBase.rem_ieee)
      end
    end
  end
  foreach(sub_ieee!, ex.args)
  esc(ex)
end
macro ieee(ex)
  sub_ieee!(ex)
end

const vIEEEFloat = Union{
  IEEEFloat,
  Vec{<:Any,<:IEEEFloat},
  VectorizationBase.VecUnroll{<:Any,<:Any,<:IEEEFloat},
}

struct Double{T<:vIEEEFloat} <: Number
  hi::T
  lo::T
end
@inline Double(x::T) where {T<:vIEEEFloat} = Double(x, zero(T))
@inline function Double(x::Vec, y::Vec)
  Double(Vec(data(x)), Vec(data(y)))
end
@inline Base.convert(::Type{Double{V}}, v::Vec) where {W,T,V<:AbstractSIMD{W,T}} =
  Double(convert(V, v), vzero(V))
@inline Base.convert(::Type{Double{V}}, v::V) where {V<:AbstractSIMD} = Double(v, vzero(V))
# @inline Base.convert(::Type{Double{V}}, m::Mask) where {V} = m
# @inline Base.convert(::Type{Double{Mask{W,U}}}, m::Mask{W,U}) where {W,U} = m
@inline Base.convert(::Type{Double{V}}, d::Double{T}) where {W,T,V<:AbstractSIMD{W,T}} =
  Double(vbroadcast(Val{W}(), d.hi), vbroadcast(Val{W}(), d.lo))
@inline Base.eltype(d::Double) = eltype(d.hi)

(::Type{T})(x::Double{T}) where {T<:vIEEEFloat} = x.hi + x.lo

Base.issubnormal(d::Double) = issubnormal(d.hi) | issubnormal(d.lo)

@inline Base.eltype(d::Double{T}) where {T<:IEEEFloat} = T
@inline Base.eltype(d::Double{S}) where {N,T,S<:Union{Vec{N,T},Vec{N,T}}} = T
@inline ifelse(u::Bool, v1::Double, v2::Double) =
  Double(ifelse(u, v1.hi, v2.hi), ifelse(u, v1.lo, v2.lo))
@inline ifelse(u::Mask, v1::Double, v2::Double) =
  Double(ifelse(u, v1.hi, v2.hi), ifelse(u, v1.lo, v2.lo))
@generated function ifelse(
  m::VecUnroll{N,W,T},
  v1::Double{V1},
  v2::Double{V2},
) where {N,W,T,V1,V2}
  q = Expr(
    :block,
    Expr(:meta, :inline),
    :(md = data(m)),
    :(v1h = v1.hi),
    :(v2h = v2.hi),
    :(v1l = v1.lo),
    :(v2l = v2.lo),
  )
  if V1 <: VecUnroll
    push!(q.args, :(v1hd = data(v1h)))
    push!(q.args, :(v1ld = data(v1l)))
  end
  if V2 <: VecUnroll
    push!(q.args, :(v2hd = data(v2h)))
    push!(q.args, :(v2ld = data(v2l)))
  end
  th = Expr(:tuple)
  tl = Expr(:tuple)
  gf = GlobalRef(Core, :getfield)
  for n ∈ 1:N+1
    ifelseₕ = Expr(:call, :ifelse, Expr(:call, gf, :md, n, false))
    ifelseₗ = Expr(:call, :ifelse, Expr(:call, gf, :md, n, false))
    if V1 <: VecUnroll
      push!(ifelseₕ.args, Expr(:call, gf, :v1hd, n, false))
      push!(ifelseₗ.args, Expr(:call, gf, :v1ld, n, false))
    else
      push!(ifelseₕ.args, :v1h)
      push!(ifelseₗ.args, :v1l)
    end
    if V2 <: VecUnroll
      push!(ifelseₕ.args, Expr(:call, gf, :v2hd, n, false))
      push!(ifelseₗ.args, Expr(:call, gf, :v2ld, n, false))
    else
      push!(ifelseₕ.args, :v2h)
      push!(ifelseₕ.args, :v2l)
    end
    push!(th.args, ifelseₕ)
    push!(tl.args, ifelseₗ)
  end
  push!(q.args, :(Double(VecUnroll($th), VecUnroll($tl))))
  q
end
@inline trunclo(x::Float64) =
  reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000) # clear lower 27 bits (leave upper 26 bits)
@inline trunclo(x::Float32) = reinterpret(Float32, reinterpret(UInt32, x) & 0xffff_f000) # clear lowest 12 bits (leave upper 12 bits)

# @inline trunclo(x::VecProduct) = trunclo(Vec(data(x)))
@inline function trunclo(x::AbstractSIMD{N,Float64}) where {N}
  reinterpret(
    Vec{N,Float64},
    reinterpret(Vec{N,UInt64}, x) & convert(Vec{N,UInt64}, 0xffff_ffff_f800_0000),
  ) # clear lower 27 bits (leave upper 26 bits)
end
@inline function trunclo(x::AbstractSIMD{N,Float32}) where {N}
  reinterpret(
    Vec{N,Float32},
    reinterpret(Vec{N,UInt32}, x) & convert(Vec{N,UInt32}, 0xffff_f000),
  ) # clear lowest 12 bits (leave upper 12 bits)
end

@inline function splitprec(x::vIEEEFloat)
  hx = trunclo(x)
  hx, x - hx
end

@inline function dnormalize(x::Double{T}) where {T}
  r = x.hi + x.lo
  Double(r, (x.hi - r) + x.lo)
end

@inline Base.flipsign(x::Double{<:vIEEEFloat}, y::vIEEEFloat) =
  Double(flipsign(x.hi, y), flipsign(x.lo, y))

@inline scale(x::Double{<:vIEEEFloat}, s::vIEEEFloat) = Double(s * x.hi, s * x.lo)

@inline Base.:(-)(x::Double{T}) where {T<:vIEEEFloat} = Double(-x.hi, -x.lo)

@inline function Base.:(<)(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat})
  x.hi < y.hi
end

@inline function Base.:(<)(x::Double{<:vIEEEFloat}, y::Union{Number,Vec})
  x.hi < y
end

@inline function Base.:(<)(x::Union{Number,Vec}, y::Double{<:vIEEEFloat})
  x < y.hi
end

# quick-two-sum x+y
@inline function dadd(x::vIEEEFloat, y::vIEEEFloat) #WARNING |x| >= |y|
  s = x + y
  Double(s, ((x - s) + y))
end

@inline function dadd(x::vIEEEFloat, y::Double{<:vIEEEFloat}) #WARNING |x| >= |y|
  s = x + y.hi
  Double(s, (((x - s) + y.hi) + y.lo))
end

@inline function dadd(x::Double{<:vIEEEFloat}, y::vIEEEFloat) #WARNING |x| >= |y|
  s = x.hi + y
  Double(s, (((x.hi - s) + y) + x.lo))
end

@inline function dadd(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat}) #WARNING |x| >= |y|
  s = x.hi + y.hi
  Double(s, ((((x.hi - s) + y.hi) + y.lo) + x.lo))
end

@inline function dsub(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat}) #WARNING |x| >= |y|
  s = x.hi - y.hi
  Double(s, ((((x.hi - s) - y.hi) - y.lo) + x.lo))
end

@inline function dsub(x::Double{<:vIEEEFloat}, y::vIEEEFloat) #WARNING |x| >= |y|
  s = x.hi - y
  Double(s, (((x.hi - s) - y) + x.lo))
end

@inline function dsub(x::vIEEEFloat, y::Double{<:vIEEEFloat}) #WARNING |x| >= |y|
  s = x - y.hi
  Double(s, (((x - s) - y.hi - y.lo)))
end

@inline function dsub(x::vIEEEFloat, y::vIEEEFloat) #WARNING |x| >= |y|
  s = x - y
  Double(s, ((x - s) - y))
end

# two-sum x+y  NO BRANCH
@inline function dadd2(x::vIEEEFloat, y::vIEEEFloat)
  s = x + y
  v = s - x
  Double(s, ((x - (s - v)) + (y - v)))
end

@inline function dadd2(x::vIEEEFloat, y::Double{<:vIEEEFloat})
  s = x + y.hi
  v = s - x
  Double(s, (x - (s - v)) + (y.hi - v) + y.lo)
end

@inline dadd2(x::Double{<:vIEEEFloat}, y::vIEEEFloat) = dadd2(y, x)

@inline function dadd2(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat})
  s = (x.hi + y.hi)
  v = (s - x.hi)
  smv = (s - v)
  yhimv = (y.hi - v)
  Double(s, ((((x.hi - smv) + yhimv) + x.lo) + y.lo))
end

@inline function dsub2(x::vIEEEFloat, y::vIEEEFloat)
  s = x - y
  v = s - x
  Double(s, ((x - (s - v)) - (y + v)))
end

@inline function dsub2(x::vIEEEFloat, y::Double{<:vIEEEFloat})
  s = (x - y.hi)
  v = (s - x)
  Double(s, (((x - (s - v)) - (y.hi + v)) - y.lo))
end

@inline function dsub2(x::Double{<:vIEEEFloat}, y::vIEEEFloat)
  s = x.hi - y
  v = s - x.hi
  Double(s, (((x.hi - (s - v)) - (y + v)) + x.lo))
end

@inline function dsub2(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat})
  s = x.hi - y.hi
  v = s - x.hi
  Double(s, ((((x.hi - (s - v)) - (y.hi + v)) + x.lo) - y.lo))
end

@inline function ifelse(
  b::Mask{N},
  x::Double{T1},
  y::Double{T2},
) where {N,T<:Union{Float32,Float64},T1<:Union{T,Vec{N,T}},T2<:Union{T,Vec{N,T}}}
  V = Vec{N,T}
  Double(ifelse(b, V(x.hi), V(y.hi)), ifelse(b, V(x.lo), V(y.lo)))
end

# two-prod-fma
@inline function dmul(x::vIEEEFloat, y::vIEEEFloat, ::True)
  z = (x * y)
  Double(z, vfmsub(x, y, z))
end
@inline function dmul(x::vIEEEFloat, y::vIEEEFloat, ::False)
  hx, lx = splitprec(x)
  hy, ly = splitprec(y)
  @ieee begin
    z = x * y
    Double(z, (((hx * hy - z) + lx * hy + hx * ly) + lx * ly))
  end
end
@inline function dmul(x::Double{<:vIEEEFloat}, y::vIEEEFloat, ::True)
  z = (x.hi * y)
  Double(z, vfmsub(x.hi, y, z) + x.lo * y)
end
@inline function dmul(x::Double{<:vIEEEFloat}, y::vIEEEFloat, ::False)
  hx, lx = splitprec(x.hi)
  hy, ly = splitprec(y)
  @ieee begin
    z = x.hi * y
    Double(z, (hx * hy - z) + lx * hy + hx * ly + lx * ly + x.lo * y)
  end
end
@inline function dmul(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat}, ::True)
  z = x.hi * y.hi
  Double(z, vfmsub(x.hi, y.hi, z) + x.hi * y.lo + x.lo * y.hi)
end
@inline function dmul(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat}, ::False)
  hx, lx = splitprec(x.hi)
  hy, ly = splitprec(y.hi)
  @ieee begin
    z = x.hi * y.hi
    Double(z, (((hx * hy - z) + lx * hy + hx * ly) + lx * ly) + x.hi * y.lo + x.lo * y.hi)
  end
end
@inline dmul(x::vIEEEFloat, y::Double{<:vIEEEFloat}) = dmul(y, x)
@inline dmul(x, y) = dmul(x, y, fma_fast())
# x^2
@inline function dsqu(x::T, ::True) where {T<:vIEEEFloat}
  z = x * x
  Double(z, vfmsub(x, x, z))
end
@inline function dsqu(x::T, ::False) where {T<:vIEEEFloat}
  hx, lx = splitprec(x)
  @ieee begin
    z = x * x
    Double(z, (hx * hx - z) + lx * (hx + hx) + lx * lx)
  end
end
@inline function dsqu(x::Double{T}, ::True) where {T<:vIEEEFloat}
  z = x.hi * x.hi
  Double(z, vfmsub(x.hi, x.hi, z) + (x.hi * (x.lo + x.lo)))
end
@inline function dsqu(x::Double{T}, ::False) where {T<:vIEEEFloat}
  hx, lx = splitprec(x.hi)
  @ieee begin
    z = x.hi * x.hi
    Double(z, (hx * hx - z) + lx * (hx + hx) + lx * lx + x.hi * (x.lo + x.lo))
  end
end
@inline dsqu(x) = dsqu(x, fma_fast())

# sqrt(x)
@inline function dsqrt(x::Double{T}, ::True) where {T<:vIEEEFloat}
  zhi = @fastmath sqrt(x.hi)
  Double(zhi, (x.lo + vfnmadd(zhi, zhi, x.hi)) / (zhi + zhi))
end
@inline function dsqrt(x::Double{T}, ::False) where {T<:vIEEEFloat}
  c = @fastmath sqrt(x.hi)
  u = dsqu(c, False())
  @ieee Double(c, (x.hi - u.hi - u.lo + x.lo) / (c + c))
end
@inline dsqrt(x) = dsqrt(x, fma_fast())

# x/y
@inline function ddiv(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat}, ::True)
  invy = inv(y.hi)
  zhi = (x.hi * invy)
  Double(zhi, ((vfnmadd(zhi, y.hi, x.hi) + vfnmadd(zhi, y.lo, x.lo)) * invy))
end
@inline function ddiv(x::Double{<:vIEEEFloat}, y::Double{<:vIEEEFloat}, ::False)
  @ieee begin
    invy = one(y.hi) / y.hi
    c = x.hi * invy
    u = dmul(c, y.hi, False())
    Double(c, ((((x.hi - u.hi) - u.lo) + x.lo) - c * y.lo) * invy)
  end
end
@inline function ddiv(x::vIEEEFloat, y::vIEEEFloat, ::True)
  ry = inv(y)
  r = (x * ry)
  Double(r, (vfnmadd(r, y, x) * ry))
end
@inline function ddiv(x::vIEEEFloat, y::vIEEEFloat, ::False)
  @ieee begin
    ry = one(y) / y
    r = x * ry
    hx, lx = splitprec(r)
    hy, ly = splitprec(y)
    Double(r, (((-hx * hy + r * y) - lx * hy - hx * ly) - lx * ly) * ry)
  end
end
@inline ddiv(x, y) = ddiv(x, y, fma_fast())
# 1/x
@inline function drec(x::vIEEEFloat, ::True)
  zhi = inv(x)
  Double(zhi, (vfnmadd(zhi, x, one(eltype(x))) * zhi))
end
@inline function drec(x::vIEEEFloat, ::False)
  @ieee begin
    c = one(x) / x
    u = dmul(c, x, False())
    Double(c, (one(eltype(u.hi)) - u.hi - u.lo) * c)
  end
end

@inline function drec(x::Double{<:vIEEEFloat}, ::True)
  zhi = inv(x.hi)
  Double(zhi, ((vfnmadd(zhi, x.hi, one(eltype(x))) - (zhi * x.lo)) * zhi))
end
@inline function drec(x::Double{<:vIEEEFloat}, ::False)
  @ieee begin
    c = inv(x.hi)
    u = dmul(c, x.hi, False())
    Double(c, (one(eltype(u.hi)) - u.hi - u.lo - c * x.lo) * c)
  end
end
@inline drec(x) = drec(x, fma_fast())


text: # Init to get started

ext: jl

max_stars_repo_name: TarrySingh/Julia

max_issues_repo_name: TarrySingh/Julia

max_forks_repo_name: TarrySingh/Julia


text:
OrdinaryDiffEq.initial_η(fpsolver::FPSolver, integrator::DDEIntegrator) = fpsolver.ηold

OrdinaryDiffEq.apply_step!(fpsolver::FPSolver, integrator::DDEIntegrator) = nothing

function DiffEqBase.postamble!(fpsolver::FPSolver, integrator::DDEIntegrator)
  integrator.destats.nfpiter += fpsolver.iter

  if OrdinaryDiffEq.nlsolvefail(fpsolver)
    integrator.destats.nfpconvfail += 1
  end
  integrator.force_stepfail = OrdinaryDiffEq.nlsolvefail(fpsolver) || integrator.force_stepfail

  nothing
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(mesh::Mesh)
  seriescolor --> :black
  legend --> false

  for elm in elements(mesh)
    @series begin
      primary --> false
      elm
    end
  end
end


text:
using StatisticalRethinkingPlots
using Test

@testset "StatisticalRethinkingPlots.jl" begin
    # Write your tests here.
    @testset "coeftab_plot.jl" begin
        include("test_coeftab_plot.jl")
    end
end

#=

Tests:

plot_density_interval       Fig3.3s.jl
grid_prob                   Fig4.4s.jl

=#


text:
using Test
using LinearAlgebra, SparseArrays
using SummationByPartsOperators

test_list = (Mattsson2012(),)

# Test consistency with constant coefficient operators
for source in test_list, acc_order in (2,4,6), T in (Float32,Float64)
    xmin = -one(T)
    xmax = 2*one(T)
    N = 101

    D2 = derivative_operator(source, 2, acc_order, xmin, xmax, N)
    D2var = try
        var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one)
    catch err
        !isa(err, ArgumentError) && throw(err)
        nothing
    end
    D2var == nothing && continue

    println(devnull, D2var)
    println(devnull, D2var.coefficients)
    @test maximum(abs, Matrix(D2) - Matrix(D2var)) < 10000*eps(T)
end

# Compare mul! with β=0 and mul! without β.
for T in (Float32, Float64), acc_order in (2,4,6)
    xmin = zero(T)
    xmax = 5*one(T)
    N = 51
    source = Mattsson2012()

    D2var_serial = var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one, Val{:serial}())
    D2var_threads= var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one, Val{:threads}())
    D2var_full   = Matrix(D2var_serial)
    D2var_sparse = sparse(D2var_serial)

    x = grid(D2var_serial)
    u = x.^5
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))

    mul!(dest1, D2var_serial, u, one(T), zero(T))
    mul!(dest2, D2var_serial, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    mul!(dest1, D2var_threads, u, one(T), zero(T))
    mul!(dest2, D2var_threads, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    mul!(dest2, D2var_serial, u)
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    mul!(dest2, D2var_full, u)
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
end


text:
# Types used in Arena.jl

const lmpoint = HyperSphere(Point(0.,0,0), 0.05)
const greenMat = MeshPhongMaterial(color=RGBA(0, 1, 0, 0.5))
const redMat = MeshPhongMaterial(color=RGBA(0, 1, 0, 0.5))

"""
	$(TYPEDEF)
Type for 2d visualization
"""
struct BotVis2
    vis::Visualizer
    # poses::Dict{Symbol, NTuple{3,Float64}}
    # landmarks::Dict{Symbol, NTuple{3,Float64}}
    cachevars::Dict{Symbol, Tuple{Symbol, Vector{Bool}, Vector{Float64}}}
end

mutable struct SubmapColorCheat
  colors
  SubmapColorCheat(;colors::Vector=
    [ RGB(0.651, 0.808, 0.890),
      RGB(0.122, 0.471, 0.706),
      RGB(0.698, 0.875, 0.541),
      RGB(0.200, 0.627, 0.172),
      RGB(0.984, 0.604, 0.600),
      RGB(0.890, 0.102, 0.110),
      RGB(0.992, 0.749, 0.043)]
  ) = new(colors)
end

mutable struct ArcPointsRangeSolve <: Function
  x1::Vector{Float64}
  x2::Vector{Float64}
  x3::Vector{Float64}
  r::Float64
  center::Vector{Float64}
  angle::Float64
  axis::Vector{Float64}
  ArcPointsRangeSolve(x1::Vector{Float64}, x2::Vector{Float64}, r::Float64) = new(x1,x2,zeros(0),r, zeros(2), 0.0, zeros(3))
  ArcPointsRangeSolve(x1::Vector{Float64}, x2::Vector{Float64}, x3::Vector{Float64}, r::Float64) = new(x1,x2,x3,r, zeros(3), 0.0, zeros(3))
end


text:
"""
    DecisionRef <: AbstractVariableRef

Holds a reference to the model, the stage the decision is taken in, and the corresponding MOI.VariableIndex.
"""
struct DecisionRef <: JuMP.AbstractVariableRef
    model::JuMP.Model
    index::MOI.VariableIndex
end
is_decision_type(::Type{DecisionRef}) = true

# Getters (model) #
# ========================== #
function get_decisions(model::JuMP.Model, s::Integer = 1)::Decisions
    !haskey(model.ext, :decisions) && return IgnoreDecisions()
    N = length(model.ext[:decisions])
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    return model.ext[:decisions][s]
end

function all_decisions(model::JuMP.Model, s::Integer = 1)
    decisions = get_decisions(model, s)::Decisions
    return all_decisions(decisions)
end

function all_known_decisions(model::JuMP.Model, s::Integer = 2)
    decisions = get_decisions(model, s)::Decisions
    return all_known_decisions(decisions)
end
"""
    all_decision_variables(model::JuMP.Model)

Returns a stage-wise list of all decisions currently in the `model`. The decisions are
ordered by creation time.
"""
function all_decision_variables(model::JuMP.Model)
    haskey(model.ext, :decisions) || error("No decisions in model.")
    N = length(model.ext[:decisions])
    return ntuple(Val{N}()) do stage
        return all_decision_variables(model, stage)
    end
end
"""
    all_decision_variables(model::JuMP.Model, stage::Integer)

Returns a list of all decisions currently in the `model` at stage `stage`. The decisions are
ordered by creation time.
"""
function all_decision_variables(model::JuMP.Model, stage::Integer)
    haskey(model.ext, :decisions) || error("No decisions in model.")
    decisions = get_decisions(model, stage)::Decisions
    return map(all_decisions(decisions)) do index
        DecisionRef(model, index)
    end
end
"""
    all_known_decision_variables(model::JuMP.Model)

Returns a stage-wise list of all known decisions currently in the `model`. The decisions are
ordered by creation time.
"""
function all_known_decision_variables(model::JuMP.Model)
    haskey(model.ext, :decisions) || error("No decisions in model.")
    N = length(model.ext[:decisions])
    return ntuple(Val{N}()) do s
        return all_known_decision_variables(model, s)
    end
end
"""
    all_known_decision_variables(model::JuMP.Model, stage::Integer)

Returns a stage-wise list of all known decisions currently in the `model` at stage `stage`. The decisions are
ordered by creation time.
"""
function all_known_decision_variables(model::JuMP.Model, stage::Integer)
    haskey(model.ext, :decisions) || error("No decisions in model.")
    decisions = get_decisions(model, stage)::Decisions
    return map(all_known_decisions(decisions)) do index
        DecisionRef(model, index)
    end
end
"""
    num_decisions(model::JuMP.Model, stage::Integer = 1)

Return the number of decisions in `model` at stage `stage`. Defaults to the first stage.
"""
function num_decisions(model::JuMP.Model, stage::Integer = 1)
    decisions = get_decisions(model, stage)::Decisions
    return num_decisions(decisions)
end
"""
    num_known_decisions(model::JuMP.Model, stage::Integer = 2)

Return the number of known decisions in `model` at stage `stage`. Defaults to the second stage.
"""
function num_known_decisions(model::JuMP.Model, stage::Integer = 2)
    stage > 1 || error("No decisions can be known in the first stage.")
    decisions = get_decisions(model, stage - 1)::Decisions
    return num_known_decisions(decisions)
end

# Getters (refs) #
# ========================== #
function stage(dref::DecisionRef)
    haskey(dref.model.ext, :stage_map) || error("No decisions in model.")
    return dref.model.ext[:stage_map][index(dref)]
end

function get_decisions(dref::DecisionRef)::Decisions
    s = stage(dref)
    s == 0 && return IgnoreDecisions()
    return get_decisions(dref.model, s)
end
"""
    decision(dref::DecisionRef)

Return the internal `Decision` associated with `dref`.
"""
function decision(dref::DecisionRef)
    decisions = get_decisions(dref)::Decisions
    return decision(decisions, index(dref))
end
"""
    state(dref::DecisionRef)

Return the `DecisionState` of `dref`.
"""
function state(dref::DecisionRef)
    return decision(dref).state
end

# Setters #
# ========================== #
function take_decisions!(model::JuMP.Model, drefs::Vector{DecisionRef}, vals::AbstractVector)
    # Check that all given decisions are in model
    map(dref -> check_belongs_to_model(dref, model), drefs)
    # Check decision length
    length(drefs) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(drefs)).")
    # Update decisions
    for (dref, val) in zip(drefs, vals)
        d = decision(dref)
        # Update state
        d.state = Taken
        # Update value
        d.value = val
        # Update
        update_decision_state!(dref, Taken)
    end
    return nothing
end

function untake_decisions!(model::JuMP.Model, drefs::Vector{DecisionRef})
    # Check that all given decisions are in model
    map(dref -> check_belongs_to_model(dref, model), drefs)
    # Update decisions
    need_update = false
    for dref in drefs
        d = decision(dref)
        if state(d) == Taken
            need_update |= true
            # Update state
            d.state = NotTaken
            update_decision_state!(dref, NotTaken)
        end
    end
    return nothing
end

# JuMP variable interface #
# ========================== #
function MOI.get(model::JuMP.Model, attr::MOI.AbstractVariableAttribute,
                 dref::DecisionRef)
    check_belongs_to_model(dref, model)
    if MOI.is_set_by_optimize(attr)
        return JuMP._moi_get_result(backend(model), attr, index(dref))
    else
        return MOI.get(backend(model), attr, index(dref))
    end
end

function MOI.set(model::Model, attr::MOI.AbstractVariableAttribute,
                 dref::DecisionRef, value)
    check_belongs_to_model(dref, model)
    MOI.set(backend(model), attr, index(dref), value)
    return nothing
end

JuMP.name(dref::DecisionRef) = MOI.get(owner_model(dref), MOI.VariableName(), dref)::String

function JuMP.set_name(dref::DecisionRef, name::String)
    return MOI.set(owner_model(dref), MOI.VariableName(), dref, name)
end

function decision_by_name(model::Model, name::String)
    index = MOI.get(backend(model), MOI.VariableIndex, name)
    if index isa Nothing
        return nothing
    else
        return DecisionRef(model, index)
    end
end

JuMP.index(dref::DecisionRef) = dref.index

function JuMP.value(dref::DecisionRef; result::Int = 1)::Float64
    if state(dref) == NotTaken
        return MOI.get(owner_model(dref), MOI.VariablePrimal(result), dref)
    end
    # If decision has been taken or is known
    # the value can be fetched directly
    return decision(dref).value
end

function JuMP.is_fixed(dref::DecisionRef)
    if state(dref) == NotTaken
        return false
    end
    return true
end

function JuMP.unfix(dref::DecisionRef)
    if state(dref) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    d = decision(dref)
    # Update state
    d.state = NotTaken
    # Update objective and constraints
    update_decision_state!(dref, NotTaken)
    return nothing
end

"""
    fix(dref::DecisionRef, val::Number)

Fix the decision associated with `dref` to `val`.

See also [`unfix`](@ref).
"""
function JuMP.fix(dref::DecisionRef, val::Number)
    d = decision(dref)
    if state(dref) == NotTaken
        # Update state
        d.state = Taken
        # Update value
        d.value = val
    else
        # Just update value
        d.value = val
    end
    # Modify decision state
    update_decision_state!(dref, d.state)
    return nothing
end

function JuMP.delete(model::JuMP.Model, dref::DecisionRef)
    if model !== owner_model(dref)
        error("The decision you are trying to delete does not " *
              "belong to the model.")
    end
    # First delete any SingleDecision constraints
    for S in [MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.EqualTo{Float64}, MOI.ZeroOne, MOI.Integer]
        ci = CI{SingleDecision,S}(index(dref).value)
        if MOI.is_valid(backend(model), ci)
            MOI.delete(backend(model), ci)
        end
    end
    # Remove SingleDecisionSet constraint
    ci = CI{MOI.SingleVariable,SingleDecisionSet{Float64}}(index(dref).value)
    MOI.delete(backend(model), ci)
    # Delete the variable corresponding to the decision
    MOI.delete(backend(model), index(dref))
    # Remove the decision
    remove_decision!(get_decisions(dref), index(dref))
    return nothing
end

function JuMP.delete(model::JuMP.Model, drefs::Vector{DecisionRef})
    isempty(drefs) && return nothing
    if any(model !== owner_model(dref) for dref in drefs)
        error("A decision you are trying to delete does not " *
              "belong to the model.")
    end
    # First delete any SingleDecision constraints
    for dref in drefs
        for S in [MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.EqualTo{Float64}, MOI.ZeroOne, MOI.Integer]
            ci = CI{SingleDecision,S}(index(dref).value)
            if MOI.is_valid(backend(model), ci)
                MOI.delete(backend(model), ci)
            end
        end
    end
    # Delete the variables corresponding to the decision
    MOI.delete(backend(model), index.(drefs))
    # Remove any MultipleDecisionSet constraint
    for ci in MOI.get(backend(model), MOI.ListOfConstraintIndices{MOI.VectorOfVariables, MultipleDecisionSet{Float64}}())
        f = MOI.get(backend(model), MOI.ConstraintFunction(), ci)::MOI.VectorOfVariables
        if all(f.variables .== index.(drefs))
            # This is the constraint
            MOI.delete(backend(model), ci)
            break
        end
    end
    # Remove any SingleDecisionSet constraints
    for dref in drefs
        ci = CI{MOI.SingleVariable,SingleDecisionSet{Float64}}(index(dref).value)
        if MOI.is_valid(backend(model), ci)
            MOI.delete(backend(model), ci)
        end
    end
    # Remove the decisions
    map(dref -> remove_decision!(get_decisions(dref), index(dref)), drefs)
    return nothing
end

JuMP.owner_model(dref::DecisionRef) = dref.model

struct DecisionNotOwned <: Exception
    dref::DecisionRef
end

function JuMP.check_belongs_to_model(dref::DecisionRef, model::AbstractModel)
    if owner_model(dref) !== model
        throw(DecisionNotOwned(dref))
    end
end

function JuMP.is_valid(model::Model, dref::DecisionRef)
    return model === owner_model(dref)
end

function JuMP.has_lower_bound(dref::DecisionRef)
    if state(dref) == Known
        return false
    end
    ci = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), ci)
end
function JuMP.LowerBoundRef(dref::DecisionRef)
    moi_lb =  MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}
    ci = moi_lb(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_lb, ScalarShape}(owner_model(dref),
                                                          ci,
                                                          ScalarShape())
end
function JuMP.set_lower_bound(dref::DecisionRef, lower::Number)
    new_set = MOI.GreaterThan(convert(Float64, lower))
    if has_lower_bound(dref)
        ci = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dref).value)
        MOI.set(backend(owner_model(dref)), MOI.ConstraintSet(), ci, new_set)
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), new_set)
    end
    return nothing
end
function JuMP.delete_lower_bound(dref::DecisionRef)
    JuMP.delete(owner_model(dref), LowerBoundRef(dref))
end
function JuMP.lower_bound(dref::DecisionRef)
    if !has_lower_bound(dref)
        error("Decision $(dref) does not have a lower bound.")
    end
    cset = MOI.get(owner_model(dref), MOI.ConstraintSet(),
                   LowerBoundRef(dref))::MOI.GreaterThan{Float64}
    return cset.lower
end

function JuMP.has_upper_bound(dref::DecisionRef)
    if state(dref) == Known
        return false
    end
    ci = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), ci)
end
function JuMP.UpperBoundRef(dref::DecisionRef)
    moi_ub =  MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}
    ci = moi_ub(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_ub, ScalarShape}(owner_model(dref),
                                                          ci,
                                                          ScalarShape())
end
function JuMP.set_upper_bound(dref::DecisionRef, lower::Number)
    new_set = MOI.LessThan(convert(Float64, lower))
    if has_upper_bound(dref)
        ci = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dref).value)
        MOI.set(backend(owner_model(dref)), MOI.ConstraintSet(), ci, new_set)
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), new_set)
    end
    return nothing
end
function JuMP.delete_upper_bound(dref::DecisionRef)
    JuMP.delete(owner_model(dref), UpperBoundRef(dref))
end
function JuMP.upper_bound(dref::DecisionRef)
    if !has_upper_bound(dref)
        error("Decision $(dref) does not have a upper bound.")
    end
    cset = MOI.get(owner_model(dref), MOI.ConstraintSet(),
                   UpperBoundRef(dref))::MOI.LessThan{Float64}
    return cset.upper
end

function JuMP.is_integer(dref::DecisionRef)
    if state(dref) == Known
        return false
    end
    ci = MOI.ConstraintIndex{SingleDecision, MOI.Integer}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), ci)
end
function JuMP.set_integer(dref::DecisionRef)
    if is_integer(dref)
        return nothing
    elseif is_binary(dref)
        error("Cannot set the decision $(dref) to integer as it is already binary.")
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), MOI.Integer())
    end
end
function JuMP.unset_integer(dref::DecisionRef)
    JuMP.delete(owner_model(dref), IntegerRef(dref))
    return nothing
end
function JuMP.IntegerRef(dref::DecisionRef)
    moi_int =  MOI.ConstraintIndex{SingleDecision, MOI.Integer}
    ci = moi_int(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_int, ScalarShape}(owner_model(dref),
                                                           ci,
                                                           ScalarShape())
end

function JuMP.is_binary(dref::DecisionRef)
    if state(dref) == Known
        return false
    end
    ci = MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), ci)
end
function JuMP.set_binary(dref::DecisionRef)
    if is_binary(dref)
        return nothing
    elseif is_integer(dref)
        error("Cannot set the decision $(dref) to binary as it is already integer.")
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), MOI.ZeroOne())
    end
end
function JuMP.unset_binary(dref::DecisionRef)
    JuMP.delete(owner_model(dref), BinaryRef(dref))
    return nothing
end
function JuMP.BinaryRef(dref::DecisionRef)
    moi_bin =  MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}
    ci = moi_bin(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_bin, ScalarShape}(owner_model(dref),
                                                           ci,
                                                           ScalarShape())
end

function JuMP.start_value(dref::DecisionRef)
    return MOI.get(owner_model(dref), MOI.VariablePrimalStart(), dref)
end
function JuMP.set_start_value(dref::DecisionRef, value::Number)
    MOI.set(owner_model(dref), MOI.VariablePrimalStart(), dref, Float64(value))
end

function Base.hash(dref::DecisionRef, h::UInt)
    return hash(objectid(owner_model(dref)), hash(dref.index, h))
end

JuMP.isequal_canonical(d::DecisionRef, other::DecisionRef) = isequal(d, other)
function Base.isequal(dref::DecisionRef, other::DecisionRef)
    return owner_model(dref) === owner_model(other) && dref.index == other.index
end

Base.iszero(::DecisionRef) = false
Base.copy(dref::DecisionRef) = DecisionRef(dref.model, dref.index)
Base.broadcastable(dref::DecisionRef) = Ref(dref)

# JuMP copy interface #
# ========================== #
function JuMP.copy_extension_data(decisions::NTuple{N,Decisions}, dest::Model, src::Model) where N
    new_decisions = ntuple(Val{N}()) do _
        Decisions()
    end
    for s in 1:N
        for dref in all_decision_variables(src, s)
            set_decision!(new_decisions[s], index(dref), decision(dref))
        end
        # TODO
        for kref in all_known_decision_variables(src, s)
            set_known_decision!(new_decisions[s], index(kref), decision(kref))
        end
    end
    return new_decisions
end

function Base.getindex(reference_map::JuMP.ReferenceMap, dref::DecisionRef)
    return DecisionRef(reference_map.model,
                       reference_map.index_map[index(dref)])
end


text:
# # BESO example: HalfMBB Beam
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`beso.ipynb`](@__NBVIEWER_ROOT_URL__/examples/beso.ipynb)
#-
# ## Commented Program
#
# What follows is a program spliced with comments.
#md # The full program, without comments, can be found in the next [section](@ref beso-plain-program).

using TopOpt

# ### Define the problem
E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0; # downward force

nels = (160, 40)
problem = HalfMBB(Val{:Linear}, nels, (1.0, 1.0), E, v, f)

# ### Define the FEA Solver and penalty functions
solver = FEASolver(Displacement, Direct, problem, xmin = 0.01, penalty = TopOpt.PowerPenalty(3.0))

# ### Define the compliance objective function and volume fraction constraint
comp = Compliance(problem, solver)
volfrac = Volume(problem, solver)
sensfilter = SensFilter(solver, rmin = 4.0)
beso = BESO(comp, volfrac, 0.5, sensfilter)

# ### Run optimization
x0 = ones(length(solver.vars))
result = beso(x0)

# ### (Optional) Visualize the result using Makie.jl
# Need to run `using Pkg; Pkg.add(Makie)` first
# ```julia
# using TopOpt.TopOptProblems.Visualization: visualize
# fig = visualize(problem; topology = result.topology)
# Makie.display(fig)
# ```

#md # ## [Plain Program](@id beso-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [beso.jl](beso.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
#========================================================================================#
#	Laboratory 200
#
# Welcome to course 200: Modelling Dynamical Systems!
#
# Author: Niall Palfreyman, 24/04/2022
#========================================================================================#
[
	Activity(
		"""
		???
		""",
		"???",
		x -> true
	),
]


text:
# Functions separated from type definitions for compiler reasons

#### Abstract Interface Functions ####

"""
Returns the constant penalty parameter value. Only required if `is_variable_dependent` returns false.
"""
function get_penalty_value(r::AbstractPenaltyParameter)::Float64
    throw(UnimplementedError("get_penalty_value is unimplemented for penalty parameter of type $(typeof(r))."))
end

"""
Returns the penalty value for the consensus variable associated with `xhid`.
"""
function get_penalty_value(r::AbstractPenaltyParameter,
                           xhid::XhatID
                           )::Float64
    throw(UnimplementedError("get_penalty_value is unimplemented for penalty parameter of type $(typeof(r))."))
end

"""
Returns true if the penalty parameter value is dependent on the initial solutions of the subproblems.
"""
function is_initial_value_dependent(p::Type{P})::Bool where P <: AbstractPenaltyParameter
    throw(UnimplementedError("is_initial_value_dependent is unimplemented for penalty parameter of type $(p)."))
end

"""
Returns true if the penalty parameter value is dependent on any data or values in the subproblems (e.g., the coefficients in the objective function).
"""
function is_subproblem_dependent(p::Type{P})::Bool where P <: AbstractPenaltyParameter
    throw(UnimplementedError("is_subproblem_dependent is unimplemented for penalty parameter of type $(p)."))
end

"""
Returns true if the penalty parameter value may differ for different consensus variables.
"""
function is_variable_dependent(p::Type{P})::Bool where P <: AbstractPenaltyParameter
    throw(UnimplementedError("is_variable_dependent is unimplemented for penalty parameter of type $(p)."))
end

"""
Returns a mapping of consensus variable ids to penalty parameter values.

Only required if `is_variable_dependent` returns false.
"""
function penalty_map(r::AbstractPenaltyParameter)::Dict{XhatID,Float64}
    throw(UnimplementedError("penalty_map is unimplemented for penalty parameter of type $(typeof(r))."))
end

"""
Performs any computations for the penalty parameter based on the initial solutions of the subproblems.

**Arguments**

*`r::AbstractPenaltyParameter` : penalty parameter struct (replace with appropriate type)
*`phd::PHData` : PH data structure used for obtaining any required variable values. See help on `PHData` for details on available functions.
"""
function process_penalty_initial_value(r::AbstractPenaltyParameter,
                                       phd::PHData,
                                       )::Nothing
    throw(UnimplementedError("process_penalty_initial_value is unimplemented for penalty parameter of type $(typeof(r))."))
end

"""
Performs any computations for the penalty parameter based data or values from the subproblems.

This function is called *before* the initial solution of subproblems. Any accessing of variable values in this function may result in undefined behavior.

**Arguments**

* `r::AbstractPenaltyParameter` : penalty parameter struct (replace with appropriate type)
* `phd::PHData` : PH data structure used for obtaining any values or information not provided by `scid` and `subproblem_dict`
* `scid::ScenarioID` : scenario id from which `subproblem_dict` comes
* `subproblem_dict::Dict{VariableID,Float64}` : Mapping specifying a value needed from a subproblem used in the computation of the penalty parameter.
"""
function process_penalty_subproblem(r::AbstractPenaltyParameter,
                                    phd::PHData,
                                    scid::ScenarioID,
                                    subproblem_dict::Dict{VariableID,Float64}
                                    )::Nothing
    throw(UnimplementedError("process_penalty_subproblem is unimplemented for penalty parameter of type $(typeof(r))."))
end

"""
Resets any mapping of consensus variable ids to penalty parameter values which `r` may have. This avoids potential bugs caused by the same instance of a penalty parameter being used for several calls to `solve`.
"""
function reset_penalty_map(r::AbstractPenaltyParameter)::Nothing
    return
end

#### Concrete Implementations ####

## Proportional Penalty Parameter ##

function ProportionalPenaltyParameter(constant::Real)
    return ProportionalPenaltyParameter(constant, Dict{XhatID,Float64}())
end

function get_penalty_value(r::ProportionalPenaltyParameter,
                           xhid::XhatID,
                           )::Float64
    return r.penalties[xhid]
end

function penalty_map(r::ProportionalPenaltyParameter)::Dict{XhatID,Float64}
    return r.penalties
end

function process_penalty_subproblem(r::ProportionalPenaltyParameter,
                                    phd::PHData,
                                    scid::ScenarioID,
                                    penalties::Dict{VariableID,Float64}
                                    )::Nothing

    for (vid, penalty) in pairs(penalties)

        xhid = convert_to_xhat_id(phd, vid)
        r_value = penalty == 0.0 ? r.constant : abs(penalty) * r.constant

        if haskey(r.penalties, xhid)
            if !isapprox(r.penalties[xhid], r_value)
                error("Penalty parameter must match across scenarios. Got $(r.penalties[xhid]) instead of $(penalty) for variable $(name(phd, xhid)) when processing penalty parameter for scenario $(scid). This error likely means the coefficients in the objective function do not match in the scenario subproblems.")
            end
        else
            r.penalties[xhid] = r_value
        end

    end

    return
end

function is_initial_value_dependent(::Type{ProportionalPenaltyParameter})::Bool
    return false
end

function is_subproblem_dependent(::Type{ProportionalPenaltyParameter})::Bool
    return true
end

function is_variable_dependent(::Type{ProportionalPenaltyParameter})::Bool
    return true
end

function reset_penalty_map(r::ProportionalPenaltyParameter)::Nothing
    for k in keys(r.penalties)
        delete!(r.penalties, k)
    end
    return
end

## Scalar Penalty Parameter ##

function get_penalty_value(r::ScalarPenaltyParameter)::Float64
    return r.value
end

function get_penalty_value(r::ScalarPenaltyParameter,
                           xhid::XhatID,
                           )::Float64
    return get_penalty_value(r)
end

function is_initial_value_dependent(::Type{ScalarPenaltyParameter})::Bool
    return false
end

function is_subproblem_dependent(::Type{ScalarPenaltyParameter})::Bool
    return false
end

function is_variable_dependent(::Type{ScalarPenaltyParameter})::Bool
    return false
end

## SEP Penalty Parameter ##

function SEPPenaltyParameter(default::Float64=1.0)
    return SEPPenaltyParameter(default, Dict{XhatID,Float64}())
end

function get_penalty_value(r::SEPPenaltyParameter,
                           xhid::XhatID,
                           )::Float64
    return r.penalties[xhid]
end

function penalty_map(r::SEPPenaltyParameter)::Dict{XhatID,Float64}
    return r.penalties
end

function is_initial_value_dependent(::Type{SEPPenaltyParameter})::Bool
    return true
end

function process_penalty_initial_value(r::SEPPenaltyParameter,
                                       phd::PHData,
                                       )::Nothing

    for (xhid, xhat) in pairs(consensus_variables(phd))

        if is_integer(xhat)

            xmin = typemax(Int)
            xmax = typemin(Int)

            for vid in variables(xhat)

                xs = branch_value(phd, vid)

                if xs < xmin
                    xmin = xs
                end

                if xs > xmax
                    xmax = xs
                end

            end

            denom = xmax - xmin + 1

        else

            denom = 0.0

            for vid in variables(xhat)
                p = probability(phd, scenario(vid))
                xs = branch_value(phd, vid)
                denom += p * abs(xs - value(xhat))
            end

            denom = max(denom, 1.0)

        end

        obj_coeff = r.penalties[xhid]
        r.penalties[xhid] = obj_coeff == 0.0 ? r.default : abs(obj_coeff)/denom

    end

    return
end

function is_subproblem_dependent(::Type{SEPPenaltyParameter})::Bool
    return true
end

function process_penalty_subproblem(r::SEPPenaltyParameter,
                                    phd::PHData,
                                    scid::ScenarioID,
                                    penalties::Dict{VariableID,Float64}
                                    )::Nothing

    for (vid, penalty) in pairs(penalties)

        xhid = convert_to_xhat_id(phd, vid)

        if haskey(r.penalties, xhid)
            if !isapprox(r.penalties[xhid], penalty)
                error("Penalty parameter must match across scenarios. Got $(r.penalties[xhid]) instead of $(penalty) for variable $(name(phd, xhid)) when processing penalty parameter for scenario $(scid). This error likely means the coefficients in the objective function do not match in the scenario subproblems.")
            end
        else
            r.penalties[xhid] = penalty
        end

    end

    return
end

function is_variable_dependent(::Type{SEPPenaltyParameter})::Bool
    return true
end

function reset_penalty_map(r::SEPPenaltyParameter)::Nothing
    for k in keys(r.penalties)
        delete!(r.penalties, k)
    end
    return
end


text:
module BasisFunctions

using Common, StaticArrays, OffsetArrays

export Gaussian, cartesianToSphericalMxs

struct Gaussian
    ζ::Float
    R::Coordinate
end

function positiveDiagonal(l::Int8, m::Int8)
    δl0 = l == 0 ? 1 : 0
    return sqrt(2 ^ δl0 * (2*l + 1) / (2 * l + 2))
end

function cartesianToSphericalMxs(lmax::Int)
    nang = (lmax+1)^2
    ncart = (lmax+1)*(lmax+2)*(lmax+3)÷6
    M = OffsetArray(Array([Matrix(zeros(Float64, ((l+1)*(l+2)÷2, 2*l+1))) for l=0:lmax]), -1)
    #  M = zeros(Float, (ncart, nang))
    for l=0:lmax #  Over common angular momentum quantum number
        iang = 0
        println()
        for m=-l:l  #  Spherical quantum number m
            iang += 1
            vm = m < 0 ? 0.5 : 0
            am = abs(m)
            vmax = Int(floor(am/2-vm))
            mfact = m == 0 ? 0.5 : 1
            N = 1 / (2^am * factorial(l)) * sqrt(2 * mfact * (factorial(l+am) * factorial(l-am)))
            #  icart = l ^ 2 #  Take into account the cartesians of all the previous ls
            icart = 0
            for k=0:l #  Cartesian quantum number k
                for j=0:l-k #  Cartesian quantum number j
                    i=l-k-j #  Cartesian quantum number i
                    #  if m == 0
                        #  println("x,y,z:",i,j,k)
                    #  end
                    icart += 1
                    for t=0:(l-am)÷2
                        kk = l - 2 * t - am
                        if kk ≠ k continue end
                        for u=0:t
                            for vv=0:vmax
                                v=vv+vm
                                ii = 2 * t + am - 2 * (u + v)
                                jj = 2 * (u + v)
                                if ii == i && jj == j && kk == k
                                    M[l][icart,iang] += (-1)^(t+vv) * 0.25 ^ t * N * binomial(l,t) * binomial(l-t,am+t) * binomial(t,u) * binomial(am,Int(2*v))
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return M
end

end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")
include("../LF_common/prec_mkl_ilu0.jl")
include("../LF_common/apply_prec_ilu0.jl")
include("build_nabla2.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  # Directly build Laplacian and ILU preconditioner from MKL
  @time nabla2 = build_nabla2(LF)
  @time prec_nabla2 = prec_mkl_ilu0( nabla2 )

  # Using Kronecker product and ILU0 preconditioner from MKL
  @time ∇2 = get_Laplacian3d_kron(LF)
  @time prec_∇2 = prec_mkl_ilu0( ∇2 )

  # Test the preconditioner
  srand(123)
  v = rand(Npoints)

  Pv_1 = apply_prec_ilu0( prec_nabla2, v )
  Pv_2 = apply_prec_ilu0( prec_∇2, v )

  println( sum(Pv_1 - Pv_2) )

end

test_main( [50, 50, 50])


text:
using MXNet

#--------------------------------------------------------------------------------
# define MLP
# the following two ways are equivalent

#-- Option 1: explicit composition
# data = mx.Variable(:data)
# fc1  = mx.FullyConnected(data = data, name=:fc1, num_hidden=128)
# act1 = mx.Activation(data = fc1, name=:relu1, act_type=:relu)
# fc2  = mx.FullyConnected(data = act1, name=:fc2, num_hidden=64)
# act2 = mx.Activation(data = fc2, name=:relu2, act_type=:relu)
# fc3  = mx.FullyConnected(data = act2, name=:fc3, num_hidden=10)
# mlp  = mx.SoftmaxOutput(data = fc3, name=:softmax)

#-- Option 2: using the mx.chain macro
# mlp = @mx.chain mx.Variable(:data)             =>
#   mx.FullyConnected(name=:fc1, num_hidden=128) =>
#   mx.Activation(name=:relu1, act_type=:relu)   =>
#   mx.FullyConnected(name=:fc2, num_hidden=64)  =>
#   mx.Activation(name=:relu2, act_type=:relu)   =>
#   mx.FullyConnected(name=:fc3, num_hidden=10)  =>
#   mx.SoftmaxOutput(name=:softmax)

#-- Option 3: using nn-factory
mlp = @mx.chain mx.Variable(:data) =>
  mx.MLP([128, 64, 10])            =>
  mx.SoftmaxOutput(name=:softmax)

# data provider
batch_size = 100
include("mnist-data.jl")
train_provider, eval_provider = get_mnist_providers(batch_size)

# setup model
model = mx.FeedForward(mlp, context=mx.cpu())

# optimizer
optimizer = mx.SGD(lr=0.1, momentum=0.9, weight_decay=0.00001)

# fit parameters
mx.fit(model, optimizer, train_provider, eval_data=eval_provider, n_epoch=20)

#--------------------------------------------------------------------------------
# Optional, demonstration of the predict API
probs = mx.predict(model, eval_provider)

# collect all labels from eval data
labels = Array[]
for batch in eval_provider
  push!(labels, copy(mx.get(eval_provider, batch, :softmax_label)))
end
labels = cat(1, labels...)

# Now we use compute the accuracy
correct = 0
for i = 1:length(labels)
  # labels are 0...9
  if indmax(probs[:,i]) == labels[i]+1
    correct += 1
  end
end
println(mx.format("Accuracy on eval set: {1:.2f}%", 100correct/length(labels)))


text:
@testset "6.zig-zag-conversion.jl" begin
    @test zig_zag_convert("PAYPALISHIRING", 4) == "PINALSIGYAHRPI"
    @test zig_zag_convert("PAYPALISHIRING", 3) == "PAHNAPLSIIGYIR"
end


text:
# conclusion shortcircuit appears faster
const CHAR0 = Char(0)

function code_unit_if(str, pos)
    res = CHAR0
    if sizeof(str) > pos
        @inbounds res = codeunit(str, pos)
    end
    return res
end

function code_unit_short_circuit(str, pos)
    @inbounds return sizeof(str) >  pos ? CHAR0 :  res = codeunit(str, pos)
end

using DataBench, BenchmarkTools, Plots
gr()
a = DataBench.gen_string_vec_var_len(10_000_000, 8);

x = @belapsed code_unit_if.(a, 8)
y = @belapsed code_unit_short_circuit.(a, 8)
bar([x,y])

x = @belapsed code_unit_if.(a, 8)
y = @belapsed code_unit_short_circuit.(a, 8)
bar([x,y])


text:
@testset "persistence helpers" begin
    n = 100
    nrepeats = 10

    @testset "sampledensity" begin
        using IntrinsicPersistence.sampledensity

        for dim in 2:nrepeats+1
            # Points
            pts = rand(dim, n)
            M   = pairwise(Euclidean(), pts)
            M2  = copy(M)
            for i in 1:n
                M2[i, i] = Inf
            end
            @test @inferred(sampledensity(M)) == maximum(minimum(M2, 1))

            # Graphs
            G = Graph(n, 2n)
            @test @inferred(sampledensity(G)) == 1
            A = sprand(n, n, 0.1)
            A += A'
            G = SimpleWeightedGraph(A)
            @test @inferred(sampledensity(G)) == maximum(A)
        end
    end

    @testset "edgetoorder" begin
        using IntrinsicPersistence.edgetoorder

        for _ in 1:nrepeats
            idxs = Tuple{Int, Int}[]
            A = spzeros(n, n)
            for i in 1:n
                j, k = 0, 0
                while j == k || A[j, k] ≠ 0
                    j, k = rand(1:n, 2)
                end
                A[j, k] = i
                A[k, j] = i
                push!(idxs, (j, k))
            end

            # Make sure this appears as a unit test.
            @test @inferred(edgetoorder(SimpleWeightedGraph(A))) ≠ nothing
            ord = edgetoorder(SimpleWeightedGraph(A))
            for i in 1:n
                @test ord[idxs[i]...] == i
            end
        end
    end

    @testset "spannergraph" begin
        using IntrinsicPersistence.spannergraph

        # Metrics that don't need a parameter.
        metrics = [Distances.Chebyshev(),
                   Distances.Cityblock(),
                   Distances.Euclidean(),
                   Distances.Hamming(),
                   Distances.HellingerDist(),
                   Distances.Jaccard(),
                   Distances.MeanAbsDeviation(),
                   Distances.RMSDeviation()]

        for dim in 2:nrepeats+1, m in metrics
            pts1 = rand(dim, n)
            M1   = pairwise(m, pts1)
            pts2 = [rand(SVector{dim, Float64}) for _ in 1:n]
            M2   = pairwise(m, reinterpret(Float64, pts2, (dim, n)))

            @test @inferred(spannergraph(pts1, m, 2maximum(M1))).weights ≈ M1
            @test @inferred(spannergraph(pts2, m, 2maximum(M2))).weights ≈ M2

            r1 = 2 * sampledensity(M1)
            r2 = 2 * sampledensity(M2)
            for i in 1:n, j in 1:n
                if M1[i, j] ≥ r1
                    M1[i, j] = 0
                    M1[j, i] = 0
                end
                if M2[i, j] ≥ r2
                    M2[i, j] = 0
                    M2[j, i] = 0
                end
            end
            @test @inferred(spannergraph(pts1, m)).weights ≈ M1
            @test @inferred(spannergraph(pts2, m)).weights ≈ M2
        end
    end
end

@testset "persistence" begin
    @testset "graphs" begin

        function squarelattice(n, m)
            G = Graph(n*m)
            L = reshape(1:n*m, n, m)
            for j in 1:m, i in 1:n
                i < n && add_edge!(G, L[i,j], L[i+1,j])
                j < m && add_edge!(G, L[i,j], L[i,j+1])
                i < n && j < m && add_edge!(G, L[i,j], L[i+1,j+1])
            end
            G
        end

        function glue!(G::Graph, r1, r2)
            length(r1) == length(r2) || throw("Edges should have the same length!")
            for i in eachindex(r1)
                add_edge!(G, r1[i], r2[i])
            end
            G
        end

        # Torus:
        nt = 51; mt = 21
        torus = glue!(squarelattice(nt,mt), 1:nt, nt*mt-nt+1:nt*mt)
        torus = glue!(torus, 1:nt:nt*mt, nt:nt:nt*mt)

        # Klein bottle:
        nk = 41; mk = 19
        klein = glue!(squarelattice(nk,mk), 1:nk, nk*mk-nk+1:nk*mk)
        klein = glue!(klein, 1:nk:nk*mk, reverse(nk:nk:nk*mk))

        # Real projective plane:
        np = 31; mp = 31
        pplane = glue!(squarelattice(np,mp), 1:np, reverse(np*mp-np+1:np*mp))
        pplane = glue!(pplane, 1:np:np*mp, reverse(np:np:np*mp))

        results = []
        for (idx, space) in enumerate([torus, klein, pplane])
            res1 = persistence(space, use_hash = true)
            res2 = persistence(space, use_hash = false)

            @test length(res1.generators) == length(res2.generators)
            for i in 1:length(res1.generators)
                @test res1.generators[i].path     == res2.generators[i].path
                @test res1.generators[i].triangle == res2.generators[i].triangle
                @test length(res1.generators[i])  == length(res2.generators[i])
            end

            push!(results, res1)
        end

        # Torus:
        @test length(results[1].generators)   == 2
        @test length(results[1].generators[1]) ≈ mt atol = 2
        @test length(results[1].generators[2]) ≈ nt atol = 2

        # Klein bottle:
        @test length(results[2].generators)   == 2
        @test length(results[2].generators[1]) ≈ mk atol = 2
        @test length(results[2].generators[2]) ≈ nk atol = 2

        # Real projective plane:
        @test length(results[3].generators)   == 1
        @test length(results[3].generators[1]) ≈ mp atol = 2

        # TODO: make a shorter path and make sure it's found.
    end

    @testset "points" begin

        circ = [sin.(linspace(0, 2π, 1000))'; cos.(linspace(0, 2π, 1000))']
        res  = persistence(circ)
        @test length(res.generators)   == 1
        @test length(res.generators[1]) ≈ 2π atol = 0.01

        rnd = rand(3, 1000)
        @test length(persistence(rnd).generators) == 0

        # TODO:
        # * randomly sampled spaces with large R
    end
end


text:
using FileIO, Compat
using FactCheck

# Stub readers---these might bork any existing readers, so don't
# run these tests while doing other things!
module TestLoadSave
import FileIO: File, @format_str
load(file::File{format"PBMText"})   = "PBMText"
load(file::File{format"PBMBinary"}) = "PBMBinary"
load(file::File{format"HDF5"})      = "HDF5"
load(file::File{format"JLD"})       = "JLD"
end

sym2loader = copy(FileIO.sym2loader)
sym2saver = copy(FileIO.sym2saver)

try
    empty!(FileIO.sym2loader)
    empty!(FileIO.sym2saver)
    file_dir = joinpath(dirname(@__FILE__), "files")
    context("Load") do

        add_loader(format"PBMText", :TestLoadSave)
        add_loader(format"PBMBinary", :TestLoadSave)
        add_loader(format"HDF5", :TestLoadSave)
        add_loader(format"JLD", :TestLoadSave)

        @fact load(joinpath(file_dir, "file1.pbm")) --> "PBMText"
        @fact load(joinpath(file_dir, "file2.pbm")) --> "PBMBinary"
        # Regular HDF5 file with magic bytes starting at position 0
        @fact load(joinpath(file_dir, "file1.h5")) --> "HDF5"
        # This one is actually a JLD file saved with an .h5 extension,
        # and the JLD magic bytes edited to prevent it from being recognized
        # as JLD.
        # JLD files are also HDF5 files, so this should be recognized as
        # HDF5. However, what makes this more interesting is that the
        # magic bytes start at position 512.
        @fact load(joinpath(file_dir, "file2.h5")) --> "HDF5"
        # JLD file saved with .jld extension
        @fact load(joinpath(file_dir, "file.jld")) --> "JLD"

        @fact_throws load("missing.fmt")
    end
finally
    merge!(FileIO.sym2loader, sym2loader)
    merge!(FileIO.sym2saver, sym2saver)
end

# A tiny but complete example
# DUMMY format is:
#  - n, a single Int64
#  - a vector of length n of UInt8s

add_format(format"DUMMY", b"DUMMY", ".dmy")

module Dummy

using FileIO, Compat

function FileIO.load(file::File{format"DUMMY"})
    open(file) do s
        skipmagic(s)
        load(s)
    end
end

function FileIO.load(s::Stream{format"DUMMY"})
    # We're already past the magic bytes
    n = read(s, Int64)
    out = Array(UInt8, n)
    read!(s, out)
    close(s)
    out
end

function FileIO.save(file::File{format"DUMMY"}, data)
    open(file, "w") do s
        write(s, magic(format"DUMMY"))  # Write the magic bytes
        write(s, convert(Int64, length(data)))
        udata = convert(Vector{UInt8}, data)
        write(s, udata)
    end
end

end

add_loader(format"DUMMY", :Dummy)
add_saver(format"DUMMY", :Dummy)

context("Save") do
    a = [0x01,0x02,0x03]
    fn = string(tempname(), ".dmy")
    save(fn, a)

    b = load(query(fn))
    @fact a --> b

    b = open(query(fn)) do s
        skipmagic(s)
        load(s)
    end
    @fact a --> b

    # low-level I/O test
    open(query(fn)) do s
        @fact position(s) --> 0
        skipmagic(s)
        @fact position(s) --> length(magic(format"DUMMY"))
        seek(s, 1)
        @fact position(s) --> 1
        seekstart(s)
        @fact position(s) --> 0
        seekend(s)
        @fact eof(s) --> true
        skip(s, -position(s)+1)
        @fact position(s) --> 1
        @fact isreadonly(s) --> true
        @fact isopen(s) --> true
        @fact read(s, 2) --> b"UM"
    end
    rm(fn)

    @fact_throws save("missing.fmt", 5)
end

del_format(format"DUMMY")

# PPM/PBM can be either binary or text. Test that the defaults work,
# and that we can force a choice.
module AmbigExt
using Compat
import FileIO: File, @format_str, Stream, stream, skipmagic

load(f::File{format"AmbigExt1"}) = open(f) do io
    skipmagic(io)
    readstring(stream(io))
end
load(f::File{format"AmbigExt2"}) = open(f) do io
    skipmagic(io)
    readstring(stream(io))
end

save(f::File{format"AmbigExt1"}, testdata) = open(f, "w") do io
    s = stream(io)
    print(s, "ambigext1")
    print(s, testdata)
end
save(f::File{format"AmbigExt2"}, testdata) = open(f, "w") do io
    s = stream(io)
    print(s, "ambigext2")
    print(s, testdata)
end
end

context("Ambiguous extension") do
    add_format(format"AmbigExt1", "ambigext1", ".aext", [:AmbigExt])
    add_format(format"AmbigExt2", "ambigext2", ".aext", [:AmbigExt])
    A = "this is a test"
    fn = string(tempname(), ".aext")
    # Test the forced version first: we wouldn't want some method in Netpbm
    # coming to the rescue here, we want to rely on FileIO's logic.
    # `save(fn, A)` will load Netpbm, which could conceivably mask a failure
    # in the next line.
    save(format"AmbigExt2", fn, A)

    B = load(fn)
    @fact B --> A
    @fact typeof(query(fn)) --> File{format"AmbigExt2"}
    rm(fn)

    save(fn, A)
    B = load(fn)
    @fact B --> A
    @fact typeof(query(fn)) --> File{format"AmbigExt1"}

    rm(fn)
end

context("Absent file") do
    @fact_throws SystemError load("nonexistent.oops")
end


text:
using CSV, DataFrames, Query
using Gadfly
using DifferentialEquations

println("Loaded libraries")

function lotka_volterra(du, u, p, t)
	du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
	du[2] = -p[3] * u[2] + p[4] * u[1] * u[2]
end

function solve_lotka_volterra(p, u0, tspan)
	prob = ODEProblem(lotka_volterra, u0, tspan, p)
	return solve(prob, RK4(), dt = 0.005, dtmin = 0.001, dtmax = 0.01)
end

function experiment(id, p, u0)
	solution = solve_lotka_volterra(p, u0, (0.0, 10.0))
	dataframe = DataFrame(
		time = solution.t,
		x = map(row -> row[1], solution.u),
		y = map(row -> row[2], solution.u),
		experiment = "exp$(id)"
	)

	CSV.write("exp$(id).csv", dataframe, header=true)
end

function gather_results(ids...)
	results = DataFrame(
		time=Float64[],
		x=Float64[],
		y=Float64[],
		experiment=String[]
	)

	for id in ids[1:end]
		df = CSV.read("exp$(id).csv")
		results = vcat(results, df)
	end

	return results
end

function plot_results(results)
	set_default_plot_size(40cm, 15cm)
	plot(results, xgroup="experiment",
		Guide.manual_color_key("Legend",
			["pred", "prey", "diff"],
			["red", "green", "lightblue"]
		),
		Geom.subplot_grid(
			layer(x="time", y="diff", Geom.line, Theme(default_color=colorant"lightblue")),
			layer(x="time", y="pred", Geom.line, Theme(default_color=colorant"red")),
			layer(x="time", y="prey", Geom.line, Theme(default_color=colorant"green"))
		),
		Guide.xlabel("time")
	)
end

function plot_phase_space(results)
	set_default_plot_size(20cm, 20cm)
	plot(results,
		x="prey", y="pred",
		color="experiment",
		Geom.path
	)
end

println("Solving the ODE")

experiment(1, [4, 4, 1, 1], [1.0; 2.0])
experiment(2, [6, 2, 1, 2], [2.0; 4.0])
experiment(3, [3, 5, 5, 6], [2.0; 1.0])
experiment(4, [2, 3, 5, 7], [3.0; 1.0])

println("ODE solved")

gathered = gather_results(1, 2, 3, 4)

summary_table = by(gathered, :experiment, df -> DataFrame(
	prey_mean = mean(df[:x]),
	pred_mean = mean(df[:y]),
	prey_min = minimum(df[:x]),
	pred_min = minimum(df[:y]),
	prey_max = maximum(df[:x]),
	pred_max = maximum(df[:y])
))

println("Summary:")
println(summary_table)

with_difference = @from r in gathered begin
	@select {
		time = r.time,
		prey = r.x,
		pred = r.y,
		diff = r.x - r.y,
		experiment = r.experiment
	}
	@orderby time
	@collect DataFrame
end

set_default_plot_size(30cm, 40cm)
vstack(
	plot_results(with_difference),
	plot_phase_space(with_difference)
)


text:
"""
Return a vector with one less dimension
with i-th item removed.

Arguments
----------
input:Array
    input vector

i:Integer
    which item to remove
"""
function smaller_vector(input::Array, i::Integer)
    return vcat(input[1:i-1],input[i+1:end])
end


text:
struct FieldNode <: FlatBuffers.Struct
    bytes::Vector{UInt8}
    pos::Base.Int
end

FlatBuffers.structsizeof(::Base.Type{FieldNode}) = 16

Base.propertynames(x::FieldNode) = (:length, :null_count)

function Base.getproperty(x::FieldNode, field::Symbol)
    if field === :length
        return FlatBuffers.get(x, FlatBuffers.pos(x), Int64)
    elseif field === :null_count
        return FlatBuffers.get(x, FlatBuffers.pos(x) + 8, Int64)
    end
    return nothing
end

function createFieldNode(b::FlatBuffers.Builder, length::Int64, nullCount::Int64)
    FlatBuffers.prep!(b, 8, 16)
    prepend!(b, nullCount)
    prepend!(b, length)
    return FlatBuffers.offset(b)
end

FlatBuffers.@scopedenum CompressionType::Int8 LZ4_FRAME ZSTD

FlatBuffers.@scopedenum BodyCompressionMethod::Int8 BUFFER

struct BodyCompression <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::BodyCompression) = (:codec, :method)

function Base.getproperty(x::BodyCompression, field::Symbol)
    if field === :codec
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), CompressionType)
        return CompressionType.LZ4_FRAME
    elseif field === :method
        o = FlatBuffers.offset(x, 6)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), BodyCompressionMethod)
        return BodyCompressionMethod.BUFFER
    end
    return nothing
end

bodyCompressionStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 2)
bodyCompressionAddCodec(b::FlatBuffers.Builder, codec::CompressionType) = FlatBuffers.prependslot!(b, 0, codec, 0)
bodyCompressionAddMethod(b::FlatBuffers.Builder, method::BodyCompressionMethod) = FlatBuffers.prependslot!(b, 1, method, 0)
bodyCompressionEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct RecordBatch <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::RecordBatch) = (:length, :nodes, :buffers, :compression)

function Base.getproperty(x::RecordBatch, field::Symbol)
    if field === :length
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int64)
    elseif field === :nodes
        o = FlatBuffers.offset(x, 6)
        if o != 0
            return FlatBuffers.Array{FieldNode}(x, o)
        end
    elseif field === :buffers
        o = FlatBuffers.offset(x, 8)
        if o != 0
            return FlatBuffers.Array{Buffer}(x, o)
        end
    elseif field === :compression
        o = FlatBuffers.offset(x, 10)
        if o != 0
            y = FlatBuffers.indirect(x, o + FlatBuffers.pos(x))
            return FlatBuffers.init(BodyCompression, FlatBuffers.bytes(x), y)
        end
    end
    return nothing
end

recordBatchStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 4)
recordBatchAddLength(b::FlatBuffers.Builder, length::Int64) = FlatBuffers.prependslot!(b, 0, length, 0)
recordBatchAddNodes(b::FlatBuffers.Builder, nodes::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, nodes, 0)
recordBatchStartNodesVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 16, numelems, 8)
recordBatchAddBuffers(b::FlatBuffers.Builder, buffers::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 2, buffers, 0)
recordBatchStartBuffersVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 16, numelems, 8)
recordBatchAddCompression(b::FlatBuffers.Builder, c::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 3, c, 0)
recordBatchEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct DictionaryBatch <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::DictionaryBatch) = (:id, :data, :isDelta)

function Base.getproperty(x::DictionaryBatch, field::Symbol)
    if field === :id
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int64)
        return Int64(0)
    elseif field === :data
        o = FlatBuffers.offset(x, 6)
        if o != 0
            y = FlatBuffers.indirect(x, o + FlatBuffers.pos(x))
            return FlatBuffers.init(RecordBatch, FlatBuffers.bytes(x), y)
        end
    elseif field === :isDelta
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Bool)
        return false
    end
    return nothing
end

dictionaryBatchStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 3)
dictionaryBatchAddId(b::FlatBuffers.Builder, id::Int64) = FlatBuffers.prependslot!(b, 0, id, 0)
dictionaryBatchAddData(b::FlatBuffers.Builder, data::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, data, 0)
dictionaryBatchAddIsDelta(b::FlatBuffers.Builder, isdelta::Base.Bool) = FlatBuffers.prependslot!(b, 2, isdelta, false)
dictionaryBatchEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

function MessageHeader(b::UInt8)
    b == 1 && return Schema
    b == 2 && return DictionaryBatch
    b == 3 && return RecordBatch
    # b == 4 && return Tensor
    # b == 5 && return SparseTensor
    return nothing
end

function MessageHeader(::Base.Type{T})::Int16 where {T}
    T == Schema && return 1
    T == DictionaryBatch && return 2
    T == RecordBatch && return 3
    # T == Tensor && return 4
    # T == SparseTensor && return 5
    return 0
end

struct Message <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Message) = (:version, :header, :bodyLength, :custom_metadata)

function Base.getproperty(x::Message, field::Symbol)
    if field === :version
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), MetadataVersion)
    elseif field === :header
        o = FlatBuffers.offset(x, 6)
        if o != 0
            T = MessageHeader(FlatBuffers.get(x, o + FlatBuffers.pos(x), UInt8))
            o = FlatBuffers.offset(x, 8)
            pos = FlatBuffers.union(x, o)
            if o != 0
                return FlatBuffers.init(T, FlatBuffers.bytes(x), pos)
            end
        end
    elseif field === :bodyLength
        o = FlatBuffers.offset(x, 10)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int64)
        return Int64(0)
    elseif field === :custom_metadata
        o = FlatBuffers.offset(x, 12)
        if o != 0
            return FlatBuffers.Array{KeyValue}(x, o)
        end
    end
    return nothing
end

messageStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 5)
messageAddVersion(b::FlatBuffers.Builder, version::MetadataVersion) = FlatBuffers.prependslot!(b, 0, version, 0)
messageAddHeaderType(b::FlatBuffers.Builder, ::Core.Type{T}) where {T} = FlatBuffers.prependslot!(b, 1, MessageHeader(T), 0)
messageAddHeader(b::FlatBuffers.Builder, header::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 2, header, 0)
messageAddBodyLength(b::FlatBuffers.Builder, bodyLength::Int64) = FlatBuffers.prependslot!(b, 3, bodyLength, 0)
messageAddCustomMetadata(b::FlatBuffers.Builder, meta::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 4, meta, 0)
messageStartCustomMetadataVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 4, numelems, 4)
messageEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


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

include("model_PublicIPAddressDnsSettings.jl")
include("model_PublicIPAddressListResult.jl")
include("model_PublicIPAddressPropertiesFormat.jl")
include("model_Resource.jl")
include("model_SubResource.jl")
include("model_IPConfiguration.jl")
include("model_PublicIPAddress.jl")


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

mutable struct IoK8sApiCoreV1EnvVar <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    value::Any # spec type: Union{ Nothing, String } # spec name: value
    valueFrom::Any # spec type: Union{ Nothing, IoK8sApiCoreV1EnvVarSource } # spec name: valueFrom

    function IoK8sApiCoreV1EnvVar(;name=nothing, value=nothing, valueFrom=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EnvVar, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1EnvVar, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(IoK8sApiCoreV1EnvVar, Symbol("valueFrom"), valueFrom)
        setfield!(o, Symbol("valueFrom"), valueFrom)
        o
    end
end # type IoK8sApiCoreV1EnvVar

const _property_map_IoK8sApiCoreV1EnvVar = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("value")=>Symbol("value"), Symbol("valueFrom")=>Symbol("valueFrom"))
const _property_types_IoK8sApiCoreV1EnvVar = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("value")=>"String", Symbol("valueFrom")=>"IoK8sApiCoreV1EnvVarSource")
Base.propertynames(::Type{ IoK8sApiCoreV1EnvVar }) = collect(keys(_property_map_IoK8sApiCoreV1EnvVar))
Swagger.property_type(::Type{ IoK8sApiCoreV1EnvVar }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EnvVar[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EnvVar }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EnvVar[property_name]

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

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


text:
Combinatorics.permutations(a, c) = permutations(a, c, sum(c))
Combinatorics.permutations(a, c, k::Integer) = MultisetPermuations(a, c, k)

struct MultisetPermuations{T,U,V}
    a::T
    c::U
    k::V
end

Base.start(p::MultisetPermuations) = begin
    s = [fill(i, p.c[i]) for i in eachindex(p.a)] |> flatten
    p.k <= length(s) ? s[1:p.k] : [length(p.a) + 1]
end

Base.next(p::MultisetPermuations, s) = begin
    permutation = [p.a[si] for si in s]
    p.k > 0 || return permutation, [length(p.a) + 1]

    s = copy(s)
    for i in length(s):-1:1
        s[i] += 1
        while s[i] <= length(p.a) && countnz(s .== s[i]) > p.c[s[i]]
            s[i] += 1
        end

        if s[i] <= length(p.a)
            t = [fill(i, p.c[i] - countnz(s .== i)) for i in eachindex(p.a)] |> flatten
            s[i+1:end] = t[1:(length(s)-i)]
            break
        end
    end

    permutation, s
end

Base.done(p::MultisetPermuations, s) = !isempty(s) && s[1] > length(p.a)

Base.eltype(p::MultisetPermuations) = eltype(typeof(p))
Base.eltype(::Type{MultisetPermuations{T,U,V}}) where {T,U,V} = Array{eltype(T),1}

Base.length(p::MultisetPermuations) = begin
    p.k == 0 && return big(1)

    csum = sum(p.c)
    p.k > csum && return big(0)
    p.k == csum && return multinomial(p.c)

    # http://www.m-hikari.com/ams/ams-2011/ams-17-20-2011/siljakAMS17-20-2011.pdf
    R, z = PowerSeriesRing(QQ, p.k + 1, "z")
    gf = prod(sum(divexact(z^j, fac(j)) for j in 0:m) for m in p.c)
    l = fac(p.k) * coeff(gf, p.k)

    @assert denominator(l) == 1
    convert(BigInt, numerator(l))
end


text:
using MultipleViewGeometry, Test
using StaticArrays, Random

# Vectors are scaled so that the last component is unity.
v = [2,2,2]
@test 𝑛(v) == [1.0, 1.0, 1.0]

v = Point2DH([2,2,2])
@test 𝑛(v) == Point2DH([1.0, 1.0, 1.0])

# Vectors which represent points at infinity are unchanged.
v = [2,2,0]
@test 𝑛(v) == [2,2,0]

v = Point2DH([2,2,0])
@test 𝑛(v) == Point2DH([2,2,0])

v = [1, 2, 3]
@test vec2antisym(v) == @SMatrix [ 0 -3  2;
                                   3  0 -1;
                                  -2  1  0]

 v = Point2DH([1, 2, 3])
 @test vec2antisym(v) == @SMatrix [ 0 -3  2;
                                    3  0 -1;
                                   -2  1  0]


text:
"""
    dropgrad(x) -> x

Drop the gradient of `x`.

    julia> gradient(2, 3) do a, b
         dropgrad(a)*b
       end
    (nothing, 2)
"""
dropgrad(x) = x
@adjoint dropgrad(x) = dropgrad(x), _ -> nothing

"""
    ignore() do
      ...
    end

Tell Zygote to ignore a block of code. Everything inside the `do` block will run
on the forward pass as normal, but Zygote won't try to differentiate it at all.
This can be useful for e.g. code that does logging of the forward pass.

Obviously, you run the risk of incorrect gradients if you use this incorrectly.
"""
ignore(f) = f()
@adjoint ignore(f) = ignore(f), _ -> nothing

"""
    hook(x̄ -> ..., x) -> x

Gradient hooks. Allows you to apply an arbitrary function to the gradient for
`x`.

    julia> gradient(2, 3) do a, b
             hook(ā -> @show(ā), a)*b
           end
    ā = 3
    (3, 2)

    julia> gradient(2, 3) do a, b
             hook(-, a)*b
           end
    (-3, 2)
"""
hook(f, x) = x

@adjoint! hook(f, x) = x, x̄ -> (nothing, f(x̄),)

"""
    @showgrad(x) -> x

Much like `@show`, but shows the gradient about to accumulate to `x`. Useful for
debugging gradients.

    julia> gradient(2, 3) do a, b
             @showgrad(a)*b
           end
    ∂(a) = 3
    (3, 2)

Note that the gradient depends on how the output of `@showgrad` is *used*, and is
not the *overall* gradient of the variable `a`. For example:

    julia> gradient(2) do a
         @showgrad(a)*a
       end
    ∂(a) = 2
    (4,)

    julia> gradient(2, 3) do a, b
             @showgrad(a) # not used, so no gradient
             a*b
           end
    ∂(a) = nothing
    (3, 2)
"""
macro showgrad(x)
  :(hook($(esc(x))) do x̄
      println($"∂($x) = ", repr(x̄))
      x̄
    end)
end

"""
    hessian(f, x)

Construct the Hessian of `f`, where `x` is a real or real array and `f(x)` is
a real.

    julia> hessian(((a, b),) -> a*b, [2, 3])
    2×2 Array{Int64,2}:
     0  1
     1  0
"""
hessian(f, x::AbstractArray) = forward_jacobian(x -> gradient(f, x)[1], x)[2]


text:
function compute_update0!(H2, dHdtau, H, min_dxy2, _dx, _dy, dx, dy)
  @inbounds @fastmath for iy = 1:size(dHdtau, 2)
    for ix = 1:size(dHdtau, 1)
      dHdtau[ix, iy] =
        (
          -((
            -((H[(ix+1)+1, iy+1] - H[ix+1, iy+1])) * _dx -
            -((H[ix+1, iy+1] - H[ix, iy+1])) * _dx
          )) * _dx -
          (
            -((H[ix+1, (iy+1)+1] - H[ix+1, iy+1])) * _dy -
            -((H[ix+1, iy+1] - H[ix+1, iy])) * _dy
          ) * _dy
        ) + 0.6 * dHdtau[ix, iy]
      H2[ix+1, iy+1] = H[ix+1, iy+1] + min_dxy2 / 4.1 * dHdtau[ix, iy]  # sets the BC as H[1]=H[end]=0
    end
  end
  return
end
function compute_update1!(H2, dHdtau, H, min_dxy2, _dx, _dy, dx, dy)
  @turbo for iy = 1:size(dHdtau, 2)
    for ix = 1:size(dHdtau, 1)
      dHdtau[ix, iy] =
        (
          -((
            -((H[(ix+1)+1, iy+1] - H[ix+1, iy+1])) * _dx -
            -((H[ix+1, iy+1] - H[ix, iy+1])) * _dx
          )) * _dx -
          (
            -((H[ix+1, (iy+1)+1] - H[ix+1, iy+1])) * _dy -
            -((H[ix+1, iy+1] - H[ix+1, iy])) * _dy
          ) * _dy
        ) + 0.6 * dHdtau[ix, iy]
      H2[ix+1, iy+1] = H[ix+1, iy+1] + min_dxy2 / 4.1 * dHdtau[ix, iy]  # sets the BC as H[1]=H[end]=0
    end
  end
  return
end

@testset "conv" begin
  lx, ly = 10.0, 10.0   # domain size
  nx, ny = 12, 12     # number of grid points
  nt = 1e3          # max number of iterations
  dx, dy = lx / nx, ly / ny
  xc, yc = LinRange(dx / 2, lx - dx / 2, nx), LinRange(dy / 2, ly - dy / 2, ny)
  dHdtau = zeros(nx - 2, ny - 2)
  dHdtau2 = zeros(nx - 2, ny - 2)
  H = exp.(.-(xc .- lx / 2) .^ 2 .- (yc .- ly / 2)' .^ 2)
  H1 = copy(H)
  H2 = copy(H)
  min_dxy2 = min(dx, dy)^2
  _dx, _dy = 1.0 / dx, 1.0 / dy
  compute_update0!(H1, dHdtau, H, min_dxy2, _dx, _dy, dx, dy)
  compute_update1!(H2, dHdtau2, H, min_dxy2, _dx, _dy, dx, dy)
  @test H1 ≈ H2
  @test dHdtau ≈ dHdtau2
end


text:
"""
    Rel

Ability to define frames in a relative fashion.

# Example
```
 Action(1:100, ground; in_global_layer=true),
 Action(1:90, :red_ball, (args...)->circ(p1, "red"), Rotation(from_rot, to_rot)),
 Action(Rel(10), :blue_ball, (args...)->circ(p2, "blue"), Rotation(2π, from_rot, :red_ball)),
 Action((video, args...)->path!(path_of_red, pos(:red_ball), "red"))
```
is the same as
```
Action(1:100, ground; in_global_layer=true),
Action(1:90, :red_ball, (args...)->circ(p1, "red"), Rotation(from_rot, to_rot)),
Action(91:100, :blue_ball, (args...)->circ(p2, "blue"), Rotation(2π, from_rot, :red_ball)),
Action(91:100, (video, args...)->path!(path_of_red, pos(:red_ball), "red"))
```

# Fields
- rel::UnitRange defines the frames in a relative fashion.
"""
struct Rel
    rel::UnitRange
end

"""
    Rel(i::Int)

Shorthand for Rel(1:i)
"""
Rel(i::Int) = Rel(1:i)


text:
using Documenter
using Stopping

makedocs(
    sitename = "Stopping.jl",
    format = Documenter.HTML(assets = ["assets/style.css"], prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [Stopping],
    pages = [
             "Home" => "index.md",
             "API" => "api.md",
             "Examples and tutorials" => "tutorial.md",
             "How to State" => "howtostate.md",
             "How to State for NLPs" => "howtostate-nlp.md",
             "How to Stop" => "howtostop.md",
             "How to Stop 2" => "howtostop-2.md",
             "How to Stop for NLPs" => "howtostop-nlp.md",
             "Solve linear algebra" => "linear-algebra.md",
             "Use a buffer function" => "buffer.md",
             "A fixed point algorithm" => "fixed-point.md",
             "Backtracking linesearch algorithm" => "backls.md",
             "Unconstrained optimization algorithm" => "uncons.md",
             "Active set algorithm" => "active-set.md",
             "Quadratic penalty algorithm" => "penalty.md",
             "Run optimization algorithms" => "run-optimsolver.md",
             "Benchmark optimization algorithms" => "benchmark.md"
            ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/Goysa2/Stopping.jl")#


text:
ENV["GKSwstype"] = "100"
using RandomLinearAlgebraSolvers
using Documenter

DocMeta.setdocmeta!(
  RandomLinearAlgebraSolvers,
  :DocTestSetup,
  :(using RandomLinearAlgebraSolvers);
  recursive = true,
)

makedocs(;
  modules = [RandomLinearAlgebraSolvers],
  authors = "Tangi Migot tangi.migot@gmail.com",
  repo = "https://github.com/tmigot/RandomLinearAlgebraSolvers.jl/blob/{commit}{path}#{line}",
  sitename = "RandomLinearAlgebraSolvers.jl",
  format = Documenter.HTML(;
    prettyurls = get(ENV, "CI", "false") == "true",
    canonical = "https://tmigot.github.io/RandomLinearAlgebraSolvers.jl",
    assets = String[],
  ),
  pages = [
    "Home" => "index.md",
    "Benchmark overdetermined" => "rectangular.md",
    "Benchmark sdp matrices" => "sdp.md",
    "Benchmark MatrixDepot" => "matrixdepot.md",
  ],
)

deploydocs(; repo = "github.com/tmigot/RandomLinearAlgebraSolvers.jl", devbranch = "main")


text:
using Documenter
using Controlz

# using Controlz before all docstrings
DocMeta.setdocmeta!(Controlz, :DocTestSetup, :(using Controlz); recursive=true)

makedocs(
    root = joinpath(dirname(pathof(Controlz)), "..", "docs"),
    modules = [Controlz],
 #     source = "src",
    doctest = true,
    sitename = "Controlz.jl",
    clean = true,
    pages = ["Controlz" => "index.md",
             "Transfer Functions" => "tfs.md",
             "Simulation" => "sim.md",
             "Visualization" => "viz.md",
             "Control systems" => "controls.md",
             "FAQ" => "faq.md"
             ],
    format = Documenter.HTML(assets = ["assets/flux.css"])
)

deploydocs(
    repo = "github.com/SimonEnsemble/Controlz.jl.git",
    devbranch = "main",
    push_preview=true
)

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text: Dict{String, Any}("title" => Any[Dict{String, Any}("value" => "Client: \"XXXX\", Job: XXXX", "type" => "string"), Dict{String, Any}("value" => "Code: XXXX", "type" => "string")])


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

"""
    length(element)

Return the length of basis (number of nodes).
"""
function length(element::AbstractElement)
    return length(element.properties)
end

"""
    size(element)

Return the size of basis (dim, nnodes).
"""
function size(element::AbstractElement)
    return size(element.properties)
end

function getindex(element::AbstractElement, field_name::String)
    return get_field(element, Symbol(field_name))
end

function setindex!(element::AbstractElement, data::T, field_name) where T<:AbstractField
    element.fields[field_name] = data
end

function setindex!(element::AbstractElement, data::Function, field_name)
    if hasmethod(data, Tuple{AbstractElement, Vector, Float64})
        # create enclosure to pass element as argument
        element.fields[field_name] = field((ip,time) -> data(element,ip,time))
    else
        element.fields[field_name] = field(data)
    end
end

function setindex!(element::AbstractElement, data, field_name)
    element.dfields[field_name] = field(data)
end

function setindex!(fieldset::Dict{Symbol, AbstractField}, field_data, field_name::String)
    setindex!(fieldset, field_data, Symbol(field_name))
end

function has_field(element, field_name::String)
    return has_field(element, Symbol(field_name))
end

#""" Return a Field object from element.
#Examples
#--------
#>>> element = Element(Seg2, [1, 2])
#>>> data = Dict(1 => 1.0, 2 => 2.0)
#>>> update!(element, "my field", data)
#>>> element("my field")
#"""
function (element::Element)(field_name::String)
    return element[field_name]
end

function size(element::AbstractElement, dim)
    return size(element)[dim]
end

""" Check existence of field. """
function haskey(element::AbstractElement, field_name)
    return has_field(element, field_name)
end

# will be deprecated
function assemble!(::Assembly, ::Problem{P}, ::AbstractElement, ::Any) where P
    @warn("One must define assemble! function for problem of type $P. " *
          "Not doing anything.")
    return nothing
end

# will be deprecated
function assemble!(assembly::Assembly, problem::Problem{P},
                   elements::Vector{Element}, time) where P
    @warn("This is default assemble! function. Decreased performance can be " *
          "expected without preallocation of memory. One should implement " *
          "`assemble_elements!(problem, assembly, elements, time)` function.")
    for element in elements
        assemble!(assembly, problem, element, time)
    end
    return nothing
end

# generally a bad idea to have functions like update! and interpolate, which
# are not explicitly giving targe?

function update!(element::AbstractElement, field_name, field_data)
    update_field!(element, field_name, field_data)
end

function update!(field::AbstractField, data)
    update_field!(field, data)
end

function update!(element::AbstractElement, field_name::String, field_data)
    update_field!(element, Symbol(field_name), field_data)
end

function update!(elements::Vector{Element}, field_name::String, field_data)
    update_field!(elements, Symbol(field_name), field_data)
end

function update!(elements::Vector{Element{T}}, field_name::String, field_data) where T
    update_field!(convert(Vector{Element}, elements), field_name, field_data)
end

function interpolate(element::AbstractElement, field_name::String, time)
    interpolate_field(element, Symbol(field_name), time)
end

function interpolate(element::AbstractElement, field_name::String, ip, time)
    interpolate(element, Symbol(field_name), ip, time)
end

# fields is now dfields
function Base.getproperty(element::Element, sym::Symbol)
    if sym == :fields
        return getfield(element, :dfields)
    else
        return getfield(element, sym)
    end
end

# getindex, when someone asks key element.fields["f"] => element.fields[:f]
function Base.getindex(fieldset::Dict{Symbol,AbstractField}, field_name::String)
    return getindex(fieldset, Symbol(field_name))
end


text:
#
# Perturbed set of parameters for next trial
#
function nextP!(p,options)
  @. p = p + (options.p0_range[2]-options.p0_range[1])*(randn()*options.rand_std)
end


text:
module Saveresult_module
    struct Saveresult_param
        alpha::Float64
        b::Float64
        dx::Float64
        s::Float64
        rarray::Vector{Float64}
        xarray::Vector{Float64}
        yarray::Vector{Float64}
        Z::Float64
    end
end