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

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

# ╔═╡ e5020bc7-a803-47c5-b38e-8879937c8a5f
using PlutoUI, CUDA, BenchmarkTools

# ╔═╡ 04b9a546-dc75-440e-bc74-9bbcc945f5f4
begin
	using Statistics
	matriz_datos = rand(Float32, 2500, 2500)
	matriz_datos_gpu = CuArray(matriz_datos)
end

# ╔═╡ f3a89703-2903-4874-9a57-5f84e401c970
using Tullio, CUDAKernels, KernelAbstractions

# ╔═╡ bbbb53fe-f599-11eb-02d3-915af1c0597f
md"""
# Uso de CUDA en Julia

El uso de CUDA en Julia se presenta mediante tres formas:
- Paralelismo implícito por drop-ins.
- Paralelismo explícito al escribir kernels.
- Interfaces con las bibliotecas y el API por medio de wrappers.
"""

# ╔═╡ 98c9d6f2-fa7a-4bb9-acee-02c4022ea03a
PlutoUI.TableOfContents(aside=true, title="CUDA")

# ╔═╡ c0a94327-49b0-426e-bb40-ad7224e80b43
with_terminal(CUDA.versioninfo)

# ╔═╡ d0c4915a-85e4-49fb-b931-c4be3ff0d142
md"""
# Programación con arreglos

Los arreglos del tipo `CuArray` están diseñados para manejar una gran parte de las operaciones típicas de arreglos en el GPU. De esta forma, el código generado en Julia puro puede ser utilizado tan solo redefiniendo los arreglos como CuArrays.

Detrás de las escenas, CUDA.jl se encarga de lo siguiente al generar un arreglo:
- Reservar la memoria en el GPU.
- Si es necesario, realizar la transferencia de memoria CPU ➡ GPU.
- Lanzar kernels para realizar la operación solicitada.
- Si el usuario lo solicita, realizar la transferencia de memoria CPU ⬅ GPU.

Por lo general, se recomienda solo utilizar operaciones vectorizadas. Algunas podrían requerir cargar el arreglo a CPU e iterarlo, y esto es lento (las denominadas operaciones escalares). CUDA.jl nos mostrará una advertencia al intentarlo para evitar malas prácticas, pero debido a que Pluto utiliza este tipo de operaciones para mostrar su salida, le diremos a CUDA.jl que permita operaciones escalares.
"""

# ╔═╡ 15050436-2551-4136-b30d-3d3346493c7b
CUDA.allowscalar(true)

# ╔═╡ 07d2caf0-f987-474d-a462-f452a42be99f
a = CuArray([1, 2, 3, 4])

# ╔═╡ ca37cd69-eb08-4002-80dd-857ef77c3689
md"""Realicemos la primera operación en el arreglo, demostrando el uso del broadcasting."""

# ╔═╡ 456b2857-c2e4-496c-8f36-c8b58fa4c497
a .+= 1

# ╔═╡ 25389b61-4e40-4c11-9d28-45d168e39214
md""" De esta forma, las operaciones diseñadas en paralelo se pueden utilizar como si se tratara de un arreglo normal de Julia, por ejemplo:"""

# ╔═╡ c5ffc0b4-d5d7-48fd-a518-7c2c98c45f82
map(x -> x^2 + 3x - 2, a)

# ╔═╡ c4b9b898-d1c2-4890-8804-c7258f0fdc94
md"""Nótese que no existe modificación alguna al arreglo original:"""

# ╔═╡ 5917412e-80a6-4d9d-b687-611b99e96dac
a

# ╔═╡ a3a61808-148b-465c-8f6f-2483cd6421bd
md"""
Ya que CUDA.jl y Julia siguen un paradigma funcional, en el cual es preferible que las estructuras creadas nunca se modifiquen.

Volviendo al tema del paralelismo, podemos comparar una implementación vectorizada con una secuencial:

**Implementación del algoritmo basado en map:**
"""

# ╔═╡ 61a34fe6-83a5-4508-80af-c449fdbba7d1
@bind tipo_algoritmo Radio(["vect" => "Map (CUDA)", "iter" => "Iterativo (CUDA)", "iter_j" => "Map (Julia)"], default="vect")

# ╔═╡ f5820a28-20ef-4fb3-a466-515ebedd0676
arr1 = collect(1:1_000_000)

# ╔═╡ 04b92cdc-f0d5-43a4-b0b6-31e044ee841b
arr1_gpu = CuArray(arr1)

# ╔═╡ 9407c96c-03bc-4066-a9ef-14ae7e5037e1
if tipo_algoritmo == "vect"
	#El uso de CUDA.@sync está recomendado para medir el tiempo real
	@benchmark CUDA.@sync map(x -> x^2 + 3x - 2, arr1_gpu)
elseif tipo_algoritmo == "iter"
	@benchmark for i in 1:1_000_000
			x = arr1_gpu[i]
			x^2 + 3x - 2
		end
else
	@benchmark map(x -> x^2 + 3x - 2, arr1)
end

# ╔═╡ 47501814-869e-4c2b-9918-323600dc3d63
md"""
Otro ejemplo:

**Implementación del algoritmo (cov):**
"""

# ╔═╡ 2680dbf0-4af4-4d18-beaf-8961203fc2cd
@bind tipo_algoritmo_2 Radio(["vect" => "CUDA", "iter" => "Julia"], default="vect")

# ╔═╡ 8c315cd4-946b-4c6e-9152-173c18efd7fa
if tipo_algoritmo_2 == "vect"
	@benchmark CUDA.@sync cov(matriz_datos_gpu)
else
	@benchmark cov(matriz_datos)
end

# ╔═╡ 61f7b6f6-698b-4f00-84c2-90b288a3c3a7
md"""
Como se puede notar, la implementación de CUDA con **operaciones fáciles de paralelizar** está en otro nivel de magnitud en cuanto a tiempo de ejecución cuando se le compara con una implementación en Julia corriendo en el CPU.
"""

# ╔═╡ 1f74f846-6d8d-48f5-9dc1-391d59cbc164
md"""
## Creación de arreglos

Como se mencionó antes, los arreglos están diseñados para ser reemplazos transparentes al usuario. Un arreglo puede iniciar su vida de la siguiente forma:
"""

# ╔═╡ 04987ac0-cf32-42c6-884f-5d19c4aafec9
arreglo_vacio = CuArray{Float32}(undef, 10)

# ╔═╡ 436a254a-124f-4686-8da9-9fcc560459eb
md"""
Como se puede notar, al estar indefinido no se realiza limpieza de la memoria que existe en la GPU, por lo que se encuentra lleno de datos basura. No obstante, es ligeramente más rápido reservarlo así que inicializándolo.
"""

# ╔═╡ a5a101bc-c5cc-445d-8cbd-181cf6d789d8
arreglo_inicializado = CuArray{Float32}([x^2 + 3x - 2 for x in 1:10])

# ╔═╡ bea5f1e1-e5ba-4f8e-848d-27313e0ae4aa
md"""
Se pueden copiar los contenidos entre arreglos:
"""

# ╔═╡ 4ed549ba-c6b9-4d75-9d2d-5fd2a35081cd
copyto!(arreglo_vacio, arreglo_inicializado)

# ╔═╡ ae014fc7-2c1a-43b9-8dac-faf4ebe14769
arreglo_vacio

# ╔═╡ 89329729-3f50-4e82-bc68-b59bcbc91835
md"""
Llenar un arreglo de un valor arbitrario:
"""

# ╔═╡ 42419f93-26a7-4c28-80d3-db7547dee130
fill!(arreglo_vacio, 7.0)

# ╔═╡ 9074daea-7744-4dd0-bba3-9c727ff93b7f
arreglo_vacio

# ╔═╡ 1b47496e-9f3f-4b4b-b19d-e7e5825869df
md"""
Y generar arreglos con las dimensiones de otro:
"""

# ╔═╡ 9337632e-7bbb-4346-acc8-d0c0eddfd1df
similar(arreglo_vacio)

# ╔═╡ e9982200-fa77-4f41-a201-493f5772fa2a
md"""
Por limitaciones de como funciona la arquitectura CUDA, algunas operaciones que no toman un tipo de arreglo en Julia no pueden ser traducidas directamente a una versión para GPU. Para estos, CUDA.jl implementa una versión propia:
"""

# ╔═╡ 4df35ae7-da44-437f-a762-03ad9b66a15a
md"""
Reemplazando a la función ones:
"""

# ╔═╡ 215d87c7-6ce8-434b-b1cc-6979b85ee2d2
arreglo_unos = CUDA.ones(length(arreglo_vacio))

# ╔═╡ 3234e223-7e03-4258-b7b1-f3fbaf7d4c49
md"""
Reemplazando a la función rand:
"""

# ╔═╡ ecee0116-16e2-41c7-971b-36ecde535260
arreglo_aleatorio = CUDA.rand(length(arreglo_vacio))

# ╔═╡ 4eab8c06-be7c-40d2-94e2-e23471b2a3f5
begin
	using Plots
	plot(arreglo_aleatorio)
	plot!(log.(accumulate(+, arreglo_unos)))
	title!("Un plot con CuArrays")
end

# ╔═╡ 90544984-f317-4e43-9f65-76a358f50e21
md"""
Reemplazando a un arreglo en uso diario:
"""

# ╔═╡ 9c8829c2-7601-48d3-aa58-071926870691
md"""
# Uso de bibliotecas de Nvidia

CUDA.jl provee acceso a las bibliotecas por medio de funciones al estilo de Julia en muchos casos. Para los restantes, da la posibilidad de llamar a la biblioteca en CUDA C directamente.
"""

# ╔═╡ 21257da0-d16b-48b3-8c66-4309cdaa6c7d
cuda_matriz_R = CuArray{Float32}(undef, (10,10))

# ╔═╡ d541be46-5341-4994-be70-cc2d95104d65
begin
	using AbstractFFTs
	fft(cuda_matriz_R)
end

# ╔═╡ fabbde47-37fa-4c5e-9dae-28a1bbc95d4a
cuda_matriz = CUDA.ones(10,10)

# ╔═╡ 32081f1d-7a00-490b-94db-a070e17a0b8a
begin
	using LinearAlgebra
	mul!(cuda_matriz_R, CUDA.rand(10,10), cuda_matriz)
end

# ╔═╡ 804c9f33-4dc0-4763-9137-125d14469bd0
begin
	using Random
	rand!(cuda_matriz)
end

# ╔═╡ 86d6ea0e-a64b-445c-86ae-1479a12519a5
md"""
**CUBLAS:**
"""

# ╔═╡ cd43ba66-5613-4501-b42c-4092fe4ce9e7
md"""
**CURAND:**
"""

# ╔═╡ 2f0edcd1-18b3-4b69-a2a1-278bfd343e3e
md"""
**CUSOLVER:**
"""

# ╔═╡ c9bbb09b-c515-4dcd-8a63-36ff47f5ed4f
qr(cuda_matriz_R)

# ╔═╡ e0261365-1cf5-415d-a556-5160fe949a7e
md"""
**CUFFT:**
"""

# ╔═╡ 985f163d-8404-4d26-ae5c-cb4422095eeb
md"""
Un ejemplo utilizando las APIs de C propias a las bibliotecas:
"""

# ╔═╡ 77c8d1b6-a90a-43c1-8110-20023450313c
begin
	out = Ref{Cint}()
	#Encontrar el índice con el valor mínimo en un arreglo dado
	CUBLAS.cublasIsamin_v2(CUBLAS.handle(), length(arreglo_inicializado),
		arreglo_inicializado, stride(arreglo_inicializado, 1), out)
	out[]
end

# ╔═╡ b4063f3e-2031-43ea-b0ad-a398c7f765d1
findmin(arreglo_inicializado)

# ╔═╡ f3457168-6da7-4daa-957e-bdf187d1aa29
md"""
Poder utilizar este tipo de llamadas a las bibliotecas permite tener la funcionalidad completa de las mismas.

# Abstracciones de alto nivel

La ventaja de CUDA.jl frente a implementaciones como las de Python (por ejemplo CuPy) es la flexibilidad de las abstracciones de Julia en conjunto a los CuArray:

"""

# ╔═╡ 0022e1db-f5e4-47f1-996e-ff12a3022e5a
begin
	arreglo_prueba = CUDA.ones(10)
	broadcast(arreglo_prueba) do x
		x+= 1
	end
end

# ╔═╡ 3c88f1f0-753c-4182-8458-ac97a7e375f3
#Equivalente a:
arreglo_prueba .+ 1

# ╔═╡ 877ca78d-d8c9-4a42-9332-a5c4979ac822
md"""
Al componer varias de estas abstracciones, es posible generar código muy poderoso y conciso:
"""

# ╔═╡ b2195cfe-043a-430f-9ee0-d6ac7f745d71
map(arreglo_prueba) do x
	x + log(x) - exp(x)^π
end

# ╔═╡ 1791ff51-f819-420a-bf03-820fc48313d7
reduce(+, arreglo_prueba)

# ╔═╡ a2843d0e-b520-4ab8-81ac-cf8b1ecf1d40
accumulate(+, arreglo_prueba)

# ╔═╡ bad0aa36-6cc1-49ff-ab3e-9db894055c64
md"""
Algo a tomar en cuenta es que CUDA.jl **no es un compilador de tensores**. Esto quiere decir que implementa las abstracciones directamente, y genera kernels escalares. Cuando se tienen operaciones múltiples, esstas no se suelen fusionar entre ellas.

Entre otras cosas, las operaciones genéricas que operan en otros arreglos no resultan directamente en un kernel para GPU:
"""

# ╔═╡ 80f29cd8-0c01-444c-8ea3-8f6655b99d05
broadcast(eachcol(cuda_matriz)) do x
	sum(x)
end

# ╔═╡ 68919e5f-18da-4bbf-ad66-2720dd556f32
md"""
Realizar la suma dentro del map resulta en lanzar un kernel distinto por columna, lo cuál es lento y no utiliza la GPU de forma eficiente.
"""

# ╔═╡ ce737b93-8409-4f2a-895c-e6dee7dcf570
@benchmark CUDA.@sync broadcast(eachcol($cuda_matriz)) do x
	sum(x)
end

# ╔═╡ f8719bda-5f64-4b2a-a5c3-93343116c0ff
md"""
Se debería entonces optar por una sola invocación:
"""

# ╔═╡ 8c7f0eb9-8db3-4e7d-a27e-a38d08f530b8
sum(cuda_matriz, dims=1)

# ╔═╡ fcadec48-a6cd-4471-95b2-4da9327eac37
@benchmark CUDA.@sync sum(cuda_matriz, dims=1)

# ╔═╡ 8e8ac21a-7310-445e-b515-d6e718f2877e
md"""
En otros casos, la operación no se compilará ya que pocas operaciones de arreglos están soportadas explícitamente en CUDA.jl:
"""

# ╔═╡ 78fab451-9329-41ea-b1df-97c57cb6c0a4
broadcast(cuda_matriz, eachcol(cuda_matriz)) do x, Y
	x + sum(Y)
end

# ╔═╡ 3370dabd-30ba-4233-a60d-0c563ffd82fe
md"""
Por suerte, existen otros paquetes desarrollados con soporte para CUDA.jl. Por ejemplo Tullio.jl, que permite generar expresiones complejas (utilizando notación de Einstein) de tal forma que se compilen como un solo kernel.
"""

# ╔═╡ 4855b2d0-3084-437a-9bff-d56bf5e4eef9
@tullio cuda_matriz_2[i,j] := cuda_matriz[i, j] + cuda_matriz[k, j]

# ╔═╡ 40ca06d9-6845-4af0-b5ed-4e3827156914
@benchmark @tullio cuda_matriz_3[i,j] := cuda_matriz[i, j] + cuda_matriz[k, j]

# ╔═╡ 3804d33f-4908-4ff8-841a-d80a9f53e6ee
md"""
Para casos en los que se requiere más control, se pueden escribir kernels utilizando Julia.

# Programación de kernels

Muchas operaciones con paralelísmo implícito se pueden ejecutar en paralelo dentro del GPU con la ayuda de CUDA.jl. En otros casos, la persona que programa puede crear un kernel para realizar la operación (y adquiere la responsabilidad que esto conlleva).

En un alto nivel, esto es sencillo. Se debe de escribir una función escalar, lanzándola en paralelo utilizando el macro @cuda y el argumento threads.
"""

# ╔═╡ 7bc4bffa-3faf-42f7-acc4-4e8c28181271
function vadd(v_c, v_a, v_b)
	i = threadIdx().x
	v_c[i] = v_a[i] + v_b[i]
	return
end

# ╔═╡ e6d23eb4-c10d-489f-8390-f344dd4b81db
let
	v_a = CuArray(1:10)
	v_b = CuArray(2:2:20)
	v_c = similar(v_a)

	@cuda threads = length(v_a) vadd(v_c, v_a, v_b)

	v_c
end

# ╔═╡ 1d2f4bc8-64bd-489e-ac68-f9d5c21f11b2
md"""
La programación de kernels rápidamente se vuelve más complicada debido a que:
- Se deben respetar las limitaciones del hardware.
- Se deben usar los recursos de hardware eficientemente.
- No todas las operaciones se pueden implementar fácilmente en un kernel escalar.
"""

# ╔═╡ 74fdc782-cba9-417a-b515-9ad8c38ad2b3
md"""
## Limitaciones de hardware

Por ejemplo, se puede exceder el número de hilos por bloque:
"""

# ╔═╡ 0b9eeec1-fbb6-4408-9b42-79f20017fb74
let
	v_a = CuArray(1:100_000)
	v_b = CuArray(2:2:200_000)
	v_c = similar(v_a)

	@cuda threads = length(v_a) vadd(v_c, v_a, v_b)
end

# ╔═╡ e40ffe89-816e-4e5e-996c-4fb2ba05d244
md"""
Podemos visualizar el número máximo de hilos por bloque:
"""

# ╔═╡ 40270d30-6259-4641-afee-d4b74571995c
CUDA.attribute(device(), CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)

# ╔═╡ 13e45f04-6626-46e8-9960-f63b0db6b01b
md"""
Podemos resolver este primer problema al utilizar más de un bloque. Para esto habrá que adaptar el cálculo de índice para que este sea válido.
"""

# ╔═╡ 4d83c990-81f2-4b94-8a8a-95905df93594
function vadd_bloques(c, a, b)
	i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
	if i <= length(a)
		c[i] = a[i] + b[i]
	end
	return
end

# ╔═╡ 08a44e3e-daac-4a07-a1d8-a1bec02c0349
let
	v_a = CuArray(1:100_000)
	v_b = CuArray(2:2:200_000)
	v_c = similar(v_a)

	@cuda threads=1024 blocks=cld(length(v_a), 1024) vadd_bloques(v_c, v_a, v_b)

	v_c
end

# ╔═╡ 5e51ce22-ec1d-47c6-ab5c-eec254b0ec72
md"""
Una forma más sencilla de calcular los recursos necesarios: **Occupancy API**
"""

# ╔═╡ 6cc5f45d-a149-4603-b5f5-2b895e37d134
let
	v_a = CuArray(1:100_000)
	v_b = CuArray(2:2:200_000)
	v_c = similar(v_a)

	kernel = @cuda launch = false vadd_bloques(v_c, v_a, v_b)

	config = CUDA.launch_configuration(kernel.fun)
end

# ╔═╡ ee8da3f3-ca78-433c-8510-40c0ea8958ba
let
	v_a = CuArray(1:100_000)
	v_b = CuArray(2:2:200_000)
	v_c = similar(v_a)

	t = min(length(v_a), 1024)
	b = cld(length(v_a), t)

	@cuda threads=t blocks=b vadd_bloques(v_c, v_a, v_b)

	(t, b, v_c)
end

# ╔═╡ 5b8de2c3-23e3-4f80-b95a-d01b3397bb2b
md"""
En este caso, el API de CUDA calculó que lo más óptimo para nuestra operación sobre esas entradas sería lanzar 40 bloques de 1024 hilos cada uno. Realmente esto no sería suficiente para nuestra operación dado el kernel actual, pero el kernel se puede modificar para que funcione de esa manera.
"""

# ╔═╡ d9689381-94a4-45a6-812a-8ca6e3027c99
md"""
## Operaciones atómicas

CUDA permite operaciones atómicas en sus kernels, pero estas siguien siendo bastante costosas.
"""

# ╔═╡ 307933d4-af70-4fcc-a0b9-7818487e88e6
function reduce_atomic(op, a, b)
	i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
	@atomic b[] = op(b[], a[i])
	return
end

# ╔═╡ 8df9cc05-c51b-4675-9441-12f9b992222e
let
	arreglo = CuArray(1:16)
	salida = CuArray([0])

	@cuda threads=length(arreglo) reduce_atomic(+, arreglo, salida)

	salida
end

# ╔═╡ 56c80ef5-d974-4a7c-9b2c-adf9bf3282b5
sum(1:16)

# ╔═╡ b6847097-a4aa-476a-b8d8-d354168553a5
md"""
Una ventaja clara de Julia se puede notar aquí. El código Julia toma una función como argumento y lo puede aplicar debido a la compilación del código a kernel.
"""

# ╔═╡ 87089c7e-db1d-42d2-bdec-57e3c42801b7
md"""
## Algunos ejemplos más avanzados

Este ejemplo muestra el cálculo de la distancia de Haversine (es decir, entre dos puntos de una esfera).

"""

# ╔═╡ bdd4e398-0049-4b1c-8114-bb7ab786d854
function haversine_gpu(lat1::Float32, lon1::Float32, lat2::Float32, lon2::Float32, radius::Float32)
    c1 = CUDA.cospi(lat1 / 180.0f0)
    c2 = CUDA.cospi(lat2 / 180.0f0)
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    d1 = CUDA.sinpi(dlat / 360.0f0)
    d2 = CUDA.sinpi(dlon / 360.0f0)
    t = d2 * d2 * c1 * c2
    a = d1 * d1 + t
    c = 2.0f0 * CUDA.asin(CUDA.min(1.0f0, CUDA.sqrt(a)))
    return radius * c
end

# ╔═╡ 834a723b-52cf-4cc8-a815-9175be6370b6
md"""
En la siguiente función, se muestran algunas de las características más avanzadas de CUDA.jl:
- Usar índices bidimensionales para bloques e hilos.
- Utilizar memoria compartida dentro del bloque en vez de solo global.
- Funciones en dispositivo (llamadas dentro de la GPU para ejecutar en la GPU).
- Sincronización de hilos dentro de un kernel.

Esta se encarga de calcular las distancias de todos los puntos a todos los puntos.
"""

# ╔═╡ 0407b18f-216f-47d2-8102-ad58326f9baa
# pairwise distance calculation kernel
function pairwise_dist_kernel(lat::CuDeviceVector{Float32}, lon::CuDeviceVector{Float32},
                              rowresult::CuDeviceMatrix{Float32}, n)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    j = (blockIdx().y-1) * blockDim().y + threadIdx().y

    if i <= n && j <= n
        # store to shared memory
        shmem = @cuDynamicSharedMem(Float32, 2*blockDim().x + 2*blockDim().y)
        if threadIdx().y == 1
            shmem[threadIdx().x] = lat[i]
            shmem[blockDim().x + threadIdx().x] = lon[i]
        end
        if threadIdx().x == 1
            shmem[2*blockDim().x + threadIdx().y] = lat[j]
            shmem[2*blockDim().x + blockDim().y + threadIdx().y] = lon[j]
        end
        sync_threads()

        # load from shared memory
        lat_i = shmem[threadIdx().x]
        lon_i = shmem[blockDim().x + threadIdx().x]
        lat_j = shmem[2*blockDim().x + threadIdx().y]
        lon_j = shmem[2*blockDim().x + blockDim().y + threadIdx().y]

        @inbounds rowresult[i, j] = haversine_gpu(lat_i, lon_i, lat_j, lon_j, 6372.8f0)
    end

    return
end

# ╔═╡ f9ec57fc-f7d1-4231-a873-0371250d48fd
let
	n = 10_000

	lat_gpu = CUDA.rand(n) .* 45
    lon_gpu = CUDA.rand(n) .* -120

    #Matriz resultante
    rowresult_gpu = CUDA.zeros(n, n)

    # calculate a 2D block size from the suggested 1D configuration
    # NOTE: we want our launch configuration to be as square as possible,
    #       because that minimizes shared memory usage
    function get_threads(threads)
        threads_x = floor(Int, sqrt(threads))
        threads_y = threads ÷ threads_x
        return (threads_x, threads_y)
    end

    # calculate the amount of dynamic shared memory for a 2D block size
    get_shmem(threads) = 2 * sum(threads) * sizeof(Float32)

    kernel = @cuda launch=false pairwise_dist_kernel(lat_gpu, lon_gpu, rowresult_gpu, n)
    config = launch_configuration(kernel.fun, shmem=threads->get_shmem(get_threads(threads)))

    # convert to 2D block size and figure out appropriate grid size
    t = get_threads(config.threads)
    b = ceil.(Int, n ./ t)
    s = get_shmem(t)

    kernel(lat_gpu, lon_gpu, rowresult_gpu, n; threads=t, blocks=b, shmem=s)

	rowresult_gpu
end

# ╔═╡ 2597a206-5396-4c33-98d9-7ea6b2dbbeca
md"""
Como se puede notar, CUDA.jl tiene un grado de flexibilidad alto. Por un lado, puede servir como un reemplazo para operaciones complejas tan solo reemplazando el tipo de arreglo que se utiliza. Por el otro lado, puede utilizar optimizaciones de bajo nivel de CUDA, enriqueciéndolas con el poder del lenguaje puro.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57"
KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

[compat]
AbstractFFTs = "~1.0.1"
BenchmarkTools = "~1.1.1"
CUDA = "~3.3.4"
CUDAKernels = "~0.3.0"
KernelAbstractions = "~0.7.0"
Plots = "~1.20.0"
PlutoUI = "~0.7.9"
Tullio = "~0.3.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[BFloat16s]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "4af69e205efc343068dc8722b8dfec1ade89254a"
uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
version = "0.1.0"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
git-tree-sha1 = "c31ebabde28d102b602bada60ce8922c266d205b"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.1.1"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.6+5"

[[CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[CUDA]]
deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "DataStructures", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "TimerOutputs"]
git-tree-sha1 = "5e696e37e51b01ae07bd9f700afe6cbd55250bce"
uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
version = "3.3.4"

[[CUDAKernels]]
deps = ["Adapt", "CUDA", "Cassette", "KernelAbstractions", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "81f76297b63c67723b1d60f5e7e002ae3393974b"
uuid = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57"
version = "0.3.0"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "e2f47f6d8337369411569fd45ae5753ca10394c6"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.0+6"

[[Cassette]]
git-tree-sha1 = "087e76b8d48c014112ba890892c33be42ad10504"
uuid = "7057c7e9-c182-5462-911a-8362d720325c"
version = "0.3.7"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.0"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
git-tree-sha1 = "ed268efe58512df8c7e224d2e170afd76dd6a417"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.13.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "344f143fa0ec67e47917848795ab19c6a455f32c"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.32.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "85d2d9e2524da988bffaf2a381864e20d2dae08d"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.2.1"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.1.5+1"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "3cc57ad0a213808473eafef4845a74766242e05f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.3.1+4"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "35895cf184ceaab11fd778b4590144034a167a2f"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.1+14"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "cbd58c9deb1d304f5a245a0b7eb841a2560cfec6"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.1+5"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+0"

[[GPUArrays]]
deps = ["AbstractFFTs", "Adapt", "LinearAlgebra", "Printf", "Random", "Serialization", "Statistics"]
git-tree-sha1 = "ececbf05f8904c92814bdbd0aafd5540b0bf2e9a"
uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
version = "7.0.1"

[[GPUCompiler]]
deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "TimerOutputs", "UUIDs"]
git-tree-sha1 = "f26f15d9c353f7091065390ea826df9e03917e58"
uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
version = "0.12.8"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "182da592436e287758ded5be6e32c406de3a2e47"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.58.1"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "d59e8320c2747553788e4fc42231489cc602fa50"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.58.1+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "7bf67e9a481712b3dbe9cb3dac852dc4b1162e02"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "44e3b40da000eab4ccb1aecdc4801c040026aeb5"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.13"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[KernelAbstractions]]
deps = ["Adapt", "Cassette", "InteractiveUtils", "MacroTools", "SpecialFunctions", "StaticArrays", "UUIDs"]
git-tree-sha1 = "5e6c70389c1b1e40adb81664ca8cea6ce8127afc"
uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
version = "0.7.0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LLVM]]
deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"]
git-tree-sha1 = "d6041ad706cf458b2c9f3e501152488a26451e9c"
uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
version = "4.2.0"

[[LLVMExtra_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "a9b1130c4728b0e462a1c28772954650039eb847"
uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
version = "0.0.7+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[LibVPX_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "12ee7e23fa4d18361e7c2cde8f8337d4c3101bc7"
uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a"
version = "1.10.0+0"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "LinearAlgebra"]
git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.2.5"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "501c20a63a34ac1d015d5304da0e645f42d91c9f"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.11"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "e39bea10478c6aff5495ab522517fae5134b40e3"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.20.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Random123]]
deps = ["Libdl", "Random", "RandomNumbers"]
git-tree-sha1 = "0e8b146557ad1c6deb1367655e052276690e71a3"
uuid = "74087812-796a-5b5d-8853-05524746bad3"
version = "1.4.2"

[[RandomNumbers]]
deps = ["Random", "Requires"]
git-tree-sha1 = "a752043df7488ca8bcbe05fa82c831b5e2c67211"
uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
version = "1.5.2"

[[RecipesBase]]
git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.1"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "2a7a2469ed5d94a98dea0e85c46fa653d76be0cd"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.3.4"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "508822dca004bf62e210609148511ad03ce8f1d8"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.0"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3fedeffc02e47d6e3eb479150c8e5cd8f15a77a0"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.10"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "fed1ec1e65749c4d96fc20dd13bea72b55457e62"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.9"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "000e168f5cc9aded17b6999a560b7c11dda69095"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.0"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TimerOutputs]]
deps = ["ExprTools", "Printf"]
git-tree-sha1 = "209a8326c4f955e2442c07b56029e88bb48299c7"
uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
version = "0.5.12"

[[Tullio]]
deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"]
git-tree-sha1 = "0288b7a395fc412952baf756fac94e4f28bfec65"
uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
version = "0.3.2"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"]
git-tree-sha1 = "2839f1c1296940218e35df0bbb220f2a79686670"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.18.0+4"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "acc685bcf777b2202a904cdcb49ad34c2fa1880c"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.14.0+4"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7a5780a0d9c6864184b3a2eeeb833a0c871f00ab"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "0.1.6+4"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d713c1ce4deac133e3334ee12f4adff07f81778f"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2020.7.14+2"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "487da2f8f2f0c8ee0e83f39d13037d6bbf0a45ab"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.0.0+3"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─bbbb53fe-f599-11eb-02d3-915af1c0597f
# ╠═e5020bc7-a803-47c5-b38e-8879937c8a5f
# ╠═98c9d6f2-fa7a-4bb9-acee-02c4022ea03a
# ╠═c0a94327-49b0-426e-bb40-ad7224e80b43
# ╟─d0c4915a-85e4-49fb-b931-c4be3ff0d142
# ╠═15050436-2551-4136-b30d-3d3346493c7b
# ╠═07d2caf0-f987-474d-a462-f452a42be99f
# ╟─ca37cd69-eb08-4002-80dd-857ef77c3689
# ╠═456b2857-c2e4-496c-8f36-c8b58fa4c497
# ╟─25389b61-4e40-4c11-9d28-45d168e39214
# ╠═c5ffc0b4-d5d7-48fd-a518-7c2c98c45f82
# ╟─c4b9b898-d1c2-4890-8804-c7258f0fdc94
# ╠═5917412e-80a6-4d9d-b687-611b99e96dac
# ╟─a3a61808-148b-465c-8f6f-2483cd6421bd
# ╟─61a34fe6-83a5-4508-80af-c449fdbba7d1
# ╠═f5820a28-20ef-4fb3-a466-515ebedd0676
# ╠═04b92cdc-f0d5-43a4-b0b6-31e044ee841b
# ╠═9407c96c-03bc-4066-a9ef-14ae7e5037e1
# ╟─47501814-869e-4c2b-9918-323600dc3d63
# ╟─2680dbf0-4af4-4d18-beaf-8961203fc2cd
# ╠═04b9a546-dc75-440e-bc74-9bbcc945f5f4
# ╠═8c315cd4-946b-4c6e-9152-173c18efd7fa
# ╟─61f7b6f6-698b-4f00-84c2-90b288a3c3a7
# ╟─1f74f846-6d8d-48f5-9dc1-391d59cbc164
# ╠═04987ac0-cf32-42c6-884f-5d19c4aafec9
# ╟─436a254a-124f-4686-8da9-9fcc560459eb
# ╠═a5a101bc-c5cc-445d-8cbd-181cf6d789d8
# ╟─bea5f1e1-e5ba-4f8e-848d-27313e0ae4aa
# ╠═4ed549ba-c6b9-4d75-9d2d-5fd2a35081cd
# ╠═ae014fc7-2c1a-43b9-8dac-faf4ebe14769
# ╟─89329729-3f50-4e82-bc68-b59bcbc91835
# ╠═42419f93-26a7-4c28-80d3-db7547dee130
# ╠═9074daea-7744-4dd0-bba3-9c727ff93b7f
# ╟─1b47496e-9f3f-4b4b-b19d-e7e5825869df
# ╠═9337632e-7bbb-4346-acc8-d0c0eddfd1df
# ╟─e9982200-fa77-4f41-a201-493f5772fa2a
# ╟─4df35ae7-da44-437f-a762-03ad9b66a15a
# ╠═215d87c7-6ce8-434b-b1cc-6979b85ee2d2
# ╟─3234e223-7e03-4258-b7b1-f3fbaf7d4c49
# ╠═ecee0116-16e2-41c7-971b-36ecde535260
# ╟─90544984-f317-4e43-9f65-76a358f50e21
# ╠═4eab8c06-be7c-40d2-94e2-e23471b2a3f5
# ╟─9c8829c2-7601-48d3-aa58-071926870691
# ╠═21257da0-d16b-48b3-8c66-4309cdaa6c7d
# ╠═fabbde47-37fa-4c5e-9dae-28a1bbc95d4a
# ╟─86d6ea0e-a64b-445c-86ae-1479a12519a5
# ╠═32081f1d-7a00-490b-94db-a070e17a0b8a
# ╟─cd43ba66-5613-4501-b42c-4092fe4ce9e7
# ╠═804c9f33-4dc0-4763-9137-125d14469bd0
# ╟─2f0edcd1-18b3-4b69-a2a1-278bfd343e3e
# ╠═c9bbb09b-c515-4dcd-8a63-36ff47f5ed4f
# ╟─e0261365-1cf5-415d-a556-5160fe949a7e
# ╠═d541be46-5341-4994-be70-cc2d95104d65
# ╟─985f163d-8404-4d26-ae5c-cb4422095eeb
# ╠═77c8d1b6-a90a-43c1-8110-20023450313c
# ╠═b4063f3e-2031-43ea-b0ad-a398c7f765d1
# ╟─f3457168-6da7-4daa-957e-bdf187d1aa29
# ╠═0022e1db-f5e4-47f1-996e-ff12a3022e5a
# ╠═3c88f1f0-753c-4182-8458-ac97a7e375f3
# ╟─877ca78d-d8c9-4a42-9332-a5c4979ac822
# ╠═b2195cfe-043a-430f-9ee0-d6ac7f745d71
# ╠═1791ff51-f819-420a-bf03-820fc48313d7
# ╠═a2843d0e-b520-4ab8-81ac-cf8b1ecf1d40
# ╟─bad0aa36-6cc1-49ff-ab3e-9db894055c64
# ╠═80f29cd8-0c01-444c-8ea3-8f6655b99d05
# ╟─68919e5f-18da-4bbf-ad66-2720dd556f32
# ╠═ce737b93-8409-4f2a-895c-e6dee7dcf570
# ╟─f8719bda-5f64-4b2a-a5c3-93343116c0ff
# ╠═8c7f0eb9-8db3-4e7d-a27e-a38d08f530b8
# ╠═fcadec48-a6cd-4471-95b2-4da9327eac37
# ╟─8e8ac21a-7310-445e-b515-d6e718f2877e
# ╠═78fab451-9329-41ea-b1df-97c57cb6c0a4
# ╟─3370dabd-30ba-4233-a60d-0c563ffd82fe
# ╠═f3a89703-2903-4874-9a57-5f84e401c970
# ╠═4855b2d0-3084-437a-9bff-d56bf5e4eef9
# ╠═40ca06d9-6845-4af0-b5ed-4e3827156914
# ╟─3804d33f-4908-4ff8-841a-d80a9f53e6ee
# ╠═7bc4bffa-3faf-42f7-acc4-4e8c28181271
# ╠═e6d23eb4-c10d-489f-8390-f344dd4b81db
# ╟─1d2f4bc8-64bd-489e-ac68-f9d5c21f11b2
# ╟─74fdc782-cba9-417a-b515-9ad8c38ad2b3
# ╠═0b9eeec1-fbb6-4408-9b42-79f20017fb74
# ╟─e40ffe89-816e-4e5e-996c-4fb2ba05d244
# ╠═40270d30-6259-4641-afee-d4b74571995c
# ╟─13e45f04-6626-46e8-9960-f63b0db6b01b
# ╠═4d83c990-81f2-4b94-8a8a-95905df93594
# ╠═08a44e3e-daac-4a07-a1d8-a1bec02c0349
# ╟─5e51ce22-ec1d-47c6-ab5c-eec254b0ec72
# ╠═6cc5f45d-a149-4603-b5f5-2b895e37d134
# ╠═ee8da3f3-ca78-433c-8510-40c0ea8958ba
# ╟─5b8de2c3-23e3-4f80-b95a-d01b3397bb2b
# ╟─d9689381-94a4-45a6-812a-8ca6e3027c99
# ╠═307933d4-af70-4fcc-a0b9-7818487e88e6
# ╠═8df9cc05-c51b-4675-9441-12f9b992222e
# ╠═56c80ef5-d974-4a7c-9b2c-adf9bf3282b5
# ╟─b6847097-a4aa-476a-b8d8-d354168553a5
# ╟─87089c7e-db1d-42d2-bdec-57e3c42801b7
# ╠═bdd4e398-0049-4b1c-8114-bb7ab786d854
# ╟─834a723b-52cf-4cc8-a815-9175be6370b6
# ╠═0407b18f-216f-47d2-8102-ad58326f9baa
# ╠═f9ec57fc-f7d1-4231-a873-0371250d48fd
# ╟─2597a206-5396-4c33-98d9-7ea6b2dbbeca
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: Curso Julia/Github - Material/src/CUDA/CUDA.jl


text:
__precompile__()

module ODEFilters

using Reexport
@reexport using DiffEqBase
using OrdinaryDiffEq
import DiffEqBase: check_error!

import Base: copy, copy!, show, size, ndims
stack(x) = copy(reduce(hcat, x)')

using LinearAlgebra
using TaylorSeries
@reexport using StructArrays
using UnPack
using RecipesBase
using ModelingToolkit
using RecursiveArrayTools
using StaticArrays

@reexport using PSDMatrices
import PSDMatrices: X_A_Xt
X_A_Xt(A, X) = X*A*X'
apply_diffusion(Q, diffusion::Diagonal) = X_A_Xt(Q, sqrt.(diffusion))
apply_diffusion(Q, diffusion::Number) = Q*diffusion

# All the Gaussian things
@reexport using GaussianDistributions
using GaussianDistributions: logpdf
const PSDGaussian{T} = Gaussian{Vector{T}, PSDMatrix{T}}
const PSDGaussianList{T} = StructArray{PSDGaussian{T}}
copy(P::Gaussian) = Gaussian(copy(P.μ), copy(P.Σ))
copy!(dst::Gaussian, src::Gaussian) = (copy!(dst.μ, src.μ); copy!(dst.Σ, src.Σ); nothing)
RecursiveArrayTools.recursivecopy(P::Gaussian) = copy(P)
show(io::IO, g::Gaussian) = print(io, "Gaussian($(g.μ), $(g.Σ))")
show(io::IO, ::MIME"text/plain", g::Gaussian{T, S}) where {T, S} =
    print(io, "Gaussian{$T,$S}($(g.μ), $(g.Σ))")
size(g::Gaussian) = size(g.μ)
ndims(g::Gaussian) = ndims(g.μ)

Base.:*(M, g::PSDGaussian) = Gaussian(M * g.μ, X_A_Xt(g.Σ, M))
GaussianDistributions.whiten(Σ::PSDMatrix, z) = Σ.L\z

import Statistics: mean, var, std
var(p::PSDGaussian{T}) where {T} = diag(p.Σ)
std(p::PSDGaussian{T}) where {T} = sqrt.(diag(p.Σ))
mean(s::PSDGaussianList{T}) where {T} = mean.(s)
var(s::PSDGaussianList{T}) where {T} = var.(s)
std(s::PSDGaussianList{T}) where {T} = std.(s)

include("priors.jl")
include("diffusions.jl")

include("algorithms.jl")
export EK0, EK1
include("alg_utils.jl")
include("caches.jl")
include("state_initialization.jl")
include("integrator_utils.jl")
include("filtering.jl")
include("perform_step.jl")
include("preconditioning.jl")
include("smoothing.jl")
include("solution.jl")

# Utils
include("jacobian.jl")
include("numerics_tricks.jl")

# Iterated Extended Kalman Smoother
include("ieks.jl")
export IEKS, solve_ieks

end

ext: jl


text:
module test_spineruntimes_lib

using Test
using SpineRuntimes: Lib

end # module test_spineruntimes_lib


text:
"""

    calNewMLD(; h, B, fric_u, Δb, f, Δt, m=0.8, n=0.2)

# Description

 The calculation follows eqn (30) of Gasper (1988). The changes are:

   1. We assume that all sunlight is absorbed at the very surface. According to Paulson and James (1988), in upper ocean
      90% of the radiation is absorbed within the first 10 meters. Since even POP has the resultion of 10m, it is reasonable
      that we ignore the structrue of light absorbtion in sea-water.

   2. Turbulent kinetic energy is added into the equation to avoid the case Δb = 0 which results in we → ∞. This follows
      the reasoning in Alexander et al. (2000)

 References:

   Gaspar, Philippe. "Modeling the seasonal cycle of the upper ocean." Journal of Physical Oceanography 18.2 (1988): 161-180.

   Paulson, Clayton A., and James J. Simpson. "Irradiance measurements in the upper ocean."
   Journal of Physical Oceanography 7.6 (1977): 952-956.

   Alexander, Michael A., James D. Scott, and Clara Deser. "Processes that influence sea surface temperature and ocean mixed
   layer depth variability in a coupled model." Journal of Geophysical Research: Oceans 105.C7 (2000): 16823-16842.

# Return values
This function returns a list with two elements. The first is a symbol. ``:we`` indicates the second value is the entrainment speed whereas ``:MLD`` indicates the second value is the diagnosed MLD.

"""
function calNewMLD(;
    h_ML   :: Float64,
    Bf     :: Float64,    # sum of sensible heat fluxes and longwave radiation flux
    J0     :: Float64,    # shortwave radiation flux (converted to buoyancy flux)
    fric_u :: Float64,
    Δb     :: Float64,
    f      :: Float64,
    Δt     :: Float64,
    ζ      :: Float64,    # decay scale length of shortwave penetration
    m::Float64 = 0.8,
    n::Float64 = 0.20,
    h_max  :: Float64,
    we_max :: Float64,
)

    if Δb < 0
        throw(ErrorException("Δb cannot be negative. Right now Δb = " * string(Δb)))
    end

    #
    # Transform the equation of entrainment into the form
    # we (h_ML × Δb + k) =  a exp(-h / λ) + b h + c S(h, ζ) := Δ
    #
    # where Δ is called "determinant"
    #       a       = 2 m u_fric^3
    #       b       = Bf
    #       c       = J0
    #       S(h, γ) = h - 2 ζ + exp(-h/ζ) (h + 2 ζ)
    #       λ       = fric_u / abs(f)
    #

    a, λ = calΔCoefficients(fric_u, f, m)
    Δ, _ = calΔ_and_∂Δ∂h(h_ML, a, Bf, J0, λ, ζ, n; need_∂Δ∂h=false)

    h_MO = solveMoninObuhkovLength_bisection(h_max, a, Bf, J0, λ, ζ, n)

    # 2019/08/05 Newton method. Not a robust way
    # h_MO = solveMoninObuhkovLength(h_max, a, Bf, J0, λ, ζ, n)

    if Δ >= 0
        k = getTKE(fric_u=fric_u)
        we = Δ / (h_ML * Δb + k)

        # First min function restrict the maximum of entrainment speed
        # Second min function restrict the maximum of MLD

        return min(h_ML + Δt * min(we, we_max), h_MO), h_MO
    else

        # h becomes diagnostic.
        # Notice that if Δ < 0, then finite positive solution of h_MO exists
        return h_MO, h_MO
    end

end

@inline function calΔCoefficients(
    u_fric::Float64,
    f::Float64,
    m::Float64,
)
    return 2.0 * m * u_fric^3.0, u_fric / abs(f)
end

function calΔ_and_∂Δ∂h(
    h::Float64,
    a::Float64,
    b::Float64,
    c::Float64,
    λ::Float64,
    ζ::Float64,
    n::Float64;
    need_∂Δ∂h=true,
)
    Bh = b * h + c * S(h, ζ)
    ( Bh < 0.0 ) && ( n = 1.0 )
    exponential = exp(-h / λ)
    return ( ( λ == 0.0 ) ? 0.0 : a * exponential ) + n * Bh,
           ( need_∂Δ∂h ) ? ( ( λ == 0.0 ) ? 0.0 : - a / λ * exponential ) + n * (b + c * ∂S∂h(h, ζ)) : 0.0
end

@inline function S(h::Float64, ζ::Float64)
    return h - 2.0 * ζ + exp(-h / ζ) * (h + 2.0 * ζ)
end

@inline function ∂S∂h(h::Float64, ζ::Float64)
    return 1 - exp(-h / ζ) * (1 + h / ζ)
end

# it is possible to have multisolution but we ignore this problem
function solveMoninObuhkovLength_bisection(
    h_max  :: Float64,
    a      :: Float64,
    b      :: Float64,
    c      :: Float64,
    λ      :: Float64,
    ζ      :: Float64,
    n      :: Float64;
    δh_max :: Float64 = 0.01,  # absolute increment threshold = 1cm
    max    :: Integer = 100,
)

    Δ, _ = calΔ_and_∂Δ∂h(h_max, a, b, c, λ, ζ, n; need_∂Δ∂h = false)

    if Δ > 0
        return h_max
    end

    bot = 0.0
    top = h_max

    local h

    for t = 1:max

        h = (bot + top) / 2.0
        Δ, _ = calΔ_and_∂Δ∂h(h, a, b, c, λ, ζ, n; need_∂Δ∂h = false)

        if Δ < 0
            top = h
        elseif Δ > 0
            bot = h
        else # Δ = 0 which is rare
            return h
        end

        if (top - bot) < δh_max
            return h
        end
    end

    throw(ErrorException("Monin-Obuhkov length iteration reached max before δh gets below δh_max."))
end

function solveMoninObuhkovLength(
    h      :: Float64,   # initial guess
    a      :: Float64,
    b      :: Float64,
    c      :: Float64,
    λ      :: Float64,
    ζ      :: Float64,
    n      :: Float64;
    η      :: Float64 = 0.01,  # relative increment threshold = 1%
    δh_max :: Float64 = 0.01,  # absolute increment threshold = 1cm
    max    :: Integer = 100,
)

    if_converge = false
    prev_δh = 0.0

    # The only contribution of negativity comes from b + c
    if b + c > 0
#        println("MO Length Inf!", h, ",", a, ",", b, ",", c, ",", λ)
        return Inf
    end

    for i=1:max
        Δ, ∂Δ∂h = calΔ_and_∂Δ∂h(h, a, b, c, λ, ζ, n)
        δh = - Δ / ∂Δ∂h
        #println(i, ":", h, "; δh=",δh, ", Δ=", Δ, "; ∂Δ∂h=", ∂Δ∂h)
        if abs((δh - prev_δh) / prev_δh) >= η && abs(δh) >= δh_max
            h, prev_δh = h + δh, δh
        else
            if_converge = true
#            if h < 10.0
#                println("Converge!", h, ",", a, ",", b, ",", c, ",", λ)
#            end
            break
        end
    end

    if if_converge
        return h
    else
        throw(ErrorException("Monin-Obuhkov length iteration cannot converge"))
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Flux
using StatsBase

struct CounterfactualData
    X::AbstractMatrix
    y::AbstractMatrix
    mutability::Union{Vector{Symbol},Nothing}
    domain::Union{Any,Nothing}
    categorical::Union{Vector{Int},Nothing}
    continuous::Union{Vector{Int},Nothing}
    standardize::Bool
    dt::StatsBase.AbstractDataTransform
    function CounterfactualData(X,y,mutability,domain,categorical,continuous,standardize,dt)
        conditions = []
        conditions = vcat(conditions..., length(size(X)) != 2 ? error("Data should be in tabular format") : true)
        conditions = vcat(conditions..., size(X)[2] != size(y)[2] ? throw(DimensionMismatch("Number of output observations is $(size(y)[2]). Expected: $(size(X)[2])")) : true)
        if all(conditions)
            new(X,y,mutability,domain,categorical,continuous,standardize,dt)
        end
    end
end

"""
    CounterfactualData(
        X::AbstractMatrix, y::AbstractMatrix;
        mutability::Union{Vector{Symbol},Nothing}=nothing,
        domain::Union{Any,Nothing}=nothing,
        categorical::Union{Vector{Int},Nothing}=nothing,
        continuous::Union{Vector{Int},Nothing}=nothing,
        standardize::Bool=false
    )

This outer constructor method prepares features `X` and labels `y` to be used with the package. Mutability and domain constraints can be added for the features. The function also accepts arguments that specify which features are categorical and which are continues. These arguments are currently not used.

# Examples

```julia-repl
using CounterfactualExplanations.Data
x, y = toy_data_linear()
X = hcat(x...)
counterfactual_data = CounterfactualData(X,y')
```

"""
function CounterfactualData(
    X::AbstractMatrix, y::AbstractMatrix;
    mutability::Union{Vector{Symbol},Nothing}=nothing,
    domain::Union{Any,Nothing}=nothing,
    categorical::Union{Vector{Int},Nothing}=nothing,
    continuous::Union{Vector{Int},Nothing}=nothing,
    standardize::Bool=false
)

    # If nothing supplied, assume all continuous:
    if isnothing(categorical) && isnothing(continuous)
        continuous = 1:size(X)[1]
    end

    # If tuple is supplied, assume it counts for all continuous variables:
    if typeof(domain) <: Tuple
        domain = [domain for i in 1:size(X)[1]]
    end

    # Data transformer:
    dt = fit(ZScoreTransform, X, dims=2)

    counterfactual_data = CounterfactualData(X, y, mutability, domain, categorical, continuous, standardize, dt)

    return counterfactual_data
end

"""
    select_factual(counterfactual_data::CounterfactualData, index::Int)

A convenience method that can be used to access the the feature matrix.
"""
select_factual(counterfactual_data::CounterfactualData, index::Int) = counterfactual_data.X[:,index]

"""
    mutability_constraints(counterfactual_data::CounterfactualData)

A convience function that returns the mutability constraints. If none were specified, it is assumed that all features are mutable in `:both` directions.
"""
mutability_constraints(counterfactual_data::CounterfactualData) = isnothing(counterfactual_data.mutability) ? [:both for i in 1:size(counterfactual_data.X)[1]] : counterfactual_data.mutability

"""
    apply_domain_constraints(counterfactual_data::CounterfactualData, x::AbstractArray)

A subroutine that is used to apply the predetermined domain constraints.
"""
function apply_domain_constraints(counterfactual_data::CounterfactualData, x::AbstractArray)

    # Continuous variables:
    if !isnothing(counterfactual_data.domain)
        for i in counterfactual_data.continuous
            x[i] = clamp(x[i], counterfactual_data.domain[i][1], counterfactual_data.domain[i][2])
        end
    end

    return x

end


text:
using Test
@testset "GpABC Full Test Suite" begin
include("kernels/abstract_kernel_tests.jl")
include("kernels/scaled_squared_distance_test.jl")
include("kernels/rbf_kernel_tests.jl")
include("kernels/matern_kernel_test.jl")
include("gp_regression/gp_regression_tests.jl")
include("abc/rejection-abc-test.jl")
include("abc/smc-abc-test.jl")
include("abc/model-selection-test.jl")
include("lna-test.jl")
end


text:
using Performance
using Test

ntests_finished = 5   # increment this counter to run the next set of tests

@testset "Performance.jl" begin
    ntests = ntests_finished
    if ntests > 0
        ntests -= 1
        include("globals.jl")
    end
    if ntests > 0
        ntests -= 1
        include("arrays.jl")
    end
    if ntests > 0
        ntests -= 1
        include("structs.jl")
    end
    if ntests > 0
        ntests -= 1
        include("type_abuse.jl")
    end
    if ntests > 0
        ntests -= 1
        include("algorithm_order.jl")
    end
end


text:
using Test

new_floor = function(direction, x)
    if direction == '('
        return x+1
    elseif direction == ')'
        return x-1
    else
        return x
    end
end

day1_part1 = function(x)
    floor = 0
    for c in x
        floor = new_floor(c, floor)
    end
    return floor
end

@testset "day1_part1" begin
    @test 0 == day1_part1("")
    @test 1 == day1_part1("(")
    @test 2 == day1_part1("((")
    @test 1 == day1_part1("(()")
    @test 0 == day1_part1("(())")
    @test 0 == day1_part1("()()")
    @test 3 == day1_part1("(((")
    @test 3 == day1_part1("(()(()(")
    @test 3 == day1_part1("))(((((")
    @test -1 == day1_part1("())")
    @test -1 == day1_part1("))(")
    @test -3 == day1_part1(")))")
    @test -3 == day1_part1(")())())")
    @test 0 == day1_part1("()\n")
end

day1_part2 = function(x)
    str = ""
    for c in x
        str *= string(c)
        if -1 == day1_part1(str)
            break
        end
    end
    return length(str)
end

@testset "day1_part2" begin
    @test 1 == day1_part2(")")
    @test 5 == day1_part2("()())")
end

const day1_input = read("1.txt", String)
println("""Day 1, Part 1: $(day1_part1(day1_input))""")
println("""Day 1, Part 2: $(day1_part2(day1_input))""")


text:
@testset "history.jl" begin

        @testset ":tally" begin
            global tallied = Variable{Int64}(data=:tally, history=true)

            @process test_process(value_durations) begin
                for (value, duration) in value_durations
                    tallied.value = value
                    work(duration)
                end
            end

            @simulation begin
                value_durations = [(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)]

                @schedule now test_process(value_durations)
                @schedule at 1.0 test_process(value_durations)
                @schedule in 3.0 test_process(value_durations)

                start_simulation()
            end

            want = [1, 1, 2, 1, 2, 3, 3, 2, 4, 4, 3, 4]
            @test tallied.history.data == want
        end

        @testset ":accumulate" begin
            global accumulated = Variable{Int64}(data=:accumulate, history=true)

            @process test_process(value_durations) begin
                for (value, duration) in value_durations
                    accumulated.value = value
                    work(duration)
                end
            end

            @simulation begin
                value_durations = [(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)]

                @schedule now test_process(value_durations)
                @schedule at 1.0 test_process(value_durations)
                @schedule in 3.0 test_process(value_durations)

                start_simulation()
            end

            want = [1, 2, 1, 2, 3, 2, 4, 3, 4]
            @test accumulated.history.data == want

            want = [2.0, 1.0, 0.0, 0.0, 2.0, 0.0, 1.0, 2.0, 3.0]
            @test accumulated.history.durations == want
        end
end


text:
using ThermodynamicIntegration
using Documenter

makedocs(;
    modules=[ThermodynamicIntegration],
    authors="Theo Galy-Fajou <theo.galyfajou@gmail.com> and contributors",
    repo="https://github.com/theogf/ThermodynamicIntegration.jl/blob/{commit}{path}#L{line}",
    sitename="ThermodynamicIntegration.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://theogf.github.io/ThermodynamicIntegration.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/theogf/ThermodynamicIntegration.jl",
)


text:
using Plots

include("QuadraturePoints.jl")
#include("element_matrices.jl")
include("create_grid.jl")
#include("global_matrices.jl")
include("exact_solution.jl")
#include("compute_norms.jl")
#include("time_integration.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

#ipoints=1: Lobatto
#ipoints=2: Legendre
#ipoints=3: Chebyshev
#ipoints=4: Equi-spaced

#{{{ Main
function main()

    #Fix Interpolation and Integration Order
    N=3
    Q=N+1
    Np=N+1
    Nq=Q+1
    ipoints=1
    qpoints=1
    space_method="CG"
    #space_method="DG"
    #ti_type="SSP"
    ti_type="LSRK"
    case=2 #1=exponential, 2=square wave
    stages=3 #1,2, or 3 only
    if (ti_type == "LSRK")
        stages=5
    end
    Ne=40
    Npoin_dg=Ne*Np
    Npoin_cg=Ne*N + 1
    dt=DFloat(0.0001)
    time_final=DFloat(2.0)
    Courant_max=0.5
    ntime = ceil(Int64, time_final / dt)
    @show(N,Q,Ne,Npoin_cg,Npoin_dg,dt,ntime,case,ti_type,stages)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq,DFloat)

    #Construct Grid and (i,e)->I map
    (coord_cg,coord_dg,intma_cg,intma_dg,periodicity_cg,periodicity_dg) = create_grid(Np,Npoin_cg,Npoin_dg,Ne,ξ,DFloat)

    #Construct Local Mass and Differentiation Matrices
#    (Me,De,Dwe,Fe) = element_matrices(ψ,dψ,Np,Nq,ωq,DFloat)

    #CG or DG?
    if (space_method == "CG")
        Npoin=Npoin_cg
        coord=copy(coord_cg)
        intma=copy(intma_cg)
        periodicity=copy(periodicity_cg)
    elseif (space_method == "DG")
        Npoin=Npoin_dg
        coord=copy(coord_dg)
        intma=copy(intma_dg)
        periodicity=copy(periodicity_dg)
    end
    @show(space_method,Npoin,Ne,N,Q)

    #Construct Global Mass and Differentiation Matrices
#    (M,D,Dw,F) = global_matrices(Me,De,Dwe,Fe,intma,coord,Npoin,Ne,Np,periodicity,DFloat)

    #Compute Initial Condition
    time=DFloat(0.0)
    (qe,u) = exact_solution(coord,Npoin,time,case,DFloat)
    q0=copy(qe)
    dx_min=coord[2]-coord[1]

    println("Done") #output

end
#}}} Main

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
import GitLab
using Base.Test

myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

myrepo = GitLab.repo_by_name("TestProject1"; headers=options)
## @show myrepo

@test get(GitLab.namefield(myrepo)) == "TestProject1"

try ## remove this later
    forks = GitLab.forks(myrepo; params=options)
    @show forks
end

yourrepo = GitLab.repo_by_name("admin-project"; headers=options)
## @show yourrepo

try ## This may fail if the fork already exists
    forked_repo = GitLab.create_fork(yourrepo; headers=options)
    @show forked_repo

    repo = GitLab.delete_fork(forked_repo; headers=options)
    @show repo
end

contributors = GitLab.contributors(myrepo, params=options)
@test GitLab.name(first(contributors)[1]["contributor"]) == "Pradeep Mudlapur"
## @show contributors

collaborators = GitLab.collaborators(myrepo, params=options)
@test GitLab.name(first(collaborators)) == "mdpradeep"
## @show collaborators

collaborator = GitLab.iscollaborator(myrepo, "mdpradeep", params=options)
@test collaborator == true
## @show collaborator

non_collaborator = GitLab.iscollaborator(myrepo, "XXXXX", params=options)
@test non_collaborator == false
## @show non_collaborator

user_data = Dict{AbstractString, Any}()
user_data["user_id"] = 3
user_data["access_level"] = 40
result = GitLab.add_collaborator(myrepo, "3", headers=options, params=user_data)
## result = GitLab.add_collaborator(myrepo, "", headers=options, params=user_data)
@test result.status <= 204

no_user_data = Dict{AbstractString, Any}()
no_user_data["user_id"] = 999
no_user_data["access_level"] = 40
try
    result = GitLab.add_collaborator(myrepo, "999", headers=options, params=no_user_data)
    @test result.status <= 204
catch e
    println("Failed (expected) to added collaborator $(e)")
    ## This is the expected behavior in this case
end

result = GitLab.remove_collaborator(myrepo, "3", headers=options, params=user_data)
@test result.status == 200
## @show result

try
    result = GitLab.remove_collaborator(myrepo, "999", headers=options, params=user_data)
    @test result.status == 200
    ## @show result
catch e
    println("Failed (expected) to remove collaborator $(e)")
    ## This is the expected behavior in this case
end

println("Repos Tests Done !!!")


text:
"""
Subroutine dstrsol

This subroutine solves the triangular systems L*x = r or L'*x = r.

MINPACK-2 Project. May 1998.
Argonne National Laboratory.
"""
function dstrsol(n, L::TronSparseMatrixCSC, r, task)
    zero = 0.0

    # Solve L*x = r and store the result in r.

    if task == 'N'
        @inbounds for j=1:n
            temp = r[j]/L.diag_vals[j]
            for k=L.colptr[j]:L.colptr[j+1]-1
                r[L.rowval[k]] = r[L.rowval[k]] - L.tril_vals[k]*temp
            end
            r[j] = temp
        end

        return
    end

    # Solve L'*x = r and store the result in r.

    if task == 'T'
        r[n] = r[n]/L.diag_vals[n]
        @inbounds for j=n-1:-1:1
            temp = zero
            for k=L.colptr[j]:L.colptr[j+1]-1
                temp = temp + L.tril_vals[k]*r[L.rowval[k]]
            end
            r[j] = (r[j] - temp)/L.diag_vals[j]
        end

        return
    end
end

function dstrsol(n, L::TronDenseMatrix, r, task)
    zero = 0.0

    # Solve L*x = r and store the result in r.

    if task == 'N'
        @inbounds for j=1:n
            temp = r[j]/L.vals[j,j]
            for k=j+1:n
                r[k] = r[k] - L.vals[k,j]*temp
            end
            r[j] = temp
        end

        return
    end

    # Solve L'*x = r and store the result in r.

    if task == 'T'
        r[n] = r[n]/L.vals[n,n]
        @inbounds for j=n-1:-1:1
            temp = zero
            for k=j+1:n
                temp = temp + L.vals[k,j]*r[k]
            end
            r[j] = (r[j] - temp)/L.vals[j,j]
        end

        return
    end
end


text:
# This file describes the methods that should be implemented for each
# new Koala transformer type. See https://github.com/ablaom/KoalaTransforms.jl
# for examples.

mutable struct SomeTransformerType <: Transformer
    param1::Float64
    param2::Int
    etc...
end

# lazy keyword constructor:
SomeTransformerType(; param1=default1, param2=default2, etc...) =
    SomeTransformerType(param1, param2, etc...)

fit(transformer::SomeTransformerType, X, parallel, verbosity;
    optional_keyword_args...) -> scheme

transform(transformer::SomeTransformerType, scheme, X) -> Xt

# optional:
inverse_transform(transformer::SomeTransformerType, scheme, Xt) -> Xt


text:
function data_driven_lyapb_quad(state0,state;horizon=1,C=1e3,ub=1e2,lb=0,tol=1e-4,numIter=1e2,postprocess=false)
    numTraj = size(state0)[2]
    dim = size(state0)[1]
    iter = 1
    gammaU = ub
    gammaL = lb
    while gammaU-gammaL > tol && iter < numIter
        iter += 1
        gamma = (gammaU + gammaL)/2
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, P[1:dim, 1:dim] in PSDCone())
        @SDconstraint(model, P >= Matrix(I,dim,dim))
        @objective(model, Max, 0)
        for i in 1:numTraj
          @constraint(model, state[:,i]'*P*state[:,i] <= gamma^(2*horizon)*state0[:,i]'*P*state0[:,i])
        end
        @SDconstraint(model, P <= C*Matrix(I,dim,dim))
        JuMP.optimize!(model)
        if termination_status(model) == MOI.OPTIMAL
          gammaU=gamma
        else
          gammaL=gamma
        end
    end
    #=
    while gammaU-gammaL > tol && iter < numIter
        iter += 1
        gamma = (gammaU + gammaL)/2
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, P[1:dim, 1:dim] in PSDCone())
        @variable(model, s>=0)
        @SDconstraint(model, P >= Matrix(I,dim,dim))
        @objective(model, Min, s)
        for i in 1:numTraj
          @constraint(model, state[:,i]'*P*state[:,i] <= gamma^(2*horizon)*state0[:,i]'*P*state0[:,i]+s)
        end
        @SDconstraint(model, P <= C*Matrix(I,dim,dim))
        JuMP.optimize!(model)
        if value.(s) < 1e-8
          gammaU=gamma
        else
          gammaL=gamma
        end
    end=#
    gamma = gammaU

    if postprocess == true
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, P[1:dim, 1:dim] in PSDCone())
        @SDconstraint(model, P >= Matrix(I,dim,dim))
        # @variable(model,t >= 0)
        # @constraint(model,t >= norm(P, 2))
        @objective(model, Min, sum(P[:].^2))
        for i in 1:numTraj
          @constraint(model, state[:,i]'*P*state[:,i] <= gamma^(2*horizon)*state0[:,i]'*P*state0[:,i])
        end
        @SDconstraint(model, P <= C*Matrix(I,dim,dim))
        JuMP.optimize!(model)
        return gamma, value.(P)
    else
        return gamma
    end

    end

function soslyap_alternating(state0,state;B,d,batchsize,K0,tol=1e-4)
    (n,m) = size(B)
    nlift = binomial(n+d-1, d)

    #transform = kronecker2veronese(n,d)
    Xscale = veroneseliftscale(n,d)

    P0 = Matrix(1.0I, nlift, nlift)
    ga0 = convergerate(state0,state,B,d,P0,K0,Xscale)
    gatemp = copy(ga0)
    Ktemp = copy(K0)
    Ptemp = copy(P0)

    dga = 1e2

    Xscale = veroneseliftscale(n,d)

    while dga > tol
        println(repeat('*', 80))
        println("Order: $d || Convergence rate: $gatemp")
        println("Minimizing P...")
        #P,ga = soslyap_alter_P(state0,state,B,d,Ktemp,gatemp)
        P_update,ga_update = soslyap_alternating_P_distributed(state0,state,batchsize,B,d,Ptemp,Ktemp,gatemp,Xscale)
        cholP = cholesky(P_update)
        println("Minimizing K...")
        K_updated,ga = soslyap_alternating_K_distributed(state0,state,batchsize,B,d,P_update,cholP,Ktemp,ga_update,Xscale)

        dga = gatemp-ga
        gatemp = copy(ga)
        Ktemp = copy(K_updated)
        Ptemp = copy(P_update)
    end
    println(repeat('*', 80))
    println("Order: $d || Convergence rate (terminated): $gatemp")

    return Ktemp, gatemp, Ptemp

end

function soslyap_alternating_K(state0,state,B,d,P,cholP,K0,ga0,Xscale)

    (n,m) = size(B)
    nlift = binomial(n+d-1, d)
    N = size(state0,2)

    V0 = zeros(N)
    for i in 1:N
        xlift0_i = veroneselift(state0[:,i],Xscale,d)
        V0[i] = transpose(xlift0_i)*P*xlift0_i
    end

    if d == 1
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, K[1:m, 1:n])
        @variable(model, ga >=0)
        @objective(model, Min, ga)
        for i in 1:N
            @SDconstraint(model, [ga*V0[i] (state[:,i]+B*K*state0[:,i])'*P;P*(state[:,i]+B*K*state0[:,i]) P] >=0)
        end
        JuMP.optimize!(model)
        return reshape(JuMP.value.(K),m,n),(sqrt(JuMP.value.(ga)))^(1/d)
    else
        solver = optimizer_with_attributes(Ipopt.Optimizer, MOI.Silent() => true)
        #solver = optimizer_with_attributes(NLopt.Optimizer)
        model = Model(solver)
        #set_optimizer_attribute(model, "algorithm", :GN_ISRES) #:GN_ISRES GN_AGS GN_ORIG_DIRECT
        #set_optimizer_attribute(model, "local_optimizer", :LD_MMA) #LD_SLSQP LD_MMA LN_COBYLA
        @variable(model, -1e1 <= K[i=1:n*m] <= 1e1)
        set_start_value.(K,reshape(K0,m*n,1))
        @variable(model, 0 <= ga <= ga0^(2*d),start = ga0^(2*d))
        @objective(model, Min, ga)

        lyapconstraint = []
        for i in 1:N
            push!(lyapconstraint,(y,x...) -> sum(sum(cholP.U[j,k]*veroneselift(state[:,i]+B*reshape(collect(x),m,n)*state0[:,i],Xscale,d)[k] for k in 1:nlift)^2 for j in 1:nlift)-y*V0[i])
        end
        for (i, f) in enumerate(lyapconstraint)
            f_sym = Symbol("f_$(i)")
            register(model, f_sym, n*m+1, f; autodiff = true)
            add_NL_constraint(model, :($(f_sym)($(ga),$(K...)) <= 0))
        end
        JuMP.optimize!(model)
        return reshape(JuMP.value.(K),m,n),(sqrt(JuMP.value.(ga)))^(1/d)

    end
end

function soslyap_alternating_K_distributed(state0,state,batchsize,B,d,P,cholP,K0,ga0,Xscale)

    N = size(state0,2)
    batchnum = ceil(Int,N/batchsize)
    partition = []
    for i in 1:batchnum
        if i<batchnum
            push!(partition,(i-1)*batchsize+1:i*batchsize)
        else
            push!(partition,(batchnum-1)*batchsize+1:N)
        end
    end

    gammaopt = copy(ga0)
    Kopt = copy(K0)

    println("K: Total number of batches: $batchnum")
    for i in 1:batchnum
        println("K: Solving batch $i...")
        K_i,gamma_i = soslyap_alternating_K(state0[:,partition[i]],state[:,partition[i]],B,d,P,cholP,K0,ga0,Xscale)
        Kopt_i,gammaopt_i = backtrackinglinesearch_K(state0,state,B,d,P,K_i,K0,Xscale)
        if gammaopt_i < gammaopt
            Kopt = copy(Kopt_i)
            gammaopt = copy(gammaopt_i)
        end
    end
    return Kopt,gammaopt
end

function backtrackinglinesearch_K(state0,state,B,d,P,Ktemp,K0,Xscale,tol=1e-2)

    gammaopt = 1e6
    Kopt = copy(K0)
    for lambda in 0:tol:1
        Klambda = (1-lambda)*K0+lambda*Ktemp
        galambda = convergerate(state0,state,B,d,P,Klambda,Xscale)
        #println(galambda)
        if galambda <= gammaopt
            gammaopt = copy(galambda)
            Kopt = copy(Klambda)
        end
    end

    return Kopt,gammaopt
end

function soslyap_alternating_P(state0,state,B,d,K,ga0,Xscale)

    n = size(B,1)
    nlift = binomial(n+d-1, d)
    N = size(state0,2)

    statelift0 = zeros(nlift,N)
    statelift = zeros(nlift,N)
    for i in 1:N
        statelift0[:,i] = veroneselift(state0[:,i],Xscale,d)
        statelift[:,i] = veroneselift(state[:,i]+B*K*state0[:,i],Xscale,d)
    end
    if d >= 2
        ga,P =  data_driven_lyapb_quad(statelift0,statelift;C=1e2,ub=ga0^d,postprocess=true)
    else
        ga,P =  data_driven_lyapb_quad(statelift0,statelift;C=1e2,ub=ga0^d,postprocess=true)
    end

    return P,ga^(1/d)
end

function soslyap_alternating_P_distributed(state0,state,batchsize,B,d,P0,K0,ga0,Xscale)
    N = size(state0,2)
    batchnum = ceil(Int,N/batchsize)
    partition = []
    for i in 1:batchnum
        if i<batchnum
            push!(partition,(i-1)*batchsize+1:i*batchsize)
        else
            push!(partition,(batchnum-1)*batchsize+1:N)
        end
    end

    gammaopt = copy(ga0)
    Popt = copy(P0)

    println("P: Total number of batches: $batchnum")
    for i in 1:batchnum
        println("P: Solving batch $i...")
        P_i, = soslyap_alternating_P(state0[:,partition[i]],state[:,partition[i]],B,d,K0,ga0,Xscale)
        Popt_i,gammaopt_i = backtrackinglinesearch_P(state0,state,B,d,P_i,P0,K0,Xscale)
        if gammaopt_i < gammaopt
            Popt = copy(Popt_i)
            gammaopt = copy(gammaopt_i)
        end
    end
    return Popt,gammaopt
end

function backtrackinglinesearch_P(state0,state,B,d,Ptemp,P0,K0,Xscale,tol=1e-2)

    gammaopt = 1e6
    Popt = copy(P0)

    for lambda in 0:tol:1
        Plambda = (1-lambda)*P0+lambda*Ptemp
        galambda = convergerate(state0,state,B,d,Plambda,K0,Xscale)
        if galambda <= gammaopt
            gammaopt = copy(galambda)
            Popt = copy(Plambda)
        end
    end

    return Popt,gammaopt
end

function convergerate(state0,state,B,d,P,K,Xscale)
    N = size(state0,2)
    gamma = 0.0
    for i in 1:N
        xlift_i = veroneselift(state[:,i]+B*K*state0[:,i],Xscale,d)
        xlift0_i = veroneselift(state0[:,i],Xscale,d)
        gamma_i = sqrt(transpose(xlift_i)*P*xlift_i/(transpose(xlift0_i)*P*xlift0_i))
        if gamma_i > gamma
            gamma = copy(gamma_i)
        end
    end
    return gamma^(1/d)
end

function LQR_alternating(state0,state;B,Q,R,kappaU,lowbound,K0,ep,tol=1e-4)

    (dim,N) = size(state0)
    xitemp = 0
    Ptemp = (Q+K0'*R*K0)*2
    for i in 1:N
        x_i = state[:,i]
        x0_i = state0[:,i]
        xi_i = sqrt((x_i+B*K0*x0_i)'*Ptemp*(x_i+B*K0*x0_i)/(x0_i'*Q*x0_i))
        if xitemp < xi_i
            xitemp = xi_i
        end
    end

    Ptemp,xitemp = LQR_alternating_xi(state0,state;B=B,Q=Q,R=R,kappaU=kappaU,Ktemp=K0,xitemp=xitemp)
    Ktemp = K0
    Ztemp = Ptemp-Q-Ktemp'*R*Ktemp
    eigZtemp = eigvals(Ztemp)
    ratio = 1-(maximum(eigZtemp)/minimum(eigZtemp))*(1-cos(delta2theta(ep,dim)))
    if ratio <0
        ratio = 0
    end
    #println("initializaiton: $Ptemp || $eigZtemp")

    dga = 1e2
    feasible = xitemp/ratio
    while dga > tol && feasible > 1
        println(repeat('*', 80))
        println("Contraction rate: $xitemp")
        println("Minimizing K...")
        K_update,xi_update = LQR_alternating_K(state0,state;B=B,Q=Q,R=R,kappaU=kappaU,lowbound=lowbound,Ptemp=Ptemp,Ktemp=Ktemp,xitemp=xitemp)

        println("Minimizing P...")
        P_update = LQR_alternating_P(state0,state;B=B,Q=Q,R=R,kappaU=kappaU,Ktemp=K_update,xitemp=xi_update)

        dga = xitemp-xi_update
        xitemp = xi_update
        Ptemp = P_update
        Ktemp = K_update

        Ztemp = Ptemp-Q-Ktemp'*R*Ktemp
        eigZtemp = eigvals(Ztemp)

        ratio = 1-(maximum(eigZtemp)/minimum(eigZtemp))*(1-cos(delta2theta(ep,dim)))
        if ratio <0
            ratio = 0
        end
        feasible = xitemp/ratio
    end
    println(repeat('*', 80))
    println("Contraction rate: $xitemp")

    return Ktemp, xitemp, Ptemp

end

function LQR_alternating_xi(state0,state;B,Q,R,kappaU,Ktemp,xitemp,tol=1e-4)

    numTraj = size(state0)[2]
    dim = size(state0)[1]
    ub = xitemp
    lb = 1
    if ub < lb
        ub = 1
    end
    stateK = zeros(dim,numTraj)

    for i in 1:numTraj
        stateK[:,i] = state[:,i]+B*Ktemp*state0[:,i]
    end

    while ub-lb > tol
        xi = (ub + lb)/2
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, P[1:dim, 1:dim] in PSDCone())
        @variable(model, v>=1e-6)
        @objective(model, Min, 0)
        for i in 1:numTraj
          @constraint(model, stateK[:,i]'*P*stateK[:,i] <= xi^2*state0[:,i]'*(P-Q-Ktemp'*R*Ktemp)*state0[:,i])
        end
        @SDconstraint(model, P >= Q+Ktemp'*R*Ktemp)
        @SDconstraint(model, P >= Q+Ktemp'*R*Ktemp+v*Matrix(I,dim,dim))
        @SDconstraint(model, P - Q - Ktemp'*R*Ktemp <=kappaU*v*Matrix(I,dim,dim))

        JuMP.optimize!(model)
        if termination_status(model) == MOI.OPTIMAL
            ub=xi
            #println(value.(v))
            #println(value.(P))
        else
            lb=xi
        end
    end

    solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
    model = Model(solver)
    @variable(model, P[1:dim, 1:dim] in PSDCone())
    @variable(model, v>=1e-6)
    @objective(model, Min, tr(P))
    for i in 1:numTraj
      @constraint(model, stateK[:,i]'*P*stateK[:,i] <= ub^2*state0[:,i]'*(P-Q-Ktemp'*R*Ktemp)*state0[:,i])
    end
    @SDconstraint(model, P >= Q+Ktemp'*R*Ktemp)
    @SDconstraint(model, P >= Q+Ktemp'*R*Ktemp+v*Matrix(I,dim,dim))
    @SDconstraint(model, P - Q - Ktemp'*R*Ktemp <=kappaU*v*Matrix(I,dim,dim))
    JuMP.optimize!(model)

    return value.(P), ub
end

function LQR_alternating_K(state0,state;B,Q,R,kappaU,lowbound,Ptemp,Ktemp,xitemp)
    numTraj = size(state0)[2]
    (n,m) = size(B)

    #=ub = xitemp^2
    lb = lowbound^2
    Kopt = zeros(m,n)

    while ub-lb > 1e-4
        xi = (ub + lb)/2
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, K[1:m, 1:n])
        @objective(model, Min, 0)
        for i in 1:numTraj
            @constraint(model, (state[:,i]+B*K*state0[:,i])'*Ptemp*(state[:,i]+B*K*state0[:,i]) <= xi*state0[:,i]'*(Ptemp-Q-K'*R*K)*state0[:,i])
        end
        @SDconstraint(model, [Ptemp-Q K';K inv(R)]>=0)
        JuMP.optimize!(model)
        if termination_status(model) == MOI.OPTIMAL
            ub=xi
            Kopt = value.(K)
        else
            lb=xi
        end
    end
    xiopt = ub=#

    solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
    model = Model(solver)
    @variable(model, K[1:m, 1:n])
    @variable(model, xitemp^2>=xi>=lowbound^2)
    #@variable(model, v>=0)
    @objective(model, Min, xi+10*(xi-xitemp^2)^2)
    for i in 1:numTraj
        @SDconstraint(model,[state0[:,i]'*(Ptemp-Q)*state0[:,i] (state[:,i]+B*K*state0[:,i])'*Ptemp state0[:,i]'*K';Ptemp*(state[:,i]+B*K*state0[:,i]) xi*Ptemp zeros(n,m);K*state0[:,i] zeros(m,n) inv(R)]>=0)
    end
    @SDconstraint(model, [Ptemp-Q K';K inv(R)]>=0)
    #@SDconstraint(model, [v*Matrix(I,m,m) K';K inv(R)]>=0)
    JuMP.optimize!(model)

    Kopt = value.(K)
    xiopt = value.(xi)

    lambdamin = 0
    for la in 0:0.01:1
        Kla = la*Ktemp+(1-la)*Kopt
        eigZ = eigvals(Ptemp-Q-Kla'*R*Kla)
        if maximum(eigZ)/minimum(eigZ) <= kappaU
            lambdamin = la
            break
        end
    end
    Knew = lambdamin*Ktemp+(1-lambdamin)*Kopt
    xinew = sqrt((1-lambdamin)*xiopt+lambdamin*xitemp^2)

    return Knew, xinew

end

function LQR_alternating_P(state0,state;B,Q,R,kappaU,Ktemp,xitemp)

    numTraj = size(state0)[2]
    (n,m) = size(B)

    stateK = zeros(dim,numTraj)

    for i in 1:numTraj
        stateK[:,i] = state[:,i]+B*Ktemp*state0[:,i]
    end

    solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
    model = Model(solver)
    @variable(model, P[1:n, 1:n] in PSDCone())
    @variable(model, v>=1e-6)
    #(model, s>=0)
    @objective(model, Min, tr(P))
    #@SDconstraint(model, P <= s*Matrix(I,n,n))
    @SDconstraint(model, P >= Q+Ktemp'*R*Ktemp)
    @SDconstraint(model, P >= Q+Ktemp'*R*Ktemp+v*Matrix(I,n,n))
    @SDconstraint(model, P - Q -Ktemp'*R*Ktemp <=kappaU*v*Matrix(I,n,n))
    for i in 1:numTraj
        @constraint(model, stateK[:,i]'*P*stateK[:,i] <= xitemp^2*state0[:,i]'*(P-Q-Ktemp'*R*Ktemp)*state0[:,i])
    end
    JuMP.optimize!(model)

    return value.(P)

end

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl