text:
using Oceananigans
using Oceananigans.Utils: prettytime
using JLD2
using GLMakie

function variable_range(file, variable)
    max_var = -Inf
    min_var = Inf
    for i in collect(parse.(Int, keys(file["timeseries/t"])))
        max_var = max(max_var, maximum(file["timeseries/$(variable)/" * string(i)][:, :, 1]))
        min_var = min(max_var, minimum(file["timeseries/$(variable)/" * string(i)][:, :, 1]))
    end
    return (min_var * 0.5, max_var * 0.5)
end

function visualize_spherical_field(filepath, variable, range = nothing)

    @show output_prefix = basename(filepath)[1:end-5]

    file = jldopen(filepath)

    iterations = parse.(Int, keys(file["timeseries/t"]))

    if range isa Nothing
        range = variable_range(file, variable)
    end
    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]

    grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                 radius = 1,
                                 latitude = (-80, 80),
                                 longitude = (-180, 180),
                                 z = (-1, 0))

    λ = xnodes(Face, grid)
    φ = ynodes(Center, grid)

    λ = repeat(reshape(λ, Nx, 1), 1, Ny)
    φ = repeat(reshape(φ, 1, Ny), Nx, 1)

    λ_azimuthal = λ .+ 180  # Convert to λ ∈ [0°, 360°]
    φ_azimuthal = 90 .- φ   # Convert to φ ∈ [0°, 180°] (0° at north pole)

    iter = Observable(0)

    plot_title = @lift "vorticity at: time = $(prettytime(file["timeseries/t/" * string($iter)]))"

    var = @lift file["timeseries/$(variable)/" * string($iter)][:, :, 1]

    # Plot on the unit sphere to align with the spherical wireframe.
    x = @. cosd(λ_azimuthal) * sind(φ_azimuthal)
    y = @. sind(λ_azimuthal) * sind(φ_azimuthal)
    z = @. cosd(φ_azimuthal)

    fig = Figure(resolution = (2000, 2000))

    fontsize_theme = Theme(fontsize = 25)
    set_theme!(fontsize_theme)

    ax = fig[1, 1] = LScene(fig)
    wireframe!(ax, Sphere(Point3f0(0), 0.99f0), show_axis=false)
    hm = surface!(ax, x, y, z, color=var, colormap=:thermal, colorrange=range)
    rotate_cam!(ax.scene, (π/5, π/6, 0))

    cb = Colorbar(fig[1, 2], hm, label=plot_title)
    cb.height = Relative(2/3)
    cb.width = 20

    record(fig, output_prefix * ".mp4", iterations, framerate=80) do i
        @info "Animating iteration $i/$(iterations[end])..."
        iter[] = i
    end

    return nothing
end

function visualize_cartesian_field(filepath, variable, range = nothing)

    @show output_prefix = basename(filepath)[1:end-5]

    file = jldopen(filepath)

    iterations = parse.(Int, keys(file["timeseries/t"]))

    if range isa Nothing
        range = variable_range(file, variable)
    end
    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]

    grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                 radius = 1,
                                 latitude = (-80, 80),
                                 longitude = (-180, 180),
                                 z = (-1, 0))

    λ = xnodes(Face, grid)
    φ = ynodes(Center, grid)

    iter = Observable(0)

    plot_title = @lift "vorticity at: time = $(prettytime(file["timeseries/t/" * string($iter)]))"

    var = @lift file["timeseries/$(variable)/" * string($iter)][:, :, 1]

    fig = Figure(resolution = (2000, 2000))

    fontsize_theme = Theme(fontsize = 25)
    set_theme!(fontsize_theme)

    ax = fig[1, 1] = LScene(fig)
    hm = heatmap!(ax, λ, φ, var, colormap=:thermal, colorrange=range)

    cb = Colorbar(fig[1, 2], hm, label=plot_title)
    cb.height = Relative(2/3)
    cb.width = 20

    record(fig, output_prefix * ".mp4", iterations, framerate=80) do i
        @info "Animating iteration $i/$(iterations[end])..."
        iter[] = i
    end

    return nothing
end

function visualize_cartesian_comparison(filepath1, filepath2, title1, title2, variable)

    @show output_prefix = "comparison"

    file1 = jldopen(filepath1)
    file2 = jldopen(filepath2)

    it1 = parse.(Int, keys(file1["timeseries/t"]))
    it2 = parse.(Int, keys(file2["timeseries/t"]))

    frames = collect(1:length(it1))

    iter = Observable(1)

    var1 = @lift file1["timeseries/$(variable)/" * string(it1[$iter])][:, :, 1]
    var2 = @lift file2["timeseries/$(variable)/" * string(it2[$iter])][:, :, 1]

    fig = Figure(resolution = (3500, 1000))

    ax = Axis(fig[1, 1], title = title1)
    hm = heatmap!(ax, var1, colormap=:balance, colorrange=(-10, 10))

    ax = Axis(fig[1, 3], title = title2)
    hm = heatmap!(ax, var2, colormap=:balance, colorrange=(-10, 10))

    record(fig, output_prefix * ".mp4", frames, framerate=10) do i
        @info "Animating iteration $i/$(frames[end])..."
        iter[] = i
    end

    return nothing
end


text:
using SafeTestsets

@safetestset "Blocks: math" begin include("test_math.jl") end
@safetestset "Blocks: nonlinear" begin include("test_nonlinear.jl") end
@safetestset "Blocks: continuous" begin include("test_continuous.jl") end
@safetestset "Analog Circuits" begin include("analog.jl") end
#@safetestset "Digital Circuits" begin include("digital.jl") end
@safetestset "RC Circuit Demo" begin include("demo.jl") end
@safetestset "Thermal Circuits" begin include("thermal.jl") end

ext: jl


text:
function readAttritionSchemes( mpSim::MPsim, catXF::XF )

    if !( XLSX.hassheet( catXF, "General" ) &&
        XLSX.hassheet( catXF, "Attrition" ) )
        @warn "Catalogue file is missing a 'General' and 'Attrition' sheet, cannot confiugre attrition schemes."
        return
    end  # if !( XLSX.hassheet( catXF, "General" ) && ...

    ws = catXF[ "General" ]
    nAttrition = ws[ "B4" ] isa Integer ? ws[ "B4" ] : 0
    ws = catXF[ "Attrition" ]
    clearSimulationAttrition!( mpSim )

    attritionSchemes = readAttritionScheme.( Ref( ws ), (1:nAttrition) .+ 1 )
    setSimulationAttrition!( mpSim, attritionSchemes )

end  # readAttritionSchemes( mpSim, catXF )

function readAttritionScheme( ws::WS, sLine::Int )

    # Get the name of the scheme.
    name = ws[ sLine, 1 ]
    name = name isa Missing ? "default" : name

    attrition = Attrition( name )

    period = ws[ sLine, 2 ]
    period = period isa Real ? period : 1.0
    setAttritionPeriod!( attrition, period )

    nPoints = ws[ sLine, 3 ]

    if ( nPoints isa Integer ) && ( nPoints > 1 )
        curve = Dict{Float64, Float64}( period * ws[ sLine, 2*ii + 2 ] =>
            ws[ sLine, 2*ii + 3 ] for ii in 1:nPoints )
        setAttritionCurve!( attrition, curve )
    end  # if ( nPoints isa Integer ) && ( nPoints > 1 )

    return attrition

end  # readAttritionScheme( ws, sLine )

function readAttritionSchemes( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    attritionPars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Attrition'" ) ) )

    if isempty( attritionPars )
        return
    end  # if isempty( attritionPars )

    attritionSchemes = map( eachindex( attritionPars[ :, :parName ] ) ) do ii
        attrition = Attrition( attritionPars[ ii, :parName ] )
        pars = split( attritionPars[ ii, :parValue ], ";" )

        setAttritionPeriod!( attrition, parse( Float64, pars[ 1 ] ) )

        curve = split.( split( pars[ 2 ][ 2:(end-1) ], "," ), ":" )
        curveDict = Dict{Float64, Float64}()

        for curvePoint in curve
            curveDict[ parse( Float64, curvePoint[ 1 ] ) ] =
                parse( Float64, curvePoint[ 2 ] )
        end  # for curvePoint in curve

        setAttritionCurve!( attrition, curveDict )

        return attrition
    end  # map( eachindex( attritionPars[ :parName ] ) ) do ii

    setSimulationAttrition!( mpSim, attritionSchemes )

end  # readAttritionSchemes( mpSim, configDB, configName )


text:
using Base.BLAS
using DataFrames
using StatsBase

immutable StackedNet{T<:FloatingPoint}
	layers::Vector{Layer{T}}
	error::ASCIIString
	error_function!::Function
	error_function_prime!::Function

	function StackedNet(units::Vector{Units}; error::ASCIIString="squared_error")
		if length(units) < 2
			return error("StackedNet units specification is too short.")
		end
		if minimum([unit.n for unit in units]) <= 0
			return error("Invalid number of units in StackedNet units specification.")
		end
		# Iterate over sequential paits of units and construct the required layers.
		layers = Array(Layer{T}, length(units)-1)
		for u = 1:length(units)-1
			units1, units2 = units[u], units[u + 1]
			layers[u] = Layer{T}(units1.n, units2.n, units2.activation)
		end
		# Set error function.
		error, error_function!, error_function_prime! = error_function_selector(error)
		# Create and return the object.
		new(layers, error, error_function!, error_function_prime!)
	end
end

# Returns the patternwise error on a specific pattern.
function error!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, Y::Matrix{T}, p::Int)
	forward!(DN, X, p)
	L = DN.layers[end]
	y = sub(Y, :, p)
	DN.error_function!(L.ACT, y, L.E)
	sum(L.E)
end

# Returns the patternwise error on a set of patterns.
function error!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, Y::Matrix{T})
	@inbounds begin
		E::T = 0.0
		for p = 1:size(X, 2)
			E += error!(DN, X, Y, p)
		end
		E = E / size(X, 2)
	end
	E
end

# Forward propagate a pattern through a StackedNet.
function forward!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, p::Int)
	@inbounds begin
		x = sub(X, :, p)
		forward!(DN.layers[1], x)
		# forward!(DN.layers[1], X, p)
		for l = 2:length(DN.layers)
			forward!(DN.layers[l], DN.layers[l - 1].ACT)
		end
	end
end

# Returns output activations in a StackedNet for a set of input patterns.
function forward!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T})
	@inbounds begin
		no = DN.layers[end].no
		np = size(X, 2)
		YH = zeros(T, (no, np))
		for p = 1:np
			forward!(DN, X, p)
			for o = 1:no
				YH[o, p] = DN.layers[end].ACT[o]
			end
		end
	end
	YH
end

# Returns output activations in a StackedNet for a set of input patterns.
function forward!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, YH::Matrix{T})
	@inbounds begin
		for p = 1:np
			forward!(DN, X, p)
			for o = 1:no
				YH[o, p] = DN.layers[end].ACT[o]
			end
		end
	end
end

# Numerically checks the analytic gradient of a StackedNet and returns a data frame of results.
function gradient_check{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, Y::Matrix{T}; step=1e-8, abs_tolerance=1e-5, rel_tolerance=1e-3)
	# Initialize data frame.
	df = DataFrame(layer=Int64[], parameter=Int64[], analytic_grad=T[], numerical_grad=T[], absolute_error=T[], ok=Bool[])
	# Iterate over layers, and parameters within each layer.
	for l = 1:length(DN.layers)
		for (M, GM) in ((DN.layers[l].W, DN.layers[l].GW), (DN.layers[l].B, DN.layers[l].GB))
			for i = 1:length(M)
				current_parameter_value = M[i]
				# Set the gradient information.
				gradient_reset!(DN)
				gradient_update!(DN, X, Y)
				# Get the analytic gradient (scale it to reflect patternwise gradient).
				analytic_gradient = GM[i] / size(X, 2)
				# Get the numerical gradient.
				M[i] = current_parameter_value - step
				EN = error!(DN, X, Y)
				M[i] = current_parameter_value + step
				EP = error!(DN, X, Y)
				numerical_gradient = (EP - EN) / (2.0 * step)
				M[i] = current_parameter_value
				# Compare gradient values.
				# if (analytic_gradient == 0.0) & (numerical_gradient == 0.0)
				# 	error = 0.0
				# else
				# 	error = abs(analytic_gradient - numerical_gradient) / (abs(analytic_gradient) + abs(numerical_gradient))
				# end
				error = abs(analytic_gradient - numerical_gradient)
				ok = error <= abs_tolerance + (rel_tolerance * abs(numerical_gradient)) ? true : false
				push!(df, (l, i, analytic_gradient, numerical_gradient, error, ok))
			end
		end
	end
	df
end

# Finds the maximum absolute gradient value.
function gradient_maxabs{T<:FloatingPoint}(DN::StackedNet{T})
	@inbounds begin
		gmax::T = -Inf
		for l = 1:length(DN.layers)
			for M in (DN.layers[l].GW, DN.layers[l].GB)
				thismaxabs = maximum(abs(M))
				if thismaxabs > gmax
					gmax = thismaxabs
				end
			end
		end
	end
	gmax
end

# Zeros out all the gradient information in a StackedNet.
function gradient_reset!{T<:FloatingPoint}(DN::StackedNet{T})
	@inbounds begin
		for l = 1:length(DN.layers)
			fill!(DN.layers[l].GW, 0.0)
			fill!(DN.layers[l].GB, 0.0)
		end
	end
end

# Increment the gradient information (GW and GB) on each layer based on a single input-output pattern.
function gradient_update!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, Y::Matrix{T}, p::Int)
	@inbounds begin
		# Forward propagate the input pattern through the network.
		forward!(DN, X, p)
		# Walk backward through the stacked network.
		nl = length(DN.layers)
		for l = nl:-1:1
			# Set references to layers.
			Lup = l == nl ? nothing : DN.layers[l + 1]
			L = DN.layers[l]
			Ldn = l == 1 ? nothing : DN.layers[l - 1]
			# Backpropagation.
			if l == nl
				y = sub(Y, :, p)
				DN.error_function_prime!(L.ACT, y, L.DE_DYH)
				backward!(L, L.DE_DYH)
			else
				backward!(L, Lup.DELTAS)
			end
			# Increment gradient information.
			axpy!(L.no, 1.0, L.DE_DNET, 1, L.GB, 1)
			if l == 1
				x = sub(X, :, p)
				ger!(1.0, X[:, p], L.DE_DNET, L.GW)
			else
				ger!(1.0, Ldn.ACT, L.DE_DNET, L.GW)
			end
		end
	end
end

# Increment the gradient information (GW and GB) on each layer based on a set of input-output pairs.
function gradient_update!{T<:FloatingPoint}(DN::StackedNet{T}, X::Matrix{T}, Y::Matrix{T})
	@inbounds begin
		for p = 1:size(X, 2)
			gradient_update!(DN, X, Y, p)
		end
	end
end

# Update the parameters of a StackedNet based on the accumulated gradient information.
function parameters_update!{T<:FloatingPoint}(DN::StackedNet{T}, learning_rate::T; reset_gradient::Bool=true)
	@inbounds begin
		#lr_use = get_learning_rate(DN, lr)
		for l = 1:length(DN.layers)
			L = DN.layers[l]
			for i = 1:length(L.B)
				L.B[i] -= learning_rate * sign(L.GB[i])
			end
			for i = 1:length(L.W)
				L.W[i] -= learning_rate * sign(L.GW[i])
			end
			# for i = 1:length(L.GB)
			# 	L.UB[i] = 0.9 * L.UB[i] + 0.1 * abs(L.GB[i])
			# 	if L.UB[i] > 0.0
			# 		L.B[i] -= (lr_use / L.UB[i]) * L.GB[i]
			# 	else
			# 		L.B[i] -= lr_use * L.GB[i]
			# 	end
			# end
			# for i = 1:length(L.GW)
			# 	L.UW[i] = 0.9 * L.UW[i] + 0.1 * abs(L.GW[i])
			# 	if L.UW[i] > 0.0
			# 		L.W[i] -= (lr_use / L.UW[i]) * L.GW[i]
			# 	else
			# 		L.W[i] -= lr_use * L.GW[i]
			# 	end
			# end
			# axpy!(length(L.W), -lr, L.GW, 1, L.W, 1)
			# axpy!(length(L.B), -lr, L.GB, 1, L.B, 1)
			if reset_gradient
				fill!(L.GW, 0.0)
				fill!(L.GB, 0.0)
			end
		end
	end
end

function get_learning_rate{T<:FloatingPoint}(DN::StackedNet{T}, max_parameter_change::T)
	gmax = gradient_maxabs(DN)
	max_parameter_change / gmax
end

function train_sgd!{T<:FloatingPoint}(DN::StackedNet{T}, X_training::Matrix{T}, Y_training::Matrix{T}; X_testing::Union(Nothing, Matrix{T})=nothing, Y_testing::Union(Nothing, Matrix{T})=nothing, custom_error::Union(Nothing, Function)=nothing, iterations::Int=1000, iterations_report::Int=100, learning_rate::T=1e-2, minibatch_size::Int=100, minibatch_replace::Bool=true)
	@inbounds begin
		num_patterns::Int64 = size(X_training, 2)
		# Minibatch size cannot be larger than number of patterns when sampling without replacement.
		if (minibatch_size > num_patterns) && !minibatch_replace
			minibatch_size = num_patterns
		end
		# Reserve space for minibatch integers.
		minibatch_ints = zeros(Int, minibatch_size)
		minibatch_domain = 1:num_patterns
		# Create results dataframe if reporting is required.
		if custom_error == nothing
			df = DataFrame(iteration=Int64[], error_training=T[], error_testing=T[])
		else
			df = DataFrame(iteration=Int64[], error_training=T[], error_testing=T[], custom_error_training=T[], custom_error_testing=T[])
		end
		# Reserve space for output predictions if custom error function specified.
		if custom_error != nothing
			YH_training = zeros(T, (DN.layers[end].no, size(X_training, 2)))
			YH_testing = zeros(T, (DN.layers[end].no, size(X_testing, 2)))
		end
		# Perform the required number of iterations of learning.
		gradient_reset!(DN)
		for iteration = 1:iterations
			# Increment the gradient information based on the minibatch.
			sort!(sample!(minibatch_domain, minibatch_ints, replace=minibatch_replace))
			for minibatch_int in minibatch_ints
				gradient_update!(DN, X_training, Y_training, minibatch_int)
			end
			# Update the parameters based on the gradient information.
			parameters_update!(DN, learning_rate, reset_gradient=true)
			# Decide whether to record performance.
			if (iteration % iterations_report) == 0
				# Construct row to add to dataframe.
				if custom_error == nothing
					row = (iteration, error!(DN, X_training, Y_training), error!(DN, X_testing, Y_testing))
				else
					forward!(DN, X_training, YH_training)
					forward!(DN, X_testing, YH_testing)
					row = (iteration, error!(DN, X_training, Y_training), error!(DN, X_testing, Y_testing), custom_error(YH_training, Y_training), custom_error(YH_testing, Y_testing))
				end
				push!(df, row)
				println("Iteration $iteration. Training error $(df[end, :error_training]). Testing error $(df[end, :error_testing]).")
			end
		end
	end
	df
end


text:
import CUDAnative
using CuArrays

# CuArrays.cufunc(::typeof(^)) = CUDAnative.pow

@diffrule CUDAnative.exp(u::Real) u CUDAnative.exp(u) * dy
@diffrule CUDAnative.pow(u::Real, v::Real) u (v * CUDAnative.pow(u, (v-1)) * dy)
@diffrule CUDAnative.pow(u::Real, v::Real) v CUDAnative.log(u) * CUDAnative.pow(u, v) * dy
@diffrule CUDAnative.log(u::Real) u dy / u
@diffrule CUDAnative.sqrt(u::Real) u (0.5f0 * CUDAnative.pow(u, -0.5f0) * dy)

# # const NON_DISPATCHED_OPS = [log, exp, sqrt, ^, ones]
# # const CUDA_NATIVE_OPS = Dict{Function,Function}(op => op for op in NON_DISPATCHED_OPS)
const CUDANATIVE_OPS = Dict{Function,Function}()

CUDANATIVE_OPS[log] = CUDAnative.log
CUDANATIVE_OPS[exp] = CUDAnative.exp
CUDANATIVE_OPS[sqrt] = CUDAnative.sqrt
CUDANATIVE_OPS[^] = CUDAnative.pow
CUDANATIVE_OPS[ones] = CUDAnative.ones

device_function(device::GPU, f::Function) = get(CUDANATIVE_OPS, f, f)

to_device(device::CPU, x::CuArray) = convert(Array, x)

function to_device(device::GPU, x)
    T = typeof(x)
    flds = fieldnames(T)
    if is_cuarray(x)
        return x
    elseif isa(x, AbstractFloat)
        return Float32(x)
    elseif isa(x, Tuple)
        return ((to_device(device, el) for el in x)...,)
    elseif isempty(flds)
        # primitive or array
        return cu(x)
    else
        # struct, recursively convert and construct type from fields
        fld_vals = [to_device(device, getfield(x, fld)) for fld in flds]
        return T(fld_vals...)
    end
end

# function cuarray_compatible_tform(tape::Tape)
#     new_tape = similar(tape)
#     changed = false
#     for op in tape
#         if op isa Call && haskey(CUDANATIVE_OPS, op.fn)
#             changed = true
#             push!(new_tape, copy_with(op, fn=CUDANATIVE_OPS[op.fn]))
#         else
#             push!(new_tape, op)
#         end
#     end
#     return new_tape, changed
# end

# """
# Transform function to CuArrays compatible.
# """
# function cuda_compatible(f, args)
#     cf = CuArrays.cufunc(f)
#     if f === cf
#         return cf
#     else
#         return transform(cuarray_compatible_tform, f, args)
#     end
# end


text:
using Revise
using JLD2
# using ViewerGL
using LarSurf

show = false
taubin = true
mask_labels = ["liver", "portalvein"]
# mask_labels = ["liver"]
# created by surface_extraction_parallel_ircad01.jl
if show
	using ViewerGL
end

for mask_label in mask_labels
	@info "mask_label: $mask_label"
	@load "ircad_$(mask_label).jld2" V FV
# @load "liver01.jld2" V FV

	t = @elapsed FVtri = LarSurf.triangulate_quads(FV)
	@info "triangulate quads time", t

	@JLD2.save "ircad_$(mask_label)_tri.jld2" V FVtri
	if show
		ViewerGL.VIEW([
		    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.9))
		    # ViewerGL.GLGrid(V,FVtri, 1)
		    # ViewerGL.GLGrid(V,EV1,ViewerGL.Point4d(.9,0,.9,0.9))
			ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
		])
	end

	if taubin
		t = @elapsed Vs = LarSurf.Smoothing.smoothing_FV_taubin(V, FV, 0.6, -0.3, 50)
		@info "smoothing taubin time", t
		# t = @elapsed FVtri = LarSurf.triangulate_quads(FV)

		if show
			ViewerGL.VIEW([
			    ViewerGL.GLGrid(Vs, FVtri, ViewerGL.Point4d(1,0,1,0.1))
				ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
			])
		end
	else
		t = @elapsed Vs = LarSurf.Smoothing.smoothing_FV(V, FVtri, 0.6, 3)
		@info "smoothing time", t
	end
	@info "Smoothing numer of Vs: $(size(Vs))"
	@JLD2.save "ircad_$(mask_label)_sm.jld2" Vs FVtri

	# tm1 = @elapsed EV1 = LarSurf.Smoothing.get_EV_quads(FVtri)
	# # tm2 = @elapsed EV2 = LarSurf.Smoothing.get_EV_quads2(FVtri)
	# # @info "time get EV" tm1 tm2
	# @info "calculate EV", t

	if show
		ViewerGL.VIEW([
		    ViewerGL.GLGrid(Vs,FVtri,ViewerGL.Point4d(1,1,1,0.1))
			ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
		])
	end

end
# ViewerGL.VIEW([
#     ViewerGL.GLGrid(Vs,FVtri,ViewerGL.Point4d(1,0,1,0.1))
# 	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
# ])


text:
"""
Function to initialize the container-type (Boostscratch) and to find the first labels to be called
Arguments:
stepno::Int -> Number of Boostingsteps which should be performed
a::Int -> Number of labels for the covariance matrix to start with
x::Int -> Number of labels which should be additionally called
path_out_myscratch::ASCIIString -> Path where the container (Boostscratch) should be stored
path_out_wantedlabels::ASCIIString -> Path where the wanted labels should be stored
"""
function boosting(stepno::Int, a::Int, x::Int, path_out_myscratch::ASCIIString)
	myscratch = Boostscratch(1, Array{Float64,1}(), Array{Float64,1}(), Unibeta(Array{Float64,1}(),Array{Int,1}(),Array{ASCIIString,1}()), Covarmat(Array{Float64,2}(), Array{ASCIIString,1}()), Array{Float64,1}(stepno), Array{ASCIIString, 1}(stepno), Array{ASCIIString,1}(), 0.1, stepno, 1, Array{Float64,1}(), a, x)
	write_boostscratch(myscratch, path_out_myscratch)
end

"""
Function to perform the boosting algorithm
Arguments:
path_boostscratch::ASCIIString -> Path where the container (Boostscratch) is stored to set the parameters needed for running
path_out_wantedlabels::ASCIIString -> Path where the wanted labels should be stored
paths_in_data::Array{ASCIIString, 1} -> Array with paths which lead to the provided covariance matrices or univariate estimators
"""
function boosting(path_boostscratch::ASCIIString, path_out_wantedlabels::ASCIIString, paths_in_data::Array{ASCIIString, 1})
	exist = true
	myscratch = Boostscratch(0, Array{Float64,1}(), Array{Float64,1}(), Unibeta(Array{Float64,1}(),Array{Int,1}(),Array{ASCIIString,1}()), Covarmat(Array{Float64,2}(), Array{ASCIIString,1}()), Array{Float64,1}(), Array{ASCIIString, 1}(), Array{ASCIIString,1}(), 0, 0, 0, Array{Float64,1}(), 0, 0)
	set_boostscratch!(myscratch, path_boostscratch)
	if(length(myscratch.pooledunibeta.unibeta) == 0)
		fetch_pool_unibeta!(myscratch, paths_in_data)
		myscratch.wantedlabels = selectionofcovs(Unibeta(myscratch.actualnom, ones(1), myscratch.pooledunibeta.labels), myscratch.x)
	else
		fetch_pool_covarmat!(myscratch, paths_in_data)
		if(myscratch.actualstepno > 1)
			reboost!(myscratch)
		end
		while(myscratch.actualstepno <= myscratch.stepno && exist == true)
			newselval = boost!(myscratch)
			myscratch.wantedlabels = getselections(myscratch, newselval)
			if(length(myscratch.wantedlabels) > 0)
				exist = false
			end
			if(exist == true)
				myscratch.actualstepno += 1
			end
		end
	end
	write_wantedlabels(myscratch, path_out_wantedlabels)
	write_boostscratch(myscratch, path_boostscratch)
	if(myscratch.actualstepno >= 2)
		println(myscratch.actualstepno," bsteps performed \n recent selected labels: ",myscratch.selections)
	end
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
# Set plot globals
ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "nul"

using Documenter, ControlSystems, Plots, LinearAlgebra, DSP
#import GR # Bug with world age in Plots.jl, see https://github.com/JuliaPlots/Plots.jl/issues/1047
gr()
default(show=false, size=(800,450))

dir = joinpath(dirname(pathof(ControlSystems)), "..")
cd(dir)

# Copied from Documenter/src/Document.jl, modified to remove # hide lines
Markdown = Documenter.Documents.Markdown
function Documenter.Documents.doctest_replace!(block::Markdown.Code)
    startswith(block.language, "jldoctest") || return false
    # suppress output for `#output`-style doctests with `output=false` kwarg
    if occursin(r"^# output$"m, block.code) && occursin(r";.*output\h*=\h*false", block.language)
        input = first(split(block.code, "# output\n", limit = 2))
        block.code = rstrip(input)
    end
    # Remove # hide lines
    block.code = Documenter.Expanders.droplines(block.code)
    # correct the language field
    block.language = occursin(r"^julia> "m, block.code) ? "julia-repl" : "julia"
    return false
end

println("Making docs")
makedocs(modules=[ControlSystems],
    format=Documenter.HTML(),
    sitename="ControlSystems.jl",
    #strict=true,
    pages=[
        "Home" => "index.md",
        "Guide" => Any[
            "Introduction" => "man/introduction.md",
            "Creating Systems" => "man/creating_systems.md",
        ],
        "Examples" => Any[
            "Design" => "examples/example.md",
        ],
        "Functions" => Any[
            "Constructors" => "lib/constructors.md",
            "Analysis" => "lib/analysis.md",
            "Synthesis" => "lib/synthesis.md",
            "Time and Frequency response" => "lib/timefreqresponse.md",
            "Plotting" => "lib/plotting.md",
        ],
    ]
)

deploydocs(repo = "github.com/JuliaControl/ControlSystems.jl.git")


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

import LightGraphs: AbstractEdge
import MetaGraphs: MetaGraph, MetaDiGraph

include("utils/graph_wrapper.jl")

export Topology, Diagram, topologies

"""
"""
struct Topology{I} <: GraphWrapper{I}
    graph::MetaGraph{I}
end

"""
"""
graph(t::Topology) = t.graph

"""
"""
function (::Type{Topology})(mdg::MetaDiGraph)
    Topology(MetaGraph(mdg))
end

"""
"""
function topologies(loops::Integer, external::Integer)
    error("Not Implemented")
end

"""
"""
function topologies(loops::Integer, incoming::Integer, outgoing::Integer)
    topologies(loops, incoming + outgoing)
end

"""
"""
function topologies(loops::Integer, inout::Pair{<: Integer, <: Integer})
    topologies(loops, first(inout), last(inout))
end

"""
"""
struct Diagram{I} <: AbstractGraph{I}
    graph::MetaDiGraph{I}
end

"""
"""
graph(t::Diagram) = t.graph


text:
using ProbNumODE
using ProbNumODE: remake_prob_with_jac
using DifferentialEquations
using LinearAlgebra
using Distributions
using StatsPlots
pyplot()

FILEDIR = dirname(@__FILE__)

# Problem and reference solution
prob = remake_prob_with_jac(lotka_volterra())
appxsol = solve(remake(prob, u0=big.(prob.u0)), abstol=1e-30, reltol=1e-30)

# Probabilistic solve and errors
sol = solve(prob, EKF1(), sigmarule=:fixedMLE, q=5,
            abstol=1e-7, reltol=1e-4, smooth=true)
errors = sol.u .- appxsol.(sol.t)

# Sampling and plotting
sp = ProbNumODE.sample(sol, 10)
sp_errs = sp .- ProbNumODE.stack(sol.u)
p = plot(
    sol.t, ProbNumODE.stack(errors), xlabel="\$t\$", legend=:bottomleft,
    label=["\$(\\hat{y}(t)-y(t))_1\$" "\$(\\hat{y}(t)-y(t))_2\$"],
)
for i in 1:size(sp)[3]
    plot!(p, sol.t, sp_errs[:, :, i], color=[1 2], label="", linewidth=0.5)
end
plot!(sol.t, ProbNumODE.stack(errors), color=:black, label="")
plot!(p, size=(400, 250))
savefig(joinpath(FILEDIR, "sampled_errors_ekf1.png"))

println("Plot saved to $FILEDIR")

max_stars_repo_path: experiments/0 Solution and Sampling/main.jl

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
module LoopVectorization

using VectorizationBase, SIMDPirates, SLEEFPirates, MacroTools, Parameters
using VectorizationBase: REGISTER_SIZE, REGISTER_COUNT, extract_data, num_vector_load_expr,
    mask, masktable, pick_vector_width_val, valmul, valrem, valmuladd, valadd, valsub
using SIMDPirates: VECTOR_SYMBOLS, evadd, evmul, vrange, reduced_add, reduced_prod
using Base.Broadcast: Broadcasted, DefaultArrayStyle
using LinearAlgebra: Adjoint, Transpose
using MacroTools: prewalk, postwalk

export LowDimArray, stridedpointer, vectorizable,
    @avx, *ˡ, ∗,
    vmap, vmap!

include("costs.jl")
include("operations.jl")
include("graphs.jl")
include("broadcast.jl")
include("determinestrategy.jl")
include("lowering.jl")
include("constructors.jl")
include("map.jl")
# include("precompile.jl")
# _precompile_()

end # module


text:
# Autogenerated wrapper script for LLVM_full_assert_jll for i686-w64-mingw32-cxx11
export clang, libclang, libclang_cpp, libllvm, liblto, llc, llvm_config, llvm_mca, mlir, opt

using Zlib_jll
JLLWrappers.@generate_wrapper_header("LLVM_full_assert")
JLLWrappers.@declare_executable_product(clang)
JLLWrappers.@declare_library_product(libclang, "libclang.dll")
JLLWrappers.@declare_library_product(libclang_cpp, "libclang-cpp.dll")
JLLWrappers.@declare_library_product(libllvm, "libLLVM-13jl.dll")
JLLWrappers.@declare_library_product(liblto, "libLTO.dll")
JLLWrappers.@declare_executable_product(llc)
JLLWrappers.@declare_executable_product(llvm_config)
JLLWrappers.@declare_executable_product(llvm_mca)
JLLWrappers.@declare_library_product(mlir, "libMLIR.dll")
JLLWrappers.@declare_executable_product(opt)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_executable_product(
        clang,
        "tools\\clang.exe",
    )

    JLLWrappers.@init_library_product(
        libclang,
        "bin\\libclang.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libclang_cpp,
        "bin\\libclang-cpp.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libllvm,
        "bin\\libLLVM-13jl.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        liblto,
        "bin\\libLTO.dll",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        llc,
        "tools\\llc.exe",
    )

    JLLWrappers.@init_executable_product(
        llvm_config,
        "tools\\llvm-config.exe",
    )

    JLLWrappers.@init_executable_product(
        llvm_mca,
        "tools\\llvm-mca.exe",
    )

    JLLWrappers.@init_library_product(
        mlir,
        "bin\\libMLIR.dll",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        opt,
        "tools\\opt.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SamplingWrapper

using JavaCall

const AbstractFieldODEStateInterpolator = @jimport org.hipparchus.ode.sampling.AbstractFieldODEStateInterpolator
const AbstractODEStateInterpolator = @jimport org.hipparchus.ode.sampling.AbstractODEStateInterpolator
const Class = @jimport java.lang.Class
const Enum = @jimport java.lang.Enum
const FieldODEFixedStepHandler = @jimport org.hipparchus.ode.sampling.FieldODEFixedStepHandler
const FieldODEStateAndDerivative = @jimport org.hipparchus.ode.FieldODEStateAndDerivative
const FieldODEStateInterpolator = @jimport org.hipparchus.ode.sampling.FieldODEStateInterpolator
const FieldODEStepHandler = @jimport org.hipparchus.ode.sampling.FieldODEStepHandler
const FieldStepNormalizer = @jimport org.hipparchus.ode.sampling.FieldStepNormalizer
const JString = @jimport java.lang.String
const ODEFixedStepHandler = @jimport org.hipparchus.ode.sampling.ODEFixedStepHandler
const ODEStateAndDerivative = @jimport org.hipparchus.ode.ODEStateAndDerivative
const ODEStateInterpolator = @jimport org.hipparchus.ode.sampling.ODEStateInterpolator
const ODEStepHandler = @jimport org.hipparchus.ode.sampling.ODEStepHandler
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const StepNormalizer = @jimport org.hipparchus.ode.sampling.StepNormalizer
const StepNormalizerBounds = @jimport org.hipparchus.ode.sampling.StepNormalizerBounds
const StepNormalizerMode = @jimport org.hipparchus.ode.sampling.StepNormalizerMode

export AbstractFieldODEStateInterpolator
export AbstractODEStateInterpolator
export FieldODEFixedStepHandler
export FieldODEStateInterpolator
export FieldODEStepHandler
export FieldStepNormalizer
export ODEFixedStepHandler
export ODEStateInterpolator
export ODEStepHandler
export StepNormalizer
export StepNormalizerBounds
export StepNormalizerMode
export compare_to
export describe_constable
export equals
export first_included
export get_class
export get_current_state
export get_declaring_class
export get_global_current_state
export get_global_previous_state
export get_interpolated_state
export get_previous_state
export handle_step
export hash_code
export init
export is_current_state_interpolated
export is_forward
export is_previous_state_interpolated
export last_included
export name
export notify
export notify_all
export ordinal
export restrict_step
export to_string
export value_of
export values
export wait

include("abstract_field_ode_state_interpolator.jl")
include("abstract_ode_state_interpolator.jl")
include("field_ode_fixed_step_handler.jl")
include("field_ode_state_interpolator.jl")
include("field_ode_step_handler.jl")
include("field_step_normalizer.jl")
include("ode_fixed_step_handler.jl")
include("ode_state_interpolator.jl")
include("ode_step_handler.jl")
include("step_normalizer.jl")
include("step_normalizer_bounds.jl")
include("step_normalizer_mode.jl")

end


text:
"""
    fillwt(itype)

Fill weight tables for all order and number of flavours.
`itype` is used to select un-polarised pdfs (1), polarised
pdfs (2) or fragmentation functions (3).

# Returns
- `nwds::Integer`: number of words used in memory.
"""
function fillwt(itype::Integer)

    itype = Ref{Int32}(itype)
    nwds = Ref{Int32}()

    # For now-obsolete integer id vars
    idnum1 = Ref{Int32}()
    idnum2 = Ref{Int32}()

    @ccall fillwt_(itype::Ref{Int32}, idnum1::Ref{Int32},
                   idnum2::Ref{Int32}, nwds::Ref{Int32})::Nothing

    nwds[]
end

"""
    dmpwgt(itype, lun, filename)

Dump weight tables of a given `itype` to `filename`.
"""
function dmpwgt(itype::Integer, lun::Integer, filename::String)

    itype = Ref{Int32}(itype)
    lun = Ref{Int32}(lun)

    @ccall dmpwgt_(itype::Ref{Int32}, lun::Ref{Int32}, filename::Cstring)::Nothing

    nothing
end

"""
    readwt(lun, filename)

Read weight tables from filename.
TODO: Fix implementation which currently crashes.
"""
function readwt(lun::Integer, filename::String)

    lun = Ref{Int32}(lun)
    ierr = Ref{Int32}()
    nwds = Ref{Int32}()

    # For now-obsolete integer id vars
    idnum1 = Ref{Int32}()
    idnum2 = Ref{Int32}()

    @ccall readwt_(lun::Ref{Int32}, filename::Cstring, idnum1::Ref{Int32},
                   idnum2::Ref{Int32}, nwds::Ref{Int32}, ierr::Ref{Int32})::Nothing

    nwds[], ierr[]

end

"""
    wtfile(itype, filename)

Maintains an up-to-date weight table in filename.
"""
function wtfile(itype::Integer, filename::String)

    itype = Ref{Int32}(itype)

    @ccall wtfile_(itype::Ref{Int32}, filename::Cstring)::Nothing

    nothing
end


text:
using HDF5, H5Sparse, MAT

function writeHDF5(filename::String, data::Tuple{String, AbstractMatrix}...)
  fid = h5open(filename, "w")
  for mat in data
    m_write_hdf5(filename, mat[1], mat[2])
  end
  close(fid)
end

function writeHDF5(filename::String; kwargs...)
  fid = h5open(filename, "w")
  for (k, v) in pairs(kwargs)
    m_write_hdf5(filename, String(k), v)
  end
  close(fid)
end

m_write_hdf5(f::String, varname::String, var::Matrix) = h5write(f, varname, var)
m_write_hdf5(f::String, varname::String, var::SparseMatrixCSC) = H5SparseMatrixCSC(f, varname, var)

function loadHDF5(filename::String)
  fid = h5open(filename)
  d = Dict()
  for key in keys(fid)
    d[key] =  m_read(fid, key, fid[key])
  end
  return d
end

function m_read(::Any, ::Any, dset::HDF5.Dataset)
  return read(dset)
end
function m_read(fid, key, dset::HDF5.Group)
  return H5SparseMatrixCSC(fid, key) |> SparseMatrixCSC
end

function writeMAT(filename::String, data::Tuple{String, AbstractMatrix}...)
  d = Dict()
  for mat in data
    d[mat[1]] = mat[2]
  end
  matwrite(filename, d)
end

function writeMAT(filename::String; kwargs...)
  d = Dict()
  for (k, v) in pairs(kwargs)
    d[String(k)] = v
  end
  matwrite(filename, d)
end

function loadMAT(filename::String)
  return matread(filename)
end

export writeHDF5, writeMAT


text:
function bob(stimulus)
    # 1) Checar si la cadena esta vacia
    # 2) Checar si la cadena tiene un signo de interrogacion (?)
    # 3) Checar si la cadena es toda en mayusculas
    # 4) Cualquier otra cosa regresa "Whatever."
    nospace = strip(stimulus)
    if isempty(nospace)
        return "Fine. Be that way!"
    elseif uppercase(stimulus) == stimulus
        if endswith(stimulus, "?")
            return "Calm down, I know what I'm doing!"
        else
            return "Whoa, chill out!"
        end
    else
        if endswith(stimulus, "?")
            return "Sure."
        end
    end
    return "Whatever."
end

max_stars_repo_name: mauricioriva/learn-julia

max_issues_repo_name: mauricioriva/learn-julia

max_forks_repo_name: mauricioriva/learn-julia


text:
###############################################################################
#
#  Localization of K[1/x] at (1/x), i.e. k_\infty(x) \subseteq k(x)
#
#
#  (C) 2021 William Hart
#
###############################################################################

export KInftyRing, KInftyElem, function_field

###############################################################################
#
#  Declaration types
#  KInftyRing / KInftyElem
#
###############################################################################

mutable struct KInftyRing{T <: FieldElement} <: Hecke.Ring
  K::Generic.RationalFunctionField{T}

  function KInftyRing{T}(K::Generic.RationalFunctionField{T}, cached::Bool) where T <: FieldElement
    return AbstractAlgebra.get_cached!(KInftyID, K, cached) do
      new{T}(K)
    end::KInftyRing{T}
  end
end

const KInftyID = Dict{Generic.RationalFunctionField, Hecke.Ring}()

mutable struct KInftyElem{T <: FieldElement} <: Hecke.RingElem
  d::Generic.Rat{T}
  parent::KInftyRing{T}
end

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

elem_type(::Type{KInftyRing{T}}) where T <: FieldElement = KInftyElem{T}

parent_type(::Type{KInftyElem{T}}) where T <: FieldElement = KInftyRing{T}

is_domain_type(::Type{KInftyElem{T}}) where {T} = true

# return the rational function field which KInfty wraps, mostly internal use
function function_field(R::KInftyRing{T}) where T <: FieldElement
  return R.K::Generic.RationalFunctionField{T}
end

parent(a::KInftyElem{T}) where T <: FieldElement = a.parent

function check_parent(a::KInftyElem{T}, b::KInftyElem{T})  where T <: FieldElement
  parent(a) != parent(b) && error("Parent objects do not match")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

data(a::KInftyElem{T}) where T <: FieldElement = a.d::Generic.Rat{T}

function numerator(a::KInftyElem{T}, canonicalise::Bool=true) where T <: FieldElement
  return numerator(data(a), canonicalise)
end

function denominator(a::KInftyElem{T}, canonicalise::Bool=true) where T <: FieldElement
  return denominator(data(a), canonicalise)
end

gen(R::KInftyRing) = R(inv(gen(R.K)))

characteristic(R::KInftyRing) = characteristic(R.K)

@doc Markdown.doc"""
     degree(a::KInftyElem)

Return the degree of the given element, i.e.
`degree(numerator) - degree(denominator)`.
"""
degree(a::KInftyElem) = degree(numerator(a, false)) - degree(denominator(a, false))

@doc Markdown.doc"""
    valuation(a::KInftyElem)

Return the degree valuation of the given element, i.e. `-degree(a)`.
"""
valuation(a::KInftyElem) = -degree(a)

zero(K::KInftyRing{T}) where T <: FieldElement = K(0)

one(K::KInftyRing{T}) where T <: FieldElement = K(1)

iszero(a::KInftyElem{T}) where T <: FieldElement = iszero(data(a))

isone(a::KInftyElem{T}) where T <: FieldElement = isone(data(a))

function is_unit(a::KInftyElem{T}) where T <: FieldElement
  return degree(numerator(data(a), false)) ==
                                            degree(denominator(data(a), false))
end

@doc Markdown.doc"""
    in(a::Generic.Rat{T}, R::KInftyRing{T}) where T <: FieldElement

Return `true` if the given element of the rational function field is an
element of `k_\infty(x)`, i.e. if `degree(numerator) <= degree(denominator)`.
"""
function in(a::Generic.Rat{T}, R::KInftyRing{T}) where T <: FieldElement
  if parent(a) != function_field(R)
    return false
  end
  return degree(numerator(a, false)) <= degree(denominator(a, false))
end

function Base.deepcopy_internal(a::KInftyElem{T}, dict::IdDict) where T <: FieldElement
  c = Base.deepcopy_internal(data(a), dict)
  parent(a)(Base.deepcopy_internal(data(a), dict))
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function AbstractAlgebra.expressify(a::KInftyElem; context = nothing)
  return AbstractAlgebra.expressify(data(a), context = context)
end

function show(io::IO, a::KInftyElem)
  print(io, AbstractAlgebra.obj_to_string(a, context = io))
end

function show(io::IO, R::KInftyRing)
  print(io, "Degree localization of ", function_field(R))
end

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(a::KInftyElem{T}) where T <: FieldElement
  parent(a)(-data(a), false)
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::KInftyElem{T}, b::KInftyElem{T})  where T <: FieldElement
  check_parent(a, b)
  return parent(a)(data(a) + data(b), false)
end

function -(a::KInftyElem{T}, b::KInftyElem{T})  where T <: FieldElement
  check_parent(a, b)
  return parent(a)(data(a) - data(b), false)
end

function *(a::KInftyElem{T}, b::KInftyElem{T})  where T <: FieldElement
  check_parent(a, b)
  return parent(a)(data(a)*data(b), false)
end

###############################################################################
#
#   Inplace
#
###############################################################################

function mul!(a::KInftyElem{T}, b::KInftyElem{T}, c::KInftyElem{T}) where {T}
  return b*c
end

function addeq!(a::KInftyElem{T}, b::KInftyElem{T}) where {T}
  return a+b
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::KInftyElem{T}, b::KInftyElem{T}) where T <: FieldElement
  check_parent(a, b)
  return data(a) == data(b)
end

###############################################################################
#
#  Inversion
#
###############################################################################

@doc Markdown.doc"""
     inv(a::KInftyElem{T}, checked::Bool = true)  where T <: FieldElement
Returns the inverse element of $a$ if $a$ is a unit.
If 'checked = false' the invertibility of $a$ is not checked and the
corresponding inverse element of the rational function field is returned.
"""
function inv(a::KInftyElem{T}, checked::Bool = true)  where T <: FieldElement
  b = inv(data(a))
  return parent(a)(b, checked)
end

###############################################################################
#
#  Exact division
#
###############################################################################

@doc Markdown.doc"""
     divides(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool = true) where T <: FieldElement

Returns tuple `(flag, c)` where `flag = true` if $b$ divides $a$ and $a = bc$,
otherwise `flag = false` and $c = 0$.
If `checked = false` the corresponding element of the rational function field
is returned and it is not checked whether it is an element of the given
localization.
"""
function divides(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool = true) where T <: FieldElement
  check_parent(a, b)
  if iszero(a)
    return true, a
  end
  if checked && degree(a) > degree(b)
    return false, parent(a)()
  else
    elem = divexact(data(a), data(b))
    return true, parent(a)(elem, false)
  end
end

@doc Markdown.doc"""
     divexact(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool = true)  where {T <: nf_elem}
Returns element 'c' of given localization such that $a = bc$ if such element
exists. If `checked = false` the corresponding element of the rational function
field is returned and it is not checked whether it is an element of the given
localization.
"""
function divexact(a::KInftyElem{T}, b::KInftyElem{T}; check::Bool = true)  where T <: FieldElement
  iszero(b) && throw(DivideError())
  d = divides(a, b, check)
  d[1] ? d[2] : error("$a not divisible by $b in the given localization")
end

###############################################################################
#
#  Euclidean division
#
###############################################################################

#= an attempt to have a canonical rep modulo the ideal
 generated by b
 - <b> = <1/t^n>, the canonical rep
 - translate via t-> 1/t into the "normal" world
 - a rep of a/b mod t^n is ((a % t^n)*invmod(b, t^n) % t^n)
 - translate back
 The "canonical" rep should be r/t^l for deg(r) < n and l = deg(r)
 In particular, for n=1, this is a scalar.

 The translation is done using reverse:
  a-> a(1/t)*t^deg(a) = reverse(a)
  so a/b -> a(1/t)/b(1/t) and this needs to be expanded by t^max/t^max to
  make both num and den polynomials
  a/b -> reverse(a, max)//reverse(b, max)
=#

function mod(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool=true) where T <: FieldElement
  check_parent(a, b)
  iszero(b) && throw(DivideError())
  iszero(a) && return deepcopy(a)

  if degree(a) > degree(b)
    n = -degree(b)
    na = numerator(a.d)
    da = denominator(a.d)
    mx = max(degree(na), degree(da))+1
    if mx != 0
      na = reverse(na, mx)
      da = reverse(da, mx)
    end
    t = gen(parent(na))
    tn = t^n
    r = mod(mod(na, tn) * invmod(da, tn), tn)
    Qt = parent(a.d)
    return parent(a)(Qt(reverse(r))//Qt(t)^degree(r))
  else
    return zero(parent(a))
  end
end

function div(a::KInftyElem{T}, b::KInftyElem{T}, check::Bool=true) where T <: FieldElement
  check_parent(a, b)
  iszero(b) && throw(DivideError())
  iszero(a) && return a
  if degree(a) > degree(b)
    return divexact(a-mod(a, b), b, check = check)
  else
    return divexact(a, b, check = check)
  end
end

function divrem(a::KInftyElem{T}, b::KInftyElem{T}, check::Bool=true) where T <: FieldElement
  check_parent(a, b)
  iszero(b) && throw(DivideError())
  iszero(a) && return a, a

  if degree(a) > degree(b)
    r = mod(a, b)
    return divexact(a-r, b, check = check), r
  else
    return divexact(a, b, check = check), parent(a)()
  end
end

###############################################################################
#
#  GCD
#
###############################################################################

function gcd(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool=true) where T <: FieldElement
  check_parent(a, b)
  t = gen(parent(a))
  iszero(a) && iszero(b) && return a
  iszero(a) && return t^-degree(b)
  iszero(b) && return t^-degree(a)

  return t^-max(degree(a), degree(b))
end

function lcm(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool=true) where T <: FieldElement
  t = gen(parent(a))
  iszero(a) && iszero(b) && return a
  iszero(a) && return t^-degree(b)
  iszero(b) && return t^-degree(a)

  return t^-min(degree(a), degree(b))
end
###############################################################################
#
#  GCDX
#
###############################################################################

function gcdx(a::KInftyElem{T}, b::KInftyElem{T}, checked::Bool=true) where T <: FieldElement
  check_parent(a, b)
  K = parent(a)
  t = gen(K)

  if iszero(a)
    iszero(b) && a, a, one(K)
    g = t^-degree(b)
    return g, zero(K), inv(canonical_unit(b))
  end

  if iszero(b) || degree(a) >= degree(b)
    g = t^-degree(a)
    return g, inv(canonical_unit(a)), zero(K)
  end

  g = t^-degree(b)
  return g, zero(K), inv(canonical_unit(b))
end

function gcdinv(a::KInftyElem{T}, b::KInftyElem{T}) where {T}
  g, q, w = gcdx(a, b)
  @assert is_unit(g)
  return one(parent(a)), q*inv(g)
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::KInftyElem{T}, b::Int) where T <: FieldElement
  return parent(a)(data(a)^b, false)
end

###############################################################################
#
#   Random functions
#
###############################################################################

RandomExtensions.maketype(R::KInftyRing, _) = elem_type(R)

function RandomExtensions.make(S::KInftyRing, vs...)
  R = function_field(S)
  if length(vs) == 1 && elem_type(R) == Random.gentype(vs[1])
    RandomExtensions.Make(S, vs[1]) # forward to default Make constructor
  else
    make(S, make(R, vs...))
  end
end

function rand(rng::AbstractRNG,
              sp::SamplerTrivial{<:Make2{<:RingElement, <:KInftyRing}})
  S, v = sp[][1:end]
  R = function_field(S)
  d = rand(rng, v)
  if iszero(d)
    return S(d, false)
  end
  if degree(numerator(d, false)) <= degree(denominator(d, false))
    return S(d, false)
  else
    return S(inv(d), false)
  end
end

rand(rng::AbstractRNG, S::KInftyRing, v...) =
  rand(rng, make(S, v...))

rand(S::KInftyRing, v...) = rand(GLOBAL_RNG, S, v...)

###############################################################################
#
#   Promotion rules
#
###############################################################################

AbstractAlgebra.promote_rule(::Type{KInftyElem{T}}, ::Type{KInftyElem{T}}) where T <: FieldElement = KInftyElem{T}

AbstractAlgebra.promote_rule(::Type{KInftyElem{T}}, ::Type{Generic.Rat{T}}) where T <: FieldElement = KInftyElem{T}

function AbstractAlgebra.promote_rule(::Type{KInftyElem{T}}, ::Type{U}) where {T <: FieldElement, U <: RingElem}
  promote_rule(T, U) == T ? KInftyElem{T} : Union{}
end

###############################################################################
#
#  Parent call overloading
#
###############################################################################

(R::KInftyRing)() = R(function_field(R)())

function (R::KInftyRing{T})(a::Generic.Rat{T}, checked::Bool=true) where T <: FieldElement
  checked && degree(numerator(a, false)) > degree(denominator(a, false)) &&
                                           error("Not an element of k_infty(x)")
  return KInftyElem{T}(a, R)
end

(R::KInftyRing)(a::RingElement) = R(function_field(R)(a), false)

function (R::KInftyRing{T})(a::KInftyElem{T}) where T <: FieldElement
  parent(a) != R && error("Cannot coerce element")
  return a
end

###############################################################################
#
#  Canonical unit
#
###############################################################################

function canonical_unit(a::KInftyElem)
  iszero(a) && return one(parent(a))
  num = numerator(a, false)
  den = denominator(a, false)
  R = parent(num)
  x = gen(R)
  return parent(a)((x^(degree(den) - degree(num))*num)//den)
end

###############################################################################
#
#  Factorisation
#
###############################################################################

function factor(d::KInftyElem)
  t = gen(parent(d))
  a = degree(d)
  return Fac(canonical_unit(d), Dict(t=>-a))
end

###############################################################################
#
#  Constructors
#
###############################################################################

function ResidueField(K::KInftyRing{T}, a::KInftyElem{T}) where {T <: FieldElement}
  F = base_ring(K.K)
  @assert degree(a) == -1
  #TODO: can be optimized, see blurb of euc. div. above
  return F, MapFromFunc(x -> leading_coefficient(numerator(mod(x, a))), y-> K(y), K, F)
end
#TODO: ResidueRing is probably "just" poly of deg < n, think about it

@doc Markdown.doc"""
    Localization(K::RationalFunctionField{T}, ::typeof(degree)) where T <: FieldElement

Return the localization of $k[1/x]$ at $(1/x)$ inside the rational function
field $k(x)$, i.e. the localization of the function field at the point at
infinity, i.e. the valuation ring for valuation $-$degree$(x)$. This is the ring
$k_\infty(x) = \{ f/g | \deg(f) \leq \deg(g)\}$.
"""
function Localization(K::Generic.RationalFunctionField{T}, ::typeof(degree); cached::Bool=true) where T <: FieldElement
  return KInftyRing{T}(K, cached)
end


text:
"""Divides a list of indices 1:N in k groups"""
function kfoldperm(N,k; type="random", breaks=[])
	n,r = divrem(N,k)
	b = collect(1:n:N+1)
	Nb = length(b)
	for i in 1:Nb
		b[i] += i > r ? r : i-1
	end
	b = sort(unique(append!(b, breaks)))
	Nbreaks = length(b) - Nb
	if type=="random"
		p = randperm(N)
	elseif type=="ordered"
		p = 1:N
	end
	[p[r] for r in [b[i]:b[i+1]-1 for i=1:k+Nbreaks]] #TODO: use RF starts and ends differently to remove PATCH in run_sim_time_iter
end

"""Calcualtes key time points that can not be missed during the simulation"""
function points_from_key_times(times ; dt)
	t = Float64[]
	for i = 1:length(times)-1
		if dt < times[i+1] - times[i]
			taux = collect(range(times[i],times[i+1];step=dt))
		else
			taux = [times[i], times[i+1]]
		end
		append!(t, taux)
	end
	t
end
"""
Uniform time-step calculation
"""
function get_uniform_times(seq,Δt;Δt_rf=1e-4)
	t, Δt = get_variable_times(seq; dt=Δt, dt_rf=Δt_rf)
end

"""
Variable time-step calculation
"""
function get_variable_times(seq; dt=1, dt_rf=1e-4)
	t = Float64[]
	ΔT = durs(seq) #Duration of sequence block
	T0 = cumsum([0; ΔT[:]]) #Start time of each block
	for i = 1:length(seq)
		s = seq[i] #Current sequence block
		t0 = T0[i]
		if is_ADC_on(s)
			ts = get_sample_times(s) .+ t0
			taux = points_from_key_times(ts; dt) # ADC sampling
			append!(t, taux)
		end
		if is_RF_on(s)
			y = s.RF[1]
			delay, T = y.delay, y.T
			t1 = t0 + delay
			t2 = t1 + sum(T)
			rf0 = t0 + get_RF_center(y) #get_RF_center includes delays
			taux = points_from_key_times([t1,rf0,t2]; dt=dt_rf) # Arbitrary RF
			append!(t, taux)
		end
		if is_GR_on(s)
			active_gradients = []
			if is_Gx_on(s) append!(active_gradients, s.GR.x) end
			if is_Gy_on(s) append!(active_gradients, s.GR.y) end
			if is_Gz_on(s) append!(active_gradients, s.GR.z) end
			for y = active_gradients
				ts = get_theo_t(y) .+ t0
				taux = points_from_key_times(ts; dt)
				append!(t, taux)
			end
		end
	end
	t = sort(unique(x -> round(x*1e8), t)) #Removing repeated points
	Δt = t[2:end] .- t[1:end-1]
	t = t[1:end-1]
	t, Δt
end

"""Calculate RF key points (start-end) to split simulation in RF and non-RF parts. """
function get_breaks_in_RF_key_points(seq::Sequence, t)
	T0 = cumsum([0; durs(seq)[:]])
	# Identifying RF key points
	key_points = Float64[]
	key_idxs = Int[]
	for (i, s) = enumerate(seq)
		if is_RF_on(s)
			t0 = T0[i] + s.RF.delay[1]	#start of RF waverform
			tf = T0[i] + s.RF.dur[1]	#end of RF waveform
			append!(key_points, [t0; tf])
			idx0 = argmin(abs.(t.-t0))
			idxf = argmin(abs.(t.-tf))
			append!(key_idxs,   [idx0; idxf])
		end
	end
	key_idxs
end


text:
const textDocument_codeAction_request_type = JSONRPC.RequestType("textDocument/codeAction", CodeActionParams, Union{Vector{Command}, Vector{CodeAction}, Nothing})
const textDocument_completion_request_type = JSONRPC.RequestType("textDocument/completion", CompletionParams, Union{CompletionList, Vector{CompletionItem}, Nothing})
const textDocument_signatureHelp_request_type = JSONRPC.RequestType("textDocument/signatureHelp", TextDocumentPositionParams, Union{SignatureHelp, Nothing})
const textDocument_definition_request_type = JSONRPC.RequestType("textDocument/definition", TextDocumentPositionParams, Union{Location, Vector{Location}, Vector{LocationLink}, Nothing})
const textDocument_formatting_request_type = JSONRPC.RequestType("textDocument/formatting", DocumentFormattingParams, Union{Vector{TextEdit}, Nothing})
const textDocument_range_formatting_request_type = JSONRPC.RequestType("textDocument/rangeFormatting", DocumentRangeFormattingParams, Union{Vector{TextEdit}, Nothing})
const textDocument_references_request_type = JSONRPC.RequestType("textDocument/references", ReferenceParams, Union{Vector{Location}, Nothing})
const textDocument_rename_request_type = JSONRPC.RequestType("textDocument/rename", RenameParams, Union{WorkspaceEdit, Nothing})
const textDocument_prepareRename_request_type = JSONRPC.RequestType("textDocument/prepareRename", PrepareRenameParams, Union{Range,Nothing})
const textDocument_documentSymbol_request_type = JSONRPC.RequestType("textDocument/documentSymbol", DocumentSymbolParams, Union{Vector{DocumentSymbol}, Vector{SymbolInformation}, Nothing})
const textDocument_documentHighlight_request_type = JSONRPC.RequestType("textDocument/documentHighlight", DocumentHighlightParams, Union{Vector{DocumentHighlight}, Nothing})
const textDocument_hover_request_type = JSONRPC.RequestType("textDocument/hover", TextDocumentPositionParams, Union{Hover, Nothing})
const textDocument_didOpen_notification_type = JSONRPC.NotificationType("textDocument/didOpen", DidOpenTextDocumentParams)
const textDocument_didClose_notification_type = JSONRPC.NotificationType("textDocument/didClose", DidCloseTextDocumentParams)
const textDocument_didChange_notification_type = JSONRPC.NotificationType("textDocument/didChange", DidChangeTextDocumentParams)
const textDocument_didSave_notification_type = JSONRPC.NotificationType("textDocument/didSave", DidSaveTextDocumentParams)
const textDocument_willSave_notification_type = JSONRPC.NotificationType("textDocument/willSave", WillSaveTextDocumentParams)
const textDocument_willSaveWaitUntil_request_type = JSONRPC.RequestType("textDocument/willSaveWaitUntil", WillSaveTextDocumentParams, Union{Vector{TextEdit}, Nothing})
const textDocument_publishDiagnostics_notification_type = JSONRPC.NotificationType("textDocument/publishDiagnostics", PublishDiagnosticsParams)
const textDocument_selectionRange_request_type = JSONRPC.RequestType("textDocument/selectionRange", SelectionRangeParams, Union{Vector{SelectionRange}, Nothing})
const textDocument_documentLink_request_type = JSONRPC.RequestType("textDocument/documentLink", DocumentLinkParams, Union{Vector{DocumentLink}, Nothing})

const workspace_executeCommand_request_type = JSONRPC.RequestType("workspace/executeCommand", ExecuteCommandParams, Any)
const workspace_symbol_request_type = JSONRPC.RequestType("workspace/symbol", WorkspaceSymbolParams, Union{Vector{SymbolInformation}, Nothing})
const workspace_didChangeWatchedFiles_notification_type = JSONRPC.NotificationType("workspace/didChangeWatchedFiles", DidChangeWatchedFilesParams)
const workspace_didChangeConfiguration_notification_type = JSONRPC.NotificationType("workspace/didChangeConfiguration", DidChangeConfigurationParams)
const workspace_didChangeWorkspaceFolders_notification_type = JSONRPC.NotificationType("workspace/didChangeWorkspaceFolders", DidChangeWorkspaceFoldersParams)
const workspace_applyEdit_request_type = JSONRPC.RequestType("workspace/applyEdit", ApplyWorkspaceEditParams, ApplyWorkspaceEditResponse)
const workspace_configuration_request_type = JSONRPC.RequestType("workspace/configuration", ConfigurationParams, Vector{Any})
const julia_activateenvironment_notification_type = JSONRPC.NotificationType("julia/activateenvironment", NamedTuple{(:envPath,),Tuple{String}})
const julia_refreshLanguageServer_notification_type = JSONRPC.NotificationType("julia/refreshLanguageServer", Nothing)

const initialize_request_type = JSONRPC.RequestType("initialize", InitializeParams, InitializeResult)
const initialized_notification_type = JSONRPC.NotificationType("initialized", InitializedParams)
const shutdown_request_type = JSONRPC.RequestType("shutdown", Nothing, Nothing)
const exit_notification_type = JSONRPC.NotificationType("exit", Nothing)
const client_registerCapability_request_type = JSONRPC.RequestType("client/registerCapability", RegistrationParams, Nothing)

const cancel_notification_type = JSONRPC.NotificationType("\$/cancelRequest", CancelParams)
const setTrace_notification_type = JSONRPC.NotificationType("\$/setTrace", SetTraceParams)
# TODO This seems to not exist in the spec?
const setTraceNotification_notification_type = JSONRPC.NotificationType("\$/setTraceNotification", Nothing)

const window_workDoneProgress_create_request_type = JSONRPC.RequestType("window/workDoneProgress/create", WorkDoneProgressCreateParams, Nothing)
const window_showMessage_notification_type = JSONRPC.NotificationType("window/showMessage", ShowMessageParams)
const window_showMessage_request_type = JSONRPC.RequestType("window/showMessageRequest", ShowMessageRequestParams, Union{MessageActionItem, Nothing})

const progress_notification_type = JSONRPC.NotificationType("\$/progress", ProgressParams)

const telemetry_event_notification_type = JSONRPC.NotificationType("telemetry/event", Any)


text:
using OVERT

a_lead = -2.0
mu = 0.0001

function acc_dynamics(x::Array{T, 1} where {T <: Real},
	                  u::Array{T, 1} where {T <: Real})
	# x1,x2,x3 are lead vehicle variables
    dx1 = x[2] # dposition
    dx2 = x[3] # dvelocity
    dx3 = -2 * x[3] + 2 * a_lead - mu * x[2]^2 # daccel
	# x4, x5, x6 are ego vehicle variables
	# which is opposite from the notation used in the paper...
	dx4 = x[5]
	dx5 = x[6]
	dx6 = -2 * x[6] + 2 * u[1] - mu * x[5]^2
    return [dx1, dx2, dx3, dx4, dx5, dx6]
end

acc_ẋ₃ = :(-2 * x3 + 2 * $a_lead - $mu * x2^2)
acc_ẋ₆ = :(-2 * x6 + 2 * u1 - $mu * x5^2)

function acc_dynamics_overt(range_dict::Dict{Symbol, Array{T, 1}} where {T <: Real},
	                              N_OVERT::Int,
							      t_idx::Union{Int, Nothing}=nothing)
	if isnothing(t_idx)
		v1 = acc_ẋ₃
		v1_oA = overapprox(v1, range_dict; N=N_OVERT, ϵ=0.5)

		v2 = acc_ẋ₆
		v2_oA = overapprox(v2, range_dict; N=N_OVERT, ϵ=0.5)
	else
		#v1 = "x2_$t_idx^2" # debug idea: try with simpler dynamics
		v1 = "-2 * x3_$t_idx + 2 * $a_lead - $mu * x2_$t_idx^2"
    	v1 = Meta.parse(v1)
		v1_oA = overapprox(v1, range_dict; N=N_OVERT, ϵ=0.5)

		v2 = "-2 * x6_$t_idx + 2 * u1_$t_idx - $mu * x5_$t_idx^2"
		v2 = Meta.parse(v2)
		v2_oA = overapprox(v2, range_dict; N=N_OVERT, ϵ=0.5)
	end
    oA_out = add_overapproximate([v1_oA, v2_oA])
    return oA_out, [v1_oA.output, v2_oA.output]
end

function acc_update_rule(input_vars, control_vars, overt_output_vars)
    integration_map = Dict(input_vars[1] => input_vars[2],
                           input_vars[2] => input_vars[3],
                           input_vars[3] => overt_output_vars[1],
                           input_vars[4] => input_vars[5],
                           input_vars[5] => input_vars[6],
                           input_vars[6] => overt_output_vars[2])
    return integration_map
end

acc_input_vars = [:x1, :x2, :x3, :x4, :x5, :x6]
acc_control_vars = [:u1]

T_gap = 1.4
acc_measurement_model = [[1, 0, 0, -1, -T_gap, 0]::Array{Float64}] # relative distance  - v_ego*Tgap
# but a measurement model can contain multiple rows of measurements!

ACC = OvertProblem(
	acc_dynamics,
	acc_dynamics_overt,
	acc_update_rule,
	acc_input_vars,
	acc_control_vars,
 	acc_measurement_model
 )

"""
Note: Compared this dynamics variable ordering in combination with controller acc_controller.nnet with controller outputs in
https://github.com/souradeep-111/sherlock/blob/master/systems_with_networks/ARCH_2019/ACC/simulation_script.m
with their different variable ordering (xxvvaa) and controller outputs match to 2 decimal places in most cases indicating successful
handling of the differing variable ordering.
Closed loop simulations will not be the same as that above script uses ode45, but it matches closely for the first few steps.
"""


text:
macro StenoGraph(ex)
    if ex.head == :block "Please use a multiline block, e.g. @StenoGraph begin ... end"
        exs = filter(x -> !isa(x, LineNumberNode), ex.args)
        exs = StenoGraphs.quote_symbols.(exs)
        exs = StenoGraphs.addition_to_vector.(exs)
        vec = Expr(:call, :vcat, exs...)
        return esc(:(vcat(unarrow.($vec)...)))
    else
        exs = StenoGraphs.quote_symbols(ex)
        exs = StenoGraphs.addition_to_vector(ex)
        return esc(:(unarrow($ex)))
    end
end


text:
@testset "Testing BF16" begin
    # Note: use `abs` to serve as an example of an activation function that cannot be fused
    # with the `innerproduct` kernel.
    activations = [identity, Flux.relu, Flux.sigmoid, abs]
    BFloat16 = OneDNN.BFloat16
    for σ in activations
        @show σ
        m = Flux.Dense(128, 256, σ)
        M32 = OneDNN.Dense(m)
        M16 = OneDNN.Dense(
            convert.(BFloat16, transpose(m.weight)),
            convert.(BFloat16, m.bias),
            σ,
        )

        x32 = randn(Float32, 128, 128)
        y = m(x32)

        X32 = OneDNN.Memory(x32)
        X16 = OneDNN.Memory(convert.(BFloat16, x32))

        @test M32.optimized_weights == false
        @test M16.optimized_weights == false

        Y32 = M32(X32)
        Y16 = M16(X16)
        # Did the automatic weight transformation take place?
        @test M32.optimized_weights
        @test M16.optimized_weights

        @test isapprox(y, OneDNN.materialize(Y32))
        @test isapprox(y, OneDNN.materialize(Y16))

        # Ensure inference works correctly
        @inferred M32(X32)
        @inferred M16(X16)

        # Now - try backprop
        y, back_jl = Zygote._pullback(m, x32)
        Y32, back_dnnl32 = Zygote._pullback(M32, X32)
        Y16, back_dnnl16 = Zygote._pullback(M16, X16)

        @test isa(Y32, OneDNN.Memory)
        @test isa(Y16, OneDNN.Memory)

        @test isapprox(y, OneDNN.materialize(Y32))
        @test isapprox(y, OneDNN.materialize(Y16))

        @inferred Zygote._pullback(M32, X32)
        @inferred Zygote._pullback(M16, X16)

        dy = Float32(0.125) * randn(Float32, size(y))
        DY32 = OneDNN.Memory(dy)
        DY16 = OneDNN.Memory(convert.(BFloat16, dy))

        grads_jl = back_jl(dy)
        grads_dnnl32 = back_dnnl32(DY32)
        grads_dnnl16 = back_dnnl16(DY16)

        @test isa(grads_jl, Tuple)
        @test isa(grads_dnnl32, Tuple)
        @test isa(grads_dnnl16, Tuple)

        # N.B.: Notice the difference in the plurality of the `weights` field ...
        # In the OneDNN code, we try to stay with the plural version since that's what the
        # oneDNN source code generally does.
        @test isapprox(
            grads_jl[1].weight, transpose(OneDNN.materialize(grads_dnnl32[1].weights))
        )
        @test isapprox(
            grads_jl[1].weight, transpose(OneDNN.materialize(grads_dnnl16[1].weights))
        )

        @test isapprox(grads_jl[1].bias, OneDNN.materialize(grads_dnnl32[1].bias))
        @test isapprox(grads_jl[1].bias, OneDNN.materialize(grads_dnnl16[1].bias))
        @test isapprox(grads_jl[2], OneDNN.materialize(grads_dnnl32[2]))
        @test isapprox(grads_jl[2], OneDNN.materialize(grads_dnnl16[2]))
        @inferred back_dnnl32(DY32)
        @inferred back_dnnl16(DY16)
    end
end


text:
@testset "generic_labelling 2D" begin

   for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))
            # Test using 8-connected neighbourhood
            labels = ImageComponentAnalysis.label_components(Generic(), test_image, trues(3,3))
            num_components = maximum(labels)
            @test num_components == test_image_results[i]
        end
    end
end

@testset "generic_labelling 3D" begin
    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
         for i = 1:20
             test_volume = eval(Symbol("test_volume_$(i)"))
             labels = ImageComponentAnalysis.label_components(Generic(), test_volume, trues(3,3,3))
             num_components = maximum(labels)
             @test num_components == test_volume_results[i]
         end
     end
end


text:
######################################################################################
## log(g(x)) & E[log(g(x))]
######################################################################################
using CpelNano

# Configuration
config = CpelNano.CpelNanoConfig()

# Parameter vector
a = 20.0; b = 100.0; c = 20.0; ϕ = [a,b,c];

# Other parameters
rs = CpelNano.RegStruct();
rs.cpg_pos = collect(100:20:3000); rs.N = length(rs.cpg_pos);
rs.chrst = 1; rs.chrend = 3000; rs.ϕhat = ϕ;
rs.ρn = rand(0.001:0.001:0.1, rs.N); rs.dn = fill(20.0, rs.N - 1);

# Re-scale model
CpelNano.rscle_grp_mod!(rs)

# Get matrices
CpelNano.get_rs_lgtrck_mats!(rs)

# Partition function
CpelNano.get_rs_logZ!(rs)

# Analysis region info
CpelNano.get_nls_reg_info!(rs,config)

# Get them all
CpelNano.get_rs_log_gs!(rs)

# Compute E[log g(X)]
CpelNano.get_rs_exp_log_g1!(rs)
CpelNano.get_rs_exp_log_g2!(rs)
rs.exps.log_g1
rs.exps.log_g2

# Compute E[X]
rs.exps.ex = CpelNano.get_E_X_log(rs.tm.log_u1, rs.tm.log_uN, rs.tm.log_Ws, rs.logZ)
rs.exps.exx = CpelNano.get_E_XX_log(rs.tm.log_u1, rs.tm.log_uN, rs.tm.log_Ws, rs.logZ)

# Compute MML
CpelNano.comp_mml!(rs)
rs.mml

# Compute NME
CpelNano.comp_nme!(rs)
rs.nme

# Compute CMD
CpelNano.comp_cmd(rs,rs)

######################################################################################
## Cross-entropy
######################################################################################
using Random
using CpelNano

# Configuration
config = CpelNano.CpelNanoConfig()

# Genomic properties
Random.seed!(123);
chrst = 1; chrend = 3000;
cpg_pos = collect(25:25:2975);
ρn = rand(0.001:0.001:0.1, length(cpg_pos));
dn = fill(25.0, length(cpg_pos) - 1)

# Parameter vectors
a1 = -10.0; b1 = -100.0; c1 = 20.0; ϕ1 = [a1,b1,c1];
a2 = 10.0; b2 = 100.0; c2 = 20.0; ϕ2 = [a2,b2,c2];

# Scale
ϕ1 = 0.5 .* ϕ1
ϕ2 = 0.5 .* ϕ2

## Region genomic properties

# Init structure 1
rs1 = CpelNano.RegStruct();
rs1.cpg_pos = cpg_pos; rs1.N = length(cpg_pos); rs1.chrst = chrst;
rs1.chrend = chrend; rs1.ϕhat = ϕ1; rs1.ρn = ρn; rs1.dn = dn;

# Init structure 2
rs2 = CpelNano.RegStruct();
rs2.cpg_pos = cpg_pos; rs2.N = length(cpg_pos); rs2.chrst = chrst;
rs2.chrend = chrend; rs2.ϕhat = ϕ2; rs2.ρn = ρn; rs2.dn = dn;

# Re-scale models
CpelNano.rscle_grp_mod!(rs1); CpelNano.rscle_grp_mod!(rs2);

# Get matrices
CpelNano.get_rs_lgtrck_mats!(rs1); CpelNano.get_rs_lgtrck_mats!(rs2);

# Partition function
CpelNano.get_rs_logZ!(rs1); CpelNano.get_rs_logZ!(rs2);

# Analysis region info
CpelNano.get_nls_reg_info!(rs1,config); CpelNano.get_nls_reg_info!(rs2, config);

# Get them all
CpelNano.get_rs_log_gs!(rs1); CpelNano.get_rs_log_gs!(rs2);

# Compute E[log g(X)]
CpelNano.get_rs_exp_log_g1!(rs1); CpelNano.get_rs_exp_log_g2!(rs1);
CpelNano.get_rs_exp_log_g1!(rs2); CpelNano.get_rs_exp_log_g2!(rs2);

# Compute E[X] and E[XX]
CpelNano.get_rs_exps!(rs1); CpelNano.get_rs_exps!(rs2);

# Mixture
mix = CpelNano.RegStruct()
mix.cpg_pos = rs1.cpg_pos; mix.N = rs1.N;
mix.ρn = ρn; mix.dn = dn; mix.nls_rgs = rs1.nls_rgs;
mix.Nl = rs1.Nl; mix.L = rs1.L; mix.ρl = rs1.ρl;
mix.dl = rs1.dl; mix.ϕhat = 0.5 .* (rs1.ϕhat + rs2.ϕhat);

# Update rest of information (CHECK WHAT'S NEEDED)
CpelNano.get_rs_lgtrck_mats!(mix)
CpelNano.get_rs_logZ!(mix)
CpelNano.get_rs_exps!(mix)
CpelNano.get_rs_log_gs!(mix)

# Mixture
exps_log_g1 = CpelNano.comp_crs_exps_log_g1(rs1, mix)
exps_log_g2 = CpelNano.comp_crs_exps_log_g2(rs2, mix)

# Cross-entropy w/ θ1 & θ2
cross_1 = CpelNano.comp_crs_ent(rs1, mix)
cross_2 = CpelNano.comp_crs_ent(rs2, mix)

######################################################################################
## CMD
######################################################################################
using Random
using CpelNano

# Configuration
config = CpelNano.CpelNanoConfig()

# Genomic properties
Random.seed!(123);
chrst = 1; chrend = 3000;
cpg_pos = collect(25:25:2975);
ρn = rand(0.001:0.001:0.1, length(cpg_pos));
dn = fill(25.0, length(cpg_pos) - 1)

# Parameter vectors
a1 = -10.0; b1 = -100.0; c1 = 20.0; ϕ1 = [a1,b1,c1];
a2 = 10.0; b2 = 100.0; c2 = 20.0; ϕ2 = [a2,b2,c2];

# Scale
ϕ1 = 0.5 .* ϕ1
ϕ2 = 0.5 .* ϕ2

## Region genomic properties

# Init structure 1
rs1 = CpelNano.RegStruct();
rs1.cpg_pos = cpg_pos; rs1.N = length(cpg_pos); rs1.chrst = chrst;
rs1.chrend = chrend; rs1.ϕhat = ϕ1; rs1.ρn = ρn; rs1.dn = dn;

# Init structure 2
rs2 = CpelNano.RegStruct();
rs2.cpg_pos = cpg_pos; rs2.N = length(cpg_pos); rs2.chrst = chrst;
rs2.chrend = chrend; rs2.ϕhat = ϕ2; rs2.ρn = ρn; rs2.dn = dn;

# Re-scale models
CpelNano.rscle_grp_mod!(rs1); CpelNano.rscle_grp_mod!(rs2);

# Get matrices
CpelNano.get_rs_lgtrck_mats!(rs1); CpelNano.get_rs_lgtrck_mats!(rs2);

# Partition function
CpelNano.get_rs_logZ!(rs1); CpelNano.get_rs_logZ!(rs2);

# Analysis region info
CpelNano.get_nls_reg_info!(rs1,config); CpelNano.get_nls_reg_info!(rs2, config);

# Get them all
CpelNano.get_rs_log_gs!(rs1); CpelNano.get_rs_log_gs!(rs2);

# Compute E[log g(X)]
CpelNano.get_rs_exp_log_g1!(rs1); CpelNano.get_rs_exp_log_g2!(rs1);
CpelNano.get_rs_exp_log_g1!(rs2); CpelNano.get_rs_exp_log_g2!(rs2);

# Compute E[X] and E[XX]
CpelNano.get_rs_exps!(rs1); CpelNano.get_rs_exps!(rs2);

# Compute MML
CpelNano.comp_mml!(rs1); CpelNano.comp_mml!(rs2);
round.([rs1.mml rs2.mml], digits=4)

# Compute NME
CpelNano.comp_nme!(rs1); CpelNano.comp_nme!(rs2);
round.([rs1.nme rs2.nme], digits=4)

# Compute CMD
cmd = CpelNano.comp_cmd(rs1, rs2)
round.(cmd, digits=4)


text:
abstract type CompiledAction <: Action end

# Utility functions
include("utils.jl")
# Compiled types
include("domain.jl")
include("state.jl")
include("action.jl")
# Helper functions for objects/accessors/formulas
include("objects.jl")
include("accessors.jl")
include("formulas.jl")
# Interface functions
include("satisfy.jl")
include("evaluate.jl")
include("initstate.jl")
include("transition.jl")
include("available.jl")
include("execute.jl")

"""
    compiled(domain, state)
    compiled(domain, problem)

Compile a `domain` and `state` and return the resulting compiled domain
and compiled state. A `problem` maybe provided instead of a state.
Note that this function must be called at the top-level in order to avoid
world-age errors.
"""
function compiled(domain::Domain, state::State)
    # Generate definitions
    domain_type, domain_typedef, domain_defs =
        generate_domain_type(domain, state)
    state_type, state_typedef, state_defs =
        generate_state_type(domain, state, domain_type)
    object_defs =
        generate_object_defs(domain, state, domain_type, state_type)
    evaluate_def =
        generate_evaluate(domain, state, domain_type, state_type)
    satisfy_def =
        generate_satisfy(domain, state, domain_type, state_type)
    action_defs =
        generate_action_defs(domain, state, domain_type, state_type)
    transition_def =
        generate_transition(domain, state, domain_type, state_type)
    # Generate return expression
    return_expr = :($domain_type(), $state_type($state))
    # Evaluate definitions
    expr = Expr(:block,
        domain_typedef, domain_defs, state_typedef, state_defs, object_defs,
        evaluate_def, satisfy_def, action_defs, transition_def, return_expr)
    return eval(expr)
end

function compiled(domain::Domain, problem::Problem)
    return compiled(domain, initstate(domain, problem))
end

"""
    compilestate(domain, state)

Return compiled version of a state compatible with the compiled `domain`.
"""
function compilestate(domain::CompiledDomain, state::State)
    S = statetype(domain)
    return S(state)
end

# Abstract a domain that is already compiled
function abstracted(domain::CompiledDomain, state::State; options...)
    absdom = abstracted(get_source(domain); options...)
    return compiled(absdom, GenericState(state))
end

# If compiled state is abstract, we can just copy construct
function abstractstate(domain::CompiledDomain, state::State)
    S = statetype(domain)
    return S(state)
end


text:
function vech_to_vec(n::S) where S <: Int

  a = zeros(Int,n,n)
  tracker = 1

  for r = 1:n

    for c = 1:n

      if c >= r
        a[r,c] = tracker
        tracker += 1
      else
        a[r,c] = a[c,r]
      end

    end

  end

  gi = vec(a)

  x = zeros(Int,n*n,round(Int,(n*(n+1)/2)))

  for r = 1:size(x,1)
    x[r,gi[r]] = 1
  end

  return x

end