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

using Markdown
using InteractiveUtils

# ╔═╡ 2f87cbf3-75a1-45f4-afef-fbd324ec351f
using PlutoUI

# ╔═╡ 790d1872-28a8-461b-889b-93165adf79f0
begin
	using Flux
	import Flux.Data: DataLoader
	import Flux.Optimise: Optimiser, WeightDecay
	import Flux: onehotbatch, onecold
	import Flux.Losses: logitcrossentropy
end

# ╔═╡ b7d682c3-60c9-46c0-9113-81b826dde629
using Statistics, Random

# ╔═╡ 9f36cfca-c197-450f-a548-4e7f6f13fd4b
using CUDA

# ╔═╡ a72d6c9e-f30d-4029-ac54-ec40c3a87a31
using Colors

# ╔═╡ f3dae82b-699c-4237-aebc-39e02e451837
md"""# Ciencia de Datos pt. 3

	- Oscar Alejando Esquivel Flores
	- Óscar Anuar Alvarado Morán
	- Mario Horacio Garrido Czacki
"""

# ╔═╡ a2a725cb-5aa5-4627-aac4-018a4a63e4b4
PlutoUI.TableOfContents(title="🔬 Ciencia de datos pt. 3 🔭", depth = 5)

# ╔═╡ a7c163f2-df19-4cf2-aee2-349c75cd28b0
md"# Aprendizaje profundo"

# ╔═╡ 449ee180-4a38-45e9-ab65-255dbe8abee8
md"## Flux.jl
- [GitHub](https://github.com/FluxML/Flux.jl)
- [fluxml](https://fluxml.ai/)"

# ╔═╡ 2173f1c7-c168-455d-b57b-0ca7c2aa2355
md"El siguiente ejemplo se obtuvo de [modelzoo](https://github.com/FluxML/model-zoo/)"

# ╔═╡ 8ca80cbc-237a-408a-a070-a716e86a272f
md"Importamos las paqueterías necesarias"

# ╔═╡ 7a747250-756c-4b8e-aaf2-63140b167049
import Logging: with_logger

# ╔═╡ 34ad3424-0c46-4ba6-a27e-f9df0646da2d
import TensorBoardLogger: TBLogger, tb_overwrite, set_step!, set_step_increment!

# ╔═╡ 34084c2e-70fa-4fb8-abf0-768af7c86297
import ProgressMeter: @showprogress

# ╔═╡ 7199b894-88bf-4b59-8dbd-4fe16091dd8b
import MLDatasets# Parece que se tiene que tener instalado previamente scipy, yo lo hice con pip3 install scipy por la instalación que tengo.

# ╔═╡ c07a4f09-42a8-4850-9874-1f14ac6a64a5
import BSON

# ╔═╡ 8978ca53-ee1f-43b2-b614-0211af245ab7
md"Definimos la arquitectura que estaremos utilizando, en este caso es una conocida como LeNet"

# ╔═╡ cb7ad0e6-f5ec-11eb-0383-2dd50c9043f2
function LeNet5(; imgsize=(28,28,1), nclasses=10)
    out_conv_size = (imgsize[1]÷4 - 3, imgsize[2]÷4 - 3, 16)

    return Chain(
            Conv((5, 5), imgsize[end]=>6, relu),
            MaxPool((2, 2)),
            Conv((5, 5), 6=>16, relu),
            MaxPool((2, 2)),
            flatten,
            Dense(prod(out_conv_size), 120, relu),
            Dense(120, 84, relu),
            Dense(84, nclasses)
          )
end

# ╔═╡ 47c47f69-2ac8-4c9f-8a04-854a66eb58f8
md"Ponemos esto porque Pluto tiene un problema al querer manipular los datos de MNIST"

# ╔═╡ d2c615e6-973b-45d4-9d5c-fb7a5157214a
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# ╔═╡ 7ad81da6-4bf4-49b4-a289-4165b9a62136
md"Mostramos un ejemplo del conjunto de datos en pantalla"

# ╔═╡ c69d92f2-db3e-4927-82cd-793da8bc8695
begin
	xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
	Gray.(xtrain[:,:,1]')
end

# ╔═╡ 3c363cce-f75b-4996-a173-9861d54f7889
md"Definimos una función para obtener los datos como dataloaders, que será la entrada de la red"

# ╔═╡ d8d474c2-704b-4ff7-85cd-70a87b49dde3
function get_data(args)
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32)

    xtrain = reshape(xtrain, 28, 28, 1, :)
    xtest = reshape(xtest, 28, 28, 1, :)

    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    train_loader = DataLoader((xtrain, ytrain), batchsize=args.batchsize, shuffle=true)
    test_loader = DataLoader((xtest, ytest),  batchsize=args.batchsize)

    return train_loader, test_loader
end

# ╔═╡ b17f27d6-8b39-4317-8fe0-dfe5e1e44a9b
md"Definimos la función de pérdidas, ésta será la entropía cruzada."

# ╔═╡ 9a432f10-efeb-4c3c-812a-aa20cd393b5f
loss(ŷ, y) = logitcrossentropy(ŷ, y)

# ╔═╡ 72f5b8e2-79a8-4cdc-bc22-baa4b9dc59e7
md"Definimos esta función que evalúa las pérdidas y obtiene la exactitud"

# ╔═╡ af218b28-7596-4611-ad02-a10a6409d9b9
md"Definimos algunas funciones que nos serán útiles después."

# ╔═╡ db1571b3-4476-42bb-a107-5955af9d1b0d
num_params(model) = sum(length, Flux.params(model))

# ╔═╡ 2dbae0a8-93e8-4aa0-8792-f610e42e847f
round4(x) = round(x, digits=4)

# ╔═╡ 7d0ea233-9a06-40c2-9305-4fd6c01019e1
function eval_loss_accuracy(loader, model, device)
    l = 0f0
    acc = 0
    ntot = 0
    for (x, y) in loader
        x, y = x |> device, y |> device
        ŷ = model(x)
        l += loss(ŷ, y) * size(x)[end]
        acc += sum(onecold(ŷ |> cpu) .== onecold(y |> cpu))
        ntot += size(x)[end]
    end
    return (loss = l/ntot |> round4, acc = acc/ntot*100 |> round4)
end

# ╔═╡ 7a977696-b174-4878-bf48-a4e3740f4e0a
md"Definimos el ciclo de entrenamiento junto con algunas configuraciones para mandar mensajes del avance de la red y sobre el dispositivo en el que nos encontremos"

# ╔═╡ bad60fef-5fc5-4f38-a904-389037a707d5
Base.@kwdef mutable struct Args
    η = 3e-4             # learning rate
    λ = 0                # L2 regularizer param, implemented as weight decay
    batchsize = 128      # batch size
    epochs = 10          # number of epochs
    seed = 0             # set seed > 0 for reproducibility
    use_cuda = true      # if true use cuda (if available)
    infotime = 1 	     # report every `infotime` epochs
    checktime = 5        # Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    tblogger = true      # log training with tensorboard
    savepath = "./modelo/"       # results path
end

# ╔═╡ d2edf222-0bf0-4a86-94c7-e65f01ca955d
function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.use_cuda && CUDA.functional()

    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ## DATA
    train_loader, test_loader = get_data(args)
    @info "Dataset MNIST: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    model = LeNet5() |> device
    @info "LeNet5 model: $(num_params(model)) trainable params"

    ps = Flux.params(model)

    opt = ADAM(args.η)
    if args.λ > 0 # add weight decay, equivalent to L2 regularization
        opt = Optimiser(WeightDecay(args.λ), opt)
    end

    ## LOGGING UTILITIES
    if args.tblogger
        tblogger = TBLogger(args.savepath, tb_overwrite)
        set_step_increment!(tblogger, 0) # 0 auto increment since we manually set_step!
        @info "TensorBoard logging at \"$(args.savepath)\""
    end

    function report(epoch)
        train = eval_loss_accuracy(train_loader, model, device)
        test = eval_loss_accuracy(test_loader, model, device)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
        if args.tblogger
            set_step!(tblogger, epoch)
            with_logger(tblogger) do
                @info "train" loss=train.loss  acc=train.acc
                @info "test"  loss=test.loss   acc=test.acc
            end
        end
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch in 1:args.epochs
        @showprogress for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                    ŷ = model(x)
                    loss(ŷ, y)
                end

            Flux.Optimise.update!(opt, ps, gs)
        end

        ## Printing and logging
        epoch % args.infotime == 0 && report(epoch)
        if args.checktime > 0 && epoch % args.checktime == 0
            !ispath(args.savepath) && mkpath(args.savepath)
            modelpath = joinpath(args.savepath, "model.bson")
            let model = cpu(model) #return model to cpu before serialization
                BSON.@save modelpath model epoch
            end
            @info "Model saved in \"$(modelpath)\""
        end
    end
	model
end

# ╔═╡ d846b086-bf2b-411a-8d50-728fbfd99645
model = train()

# ╔═╡ 15b2ed6a-2dd7-42fc-b66f-df8e185049cf
md"Veamos cómo se comporta en el conjunto de pruebas"

# ╔═╡ 8a256b50-b395-4818-b7ca-0896d7179828
begin
	xtest, ytest = MLDatasets.MNIST.testdata(Float32)
	xtest = reshape(xtest, 28, 28, 1, :)
	test_loader = DataLoader((xtest, ytest),  batchsize=128)
	0
end

# ╔═╡ c425e0c1-50b8-4320-bfe5-b001496a97a2
begin
	pred = 0
	model_cpu = model |> cpu
	for (x,y) ∈ test_loader
		pred = model_cpu(x)
		break
	end
end

# ╔═╡ 05290564-2ca7-4a83-8367-68a069a0e612
var = 127

# ╔═╡ af215c7e-fd57-4ba7-ab98-82e86cbab88c
begin
	f1(x) = getindex.(argmax(x, dims=1), 1)
	preds = f1(pred) .- 1;
	preds[var]
end

# ╔═╡ 71701f5f-fabc-4db0-bde2-9a6b3e65ed6a
ytest[var]

# ╔═╡ e478c1d4-a838-49fc-a1f3-9bd02189ccae
Gray.(xtest[:,:,1,var]')

# ╔═╡ 322b2033-adf1-489b-bf3a-957c942fe783
md"# Big Data"

# ╔═╡ a77eefa6-c3b1-47a9-b45a-2045001407f5
md"## Spark.jl
- [GitHub](https://github.com/dfdx/Spark.jl)"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
MLDatasets = "eb30cadb-4394-5ae3-aed4-317e484a6458"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
TensorBoardLogger = "899adc3e-224a-11e9-021f-63837185c80f"

[compat]
BSON = "~0.3.3"
CUDA = "~3.3.4"
Colors = "~0.12.8"
Flux = "~0.12.6"
MLDatasets = "~0.5.9"
PlutoUI = "~0.7.9"
ProgressMeter = "~1.7.1"
TensorBoardLogger = "~0.1.18"
"""

# ╔═╡ 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"

[[AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

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

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

[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "2e004e61f76874d153979effc832ae53b56c20ee"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.22"

[[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"

[[BSON]]
git-tree-sha1 = "92b8a8479128367aaab2620b8e73dff632f5ae69"
uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
version = "0.3.3"

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

[[BinDeps]]
deps = ["Libdl", "Pkg", "SHA", "URIParser", "Unicode"]
git-tree-sha1 = "1289b57e8cf019aede076edab0587eb9644175bd"
uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
version = "1.0.2"

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[Blosc]]
deps = ["Blosc_jll"]
git-tree-sha1 = "84cf7d0f8fd46ca6f1b3e0305b4b4a37afe50fd6"
uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
version = "0.7.0"

[[Blosc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "e747dac84f39c62aff6956651ec359686490134e"
uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
version = "1.21.0+0"

[[BufferedStreams]]
deps = ["Compat", "Test"]
git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f"
uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d"
version = "1.0.0"

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

[[CRC32c]]
uuid = "8bf52ea8-c179-5cab-976a-9e18b702a9bc"

[[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"

[[ChainRules]]
deps = ["ChainRulesCore", "Compat", "LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "11567f2471013449c2fcf119f674c681484a130e"
uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2"
version = "1.5.1"

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

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.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"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[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"

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

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

[[DataDeps]]
deps = ["BinaryProvider", "HTTP", "Libdl", "Reexport", "SHA", "p7zip_jll"]
git-tree-sha1 = "4f0e41ff461d42cfc62ff0de4f1cd44c6e6b3771"
uuid = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
version = "0.7.7"

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

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

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

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[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"

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

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.10.1"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8c8eac2af06ce35973c3eadb4ab3243076a408e7"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.1"

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

[[Flux]]
deps = ["AbstractTrees", "Adapt", "ArrayInterface", "CUDA", "CodecZlib", "Colors", "DelimitedFiles", "Functors", "Juno", "LinearAlgebra", "MacroTools", "NNlib", "NNlibCUDA", "Pkg", "Printf", "Random", "Reexport", "SHA", "Statistics", "StatsBase", "Test", "ZipFile", "Zygote"]
git-tree-sha1 = "1286e5dd0b4c306108747356a7a5d39a11dc4080"
uuid = "587475ba-b771-5e3f-ad9e-33799f191a9c"
version = "0.12.6"

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "b5e930ac60b613ef3406da6d4f42c35d8dc51419"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.19"

[[Functors]]
deps = ["MacroTools"]
git-tree-sha1 = "4cd9e70bf8fce05114598b663ad79dfe9ae432b3"
uuid = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
version = "0.2.3"

[[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"

[[GZip]]
deps = ["Libdl"]
git-tree-sha1 = "039be665faf0b8ae36e089cd694233f5dee3f7d6"
uuid = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
version = "0.5.1"

[[Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "2c1cf4df419938ece72de17f368a021ee162762e"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.0"

[[HDF5]]
deps = ["Blosc", "Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires"]
git-tree-sha1 = "83173193dc242ce4b037f0263a7cc45afb5a0b85"
uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
version = "0.15.6"

[[HDF5_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "fd83fa0bde42e01952757f01149dd968c06c4dba"
uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
version = "1.12.0+1"

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

[[IRTools]]
deps = ["InteractiveUtils", "MacroTools", "Test"]
git-tree-sha1 = "95215cd0076a150ef46ff7928892bc341864c73c"
uuid = "7869d1d1-7146-5819-86e3-90919afe41df"
version = "0.4.3"

[[IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[ImageCore]]
deps = ["AbstractFFTs", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "db645f20b59f060d8cfae696bc9538d13fd86416"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.8.22"

[[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"

[[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"

[[Juno]]
deps = ["Base64", "Logging", "Media", "Profile"]
git-tree-sha1 = "07cb43290a840908a771552911a6274bc6c072c7"
uuid = "e5e0dc1b-0480-54bc-9374-aad01c23163d"
version = "0.8.4"

[[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"

[[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"

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

[[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"

[[Lz4_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
version = "1.9.3+0"

[[MAT]]
deps = ["BufferedStreams", "CodecZlib", "HDF5", "SparseArrays"]
git-tree-sha1 = "5c62992f3d46b8dce69bdd234279bb5a369db7d5"
uuid = "23992714-dd62-5051-b70f-ba57cb901cac"
version = "0.10.1"

[[MLDatasets]]
deps = ["BinDeps", "ColorTypes", "DataDeps", "DelimitedFiles", "FixedPointNumbers", "GZip", "MAT", "PyCall", "Requires"]
git-tree-sha1 = "65cb0a663d65d0b782ba74bfc3982ba51eb85485"
uuid = "eb30cadb-4394-5ae3-aed4-317e484a6458"
version = "0.5.9"

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

[[MappedArrays]]
git-tree-sha1 = "18d3584eebc861e311a552cbb67723af8edff5de"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.0"

[[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"

[[Media]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "75a54abd10709c01f1b86b84ec225d26e840ed58"
uuid = "e89f7d12-3494-54d1-8411-f7d8b9ae1f27"
version = "0.5.0"

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

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

[[MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

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

[[NNlib]]
deps = ["Adapt", "ChainRulesCore", "Compat", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "16520143f067928bb69eee59ac8bca06be1e43b8"
uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
version = "0.7.27"

[[NNlibCUDA]]
deps = ["CUDA", "LinearAlgebra", "NNlib", "Random", "Statistics"]
git-tree-sha1 = "a7de026dc0ff9f47551a16ad9a710da66881b953"
uuid = "a00861dc-f156-4864-bf3c-e6376f28a68d"
version = "0.1.7"

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

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

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "5cc97a6f806ba1b36bac7078b866d4297ae8c463"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.4"

[[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"

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

[[PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "0fa5e78929aebc3f6b56e1a88cf505bb00a354c4"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.8"

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

[[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"

[[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"

[[Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[ProtoBuf]]
deps = ["Logging", "protoc_jll"]
git-tree-sha1 = "3142e13697a7f6363798d1241daf3e6cead83008"
uuid = "3349acd9-ac6a-5e09-bcdb-63829b23a429"
version = "0.11.2"

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

[[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 = "043da614cc7e95c703498a491e2c21f58a2b8111"
uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
version = "1.5.3"

[[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"

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

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

[[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"

[[StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[Static]]
deps = ["IfElse"]
git-tree-sha1 = "62701892d172a2fa41a1f829f66d2b0db94a9a63"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.3.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"

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

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

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

[[TensorBoardLogger]]
deps = ["CRC32c", "FileIO", "ImageCore", "ProtoBuf", "Requires", "StatsBase"]
git-tree-sha1 = "96d160e85038f6d89e6c9b91492466f9a7d454f2"
uuid = "899adc3e-224a-11e9-021f-63837185c80f"
version = "0.1.18"

[[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"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.5"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[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"

[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

[[ZipFile]]
deps = ["Libdl", "Printf", "Zlib_jll"]
git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7"
uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
version = "0.9.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"

[[Zygote]]
deps = ["AbstractFFTs", "ChainRules", "ChainRulesCore", "DiffRules", "Distributed", "FillArrays", "ForwardDiff", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NaNMath", "Random", "Requires", "SpecialFunctions", "Statistics", "ZygoteRules"]
git-tree-sha1 = "f01bac579bb397ab138aed7e9e3f80ef76d055f7"
uuid = "e88e6eb3-aa80-5325-afca-941959d7151f"
version = "0.6.19"

[[ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.1"

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

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

[[protoc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "89b92b537ffde09cab61ad20636da135d0791007"
uuid = "c7845625-083e-5bbe-8504-b32d602b7110"
version = "3.15.6+0"
"""

# ╔═╡ Cell order:
# ╟─f3dae82b-699c-4237-aebc-39e02e451837
# ╠═2f87cbf3-75a1-45f4-afef-fbd324ec351f
# ╟─a2a725cb-5aa5-4627-aac4-018a4a63e4b4
# ╟─a7c163f2-df19-4cf2-aee2-349c75cd28b0
# ╟─449ee180-4a38-45e9-ab65-255dbe8abee8
# ╟─2173f1c7-c168-455d-b57b-0ca7c2aa2355
# ╟─8ca80cbc-237a-408a-a070-a716e86a272f
# ╠═790d1872-28a8-461b-889b-93165adf79f0
# ╠═b7d682c3-60c9-46c0-9113-81b826dde629
# ╠═7a747250-756c-4b8e-aaf2-63140b167049
# ╠═34ad3424-0c46-4ba6-a27e-f9df0646da2d
# ╠═34084c2e-70fa-4fb8-abf0-768af7c86297
# ╠═7199b894-88bf-4b59-8dbd-4fe16091dd8b
# ╠═c07a4f09-42a8-4850-9874-1f14ac6a64a5
# ╠═9f36cfca-c197-450f-a548-4e7f6f13fd4b
# ╠═a72d6c9e-f30d-4029-ac54-ec40c3a87a31
# ╟─8978ca53-ee1f-43b2-b614-0211af245ab7
# ╠═cb7ad0e6-f5ec-11eb-0383-2dd50c9043f2
# ╠═47c47f69-2ac8-4c9f-8a04-854a66eb58f8
# ╠═d2c615e6-973b-45d4-9d5c-fb7a5157214a
# ╟─7ad81da6-4bf4-49b4-a289-4165b9a62136
# ╠═c69d92f2-db3e-4927-82cd-793da8bc8695
# ╟─3c363cce-f75b-4996-a173-9861d54f7889
# ╠═d8d474c2-704b-4ff7-85cd-70a87b49dde3
# ╟─b17f27d6-8b39-4317-8fe0-dfe5e1e44a9b
# ╠═9a432f10-efeb-4c3c-812a-aa20cd393b5f
# ╟─72f5b8e2-79a8-4cdc-bc22-baa4b9dc59e7
# ╠═7d0ea233-9a06-40c2-9305-4fd6c01019e1
# ╟─af218b28-7596-4611-ad02-a10a6409d9b9
# ╠═db1571b3-4476-42bb-a107-5955af9d1b0d
# ╠═2dbae0a8-93e8-4aa0-8792-f610e42e847f
# ╟─7a977696-b174-4878-bf48-a4e3740f4e0a
# ╠═bad60fef-5fc5-4f38-a904-389037a707d5
# ╠═d2edf222-0bf0-4a86-94c7-e65f01ca955d
# ╠═d846b086-bf2b-411a-8d50-728fbfd99645
# ╟─15b2ed6a-2dd7-42fc-b66f-df8e185049cf
# ╠═8a256b50-b395-4818-b7ca-0896d7179828
# ╠═c425e0c1-50b8-4320-bfe5-b001496a97a2
# ╠═05290564-2ca7-4a83-8367-68a069a0e612
# ╠═af215c7e-fd57-4ba7-ab98-82e86cbab88c
# ╠═71701f5f-fabc-4db0-bde2-9a6b3e65ed6a
# ╠═e478c1d4-a838-49fc-a1f3-9bd02189ccae
# ╟─322b2033-adf1-489b-bf3a-957c942fe783
# ╟─a77eefa6-c3b1-47a9-b45a-2045001407f5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: Curso Julia/Github - Material/src/CienciaDeDatos/AprendizajeProfundo.jl


text:
using .Origami

const evaluate = Origami.evaluate
const P = Point{Float64}

@testset "Origami" begin
    @test isempty(Origami.normal!([P(0,0), P(1,1)]))
    @test evaluate(P[],[], P(0,0)) == 0
    shape = [P(0,0), P(0,1000), P(1000,1000), P(1000,0), P(0,0)]
    @test evaluate(shape, [], P(1,1)) == 1
    folds = [(P(-5,-5), P(10,10)), (P(10,750), P(0,750))]
    probes = [P(100,600), P(800,600), P(300,400), P(100,100), P(500,500), P(200,500)]
    result = [4, 2, 2, 0, 0, 2]
    for (p, r) in zip(probes, result)
        @test evaluate(shape, folds, p) == r
    end
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using StatsPlots, Distributions, CSV, DataFrames, Measures
	include(datadir("mvParams.jl"))
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.25"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
realData = CSV.read(datadir("temperatures.csv"));

# ╔═╡ 14d9ae82-2510-11eb-1223-13889002f669
begin
	N = 10^5
	biNorm = MvNormal(meanVect,covMat)
	syntheticData = DataFrame(rand(MvNormal(meanVect,covMat),N)')
	rename!(syntheticData,[:x1=>:Brisbane, :x2 => :GoldCoast])
end

# ╔═╡ 14d9e104-2510-11eb-3fe6-d3ca0f565047
begin
	default(c=cgrad([:blue, :red]),
		xlabel="Brisbane Temperature",
		ylabel="Gold Coast Temperature")

	p1 = marginalhist(realData.Brisbane, realData.GoldCoast, bins=10:45)
	p2 = marginalhist(syntheticData.Brisbane, syntheticData.GoldCoast, bins=10:.5:45)

	plot(p1,p2, size = (1000,500), margin = 10mm)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.25"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═14d9ae82-2510-11eb-1223-13889002f669
# ╠═14d9e104-2510-11eb-3fe6-d3ca0f565047
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
@testset "parents and children" begin
    cell = Cell(SVector(-1., -2, -3), SVector(2., 4, 6), 0)
    @test cell.boundary.origin == SVector(-1., -2, -3)
    @test cell.boundary.widths == SVector(2., 4, 6)
    @test cell.children === nothing
    @test isleaf(cell)

    split!(cell, [i for i in 1:8])
    @test !isleaf(cell)
    @test cell[1,1,1].boundary.origin == cell.boundary.origin
    @test cell[2,2,2].boundary.origin + cell[2,2,2].boundary.widths == cell.boundary.origin + cell.boundary.widths
    @test cell[1,1,1].data == 1
    @test cell[2,1,1].data == 2
    for child in children(cell)
        @test parent(child) === cell
    end
end

@testset "type promotions" begin
    c1 = RegionTrees.Cell(SVector(0.0, 0.0), SVector(1.0, 1.0))
    @test typeof(c1.boundary.origin) == SVector{2, Float64}
    c2 = RegionTrees.Cell(SVector(0, 0), SVector(1, 1))
    @test typeof(c2.boundary.origin) == SVector{2, Float64}

    @inferred RegionTrees.Cell(SVector(0, 0), SVector(1, 1))
    @inferred RegionTrees.Cell(SVector(0.0, 0.0), SVector(1.0, 1.0))
end

@testset "find leaf" begin
    cell = Cell(SVector(-1., -2, -3), SVector(2., 4, 6), 0)
    @test findleaf(cell, SVector(-0.01, -0.01, -0.01)) === cell

    split!(cell, [0 for i in 1:8])
    @test findleaf(cell, SVector(-0.01, -0.01, -0.01)) === cell[1,1,1]
    @test findleaf(cell, SVector(0.01, -0.01, -0.01)) === cell[2,1,1]
end


text:
__precompile__()

module Plexifications

import Base: +, -, *, /, \, real, isreal, inv, abs2, norm, show, zero, one, iszero, ==, exp
# import Unreal: unreal

"Error when finding the inverse of a zero divisor."
const ZeroDivisorInverse = "inverse of zero divisor"

"Error when the denominator in a quotient is zero."
const ZeroDenominator = "denominator is zero"

include("Plexification.jl")
include("Complexification.jl")
include("Perplexification.jl")
include("Nilplexification.jl")
include("Cmplex.jl")
include("Prplex.jl")
include("Nlplex.jl")
include("BiComplex.jl")
include("BiPerplex.jl")
include("BiNilplex.jl")
include("DualComplex.jl")
include("DualPerplex.jl")
include("TriComplex.jl")
include("TriPerplex.jl")
include("TriNilplex.jl")
include("HyperComplex.jl")
include("HyperPerplex.jl")
include("random.jl")

# Types
export Plexification, Complexification, Perplexification, Nilplexification,
       Cmplex, Prplex, Nlplex,
       BiComplex, BiPerplex, BiNilplex, DualComplex, DualPerplex,
       TriComplex, TriPerplex, TriNilplex, HyperComplex, HyperPerplex

# Methods
export unreal, array, iszerodivisor, random, quadrance,
       barH, barR, barA,
       barG, barQ, barB,
       barF, barP, barC,
       commutator, associator, jacobiator, alternatorL, alternatorR, flexator,
       lozenge, anti_lozenge

end


text:
# SteinBanana
# The banana posterior distribution p(theta | y_1,...,y_L) induced by the model
#
#   y_l | theta ~iid N(theta1 + theta2^2, sigma2y)
#   theta_1, theta_2 ~iid N(0, sigma2theta)
#
# This distribution is analyzed in the papers
#
#   http://arxiv.org/pdf/1211.3759v2.pdf
#   http://jmlr.org/proceedings/papers/v33/dubois14.pdf

using Distributions: Normal, logpdf
using SteinDistributions

type SteinBanana <: SteinPosterior
    y::Array{Float64, 1}
    sigma2y::Float64
    sigma2theta::Float64
    c1::Array{Float64, 1}
    c2::Array{Float64, 1}
    c3::Array{Float64, 1}
end

# Constructor with uniform setting of (c1,c2,c3) and
# default setting of sigma2y and sigma2theta
defaultsigma2y = 4.0
SteinBanana(y; sigma2y = defaultsigma2y, sigma2theta = 1.0) =
SteinBanana(y, sigma2y, sigma2theta,
            [1.0 for i=1:2],[1.0 for i=1:2],[1.0 for i=1:2])

# Returns lower bound of support
function supportlowerbound(d::SteinBanana, j::Int64)
    -Inf
end

# Returns upper bound of support
function supportupperbound(d::SteinBanana, j::Int64)
    Inf
end

# Draws n datapoints from the likelihood distribution
#   N(theta1 + theta2^2, sigma2y)
#
# Args:
#   mu = theta1 + theta2^2
function randbanana(mu::Float64, n::Int64)
    datanormal = Normal(mu, sqrt(defaultsigma2y))
    rand(datanormal, n)
end

# Computes the log prior for N(theta1 + theta2^2, sigma2y)
function logprior(d::SteinBanana,
                  theta::Array{Float64, 1})
    priornormal = Normal(0.0, sqrt(d.sigma2theta))
    sum(logpdf(priornormal, theta))
end

# Returns log p(theta | d.y[idx]), the log density of the posterior distribution
# evaluated at theta, conditioned only on those datapoints indexed by idx
function loglikelihood(d::SteinBanana,
                       theta::Array{Float64, 1};
                       idx=1:length(d.y))
    y = d.y[idx]
    theta0 = theta[1] + theta[2]^2

    datanormal = Normal(theta0, sqrt(d.sigma2y))
    sum(logpdf(datanormal, y))
end

# Returns grad_{theta} log p(d.y[idx] | theta)
function gradloglikelihood(d::SteinBanana,
                           theta::Array{Float64, 1};
                           idx=1:length(d.y))
    y = d.y[idx]
    theta0 = theta[1] + theta[2]^2

    gradtheta1 = sum(y - theta0) / d.sigma2y
    gradtheta2 = 2 * theta[2] * gradtheta1

    [gradtheta1, gradtheta2]
end

# Returns grad_{theta} log p(theta)
function gradlogprior(d::SteinBanana,
                      theta::Array{Float64, 1})
    -(1/d.sigma2theta) .* theta
end

# diaghessian(x) = (grad_1 b1(x), grad_2 b2(x), ....) where b = grad log p.
function diaghessian(d::SteinBanana, theta::Array{Float64,1}; idx=1:length(d.y))
    m = length(idx)
    priorterm = -1.0/d.sigma2theta
    ysum = sum(d.y[idx])
    likelihoodterm = [-m, -2*ysum + 2*m*theta[1] + 6*m*theta[2]^2] ./ d.sigma2y
    priorterm .+ likelihoodterm
end

function numdatapoints(d::SteinBanana)
    length(d.y)
end


text:
using MyProject
using Test

@testset "foo checks" begin
    @test foo(0) < 1E-16
end


text:
mutable struct SingleDUParams
    d::Float64
    t::Float64
    z::Float64
    dz::Float64
    ϕ::Float64
    t₀::Float64
end

struct ROyFit
    sdp::SingleDUParams
    sdp_initial::SingleDUParams
    Q::Float64
    selected_hits::Vector{CalibratedHit}
    model::Model
end

"""
    function svdfit(hits::Vector{CalibratedHit})

Uses SVD to do a fast and dirty track prefit. Provide hits with a multiplicity
of at least 2.
"""
function svdfit(hits::Vector{CalibratedHit})
    t₀ = hits[div(length(hits), 2)].t
    pos, dir = svdfit(matrix([h.pos for h in hits]))
    if (last(hits).pos - first(hits).pos) ⋅ dir  < 0.0
        dir *= -1
    end
    return Track(dir, pos, t₀)
end

"""
    function prefit(hits::Vector{CalibratedHit})

Performs the prefit algorithm which was used in DUMAND II.
"""
function prefit(hits::Vector{CalibratedHit})
    N = length(hits)
    D = 0.0
    dir = [0.0, 0.0, 0.0]
    pos = [0.0, 0.0, 0.0]
    pes = [max(1, (h.tot - 26.3) / 4.5) for h in hits]
    # pes = [h.tot for h in hits]
    # pes = [h.multiplicity.count for h in hits]
    @inbounds for i ∈ 1:N
        @inbounds for k ∈ 1:N
            if i == k
                continue
            end
            t_i = hits[i].t
            t_k = hits[k].t
            q_ik = pes[i] * pes[k]
            t_ik = t_i * t_k
            pos += q_ik * (hits[i].pos*(t_k^2 - t_ik) + hits[k].pos*(t_i^2 - t_ik))
            dir += q_ik * (hits[i].pos - hits[k].pos) * (t_i - t_k)
            D += q_ik * (t_i - t_k)^2
        end
    end
    # dir = normalize(dir/D)
    return NeRCA.Track(dir/D, pos/D, 0)
end

"""
    function make_cherenkov_calculator(track::Track; v=2.99792458e8)

Returns a function which calculates the arrival time of a Cherenkov photon
at a given position.
"""
function make_cherenkov_calculator(track::Track; v=2.99792458e8, n=N_SEAWATER)
    c_medium = c_0/n
    β = v/c_0
    θ = acos(min(1/(β*n), 1))
    θ′ = π - θ
    track_dir = normalize(track.dir)
    t₀ = track.time

    pos::Position -> begin

        # minimal distance between hit and track
        track_distance = pld3(track.pos, pos, track_dir)
        # the path the Cherenkov photon has to travel
        cherenkov_path_length = track_distance / sin(θ)
        distance_vector = pos - track.pos
        # distance between particle position and hit
        distance = norm(distance_vector)

        t_cherenkov = cherenkov_path_length / c_medium*1e9

        η = angle_between(track_dir, distance_vector)
        χ = θ - η
        particle_travel_path_length = distance / sin(θ′) * sin(χ)
        t = t₀ + particle_travel_path_length / v*1e9 + t_cherenkov
        return t
    end
end

"""
    function make_cherenkov_calculator(track::Track, event_info::Union{MCEventInfo,DAQEventInfo}; v=2.99792458e8)

Returns a function which calculates the arrival time of a Cherenkov photon
at a given position.
"""
function make_cherenkov_calculator(track::Track, event_info::Union{MCEventInfo,DAQEventInfo}; v=2.99792458e8)
    jte_time = make_mc_time_converter(event_info)(track.time)
    make_cherenkov_calculator(Track(track.dir, track.pos, jte_time), v=v)
end

"""
    function make_cherenkov_calculator(d_closest, t_closest, z_closest, dir_z, t₀)

Returns a function which calculates the arrival time of a Cherenkov photon
at a given position.
"""
function make_cherenkov_calculator(d_closest, t_closest, z_closest, dir_z, t₀; n=N_SEAWATER)
    c_ns = c_0 / 1e9
    d_γ(z) = n/√(n^2 - 1) * √(d_closest^2 + (z-z_closest)^2 * (1 - dir_z^2))
    t(z) = (t₀) + 1/c_ns * ((z - z_closest)*dir_z + (n^2 - 1)/n * d_γ(z))
    d_γ, t
end

"""
    function make_cherenkov_calculator(sdp::SingleDUParams)

Returns a function which calculates the arrival time of a Cherenkov photon
at a given position.
"""
function make_cherenkov_calculator(sdp::SingleDUParams; n=N_SEAWATER)
    c_ns = c_0 / 1e9
    d_γ(z) = n/√(n^2 - 1) * √(sdp.d^2 + (z-sdp.z)^2 * (1 - sdp.dz^2))
    t(z) = (sdp.t₀) + 1/c_ns * ((z - sdp.z)*sdp.dz + (n^2 - 1)/n * d_γ(z))
    d_γ, t
end

"""
    function single_du_params(track::NeRCA.Track)

Calculates five parameters to describe a track for a single DU case.
"""
function single_du_params(track::NeRCA.Track)
    c_ns = c_0 / 1e9
    pos = track.pos
    dir = Direction(normalize(track.dir))
    t₀ = track.time
    proj = pos ⋅ dir
    z_closest = (pos.z - dir.z*proj) / (1 - dir.z^2)
    t_closest = t₀ + (z_closest * dir.z - proj)/c_ns
    p_t_closest = pos + c_ns * (t_closest - t₀) * dir
    d_closest = √(p_t_closest[1]^2 + p_t_closest[2]^2)

    d_closest, t_closest, z_closest, dir.z, t₀
end

"""
    function single_du_params(t::NeRCA.MCTrack)

Calculates five parameters to describe a MC track for a single DU case.
"""
function single_du_params(t::NeRCA.MCTrack, time)
    single_du_params(NeRCA.Track([t.dir_x, t.dir_y, t.dir_z], [t.pos_x, t.pos_y, t.pos_z], time))
end

# TODO: delete this function
function cherenkov_plausible(Δt, Δz; time_extra=10, n=N_SEAWATER)
    Δt < Δz * n / NeRCA.c_0*1e9 + time_extra
end

function create_hit_pool(hits::T) where T<:Vector{NeRCA.CalibratedHit}
    hit_pool = Dict{Int, T}()
    for hit in hits
        if !haskey(hit_pool, hit.floor)
            hit_pool[hit.floor] = T()
            sizehint!(hit_pool[hit.floor], length(hits))
        end
        push!(hit_pool[hit.floor], hit)
    end
    hit_pool
end

"""
    function select_hits(du_hits, hit_pool; Δt₋=10, Δz=9, new_hits=nothing)

Returns the seed hits suited for a Cherenkov hit time residual based
reconstruction algorithm.

The `du_hits` should only contain hits for a single DU. The hit_pool holds
all other hit candidates (e.g. created by `create_hit_pool()`). `Δt₋` is
the allowed negative time error for the arrival time, `Δz` distance
between two floors.
"""
function select_hits(du_hits, hit_pool; Δt₋=10, Δz=9, new_hits=nothing)
    intervals = Dict{Int16, Interval{Float64}}()
    hits₀ = NeRCA.create_hit_pool(du_hits)

    extended = Vector{CalibratedHit}()

    function time_interval(t₀, tₜ, Δfloor)
        t₋ = tₜ - Δfloor*Δt₋
        t₊ = max(t₀ + Δfloor*Δz*NeRCA.N_SEAWATER/NeRCA.c_0*1e9 + Δt₋, tₜ + Δfloor*Δt₋)
        @interval(t₋, t₊)
    end

    function add_if_in_interval(floor, interval)
        # println("    searching for hit candidate on floor $floor")
        if haskey(hit_pool, floor)
            floor_hits = hit_pool[floor]
            # println("       $(length(floor_hits)) hits to check")
            for (iₕ, hit) ∈ enumerate(floor_hits)
                # println("           -> hit time: $(hit.t)")
                if hit.t ∈ interval
                    # println("           --> in interval!")
                    if  haskey(hits₀, floor) && hit.t > hits₀[floor][1].t
                        # println("             !! got a hit already which is earlier")
                        # println("             -- removing consecutive hits")
                        deleteat!(floor_hits, iₕ:length(floor_hits))
                        return
                    end
                    # println("             ++ adding hit")
                    push!(extended, hit)
                    # println("             -- removing consecutive hits")
                    deleteat!(floor_hits, iₕ:length(floor_hits))
                    return
                end
            end
        end
    end

    hit_seeds = new_hits == nothing ? du_hits : new_hits

    for hit in hit_seeds
        t₀ = hit.t
        floor = hit.floor
        # println("Checking floor $floor")
        for (floor₋, floor₋₋, floor₊, floor₊₊) in [(floor-1, floor-2, floor+1, floor+2), (floor+1, floor+2, floor-1, floor-2)]
            if haskey(hits₀, floor₋)
                t₋ = first(hits₀[floor₋]).t
                if haskey(hits₀, floor₋₋)
                    # println(" => Got hit on $(floor₋₋) and $(floor₋)")
                    t₋₋ = first(hits₀[floor₋₋]).t
                    t₊ = t₋ + 2(t₋ - t₋₋)
                    t₊₊ = t₋ + 3(t₋ - t₋₋)
                    interval₊ = time_interval(t₀, t₊, 1)
                    interval₊₊ = time_interval(t₀, t₊₊, 2)
                    # println("    interval+: $interval₊")
                    # println("    interval++: $interval₊₊")
                    #interval₊ = @interval(minimum([t₀, t₋, t₋₋]) - Δt, maximum([t₀, t₋, t₋₋]) + Δt)
                    #interval₊₊ = @interval(minimum([t₀, t₋, t₋₋]) - 2Δt, maximum([t₀, t₋, t₋₋]) + 2Δt)
                    add_if_in_interval(floor₊, interval₊)
                    add_if_in_interval(floor₊₊, interval₊₊)
                else
                    t₊ = t₀ + (t₀ - t₋)
                    t₊₊ = t₀ + 2(t₀ - t₋)
                    interval₊ = time_interval(t₀, t₊, 1)
                    interval₊₊ = time_interval(t₀, t₊₊, 2)
                    add_if_in_interval(floor₊, interval₊)
                    add_if_in_interval(floor₊₊, interval₊₊)
                end
            elseif haskey(hits₀, floor₋₋)
                t₋₋ = first(hits₀[floor₋₋]).t
                t₊ = t₀ + (t₀ - t₋₋)/2
                t₊₊ = t₀ + 2(t₀ - t₋₋)/2
                interval₊ = time_interval(t₀, t₊, 1)
                interval₊₊ = time_interval(t₀, t₊₊, 2)
                add_if_in_interval(floor₊, interval₊)
                add_if_in_interval(floor₊₊, interval₊₊)
            end
        end

    end

    if length(extended) == 0
        return du_hits
    end
    extension = unique(h->h.floor, sort(vcat(du_hits, extended), by=h->h.t))
    return select_hits(extension, hit_pool; Δt₋=Δt₋, Δz=Δz, new_hits=extended)
end

struct SingleDUMinimiser <: Function
    z_positions::Vector{Float64}
    times::Vector{Float64}
    pmt_directions::Vector{Direction}
    tots::Vector{Int16}
    multiplicities::Vector{Int}
    max_multiplicity::Float64
    nphes::Vector{Float64}
    average_coinc_tots::Vector{Float64}
end

function SingleDUMinimiser(hits::Vector{CalibratedHit}, triggered_hits::Vector{CalibratedHit})
    n = length(hits)
    n_triggered = length(triggered_hits)
    z_positions = Vector{Float64}()
    times = Vector{Float64}()
    multiplicities = Vector{Int32}()
    pmt_directions = Vector{Direction}()
    tots = Vector{Int16}()
    sizehint!(z_positions, n)
    sizehint!(times, n)
    sizehint!(multiplicities, n)
    sizehint!(pmt_directions, n_triggered)
    @inbounds for i ∈ 1:n
        hit = hits[i]
        push!(z_positions, hit.pos.z)
        push!(tots, hit.tot)
        push!(times, hit.t)
        push!(multiplicities, hit.multiplicity.count)
        push!(pmt_directions, hit.dir)
    end
    max_multiplicity = maximum(multiplicities)
    nphe = Vector{Float32}()
    average_coinc_tots = Vector{Float32}()
    for hit in hits
        coinc_hits = filter(h->h.multiplicity.id == hit.multiplicity.id, hits)
        average_coinc_tot = mean([h.tot for h in coinc_hits])
        push!(average_coinc_tots, average_coinc_tot)
        estimated_nphes = sum([nphes(h.tot) for h in coinc_hits])
        push!(nphe, estimated_nphes)
    end
    SingleDUMinimiser(z_positions, times, pmt_directions, tots, multiplicities, max_multiplicity, nphe, average_coinc_tots)
end

"""
    function (s::SingleDUMinimiser)(d_closest, t_closest, z_closest, dir_z, ϕ, t₀)

The quality function to be minimised when performing the single DU fit.
`d_closest` is the closest distance between the track (starting at `t₀` and
the DU at time `t_closest` with the z-coordinate `z_closest`. The direction
is `dir_z` and `ϕ` (azimuth).

"""
function (s::SingleDUMinimiser)(d_closest, t_closest, z_closest, dir_z, ϕ, t₀)
    n = length(s.times)

    d_γ, ccalc = make_cherenkov_calculator(d_closest, t_closest, z_closest, dir_z, t₀)

    Q = 0.0
    # println("===============")
    @inbounds for i ∈ 1:n
        t = s.times[i] - slew(s.tots[i])
        z = s.z_positions[i]
        m = s.multiplicities[i]
        t_exp = ccalc(z)
        Δt = abs(t - t_exp)
        zenith_acceptance = 1 - NeRCA.zenith(Direction(0,0,1))/π
        photon_distance = d_γ(z)

        timing = Δt^2

        distance_term = (s.nphes[i] - 72.0 / photon_distance * zenith_acceptance)^2

        pmtᵩ = azimuth(-s.pmt_directions[i])
        ξ = (ϕ - pmtᵩ)^2

        Q += timing + distance_term + ξ
        # println((timing, distance_term, ξ))
    end
    # println("---------------")

    Q
end

struct MultiDUMinimiser <: Function
    hits::Vector{CalibratedHit}
end

function (m::MultiDUMinimiser)(x, y, z, θ, ϕ, t₀)
    pos = Position(x, y, z)
    # dir = Direction(cos(θ)*cos(ϕ), cos(θ)*sin(ϕ), sin(θ))
    dir = Direction(ϕ, θ)
    track = Track(dir, pos, t₀)

    # ccalc = make_cherenkov_calculator(track, v=v*1e9)
    ccalc = make_cherenkov_calculator(track)

    Q = 0.0
    for hit in m.hits
        t_exp = ccalc(hit.pos)
        Δt = abs(hit.t - t_exp)
        # println("expected: $(t_exp), got: $(hit.t), delta: $(Δt)")
        Q += (Δt - 2)^2
    end
    Q
end

function multi_du_fit(prefit, hits; print_level=0)
    ϕ_start = NeRCA.azimuth(prefit.dir)
    θ_start = NeRCA.zenith(prefit.dir)
    t₀_start = prefit.time
    pos = prefit.pos

    m = NeRCA.MultiDUMinimiser(hits)

    model = Model(
        with_optimizer(Ipopt.Optimizer,
                       print_level=print_level,
                       max_iter=200,
                       tol=1e-3)
       )
    register(model, :qfunc, 6, m, autodiff=true)

    Δpos = 50

    @variable(model, pos.x - Δpos <= x <= pos.x + Δpos, start=pos.x)
    @variable(model, pos.y - Δpos <= y <= pos.y + Δpos, start=pos.y)
    @variable(model, pos.z - Δpos <= z <= pos.z + Δpos, start=pos.z)
    @variable(model, -1.0 <= θ <= 1.0, start=θ_start)
    @variable(model, -3*2π <= ϕ <= 3*2π, start=ϕ_start)
    @variable(model, t₀, start=t₀_start)

    @NLobjective(model, Min, qfunc(x, y, z, θ, ϕ, t₀))

    optimize!(model);

    return Track(Direction(value(ϕ), value(θ)), Position(value(x), value(y), value(z)), value(t₀))
end

@with_kw struct SingleDURecoParams
    floor_distance::Int = 9
    Δt₋::Int = 10
    multiplicity::Int = 3
    min_hits::Int = 3
    Δt::Float64 = 10
end

DrWatson.default_prefix(s::SingleDURecoParams) = "SingleDUReco"

function startparams(SingleDUParams, du_hits::Vector{NeRCA.CalibratedHit})
    brightest_floor = NeRCA.most_frequent(h -> h.floor, du_hits)

    hits_on_brightest_floor = filter(h -> h.floor == brightest_floor, du_hits)
    thits_on_brightest_floor = filter(h -> h.triggered, hits_on_brightest_floor)

    if length(thits_on_brightest_floor) == 0
        z_closest = hits_on_brightest_floor[1].pos.z
        hit_time = mean([h.t for h in hits_on_brightest_floor])
    else
        closest_hit = thits_on_brightest_floor[1]
        z_closest = closest_hit.pos.z
        hit_time = closest_hit.t
    end

    # ϕ = azimuth(sum([-h.dir for h in du_hits]) ./ length(du_hits))

    SingleDUParams(10.0, hit_time, z_closest, -0.9, π, hit_time)
end

function single_du_fit(du_hits::Vector{NeRCA.CalibratedHit}, par::SingleDURecoParams; print_level=0)
    sort!(du_hits, by = h -> h.t)
    sort!(du_hits, by=h->h.dom_id)
    count_multiplicities!(du_hits)

    hit_pool = create_hit_pool(du_hits)

    hits₀ = unique(h->h.floor, nfoldhits(du_hits, par.Δt, par.multiplicity))
    if length(hits₀) < par.min_hits
        hits₀ = unique(h->h.floor, triggered(du_hits))
    end
    shits = select_hits(hits₀, hit_pool; Δt₋ = par.Δt₋, Δz = par.floor_distance)

    qfunc = SingleDUMinimiser(shits, filter(h->h.triggered, du_hits))

    model = Model(with_optimizer(Ipopt.Optimizer, print_level=print_level, tol=1e-3))

    register(model, :qfunc, 6, qfunc, autodiff=true)

    max_z = maximum([h.pos.z for h in shits]) + par.floor_distance
    min_z = minimum([h.pos.z for h in shits]) - par.floor_distance

    sdp₀ = startparams(SingleDUParams, du_hits)

    @variable(model, 1 <= d_closest <= 100, start=sdp₀.d)
    @variable(model, sdp₀.t - 1000 <= t_closest <= sdp₀.t + 1000, start=sdp₀.t)
    @variable(model, min_z <= z_closest <= max_z, start=sdp₀.z)
    @variable(model, -0.999 <= dir_z <= .999, start=sdp₀.dz)
    @variable(model, -4π <= ϕ₀ <= 4π, start=sdp₀.ϕ)
    @variable(model, t₀, start=sdp₀.t₀)

    @NLobjective(model, Min, qfunc(d_closest, t_closest, z_closest, dir_z, ϕ₀, t₀))

    optimize!(model);
    values = map(value, [d_closest, t_closest, z_closest, dir_z, ϕ₀, t₀])
    sdp = SingleDUParams(values...)
    ϕ = estimate_azimuth(sdp, shits, create_hit_pool(du_hits))
    sdp.ϕ = ϕ

    return ROyFit(sdp, sdp₀, qfunc(values...)/length(shits), shits, model)
end

function estimate_azimuth(
    sdp::SingleDUParams,
    direct_hits::Vector{CalibratedHit},
    hit_pool::Dict{Int, Vector{CalibratedHit}};
    Δt=20, n=N_SEAWATER
)
    ϕ = azimuth(sum([-h.dir for h in direct_hits]) ./ length(direct_hits))

    dᵧ, ccalc = make_cherenkov_calculator(sdp)

    late_hits = Vector{CalibratedHit}()
    for (floor, hits) in hit_pool
        t₀ = ccalc(first(hits).pos.z)
        for hit in hits
            if hit.t - t₀ > Δt
                push!(late_hits, hit)
            end
        end
    end

    if length(late_hits) == 0
        return ϕ
    end

    ϕₗ = azimuth(sum([-h.dir for h in late_hits]) ./ length(late_hits))
    if ϕ - ϕₗ < 0
        return ϕ - cos(1/n)
    else
        return ϕ + cos(1/n)
    end
end

function royfit_rbr(filename::AbstractString, detx::AbstractString, sparams::SingleDURecoParams)
end


text:
using MittagLeffler
using Base.Test

myapp(x,y) = abs(x-y) < 1e-10

@test myapp(mittleff(.5,.5,.5),1.5403698281390346)
@test myapp(mittleff(1.5,.5,.5),1.1448466286155243)
@test myapp(mittleff(2.3, .7 + 2. * im), 1.201890136368392 + 0.7895394560075035 * im)
@test myapp(mittleff(2.3, .7 + 0.2 * im) , (1.268233154873853 + 0.07914994421659409im))
@test mittleff(big".3",100.0) == big"8.721285946907744692995882256235296113802695745418015206361825134909144332670706e+2015816"

# allow alpha == beta == 1  --> exp(z)
@test myapp(mittleff(1,-2), 0.1353352832366127)

myder(f,x,h) = (f(x+h/2)-f(x-h/2))/h

@test myapp(myder(z -> mittleff(.4,z),.4,1e-5), mittleffderiv(.4,.4))
@test abs(myder(z -> mittleff(big".4",z),big".4",BigFloat(1//10^16)) - mittleffderiv(big".4",big".4")) < 1e-32

nothing


text:
using DimensionalAnalysis
using Base.Test

@test typeof( dimensionless()) == PhysicalQuantity{Int,0,0,0,0,0,0,0}
@test typeof( dimensionless())   <: Dimensionless
@test typeof( 3*dimensionless()) <: Dimensionless
@test typeof( dimensionless()/4 ) <: Dimensionless

@test dimensionless( 2 ) > dimensionless( 1 )
@test !( dimensionless( 1 ) > dimensionless( 2 ))
@test !( dimensionless( 2 ) > dimensionless( 2 ))

@test dimensionless( 2 ) >= dimensionless( 1 )
@test dimensionless( 2 ) >= dimensionless( 2 )
@test !( dimensionless( 1 ) >= dimensionless( 2 ))

@test dimensionless( 2 ) == dimensionless( 2 )
@test !( dimensionless( 1 ) == dimensionless( 2 ))
@test !( dimensionless( 2 ) == dimensionless( 1 ))

@test dimensionless( 1 ) <= dimensionless( 2 )
@test dimensionless( 1 ) <= dimensionless( 1 )
@test !( dimensionless( 2 ) <= dimensionless( 1 ))

@test dimensionless( 1 ) < dimensionless( 2 )
@test !( dimensionless( 1 ) < dimensionless( 1 ))
@test !( dimensionless( 2 ) < dimensionless( 1 ))

@test dimensionless( 2 ) > 1
@test 2 > dimensionless( 1 )
@test !( dimensionless( 1 ) > 2)
@test !( 1 > dimensionless( 2 ))
@test !( dimensionless( 1 ) > 1 )
@test !( 1 > dimensionless( 1 ))

@test dimensionless( 2 ) >= 1
@test !( 1 >= dimensionless( 2 ))

@test typeof( mass()) == PhysicalQuantity{Int,1,0,0,0,0,0,0}
@test typeof( mass()) <: Mass

@test typeof( length_()) == PhysicalQuantity{Int,0,1,0,0,0,0,0}
@test typeof( length_()) <: Length

@test typeof( time_()) == PhysicalQuantity{Int,0,0,1,0,0,0,0}
@test typeof( time_()) <: Time

@test typeof( temperature()) == PhysicalQuantity{Int,0,0,0,1,0,0,0}
@test typeof( temperature()) <: Temperature

@test typeof( intensity()) == PhysicalQuantity{Int,0,0,0,0,1,0,0}
@test typeof( intensity()) <: Intensity

@test typeof( current()) == PhysicalQuantity{Int,0,0,0,0,0,1,0}
@test typeof( current()) <: Current

@test typeof( amount()) == PhysicalQuantity{Int,0,0,0,0,0,0,1}
@test typeof( amount()) <: Amount

@test typeof( length_()*length_() ) <: Area
@test typeof( length_()*length_()*length_()) <: Volume
@test typeof( length_()*area()) <: Volume

@test typeof( density()) <: Density
@test typeof( mass()/volume() ) <: Density

@test typeof( speed()) <: Speed
@test typeof( length_()/time_()) <: Speed

@test typeof( acceleration())  <: Acceleration
@test typeof( speed()/time_()) <: Acceleration

@test typeof( force()) <: Force
@test typeof( mass()*acceleration()) <: Force

@test typeof( energy()) <: Energy
@test typeof( force()*length_()) <: Energy

@test typeof( power()) <: Power
@test typeof( energy()/time_()) <: Power

@test typeof( stress()) <: Stress
@test typeof( force()/area()) <: Stress

@test typeof( pressure()) <: Pressure
@test typeof( force()/area()) <: Pressure

@test typeof( diffusivity()) <: Diffusivity
@test typeof( area()/time_()) <: Diffusivity

@test typeof( kinematicViscosity()) <: KinematicViscosity
@test typeof( kinematicViscosity()) <: Diffusivity

@test typeof( dynamicViscosity()) <: DynamicViscosity
@test typeof( density()*diffusivity()) <: DynamicViscosity

@test typeof( energy()/area()) <: SurfaceEnergy
@test typeof( surfaceEnergy()) <: SurfaceEnergy

@test typeof( density()*speed()*length_()/dynamicViscosity()) <: Dimensionless
@test typeof( speed()*length_()/kinematicViscosity()) <: Dimensionless

let ρ = density(), v = speed(), d = length_(), μ = dynamicViscosity()
    Re = ReynoldsNumber( ρ, v, d, μ )
    @test typeof( Re ) <: Dimensionless
end

let v = speed(), d = length_(), ν = kinematicViscosity()
    Re = ReynoldsNumber( v, d, ν )
    @test typeof( Re ) <: Dimensionless
end

let ρ = density(), v = speed(), d = length_(), σ = surfaceEnergy()
    We = WeberNumber( ρ, v, d, σ )
    @test typeof( We ) <: Dimensionless
end

let ρ = density(),
    v = speed(),
    d = length_(),
    σ = surfaceEnergy(),
    We = WeberNumber( ρ, v, d, σ )
    We_m = ModifiedWeberNumber( ρ, v, d, σ )
    @test typeof( ρ )    <: Density
    @test typeof( d )    <: Length
    @test typeof( σ )    <: SurfaceEnergy
    @test typeof( We_m ) <: Dimensionless
    @test We/We_m == dimensionless( 12 )
end

max_stars_repo_name: sabjohnso/DimensionalAnalysis

max_issues_repo_name: sabjohnso/DimensionalAnalysis

max_forks_repo_name: sabjohnso/DimensionalAnalysis


text:
"""
    ArrayAllocators.ByteCalculators

Defines calculators for computing the number of bytes needed to allocate an array while detecting integer overflow.

# Examples

```julia
using ArrayAllocators.ByteCalculators

bc = CheckedMulByteCalculator{UInt8}(1024, 2048)
elsize(bc)
nbytes(bc)
```
"""
module ByteCalculators

export nbytes, elsize
export UnsafeByteCalculator, CheckedMulByteCalculator, WideningByteCalculator

"""
    AbstractByteCalculator

Parent abstract type for byte calculators, which calculate the total number of bytes of memory to allocate.
"""
abstract type AbstractByteCalculator{T} end

# Allow dimensions to be passed as individual arguments
function (::Type{B})(dims::Int...) where {T, B <: AbstractByteCalculator{T}}
    return B(dims)
end

function (::Type{B})(dims::AbstractUnitRange...) where {T, B <: AbstractByteCalculator{T}}
    return B(length.(dims))
end
function (::Type{B})(dims::NTuple{N, AbstractUnitRange}) where {N, T, B <: AbstractByteCalculator{T}}
    return B(length.(dims))
end
elsize(::AbstractByteCalculator{T}) where T = isbitstype(T) ? sizeof(T) : sizeof(Ptr)
nbytes(b::AbstractByteCalculator{T}) where T = elsize(b) * length(b)

"""
    UnsafeByteCalculator

Calculate number of bytes to allocate for an array without any integer overflow checking.
"""
struct UnsafeByteCalculator{T} <: AbstractByteCalculator{T}
    dims::Dims
    UnsafeByteCalculator{T}(dims::Dims) where T = new{T}(dims)
end
Base.length(b::UnsafeByteCalculator{T}) where T = prod(b.dims)

"""
    CheckedMulByteCalculator

Calculate the number of bytes by using `Base.checked_mul` to check if the product of the dimensions (length)
or the product of the length and the element size will cause an integer overflow.
"""
struct CheckedMulByteCalculator{T} <: AbstractByteCalculator{T}
    dims::Dims
    CheckedMulByteCalculator{T}(dims::Dims) where T = new{T}(dims)
end
function Base.length(b::CheckedMulByteCalculator)
    try
        return reduce(Base.checked_mul, b.dims)
    catch err
        if err isa OverflowError
            rethrow(OverflowError("The product of the dimensions results in integer overflow."))
        else
            rethrow(err)
        end
    end
end
function nbytes(b::CheckedMulByteCalculator{T}) where T
    element_size = sizeof(T)
    len = length(b)
    if len > typemax(typeof(element_size)) ÷ element_size
        throw(OverflowError("The product of array length and element size will cause an overflow."))
    end
    return element_size * len
end

"""
    WideningByteCalculator

Widens `eltype(Dims)`, Int in order to catch integer overflow.
"""
struct WideningByteCalculator{T} <: AbstractByteCalculator{T}
    dims::Dims
    WideningByteCalculator{T}(dims::Dims) where T = new{T}(dims)
end

function Base.length(bc::WideningByteCalculator)
    D = eltype(Dims)
    M = typemax(D)
    W = widen(D)

    len = reduce(bc.dims) do a,b
        c = W(a) * W(b)
        if c > M
            throw(OverflowError("The product of the dimensions results in integer overflow."))
        end
        return D(c)
    end

    return len
end
function nbytes(b::WideningByteCalculator{T}) where T
    nb = widen(sizeof(T)) * widen(length(b))
    if nb > typemax(eltype(Dims))
        throw(OverflowError("The product of array length and element size will cause an overflow."))
    end
    return nb
end

end # module ByteCalculators


text:
# By year

## ----------- Status by year (in college quality X; drop; grad)

raw_status_year(ds :: DataSettings; momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :progress, momentType,
        "q_dist_by_y.dat", ds);

function status_by_year(ds :: DataSettings)
    load_fct = mt -> begin
        # Read into DataFrameXY
        d = read_by_xy(raw_status_year(ds; momentType = mt));
        nr, nc = size(d.df);
        # Drop header column. Keep the rest.
        Matrix{Double}(d.df[1:nr, 2:nc]);
    end
    m, ses, cnts = choice_prob_from_xy(load_fct);
    @assert check_float_array(m, 0.0, 1.0);
    @assert size(m, 2) == n_colleges(ds) + 2  "Expecting columns for quality, drop, grad";
    return m, ses, cnts
end

## ----------  Fraction dropping out of 4y colleges by year
# Out of total entrants. First 3 years only.

# Return the fractions; std errors are not exactly right
function frac_drop_4y_by_year(ds :: DataSettings)
    # Rows are years, cols are colleges + dropouts + grads
    m, _, cnts = load_moment(ds, :statusByYear);
    @assert size(m, 2) == (n_colleges(ds) + 2)  "Invalid size: $(size(m))"

    # Keep 4y colleges
    idx4yV = (n_2year(ds) + 1) : n_colleges(ds);
    mass_tV = sum(m[:, idx4yV]; dims = 2);
    cnts_tV = sum(cnts[:, idx4yV]; dims = 2);
    # m = m[:, idx4yV];
    # cnts = cnts[:, idx4yV];

    # Only constructed for first 3 years
    T = 3;
    massEnter = mass_tV[1];     # sum(m[1, :]);
    @assert 0.3 < massEnter < 0.7
    fracDrop_tV = zeros(T);
    cnts_tV = cnts_tV[1 : T];
    # cnts_tV = zeros(Int, T);
    for t = 1 : T
        # Number in 4y in t
        # cnts_tV[t] = sum(cnts[t, :]);
        # There are no graduates, so everyone who leaves is a dropout.
        massDrop = mass_tV[t] - mass_tV[t+1];
        # massDrop = sum(m[t, :]) - sum(m[t+1, :]);
        @assert 0.0 < massDrop < 0.3;
        fracDrop_tV[t] = massDrop / massEnter;
    end
    @assert check_float_array(fracDrop_tV, 0.02, 0.15);

    # This is not exactly right, but the best we can do
    ses_tV = fracDrop_tV .* (1.0 .- fracDrop_tV) ./ (cnts_tV[1 : T] .^ 0.5);
    return fracDrop_tV, ses_tV, cnts_tV
end

## ----------  90th percentile of cumulative loans by year.
# No counts or std errors.
function cum_loans90_year(ds :: DataSettings)
    outV = cum_loans_year(ds; percentile = 90);
    cntV = zeros(Int, size(outV));
    sesV = zeros(Double, size(outV));
    return outV, sesV, cntV
end

function cum_loans_year(ds :: DataSettings; percentile = nothing)
    T = ds.Tmax;
    outV = zeros(Double, T);
    for t = 1 : T
        outV[t] = cum_loans_year(ds, t; percentile = percentile);
    end
    @assert all(outV .>= 0.0)  "Negative loans";
    return outV
end

"""
	$(SIGNATURES)

Read cumulative loans for one year. Optionally a percentile. Otherwise the mean.
"""
function cum_loans_year(ds :: DataSettings, t :: Integer; percentile = nothing)
    rf = raw_cum_loans_qual_year(ds, t; momentType = :mean, percentile = percentile);
    return Double(read_total(rf));
end

# -------------


text:
module InnerLoaders

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export innerloader

register_normal_element("q__inner__loading", context = @__MODULE__)

function innerloader(label::String = "",
                          state::Union{Symbol,Nothing} = nothing,
                          args...;
                          kwargs...)

  q__inner__loading(  args...; state = state, label = label, kw(kwargs)...)
end

end


text:
using HDF5
using Unitful

include("line.jl")
include("atmosphere.jl")
include("voronoi_utils.jl")

"""
    function write_arrays(x::AbstractArray, y::AbstractArray, z::AbstractArray,
                          fname::String)

Writes the arrays z, x, and y to a file with filename fname.
Arrays are written in columns [ row number ] [ x ] [ y ] [ z ]
"""
function write_arrays(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64},
                      fname::String)

    if length(z) != length(y) || length(y) != length(x)
        println("Wrong shapes of input data")
        exit()
    end

    open(fname, "w") do io
        for i in 1:length(z)
            println(io, "$i\t$(x[i])\t$(y[i])\t$(z[i])")
        end
    end
end

function write_arrays(x::Vector{<:Unitful.Length}, y::Vector{<:Unitful.Length}, z::Vector{<:Unitful.Length},
                      fname::String)

    x = ustrip.(x)
    y = ustrip.(y)
    z = ustrip.(z)

    if length(z) != length(y) || length(y) != length(x)
        println("Wrong shapes of input data")
        exit()
    end

    open(fname, "w") do io
        for i in 1:length(z)
            println(io, "$i\t$(x[i])\t$(y[i])\t$(z[i])")
        end
    end
end

function write_boundaries(z_min, z_max, x_min, x_max, y_min, y_max, fname::String)
    open(fname, "w") do io
        println(io, "z_min = $z_min")
        println(io, "z_max = $z_max")
        println(io, "x_min = $x_min")
        println(io, "x_max = $x_max")
        println(io, "y_min = $y_min")
        println(io, "y_max = $y_max")
    end
end

"""
    write_to_file(Array,
                  output_path::String)
Write 'Array' to the output file.
"""
function write_to_file(S_λ::Array{<:UnitsIntensity_λ, 4},
                       output_path::String)
    h5open(output_path, "r+") do file
        file["source_function"][:,:,:,:] = ustrip.(S_λ .|> u"kW*m^-2*nm^-1")
    end
end
function write_to_file(S_λ::Matrix{<:UnitsIntensity_λ},
                       output_path::String)
    h5open(output_path, "r+") do file
        file["source_function"][:,:] = ustrip.(S_λ .|> u"kW*m^-2*nm^-1")
    end
end

function write_to_file(populations::Array{<:NumberDensity, 4},
                       output_path::String)
    h5open(output_path, "r+") do file
        file["populations"][:,:,:,:] = ustrip.(populations .|> u"m^-3")
    end
end
function write_to_file(populations::Matrix{<:NumberDensity},
                       output_path::String)
    h5open(output_path, "r+") do file
        file["populations"][:,:] = ustrip.(populations .|> u"m^-3")
    end
end

function write_to_file(atmos::Atmosphere,
                       output_path::String;
                       ghost_cells=false)
    h5open(output_path, "r+") do file
        for field in fieldnames(Atmosphere)
            if String(field) in ["z", "x", "y"]
                if ghost_cells && String(field) in ["x", "y"]
                    file[String(field)][:] = ustrip.(getfield(atmos, field))[2:end-1]
                else
                    file[String(field)][:] = ustrip.(getfield(atmos, field))
                end
            else
                if ghost_cells
                    file[String(field)][:, :, :] = ustrip.(getfield(atmos, field))[:, 2:end-1, 2:end-1]
                else
                    file[String(field)][:, :, :] = ustrip.(getfield(atmos, field))
                end
            end
        end
    end
end

function write_to_file(sites::VoronoiSites,
                       output_path::String)

    h5open(output_path, "r+") do file
        fieldnames = ["positions", "temperature", "electron_density",
                      "hydrogen_populations", "velocity_z", "velocity_x", "velocity_y"]

        for field in fieldnames
            if field == "positions"
                file[field][:, :] = ustrip.(getfield(sites, Symbol(field)))
            else
                file[field][:] = ustrip.(getfield(sites, Symbol(field)))
            end
        end

        file["boundaries"][:] = ustrip.([sites.z_min, sites.z_max,
                                         sites.x_min, sites.x_max,
                                         sites.y_min, sites.y_max])

    end

end

function write_to_file(difference::Float64,
                       iteration::Int,
                       output_path::String)
    h5open(output_path, "r+") do file
       file["convergence"][iteration] = difference
    end
end

function write_to_file(n::Int,
                       field::String,
                       output_path::String)
    h5open(output_path, "r+") do file
       file[field][:] = n
    end
end

function write_to_file(line::HydrogenicLine,
                       output_path::String)
    h5open(output_path, "r+") do file
        file["wavelength"][:] = ustrip.(line.λ .|> u"nm")
        file["line_center"][:] = ustrip(line.λ0 |> u"nm")
    end
end

"""
    create_output_file(output_path::String, nλ::Int64, atmosphere_size::Tuple)
Initialise all output variables for the full atom mode.
"""
function create_output_file(output_path::String, nλ::Int, atmosphere_size::Tuple,
                            maxiter::Int)

    nz, nx, ny = atmosphere_size

    h5open(output_path, "w") do file
        write(file, "source_function", Array{Float64, 4}(undef, (nλ, nz, nx, ny)))
        write(file, "populations", Array{Float64, 4}(undef, (nz, nx, ny, 3)))

        write(file, "z", Vector{Float64}(undef, nz))
        write(file, "x", Vector{Float64}(undef, nx))
        write(file, "y", Vector{Float64}(undef, ny))

        write(file, "temperature", Array{Float64, 3}(undef, (nz, nx, ny)))
        write(file, "hydrogen_populations", Array{Float64, 3}(undef, (nz, nx, ny)))
        write(file, "electron_density", Array{Float64, 3}(undef, (nz, nx, ny)))

        write(file, "velocity_z", Array{Float64, 3}(undef, (nz, nx, ny)))
        write(file, "velocity_x", Array{Float64, 3}(undef, (nz, nx, ny)))
        write(file, "velocity_y", Array{Float64, 3}(undef, (nz, nx, ny)))

        write(file, "convergence", zeros(Float64, maxiter+1))

        write(file, "n_bb", Vector{Int}(undef, 1))
        write(file, "n_bf", Vector{Int}(undef, 1))

        write(file, "wavelength", Vector{Float64}(undef, nλ))
        write(file, "line_center", Vector{Float64}(undef, 1))

        write(file, "time", Vector{Float64}(undef, 1))
    end
end

"""
    create_output_file(output_path::String, nλ::Int64, atmosphere_size::Tuple)
Initialise all output variables for the full atom mode.
"""
function create_output_file(output_path::String, nλ::Int, n_sites::Int,
                            maxiter::Int)

    h5open(output_path, "w") do file
        write(file, "source_function", Matrix{Float64}(undef, (nλ, n_sites)))
        write(file, "populations", Matrix{Float64}(undef, (n_sites, 3)))

        write(file, "positions", Matrix{Float64}(undef, (3, n_sites)))

        write(file, "temperature", Vector{Float64}(undef, n_sites))
        write(file, "hydrogen_populations", Vector{Float64}(undef, n_sites))
        write(file, "electron_density", Vector{Float64}(undef, n_sites))

        write(file, "velocity_z", Vector{Float64}(undef, n_sites))
        write(file, "velocity_x", Vector{Float64}(undef, n_sites))
        write(file, "velocity_y", Vector{Float64}(undef, n_sites))

        write(file, "boundaries", Vector{Float64}(undef, 6))

        write(file, "convergence", zeros(Float64, maxiter+1))

        write(file, "n_bb", Vector{Int}(undef, 1))
        write(file, "n_bf", Vector{Int}(undef, 1))

        write(file, "wavelength", Vector{Float64}(undef, nλ))
        write(file, "line_center", Vector{Float64}(undef, 1))

        write(file, "time", Vector{Float64}(undef, 1))
    end
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
struct TapirEx{
    Basesize<:Union{Integer,Nothing},
    SIMD,
    Stoppable<:Union{Bool,Nothing},
    TaskGroup,
} <: FoldsBase.Executor
    basesize::Basesize
    simd::Val{SIMD}
    stoppable::Stoppable
    taskgroup::TaskGroup

    TapirEx(
        basesize::Basesize,
        simd::Val{SIMD},
        stoppable::Stoppable,
        taskgroup::TaskGroup,
    ) where {Basesize,SIMD,Stoppable,TaskGroup} =
        new{Basesize,SIMD,Stoppable,Core.Typeof(taskgroup)}(
            basesize,
            simd,
            stoppable,
            taskgroup,
        )
end

TapirEx(;
    basesize::Union{Integer,Nothing} = nothing,
    simd::SIMDFlag = Val(false),
    stoppable::Union{Bool,Nothing} = nothing,
    taskgroup = Tapir.taskgroup,
) = TapirEx(basesize, asval(Bool, simd), stoppable, taskgroup)

Transducers.transduce(xf, rf, init, xs, ex::TapirEx) = _transduce(xf, rf, init, xs, ex)

function _transduce(xform::Transducer, step::F, init, coll0, ex::TapirEx) where {F}
    (; basesize, simd, stoppable, taskgroup) = ex
    rf0 = _reducingfunction(xform, step; init = init)
    rf, coll = retransform(rf0, coll0)
    if stoppable === nothing
        stoppable = _might_return_reduced(rf, init, coll)
    end
    acc = @return_if_reduced _reduce(
        stoppable ? CancellableDACContext() : NoopDACContext(),
        maybe_usesimd(rf, simd),
        init,
        SizedReducible(
            coll,
            basesize === nothing ? amount(coll) ÷ Threads.nthreads() : basesize,
        ),
        taskgroup,
    )
    result = complete(rf, acc)
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
    end
    return result
end

function _reduce(ctx, rf::R, init::I, reducible::Reducible, taskgroup) where {R,I}
    if should_abort(ctx)
        return init
    end
    if issmall(reducible)
        acc = _reduce_basecase(rf, init, reducible)
        if acc isa Reduced
            cancel!(ctx)
        end
        return acc
    else
        left, right = _halve(reducible)
        fg, bg = splitcontext(ctx)
        @static if USE_TAPIR_OUTPUT
            Tapir.@output a0 b0
            Tapir.@sync taskgroup() begin
                Tapir.@spawn b0 = _reduce(bg, rf, init, right, taskgroup)
                a0 = _reduce(fg, rf, init, left, taskgroup)
            end
        else
            Tapir.@sync taskgroup() begin
                Tapir.@spawn $b0 = _reduce(bg, rf, init, right, taskgroup)
                $a0 = _reduce(fg, rf, init, left, taskgroup)
            end
        end
        a = @return_if_reduced a0
        should_abort(ctx) && return a  # slight optimization
        b0 isa Reduced && return combine_right_reduced(rf, a, b0)
        return combine(rf, a, b0)
    end
end


text:
const ad = AutoHOOT.autodiff

function construct_gradient!(net_sum::NetworkSum, innodes::Array, feed_dict::Dict)
  for net in net_sum.nodes
    inputs = ad.get_all_inputs(net)
    vars = [n for n in inputs if n in innodes]
    grads = ad.gradients(net, vars)
    for (grad, var) in zip(grads, vars)
      push!(feed_dict[var], grad)
    end
  end
end

# vector-jacobian product
function vjps(executor::Executor, vars, vector)
  node_dict = copy(executor.feed_dict)
  for t in vector
    update_dict!(node_dict, t)
  end
  # add vector to the executor
  @assert(length(vector) == length(executor))
  vec_nodes = [retrieve_key(node_dict, t) for t in vector]
  net_sums = [
    inner(net_sum, vnode) for (net_sum, vnode) in zip(executor.net_sums, vec_nodes)
  ]
  # get the gradient graph
  innodes = [retrieve_key(node_dict, t) for t in vars]
  network_sum_dict = Dict()
  for n in innodes
    network_sum_dict[n] = NetworkSum()
  end
  for net_sum in net_sums
    construct_gradient!(net_sum, innodes, network_sum_dict)
  end
  return Executor([network_sum_dict[n] for n in innodes], node_dict)
end

@non_differentiable vjps(executor::Executor, vars, vector)

batch_tensor_contraction(executor::Executor, vars...; kwargs...) = run(executor; kwargs...)

function ChainRulesCore.rrule(
  ::typeof(batch_tensor_contraction), executor::Executor, vars...; kwargs...
)
  function pullback(v)
    output = batch_tensor_contraction(vjps(executor, vars, v), vars...; kwargs...)
    return (NoTangent(), NoTangent(), Tuple(output)...)
  end
  return batch_tensor_contraction(executor, vars...; kwargs...), pullback
end

"""Perform a batch of tensor contractions, each one defined by a tensor network.
Parameters
----------
networks: An array of networks. Each network is represented by an array of tensors
vars: the tensors to take derivative of
Returns
-------
A list of tensors representing the contraction outputs of each network.
"""
function batch_tensor_contraction(
  networks::Vector{<:Vector{<:AbstractTensor}}, vars...; kwargs...
)
  return batch_tensor_contraction(Executor(networks), vars...; kwargs...)
end

function batch_tensor_contraction(trees::Vector{SubNetwork}, vars...; kwargs...)
  return batch_tensor_contraction(Executor(trees), vars...; kwargs...)
end

function batch_tensor_contraction(
  networks::Vector{<:Vector{<:AbstractTensor}}, cache::NetworkCache, vars...; kwargs...
)
  return batch_tensor_contraction(Executor(networks, cache), vars...; kwargs...)
end

function batch_tensor_contraction(
  trees::Vector{SubNetwork}, cache::NetworkCache, vars...; kwargs...
)
  return batch_tensor_contraction(Executor(trees, cache), vars...; kwargs...)
end

function batch_tensor_contraction(
  tensortype, networks::Vector{Vector{ITensor}}, vars...; kwargs...
)
  networks, vars = abstract_network(tensortype, networks, vars; kwargs...)
  out = batch_tensor_contraction(Executor(networks), vars...; kwargs...)
  return out
end

function batch_tensor_contraction(
  tensortype, trees::Vector{SubNetwork}, vars...; optimize=true, kwargs...
)
  trees, vars = abstract_network(tensortype, trees, vars; kwargs...)
  out = batch_tensor_contraction(Executor(trees; optimize=optimize), vars...; kwargs...)
  return out
end


text:
using LLVM.Interop
using Core: LLVMPtr

function create_compute_context()
    ctx    = create_some_context()
    device = first(devices(ctx))
    queue  = CmdQueue(ctx)
    return (device, ctx, queue)
end

opencl_version(obj :: CLObject) = api.parse_version(obj[:version])
opencl_version(c :: Context)  = opencl_version(first(devices(c)))
opencl_version(q :: CmdQueue) = opencl_version(q[:device])

const _versionDict = Dict{Ptr{Nothing}, VersionNumber}()

_deletecached!(obj :: CLObject) = delete!(_versionDict, pointer(obj))

function check_version(obj :: CLObject, version :: VersionNumber)
    version <= get!(_versionDict, pointer(obj)) do
        opencl_version(obj)
    end
end

"""
Format string using dict-like variables, replacing all accurancies of
`%(key)` with `value`.

Example:
    s = "Hello, %(name)"
    format(s, name="Tom")  ==> "Hello, Tom"
"""
function format(s::String; vars...)
    for (k, v) in vars
        s = replace(s, "%($k)" => v)
    end
    s
end

function build_kernel(ctx::Context, program::String,
                      kernel_name::String; vars...)
    src = format(program; vars...)
    p = Program(ctx, source=src)
    build!(p)
    return Kernel(p, kernel_name)
end

# cache for kernels; dict of form `(program_file, kernel_name, vars) -> kernel`
const CACHED_KERNELS = Dict{Tuple{String, String, Dict}, Kernel}()

function get_kernel(ctx::Context, program_file::String,
                    kernel_name::String; vars...)
    key = (program_file, kernel_name, Dict(vars))
    if in(key, keys(CACHED_KERNELS))
        return CACHED_KERNELS[key]
    else
        kernel = build_kernel(ctx, Base.read(program_file, String), kernel_name; vars...)
        CACHED_KERNELS[key] = kernel
        return kernel
    end
end

min_v11(elem) = check_version(elem, v"1.1")
min_v12(elem) = check_version(elem, v"1.2")
min_v20(elem) = check_version(elem, v"2.0")
min_v21(elem) = check_version(elem, v"2.1")
min_v22(elem) = check_version(elem, v"2.2")

const opencl_builtins = String["printf"]

# OpenCL functions need to be mangled according to the C++ Itanium spec. We implement a very
# limited version of that spec here, just enough to support OpenCL built-ins.
#
# This macro also keeps track of called builtins, generating `ccall("extern...", llvmcall)`
# expressions for them (so that we can exclude them during IR verification).
macro builtin_ccall(name, ret, argtypes, args...)
    @assert Meta.isexpr(argtypes, :tuple)
    argtypes = argtypes.args

    function mangle(T::Type)
        if T == Cint
            "i"
        elseif T == Cuint
            "j"
        elseif T == Clong
            "l"
        elseif T == Culong
            "m"
        elseif T == Cshort
            "s"
        elseif T == Cushort
            "t"
        elseif T == Cchar
            "c"
        elseif T == Cuchar
            "h"
        elseif T == Cfloat
            "f"
        elseif T == Cdouble
            "d"
        elseif T <: LLVMPtr
            elt, as = T.parameters

            # mangle address space
            ASstr = if as == AS.Global
                "CLglobal"
            #elseif as == AS.Global_device
            #    "CLdevice"
            #elseif as == AS.Global_host
            #    "CLhost"
            elseif as == AS.Local
                "CLlocal"
            elseif as == AS.Constant
                "CLconstant"
            elseif as == AS.Private
                "CLprivate"
            elseif as == AS.Generic
                "CLgeneric"
            else
                error("Unknown address space $AS")
            end

            # encode as vendor qualifier
            ASstr = "U" * string(length(ASstr)) * ASstr

            # XXX: where does the V come from?
            "P" * ASstr * "V" * mangle(elt)
        else
            error("Unknown type $t")
        end
    end

    # C++-style mangling; very limited to just support these intrinsics
    # TODO: generalize for use with other intrinsics? do we need to mangle those?
    mangled = "_Z$(length(name))$name"
    for t in argtypes
        # with `@eval @builtin_ccall`, we get actual types in the ast, otherwise symbols
        t = (isa(t, Symbol) || isa(t, Expr)) ? eval(t) : t
        mangled *= mangle(t)
    end

    push!(opencl_builtins, mangled)
    esc(quote
        @typed_ccall($mangled, llvmcall, $ret, ($(argtypes...),), $(args...))
    end)
end

## device overrides

# local method table for device functions
@static if isdefined(Base.Experimental, Symbol("@overlay"))
Base.Experimental.@MethodTable(method_table)
else
const method_table = nothing
end

# list of overrides (only for Julia 1.6)
const overrides = quote end

macro device_override(ex)
    code = quote
        $GPUCompiler.@override($method_table, $ex)
    end
    if isdefined(Base.Experimental, Symbol("@overlay"))
        return esc(code)
    else
        push!(overrides.args, code)
        return
    end
end

macro device_function(ex)
    ex = macroexpand(__module__, ex)
    def = ExprTools.splitdef(ex)

    # generate a function that errors
    def[:body] = quote
        error("This function is not intended for use on the CPU")
    end

    esc(quote
        $(ExprTools.combinedef(def))
        @device_override $ex
    end)
end