text:
#=
  to maintain the package intent correctly
     implicit construction of a Float64 does not require the target become a Float32

  to maintain the package intent correctly
    some primitive operations must be taken with respect to Float32

  to maintain the package intent correctly
     explicit construction of T requires the target become a Float32
=#

"""
    Robust32s

A module for robust Float32 computation.

Exports: `Robust32`, `ComplexR32`
""" Robust32s

module Robust32s

export Robust32, ComplexR32

import Base: convert, promote_rule, show, string,
             Float64, Float32, Float16, prevfloat, nextfloat,
             ==, !=, <, <=, >, >=, isless, isequal, +, -, *, \, /, ^, fma, muladd,
             signbit, precision, significand, exponent, sign, eps, inv, sqrt, cbrt, hypot, clamp, clamp!,
             min, max, minmax, frexp, ldexp, abs, copysign, flipsign, zero, one, iszero, isone,
             isfinite, issubnormal, isinf, isnan, float, floatmin, floatmax, maxintfloat, typemax, typemin,
             mod, rem, div, fld, cld, divrem, fldmod,
             floor, ceil, trunc, round,
             rand, randn, evalpoly,
             sum, prod, cumsum, cumprod, accumulate, sum!, prod!, cumsum!, cumprod!, accumulate!

import Base.Math: abs2, acos, acosd, acosh, acot, acotd, acoth, acsc, acscd, acsch, asec, asecd, asech,
                  asin, asind, asinh, atan, atand, atanh, cos, cosc, cosd, cosh, cospi, cot, cotd, coth,
                  csc, cscd, csch, deg2rad, evalpoly, exp, exp10, exp2, expm1,
                  log, log10, log1p, log2, mod2pi, modf, rad2deg, rem2pi, sec, secd, sech,
                  sin, sinc, sincos, sincosd, sind, sinh, sinpi, tan, tand, tan #  sincospi

using Random
using LinearAlgebra
# using Gaius

struct As64 end # internal use only

struct Robust32 <: AbstractFloat
    val::Float64

    Robust32(x::Float64) = new(Float64(Float32(x)))
    Robust32(::Type{As64}, x::Float64) = new(x)
end

Robust32(x::Robust32) = x # idempotency

value64(x::Robust32) = x.val
value32(x::Robust32) = Float32(x.val)

const ComplexR32 = Complex{Robust32}

value64(x::ComplexR32) = (value64(x.re), value64(x.im))
value32(x::ComplexR32) = (value32(x.re), value32(x.im))

show(io::IO, x::Robust32) = show(io, value32(x))
string(x::Robust32) = string(value32(x))

# internal use only
Rob32(x::Real) = Robust32(As64, Float64(x))
Rob32(x::Float64) = Robust32(As64, x)
Rob32(x::Robust32) = x
Rob32(x::Float32) = Rob32(Meta.parse(string(x)))
Rob32(x::Float16) = Rob32(Meta.parse(string(x)))

# process `divrem`, `fldmod` automatically
Rob32(x::Tuple{Float64,Float64}) = (Robust32(As64, x[1]), Robust32(As64, x[2]))

Float64(x::Robust32) = Float64(value32(x))
convert(::Type{Float64}, x::Robust32) = Float64(x)
promote_rule(::Type{Robust32}, ::Type{Float64}) = Robust32
convert(::Type{Robust32}, x::Float64) = Rob32(x)
# Robust32(x::Float64) defined as an internal constructor

Float32(x::Robust32) = value32(x)
Robust32(x::Float32) = Rob32(Float64(x))
convert(::Type{Float32}, x::Robust32) = Float32(x)
promote_rule(::Type{Robust32}, ::Type{Float32}) = Robust32
convert(::Type{Robust32}, x::Float32) = Robust32(x)

for T in (:BigFloat, :Float16)
  @eval begin
    Base.$T(x::Robust32) = $T(value64(x))
    Robust32(x::$T) = Rob32(Float64(x))
    convert(::Type{$T}, x::Robust32) = $T(x)
    promote_rule(::Type{Robust32}, ::Type{$T}) = Robust32
    convert(::Type{Robust32}, x::$T) = Robust32(x)
  end
end

for T in (:BigInt, :Int128, :Int64, :Int32, :Int16, :Int8,
                   :UInt128, :UInt64, :UInt32, :UInt16, :UInt8)
  @eval begin
    Base.$T(x::Robust32) = $T(value32(x))
    Robust32(x::$T) = Rob32(Float64(Float32(x)))
    convert(::Type{$T}, x::Robust32) = $T(x)
    promote_rule(::Type{Robust32}, ::Type{$T}) = Robust32
    convert(::Type{Robust32}, x::$T) = Robust32(x)
  end
end

convert(::Type{Rational{T}}, x::Robust32) where {T} = convert(Rational{T}, value32(x))
convert(::Type{Rational}, x::Robust32) = convert(Rational{Int64}, x)
promote_rule(::Type{Robust32}, ::Type{Rational}) = Robust32

Robust32(x::T) where {T<:Irrational} = Robust32(Float64(x))
convert(::Type{Robust32}, x::T) where {T<:Irrational} = Rob32(Float64(x))
promote_rule(::Type{Robust32}, ::Type{Irrational})  = Robust32

const Robust32_0 = Rob32(0.0)
const Robust32_1 = Rob32(1.0)
const Robust32_2 = Rob32(2.0)

Robust32(x::Bool) = x ? Robust32_1 : Robust32_0
Bool(x::Robust32) = iszero(x) ? false : true
Base.hash(x::Robust32, h::UInt64) = Base.hash(value32(x), h)

Base.decompose(x::Robust32) = Base.decompose(value32(x))
Base.rtoldefault(x::Robust32) = Base.rtoldefault(Float32(x))

precision(::Type{Robust32}) = precision(Float32)

for F in (:floatmin, :floatmax, :maxintfloat, :typemax, :typemin)
  @eval $F(::Type{Robust32}) = Robust32($F(Float32))
end

eps(::Type{Robust32}) = Robust32(eps(Float32))
eps(x::Robust32) = Robust32(eps(value32(x)))

significand(x::Robust32) = significand(value32(x))
exponent(x::Robust32) = exponent(value32(x))
sign(x::Robust32) = sign(value32(x))

iszero(x::Robust32) = iszero(value32(x))
isone(x::Robust32) = isone(value32(x))
isfinite(x::Robust32) = isfinite(value32(x))
issubnormal(x::Robust32) = issubnormal(value32(x))
isinf(x::Robust32) = isinf(value64(x))
isnan(x::Robust32) = isnan(value64(x))

signbit(x::Robust32) = signbit(value32(x))

zero(::Type{Robust32}) = Robust32_0
one(::Type{Robust32}) = Robust32_1
two(::Type{Robust32}) = Robust32_2
zero(x::Robust32) = zero(Robust32)
one(x::Robust32) = one(Robust32)
two(x::Robust32) = two(Robust32)

nextfloat(x::Robust32, n::Int=1) = Rob32(nextfloat(value64(x),n))
prevfloat(x::Robust32, n::Int=1) = Rob32(prevfloat(value64(x),n))

#=
   >>> Important Implementation Note <<<

   Base.frexp(x::Robust32) = frexp64(x)
   Base.ldexp(fr::Robust32, xp::Int) = ldexp64(fr, xp)

   -----

   frexp32(::Robust32)
   ldexp32(::Robust32, ::Int)

   This pair is to be used with input data (`ldexp32`)
   and output values (`frexp32`) should the need arise.

   For example, if one prefer to generate
   initial input data values by combining
   (significand, exponent) pairings,
   use `ldexp32.(significands, exponents)`.

   Similarly, if one prefer to store results
   in the form of (significand, exponent) pairings
   use `frexp32.(results)`.

   ----

   frexp64(::Robust32)
   ldexp64(::Robust32, ::Int)

   This pair works with 64-bit values directly.
   That design decision makes the following hold:
      `ldexp64(frexp64(x::Robust32)) == x`,

   _otherwise it would not hold much of the time_
   _consider_ `x = sqrt(Robust32(2.0))`.
=#

frexp(x::Robust32) = frexp64(x)
ldexp(fr::Robust32, xp::Int) = ldexp64(fr, xp)

function frexp64(x::Robust32)
  fr64, xp = frexp(value64(x))
  fr = Rob32(fr64)
  return (fr, xp)
end

function ldexp64(fr::Robust32, ex::Int)
  fr64 = value64(fr)
  return Rob32(ldexp(fr64, ex))
end
ldexp64(x::Tuple{Robust32, Int}) = ldexp64(x[1], x[2])

function frexp32(x::Robust32)
  fr64, xp = frexp(Float64(x))
  fr = Rob32(fr64)
  return (fr, xp)
end

function ldexp32(fr::Robust32, ex::Int)
  fr64 = Float64(fr)
  return Rob32(ldexp(fr64, ex))
end
ldexp32(x::Tuple{Robust32, Int}) = ldexp32(x[1], x[2])

for F in (:-, :abs, :inv, :sqrt, :cbrt)
  @eval $F(x::Robust32) = Rob32($F(value64(x)))
end

for F in (:(==), :(!=), :(<), :(<=), :(>), :(>=), :isless, :isequal)
  @eval begin
    $F(x::Robust32, y::Robust32) = $F(value32(x), value32(y))
    $F(x::Robust32, y::Float64) = $F(value32(x), y)
    $F(x::Robust32, y::Float32) = $F(value32(x), y)
    $F(x::Float64, y::Robust32) = $F(x, value32(y))
    $F(x::Float32, y::Robust32) = $F(x, value32(y))
    $F(x::AbstractFloat, y::Robust32) = $F(promote(x,y)...)
    $F(x::Robust32, y::AbstractFloat) = $F(promote(x,y)...)
    $F(x::Rational, y::Robust32) = $F(promote(x,y)...)
    $F(x::Robust32, y::Rational) = $F(promote(x,y)...)
  end
end

^(x::Robust32, y::Robust32) = Rob32(value64(x)^value64(y))
^(x::Robust32, y::Integer) = Rob32(value64(x)^y)
^(x::Real, y::Robust32) = Rob32(Float64(x^value64(y)))

for F in (:+, :-, :*, :/, :\, :hypot, :copysign, :flipsign,
          :mod, :rem, :div, :fld, :cld)
  @eval begin
    $F(x::Robust32, y::Robust32) = Rob32($F(value64(x), value64(y)))
    $F(x::Robust32, y::Float32) = $F(x, Rob32(y))
    $F(x::Float32, y::Robust32) = $F(Rob32(x), y)
    $F(x::Robust32, y::Float16) = $F(x, Rob32(y))
    $F(x::Float16, y::Robust32) = $F(Rob32(x), y)
    $F(x::Robust32, y::T) where {T<:Real} = $F(x, Rob32(y))
    $F(x::T, y::Robust32) where {T<:Real} = $F(Rob32(x), y)
  end
end

# divrem(x, y) = divrem(x, y, RoundToZero)
# fldmod(x,y) = divrem(x, y, RoundDown)
divrem(x::Robust32, y::Robust32) = Rob32(divrem(value64(x), value64(y), RoundToZero))
divrem(x::Robust32, y::Float64) = Rob32(divrem(value64(x), y, RoundToZero))
divrem(x::Robust32, y::Float32) = Rob32(divrem(value64(x), Float64(y), RoundToZero))
divrem(x::Float64, y::Robust32) = Rob32(divrem(x, value64(y), RoundToZero))
divrem(x::Float32, y::Robust32) = Rob32(divrem(Float64(x), value64(y), RoundToZero))

fldmod(x::Robust32, y::Robust32) = Rob32(fldmod(value64(x), value64(y), RoundDown))
fldmod(x::Robust32, y::Float64) = Rob32(fldmod(value64(x), y, RoundDown))
fldmod(x::Robust32, y::Real) = Rob32(fldmod(value64(x), Float64(y), RoundDown))
fldmod(x::Float64, y::Robust32) = Rob32(fldmod(x, value64(y), RoundDown))
fldmod(x::Real, y::Robust32) = Rob32(fldmod(Float64(x), value64(y), RoundDown))

for F in (:hypot, :clamp)
  @eval begin
    $F(x::Robust32, y::Robust32, z::Robust32) = Robust32($F(value64(x), value64(y), value64(z)))
    $F(x::Robust32, y::Robust32, z::T) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::T, y::Robust32, z::Robust32) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::Robust32, y::T, z::Robust32) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::Robust32, y::T, z::T) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::T, y::Robust32, z::T) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::T, y::T, z::Robust32) where {T<:Real} = $F(promote(x,y,z)...)
  end
end

for F in (:floor, :ceil, :trunc)
  @eval Base.$F(x::Robust32) = Rob32($F(value64(x)))
    for T in (:BigInt, :Int128, :Int64, :Int32, :Int16, :Int8)
      @eval Base.$F(::Type{$T}, x::Robust32) = $T(Float32($F(value64(x))))
    end
end

Base.round(x::Robust32; digits=0, sigdigits=0, base=10) =
    iszero(sigdigits) ? Rob32(round(value64(x), digits=digits, base=base)) : Rob32(round(value64(x), sigdigits=sigdigits, base=base))

for T in (:BigInt, :Int128, :Int64, :Int32, :Int16, :Int8)
  @eval Base.round(::Type{$T}, x::Robust32) = round($T, Float64(x))
end

include("minmax.jl")

for F in (:abs2, :acos, :acosd, :acosh, :acot, :acotd, :acoth, :acsc, :acscd, :acsch, :asec,
          :asecd, :asech, :asin, :asind, :asinh, :atan, :atand, :atanh, :cos, :cosc,
          :cosd, :cosh, :cospi, :cot, :cotd, :coth, :csc, :cscd, :csch, :deg2rad,
          :exp, :exp10, :exp2, :expm1, :log, :log10, :log1p, :log2, :mod2pi,
          :rad2deg, :rem2pi, :sec, :secd, :sech, :sin, :sinc, :sind, :sinh,
          :sinpi, :tan, :tand, :tanh)
    @eval $F(x::Robust32) = Rob32($F(value64(x)))
end

atan(x::Robust32, y::Robust32) = Rob32(atan(value64(x), value64(y)))

for F in (:modf, :sincos, :sincosd) # , :sincospi)
  @eval function $F(x::Robust32)
            s, c = $F(value64(x))
            return Rob32(s), Rob32(c)
         end
end

for F in (:fma, :muladd)
  @eval begin
    $F(x::Robust32, y::Robust32, z::Robust32) = Robust32($F(value64(x), value64(y), value64(z)))
    $F(x::Robust32, y::Robust32, z::T) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::T, y::Robust32, z::Robust32) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::Robust32, y::T, z::Robust32) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::Robust32, y::T, z::T) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::T, y::Robust32, z::T) where {T<:Real} = $F(promote(x,y,z)...)
    $F(x::T, y::T, z::Robust32) where {T<:Real} = $F(promote(x,y,z)...)
  end
end

function evalpoly(x::Robust32, p::NTuple{N, Robust32}) where {N}
    Rob32(evalpoly(value64(x), map(value64, p)))
end
function evalpoly(x::T, p::NTuple{N, Robust32}) where {T,N}
    Rob32(evalpoly(Float64(x), map(value64, p)))
end
function evalpoly(x::Robust32, p::NTuple{N, T}) where {T,N}
    Rob32(evalpoly(value64(x), p))
end

#=
     provide(x) uses reinterpret

     provide(x::Vector{Robust32}), provide(x::Matrix{Robust32})
        presents an Array of Robust32s as an Array of Float64s
     provide(x::Vector{Float64}), provide(x::Matrix{Float64})
        presents an Array of Float64s as an Array of Robust32s
     rewrap(x) uses unsafe_wrap

     rewrap(x::Vector{Robust32}, rewrap(x::Matrix{Robust32})
        obtains an Array of Robust32s as an Array of Float64s
     rewrap(x::Vector{Float64}, rewrap(x::Matrix{Float64})
        obtains an Array of Float64s  as an Array of Robust32s
=#

for (R,F) in ((:Robust32, :Float64), (:ComplexR32, :ComplexF64))
  for A in (:Vector, :Matrix)
    @eval begin
        provide(x::$A{$R}) = reinterpret($F, x)
        provide(x::$A{$F}) = reinterpret($R, x)

        rewrap(x::$A{$R}) = unsafe_wrap($A{$F}, cvtptr($F, x), size(x))
        rewrap(x::$A{$F}) = unsafe_wrap($A{$R}, cvtptr($R, x), size(x))
    end
  end
end

@inline cvtptr(::Type{T}, m::Array{S,N}) where {N,T,S} =
    convert(Ptr{T}, pointer(m,1))

include("randnums.jl")
include("linearalgebra.jl")

end  # Robust32s


text:
@params mutable struct GESOResult{T}
    topology::AbstractVector{T}
    objval::T
    change::T
    converged::Bool
    fevals::Int
end
@params struct GESO{T} <: TopOptAlgorithm
    obj
    constr
    vars::AbstractVector{T}
    topology::AbstractVector{T}
    Pcmin::T
    Pcmax::T
    Pmmin::T
    Pmmax::T
    Pen::T
    string_length::Int
    var_volumes::AbstractVector{T}
    cum_var_volumes::AbstractVector{T}
    order::AbstractVector{Int}
    genotypes::BitArray{2}
    children::BitArray{2}
    var_black::BitVector
    maxiter::Int
    penalty
    sens::AbstractVector{T}
    old_sens::AbstractVector{T}
    obj_trace::MVector{10, T}
    tol::T
    sens_tol::T
    result::GESOResult{T}
end

function GESO(obj::Objective{<:ComplianceFunction}, constr::Constraint{<:VolumeFunction}; maxiter = 1000, tol = 0.001, p = 3., Pcmin = 0.6, Pcmax = 1., Pmmin = 0.5, Pmmax = 1., Pen = 3., sens_tol = tol/100, string_length = 4)
    penalty = obj.solver.penalty
    penalty = @set penalty.p = p
    T = typeof(obj.comp)
    nel = getncells(obj.problem.ch.dh.grid)
    black = obj.problem.black
    white = obj.problem.white
    nvars = nel - sum(black) - sum(white)
    vars = zeros(T, nvars)

    topology = zeros(T, nel)
    result = GESOResult(topology, T(NaN), T(NaN), false, 0)
    sens = zeros(T, nvars)
    old_sens = zeros(T, nvars)
    obj_trace = zeros(MVector{10, T})
    var_volumes = constr.cellvolumes[.!black .& .!white]
    cum_var_volumes = zeros(T, nvars)
    order = zeros(Int, nvars)
    genotypes = trues(string_length, nvars)
    children = trues(string_length, nvars)
    var_black = trues(nvars)

    return GESO(obj, constr, vars, topology, Pcmin, Pcmax, Pmmin, Pmmax, Pen, string_length, var_volumes, cum_var_volumes, order, genotypes, children, var_black, maxiter, penalty, sens, old_sens, obj_trace, tol, sens_tol, result)
end

update_penalty!(b::GESO, p::Number) = (b.penalty.p = p)

get_progress(current_volume, total_volume, design_volume) = clamp(min((total_volume - current_volume) / (total_volume - design_volume), current_volume / design_volume), 0, 1)

get_probs(b::GESO, Prg) = (b.Pcmin + (b.Pcmax - b.Pcmin)*Prg^b.Pen, b.Pmmin + (b.Pmmax - b.Pmmin)*Prg^b.Pen)

function crossover!(children, genotypes, i, j)
    for k in 1:size(genotypes, 1)
        r = rand()
        if r < 0.5
            children[k,i] = genotypes[k,i]
        else
            children[k,i] = genotypes[k,j]
        end
    end
    return
end

function update!(var_black, children, genotypes, Pc, Pm, high_class, mid_class, low_class)
    topology_changed = false
    while !topology_changed
        for i in high_class
            r = rand()
            j = i
            if length(high_class) > 1
                if r < Pc
                    while i == j
                        j = rand(high_class)
                    end
                elseif r < 0.5 + 0.5*Pc
                    j = rand(mid_class)
                else
                    j = rand(low_class)
                end
            else
                if r < 0.5
                    j = rand(mid_class)
                else
                    j = rand(low_class)
                end
            end
            crossover!(children, genotypes, i, j)
        end
        for i in mid_class
            r = rand()
            j = i
            if length(mid_class) > 1
                if r < Pc
                    while i == j
                        j = rand(mid_class)
                    end
                elseif r < 0.5 + 0.5*Pc
                    j = rand(high_class)
                else
                    j = rand(low_class)
                end
            else
                if r < 0.5 + 0.5*Pc
                    j = rand(high_class)
                else
                    j = rand(low_class)
                end
            end
            crossover!(children, genotypes, i, j)
        end
        for i in low_class
            r = rand()
            j = i
            if length(low_class) > 1
                if r < Pc
                    while i == j
                        j = rand(low_class)
                    end
                elseif r < 0.5 + 0.5*Pc
                    j = rand(mid_class)
                else
                    j = rand(high_class)
                end
            else
                if r < 0.5
                    j = rand(mid_class)
                else
                    j = rand(high_class)
                end
            end
            crossover!(children, genotypes, i, j)
        end
        genotypes .= children

        for i in high_class
            for j in 1:size(genotypes, 1)
                r = rand()
                if r < Pm && !genotypes[j,i]
                    genotypes[j,i] = !genotypes[j,i]
                end
            end
            if any(@view genotypes[:,i]) != var_black[i]
                var_black[i] = !var_black[i]
                topology_changed = true
            end
        end
        for i in mid_class
            for j in 1:size(genotypes, 1)
                r = rand()
                if r < Pm && genotypes[j,i]
                    genotypes[j,i] = !genotypes[j,i]
                end
            end
            if any(@view genotypes[:,i]) != var_black[i]
                var_black[i] = !var_black[i]
                topology_changed = true
            end
        end
        for i in low_class
            for j in 1:size(genotypes, 1)
                r = rand()
                if r < Pm && genotypes[j,i]
                    genotypes[j,i] = !genotypes[j,i]
                end
            end
            if any(@view genotypes[:,i]) != var_black[i]
                var_black[i] = !var_black[i]
                topology_changed = true
            end
        end
    end

    return var_black
end

function (b::GESO{TO, TC, T})(x0=copy(b.obj.solver.vars); seed=NaN) where {TO<:Objective{<:ComplianceFunction}, TC<:Constraint{<:VolumeFunction}, T}
    @unpack sens, old_sens, tol, maxiter = b
    @unpack obj_trace, topology, sens_tol, vars = b
    @unpack Pcmin, Pcmax, Pmmin, Pmmax, Pen = b
    @unpack string_length, genotypes, children, var_black = b
    @unpack cum_var_volumes, var_volumes, order = b
    @unpack varind, black, white = b.obj.f.problem
    @unpack total_volume, cellvolumes, fixed_volume = b.constr.f
    V = b.constr.s
    design_volume = V * total_volume

    nel = length(x0)
    nvars = length(vars)

    # Set seed
    isnan(seed) || Random.seed!(seed)

    # Initialize the topology
    for i in 1:length(topology)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = round(x0[varind[i]])
            vars[varind[i]] = topology[i]
        end
    end

    check(x) = x > design_volume - fixed_volume
    #rrmax = clamp(1 - design_volume/current_volume, 0, 1)
    current_volume = dot(vars, var_volumes) + fixed_volume
    vol = current_volume/total_volume
    # Main loop
    change = T(1)
    iter = 0
    while (change > tol || vol > V) && iter < maxiter
        iter += 1
        if iter > 1
            old_sens .= sens
        end
        for j in max(2, 10-iter+2):10
            obj_trace[j-1] = obj_trace[j]
        end
        obj_trace[10] = b.obj(vars, sens)
        rmul!(sens, -1)
        if iter > 1
            @. sens = (sens + old_sens) / 2
        end

        # Classify the cells by their sensitivities
        sortperm!(order, sens, rev=true)
        accumulate!(+, cum_var_volumes, view(var_volumes, order))
        N1 = findfirst(check, cum_var_volumes) - 1
        N2 = (nel - N1) ÷ 2
        N3 = nvars - N1 - N2
        high_class = @view order[1:N1]
        mid_class = @view order[N1+1:N1+N2]
        low_class = @view order[N1+N2+1:end]

        # Crossover and mutation
        Prg = get_progress(current_volume, total_volume, design_volume)
        Pc, Pm = get_probs(b, Prg)
        vars .= update!(var_black, children, genotypes, Pc, Pm, high_class, mid_class, low_class)

        # Update crossover and mutation probabilities
        current_volume = dot(vars, var_volumes) + fixed_volume
        vol = current_volume/total_volume

        if iter >= 10
            l = sum(@view obj_trace[1:5])
            h = sum(@view obj_trace[6:10])
            change = abs(l-h)/h
        end
    end

    for i in 1:length(topology)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = vars[varind[i]]
        end
    end

    b.result.objval = obj_trace[10]
    b.result.change = change
    b.result.converged = change <= tol
    b.result.fevals = iter

    return b.result
end

ext: jl


text:
struct Hundsdorfer <: FdmSchemeDescType end
struct Douglas <: FdmSchemeDescType end

mutable struct FdmSchemeDesc{F <: FdmSchemeDescType}
  schemeType::F
  theta::Float64
  mu::Float64
end

FdmSchemeDesc(t::Hundsdorfer) = FdmSchemeDesc{Hundsdorfer}(t, 0.5 + sqrt(3.0) / 6.0, 0.5)
FdmSchemeDesc(t::Douglas) = FdmSchemeDesc{Douglas}(t, 0.5, 0.0)

## Schemes ##
mutable struct HundsdorferScheme <: FdScheme
  theta::Float64
  mu::Float64
  map::FdmLinearOpComposite
  bcSet::FdmBoundaryConditionSet
  dt::Float64
end

HundsdorferScheme(theta::Float64, mu::Float64, map::FdmLinearOpComposite, bcSet::FdmBoundaryConditionSet) = HundsdorferScheme(theta, mu, map, bcSet, 0.0)

mutable struct DouglasScheme <: FdScheme
  theta::Float64
  map::FdmLinearOpComposite
  bcSet::FdmBoundaryConditionSet
  dt::Float64
end

DouglasScheme(theta::Float64, map::FdmLinearOpComposite, bcSet::FdmBoundaryConditionSet) = DouglasScheme(theta, map, bcSet, 0.0)

mutable struct MixedScheme <: FdScheme
  L::TridiagonalOperator
  Ident::TridiagonalOperator
  dt::Float64
  theta::Float64
  bcSet::FdmBoundaryConditionSet
  explicitPart::TridiagonalOperator
  implicitPart::TridiagonalOperator

  function MixedScheme(L::TridiagonalOperator, theta::Float64, bcSet::FdmBoundaryConditionSet)
    Ident = TridiagIdentity(L.n)
    dt = 0.0
    return new(L, Ident, dt, theta, bcSet)
  end
end

# add'l constructors
CrankNelson(L::TridiagonalOperator, bcs::FdmBoundaryConditionSet) = MixedScheme(L, 0.5, bcs)

function set_step!(evolver::MixedScheme, dt::Float64)
  evolver.dt = dt
  if evolver.theta != 1.0 # there is an explicit part
    evolver.explicitPart = evolver.Ident - ((1.0 - evolver.theta) * evolver.dt) * evolver.L
  end

  if evolver.theta != 0.0 # there is an implicit part
    evolver.implicitPart = evolver.Ident + (evolver.theta * evolver.dt) * evolver.L
  end

  return evolver
end

function step!(evolver::MixedScheme, a::Vector{Float64}, t::Float64)
  for i in eachindex(evolver.bcSet.conditions)
    set_time!(evolver.bcSet.conditions[i], t)
  end

  if evolver.theta != 1.0 # there is an explicit part
    # TODO Tridiagonal time dependent
    for i in eachindex(evolver.bcSet.conditions)
      apply_before_applying!(evolver.bcSet.conditions[i], evolver.explicitPart)
    end
    a[:] = apply_to(evolver.explicitPart, a)
    for i in eachindex(evolver.bcSet.conditions)
      apply_after_applying!(evolver.bcSet.conditions[i], a)
    end
  end

  if evolver.theta != 0.0 # there is an implicit part
    # TODO Tridiagonal time dependent
    for i in eachindex(evolver.bcSet.conditions)
      apply_before_solving!(evolver.bcSet.conditions[i], evolver.implicitPart, a)
    end
    a[:] = solve_for(evolver.implicitPart, a)
    for i in eachindex(evolver.bcSet.conditions)
      apply_after_solving!(evolver.bcSet.conditions[i], a)
    end
  end

  return a, evolver
end

mutable struct ParallelEvolver{E <: FdScheme} <: FdScheme
  evolvers::Vector{E}
end

function ParallelEvolver(L::Vector{TridiagonalOperator}, bcs::Matrix{BC}, evolverFunc::Function) where {BC <: BoundaryCondition}
  # build first for Vector
  evolv1 = evolverFunc(L[1], FdmBoundaryConditionSet(bcs[:, 1]))
  evolvers = Vector{typeof(evolv1)}(undef, length(L))
  evolvers[1] = evolv1

  @simd for i = 2:length(L)
    @inbounds evolvers[i] = evolverFunc(L[i], FdmBoundaryConditionSet(bcs[:, i]))
  end

  return ParallelEvolver(evolvers)
end

function step!(evolver::ParallelEvolver, a::Vector{Vector{Float64}}, t::Float64)
  for i in eachindex(evolver.evolvers)
    step!(evolver.evolvers[i], a[i], t)
  end

  return a, evolver
end

function set_step!(evolver::ParallelEvolver, dt::Float64)
  for i in eachindex(evolver.evolvers)
    set_step!(evolver.evolvers[i], dt)
  end

  return evolver
end

set_step!(evolver::FdScheme, dt::Float64) = evolver.dt = dt

function step!(evolver::HundsdorferScheme, a::Vector{Float64}, t::Float64)
  t - evolver.dt > -1e-8 || error("a step towards negative time given")

  set_time!(evolver.map, max(0.0, t - evolver.dt), t)
  set_time!(evolver.bcSet, max(0.0, t - evolver.dt))

  apply_before_applying!(evolver.bcSet, evolver.map)
  y = a + evolver.dt * apply(evolver.map, a)
  apply_after_applying!(evolver.bcSet, y)

  y0 = copy(y)

  @simd for i = 1:get_size(evolver.map)
    @inbounds rhs = y - evolver.theta * evolver.dt * apply_direction(evolver.map, i, a)
    @inbounds y = solve_splitting(evolver.map, i, rhs, -evolver.theta * evolver.dt)
  end

  apply_before_applying!(evolver.bcSet, evolver.map)
  yt = y0 + evolver.mu * evolver.dt * apply(evolver.map, y - a)
  apply_after_applying!(evolver.bcSet, yt)

  @simd for i = 1:get_size(evolver.map)
    @inbounds rhs = yt - evolver.theta * evolver.dt * apply_direction(evolver.map, i, y)
    @inbounds yt = solve_splitting(evolver.map, i, rhs, -evolver.theta * evolver.dt)
  end

  apply_after_solving!(evolver.bcSet, yt)

  a[:] = yt

  return a, evolver
end

function step!(evolver::DouglasScheme, a::Vector{Float64}, t::Float64)
  t - evolver.dt > -1e-8 || error("a step towards negative time given")

  set_time!(evolver.map, max(0.0, t - evolver.dt), t)
  set_time!(evolver.bcSet, max(0.0, t - evolver.dt))

  apply_before_applying!(evolver.bcSet, evolver.map)
  y = a + evolver.dt * apply(evolver.map, a)
  apply_after_applying!(evolver.bcSet, y)

  @simd for i = 1:get_size(evolver.map)
    @inbounds rhs = y - evolver.theta * evolver.dt * apply_direction(evolver.map, i, a)
    @inbounds y = solve_splitting(evolver.map, i, rhs, -evolver.theta * evolver.dt)
  end

  apply_after_solving!(evolver.bcSet, y)

  a[:] = y

  return a, evolver
end


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

@doc raw"""ExternalMetricSource indicates how to scale on a metric not associated with any Kubernetes object (for example length of queue in cloud messaging service, or QPS from loadbalancer running outside of cluster). Exactly one \"target\" type should be set.

    IoK8sApiAutoscalingV2beta1ExternalMetricSource(;
        metricName=nothing,
        metricSelector=nothing,
        targetAverageValue=nothing,
        targetValue=nothing,
    )

    - metricName::String : metricName is the name of the metric in question.
    - metricSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : metricSelector is used to identify a specific time series within a given metric.
    - targetAverageValue::IoK8sApimachineryPkgApiResourceQuantity : targetAverageValue is the target per-pod value of global metric (as a quantity). Mutually exclusive with TargetValue.
    - targetValue::IoK8sApimachineryPkgApiResourceQuantity : targetValue is the target value of the metric (as a quantity). Mutually exclusive with TargetAverageValue.
"""
mutable struct IoK8sApiAutoscalingV2beta1ExternalMetricSource <: SwaggerModel
    metricName::Any # spec type: Union{ Nothing, String } # spec name: metricName
    metricSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: metricSelector
    targetAverageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: targetAverageValue
    targetValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: targetValue

    function IoK8sApiAutoscalingV2beta1ExternalMetricSource(;metricName=nothing, metricSelector=nothing, targetAverageValue=nothing, targetValue=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricSource, Symbol("metricName"), metricName)
        setfield!(o, Symbol("metricName"), metricName)
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricSource, Symbol("metricSelector"), metricSelector)
        setfield!(o, Symbol("metricSelector"), metricSelector)
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricSource, Symbol("targetAverageValue"), targetAverageValue)
        setfield!(o, Symbol("targetAverageValue"), targetAverageValue)
        validate_property(IoK8sApiAutoscalingV2beta1ExternalMetricSource, Symbol("targetValue"), targetValue)
        setfield!(o, Symbol("targetValue"), targetValue)
        o
    end
end # type IoK8sApiAutoscalingV2beta1ExternalMetricSource

const _property_map_IoK8sApiAutoscalingV2beta1ExternalMetricSource = Dict{Symbol,Symbol}(Symbol("metricName")=>Symbol("metricName"), Symbol("metricSelector")=>Symbol("metricSelector"), Symbol("targetAverageValue")=>Symbol("targetAverageValue"), Symbol("targetValue")=>Symbol("targetValue"))
const _property_types_IoK8sApiAutoscalingV2beta1ExternalMetricSource = Dict{Symbol,String}(Symbol("metricName")=>"String", Symbol("metricSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("targetAverageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("targetValue")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricSource }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1ExternalMetricSource))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1ExternalMetricSource[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1ExternalMetricSource }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1ExternalMetricSource[property_name]

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

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


text:
using ChemicalEngineeringDesign
using Test

@testset "ChemicalEngineeringDesign.jl" begin
    # Write your own tests here.
end


text:
using CoupledDipoles, Revise
using Random
Random.seed!(1134)

using Statistics: mean, var
using StatsBase: fit, normalize, Histogram
using ProgressMeter, Plots, LinearAlgebra

### ------------ ATOMS SPECS ---------------------
ρλ³ = 45
kL = 32.4;
ρ_over_k₀³ = ρλ³ / (2π)^3
N = floor(Int, ρ_over_k₀³ * kL^3)

### ------------ LASER SPECS ---------------------
Δ = 1.0
s = 1e-6

### -------- PRODUCE INTENSITIES -----------------
sensors = get_sensors_ring(; num_pts=360, kR=1.5kL, θ=5π / 12)
many_intensities = Float64[]
maxRep = 40

@showprogress for rep in 1:maxRep
    atoms = Atom(Cube(), N, kL)
    laser = Laser(Gaussian3D(kL/8), s, Δ)
    simulation = LinearOptics(Scalar(), atoms, laser)

    βₙ = get_steady_state(simulation)
    intensities = get_intensities_over_sensors(simulation, βₙ, sensors)

    for i in intensities
        push!(many_intensities, copy(i))
    end
end
all_intensities_over_mean = many_intensities ./ mean(many_intensities)

### ------------ CREATE HISTOGRAM ---------------------
bins = 10.0 .^ range(log10(1e-10), log10(100); length=100)
h = fit(Histogram, all_intensities_over_mean, bins)

h_norm = normalize(h; mode=:pdf)
bins_edges = collect(h_norm.edges[1])
bins_centers = [
    sqrt(bins_edges[i] * bins_edges[i + 1]) for i in 1:(length(bins_edges) - 1)
]

### ------------ FINAL PLOT ---------------------
## theoretical curve
x_ray = range(0.01, 100; step=0.15)
y_ray = exp.(-x_ray)
plot(x_ray, y_ray; linestyle=:dash, label="Rayleigh", lw=4)

scatter!(
    bins_centers, h_norm.weights; label="N=$(N), kL=$(kL)", markershape=:circle, markersize=5
)
plot!(; guidefont=15, tickfont=15, legendfont=10, size=(1000, 500))
plot!(;
    left_margin=5Plots.mm,
    right_margin=5Plots.mm,
    top_margin=5Plots.mm,
    bottom_margin=5Plots.mm,
    c=:blue,
    gridalpha=0.8,
    minorgrid=true,
    minorgridalpha=0.5,
    minorgridstyle=:dashdot,
    yscale=:log10,
    xscale=:log10,
    xlims=(10^-1, 10^2),
    xticks=[10^-1, 10^0, 10^1, 10^2],
    ylims=(10^-6, 10^1),
    yticks=[10^1, 10^0, 10^-2, 10^-4, 10^-6],
)

plot!(; ylims=(1e-6, 10), xlims=(1e-1, 100), scale=:log10)
xlabel!("I")
ylabel!("P(I)")
title!("variance = $(  round(var(all_intensities_over_mean),digits=3 ))")


text:
export SVMLearner

type SVMLearner <: Learner
    lambda::Float64
    p::Int
    coefs::Vector{Float64}
    gr::Vector{Float64}
    optimizer::AbstractSGD
    initialized::Bool
    function SVMLearner(optimizer::AbstractSGD, lambda::Float64)
        lambda >= 0.0 || error("lambda should be non-negative")
        obj = new()
        obj.optimizer = optimizer
        obj.lambda = lambda
        obj.initialized = false
        obj
    end
end

SVMLearner(optimizer::AbstractSGD, lambda::Float64) = SVMLearner{Float64}(optimizer, lambda)

function Base.show(io::IO, obj::SVMLearner)
    print(io, "SVMLearner lambda: $(obj.lambda)")
    print(io, "\nOptimizer: ")
    show(io, obj.optimizer)
    if obj.initialized
        print(io, "\nCoefficients: ")
        show(io, obj.coefs)
    end
end

function init!(obj::SVMLearner, p)
    obj.initialized && error("already initialized")
    obj.coefs = zeros(p)
    obj.gr = Array(Float64, p)
    obj.initialized = true
end

function grad!(obj::SVMLearner, x, y)
   error()
end

function update!(obj::SVMLearner, x::DSMat{Float64}, y::Vector{Float64})
    obj.initialized || init!(obj, size(x, 2))
    grad!(obj, x, y)
    update!(obj.optimizer, obj.coefs, obj.gr)
    obj
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

module ComputabilityTheory

using Printf: @printf

include("utils.jl")

# Abstract Types
export Machine, TuringMachine, MachineComponent, Programme,
        ProgrammeComponent

export ∸, pair_tuple, algebraic, π, cℤ, cℤ⁻¹

export Tape, Left, Stay, Right, MachineState, Rule, TMProgramme,
        run_turing_machine, RegisterMachine, run_goto_programme

export Sequence, Instruction, GoToProgramme, increment, decrement,
        goto, ifzero_goto, halt, show_programme, rand

include("abstract_types.jl")
include("coding.jl")
include("goto.jl")
include("machines.jl")

end # end module


text:
# Below url is pinned to a specific commit -- if xiaoliangstd updates the repo
# with changes, you may want to investigate them.
mini_cheetah_url = "https://raw.githubusercontent.com/xiaoliangstd/MiniCheetah_description/f8e30e3d9f1077d096483010ddfbaa16e928a37c/cheetah_description/"

urdfpath = "xacro/mini_cheetah.urdf"

meshpaths = [
    "meshes/mini_body.obj";
    "meshes/mini_abad.obj";
    "meshes/mini_upper_link.obj";
    "meshes/mini_lower_link.obj"
    ]

datadir = "cheetah_description"
ispath(datadir) || mkpath(datadir)
download(mini_cheetah_url * urdfpath, joinpath(datadir, "mini_cheetah.urdf"))

for meshpath in meshpaths
    meshdir, meshfilename = splitdir(meshpath)
    meshdir = joinpath(datadir, meshdir)
    ispath(meshdir) || mkpath(meshdir)
    download(mini_cheetah_url * meshpath, joinpath(datadir, meshpath))
end


text: ## run through weave...


text:
@testset "Testing Simple Ops" begin
    @testset "Testing Eltwise" begin
        # Testing normal eltwise functions

        ## Linear
        f = (x, α, β) -> (α .* x) .+ β
        x = rand(Float32, 10, 10)
        X = OneDNN.Memory(x)
        linear = OneDNN.Linear(1, 2)
        Y = OneDNN.eltwise(linear, X)
        Z = OneDNN.materialize(Y)
        @test isapprox(Z, f(x, 1, 2))

        # Try mutating version
        Z .= zero(eltype(Z))
        @test all(iszero, Z)
        OneDNN.eltwise!(OneDNN.Memory(Z), X, OneDNN.forward_expand(linear)...)
        @test isapprox(Z, f(x, 1, 2))

        # Backprop
        ff(x) = 0.5 .* x .+ 2
        y, back = Zygote._pullback(ff, x)
        z = ones(Float32, size(y))
        Z = OneDNN.Memory(z)
        result = OneDNN.eltwise_backward(OneDNN.Linear(0.5, 2), Z, X)

        expected = back(z)[2]
        @test isapprox(OneDNN.materialize(result), expected)

        #####
        ##### Aliased Functions
        #####

        fns = [abs, Flux.sigmoid, sqrt, Flux.relu, log]
        x = rand(Float32, 100, 100)
        X = OneDNN.Memory(x)
        for f in fns
            @show f
            F(x) = f.(x)
            y, back_expected = Zygote._pullback(F, x)
            Y, back_result = Zygote._pullback(OneDNN.eltwise, f, X)

            @test isapprox(OneDNN.materialize(Y), y)
            @inferred Zygote._pullback(OneDNN.eltwise, f, X)

            dy = rand(Float32, size(y))
            dY = OneDNN.Memory(dy)

            dX = back_result(dY)[3]
            dx = back_expected(dy)[2]
            @test isapprox(dx, OneDNN.materialize(dX))
            @inferred back_result(dY)

            # Make sure the alias rules work as well.
            YY, back_result_2 = Zygote._pullback(f, X)
            @test isapprox(OneDNN.materialize(YY), y)
            dXX = back_result_2(dY)[2]
            @test isapprox(dx, OneDNN.materialize(dXX))
            @inferred back_result_2(dY)
        end
    end

    @testset "Testing Binary" begin
        x = randn(Float32, 100, 100)
        y = randn(Float32, 100, 100)

        X = OneDNN.Memory(x)
        Y = OneDNN.Memory(y)

        fns = [+, -, *, /, max, min]
        for f in fns
            @show f
            z = f.(x, y)
            Z = f.(X, Y)
            @test isa(Z, OneDNN.Memory)
            @test isapprox(z, OneDNN.materialize(Z))
        end
    end
end


text:
# module CodeAnalysis

# # TODO: Make OS-independent
# const newline = "\n"
# struct JuliaModule
#   modulename::AbstractString
#   code::AbstractString
#   submodules::Array{JuliaModule}
# end
const DEBUG = false

# TODO: add filter for comments
has_a_module(filename) = occursin("\nmodule ", open(f->read(f, String), filename))
has_an_include(filename) = occursin("\ninclude(", open(f->read(f, String), filename))

"""
    get_single_src(filename)

Returns `contents`, a string containing the
entire source code from `filename`
recursively through `include` statements.
"""
function get_single_src(filename)
  contents = open(f->read(f, String), filename)
  if has_an_include(filename)
    # TODO: Avoid commented `include`'s (requires greedy check of #-\n and scope-check on #= =#)
    incl = "include("
    contents_split = split(contents, "\n"*incl)
    subcontents = Dict()
    for c in contents_split[2:end]
      DEBUG && println("-------")
      i = findfirst(')', c)
      DEBUG && @show filename
      DEBUG && @show i
      # @show c
      jp = "joinpath("
      k_raw = c[1:i-1]
      k_raw = replace(k_raw, jp => "")
      DEBUG && @show k_raw

      k = k_raw
      DEBUG && @show k

      k_processed = k
      if occursin(",", k)
        k_processed = replace(k_processed, "\""=> "")
        k_processed = split(k_processed, ",")
        filter!(x->!(x=="\n"), k_processed)
        filter!(x->!(x==""), k_processed)
        k_processed = strip.(k_processed)
        k_processed = convert.(String, k_processed)
        k_processed = joinpath(k_processed...)
        s_to_replace = "$(jp)$(k_raw))" # or s_to_replace = "$(jp)$(k_raw)\n)"
      else
        k_processed = replace(k_processed, "\""=> "")
        s_to_replace = "$(incl)$(k_raw))" # or s_to_replace = "$(jp)$(k_raw)\n)"
      end
      i = findfirst(s_to_replace, contents)
      DEBUG && @show s_to_replace
      DEBUG && @show i
      DEBUG && println("-------")
      subcontents[k_processed] = get_single_src(joinpath(dirname(filename), k_processed))
      contents = replace(contents, s_to_replace=>subcontents[k_processed])
    end
  end
  return contents
end

# contents = get_single_src(joinpath(@__DIR__, "..", "src", "Common", "MoistThermodynamics", "MoistThermodynamics.jl")); # works
# contents = get_single_src(joinpath(@__DIR__, "..", "src", "Atmos", "Model", "AtmosModel.jl")); # works
contents = get_single_src(joinpath(@__DIR__, "..", "src", "CLIMA.jl"));
# for line in split(contents, "\n")
#   @show line
# end
nothing

# function assemble_modules(src_dir, package_name)
#   all_files = [joinpath(root, f) for (root, dirs, files) in Base.Filesystem.walkdir(src_dir) for f in files]
#   filter!(x->endswith(x, ".jl"), all_files)
#   filter!(x->!endswith(x, ".jll"), all_files)
#   root_src = joinpath(src_dir, package_name*".jl")

#   contents = open(f->read(f, String), root_src)
#   contents = open(f->read(f, String), root_src)
#   if has_an_include(root_src)
#   end
#   for line in split(contents, newline)
#     if occursin(line)

#     end
#   end

#   return JuliaModule(package_name, code, submodules)
# end

# function explicit_import(filename, _usingmodule)

#   parameters_found = []
#   contents = open(f->read(f, String), filename)
#   for pp in planet_parameters
#     if occursin(pp, contents)
#       push!(parameters_found, pp)
#     end
#   end

#   i = findfirst("$(_usingmodule)\n", contents)
#   if i ≠ nothing && !isempty(parameters_found)
#     contents = replace(contents, "$(_usingmodule)\n" => "$(_usingmodule): "*join(parameters_found, ", ")*"\n")
#   end
#   if isempty(parameters_found)
#     contents = replace(contents, "$(_usingmodule)\n" => "")
#   end
#   open(filename,"w") do io
#     print(io, contents)
#   end
# end

# src_dir = joinpath(@__DIR__, "..")
# all_files = [joinpath(root, f) for (root, dirs, files) in Base.Filesystem.walkdir(src_dir) for f in files]
# _usingmodule = "using CLIMA.PlanetParameters"

# filter!(x->endswith(x, ".jl"), all_files)
# filter!(x->!endswith(x, "runtests.jl"), all_files)
# filter!(x->occursin("$(_usingmodule)", open(f->read(f, String), x)), all_files)
# filter!(x->!occursin("$(_usingmodule):", open(f->read(f, String), x)), all_files)

# for driver in all_files
#   explicit_import(driver, _usingmodule)
# end

# end


text:
int_rules_1_2_1_8 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.1.8*P(x)*(a+b*x+c*x^2)^p =#
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand(~Pq * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, ~x) && IGtQ(~p, -2))
    @apply_utils Antiderivative(~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => Antiderivative(~x * PolynomialQuotient(~Pq, ~x, ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(Coeff(~Pq, ~x, 0), 0) && Not(MatchQ(~Pq, (~x) ^ ~(m') * ~(u') <-- IntegerQ(m)))))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * ~x) ^ (2 * FracPart(~p)))) * Antiderivative(~Pq * (~b + 2 * ~c * ~x) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 1)], (((~b * f - 2 * ~a * g) + (2 * ~c * f - ~b * g) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * Q - (2 * ~p + 3) * (2 * ~c * f - ~b * g), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~p, -1)))
    @apply_utils Antiderivative(~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x), e = Coeff(~Pq, ~x, Expon(~Pq, ~x))], (e * (~x) ^ (q - 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * (q + 2 * ~p + 1)) + (1 / (~c * (q + 2 * ~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * ExpandToSum(((~c * (q + 2 * ~p + 1) * ~Pq - ~a * e * (q - 1) * (~x) ^ (q - 2)) - ~b * e * (q + ~p) * (~x) ^ (q - 1)) - ~c * e * (q + 2 * ~p + 1) * (~x) ^ q, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(LeQ(~p, -1))))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/1.2.1.8 P(x) (a+b x+c x^2)^p.jl


text:
@testset "Input Widgets" begin

function Item(name::String, child::Widget)
    value = map(x -> "Value: $x", child.value)
    child["margin-top"] = 10
    return Expander(name, expanded = true,
        Box(:v, spacing = 10,
            child,
            Label("", label = value, halign = Center)
        ) # Box
    ) # Expander
end

inputs = @map {
    Slider       = Slider(0:100, hexpand = true),
    SpinButton   = SpinButton(-50:50),
    TextField    = TextField(text = "Lorem ipsum dolor"),
    ColorButton  = ColorButton("#f1a"),
    ToggleButton = ToggleButton("Click-me"),
    CheckBox     = CheckBox(true),
    Switch       = Switch(false, halign = Center)
}

widgets = map(collect(inputs)) do (name, widget)
    on(widget) do x
        @info "Value of $name widget changed to $x"
    end
    Item(string(name), widget)
end

app = Window("Input widgets test", 300, 800,
    Scrolled(
        Box(:v, margin = 20, spacing = 10,
            widgets...
        )
    ) #Scrolled
) # Window

@test value(inputs[:Slider])        == 50
@test value(inputs[:SpinButton])    == 0
@test value(inputs[:TextField])     == "Lorem ipsum dolor"
@test value(inputs[:ColorButton])   == parse(Colorant, "#f1a")
@test value(inputs[:ToggleButton])  == false
@test value(inputs[:CheckBox])      == true
@test value(inputs[:Switch])        == false

@showall app
# @waitfor app.destroy
Lynx.destroy(app)

end # testset


text:
export add!

function add! end

function (sim::Simulation)(::typeof(add!),m::AbstractArray,valueF::Function,maskF = EVERYWHERE ,gridtype::GridType = p̂)
	mask = maskF.(Coordinates(sim.grid,gridtype))
	vals = valueF.(Coordinates(sim.grid,gridtype));
	m[mask] .+= vals[mask];
end

function (sim::Simulation)(::typeof(add!),m::AbstractArray,vals::AbstractArray,maskF = EVERYWHERE ,gridtype::GridType = p̂)
	mask = maskF.(Coordinates(sim.grid,gridtype))
	m[mask] .+= vals[mask];
end


text:
module Zygote

using LinearAlgebra, Statistics
using LinearAlgebra: copytri!, AbstractTriangular

import ZygoteRules: @adjoint, @adjoint!, AContext, adjoint, _pullback, pullback,
  literal_getproperty, literal_getfield

using ChainRules: ChainRules, rrule, unthunk
using IRTools
using MacroTools, Requires
using MacroTools: @forward

import Distributed: pmap, CachingPool, workers
export Params, gradient, jacobian, hessian, pullback, pushforward, @code_adjoint

const Numeric{T<:Number} = Union{T, AbstractArray{<:T}}

include("tools/idset.jl")
include("tools/buffer.jl")
include("tools/builtins.jl")

include("forward/Forward.jl")
using .Forward

include("compiler/reverse.jl")
include("compiler/emit.jl")
include("compiler/chainrules.jl")
include("compiler/interface.jl")
include("compiler/show.jl")

include("lib/grad.jl")
include("lib/lib.jl")
include("lib/number.jl")
include("lib/base.jl")
include("lib/array.jl")
include("lib/buffer.jl")
include("lib/broadcast.jl")
include("lib/forward.jl")
include("lib/utils.jl")
include("lib/range.jl")
@init @require Distances="b4f34e82-e78d-54a5-968a-f98e89d6e8f7" include("lib/distances.jl")
@init @require StatsFuns="4c63d2b9-4356-54db-8cca-17b64c39e42c" include("lib/statsfuns.jl")

# we need to define this late, so that the genfuncs see lib.jl
# Move using statements out of this file to help with sysimage building
using IRTools: varargs!, inlineable!, pis!, slots!
using IRTools.Inner: argnames!, update!
include("compiler/interface2.jl")

include("profiler/Profile.jl")

@init @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
  include("flux.jl")
end

@init @require Colors="5ae59095-9a9b-59fe-a467-6f913c188581" begin
  @nograd Colors.ColorTypes._parameter_upper_bound
end

using InteractiveUtils
precompile() = Requires.@include("precompile.jl")

# helps to work around 265-y issues
function refresh()
  Requires.@include("compiler/interface2.jl")
  precompile()
  return
end

macro profile(ex)
  @capture(ex, f_(x__)) || error("@profile f(args...)")
  quote
    _, back = _pullback($(esc(f)), $(esc.(x)...))
    Profile.juno(Profile.profile(back))
  end
end

end # module


text:
using StructTypes
using JSON3

struct SaveData
  input :: MDInput
  traj :: MDTraj
end

StructTypes.StructType(::Type{MDInput}) = StructTypes.Struct()
StructTypes.StructType(::Type{MDTraj}) = StructTypes.Struct()
StructTypes.StructType(::Type{SaveData}) = StructTypes.Struct()
StructTypes.StructType(::Type{Atoms}) = StructTypes.Struct()

function write(input :: MDInput, traj :: MDTraj, filename :: String)
  S = SaveData(input,traj)
  f = open(filename,"w")
  JSON3.write(f,S)
  close(f)
end

function read(filename :: String)
  input, traj = read_wg(filename)
  GC.gc()
  return input, traj
end

function read_wg(filename :: String)
  f = open(filename,"r")
  S = JSON3.read(f,SaveData)
  traj = MDTraj(S.input.n,S.input.nsave)
  # Need to convert the MDTraj structure because the coordinates
  # are read in a single column, instead of the [n,2] array we want.
  for i in 1:S.input.nsave
    for iat in 1:S.input.n
      traj.atoms[i].x[iat,1] = S.traj.atoms[i].x[iat]
      traj.atoms[i].x[iat,2] = S.traj.atoms[i].x[iat+S.input.n]
      traj.atoms[i].status[iat] = S.traj.atoms[i].status[iat]
    end
    traj.potential[i] = S.traj.potential[i]
    traj.kinetic[i] = S.traj.kinetic[i]
    traj.total[i] = S.traj.total[i]
    traj.time[i] = S.traj.time[i]
    traj.nenc[i] = S.traj.nenc[i]
    traj.U[i] = S.traj.U[i]
    traj.S[i] = S.traj.S[i]
    traj.D[i] = S.traj.D[i]
    traj.I[i] = S.traj.I[i]
  end
  input = S.input
  S = nothing
  return input, traj
end


text:
function process_pen_kwargs(kwargs)
    global _DEFAULT_PEN_ARGS
    D = OrderedDict{Symbol,Any}(kwargs)
    penkwargs = Any[]
    for k in keys(_DEFAULT_PEN_ARGS)
        if k in keys(D)
            push!(penkwargs,(k,pop!(D,k)))
        end
    end
    if length(penkwargs) > 0 && :pen in keys(D)
        error("Both pen and $penkwargs keyword arguments found")
    end
    if length(penkwargs) > 0
        D[:pen] = Pen(;penkwargs...)
    end
    collect(D)
end

function kwargstring(P,D::OrderedDict;semicolon=true)
    result = filterjoin([(x=getfield(P,s);
                    x == D[s] ? "" : "$(string(s))=$(string(x))")
                    for s in keys(D)]...)
    (result == "" ? "" : """$(semicolon ? ";" : "")$result""")
end

function updatedvals(defOptions::OrderedDict,newOptions)
    newOptionsDict = Dict(newOptions)
    for (k,v) in newOptionsDict
        if !(k in keys(defOptions))
            error("Unknown option $k")
        end
    end
    D = deepcopy(defOptions)
    merge!(D,newOptionsDict)
    values(D)
end

"""
    splitkwargs(pooled_kwargs,options_to_separate)

Split a pooled list of keyword arguments into kwargs
in `options_to_separate` and the rest
"""
function splitkwargs(pooled_kwargs,options_to_separate)
    remaining_kwargs = Dict(pooled_kwargs)
    separated_kwargs = Dict{Symbol,Any}()
    for (s,v) in remaining_kwargs
        if s in keys(options_to_separate)
            separated_kwargs[s] = pop!(remaining_kwargs,s)
        end
    end
    map(collect,(separated_kwargs,remaining_kwargs))
end

"""
    filterjoin(args...)

Concatenate nonempty string representations of `args`,
separated by commas
"""
filterjoin(args...) = join(Iterators.filter(s->length(s)>0,map(string,args)),",")

"""
    enclosequote(s)

Enclose `s` in quotation marks, unless it starts
with "Label"
"""
enclosequote(s) = startswith(s,"Label") ? s : "\"$s\""


text:
module ExtendedVariantModeExtendedVariantTrigger

using ..Ahorn, Maple

# obsolete
@mapdef Trigger "ExtendedVariantMode/ExtendedVariantTrigger" ExtendedVariantTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="Gravity", enable::Bool=true, newValue::Integer=10, revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/BooleanExtendedVariantTrigger" BooleanExtendedVariantTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="DisableNeutralJumping", newValue::Bool=true, revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	coversScreen::Bool=false, flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/IntegerExtendedVariantTrigger" IntegerExtendedVariantTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="DashCount", enable::Bool=true, newValue::Int=1, revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	coversScreen::Bool=false, flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/FloatExtendedVariantTrigger" FloatExtendedVariantTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="Gravity", enable::Bool=true, newValue::Number=1.0, revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	coversScreen::Bool=false, flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/BadelineAttackPatternTrigger" BadelineAttackPatternTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="BadelineAttackPattern", enable::Bool=true, newValue::Int=1, revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	coversScreen::Bool=false, flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/DashDirectionTrigger" DashDirectionTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	topLeft::Bool=true, top::Bool=true, topRight::Bool=true, left::Bool=true, right::Bool=true, bottomLeft::Bool=true, bottom::Bool=true, bottomRight::Bool=true,
	coversScreen::Bool=false, flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/DisplaySpeedometerTrigger" DisplaySpeedometerTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	newValue::String="BOTH", revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false, coversScreen::Bool=false,
	flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/DontRefillDashOnGroundTrigger" DontRefillDashOnGroundTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	newValue::String="ON", revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false, coversScreen::Bool=false,
	flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/JumpCountTrigger" JumpCountTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="JumpCount", newValue::Int=1, infinite::Bool=false, revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false,
	coversScreen::Bool=false, flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/MadelineBackpackModeTrigger" MadelineBackpackModeTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	newValue::String="NoBackpack", revertOnLeave::Bool=false, revertOnDeath::Bool=true, delayRevertOnDeath::Bool=false, withTeleport::Bool=false, coversScreen::Bool=false,
	flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/WindEverywhereTrigger" WindEverywhereTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	newValue::String="Left", revertOnLeave::Bool=false, revertOnDeath::Bool=true, withTeleport::Bool=false, delayRevertOnDeath::Bool=false, coversScreen::Bool=false,
	flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

@mapdef Trigger "ExtendedVariantMode/ColorGradeTrigger" ColorGradeTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	colorGrade::String="none", revertOnLeave::Bool=false, revertOnDeath::Bool=true, withTeleport::Bool=false, delayRevertOnDeath::Bool=false, coversScreen::Bool=false,
	flag::String="", flagInverted::Bool=false, onlyOnce::Bool=false)

const placements = Ahorn.PlacementDict(
	"Extended Variant Trigger (Boolean) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		BooleanExtendedVariantTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Integer) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		IntegerExtendedVariantTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Float) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		FloatExtendedVariantTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Badeline Attack Pattern) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		BadelineAttackPatternTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Dash Direction) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		DashDirectionTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Display Speedometer) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		DisplaySpeedometerTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Don't Refill Dash On Ground) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		DontRefillDashOnGroundTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Jump Count) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		JumpCountTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Madeline Backpack Mode) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		MadelineBackpackModeTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Wind Everywhere) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		WindEverywhereTrigger,
		"rectangle"
	),
	"Extended Variant Trigger (Color Grading) (Extended Variant Mode)" => Ahorn.EntityPlacement(
		ColorGradeTrigger,
		"rectangle"
	),
)

Ahorn.editingOptions(trigger::BooleanExtendedVariantTrigger) = Dict{String, Any}(
	"variantChange" => Ahorn.ExtendedVariantDictionary.BooleanVariants
)
Ahorn.editingOptions(trigger::IntegerExtendedVariantTrigger) = Dict{String, Any}(
	"variantChange" => Ahorn.ExtendedVariantDictionary.IntegerVariants
)
Ahorn.editingOptions(trigger::FloatExtendedVariantTrigger) = Dict{String, Any}(
	"variantChange" => Ahorn.ExtendedVariantDictionary.FloatVariants
)
Ahorn.editingOptions(trigger::BadelineAttackPatternTrigger) = Dict{String, Any}(
	"newValue" => Dict{String, Int}(
		"Slow Shots" => 1,
		"Beam > Shot" => 2,
		"Fast double-shots" => 3,
		"Double-shots > Beam" => 4,
		"Beam > Double-shots" => 5,
		"Slow double-shots" => 9,
		"Nothing" => 10,
		"Fast beam" => 14,
		"Slow beam" => 15
	)
)
Ahorn.editingOptions(trigger::DisplaySpeedometerTrigger) = Dict{String, Any}(
	"newValue" => Dict{String, String}(
		"Disabled" => "DISABLED",
		"Horizontal" => "HORIZONTAL",
		"Vertical" => "VERTICAL",
		"Both" => "BOTH"
	)
)
Ahorn.editingOptions(trigger::DontRefillDashOnGroundTrigger) = Dict{String, Any}(
	"newValue" => Dict{String, String}(
		"Default" => "DEFAULT",
		"On" => "ON",
		"Off" => "OFF"
	)
)
Ahorn.editingOptions(trigger::MadelineBackpackModeTrigger) = Dict{String, Any}(
	"newValue" => Dict{String, String}(
		"Default" => "Default",
		"No Backpack" => "NoBackpack",
		"With Backpack" => "Backpack"
	)
)
Ahorn.editingOptions(trigger::WindEverywhereTrigger) = Dict{String, Any}(
	"newValue" => Dict{String, String}(
		"Default" => "Default",
		"Left" => "Left",
		"Right" => "Right",
		"Left Strong" => "LeftStrong",
		"Right Strong" => "RightStrong",
		"Right Crazy" => "RightCrazy",
		"Left On-Off" => "LeftOnOff",
		"Right On-Off" => "RightOnOff",
		"Alternating" => "Alternating",
		"Left On-Off Fast" => "LeftOnOffFast",
		"Right On-Off Fast" => "RightOnOffFast",
		"Down" => "Down",
		"Up" => "Up",
		"Random" => "Random"
	)
)
Ahorn.editingOptions(trigger::ColorGradeTrigger) = Dict{String, Any}(
	"colorGrade" => [
		"none",
		"oldsite",
		"panicattack",
		"templevoid",
		"reflection",
		"credits",
		"cold",
		"hot",
		"feelingdown",
		"golden",
		"max480/extendedvariants/celsius/tetris",
		"max480/extendedvariants/greyscale",
		"max480/extendedvariants/sepia",
		"max480/extendedvariants/inverted",
		"max480/extendedvariants/rgbshift1",
		"max480/extendedvariants/rgbshift2",
		"max480/extendedvariants/hollys_randomnoise"
	]
)

Ahorn.editingOrder(trigger::DashDirectionTrigger) = String["x", "y", "width", "height", "topLeft", "top", "topRight", "left", "right", "bottomLeft", "bottom", "bottomRight"]

# those are internally Integer Extended Variant Triggers that render differently: Jump Count has the "infinite" checkbox, and Badeline Attack Pattern has a dropdown.
Ahorn.editingIgnored(trigger::JumpCountTrigger, multiple::Bool=false) = String["variantChange"]
Ahorn.editingIgnored(trigger::BadelineAttackPatternTrigger, multiple::Bool=false) = String["variantChange"]

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -9.06497171739458 <= q <= 7.221276301954956)

                     add_NL_constraint(m, :((0.6131370972312755 + 0.7644530725075485*$(x[1]) + 0.7481034782956382*$(x[2]) + 0.43160377180111187*$(x[3]))/(1 + exp(-(0.6131370972312755 + 0.7644530725075485*$(x[1]) + 0.7481034782956382*$(x[2]) + 0.43160377180111187*$(x[3])))) + (-0.2235485813631728 + -0.4821458309603419*$(x[1]) + -0.0034703029219445014*$(x[2]) + 0.5577106220080164*$(x[3]))/(1 + exp(-(-0.2235485813631728 + -0.4821458309603419*$(x[1]) + -0.0034703029219445014*$(x[2]) + 0.5577106220080164*$(x[3])))) + (0.3888102587500133 + 0.8414638325842052*$(x[1]) + 0.671097918031764*$(x[2]) + -0.05238562929678103*$(x[3]))/(1 + exp(-(0.3888102587500133 + 0.8414638325842052*$(x[1]) + 0.671097918031764*$(x[2]) + -0.05238562929678103*$(x[3])))) + (-0.9543429273497122 + 0.7406334788256062*$(x[1]) + 0.9764525500434011*$(x[2]) + 0.3990593176587991*$(x[3]))/(1 + exp(-(-0.9543429273497122 + 0.7406334788256062*$(x[1]) + 0.9764525500434011*$(x[2]) + 0.3990593176587991*$(x[3])))) + (-0.7459035549882156 + 0.4045099462664137*$(x[1]) + 0.9504581535443215*$(x[2]) + -0.1195761049288735*$(x[3]))/(1 + exp(-(-0.7459035549882156 + 0.4045099462664137*$(x[1]) + 0.9504581535443215*$(x[2]) + -0.1195761049288735*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Sekira

t = @summary function t(meta, data)
    nt = meta.n_target
    μ_shift, σ_shift = mean_and_se(data.shift)
    μ_norm, σ_norm = mean_and_se(data.norm)
    μ_proj, _, _, σ_proj_l, σ_proj_u = med_and_errs(ratio_of_means(data.hproj, data.vproj))
    μ_coherence, σ_coherence = mean_and_se(data.coherence)
    μ_ref_coherence, σ_ref_coherence = mean_and_se(data.single_coherence)
    μ_loneliness, σ_loneliness = mean_and_se(data.loneliness)
end

for dir in ARGS
    @info "Processing `$dir`"
    RimuIO.save_df(joinpath(dir, "results", "summary.arrow"), summarize(t, dir; skip=50_000))
end


text:
@inline colnames(table)::Vector{Symbol} = collect(propertynames(columns(table)))
@inline getcol(table, x) = getproperty(columns(table), x)

function subset_cols(table, cols::Tuple{Vararg{Symbol}})
    nt = columntable(table)
    newcols = NamedTuple{cols}(nt)
    materializer(table)(newcols)
end
function subset_cols(table::NamedTuple, cols::Tuple{Vararg{Symbol}})
    NamedTuple{cols}(table)
end
subset_cols(table, cols::Symbol...) = subset_cols(table, cols)

subset_cols(table, cols::Tuple{Vararg{Symbol,0}}) = table

ensure_tuple(x) = tuple(x)
ensure_tuple(x::Tuple) = x
ensure_tuple(x::AbstractArray) = Tuple(x)


text:
module CellFieldsTests

using Test
using FillArrays
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.CellData
#using Gridap.FESpaces
using Random

domain = (0,1,0,1)
cells = (3,3)
model = CartesianDiscreteModel(domain,cells)

trian = Triangulation(model)
trian_N =BoundaryTriangulation(model)
trian_D =BoundaryTriangulation(model,tags="tag_8")
trian_S =SkeletonTriangulation(model)
trian_0 =Triangulation(model,Int32[])

trian_sv = view(trian_S,[1,5,4,2])
nv = get_normal_vector(trian_sv)
@test isa(nv,SkeletonPair)
@test isa(nv.plus,CellField)
@test isa(nv.minus,CellField)

ϕ = GenericCellField(get_cell_map(trian),trian,ReferenceDomain())
ϕinv = GenericCellField(lazy_map(inverse_map,get_cell_map(trian)),trian,PhysicalDomain())
xref = get_cell_points(trian)
xphy = CellPoint(get_array(xref),trian,PhysicalDomain())
test_array(ϕ(xref),collect1d(ϕ(xphy)),(i,j)->all(map(≈,i,j)))
test_array(ϕinv(xref),collect1d(ϕinv(xphy)),(i,j)->all(map(≈,i,j)))

x = get_cell_points(trian)
@test DomainStyle(x) == ReferenceDomain()
@test get_array(x) == get_cell_coordinates(trian)
@test get_data(x) == get_cell_ref_coordinates(trian)

px = get_physical_coordinate(trian)
test_array(px(x),collect1d(get_array(x)))

_x = change_domain(x,PhysicalDomain())
@test DomainStyle(_x) == PhysicalDomain()
@test get_array(_x) == get_cell_coordinates(trian)
@test get_data(_x) == get_cell_coordinates(trian)

_x = change_domain(x,ReferenceDomain())
@test DomainStyle(_x) == ReferenceDomain()
@test get_array(x) == get_cell_coordinates(trian)
@test get_data(x) == get_cell_ref_coordinates(trian)

ffun(x) = 2*x[1]
f = CellField(ffun,trian)
fx = f(x)
r = map(xs->ffun.(xs),get_array(x))
r = reshape(r,length(r))
test_array(fx,r,≈)

x_0 = get_cell_points(trian_0)
fx_0 = f(x_0)
test_array(fx_0,collect(fx_0))

n_S = get_normal_vector(trian_S)
x_S = get_cell_points(trian_S)

ts = get_triangulation(∇(f))
tt = get_triangulation(n_S.plus)

@test is_change_possible(ts,tt) == true

nf_S = n_S⋅∇(f)

jnf_S = jump(n_S⋅∇(f))
jnfx_S = jnf_S(x_S)
test_array(jnfx_S,0*collect(jnfx_S))

h = CellField(rand(num_cells(trian_S)),trian_S)*jump(∇(f))
hx_S = h(x_S)
test_array(hx_S,collect(hx_S))

h = 3*mean(f)⋅jump(n_S⋅∇(f))
hx_S = h(x_S)
test_array(hx_S,0*collect(hx_S))

aa(f) = 4*f
test_array((aa∘f)(x),4*r,≈)

f1 = f
f2 = 2*f
b(f1,f2) = f1+f2
test_array((b∘(f1,f2))(x),3*r,≈)

f2 = ones(num_cells(trian))
test_array((b∘(f1,f2))(x),map(i->i.+1,r),≈)

v = GenericCellField(get_cell_shapefuns(trian),trian,ReferenceDomain())
vx = v(x)
test_array(vx,collect(vx))

u = GenericCellField(lazy_map(transpose,get_data(v)),v.trian,v.domain_style)
m = v*u
test_array(m(x),collect(m(x)))
m = ∇(v)⋅∇(u)
test_array(m(x),collect(m(x)))

∇vx = ∇(v)(x)
test_array(∇vx,collect(∇vx))

∇fx = ∇(f)(x)
test_array(∇fx,collect(∇fx))

k = VectorValue(1.0,2.0)
∇kfx = ((∇+k)(f))(x)
test_array(∇kfx,collect(∇kfx))

∇kvx = ((∇+k)(v))(x)
test_array(∇kvx,collect(∇kvx))

β(x) = 2*x[1]
α = CellField(x->2*x,trian)
ax = ((∇+k)(β*α))(x)
test_array(ax,collect(ax))

ν = CellField(x->2*x,trian)
ax =((∇-k)⋅ν)(x)
test_array(ax,collect(ax))

ax =((∇-k)×ν)(x)
test_array(ax,collect(ax))

ax =((∇-k)⊗ν)(x)
test_array(ax,collect(ax))

ax =(∇.*ν)(x)
test_array(ax,collect(ax))

ax =(ν.*ν)(x)
test_array(ax,collect(ax))

ax =((∇-k).*ν)(x)
test_array(ax,collect(ax))

ax =(ν⊗(∇-k))(x)
test_array(ax,collect(ax))

σ(x) = diagonal_tensor(VectorValue(1*x[1],2*x[2]))
Fields.gradient(::typeof(σ)) = x-> ThirdOrderTensorValue{2,2,2,Float64}(1,0,0,0,0,0,0,2)
ax = ((∇+k)(σ⋅α))(x)
test_array(ax,collect(ax))

h = Operation(*)(2,f)
hx = h(x)
test_array(hx,2*fx)

a = fill(2,num_cells(trian))
h = Operation(*)(a,f)
hx = h(x)
test_array(hx,2*fx)

fx = evaluate(ffun,x)
test_array(fx,r,≈)

f_N = CellField(ffun,trian_N)
x_N = get_cell_points(trian_N)
fx_N = f_N(x_N)
test_array(fx_N,collect(fx_N))

n_N = get_normal_vector(trian_N)
nx_N = n_N(x_N)
test_array(nx_N,collect(nx_N))

h = f*n_N
hx = h(x_N)
test_array(hx,collect(hx))

gfun(x) = 3*x
g = CellField(gfun,trian)

h = Operation(*)(f,g)
gx = g(x)
hx = h(x)
r = map((i,j)->broadcast(*,i,j),fx,gx)
test_array(hx,r)

h_N = Operation(*)(f_N,g)
gx_N = g(x_N)
r = map((i,j)->broadcast(*,i,j),fx_N,gx_N)
hx_N = h_N(x_N)
test_array(hx_N,r)

g_D = CellField(gfun,trian_D)

cell_h = rand(num_cells(trian))
h = CellField(cell_h,trian)
test_array(h(x),collect(h(x)))
test_array(h(x_N),collect(h(x_N)))

h_N = (2*f_N+g)⋅g
hx_N = h_N(x_N)
test_array(hx_N,collect(hx_N))

end # module


text:
function SAC_target(π)
    (π⁻, 𝒫, 𝒟, γ; kwargs...) -> begin
        ap, logprob = exploration(actor(π), 𝒟[:sp])
        y = 𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* (min.(value(π⁻, 𝒟[:sp], ap)...) .- exp(𝒫[:SAC_log_α][1]).*logprob)
    end
end

function SAC_deterministic_target(π)
    (π⁻, 𝒫, 𝒟, γ; kwargs...) -> begin
        y = 𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* min.(value(π⁻, 𝒟[:sp], action(actor(π), 𝒟[:sp]))...)
    end
end

function SAC_max_Q_target(π)
    (π⁻, 𝒫, 𝒟, γ; kwargs...) -> begin
        error("not implemented")
        #TODO: Sample some number of actions and then choose the max
    end
end

function SAC_actor_loss(π, 𝒫, 𝒟; info = Dict())
    a, logprob = exploration(π.A, 𝒟[:s])
    ignore() do
        info["entropy"] = -mean(logprob)
    end
    mean(exp(𝒫[:SAC_log_α][1]).*logprob .- min.(value(π, 𝒟[:s], a)...))
end

function SAC_temp_loss(π, 𝒫, 𝒟; info = Dict())
    ignore() do
        info["SAC alpha"] = exp(𝒫[:SAC_log_α][1])
    end
    _, logprob = exploration(π.A, 𝒟[:s])
    target_α = logprob .+ 𝒫[:SAC_H_target]
    -mean(exp(𝒫[:SAC_log_α][1]) .* target_α)
end

function SAC(;π::ActorCritic{T, DoubleNetwork{ContinuousNetwork, ContinuousNetwork}},
              ΔN=50,
              SAC_α::Float32=1f0,
              SAC_H_target::Float32 = Float32(-prod(dim(action_space(π)))),
              π_explore=GaussianNoiseExplorationPolicy(0.1f0),
              SAC_α_opt::NamedTuple=(;),
              a_opt::NamedTuple=(;),
              c_opt::NamedTuple=(;),
              a_loss=SAC_actor_loss,
              c_loss=double_Q_loss(),
              target_fn=SAC_target(π),
              prefix="",
              log::NamedTuple=(;),
              𝒫::NamedTuple=(;),
              param_optimizers=Dict(),
              kwargs...) where T

    𝒫 = (SAC_log_α=[Base.log(SAC_α)], SAC_H_target=SAC_H_target, 𝒫...)
    OffPolicySolver(;agent=PolicyParams(π=π, π_explore=π_explore, π⁻=deepcopy(π)),
                     ΔN=ΔN,
                     𝒫=𝒫,
                     log=LoggerParams(;dir = "log/sac", log...),
                     param_optimizers=Dict(Flux.params(𝒫[:SAC_log_α]) => TrainingParams(;loss=SAC_temp_loss, name="temp_", SAC_α_opt...), param_optimizers...),
                     a_opt=TrainingParams(;loss=a_loss, name=string(prefix, "actor_"), a_opt...),
                     c_opt=TrainingParams(;loss=c_loss, name=string(prefix, "critic_"), epochs=ΔN, c_opt...),
                     target_fn=target_fn,
                     kwargs...)
end


text:
# example of CUE GMM: draws from N(0,1)
using Econometrics, Statistics, LinearAlgebra
y = randn(100)
# 3 moment conditions
moments = theta -> [(y .-theta[1]) ((y.^2.0) .- theta[2]) ((y .- theta[1]).^3.0)]
m = theta -> mean(moments(theta), dims=1)[:]
weight = theta -> inv(NeweyWest(moments(theta)))
obj = theta -> dot(m(theta), weight(theta), m(theta))
theta = [0.0, 1.0]
results = fminunc(obj, theta)

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module BootlegCassette

using IRTools: IRTools, @code_ir, @dynamo, argument!, IR, xcall, isexpr, var

abstract type Context end

macro context(_ctx::Symbol)
    ctx = esc(_ctx)
    quote
        struct $ctx{T} <: Context
            metadata::T
        end
        $ctx(;metadata=nothing) = $ctx(metadata)
        $ctx{T}(;metadata=nothing) where {T} = $ctx{T}(metadata)
    end
end

function similarcontext(ctx::Ctx; metadata=nothing) where {Ctx<: Context}
    Ctx(metadata=metadata)
end

function default_is_blacklisted(x::GlobalRef, ::Type{<:Context})
    x.mod == Core || x ∈ (
        GlobalRef(Base, :isdispatchtuple),
        GlobalRef(Base, :eltype),
        GlobalRef(Base, :convert),
        GlobalRef(Base, :getproperty),
        GlobalRef(Base, :throw),
    )
end
default_is_blacklisted(::Any, ::Type{<:Context}) = false

is_blacklisted(x, y) = default_is_blacklisted(x, y)

@dynamo function overdub(ctx::C, f, args...) where {C <: Context}
    ir = IR(f, args...)
    ir == nothing && return
    ctx_arg = argument!(ir, at=1)
    for (x, st) in ir
        isexpr(st.expr, :call) || continue
        is_blacklisted(st.expr.args[1], ctx) && continue
        ir[x] = xcall(overdub_pass, ctx_arg, st.expr.args...)
    end
    return ir
end

@inline function overdub_pass(ctx, f, args...)
    prehook(ctx, f, args...)
    res = overdub(ctx, f, args...)
    posthook(ctx, res, f, args...)
    res
end

prehook(args...)  = nothing
@inline overdub(ctx::Context, f::Core.Builtin, args...) = f(args...)
posthook(args...) = nothing

recurse(args...) = invoke(overdub, Tuple{Vararg}, args...)
canrecurse(ctx, args...) = IR(typeof.(args)...) !== nothing

end # module


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

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion <: SwaggerModel
    additionalPrinterColumns::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition} } # spec name: additionalPrinterColumns
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    schema::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation } # spec name: schema
    served::Any # spec type: Union{ Nothing, Bool } # spec name: served
    storage::Any # spec type: Union{ Nothing, Bool } # spec name: storage
    subresources::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresources } # spec name: subresources

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion(;additionalPrinterColumns=nothing, name=nothing, schema=nothing, served=nothing, storage=nothing, subresources=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("additionalPrinterColumns"), additionalPrinterColumns)
        setfield!(o, Symbol("additionalPrinterColumns"), additionalPrinterColumns)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("schema"), schema)
        setfield!(o, Symbol("schema"), schema)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("served"), served)
        setfield!(o, Symbol("served"), served)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("storage"), storage)
        setfield!(o, Symbol("storage"), storage)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("subresources"), subresources)
        setfield!(o, Symbol("subresources"), subresources)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion = Dict{Symbol,Symbol}(Symbol("additionalPrinterColumns")=>Symbol("additionalPrinterColumns"), Symbol("name")=>Symbol("name"), Symbol("schema")=>Symbol("schema"), Symbol("served")=>Symbol("served"), Symbol("storage")=>Symbol("storage"), Symbol("subresources")=>Symbol("subresources"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion = Dict{Symbol,String}(Symbol("additionalPrinterColumns")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition}", Symbol("name")=>"String", Symbol("schema")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation", Symbol("served")=>"Bool", Symbol("storage")=>"Bool", Symbol("subresources")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresources")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion[property_name]

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

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


text:
Dict([
    "opts" => Dict(:T=>10.0, :alg=>TMJets(orderT=5, orderQ=2, abstol=1e-10), :disjointness_method=>BoxEnclosure(), :intersect_source_invariant=>false),
    "X0"   => [(2, Hyperrectangle(low=[4.9, -0.2], high=[5.1, 0.0]))],
    "info" => Dict(
        "name" => "bouncing_ball_nonlinear",
        "dim" => 2,
        "linear" => false,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using ThreadedSparseArrays
using SparseArrays
using LinearAlgebra
using Test

@testset "ThreadedSparseArrays.jl" begin
    N = 4000
    n = 200
    T = ComplexF64

    C = sprand(T, N, n, 0.05)
    @testset "$(Mat)_R" for Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        eye = Matrix(one(T)*I, N, N)
        out = zeros(T, N, n)
        LinearAlgebra.mul!(out, eye, Ct)
        ref = eye*C
        @test norm(ref-out) == 0
    end

    @testset "$(Mat)_L" for Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        eye = Matrix(one(T)*I, n, n)
        out = zeros(T, N, n)
        LinearAlgebra.mul!(out, Ct, eye)
        ref = C*eye
        @test norm(ref-out) == 0
    end

    @testset "$(Mat)_L_$(op)" for op in [adjoint,transpose], Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        eye = Matrix(one(T)*I, N, N)
        out = zeros(T, n, N)
        LinearAlgebra.mul!(out, op(Ct), eye)
        ref = op(C)*eye
        @test norm(ref-out) == 0
    end

    x = rand(0:1,N)
    @testset "$(Mat)_L_$(op)_vec" for op in [adjoint,transpose], Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        out = zeros(T, n)
        LinearAlgebra.mul!(out, op(Ct), x)
        ref = op(C)*x
        @test norm(ref-out) == 0
    end

    # These test below are here to ensure we don't hit ambiguity warnings.
    # The implementations are not (currently) threaded.
    sx = sprand(Bool,n,0.05)
    @testset "$(Mat)_L_sparsevec" for Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        out = Ct*sx
        ref = C*sx
        @test norm(ref-out) == 0
        @test typeof(ref)==typeof(out)
    end

    sx = sprand(Bool,N,0.05)
    @testset "$(Mat)_L_$(op)_sparsevec" for op in [adjoint,transpose], Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        out = op(Ct)*sx
        ref = op(C)*sx
        @test norm(ref-out) == 0
        @test typeof(ref)==typeof(out)
    end

    sx = sparse(rand(1:n,10),1:10,true,n,10)
    @testset "$(Mat)_L_sparse" for Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        out = Ct*sx
        ref = C*sx
        @test norm(ref-out) == 0
        @test typeof(ref)==typeof(out)
    end

    sx = sparse(rand(1:N,10),1:10,true,N,10)
    @testset "$(Mat)_L_$(op)_sparse" for op in [adjoint,transpose], Mat in [ThreadedSparseMatrixCSC]
        Ct = Mat(C)

        out = op(Ct)*sx
        ref = op(C)*sx
        @test norm(ref-out) == 0
        @test typeof(ref)==typeof(out)
    end

end


text:
x, y = 1, 2
z = x + y


text:
function BLPConstant(handle::Ptr{Cvoid})
    ptr_check(handle, "Failed to create BLPConstant")
    name = BLPName(blpapi_Constant_name(handle))
    description = unsafe_string(blpapi_Constant_description(handle))
    status = SchemaStatus(blpapi_Constant_status(handle))
    datatype = BLPDataType(blpapi_Constant_datatype(handle))
    value = parse_blpapi_const_value(handle, datatype)
    return BLPConstant{datatype, typeof(value)}(name, description, status, datatype, value)
end

function BLPConstantList(handle::Ptr{Cvoid})
    ptr_check(handle, "Failed to create BLPConstantList")
    name = BLPName(blpapi_ConstantList_name(handle))
    description = unsafe_string(blpapi_ConstantList_description(handle))
    status = SchemaStatus(blpapi_ConstantList_status(handle))
    datatype = BLPDataType(blpapi_ConstantList_datatype(handle))

    num_constants = blpapi_ConstantList_numConstants(handle)
    if num_constants == 0
        list = Vector{BLPConstant}()
    else
        @assert num_constants > 0
        list = [ BLPConstant(blpapi_ConstantList_getConstantAt(handle, i)) for i in 0:(num_constants-1) ]
    end

    return BLPConstantList{datatype}(name, description, status, datatype, list)
end

function parse_blpapi_const_value(handle::Ptr{Cvoid}, datatype::BLPDataType)

    _error_check(code) = error_check(code, "Error parsing $datatype")

    if datatype == BLPAPI_DATATYPE_CHAR
        buffer = Ref{Cchar}()
        err = blpapi_Constant_getValueAsChar(handle, buffer)
        _error_check(err)
        return Char(buffer[])

    elseif datatype == BLPAPI_DATATYPE_INT32
        buffer = Ref{Int32}()
        err = blpapi_Constant_getValueAsInt32(handle, buffer)
        _error_check(err)
        return buffer[]

    elseif datatype == BLPAPI_DATATYPE_INT64
        buffer = Ref{Int64}()
        err = blpapi_Constant_getValueAsInt64(handle, buffer)
        _error_check(err)
        return buffer[]

    elseif datatype == BLPAPI_DATATYPE_FLOAT32
        buffer = Ref{Float32}()
        err = blpapi_Constant_getValueAsFloat32(handle, buffer)
        _error_check(err)
        return buffer[]

    elseif datatype == BLPAPI_DATATYPE_FLOAT64
        buffer = Ref{Float64}()
        err = blpapi_Constant_getValueAsFloat64(handle, buffer)
        _error_check(err)
        return buffer[]

    elseif datatype == BLPAPI_DATATYPE_DATETIME
        buffer = Ref{BLPDateTime}()
        err = blpapi_Constant_getValueAsDatetime(handle, buffer)
        _error_check(err)
        return buffer[]

    elseif datatype == BLPAPI_DATATYPE_STRING
        buffer = Ref{Cstring}()
        err = blpapi_Constant_getValueAsString(handle, buffer)
        return unsafe_string(buffer[])
    end

    error("Unsupported datatype $datatype.")
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
function State2XY(DP::DrivePOMDP, st::CarSt)
    Pos = GetGlobalPosition(DP.Routes[st.r], st.s)
    return (Pos.x, Pos.y)
end

function Junction(halflen::Float64, gauge::Float64, θ::Float64; style::Symbol) # gauge = Spurweite
    Ycenter=Matrix{Float64}
    if style == :Crossroad
        plot(xlim=(-halflen, halflen), ylim=(-halflen*sin(θ), halflen), size=(650, 650*sin(θ)), grid=false, legend=false, showaxis=false)
        Ycenter = [0.0 0.0; -halflen halflen]
    elseif style == :TJunction
        plot(xlim=(-halflen, halflen), ylim=(-0.3*halflen, halflen), size=(650, 650/2*1.3), grid=false, legend=false, showaxis=false)
        plot!([-halflen, halflen], [-gauge, -gauge], c=:black)
        Ycenter = [0.0 0.0; 0.0 halflen]
    else
        error("style $style not supported, please choose :Crossroad or :TJunction and try again.")
    end
    plot!([-halflen, halflen], [0.0, 0.0], c=:lightgray, linestyle=:dash) # x axis; center line
    α = π/2 - θ
    YcenterR = Rmat(α)*Ycenter
    plot!(YcenterR[1,:], YcenterR[2,:], c=:lightgray, linestyle=:dash) # y axis; center line

    rr = 1.0
    # circular fillet left
    pts_c1 = Plots.partialcircle(-π/2, π/2-θ, 100, rr) # n=100 steps
    x_c1, y_c1 = Plots.unzip(pts_c1)
    x_c1 .-= gauge*(tan(α)+1/sin(θ)) + rr/tan(θ/2)
    y_c1 .+= gauge + rr
    pushfirst!(x_c1, -halflen)
    pushfirst!(y_c1, gauge)
    push!(x_c1, -gauge/sin(θ)-halflen*cos(θ))
    push!(y_c1, halflen*sin(θ))
    plot!(x_c1, y_c1, c=:black)

    # circular fillet right
    pts_c2 = Plots.partialcircle(-π/2-θ, -π/2, 100, rr) # n=100 steps
    x_c2, y_c2 = Plots.unzip(pts_c2)
    x_c2 .+= -(gauge+rr*(1-cos(θ)))/tan(θ) + gauge/sin(θ) + rr*sin(θ)
    y_c2 .+= gauge+rr
    pushfirst!(x_c2, -halflen*cos(θ)+gauge/sin(θ))
    pushfirst!(y_c2, halflen*sin(θ))
    push!(x_c2, halflen)
    push!(y_c2, gauge)
    plt = plot!(x_c2, y_c2, c=:black)

    if style == :Crossroad
        pts_c1_mirror = Rmat(π)*[transpose(x_c1); transpose(y_c1)]
        plt = plot!(plt, pts_c1_mirror[1,:], pts_c1_mirror[2,:], c=:black)
        pts_c2_mirror = Rmat(π)*[transpose(x_c2); transpose(y_c2)]
        plt = plot!(plt, pts_c2_mirror[1,:], pts_c2_mirror[2,:], c=:black)
    end

    return plt
end

function DrawStopline!(plt::Plots.Plot{Plots.GRBackend}, R::Route, StoplineLength::Float64, gauge::Float64)
    StopPoint, θ = GetGlobalPosAngle(R, StoplineLength)
    Δx = gauge/2 * sin(θ)
    Δy = gauge/2 * cos(θ)
    X = [StopPoint.x - Δx, StopPoint.x + Δx]
    Y = [StopPoint.y + Δy, StopPoint.y - Δy]
    plot!(plt, X, Y, c=:blue, linewidth=4)
end

function PlotRoutes(DP::DrivePOMDP, plt::Plots.Plot{Plots.GRBackend})
    for i in 1:length(DP.Routes)
        if i == 1
            for G in DP.Routes[i].Geos
                PlotGeometry(G, plt; color=:blue)
            end
        else
            for G in DP.Routes[i].Geos
                PlotGeometry(G, plt)
            end
        end
    end
    return plt
end

function PlotGeometry(Geo::Geometry, plt::Plots.Plot{Plots.GRBackend}; color::Symbol=:red)
    if Geo.GType == "arc"
        if Geo.Ep.θ > Geo.Sp.θ
            δ = -π/2
        else
            δ = π/2
        end
        pts = Plots.partialcircle(Geo.Sp.θ+δ, Geo.Ep.θ+δ, 100, abs(1/Geo.Sp.cur))
        Center = ArcCenter(Geo)
        x, y = Plots.unzip(pts)
        x .+= Center.x
        y .+= Center.y
        plot!(plt, x, y, c=color, linestyle=:dot)
    elseif Geo.GType == "line"
        plot!(plt, [Geo.Sp.x, Geo.Ep.x], [Geo.Sp.y, Geo.Ep.y], c=color, linestyle=:dot)
    else
        error("Can not print this type of geomatry: $(Geo.GType)")
    end
end

function SVAInit(DP::DrivePOMDP, StsArray::Vector{Sts})
    t = (length(StsArray)-1)*DP.Δt
    Splt = plot(xlim=(0.0, t*1.2), ylim=(0.0, DP.sV[end]*1.2), size=(200, 150), grid=true, legend=false, ylabel="S(m)", left_margin=-2mm, titleposition=:left, showaxis=true, xlabel="t/s", xticks =0:1.0:UInt16(floor(t*1.2)), yticks =0:4*DP.Δs:DP.sV[end])
    Vplt = plot(xlim=(0.0, t*1.2), ylim=(0.0, DP.vV[end]*1.2), size=(200, 150), grid=true, legend=false, ylabel="V(m/s)", left_margin=-2mm, titleposition=:left, showaxis=true, xlabel="t/s", xticks =0:1.0:UInt16(floor(t*1.2)), yticks =DP.vV)
    Aplt = plot(xlim=(0.0, t*1.2), ylim=(DP.Aset.min-0.5, DP.Aset.max+0.5), size=(200, 150), grid=true, legend=false, ylabel="A(m/s²)", left_margin=-2mm, titleposition=:left, showaxis=true, xlabel="t/s", xticks =0:1.0:UInt16(floor(t*1.2)), yticks =DP.Aset.min:1.0:DP.Aset.max)
    return Splt, Vplt, Aplt
end

function PlotSVA!(Plts::NTuple{3, Plots.Plot{Plots.GRBackend}}, DP::DrivePOMDP, Step::Int64, StPair::NTuple{2, CarSt}, AvecPair::NTuple{2, Vector{Float64}}) # should plot standalone and the returned plt should be used as subplot
    plot!(Plts[1], [(Step-1)*DP.Δt, Step*DP.Δt], [StPair[1].s, StPair[2].s], c=:black)
    plot!(Plts[2], [(Step-1)*DP.Δt, Step*DP.Δt], [StPair[1].v, StPair[2].v], c=:black)
    #plot!(Plts[3], [(Step-1)*DP.Δt, (Step-1)*DP.Δt], [APair[1], APair[2]], c=:red, linestyle = :dot)
    Astep = Int64(DP.Δt/0.1)
    for i in 1:Astep+1
        if Step > 1 && i == 1
            plot!(Plts[3], [(Step-1)*DP.Δt+0.1*(i-1), (Step-1)*DP.Δt+0.1*(i-1)], [AvecPair[1][end], AvecPair[2][i]], c=:black)
        elseif i >= 2
            plot!(Plts[3], [(Step-1)*DP.Δt+0.1*(i-2), (Step-1)*DP.Δt+0.1*(i-1)], [AvecPair[2][i-1], AvecPair[2][i]], c=:black)
        end
    end
    return nothing
end

function Visualisation(DP::DrivePOMDP, StsArray::Vector{Sts}, ObsArray::Vector{CarOb}, ActVec::Vector{Symbol}, AccMat::Vector{Vector{Float64}}, halflen::Float64, gauge::Float64, θ::Float64, style::Symbol, picname::String)
    clibrary(:Plots)
    pltJunction = Junction(halflen, gauge, θ; style=style)
    pltMap = PlotRoutes(DP, pltJunction)
    DrawStopline!(pltMap, DP.Routes[DP.SsInit.Ego.r], DP.Stopline, gauge)
    pltSVA_Ego = SVAInit(DP, StsArray)
    plt_Ego = plot(pltSVA_Ego..., layout=(3,1), size=(250, 650))

    anim = @animate for i in 1:length(StsArray)
        if i > 1
            plot!(pltMap, [ObsArray[i-1].x], [ObsArray[i-1].y], seriestype=:scatter, markersize = 4, c=:violet, markerstrokecolor=:black)
        end
        if i == 2
            StPair = (StsArray[i-1].Ego, StsArray[i].Ego)
            PlotSVA!(pltSVA_Ego, DP, i-1, StPair, (AccMat[i-1], AccMat[i-1]))
        end
        if i > 2
            StPair = (StsArray[i-1].Ego, StsArray[i].Ego)
            PlotSVA!(pltSVA_Ego, DP, i-1, StPair, (AccMat[i-2], AccMat[i-1]))
        end

        if ActVec[i] == :takeover
            plot!(pltMap, State2XY(DP, StsArray[i].Ego), seriestype=:scatter, markersize = 5, c=:green, markerstrokecolor=:black)
        else
            plot!(pltMap, State2XY(DP, StsArray[i].Ego), seriestype=:scatter, markersize = 5, c=:orange, markerstrokecolor=:black)
        end
        plot!(pltMap, State2XY(DP, StsArray[i].Other), seriestype=:scatter, markersize = 5, c=:red, markerstrokecolor=:black)

        if style == :TJunction
            column1 = @layout [a{0.2h}; b{0.6h}; c{0.2h}]
            plt_null = plot(grid=false, legend=false, showaxis=false)
            plt_col1 = plot(plt_null, pltMap, plt_null, layout=column1)
            row1 = @layout [a{0.7w} b{0.3w}]
            plot(plt_col1, plt_Ego, layout=row1, size=(900, 650))
        elseif style == :Crossroad
            plot(pltMap, plt_Ego, layout=grid(1, 2, widths=[0.7, 0.3]), size=(900, 650))
        else
            error("Unsupported Style: $style. Please choose either :Crossroad or :TJunction.")
        end
        savefig("output/$picname.frame$i.pdf")
    end
    gif(anim, "output/$picname.gif", fps = Int(floor(1/DP.Δt)))
end


text:
using Documenter, FunManifolds, Manifolds, ManifoldsBase

#cd("docs")
makedocs(
    modules=[FunManifolds, Manifolds, ManifoldsBase],
    authors="Mateusz Baran",
    format=Documenter.HTML(prettyurls=false),
    sitename="FunManifolds.jl",
    pages=[
        "Home" => "index.md",
        "Manual" => Any["man/examples.md", "man/manifolds.md"],
        "Library" => Any["lib/public.md", "lib/internals.md"],
    ],
)


text:
import GR

j = [0.19454, 0.1194, 0.19454, 0.41129, 0.10807, 0.3875, 0.10807, 0.07711, 0.10807, 0.06075, 0.10681, 0.04973, 0.1043, 0.04407, 0.10178, 0.0384, 0.09813, 0.03557, 0.09335, 0.03557, 0.09108, 0.03557, 0.08856, 0.03645, 0.08579, 0.03822, 0.08303, 0.03998, 0.07963, 0.0435, 0.0756, 0.04879, 0.07207, 0.05357, 0.06761, 0.05842, 0.06219, 0.06333, 0.05678, 0.06824, 0.04954, 0.07069, 0.04048, 0.07069, 0.0284, 0.07069, 0.01877, 0.06767, 0.0116, 0.06163, 0.00442, 0.05559, 0.00083, 0.04816, 0.00083, 0.03935, 0.00083, 0.02878, 0.00725, 0.01971, 0.02009, 0.01216, 0.03293, 0.00461, 0.05143, 0.00083, 0.0756, 0.00083, 0.09347, 0.00083, 0.10971, 0.00228, 0.12431, 0.00518, 0.13891, 0.00807, 0.15143, 0.0138, 0.16188, 0.02236, 0.17233, 0.03092, 0.18038, 0.04294, 0.18605, 0.05842, 0.19171, 0.0739, 0.19454, 0.09423, 0.19454, 0.1194]

u = [0.31464, 0.40298, 0.22855, 0.40298, 0.22855, 0.20323, 0.22855, 0.1909, 0.23113, 0.17932, 0.23629, 0.16849, 0.24145, 0.15767, 0.24856, 0.14823, 0.25763, 0.14017, 0.26669, 0.13212, 0.27726, 0.12576, 0.28935, 0.1211, 0.30143, 0.11644, 0.31452, 0.11412, 0.32862, 0.11412, 0.3407, 0.11412, 0.35316, 0.11701, 0.366, 0.1228, 0.37884, 0.12859, 0.39092, 0.13627, 0.40225, 0.14583, 0.40225, 0.1194, 0.48834, 0.1194, 0.48834, 0.40298, 0.40225, 0.40298, 0.40225, 0.19832, 0.3957, 0.18825, 0.38803, 0.17988, 0.37922, 0.17321, 0.3704, 0.16654, 0.36273, 0.1632, 0.35618, 0.1632, 0.35039, 0.1632, 0.34498, 0.16427, 0.33994, 0.16641, 0.33491, 0.16855, 0.3305, 0.17139, 0.32673, 0.17491, 0.32295, 0.17843, 0.31999, 0.18265, 0.31785, 0.18756, 0.31571, 0.19247, 0.31464, 0.19769, 0.31464, 0.20323]

l = [0.60797, 0.1194, 0.60797, 0.54081, 0.52225, 0.51702, 0.52225, 0.1194]

i = [0.64226, 0.3875, 0.64226, 0.1194, 0.72835, 0.1194, 0.72835, 0.41129]

a = [0.91001, 0.26591, 0.91001, 0.18397, 0.9012, 0.17743, 0.89321, 0.17201, 0.88604, 0.16773, 0.87886, 0.16345, 0.87175, 0.16132, 0.8647, 0.16132, 0.86118, 0.16132, 0.8579, 0.16238, 0.85488, 0.16452, 0.85186, 0.16666, 0.84909, 0.1695, 0.84658, 0.17302, 0.84406, 0.17655, 0.84211, 0.18076, 0.84072, 0.18567, 0.83934, 0.19058, 0.83865, 0.19568, 0.83865, 0.20096, 0.83865, 0.20776, 0.84079, 0.21431, 0.84507, 0.2206, 0.84934, 0.22689, 0.85495, 0.23281, 0.86187, 0.23835, 0.86879, 0.24388, 0.87647, 0.24898, 0.8849, 0.25364, 0.89334, 0.2583, 0.90171, 0.26239, 0.91001, 0.26591, 0.99648, 0.1194, 0.99648, 0.33199, 0.99648, 0.34382, 0.99422, 0.35452, 0.98969, 0.36409, 0.98516, 0.37365, 0.97811, 0.38184, 0.96854, 0.38863, 0.95898, 0.39543, 0.94683, 0.40065, 0.9321, 0.4043, 0.91738, 0.40795, 0.89982, 0.40978, 0.87943, 0.40978, 0.86281, 0.40978, 0.84708, 0.40802, 0.83223, 0.40449, 0.81737, 0.40097, 0.80428, 0.39618, 0.79296, 0.39014, 0.78163, 0.3841, 0.77263, 0.37686, 0.76596, 0.36843, 0.75929, 0.36, 0.75595, 0.35087, 0.75595, 0.34105, 0.75595, 0.33048, 0.75973, 0.32173, 0.76728, 0.31481, 0.77483, 0.30789, 0.78465, 0.30443, 0.79673, 0.30443, 0.80454, 0.30443, 0.81102, 0.30556, 0.81618, 0.30783, 0.82134, 0.31009, 0.8253, 0.31311, 0.82807, 0.31689, 0.83084, 0.32066, 0.83279, 0.32507, 0.83393, 0.3301, 0.83506, 0.33514, 0.83563, 0.3403, 0.83563, 0.34559, 0.83563, 0.3549, 0.83827, 0.36283, 0.84356, 0.36937, 0.84884, 0.37592, 0.85828, 0.37919, 0.87188, 0.37919, 0.88346, 0.37919, 0.89271, 0.37542, 0.89963, 0.36786, 0.90655, 0.36031, 0.91001, 0.34747, 0.91001, 0.32935, 0.91001, 0.30141, 0.90057, 0.29914, 0.88899, 0.29562, 0.87729, 0.2919, 0.86546, 0.288, 0.85362, 0.2841, 0.84223, 0.27982, 0.83128, 0.27516, 0.82033, 0.27051, 0.81007, 0.26541, 0.80051, 0.25987, 0.79094, 0.25433, 0.78257, 0.24804, 0.7754, 0.24099, 0.76822, 0.23394, 0.76256, 0.22601, 0.75841, 0.2172, 0.75425, 0.20839, 0.75217, 0.19857, 0.75217, 0.18775, 0.75217, 0.17717, 0.75413, 0.16742, 0.75803, 0.15848, 0.76193, 0.14955, 0.76753, 0.14181, 0.77483, 0.13526, 0.78213, 0.12872, 0.791, 0.12355, 0.80145, 0.11978, 0.8119, 0.116, 0.82367, 0.11411, 0.83676, 0.11411, 0.84632, 0.11411, 0.85463, 0.11481, 0.86168, 0.11619, 0.86873, 0.11758, 0.87508, 0.11953, 0.88075, 0.12204, 0.88641, 0.12456, 0.89151, 0.12752, 0.89604, 0.13092, 0.90057, 0.13432, 0.90523, 0.13803, 0.91001, 0.14206, 0.91001, 0.1194]

function draw_path(p, codes, fill, border)
  GR.setfillcolorind(fill)
  GR.setbordercolorind(border)
  op = border == 0 ? 'f' : 'F'
  GR.path(p[1:2:end], p[2:2:end], codes)
end

function circle(xmin, ymin, xmax, ymax, fill, border)
  GR.setfillintstyle(4)
  GR.setfillcolorind(fill)
  GR.setbordercolorind(border)
  GR.fillarc(xmin, xmax, ymin, ymax, 0, 360)
end

GR.setcolorrep(1, 0.14, 0.14, 0.14)
GR.setcolorrep(2, 0.251, 0.388, 0.847) # darker blue
GR.setcolorrep(3, 0.796, 0.235, 0.2)   # darker red
GR.setcolorrep(4, 0.584, 0.345, 0.698) # darker purple
GR.setcolorrep(5, 0.22, 0.596, 0.149)  # darker green
GR.setcolorrep(6, 0.4, 0.51, 0.878)    # lighter blue
GR.setcolorrep(7, 0.835, 0.388, 0.361) # lighter red
GR.setcolorrep(8, 0.667, 0.475, 0.757) # lighter purple
GR.setcolorrep(9, 0.376, 0.678, 0.318) # lighter green

GR.setviewport(0, 1, 0, 1)
GR.setwindow(0, 1, 0, 1)
GR.updatews()

for s in 0.1:0.01:0.5
  GR.clearws()
  GR.setviewport(0.5 - s, 0.5 + s, 0.5 - s, 0.5 + s)
  GR.setborderwidth(s * 8)

  draw_path(j, "MLLLCCCCCCCCCCCCCCf", 1, 0)
  draw_path(u, "MLLCCCCCCLLLLLCCCCCCf", 1, 0)
  draw_path(l, "MLLLf", 1, 0)
  draw_path(i, "MLLLf", 1, 0)
  draw_path(a, "MLCCCCCCCCCCZMLCCCCCCCCCCCCCCCCCCLLCCCCCCCCCCCCCCLf", 1, 0)

  circle(0.09698, 0.43274, 0.20364, 0.5394, 2, 6)
  circle(0.63794, 0.43274, 0.7446, 0.5394, 3, 7)
  circle(0.70794, 0.55274, 0.8146, 0.6594, 4, 8)
  circle(0.77794, 0.43274, 0.8846, 0.5394, 5, 9)

  GR.updatews()

  sleep(0.01)
end