text:
using Distributed, DataFrames, CSV
using MINLPLib, JuMP, SCIP

new_lib = "ImprovedCompMult"

# build library
build_lib = false
instance_names = ["bearing", "ex6_2_10", "ex6_2_10", "ex6_2_11", "ex6_2_12", "ex6_2_13", "ex7_2_4", "ex7_3_1", "ex7_3_2", "ex14_1_8", "ex14_1_9"]
if build_lib
    source_lib = "global"
    for n in instance_names
        MINLPLib.add_to_lib(new_lib, source_lib, n)
    end
end

# solution handler
struct TrialSummary
    was_error::Bool
    solver_name::String
    instance_name::String
    term_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    dual_status::MOI.ResultStatusCode
    obj_bound::Float64
    obj_value::Float64
    solve_time::Float64
end
TrialSummary() = TrialSummary(true, "", "", MOI.OTHER_ERROR, MOI.OTHER_RESULT_STATUS, MOI.OTHER_RESULT_STATUS, -Inf, -Inf, 0.0)

function benchmark_problem(d, lib, sname::String)
    m = fetch_model(lib, sname)
    set_optimizer(m, d)

    set_time_limit_sec(m, 1000.0)
    set_silent(m)

    JuMP.optimize!(m)
    p = primal_status(m)
    ds = dual_status(m)
    t = termination_status(m)
    v = objective_value(m)
    b = objective_bound(m)
    n = solve_name(m)
    s = solve_time(m)
    return TrialSummary(false,n,sname,t,p,ds,v,b,s)
end

trial_summaries = TrialSummary[]

# Need to assign tolerance here (since tolerances aren't standardized among solvers)
function build_scip()
    m = SCIP.Optimizer(limits_gap=1E-3,    # absolute tolerance
                       limits_absgap=1E-3  # relative tolerance
                       )
    return m
end
append!(trial_summaries, pmap(x->benchmark_problem(build_scip, new_lib, x), instance_names; on_error=ex->TrialSummary()))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


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

using Markdown
using InteractiveUtils

# ╔═╡ b107d7f9-90c3-4a9f-9b55-0e1dc323ee56
using PlutoUI

# ╔═╡ ad6469ef-68a5-4819-bf98-85d06d37b960
using CSV

# ╔═╡ f31b4062-c763-4647-84ef-b3db25406fc8
using HTTP

# ╔═╡ 6cab610d-e713-4b08-8a5b-5eab8820732f
using DataFrames

# ╔═╡ 6fc301b3-18f2-44e8-91d2-268e10c3a0e6
using BenchmarkTools

# ╔═╡ 6b4a8fd1-60f5-4760-a453-86a75a201a5f
md"""# Notebok sobre Ciencia de Datos

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

# ╔═╡ ef6340c4-f4f7-11eb-30f9-599dcebb98be
PlutoUI.TableOfContents(title="🔬 Ciencia de datos 🔭")

# ╔═╡ d24051dc-5efa-48c8-b65b-c756ba43b864
md"# Manejo de datos"

# ╔═╡ b1a67905-09f3-41aa-9a87-362d051fc69f
md"## CSV.jl
- [GitHub](https://github.com/JuliaData/CSV.jl)
- [JuliaData](https://csv.juliadata.org/stable/)"

# ╔═╡ 6069ba43-4a4d-4c59-a465-098385eb6e80
md"### Leyendo csv
Usemos la función `CSV.File` para leer archivos `.csv`. Podemos hacerlo desde archivos locales:"

# ╔═╡ a3c21e5d-a030-4cf3-aed3-e7100a885d37
hotel = CSV.File("Julia-IIMAS-v1.0/datos/hotel_bookings.csv")

# ╔═╡ 6baca64d-6dcb-44d7-b03d-69c0322631fe
md"O también podemos leer desde una página con la paquetería `HTTP.jl`"

# ╔═╡ e32d05c9-5896-4435-8f4c-15f887d6b256
begin
	url = "https://raw.githubusercontent.com/OscarAlvaradoM/Pokemon/master/datos/Disp8_completo.csv"
	pokemon = CSV.File(HTTP.get(url).body)
end

# ╔═╡ bf9cd655-eafe-41d2-ab06-f0a53cde9455
pokemon.names

# ╔═╡ f653a134-d37a-4959-bf81-40c2d20d4a87
pokemon.cols

# ╔═╡ 7ea3c9a3-a369-4ee9-908d-ed8b7d1ae3cb
md"### Delimitador"

# ╔═╡ 36b56e8b-dde9-4fcb-a3f9-a93fb3101cdf
md"""
	col1::col2
	1::2
	3::4
```Julia
CSV.File(file; delim="::")
```
"""

# ╔═╡ fcd2724d-aa29-4e7b-8d39-818463e7ceb5
md"### Sin encabezado"

# ╔═╡ 6165b59f-2c9c-44d8-b914-a307f99a5a48
md"""
	1,2,3
	4,5,6
	7,8,9
```Julia
CSV.File(file; header=false)
CSV.File(file; header=["col1", "col2", "col3"])
CSV.File(file; header=[:col1, :col2, :col3])
```
"""

# ╔═╡ c36c1741-21fe-4758-a82d-dbf3439830cd
md"### Datarow"

# ╔═╡ 51e1ee64-9541-4ee2-9499-748811dacac9
md"""
	col1,col2,col3
	metadata1,metadata2,metadata3
	extra1,extra2,extra3
	1,2,3
	4,5,6
	7,8,9

```Julia
CSV.File(file; datarow=4)
CSV.File(file; skipto=4)
```
"""

# ╔═╡ e2da47f5-d1bf-4060-8113-5fa2a7de56ff
md"### Leyendo lotes"

# ╔═╡ 4ece0521-051e-475c-bec9-10f38beb009d
md"""
	col1,col2,col3
	1,2,3
	4,5,6
	7,8,9
	10,11,12
	13,14,15
	16,17,18
	19,20,21

```Julia
CSV.File(file; limit=3)
CSV.File(file; skipto=4, limit=1)
CSV.File(file; skipto=7, footerskip=1)
```
"""

# ╔═╡ 236f6a7d-86a5-4bbc-a0d8-49f06e1714bf
md"### Reenglones comentados"

# ╔═╡ 66e15e9e-f283-487d-8e73-ae5ec23edb54
md"""
	col1,col2,col3
	# this row is commented and we'd like to ignore it while parsing
	1,2,3
	4,5,6

```Julia
CSV.File(file; comment="#")
CSV.File(file; datarow=3)
```
"""

# ╔═╡ 72e80a6f-8d30-4a85-a154-d2fbfc653421
md"### Valores nulos"

# ╔═╡ 0b5f90c0-7b51-484a-a023-0f533167bf7b
md"""
	code,age,score
	0,21,3.42
	1,42,6.55
	-999,81,NA
	-999,83,NA

```Julia
CSV.File(file; missingstring="-999")
CSV.File(file; missingstrings=["-999", "NA"])
```
"""

# ╔═╡ 720f5fe8-0616-4480-be24-b81764c4361f
md"### Fechas"

# ╔═╡ 7cbe09ae-5dcc-4c42-b824-bd401ed7f18b
md"""
	code,date
	0,2019/01/01
	1,2019/01/02

```Julia
CSV.File(file; dateformat="yyyy/mm/dd")
```
"""

# ╔═╡ 6b8336a6-fac6-4042-84e5-2cf4a0bf83b7
md"### Asignando tipos"

# ╔═╡ 9d3653a3-a617-4694-bc87-72d2606245c9
md"""
	col1,col2,col3
	1,2,3
	4,5,invalid
	6,7,8

```Julia
CSV.File(file; types=Dict(3 => Int))
CSV.File(file; types=Dict(:col3 => Int))
CSV.File(file; types=Dict("col3" => Int))
CSV.File(file; types=[Int, Int, Int])
CSV.File(file; types=[Int, Int, Int], silencewarnings=true)
CSV.File(file; types=[Int, Int, Int], strict=true)
```
"""

# ╔═╡ 7f43e8c7-7ac8-4aa2-acae-87e0c3814da1
md"### Seleccionar o eliminar columnas de un archivo"

# ╔═╡ ffb89c6b-cdd1-4a64-8ad7-b512373dfa1a
md"""
	a,b,c
	1,2,3
	4,5,6
	7,8,9

```Julia
# select
CSV.File(file; select=[1, 3])
CSV.File(file; select=[:a, :c])
CSV.File(file; select=["a", "c"])
CSV.File(file; select=[true, false, true])
CSV.File(file; select=(i, nm) -> i in (1, 3))
# drop
CSV.File(file; drop=[2])
CSV.File(file; drop=[:b])
CSV.File(file; drop=["b"])
CSV.File(file; drop=[false, true, false])
CSV.File(file; drop=(i, nm) -> i == 2)
```
"""

# ╔═╡ acf3fc88-d7cd-4bf9-ba4b-dd166984d7aa
md"## DataFrames.jl
- [GitHub](https://github.com/JuliaData/DataFrames.jl)
- [JuliaData](https://dataframes.juliadata.org/stable/#)"

# ╔═╡ 92056bbd-8dde-4b38-8780-02e248c6cfd3
md"### Transformando de un CSV.File a DataFrame"

# ╔═╡ fdbd56bb-4a02-436a-9ee3-552b373547c7
df_poke = CSV.read(HTTP.get("https://raw.githubusercontent.com/OscarAlvaradoM/Pokemon/master/datos/Disp8_completo.csv").body, DataFrame)

# ╔═╡ ccd0662b-455d-453f-be8c-d0a998910a7f
df_pokemon = DataFrame(pokemon)

# ╔═╡ 97d67ce7-cfb5-47fd-9ea3-272ded84e9fd
methodswith(typeof(df_pokemon))

# ╔═╡ 21f1d7dd-a865-4e8b-ab03-cd81b91f5f29
names(df_pokemon)

# ╔═╡ 6dc1461e-ad0e-4020-825a-a696de245ca8
nrow(df_pokemon)

# ╔═╡ 7dd8b809-0870-4945-97d2-2a655d0fcc63
md"### Creando DataFrames desde 0"

# ╔═╡ fd76cca7-48c0-474a-b826-3be1d7c54c14
DataFrame(A=1:3, B=5:7, fixed=1)

# ╔═╡ a5f4bad2-2919-4241-9178-e944e42d69a5
DataFrame("Edad" => ["?", 23, 24], "Nombre" => ["Oscar", "Mario", "Anuar"])

# ╔═╡ b2fbf649-32b6-4dd8-894c-9ca2e2fa1e9c
begin
	dict_edades1 = Dict("Edad" => ["?", 23, 24], "Nombre" => ["Oscar", "Mario", "Anuar"])
	DataFrame(dict_edades1)
end

# ╔═╡ 8b1d8c99-7fa1-4054-ae7d-52b01a97b0b3
begin
	dict_edades2 = Dict(:Edad => ["?", 23, 24], :Nombre => ["Oscar", "Mario", "Anuar"])
	DataFrame(dict_edades2)
end

# ╔═╡ 76dbe6fa-7426-4ef1-874a-7463602e6c7c
DataFrame((a=[1, 2], b=[3, 4]))

# ╔═╡ bbae294a-01f4-47a3-9cfb-99e8e0978796
DataFrame([(a=1, b=0), (a=2, b=0)])

# ╔═╡ 064bea34-0890-4d27-8c5c-63d7f66e10ac
DataFrame([1 0; 2 0], ["A", "B"])

# ╔═╡ 35fef6ab-eaf7-4e39-8424-1e77ef29a4f1
DataFrame([1 0; 2 0], :auto)

# ╔═╡ e13caf23-fe50-456a-a7f0-1df74580e6c6
md"### Operaciones básicas sobre DataFrames"

# ╔═╡ 2d0dd71a-d577-4d56-bf42-31ce47d0bebe
names(df_pokemon)

# ╔═╡ 41424723-6627-4462-9e31-36ef20933230
md"Esto no es una copia"

# ╔═╡ 816b06b0-82b9-45ce-8cbb-22061363c936
df_pokemon.Nombre

# ╔═╡ 8784abe2-688d-401f-bce0-0ad02b4b7e76
md"Esto es una copia:"

# ╔═╡ efca9fcb-2378-49ad-b82b-a528674944da
df_pokemon[:, :Nombre]

# ╔═╡ 0cd53836-8853-46da-86f4-2a11bee5bed8
md"Y esto ¿qué es?"

# ╔═╡ 0eb90afd-4e1c-4bf2-b8a4-b7241ed7d042
df_pokemon[!, :Nombre]

# ╔═╡ 7c787df1-3f7e-4057-91e2-db9b0a42b11f
df_pokemon.Nombre === df_pokemon[!, :Nombre]

# ╔═╡ e426364d-dc2d-45a7-8e3f-a7d63424433e
df_pokemon.Nombre === df_pokemon[:, :Nombre]

# ╔═╡ d96f14f8-dc5c-4adf-899e-5e9497e9bba1
eltype.(eachcol(df_pokemon))

# ╔═╡ 79bf46e7-96c3-488d-9ced-e35798d6ed09
size(df_pokemon)

# ╔═╡ d6b8b7f2-2eff-4201-b53f-c10848942f9b
nrow(df_pokemon)

# ╔═╡ dfe4c5c0-d51d-4aae-abe7-604619a16e17
ncol(df_pokemon)

# ╔═╡ 22b933e1-9763-4922-b72d-a9f59477aaaf
describe(df_pokemon, cols=1:10)

# ╔═╡ 4e6a6592-72dd-4e59-9cb8-e4aa75b82c3f
mapcols(Altura_m -> Altura_m.^2, df_pokemon)

# ╔═╡ dc16d155-379f-4d2a-8e01-825c384c117c
df_pokemon.Altura_m

# ╔═╡ f413974c-b210-4d9e-bc3f-c5b0d3f28e69
first(df_pokemon, 8)

# ╔═╡ 5f96e230-7e61-4a1e-962f-8ca5f0d0be60
last(df_pokemon, 5)

# ╔═╡ 466e107d-1c7a-4a15-a08b-c9120b84a917
df_pokemon[1:5, [:Nombre, :Clasificación]]

# ╔═╡ b7935362-ac6b-47a8-9106-cc2beaac70a8
df_pokemon[[1, 6, end], :]

# ╔═╡ 8d0b9978-b04a-436d-bf13-6946b043e48d
@view df_pokemon[end:-1:1, [:Nombre, :Generación]]

# ╔═╡ 38452d24-76bd-4d75-9f9b-39d5d0d2dae8
@btime @view $df_pokemon[1:end-1, 1:end-1]

# ╔═╡ 68733a99-4d7e-4256-a168-2032de50ab80
@btime $df_pokemon[1:end-1, 1:end-1];

# ╔═╡ e0968f67-1995-45d1-81e5-f6fb875312d5
df_pokemon.Contra_acero .= 1.

# ╔═╡ e51c2334-fd90-4847-9877-05ed2390d45a
insertcols!(df_pokemon, 1, :Juego => "Pokemon")

# ╔═╡ 432e6d0a-7671-4b3c-a8c4-0508a418cee7
df_pokemon[:, Not(:Evolución)]

# ╔═╡ 7e693874-9155-4c2f-8b4c-3543454fbb56
df_pokemon[:, Between(:Juego, :Generación)]

# ╔═╡ 7cc17b78-8f2a-49de-a7a6-99ab6d89b0e0
df_pokemon[:, Cols(:Altura_m, Between(:Juego, :Generación))]

# ╔═╡ ab72fb5f-90ac-446c-88d6-938852f3816e
md"Haciendo algunos joins"

# ╔═╡ 9c851f5a-4c8f-4388-86a1-ff6a3d83c8fe
people = DataFrame(ID = [20, 40], Name = ["John Doe", "Jane Doe"])

# ╔═╡ 6f6459c8-8242-4340-a588-0512eff5d5c3
jobs = DataFrame(ID = [20, 60], Job = ["Lawyer", "Doctor"])

# ╔═╡ db565011-1563-41cc-84d0-e14c8237aa8a
innerjoin(people, jobs, on = :ID)

# ╔═╡ f518b097-a6f4-4bf5-ba36-0088f5fc6b94
leftjoin(people, jobs, on = :ID)

# ╔═╡ 98693a39-71de-4e60-bf66-42e30618e254
rightjoin(people, jobs, on = :ID)

# ╔═╡ f1ea0441-6910-4704-906a-3e6cf05b4105
df = DataFrame(i = 1:5,x = [missing, 4, missing, 2, 1],y = [missing, missing, "c", "d", "e"])

# ╔═╡ a0d7c872-8c54-4eed-9a01-1ca92c2d6a28
dropmissing(df)

# ╔═╡ 65729744-b221-488d-b1ad-85981489916c
dropmissing(df, :x)

# ╔═╡ c0b3f6a7-341c-435a-b528-3a190f5c500c
md"Una buena comparación con Pandas y R:

[DataFrames (JuliaData)](https://dataframes.juliadata.org/stable/man/comparisons/)"

# ╔═╡ 185015b8-d252-4beb-96f0-4d2483f79344
md"# Tarea
Buscar un conjunto de datos de tu interés con el cuál trabajar para aprendizaje automatizado, leerlo y analizarlo con julia"

# ╔═╡ e9fcbd2a-6c9b-4e07-bcae-b990ebc1a3da
md"# Aprendizaje Automatizado"

# ╔═╡ 042c51dc-80d8-4c19-9c87-4e3e83de5e19
md"## MLJ.jl
- [GitHub](https://github.com/alan-turing-institute/MLJ.jl)
- [Alan Turing Institute](https://alan-turing-institute.github.io/MLJ.jl/dev/)"

# ╔═╡ 44899540-8148-4b83-a88a-9b8b4d73ce82

# ╔═╡ 18624127-8d76-4b10-ab15-5e7cd4caea79
md"## ScikitLearn.jl
- [Docs](https://scikitlearnjl.readthedocs.io/en/latest/)"

# ╔═╡ caa805f8-7ece-429b-b63b-0542dbe159c1

# ╔═╡ 4247ed0b-7722-4ea9-a243-a48644bc65ee
md"# Aprendizaje profundo"

# ╔═╡ b0bc0de6-a0ab-49c6-a425-6df330e877d4
md"## Flux.jl
- [GitHub](https://github.com/FluxML/Flux.jl)
- [fluxml](https://fluxml.ai/)"

# ╔═╡ 798ee25d-3767-462d-b7d8-522f20754f78

# ╔═╡ 4ec0cfec-f200-4063-9cf0-5979cd369af9
md"# Big Data"

# ╔═╡ 03a434ee-d037-4607-b346-94b31e5a0fa5
md"## Spark.jl
- [GitHub](https://github.com/dfdx/Spark.jl)"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
BenchmarkTools = "~1.1.1"
CSV = "~0.8.5"
DataFrames = "~1.2.2"
HTTP = "~0.9.13"
PlutoUI = "~0.7.9"
"""

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

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

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

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

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

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

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

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

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

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

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

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

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

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

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

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

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

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

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

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

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

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

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

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

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

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

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

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

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

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

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

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "a3a337914a035b2d59c9cbe7f1a38aaba1265b02"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.6"

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─6b4a8fd1-60f5-4760-a453-86a75a201a5f
# ╠═b107d7f9-90c3-4a9f-9b55-0e1dc323ee56
# ╠═ef6340c4-f4f7-11eb-30f9-599dcebb98be
# ╟─d24051dc-5efa-48c8-b65b-c756ba43b864
# ╟─b1a67905-09f3-41aa-9a87-362d051fc69f
# ╠═ad6469ef-68a5-4819-bf98-85d06d37b960
# ╟─6069ba43-4a4d-4c59-a465-098385eb6e80
# ╠═a3c21e5d-a030-4cf3-aed3-e7100a885d37
# ╟─6baca64d-6dcb-44d7-b03d-69c0322631fe
# ╠═f31b4062-c763-4647-84ef-b3db25406fc8
# ╠═e32d05c9-5896-4435-8f4c-15f887d6b256
# ╠═bf9cd655-eafe-41d2-ab06-f0a53cde9455
# ╠═f653a134-d37a-4959-bf81-40c2d20d4a87
# ╟─7ea3c9a3-a369-4ee9-908d-ed8b7d1ae3cb
# ╟─36b56e8b-dde9-4fcb-a3f9-a93fb3101cdf
# ╟─fcd2724d-aa29-4e7b-8d39-818463e7ceb5
# ╟─6165b59f-2c9c-44d8-b914-a307f99a5a48
# ╟─c36c1741-21fe-4758-a82d-dbf3439830cd
# ╟─51e1ee64-9541-4ee2-9499-748811dacac9
# ╟─e2da47f5-d1bf-4060-8113-5fa2a7de56ff
# ╟─4ece0521-051e-475c-bec9-10f38beb009d
# ╟─236f6a7d-86a5-4bbc-a0d8-49f06e1714bf
# ╟─66e15e9e-f283-487d-8e73-ae5ec23edb54
# ╟─72e80a6f-8d30-4a85-a154-d2fbfc653421
# ╟─0b5f90c0-7b51-484a-a023-0f533167bf7b
# ╟─720f5fe8-0616-4480-be24-b81764c4361f
# ╟─7cbe09ae-5dcc-4c42-b824-bd401ed7f18b
# ╟─6b8336a6-fac6-4042-84e5-2cf4a0bf83b7
# ╟─9d3653a3-a617-4694-bc87-72d2606245c9
# ╟─7f43e8c7-7ac8-4aa2-acae-87e0c3814da1
# ╟─ffb89c6b-cdd1-4a64-8ad7-b512373dfa1a
# ╟─acf3fc88-d7cd-4bf9-ba4b-dd166984d7aa
# ╠═6cab610d-e713-4b08-8a5b-5eab8820732f
# ╟─92056bbd-8dde-4b38-8780-02e248c6cfd3
# ╠═fdbd56bb-4a02-436a-9ee3-552b373547c7
# ╠═ccd0662b-455d-453f-be8c-d0a998910a7f
# ╠═97d67ce7-cfb5-47fd-9ea3-272ded84e9fd
# ╠═21f1d7dd-a865-4e8b-ab03-cd81b91f5f29
# ╠═6dc1461e-ad0e-4020-825a-a696de245ca8
# ╟─7dd8b809-0870-4945-97d2-2a655d0fcc63
# ╠═fd76cca7-48c0-474a-b826-3be1d7c54c14
# ╠═a5f4bad2-2919-4241-9178-e944e42d69a5
# ╠═b2fbf649-32b6-4dd8-894c-9ca2e2fa1e9c
# ╠═8b1d8c99-7fa1-4054-ae7d-52b01a97b0b3
# ╠═76dbe6fa-7426-4ef1-874a-7463602e6c7c
# ╠═bbae294a-01f4-47a3-9cfb-99e8e0978796
# ╠═064bea34-0890-4d27-8c5c-63d7f66e10ac
# ╠═35fef6ab-eaf7-4e39-8424-1e77ef29a4f1
# ╟─e13caf23-fe50-456a-a7f0-1df74580e6c6
# ╠═2d0dd71a-d577-4d56-bf42-31ce47d0bebe
# ╟─41424723-6627-4462-9e31-36ef20933230
# ╠═816b06b0-82b9-45ce-8cbb-22061363c936
# ╟─8784abe2-688d-401f-bce0-0ad02b4b7e76
# ╠═efca9fcb-2378-49ad-b82b-a528674944da
# ╟─0cd53836-8853-46da-86f4-2a11bee5bed8
# ╠═0eb90afd-4e1c-4bf2-b8a4-b7241ed7d042
# ╠═7c787df1-3f7e-4057-91e2-db9b0a42b11f
# ╠═e426364d-dc2d-45a7-8e3f-a7d63424433e
# ╠═d96f14f8-dc5c-4adf-899e-5e9497e9bba1
# ╠═79bf46e7-96c3-488d-9ced-e35798d6ed09
# ╠═d6b8b7f2-2eff-4201-b53f-c10848942f9b
# ╠═dfe4c5c0-d51d-4aae-abe7-604619a16e17
# ╠═22b933e1-9763-4922-b72d-a9f59477aaaf
# ╠═4e6a6592-72dd-4e59-9cb8-e4aa75b82c3f
# ╠═dc16d155-379f-4d2a-8e01-825c384c117c
# ╠═f413974c-b210-4d9e-bc3f-c5b0d3f28e69
# ╠═5f96e230-7e61-4a1e-962f-8ca5f0d0be60
# ╠═466e107d-1c7a-4a15-a08b-c9120b84a917
# ╠═b7935362-ac6b-47a8-9106-cc2beaac70a8
# ╠═8d0b9978-b04a-436d-bf13-6946b043e48d
# ╠═6fc301b3-18f2-44e8-91d2-268e10c3a0e6
# ╠═38452d24-76bd-4d75-9f9b-39d5d0d2dae8
# ╠═68733a99-4d7e-4256-a168-2032de50ab80
# ╠═e0968f67-1995-45d1-81e5-f6fb875312d5
# ╠═e51c2334-fd90-4847-9877-05ed2390d45a
# ╠═432e6d0a-7671-4b3c-a8c4-0508a418cee7
# ╠═7e693874-9155-4c2f-8b4c-3543454fbb56
# ╠═7cc17b78-8f2a-49de-a7a6-99ab6d89b0e0
# ╟─ab72fb5f-90ac-446c-88d6-938852f3816e
# ╠═9c851f5a-4c8f-4388-86a1-ff6a3d83c8fe
# ╠═6f6459c8-8242-4340-a588-0512eff5d5c3
# ╠═db565011-1563-41cc-84d0-e14c8237aa8a
# ╠═f518b097-a6f4-4bf5-ba36-0088f5fc6b94
# ╠═98693a39-71de-4e60-bf66-42e30618e254
# ╠═f1ea0441-6910-4704-906a-3e6cf05b4105
# ╠═a0d7c872-8c54-4eed-9a01-1ca92c2d6a28
# ╠═65729744-b221-488d-b1ad-85981489916c
# ╟─c0b3f6a7-341c-435a-b528-3a190f5c500c
# ╠═185015b8-d252-4beb-96f0-4d2483f79344
# ╠═e9fcbd2a-6c9b-4e07-bcae-b990ebc1a3da
# ╟─042c51dc-80d8-4c19-9c87-4e3e83de5e19
# ╠═44899540-8148-4b83-a88a-9b8b4d73ce82
# ╠═18624127-8d76-4b10-ab15-5e7cd4caea79
# ╠═caa805f8-7ece-429b-b63b-0542dbe159c1
# ╟─4247ed0b-7722-4ea9-a243-a48644bc65ee
# ╟─b0bc0de6-a0ab-49c6-a425-6df330e877d4
# ╠═798ee25d-3767-462d-b7d8-522f20754f78
# ╟─4ec0cfec-f200-4063-9cf0-5979cd369af9
# ╟─03a434ee-d037-4607-b346-94b31e5a0fa5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_path: Curso Julia/3 Ciencia de datos/CienciaDeDatos.jl


text:
n = 30
m = 20
p = 5

@testset "Testing $(MLK.kernel)" begin
    for T in (Float32, Float64)
        x = rand(T,p)
        y = rand(T,p)

        x_alt = rand(T == Float32 ? Float64 : Float32, p)

        for f in kernel_functions
            P = (get(kernel_functions_base, f, SquaredEuclidean))()
            F = convert(f{T}, (f)())

            @test isapprox(MLK.kernel(F, x[1], y[1]), MLK.kappa(F, MLK.base_evaluate(P, x[1], y[1])))
            @test isapprox(MLK.kernel(F, x, y),       MLK.kappa(F, MLK.base_evaluate(P, x, y)))

            z = MLK.kernel(F, x_alt[1], y[1])
            @test typeof(z) == T

            z = MLK.kernel(F, x_alt, y)
            @test typeof(z) == T
        end
    end
end

@testset "Testing $(MLK.kappamatrix!)" begin
    for T in (Float32, Float64)
        X = rand(T,n,m)

        for f in kernel_functions
            F = convert(f{T}, (f)())

            K_tmp = [MLK.kappa(F,X[i]) for i in Base.Cartesian.CartesianIndices(size(X))]
            K_tst = MLK.kappamatrix!(F, copy(X))

            @test isapprox(K_tmp, K_tst)
        end
    end
end

@testset "Testing $(MLK.symmetric_kappamatrix!)" begin
    for T in (Float32, Float64)
        X = LinearAlgebra.copytri!(rand(T,n,n), 'U')

        for f in kernel_functions
            F = convert(f{T}, (f)())

            K_tmp = [MLK.kappa(F,X[i]) for i in Base.Cartesian.CartesianIndices(size(X))]
            K_tst = MLK.symmetric_kappamatrix!(F, copy(X), true)

            @test isapprox(K_tmp, K_tst)
            @test_throws DimensionMismatch MLK.symmetric_kappamatrix!(F, rand(T, p, p+1), true)
        end
    end
end

@testset "Testing $(MLK.kernelmatrix!)" begin
    for T in (Float32, Float64)
        X_set = [rand(T,p) for i = 1:n]
        Y_set = [rand(T,p) for i = 1:m]

        K_tst_nn = Array{T}(undef, n, n)
        K_tst_nm = Array{T}(undef, n, m)

        for layout in (Val(:row), Val(:col))
            X = layout == Val(:row) ? permutedims(hcat(X_set...)) : hcat(X_set...)
            Y = layout == Val(:row) ? permutedims(hcat(Y_set...)) : hcat(Y_set...)

            for f in kernel_functions
                F = convert(f{T}, (f)())

                K = [MLK.kernel(F,x,y) for x in X_set, y in X_set]
                @test isapprox(K, MLK.kernelmatrix!(layout, K_tst_nn, F, X, true))

                K = [MLK.kernel(F,x,y) for x in X_set, y in Y_set]
                @test isapprox(K, MLK.kernelmatrix!(layout, K_tst_nm, F, X, Y))
            end
        end
    end
end

@testset "Testing $(MLK.kernelmatrix)" begin
    for T in (Float32, Float64)
        X_set = [rand(Float32,p) for i = 1:n]
        Y_set = [rand(Float32,p) for i = 1:m]

        K_tst_nn = Array{T}(undef, n, n)
        K_tst_nm = Array{T}(undef, n, m)

        for layout in (Val(:row), Val(:col))
            isrowmajor = layout == Val(:row)
            X = convert(Array{T}, isrowmajor ? permutedims(hcat(X_set...)) : hcat(X_set...))
            Y = convert(Array{T}, isrowmajor ? permutedims(hcat(Y_set...)) : hcat(Y_set...))

            X_alt = convert(Array{T == Float32 ? Float64 : Float32}, X)
            Y_alt = convert(Array{T == Float32 ? Float64 : Float32}, Y)

            for f in kernel_functions
                F = convert(f{T}, (f)())

                K_tmp = MLK.kernelmatrix!(layout, K_tst_nn, F, X, true)

                K_tst = MLK.kernelmatrix(layout, F, X, true)
                @test isapprox(K_tmp, K_tst)
                @test eltype(K_tst) == T

                K_tst = MLK.kernelmatrix(layout, F, X_alt)
                @test isapprox(K_tmp, K_tst)
                @test eltype(K_tst) == T

                if isrowmajor
                    K_tst = MLK.kernelmatrix(F, X_alt)
                    @test isapprox(K_tmp, K_tst)
                    @test eltype(K_tst) == T
                end

                K_tmp = MLK.kernelmatrix!(layout, K_tst_nm, F, X, Y)

                K_tst = MLK.kernelmatrix(layout, F, X, Y)
                @test isapprox(K_tmp, K_tst)
                @test eltype(K_tst) == T

                K_tst = MLK.kernelmatrix(layout, F, X_alt, Y_alt)
                @test isapprox(K_tmp, K_tst)
                @test eltype(K_tst) == T

                if isrowmajor
                    K_tst = MLK.kernelmatrix(F, X_alt, Y_alt)
                    @test isapprox(K_tmp, K_tst)
                    @test eltype(K_tst) == T
                end
            end
        end
    end
end

@testset "Testing $(MLK.centerkernelmatrix!)" begin
    for T in FloatingPointTypes
        X = rand(T, n, p)
        Y = rand(T, m, p)

        K = X*permutedims(X)
        MLK.centerkernelmatrix!(K)

        Xc = X .- Statistics.mean(X, dims = 1)
        Kc = Xc*permutedims(Xc)

        @test isapprox(K, Kc)

        K = X*permutedims(Y)
        MLK.centerkernelmatrix!(K)

        Yc = Y .- Statistics.mean(Y, dims = 1)
        Kc = Xc*permutedims(Yc)

        @test isapprox(K, Kc)
    end
end


text:
colorbuffer(screen) = error("Color buffer retrieval not implemented for $(typeof(screen))")

"""
    scene2image(scene::Scene)

Buffers the `scene` in an image buffer.
"""
function scene2image(scene::Scene)
    scene.updated[] = true
    d = global_gl_screen()
    force_update!()
    yield()
    display(d, scene)
    colorbuffer(d)
end

"""
    save(path::String, scene::Scene)

Saves an image of the `scene` at the specified `path`.
"""
function save(path::String, scene::Scene)
    img = scene2image(scene)
    if img != nothing
        try
            FileIO.save(path, img)
        catch e
            # TODO print error?? But it's super long if its JuliaPlots/Makie.jl#138
            error("Failed to save Image. You likely need to install ImageMagick with `]add ImageMagick`")
        end
    else
        # TODO create a screen
        error("Scene isn't displayed on a screen")
    end
end

# Base.showable(::MIME"text/html", scene::VideoStream) = true
Base.showable(::MIME"image/png", scene::Scene) = true

# Let's not lie: with only the OpenGL backend, we don't really support any svg/html mime

# function show(io::IO, mime::MIME"application/javascript", scene::Scene)
#     #TODO use WebIO
#     print(io, scene2javascript(scene))
# end

# function show(io::IO, mime::MIME"text/html", scene::Scene)
#     print(io, "<img src=\"data:image/png;base64,")
#     b64pipe = Base64EncodePipe(io)
#     img = scene2image(scene)
#     FileIO.save(FileIO.Stream(FileIO.format"PNG", b64pipe), img)
#     print(io, "\">")
# end

# function svg(scene::Scene, path::Union{String, IO})
#     cs = CairoBackend.CairoScreen(scene, path)
#     CairoBackend.draw_all(cs, scene)
# end

# function show(io::IO, m::MIME"image/svg+xml", scene::Scene)
#     if false#AbstractPlotting.is2d(scene)
#         svg(scene, io)
#     else
#         show(io, MIME"text/html"(), scene)
#     end
# end

function show(io::IO, m::MIME"image/png", scene::Scene)
    img = scene2image(scene)
    FileIO.save(FileIO.Stream(FileIO.format"PNG", io), img)
end

export VideoStream, recordframe!, finish, record

# Stepper for generating progressive plot examples
mutable struct Stepper
    scene::Scene
    folder::String
    step::Int
end

function Stepper(scene, path)
    ispath(path) || mkpath(path)
    Stepper(scene, path, 1)
end

"""
    step!(s::Stepper)
steps through a `Makie.Stepper` and outputs a file with filename `filename-step.jpg`.
This is useful for generating progressive plot examples.
"""
function step!(s::Stepper)
    Makie.save(joinpath(s.folder, basename(s.folder) * "-$(s.step).jpg"), s.scene)
    s.step += 1
    return s
end

export Stepper, step!


text:
using StaticArrays
@testset "Correlated Wiener Process" begin

using DiffEqNoiseProcess, DiffEqBase, Test, Statistics

# oop
ρ = 0.3
Γ = [1 ρ;ρ 1]
W = CorrelatedWienerProcess(Γ, 0., zeros(2), zeros(2))

dt = 0.1
calculate_step!(W,dt,nothing,nothing)
for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

# inplace
W = CorrelatedWienerProcess!(Γ, 0., zeros(2), zeros(2))

dt = 0.1
calculate_step!(W,dt,nothing,nothing)
for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.01)

output_func = (sol,i) -> (sol.dW,false)
ensemble_prob = EnsembleProblem(prob,output_func=output_func)
@time sol = Array(solve(ensemble_prob,dt=dt,trajectories=1_000_000))

@test zero(mean(sol,dims=2)[:]) ≈ mean(sol,dims=2)[:] atol=1e-2
@test Γ ≈ cov(sol,dims=2)/dt rtol=1e-2

# with StaticArrays
Γ = @SMatrix [1.0 ρ
              ρ 1.0]
W = CorrelatedWienerProcess(Γ, 0., @SVector(zeros(2)), @SVector(zeros(2)))

dt = 0.1
calculate_step!(W,dt,nothing,nothing)
for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

@test typeof(sol[1]) <: SArray
@test typeof(sol[end]) <: SArray

end


text:
module ExtremeLearningMachines

using JLD
using LinearAlgebra
using Quaternions

export ELM
export Regression, Classification
export train!, predict
export Sigmoid, HardLimit, ReLU, Cosine, Qubit, HyperbolicTangent
export Batch, BatchCM, BatchL2, BatchCML2
export Uniform
export save, load

include("common.jl")
include("mapping.jl")
include("initializer.jl")
include("training.jl")
include("core.jl")
include("io.jl")

end # module


text:
using  Languages, TextAnalysis, Flux, PyPlot, Statistics, MLDataUtils, Embeddings, MLLabelUtils
using Flux: logitcrossentropy, normalise, onecold, onehotbatch
using Parameters: @with_kw
using BSON: @save,@load
@with_kw mutable struct Args
    lr::Float64 = 0.01 # learning rate
    τ::Int = 110 # max iterations
    path::String = "Rango/data/sample.txt" # path to data
    max_length::Int = 4 # Max length of the doc to be considered
    max_features::Int = 8 #Number of features per word to learn for training
end
# Padding the corpus wrt the longest document
function pad_corpus(c, pad_size, word_dict)
    M = []
    for doc in 1:length(c)
        tks = tokens(c[doc])
        if length(tks) >= pad_size
            tk_indexes = [tk_idx(w,word_dict) for w in tks[1:pad_size]]
        end
        if length(tks) < pad_size
            tk_indexes = zeros(Int64,pad_size - length(tks))
            tk_indexes = vcat(tk_indexes, [tk_idx(w,word_dict) for w in tks])
        end
        doc == 1 ? M = tk_indexes' : M = vcat(M, tk_indexes')
    end
    return M
end
# function to return the index of the word in the word dictionary
tk_idx(s,word_dict) = haskey(word_dict, s) ? i = word_dict[s] : i = 0
function get_processed_data(args)
    f = open(args.path)
    doc_array = readlines(f)
    labels, texts = [], []
    for doc in doc_array
        content = split(doc)
        push!(labels,content[1])
        push!(texts,join(content[2:end]," "))
    end
    # pushing the text from the files to the string documents
    docs = []
    for i in 1:length(texts)
        push!(docs, StringDocument(texts[i]))
    end
    # building a Corpus
    corpus = Corpus(docs)
    # updating the lexicon and creating the word dict
    update_lexicon!(corpus)
    doc_term_matrix = DocumentTermMatrix(corpus)
    word_dict = doc_term_matrix.column_indices
    # splitting words in the document
    word_docs = map(s -> split(s,r"[,. ]",keepempty=false),texts)
    # pad size is the number of words in the maximum word document
    # Can set a fixed length or the max doc length
    pad_size = args.max_length
    # padding the docs
    padded_docs = pad_corpus(corpus, pad_size,word_dict)
    # forming the data with the labels
    X = padded_docs'
    (X_train, y_train), (X_test, y_test) = splitobs((X, labels); at = 0.67)
    klasses = sort(unique(labels))
    y_train = convert(Array{Bool}, y_train .== "True")
    y_test = convert(Array{Bool}, y_test .== "True")
    y_train = reshape(y_train,1,size(y_train,1))
    y_test = reshape(y_test,1,size(y_test,1))
    train_data = [(X_train, y_train)]
    test_data = [(X_test, y_test)]
    return train_data, test_data, doc_term_matrix, klasses
end
# Function to build confusion matrix
function confusion_matrix(X, y, model)
    ŷ = onehotbatch(onecold(model(X)), 1:2)
    y * transpose(ŷ)
end
function getArgs(; kws...)
    # Initialize hyperparameter arguments
    global args = Args(; kws...)
end
function test(model, test)
    # Testing model performance on test data
    X_test, y_test = test[1]
    accuracy_score = accuracy(X_test, y_test, model)
    println("\nAccuracy: $accuracy_score")
    # Sanity check. TO be done in testing
    # @assert accuracy_score > 0.6
    println("\nConfusion Matrix:\n")
    display(confusion_matrix(X_test, y_test, model))
end
getArgs()
#Loading processed data
train_data, test_data, doc_term_matrix, klasses = get_processed_data(args)
word_dict = doc_term_matrix.column_indices
X_train, y_train = train_data[1]
max_features = args.max_features
pad_size = args.max_length
# Building Flux Embeddings
N = size(X_train,2)  #Number of documents
# features per word to learn, depends on the size of the corpus, larger corpus will probably need a higher dimension
# number of words in the vocabulary, should always be higher than the maximum index in our dictionary.
vocab_size = maximum(word_dict)[2] + 1
# Embedding layer for Flux model
# glorot_normal returns an Array of size dims containing random variables taken from a normal distribution with mean 0 and standard deviation (2 / sum(dims)).
struct EmbeddingLayer
   W
   EmbeddingLayer(mf, vs) = new(Flux.glorot_normal(mf, vs))
end
# Enabling Flux to learn
@Flux.treelike EmbeddingLayer
(m::EmbeddingLayer)(x) = m.W * Flux.onehotbatch(reshape(x, pad_size*N), 0:vocab_size-1)
m = Chain(EmbeddingLayer(max_features, vocab_size),
          x -> reshape(x, max_features, pad_size, N),
          x -> sum(x, dims=2),
          x -> reshape(x, max_features, N),
          Dense(max_features, 1, σ)
)
# Defining loss function to be used in training
# Training
# Gradient descent optimiser with learning rate `args.lr`
loss_h = []
accuracy_train = []
accuracy(x, y) = mean(x .== y)
loss(x, y) = sum(Flux.binarycrossentropy.(m(x), y))
optimizer = Flux.Descent(args.lr)
println("Starting training.")
for epoch in 1:args.τ
    Flux.train!(loss, Flux.params(m), train_data, optimizer)
    println(loss(X_train, y_train), " ", accuracy(m(X_train).>0.5,y_train))
    push!(loss_h, loss(X_train, y_train))
    push!(accuracy_train, accuracy(m(X_train).>0.5,y_train))
end

max_stars_repo_name: akaysh/Rango

max_issues_repo_name: akaysh/Rango

max_forks_repo_name: akaysh/Rango


text:
using GNSSTools
using JLD
using PyPlot
pygui(true)

function vsm(trackresults, K)
    ZP = trackresults.ZP
    M = floor(Int, length(ZP)/K)
    CN0_est = Array{Float64}(undef, M)
    t_est = Array{Float64}(undef, M)
    for i in 1:M
        t_est[i] = trackresults.t[i*K]
        zp = ZP[(i-1)*K+1:i*K]
        Zₖ = real.(zp).^2 + imag.(zp).^2
        Z̄ = mean(Zₖ)
        σ² = var(Zₖ)
        A = sqrt(Z̄^2 - σ²)
        σ²_iq = (Z̄ - A)/2
        CN0 = 10log10(A/(2*σ²_iq*trackresults.T))
        CN0_est[i] = CN0
    end
    return (CN0_est, t_est)
end

function snr_est(trackresults, K)
    ZP = trackresults.ZP
    M = floor(Int, length(ZP)/K)
    SNR_est = Array{Float64}(undef, M)
    t_est = Array{Float64}(undef, M)
    for i in 1:M
        t_est[i] = trackresults.t[i*K]
        zp = ZP[(i-1)*K+1:i*K]
        I² = mean(real.(zp))^2
        Q² = mean(imag.(zp))^2
        SNR_est[i] = I² / (2*Q²)
    end
    return (SNR_est, t_est)
end

h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # Rakon IT5300B TCXO

#------------------- Lab Data Sample for L1 C/A Signal ------------------------#
# file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/known_oscillator_data/20210327_013800_000005_Aero-North-TCXO_1575.42_25.sc16"
# file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/20201217_224100_GPSL1_1575.42_25.0.sc8"
file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/20201211_211745_GPSL1_1575.42_25.0.sc4"
data = loaddata(:sc4, file, 25e6, L1_freq, L1_freq, 3, skip_to=1)

replica = definereplica(define_l1ca_code_type(), 25e6, 2e-3)
detectable_prns = []
for i in 1:32
    fd_est, n0_est, SNR_est,
    P_d, above_threshold = courseacquisition(data, replica, i; M=20,
                                             Δfd=1/(1*replica.t_length))
    if (P_d >= 0.9)
        detectable = true
        push!(detectable_prns, i)
    else
        detectable = false
    end
    println("PRN $i\t$SNR_est\t$(Float64(round(P_d, digits=2)))\t$detectable")
end
println("$(length(detectable_prns)) Detected PRNs: $detectable_prns")

fig = figure()
acquisition_T = 2e-3
fine_acq_T = 10e-3
tracking_T = 1e-3
M = 10
for i in 1:32
    ax = fig.add_subplot(8, 4, i)
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $i")
    acqresults, trackresults = process(data, define_l1ca_code_type(), i;
                                       show_plot=false,
                                       fine_acq_method=:fft, state_num=2,
                                       dynamickf=true,
                                       # h₀=h_parms[3],
                                       # h₋₂=h_parms[1],
                                       dll_b=0.01,
                                       acquisition_T=acquisition_T,
                                       fine_acq_T=fine_acq_T,
                                       tracking_T=tracking_T,
                                       M=M)
    ax.plot(trackresults.t, real.(trackresults.ZP))
    ax.plot(trackresults.t, imag.(trackresults.ZP))
end
subplots_adjust(hspace=0.4, wspace=0.4,  top=0.93, left=0.08, right=0.93,
                bottom=0.01)
plot_title = string("Course Acq. T=$(floor(Int, acquisition_T*1000))ms;    ",
                    "Fine Acq. T=$(floor(Int, fine_acq_T*1000))ms;    ",
                    "Tracking T=$(floor(Int, tracking_T*1000))ms;    ",
                    "Noncoherent Integrations=$(M)")
suptitle(plot_title)

acqresults, trackresults, corr_result, SNR_est,
P_d, above_threshold = process(data,
                               define_l1ca_code_type(),
                               3;
                               fine_acq_method=:fft,
                               return_corrresult=true,
                               return_Pd=true,
                               state_num=3,
                               # h₀=h_parms[3],
                               # h₋₂=h_parms[1],
                               q_a=1,
                               dll_b=1,
                               acquisition_T=2e-3,
                               fine_acq_T=10e-3,
                               tracking_T=1e-3,
                               M=10,
                               use_fine_acq=true,
                               show_plot=false,
                               σω=1);

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

#--------------------- Lab Data Sample for L5 Signal --------------------------#
# file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/20201211_211300_GPSL5_1176.45_25.0.sc4"
file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/20200622_203230_000100_g5_ref.jld"
# data = loaddata(:sc4, file, 25e6, L5_freq, L5_freq, 1, skip_to=0)
data = load(file, "data")

sigtype = definesignaltype(definecodetype([l5q_codes, nh20],
                                          [L5_chipping_rate, nh_chipping_rate]),
                                          L5_freq, "Q")
replica = definereplica(sigtype, 25e6, 20e-3)
detectable_prns = []
for i in 1:32
    fd_est, n0_est, SNR_est,
    P_d, above_threshold = courseacquisition(data, replica, i; M=1,
                                             Δfd=1/(1*replica.t_length))
    if above_threshold && (P_d >= 0.9)
        push!(detectable_prns, i)
    end
    println("PRN $i\t$SNR_est\t$(Float64(round(P_d, digits=2)))\t$above_threshold")
end
println("Detected PRNs: $detectable_prns")

fig = figure()
acquisition_T = 1e-3
fine_acq_T = 100e-3
tracking_T = 1e-3
M = 20
for i in 1:32
    ax = fig.add_subplot(8, 4, i)
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $i")
    acqresults, trackresults = process(data, sigtype, i, "Q";
                                       show_plot=false,
                                       fine_acq_method=:fft, state_num=2,
                                       dynamickf=true, h₀=h_parms[3],
                                       h₋₂=h_parms[1], dll_b=0.01,
                                       acquisition_T=acquisition_T,
                                       fine_acq_T=fine_acq_T,
                                       tracking_T=tracking_T,
                                       M=M)
    ax.plot(trackresults.t, real.(trackresults.ZP))
    ax.plot(trackresults.t, imag.(trackresults.ZP))
end
subplots_adjust(hspace=0.4, wspace=0.4,  top=0.93, left=0.08, right=0.93,
                bottom=0.01)
plot_title = string("Course Acq. T=$(floor(Int, acquisition_T*1000))ms;    ",
                    "Fine Acq. T=$(floor(Int, fine_acq_T*1000))ms;    ",
                    "Tracking T=$(floor(Int, tracking_T*1000))ms;    ",
                    "Noncoherent Integrations=$(M)")
suptitle(plot_title)

acqresults, trackresults, corr_result, SNR_est,
P_d, above_threshold = process(data,
                               sigtype,
                               6, "Q";
                               fine_acq_method=:fft,
                               return_corrresult=true,
                               return_Pd=true,
                               state_num=2,
                               dynamickf=true,
                            #    h₀=h_parms[3],
                            #    h₋₂=h_parms[1],
                               q_a=1,
                               dll_b=0.1,
                               acquisition_T=1e-3,
                               fine_acq_T=100e-3,
                               tracking_T=20e-3,
                               M=20,
                               show_plot=true,
                               use_fine_acq=true);

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


text:
function build_metabolite_keggid_matching_table()::Some

    try

        # initialize -
        metabolite_matching_table = Dict{String,String}()

        # build the url -
        url_string = _KEGG_LIST_COMPOUND_URL

        # make the call -
        http_body = http_get_call_with_url(url_string) |> check

        # ok, so no body?
        if (length(http_body) == 0)
            return Some(nothing)
        end

        # populate the table -
        records_array = split(http_body, "\n")
        for record in records_array

            # records take the form:
            # kegg_compound_id\t{Names;*}

            if (record != "")
                # split around the \t -
                tmp_compound_name = string(split(record, "\t")[1])
                kegg_compound_id = string(split(tmp_compound_name, ":")[2])
                compound_name_array = split(record, "\t")[2]
                if (occursin(";", compound_name_array) == true)

                    # ok, so we have multiple possible names -
                    synomym_array = split(compound_name_array, ";")
                    for synonum in synomym_array
                        metabolite_matching_table[lowercase(lstrip(string(synonum)))] = kegg_compound_id
                    end
                else
                    metabolite_matching_table[lowercase(lstrip(string(compound_name_array)))] = kegg_compound_id
                end
            end
        end

        # return -
        return Some(metabolite_matching_table)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_compound_records(compound_array::Array{String,1})::Some

    try

        # initialize -
        kegg_object_array = Array{KEGGCompound,1}()

        # main loop -
        for compound in compound_array

            # build compound object -
            object = get_compound_records(compound) |> check
            if (isnothing(object) == false)
                # push -
                push!(kegg_object_array, object)
            end
        end

        # return -
        return Some(kegg_object_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_compound_records(compound::String)::Some

    try

        # do we have a bad compound code?
        if (contains(compound, " ") == true || compound[1] != 'C' ||
            contains(compound, "(") == true)
            return Some(nothing)
        end

        # build url -
        url_string = "$(_KEGG_GET_URL)/$(compound)"

        # make the call -
        http_body = http_get_call_with_url(url_string) |> check
        if (isempty(http_body) == true)
            return Some(nothing)
        end

        # initailize -
        compound_object = KEGGCompound()

        # ok, so we need to split this around the \n
        record_components = string.(split(http_body, "\n"))
        compound_object.compound_id = compound

        # NAME is always the [2] -
        tmp_name = split(record_components[2], repeat(" ", 5))[2] |> lstrip |> rstrip |> lowercase
        tmp_compound_name = (last(tmp_name) == ';') ? tmp_name[1:end-1] : tmp_name
        compound_object.compound_name = replace(tmp_compound_name, "+" => "")

        # FORMULA section -
        formula_section = extract_db_file_section(record_components, "FORMULA")
        if (isnothing(formula_section) == true)
            return Some(nothing) # skip that compound -
        else
            tmp_compound_formula = string(split(formula_section, repeat(" ", 5))[2] |> lstrip |> rstrip)
            compound_object.compound_formula = tmp_compound_formula
            compound_object.atom_dictionary = extract_atom_dictionary(tmp_compound_formula)
        end

        # MOL_WEIGHT -> do we have a molecular weight?
        mw_section = extract_db_file_section(record_components, "MOL_WEIGHT")
        if (isnothing(mw_section) == true)
            compound_object.compound_mw = 0.0
        else
            compound_object.compound_mw = parse(Float64, (split(mw_section, repeat(" ", 2))[2] |> lstrip |> rstrip))
        end

        # return -
        return Some(compound_object)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_compound_records_for_reaction(reaction_code::String)::Some

    try

        # get the reaction object with metabolites in KEGG markup (not english names) -
        kegg_reaction_object = get_reaction_for_rn_number(reaction_code) |> check
        if (isnothing(kegg_reaction_object) == true)
            return Some(nothing)
        end

        # initialize -
        reaction_compound_record_array = Array{KEGGCompound,1}()

        # pool the forward and reverse reaction strings -
        tmp_reaction_phrase_pool = Array{String,1}()
        push!(tmp_reaction_phrase_pool, kegg_reaction_object.reaction_forward)
        push!(tmp_reaction_phrase_pool, kegg_reaction_object.reaction_reverse)

        # build a local metabolite code pool -
        local_metabolite_code_pool = Array{String,1}()
        for reaction_phrase in tmp_reaction_phrase_pool

            # split around the + and strip space -
            compound_code_array = split(reaction_phrase, "+")
            for raw_compound_code in compound_code_array

                # turn code into string, strip spaces -
                tmp_code = string(raw_compound_code) |> lstrip |> rstrip

                # cache -
                value = "cpd:$(tmp_code)"
                push!(local_metabolite_code_pool, value)
            end
        end

        for metabolite_code in local_metabolite_code_pool
            compound_record = get_compound_records(metabolite_code) |> check
            if (isnothing(compound_record) == true)
                return Some(nothing)
            end
            push!(reaction_compound_record_array, compound_record)
        end

        return Some(reaction_compound_record_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end


text:
"""
BinaryVectorProbability

Includes types and methods for representing probability distributions on the set of binary vectors.
"""
module BinaryVectorProbability

using JLD
using Optim

import Base: show, ==, isequal, hash
include("metadata.jl")

################################################################################
#### Abstract type to use as placeholder/to get all docstrings in one place
################################################################################
"""
    AbstractBinaryVectorDistribution

Abstract type used as placeholder.
"""
abstract type AbstractBinaryVectorDistribution end

################################################################################
#### Common functions
################################################################################

_NI(m) = error("Not implemented: $m")

# override this definition for specific types of distributions. This should represent
# equality in the sense of P(x) = Q(x) for all x; other attributes shouldn't matter. on the
# other hand, isequal is a more stringent check used for hashing.
==(P::AbstractBinaryVectorDistribution, Q::AbstractBinaryVectorDistribution) = false

"""
    n_bits(P::AbstractBinaryVectorDistribution)

Size of the sample space
"""
n_bits(P::AbstractBinaryVectorDistribution) = _NI("n_bits($(typeof(P)))")

"""
    _binary_to_int(x::Union{BitVector, Vector{Bool}})

Returns the integer that `x` represents. Will use `x.chunks` for `BitVector`s
and the "dot product method" for `Vector{Bool}`s. Note that this will really
mess you up if `length(x) > 63`

"""
_binary_to_int(x::BitVector) = Int(x.chunks[1])
_binary_to_int(x::AbstractVector{Bool}) = dot([2^i for i = 0:(length(x) - 1)], x)

"""
    pdf(P)

Returns the function `x -> pdf(P, x)`
"""
pdf(P::AbstractBinaryVectorDistribution) = x -> pdf(P, x)

"""
    get_pdf(P)

Returns the pdf of `P` as a `Vector{Float64}`, where `pdf[i]` is the probability
of `digits(Bool, i-1, 2, n_bits(P))`

"""
function get_pdf(P::AbstractBinaryVectorDistribution)
    if get(P.metadata, :pdf_computed, false)
        return full(P.cache[:pdf])
    else
        P.metadata[:pdf_computed] = true
        return [pdf(P, digits(Bool,x,2,n_bits(P))) for x in 0:(2^n_bits(P) - 1)]
    end
end

"""
    get_cdf(P)

Returns the "cdf" of `P` as a `Vector{Float64}`, meaning `cumsum(get_pdf(P))`.
(I mean, it's a "cdf" in the sense that you can put a total order {0,1}^n, but
really it's just used for sampling from the distribution.)

"""
function get_cdf(P::AbstractBinaryVectorDistribution)
    if haskey(P.cache, :cdf)
        cdf = P.cache[:cdf]
    else
        cdf = cumsum_kbn(get_pdf(P))
        P.cache[:cdf] = cdf
    end
    return cdf
end

"""
    random(P, n_sample)

Naive sampling method to draw random vectors from a distribution. Relies on an
implementation of `get_cdf` for the distribution (so, don't use with distributions on more
than ~20 neurons). Numerical accuracy not guaranteed. Returns a `BitMatrix`.

"""
function random(P::AbstractBinaryVectorDistribution, n_samples=1)
    X = falses(n_bits(P), n_samples)
    r = rand(n_samples)
    cdf = get_cdf(P)
    for s = 1:n_samples
        k = searchsortedfirst(cdf, r[s])
        X[:,s] = digits(Bool, k - 1, 2, n_bits(P))
    end
    return X
end

"""
    expectation_matrix(P)

Returns a matrix with the expected values of bits and pairs of bits. Default
behavior is to call `get_pdf(P)` and then loop through those values.

"""
function expectation_matrix(P::AbstractBinaryVectorDistribution)
    em = zeros(n_bits(P), n_bits(P))
    p = get_pdf(P)
    for k in sortperm(p)
        x = digits(Bool, k-1, 2, n_bits(P))
        em += p[k] * (x * x')
    end
    return em
end

"""
    entropy(P)

Entropy, in nats (i.e. natural log is used) of the given distribution. Sets a
metadata value.

"""
function entropy(P::AbstractBinaryVectorDistribution)
    if !haskey(P.metadata, :entropy)
        #MAYBEDO figure out how the get! function works. I assume since it's a
        #function call, it'll evalute all arguments first. But, if that's not
        #the case this could be cleaned up, I suppose.
        P.metadata[:entropy] = -sum_kbn([p * log(p) for p in get_pdf(P) if 0.0 < p < 1.0])
    end
    return P.metadata[:entropy]
end

"""
    entropy2(P)

Entropy, in bits (i.e. log base 2 is used) of the given distribution. Sets a
metadata value.

"""
function entropy2(P::AbstractBinaryVectorDistribution)
    if !haskey(P.metadata, :entropy2)
        P.metadata[:entropy2] = -sum_kbn([p * log2(p) for p in get_pdf(P) if 0.0 < p < 1.0])
    end
    return P.metadata[:entropy2]
end

################################################################################
#### Constants
################################################################################
"""
    ISING_METHOD_THRESHOLD

If N_neurons is greater than this number, switch from exact sampling to Gibbs
sampling and use MPF instead of LogLikelihood to for the parameters.
"""
const ISING_METHOD_THRESHOLD = 20
const METH_THRESH_WARN_KEY = 1

################################################################################
#### Include definitions of concrete distributions, other functions
################################################################################
include("DataDistribution.jl")
include("BernoulliCodeDistribution.jl")
include("IsingDistribution.jl")

include("fileio.jl")
include("objective_functions.jl")
include("optimizers.jl")
include("fitters.jl")

export AbstractBinaryVectorDistribution,
       show, ==, isequal, hash,
       n_bits, random, expectation_matrix, entropy, entropy2, pdf, get_pdf, get_cdf,
       savedistribution, loaddistribution,
       DataDistribution, BernoulliCodeDistribution, IsingDistribution,
       first_order_model, second_order_model, data_model,
       loglikelihood, MPF_objective,
       negloglikelihood, dnegloglikelihood!, K_MPF, dK_MPF!

end

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbcirc(radius::F=1.0,center::Vector{F}=[0.0;0.0;0.0],
            sang::F=0.0,eang::F=2.0*π)::NURBS1D where {F<:AbstractFloat}

Constructs NURBS data structure for a circular arc in the x-y plane. If no
arguments are supplied a unit circle with center `[0.0;0.0]` is constructed.

Angles are defined as positive in the anti-clockwise direction.

# Arguments:
- `radius`: Radius of the circle, default `1.0`
- `center`: Center of the circle, default `[0.0;0.0;0.0]`
- `sang`: Start angle, default `0` degrees
- `eang`: End angle `360` degrees

# Output:
- `crv`: NURBS curve for a circular arc.

# Examples:
```julia
julia> crv = nrbcirc()
julia> crv = nrbcirc(radius)
julia> crv = nrbcirc(radius,center)
julia> crv = nrbcirc(radius,center,sang,eang)
```
"""
function nrbcirc(radius::F=1.0,center::Vector{F}=[0.0;0.0;0.0],
                 sang::F=0.0,eang::F=2.0*π)::NURBS1D where {F<:AbstractFloat}
# TODO: maybe revisit this function and split it up
# TODO: optional arguments in julia are always at the end and providing [] as
#       argument input is a bad habbit i changed the argument behaviour for
#       center > examine if that did change the behaviour in this toolbox

# sweep angle of arc
sweep = eang - sang;
if sweep < 0
  sweep = 2 * π + sweep;
end # if

if abs(sweep) <= π/2
  narcs = 1;
  knots = vec([0.0 0.0 0.0 1.0 1.0 1.0]);
elseif abs(sweep) <= π
  narcs = 2;
  knots = vec([0.0 0.0 0.0 1/2 1/2 1.0 1.0 1.0]);
elseif abs(sweep) <= π*3/2
  narcs = 3;
  knots = vec([0.0 0.0 0.0 1/3 1/3 2/3 2/3 1.0 1.0 1.0]);
else
  narcs = 4;
  knots = vec([0.0 0.0 0.0 1/4 1/4 1/2 1/2 3/4 3/4 1.0 1.0 1.0]);
end # if

# arc segment sweep angle/2
dsweep = sweep / (2 * narcs);

# determine middle control point and weight
wm = cos(dsweep);
x  = radius * wm;
y  = radius * sin(dsweep);
xm = x + y * tan(dsweep);

# arc segment control points
ctrlpt = [ x wm*xm x;  # w*x - coordinate
          -y   0   y;  # w*y - coordinate
           0   0   0;  # w*z - coordinate
           1   wm  1]; # w   - coordinate

# build up complete arc from rotated segments
coefs = zeros(F,(4,2*narcs+1)); # nurbs control points of arc
xx = vecrotz(sang + dsweep);
coefs[:,1:3] = xx * ctrlpt;     # rotate to start angle
xx = vecrotz(2 * dsweep);
for n in 2:narcs
  m = (0:1) .+ (2*n);
  coefs[:,m] = xx * coefs[:,m .- 2];
end # if

# vectrans arc if necessary
if !(all(i->i==0,center))
  xx = vectrans(center);
  coefs = xx * coefs;
end # if

return nrbmak(coefs,knots);
end # nrbcirc


text:
using DataFrames, CSV, Statistics, StatsBase
data = CSV.read("../externals/Core.Math.Data/data/Pejcic_318.csv", copycols = true)

# println("Sample Mean: ", mean(data[:ATT]))
# println("Harmonic <= Geometric <= Arithmetic ",(harmmean(data), geomean(data), mean(data)))
# println("Sample Variance: ",var(data))
# println("Sample Standard Deviation: ",std(data))
# println("Minimum: ", minimum(data))
# println("Maximum: ", maximum(data))
# println("Median: ", median(data))
# println("95th percentile: ", percentile(data, 95))
# println("0.95 quantile: ", quantile(data, 0.95))
# println("Interquartile range: ", iqr(data),"\n")

for col in eachcol(data,true);
    if (eltype(col[2]) <: Real) println(mean(col[2]));
    end;
end;


text:
import Base: one, size, zero

struct SizedType{S,T}
    sz::S
    SizedType(::Type{T}, sz::S) where {S,T} = new{S,T}(sz)
end

"""
    sizedtype(x) -> T

Return an extended type which preserves `size` information. Makes `one(T)` and `zero(T)`
for vectors possible.
"""
sizedtype(x::T) where {T} =  SizedType{T}(())
sizedtype(x::Array{T}) where {T} = SizedType(T, size(x))
size(st::SizedType) = st.sz
basetype(st::SizedType{S,T}) where {S,T} = T

one(st::SizedType) = ones(basetype(st), size(st))
one(st::SizedType{Tuple{}}) = one(basetype(st))

zero(st::SizedType) = zeros(basetype(st), size(st))
zero(st::SizedType{Tuple{}}) = zero(basetype(st))


text:
module EnhancedLogging

using Printf, Logging, Distributed
import Logging: handle_message, shouldlog, min_enabled_level, catch_exceptions, LogLevel

export ProgressLevel, EnhancedConsoleLogger, WorkerLogger

const ProgressLevel = LogLevel(-1)

include("EnhancedConsoleLogger.jl")
include("WorkerLogger.jl")

end # module


text:
function test_func(x)
    function max_missing(x)
        """Finds the maximum value of a array where missing are present."""
            y = replace(copy(x), missing => 0)
            return(maximum(y))
        end
    max = max_missing(x)
    return x.+ max
end

@testset "cross_apply" begin
    for i in 1:10
        x = rand(10:30)
        y = rand(10:30)
        z = rand(10:30)
        rad = rand(20:100.0, x, y, z)
        @test size(cross_apply(test_func, rad)) == size(rad)
    end
end

@testset "long_apply" begin
    for i in 1:10
        x = rand(10:30)
        y = rand(10:30)
        z = rand(10:30)
        rad = rand(20:100.0, x, y, z)
        mask = rand(0:1, x, y)
        @test size(long_apply(test_func, rad, mask)) == size(rad)
    end
end


text:
# ---
# title: 1486. XOR Operation in an Array
# id: problem1486
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Bit Manipulation
# link: <https://leetcode.com/problems/xor-operation-in-an-array/description/>
# hidden: true
# ---
#
# Given an integer `n` and an integer `start`.
#
# Define an array `nums` where `nums[i] = start + 2*i` (0-indexed) and `n ==
# nums.length`.
#
# Return the bitwise XOR of all elements of `nums`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 5, start = 0
#     Output: 8
#     Explanation: Array nums is equal to [0, 2, 4, 6, 8] where (0 ^ 2 ^ 4 ^ 6 ^ 8) = 8.
#     Where "^" corresponds to bitwise XOR operator.
#
#
# **Example 2:**
#
#
#
#     Input: n = 4, start = 3
#     Output: 8
#     Explanation: Array nums is equal to [3, 5, 7, 9] where (3 ^ 5 ^ 7 ^ 9) = 8.
#
# **Example 3:**
#
#
#
#     Input: n = 1, start = 7
#     Output: 7
#
#
# **Example 4:**
#
#
#
#     Input: n = 10, start = 5
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 1000`
#   * `0 <= start <= 1000`
#   * `n == nums.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Cambrian
using NSGAII

function pop3D()
    liste=Array{Array{Int64,1}}(undef,0)
    n=1
    for i in 1:4
        push!(liste,[n,n,n])
        push!(liste,[n+1,n,n])
        push!(liste,[n,n+1,n])
        push!(liste,[n,n,n+1])
        n+=1
    end
    r=reverse([1,1,1,2,3,3,3,4,5,5,5,6,7,7,7,8])
    return liste,r
end

function pop1D()
    liste=[[i] for i in 1:20]
    ranks=reverse([i for i in 1:20])
    return liste,ranks
end

function subPop1D()
    liste=[[1],[1],[1]]
    return liste
end

function subPop3D()
    liste=[[5,3,4],[3,5,4],[2,4,5],[4,4,4],[4,4,4]]
    return liste
end

function subPopEqual3D()
    liste=[[5,5,5] for i in 1:8]
    return liste
end

function maxPop(pop)
    n=length(pop[1].fitness)
    max=[]
    for i in 1:n
        sort!(pop,by=x->x.fitness[i])
        push!(max,pop[end])
    end
    return max
end

function minPop(pop)
    n=length(pop[1].fitness)
    min=[]
    for i in 1:n
        sort!(pop,by=x->x.fitness[i])
        push!(min,pop[1])
    end
    return min
end


text:
module TreeShaker

# Deps
using SnoopCompile # for snooping
using Random # for snoop log names
using Pkg # context/manifest info
using Suppressor # for optionally suppressing build & test output

# Includes
include("trunk.jl")

# Exports
export shake

end # module


text:
# If you want to write expressions of the form
#     @test_formany x::MyType, f(x)
# without going through the steps of writing
#     @test_formany 1 < x < 10, 1 < y < 10, f(MyType(x,y))...
# You can define a custom generator. An example is below.

using Checkers

type ExampleType
	field1::Float64
end

## Fully qualify 'custom_generator' if not explicitly importing from Checkers
#TODO: Let's explain more how this works and also show a more complicated
#example that doesn't just wrap a single number.
function Checkers.custom_generator(a::Number, b::Number, type_a::Bool, type_b::Bool, mytype)
    (size) -> ExampleType(a+(b-a)*rand(Float64))
end

@testset "Testing custom generator for ExampleType" begin
    #This should pass
	@test_formany 1 < x::ExampleType < 10, 1 < x.field1 < 10
    #This should fail because we did not design the generator
    #to choose only 1 as the value of field1.
	@test_formany 1 < x::ExampleType < 10, x.field1 == 1
end

# elementary example with complex types
type multifield
	field1::Float64
	field2::Float64
	field3::Float64
end

function Checkers.custom_generator{T<:multifield}(var_type::Type{T}, p::Number)
	return (size) -> multifield(p*rand()/size, p*rand(), p*rand()*size)
end

@test_exists mf::multifield, 10, mf.field1<.1 && 1< mf.field2 <10 && mf.field3>100


text:
# Afft.jl

using MIRT: Afft

using LinearMapsAA: LinearMapAM, LinearMapAO
using Test: @test

samp = trues(3,2); samp[2] = false

@testset "AM" begin
	A = Afft(samp ; operator=false)
	@test A isa LinearMapAM
	@test Matrix(A)' ≈ Matrix(A')
	@test A * vec(ones(3,2)) == [6, 0, 0, 0, 0]
	@test A' * [1, 0, 0, 0, 0] == vec(ones(3,2))
	@test A.name == "fft"
	@test eltype(A) == ComplexF32
end

@testset "AO" begin
	A = Afft(samp)
	@test A isa LinearMapAO
	@test A * ones(3,2) == [6, 0, 0, 0, 0]
	@test A' * [1, 0, 0, 0, 0] == ones(3,2)
	@test Matrix(A)' ≈ Matrix(A')
end


text:
## test_gaussian_random_fields.jl : test implementation of GRF constructors

@testset "Gaussian random fields   " begin

    # vectors instead of linspaces
    cov = CovarianceFunction(2,Exponential(0.1))
    pts1 = collect(0:0.05:1)
    pts2 = collect(0:0.01:1)
    grf = GaussianRandomField(cov,KarhunenLoeve(500),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Exponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{KarhunenLoeve{500}})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    # test non-equidistant grid
    cov = CovarianceFunction(2,Matern(0.3,1.))
    pts1 = 0.5 * (cos.((0:0.05:1)*pi) .+ 1)
    pts2 = 0.5 * (cos.((0:0.05:1)*pi) .+ 1)
    grf = GaussianRandomField(cov,Cholesky(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    # test wrong dimension of points
    cov = CovarianceFunction(2,Exponential(0.3))
    pts = 0:0.1:1
    @test_throws MethodError GaussianRandomField(cov, Cholesky(), pts)

    # test sample with random number vector (and wrong length)
    cov = CovarianceFunction(1,Matern(1.0,2.,σ=1.))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,Cholesky(),pts)
    @test length(sample(grf,xi=rand(length(pts)))) == length(pts)
    @test_throws DimensionMismatch sample(grf,xi=rand(length(pts)+1))

    # test supplying a mean value for the field
    cov = CovarianceFunction(2,Exponential(1))
    pts1 = 0:0.05:1
    pts2 = 0:0.05:1
    grf = GaussianRandomField(1,cov,Cholesky(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Exponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)
    @test size(grf.mean,1) == length(pts1)
    @test size(grf.mean,2) == length(pts2)
    @test all(grf.mean.==1)

    grf = GaussianRandomField(2.0*ones(length(pts1),length(pts2)),cov,Cholesky(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Exponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)
    @test size(grf.mean,1) == length(pts1)
    @test size(grf.mean,2) == length(pts2)
    @test all(grf.mean.==2.)

    @test_throws DimensionMismatch GaussianRandomField(2.0*ones(5,10),cov,Cholesky(),pts1,pts2)

end


text:
using CanvasWebIO, WebIO, InteractNext, Mux, Plots

inspectdr(show=false)
port = 8000

height = 800
width = 800
linesy = 8
linesx = 8
scaling = 1./400
style = Dict(:stroke=>"black", :strokeWidth=>"2")
svggrid = [[dom"svg:line[x1=0, y1=$(i*height/linesy), x2=$width, y2=$(i*height/linesy)]"(style=style)
            for i in 0:linesy]...,
           [dom"svg:line[x1=$(i*width/linesx), y1=0, x2=$(i*width/linesx), y2=$height]"(style=style)
            for i in 1:linesx]...,
           dom"svg:text[x = $(width/2-10), y=$(height+15)]"("1.0"),
           dom"svg:text[x = $(width-10), y=$(height+15)]"("2.0"),
           dom"svg:text[x = 0, y=$(height/2-5)]"("1.0"),
           dom"svg:text[x = 0, y=$(15)]"("2.0"),
          ]

canvas = Canvas([width+15, height+15])
marker1 = dom"svg:circle[id=marker1, cy=153, cx=273, r=10]"()
marker2 = dom"svg:rect[id=marker2, height=20, width=20, x=210, y=210]"()
addmovable!(canvas, marker1)
addmovable!(canvas, marker2)
for i in svggrid
    addstatic!(canvas, i)
end
ui = @manipulate for r in 0.1:0.1:10
    t = 0:0.02:2pi

    xf1 = canvas["marker1"][][1]*scaling
    yf1 = (height-canvas["marker1"][][2])*scaling
    plot(xf1*cos.(t), yf1*sin.(t), xlim = [-r,r], ylim=[-r,r])

    xf2 = canvas["marker2"][][1]*scaling
    yf2 = (height-canvas["marker2"][][2])*scaling
    plot(xf1*cos.(t), yf1*sin.(t), xlim = [-r,r], ylim=[-r,r])
    plot!([-xf2, -xf2, xf2, xf2, -xf2], [-yf2, yf2, yf2, -yf2, -yf2])
end

on(canvas["marker1"]) do val
    observe(r)[] = observe(r)[]
end

on(canvas["marker2"]) do val
    observe(r)[] = observe(r)[]
end

style = Dict(:display => "inline-table",
             :verticalAlign => "top",
             :width => "50%")

webio_serve(page("/",
                 req -> Node(:div,
                             Node(:div, canvas(), style=style),
                             Node(:div, ui, style=style))),
                 port)


text:
using CatViews

"""
    gcmvector{T, N}

gcmvector data structure that can be used for
  subsetting and indexing into a gcmarray.

```
gcmvector{T,N}(grid::gcmgrid,f::Array{Array{T,1},N},
         fSize::Array{NTuple{N, Int}},fIndex::Array{Int,1})

```
"""
struct gcmvector{T,N} <: AbstractMeshArray{T,N}
   grid::gcmgrid
   f::Array{Array{T,1},N}
   fSize::Array{NTuple{2, Int}}
   fIndex::Array{Int,1}
end

## Customize standard functions

import Base: findall, size, show, getindex, setindex!, copyto!

Base.size(A::gcmvector) = size(A.f)
Base.size(A::gcmvector, dim::Integer) = size(A)[dim]

function Base.getindex(A::gcmvector{T,N}, I::Vararg{Union{Int,Array{Int},AbstractUnitRange,Colon}, N}) where {T,N}
  return A.f[I...]
end

function Base.setindex!(A::gcmvector{T,N}, v, I::Vararg{Int, N}) where {T,N}
  return (A.f[I...] = v)
end

function Base.show(io::IO, z::gcmvector{T,N}) where {T,N}
    printstyled(io, " gcmvector \n",color=:normal)
    printstyled(io, "  grid type   = ",color=:normal)
    printstyled(io, "$(z.grid.class)\n",color=:blue)
    printstyled(io, "  data type   = ",color=:normal)
    printstyled(io, "$(eltype(z))\n",color=:blue)
    printstyled(io, "  tile array  = ",color=:normal)
    printstyled(io, "$(size(z))\n",color=:blue)
    printstyled(io, "  tile sizes  = ",color=:normal)
    printstyled(io, "$(size(z.f[1]))\n",color=:blue)
    for iFace=2:length(z.fIndex)
      printstyled(io, "                ",color=:normal)
      printstyled(io, "$(size(z.f[iFace]))\n",color=:blue)
    end
  return
end

function Base.view(A::gcmvector{T, N}, I::Vararg{Union{Int,AbstractUnitRange,Colon}, N}) where {T,N}
  tmpf=view(A.f,I...)
  return gcmvector{eltype(A),ndims(tmpf)}(A.grid,tmpf,A.fSize,A.fIndex)
end

## Customize broadcasting

Base.BroadcastStyle(::Type{<:gcmvector}) = Broadcast.ArrayStyle{gcmvector}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{gcmvector}}, ::Type{ElType}) where ElType
  # Scan the inputs for the gcmarray:
  A = find_gcmvector(bc)
  # Create the gcmvector output:
  return gcmvector{ElType,ndims(A)}(A.grid,similar(A.f),A.fSize,A.fIndex)
end

find_gcmvector(bc::Base.Broadcast.Broadcasted) = find_gcmvector(bc.args)
find_gcmvector(args::Tuple) = find_gcmvector(find_gcmvector(args[1]), Base.tail(args))
find_gcmvector(x) = x
find_gcmvector(a::gcmvector, rest) = a
find_gcmvector(::Any, rest) = find_gcmvector(rest)

# Specialize this method if all you want to do is specialize on typeof(dest)
@inline function copyto!(dest::gcmvector, bc::Broadcast.Broadcasted{Nothing})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && bc.args isa Tuple{AbstractArray} # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    bc′ = Broadcast.preprocess(dest, bc)
    @simd for I in eachindex(bc′)
        #@inbounds dest[I] = bc′[I]
        @inbounds dest[I] = gcmarray_getindex_evalf(bc′,I)
    end
    return dest
end

##

"""
    findall(A::gcmarray{Bool})

Return a gcmvector of the true indices in A. This allows:
```
findall(A.<0) #gcmvector of CartesianIndex{2}
A[findall(A.<0)] #gcmvector of eltype(A)
view(A,findall(A.<0)) #CatView of eltype(A)

A[findall(A.<0)]=B[findall(A.<0)]
A[findall(A.<0)].=view(B,findall(A.<0))
A[findall(A.<0)].=NaN
```
"""
function findall(A::gcmarray{Bool,N}) where {N}
  tmpOut=Array{Array{CartesianIndex{2},1},N}(undef,size(A))
  for a in eachindex(A)
    tmpOut[a]=findall(A[a])
  end
  return gcmvector{CartesianIndex{2},N}(A.grid,tmpOut,A.fSize,A.fIndex)
end

function Base.getindex(A::gcmarray{T,N}, B::gcmvector{CartesianIndex{2},N}) where {T,N}
  tmpOut=Array{Array{T,1},N}(undef,size(A))
  for a in eachindex(A)
    tmpOut[a]=A[a][B[a]]
  end
  return gcmvector{T,N}(A.grid,tmpOut,A.fSize,A.fIndex)
end

function Base.view(A::gcmarray{T,N}, B::gcmvector{CartesianIndex{2},N}) where {T,N}
   tmpOut=missing
   for a in eachindex(A)
     ismissing(tmpOut) ? tmpOut=view(A[a],B[a]) : tmpOut=CatView(tmpOut,view(A[a],B[a]))
   end
   return tmpOut
end

function Base.setindex!(A::gcmarray{T,N}, B::gcmvector{T,N}, C::gcmvector{CartesianIndex{2},N}) where {T,N}
  D=A
  for a in eachindex(A)
    D[a][C[a]]=B[a]
  end
  return (A = D)
end

function Base.setindex!(A::gcmarray{T,N}, v::Number, C::gcmvector{CartesianIndex{2},N}) where {T,N}
  D=A
  for a in eachindex(A)
    D[a][C[a]].=v
  end
  return (A = D)
end


text:
@testset "Chebyshev accuracy" begin
    for i in 1:100
        @test accuracy((θ,x)->evaluate(Chebyshev(length(θ)), θ, x), rand_θ(10),
                       linspace(-1,1,1001), [0.1])[1] ≥ 14
    end
    for i in 1:100
        @test accuracy((θ,x)->evaluate(Chebyshev(length(θ)), θ, Partial(x)), rand_θ(10),
                       linspace(-1,1,1001), [0.1])[1] ≥ 13
    end
end


text:
using QMTK
using Compat.Test
using Compat.Distributed

@testset "Metropolis-Hasting" begin

space = RealSpace(min=-4, max=4)
sampler = MHSampler(space)

@everywhere function normal(x::T) where T <: Real
    factor = 1 / sqrt(2 * pi)
    return factor * exp(- x^2/2)
end

plan = SamplePlan(normal, sampler)
data = sample!(plan)

@test length(data) == 1000

end


text:
#########
# Binning / elevation-bands
#########
x = -100.3:0.56:-12
y = 1001:2.3:1100
ele = x .+ y'
mask = trues(size(ele))
mask[1:10,:] = false
mask[[1,end],:] = false
mask[:,[1,end]] = false
g = Gridded(x,y,ele)
# negative step
bands, bandi = bin_grid(g, -10.0, mask)
@test length(bands)==length(bandi)
binmat = VAWTools.bins2matrix(g, bands, bandi)
@test all(binmat[1:10,:].==0)
@test all(binmat[:,1].==0)
# positive step
bands, bandi = bin_grid(g, 10.0, mask)
@test length(bands)==length(bandi)
binmat = VAWTools.bins2matrix(g, bands, bandi)
@test all(binmat[1:10,:].==0)
@test all(binmat[:,1].==0)
# adjusting step
bands, bandi = bin_grid(g, 200.0, mask)
@test length(bands)==length(bandi)
binmat = VAWTools.bins2matrix(g, bands, bandi)
@test all(binmat[1:10,:].==0)
@test all(binmat[:,1].==0)
# range steps:
bands_, bandi_ = bin_grid(g, bands, mask)
@test all(bands_.==bands)
@test all(bandi_.==bandi)
# non range steps:
@test_throws MethodError bands_, bandi_ = bin_grid(g, collect(bands), mask)
# @test all(bands_.==bands)
# @test all(bandi_.==bandi)

bands, bandi = bin_grid(g, 10.0, mask)
binmat = VAWTools.bins2matrix(g, bands, bandi)
xx,yy = -50:11:23, 900.0:7:1200
othergrid = Gridded(xx, yy, rand(length(xx),length(yy)))
bandii = VAWTools.bandi_for_other_grid(bands, bandi, binmat, g, othergrid)
bandii2 = VAWTools.bandi_for_other_grid(bands, bandi, g, othergrid)
@test bandii==bandii2

binmatt = VAWTools.bins2matrix(othergrid, bands, bandii)
@test binmatt == [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  6   7   7   8   9  10  10  11  12  12  13  14  14  15  0  0  0  0  0  0  0  0  0  0  0  0  0  0
                   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  7   8   9   9  10  11  11  12  13  13  14  15  15  16  0  0  0  0  0  0  0  0  0  0  0  0  0  0
                   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  8   9  10  10  11  12  12  13  14  14  15  16  17  17  0  0  0  0  0  0  0  0  0  0  0  0  0  0
                   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  9  10  11  11  12  13  14  14  15  16  16  17  18  18  0  0  0  0  0  0  0  0  0  0  0  0  0  0
                   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
                   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
                   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]

# 1D glacier
bands_g, bandi_g, malphas, areas, lengths, widths, x, xmid, dem, alpha2d =
    make_1Dglacier(g, 10.0, mask)
@test bands_g==bands
@test bandi_g==bandi
bands_g, bandi_g, malphas, areas, lengths, widths, x, xmid, dem, alpha2d =
    make_1Dglacier(g, bands_g, mask)
@test all(bands_g.==bands)
@test all(bandi_g.==bandi)

###############
# Fluxes
###############
tmp       = [14 14 15 15 16 17 18 18 19 0 0 0;
             14 14 15 15 16 17 18 19 19 0 0 0;
             14 15 15 16 16 17 18 19 19 0 0 0;
             14 15 15 16 17 18 18 19 0 0 0 0;
             15 15 16 16 17 18 19 19 0 0 0 0;
             15 16 16 17 17 18 19 20 0 0 0 0;
             15 16 16 17 18 17 20 21 0 0 0 0;
             16 16 17 18 18 17 21 22 23 0 0 0;
             16 17 17 18 19 20 21 23 24 0 0 0;
             16 17 18 19 19 20 22 23 0 0 24 24;
             17 18 19 19 20 21 22 0 23 24 24 24]-13;
elevation = zeros(size(tmp,1)+2,size(tmp,2)+2)-13
elevation[2:end-1,2:end-1] = tmp
mask = elevation.>0
bands, bandi = VAWTools.bin_grid(elevation,1,mask)
binmat = VAWTools.bins2matrix(elevation,bands,bandi)
# using PyPlot
# bm = binmat*1.0
# bm[bm.==0] = NaN
# matshow(bm',cmap="flag", origin="lower", extent=(1,size(bm,1)+1,1,size(bm,2)+1)); colorbar()
# xticks(1:size(bm,1))
# yticks(1:size(bm,2))
# grid()
@test all(elevation[mask].==binmat[mask])
@test all(elevation[(!).(mask)].==-13)
@test all(binmat[(!).(mask)].==0)
edges_at_boundaries = VAWTools.get_cells_on_boundary(bands, bandi, binmat)
@test sort(collect(keys(edges_at_boundaries)))== Tuple{Int64,Int64}[(1,0),(1,2),(2,0),(2,1),(2,3),(3,0),(3,2),(3,4),(4,0),(4,3),(4,5),(4,7),(4,8),(5,0),(5,4),(5,6),(6,0),(6,5),(6,7),(7,0),(7,4),(7,6),(7,8),(7,9),(8,0),(8,4),(8,7),(8,9),(8,10),(9,0),(9,7),(9,8),(9,10),(10,0),(10,8),(10,9),(10,11),(11,0),(11,10)]

boundaries = VAWTools.calc_boundaries(bands, bandi, binmat)
# symmetry
for i=1:length(bands), j=1:length(bands)
    if haskey(boundaries[j],i)
        @test length(boundaries[i][j])==length(boundaries[j][i])
    end
end
@test length(boundaries[2-1][3-1])==1
@test length(boundaries[3-1][2-1])==1
@test length(boundaries[3-1][4-1])==1
@test length(boundaries[4-1][5-1])==1
@test length(boundaries[4][3])==1
@test length(boundaries[4][5])==2
@test length(boundaries[4][7])==2
@test length(boundaries[4][8])==1
@test length(boundaries[5][4])==2
@test length(boundaries[5][6])==2
@test length(boundaries[6][5])==2
@test length(boundaries[8-1][7-1])==2
@test length(boundaries[8-1][9-1])==3
@test length(boundaries[9-1][8-1])==3
@test length(boundaries[9-1][10-1])==3
@test length(boundaries[9-1][11-1])==1
@test length(boundaries[10-1][9-1])==3
@test length(boundaries[10-1][11-1])==2
@test length(boundaries[10-1][1-1])==1
@test length(boundaries[11-1][1-1])==3
@test length(boundaries[11-1][10-1])==2
@test length(boundaries[11-1][12-1])==2
@test length(boundaries[12-1][11-1])==2
@test length(boundaries[12-1][1-1])==2

q1d = bands+1
u_trial = ones(elevation)
thick = sqrt.(elevation+13)
ux,uy = (-).(gradient3by3(1:13,1:14,elevation,mask))
dx = 1.0
lengths = rand(length(bands))
flux_dir_window = 4
window_frac = 0.4
plotyes=false
x=nothing
y=nothing
ubar, mask_ubar, facs = VAWTools._calc_u(q1d, boundaries, u_trial, thick,
                                                ux, uy, dx, mask, bands,
                                                flux_dir_window, # in [m]
                                                plotyes,x,y)

u2d, u2d_at_bands, scaling_factors_1d, mask_u2d_at_bands = VAWTools.calc_u(q1d, boundaries, u_trial, thick,
                                                                           ux, uy, dx, mask, bands, lengths,
                                                                           flux_dir_window, # in [m]
                                                                           window_frac)

# more tests would be good...
@test isequal(u2d_at_bands, ubar)

### Round trip 1D -> 2D -> 1D
#
# Note that 2D -> 1D -> 2D will not round-trip as information
# is lost in the first ->

x = -100.3:0.56:-12
y = 1001:2.3:1100
ele = 0.01*x.^2 .+ y';
ele[:, end÷2:end] += 120; # make a step, just to make it awkward
# contour(x,y,ele')
dem = Gridded(x,y,ele)
mask = sqrt.((x.+60).^2 .+ (y'.-1050).^2) .<45 #; ele[.!mask] = NaN; heatmap(ele)
maxs = [maximum(ele),maximum(ele[mask])]
mins = [minimum(ele),minimum(ele[mask])]

bands = reverse(mean(mins):5:maxs[1]+20)
vals1d = 0.4 * [149, 957, -777, -732, -260, 168, 491, 424, -709, 505, -52, 252, -6, 388, -98, -378, -568, -363, -885, 705, 761, 485, -346, -803, 950, -245, -253, 471, -303, 567, -108,
 773, -545, 23, -57, 524, -730, 521, -232, 600, 311, -862, -284, 151, 694, 541, 565, 707, 184, -189, 192, -471, -283, 293, -669, -244, -808, 155, -553, -392, -769, -940,
 828, -489, 863, 126]
# test that it's good input data
@test length(bands) == length(vals1d)
@test length(unique(vals1d))==length(vals1d)

bands_, bandi = VAWTools.bin_grid(dem, bands, mask, min_bin_number_ends=5);
@test length(bands_)==53
bands_, bandi = VAWTools.bin_grid(dem, bands, mask);
@test bands==bands_

# check that mask works:
@test !any(mask[vcat(bandi...)].==false)
mm = trues(mask)
mm[vcat(bandi...)] = false
@test !any(mask[mm])

vals2d = VAWTools.map_back_to_2D(size(dem.v), bandi, vals1d)

# test that all values are the same in each band
for bi in bandi
    for b in bi
        @test vals2d[bi[1]]==vals2d[b]
    end
end

vals1d_ = VAWTools.map_onto_bands(bandi, vals2d)
@test all(vals1d_[length.(bandi).>0] .≈ vals1d[length.(bandi).>0])
@test all(isnan.(vals1d_[length.(bandi).==0]))


text:
function transform(entity::HomogeneousCoordinates, coordinate_system::CanonicalToHartley, ℳ::Vector{<:AbstractArray})
    𝒪, 𝐓 = hartley_normalization(ℳ)
end

function transform(entity::CovarianceMatrices, coordinate_system::CanonicalToHartley, Λ::Vector{<:AbstractArray}, 𝐓::AbstractArray)

    Λ₂ = map(Λ) do 𝚲ₙ
         (𝐓 * 𝚲ₙ * 𝐓')
    end

end


text:
"""
        P16(X_cr::T=0.05, η_max::T=0.5) where T

Constant efficiency as in Pfrommer+ 2016, doi: 10.1093/mnras/stw2941
"""
struct P16{T} <: AbstractShockAccelerationEfficiency
    X_cr::T
    η_max::T

    P16(X_cr::T = 0.05, η_max::T = 0.5) where {T} = new{T}(X_cr, η_max)
end

"""
    η_Ms_acc(η_model::P16, M::T) where T

Constant efficiency as in Pfrommer+ 2016, doi: 10.1093/mnras/stw2941
"""
function η_Ms_acc(η_model::P16, M::Real)
    return 0.5
end

"""
    η_Ms_reacc(η_model::P16, M::T) where T

Constant efficiency as in Pfrommer+ 2016, doi: 10.1093/mnras/stw2941
"""
function η_Ms_reacc(η_model::P16, M::Real)
    return 0.5
end


text:
"""model definition of point2d"""

struct StatePoint2D <: StatePoint
    x::Float64
    y::Float64
end

function get_mid_status(p::StatePoint2D, q::StatePoint2D)::StatePoint2D
    return StatePoint2D((p.x + q.x) / 2, (p.y + q.y) / 2)
end

function dist(q1::StatePoint2D, q2::StatePoint2D)::Float64
    return norm([q1.x - q2.x, q1.y - q2.y])
end

function dist(q_from::StatePoint2D, q_to::StatePoint2D, o::CircleObstacle2D)::Float64
    return dist([q_from.x, q_from.y], [q_to.x, q_to.y], [o.x, o.y])
end

function dist(q_from::StatePoint2D, q_to::StatePoint2D, q_static::StatePoint2D)::Float64
    return dist([q_from.x, q_from.y], [q_to.x, q_to.y], [q_static.x, q_static.y])
end

function dist(q::StatePoint2D, o::CircleObstacle2D)::Float64
    return norm([q.x - o.x, q.y - o.y])
end

function dist(
    a_from::StatePoint2D,
    a_to::StatePoint2D,
    b_from::StatePoint2D,
    b_to::StatePoint2D,
)::Float64
    return dist(
        [a_from.x, a_from.y],
        [a_to.x, a_to.y],
        [b_from.x, b_from.y],
        [b_to.x, b_to.y],
    )
end

function gen_connect(
    q::StatePoint2D,  # to identify type
    obstacles::Vector{CircleObstacle2D},
    rads::Vector{Float64},
)::Function

    # check: q \in C_free
    f(q::StatePoint2D, i::Int64)::Bool = begin
        # within field
        any(x -> (x < rads[i] || 1 - rads[i] < x), [q.x, q.y]) && return false

        # obstacles
        any([dist(q, o) < o.r + rads[i] for o in obstacles]) && return false

        return true
    end

    f(q_from::StatePoint2D, q_to::StatePoint2D, i::Int64)::Bool = begin
        any(x -> (x < rads[i] || 1 - rads[i] < x), [q_to.x, q_to.y]) && return false

        # check: collisions with static obstacles
        any([dist(q_from, q_to, o) < o.r + rads[i] for o in obstacles]) && return false

        return true
    end

    return f
end

function gen_uniform_sampling(q::StatePoint2D)::Function
    () -> begin
        return StatePoint2D(rand(2)...)
    end
end

function plot_motion!(q_from::StatePoint2D, q_to::StatePoint2D, rad::Float64, params)
    plot!([q_from.x, q_to.x], [q_from.y, q_to.y]; params...)
end

function plot_start_goal!(q_init::StatePoint2D, q_goal::StatePoint2D, rad::Float64, params)
    plot!([q_init.x], [q_init.y]; markershape = :hex, params...)
    plot!([q_goal.x], [q_goal.y]; markershape = :star, params...)
end

function plot_agent!(q::StatePoint2D, rad::Float64, color::String)
    plot_circle!(q.x, q.y, rad, color, 3.0, fillalpha = 0.1)
end

function gen_random_instance_StatePoint2D(; params...)
    return gen_random_instance(StatePoint2D(0, 0); params...)
end


text:
# Author: Lukasz Bratos

# Funkcja znajduje minimalna wartosc w przedziale [1, 2],
# ktora nie spelna rownania x * 1/x == 1 w arytmetyce zmiennopozycyjnej.
function betweenOneAndTwo()
    x = one(Float64)
    while nextfloat(x) * (one(Float64) / nextfloat(x)) == one(Float64) && x < 2.0
        x = nextfloat(x)
    end
    return nextfloat(x)
end

# Funkcja do znajdowania najmniejszej dodatniej liczby
# ktora nie spelnia rownania x * 1/x == 1
function theSmallestOne()
    x = zero(Float64)
    while nextfloat(x) * (one(Float64) / nextfloat(x)) == one(Float64)
        x = nextfloat(x)
    end
    return nextfloat(x)
end

println("The smallest number in [1, 2] that doesn't satisfy the quation x * 1/x = 1: ", betweenOneAndTwo())
#println("The smallest number in [-inf, +inf] that doesn't satisfy the quation x * 1/x = 1: ", theSmallestOne())


text:
# Tests if RiskRegion and BruteForceRiskRegion define (approximately) the same region by
# testing for membership of both for many randomly generated points.

using LinearAlgebra: I
using Distributions

function verify_RiskRegion(dist::Sampleable, Ω::AbstractRiskRegion, Ω_bf::TailRiskScenGen.BruteForceRiskRegion, n::Int64)
    for i in 1:n
        y = rand(dist)
        ellipse_res = y ∈ Ω
        bf_res = y ∈ Ω_bf
        if  ellipse_res != bf_res
            println("\tERROR: failed at point ", i, " = ", y)
            println("\tEllipse results = ", ellipse_res)
            println("\tTransformed size = ", transformed_size(Ω, y))
            println("\tΩ.α = ", Ω.α)
            println("\tBrute Force results = ", bf_res)
            break
        end
    end
end

using Random
Random.seed!(1)

# Test RiskRegion for multivariate normal distribution

print("\tTest 1 - MvNormal...")
dim = 3
num_points = 1000
lattice_width = 100
β = 0.99

Z = Normal()
μ = rand(Z, dim)
A = rand(Z, dim, dim)
Σ = A'A
## μ = fill(0.0, dim)
## Σ = eye(dim)
dist = MvNormal(μ, Σ)
K = FiniteCone(Array{Float64}(I, dim, dim))

Ω_ellipse = EllipticalRiskRegion(dist, K, β)
Ω_bf = TailRiskScenGen.BruteForceRiskRegion(dist, K, β, lattice_width)

verify_RiskRegion(dist, Ω_ellipse, Ω_bf, num_points)

print("success\n")

# Test RiskRegion for multivariate T distribution

print("\tTest 2 - MvTDist...")
dim = 3
num_points = 1000
lattice_width = 100
β = 0.95

μ = randn(dim)
A = randn(dim, dim)/sqrt(dim)
Σ = A'A
df = 4.0
## μ = fill(0.0, dim)
## Σ = eye(dim)
dist = MvTDist(df, μ, Σ)
K = FiniteCone(Array{Float64}(I, dim, dim))
Ω_ellipse = EllipticalRiskRegion(dist, K, β)
Ω_bf = TailRiskScenGen.BruteForceRiskRegion(dist, K, β, lattice_width)

verify_RiskRegion(dist, Ω_ellipse, Ω_bf, num_points)
print("success\n")
# println("μ = ", μ)
# println("Σ = \n", Σ)