text:
"""
General time profile type (gradient sequence).

A `TimeProfile` should implement a calling method returning the value of the
time profile at time `t`. It can optionally overwrite the `echotime` `integral`,
and `bvalue_no_q` methods, the first providing a default value and the latter
two computing numerical approximation to the integral quantities if not
overwritten.
"""
abstract type TimeProfile{T<:Real} <: Function end

"""
    f = PGSE(δ, Δ)

Pulsed Gradient Spin Echo sequence with pulse duration `δ` and time between
pulses `Δ-δ`.
"""
struct PGSE{T} <: TimeProfile{T}
    δ::T
    Δ::T
    function PGSE(δ::T, Δ::T) where {T}
        @assert δ ≤ Δ
        new{T}(δ, Δ)
    end
end
PGSE(δ::Int, Δ::Int) = PGSE(float(δ), float(Δ))
PGSE(δ::Int, Δ) = PGSE(promote(δ, Δ)...)
PGSE(δ, Δ::Int) = PGSE(promote(δ, Δ)...)

"""
    f = CosOGSE(δ, Δ, nperiod)

Oscillating Gradient Spin Echo sequence with two cos-pulses of duration `δ`
separated by a pause of duration `Δ-δ` for `nperiod` periods per pulse.
"""
struct CosOGSE{T} <: TimeProfile{T}
    δ::T
    Δ::T
    nperiod::Int
    function CosOGSE(δ::T, Δ::T, n) where {T}
        @assert δ ≤ Δ
        @assert n > zero(n)
        new{T}(δ, Δ, n)
    end
end

"""
    f = SinOGSE(δ, Δ, nperiod)

Oscillating Gradient Spin Echo sequence with two sin-pulses of duration `δ`
separated by a pause of duration `Δ-δ` for `nperiod` periods per pulse.
"""
struct SinOGSE{T} <: TimeProfile{T}
    δ::T
    Δ::T
    nperiod::Int
    function SinOGSE(δ::T, Δ::T, n) where {T}
        @assert δ ≤ Δ
        @assert n > zero(n)
        new{T}(δ, Δ, n)
    end
end

"""
    f = DoublePGSE(δ, Δ)

Double Pulsed Gradient Spin Echo sequence with four pulses of duration `δ`,
separated by pauses of duration `Δ-δ`, `0`, and `Δ-δ` repsectively.
"""
struct DoublePGSE{T} <: TimeProfile{T}
    δ::T
    Δ::T
    tpause::T
    function DoublePGSE(δ::T, Δ::T, tpause::T) where {T}
        @assert δ ≤ Δ
        @assert tpause ≥ zero(tpause)
        new{T}(δ, Δ, tpause)
    end
end

"""
    ft = (f::TimeProfile)(t)

Return the time profile value at time `t`.
"""
function (f::TimeProfile)(t)
    error("Not implemented")
end

function (f::PGSE)(t)
    (zero(t) ≤ t < f.δ) - (f.Δ ≤ t < f.Δ + f.δ)
end

function (f::DoublePGSE)(t)
    δ, Δ, p = f.δ, f.Δ, f.tpause
    (
        (zero(t) ≤ t < δ) - (Δ ≤ t < Δ + δ) + (p + Δ + δ ≤ t < p + Δ + 2δ) -
        (p + 2Δ + δ ≤ t < p + 2Δ + 2δ)
    )
end

function (f::CosOGSE)(t)
    δ, Δ, n = f.δ, f.Δ, f.nperiod
    (t < δ) * cos(2π * nperiod * t / δ) - (Δ ≤ t) * cos(2π * nperiod * (t - Δ) / δ)
end

function (f::SinOGSE)(t)
    δ, Δ, n = f.δ, f.Δ, f.nperiod
    (t < δ) * sin(2π * nperiod * t / δ) - (Δ ≤ t) * sin(2π * nperiod * (t - Δ) / δ)
end

"""
    integral(f, t = echotime(f))

Integral of time profile `f` between `0` and `t`. Unless specified, the echotime
is used as the upper integral limit.

For the `PGSE`, `SinOGSE`, `CosOGSE` and `DoublePGSE` sequences, analytical
expressions are available. Otherwise a numerical integral is computed.
"""
function integral(f::TimeProfile, t = echotime(f))
    quadgk(f, zero(t), t)
end

function integral(f::PGSE, t = echotime(f))
    δ, Δ = f.δ, f.Δ
    (zero(t) ≤ t < δ) * t + (δ ≤ t < Δ + δ) * δ - (Δ ≤ t < Δ + δ) * (t - Δ)
end

function integral(f::DoublePGSE, t = echotime(f))
    δ, Δ, p = f.δ, f.Δ, f.tpause
    tmid = p + Δ + δ
    if zero(t) ≤ t < δ
        t
    elseif δ ≤ t < Δ
        δ
    elseif Δ ≤ t < Δ + δ
        δ - (t - Δ)
    elseif tmid ≤ t < tmid + δ
        t - tmid
    elseif tmid + δ ≤ t < tmid + Δ
        δ
    elseif tmid + Δ ≤ t < tmid + Δ + δ
        δ - (t - Δ - tmid)
    else
        zero(t)
    end
end

function integral(f::CosOGSE, t = echotime(f))
    δ, Δ, n = f.δ, f.Δ, f.nperiod
    ((t < δ) * sin(2π * n * t / δ) - (Δ ≤ t) * sin(2π * n * (t - Δ) / δ)) * δ / (2π * n)
end

function integral(f::SinOGSE, t = echotime(f))
    δ, Δ, n = f.δ, f.Δ, f.nperiod
    ((t < δ) * (1 - cos(2π * n * t / δ)) - (Δ ≤ t) * (1 - cos(2π * n * (t - Δ) / δ))) * δ /
    (2π * n)
end

"""
    b = bvalue_no_q(f)

Compute the time profile contribution to the b-value. To obtain the b-value,
multiply the result by `q^2 = (γg)^2`, where `γ` is the gyromagnetic ratio of the
water proton, `g` is the gradient amplitude, and `b = q^2 * bvalue_no_q(f)`.
"""
function bvalue_no_q(f::TimeProfile{T}) where {T}
    quadgk(s -> integral(f, s)^2, zero(T), echotime(f))
end

function bvalue_no_q(f::PGSE)
    f.δ^2 * (f.Δ - f.δ / 3)
end

function bvalue_no_q(f::DoublePGSE)
    2f.δ^2 * (f.Δ - f.δ / 3)
end

function bvalue_no_q(f::CosOGSE)
    f.δ^3 / (2π * f.nperiod)^2
end

function bvalue_no_q(f::SinOGSE)
    3f.δ^3 / (2π * f.nperiod)^2
end

"""
    intervals(f)

Get characteristic intervals of the time profile `f`.
"""
function intervals(f::TimeProfile)
    error("Not implemented")
end

function intervals(f::Union{PGSE,CosOGSE,SinOGSE})
    δ, Δ = f.δ, f.Δ
    if δ < Δ
        [zero(δ), δ, Δ, Δ + δ]
    else
        # No pause between pulses
        [zero(δ), δ, 2f.δ]
    end
end

function intervals(f::DoublePGSE)
    δ, Δ, p = f.δ, f.Δ, f.tpause
    if δ < Δ
        i = [zero(δ), δ, Δ, Δ + δ]
    else
        # No pause between pulses
        i = [zero(δ), δ, 2f.δ, 3f.δ, 4f.δ]
    end
    if p > zero(p)
        [i; i .+ p .+ Δ .+ δ]
    else
        [i; i[2:end] .+ p .+ Δ .+ δ]
    end
end

"""
    TE = echotime(f)

Get echo time `TE` of the time profile `f`, which is the end of the last characteristic
interval.
"""
function echotime(f::TimeProfile)
    intervals(f)[end]
end

"""
    get_interval(f, t)

Get the characteristic interval number of time profile `f` containing time `t`.
If `intervals[i] ≤ t < intervals[i+1]`, `i` is returned.
"""
function get_interval(f::TimeProfile, t)
    i = intervals(f)
    findfirst(i[1:end-1] .≤ t .< i[2:end])
end

"""
    markers = constant_intervals(f)

Get markers for constant intervals. Default is false on all intervals.
"""
function constant_intervals(f::TimeProfile)
    fill(false, length(intervals(f)) - 1)
end

function constant_intervals(f::Union{PGSE,DoublePGSE})
    fill(true, length(intervals(f)) - 1)
end

function constant_intervals(f::Union{CosOGSE,SinOGSE})
    if f.δ < f.Δ
        # OGSE is constant between the pulses
        c = [false, true, false]
    else
        # No pause between pulses
        c = [false, false]
    end
    c
end

"""
    is_constant(f, t)

Return `true` if time profile `f` is constant on the characteristic interval containing or
starting at `t`.
"""
function is_constant(f::TimeProfile, t)
    constant_intervals(f)[get_interval(f, t)]
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ b5a8675e-0e2a-11eb-1787-55285590d7de
using Plots

# ╔═╡ bcd20d70-0e2a-11eb-34c1-5f49700df9cd
using LaTeXStrings

# ╔═╡ 55991350-0e2b-11eb-3419-0f9f9ac0293f
using PlutoUI

# ╔═╡ 0cbfa5d0-2506-11eb-1485-7733b85a975b
using SymEngine
#hay muchos paquete simbólicos, el más completo (límites integrales, sustiuciones...) es SymPy, Symbolics es nuevo y en desarrollo y SymEngine es simple y suficiente para lo que vamos a hacer

# ╔═╡ 5aded872-b8b4-11eb-0277-93d61bffe2fd
md"# Tema 5 Control en modo deslizante"

# ╔═╡ f8105520-250b-11eb-1a6b-0129a812bcee
md"""# Vamos a crear nuetro propio simulador con *Euler*
¿En serio? ¿Euler no era un método malo?
No es el mejor método pero... los otros tampoco funcionan bien si la función es discontinua y a nuestro *le podemos hacer perrerías* (usar funciones discontinuas, introducir ruido, etc...)
"""

# ╔═╡ 96dbe960-2509-11eb-1248-b7c75a8974dd
#uso euler por que en un sistema discontinuo el resto de solvers no ofrece mucha ventaja y así tenemos más control sobre lo que está pasando

function simular(x0,planta,control,parametros,tf,Δt)
	t=0:Δt:tf
	estado=copy(x0);
	x=[]
	push!(x,estado)
	u=Float64[control(estado,t[1],parametros)[1]]

	for n=1:(length(t)-1)
		señal_de_control=control(estado,t[n],parametros)[1] #el [1] se entenderá luego
		derivada=planta(estado,señal_de_control,t[n],parametros)
		estado=estado+derivada*Δt #euler
		push!(x,estado)
		push!(u,señal_de_control)
	end

	return (x,u,t)
end

# ╔═╡ ba989850-0a31-11eb-1f3a-e363a7e0207e
md"""# Primer ejemplo (primer orden):

$\dot x=sin(t)+acos(x)+u=f(x,t)+u$

$a \in [0.5,0.9]$

Luego sabemos que

$|f(x,t)|<|sin(t)|+|acos(x)| \le 1+|a| \le 1+0.9=1.9=\rho$

## Hay que aplicar un $u$ mayor que $\rho$ y con signo opuesto a $x$... luego haremos esto:

$u=-sign(x)(\rho+\epsilon)=-2sign(x)$

Donde hemos tomado $\epsilon=0.1$
"""

# ╔═╡ 94294500-0e2a-11eb-3a25-617242dfe932
function planta_ejemplo1(x,u,t,p)
	a=p[1]#parámetro desconocido
	dx=sin(t)+a*cos(x)+u
end

# ╔═╡ 1d402de0-0e2b-11eb-10e8-79bdc072ee9b
control_1a(x,t,p)=-2*sign(x)

# ╔═╡ 3b89eed0-0e2b-11eb-3cad-8d86b96accc5
@bind x1_0 Slider(-1:0.1:1)

# ╔═╡ 6b48fc60-0e2b-11eb-1ad1-075f865e4c78
@bind a1 Slider(0.5:0.01:0.9)

# ╔═╡ f77e1fc0-0e2c-11eb-0051-d7eee614b171
md"""# Hemos hecho *trampa*, ¿Dónde está el control equivalente?
El control equivalente es el que mantiene $s=x$ en cero es decir:

$x(t)=0 \to  \dot x = 0 =sin(t)+acos(x)+u=sin(t)+acos(0)+u$
$0=sin(t)+a+u \to u=-a-sin(t)$

No está ni hace falta, ¡Nunca se ejecuta! Vamos a verlo
"""

# ╔═╡ e3723e20-0efa-11eb-1331-e710fbcd289f
function control_1b(x,t,p)
	if x>0
		u=-2
	elseif x<0
		u=2
	else
		u=1000000 #NO es el control equivalente es para ver que sucede ;)
	end
	return u
end

# ╔═╡ 7cc77a40-0efb-11eb-145d-4b535e701cf8
md"""no ha pasado nada, u nunca vale 1000000
... a no ser que x sea exactamente cero lo que es muy raro a no ser que lo forcemos
(poniendo $x_0=0$) pero ni así tiene mucha improtancia

## Significado del control equivalente

Aunque no lo veamos el control equivalente es "en promedio" igual que el control que estamos aplicando, vamos a verlo
"""

# ╔═╡ 29de3480-25c7-11eb-1d21-2107f4712127
media_móvil(v,n) = [sum(@view v[i:(i+n-1)])/n for i in 1:(length(v)-(n-1))]

# ╔═╡ 1ce32ff0-25c8-11eb-1371-e360c5c3ed2b
@bind N Slider(1:30)

# ╔═╡ 674476e0-2810-11eb-000b-0742e1bce381
md"Ver real $(@bind ver CheckBox())"

# ╔═╡ e32619f0-0efa-11eb-020d-c9743701e335
let
	Δt=0.001
	tf=2
	x0=x1_0
	parametros=[a1]
	x,u,t=simular(x0,planta_ejemplo1,control_1a,parametros,tf,Δt)

	u_equivalente=-a1.-sin.(t)

	u_promedio=media_móvil(u,N)
	t_promedio=media_móvil(t,N)

	fig1=plot()
	if ver==true
	   fig1=plot!(fig1,t,u, xaxis=L"t",yaxis="ues",label="u_real")
	end
	fig1=plot!(fig1,t,u_equivalente,label="u_equivalente")
	fig1=plot!(fig1,t_promedio,u_promedio,label="u_promedio")
end

# ╔═╡ 14080630-0f07-11eb-068d-09c4055d081b
md"""# Regularización y control equivalente
En vez de usar el signo que es discontinuo podemos usar una versión continua del mismo
"""

# ╔═╡ f47ce110-24fe-11eb-036a-e5d8bad1955f
@bind ϵ Slider(0.01:0.01:0.2)   #\epsilon + TAB

# ╔═╡ 1f38d210-24ff-11eb-3cbe-9b1e70108439
let
	ordenadas=-1:0.01:1
	fig1=plot(ordenadas,sign.(ordenadas), xaxis=L"x",label="sign(x)")
	fig1=plot!(fig1,ordenadas,tanh.(ordenadas/ϵ), label="tahn(x/$ϵ)")

end

# ╔═╡ e71e9160-2500-11eb-1b6c-d59e43143b66
md"Veamos que efecto tiene camniar el signo por su vesion continua"

# ╔═╡ 24a65e50-2501-11eb-2363-67d7e5bcb405
control_1c(x,t,p)=-2*tanh(x/p[2])

# ╔═╡ e01fcbe0-2500-11eb-3b94-a7662624c5ed
md"""## Entonces lo mejor es disminuir $\epsilon$ todo lo posible ¿No?

*No exactamente*, en la práctica tiene sus problemas, por ejemplo si la medida del estado es ruidosa...
Veámoslo:
"""

# ╔═╡ 16e8dac2-25cb-11eb-33d1-a95dc487b34d
function control_1d(x,t,p)
	_,ϵ,nivel_ruido=p
	ruido=nivel_ruido*2*(rand()-0.5)
	x_medido=x+ruido
	u=-2*tanh(x_medido/ϵ)
end

# ╔═╡ 612d43c0-2506-11eb-3a14-8145b79e349f
@bind ϵ2 Slider(0.01:0.01:0.2)

# ╔═╡ d8e35250-25ca-11eb-3d59-b3c6396a4fb1
@bind nivel_ruido Slider(0:0.01:0.2)

# ╔═╡ f694c060-2505-11eb-26ef-1fc536e2260d
md"""# Ejemplo 2 vamos al tracking directamente ☺
La estabilización es como el tracking pero con $r=0$, como ya lo vimos en el tema anterior podemos ir más rápido en este.

$\dot x_1=sin(x_1)+x_2$

$\dot x_2=ax_1+bx_2^2+cu$

$ y=x_1 $

con

$ a \in [-1,7] $,  $ b \in [1,2] $ y $ c \in [1.5,2.8] $

Queremos seguir la referenica $ y_r=sin(t) $

"""

# ╔═╡ 595542f0-25cf-11eb-308d-bb524a624b04
md"Defimos las función igual que en el tema 4"

# ╔═╡ 41d2207e-25cf-11eb-1cbf-a599e3a5d2fc
begin
	function derivada(funcion,variables,derivadas)
		d=0
		for i=1:length(variables)
			d=d+diff(funcion,variables[i])*derivadas[i]
		end
		return d
	end

	function derivada(f,x,dx,n) #multiple distpach
		d=f
		for i=1:n
			d=derivada(d,x,dx)
		end
		return d
	end
end

# ╔═╡ 32f2cfd0-2505-11eb-36a3-cdaec55a91cc
x1,x2,u,a,b,c=symbols("x1,x2,u,a,b,c") #o @vars x1... pero inicializando antes

# ╔═╡ ff0eafe0-0e2a-11eb-3e1c-5dd70c43c3c1
let
	Δt=0.001
	tf=2
	x0=x1_0
	parametros=[a1]
	x,u,t=simular(x0,planta_ejemplo1,control_1a,parametros,tf,Δt)

	fig1=plot(t,x, xaxis=L"t",yaxis=L"x(t)",label=:none)
	title!(fig1,"a=$a1")
	fig2=plot(t,u, xaxis=L"t",yaxis=L"u(t)",legend=false)
	l = @layout [a ; b]
	plot(fig1,fig2, layout=l)
end

# ╔═╡ 155f6de0-25c6-11eb-3a8e-851538372b3e
let
	Δt=0.001
	tf=2
	x0=x1_0
	parametros=[a1]
	x,u,t=simular(x0,planta_ejemplo1,control_1b,parametros,tf,Δt)

	fig1=plot(t,x, xaxis=L"t",yaxis=L"x(t)",label=:none)
	fig1=hline!([0],line=:dot,label=:none)
	title!(fig1,"a=$a1")
	fig2=plot(t,u, xaxis=L"t",yaxis=L"u(t)",legend=false)
	l = @layout [a ; b]
	plot(fig1,fig2, layout=l)
end

# ╔═╡ e023ea92-2500-11eb-00ff-4f27947f46d4
let
	Δt=0.001
	tf=2
	x0=x1_0
	parametros=[a1,ϵ]
	x,u,t=simular(x0,planta_ejemplo1,control_1c,parametros,tf,Δt)

	u_equivalente=-a1.-sin.(t)

	fig1=plot(t,x, xaxis=L"t",yaxis=L"x(t)",label=:none)
	title!(fig1,"a=$a1 ϵ=$ϵ")
	fig1=hline!([0],line=:dot,label=:none)

	fig2=plot(t,u, xaxis=L"t",yaxis=L"u_{aplicado}",legend=false)
	fig2=plot!(fig2,t,u_equivalente,
		    xaxis=L"t",yaxis=L"u_{es}",label=L"u_{equivalente}")

	l = @layout [a ; b]
	plot(fig1,fig2, layout=l)
end

# ╔═╡ 5a31c47e-2504-11eb-37af-07b3b88360eb
let
	Δt=0.001
	tf=2
	x0=x1_0
	parametros=[a1,ϵ2,nivel_ruido]
	x,u,t=simular(x0,planta_ejemplo1,control_1d,parametros,tf,Δt)

	u_equivalente=-a1.-sin.(t)

	fig1=plot(t,x, xaxis=L"t",yaxis=L"x(t)",label=:none)
	title!(fig1,"a=$a1, ϵ=$ϵ, nivel del ruido=$nivel_ruido")
	fig1=hline!([0],line=:dot,label=:none)

	fig2=plot(t,u, xaxis=L"t",yaxis=L"u_{aplicado}",legend=false)
	fig2=plot!(fig2,t,u_equivalente,
		    xaxis=L"t",yaxis=L"u_{es}",label=L"u_{equivalente}")

	l = @layout [a ; b]
	plot(fig1,fig2, layout=l)
end

# ╔═╡ b1c43a50-25ce-11eb-3478-f13e9174ccc9
begin
    dx1=sin(x1)+x2
    dx2=a*x1+b*x2^2+c*u
	dx=[dx1 dx2]
	y=x1
end

# ╔═╡ 992e79f0-25cf-11eb-0885-9f431d9e2ed6
dy=derivada(y,[x1 x2],[dx1 dx2],1)

# ╔═╡ 806542f0-25cf-11eb-39c8-db55ed5e86a7
ddy=derivada(y,[x1 x2],[dx1 dx2],2)

# ╔═╡ c5f21a70-25d2-11eb-2af8-395671e44537
md"Ok tiene el grado correcto podemos usar la técnica"

# ╔═╡ 2cd935f0-25d0-11eb-1b07-a7002584aa16
md"""## Defimos ahora el error y su dinámica

Recordemos $e_1 = y - y_r$  y ahora a calcular sus derivadas:
"""

# ╔═╡ 77f468f0-2816-11eb-2dc5-25eb17f3a834
yr,dyr,ddyr=symbols("yr,dyr,ddyr")

# ╔═╡ e1f2879e-25d2-11eb-04d4-b53a77d90c5f
e1=yr-y

# ╔═╡ 14e23520-2817-11eb-2f93-0779085f77a5
md"""La teoría es que $S=(\frac{d}{dt}+k)^{n-1}e_1$

Para el caso n=2

$S=\dot e_1 +ke_1$"""

# ╔═╡ 37fe65b0-2817-11eb-2cdb-2921b513978c
k=symbols("k")

# ╔═╡ 9ac55fb0-2816-11eb-12da-a9ddac1c913c
de1=derivada(e1,[x1 x2 yr dyr],[dx1 dx2 dyr ddyr])

# ╔═╡ 0ed1a02e-2817-11eb-0b92-e5feb0f17c63
S=de1+k*e1

# ╔═╡ a3caae20-2817-11eb-09e6-59d6464cef2a
dS=derivada(S,[x1 x2 yr dyr],[dx1 dx2 dyr ddyr])

# ╔═╡ aff35570-25cf-11eb-1627-43fe902b9db5
md"Identificamos términos $\dot S=-f-gu$"

# ╔═╡ a8a18df2-25cf-11eb-38b7-e128e53ba029
g=-diff(dS,u)

# ╔═╡ 081bb16e-25d0-11eb-28e6-c3e5a042fce2
f=-(dS+g*u) #el resto vamos.... OJO A LOS SIGNOS

# ╔═╡ f77c9930-25d2-11eb-287f-655c63ea721e
md"""Acotar $|f/g|$ un proceso **artesanal** se puede hacer **de muchas formas**
recordemos que $ a \in [-1,7] $,  $ b \in [1,2] $ y $ c \in [1.5,2.8] $
"""

# ╔═╡ cf8b68e2-b8ba-11eb-0177-8fa5c7f332bd
expand(f/g)

# ╔═╡ 097834c0-b8bb-11eb-3e52-c1d543ef9fb6
md"Y ahora vamos a acotar"

# ╔═╡ 9c2a7480-25d2-11eb-0479-751c0dc8b88e
rho=(abs(ddyr) + 7*abs(x1) + 2*x2^2 + k*abs(dyr)+ k*abs(x2) + k + abs(x2) + 1)/1.5

# ╔═╡ abc0604e-2823-11eb-1708-83268e22b7bc
md"## Todo junto"

# ╔═╡ 2b9f72b0-281b-11eb-1fbc-15b3e67b017d
referencia(t)=0.0.*(sin(t),cos(t),-sin(t))

# ╔═╡ 42884660-281a-11eb-015e-29ff1b21dfe0
function control_seguimiento(x,t,p)
	ϵ=0.01;
    x1,x2=x
	ε=0.1;#el epsilon bonito es varepsilon
    k=1;#una elección arbitraria
    yr,dyr,ddyr=referencia(t)

	y=x1; #mucho cuidado si no se define puede usar las globales y se cuelga
	e1=yr-y
	de1=dyr - (x2 + sin(x1))

	S=de1+k*e1

	#RELLENAR ESTO
	u=tanh(S / ϵ)*((abs(ddyr) + 7*abs(x1) + 2*x2^2 + k*abs(dyr)+ k*abs(x2) + k + abs(x2) + 1)/1.5 + ε)

	return (u,S,e1,de1)  #Saco todo esto para pintarlo luego fácilmente
	#por eso la función simular se queda SOLO CON u
end

# ╔═╡ 9acf5880-281b-11eb-3aa5-9b6df9ac7501
function planta_ejemplo2(x,u,t,p)
	a,b,c=p
	x1,x2=x #OJO a las variables, si por ejemplo no las ponemos...
    dx1=sin(x1)+x2
    dx2=a*x1+b*x2^2+c*u
	dx=[dx1 dx2]
end

# ╔═╡ e8ee4712-281b-11eb-0fbf-7783935a0075
@bind a2 Slider(-1:0.1:7)

# ╔═╡ 0e04e810-281c-11eb-3e72-c58bf8d7ce0d
@bind b2 Slider(1:0.1:2)

# ╔═╡ 1e143c10-281c-11eb-1117-29f72834203e
@bind c2 Slider(1.5:0.1:2.8)

# ╔═╡ 756a3570-2833-11eb-2596-c37768e5a75b
@bind x1_0b Slider(-1:0.1:1)

# ╔═╡ be9d11e0-2833-11eb-1e9f-d3bfc428c893
@bind variables Select(["estados"=>"Ver estados", "errores"=>"Ver errores"])

# ╔═╡ 8bdbadb0-281b-11eb-1786-3bef9aee1d39
let
	Δt=0.001
	tf=10
	x0=[x1_0b 0]
	parametros=[a2,b2,c2]
	estados,u,t=simular(x0,planta_ejemplo2,control_seguimiento,parametros,tf,Δt)

	#hay que sacar los datos de los estados, es decir parsar de un vecror de vectores a vectores sencillos
	x1=[x[1] for x in estados]
	x2=[x[2] for x in estados]
	y=x1

	refs=referencia.(t)
	yr=[ref[1] for ref in refs]
	dyr=[ref[2] for ref in refs]

	#truco para pintar más cosas sin tener que reescribirlas
	Ss=[]
	e1s=[]
	de1s=[]
	for i=1:length(t)
		_,S,e1,de1=control_seguimiento(estados[i],t[i],parametros)
		push!(Ss,S)
		push!(e1s,e1)
        push!(de1s,de1)
	end

	if variables=="estados"
		plot(t,x1, xaxis=L"t",yaxis=L"x(t)",label=L"y=x_1")
	    plot!(t,x2, xaxis=L"t",label=L"x_2")
	    title!("a=$a2, b=$a2, c=$a2")
	    fig1=plot!(t,yr,line=:dot,label=L"y_r")
	else
		plot(t,e1s, xaxis=L"t",yaxis=L"e(t)",label=L"e_1")
	    plot!(t,de1s, xaxis=L"t",label=L"e_2=\dot e_1")
        title!("a=$a2, b=$a2, c=$a2")
		fig1=hline!([0],line=:dot,label=:none)
	end

	fig2=plot(t,Ss, xaxis=L"t",yaxis=L"s",legend=false)

	fig3=plot(t,u, xaxis=L"t",yaxis=L"u",legend=false)

    if variables=="estados"
		fig4=plot(x1,x2, xaxis=L"x_1", yaxis=L"x_2",label=:none)
	else
		fig4=plot(e1s,de1s, xaxis=L"e_1", yaxis=L"e_2",label=:none)
	end

	l = @layout [a ; b ;c ; d]
	plot(fig1,fig2,fig3, fig4, layout=l, size=(600,800))

end

# ╔═╡ 28d56f3e-b8c0-11eb-1329-21f533ccddd9

# ╔═╡ Cell order:
# ╟─5aded872-b8b4-11eb-0277-93d61bffe2fd
# ╟─f8105520-250b-11eb-1a6b-0129a812bcee
# ╠═96dbe960-2509-11eb-1248-b7c75a8974dd
# ╟─ba989850-0a31-11eb-1f3a-e363a7e0207e
# ╠═b5a8675e-0e2a-11eb-1787-55285590d7de
# ╠═bcd20d70-0e2a-11eb-34c1-5f49700df9cd
# ╠═94294500-0e2a-11eb-3a25-617242dfe932
# ╠═1d402de0-0e2b-11eb-10e8-79bdc072ee9b
# ╠═55991350-0e2b-11eb-3419-0f9f9ac0293f
# ╠═3b89eed0-0e2b-11eb-3cad-8d86b96accc5
# ╠═6b48fc60-0e2b-11eb-1ad1-075f865e4c78
# ╠═ff0eafe0-0e2a-11eb-3e1c-5dd70c43c3c1
# ╟─f77e1fc0-0e2c-11eb-0051-d7eee614b171
# ╠═e3723e20-0efa-11eb-1331-e710fbcd289f
# ╠═155f6de0-25c6-11eb-3a8e-851538372b3e
# ╟─7cc77a40-0efb-11eb-145d-4b535e701cf8
# ╠═29de3480-25c7-11eb-1d21-2107f4712127
# ╠═1ce32ff0-25c8-11eb-1371-e360c5c3ed2b
# ╟─674476e0-2810-11eb-000b-0742e1bce381
# ╠═e32619f0-0efa-11eb-020d-c9743701e335
# ╟─14080630-0f07-11eb-068d-09c4055d081b
# ╟─1f38d210-24ff-11eb-3cbe-9b1e70108439
# ╠═f47ce110-24fe-11eb-036a-e5d8bad1955f
# ╠═e71e9160-2500-11eb-1b6c-d59e43143b66
# ╠═24a65e50-2501-11eb-2363-67d7e5bcb405
# ╟─e023ea92-2500-11eb-00ff-4f27947f46d4
# ╟─e01fcbe0-2500-11eb-3b94-a7662624c5ed
# ╠═16e8dac2-25cb-11eb-33d1-a95dc487b34d
# ╟─5a31c47e-2504-11eb-37af-07b3b88360eb
# ╠═612d43c0-2506-11eb-3a14-8145b79e349f
# ╠═d8e35250-25ca-11eb-3d59-b3c6396a4fb1
# ╟─f694c060-2505-11eb-26ef-1fc536e2260d
# ╠═0cbfa5d0-2506-11eb-1485-7733b85a975b
# ╟─595542f0-25cf-11eb-308d-bb524a624b04
# ╟─41d2207e-25cf-11eb-1cbf-a599e3a5d2fc
# ╠═32f2cfd0-2505-11eb-36a3-cdaec55a91cc
# ╠═b1c43a50-25ce-11eb-3478-f13e9174ccc9
# ╠═992e79f0-25cf-11eb-0885-9f431d9e2ed6
# ╠═806542f0-25cf-11eb-39c8-db55ed5e86a7
# ╟─c5f21a70-25d2-11eb-2af8-395671e44537
# ╟─2cd935f0-25d0-11eb-1b07-a7002584aa16
# ╠═77f468f0-2816-11eb-2dc5-25eb17f3a834
# ╠═e1f2879e-25d2-11eb-04d4-b53a77d90c5f
# ╟─14e23520-2817-11eb-2f93-0779085f77a5
# ╠═37fe65b0-2817-11eb-2cdb-2921b513978c
# ╠═9ac55fb0-2816-11eb-12da-a9ddac1c913c
# ╠═0ed1a02e-2817-11eb-0b92-e5feb0f17c63
# ╠═a3caae20-2817-11eb-09e6-59d6464cef2a
# ╟─aff35570-25cf-11eb-1627-43fe902b9db5
# ╠═a8a18df2-25cf-11eb-38b7-e128e53ba029
# ╠═081bb16e-25d0-11eb-28e6-c3e5a042fce2
# ╟─f77c9930-25d2-11eb-287f-655c63ea721e
# ╠═cf8b68e2-b8ba-11eb-0177-8fa5c7f332bd
# ╟─097834c0-b8bb-11eb-3e52-c1d543ef9fb6
# ╠═9c2a7480-25d2-11eb-0479-751c0dc8b88e
# ╟─abc0604e-2823-11eb-1708-83268e22b7bc
# ╠═2b9f72b0-281b-11eb-1fbc-15b3e67b017d
# ╠═42884660-281a-11eb-015e-29ff1b21dfe0
# ╠═9acf5880-281b-11eb-3aa5-9b6df9ac7501
# ╠═e8ee4712-281b-11eb-0fbf-7783935a0075
# ╠═0e04e810-281c-11eb-3e72-c58bf8d7ce0d
# ╠═1e143c10-281c-11eb-1117-29f72834203e
# ╠═756a3570-2833-11eb-2596-c37768e5a75b
# ╠═be9d11e0-2833-11eb-1e9f-d3bfc428c893
# ╠═8bdbadb0-281b-11eb-1786-3bef9aee1d39
# ╠═28d56f3e-b8c0-11eb-1329-21f533ccddd9

ext: jl

max_stars_repo_path: Tema 5.jl

max_stars_repo_name: Dictino/CNL

max_issues_repo_path: Tema 5.jl

max_issues_repo_name: Dictino/CNL

max_forks_repo_path: Tema 5.jl

max_forks_repo_name: Dictino/CNL


text:
α = U1(rand() * 2pi)
β = U1(rand() * 2pi)
a = SU2(α, β)
b = SU2(U1(rand() * 2pi), U1(rand() * 2pi))
M = Matrix(a)
N = [α -conj(β); β conj(α)]

@test a == SU2(M)
@test size(Matrix(a)) == (2, 2)
@test a + b isa Array{Complex,2}
@test a + M isa Array{Complex,2}
@test M + a isa Array{Complex,2}
@test a + N isa Array{Complex,2}
@test N + a isa Array{Complex,2}
@test a * b isa SU2
@test a * M isa SU2
@test M * a isa SU2
@test a * N isa Array{Complex,2}
@test N * a isa Array{Complex,2}
@test a .* α isa SU2
@test α .* a isa SU2


text:
# ## Pie

#md # ![FILE_NAME.png](images/FILE_NAME.png)

## ref https://docs.juliaplots.org/latest/generated/pyplot/
x = ["Nerds","Hackers","Scientists"]
y = [0.4,0.35,0.25]
pie(x, y, title="The Julia Community",l = 0.5)


text:
rng = StableRNG(614)

# convert a Binary vector into vector of +1 or -1 values
# (for testing only):
pm1(y) = Int8(2) .* (Int8.(MLJBase.int(y))) .- Int8(3)

const MARGIN_LOSSES = MLJBase.MARGIN_LOSSES
const DISTANCE_LOSSES = MLJBase.DISTANCE_LOSSES

# using `WeightedSum` instead of `WeightedMean`; see
# https://github.com/JuliaML/LossFunctions.jl/issues/149
WeightedSum(w) = LossFunctions.AggMode.WeightedMean(w, normalize=false)

@testset "naked" begin
    @test MLJBase.naked(MLJBase.LossFunctions.PeriodicLoss{Float64}) ==
        :PeriodicLoss
    @test MLJBase.naked(MLJBase.LossFunctions.PeriodicLoss) ==
        :PeriodicLoss
end

@testset "LossFunctions.jl - binary" begin
    y = categorical(["yes", "yes", "no", "yes"])
    yes, no = y[1], y[3]
    dyes = MLJBase.UnivariateFinite([yes, no], [0.6, 0.4])
    dno =  MLJBase.UnivariateFinite([yes, no], [0.3, 0.7])
    yhat = [dno, dno, dyes, dyes]
    w = [1, 2, 3, 4]

    @test MLJBase.ZeroOneLoss()(yhat, y) ≈ [1, 1, 1, 0]
    @test MLJBase.zero_one_loss(yhat,y, w) ≈ [1, 2, 3, 0]

    N = 10
    y = categorical(rand(rng, ["yes", "no"], N), ordered=true)
    levels!(y, ["no", "yes"])
    no, yes = MLJBase.classes(y[1])
    @test pm1([yes, no]) in [[+1, -1], [-1, +1]]
    ym = pm1(y) # observations for raw LossFunctions measure
    p_vec = rand(N)
    yhat = MLJBase.UnivariateFinite([no, yes], p_vec, augment=true)
    yhatm = MLJBase._scale.(p_vec) # predictions for raw LossFunctions measure
    w = rand(rng, N)

    for M_ex in MARGIN_LOSSES
        m = eval(:(MLJBase.$M_ex()))
        @test m(yhat, y) ≈ LossFunctions.value(getfield(m, :loss), ym, yhatm)
        @test MLJBase.Mean()(m(yhat, y, w)) ≈
            LossFunctions.value(getfield(m, :loss),
                                ym,
                                yhatm,
                                WeightedSum(w))/N
    end
end

@testset "LossFunctions.jl - continuous" begin
    # losses for continuous targets:
    N    = 10
    y    = randn(rng, N)
    yhat = randn(rng, N)
    X    = nothing
    w    = rand(rng, N)

    for M_ex in DISTANCE_LOSSES
        m = eval(:(MLJBase.$M_ex()))
        m_ex = MLJBase.snakecase(M_ex)
        @test m == eval(:(MLJBase.$m_ex))
        @test m(yhat, y) ≈
            LossFunctions.value(getfield(m, :loss), y, yhat)
        @test mean(m(yhat ,y, w)) ≈
            LossFunctions.value(getfield(m, :loss), y, yhat,
                                WeightedSum(w))/N

    end
end


text:
using OmegaCore
using Distributions
using Test

@testset "Pointwie" begin
  x = 1 ~ Normal(0, 1)
  y = abs(x)
  y_, x_ = randsample((y, x))
  @test y_ == abs(x_)
  @test isinferred(randsample, x)
end

@testset "Pointwise" begin
  x = 1 ~ Normal(0, 1)
  y = abs(x)
  y_, x_ = randsample((y, x))
  @test y_ == abs(x_)
  @test isinferred(randsample, x)
end


text:
# compare output from KernelDensity.jl with SINC kernel
# to what Comte-Butucea yields. (should be the same)
# also look at how DeLaValleePoussin performs

using EmpiricalBayes
using Test
using Distributions
using Random

true_dist = MixtureModel([ Normal(-0.3,.5), Normal(1.05,.5)])
marginal_grid = collect(range(-6, stop=6, length=1001));

marginal_h = marginal_grid[2] - marginal_grid[1] # should really replace by range

prior_grid = collect(range(-3, stop=3, length=121));
d_true = NormalConvolutionProblem(true_dist, marginal_grid);

ds = MixingNormalConvolutionProblem(Normal, 0.2, prior_grid, marginal_grid);

m = 10_000

Random.seed!(2)
Xs = rand(d_true, m)

#----------------------------------------
#---- Test with actual Sinc Kernel-------
#----------------------------------------
f_sinc = sinc_kde(Xs, marginal_grid, SincKernel)

x_tst_1 = f_sinc.x[700]
f_tst_1 = f_sinc.density[700]

# Let us compare to the result from comte_butucea
f_cbt_1 = estimate(Xs, ComteButucea,
  MarginalDensityTarget(x_tst_1), marginal_grid)

@test f_cbt_1 ≈ f_tst_1 atol=0.001

#----- Check default dispatch ---------------------------------
dv_kde = sinc_kde(Xs, marginal_grid, DeLaValleePoussinKernel)
dv_kde_auto_dispatch = sinc_kde(Xs, marginal_grid)
@test dv_kde.density == dv_kde_auto_dispatch.density

#--------------------------------------------------------------
#---- Tests with both Sinc and DeLaValleePoussin Kernels-------
#--- Are our bands longer than true bands?
#--------------------------------------------------------------

f_true = pdf.(Ref(d_true), marginal_grid)

for T in [SincKernel, DeLaValleePoussinKernel]
  f_kde= sinc_kde(Xs, marginal_grid, T)
  true_C = maximum(abs.(f_true .- f_kde.density))

  f_nb = fit(BinnedMarginalDensityNeighborhood, Xs, marginal_grid, T)
  f_nb.C_std

  f_nb_ds = fit(BinnedMarginalDensityNeighborhood, Xs, ds, T)
  @show true_C, f_nb_ds.C_std, f_nb_ds.C_bias, f_nb_ds.C_std + f_nb_ds.C_bias
  @test true_C <= f_nb_ds.C_std + f_nb_ds.C_bias
end

# What are these doing here?
Xs_test = rand(d_true, m)
full_CI = CEB_ci(Xs, Xs_test, ds, PosteriorTarget(PosteriorMeanNumerator(2.0)))

full_CI[1].ci_left, full_CI[1].ci_right

# Use output from first run to speed up 2nd run by factor 4!
full_CI_fast_comp = CEB_ci(Xs, Xs_test, ds, PosteriorTarget(PosteriorMeanNumerator(2.0)),
                            full_CI[2], full_CI[3], full_CI[4].C_bias)

full_CI_fast_comp[1].ci_left, full_CI_fast_comp[1].ci_right

@test full_CI[1].ci_left ≈ full_CI_fast_comp[1].ci_left atol=0.005
@test full_CI[1].ci_right ≈ full_CI_fast_comp[1].ci_right atol=0.005

full_CI_inf = CEB_ci(Xs, Xs_test, ds, PosteriorTarget(PosteriorMeanNumerator(2.0)); C=Inf)

full_CI_inf[1].ci_left, full_CI_inf[1].ci_right

posterior_stats(d_true,  PosteriorTarget(PosteriorMeanNumerator(2.0)))


text:
function round2sonum8(x::Float64)
    return signbit(x) ? -sonum8[findFirstSmaller(-x,bounds8)-2)] : sonum8[findFirstSmaller(x,bounds8)-2)]
end

function round2sonum16(x::Float64)
    return signbit(x) ? -sonum16[findFirstSmaller(-x,bounds16)-2)] : sonum16[findFirstSmaller(-x,bounds16)-2)]
end


text:
abstract type PMUType end

pmutype(::T) where {T} = error("`pmutype` not defined for arguments of type $T")

# Defaults
_unitstatus(x::PMUType, i...) = error("`unitstatus` undefined for $(typeof(x))")
_unitcontrol(x::PMUType, i...) = error("`unitcontrol` undefined for $(typeof(x))")
_counter(x::PMUType, i...) = error("`counter` undefined for $(typeof(x))")
_control(x::PMUType, i...) = error("`control` undefined for $(typeof(x))")
_extras(x::PMUType, i...) = error("`extras` undefined for $(typeof(x))")

writetype(::PMUType) = UInt32

numcounters(x::PMUType) = error("`numcounters` undefined for $(typeof(x))")
numbytes(x::PMUType) = sizeof(UInt64) * numcounters(x)

unpack(_) = ()
unitstatus(x, i...) = _unitstatus(pmutype(x), indexzero.((unpack(x)..., i...))...)
unitcontrol(x, i...) = _unitcontrol(pmutype(x), indexzero.((unpack(x)..., i...))...)
counter(x, i...) = _counter(pmutype(x), indexzero.((unpack(x)..., i...))...)
control(x, i...) = _control(pmutype(x), indexzero.((unpack(x)..., i...))...)
extras(x, i...) = _extras(pmutype(x), indexzero.((unpack(x)..., i...))...)
writetype(x) = writetype(pmutype(x))

numcounters(x) = numcounters(pmutype(x))
numbytes(x) = numbytes(pmutype(x))

### Integrated Memory Controller
struct IMC{T<:AbstractCPU} <: PMUType end
_unitstatus(::IMC{SkylakeServer}) = IndexZero(0xF8)
_unitcontrol(::IMC{SkylakeServer}) = IndexZero(0xF4)
_counter(::IMC{SkylakeServer}, i::IndexZero) = IndexZero(0xA0 + value(i) * 0x8)
_control(::IMC{SkylakeServer}, i::IndexZero) = IndexZero(0xD8 + value(i) * 0x4)
numcounters(::IMC) = 4

# For now, only read the fixed counters for IceLake servers.
# There are 4 such counters, starting at address 0x2290 and they are
# DRAM Read, DRAM Write, PM Read, and PM Write respectively
_counter(::IMC{IcelakeServer}, i) = IndexZero(0x2290 + value(i) * 0x8)

### CHA Counters
struct CHA <: PMUType end
_unitstatus(::CHA, i) = IndexZero(0xE07 + value(i) * 0x10)
_unitcontrol(::CHA, i) = IndexZero(0xE00 + value(i) * 0x10)
_counter(::CHA, cha, i) = IndexZero(0xE08 + value(cha) * 0x10 + value(i))
_control(::CHA, cha, i) = IndexZero(0xE01 + value(cha) * 0x10 + value(i))
_extras(::CHA, cha, i) = IndexZero(0xE05 + value(cha) * 0x10 + value(i))
writetype(::CHA) = UInt64
numcounters(::CHA) = 4

# Customize for various types Specialize
abstract type AbstractUncorePMU end

##### IMC Uncore PMU
# PMU implementation for monitoring the integrated memory controller
struct IMCUncorePMU <: AbstractUncorePMU
    # A handle to the underlying
    handle::Handle
end
unwrap(x::IMCUncorePMU) = x.handle
pmutype(::IMCUncorePMU) = IMC{SkylakeServer}()
Base.close(x::IMCUncorePMU) = close(x.handle)

# IceLake IMC PMU
# For now - only return the free-running counters
struct IMCUncoreICX <: AbstractUncorePMU
    mmio::MMIO
end
unwrap(x::IMCUncoreICX) = x.mmio
pmutype(::IMCUncoreICX) = IMC{IcelakeServer}()
Base.close(::IMCUncoreICX) = nothing

##### CHA Uncore PMU
# PMU implementation for monitoring the CHA
struct CHAUncorePMU <: AbstractUncorePMU
    # We hold on to a single handle for the MSR path, shared by all PMUs
    handle::Handle
    # The number of this CHA
    cha::IndexZero{Int}
    buffer::Vector{UInt8}

    # Allow passing a buffer, or manually create one
    function CHAUncorePMU(handle::Handle, cha, buffer = zeros(UInt8, numbytes(CHA())))
        resize!(buffer, numbytes(CHA()))

        return new(handle, indexzero(cha), buffer)
    end
end

unwrap(x::CHAUncorePMU) = x.handle
pmutype(::CHAUncorePMU) = CHA()
unpack(x::CHAUncorePMU) = (x.cha,)
Base.close(x::CHAUncorePMU) = close(x.handle)

#####
##### Low level accessing functions
#####

function setunitstatus!(U::AbstractUncorePMU, v)
    write(unwrap(U), convert(writetype(U), v), unitstatus(U))
end

function getunitstatus(U::AbstractUncorePMU)
    return read(unwrap(U), UInt32, unitstatus(U))
end

function setunitcontrol!(U::AbstractUncorePMU, v)
    write(unwrap(U), convert(writetype(U), v), unitcontrol(U))
end

function getunitcontrol(U::AbstractUncorePMU)
    return read(unwrap(U), UInt32, unitcontrol(U))
end

function setcontrol!(U::AbstractUncorePMU, counter, v)
    return write(unwrap(U), convert(writetype(U), v), control(U, counter))
end

function getcontrol(U::AbstractUncorePMU, i)
    return read(unwrap(U), UInt32, control(U, i))
end

function getcounter(U::AbstractUncorePMU, i)
    return CounterValue(read(unwrap(U), UInt64, counter(U, i)))
end

function setextra!(U::AbstractUncorePMU, i, v)
    write(unwrap(U), convert(writetype(U), v), extras(U, i))
end

function getextra(U::AbstractUncorePMU, i)
    return read(unwrap(U), UInt32, extras(U, i))
end

#####
##### Some higher level functions
#####

function getallcounters(U::AbstractUncorePMU)
    return ntuple(i -> getcounter(U, i), Val(numcounters(U)))
end

function reset!(U::AbstractUncorePMU)
    # Write to the unit control to clear all counters and control registers
    val = setbits(zero(writetype(U)), (0, 1, 8, 16, 17))
    setunitcontrol!(U, val)
end

function enable!(U::AbstractUncorePMU)
    val = setbits(zero(writetype(U)), (16, 17))
    setunitcontrol!(U, val)
end


text:
__precompile__(true)
module PkgMirrors

include("types.jl")
const PkgMirror = Types.PkgMirror
export PkgMirror

include("utils.jl")
const current = Utils.current
export current

const MIRRORS = Dict{String, String}()

export PKG
include("pkg.jl")

include("api.jl")

const availables = API.availables
const setmirror = API.setmirror
const activate = API.activate
const deactivate = API.deactivate
const clear = API.clear
const activated = API.activated
export availables, setmirror, activate, deactivate, clear, activated

function __init__()
    open(joinpath(@__DIR__, "../data/mirror_list.txt")) do fi
        for line = readlines(fi)
            tmp = split(line, ' ')
            MIRRORS[tmp[1]] = tmp[2]
        end
    end
    API.activate(first=false)
    finalizer(Types.FINALIZER) do _
        activated() && API.deactivate()
    end
end

end


text:
"""
    AbstractNLPEvaluator

AbstractNLPEvaluator implements the bridge between the
problem formulation (see `AbstractFormulation`) and the optimization
solver. Once the problem formulation bridged, the evaluator allows
to evaluate in a straightfoward fashion the objective and the different
constraints, but also the corresponding gradient and Jacobian objects.

"""
abstract type AbstractNLPEvaluator end

abstract type AbstractNLPAttribute end

"""
    Variables <: AbstractNLPAttribute end

Attribute corresponding to the optimization variables attached
to a given `AbstractNLPEvaluator`.
"""
struct Variables <: AbstractNLPAttribute end

"""
    Constraints <: AbstractNLPAttribute end

Attribute corresponding to the constraints  attached
to a given `AbstractNLPEvaluator`.
"""
struct Constraints <: AbstractNLPAttribute end

"""
    n_variables(nlp::AbstractNLPEvaluator)
Get the number of variables in the problem.
"""
function n_variables end

"""
    n_constraints(nlp::AbstractNLPEvaluator)
Get the number of constraints in the problem.
"""
function n_constraints end

# Callbacks
"""
    objective(nlp::AbstractNLPEvaluator, u)::Float64

Evaluate the objective at point `u`.
"""
function objective end

"""
    gradient!(nlp::AbstractNLPEvaluator, g, u)

Evaluate the gradient of the objective at point `u`. Store
the result inplace in the vector `g`, which should have the same
dimension as `u`.

"""
function gradient! end

"""
    constraint!(nlp::AbstractNLPEvaluator, cons, u)

Evaluate the constraints of the problem at point `u`. Store
the result inplace, in the vector `cons`.

## Note
The vector `cons` should have the same dimension as the result
returned by `n_constraints(nlp)`.

"""
function constraint! end

"""
    jacobian_structure!(nlp::AbstractNLPEvaluator, rows, cols)

Return the sparsity pattern of the Jacobian matrix.
"""
function jacobian_structure! end

"""
    jacobian!(nlp::ReducedSpaceEvaluator, jac, u)

Evaluate the Jacobian of the constraints at position `u`. Store
the result inplace, in the array `jac`.
"""
function jacobian! end

"""
    jtprod!(nlp::ReducedSpaceEvaluator, jv, u, v)

Evaluate the transpose Jacobian-vector product ``J^{T} v``.
The vector `jv` is modified inplace.

Let `(n, m) = n_variables(nlp), n_constraints(nlp)`.

* `u` is a vector with dimension `n`
* `v` is a vector with dimension `m`
* `jv` is a vector with dimension `n`

"""
function jtprod! end

"""
    hessian!(nlp::AbstractNLPEvaluator, hess, u)

Evaluate the Hessian of the problem at point `u`. Store
the result inplace, in the vector `hess`.

"""
function hessian! end

# Utilities
"""
    primal_infeasibility(nlp::AbstractNLPEvaluator, u)

Return primal infeasibility associated to current model `nlp` evaluated
at control `u`.

"""
function primal_infeasibility end

"""
    primal_infeasibility!(nlp::AbstractNLPEvaluator, cons, u)

Return primal infeasibility associated to current model `nlp` evaluated
at control `u`. Modify vector `cons` inplace.

"""
function primal_infeasibility! end

"""
    reset!(nlp::AbstractNLPEvaluator)

Reset evaluator `nlp` to default configuration.

"""
function reset! end

include("common.jl")
include("reduced_evaluator.jl")
include("penalty.jl")
include("auglag.jl")
include("MOI_wrapper.jl")


text:
##### ArbUInt

export ArbUInt

struct ArbUInt <: Unsigned
    data::Vector{ArbDigit}
    ArbUInt(data::AbstractVector{ArbDigit}) = normalize!(new(data))
end
ArbUInt(data::AbstractVector) = ArbUInt(convert(Vector{ArbDigit}, data))
ArbUInt() = zero(ArbUInt)
ArbUInt(x) = ArbUInt(ArbDigit[x])
ArbUInt(a::ArbUInt) = deepcopy(a)
ArbUInt(x::DoubleArbDigit) = ArbUInt(reverse!([fromDoubleArbDigit(x)...]))
if DoubleArbDigit === UInt64
    ArbUInt(x::UInt128) = ArbUInt(reverse!([u32_from_u128(x)...]))
end

Base.promote_type(::Type{<:Integer}, ::Type{ArbUInt}) = ArbUInt

Base.show(io::IO, m::MIME"text/plain", a::ArbUInt) = print(io, "0x", string(a; base=16))

Base.deepcopy(a::ArbUInt) = ArbUInt(deepcopy(a.data))

Base.:(==)(a::ArbUInt, b::ArbUInt) = length(a.data) == length(b.data) && a.data == b.data
Base.hash(a::ArbUInt, h::UInt) = hash(a.data, h)
Base.:(<=)(a::ArbUInt, b::ArbUInt) = a == b || a < b

Base.:(<)(a::ArbUInt, b::ArbUInt) = return _less(<, a.data, b.data)
function _less(op, a, b)
    length(a) < length(b) && return true
    length(a) > length(b) && return false
    idx = lastindex(a) # we only hit this method when length(a) === length(b)
    while idx > 0 && a[idx] == b[idx]
        idx = prevind(a, idx)
    end
    idx != 0 && op(a[idx], b[idx])
end

Base.zero(::Type{ArbUInt}) = ArbUInt(ArbDigit[])
Base.one(::Type{ArbUInt}) = ArbUInt([one(ArbDigit)])
Base.isone(a::ArbUInt) = isone(length(a.data)) && isone(a.data[end])
Base.iszero(a::ArbUInt) = isempty(a.data)

function normalize!(a::ArbUInt)
    isempty(a.data) && return a
    idx = lastindex(a.data)
    while checkbounds(Bool, a.data, idx) && iszero(a.data[idx])
        idx = prevind(a.data, idx)
    end
    resize!(a.data, idx)
    a
end

function set_one!(a::ArbUInt)
    resize!(a.data, 1)
    a.data[begin] = one(ArbDigit)
    a
end

function set_zero!(a::ArbUInt)
    resize!(a.data, 0)
    a
end

function Base.trailing_zeros(a::ArbUInt)
    i = findfirst(!iszero, a.data)
    i === nothing && return length(a.data) * BITS
    zeros = trailing_zeros(a.data[i])
    (i-1) * BITS + zeros
end

function Base.trailing_ones(a::ArbUInt)
    i = findfirst(!iszero, a.data)
    i === nothing && return 0
    ones = trailing_ones(a.data[i])
    (i-1) * BITS + ones
end

Base.leading_ones(a::ArbUInt) = leading_ones(first(a.data))
function Base.leading_zeros(a::ArbUInt)
    i = findlast(!iszero, a.data)
    i === nothing && return length(a.data) * BITS
    zeros = leading_zeros(a.data[i])
    (length(a.data) - i) * BITS + zeros
end

Base.count_ones(a::ArbUInt) = isempty(a.data) ? 0 : mapreduce(count_ones, +, a.data)
function Base.count_zeros(a::ArbUInt)
    iszero(a) || isone(a) && return 0
    nzeros = mapreduce(count_zeros, +, a.data)
    nzeros - leading_zeros(a)
end

function is_bit_set(a::ArbUInt, bit)
    bits_per_digit = UInt(BITS)
    digit_index = bit ÷ bits_per_digit
    iszero(digit_index) && return false
    digit = a.data[digit_index]
    mask = one(ArbDigit) << (bit % bits_per_digit)
    return (digit & mask) != 0
end

function set_bit!(a::ArbUInt, bit, value)
    bits_per_digit = UInt(BITS)
    digit_index = bits_per_digit > bit ? typemax(UInt) : bit ÷ bits_per_digit
    mask = one(ArbDigit) << (bit % bits_per_digit)
    if value
        if digit_index > length(a.data)
            resize!(a.data, digit_index)
        end
        a.data[digit_index] |= mask
    elseif digit_index <= length(a.data)
        a.data[digit_index] &= ~mask
        normalize!(a.data)
    end
    a
end

include("rng.jl")
include("strings.jl")
include("addition.jl")
include("subtraction.jl")
include("bits.jl")
include("shifts.jl")
include("multiplication.jl")

### SUBTRACTION ###

# TODO: Implement signed numbers

### DIVISION ###

# TODO: Implement division

### POWER ###

# TODO: Implement power function

### MODULO ###

# TODO: Implement modulo arithmetic


text:
using BayesMultiFit
BAMF = BayesMultiFit

"""
SLIVER1DXMeasType(images::Int32, inttime::Float32, invx::Float32, invy::Float32)

Contains the information necessary to create a SLIVER measurement, the integration time, inttime, the xy coordinates of the inversion point, invx and invy, and the static number of images a SLIVER measurement creates, images.

THIS IS AN EXAMPLE OF AN ADAPTIVE MEASUREMENT TYPE
"""
struct SLIVER1DXMeasType <: BAMF.MeasType
    images::Int32
    inttime::Float32
    invx::Float32
end

function SLIVER1DXMeasType(inttime::Float32, invx::Float32)
    return SLIVER1DXMeasType(Int32(2), inttime, invx)
end

function SLIVER1DXMeasType(info::Tuple, inttime::Float32)
    invx, = info
    return SLIVER1DXMeasType(inttime, invx)
end

function BAMF.pickMeasType(::Val{Int32(3)}, invx::Float32, invy::Float32)
    return (SLIVER1DXMeasType, (invx,))
end

function BAMF.genimage(s::BAMF.BAMFState, sz::Int32, psf::BAMF.PSF, meas::SLIVER1DXMeasType)
    image = zeros(Float32, sz, sz, 2)
    cntr = Int32(ceil((sz + 1) / 2))
    for n_emit in 1:s.n, ii in 1:sz, jj in 1:sz
        r1 = sqrt((s.y[n_emit] - ii)^2 + (s.x[n_emit] - meas.invx - jj + cntr)^2)
        a1 = sqrt(1 / 2f0 * meas.inttime * s.photons[n_emit]) * BAMF.spread(psf, r1)
        r2 = sqrt((-s.y[n_emit] + ii)^2 + (-s.x[n_emit] + meas.invx - jj + cntr)^2)
        a2 = sqrt(1 / 2f0 * meas.inttime * s.photons[n_emit]) * BAMF.spread(psf, r2)
        image[ii, jj, 1] += 0.5f0*abs(a1 + a2)^2
        image[ii, jj, 2] += 0.5f0*abs(a1 - a2)^2
    end
    return image
end


text:
using Pkg
Pkg.activate(pwd())

# # Variables

x = 2
typeof(x)

y = x + 1
typeof(y)

x = 4
x = 3.1415

Float64

# ### Exercise:
# Create the following three variables:
# 1. Variable `x` with value `1.234`.
# 2. Variable `y` with value `1//2`.
# 3. Variable `z` with value `x + y*im`.
# What are the types of these three variables?
#
# ---
# ### Solution:

# ---
#
# ## Variable Names

I_am_float = 3.1415
CALL_ME_RATIONAL = 1//3
MyString = "MyVariable"

#+

α = 1
δ = 1
🍕 = "It's time for pizza!!!"

#+

π = 2
π

#+

ℯ
ℯ = 2

#+

struct = 3


text:
using IRTracker
using DynamicPPL
using Distributions
using DataStructures: SortedDict

export strip_dependencies, makegraph, showgraph

try_getindex(expr::TapeReference) = expr[]
try_getindex(expr::TapeConstant) = getvalue(expr)

const CallingNode{F} = Union{PrimitiveCallNode{<:Any, F}, NestedCallNode{<:Any, F}}

istilde(::CallingNode{typeof(DynamicPPL.tilde_assume)}) = true
istilde(::CallingNode{typeof(DynamicPPL.tilde_observe)}) = true
istilde(::CallingNode{typeof(DynamicPPL.dot_tilde_assume)}) = true
istilde(::CallingNode{typeof(DynamicPPL.dot_tilde_observe)}) = true
istilde(::AbstractNode) = false

"""
    tilde_parameters(node)

Extract a tuple of variablen name, distribution, and value argument of any tilde call.

The variable name can be `nothing` for constant observations.  The whole value is `nothing` for
nodes that are not tilde calls.
"""
tilde_parameters(node::AbstractNode) = nothing

function tilde_parameters(node::CallingNode{typeof(DynamicPPL.tilde_assume)})
    args = getarguments(node) # rng, ctx, sampler, right, vn, inds, vi
    vn, dist, value, vi = args[5], args[4], getvalue(node), args[7]
    return vn, dist, TapeConstant(value), vi
end

function tilde_parameters(node::CallingNode{typeof(DynamicPPL.tilde_observe)})
    args = getarguments(node)
    if length(args) == 7
        # ctx, sampler, right, left, vname, vinds, vi
        vn, dist, value, vi = args[5], args[3], args[4], args[7]
        return vn, dist, value, vi
    elseif length(args) == 5
        # ctx, sampler, right, left, vi
        dist, value, vi = args[3], args[4], args[5]
        return nothing, dist, value, vi
    else
        throw(ArgumentError("$node has unknown argument structure, you should not have reached this point!"))
    end
end

function tilde_parameters(node::CallingNode{typeof(DynamicPPL.dot_tilde_assume)})
    args = getarguments(node) # rng, ctx, sampler, right, left, vn, inds, vi
    vn, dist, value, vi = args[6], args[4], args[5], args[8]
    return vn, dist, value, vi
end

function tilde_parameters(node::CallingNode{typeof(DynamicPPL.dot_tilde_observe)})
    args = getarguments(node)
    if length(args) == 7
        # ctx, sampler, right, left, vn, inds, vi
        vn, dist, value, vi = args[5], args[3], args[4], args[7]
        return vn, dist, value, vi
    elseif length(args) == 5
        # ctx, sampler, right, left, vi
        dist, value, vi = args[3], args[4], args[5]
        return nothing, dist, value, vi
    else
        throw(ArgumentError("$node has unknown argument structure, you should not have reached this point!"))
    end
end

"""
    model_argument_nodes(root)

Extract from a stripped model trace all model argument `getindex` calls.
"""
function model_argument_nodes(root)
    # From the beginning of the trace,
    # ```
    # @1: [Arg:§1:%1] @11#1 → ##evaluator#465
    # @2: [Arg:§1:%2] @11#2 → Random._GLOBAL_RNG()
    # @3: [Arg:§1:%3] @11#3 → Model{...}
    # @4: [Arg:§1:%4] @11#4 → VarInfo (0 variables, dimension 0; logp: 0.0)
    # @5: [Arg:§1:%5] @11#5 → SampleFromPrior()
    # @6: [Arg:§1:%6] @11#6 → DefaultContext()
    # @7: [§1:%7] ⟨getproperty⟩(@3, ⟨:args⟩) → (x = [0.1, 0.05, 1.0],)
    # @8: [§1:%8] ⟨getproperty⟩(@7, ⟨:x⟩) → [0.1, 0.05, 1.0]
    # ```
    # extract only the `getproperty(@n, ⟨:varname⟩)` line for each argument.

    argument_nodes = Vector{AbstractNode}()
    model_node = getchildren(root)[3]
    modelargs_nodes = Vector{AbstractNode}()

    for child in getchildren(root)
        if child isa CallingNode{typeof(getproperty)}
            value, property_name = getarguments(child)
            if getvalue(property_name) == :args && try_getindex(value) == model_node
                push!(modelargs_nodes, child)
            elseif try_getindex(value) ∈ modelargs_nodes
                # a `getproperty(args, symbol)` whose parent is a `args = getproperty(model, :args)`
                # is a model argument, for sure
                push!(argument_nodes, child)
            end
        end
    end

    return argument_nodes
end

"""
    add_candidates!(dependencies, candidates, mutants, node)

Traverse refernces of `node` backwards and turn the reached `candidates` into `dependencies`.
"""
function add_candidates!(dependencies, candidates, mutants, node)
    current_refs = Vector{AbstractNode}(referenced(node))
    push!(dependencies, node)

    while !isempty(current_refs)
        node = pop!(current_refs)
        node ∈ dependencies && continue

        if node ∈ candidates
            push!(dependencies, node)
            union!(current_refs, referenced(node))
        end

        if haskey(mutants, node)
            push!(dependencies, node)
            union!(current_refs, mutants[node])
        end
    end
end

"""
    strip_dependencies(root)

Slice out only those nodes that are between tildes or model arguments, taking some care to get
mutated things (vector + `setindex!` etc.) right.
"""
function strip_dependencies(root)
    # "dependencies" are nodes that are in the dependency graph for sure (i.e., tildes).
    # "candidates" are nodes that follow a dependency node; are included as dependencies, if the
    # chain ends up in a dependency node again (e.g., x ~ D, y = f(x), z ~ D2(y)) -- this is
    # what `add_candidates!` does by working backwards.
    # "mutants" are backedges (from mutated to mutators) in the dependency graph, resulting from
    # later operations mutating earlier values(e.g, z = zeros(), z[1] ~ D, x ~ D2(z[1])).

    dependencies = model_argument_nodes(root)
    candidates = Set{AbstractNode}()
    mutants = Dict{AbstractNode, Vector{AbstractNode}}()

    for child in getchildren(root)
        params = tilde_parameters(child)

        if !isnothing(params)
            # tilde node: is a dependency for sure
            push!(dependencies, child)
        else
            # non-tilde node: make candidate, if any parent is a candidate or dependency

            # this `t` is critical -- exploding type inference!
            t = any(r ∈ candidates || r ∈ dependencies for r in referenced(child))
            if t
                push!(candidates, child)

                if child isa CallingNode{typeof(setindex!)}
                    # if the candidate is mutating an array, also make the mutated object a candidate
                    # and record the backedge from the mutated one to:
                    # `n = setindex!(x, v, i)` makes `x` a candidate, and a backedge from `x` to `n`.

                    mutated = getargument(child, 1)
                    if mutated isa TapeReference
                        push!(candidates, mutated[])
                        push!(get!(Vector{AbstractNode}, mutants, mutated[]), child)
                    end
                end
            end
        end
    end

    # Go over nodes once more and turn candidates between dependencies into dependencies, too.
    for child in getchildren(root)
        params = tilde_parameters(child)
        if !isnothing(params)
            vn, dist, value = params
            vn isa TapeReference && add_candidates!(dependencies, candidates, mutants, vn[])
            dist isa TapeReference && add_candidates!(dependencies, candidates, mutants, dist[])
            value isa TapeReference && add_candidates!(dependencies, candidates, mutants, value[])
        end
    end

    return sort(dependencies; by=getposition) |> unique
end

struct Reference
    number::Int
end

Base.show(io::IO, r::Reference) = print(io, "⟨", r.number, "⟩")
Base.isless(q::Reference, r::Reference) = isless(q.number, r.number)
Base.hash(r::Reference, h::UInt) = hash(r.number, h)
Base.:(==)(q::Reference, r::Reference) = q.number == r.number

abstract type Statement end

struct Assumption{dotted, TN<:VarName, TVal} <: Statement
    vn::TN
    dist_ref::Reference
    value::TVal
    logp::Float64
end

Assumption{dotted}(vn, dist_ref, value, logp) where {dotted} =
    Assumption{dotted, typeof(vn), typeof(value)}(vn, dist_ref, value, logp)

struct Observation{dotted, TN<:Union{Nothing, VarName}, TVal} <: Statement
    vn::TN
    dist_ref::Reference
    value::TVal
    logp::Float64
end

Observation{dotted}(vn, dist_ref, value, logp) where {dotted} =
    Observation{dotted, typeof(vn), typeof(value)}(vn, dist_ref, value, logp)

struct Call{TD<:Union{Nothing, Tuple{VarName, Reference}}, TF, TArgs<:Tuple, TVal} <: Statement
    definition::TD
    f::TF
    args::TArgs
    value::TVal
end

Call(f, args, value) = Call(nothing, f, args, value)

struct Constant{TVal} <: Statement
    value::TVal
end

const Tilde{dotted} = Union{Assumption{dotted}, Observation{dotted}}

_shortname(d::Type{<:Distribution}) = string(nameof(d))
_shortname(other) = string(other)

function Base.show(io::IO, stmt::Assumption)
    print(io, stmt.vn, (isdotted(stmt) ? " .~ " : " ~ "), stmt.dist_ref, " → ", stmt.value)
end
function Base.show(io::IO, stmt::Observation)
    print(io, (isnothing(stmt.vn) ? stmt.value : stmt.vn))
    print(io, (isdotted(stmt) ? " .⩪ " : " ⩪ "), stmt.dist_ref,  " ← ", stmt.value)
end
function Base.show(io::IO, stmt::Call)
    if !isnothing(stmt.definition)
        vn = stmt.definition[1]
        print(io, vn, " = ")
    end
    print(io, stmt.f)
    print(io, "(")
    join(io, stmt.args, ", ")
    print(io, ") → ", stmt.value)
end
Base.show(io::IO, stmt::Constant) = print(io, stmt.value)

IRTracker.getvalue(stmt::Statement) = stmt.value

isdotted(::Type{<:Assumption{dotted}}) where {dotted} = dotted
isdotted(::Type{<:Observation{dotted}}) where {dotted} = dotted
isdotted(tilde::Tilde) = isdotted(typeof(tilde))

getvn(stmt::Tilde) = stmt.vn
getvn(stmt::Call) = isnothing(stmt.definition) ? nothing : stmt.definition[1]
getvn(::Constant) = nothing

"""
    parent_variables(graph, stmt)

Return all `Assumption`s (potentially with indexing) that the tilde `stmt` depends on directly.
"""
function parent_variables(graph, stmt::Tilde)
    result = Set{Pair{VarName, Assumption}}()
    parent_variables!(result, graph, graph[stmt.dist_ref])
    stmt.value isa Reference && parent_variables!(result, graph, stmt.value)
    return result
end

parent_variables!(result, graph, stmt) = result
function parent_variables!(result, graph, ref::Reference)
    # follow back implicit dependencies on unmarked locations without indexing, e.g.
    # ⟨2⟩ = [0.1, -0.2]
    # ⟨3⟩ = Array{Float64,1}(array initializer with undefined values, 2) → ...
    # ⟨6⟩ = m[1] ~ Normal() → ...
    # ⟨7⟩ = m[1] = getindex(⟨3⟩, 1) → ...
    # ⟨10⟩ = m[2] ~ Normal() → ...
    # ⟨11⟩ = m[2] = getindex(⟨3⟩, 2) → ...
    # ⟨13⟩ = x ⩪ MvNormal(⟨3⟩) ← ⟨2⟩
    # `x` depends on the whole of `m`, stored in `⟨3⟩`, and thus implicitly
    # on `m[1]` and `m[2]` (cf. reverse_deps test case).

    if haskey(graph.rv_locations, ref)
        for (ix, loc) in graph.rv_locations[ref]
            loc == ref && continue # skip one-line array tildes, e.g., z ~ filldist(...)
            parent_variables!(result, graph, loc)
        end
    end

    return parent_variables!(result, graph, graph[ref])
end
function parent_variables!(result, graph, stmt::Assumption)
    # direct dependency:
    # ⟨31⟩ = w ~ Dirichlet()
    # ⟨42⟩ = z[1] ~ DiscreteNonParametric(⟨31⟩)
    return push!(result, stmt.vn => stmt)
end
function parent_variables!(result, graph, stmt::Observation)
    @warn "The parent statement is an observation ($(stmt)), something weird has happened..."
    return result
end
function parent_variables!(result, graph, stmt::Call{<:Nothing})
    # go through all argumens of an unmarked location:
    # ⟨85⟩ = getindex(⟨2⟩, ⟨72⟩)
    # ⟨86⟩ = x[3] ⩪ Normal(⟨80⟩, 1.0) ← ⟨85⟩

    for arg in stmt.args
        parent_variables!(result, graph, arg)
    end

    return result
end
function parent_variables!(result, graph, stmt::Call{<:Tuple, typeof(getindex)})
    # special case: don't add whole array (`z`) for indexing calls
    # ⟨12⟩ = z ~ filldist(⟨10⟩, ⟨5⟩)
    # ⟨30⟩ = z[2] = getindex(⟨12⟩, ⟨28⟩)
    # ⟨31⟩ = μ[2] = getindex(⟨7⟩, ⟨30⟩)
    for arg in Base.tail(stmt.args)
        parent_variables!(result, graph, arg)
    end
    vn, location = stmt.definition
    push!(result, vn => graph[location])
    return result
end

struct Graph
    """Entries in graph: mapping from references to statements."""
    statements::SortedDict{Reference, Statement}

    """Association which nodes in the trace are associated with which references in the graph"""
    reference_mapping::Dict{AbstractNode, Reference}

    """
    Remembers in which references random variables are actually stored; this is either the location
    of the tilde statement (`beta ~ Normal()`, `z ~ filldist(...)`), or the array position it is
    assigned to (`z[i] ~ Normal()`).

    The container is first ordered by reference of storage, then by index for `subsumes` testing.
    """
    rv_locations::Dict{Reference, Vector{Pair{Tuple, Reference}}}
end

Graph() = Graph(SortedDict{Reference, Statement}(),
                Dict{AbstractNode, Reference}(),
                Dict{Reference, Vector{Pair{Tuple, Reference}}}())

Base.IteratorSize(::Type{Graph}) = Base.HasLength()
Base.length(graph::Graph) = length(graph.statements)
Base.IteratorEltype(::Type{Graph}) = Base.HasEltype()
Base.eltype(graph::Graph) = eltype(graph.statements)
Base.getindex(graph::Graph, ref) = graph.statements[ref]
Base.setindex!(graph::Graph, stmt, ref) = graph.statements[ref] = stmt
Base.firstindex(graph::Graph) = firstindex(graph.statements)
Base.lastindex(graph::Graph) = lastindex(graph.statements)
Base.get(graph::Graph, ref, default) = get(graph.statements, ref, default)
Base.haskey(graph::Graph, ref) = haskey(graph.statements, ref)
Base.delete!(graph::Graph, ref) = delete!(graph.statements, ref)
Base.keys(graph::Graph) = keys(graph.statements)
Base.values(graph::Graph) = values(graph.statements)
Base.iterate(graph::Graph) = iterate(graph.statements)
Base.iterate(graph::Graph, state) = iterate(graph.statements, state)

# just intended for debugging
Base.getindex(graph::Graph, i::Int) = graph[Reference(i)]

# to convert an atom to a runtime value: either the value of the referenced node, or a literal
tovalue(graph, ref::Reference) = getvalue(graph[ref])
tovalue(graph, other) = other

function Base.mapreduce(f, op, graph::Graph; init)
    for kv in graph
        init = op(init, f(kv))
    end

    return init
end

function Base.filter(f, graph::Graph; init=Vector{eltype(graph)}())
    for kv in graph
        f(kv) && push!(init, kv)
    end

    return init
end

Base.map(f, graph::Graph; init=Vector{eltype(graph)}()) = mapreduce(f, push!; init=init)

getstatements(graph::Graph) = graph.statements

hasmapping(graph::Graph, node) = haskey(graph.reference_mapping, node)
function getmapping(graph::Graph, node)
    ref = graph.reference_mapping[node]
    return ref
end
function getmapping(graph::Graph, node, default)
    if hasmapping(graph, node)
        getmapping(graph, node)
    else
        return default
    end
end
setmapping!(graph::Graph, (node, ref)::Pair) = graph.reference_mapping[node] = ref

function deletemapping!(graph::Graph, node)
    if hasmapping(graph, node)
        delete!(graph, getmapping(graph, node))
    end
    return graph
end

"""Look up the matching location of an RV in graph.rv_locations[ref] for given index tuple."""
function _find_indexing(indexings, indexing)
    for (ix, r) in indexings
        if DynamicPPL.subsumes(ix, indexing)
            return r
        end
    end
    return nothing
end

function getrv(graph::Graph, ref::Reference, indexing)
    _find_indexing(graph.rv_locations[ref], indexing)
end

function setrv!(graph, ((ref, indexing), vn))
    push!(get!(graph.rv_locations, ref, valtype(graph.rv_locations)()),
          indexing => vn)
end

function hasrv(graph, ref, indexing)
    indexings = get(graph.rv_locations, ref, nothing)
    return !isnothing(indexings) && !isnothing(_find_indexing(indexings, indexing))
end

function makereference!(graph, node)
    newref = Reference(length(graph.reference_mapping) + 1)
    setmapping!(graph, node => newref)
    return newref
end

function Base.show(io::IO, graph::Graph)
    for (ref, stmt) in graph
        println(io, ref, " = ", stmt)
    end
end

function convertcall(graph, node::CallingNode)
    f = convertvalue(graph, getfunction(node))
    args = convertvalue.(Ref(graph), getarguments(node))
    return Call(f, args, getvalue(node))
end

convertvalue(graph, expr::TapeReference) = getmapping(graph, expr[], getvalue(expr))
convertvalue(graph, expr::TapeConstant) = getvalue(expr)

convertdist!(graph, dist_expr::TapeConstant) = Constant(getvalue(dist))
function convertdist!(graph, dist_expr::TapeReference)
    ref = getmapping(graph, dist_expr[])
    return ref
end

convertvn!(graph, vn_expr::Nothing) = nothing
convertvn!(graph, vn_expr::TapeConstant) = getvalue(vn_expr)
function convertvn!(graph, vn_expr::TapeReference)
    # @108: ⟨tuple⟩(@82) = (2,)
    # @109: ⟨tuple⟩(@108) = ((2,),)
    # @110: ⟨VarName⟩(⟨:x⟩, @109) = x[2]

    # extract the nodes that compose the indices of a varname
    vn_node = vn_expr[] # @110[]
    vn_name = getvalue(getargument(vn_node, 1))
    indices_node = try_getindex(getargument(vn_node, 2)) # @109[]
    index_element_nodes = getindex.(getarguments(indices_node)) # (@108[],)
    index_nodes = Tuple(try_getindex.(getarguments(index)) for index in index_element_nodes) # ((@82[],),)
    index_refs = map(ix -> getmapping.(Ref(graph), ix, ix), index_nodes)  # ((getmapping(graph, @82[]),),)

    # delete statements of all nodes that were involved in the construction (ie., @108, @109, @110)
    deletemapping!(graph, vn_node)
    deletemapping!(graph, indices_node)
    deletemapping!.(Ref(graph), index_element_nodes)

    return VarName(vn_name, map.(ix -> tovalue(graph, ix), index_refs))
end

function pushtilde!(graph, callingnode, tilde_constructor)
    vn_expr, dist_expr, value_expr, vi_expr = tilde_parameters(callingnode)
    vn = convertvn!(graph, vn_expr)
    dist_ref = convertdist!(graph, dist_expr)
    value = convertvalue(graph, value_expr)

    ref = makereference!(graph, callingnode)
    graph[ref] = tilde_constructor(vn, dist_ref, value, -Inf)

    if !isnothing(vn)
        indexing_values = tovalue.(Ref(graph), DynamicPPL.getindexing(vn))
        if isdotted(tilde_constructor)
            @warn "Broadcasted tildes ($(graph[ref])) are not fully supported!"
            setrv!(graph, (value, indexing_values) => ref)
        elseif indexing_values == ()
            # this is the case where an unindexed variable is sampled -- the value is equivalent to
            # the current reference; otherwise, a `setindex!` should follow in the next line.
            setrv!(graph, (ref, indexing_values) => ref)
        end
    end

    return graph
end

pushnode!(graph, node::CallingNode{typeof(DynamicPPL.tilde_assume)}) =
    pushtilde!(graph, node, Assumption{false})
pushnode!(graph, node::CallingNode{typeof(DynamicPPL.dot_tilde_assume)}) =
    pushtilde!(graph, node, Assumption{true})
pushnode!(graph, node::CallingNode{typeof(DynamicPPL.tilde_observe)}) =
    pushtilde!(graph, node, Observation{false})
pushnode!(graph, node::CallingNode{typeof(DynamicPPL.dot_tilde_observe)}) =
    pushtilde!(graph, node, Observation{true})

function pushnode!(graph, node::CallingNode)
    ref = makereference!(graph, node)
    graph[ref] = convertcall(graph, node)
    return graph
end

function pushnode!(graph, node::CallingNode{typeof(DynamicPPL.matchingvalue)})
    # special handling for model arguments
    # @7: ⟨getproperty⟩(@6, ⟨:x⟩) = [0.5, 1.1]
    # @8: ⟨DynamicPPL.matchingvalue⟩(@4, @3, @7) = [0.5, 1.1]

    value_expr = getargument(node, 3) # @7 above
    if value_expr isa TapeReference
        getproperty_node = value_expr[]
        delete!(graph, getmapping(graph, getproperty_node))
    else # value vas a literal, not a function argument
        @warn "this probably shouldn't have happened..."
    end

    ref = makereference!(graph, node)
    graph[ref] = Constant(getvalue(node))
    return graph
end

function pushnode!(graph, node::CallingNode{typeof(setindex!)})
    # @32: [§5:%34] ⟨DynamicPPL.tilde_assume⟩(@5, @4, @20, @29, @31, @3) = -0.031672055938280076
    # @33: [§5:%35] ⟨setindex!⟩(@19, @32, ⟨1⟩) = [...]
    argument_exprs = try_getindex.(getarguments(node))
    arguments = getmapping.(Ref(graph), argument_exprs, argument_exprs)
    mutated, value, indexing = arguments[1], arguments[2], arguments[3:end]
    indexing_values = (tovalue.(Ref(graph), indexing),)

    if value isa Reference && graph[value] isa Tilde
        # this is the case where in the previous line, an indexed RV was sampled; we need
        # to associate the value with the array position it is assigned to here.
        setrv!(graph, (mutated, indexing_values) => value)

        # # we can also replace the `setindex!` call that followed the tilde node
        # # by a `getindex` on the same value, to preserve the information of the dependency of
        # # between the variable, the array, and the index
        # # PROBABLY UNNEEDED? see reverse_deps test case in test_conditionals
        definition = (graph[value].vn, value)
        ref = makereference!(graph, node)
        graph[ref] = Call(definition, getindex, (mutated, indexing...), tovalue(graph, value))
    else
        invoke(pushnode!, Tuple{Graph, CallingNode}, graph, node)
    end

    return graph
end

function pushnode!(graph, node::CallingNode{typeof(getindex)})
    # @180: ⟨getindex⟩(@9, @135) → -0.05
    argument_exprs = try_getindex.(getarguments(node))
    arguments = getmapping.(Ref(graph), argument_exprs, argument_exprs)
    array, indexing = arguments[1], arguments[2:end]
    indexing_values = (tovalue.(Ref(graph), indexing),)

    if hasrv(graph, array, indexing_values)
        rv = getrv(graph, array, indexing_values)
        ref = makereference!(graph, node)
        definition = (VarName(graph[rv].vn, indexing_values), rv)
        graph[ref] = Call(definition, getindex, (array, indexing...), getvalue(node))
    else
        invoke(pushnode!, Tuple{Graph, CallingNode}, graph, node)
    end

    return graph
end

function pushnode!(graph, node::ConstantNode)
    ref = makereference!(graph, node)
    graph[ref] = Constant(getvalue(node))
    return graph
end

function pushnode!(graph, node::ArgumentNode)
    if !isnothing(node.branch_node)
        # don't add, but remember branch arguments
        original_ref = getmapping(graph, referenced(node)[1])
        setmapping!(graph, node => original_ref)
    else
        # handle function argument nodes like constants
        ref = makereference!(graph, node)
        graph[ref] = Constant(getvalue(node))
    end

    return graph
end

function makegraph(slice::Vector{<:AbstractNode})
    return foldl(slice, init=Graph()) do graph, node
        pushnode!(graph, node)
    end
end


text:
"""
	aMK(x,columns)

Function to estimate the asymptotic value of α(x).

# Arguments
 - `sfs::Matrix{Float64}`: SFS array.
 - `divergence::Matrix{Float64}`: SFS array.
# Returns
 - `Array{Float64,2}`: Array of array containing asymptotic values, lower confidence interval, higher confidence interval
"""
function aMK(α::Vector{Float64})

	# Model
	model(x,p) = @. p[1] + p[2]*exp(-x*p[3])

	# Fit values
	fitted1    = curve_fit(model,collect(1:size(α,1)),α,[-1.0,-1.0,1.0];lower=[-1.0,-1.0,1.0],upper=[1.0, 1.0, 10.0])
	fitted2    = curve_fit(model,collect(1:size(α,1)),α,fitted1.param)
	asymp      = model(size(α,1),fitted2.param)

	ci   = try
		[confidence_interval(fitted2)[1][1],confidence_interval(fitted2)[1][2]]
	catch
		[0.0,0.0]
	end

	return(asymp,ci,fitted2.param)
end

"""
	imputedMK(sfs,divergence,m,cutoff)

Function to estimate the imputedMK

# Arguments
 - `sfs::Matrix{Float64}`: SFS array
 - `divergence::Array`: divegence count array
 - `m::Union{Nothing,Array}`: non-synonymous and synonymous sites
 - `cutoff{Float64}`: frequency cutoff to perform imputedMK
# Returns
 - `Dict: Dictionary containing results
"""
function imputedMK(sfs::Matrix{Float64},divergence::Matrix{Int64},cutoff::Float64=0.15;m::T=nothing) where {T<:Union{Nothing,Array}}

	output = OrderedDict{String,Float64}()

	pn = sum(sfs[:,2])
	ps = sum(sfs[:,3])
	dn = divergence[1]
	ds = divergence[2]

	deleterious = 0
	### Estimating slightly deleterious with pn/ps ratio
	flt_low = (sfs[:, 1] .<= cutoff)
	pn_low   = sum(sfs[flt_low,2])
	ps_low   = sum(sfs[flt_low,3])

	flt_inter = (sfs[:, 1] .>= cutoff) .& (sfs[:, 1] .<= 1)
	pn_inter = sum(sfs[flt_inter,2])
	ps_inter = sum(sfs[flt_inter,3])

	ratio_ps       = ps_low / ps_inter
	deleterious   = pn_low - (pn_inter * ratio_ps)

	if (deleterious > pn) || (deleterious < 0)
		deleterious = 0
		pn_neutral = round(pn - deleterious,digits=3)
	else
		pn_neutral = round(pn - deleterious,digits=3)
	end

	output["alpha"] = round(1 - ((pn_neutral/ps) * (ds / dn)),digits=3)

	#  method = :minlike same results R, python two.sides
	output["pvalue"] = pvalue(FisherExactTest(Int(ps),Int(ceil(pn_neutral)),Int(ds),Int(dn)))

	if (!isnothing(m))

		mn = m[1]; ms = m[2]
		# ## Estimation of b: weakly deleterious
		output["b"] = (deleterious / ps) * (ms / mn)

		## Estimation of f: neutral sites
		output["f"] = (ms * pn_neutral) / (mn * ps)

		## Estimation of d, strongly deleterious sites
		output["d"] = 1 - (output["f"] + output["b"])

		ka      = dn / mn
		ks       = ds / ms
		output["omega"]    = ka/ ks

		# Omega A and Omega D
		output["omegaA"] = output["omega"] * output["alpha"]
		output["omegaD"] = output["omega"] - output["omegaA"]
	end

	return output
end

"""
	fwwMK(sfs,divergence,m,cutoff)

Function to estimate the Fay, Wycoff and Wu's MK extension.

# Arguments
 - `sfs::Matrix{Float64}`: SFS array
 - `divergence::Array`: divegence count array
 - `m::Union{Nothing,Array}`: non-synonymous and synonymous sites
 - `cutoff{Float64}`: frequency cutoff to perform fwwMK
# Returns
 - `Dict: Dictionary containing results
"""
function fwwMK(;sfs::Array,divergence::Array,m::T,cutoff::Float64=0.15) where {T<:Union{Nothing,Array}}

	output = OrderedDict{String,Float64}()

	ps = sum(sfs[:,2])
	pn = sum(sfs[:,3])
	dn = divergence[1]
	ds = divergence[2]

	deleterious = 0
	### Estimating slightly deleterious with pn/ps ratio
	flt_inter = (sfs[:, 1] .>= cutoff) .& (sfs[:, 1] .<= 1)
	pn_inter = sum(sfs[flt_inter,2])
	ps_inter = sum(sfs[flt_inter,3])

	# output['alpha'] = 1 - (((pn - deleterious) / ps) * (ds / dn))
	output["alpha"] = round(1 - ((pn_inter/ps_inter) * (ds / dn)),digits=3)
	#  method = :minlike same results R, python two.sides
	output["pvalue"] = pvalue(FisherExactTest(Int(ps_inter),Int(ceil(pn_inter)),Int(ds),Int(dn)))

	if (!isnothing(m))

		mn = m[1]; ms = m[2]
		ka      = dn / mn
		ks       = ds / ms
		output["omega"]    = ka/ ks

		# Omega A and Omega D
		output["omegaA"] = output["omega"] * output["alpha"]
		output["omegaD"] = output["omega"] - output["omegaA"]
	end

	return output
end

"""
	standardMK(x,columns)

Function to estimate the original α value

# Arguments
 - `sfs::Matrix{Float64}`: SFS array
 - `divergence::Array`: divegence count array
 - `m::Union{Nothing,Array}`: non-synonymous and synonymous sites# Returns

# Output
 - `Dict: Dictionary containing results
"""
function standardMK(;sfs::Array,divergence::Array,m::T=nothing) where {T<:Union{Nothing,Array}}

	output = OrderedDict{String,Float64}()

	pn = sum(sfs[:,2])
	ps = sum(sfs[:,3])
	dn = divergence[1]
	ds = divergence[2]

	output["alpha"] = round(1 - ((pn/ps) * (ds / dn)),digits=3)
	#  method = :mnnlike same results R, python two.sides
	output["pvalue"] = pvalue(FisherExactTest(Int(ps),Int(ceil(pn)),Int(ds),Int(dn)))

	if (!isnothing(m))

		mn = m[1]; ms = m[2]

		ka      = dn / mn
		ks       = ds / ms
		output["omega"]    = ka/ ks

		# Omega A and Omega D
		output["omegaA"] = output["omega"] * output["alpha"]
		output["omegaD"] = output["omega"] - output["omegaA"]
	end

	return output
end


text:
function _precompile_()
  ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
  # for T in (Bool, Int, Float32, Float64)
  #   for A in (Vector, Matrix)
  #     precompile(stridedpointer, (A{T},))
  #     precompile(stridedpointer, (LinearAlgebra.Adjoint{T,A{T}},))
  #   end
  # end

  precompile(offset_ptr, (Symbol, Symbol, Char, Int, Int, Int, Int, Int, Bool, Int))
  precompile(vload_quote_llvmcall_core, (Symbol, Symbol, Symbol, Int, Int, Int, Int, Bool, Bool, Int))
  precompile(vstore_quote, (Symbol, Symbol, Symbol, Int, Int, Int, Int, Bool, Bool, Bool, Bool, Int))

  precompile(Tuple{typeof(transpose_vecunroll_quote_W_smaller),Int,Int})   # time: 0.02420761

  precompile(Tuple{typeof(horizontal_reduce_store_expr),Int,Int,NTuple{4, Int},Symbol,Symbol,Bool,Int,Bool})   # time: 0.02125804
  precompile(Tuple{typeof(transpose_vecunroll_quote_W_larger),Int,Int})   # time: 0.01755242
  precompile(Tuple{typeof(shufflevector_instrs),Int,Type,Vector{String},Int})   # time: 0.0159487
  precompile(Tuple{typeof(transpose_vecunroll_quote),Int})   # time: 0.014891806
  precompile(Tuple{typeof(align),Int,Int})   # time: 0.013784537
  precompile(Tuple{typeof(align),Int})   # time: 0.013609074

  precompile(Tuple{typeof(vstore_transpose_quote),Int64,Int64,Int64,Int64,Int64,Int64,Int64,Bool,Bool,Bool,Int64,Int64,Symbol,UInt64,Bool})   # time: 0.006213663
  precompile(Tuple{typeof(vstore_unroll_i_quote),Int64,Int64,Int64,Bool,Bool,Bool,Int64,Bool})   # time: 0.002936335

  precompile(Tuple{typeof(_shuffle_load_quote), Symbol, Int, NTuple{9,Int}, Symbol, Symbol, Int, Int, Bool, Int, UInt})
  precompile(Tuple{typeof(_shuffle_store_quote), Symbol, Int, NTuple{9,Int}, Symbol, Symbol, Int, Int, Bool, Bool, Bool, Int, Bool})

  precompile(Tuple{typeof(collapse_expr),Int64,Symbol,Int64})   # time: 0.003906299

  # precompile(_pick_vector_width, (Type, Vararg{Type,100}))
    # the `"NATIVE_PRECOMPILE_VECTORIZATIONBASE" ∈ keys(ENV)` isn't respected, seems
    # like it gets precompiled anyway given that the first condition is `true`.
    # if VERSION ≥ v"1.7.0-DEV.346" && "NATIVE_PRECOMPILE_VECTORIZATIONBASE" ∈ keys(ENV)
    #     set_features!()
    #     for T ∈ (Float32, Float64)
    #         W = pick_vector_width(T)
    #         precompile(>=, (Int, MM{W, 1, Int}))
    #         for op ∈ (-, Base.FastMath.sub_fast)
    #             precompile(op, (Vec{W, T}, ))
    #         end
    #         for op ∈ (+, -, *, Base.FastMath.add_fast, Base.FastMath.sub_fast, Base.FastMath.mul_fast)
    #             precompile(op, (Vec{W, T}, Vec{W, T}))
    #         end
    #         for op ∈ (VectorizationBase.vfmadd, VectorizationBase.vfmadd_fast)
    #             precompile(op, (Vec{W, T}, Vec{W, T}, Vec{W, T}))
    #         end
    #     end
    # end
end


text:
export Problems

module Problems

struct Problem
    solutions::Dict{String, Function}
    Problem(problems::Dict{String, Function}) = new(problems)
    Problem(func::Function) = new(Dict("Solution" => func))
end

function benchmark(p::Problem, x...)
    key = collect(keys(p.solutions))[1]
    sol = 0
    # Initial run with default small value
    for k in collect(keys(p.solutions))
        sol = p.solutions[k]()
    end
    for problem in p.solutions
        println("  ", problem[1], ":")
        print("  ")
        @time sol = problem[2](x...)
    end
    println("  Answer: ", sol)
end

end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

const PULSE_INTERVAL_MILLIS = 5000

mutable struct  CoreSession <: AbstractCoreSession
    client::CoreClient
    database::CoreDatabase
    sessionID::Bytes
    transactions::Dict{UUID,T} where {T<:Union{Nothing,<:AbstractCoreTransaction}}
    type::Int
    accessLock::ReentrantLock
    options::TypeDBOptions
    isOpen::Bool
    networkLatencyMillis::Int
    timer::Optional{Controller}
    request_timeout::Real
end

Base.show(io::IO, session::T) where {T<:AbstractCoreSession} = print(io, "Session(ID: $(bytes2hex(session.sessionID)))")

function CoreSession(client::T,
                     database::String,
                     type::Int32,
                     options::TypeDBOptions = TypeDBOptions();
                     request_timout::Real=6) where {T<:AbstractCoreClient}
    try
        #building open_request
        open_req = SessionRequestBuilder.open_req(
            database, type , copy_to_proto(options, Proto.Options)
        )
        # open the session
        startTime = now()
        grpc_controller = gRPCController()
        req_result, status  = Proto.session_open(client.core_stub.blockingStub, grpc_controller, open_req)
        res_id = grpc_result_or_error(req_result, status, result->result.session_id)
        endTime = now()

        database = CoreDatabase(database)
        networkLatencyMillis = (endTime - startTime).value
        session_id =  res_id
        transactions = Dict{UUID,AbstractCoreTransaction}()
        is_open = true

        t = Controller(true, 0.5)

        result = CoreSession(client, database, session_id, transactions, type, ReentrantLock() ,options, is_open, networkLatencyMillis, t, request_timout)

        # initialize the ongoing pulse request
        make_pulse_request(result, t)

        # add the session to the client
        client.sessions[session_id] = result

        # the following transaction is useless for the construction of the the Session
        # and is only to initialize the first compilation of the ProtoBuf machinary
        # underneath the transaction. In the environment here it is possible to control
        # the tick of the pulse request a close as can be. If somebody find a way to
        # avoid the long starting time of the transaction this can be removed.
        trans = nothing
        try
            trans = transaction(result, Proto.Transaction_Type[:READ])
        catch ex
            @debug "First attempt for transaction done and not successful"
        end
        trans !== nothing && close(trans)

        t.duration_in_seconds =  (PULSE_INTERVAL_MILLIS / 1000) - 0.5

        return result
    catch ex
        throw(TypeDBClientException("Error construct a CoreSession",ex))
    end
end

"""
    make_pulse_request(session::T) where {T<:AbstractCoreSession}

This function make a pulse request to keep the session alive.
"""
function make_pulse_request(session::AbstractCoreSession, controller::Controller)
    @async begin
        while controller.running
            try
                pulsreq = SessionRequestBuilder.pulse_req(session.sessionID)
                req_result, status = Proto.session_pulse(session.client.core_stub.blockingStub, gRPCController() , pulsreq)
                result = grpc_result_or_error(req_result,status, result->result)
                @debug "Time: $(Dates.now())"

                if result.alive === false
                    close(session)
                    @info "$session is closed"
                end
            catch ex
                close(session)
                throw(ex)
            end
            sleep(controller.duration_in_seconds)
        end
    end
end

transaction(session::AbstractCoreSession, type::EnumType) = transaction(session, type, typedb_options_core())
function transaction(session::AbstractCoreSession, type::EnumType, options::TypeDBOptions)
    try
        lock(session.accessLock)
        if !session.isOpen
            throw(TypeDBClientException(CLIENT_SESSION_CLOSED, bytes2hex(session.sessionID)))
        end

        transactionRPC = CoreTransaction(session, session.sessionID, type, options)
        session.transactions[transactionRPC.transaction_id] = transactionRPC

        return transactionRPC
    finally
       unlock(session.accessLock)
    end
end

function close(session::AbstractCoreSession)
    try
        lock(session.accessLock)
        if session.isOpen
            for (_,trans) in session.transactions
                safe_close(trans)
            end
            remove_session(session.client, session)
            safe_close(session.timer)

            req = SessionRequestBuilder.close_req(session.sessionID)
            stub = session.client.core_stub.blockingStub
            Proto.session_close(stub, gRPCController(), req )

            session.isOpen = false
        end
    catch ex
        if ex isa gRPCServiceCallException
            @info "Session with id: $(session.sessionID) isn't open on the server"
        else
            @error TypeDBClientException("Unexpected error while closing session ID: $(session.sessionID)",ex)
        end
    finally
        unlock(session.accessLock)
    end
    return true
end

function Base.delete!(session::AbstractCoreSession, trans_id::UUID)
    try
        lock(session.accessLock)
        delete!(session.transactions, trans_id)
    finally
        unlock(session.accessLock)
    end
end

function is_open(session::AbstractCoreSession)
    return session.isOpen
end


text:
IndentStringCache = Dict{Int, ASCIIString}()

indentstr(indent) = get!(IndentStringCache, indent, join([" " for i in 1:indent]))

iprint(s, indent) = print(indentstr(indent), s)

function pp(d, indent = 0)
	print(indentstr(indent))
	show(d)
end

pp(s::String, indent = 0) = print(indentstr(indent), s)

# Pretty-printing of Julia Expr's
function pp(ex::Expr, indent = 0)
	# Skip line nodes unless told to print them
	if ex.head == :line && includeLineNodes == false
		return()
	end

	pp("head: $(ex.head)\n", indent)
	pp("args:\n", indent+1)
	for arg in ex.args
		pp(arg, indent+2)
		println()
	end
end


text:
# distritbutions
@rule typeof(+)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_var(m_in1)
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(mout - min1, vout + vin1)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in1::MultivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_cov(m_in1)
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(mout - min1, vout + vin1)
end

# PointMass
@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::PointMass) = PointMass(mean(m_out) - mean(m_in1))

@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::NormalMeanPrecision) = begin
    min1, win1 = mean_precision(m_in1)
    return NormalMeanPrecision(mean(m_out) - min1, win1)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::NormalMeanPrecision, m_in1::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return NormalMeanPrecision(mout - mean(m_in1), wout)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::MvNormalMeanPrecision) = begin
    min1, win1 = mean_precision(m_in1)
    return MvNormalMeanPrecision(mean(m_out) - min1, win1)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::MvNormalMeanPrecision, m_in1::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return MvNormalMeanPrecision(mout - mean(m_in1), wout)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::UnivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_var(m_in1)
    return NormalMeanVariance(mean(m_out) - min1, vin1)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass) = begin
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(mout - mean(m_in1), vout)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::MultivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_cov(m_in1)
    return MvNormalMeanCovariance(mean(m_out) - min1, vin1)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in1::PointMass) = begin
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(mout - mean(m_in1), vout)
end