text:
using BytePairEncoding, JSON, HTTP
using Flux
using StatsBase
using Transformers
using Transformers.Basic
using Transformers.GenerativePreTrain
using Transformers.Pretrain
using Transformers.HuggingFace

isdir("model_vocab")||mkdir("model_vocab")

# This is a temporary fix, will be updated once Tokenizer API is ready
isfile("model_vocab/bpe.out") || HTTP.download("https://s3.amazonaws.com/models.huggingface.co/bert/gpt2-merges.txt", "./model_vocab/bpe.out")
isfile("model_vocab/vocab.json") || HTTP.download("https://s3.amazonaws.com/models.huggingface.co/bert/gpt2-vocab.json", "./model_vocab/vocab.json")

labels = map(x->x[1], sort!(collect(pairs(JSON.parsefile("./model_vocab/vocab.json"))), by=x->x[2]))
encoder = Vocabulary(labels, "<|endoftext|>")
bpe = ByteLevelBPE("./model_vocab/bpe.out")

model = hgf"gpt2:lmheadmodel"

function encode_text(text)
  xs = ["<|endoftext|>"; bpe(text)]
  return encoder(xs)
end

function temp_softmax(logits; temperature=1.2)
  return softmax(logits ./ temperature)
end

function top_k_sample(probs; k=10)
  sorted = sort(probs, rev = true)
  indexes = partialsortperm(probs, 1:k, rev=true)
  index = sample(indexes, ProbabilityWeights(sorted[1:k]), 1)
  return index
end

function generate_text(;context="", max_length=50)
  input_ = encode_text(context)
  for i in 1:max_length
    input_ids = reshape(@view(input_[:]), :, 1)
    outputs = model(input_ids; output_attentions=false,
                    output_hidden_states=false,
                    use_cache=false)
    logits = @view outputs.logits[:, end, 1]
    probs = temp_softmax(logits)
    new_token = top_k_sample(probs)[1]
    push!(input_, new_token)
  end
  return input_
end

function decode_text(text_token_ids)
  text_list = decode(encoder, text_token_ids)
  text = BytePairEncoding.UnMap(bpe.codemap)(join(text_list))
  return text
end

function generate(prompt, max_length)
  text_token_ids = generate_text(context = prompt; max_length=max_length)
  gen_text = decode_text(text_token_ids)
  print("\n\nGenerated Text: ")
  println(gen_text)
end

generate( "Fruits are very good for ", 100)


text:
# # Arithmetical terms for rational numbers
# Parsing is reading and transforming a sequence of characters.
# This example reads and evaluates arithmetical terms for rational numbers.
# Subterms can use algebraic operators `+-*/` that will be evaluated with
function evaluate( (start, operation_values) )
    aggregated_value::Rational{Int} = start
    for (op,val) in operation_values
        aggregated_value = eval( Expr(:call, Symbol(op),
			              aggregated_value, val
			              ))
    end
    return aggregated_value
end
evaluate( (0, [ ('+',1), ('+',1) ]) )
evaluate( (1, [ ('*',2), ('*',3) ]) )
#
using CombinedParsers
# `CombinedParsers` provides constructors to combine parsers and transform (sub-)parsings arbitrarily with julia syntax.
# Combinator constructors are discussed in the [user guide](user.md).
using TextParse
# Term expressions are sequences of subterms interleaved with operators.
# Sub terms are [`Either`](@ref) fast `TextParse.Numeric(Int)` integer numbers, converted to `Rational{Int}`,
@syntax subterm = Either{Rational{Int}}(Any[TextParse.Numeric(Int)]);
# A subterm can also be a nested term in parentheses
@syntax for parenthesis in subterm
    mult = evaluate |> join(subterm, CharIn("*/"), infix=:prefix )
    adds = evaluate |> join(mult,    CharIn("+-"), infix=:prefix )
    Sequence(2,'(',adds,')')
end;
# This `CombinedParser` definition in 5,5 lines is sufficient for doing arithmetics:
# [`Base.join`](@ref)(x,infix; infix=:prefix) is shorthand for `x `[`(*)`](@ref)` `[`Repeat`](@ref)`( infix * x  )`,
# and `f |> parser` is shorthand for [`map`](@ref)(f,parser)`.
@syntax term = adds;
# registers a `@term_string` macro for parsing and transforming.
term"(1+2)/5"

# The defined `CombinedParser` `term` can be used as a function for colorful logging of the parsing process.
term("1/((1+2)*4+3*(5*2))",log = [:parenthesis])
# [Is every rational answer ultimately the inverse of a universal question in life?](https://en.wikipedia.org/wiki/Phrases_from_The_Hitchhiker%27s_Guide_to_the_Galaxy#Answer_to_the_Ultimate_Question_of_Life,_the_Universe,_and_Everything_(42))
#
# The parser representation can be printed as a tree
term

# ### Benchmarks
# Parsing times for Int, operators, brackets are
using BenchmarkTools
@benchmark match(term,"(1+2)/5")

# in unfair benchmark-comparison with the more expressive Julia syntax parser
@benchmark Meta.parse("(1+2)/5")

# Parsing and transforming (here `eval`)
@benchmark term("(1+2)/5")

# compared to Julia
@benchmark eval(Meta.parse("(1+2)/5"))

ext: jl


text:
# walk.jl
# implement walk method to locate a point

# TODO: Walk function
# INPUT: index of starting tet, point to locate, Array of points, Array of tets
# OUTPUT: index of tet containing point
#
# While we haven't found the containing tet (new index does not match last index)
#   update the current tet with STEP or the current index
#

# TODO: Step function
# INPUT: Current tet, points of current tet, point to locate
# OUTPUT: Index of next tet
#
# Iterate through the faces
#   if orient of the face and the target point is negative
#       return the neighbor index associated with that face
# if it falls through, return nothing
#


text:
function _JBox(
    cornerpoint1::Point,
    cornerpoint2::Point,
    color,
    action::Symbol,
    vertices::Bool,
)
    sethue(color)
    verts = box(cornerpoint1, cornerpoint2, action, vertices = vertices)
    return verts[2]
end
function _JBox(points::Array, color, action::Symbol, vertices::Bool)
    sethue(color)
    verts = box(points, action, vertices = vertices)
    return verts[2]
end
function _JBox(pt::Point, width::Real, height::Real, color, action::Symbol, vertices::Bool)
    sethue(color)
    box(pt, width, height, action, vertices = vertices)
    return Point(pt.x - width / 2, pt.y + height / 2)
end
function _JBox(
    pt::Point,
    width::Real,
    height::Real,
    cornerradius::Float64,
    color,
    action::Symbol,
)
    sethue(color)
    box(pt, width, height, cornerradius, action)
    return Point(pt.x - width / 2, pt.y + height / 2)
end

"""
    JBox(cornerpoint1::Point, cornerpoint2::Point; color="black", action=:stroke, vertices=false)

Create a box (rectangle) between two points and do an action.
Returns the top left corner point of the box.
"""
JBox(
    cornerpoint1::Point,
    cornerpoint2::Point;
    color = "black",
    action = :stroke,
    vertices = false,
) =
    (
        args...;
        cornerpoint1 = cornerpoint1,
        cornerpoint2 = cornerpoint2,
        color = color,
        action = action,
        vertices = vertices,
    ) -> _JBox(cornerpoint1, cornerpoint2, color, action, vertices)

"""
    JBox(points::Array; color="black", action=:stroke, vertices=false)

Create a box/rectangle using the first two points of an array of Points to defined opposite corners.
Returns the top left corner point of the box.
"""
JBox(points::Array; color = "black", action = :stroke, vertices = false) =
    (args...; points = points, color = color, action = action, vertices = vertices) ->
        _JBox(points, color, action, vertices)

"""
    JBox(pt::Point, width::Real, height::Real; color="black", action=:stroke, vertices=false)

Create a box/rectangle centered at point pt with width and height. Use vertices=true to return an array of the four corner points rather than draw the box.
"""
JBox(
    pt::Point,
    width::Real,
    height::Real;
    color = "black",
    action = :stroke,
    vertices = false,
) =
    (
        args...;
        pt = pt,
        width = width,
        height = height,
        color = color,
        action = action,
        vertices = vertices,
    ) -> _JBox(pt, width, height, color, action, vertices)

"""
    JBox(x::Int64, y::Int64, width::Real, height::Real; color="black", action=:stroke)

Create a box/rectangle centered at point x/y with width and height.
"""
JBox(x::Int64, y::Int64, width::Real, height::Real; color = "black", action = :stroke) =
    JBox(Point(x, y), width, height, color = color, action = action)

"""
    JBox(pt::Point, width::Real, height::Real, cornerradius::Float64; color="black", action=:stroke)

Draw a box/rectangle centered at point pt with width and height and round each corner by cornerradius.
"""
JBox(
    pt::Point,
    width::Real,
    height::Real,
    cornerradius::Float64;
    color = "black",
    action = :stroke,
) =
    (
        args...;
        pt = pt,
        width = width,
        height = height,
        cornerradius = cornerradius,
        color = color,
        action = action,
    ) -> _JBox(pt, width, height, cornerradius, color, action)


text:
export ft_getdeviceinfolist
export FT_FLAGS_OPENED

include("FT_DEVICE.jl")  # load constants

const handleis32 = Culong == UInt32

type FtDeviceListInfoNode
  flags         :: UInt32
  devicetype    :: UInt32
  id            :: UInt32
  locid         :: UInt32
  serialnumber  :: ASCIIString
  description   :: ASCIIString
  handle        :: Culong
end

function Base.show(io::IO, i::FtDeviceListInfoNode)
  println(typeof(i),":")
  f()=@printf(io,"flags = 0x%08x \n",i.flags);f()
  f()=@printf(io,"devicetype = 0x%08x \n",i.devicetype);f()
  f()=@printf(io,"id = 0x%08x \n",i.id);f()
  f()=@printf(io,"locid = 0x%08x \n",i.locid);f()
  println(io,"serialnumber = ",i.serialnumber)
  println(io,"description = ",i.description)
  if handleis32
    f()=@printf(io,"handle = 0x%08x \n",i.handle);f()
  else
    f()=@printf(io,"handle = 0x%016x \n",i.handle);f()
  end
end

function Base.showcompact(io::IO, i::FtDeviceListInfoNode)
  println(io,"sn=",rpad(i.serialnumber,17,' ')," description=",rpad(i.description,40,' '))
  if handleis32
    f()=@printf(io,"f=0x%08x, dt=0x%08x, id=0x%08x, lid=0x%08x, h=0x%08x \n",i.flags, i.devicetype, i.id, i.locid, i.handle);f()
  else
    f()=@printf(io,"f=0x%08x, dt=0x%08x, id=0x%08x, lid=0x%08x, h=0x%016x \n",i.flags, i.devicetype, i.id, i.locid, i.handle);f()
  end
end

function Base.show(io::IO, a::Array{FTD2XX.FtDeviceListInfoNode,1})
  println(io,length(a),"-element ",typeof(a),":")
  for i in a
    showcompact(io,i)
  end
end

# FT_FLAGS (see FT_DEVICE_LIST_INFO_NODE)
const FT_FLAGS_OPENED = 0x00000001

# FT_DEVICE_LIST_INFO_NODE (see ft_getdeviceinfolist and ft_getdeviceinfoDetail)
immutable FtDeviceListFtDeviceListInfoNode_C
  Flags :: Cuint
  Type  :: Cuint
  ID    :: Cuint
  LocId :: Cuint
  sn1  :: UInt8; sn2  :: UInt8; sn3  :: UInt8; sn4  :: UInt8  # SerialNumber
  sn5  :: UInt8; sn6  :: UInt8; sn7  :: UInt8; sn8  :: UInt8
  sn9  :: UInt8; sn10 :: UInt8; sn11 :: UInt8; sn12 :: UInt8
  sn13 :: UInt8; sn14 :: UInt8; sn15 :: UInt8; sn16 :: UInt8
  d1   :: UInt8; d2   :: UInt8; d3   :: UInt8; d4   :: UInt8  # Description
  d5   :: UInt8; d6   :: UInt8; d7   :: UInt8; d8   :: UInt8
  d9   :: UInt8; d10  :: UInt8; d11  :: UInt8; d12  :: UInt8
  d13  :: UInt8; d14  :: UInt8; d15  :: UInt8; d16  :: UInt8
  d17  :: UInt8; d18  :: UInt8; d19  :: UInt8; d20  :: UInt8
  d21  :: UInt8; d22  :: UInt8; d23  :: UInt8; d24  :: UInt8
  d25  :: UInt8; d26  :: UInt8; d27  :: UInt8; d28  :: UInt8
  d29  :: UInt8; d30  :: UInt8; d31  :: UInt8; d32  :: UInt8
  d33  :: UInt8; d34  :: UInt8; d35  :: UInt8; d36  :: UInt8
  d37  :: UInt8; d38  :: UInt8; d39  :: UInt8; d40  :: UInt8
  d41  :: UInt8; d42  :: UInt8; d43  :: UInt8; d44  :: UInt8
  d45  :: UInt8; d46  :: UInt8; d47  :: UInt8; d48  :: UInt8
  d49  :: UInt8; d50  :: UInt8; d51  :: UInt8; d52  :: UInt8
  d53  :: UInt8; d54  :: UInt8; d55  :: UInt8; d56  :: UInt8
  d57  :: UInt8; d58  :: UInt8; d59  :: UInt8; d60  :: UInt8
  d61  :: UInt8; d62  :: UInt8; d63  :: UInt8; d64  :: UInt8
  FT_HANDLE     :: Cuint
  FT_HANDLE_MSB :: Cuint
end

function ft_getdeviceinfolist(lpdwNumDevs::Integer = ft_createdeviceinfolist())
  ftdeviceinfolist = Array(FtDeviceListFtDeviceListInfoNode_C,lpdwNumDevs)
  n = Ref{Cuint}(lpdwNumDevs)
  ft_status = ccall((:FT_GetDeviceInfoList, d2xx),
                      Cuint,
                      (Ptr{FtDeviceListFtDeviceListInfoNode_C},Ref{Cuint}),
                      ftdeviceinfolist,n)
  checkstatus(ft_status)
  infonodearray = FtDeviceListInfoNode[]
  sizehint!(infonodearray,lpdwNumDevs)
  for node in ftdeviceinfolist
    flags = node.Flags
    devicetype = node.Type
    id = node.ID
    locid = node.LocId
    sn = Array(UInt8,16)
    for i in 1:16 # 5:20
      sn[i] = node.(i+4)
    end
    endofstring = findfirst(sn,0)-1
    serialnumber = convert(ASCIIString,sn[1:endofstring])
    d = Array(UInt8,64)
    for i in 1:64 # 21:85
      d[i] = node.(i+20)
    end
    endofstring = findfirst(d,0)-1
    description = convert(ASCIIString,d[1:endofstring])
    # this seems wrong, but I don't have a better idea.
    if Culong == UInt32
      handle = node.FT_HANDLE
    else
      handle = UInt64(node.FT_HANDLE)
      handle = handle | (node.FT_HANDLE_MSB <<32)
    end
    push!(infonodearray,FtDeviceListInfoNode(flags,devicetype,id,locid,serialnumber,
                                 description,handle))
  end
  return infonodearray
end


text:
version = v"1.0"

# Collection of sources required to build this package
sources = [
    GitSource("https://github.com/oneapi-src/level-zero.git",
              "fcc7b7aceacf3cbfabaf3c0952ae0cc02d083592"),
]

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
FractalLattice = Union{Classical2dFractalLattice,Quantum2dFractalLattice}

function getTensorP(lattice::FractalLattice, which::Int = 0)
	if which == 0
		return lattice.tensorP
	else
		return lattice.tensorP[which]
	end
end

function getTensorQ(lattice::FractalLattice)
	return lattice.tensorQ
end

function setTensorP!(lattice::FractalLattice, which::Int, tensorP::Array)
	lattice.tensorP[which] = tensorP
end

function setTensorQ!(lattice::FractalLattice, tensorQ::Array)
	lattice.tensorQ = tensorQ
end

function getLegextension(lattice::FractalLattice)
	return lattice.legextension
end

function getHausdorffDim(lattice::FractalLattice)
	L = getLegextension(lattice)
	d = log(4.0+8.0*L)/log(2.0+2.0*L)
	# d = log(convert(T,4+8*L))/log(convert(T,2+2*L))
	return d
end

function getFractalDim(lattice::FractalLattice)
	L = getLegextension(lattice)
	d = 1.0 + log(2.0)/log(2.0+2.0*L)
	return d
end


text:
using YaoBase
export Concentrator, concentrate

"""
    Concentrator{N, T, BT <: AbstractBlock} <: AbstractContainer{BT, N, T}

concentrates serveral lines together in the circuit, and expose
it to other blocks.
"""
struct Concentrator{N, BT <: AbstractBlock, C} <: AbstractContainer{BT, N}
    content::BT
    locations::NTuple{C, Int}
end

function Concentrator{N}(block::BT, locations::NTuple{C, Int}) where {N, M, C, BT<:AbstractBlock{M}}
    if !(length(locations) == M && N>=M)
        throw(LocationConflictError("length of locations must be equal to the size of block, and smaller than size of itself."))
    end
    return Concentrator{N, BT, C}(block, locations)
end

"""
    concentrate(n, block, locs)

Create a [`Concentrator`](@ref) block with total number of current active qubits `n`,
which concentrates given wire location together to `length(locs)` active qubits,
and relax the concentration afterwards.

# Example

Concentrator is equivalent to [`put`](@ref) a block on given position mathematically, but more efficient
and convenient for large blocks.

```jldoctest
julia> r = rand_state(3)
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 3/3

julia> apply!(copy(r), concentrate(X, 1)) ≈ apply!(copy(r), put(1=>X))
true
```

It works for in-contigious locations as well

```jldoctest
julia> r = rand_state(4)
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 4/4

julia> cc = concentrate(4, kron(X, Y), (1, 3))
nqubits: 4
Concentrator: (1, 3)
└─ kron
   ├─ 1=>X gate
   └─ 2=>Y gate

julia> pp = chain(4, put(1=>X), put(3=>Y))
nqubits: 4
chain
├─ put on (1)
│  └─ X gate
└─ put on (3)
   └─ Y gate

julia> apply!(copy(r), cc) ≈ apply!(copy(r), pp)
true
```
"""
function concentrate(n::Int, block::AbstractBlock, locs)
    return Concentrator{n}(block, Tuple(locs))
end

# support lazy qubits
concentrate(n::Int, block::Function, locs) = concentrate(n, parse_block(length(locs), block), locs)

"""
    concentrate(block, locs) -> f(n)

Lazy curried version of [`concentrate`](@ref).
"""
concentrate(block::AbstractBlock, locs) = @λ(n->concentrate(n, block, locs))
concentrate(block::Function, locs) = @λ(n->concentrate(n, block, locs))

occupied_locs(c::Concentrator) = c.locations
chsubblocks(pb::Concentrator{N}, blk::AbstractBlock) where N =
    Concentrator{N}(blk, occupied_locs(pb))
PreserveStyle(::Concentrator) = PreserveAll()

function apply!(r::AbstractRegister, c::Concentrator)
    _check_size(r, c)
    focus!(r, occupied_locs(c))
    apply!(r, c.content)
    relax!(r, occupied_locs(c), to_nactive=nqubits(c))
    return r
end

function mat(::Type{T}, c::Concentrator{N, <:AbstractBlock}) where {N, T}
    mat(T, PutBlock{N}(c.content, c.locations))
end

Base.adjoint(blk::Concentrator{N}) where N =
    Concentrator{N}(adjoint(blk.content), occupied_locs(blk))

function Base.:(==)(a::Concentrator{N, BT}, b::Concentrator{N, BT}) where {N, BT}
    return a.content == b.content && a.locations == b.locations
end

YaoBase.nqubits(::Concentrator{N}) where N = N
YaoBase.nactive(c::Concentrator) = length(c.locations)

function YaoBase.iscommute(x::Concentrator{N}, y::Concentrator{N}) where N
    if occupied_locs(x) == occupied_locs(y)
        return iscommute(x.content, y.content)
    else
        return iscommute_fallback(x, y)
    end
end


text:
using LinearAlgebra
using OpenCL

const sum_kernel = "
   __kernel void sum(__global float *a,
                     __global const float *b)
    {
      int gid = get_global_id(0);
      a[gid] = a[gid] + b[gid];
    }
"
a = zeros(Float32, 50_000)
b = ones(Float32, 50_000)

device, ctx, queue = cl.create_compute_context()

a_buff = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=a)
b_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=b)

p = cl.Program(ctx, source=sum_kernel) |> cl.build!
k = cl.Kernel(p, "sum")

queue(k, size(a), nothing, a_buff, b_buff)

r = cl.read(queue, a_buff)

if isapprox(r, fill(Float32(1), 50000))
    @info "Success!"
else
    @error "Norm should be 0.0f"
end

max_stars_repo_name: zazbone/ComputationalPhysics

max_issues_repo_name: zazbone/ComputationalPhysics

max_forks_repo_name: zazbone/ComputationalPhysics


text:
using Kinetic, Plots, Flux
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress

cd(@__DIR__)
#@load "../nn_scalar.jld2" nn
@load "/home2/vavrines/Coding/Flowmachine/src/1d/sampler/nn_rif.jld2" nn

set = Setup(case = "sod", space = "1d2f1v", maxTime = 0.15)
ps = PSpace1D(0.0, 1.0, 200, 1)
vs = VSpace1D(-5.0, 5.0, 100)
gas = Gas(Kn = 1e-4)
ib = IB2F(ib_sod(set, ps, vs, gas)...)
ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks, ks.ps, :static_array; structarray = true)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ctr[1].w)

@showprogress for iter = 1:nt
    reconstruct!(ks, ctr)

    for i = 1:ks.ps.nx+1
        w = (ctr[i-1].w .+ ctr[i].w) ./ 2
        prim = (ctr[i-1].prim .+ ctr[i].prim) ./ 2
        sw = (ctr[i].w .- ctr[i-1].w) / ks.ps.dx[i]

        #L = abs(ctr[i].w[1] / sw[1])
        #ℓ = (1/prim[end])^ks.gas.ω / prim[1] * sqrt(prim[end]) * ks.gas.Kn
        #KnGLL = ℓ / L
        #isNS = ifelse(KnGLL > 0.05, false, true)

        #h = (ctr[i-1].h .+ ctr[i].h) ./ 2
        #x, y = regime_data(ks, w, prim, sw, h)
        #isNS = ifelse(onecold(y) == 1, true, false)

        tau = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)
        regime = Int(round(nn([w; sw; tau])[1]))

        if regime == 0
            flux_gks!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[i-1].w .+ ctr[i-1].sw .* ks.ps.dx[i-1] / 2,
                ctr[i].w .- ctr[i].sw .* ks.ps.dx[i] / 2,
                ks.vs.u,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1] / 2,
                ks.ps.dx[i] / 2,
                ctr[i-1].sw,
                ctr[i].sw,
            )
        else
            flux_kfvs!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[i-1].h .+ 0.5 .* ctr[i-1].sh .* ks.ps.dx[i-1],
                ctr[i-1].b .+ 0.5 .* ctr[i-1].sb .* ks.ps.dx[i-1],
                ctr[i].h .- 0.5 .* ctr[i].sh .* ks.ps.dx[i],
                ctr[i].b .- 0.5 .* ctr[i].sb .* ks.ps.dx[i],
                ks.vs.u,
                ks.vs.weights,
                dt,
                ctr[i-1].sh,
                ctr[i-1].sb,
                ctr[i].sh,
                ctr[i].sb,
            )
        end
    end

    update!(ks, ctr, face, dt, res; coll = :bgk, bc = [:fix, :fix])

    t += dt
end

plot(ks, ctr)

cd(@__DIR__)
#@save "pure_kinetic.jld2" ks ctr
#@save "pure_ns.jld2" ks ctr
#@save "nn.jld2" ks ctr
#@save "kngll.jld2" ks ctr

begin
    #regime = zeros(Int, ks.ps.nx)
    #for i = 1:ks.ps.nx
    #    sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
    #    x, y = regime_data(ks, ctr[i].w, sw, ctr[i].h)
    #    regime[i] = Int(round(y))
    #end
    regime = zeros(Int, ks.ps.nx)
    for i = 1:ks.ps.nx
        sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
        tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
        x, y = regime_data(ks, ctr[i].w, sw, ctr[i].h)

        #regime[i] = nn(x) |> onecold
        regime[i] = round(nn(x)[1])
    end

    regime1 = zeros(Int, ks.ps.nx)
    KnGLL = zeros(ks.ps.nx)
    for i = 1:ks.ps.nx
        sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
        L = abs(ctr[i].w[1] / sw[1])
        ℓ = (1/ctr[i].prim[end])^ks.gas.ω / ctr[i].prim[1] * sqrt(ctr[i].prim[end]) * ks.gas.Kn

        KnGLL[i] = ℓ / L
        regime1[i] = ifelse(KnGLL[i] >= 0.05, 1, 0)
    end
end

plot(ks.ps.x[1:ks.ps.nx], regime)
plot!(ks.ps.x[1:ks.ps.nx], regime1)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using JLD2, Distributions, Serialization
import GEAS:sim_pt_QTL, breeding_program, top_QTL, QTL

"""
    function run_2021_06
---
This test is the first formal test on editint.
I will use the real genotype first.
This is because that the genotypes from `QMSimu` and `ms` still look wierd to me.

## Scenario to be considered
- Number of QTL, 100:450:1000
- Number of full sibs per family, 60:20:100.  Half of them are to be challenged
- number of QTL to be edited, 3, or 5 biggest in the base
  - I send a subset of the simulated QTL for editing
  - The subset is determined in the breeding program, not here.
- Correlation `r` between the QTL of the 2 traits.
- Number of sibs per full sib family
  - Always challenge half of them
- Editing vs. non editing
- Number of repeats will be as number of instances to be submitted on `Saga`.
"""
function run_2021_06()
    @load "dat/run/base.jld2" base
    rst = "dat/editing"
    isdir(rst) && rm(rst, force=true, recursive=true)
    isdir(rst) || mkpath(rst)

    for nqtl in 100:450:1000
        for r in -.3:.15:.3
            d = MvNormal(zeros(2), [1 r; r 1])
            qtl = sim_pt_QTL(base, nqtl, d)
            RR = Int(r*100)
            serialize("$rst/qtl-$nqtl-$RR.ser", qtl)
            for nsib in 30:10:50
                par = begin
                    Parameters = Dict(
                        :nSire=> 100,
                        :nDam => 200, # male:female = 1:2
                        :nSib => 2nsib,
                        :nC7e => nsib, # number of sib to be challenged
                        :nG8n => 10,   # number of generations
                        :nQTL => [nqtl, nqtl],
                        :h²   => [.5, .5],
                        :p8e  => .5, # percentage of affected in the binary trait
                        :e19e => 1., # edit_successsful_rate
                        :w4t  => 1, # weight on the binary trait EBV
                        :nk3n => 3  # number of known QTL
                    )
                    (; Parameters...) # named tuple.  contents as above
                end
                prd1, snp1 = breeding_program(base, par, qtl)
                serialize("$rst/prd-$nqtl-$RR-$nsib-n.ser", prd1)
                serialize("$rst/snp-$nqtl-$RR-$nsib-n.ser", snp1)
                prd2, snp2 = breeding_program(base, par, qtl, edit=true)
                serialize("$rst/prd-$nqtl-$RR-$nsib-y.ser", prd2)
                serialize("$rst/snp-$nqtl-$RR-$nsib-y.ser", snp2)
            end
        end
    end
end

function my_test()
    @load "dat/run/base.jld2" base
    r = .5
    d = MvNormal(zeros(2), [1 r; r 1])
    qtl = sim_pt_QTL(base, 100, d)
    q = QTL(qtl[2].pos[1:3], qtl[2].effect[1:3], 0, 0)
end


text:
#=
This file contains explicit SIMD instructions for tensors.
Many of the methods defined outside this file will use SIMD-instructions
if julia is ran with -O3. Even if -O3 is enabled, the compiler is sometimes
thrown off guard, and therefore, explicit SIMD routines are
defined. This will enable SIMD-instructions even if julia is ran with
the default -O2.

The functions here are only defined for tensors of the same
element type. Otherwise it does work. Promotion should take
care of this before the tensors enter the functions here.

The file is organized as follows:
(1): + and - between tensors
(2): * and / between tensor and number
(3): dot
(4): dcontract
(5): otimes
(6): norm
=#

import SIMD
const SVec{N, T} = SIMD.Vec{N, T}

const SIMDTypes = Union{Float16, Float32, Float64}

const AllSIMDTensors{T <: SIMDTypes} = Union{Tensor{1, 1, T, 1}, Tensor{1, 2, T, 2}, Tensor{1, 3, T, 3},
                                             Tensor{2, 1, T, 1}, Tensor{2, 2, T, 4}, Tensor{2, 3, T, 9},
                                             Tensor{4, 1, T, 1}, Tensor{4, 2, T, 16}, #=Tensor{4, 3, T, 81},=#
                                             SymmetricTensor{2, 1, T, 1}, SymmetricTensor{2, 2, T, 3}, SymmetricTensor{2, 3, T, 6},
                                             SymmetricTensor{4, 1, T, 1}, SymmetricTensor{4, 2, T, 9}, SymmetricTensor{4, 3, T, 36}}

# SIMD sizes accepted by LLVM between 1 and 100
const SIMD_CHUNKS = (1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 16, 17, 18, 20, 24, 32, 33, 34, 36, 40, 48, 64, 65, 66, 68, 72, 80, 96)

# factors for the symmetric tensors, return a quote
function symmetric_factors(order, dim, T)
    if order == 2
        dim == 1 && return :(SVec{1, $T}(($T(1),)))
        dim == 2 && return :(SVec{3, $T}(($T(1),$T(2),$T(1))))
        dim == 3 && return :(SVec{6, $T}(($T(1),$T(2),$T(2),$T(1),$T(2),$T(1))))
    elseif order == 4
        dim == 1 && return :(SVec{1, $T}(($T(1),)))
        dim == 2 && return :(SVec{9, $T}(($T(1),$T(2),$T(1),$T(2),$T(4),$T(2),$T(1),$T(2),$T(1))))
        dim == 3 && return :(SVec{36,$T}(($T(1),$T(2),$T(2),$T(1),$T(2),$T(1),$T(2),$T(4),$T(4),$T(2),$T(4),$T(2),
                                          $T(2),$T(4),$T(4),$T(2),$T(4),$T(2),$T(1),$T(2),$T(2),$T(1),$T(2),$T(1),
                                          $T(2),$T(4),$T(4),$T(2),$T(4),$T(2),$T(1),$T(2),$T(2),$T(1),$T(2),$T(1))))
    end
end

# Tensor from SVec (Note: This needs to be two separate methods in order to dispatch correctly)
@generated function (::Type{Tensor{order, dim}})(r::SVec{N, T}) where {order, dim, N, T}
    return quote
        $(Expr(:meta, :inline))
        Tensor{$order, $dim}($(Expr(:tuple, [:(r[$i]) for i in 1:N]...)))
    end
end
@generated function (::Type{SymmetricTensor{order, dim}})(r::SVec{N, T}) where {order, dim, N, T}
    return quote
        $(Expr(:meta, :inline))
        SymmetricTensor{$order, $dim}($(Expr(:tuple, [:(r[$i]) for i in 1:N]...)))
    end
end
# constructor from several SVecs which happens in dot, dcontract and otimes
@generated function (::Type{Tensor{order, dim}})(r::NTuple{M, SVec{N, T}}) where {order, dim, M, N, T}
    return quote
        $(Expr(:meta, :inline))
        @inbounds return Tensor{$order, $dim}($(Expr(:tuple, [:(r[$j][$i]) for i in 1:N, j in 1:M]...)))
    end
end
@generated function (::Type{SymmetricTensor{order, dim}})(r::NTuple{M, SVec{N, T}}) where {order, dim, M, N, T}
    return quote
        $(Expr(:meta, :inline))
        @inbounds return SymmetricTensor{$order, $dim}($(Expr(:tuple, [:(r[$j][$i]) for i in 1:N, j in 1:M]...)))
    end
end

# load Tensor into SVec
@inline tosimd(A::NTuple{N, T}) where {N, T} = SVec{N, T}(A)

# load given range of linear indices into SVec
@generated function tosimd(D::NTuple{N, T}, ::Type{Val{strt}}, ::Type{Val{stp}}) where {N, T, strt, stp}
    expr = Expr(:tuple, [:(D[$i]) for i in strt:stp]...)
    M = length(expr.args)
    return quote
        $(Expr(:meta, :inline))
        @inbounds return SVec{$M, T}($expr)
    end
end

################################
# (1): + and - between tensors #
################################
@inline function Base.:+(S1::TT, S2::TT) where {TT <: AllSIMDTensors}
    @inbounds begin
        D1 = get_data(S1); SV1 = tosimd(D1)
        D2 = get_data(S2); SV2 = tosimd(D2)
        r = SV1 + SV2
        return get_base(TT)(r)
    end
end
@generated function Base.:+(S1::Tensor{4, 3, T}, S2::Tensor{4, 3, T}) where {T <: SIMDTypes}
    return quote
        $(Expr(:meta, :inline))
        @inbounds begin
            D1 = get_data(S1); SV180 = tosimd(D1, Val{1}, Val{80})
            D2 = get_data(S2); SV280 = tosimd(D2, Val{1}, Val{80})
            r = SV180 + SV280
            r81 = D1[81] + D2[81]
            return Tensor{4, 3}($(Expr(:tuple, [:(r[$i]) for i in 1:80]..., :(r81))))
        end
    end
end
@inline function Base.:-(S1::TT, S2::TT) where {TT <: AllSIMDTensors}
    @inbounds begin
        D1 = get_data(S1); SV1 = tosimd(D1)
        D2 = get_data(S2); SV2 = tosimd(D2)
        r = SV1 - SV2
        return get_base(TT)(r)
    end
end
@generated function Base.:-(S1::Tensor{4, 3, T}, S2::Tensor{4, 3, T}) where {T <: SIMDTypes}
    return quote
        $(Expr(:meta, :inline))
        @inbounds begin
            D1 = get_data(S1); SV180 = tosimd(D1, Val{1}, Val{80})
            D2 = get_data(S2); SV280 = tosimd(D2, Val{1}, Val{80})
            r = SV180 - SV280
            r81 = D1[81] - D2[81]
            return Tensor{4, 3}($(Expr(:tuple, [:(r[$i]) for i in 1:80]..., :(r81))))
        end
    end
end

##########################################
# (2): * and / between tensor and number #
##########################################
@inline function Base.:*(n::T, S::AllSIMDTensors{T}) where {T <: SIMDTypes}
    @inbounds begin
        D = get_data(S); SV = tosimd(D)
        r = n * SV
        return get_base(typeof(S))(r)
    end
end
@inline function Base.:*(S::AllSIMDTensors{T}, n::T) where {T <: SIMDTypes}
    @inbounds begin
        D = get_data(S); SV = tosimd(D)
        r = SV * n
        return get_base(typeof(S))(r)
    end
end
@inline function Base.:/(S::AllSIMDTensors{T}, n::T) where {T <: SIMDTypes}
    @inbounds begin
        D = get_data(S); SV = tosimd(D)
        r = SV / n
        return get_base(typeof(S))(r)
    end
end
@generated function Base.:*(n::T, S::Tensor{4, 3, T}) where {T <: SIMDTypes}
    return quote
        $(Expr(:meta, :inline))
        @inbounds begin
            D = get_data(S); SV80 = tosimd(D, Val{1}, Val{80})
            r = n * SV80; r81 = n * D[81]
            return Tensor{4, 3}($(Expr(:tuple, [:(r[$i]) for i in 1:80]..., :(r81))))
        end
    end
end
@generated function Base.:*(S::Tensor{4, 3, T}, n::T) where {T <: SIMDTypes}
    return quote
        $(Expr(:meta, :inline))
        @inbounds begin
            D = get_data(S); SV80 = tosimd(D, Val{1}, Val{80})
            r = SV80 * n; r81 = D[81] * n
            return Tensor{4, 3}($(Expr(:tuple, [:(r[$i]) for i in 1:80]..., :(r81))))
        end
    end
end
@generated function Base.:/(S::Tensor{4, 3, T}, n::T) where {T <: SIMDTypes}
    return quote
        $(Expr(:meta, :inline))
        @inbounds begin
            D = get_data(S); SV80 = tosimd(D, Val{1}, Val{80})
            r = SV80 / n; r81 = D[81] / n
            return Tensor{4, 3}($(Expr(:tuple, [:(r[$i]) for i in 1:80]..., :(r81))))
        end
    end
end

############
# (3): dot #
############
# 2s-1
@inline LinearAlgebra.dot(S1::SymmetricTensor{2, dim, T}, S2::Vec{dim, T}) where {dim, T <: SIMDTypes} = dot(promote(S1), S2)
# 2-1
@inline function LinearAlgebra.dot(S1::Tensor{2, 2, T}, S2::Vec{2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1}, Val{2})
        SV12 = tosimd(D1, Val{3}, Val{4})
        r = muladd(SV12, D2[2], SV11 * D2[1])
        return Tensor{1, 2}(r)
    end
end
@inline function LinearAlgebra.dot(S1::Tensor{2, 3, T}, S2::Vec{3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1}, Val{3})
        SV12 = tosimd(D1, Val{4}, Val{6})
        SV13 = tosimd(D1, Val{7}, Val{9})
        r = muladd(SV13, D2[3], muladd(SV12, D2[2], SV11 * D2[1]))
        return Tensor{1, 3}(r)
    end
end

# 2s-2 / 2-2s
@inline function LinearAlgebra.dot(S1::AbstractTensor{2, dim, T}, S2::AbstractTensor{2, dim, T}) where {dim, T <: SIMDTypes}
    SS1, SS2 = promote_base(S1, S2); dot(SS1, SS2)
end
# 2-2
@inline function LinearAlgebra.dot(S1::Tensor{2, 2, T}, S2::Tensor{2, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1}, Val{2})
        SV12 = tosimd(D1, Val{3}, Val{4})
        r1 = muladd(SV12, D2[2], SV11 * D2[1])
        r2 = muladd(SV12, D2[4], SV11 * D2[3])
        return Tensor{2, 2}((r1, r2))
    end
end
@inline function LinearAlgebra.dot(S1::Tensor{2, 3, T}, S2::Tensor{2, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1}, Val{3})
        SV12 = tosimd(D1, Val{4}, Val{6})
        SV13 = tosimd(D1, Val{7}, Val{9})
        r1 = muladd(SV13, D2[3], muladd(SV12, D2[2], SV11 * D2[1]))
        r2 = muladd(SV13, D2[6], muladd(SV12, D2[5], SV11 * D2[4]))
        r3 = muladd(SV13, D2[9], muladd(SV12, D2[8], SV11 * D2[7]))
        return Tensor{2, 3}((r1, r2, r3))
    end
end

##################
# (4): dcontract #
##################
# 2s-2 / 2-2s
@inline function dcontract(S1::AbstractTensor{2, dim, T}, S2::AbstractTensor{2, dim, T}) where {dim, T <: SIMDTypes}
    SS1, SS2 = promote_base(S1, S2); dcontract(SS1, SS2)
end
# 2-2
@inline function dcontract(S1::Tensor{2, dim, T}, S2::Tensor{2, dim, T}) where {dim, T <: SIMDTypes}
    SV1 = tosimd(get_data(S1))
    SV2 = tosimd(get_data(S2))
    r = SV1 * SV2
    return sum(r)
end
# 2s-2s
@generated function dcontract(S1::SymmetricTensor{2, dim, T}, S2::SymmetricTensor{2, dim, T}) where {dim, T <: SIMDTypes}
    F = symmetric_factors(2, dim, T)
    return quote
        $(Expr(:meta, :inline))
        F = $F
        SV1 = tosimd(get_data(S1))
        SV2 = tosimd(get_data(S2))
        SV1SV2 = SV1 * SV2; FSV1SV2 = F * SV1SV2
        return sum(FSV1SV2)
    end
end

# 4-2s
@inline dcontract(S1::Tensor{4, 3, T}, S2::SymmetricTensor{2, 3, T}) where {T <: SIMDTypes} = dcontract(S1, promote(S2))
# 4-2
@inline function dcontract(S1::Tensor{4, 2, T}, S2::Tensor{2, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1}, Val{4})
        SV12 = tosimd(D1, Val{5}, Val{8})
        SV13 = tosimd(D1, Val{9}, Val{12})
        SV14 = tosimd(D1, Val{13}, Val{16})
        r = muladd(SV14, D2[4], muladd(SV13, D2[3], muladd(SV12, D2[2], SV11 * D2[1])))
        return Tensor{2, 2}(r)
    end
end
@inline function dcontract(S1::Tensor{4, 3, T}, S2::Tensor{2, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{9})
        SV12 = tosimd(D1, Val{10}, Val{18})
        SV13 = tosimd(D1, Val{19}, Val{27})
        SV14 = tosimd(D1, Val{28}, Val{36})
        SV15 = tosimd(D1, Val{37}, Val{45})
        SV16 = tosimd(D1, Val{46}, Val{54})
        SV17 = tosimd(D1, Val{55}, Val{63})
        SV18 = tosimd(D1, Val{64}, Val{72})
        SV19 = tosimd(D1, Val{73}, Val{81})
        r = muladd(SV19, D2[9], muladd(SV18, D2[8], muladd(SV17, D2[7], muladd(SV16, D2[6], muladd(SV15, D2[5], muladd(SV14, D2[4], muladd(SV13, D2[3], muladd(SV12, D2[2], SV11 * D2[1]))))))))
        return Tensor{2, 3}(r)
    end
end
@inline function dcontract(S1::SymmetricTensor{4, 3, T}, S2::SymmetricTensor{2, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{6})
        SV12 = tosimd(D1, Val{7},  Val{12})
        SV13 = tosimd(D1, Val{13}, Val{18})
        SV14 = tosimd(D1, Val{19}, Val{24})
        SV15 = tosimd(D1, Val{25}, Val{30})
        SV16 = tosimd(D1, Val{31}, Val{36})
        D21 = D2[1]; D22 = D2[2] * T(2); D23 = D2[3] * T(2)
        D24 = D2[4]; D25 = D2[5] * T(2); D26 = D2[6]
        r = muladd(SV16, D26, muladd(SV15, D25, muladd(SV14, D24, muladd(SV13, D23, muladd(SV12, D22, SV11 * D21)))))
        return SymmetricTensor{2, 3}(r)
    end
end

# 4s-4 / 4-4s
@inline function dcontract(S1::AbstractTensor{4, dim, T}, S2::AbstractTensor{4, dim, T}) where {dim, T <: SIMDTypes}
    SS1, SS2 = promote_base(S1, S2); dcontract(SS1, SS2)
end
# 4-4
@inline function dcontract(S1::Tensor{4, 2, T}, S2::Tensor{4, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{4})
        SV12 = tosimd(D1, Val{5},  Val{8})
        SV13 = tosimd(D1, Val{9},  Val{12})
        SV14 = tosimd(D1, Val{13}, Val{16})
        r1 = muladd(SV14, D2[4],  muladd(SV13, D2[3],  muladd(SV12, D2[2],  SV11 * D2[1])))
        r2 = muladd(SV14, D2[8],  muladd(SV13, D2[7],  muladd(SV12, D2[6],  SV11 * D2[5])))
        r3 = muladd(SV14, D2[12], muladd(SV13, D2[11], muladd(SV12, D2[10], SV11 * D2[9])))
        r4 = muladd(SV14, D2[16], muladd(SV13, D2[15], muladd(SV12, D2[14], SV11 * D2[13])))
        return Tensor{4, 2}((r1, r2, r3, r4))
    end
end
function dcontract(S1::Tensor{4, 3, T}, S2::Tensor{4, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{9})
        SV12 = tosimd(D1, Val{10}, Val{18})
        SV13 = tosimd(D1, Val{19}, Val{27})
        SV14 = tosimd(D1, Val{28}, Val{36})
        SV15 = tosimd(D1, Val{37}, Val{45})
        SV16 = tosimd(D1, Val{46}, Val{54})
        SV17 = tosimd(D1, Val{55}, Val{63})
        SV18 = tosimd(D1, Val{64}, Val{72})
        SV19 = tosimd(D1, Val{73}, Val{81})
        r1 = muladd(SV19, D2[9],  muladd(SV18, D2[8],  muladd(SV17, D2[7],  muladd(SV16, D2[6],  muladd(SV15, D2[5],  muladd(SV14, D2[4],  muladd(SV13, D2[3],  muladd(SV12, D2[2],  SV11 * D2[1] ))))))))
        r2 = muladd(SV19, D2[18], muladd(SV18, D2[17], muladd(SV17, D2[16], muladd(SV16, D2[15], muladd(SV15, D2[14], muladd(SV14, D2[13], muladd(SV13, D2[12], muladd(SV12, D2[11], SV11 * D2[10]))))))))
        r3 = muladd(SV19, D2[27], muladd(SV18, D2[26], muladd(SV17, D2[25], muladd(SV16, D2[24], muladd(SV15, D2[23], muladd(SV14, D2[22], muladd(SV13, D2[21], muladd(SV12, D2[20], SV11 * D2[19]))))))))
        r4 = muladd(SV19, D2[36], muladd(SV18, D2[35], muladd(SV17, D2[34], muladd(SV16, D2[33], muladd(SV15, D2[32], muladd(SV14, D2[31], muladd(SV13, D2[30], muladd(SV12, D2[29], SV11 * D2[28]))))))))
        r5 = muladd(SV19, D2[45], muladd(SV18, D2[44], muladd(SV17, D2[43], muladd(SV16, D2[42], muladd(SV15, D2[41], muladd(SV14, D2[40], muladd(SV13, D2[39], muladd(SV12, D2[38], SV11 * D2[37]))))))))
        r6 = muladd(SV19, D2[54], muladd(SV18, D2[53], muladd(SV17, D2[52], muladd(SV16, D2[51], muladd(SV15, D2[50], muladd(SV14, D2[49], muladd(SV13, D2[48], muladd(SV12, D2[47], SV11 * D2[46]))))))))
        r7 = muladd(SV19, D2[63], muladd(SV18, D2[62], muladd(SV17, D2[61], muladd(SV16, D2[60], muladd(SV15, D2[59], muladd(SV14, D2[58], muladd(SV13, D2[57], muladd(SV12, D2[56], SV11 * D2[55]))))))))
        r8 = muladd(SV19, D2[72], muladd(SV18, D2[71], muladd(SV17, D2[70], muladd(SV16, D2[69], muladd(SV15, D2[68], muladd(SV14, D2[67], muladd(SV13, D2[66], muladd(SV12, D2[65], SV11 * D2[64]))))))))
        r9 = muladd(SV19, D2[81], muladd(SV18, D2[80], muladd(SV17, D2[79], muladd(SV16, D2[78], muladd(SV15, D2[77], muladd(SV14, D2[76], muladd(SV13, D2[75], muladd(SV12, D2[74], SV11 * D2[73]))))))))
        return Tensor{4, 3}((r1, r2, r3, r4, r5, r6, r7, r8, r9))
    end
end
@inline function dcontract(S1::SymmetricTensor{4, 3, T}, S2::Tensor{2, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{6})
        SV12 = tosimd(D1, Val{7}, Val{12})
        SV13 = tosimd(D1, Val{13}, Val{18})
        SV14 = tosimd(D1, Val{19}, Val{24})
        SV15 = tosimd(D1, Val{25}, Val{30})
        SV16 = tosimd(D1, Val{31}, Val{36})
        D21 = D2[1]; D22 = D2[2] + D2[4]; D23 = D2[3] + D2[7]
        D24 = D2[5]; D25 = D2[6] + D2[8]; D26 = D2[9]
        r = muladd(SV16, D26, muladd(SV15, D25, muladd(SV14, D24, muladd(SV13, D23, muladd(SV12, D22, SV11 * D21)))))
        return SymmetricTensor{2, 3}(r)
    end
end
@inline function dcontract(S1::SymmetricTensor{4, 2, T}, S2::SymmetricTensor{4, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{3})
        SV12 = tosimd(D1, Val{4},  Val{6})
        SV13 = tosimd(D1, Val{7},  Val{9})
        D21  = D2[1]; D22 = D2[2] * T(2); D23 = D2[3]
        D24  = D2[4]; D25 = D2[5] * T(2); D26 = D2[6]
        D27  = D2[7]; D28 = D2[8] * T(2); D29 = D2[9]
        r1 = muladd(SV13, D23, muladd(SV12, D22, SV11 * D21))
        r2 = muladd(SV13, D26, muladd(SV12, D25, SV11 * D24))
        r3 = muladd(SV13, D29, muladd(SV12, D28, SV11 * D27))
        return SymmetricTensor{4, 2}((r1, r2, r3))
    end
end
function dcontract(S1::SymmetricTensor{4, 3, T}, S2::SymmetricTensor{4, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV11 = tosimd(D1, Val{1},  Val{6})
        SV12 = tosimd(D1, Val{7},  Val{12})
        SV13 = tosimd(D1, Val{13}, Val{18})
        SV14 = tosimd(D1, Val{19}, Val{24})
        SV15 = tosimd(D1, Val{25}, Val{30})
        SV16 = tosimd(D1, Val{31}, Val{36})
        D21  = D2[1];  D22  = D2[2]  * T(2); D23  = D2[3]  * T(2); D24  = D2[4];  D25  = D2[5]  * T(2); D26  = D2[6]
        D27  = D2[7];  D28  = D2[8]  * T(2); D29  = D2[9]  * T(2); D210 = D2[10]; D211 = D2[11] * T(2); D212 = D2[12]
        D213 = D2[13]; D214 = D2[14] * T(2); D215 = D2[15] * T(2); D216 = D2[16]; D217 = D2[17] * T(2); D218 = D2[18]
        D219 = D2[19]; D220 = D2[20] * T(2); D221 = D2[21] * T(2); D222 = D2[22]; D223 = D2[23] * T(2); D224 = D2[24]
        D225 = D2[25]; D226 = D2[26] * T(2); D227 = D2[27] * T(2); D228 = D2[28]; D229 = D2[29] * T(2); D230 = D2[30]
        D231 = D2[31]; D232 = D2[32] * T(2); D233 = D2[33] * T(2); D234 = D2[34]; D235 = D2[35] * T(2); D236 = D2[36]
        r1 = muladd(SV16, D26,  muladd(SV15, D25,  muladd(SV14, D24,  muladd(SV13, D23,  muladd(SV12, D22,  SV11 * D21 )))))
        r2 = muladd(SV16, D212, muladd(SV15, D211, muladd(SV14, D210, muladd(SV13, D29,  muladd(SV12, D28,  SV11 * D27 )))))
        r3 = muladd(SV16, D218, muladd(SV15, D217, muladd(SV14, D216, muladd(SV13, D215, muladd(SV12, D214, SV11 * D213)))))
        r4 = muladd(SV16, D224, muladd(SV15, D223, muladd(SV14, D222, muladd(SV13, D221, muladd(SV12, D220, SV11 * D219)))))
        r5 = muladd(SV16, D230, muladd(SV15, D229, muladd(SV14, D228, muladd(SV13, D227, muladd(SV12, D226, SV11 * D225)))))
        r6 = muladd(SV16, D236, muladd(SV15, D235, muladd(SV14, D234, muladd(SV13, D233, muladd(SV12, D232, SV11 * D231)))))
        return SymmetricTensor{4, 3}((r1, r2, r3, r4, r5, r6))
    end
end

###############
# (5): otimes #
###############
# 1-1
@inline function otimes(S1::Vec{2, T}, S2::Vec{2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV1 = tosimd(D1)
        r1 = SV1 * D2[1]
        r2 = SV1 * D2[2]
        return Tensor{2, 2}((r1, r2))
    end
end
@inline function otimes(S1::Vec{3, T}, S2::Vec{3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV1 = tosimd(D1)
        r1 = SV1 * D2[1]; r2 = SV1 * D2[2]; r3 = SV1 * D2[3]
        return Tensor{2, 3}((r1, r2, r3))
    end
end

# 2s-2 / 2-2s
@inline function otimes(S1::AbstractTensor{2, dim, T}, S2::AbstractTensor{2, dim, T}) where {dim, T <: SIMDTypes}
    SS1, SS2 = promote_base(S1, S2); otimes(SS1, SS2)
end

# 2-2
@inline function otimes(S1::Tensor{2, 2, T}, S2::Tensor{2, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV1 = tosimd(D1)
        r1 = SV1 * D2[1]; r2 = SV1 * D2[2]; r3 = SV1 * D2[3]; r4 = SV1 * D2[4]
        return Tensor{4, 2}((r1, r2, r3, r4))
    end
end
@inline function otimes(S1::Tensor{2, 3, T}, S2::Tensor{2, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV1 = tosimd(D1)
        r1 = SV1 * D2[1]; r2 = SV1 * D2[2]; r3 = SV1 * D2[3]
        r4 = SV1 * D2[4]; r5 = SV1 * D2[5]; r6 = SV1 * D2[6]
        r7 = SV1 * D2[7]; r8 = SV1 * D2[8]; r9 = SV1 * D2[9]
        return Tensor{4, 3}((r1, r2, r3, r4, r5, r6, r7, r8, r9))
    end
end

# 2s-2s
@inline function otimes(S1::SymmetricTensor{2, 2, T}, S2::SymmetricTensor{2, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV1 = tosimd(D1)
        r1 = SV1 * D2[1]; r2 = SV1 * D2[2]; r3 = SV1 * D2[3]
        return SymmetricTensor{4, 2}((r1, r2, r3))
    end
end
@inline function otimes(S1::SymmetricTensor{2, 3, T}, S2::SymmetricTensor{2, 3, T}) where {T <: SIMDTypes}
    @inbounds begin
        D1 = get_data(S1); D2 = get_data(S2)
        SV1 = tosimd(D1)
        r1 = SV1 * D2[1]; r2 = SV1 * D2[2]; r3 = SV1 * D2[3]
        r4 = SV1 * D2[4]; r5 = SV1 * D2[5]; r6 = SV1 * D2[6]
        return SymmetricTensor{4, 3}((r1, r2, r3, r4, r5, r6))
    end
end

#############
# (6): norm #
#############
# order 1 and order 2 norms rely on dot and dcontract respectively
@inline function LinearAlgebra.norm(S::Tensor{4, 2, T}) where {T <: SIMDTypes}
    @inbounds begin
        SV = tosimd(get_data(S))
        SVSV = SV * SV
        r = sum(SVSV)
        return sqrt(r)
    end
end
@generated function LinearAlgebra.norm(S::Tensor{4, 3, T}) where {T <: SIMDTypes}
    return quote
        $(Expr(:meta, :inline))
        @inbounds begin
            D = get_data(S)
            SV80 = tosimd(D, Val{1}, Val{80})
            SV80SV80 = SV80 * SV80
            r80 = sum(SV80SV80)
            r = r80 + D[81] * D[81]
            return sqrt(r)
        end
    end
end
@generated function LinearAlgebra.norm(S::SymmetricTensor{4, dim, T}) where {dim, T <: SIMDTypes}
    F = symmetric_factors(4, dim, T)
    return quote
        $(Expr(:meta, :inline))
        F = $F
        SV = tosimd(get_data(S))
        SVSV = SV * SV; FSVSV = F * SVSV; r = sum(FSVSV)
        return sqrt(r)
    end
end


text:
export graph_bbc_exp

"""
    (graph,cref)=graph_bbc_exp(k;T=Float64)

Computes a polynomial evaluation approximating the exponential using `k` matrix
multiplications following the procedure in the reference. The coefficients are
directly copied from the paper.

The evaluation is embedded in the `degopt`-format, see [`graph_degopt`](@ref),
and for `k< 3`, the evaluation is using the Paterson–Stockmeyer method.

**Reference**

1. P. Bader, S. Blanes, and F. Casas. "Computing the matrix exponential with an
   optimized Taylor polynomial approximation". Mathematics, 7(12), 2019.
   DOI: [10.3390/math7121174](https://doi.org/10.3390/math7121174)
"""
function graph_bbc_exp(k; T = Float64)
    if (k < 3)
        if k == 2
            deg = 4
        elseif k == 1
            deg = 2
        else
            deg = 1
        end
        return graph_ps_degopt(convert.(T, 1 ./ factorial.(0:deg)))

    elseif (k == 3)
        # Equation 13
        value177 = convert(T, 177) # Lessen roundoff errors
        y0 = 1
        y1 = 1
        y2 = (857 - 58 * sqrt(value177)) / 630
        x3 = 2 / 3
        x1 = x3 * (1 + sqrt(value177)) / 88
        x2 = x3 * (1 + sqrt(value177)) / 352
        x4 = (-271 + 29 * sqrt(value177)) / (315 * x3)
        x5 = 11 * (-1 + sqrt(value177)) / (1260 * x3)
        x6 = 11 * (-9 + sqrt(value177)) / (5040 * x3)
        x7 = (89 - sqrt(value177)) / (5040 * x3^2)

        # A2=A^2
        v1a = [0; 1.0]
        v1b = [0; 1.0]

        # A4= A2 * (x1*A+x2*A2)
        v2a = [0; 0; 1.0]
        v2b = [0; x1; x2]

        # A8=(x3*A2+A4)*(x4*I+x5*A+x6*A2+x7*A4)
        v3a = [0; 0; x3; 1.0]
        v3b = [x4; x5; x6; x7]

        y = [y0; y1; y2; 0; 1.0]

        xv = [(v1a, v1b); (v2a, v2b); (v3a, v3b)]
        (graph, cref) = graph_degopt(xv, y)

    elseif (k == 4)
        a01 = -0.01860232051462055322
        a02 = 4.60000000000000000000
        a03 = 0.21169311829980944294
        a04 = 0
        a11 = -0.00500702322573317730
        a12 = 0.99287510353848683614
        a13 = 0.15822438471572672537
        a14 = -0.13181061013830184015
        a21 = -0.57342012296052226390
        a22 = -0.13244556105279963884
        a23 = 0.16563516943672741501
        a24 = -0.02027855540589259079
        a31 = -0.13339969394389205970
        a32 = 0.00172990000000000000
        a33 = 0.01078627793157924250
        a34 = -0.00675951846863086359

        v1a = [0; 1.0]
        v1b = [0; 1.0]
        v2a = [0; 0; 1.0]
        v2b = [0; 1; 0.0]

        # B4^2
        v3a = [a04; a14; a24; a34]
        v3b = [a04; a14; a24; a34]

        # (B2+A6)*A6
        v4b = [a03; a13; a23; a33; 1] # A6
        v4a = v4b + [a02; a12; a22; a32; 0] # B2+A6

        y = [a01; a11; a21; a31; 0; 1]

        #xv=[v1,v2,v3,v4];
        xv = [(v1a, v1b); (v2a, v2b); (v3a, v3b); (v4a, v4b)]
        (graph, cref) = graph_degopt(xv, y)

    elseif (k == 5)
        a01 = 0
        a11 = -0.10036558103014462001
        a21 = -0.00802924648241156960
        a31 = -0.00089213849804572995

        b01 = 0
        b11 = 0.39784974949964507614
        b21 = 1.36783778460411719922
        b31 = 0.49828962252538267755
        b61 = -0.00063789819459472330

        b02 = -10.9676396052962062593
        b12 = 1.68015813878906197182
        b22 = 0.05717798464788655127
        b32 = -0.00698210122488052084
        b62 = 0.00003349750170860705

        b03 = -0.09043168323908105619
        b13 = -0.06764045190713819075
        b23 = 0.06759613017704596460
        b33 = 0.02955525704293155274
        b63 = -0.00001391802575160607

        b04 = 0
        b14 = 0
        b24 = -0.09233646193671185927
        b34 = -0.01693649390020817171
        b64 = -0.00001400867981820361

        v1a = [0; 1.0]
        v1b = [0; 1.0]
        v2a = [0; 1.0; 0.0]
        v2b = [0; 0; 1.0]
        v3a = [0; 0; 0; 1.0]
        v3b = [0; 0; 0; 1.0]

        # B1B5
        v4a = [a01; a11; a21; a31; 0]
        v4b = [b04; b14; b24; b34; b64]
        # (B3+A9) * A9
        v5b = [b03; b13; b23; b33; b63; 1]  # A9
        v5a = v5b + [b02; b12; b22; b32; b62; 0]  # A9+B3
        v5 = vcat(v5a, v5b)

        # T18: B2+(B3+A9)*A9
        y = [b01; b11; b21; b31; b61; 0; 1]

        xv = [(v1a, v1b); (v2a, v2b); (v3a, v3b); (v4a, v4b); (v5a, v5b)]
        (graph, cref) = graph_degopt(xv, y)

    elseif (k == 6)
        # Coefficients unknown
        error("k=6 for BBC the coefficients are not documented.")
    else
        error("Incorrect k")
    end
    return (graph, cref)
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10]);F = [4900*x10^4-10640*x10^2*x3^2+9*x10^2*x4^2+2520*x10^2*x5*x9-13020*x10^2*x7*x9+4060*x10^2*x8^2+354*x10*x2*x4*x5-168*x10*x3*x4*x8+90*x10*x4^3+498*x10*x4^2*x6+3481*x2^2*x5^2-3304*x2*x3*x5*x8+1770*x2*x4^2*x5+9794*x2*x4*x5*x6+5776*x3^4-2736*x3^2*x5*x9+14136*x3^2*x7*x9-3624*x3^2*x8^2-840*x3*x4^2*x8-4648*x3*x4*x6*x8+225*x4^4+2490*x4^3*x6+6889*x4^2*x6^2+324*x5^2*x9^2-3348*x5*x7*x9^2+1044*x5*x8^2*x9+8649*x7^2*x9^2-5394*x7*x8^2*x9+841*x8^4+980*x10^2*x4+318*x10*x4*x5+6254*x2*x5^2-1064*x3^2*x4-2968*x3*x5*x8+1590*x4^2*x5+8798*x4*x5*x6+252*x4*x5*x9-1302*x4*x7*x9+406*x4*x8^2+49*x4^2+2809*x5^2, 354*x10*x4*x5+6962*x2*x5^2-3304*x3*x5*x8+1770*x4^2*x5+9794*x4*x5*x6+6254*x5^2, -21280*x10^2*x3-168*x10*x4*x8-3304*x2*x5*x8+23104*x3^3-5472*x3*x5*x9+28272*x3*x7*x9-7248*x3*x8^2-840*x4^2*x8-4648*x4*x6*x8-2128*x3*x4-2968*x5*x8, 18*x10^2*x4+354*x10*x2*x5-168*x10*x3*x8+270*x10*x4^2+996*x10*x4*x6+3540*x2*x4*x5+9794*x2*x5*x6-1680*x3*x4*x8-4648*x3*x6*x8+900*x4^3+7470*x4^2*x6+13778*x4*x6^2+980*x10^2+318*x10*x5-1064*x3^2+3180*x4*x5+8798*x5*x6+252*x5*x9-1302*x7*x9+406*x8^2+98*x4, 2520*x10^2*x9+354*x10*x2*x4+6962*x2^2*x5-3304*x2*x3*x8+1770*x2*x4^2+9794*x2*x4*x6-2736*x3^2*x9+648*x5*x9^2-3348*x7*x9^2+1044*x8^2*x9+318*x10*x4+12508*x2*x5-2968*x3*x8+1590*x4^2+8798*x4*x6+252*x4*x9+5618*x5, 498*x10*x4^2+9794*x2*x4*x5-4648*x3*x4*x8+2490*x4^3+13778*x4^2*x6+8798*x4*x5, -13020*x10^2*x9+14136*x3^2*x9-3348*x5*x9^2+17298*x7*x9^2-5394*x8^2*x9-1302*x4*x9, 8120*x10^2*x8-168*x10*x3*x4-3304*x2*x3*x5-7248*x3^2*x8-840*x3*x4^2-4648*x3*x4*x6+2088*x5*x8*x9-10788*x7*x8*x9+3364*x8^3-2968*x3*x5+812*x4*x8, 2520*x10^2*x5-13020*x10^2*x7-2736*x3^2*x5+14136*x3^2*x7+648*x5^2*x9-6696*x5*x7*x9+1044*x5*x8^2+17298*x7^2*x9-5394*x7*x8^2+252*x4*x5-1302*x4*x7, 19600*x10^3-21280*x10*x3^2+18*x10*x4^2+5040*x10*x5*x9-26040*x10*x7*x9+8120*x10*x8^2+354*x2*x4*x5-168*x3*x4*x8+90*x4^3+498*x4^2*x6+1960*x10*x4+318*x4*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
module RealContinuedFractions

export ContinuedFraction, contfrac, fromcontfrac, convergents, convergent

"""
    ContinuedFraction(q::Q)

Type representing a continued fraction, storing the terms with the type `Q`.
"""
struct ContinuedFraction{Q}
    q::Q # quotients
end

import Base: ==, eltype, length

==(a::ContinuedFraction, b::ContinuedFraction) = a.q == b.q

eltype(cf::ContinuedFraction) = eltype(cf.q)
length(cf::ContinuedFraction) = length(cf.q)

"""
    contfrac(x::Real)
    contfrac(x::Real, n::Integer)
    contfrac(T::Type, x::Real)
    contfrac(T::Type, x::Real, n::Integer)

Compute the first `n` terms of the continued fraction of `x`, representing it with type `T`
(defaults to `Int`).
"""
function contfrac end

contfrac(x::Real) = contfrac(Int, x)
contfrac(x::Real, n::Integer) = contfrac(Int, x, n)

function contfrac(T::Type, x::Real)
    q = T[]
    while true
        fpart, ipart = modf(x)
        push!(q, T(ipart))
        fpart == 0 && return ContinuedFraction(q)
        x = inv(fpart)
    end
end

function contfrac(T::Type, x::Real, n::Integer)
    q = T[]
    for _ in 1:n
        fpart, ipart = modf(x)
        push!(q, T(ipart))
        fpart == 0 && return ContinuedFraction(q)
        x = inv(fpart)
    end
    ContinuedFraction(q)
end

"""
    fromcontfrac(cf::ContinuedFraction)
    fromcontfrac(T::Type, cf::ContinuedFraction)

Evaluate the continued fraction `cf` using the type `T` (defaults to the rational type
associated with `eltype(cf)`).
"""
function fromcontfrac end

fromcontfrac(cf::ContinuedFraction) = fromcontfrac(_rational(eltype(cf)), cf)

function fromcontfrac(T::Type, cf::ContinuedFraction)
    T = _rational(T)
    # return foldr((x, y) -> T(x) + inv(y), cf.q[1:end-1]; init=T(cf.q[end]))
    indices = reverse(eachindex(cf.q))
    x = T(cf.q[indices[1]])
    for i in indices[2:end]
        x = T(cf.q[i]) + inv(x)
    end
    x
end

"""
    convergents(cf::ContinuedFraction)
    convergents(T::Type, cf::ContinuedFraction)

Compute the convergents of the continued fraction.
"""
function convergents end

convergents(cf::ContinuedFraction) = convergents(eltype(cf), cf)

convergents(::Type{<:Rational{T}}, cf::ContinuedFraction) where T = convergents(T, cf)

function convergents(T::Type, cf::ContinuedFraction)
    na, da, nb, db = T(1), T(0), T(cf.q[1]), T(1)
    v = Rational{T}[nb//db]
    for x in cf.q[2:end]
        na, da, nb, db = nb, db, na+x*nb, da+x*db
        push!(v, nb//db)
    end
    v
end

"""
    convergents(cf::ContinuedFraction)
    convergents(T::Type, cf::ContinuedFraction)

Compute the last convergent of the continued fraction.

These are almost equivalent to

    fromcontfrac(cf)
    fromcontfrac(Rational{T}, cf)

but perform the computation in the opposite order.

Moreover, it does not use the type `Rational` internally, so it does not check for overflow.
"""
function convergent end

convergent(cf::ContinuedFraction) = convergent(eltype(cf), cf)

convergent(::Type{<:Rational{T}}, cf::ContinuedFraction) where T = convergent(T, cf)

function convergent(T::Type, cf::ContinuedFraction)
    na, da, nb, db = T(1), T(0), T(cf.q[1]), T(1)
    for x in cf.q[2:end]
        na, da, nb, db = nb, db, na+x*nb, da+x*db
    end
    nb//db
end

"""
    _rational(::Type)

Promote integral types to rational and leave others as they are.
"""
function _rational end

_rational(T::Type{<:Integer}) = Rational{T}
# _rational(T::Type{<:Rational}) = T
_rational(T::Type) = T

end # module

max_stars_repo_name: FedericoStra/RealContinuedFractions

max_issues_repo_name: FedericoStra/RealContinuedFractions

max_forks_repo_name: FedericoStra/RealContinuedFractions


text:
using Test
using ChirpSignal

@test begin
    fs = 16000.0;
    fl = 500.0;
    fh = 8000.0;
    T  = 5.0;

    wav = chirp(T, fs, fl, fh; method="linear");
    println("linear chirp ...");

    wav = chirp(T, fs, fl, fh; method="quadratic");
    println("quadratic chirp ...");

    wav = chirp(T, fs, fl, fh; method="logarithmic");
    println("logarithmic chirp ...");

    wav = chirp(T, fs, fl, fh; method="exponential");
    println("exponential chirp ...");

    fun(x) = 2000*x
    wav = chirp(2.0, 16000.0, fun);
    wav = chirp(1.0, 16000.0, x -> 2000*x^2);
    println("customizable chirp ...");
    true
end

max_stars_repo_name: sonosole/ChirpSignal

max_issues_repo_name: sonosole/ChirpSignal

max_forks_repo_name: sonosole/ChirpSignal


text:
### Filters ###

"""
`AbstractFilter` is an abstract type to encapsulate different kinds
of discrete gaussian filters
"""
abstract type AbstractFilter end

"""
    KalmanFilter(d::LinearDynamicsModel,o::LinearObservationModel)

Construct Kalman filter with LinearDynamicsModel d and
LinearObservationModel o.
"""
struct KalmanFilter <: AbstractFilter
    d::LinearDynamicsModel
    o::LinearObservationModel
end

"""
    ExtendedKalmanFilter(d::DynamicsModel,o::ObservationModel)
    KalmanFilter(d::DynamicsModel,o::ObservationModel,λ::Number,
        α::Float,β::Float)

Construct Extended Kalman filter with DynamicsModel d and
ObservationModel o.
"""
struct ExtendedKalmanFilter <: AbstractFilter
    d::DynamicsModel
    o::ObservationModel
    ExtendedKalmanFilter(d,o) =
        (d isa LinearDynamicsModel && o isa LinearObservationModel) ?
        error("linear models: use kf over ekf") : new(d,o)
end

"""
    UnscentedKalmanFilter(d::DynamicsModel,o::ObservationModel,λ::Number,
        α::Float,β::Float)
    UnscentedKalmanFilter(d::DynamicsModel,o::ObservationModel,λ::Number)
    UnscentedKalmanFilter(d::DynamicsModel,o::ObservationModel)

Construct Unscented Kalman filter with DynamicsModel d, ObservationModel o,
and UKF parameters λ, α, and β. Default constructor uses α/β formulation from
Probabilistic Robotics, second constructor reduces complexity, third
constructor defaults λ to 2, as is commonly done.
"""
struct UnscentedKalmanFilter{a<:Number,b<:Number,c<:Number} <: AbstractFilter
    d::DynamicsModel
    o::ObservationModel
    λ::a
    α::b
    β::c
    #=
    UnscentedKalmanFilter(d,o,λ,α,β) = (d isa LinearDynamicsModel && o isa LinearObservationModel) ?
        error("linear models: use kf over ukf") : new(d,o,λ,α,β)
    =#
end

function UnscentedKalmanFilter(d::DynamicsModel,o::ObservationModel,λ::a) where a<:Number
    return UnscentedKalmanFilter{a,Int8,Int8}(d,o,λ,1,0)
end

function UnscentedKalmanFilter(d::DynamicsModel,o::ObservationModel)
    return UnscentedKalmanFilter{Int8,Int8,Int8}(d,o,2,1,0)
end

### Belief States ###

"""
    GaussianBelief(μ::AbstractVector,Σ::Symmetric)
    GaussianBelief(μ::AbstractVector,Σ::AbstractMatrix)

Construct a gaussian belief, consisting of mean vector μ
and symmetric covariance matrix Σ
"""
struct GaussianBelief{a<:Number,b<:Number}
    μ::AbstractVector{a}
    Σ::Symmetric{b}
    #=
    function GaussianBelief{a,b}(μ::Vector{a},
        Σ::Symmetric{b}) where {a<:Number,b<:Number}
        (size(μ,1) == size(Σ,1)) ?
        new(μ,Σ) : error("first dimensions not matching")
    end
    =#
end

function GaussianBelief(μ::AbstractVector,Σ::AbstractMatrix)
    return GaussianBelief(μ,Symmetric(Σ))
end

function Base.rand(rng::AbstractRNG, b::GaussianBelief)
    return b.μ + cholesky(b.Σ).L * randn(size(b.Σ,1))
end


text:
export
reset_default_graph,
get_collection,
enable_eager_execution,
control_dependencies,
has_gpu,
while_loop,
if_else,
stop_gradient,
independent,
tensor,
tensorname,
has_mpi,
get_mpi,
get_mpirun,
@cpu,
@gpu

# only for eager eager execution
enable_eager_execution() = tf.enable_eager_execution()
Base.:values(o::PyObject) = o.numpy()

"""
    reset_default_graph()

Resets the graph by removing all the operators.
"""
function reset_default_graph()
    global STORAGE
    tf.compat.v1.reset_default_graph()
    STORAGE = Dict{String, Any}()
    nothing
end
"""
    get_collection(name::Union{String, Missing})

Returns the collection with name `name`. If `name` is `missing`, returns all the trainable variables.
"""
function get_collection(name::Union{String, Missing}=missing)
    if !ismissing(name) && (name in [TRAINABLE_VARIABLES, UPDATE_OPS])
        return tf.compat.v1.get_collection(name)
    end
    if ismissing(name)
        res = tf.compat.v1.get_collection(TRAINABLE_VARIABLES)
    else
        res = []
        vs = tf.compat.v1.get_collection(TRAINABLE_VARIABLES)
        rname = @eval @r_str $name
        for v in vs
            if occursin(rname, v.name)
                push!(res, v)
            end
        end
    end
    return unique(res)
end

"""
    tensor(s::String)

Returns the tensor with name `s`. See [`tensorname`](@ref).
"""
function tensor(s::String)
    tf.get_default_graph().get_tensor_by_name(s)
end

"""
    tensorname(o::PyObject)

Returns the name of the tensor. See [`tensor`](@ref).
"""
function tensorname(o::PyObject)
    o.name
end

function jlargs(kwargs)
    kwargs = Dict{Any, Any}(kwargs)
    if :axis in keys(kwargs)
        @error("axis is not a valid keyword, using dims instead (base = 1)")
    end
    if :dtype in keys(kwargs)
        kwargs[:dtype] = DTYPE[kwargs[:dtype]]
    end
    if :dims in keys(kwargs)
        kwargs[:axis] = kwargs[:dims] .- 1
        if isa(kwargs[:axis], Array)
            kwargs[:axis] = Tuple(kwargs[:axis])
        end
        delete!(kwargs, :dims)
    end
    if :colocate in keys(kwargs)
        kwargs[:colocate_gradients_with_ops] = kwargs[:colocate]
        delete!(kwargs, :colocate)
    end
    kwargs
end

# control_dependencies can be used to fix the memory problem
# https://stackoverflow.com/questions/39350164/tensorflow-parallel-for-loop-results-in-out-of-memory

"""
    control_dependencies(f, ops::Union{Array{PyObject}, PyObject})

Executes all operations in `ops` before any operations _created_ inside the block.
```julia
op1 = tf.print("print op1")
op3 = tf.print("print op3")
control_dependencies(op1) do
    global op2 = tf.print("print op2")
end
run(sess, [op2,op3])
```
In this example, `op1` must be executed before `op2`. But there is no guarantee when `op3` will be executed.
There are several possible outputs of the program such as

```julia-repl
print op3
print op1
print op2
```
or
```
print op1
print op3
print op2
```
"""
function control_dependencies(f, ops::Union{Tuple, Array{PyObject}, PyObject})
    if isa(ops, PyObject)
        ops = [ops]
    end
    @pywith tf.control_dependencies(ops) begin
        f()
    end
end

"""
    bind(op::PyObject, ops...)

Adding operations `ops` to the dependencies of `op`. `ops` are guaranteed to be executed **before** `op`.
The function is useful when we want to execute `ops` but `ops` is not
in the dependency of the final output. For example, if we want to print `i` each time `i` is evaluated
```julia
i = constant(1.0)
op = tf.print(i)
i = bind(i, op)
```
"""
function Base.:bind(op::PyObject, ops...)
    local op1
    control_dependencies(ops) do
        op1 = tf.identity(op)
    end
    return op1
end

@doc raw"""
    while_loop(condition::Union{PyObject,Function}, body::Function, loop_vars::Union{PyObject, Array{Any}, Array{PyObject}};
        parallel_iterations::Int64=10, kwargs...)

Loops over `loop_vars` while `condition` is true. This operator only creates one extra node to mark the loops in the computational graph.

# Example

The following script computes
```math
\sum_{i=1}^{10} i
```

```julia
function condition(i, ta)
    i <= 10
end
function body(i, ta)
    u = read(ta, i-1)
    ta = write(ta, i, u+1)
    i+1, ta
end
ta = TensorArray(10)
ta = write(ta, 1, constant(1.0))
i = constant(2, dtype=Int32)
_, out = while_loop(condition, body, [i, ta])
summation = stack(out)[10]
```
"""
function while_loop(condition::Union{PyObject,Function}, body::Function, loop_vars::Union{PyObject, Array{Any}, Array{PyObject}};
        parallel_iterations::Int64=10, kwargs...)
    # @warn "TensorArray must be initialized (writedown at index 1) outside" maxlog=1
    if isa(loop_vars, PyObject)
        lv = [loop_vars]
    else
        lv = loop_vars
    end
    if get_dtype(loop_vars[1])!=Int32
        error("Loop index must be Int32, got $(get_dtype(loop_vars[1]))")
    end

    res = tf.while_loop(condition, body, loop_vars=lv; parallel_iterations=parallel_iterations, kwargs...)

    if isa(loop_vars, PyObject)
        return res[1]
    else
        return res
    end
end

function if_else_v1(condition::Union{PyObject}, fn1, fn2, args...;kwargs...)
    fn1_ = ifelse(isa(fn1, Function), fn1, ()->fn1)
    fn2_ = ifelse(isa(fn2, Function), fn2, ()->fn2)
    tf.cond(condition, fn1_, fn2_, args...;kwargs...)
end

function if_else_v2(condition::PyObject, fn1::Union{Nothing, PyObject, Array},
        fn2::Union{Nothing, PyObject, Array})
    fn1 = convert_to_tensor(fn1)
    fn2 = convert_to_tensor(fn2)
    tf.compat.v2.where(condition, fn1, fn2)
end

"""
    if_else(condition::Union{PyObject,Array,Bool}, fn1, fn2, args...;kwargs...)

- If `condition` is a scalar boolean, it outputs `fn1` or `fn2` (a function with no input argument or a tensor) based on whether `condition` is true or false.
- If `condition` is a boolean array, if returns `condition .* fn1 + (1 - condition) .* fn2`

!!! info
    If you encounter an error like this:
    ```
    tensorflow.python.framework.errors_impl.InvalidArgumentError: Retval[0] does not have value
    ```
    It's probably that your code within `if_else` is not valid.
"""
function if_else(condition::Union{PyObject,Array,Bool}, fn1, fn2, args...;kwargs...)
    if isa(condition, Array) || isa(condition, Bool)
        condition = convert_to_tensor(condition)
    end
    if isa(condition, Function) || (eltype(condition)<:Bool && length(size(condition))==0)
        if_else_v1(condition, fn1, fn2, args...;kwargs...)
    else
        if_else_v2(condition, fn1, fn2)
    end
end

"""
    has_gpu()

Checks if GPU is available.

!!! note
    ADCME will use GPU automatically if GPU is available. To disable GPU, set the environment variable `ENV["CUDA_VISIBLE_DEVICES"]=""` before importing ADCME
"""
function has_gpu()
    s = tf.test.gpu_device_name()
    if length(s)==0
        return false
    else
        return true
    end
end

"""
    has_mpi(verbose::Bool = true)

Determines whether MPI is installed.
"""
function has_mpi(verbose::Bool = true)
    if Sys.iswindows()
        if haskey(ENV, "MSMPI_INC") && haskey(ENV, "MSMPI_LIB64")
            return true
        else
            return false
        end
    end
    if haskey(ENV, "MPI_INCLUDE_PATH") && haskey(ENV, "MPI_C_LIBRARIES")
        if !(isdir(ENV["MPI_INCLUDE_PATH"]) && "mpi.h" in readdir(ENV["MPI_INCLUDE_PATH"]))
            error("mpi.h is not found in ENV[\"MPI_INCLUDE_PATH\"] = $(ENV["MPI_INCLUDE_PATH"])")
        end
        if !isfile(ENV["MPI_C_LIBRARIES"])
            error("ENV[\"MPI_C_LIBRARIES\"]=$(ENV["MPI_C_LIBRARIES"]) does not exists.")
        end
        verbose && (@info "Use MPI libraries: $(ENV["MPI_C_LIBRARIES"])")
        return true
    end
    if isfile(get_library(joinpath(ADCME.LIBDIR, "mpi"))) && isfile(joinpath(ADCME.INCDIR, "mpi.h"))
        verbose && (@info "Use default MPI library (OpenMPI)")
        return true
    end
    return false
end

"""
    get_mpi()

Returns the MPI include directory and shared library.
"""
function get_mpi()
    if Sys.iswindows()
        if haskey(ENV, "MSMPI_INC") && haskey(ENV, "MSMPI_LIB64")
            return ENV["MSMPI_INC"], joinpath(ENV["MSMPI_LIB64"], "msmpi.lib")
        else
            return false
        end
    end
    if haskey(ENV, "MPI_INCLUDE_PATH") && haskey(ENV, "MPI_C_LIBRARIES")
        return ENV["MPI_INCLUDE_PATH"], ENV["MPI_C_LIBRARIES"]
    end
    if isfile(get_library(joinpath(ADCME.LIBDIR, "mpi"))) && isfile(joinpath(ADCME.INCDIR, "mpi.h"))
        return ADCME.INCDIR, get_library(joinpath(ADCME.LIBDIR, "mpi"))
    end
    error("""MPI Library is not found.
- On Windows, you can download Microsoft MPI (https://docs.microsoft.com/en-us/message-passing-interface/microsoft-mpi)
- On Unix, you can install OpenMPI via `install_openmpi()`
""")
end

"""
    get_mpirun()

Returns the **default** mpirun executable.
"""
function get_mpirun()
    if !has_mpi(false)
        error("""MPI Library is not found.
- On Windows, you can download Microsoft MPI (https://docs.microsoft.com/en-us/message-passing-interface/microsoft-mpi)
- On Unix, you can install OpenMPI via `install_openmpi()`
""")
    end
    if Sys.iswindows()
        return joinpath(ENV["MPI_BIN"], "mpiexec.exe")
    else
        if haskey(ENV, "MPI_INCLUDE_PATH") && haskey(ENV, "MPI_C_LIBRARIES")
            @warn("You are not using the default MPI. Trying to detect the executable...")
            mpirun = abspath(joinpath(ENV["MPI_INCLUDE_PATH"], "..", "bin", "mpirun"))
            if !isfile(mpirun)
                error("Failed.")
            else
                return mpirun
            end
        end

        mpirun =  joinpath(ADCME.BINDIR, "mpirun")
        if !isfile(mpirun)
            error("Failed.")
        end
        return mpirun
    end
end

"""
    independent(o::PyObject, args...; kwargs...)

Returns `o` but when computing the gradients, the top gradients will not be back-propagated into dependent variables of `o`.
"""
independent(o::PyObject, args...; kwargs...) = stop_gradient(o, args...; kwargs...)

@deprecate stop_gradient independent
function stop_gradient(o::PyObject, args...;kwargs...)
    tf.stop_gradient(o, args...;kwargs...)
end

macro cpu(device_id, expr)
    device = "/cpu:"*string(device_id)
    quote
        @pywith tf.device($device) begin
            $(esc(expr))
        end
    end
end

macro cpu(expr)
    device = "/cpu:0"
    quote
        @pywith tf.device($device) begin
            $(esc(expr))
        end
    end
end

macro gpu(device_id, expr)
    device = "/gpu:"*string(device_id)
    quote
        @pywith tf.device($device) begin
            $(esc(expr))
        end
    end
end

macro gpu(expr)
    device = "/gpu:0"
    quote
        @pywith tf.device($device) begin
            $(esc(expr))
        end
    end
end


text:
@testset "1601.maximum-number-of-achievable-transfer-requests.jl" begin
    @test maximum_requests(
        20,
        [
            [7, 11],
            [6, 4],
            [12, 15],
            [8, 3],
            [13, 12],
            [14, 19],
            [5, 16],
            [5, 6],
            [15, 0],
            [6, 5],
            [6, 0],
            [12, 15],
            [18, 8],
            [15, 5],
            [3, 16],
            [18, 13],
        ],
    ) == 2
    @test maximum_requests(5, [[0, 1], [1, 0], [0, 1], [1, 2], [2, 0], [3, 4]]) == 5
    @test maximum_requests(3, [[0, 0], [1, 2], [2, 1]]) == 3
    @test maximum_requests(4, [[0, 3], [3, 1], [1, 2], [2, 0]]) == 4
    @test maximum_requests(
        20,
        [
            [0, 7],
            [2, 0],
            [2, 15],
            [5, 6],
            [6, 16],
            [6, 16],
            [6, 17],
            [7, 14],
            [11, 2],
            [13, 2],
            [14, 11],
            [15, 6],
            [16, 5],
            [16, 6],
            [17, 13],
            [19, 10],
        ],
    ) == 15
end


text:
"""
    FerriteViz.update!(plotter::MakiePlotter, u::Vector)
Updates the Observable `plotter.u` and thereby, triggers the plot to update.
"""
function update!(plotter::MakiePlotter, u::Vector)
    @assert length(plotter.u[]) == length(u)
    plotter.u[] .= u
    Makie.notify(plotter.u)
end

function Makie.convert_arguments(P::Type{<:Makie.Mesh}, plotter::MakiePlotter)
    return Makie.convert_arguments(P,plotter.physical_coords,plotter.triangles)
end

"""
    solutionplot(plotter::MakiePlotter; kwargs...)
    solutionplot(dh::AbstractDofHandler, u::Vector; kwargs...)
    solutionplot!(plotter::MakiePlotter; kwargs...)
    solutionplot!(dh::AbstractDofHandler, u::Vector; kwargs...)
Solutionplot produces the classical contour plot onto the finite element mesh. Most important
keyword arguments are:

- `field::Symbol=:default` representing the field which gets plotted, defaults to the first field in the `dh`.
- `deformation_field::Symbol=:default` field that transforms the mesh by the given deformation, defaults to no deformation
- `process::Function=postprocess` function to construct nodal scalar values from a vector valued problem
- `colormap::Symbol=:cividis`
- `deformation_scale=1.0`
- `shading=false`
- `scale_plot=false`
- `transparent=false`
"""
@recipe(SolutionPlot) do scene
    Attributes(
    scale_plot=false,
    shading=false,
    field=:default,
    deformation_field=:default,
    process=postprocess,
    colormap=:cividis,
    colorrange=(0,1),
    transparent=false,
    deformation_scale = 1.0,
    )
end

function Makie.plot!(SP::SolutionPlot{<:Tuple{<:MakiePlotter}})
    plotter = SP[1][]
    solution = @lift($(SP[:field])===:default ? reshape(transfer_solution(plotter,$(plotter.u); field_idx=1, process=$(SP[:process])), num_vertices(plotter)) : reshape(transfer_solution(plotter,$(plotter.u); field_idx=Ferrite.find_field(plotter.dh,$(SP[:field])), process=$(SP[:process])), num_vertices(plotter)))
    u_matrix = @lift($(SP[:deformation_field])===:default ? zeros(0,3) : transfer_solution(plotter,$(plotter.u); field_idx=Ferrite.find_field(plotter.dh,$(SP[:deformation_field])), process=identity))
    coords = @lift($(SP[:deformation_field])===:default ? plotter.physical_coords : plotter.physical_coords .+ ($(SP[:deformation_scale]) .* $(u_matrix)))
    mins = @lift(minimum($solution))
    maxs = @lift(maximum($solution))
    SP[:colorrange] = @lift(isapprox($mins,$maxs) ? (0,1e-8) : ($mins,$maxs))
    return Makie.mesh!(SP, coords, plotter.triangles, color=solution, shading=SP[:shading], scale_plot=SP[:scale_plot], colormap=SP[:colormap], transparent=SP[:transparent])
end

"""
    cellplot(plotter::MakiePlotter,σ::Vector{T}; kwargs...) where T
    cellplot!(plotter::MakiePlotter,σ::Vector{T}; kwargs...) where T
`cellplot` plots constant scalar data on the cells of the finite element mesh. If `T` is not a number, the keyword argument `process`
can be passed in order to reduce the elements of `σ` to a scalar.

keyword arguments are:

- `deformation_field::Symbol=:default` field that transforms the mesh by the given deformation, defaults to no deformation
- `process::Function=identity` function to construct cell scalar values. Defaults to `identity`, i.e. scalar values.
- `colormap::Symbol=:cividis`
- `deformation_scale=1.0`
- `shading=false`
- `scale_plot=false`
- `transparent=false`
"""
@recipe(CellPlot) do scene
    Attributes(
    scale_plot=false,
    shading=false,
    deformation_field=:default,
    process=identity,
    colormap=:cividis,
    colorrange=(0,1),
    transparent=false,
    deformation_scale = 1.0,
    )
end

function Makie.plot!(CP::CellPlot{<:Tuple{<:MakiePlotter{dim},Vector}}) where dim
    plotter = CP[1][]
    qp_values = CP[2][]
    u_matrix = @lift($(CP[:deformation_field])===:default ? zeros(0,3) : transfer_solution(plotter,$(plotter.u); field_idx=Ferrite.find_field(plotter.dh,$(CP[:deformation_field])), process=identity))
    coords = @lift($(CP[:deformation_field])===:default ? plotter.physical_coords : plotter.physical_coords .+ ($(CP[:deformation_scale]) .* $(u_matrix)))
    mins = minimum(qp_values)
    maxs = maximum(qp_values)
    CP[:colorrange] = isapprox(mins,maxs) ? (0,1e-8) : (mins,maxs)
    solution =  @lift(reshape(transfer_scalar_celldata(plotter, qp_values; process=$(CP[:process])), num_vertices(plotter)))
    return Makie.mesh!(CP, coords, plotter.triangles, color=solution, shading=CP[:shading], scale_plot=CP[:scale_plot], colormap=CP[:colormap], transparent=CP[:transparent])
end

"""
    wireframe(plotter::MakiePlotter; kwargs...)
    wireframe(dh::AbstractDofHandler, u::Vector; kwargs...)
    wireframe(grid::AbstractGrid; kwargs...)
    wireframe!(plotter::MakiePlotter; kwargs...)
    wireframe!(dh::AbstractDofHandler, u::Vector; kwargs...)
    wireframe!(grid::AbstractGrid; kwargs...)
Plots the finite element mesh, optionally labels it and transforms it if a suitable `deformation_field` is given.

- `plotnodes::Bool=true` plots the nodes as circles/spheres
- `strokewidth::Int=2` how thick faces/edges are drawn
- `color::Symbol=theme(scene,:linecolor)` color of the faces/edges and nodes
- `markersize::Int=30` size of the nodes
- `deformation_field::Symbol=:default` field that transforms the mesh by the given deformation, defaults to no deformation
- `deformation_scale::Number=1.0` scaling of the deformation
- `cellsets=false` Color cells based on their cellset association. If no cellset is found for a cell, the cell is marked blue.
- `nodelables=false` global node id labels
- `nodelabelcolor=:darkblue`
- `celllabels=false` global cell id labels
- `celllabelcolor=:darkred`
- `textsize::Int=15` size of the label's text
- `visible=true`
"""
@recipe(Wireframe) do scene
    Attributes(
    plotnodes=true,
    strokewidth=2,
    color=theme(scene, :linecolor),
    markersize=30,
    deformation_field=:default,
    visible=true,
    deformation_scale=1,
    textsize=15,
    offset=(0.0,0.0),
    nodelabels=false,
    nodelabelcolor=:darkblue,
    celllabels=false,
    celllabelcolor=:darkred,
    cellsets=false,
    )
end

function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
    plotter = WF[1][]
    #triangle representation
    # can't use triangle representation, since we don't know by this information which edges to draw
    # further it would be in the incompressible example 2600 nodes vs 15000 in triangle representation
    # u_matrix = @lift($(WF[:deformation_field])===:default ? zeros(0,3) : transfer_solution(plotter; field_idx=Ferrite.find_field(plotter.dh,$(WF[:deformation_field])), process=identity))
    # coords = @lift($(WF[:deformation_field])===:default ? plotter.physical_coords : plotter.physical_coords .+ ($(WF[:scale]) .* $(u_matrix)))
    #original representation
    nodal_u_matrix = @lift($(WF[:deformation_field])===:default ? zeros(0,3) : dof_to_node(plotter.dh, $(WF[1][].u); field=Ferrite.find_field(plotter.dh,$(WF[:deformation_field])), process=identity))
    gridnodes = @lift($(WF[:deformation_field])===:default ? plotter.gridnodes : plotter.gridnodes .+ ($(WF[:deformation_scale]) .* $(nodal_u_matrix)))
    lines = @lift begin
        dim > 2 ? (lines = Point3f[]) : (lines = Point2f[])
        for cell in Ferrite.getcells(plotter.dh.grid)
            boundaryentities = dim < 3 ? Ferrite.faces(cell) : Ferrite.edges(cell)
            append!(lines, [$gridnodes[e,:] for boundary in boundaryentities for e in boundary])
        end
        lines
    end
    nodes = @lift($(WF[:plotnodes]) ? $(gridnodes) : zeros(Float32,0,3))
    #plot cellsets
    cellsets = plotter.dh.grid.cellsets
    cellset_to_value = Dict{String,Int}()
    for (cellsetidx,(cellsetname,cellset)) in enumerate(cellsets)
        cellset_to_value[cellsetname] = cellsetidx
    end
    cellset_u = zeros(Ferrite.getncells(plotter.dh.grid))
    for (cellidx,cell) in enumerate(Ferrite.getcells(plotter.dh.grid))
        for (cellsetname,cellsetvalue) in cellset_to_value
            if cellidx in cellsets[cellsetname]
                cellset_u[cellidx] = cellsetvalue
            end
        end
    end
    u_matrix = @lift($(WF[:deformation_field])===:default ? zeros(0,3) : transfer_solution(plotter,$(plotter.u); field_idx=Ferrite.find_field(plotter.dh,$(WF[:deformation_field])), process=identity))
    coords = @lift($(WF[:deformation_field])===:default ? plotter.physical_coords : plotter.physical_coords .+ ($(WF[:deformation_scale]) .* $(u_matrix)))
    colorrange = isempty(cellset_to_value) ? (0,1) : (0,maximum(values(cellset_to_value)))
    cellset_u =  reshape(transfer_scalar_celldata(plotter, cellset_u; process=identity), num_vertices(plotter))
    Makie.mesh!(WF, coords, plotter.triangles, color=cellset_u, shading=false, scale_plot=false, colormap=:darktest, visible=WF[:cellsets])
    #plot the nodes
    Makie.scatter!(WF,gridnodes,markersize=WF[:markersize], color=WF[:color], visible=WF[:visible])
    #set up nodelabels
    nodelabels = @lift $(WF[:nodelabels]) ? ["$i" for i in 1:size($gridnodes,1)] : [""]
    nodepositions = @lift $(WF[:nodelabels]) ? [dim < 3 ? Point2f(row) : Point3f(row) for row in eachrow($gridnodes)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
    #set up celllabels
    celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:Ferrite.getncells(plotter.dh.grid)] : [""]
    cellpositions = @lift $(WF[:celllabels]) ? [midpoint(cell,$gridnodes) for cell in Ferrite.getcells(plotter.dh.grid)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
    Makie.text!(WF,nodelabels, position=nodepositions, textsize=WF[:textsize], offset=WF[:offset],color=WF[:nodelabelcolor])
    Makie.text!(WF,celllabels, position=cellpositions, textsize=WF[:textsize], color=WF[:celllabelcolor], align=(:center,:center))
    #plot edges (3D) /faces (2D) of the mesh
    Makie.linesegments!(WF,lines,color=WF[:color], linewidth=WF[:strokewidth], visible=WF[:visible])
end

function Makie.plot!(WF::Wireframe{<:Tuple{<:Ferrite.AbstractGrid{dim}}}) where dim
    grid = WF[1][]
    coords = [Ferrite.getcoordinates(node)[i] for node in Ferrite.getnodes(grid), i in 1:dim]
    dim > 2 ? (lines = Point3f[]) : (lines = Point2f[])
    for cell in Ferrite.getcells(grid)
        boundaryentities = dim < 3 ? Ferrite.faces(cell) : Ferrite.edges(cell)
        append!(lines, [coords[e,:] for boundary in boundaryentities for e in boundary])
    end
    nodes = @lift($(WF[:plotnodes]) ? coords : zeros(Float32,0,3))
    Makie.scatter!(WF,nodes,markersize=WF[:markersize], color=WF[:color])
    nodelabels = @lift $(WF[:nodelabels]) ? ["$i" for i in 1:size(coords,1)] : [""]
    nodepositions = @lift $(WF[:nodelabels]) ? [dim < 3 ? Point2f(row) : Point3f(row) for row in eachrow(coords)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
    celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:Ferrite.getncells(grid)] : [""]
    cellpositions = @lift $(WF[:celllabels]) ? [midpoint(cell,coords) for cell in Ferrite.getcells(grid)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
    #cellsetsplot
    dh = Ferrite.DofHandler(grid)
    cellsets = grid.cellsets
    cellset_to_value = Dict{String,Int}()
    for (cellsetidx,(cellsetname,cellset)) in enumerate(cellsets)
        cellset_to_value[cellsetname] = cellsetidx
    end
    cellset_u = zeros(Ferrite.getncells(grid))
    for (cellidx,cell) in enumerate(Ferrite.getcells(grid))
        for (cellsetname,cellsetvalue) in cellset_to_value
            if cellidx in cellsets[cellsetname]
                cellset_u[cellidx] = cellsetvalue
            end
        end
    end
    plotter = MakiePlotter(dh,cellset_u)
    cellset_u =  reshape(transfer_scalar_celldata(plotter, cellset_u; process=identity), num_vertices(plotter))
    colorrange = isempty(cellset_to_value) ? (0,1) : (0,maximum(values(cellset_to_value)))
    Makie.mesh!(WF, plotter.physical_coords, plotter.triangles, color=cellset_u, shading=false, scale_plot=false, colormap=:darktest, visible=WF[:cellsets])
    Makie.text!(WF,nodelabels, position=nodepositions, textsize=WF[:textsize], offset=WF[:offset],color=WF[:nodelabelcolor])
    Makie.text!(WF,celllabels, position=cellpositions, textsize=WF[:textsize], color=WF[:celllabelcolor], align=(:center,:center))
    Makie.linesegments!(WF,lines,color=WF[:color], strokewidth=WF[:strokewidth])
end

"""
    surface(plotter::MakiePlotter; kwargs...)
    surface(dh::AbstractDofHandler, u::Vector; kwargs...)
    surface!(plotter::MakiePlotter; kwargs...)
    surface!(dh::AbstractDofHandler, u::Vector; kwargs...)
Uses the given `field` and plots the scalar values as a surface. If it's a vector valued problem, the nodal vector
values are transformed to a scalar based on `process` which defaults to the magnitude. Only availble in `dim=2`.

- `field = :default`
- `process = postprocess`
- `scale_plot = false`
- `shading = false`
- `colormap = :cividis`
"""
@recipe(Surface) do scene
    Attributes(
    field = :default,
    process = postprocess,
    scale_plot = false,
    shading = false,
    colormap = :cividis,
    )
end

function Makie.plot!(SF::Surface{<:Tuple{<:MakiePlotter{2}}})
    plotter = SF[1][]
    field = @lift($(SF[:field])===:default ? 1 : Ferrite.find_field(plotter.dh,$(SF[:field])))
    solution = @lift(reshape(transfer_solution(plotter, $(plotter.u); field_idx=$(field), process=$(SF[:process])), num_vertices(plotter)))
    points = @lift([Point3f(coord[1], coord[2], $(solution)[idx]) for (idx, coord) in enumerate(eachrow(plotter.physical_coords))])
    return Makie.mesh!(SF,points, plotter.triangles, color=solution, scale_plot=SF[:scale_plot], shading=SF[:shading], colormap=SF[:colormap])
end

"""
    arrows(plotter::MakiePlotter; kwargs...)
    arrows(dh::AbstractDofHandler, u::Vector; kwargs...)
    arrows!(plotter::MakiePlotter; kwargs...)
    arrows!(dh::AbstractDofHandler, u::Vector; kwargs...)
At every node position a arrows is drawn, where the arrow tip ends at the node. Only works in `dim >=2`. If a `color` is specified
the arrows are unicolored. Otherwise the color corresponds to the magnitude, or any other scalar value based on the `process` function.

- `arrowsize = 0.08`
- `normalize = true`
- `field = :default`
- `color = :default`
- `colormap = :cividis`
- `process=postprocess`
- `lengthscale = 1f0`
"""
@recipe(Arrows) do scene
    Attributes(
    arrowsize = 0.08,
    normalize = true,
    field = :default,
    color = :default,
    colormap = :cividis,
    process=postprocess,
    lengthscale = 1f0,
    )
end

function Makie.plot!(AR::Arrows{<:Tuple{<:MakiePlotter{dim}}}) where dim
    plotter = AR[1][]
    field = @lift($(AR[:field])===:default ? 1 : Ferrite.find_field(plotter.dh,$(AR[:field])))
    @assert Ferrite.getfielddim(plotter.dh,field[]) > 1
    solution = @lift(transfer_solution(plotter, $(plotter.u); field_idx=$(field), process=identity))
    if dim  == 2
        ps = [Point2f(i) for i in eachrow(plotter.physical_coords)]
        ns = @lift([Vec2f(i) for i in eachrow($(solution))])
        lengths = @lift($(AR[:color])===:default ? $(AR[:process]).($(ns)) : ones(length($(ns)))*$(AR[:color]))
    elseif dim  == 3
        ps = [Point3f(i) for i in eachrow(plotter.physical_coords)]
        ns = @lift([Vec3f(i) for i in eachrow($(solution))])
        lengths = @lift($(AR[:color])===:default ? $(AR[:process]).($(ns)) : ones(length($(ns)))*$(AR[:color]))
    else
        error("Arrows plots are only available in dim ≥ 2")
    end
    Makie.arrows!(AR, ps, ns, arrowsize=AR[:arrowsize], normalize=AR[:normalize], colormap=AR[:colormap], color=lengths, lengthscale=AR[:lengthscale])
end

"""
    ferriteviewer(plotter::MakiePlotter)
    ferriteviewer(plotter::MakiePlotter, u_history::Vector{Vector{T}}})
Constructs a viewer with a `solutionplot`, `Colorbar` as well as sliders,toggles and menus to change the current view.
If the second dispatch is called a timeslider is added, in order to step through a set of solutions obtained from a simulation.
"""
function ferriteviewer(plotter::MakiePlotter{dim}) where dim
    #set up figure and axis, axis either LScene in 3D or Axis if dim < 2
    fig = Figure()
    dim > 2 ? (ax = LScene(fig[1,1])) : (ax = Axis(fig[1,1]))
    #toggles and their labels for switching plot types on/off
    toggles = [Toggle(fig, active=active) for active in [true,true,false]]
    labels = [Label(fig,label) for label in ["mesh", "deformation", "labels"]]
    #setup the deformation_field as Node (Observable)
    deformation_field = Node(Ferrite.getfieldnames(plotter.dh)[1])
    #solutionplot main plot of the viewer
    solutionp = solutionplot!(plotter,colormap=:cividis,deformation_field=@lift $(toggles[2].active) ? $(deformation_field) : :default)

    #setting up various sliders
    markerslider = Slider(fig, range = 0:1:100, startvalue=5)
    strokewidthslider = Slider(fig, range = 0:1:10, startvalue=1)
    markersize = lift(x->x,markerslider.value)
    strokewidth = lift(x->x,strokewidthslider.value)

    #plot the fe-mesh
    wireframep = wireframe!(plotter,markersize=markersize,strokewidth=strokewidth,deformation_field= @lift $(toggles[2].active) ? $(deformation_field) : :default)
    #connect fe-mesh plot to the toggle
    connect!(wireframep.visible,toggles[1].active)
    connect!(wireframep.nodelabels,toggles[3].active)
    connect!(wireframep.celllabels,toggles[3].active)

    #set up dropdown menus for colormap, field, deformation field and processing function
    menu_cm = Menu(fig, options=["cividis", "inferno", "thermal"],label="colormap", direction=:up)
    menu_field = Menu(fig, options=Ferrite.getfieldnames(plotter.dh))
    menu_deformation_field = Menu(fig, options=Ferrite.getfieldnames(plotter.dh))
    menu_process = Menu(fig, options=[x₁,x₂,x₃,("magnitude",l2),("manhatten norm",l1),identity])
    #align all menus as a vgrid under each other
    fig[1,3] = vgrid!(grid!(hcat(toggles,labels), tellheight=false),
                      Label(fig,"nodesize",width=nothing), markerslider,
                      Label(fig,"strokewidth",width=nothing), strokewidthslider,
                      Label(fig,"processing function",width=nothing), menu_process,
                      Label(fig,"field",width=nothing), menu_field,
                      Label(fig, "deformation field",width=nothing),menu_deformation_field,
                      Label(fig, "colormap",width=nothing),menu_cm)
    #add colorbar
    cb = Colorbar(fig[1,2], solutionp)

    #event handling for selecting stuff from the menus
    on(menu_cm.selection) do s
        cb.colormap = s
        solutionp.colormap = s
    end

    on(menu_field.selection) do field
        solutionp.field = field
    end

    on(menu_deformation_field.selection) do field
        solutionp.deformation_field = field
        wireframep.deformation_field = field
    end

    on(menu_process.selection) do process_function
        solutionp.process=process_function
    end

    return fig
end

function ferriteviewer(plotter::MakiePlotter, data::Vector{Vector{T}}) where T
    fig = ferriteviewer(plotter)
    timeslider = labelslider!(fig, "timestep n:", 1:length(data); format = x->"$x", sliderkw = Dict(:snap=>false))
    fig[2,1] = timeslider.layout
    @lift(FerriteViz.update!(plotter,data[$(timeslider.slider.value)]))
    display(fig)
end

####### One Shot Methods #######
const FerriteVizPlots = Union{Type{<:Wireframe},Type{<:SolutionPlot},Type{<:Arrows}}

function Makie.convert_arguments(P::FerriteVizPlots, dh::Ferrite.AbstractDofHandler, u::Vector)
    return (MakiePlotter(dh,u),)
end


text:
""" Abstract supertype for all RBMs """
abstract type AbstractRBM end

"""
Abstract super type for RBMs with binary and Bernoulli distributed hidden nodes.
"""
abstract type AbstractXBernoulliRBM <: AbstractRBM end

"""
    BernoulliRBM(weights, visbias, hidbias)
Encapsulates the parameters of an RBM with Bernoulli distributed nodes.
* `weights`: matrix of weights with size
  (number of visible nodes, number of hidden nodes)
* `visbias`: bias vector for visible nodes
* `hidbias`: bias vector for hidden nodes
"""
struct BernoulliRBM <: AbstractXBernoulliRBM
   weights::Array{Float64,2}
   visbias::Array{Float64,1}
   hidbias::Array{Float64,1}
end

"""
    GaussianBernoulliRBM(weights, visbias, hidbias, sd)
Encapsulates the parameters of an RBM with Gaussian distributed visible nodes
and Bernoulli distributed hidden nodes.
"""
struct GaussianBernoulliRBM <: AbstractXBernoulliRBM
   weights::Array{Float64,2}
   visbias::Array{Float64,1}
   hidbias::Array{Float64,1}
   sd::Array{Float64,1}
end

"""
    GaussianBernoulliRBM2(weights, visbias, hidbias, sd)
Encapsulates the parameters of an RBM with Gaussian distributed visible nodes
and Bernoulli distributed hidden nodes with the alternative energy formula
proposed by KyungHyun Cho.
"""
struct GaussianBernoulliRBM2 <: AbstractXBernoulliRBM
   weights::Array{Float64,2}
   visbias::Array{Float64,1}
   hidbias::Array{Float64,1}
   sd::Array{Float64,1}
end

"""
    BernoulliGaussianRBM(weights, visbias, hidbias)
Encapsulates the parameters of an RBM with Bernoulli distributed visible nodes
and Gaussian distributed hidden nodes.
The standard deviation of the Gaussian distribution is 1.
"""
struct BernoulliGaussianRBM <: AbstractRBM
   weights::Array{Float64,2}
   visbias::Array{Float64,1}
   hidbias::Array{Float64,1}
end

"""
    Binomial2BernoulliRBM(weights, visbias, hidbias)
Encapsulates the parameters of an RBM with 0/1/2-valued, Binomial (n=2)
distributed visible nodes, and Bernoulli distributed hidden nodes.
This model is equivalent to a BernoulliRBM in which every two visible nodes are
connected with the same weights to each hidden node.
The states (0,0) / (1,0) / (0,1) / (1,1) of the visible nodes connected with
with the same weights translate as states 0 / 1 / 1 / 2 in the
Binomial2BernoulliRBM.
"""
struct Binomial2BernoulliRBM <: AbstractXBernoulliRBM
   weights::Matrix{Float64}
   visbias::Vector{Float64}
   hidbias::Vector{Float64}
end

"""
NegativeBinomialBernoulliRBM(weights, visbias, hidbias)
Encapsulates the parameters of an RBM with Negative Binomial distributed visible nodes
and Bernoulli distributed hidden nodes.
"""
struct NegativeBinomialBernoulliRBM <: AbstractXBernoulliRBM
   weights::Array{Float64,2}
   visbias::Array{Float64,1}
   hidbias::Array{Float64,1}
   inversedispersion::AbstractArray
   dropoutmid::Array{Float64,1}
   dropoutshape::Array{Float64,1}
end

"""
    ranges(numbers)
Returns a vector of consecutive integer ranges, the first starting with 1.
The i'th such range spans over `numbers[i]` items.
"""
function ranges(numbers::Vector{Int})
   ranges = Vector{UnitRange{Int}}(undef, length(numbers))
   offset = 0
   for i in eachindex(numbers)
      ranges[i] = offset .+ (1:numbers[i])
      offset += numbers[i]
   end
   ranges
end

"""
    PartitionedRBM(rbms)
Encapsulates several (parallel) AbstractRBMs that form one partitioned RBM.
The nodes of the parallel RBMs are not connected between the RBMs.
"""
struct PartitionedRBM{R<:AbstractRBM} <: AbstractRBM
   rbms::Vector{R}
   visranges::Vector{UnitRange{Int}}
   hidranges::Vector{UnitRange{Int}}

   function PartitionedRBM{R}(rbms::Vector{R}) where R
      visranges = ranges([length(rbm.visbias) for rbm in rbms])
      hidranges = ranges([length(rbm.hidbias) for rbm in rbms])
      new(rbms, visranges, hidranges)
   end
end

""" Singleton-Placeholder for `AbstractRBM`s """
struct NoRBM <: AbstractRBM
end

const BasicDBM = Vector{BernoulliRBM}

"A DBM with only Bernoulli distributed nodes which may contain partitioned layers."
const PartitionedBernoulliDBM =
      Vector{<:Union{BernoulliRBM, PartitionedRBM{BernoulliRBM}}}

const MultimodalDBM = Vector{<:AbstractRBM}

const AbstractBM = Union{MultimodalDBM, AbstractRBM}


text:
import ReinforcementLearningBase: reset!
function testlearn()
    mdp = MDP()
    learner = Sarsa()
    x = RLSetup(
        learner = Sarsa(),
        environment = mdp,
        callbacks = [TotalReward(), RecordAll()],
        stoppingcriterion = ConstantNumberSteps(10),
    )
    seed!(13452)
    reset!(mdp)
    learn!(x)
    learn!(x)

    x2 = RLSetup(
        learner = Sarsa(),
        environment = mdp,
        callbacks = [TotalReward(), RecordAll()],
        stoppingcriterion = ConstantNumberSteps(20),
    )
    seed!(13452)
    reset!(mdp)
    learn!(x2)
    @test x.learner.params == x2.learner.params
end
testlearn()


text:
using NestedSamplers.Bounds: Ellipsoid, MultiEllipsoid, fit, scale!, decompose, volume, volume_prefactor

const NMAX = 20

@testset "ndims=$N" for N in 1:NMAX
    @testset "Spheres" begin
        scale = 5
        center = 2scale .* ones(N)
        A = Matrix(I, N, N) ./ scale^2
        ell = Ellipsoid(center, A)
        @test volume(ell) ≈ volume_prefactor(N) * scale^N
        axs, axlens = decompose(ell)
        @test axlens ≈ fill(scale, N)
        @test axs ≈ Bounds.axes(ell) ≈ scale * Matrix(I, N, N)
    end

    @testset "Scaling" begin
        scale = 1.5
        center = zeros(N)
        A = diagm(0 => rand(rng, N))
        ell = Ellipsoid(center, A)

        ell2 = Ellipsoid(center, A ./ scale^2)

        scale!(ell, scale^N)

        @test volume(ell) ≈ volume(ell2)
        @test ell.A ≈ ell2.A
        @test all(decompose(ell) .≈ decompose(ell2))
    end

    @testset "Contains" begin
        E = 1e-7
        ell = Ellipsoid(N)

        # Point just outside unit n-Spheres
        pt = (1 / √N + E) .* ones(N)
        @test pt ∉ ell

        # point just inside
        pt = (1 / √N - E) .* ones(N)
        @test pt ∈ ell

        A = diagm(0 => rand(rng, N))
        ell = Ellipsoid(zeros(N), A)

        for i in 1:N
            axlen = 1 / sqrt(A[i, i])
            pt = zeros(N)
            pt[i] = axlen + E
            @test pt ∉ ell
            pt[i] = axlen - E
            @test pt ∈ ell
        end
    end
end # testset