text:
using LoopVectorization.VectorizationBase: REGISTER_SIZE

pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
include(joinpath(LOOPVECBENCHDIR, "looptests.jl"))

const LIBCTEST = joinpath(LOOPVECBENCHDIR, "libctests.so")
const LIBFTEST = joinpath(LOOPVECBENCHDIR, "libftests.so")

# requires Clang with polly to build
cfile = joinpath(LOOPVECBENCHDIR, "looptests.c")
if !isfile(LIBCTEST) || mtime(cfile) > mtime(LIBCTEST)
    run(`clang -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -lm -mllvm -polly -mllvm -polly-vectorizer=stripmine -shared -fPIC $cfile -o $LIBCTEST`)
end
ffile = joinpath(LOOPVECBENCHDIR, "looptests.f90")
if !isfile(LIBFTEST) || mtime(ffile) > mtime(LIBFTEST)
    # --param max-unroll-times defaults to ≥8, which is generally excessive
    run(`gfortran -Ofast -march=native -funroll-loops --param max-unroll-times=4 -floop-nest-optimize -mprefer-vector-width=$(8REGISTER_SIZE) -shared -fPIC $ffile -o $LIBFTEST`)
end

for order ∈ (:kmn, :knm, :mkn, :mnk, :nkm, :nmk)
    gemm = Symbol(:gemm_, order)
    @eval function $(Symbol(:c, gemm, :!))(C, A, B)
        M, N = size(C); K = size(B, 1)
        ccall(
            ($(QuoteNode(gemm)), LIBCTEST), Cvoid,
            (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
            C, A, B, M, K, N
        )
    end
    @eval function $(Symbol(:f, gemm, :!))(C, A, B)
        M, N = size(C); K = size(B, 1)
        ccall(
            ($(QuoteNode(gemm)), LIBFTEST), Cvoid,
            (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
            C, A, B, Ref(M), Ref(K), Ref(N)
        )
    end
end

function fgemm_builtin!(C, A, B)
    M, N = size(C); K = size(B, 1)
    ccall(
        (:gemmbuiltin, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
        C, A, B, Ref(M), Ref(K), Ref(N)
    )
end
function cAtmulB!(C, A, B)
    M, N = size(C); K = size(B, 1)
    ccall(
        (:AtmulB, LIBCTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C, A, B, M, K, N
    )
end
function fAtmulB!(C, A, B)
    M, N = size(C); K = size(B, 1)
    ccall(
        (:AtmulB, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
        C, A, B, Ref(M), Ref(K), Ref(N)
    )
end
function fAtmulB_builtin!(C, A, B)
    M, N = size(C); K = size(B, 1)
    ccall(
        (:AtmulBbuiltin, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
        C, A, B, Ref(M), Ref(K), Ref(N)
    )
end

function cdot(a, b)
    N = length(a)
    ccall(
        (:dot, LIBCTEST), Float64,
        (Ptr{Float64}, Ptr{Float64}, Clong),
        a, b, N
    )
end
function fdot(a, b)
    N = length(a)
    d = Ref{Float64}()
    ccall(
        (:dot, LIBFTEST), Cvoid,
        (Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}),
        d, a, b, Ref(N)
    )
    d[]
end
function cselfdot(a)
    N = length(a)
    ccall(
        (:selfdot, LIBCTEST), Float64,
        (Ptr{Float64}, Clong),
        a, N
    )
end
function fselfdot(a)
    N = length(a)
    d = Ref{Float64}()
    ccall(
        (:selfdot, LIBFTEST), Cvoid,
        (Ref{Float64}, Ptr{Float64}, Ref{Clong}),
        d, a, Ref(N)
    )
    d[]
end
function cdot3(x, A, y)
    M, N = size(A)
    ccall(
        (:dot3, LIBCTEST), Float64,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        x, A, y, M, N
    )
end
function fdot3(x, A, y)
    M, N = size(A)
    d = Ref{Float64}()
    ccall(
        (:dot3, LIBFTEST), Cvoid,
        (Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
        d, x, A, y, Ref(M), Ref(N)
    )
    d[]
end

function cgemv!(y, A, x)
    M, K = size(A)
    ccall(
        (:gemv, LIBCTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y, A, x, M, K
    )
end
function fgemv!(y, A, x)
    M, K = size(A)
    ccall(
        (:gemv, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
        y, A, x, Ref(M), Ref(K)
    )
end
function fgemv_builtin!(y, A, x)
    M, K = size(A)
    ccall(
        (:gemvbuiltin, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
        y, A, x, Ref(M), Ref(K)
    )
end

function caplusBc!(D, a, B, c)
    M, K = size(B)
    ccall(
        (:aplusBc, LIBCTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        D, a, B, c, M, K
    )
end
function faplusBc!(D, a, B, c)
    M, K = size(B)
    ccall(
        (:aplusBc, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
        D, a, B, c, Ref(M), Ref(K)
    )
end
function cOLSlp(y, X, β)
    N, P = size(X)
    ccall(
        (:OLSlp, LIBCTEST), Float64,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y, X, β, N, P
    )
end
function fOLSlp(y, X, β)
    N, P = size(X)
    lp = Ref{Float64}()
    ccall(
        (:OLSlp, LIBFTEST), Cvoid,
        (Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
        lp, y, X, β, Ref(N), Ref(P)
    )
    lp[]
end
function cvexp!(b, a)
    N = length(b)
    ccall(
        (:vexp, LIBCTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Clong),
        b, a, N
    )
end
function fvexp!(b, a)
    N = length(b)
    ccall(
        (:vexp, LIBFTEST), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ref{Clong}),
        b, a, Ref(N)
    )
end
function fvexpsum(a)
    N = length(a)
    s = Ref{Float64}()
    ccall(
        (:svexp, LIBFTEST), Cvoid,
        (Ref{Float64}, Ptr{Float64}, Ref{Clong}),
        s, a, Ref(N)
    )
    s[]
end


text:
import ..BifurcationsBase

dim_ds_state(sweep::Codim1LCSweep) = length(sweep.eigvals[1])  # TODO: don't

abstract type AbstractLimitCycleData end

struct LimitCycleData{M, R, V, P} <: AbstractLimitCycleData
    state::M
    period::R
    param_value::V
    prob::P
end

BifurcationsBase.problem_of(lc::LimitCycleData) = lc.prob
BifurcationsBase.contkind(lc::LimitCycleData) = contkind(lc.prob)  # TODO: don't

LimitCycleData(u::AbstractVector, n::Integer, prob) =
    LimitCycleData(
        reshape((@view u[1:end-2]), n, :),
        u[end - 1],  # period
        u[end],      # param_value
        prob,
    )

function limitcycles(sweep::Codim1LCSweep)
    n = dim_ds_state(sweep)
    u_list = as(sweep, ContinuationSweep).u
    sol = as(sweep, ContinuationSweep).sol.value  # TODO: don't
    if sol !== nothing
        prob = sol.prob
    else
        prob = nothing
    end
    return LimitCycleData.(u_list, (n,), (prob,))
end

limitcycles(sol::Codim1LCSolution) = vcat(limitcycles.(sol.sweeps)...)
limitcycles(solver::Codim1LCSolver) = limitcycles(solver.sol)

function BifurcationsBase.measure(lc::AbstractLimitCycleData, i::Integer)
    state = @view lc.state[i, :]
    xs = similar(state, length(state) + 1)
    copyto!(xs, state)
    xs[end] = state[1]
    return xs
end

# TODO: make it more type-based
function BifurcationsBase.measure(lc::LimitCycleData, key::Symbol)
    if key in (:p1, :parameter)
        return lc.param_value
    elseif key == :period
        return lc.period
    end
    error("Unsupported key: $key")
end

abstract type LCMeasurement end

struct GenericMeasurement{F} <: LCMeasurement
    f
end

(key::GenericMeasurement)(lc::LimitCycleData) = key.f(lc)

"""
Coordinate-wise measurement of the state.
"""
struct CWStateMeasurement{F} <: LCMeasurement
    index::Int
    f::F
end

(key::CWStateMeasurement)(lc::LimitCycleData) =
    key.f(@view lc.state[key.index, :])

function BifurcationsBase.measure(sweep::Codim1LCSweep, key::LCMeasurement)
    return key.(limitcycles(sweep))
end

# TODO: make it more type-based
function BifurcationsBase.measure(sweep::Codim1LCSweep, key::Symbol)
    u_list = as(sweep, ContinuationSweep).u
    if key in (:p1, :parameter)
        return [u[end] for u in u_list]
    elseif key == :period
        return [u[end - 1] for u in u_list]
    end
    error("Unsupported key: $key")
end

#=
period_extrema(ctx) = nan_(extrema, lc.period for lc in limitcycles(ctx))
=#

ext: jl


text: mango(quantity, price) = (quantity - quantity ÷ 3) * price

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


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

using Markdown
using InteractiveUtils

# ╔═╡ 63b9e832-3b3b-11eb-0d45-b54bff3a4c10
using Pkg; Pkg.activate("."); using Plots

# ╔═╡ a79e42e0-3a31-11eb-3791-2dde6b25172b
md"
# Proyecto 1

## Descripción
Se escribirá un código para visualizar diferentes familias del conjunto de Julia y Mandelbrot, en el proceso se modelarán diversas formas biológicas. Ambos conjuntos se forman a partir del estudio de ciertas funciones, $f_c$, sobre las cuales se iteran los números complejos. La notación para la función, $f_c$, denota que ésta depende de un parámetro $c \in C$.

El conjunto de Julia $J_{f_c}$, para una constante $c$ particular, está formado por todos aquellos puntos del plano complejo en el cual la sucesión siguiente es acotada:

$$z_0=z$$
$$z_{n+1}=f_c(z_n)$$

El conjunto de Mandelbrot $M_{f_c}$, de forma similar, está formado por el conjunto de puntos del plano complejo, en el cual $c$ pertenece al conjunto si la siguiente sucesión es acotada:

$$z_0=0$$
$$z_{n+1}=f_c(z_n)$$

El conjunto de puntos que resulta de la sucesión, para cada $z$ implicado, se conoce como la órbita de $z$ por $f_c$. Para determinar los elementos de una familia del conjunto de $J_{f_c}$ y $M_{f_c}$, se utiliza un método iterativo, verificando el criterio convergencia siguiente, en cada paso de la iteración donde se obtiene un elemento de la órbita:

$$|z_n|<2$$

Si en algún punto de la iteración no se cumple el criterio, el punto del plano complejo analizado no pertenece al conjunto.

El conjunto de Julia y Mandelbrot están dados por una familia muy popular de sistemas dinámicos complejos, esto es la familia de polinomios cuadráticos complejos:

$$f_c=z^2+c$$

Formas biológicas han sido representadas mediante procesos iterativos, como el descrito anteriormente. Pero, empleando un conjunto de funciones $(f_c)$ más complejas que las estudiadas en los conjuntos de Julia y Mandelbrot. A continuación, se listan algunas de las funciones usadas para reproducir formas biológicas:

$$f_1(z,c)=\sin(z)+z^2+c$$
$$f_2(z,c)=z^z+z^6+c$$
$$f_3(z,c)=z^z+z^6+c$$
$$f_4(z,c)=z^5+c$$
$$f_5(z,c)=z^3+c$$

La iteración es realizada hasta que la magnitud de la componente real e imaginaria de los elementos de la órbita sobrepasa un determinado valor $\tau$ (un valor puede ser $\tau=100$), esto significa que el proceso iterativo se mentiene mientras se cumpla que la parte real o imaginaria satisfacen la condición:

$$|Re(z_n)|<\tau$$

$$|Im(z_n)|<\tau$$

Una manera de graficar los conjuntos de Julia, Mandelbrot y las formas biológicas, es graficar solamente los puntos que pertenecen al conjunto con un color determinado (monocromático).

Para agregar colores a los gráficos, generalmente se utiliza el algoritmo de: Tiempo de Escape, aunque no es el único. El algoritmo consiste en asignar un color dependiendo del número de iteraciones realizadas antes que el elemento de la órbita no satisfaga el criterio de convergencia. De esta manera se visualiza aquellos elementos que pertenecen al conjunto y los que no, además, revela la velocidad de divergencia de la órbita en cada punto.

## Requisitos
Para el desarrollo del proyecto deberá utilizar solamente los elementos abarcados en la primera semana del workshop, esto es, el uso de: estructuras de datos y control, tipos de datos primitivos del lenguaje. Por lo que **no utilizará** ningún paquete adicional, creado o desarrollado por la comunidad de Julia, a excepción del paquete Plots.

Para completar este proyecto deberá realizar los siguientes items:
1. Graficar un conjunto de Julia, para un $c$ arbitrario (monocromático).
2. Graficar el conjunto de Mandelbrot (monocromático).
3. Graficar una forma biológica, par un c arbitrario (monocromático).
4. Hacer una gráfica a color del conjunto de Mandelbrot usando el algoritmo de Tiempo de Escape.
5. Realizar una gráfica a color para un conjunto de Julia, para un $c$ arbitrario .
6. Realizar una gráfica a color para una forma biológica.

Para realizar las gráficas, deberá tomar en cuenta el uso de 'structs' para representar estructuras de datos particulares y también el uso de 'multiple dispatch'.

## Criterios de Evaluación
Durante la evaluación del proyecto se tomarán en cuenta los siguientes aspectos:
- Cumplimiento de los requerimientos mencionados.
- Claridad y orden el código presentado.
- Completitud del proyecto, haciendo todas las partes que lo componen.
- Analizar y comentar el trabajo de otros grupos (Mediante el peer review).
"

# ╔═╡ 231273f6-3b3b-11eb-2b3b-2308e244a358
md"
---------
### AYUDA
---------
#### Gráficas
Para realizar las gráficas se utilizará el paquete Plots, por lo tanto el paquete debe ser importado:
"

# ╔═╡ 16040e00-539c-11eb-0f12-db5ae4d9970f
Plots.default(size = (650,450))

# ╔═╡ cd1023fa-3b3b-11eb-22ec-2724be9b2af8
md"
A continuación se mostrarán ejemplos para realizar gráficas con el paquete, los cuales pueden usarse para mostrar las formas biológicas y los conjuntos de Julia y Mandelbrot, una vez se hayan obtenido los datos necesarios.

Gráfico de números complejos:
"

# ╔═╡ 4070b426-3b3d-11eb-1ef2-73f2f10c9ce6
begin
	r = 2
	zcomplejo = [r*cos(θ)+r*sin(θ)*im for θ in 0:0.3:2*pi]
	scatter(zcomplejo, seriescolor=:white,
		    markerstrokecolor=:blue,
		    aspectratio=1,
		    title="Gráfica de ejemplo")
end

# ╔═╡ 5a70e2f4-3b3f-11eb-1dfe-a777db142c6b
md"
Para colorear los gráficos del conjunto de Julia, Mandelbrot y las formas biológicas, utilizando el algoritmo de Tiempo de Escape,  se pueden utilizar mapas de calor.
"

# ╔═╡ 48238fda-3b3e-11eb-00c6-278430e8d319
begin
	x = -2:0.3:2
	y = -2:0.3:2
	colores = [i₁*i₂ for i₁ in x, i₂ in y]
	heatmap(x, y, colores, color=cgrad([:black,:blue,:white]),
		    title="Ejemplo Mapa de Calor", xlabel="Re(z)", ylabel="Im(z)")
end

# ╔═╡ 0b52804a-3b41-11eb-3178-c5d7e9b3dca7
md"
Para el mapa de calor, se necesita hacer una malla (las listas 'x' y 'y' del ejemplo anterior) y una matriz donde se guardan los valores que definen el color para cada punto de la malla.

#### Pseudocódigo
Para clarificar los pasos a realizar, se muestra el siguiente algoritmo, cuando se desee obtener un gráfico con colores:

*NOTAS:*
- Se utilizan algunas funciones que pertenecen al core y base de Julia, si no conoce o intuye su uso averigue sobre la función mediante `?<nombre_función>` o `methods(<nombre_función>)`.
- Usar `test_convergencia = b` para la creación de formas biológicas.
- Las variables rz, iz son las componentes real imaginaria del número complejo respectivamente, que define un punto de la región del plano complejo por analizar.
- La variable i es el contador del número de iteraciones.
- c y z representan números complejos.

~~~julia
1  Para rz = -2:0.3:2 Hacer #Elección del rango puede ser diferente
2  	  Para iz = -2:0.3:2 Hacer #Elección del rango puede ser diferente
3		z = complex(rz,iz)
4		#La precisión de la imagen es mejor si incrementa las iteraciones
5		Para i = 1:30
6			z = f(z) + c
7			#Se verfica que z = complex(rz,iz) esté en el conjunto
8			#T depende del conjunto que se desea graficar
9			Si abs(z) > T Entonces break Fin
10	   	Fin
11		color = i #Indica la velocidad de divergencia de la orbita
12		Si test_convergencia = a Entonces
13			Graficar punto en (rz,iz,color)
14		Fin
15		Si test_convergencia = b Entonces
16			Si real(z) < T || imag(z) < T Entonces
17				Graficar punto en (rz,iz,color)
18			Fin
19		Fin
20	  Fin
21  Fin
~~~
"

# ╔═╡ c647eeb2-3b44-11eb-2637-9197f59f6490
md"
# Solución
"

# ╔═╡ ce1cb762-4629-11eb-0598-15119d777749
md"
La solución del problema planteado puede seguir un camino diferente, pero siempre basado en el algoritmo proporcionado. Para completar la parte obligatoria de este proyecto rellene las partes que están incompletas de este documento, tome en cuenta la signatura de las funciones descritas y lo que esta debe hacer. A medida avance en el documento haga uso de las definiciones previas donde sea conveniente.
"

# ╔═╡ b667cda0-46e0-11eb-370d-41d8d7cdf7de
md"
## Definición de objetos y funciones
"

# ╔═╡ 56c3949c-462c-11eb-058f-ef53a7395810
md"
+ A continuación construya un struct `Grid` (Grilla), donde se defina la región del plano complejo que debe ser analizado (límites en las abcisas y ordenadas) y el espacio entre cada punto de la grilla.
"

# ╔═╡ f915f344-3a00-11eb-16a9-6be5b192810c
"""*
 Representa un objeto grilla
"""
struct Grid{T<:Real}

end

# ╔═╡ f1207c62-4631-11eb-30cb-691cbcca8574
md"
+ Defina una función que cree una matriz de rangos para los ejes de las abcisas y ordenas.
"

# ╔═╡ 1dc63efc-39d6-11eb-3b8d-15da426e82a1
"""
	makeGrid(g::Grid)
Crear una matriz de rangos [rango_abcisas,rango_ordenadas]
"""
function makeGrid(g::Grid)

end

# ╔═╡ 774af954-4635-11eb-31dd-17a273bf9b15
md"
+ Defina funciones de una línea para las formas biológicas y la familia de polinomios cuadráticos complejos
"

# ╔═╡ acda722e-4636-11eb-0a76-29912beacb1c
begin #*

end

# ╔═╡ 36bed290-4633-11eb-26fe-d15d7d6f03f8
md"
+ Defina una función que retorne `true` si el criterio de convergencia para los conjuntos $J_{_c}$ y $M_{f_c}$ se cumple, `false` en caso contrario.
"

# ╔═╡ 60604b20-4637-11eb-1874-49e2756715ac
""" *
	testJM(z::Complex)
Comprobar el criterio de convergencia para los conjuntos de Julia y Mandelbrot
"""
function testJM(z::Complex)

end

# ╔═╡ 5466484a-4637-11eb-115f-ab417bd2b802
md"
+ Defina una función que `true` si el criterio de convergencia para las formas biológicas se cumple, `false` en caso contrario.
"

# ╔═╡ 6b6e04f0-4638-11eb-1401-6db76349254b
"""
	testbiomorph(z::Complex,τ::Real)
Comprobar el criterio de convergencia para los conjuntos de Julia y Mandelbrot
"""
function testbiomorph(z::Complex,τ::Real)

end

# ╔═╡ 71f98b8a-4641-11eb-2c9c-7f49f8577365
md"
+ Defina una función que itere un número sobre una función $f_c$ y devuelva true si se cumple el criterio de divergencia para $J_{f_c}$, $M_{f_c}$ y formas biológicas
"

# ╔═╡ c84bee02-4640-11eb-0dc2-57337e999b50
"""
	iterate(test::Function,f::Function,z::Complex,iter::Integer)
Interar z sobre una funcion f
"""
function iterate(test::Function,f::Function,z::Complex,iter::Integer)

end

# ╔═╡ fb3a891a-4634-11eb-2228-67741e0f82e2
md"
+ Defina un función que devuelva el número de iteraciones realizadas para un dado z (número complejo), hasta que se no se cumple el criterio de divergencia.
"

# ╔═╡ 828790ba-39d9-11eb-03bf-cbf806cba31a
"""*
	colormap(f::Function,test::Function,z::Complex,c::Complex,iter::Integer)
Retornar el número de iteraciones para un valor dado de z minetras un criterio de convergencia sea válido
"""
function colormap(f::Function,test::Function,z::Complex,c::Complex,iter::Integer)

end

# ╔═╡ 393b931a-46a0-11eb-361d-8da678672d65
md"+ Usando multiple dispatch, defina nuevamente la función colormap agregando un nuevo parámetro de entrada $\tau$, de modo que se pueda usar el test de convergencia para la forma biológica"

# ╔═╡ 1c510af2-46a0-11eb-08ac-a1d54104fd11
function colormap(f::Function,test::Function,z::Complex,c::Complex,iter::Integer,τ::Integer)

end

# ╔═╡ 24bba08a-46e2-11eb-0e42-412ca935f3e7
md"
## Solución a los incisos
"

# ╔═╡ 94d3f144-46e2-11eb-1f3b-7f0ba366c99a
md"
Usando las definiciones anteriores realice lo indicado.
"

# ╔═╡ c7030f70-469e-11eb-1fc0-03060512aa15
md"
##### 1) Graficar un conjunto de Julia, para un $c$ arbitrario (monocromático).
"

# ╔═╡ 3f322d18-463b-11eb-1560-0dfa38f7f644
md"
+ Defina una funcion que retorne el conjunto de puntos que pertenecen al conjunto de Julia.
"

# ╔═╡ 73644592-463b-11eb-09c0-952848e159e3
"""
	setjulia()
Construir un array con los puntos del plano complejo que pertenecen al conjunto, dada la región de análisis del plano complejo
Entradas:
- f::Function
- test::Function
- grid::Array{T,2} where T
- c::Complex
- iter::Integer
"""
function setjulia(
		f::Function,
		test::Function,
		grid::Array{T,2} where T,
		c::Complex,
		iter::Integer)

end

# ╔═╡ 553b382e-469e-11eb-3b55-3d189b08538c
begin

end

# ╔═╡ 639f176a-46e3-11eb-0d26-21e096ca0863
md"
##### 2) Graficar el conjunto de Mandelbrot (monocromático).
"

# ╔═╡ 687d2356-463b-11eb-33a6-53321e36f77a
md"
+ Defina una funcion que retorne el conjunto de puntos que pertenecen al conjunto de Mandelbrot
"

# ╔═╡ ce8536ce-469d-11eb-3f57-59b7ce9ac946
"""
	setmandelbrot()
Construir un array con los puntos del plano complejo que pertenecen al conjunto, dada la región de análisis del plano complejo
Entradas:
- f::Function
- test::Function
- grid::Array{T,2} where T
- iter::Integer
"""
function setmandelbrot(
		f::Function,
		test::Function,
		grid::Array{T,2} where T,
		iter::Integer)

end

# ╔═╡ fc4feb92-3b5d-11eb-28ff-630169bd8fc2
begin

end

# ╔═╡ 847b1e16-46e3-11eb-1581-95c456a7482f
md"
##### 3) Graficar una forma biológica, par un c arbitrario (monocromático).
"

# ╔═╡ e270a720-469f-11eb-1608-57230535b389
md"
+ Defina una funcion que retorne el conjunto de puntos que pertenecen al conjunto de una forma biológica.
"

# ╔═╡ f4f08ecc-469f-11eb-0d11-91bbe2f63dc8
"""
	setbiomorph()
Construir un array con los puntos del plano complejo que pertenecen al conjunto, dada la región de análisis del plano complejo
Entradas:
- f::Function
- test::Function
- grid::Array{T,2} where T
- c::Complex
- iter::Integer
- τ::Integer
"""
function setbiomorph(
		f::Function,
		test::Function,
		grid::Array{T,2} where T,
		c::Complex,
		iter::Integer,
		τ::Integer)

end

# ╔═╡ 9772cb98-469f-11eb-00b0-312e5a3b75dd
begin

end

# ╔═╡ b27d4370-46e3-11eb-0b7d-fbe40c187c17
md"
##### 4) Hacer una gráfica a color del conjunto de Mandelbrot

Usando el algoritmo de Tiempo de Escape y multiple dispatch, defina nuevamente la función setmandelbrot, para obtener los datos necesarios (matriz) que defina los colores para cada punto de una región arbitraria.
"

# ╔═╡ 5af85e8e-46c8-11eb-3de2-cf62fcb0614f
"""
	setmandelbrot(f::Function,test::Function,grid::Array{T,1} where T,iter::Integer)
Crear una matriz de enteros correspondiente a los colores
"""
function setmandelbrot(
		f::Function,
		test::Function,
		grid::Array{T,1} where T,
		iter::Integer
	)

end

# ╔═╡ e4797f76-40c9-11eb-02d6-35073745bec0
begin

end

# ╔═╡ e531a9be-46e3-11eb-0219-0bb7ff3871c2
md"
##### 5) Realizar una gráfica a color para un conjunto de Julia, para un $c$ arbitrario.

Usando el algoritmo de Tiempo de Escape y multiple dispatch, defina nuevamente la función setjulia, para obtener los datos necesarios (matriz) que defina los colores para cada punto de una región arbitraria.
"

# ╔═╡ af6e5546-46d6-11eb-21cf-ef7281e43d44
begin

end

# ╔═╡ 27c6bc88-46e4-11eb-2f01-9531af826324
md"
##### 6) Realizar una gráfica a color para una forma biológica.

Usando el algoritmo de Tiempo de Escape y multiple dispatch, defina nuevamente la función setbiomorph, para obtener los datos necesarios (matriz) que defina los colores para cada punto de una región arbitraria.
"

# ╔═╡ 278afcae-46dd-11eb-2551-b7e85c38f1f2
begin

end

# ╔═╡ c05869d8-46e4-11eb-1a10-e71db604ca8c
md"
##### 7) Completar un ejercicio que forma parte de los  challenges problems.
"

# ╔═╡ 3ad65006-508d-11eb-0e27-b1bbeffc47aa

# ╔═╡ Cell order:
# ╟─a79e42e0-3a31-11eb-3791-2dde6b25172b
# ╟─231273f6-3b3b-11eb-2b3b-2308e244a358
# ╠═63b9e832-3b3b-11eb-0d45-b54bff3a4c10
# ╠═16040e00-539c-11eb-0f12-db5ae4d9970f
# ╟─cd1023fa-3b3b-11eb-22ec-2724be9b2af8
# ╠═4070b426-3b3d-11eb-1ef2-73f2f10c9ce6
# ╟─5a70e2f4-3b3f-11eb-1dfe-a777db142c6b
# ╠═48238fda-3b3e-11eb-00c6-278430e8d319
# ╟─0b52804a-3b41-11eb-3178-c5d7e9b3dca7
# ╟─c647eeb2-3b44-11eb-2637-9197f59f6490
# ╟─ce1cb762-4629-11eb-0598-15119d777749
# ╟─b667cda0-46e0-11eb-370d-41d8d7cdf7de
# ╟─56c3949c-462c-11eb-058f-ef53a7395810
# ╠═f915f344-3a00-11eb-16a9-6be5b192810c
# ╟─f1207c62-4631-11eb-30cb-691cbcca8574
# ╠═1dc63efc-39d6-11eb-3b8d-15da426e82a1
# ╟─774af954-4635-11eb-31dd-17a273bf9b15
# ╠═acda722e-4636-11eb-0a76-29912beacb1c
# ╟─36bed290-4633-11eb-26fe-d15d7d6f03f8
# ╠═60604b20-4637-11eb-1874-49e2756715ac
# ╟─5466484a-4637-11eb-115f-ab417bd2b802
# ╠═6b6e04f0-4638-11eb-1401-6db76349254b
# ╟─71f98b8a-4641-11eb-2c9c-7f49f8577365
# ╠═c84bee02-4640-11eb-0dc2-57337e999b50
# ╟─fb3a891a-4634-11eb-2228-67741e0f82e2
# ╠═828790ba-39d9-11eb-03bf-cbf806cba31a
# ╟─393b931a-46a0-11eb-361d-8da678672d65
# ╠═1c510af2-46a0-11eb-08ac-a1d54104fd11
# ╟─24bba08a-46e2-11eb-0e42-412ca935f3e7
# ╟─94d3f144-46e2-11eb-1f3b-7f0ba366c99a
# ╟─c7030f70-469e-11eb-1fc0-03060512aa15
# ╟─3f322d18-463b-11eb-1560-0dfa38f7f644
# ╠═73644592-463b-11eb-09c0-952848e159e3
# ╠═553b382e-469e-11eb-3b55-3d189b08538c
# ╟─639f176a-46e3-11eb-0d26-21e096ca0863
# ╟─687d2356-463b-11eb-33a6-53321e36f77a
# ╠═ce8536ce-469d-11eb-3f57-59b7ce9ac946
# ╠═fc4feb92-3b5d-11eb-28ff-630169bd8fc2
# ╟─847b1e16-46e3-11eb-1581-95c456a7482f
# ╟─e270a720-469f-11eb-1608-57230535b389
# ╠═f4f08ecc-469f-11eb-0d11-91bbe2f63dc8
# ╠═9772cb98-469f-11eb-00b0-312e5a3b75dd
# ╟─b27d4370-46e3-11eb-0b7d-fbe40c187c17
# ╠═5af85e8e-46c8-11eb-3de2-cf62fcb0614f
# ╠═e4797f76-40c9-11eb-02d6-35073745bec0
# ╟─e531a9be-46e3-11eb-0219-0bb7ff3871c2
# ╠═af6e5546-46d6-11eb-21cf-ef7281e43d44
# ╟─27c6bc88-46e4-11eb-2f01-9531af826324
# ╠═278afcae-46dd-11eb-2551-b7e85c38f1f2
# ╟─c05869d8-46e4-11eb-1a10-e71db604ca8c
# ╠═3ad65006-508d-11eb-0e27-b1bbeffc47aa


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

abstract type AbstractTerminationCriterion end
abstract type AbstractTermination end

terminate(quasigradient::AbstractQuasiGradient, k::Integer, f::Float64, x::AbstractVector, ∇f::AbstractVector) = terminate(quasigradient.criterion, k, f, x, ∇f)

"""
    RawTerminationParameter

An optimizer attribute used for raw parameters of the termination criterion. Defers to `RawParameter`.
"""
struct RawTerminationParameter <: TerminationParameter
    name::Any
end

function MOI.get(termination::AbstractTermination, param::RawTerminationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(termination.parameters)))
        error("Unrecognized parameter name: $(name) for termination $(typeof(termination)).")
    end
    return getfield(termination.parameters, name)
end

function MOI.set(termination::AbstractTermination, param::RawTerminationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(termination.parameters)))
        error("Unrecognized parameter name: $(name) for termination $(typeof(termination)).")
    end
    setfield!(termination.parameters, name, value)
    return nothing
end

include("iteration.jl")
include("objective.jl")
include("gradient.jl")


text:
include("abstract_learner.jl")
include("tabular_learner.jl")


text:
function koch(t, x)
    if x < 3
        forward(t, x)
    else
        koch(t, x / 3)
        turn(t, -60)
        koch(t, x / 3)
        turn(t, 120)
        koch(t, x / 3)
        turn(t, -60)
        koch(t, x / 3)
    end
end


text:
module SelectionUI

import REPL
import REPL.TerminalMenus
using REPL.TerminalMenus: request  # Called from parent package

using ..MemoryInspector: SelectionOptions

mutable struct InspectMenu <: TerminalMenus.AbstractMenu
    options::Vector{String}
    pagesize::Int
    pageoffset::Int
    selected::Int
    toggle::Union{Nothing, Symbol}
    selected_command::Union{Nothing, SelectionOptions.Option}
    scroll_horizontal::Int
end

function show_as_line(el)
    reduced_displaysize = displaysize(stdout) .- (0, 3)
    buf = IOBuffer()
    show(IOContext(buf, :limit=>true, :displaysize=>reduced_displaysize), el)
    String(take!(buf))
end

function InspectMenu(options; pagesize::Int=10)
    #options = vcat(map(show_as_line, callsites), ["↩"])
    length(options) < 1 && error("InspectMenu must have at least one option")

    # if pagesize is -1, use automatic paging
    pagesize = pagesize == -1 ? length(options) : pagesize
    # pagesize shouldn't be bigger than options
    pagesize = min(length(options), pagesize)
    # after other checks, pagesize must be greater than 1
    pagesize < 1 && error("pagesize must be >= 1")

    pageoffset = 0
    selected = -1 # none

    scroll_horizontal = 0

    InspectMenu(options, pagesize, pageoffset, selected, nothing, nothing, scroll_horizontal)
end

TerminalMenus.options(m::InspectMenu) = m.options
TerminalMenus.cancel(m::InspectMenu) = m.selected = -1

# We don't use header() since we're manually printing a header before invoking the menu.
TerminalMenus.header(m::InspectMenu) = ""

function TerminalMenus.keypress(m::InspectMenu, key::UInt32)
    if key == Int(TerminalMenus.ARROW_RIGHT)
        m.scroll_horizontal -= 1
    elseif key == Int(TerminalMenus.ARROW_LEFT)
        m.scroll_horizontal += 1
    elseif key == Int('J')
        m.selected_command = SelectionOptions.JUMP
        return true
    elseif key == UInt32('v')
        m.selected_command = SelectionOptions.TOGGLE_VALUE
        return true
    end
    return false
end

function TerminalMenus.pick(menu::InspectMenu, cursor::Int)
    menu.selected = cursor
    return true #break out of the menu
end

function TerminalMenus.writeLine(buf::IOBuffer, menu::InspectMenu, idx::Int, cursor::Bool)
    term_width = REPL.Terminals.width(TerminalMenus.terminal)

    cursor_len = length(TerminalMenus.CONFIG[:cursor])
    # print a ">" on the selected entry
    cursor ? print(buf, TerminalMenus.CONFIG[:cursor]) : print(buf, repeat(" ", cursor_len))
    print(buf, " ") # Space between cursor and text

    line = replace(menu.options[idx], "\n" => "\\n")
    line,menu.scroll_horizontal = _custom_trimWidth(line, term_width, cursor, cursor_len, menu.scroll_horizontal)

    print(buf, line)
end

function _custom_trimWidth(str::String, term_width::Int, highlighted=false, pad::Int=2, scroll::Int=0)
    max_str_len = term_width - pad - 5
    str_len = length(str)
    if highlighted
        scroll = min(max(scroll, 0), min(str_len, str_len-max_str_len-pad))
    end
    if str_len <= max_str_len || str_len < 6
        return str, scroll
    end
    if !highlighted
        return string(str[1:max_str_len], ".."), scroll
    else
        if scroll > 0
            return string(str[max(1,end-max_str_len-scroll):max(max_str_len,end-scroll)], ".."), scroll
        else
            return string("..", str[length(str) - max_str_len-scroll:end]), scroll
        end
    end
end

end


text:
#
# Plural constants.
#

CollectionQuery(idom, odom, set::AbstractVector) =
    Query(
        CollectionSig(set),
        Input(convert(Domain, idom)),
        Output(convert(Domain, odom)) |> setoptional() |> setplural())

CollectionQuery(odom, set::AbstractVector) =
    CollectionQuery(Void, odom, set)

CollectionQuery{T}(set::AbstractVector{T}) =
    CollectionQuery(Void, T, set)

immutable CollectionSig <: AbstractPrimitive
    set::AbstractVector
end

ev(sig::CollectionSig, ::Input, oty::Output, iflow::InputFlow) =
    OutputFlow(oty, collection_impl(sig.set, length(iflow)))

function collection_impl(set::AbstractVector, len::Int)
    card = length(set)
    if len == 1
        return PluralColumn(Int[1, card+1], set)
    end
    offs = Int[1 + i*card for i = 0:len]
    idxs = len > 0 ? Int[j for i = 1:len for j = 1:card] : Int[]
    vals = set[idxs]
    return PluralColumn(offs, vals)
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
    InitialValueProblem <: AbstractInitialValueProblem

A composite type for an [`AbstractInitialValueProblem`](@ref).

# Constructors
```julia
InitialValueProblem(rhs, u0, tspan)
InitialValueProblem(rhs, u0, t0, tN)
IVP(args...; kwargs...)
```

# Arguments
- `rhs :: Union{Function, AbstractRightHandSideFunction}` : right-hand side derivative.
- `u0 :: Union{Number, AbstractVector{Number}}` : initial condition.
- `tspan :: Tuple{Real, Real}` : time domain.
- `t0 :: Real` : initial time.
- `tN :: Real` : final time.

# Functions
- [`show`](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
mutable struct InitialValueProblem{rhs_T, u0_T, tspan_T} <: AbstractInitialValueProblem
    rhs::rhs_T
    u0::u0_T
    tspan::tspan_T
end

InitialValueProblem(rhs::AbstractRightHandSideFunction, u0::Number, tspan) = InitialValueProblem(rhs, [u0], tspan)
InitialValueProblem(f::Function, u0, tspan) = InitialValueProblem(RHS(f), u0, tspan)
InitialValueProblem(rhs::AbstractRightHandSideFunction, u0, t0::Real, tN::Real) = InitialValueProblem(rhs, u0, (t0, tN))
InitialValueProblem(f::Function, u0, t0, tN) = InitialValueProblem(RHS(f), u0, t0, tN)
@doc (@doc InitialValueProblem) IVP(args...; kwargs...) = InitialValueProblem(args...; kwargs...)

############################################################################################
#                                         PRINTING                                         #
############################################################################################

"""
    show(io::IO, problem::InitialValueProblem)

prints a full description of `problem` and its contents to a stream `io`.
"""
Base.show(io::IO, problem::InitialValueProblem) = _show(io, problem)

"""
    summary(io::IO, problem::InitialValueProblem)

prints a brief description of `problem` to a stream `io`.
"""
Base.summary(io::IO, problem::InitialValueProblem) = _summary(io, problem)


text:
# # ==================== born forward ===============
export add_virtual_source,
       sourceside_reconstruct_forward,
       born_approximation_forward!,
# # ==================== born adjoint ===============
       apply_image_condition!,
       sourceside_reconstruct_backward,
       born_approximation_adjoint,
# # ==================== born =======================
       born_approximation_forward,
       born_approximation_adjoint,
       born_approximation,
       recordings_axpby!,
       image_axpby!,
       recordings_norm,
# # ==================== get reflections ============
       get_shotgather_adaptive,
       get_shotgather,
       get_reflections,
       get_wavefield_bound,
       initialize_lsrtm,
# # ==================== filters for image ==========
       laplace_filter,
       trace_interpolation

include("born/forward.jl")
include("born/adjoint.jl")
include("born/born.jl")
include("get_recordings.jl")
include("filter.jl")
include("upsampling.jl")


text:
function getNodeProperty(G::AbstractArray{Bool,3})
	Nx,Ny,Nz = size(G,1),size(G,2),size(G,3)
	I = Array{UInt8}(Nx-2,Ny-2,Nz-2)
	for l = 2:Nx-1, m = 2:Ny-1, i = 2:Nz-1
		if G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		   G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		   I[l-1,m-1,i-1] = 1

			#walls
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 2 #wall l = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 3 #wall l = Nx
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 0 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 4 #wall m = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 0 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 5 #wall m = Ny
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 6 #wall i = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 7 #wall i = Nz

			#edges
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 0 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 8 #edge l = 1, m = 1
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 0 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 9 #edge l = 1, m = Ny
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 1 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 10 #edge l = 1, i = 1
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 11 #edge l = 1, i = Nz

		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 0 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 12 #edge l = Nx, m = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 0 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 13 #edge l = Nx, m = Ny
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 14 #edge l = Nx, i = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 1 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 15 #edge l = Nx, i = Nz

		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 0 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 16 #edge m = 1, i = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 0 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 17 #edge m = 1, i = Nz
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 0 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 18 #edge m = Ny, i = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 0 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 19 #edge m = Ny, i = Nz

			#corners
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 0 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 20 #edge l = 1, m = 1, i = 1
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 0 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 1 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 21 #edge l = 1, m = 1, i = Nz
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 1 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 0 && G[l,m,i+1] == 1
		       I[l-1,m-1,i-1] = 22 #edge l = 1, m = Ny, i = 1
		elseif G[l-1,m,i] == 0 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 1 && G[l,m+1,i] == 0 && G[l,m,i+1] == 0
		       I[l-1,m-1,i-1] = 23 #edge l = 1, m = Ny, i = Nz

		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 0 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 1 && G[l,m,i+1] == 1
			I[l-1,m-1,i-1] = 24 #edge l = Nx, m = 1, i = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 0 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 1 && G[l,m,i+1] == 0
			I[l-1,m-1,i-1] = 25 #edge l = Nx, m = 1, i = Nz
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 0 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 0 && G[l,m,i+1] == 1
			I[l-1,m-1,i-1] = 26 #edge l = Nx, m = Ny, i = 1
		elseif G[l-1,m,i] == 1 && G[l,m-1,i] == 1 && G[l,m,i-1] == 1 &&
		       G[l+1,m,i] == 0 && G[l,m+1,i] == 0 && G[l,m,i+1] == 0
			I[l-1,m-1,i-1] = 27 #edge l = Nx, m = Ny, i = Nz
		end
	end
	return I
end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
using CSFML.LibCSFML

mode = sfVideoMode(1280, 720, 32)

window = sfRenderWindow_create(mode, "SFML window", sfResize | sfClose, C_NULL)
@assert window != C_NULL

texture = sfTexture_createFromFile(joinpath(@__DIR__, "julia-tan.png"), C_NULL)
@assert texture != C_NULL

sprite = sfSprite_create()
sfSprite_setTexture(sprite, texture, sfTrue)

font = sfFont_createFromFile(joinpath(@__DIR__, "Roboto-Bold.ttf"))
@assert font != C_NULL

text = sfText_create()
sfText_setString(text, "Hello SFML")
sfText_setFont(text, font)
sfText_setCharacterSize(text, 50)

music = sfMusic_createFromFile(joinpath(@__DIR__, "Chrono_Trigger.ogg"))
@assert music != C_NULL

sfMusic_play(music)

try
    event_ref = Ref{sfEvent}()
    while Bool(sfRenderWindow_isOpen(window))
        # process events
        while Bool(sfRenderWindow_pollEvent(window, event_ref))
            event_ptr = Base.unsafe_convert(Ptr{sfEvent}, event_ref)
            # close window : exit
            GC.@preserve event_ref begin
                ty = unsafe_load(event_ptr.type)
                ty == sfEvtClosed && sfRenderWindow_close(window)
                ty == sfEvtResized && println("Trigger sfEvtResized.")
                ty == sfEvtLostFocus && println("Trigger sfEvtLostFocus.")
                ty == sfEvtGainedFocus && println("Trigger sfEvtGainedFocus.")
                ty == sfEvtTextEntered && println("Trigger sfEvtTextEntered: $(unsafe_load(event_ptr.text).unicode)")
                ty == sfEvtKeyPressed && println("Trigger sfEvtKeyPressed: $(unsafe_load(event_ptr.key).code)")
                ty == sfEvtKeyReleased && println("Trigger sfEvtKeyReleased: $(unsafe_load(event_ptr.key).code)")
                ty == sfEvtMouseWheelMoved && println("Trigger sfEvtMouseWheelMoved: $(unsafe_load(event_ptr.mouseWheel).x), $(unsafe_load(event_ptr.mouseWheel).y)")
                ty == sfEvtMouseWheelScrolled && println("Trigger sfEvtMouseWheelScrolled: $(unsafe_load(event_ptr.mouseWheelScroll).wheel)")
                ty == sfEvtMouseButtonPressed && println("Trigger sfEvtMouseButtonPressed: $(unsafe_load(event_ptr.mouseButton).button)")
                ty == sfEvtMouseButtonReleased && println("Trigger sfEvtMouseButtonReleased: $(unsafe_load(event_ptr.mouseButton).x), $(unsafe_load(event_ptr.mouseButton).y)")
                ty == sfEvtMouseMoved && println("Trigger sfEvtMouseMoved: $(unsafe_load(event_ptr.mouseMove).x), $(unsafe_load(event_ptr.mouseMove).y)")
            end
        end
        # clear the screen
        sfRenderWindow_clear(window, sfColor_fromRGBA(0,0,0,1))
        # draw the sprite
        sfRenderWindow_drawSprite(window, sprite, C_NULL)
        # draw the text
        sfRenderWindow_drawText(window, text, C_NULL)
        # update the window
        sfRenderWindow_display(window)
    end
catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    sfMusic_destroy(music)
    sfText_destroy(text)
    sfFont_destroy(font)
    sfSprite_destroy(sprite)
    sfTexture_destroy(texture)
    sfRenderWindow_destroy(window)
end


text:
#####################################################
#Weight update functions

function backprop!(net::AE, learningrate, weight_decay;
				   nonlinearity_diff = Array{Function, 1},
					 BPuptolayer = net.nl)
	for i in net.nl:-1:2
		nonlinearity_diff[i](net.ax[i], net.e[i])
		BLAS.gemv!('T', 1., net.w[i], net.e[i], 0., net.e[i-1])
	end
	nonlinearity_diff[1](net.ax[1], net.e[1])
	for i in 1:BPuptolayer
		BLAS.ger!(learningrate[i], net.e[i], net.x[i], net.w[i])
		BLAS.axpy!(learningrate[i], net.e[i], net.b[i])
		if weight_decay != 0.
			scale!(net.w[i],1. - weight_decay)
			scale!(net.b[i],1. - weight_decay)
		end
	end
end

function _normalize_inputweights!(weights)
	for j in 1:size(weights)[1]
		weights[j,:] .*= 1. / norm(weights[j,:])
	end
end

#Algorithm for sparse coding as proposed in:
#Brito CSN, GerstnerW(2016) PLoS Comput Biol 12(9).
#Implements the second part of algorithm: weight update
#omitnonlocal: boolian if non local term should be omitted:
#true: will be omitted
#false: will be taken into account
function _backprop_sparse_Brito!(net::AE_sparse, learningrate; omitnonlocal = true, patches = false, w_connect = [])
	if net.nl > 1
		error("Network/AE is deeper than expected for a sparse network (more than one hidden layer)!\n")
	end
	if omitnonlocal
		for i in 1:net.nl
			#BLAS.ger!(learningrate[i],net.ax[i],net.x[i],net.w[i])
			#Above: might be useful for ordinary sparse coding à la mairal which uses .ax and .x differently!
			BLAS.ger!(learningrate[i],net.x[i+1],net.x[i],net.w[i])
			#BLAS.ger!(learningrate[i],net.x[i+1].*(net.x[i+1]-0.1),net.x[i],net.w[i]) #quadratic Hebbian
			_normalize_inputweights!(net.w[i])
	  end
	else
		for i in 1:net.nl
			#BLAS.ger!(learningrate[i],net.ax[i],(net.x[i]-BLAS.gemv('T',net.w[i],net.ax[i])),net.w[i])
			BLAS.ger!(learningrate[i],net.x[i+1],(net.x[i]-BLAS.gemv('T',net.w[i],net.x[i+1])),net.w[i])
			_normalize_inputweights!(net.w[i])
		end
	end
	if patches
		net.w[1] .*= w_connect
	end
end
function _backprop_sparse_Földiak!(net::AE_sparse;
		p = 0.1, #2. / length(net.ax[1]), #0.4
		learningrate_v = 1e-2, learningrate_w = 1e-3, learningrate_thr = 1e-2,
		patches = false, w_connect = [])

	#BLAS.ger!(learningrate_v, Float64.(net.x[2] .> net.b[1]), Float64.(net.x[2] .> net.b[1]), net.A[1])
	BLAS.ger!(learningrate_v,net.x[2],net.x[2],net.A[1])
	net.A[1] .+= - learningrate_v * p^2
	for j in 1:size(net.A[1])[1]
		net.A[1][j,j] = 0. #no self-inhibition
	end
	clamp!(net.A[1],0.,Inf64) #Dale's law
	net.w[1] = Diagonal(1 .- learningrate_w * net.x[2]) * net.w[1]
	BLAS.ger!(learningrate_w,net.x[2],net.x[1],net.w[1])
	#BLAS.axpy!(learningrate_thr,Float64.(net.x[2] .> net.b[1]) .- p, net.b[1])
	BLAS.axpy!(learningrate_thr,net.x[2] .- p, net.b[1]) # ...,Int.(net.x[2] .> net.b[1]) .- p,...
	if patches
		net.w[1] .*= w_connect
	end
end
function backprop_sparse_Brito!(net::AE_sparse, learningrate; omitnonlocal = true, patches = false, w_connect = [], Földiak_model = false)
	if Földiak_model
		_backprop_sparse_Földiak!(net; patches = patches)
	else
		_backprop_sparse_Brito!(net, learningrate; omitnonlocal = omitnonlocal, patches = patches, w_connect = w_connect)
	end
end

function backprop_sparse_Brito!(net::AE_sparse_patchy, learningrate; omitnonlocal = true, Földiak_model = false)
	for AE_sparse_patch in net.AE_sparse_patches
		backprop_sparse_Brito!(AE_sparse_patch, learningrate;
			omitnonlocal = omitnonlocal, patches = false, Földiak_model = Földiak_model)
	end
end


text:
using LinearAlgebra

## Iterate Riccati difference equation
function Riccati(A_tld,B_tld,R,Q_tld,tol)

    P = ones(size(A_tld)[2],size(A_tld)[1]);
    F = ones(size(R)[1],size(W)[1]);

    pDiff = 10;
    fDiff = 10;
    iter = 0;

    while pDiff >= tol.*opnorm(P,1) && fDiff > tol.*opnorm(F,1)
        A = R+B_tld'*P*B_tld;

        P1 = Q_tld+A_tld'*P*A_tld-A_tld'*P*B_tld*(A\B_tld')*P*A_tld;
        pDiff = opnorm(P1-P,1);

        F1 = A\B_tld'*P*A_tld;
        fDiff = opnorm(F1-F,1);

        iter = iter + 1;

        P = copy(P1)
        F = copy(F1)
    end

    return P,F

end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function drawtriangle(points, degree)
    sethue(cols[degree])
    poly(points, :fill)
end

function sierpinskitriangle(points, degree)
    drawtriangle(points, degree)
    if degree > 1
        p1, p2, p3 = points
        sierpinskitriangle([p1, midpoint(p1, p2),
                                midpoint(p1, p3)], degree-1)
        sierpinskitriangle([p2, midpoint(p1, p2),
                                midpoint(p2, p3)], degree-1)
        sierpinskitriangle([p3, midpoint(p3, p2),
                                midpoint(p1, p3)], degree-1)
    end
end

function drawsierpinski(fname, n)
  Drawing(200, 200, fname)
  origin()
  background("ivory")
  circle(O, 75, :clip)
  my_points = ngon(O, 150, 3, -pi/2, vertices=true)
  depth = 8 # 12 is ok, 20 is right out
  sierpinskitriangle(my_points, n)
  @test finish() == true
end

cols = distinguishable_colors(8)
fname = "sierpinskitriangle.svg"
drawsierpinski(fname, 8)
println("...finished test: output in $(fname)")


text:
# TSP instance for using BRKGA using tsplib
module TSP

using Distances, BRKGA, Logging

export TSPProblem
export readTSP, decoder
export npairs, pair, param, nalleles

type TSPProblem <: BRKGAProblem
    pairs::Vector{Vector{Float64}}
    param::Parameters
end

# Funtions for TSPProblem
npairs(tsp::TSPProblem)         = length(tsp.pairs)
pair(tsp::TSPProblem, idx::Int) = tsp.pairs[idx]
BRKGA.nalleles(tsp::TSPProblem) = npairs(tsp)

"""

Read the TSPlib file
"""
function readTSP(file::String)
    pairs = Vector{Vector{Float64}}
    open(file,"r") do f
        readline(f) # NAME:
        readline(f) # COMMENT:
        readline(f) # TYPE:

        line = readline(f) # DIMENSION:
        items = split(line, " ")
        dim = parse(Int, items[2])
        readline(f) # EDGE_WEIGHT_TYPE
        readline(f) #NODE_COORD_SECTION

        pairs = Array(Vector{Float64}, dim)
        for i = 1:dim
            line = readline(f)
            items = matchall(r"[0-9]+", line)
            x = parse(Int64, items[end-1])
            y = parse(Int64, items[end])
            pairs[i] = [x, y]
        end
    end
    return pairs
end

"""

Evaluate distance between two nodes
"""
function distance(p1::Vector{Float64}, p2::Vector{Float64})
    return floor(Int64, euclidean(p1, p2)+0.5)
end

"""

Decoder for TSP
"""
function BRKGA.decoder(tsp::TSPProblem, chromosome::Vector{Float64})

    idx = sortperm(chromosome)

    # Evaluate the tour distance
    dist = 0
    for i = 2:npairs(tsp)
        p1 = pair(tsp, idx[i-1])
        p2 = pair(tsp, idx[i])
        dist += distance(p1, p2)
    end
    firt = pair(tsp, idx[1])
    last = pair(tsp, idx[end])
    dist += distance(firt, last)
    return dist
end

end #TSP


text:
# fields.jl

using CUDA
using OffsetArrays
using Adapt

Adapt.adapt_structure(to, x::OffsetArray) = OffsetArray(adapt(to, parent(x)), x.offsets)

"""Initialize a 3D array with the correct dimensions, possible with boundary offset, possibly as CUDA array."""
function initField(inival, nx, ny, nz, useCUDA, useOffset, myFloat)
    # Initialize a field/array

    # Convert inival to correct Float type
    inival = myFloat(inival)

    # If it is an offset array, increase size by 2
    if useOffset
        var = fill(inival, (nx + 2, ny + 2, nz + 2))
    else
        var = fill(inival, (nx, ny, nz))
    end

    # CUDA array y/n
    if useCUDA
        # Convert to CUDA Arrays
        var = CuArray(var)
    end

    # When offset array, shift start index to 0
    if useOffset
        var = OffsetArray(var, (0:nx + 1, 0:ny + 1, 0:nz + 1))
    end

    return var
end

# fields.jl

using CUDA
using OffsetArrays

"""Initialize a 1D vector with the correct dimensions, possible with boundary offset, possibly as CUDA array."""
function initField1D(inival, nx, useCUDA, useOffset, myFloat)
    # Initialize a field/array

    # Convert inival to correct Float type
    inival = myFloat(inival)

    # If it is an offset array, increase size by 2
    if useOffset
        var = fill(inival, nx + 2)
    else
        var = fill(inival, nx)
    end

    # CUDA array y/n
    if useCUDA
        # Convert to CUDA Arrays
        var = CuArray(var)
    end

    # When offset array, shift start index to 0
    if useOffset
        var = OffsetArray(var, 0:nx + 1)
    end

    return var
end


text:
julia> arr = [1, 2, nothing, 3]
4-element Array{Union{Nothing, Int64},1}:
 1
 2
  nothing
 3

julia> x = arr .+ 5
ERROR: MethodError: no method matching +(::Nothing, ::Int64)
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:502
  +(::Complex{Bool}, ::Real) at complex.jl:292
  +(::Missing, ::Number) at missing.jl:93
  ...

julia> arr = [1, 2, missing, 3]
4-element Array{Union{Missing, Int64},1}:
 1
 2
  missing
 3

julia> x = arr .+ 5
4-element Array{Union{Missing, Int64},1}:
 6
 7
  missing
 8

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestFeatureSelectors

using Random
using Test
using CSV
using AutoMLPipeline
using AutoMLPipeline.DecisionTreeLearners
using AutoMLPipeline.AbsTypes
using AutoMLPipeline.SKLearners
using AutoMLPipeline.SKPreprocessors
using AutoMLPipeline.Utils
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.Pipelines
using AutoMLPipeline.EnsembleMethods
using AutoMLPipeline.FeatureSelectors
using Statistics
using DataFrames

function iris_test()
    data = getiris()
    X = data[:,1:5]
    catfeat = FeatureSelector([5])
    numfeat = FeatureSelector([1,2,3,4])
    autocat = CatFeatureSelector()
    autonum = NumFeatureSelector()
    @test (fit_transform!(catfeat,X) .== X[:,5]) |> Matrix |> sum == 150
    @test (fit_transform!(numfeat,X) .== X[:,1:4]) |> Matrix |> sum == 600
    @test (fit_transform!(autocat,X) .== X[:,5]) |> Matrix |> sum == 150
    @test (fit_transform!(autonum,X) .== X[:,1:4]) |> Matrix |> sum == 600
    catnumdata = hcat(X,repeat([1,2,3,4,5],30))
    catnum = CatNumDiscriminator()
    res = fit_transform!(catnum,catnumdata)
    @test eltype(catnumdata[:,6]) <: Number
    @test eltype(catnumdata[:,2]) <: Number
    @test eltype(catnumdata[:,4]) <: Number
    @test eltype(res[:,6]) <: String
    @test eltype(res[:,2]) <: String
    @test eltype(res[:,4]) <: String
    catnumdata = hcat(X,repeat([1,2,3,4,5],30))
    catnum = CatNumDiscriminator(0)
    res = fit_transform!(catnum,catnumdata)
    @test eltype(res[:,6]) <: Number
    catnumdata = hcat(X,repeat([1,2,3,4,5],30))
    catnum = CatNumDiscriminator(5)
    ppp=@pipeline catnum |>  ((autocat |> OneHotEncoder()) + (autonum |> SKPreprocessor("PCA")))
    res=fit_transform!(ppp,catnumdata)
    @test ncol(res) == 12
end
@testset "Feature Selectors: Iris" begin
    Random.seed!(123)
    iris_test()
end

function diabetes_test()
    Random.seed!(123)
    diabetesdf = CSV.read(joinpath(dirname(pathof(AutoMLPipeline)),"../data/diabetes.csv"))
    X = diabetesdf[:,1:end-1]
    Y = diabetesdf[:,end] |> Vector
    pca = SKPreprocessor("PCA")
    ica = SKPreprocessor("FastICA")
    dt = SKLearner("DecisionTreeClassifier")
    rf = SKLearner("RandomForestClassifier")
    rbs = SKPreprocessor("RobustScaler")
    jrf = RandomForest()
    lsvc = SKLearner("LinearSVC")
    ohe = OneHotEncoder()
    catf = CatFeatureSelector()
    numf = NumFeatureSelector()

    disc = CatNumDiscriminator(0)
    pl = @pipeline disc |> ((numf |>  pca) + (catf |> ohe)) |> jrf
    @test crossvalidate(pl,X,Y,"accuracy_score",10,false).mean > 0.60

    pl = @pipeline disc |> ((numf |> rbs |>  pca) + (catf |> ohe)) |> lsvc
    @test crossvalidate(pl,X,Y,"accuracy_score",10,false).mean > 0.60

    pl = @pipeline disc |> ((numf |> rbs |>  ica) + (catf |> ohe)) |> rf
    @test crossvalidate(pl,X,Y,"accuracy_score",10,false).mean > 0.60

    disc = CatNumDiscriminator(20)
    pl = @pipeline disc |> ( (catf |> ohe)) |> jrf
    @test crossvalidate(pl,X,Y,"accuracy_score",20,false).mean > 0.60

    disc = CatNumDiscriminator(50)
    pl = @pipeline disc |> ((numf |> rbs |>  pca) + (catf |> ohe)) |> lsvc
    @test crossvalidate(pl,X,Y,"accuracy_score",20,false).mean > 0.60

    disc = CatNumDiscriminator(100)
    pl = @pipeline disc |> ((numf |> rbs |>  ica) + (catf |> ohe)) |> rf
    @test crossvalidate(pl,X,Y,"accuracy_score",20,false).mean > 0.60
end
@testset "Feature Selectors: Diabetes" begin
    Random.seed!(123)
    diabetes_test()
end

end


text:
# communications
function hopen(spec::AbstractString)
    parts = split(spec, ":"; limit=3)
    n = length(parts)
    if n == 1
        port = parse(Int, parts[1])
        return hopen(port)
    end
    host = parts[1]
    port = parse(Int, parts[2])
    if n == 2
        return hopen(host, port)
    else
        user = parts[3]
        return hopen(host, port, user)
    end
end
hopen(host::String, port::Integer, user::String, timeout::Integer) =
    khpun(host, port, user, timeout)
hopen(host::String, port::Integer, user::String) = khpu(host, port, user)
hopen(host::String, port::Integer) = khp(host, port)
hopen(port::Integer) = hopen("localhost", port)
function hopen(;host="localhost", port=-1, user="", timeout=-1)
    if port == -1
        throw(ArgumentError("A port value must be specified"))
    end
    if timeout == -1
        if user == ""
            return khp(host, port)
        else
            return khpu(host, port, user)
        end
    else
        return khpun(host, port, user, timeout)
    end
end

function hopen(f::Function; host="localhost", port=-1, user="", timeout=-1)
    h = hopen(;host=host, port=port, user=user, timeout=timeout)
    try
        return f(h)
    finally
        kclose(h)
    end
end

function hopen(f::Function, host::String, port::Integer)
    h = hopen(host, port)
    try
        return f(h)
    finally
        kclose(h)
    end
end

hopen(f::Function, port::Integer) = hopen(f, "", port)

const hclose = kclose

hget(h::Integer, m::String) = _E(k(h, m))
hget(h::Integer, m::String, x...) = _E(k(h, m, map(K_new, x)...))

function hget(h::Tuple{String,Integer}, m)
   h = hopen(h...)
   try
       return hget(h, m)
   finally
       kclose(h)
   end
end

function hget(h::Tuple{String,Integer}, m, x...)
   h = hopen(h...)
   try
       return hget(h, m, x...)
   finally
       kclose(h)
   end
end

function open_default_kdb_handle(msg_io=STDERR)
    if KDB_HANDLE[] > 0
        return KDB_HANDLE[]
    end
    server_spec = get(ENV, "KDB", "")
    if isempty(server_spec)
        port, process = Q.Kdb.start()
        handle = hopen(port)
        atexit() do
            rc = Q.Kdb.stop(handle, process)
            info("Slave kdb+ exited with code $rc.")
        end
        info(msg_io, "Connected to $port.")
    else
        try
            handle = hopen(server_spec)
            info(msg_io, "Connected to $server_spec.")
        catch error
            warn(msg_io,
                 "Could not connect to $server_spec. $error")
            handle = -1
        end
    end
    KDB_HANDLE[] = handle
end
# Client-side q`..`
struct _Q
end

export @q_cmd
struct Qcmd
    cmd::String
end

macro q_cmd(cmd) Qcmd(cmd) end

Base.show(io::IO, x::Qcmd) = print(io, "q`", x.cmd, "`")

function (x::Qcmd)(args...)
    K(k(KDB_HANDLE[], x.cmd, map(K_new, args)...))
end
# Initialise memory without making a connection
khp("", -1)


text:
module OneClassActiveLearning

using Dates
using DelimitedFiles
using Memento
using Printf
using Random
using Reexport
using Statistics
using SVDD
using DataStructures

import DataStructures
import Formatting
import JSON
import JuMP
import MLDataUtils
import MLKernels
import MLLabelUtils
import ROCAnalysis
import Unmarshal
import ValueHistories

import Base.show
import StatsBase: countmap

include("data_util.jl")
include("evaluate.jl")

include("QueryStrategies/QueryStrategies.jl")
@reexport using .QueryStrategies

include("Oracles/Oracles.jl")
@reexport using .Oracles

include("result.jl")
include("serialize.jl")
include("run.jl")
include("log_util.jl")

const LOGGER = getlogger(@__MODULE__)

function __init__()
    Memento.register(LOGGER)
end

export
    load_data,
    SplitStrategy,
    FullSplitStrat, UnlabeledSplitStrat, UnlabeledAndLabeledInlierSplitStrat,
    LabeledSplitStrat, LabeledInlierSplitStrat, LabeledOutlierSplitStrat,
    DataSplits,
    get_train, get_test, get_query, calc_mask,
    get_splits_and_init_pools, get_initial_pools,

    ConfusionMatrix,
    cohens_kappa,
    matthews_corr,
    roc_auc,
    get_positive,
    get_negative,
    tpr,
    fpr,
    tnr,
    sensitivity,
    specificity,
    recall,
    get_n,
    f1_score,

    DataStats, Result
end


text:
# TODO: add more target distributions and make them iteratable

# Dimension of testing distribution
const D = 5

# Deterministic tolerance
const DETATOL = 1e-3 * D
# Random tolerance
const RNDATOL = 5e-2 * D

using Distributions: logpdf, MvNormal
using ForwardDiff: gradient

logπ(θ::AbstractVector{T}) where {T<:Real} = logpdf(MvNormal(zeros(D), ones(D)), θ)
∂logπ∂θ = θ -> gradient(logπ, θ)


text:
module Forcings

export Forcing, ContinuousForcing, DiscreteForcing, Relaxation, GaussianMask, LinearTarget, AdvectiveForcing

using Oceananigans.Fields

include("continuous_forcing.jl")
include("discrete_forcing.jl")
include("relaxation.jl")
include("advective_forcing.jl")
include("forcing.jl")
include("model_forcing.jl")
include("multiple_forcings.jl")

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.3773826343471955 <= q <= 0.8255956851677421)

                     add_NL_constraint(m, :(log(1 + exp(-0.5817410001457377 + -0.38051045291676777*log(1 + exp(-0.7752766980375418 + 0.20887750352529233*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + -0.43558326218422394*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + -0.5061279411570236*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))) + -0.828517562717189*log(1 + exp(-0.31567568822764036 + 0.26872838410166544*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + 0.5533650825489542*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + 0.11778701597368801*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))) + -0.13996198025628548*log(1 + exp(0.6085123540494801 + -0.012952465394535206*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + 0.6228552648202648*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + 0.7276708970518393*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))))) + log(1 + exp(-0.08252400134510385 + 0.3442818429504264*log(1 + exp(-0.7752766980375418 + 0.20887750352529233*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + -0.43558326218422394*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + -0.5061279411570236*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))) + -0.7813134096358496*log(1 + exp(-0.31567568822764036 + 0.26872838410166544*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + 0.5533650825489542*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + 0.11778701597368801*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))) + -0.5888247957791601*log(1 + exp(0.6085123540494801 + -0.012952465394535206*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + 0.6228552648202648*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + 0.7276708970518393*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))))) + log(1 + exp(-0.2140493400207517 + 0.6401317985056014*log(1 + exp(-0.7752766980375418 + 0.20887750352529233*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + -0.43558326218422394*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + -0.5061279411570236*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))) + -0.4453350709432602*log(1 + exp(-0.31567568822764036 + 0.26872838410166544*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + 0.5533650825489542*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + 0.11778701597368801*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))) + 0.9148591417614278*log(1 + exp(0.6085123540494801 + -0.012952465394535206*log(1 + exp(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2]))) + 0.6228552648202648*log(1 + exp(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2]))) + 0.7276708970518393*log(1 + exp(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "ArrayIndices" begin
    @test ArrayIndices() isa ArrayIndices{Any}

    inds = ArrayIndices{Int64}()
    @test ArrayIndices{Int64, Vector{Int64}}() == inds

    @test isinsertable(inds)
    @test length(inds) == 0
    @test keys(inds) === inds
    @test isempty(inds)
    @test isequal(copy(inds), inds)
    @test_throws IndexError inds[10]
    @test length(unset!(inds, 10)) == 0
    io = IOBuffer(); print(io, inds); @test String(take!(io)) == "{}"
    if VERSION < v"1.6-"
        io = IOBuffer(); show(io, MIME"text/plain"(), inds); @test String(take!(io)) == "0-element ArrayIndices{Int64,Array{Int64,1}}"
    else
        io = IOBuffer(); show(io, MIME"text/plain"(), inds); @test String(take!(io)) == "0-element ArrayIndices{Int64, Vector{Int64}}"
    end
    @test_throws IndexError delete!(inds, 10)

    insert!(inds, 10)

    @test length(inds) == 1
    @test keys(inds) === inds
    @test !isempty(inds)
    @test isequal(copy(inds), inds)
    @test inds[10] == 10
    @test_throws IndexError insert!(inds, 10)
    @test length(set!(inds, 10)) == 1
    @test_throws IndexError insert!(inds, 10)
    io = IOBuffer(); print(io, inds); @test String(take!(io)) == "{10}"
    if VERSION < v"1.6-"
        io = IOBuffer(); show(io, MIME"text/plain"(), inds); @test String(take!(io)) == "1-element ArrayIndices{Int64,Array{Int64,1}}\n 10"
    else
        io = IOBuffer(); show(io, MIME"text/plain"(), inds); @test String(take!(io)) == "1-element ArrayIndices{Int64, Vector{Int64}}\n 10"
    end
    @test !isequal(inds, empty(inds))
    @test isequal(inds, copy(inds))
    @test isempty(empty(inds))

    delete!(inds, 10)

    @test isequal(inds, ArrayIndices{Int64}())
    @test union!(inds, Indices([1,2,3,4])) == Indices([1, 2, 3, 4])
    @test filter!(iseven, inds) == Indices([2, 4])

    empty!(inds)
    for i = 2:2:100
        insert!(inds, i)
    end
    @test issetequal(inds, Indices(2:2:100))
    @test all(in(i, inds) == iseven(i) for i in 2:100)
    @test isempty(empty!(inds))

    @testset "covert" begin
        i = Indices{Int32}([1,2,3])
        ai = ArrayIndices{Int32}([1,2,3])

        @test convert(AbstractIndices{Int32}, ai) === ai
        @test convert(AbstractIndices{Int64}, ai)::ArrayIndices{Int64} == i

        @test convert(ArrayIndices{Int32}, ai) === ai
        @test convert(ArrayIndices{Int64}, ai)::ArrayIndices{Int64} == i

        @test convert(ArrayIndices{Int32}, i)::ArrayIndices{Int32} == i
        @test convert(ArrayIndices{Int64}, i)::ArrayIndices{Int64} == i
    end

    @testset "rand" begin
        inds = ArrayIndices(["a", "b", "c", "d", "e"])
        for i = 1:100
            @test rand(inds) in ["a", "b", "c", "d", "e"]
        end
    end
end


text:
module ABCDistances

using StatsBase: Weights, cov, sample, mean
using Distributions
using ProgressMeter
import Base.show, Base.copy

export
  ABCDistance, Lp, Euclidean, Logdist, WeightedEuclidean, MahalanobisEmp, ##ABCDistance types
  init, evaldist, ##ABCDistance methods
  ABCInput, RefTable, ABCOutput, ABCRejOutput, ABCPMCOutput, ##General ABC types
  sortABCOutput!, show, copy, parameter_means, parameter_vars, parameter_covs, ##General ABC methods
  abcRejection, abcPMC2, abcPMC3, abcPMC3V, abcPMC4, abcPMC5, ##ABC algorithms
  rgk, rgk_os, ##g&k methods
  Stoichiometry, gillespie_partial_sim, gillespie_sim ##Gillespie algorithm

include("distances.jl")
include("abc.jl")
include("abc_rejection.jl")
include("abc_pmc.jl")
include("abc_pmc_interface.jl")
include("gk.jl")
include("gillespie.jl")
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# abstract types
abstract type AbstractUGF end
abstract type AbstractInfo end
abstract type AbstractDistribution{N,R} end
abstract type AbstractSTD{T} <: _LG.AbstractGraph{T} end
abstract type AbstractNetwork{T} <: _LG.AbstractGraph{T} end
abstract type AbstractStochasticProcess end
abstract type AbstractMarkovProcess <: AbstractStochasticProcess end

# union types
UIE = Union{Int,_LG.AbstractEdge}
Single = Union{Bool,Number,String,Symbol,AbstractSTD,Tuple}

# dict types
LibDict = Dict{UIE,Vector{Int}}
PropDict = Dict{Symbol,Any}

# broadcastable
Broadcast.broadcastable(dst::AbstractDistribution) = Ref(dst)


text:
using JLD2
@load "plotting_data/plotting_data_30.jld2" thirty
@load "plotting_data/plotting_data_60.jld2" sixty
@load "plotting_data/plotting_data_100.jld2" hundred

mat"
idx = 5
%(1:idx:end)
figure
hold on
subplot(3,1,1)
hold on
h(4)=plot($thirty.t_hist(1:idx:end),$thirty.a_hist(1:idx:end)/1000,'linewidth',2)
h(3)=plot($sixty.t_hist(1:idx:end),$sixty.a_hist(1:idx:end)/1000,'linewidth',2)
h(2)=plot($hundred.t_hist(1:idx:end),$hundred.a_hist(1:idx:end)/1000,'linewidth',2)
h(1) = plot([0 $hundred.t_hist(end)], 42e3*ones(2,1),'Color','#9e9e9e','linewidth',2);
legend([h(1)],'GEO SMA','Location','SouthEast')
ylabel('Semi-Major Axis (km)')
xlabel('Time (days)')
xlim([0,$hundred.t_hist(end)])
ylim([20e3,46e3])
subplot(3,1,2)
hold on
plot($thirty.t_hist(1:idx:end),$thirty.e_hist(1:idx:end),'linewidth',2)
plot($sixty.t_hist(1:idx:end),$sixty.e_hist(1:idx:end),'linewidth',2)
plot($hundred.t_hist(1:idx:end),$hundred.e_hist(1:idx:end),'linewidth',2)
ylabel('Eccentricity')
xlabel('Time (days)')
xlim([0,$hundred.t_hist(end)])
subplot(3,1,3)
hold on
plot($thirty.t_hist(1:idx:end),$thirty.i_hist(1:idx:end),'linewidth',2)
plot($sixty.t_hist(1:idx:end),$sixty.i_hist(1:idx:end),'linewidth',2)
plot($hundred.t_hist(1:idx:end),$hundred.i_hist(1:idx:end),'linewidth',2)
legend('30 day','60 day','100 day','location','northeast')
ylabel('Inclination (deg)')
xlabel('Time (days)')
xlim([0,$hundred.t_hist(end)])
hold off
%addpath('matlab2tikz')
addpath('/Users/kevintracy/devel/Low-Thrust-TrajOpt/matlab2tikz/src')
disp('pre')
%matlab2tikz('mytest.tikz')
disp('post')
"

# mat"
# figure
# hold on
# subplot(2,2,1)
# hold on
# plot($thirty.t_hist,$thirty.a_hist/1000)
# plot($sixty.t_hist,$sixty.a_hist/1000)
# plot($hundred.t_hist,$hundred.a_hist/1000)
# h(1) = plot($hundred.t_hist, 42e3*ones(length($hundred.t_hist),1),'k--')
# legend([h(1)],'GEO SMA','Location','SouthEast')
# ylabel('Semi-Major Axis (km)')
# xlabel('Time (days)')
# xlim([0,$hundred.t_hist(end)])
# ylim([20e3,46e3])
# subplot(2,2,2)
# hold on
# plot($thirty.t_hist,$thirty.e_hist)
# plot($sixty.t_hist,$sixty.e_hist)
# plot($hundred.t_hist,$hundred.e_hist)
# ylabel('Eccentricity')
# xlabel('Time (days)')
# xlim([0,$hundred.t_hist(end)])
# subplot(2,2,3)
# hold on
# h1=plot($thirty.t_hist,$thirty.i_hist)
# h2=plot($sixty.t_hist,$sixty.i_hist)
# h3=plot($hundred.t_hist,$hundred.i_hist)
# %legend('30 day','60 day','100 day','location','southoutside')
# ylabel('Inclination (deg)')
# xlabel('Time (days)')
# xlim([0,$hundred.t_hist(end)])
# hold off
#
# hL = subplot(2,2,4);
# poshL = get(hL,'position');     % Getting its position
#
# lgd = legend(hL,[h1;h2;h3],'RandomPlot1','RandomPlot2','RandomPlot3');
# set(lgd,'position',poshL);      % Adjusting legend's position
# axis(hL,'off');
# set(hL,'position',[0 0 0.1 0.1])
# saveas(gcf,'test.png')
# "

mat"
close all
close all hidden
idx = 3
st = 3
figure
hold on
plot($thirty.t_hist(st:idx:end),.1*$thirty.Unorm(st:idx:end),'linewidth',1.5)
plot($sixty.t_hist(st:idx:end),.1*$sixty.Unorm(st:idx:end),'linewidth',1.5)
plot($hundred.t_hist(st:idx:end),.1*$hundred.Unorm(st:idx:end),'linewidth',1.5)
legend('30-Day Transfer','60-Day Transfer','100-Day Transfer')
xlabel('Time (days)')
ylabel('Thrust (N)')
addpath('/Users/kevintracy/devel/Low-Thrust-TrajOpt/matlab2tikz/src')
disp('pre')
%matlab2tikz('unorm.tex')
disp('post')
hold off
"

# plot the thrust angles for all trajectories

mat"
figure
hold on

subplot(3,1,1)
plot($thirty.t_hist,rad2deg($thirty.angles'))
xlabel('Time (days)')
ylabel('Thrust Angle (deg)')
legend('In-Plane','Out-of-Plane','Location','SouthWest')
ylim([-180 180])
yticks([-180 0 180])
xlim([0 $thirty.t_hist(end)])
hold off

subplot(3,1,2)
plot($sixty.t_hist,rad2deg($sixty.angles'))
xlabel('Time (days)')
ylabel('Thrust Angle (deg)')
legend('In-Plane','Out-of-Plane','Location','SouthWest')
ylim([-180 180])
yticks([-180 0 180])
xlim([0 $sixty.t_hist(end)])
hold off

subplot(3,1,3)
plot($hundred.t_hist,rad2deg($hundred.angles'))
xlabel('Time (days)')
ylabel('Thrust Angle (deg)')
legend('In-Plane','Out-of-Plane','Location','SouthWest')
ylim([-180 180])
yticks([-180 0 180])
xlim([0 $hundred.t_hist(end)])
hold off

"

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
struct FDHeatProblem{TBCL, TBCR,TTS}
    a # (a,b) define the domain
    b
    Δx # mesh spacing
    nx # number of points in the interior, excluding a and b
    x # mesh on which we solve the heat equation
    Δt # time step
    nΔt # number of time steps
    tmax # final time
    u0 # Initial condition
    A # the matrix (default to sparse)
    rhs # right hand side of our system
    left_bc::TBCL
    right_bc::TBCR
    time_step::TTS
end

struct DirichletBC
    u_bc
end

struct NeumannBC
    v_bc
end

struct PeriodicBC
end

struct BEuler
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#display lenslet systems: emitter, paraxial lenslet, paraxial eye lens, retina detector.
#start with eyebox plane and hexagonal tiling. Project onto display surface, generate lenslets automatically.

using OpticSim.Geometry:Transform,world2local
using OpticSim:plane_from_points,surfaceintersection,closestintersection,Ray,Plane,ConvexPolygon,Sphere
using OpticSim
using OpticSim.Repeat:tilevertices,HexBasis1,tilesinside
using Unitful:upreferred
using Unitful.DefaultSymbols:°

# using OpticSim.Vis:drawcells

"""compute the mean of the columns of `a`. If `a` is an `SMatrix` this is very fast and will not allocate."""
columncentroid(a::AbstractMatrix) = sum(eachcol(a))/size(a)[2] #works except for the case of zero dimensional matrix.
export columncentroid

function project(point::AbstractVector{T},projectionvector::AbstractVector{T},surface::OpticSim.Surface{T}) where{T}
    ray = Ray(point, projectionvector)
    intsct = surfaceintersection(surface, ray)
    pointintsct = closestintersection(intsct,false)
    if pointintsct === nothing #point didn't project to a point on the surface
        return nothing
    else
        return OpticSim.point(pointintsct)
    end
end
export project

"""project the vertices of a polygon represented by `vertices` onto `surface` using the point as the origin and `projectionvector` as the projection direction. Return nothing if any of the projected points do not intersect the surface. The projected vertices are not guaranteed to be coplanar."""
function project(vertices::AbstractMatrix{T}, projectionvector::AbstractVector{T}, surface::OpticSim.Surface{T}) where {T}
    result = similar(vertices)

    for i in 1:size(vertices)[2]
        origin = vertices[:, i]
        pt = project(origin,projectionvector,surface)

        if pt === nothing #one of the polygon points didn't project to a point on the surface so reject the polygon
            return nothing
        else
            result[:, i] = pt
        end
    end
    return result
end

"""Finds the best fit plane to `vertices` then projects `vertices` onto this plane by transforming from the global to the local coordinate frame. The projected points are represented in the local coordinate frame of the plane."""
function projectonbestfitplane(vertices::AbstractMatrix{T}) where{T}
    @assert size(vertices)[1] == 3 "projection only works for 3D points"

    center, _, localrotation  = plane_from_points(vertices)
    toworld = Transform(localrotation,center) #compute local to world transformation
    tolocal = world2local(toworld)
    numpts = size(vertices)[2]
    result = tolocal * SMatrix{3,numpts}(vertices)

    return result,toworld,tolocal #project onto best fit plane by setting z coordinate to zero.
end
export projectonbestfitplane

function testproject()
    normal = SVector(0.0, 0, 1)
    hex = HexBasis1()
    verts =vcat(tilevertices((0, 0), hex), [0.0 0 0 0 0 0])
    surf = Plane(normal, SVector(0.0, 0, 10))

    project(verts, normal, surf)
end
export testproject

function testprojectonplane()
    verts = tilevertices(HexBasis1())
    verts = vcat(verts,[0 for _ in 1:6]')
    projectonbestfitplane(SMatrix{size(verts)...}(verts))
end
export testprojectonplane

"""projects convex polygon, represented by `vertices`, onto `surface` along vector `normal`. Assumes original polygon is convex and that the projection will be convex. No guarantee that this will be true but for smoothly curved surfaces that are not varying too quickly relative to the size of the polygon it should be true."""
function planarpoly(projectedpoints::AbstractMatrix{T}) where{T}
    planarpoints,toworld,_ = projectonbestfitplane(projectedpoints)
    numpts = size(planarpoints)[2]
    temp = SMatrix{2,numpts}(planarpoints[1:2,:])
    vecofpts = collect(reinterpret(reshape,SVector{2,T},temp))
    #this code is inefficient because of different data structures used by convex polygon and most other code for representing lists of vertices.
     return ConvexPolygon(toworld,vecofpts)
end

spherepoint(radius,θ,ϕ) = radius .* SVector(cos(θ)sin(ϕ),sin(θ),cos(θ)cos(ϕ))
export spherepoint

"""Computes points on the edges of the spherical rectangle defined by the range of θ,ϕ. This is used to determine lattice boundaries on the eyebox surface."""
function spherepoints(radius, θmin,θmax,ϕmin,ϕmax)
    θedges =  [spherepoint(radius,ϕ,θ) for θ in θmin:.5:θmax, ϕ in (ϕmin,ϕmax)]
    ϕedges =  [spherepoint(radius,ϕ,θ) for ϕ in ϕmin:.5:ϕmax, θ in (θmin,θmax)]
    allpoints = vcat(reshape(θedges,reduce(*,size(θedges))),reshape(ϕedges,reduce(*,size(ϕedges))))

    reshape(reinterpret(Float64,allpoints),3,length(allpoints)) #return points as 3xn matrix with points as columns
end
export spherepoints

"""given a total fov in θ  and ϕ compute sample points on the edges of the spherical rectangle."""
spherepoints(radius,θ,ϕ) = spherepoints(radius,-θ/2,θ/2,-ϕ/2,ϕ/2)

function testprojection()
    pts = spherepoints(1.0,-.2,-.2,1.0,1.1)
    surf = Plane(0.0,0.0,-1.0,0.0,0.0,0.0)
    dir = [0.0,0.0,-1.0]
    project(pts,dir,surf)
end
export testprojection

function bounds(pts::AbstractMatrix{T}) where{T}
    return [extrema(row) for row in eachrow(pts)]
end
export bounds

function eyeboxbounds(eyebox::OpticSim.Plane,dir::AbstractVector, radius,fovθ,fovϕ)
    pts = spherepoints(radius,fovθ,fovϕ)
    projectedpts = project(pts,dir,eyebox)
    return bounds(projectedpts)
end
export eyeboxbounds

function boxtiles(bbox,lattice)
    tiles = tilesinside(bbox[1][1],bbox[2][1],bbox[1][2],bbox[2][2],lattice)
end
export boxtiles

eyeboxtiles(eyebox,dir,radius,fovθ,fovϕ,lattice) = boxtiles(eyeboxbounds(eyebox,dir,radius,fovθ,fovϕ),lattice)
export eyeboxtiles

function spherepolygon(vertices,projectiondirection,sph::Sphere)
    @assert typeof(vertices) <: SMatrix
    projectedpts = project(vertices,projectiondirection,sph)
    return planarpoly(projectedpts) #polygon on best fit plane to projectedpts
end

#TODO need to figure out what to use as the normal (eventually this will need to take into account the part of the eyebox the lenslet should cover)
#TODO write test for planarpoly and generation of Paraxial lens system using planar hexagon as lenslet.
function spherepolygons(eyebox::Plane{T,N},dir,radius,fovθ,fovϕ,lattice) where{T,N}
    #if fovθ, fovϕ are in degrees convert to radians. If they are unitless then the assumption is that the represent radians
    eyeboxz = eyebox.pointonplane[3]
    θ = upreferred(fovθ) #converts to radians if in degrees
    ϕ = upreferred(fovϕ) #converts to radians if in degrees
    tiles = eyeboxtiles(eyebox,dir,radius,θ,ϕ,lattice)
    shapes = Vector{ConvexPolygon{6,T}}(undef,0)
    coordinates = Vector{Tuple{Int64,Int64}}(undef,0)
    for coords in eachcol(tiles)
        twodverts = tilevertices(coords,lattice)
        numverts = size(twodverts)[2]
        temp = MMatrix{1,numverts,T}(undef)
        fill!(temp,eyeboxz)
        zerorow = SMatrix{1,numverts,T}(temp)
        vertices = SMatrix{N,numverts,T}(vcat(twodverts,zerorow))
        @assert typeof(vertices) <: SMatrix
        push!(shapes,spherepolygon(vertices,-dir,Sphere(radius)))
        push!(coordinates,Tuple{T,T}(coords))
    end
    shapes,coordinates
end

function spherelenslets(eyebox::Plane{T,N},dir,radius,fovθ,fovϕ,lattice) where{T,N}
    lenspolys,tilecoords = spherepolygons(eyebox,dir,radius,fovθ,fovϕ,lattice)
    result = Vector{ParaxialLens{T}}(undef,length(lenspolys))
    empty!(result)
    for poly in lenspolys
        lenslet = ParaxialLensConvexPoly(2.0,poly,SVector{2,T}(T(0),T(0)))
        push!(result,lenslet)
    end
    return result,tilecoords
end
export spherelenslets

function testspherelenslets()
    spherelenslets(Plane(0.0,0.0,1.0,0.0,0.0,12.0),[0.0,0.0,-1.0],30.0,55°,45°,HexBasis1())
end
export testspherelenslets

function testspherepolygons()
    eyebox = Plane(0.0,0.0,1.0,0.0,0.0,12.0)
    dir = [0.0,0.0,-1.0]
    radius = 30.0
    fovθ = 55°
    fovϕ = 45°
    lattice = HexBasis1()

    return spherepolygons(eyebox,dir,radius,fovθ,fovϕ,lattice)
end
export testspherepolygons

# function testeyeboxtiles()
#     tiles = eyeboxtiles(Plane(0.0,0.0,1.0,0.0,0.0,12.0),[0.0,0.0,-1.0],30,deg2rad(55),deg2rad(45),HexBasis1())
#     drawcells(HexBasis1(),10,tiles)
#     tiles = eyeboxtiles(Plane(0.0,0.0,1.0,0.0,0.0,12.0),[0.0,0.0,-1.0],30,deg2rad(25),deg2rad(45),HexBasis1())
#     drawcells(HexBasis1(),10,tiles)
# end
# export testeyeboxtiles


text:
@testset "Citations" begin
    product_years = [2009, 2013, 2018]
    @testset "WOA$(WOA.my_product_year(y))" for y in product_years
        tracers = ["t", "s", "I", "C", "o", "O", "A", "i", "p", "n"]
        @testset "$(WOA.my_varname(t))" for t in tracers
            citation_str = WOA.citation(y, t)
            @test citation_str isa String
            println("Citation for $(WOA.my_varname(t)) from WOA$(WOA.my_product_year(y))")
            println(citation_str * "\n")
        end
    end
end


text:
"""
    bootstrap!{T}(r::Matrix{T}, m::LinearMixedModel{T}, f!::Function;
        β=fixef(m), σ=sdest(m), θ=getθ(m))

Overwrite columns of `r` with the results of applying the mutating extractor `f!`
to parametric bootstrap replications of model `m`.

The signature of `f!` should be
    f!{T}(v::AbstractVector{T}, m::LinearMixedModel{T})

# Named Arguments

`β::Vector{T}`, `σ::T`, and `θ::Vector{T}` are the values of the parameters in `m`
for simulation of the responses.
"""
function bootstrap!{T}(r::Matrix{T}, m::LinearMixedModel{T}, f!::Function;
    β=fixef(m), σ=sdest(m), θ=getθ(m))
    y₀ = copy(model_response(m)) # to restore original state of m
    for i in 1 : size(r, 2)
        f!(Compat.view(r, :, i), refit!(simulate!(m, β = β, σ = σ, θ = θ)))
    end
    refit!(m, y₀)               # restore original state of m
    r
end

"""
    reevaluateAend!(m::LinearMixedModel)

Reevaluate the last column of `m.A` from `m.trms`.  This function should be called
after updating the response, `m.trms[end]`.
"""
function reevaluateAend!(m::LinearMixedModel)
    A, trms, sqrtwts, wttrms = m.A, m.trms, m.sqrtwts, m.wttrms
    wttrmn = wttrms[end]
    if !isempty(sqrtwts)
        A_mul_B!(wttrmn, sqrtwts, trms[end])
    end
    for i in eachindex(wttrms)
        Ac_mul_B!(A[i, end], wttrms[i], wttrmn)
    end
    m
end

"""
    refit!{T}(m::LinearMixedModel{T}[, y::Vector{T}])

Refit the model `m` after installing response `y`.

If `y` is omitted the current response vector is used.
"""
refit!(m::LinearMixedModel) = fit!(cfactor!(resetθ!(reevaluateAend!(m))))
function refit!{T}(m::LinearMixedModel{T}, y)
    resp = m.trms[end]
    if length(y) ≠ size(resp, 1)
        throw(DimensionMismatch("length(y) = $(length(y)), should be $(size(resp, 1))"))
    end
    copy!(resp, y)
    refit!(m)
end

"""
    resetθ!(m::LinearMixedModel)

Reset the value of `m.θ` to the initial values and mark the model as not having been fit
"""
function resetθ!(m::LinearMixedModel)
    opt = m.opt
    opt.feval = -1
    opt.fmin = Inf
    setθ!(m, opt.initial) |> cfactor!
end

"""
    unscaledre!{T}(y::Vector{T}, M::ReMat{T}, b::Matrix{T})

Add unscaled random effects defined by `M` and `b` to `y`.
"""
function unscaledre!{T<:AbstractFloat,S,R<:Integer}(y::Vector{T}, M::ScalarReMat{T,S,R}, b::Matrix{T})
    z = M.z
    if length(y) ≠ length(z) || size(b, 1) ≠ 1
        throw(DimensionMismatch())
    end
    inds = M.f.refs
    @inbounds for i in eachindex(y)
        y[i] += b[inds[i]] * z[i]
    end
    y
end

"""
    unscaledre!{T}(y::AbstractVector{T}, M::ReMat{T}, L::LowerTriangular{T})

Add unscaled random effects defined by `M` and `L * randn(1, length(M.f.pool))` to `y`.
"""
function unscaledre!{T}(y::AbstractVector{T}, M::ScalarReMat{T}, L::LowerTriangular{T})
    unscaledre!(y, M, A_mul_B!(L, randn(1, length(M.f.pool))))
end

function unscaledre!{T}(y::AbstractVector{T}, M::VectorReMat{T}, b::DenseMatrix{T})
    Z = M.z
    k, n = size(Z)
    l = length(M.f.pool)
    if length(y) ≠ n || size(b) ≠ (k, l)
        throw(DimensionMismatch("length(y) = $(length(y)), size(M) = $(size(M)), size(b) = $(size(b))"))
    end
    inds = M.f.refs
    for i in eachindex(y)
        ii = inds[i]
        for j in 1:k
            y[i] += Z[j,i] * b[j, ii]
        end
    end
    y
end

unscaledre!(y::AbstractVector, M::VectorReMat, L::LowerTriangular) =
    unscaledre!(y, M, A_mul_B!(L, randn(size(M.z, 1), length(M.f.pool))))

"""
    simulate!(m::LinearMixedModel; β=fixef(m), σ=sdest(m), θ=getΘ(m))

Overwrite the response (i.e. `m.trms[end]`) with a simulated response vector from model `m`.
"""
function simulate!{T}(m::LinearMixedModel{T}; β = coef(m), σ = sdest(m), θ = T[])
    if !isempty(θ)
        setθ!(m, θ)
    end
    trms, Λ = m.trms, m.Λ
    y = randn!(model_response(m)) # initialize to standard normal noise
    for j in eachindex(Λ)         # add the unscaled random effects
        unscaledre!(y, trms[j], Λ[j])
    end
                                  # scale by σ and add fixed-effects contribution
    BLAS.gemv!('N', 1.0, trms[end - 1], β, σ, y)
    m
end

StatsBase.model_response(m::LinearMixedModel) = vec(m.trms[end])


text:
function process_letter(input::Char)
    input_n = Int(input)

    if 97 <= input_n <= 122
        return Char(219 - input_n)
    elseif 48 <= input_n <= 57
        return input
    else
        error("Input string not fully stripped.")
    end
end

function insert_space(input::String)
    input_array = collect(input)
    len = length(input_array)

    new_array = Char[]

    added_len = len + div(len-1, 5)

    for i in 1:added_len
        if i % 6 == 0
            push!(new_array, ' ')
        else
            push!(new_array, input_array[i - div(i, 6)])
        end
    end

    return String(new_array)
end

function strip_string(input::String)
    lc_array = collect(lowercase(input))
    stripped = [ch for ch in lc_array if (97 <= Int(ch) <= 122) | (48 <= Int(ch) <= 57)]
    return stripped
end

function encode(input::String)

    stripped = strip_string(input)
    encoded = String([process_letter(ch) for ch in stripped])

    return insert_space(encoded)
end

function decode(input::String)
    stripped = strip_string(input)

    return String([process_letter(ch) for ch in stripped])
end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
using LinearAlgebra, Distributions, JLD2, DataDeps
using OceanBoundaryLayerParameterEstimation

using ParameterEstimocean
using ParameterEstimocean.Parameters: closure_with_parameters
using ParameterEstimocean.Observations: set!
using ParameterEstimocean.PseudoSteppingSchemes

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputWriters: TimeInterval, JLD2OutputWriter
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity
using Oceananigans: run!

regenerate_synthetic_observations = true
architecture = CPU()
# Δt = 5.0

# parameter_names = (:CᵂwΔ,  :Cᵂu★, :Cᴰ,
#                    :Cˢc,   :Cˢu,  :Cˢe,
#                    :Cᵇc,   :Cᵇu,  :Cᵇe,
#                    :Cᴷc⁻,  :Cᴷu⁻, :Cᴷe⁻,
#                    :Cᴷcʳ,  :Cᴷuʳ, :Cᴷeʳ,
#                    :CᴷRiᶜ, :CᴷRiʷ)

# parameter_set = ParameterSet{CATKEVerticalDiffusivity}(Set(parameter_names),
#                              nullify = Set([:Cᴬu, :Cᴬc, :Cᴬe]))

# closure = closure_with_parameters(CATKEVerticalDiffusivity(Float64;), parameter_set.settings)

# true_parameters = (CᵂwΔ = 4.46,    Cᵂu★ = 4.56,   Cᴰ = 2.91,
#                    Cˢc = 0.426,    Cˢu = 0.628,   Cˢe = 0.711,
#                    Cᵇc = 0.0723,   Cᵇu = 0.596,   Cᵇe = 0.637,
#                    Cᴷc⁻ = 0.343,   Cᴷu⁻ = 0.343,  Cᴷe⁻ = 1.42,
#                    Cᴷcʳ = -0.891,  Cᴷuʳ = -0.721, Cᴷeʳ = 1.50,
#                    CᴷRiᶜ = 2.03,   CᴷRiʷ = 0.101)

# true_closure = closure_with_parameters(closure, true_parameters)

# directory = "lesbrary_catke_perfect_model_6_days"
# isdir(directory) || mkpath(directory)

Δt = 10minutes

parameter_set = RiBasedParameterSet

# closure = closure_with_parameters(RiBasedVerticalDiffusivity(), parameter_set.settings)
closure = RiBasedVerticalDiffusivity()

true_parameters = parameter_set.settings
true_closure = closure

directory = "lesbrary_ri_based_perfect_model_6_days"
isdir(directory) || mkpath(directory)

###
### Generate data at high resolution using 6-day LESbrary as a template for ICs and BCs but CATKE "perfect model" parameters.
###

function run_synthetic_single_column_simulation!(observation, prefix;
                                                Nensemble = 30,
                                                architecture = CPU(),
                                                closure = ConvectiveAdjustmentVerticalDiffusivity(),
                                                Δt = 10.0,
                                                tracers = (:b, :e)
                                            )

    grid = RectilinearGrid(architecture,
                            size = observation.grid.Nz,
                            topology = (Flat, Flat, Bounded),
                            z = (-observation.grid.Lz, 0))

    p = observation.metadata.parameters
    Qᵘ = p.boundary_condition_u_top
    Qᵇ = p.boundary_condition_θ_top
    N² = p.N²_deep

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(; grid, tracers, closure,
                                          buoyancy = BuoyancyTracer(),
                                          boundary_conditions = (u=u_bcs, b=b_bcs),
                                          coriolis = FPlane(f = observation.metadata.coriolis.f))

    set!(model, observation, 1)

    init_with_parameters(file, model) = file["parameters"] = (; Qᵇ, Qᵘ, Δt, N², tracers)

    simulation = Simulation(model; Δt, stop_time=observation.times[end])

    model = simulation.model

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(0.5days),
                                                          array_type = Array{Float64},
                                                          prefix = prefix,
                                                          field_slicer = nothing,
                                                          init = init_with_parameters,
                                                          force = true)

    run!(simulation)
end

fields_by_case = Dict(
   "free_convection" => (:b, :e),
   "weak_wind_strong_cooling" => (:b, :u, :v, :e),
   "strong_wind_weak_cooling" => (:b, :u, :v, :e),
   "strong_wind" => (:b, :u, :v, :e),
   "strong_wind_no_rotation" => (:b, :u, :e)
)

regenerate_synthetic_observations && begin

    transformation = (b = Transformation(normalization=ZScore()),
                    u = Transformation(normalization=ZScore()),
                    v = Transformation(normalization=ZScore()),
                    e = Transformation(normalization=RescaledZScore(1e-1)))

    # lesbrary_observations = SixDaySuite(; transformation=transformation, times=nothing, Nz=128)

    field_names = (:b, :u, :v, :e)

    six_day_suite_path(case) = "six_day_suite_2m/$(case)_instantaneous_statistics.jld2"

    for (case, forward_map_names) in zip(keys(fields_by_case), values(fields_by_case))

        lesbrary_data_path = @datadep_str six_day_suite_path(case)

        lesbrary_observation = SyntheticObservations(lesbrary_data_path; transformation,
                                                times=nothing,
                                                field_names=(:b, :u, :v, :e),
                                                forward_map_names)

        prefix = joinpath(directory, case)
        run_synthetic_single_column_simulation!(lesbrary_observation, prefix;
                                                Nensemble=1,
                                                architecture,
                                                closure=true_closure,
                                                Δt)
    end
end

###
### Build coarse-grained SyntheticObservations from data.
###

Nz = 128
Nensemble = 20

training_times = [1.0day, 1.5days, 2.0days, 2.5days, 3.0days]
validation_times = [3.0days, 3.5days, 4.0days]
testing_times = [4.0days, 4.5days, 5.0days, 5.5days, 6.0days]

path_fn(case) = joinpath(directory, case) * ".jld2"

begin
    file = jldopen(path_fn("free_convection"))
    @show iterations = parse.(Int, keys(file["timeseries/t"]))
    @show times = [file["timeseries/t/$i"] for i in iterations]
    close(file)
end

build_prior(name) = ScaledLogitNormal(bounds=bounds(name, parameter_set))
free_parameters = FreeParameters(named_tuple_map(parameter_set.names, build_prior))

function inverse_problem(Nensemble, times)
    observations = SyntheticObservationsBatch(path_fn, transformation, times, Nz; datadep=false, architecture)
    simulation = lesbrary_ensemble_simulation(observations; Nensemble, architecture, closure, Δt)

    # simulation = ensemble_column_model_simulation(observations;
    #                                               Nensemble,
    #                                               architecture,
    #                                               tracers = (:b, :e),
    #                                               closure)

    simulation.Δt = Δt

    inverse_problem = InverseProblem(observations, simulation, free_parameters; output_map = ConcatenatedOutputMap())
    return inverse_problem
end

###
### Generate preliminary inverse problems.
###

training = inverse_problem(Nensemble, training_times)
validation = inverse_problem(Nensemble, validation_times)
testing = inverse_problem(Nensemble, testing_times)

resampler = Resampler(resample_failure_fraction=0.5, acceptable_failure_fraction=1.0)
eki = EnsembleKalmanInversion(training; resampler, convergence_rate=0.8)

stepping_scheme = ConstantConvergence(convergence_ratio = 0.7)
iterate!(eki; iterations = 1, pseudo_stepping = stepping_scheme)

visualize!(training, true_parameters;
                    field_names = [:u, :v, :b, :e],
                    directory,
                    filename = "true_parameter_training_realizations.png"
                    )


text:
# Circular Manifold definition and associated arithmetic

function get2DMu(mus, Lambdas; diffop::Function=-, periodicmanifold::Function=(x)->x, initrange::Tuple{Float64, Float64}=(-1e-5,1e-5) )::Float64
  # TODO: can be solved as the null space basis, but requires proper scaling
  gg = (res, x) -> solveresid2DLinear!(res, x, mus, Lambdas, diffop=diffop)
  initr = initrange[2]-initrange[1]
  x0 = [initr*rand()+initrange[1]]
  r = NLsolve.nlsolve(gg, x0)
  xs = sign(r.zero[1] - x0[1])*1e-3 .+ r.zero
  r = NLsolve.nlsolve(gg, xs)
  return periodicmanifold(r.zero[1])
end

function get2DMuMin(mus, Lambdas; diffop::Function=-, periodicmanifold::Function=(x)->x, initrange::Tuple{Float64, Float64}=(-1e-5,1e-5), method=Optim.Newton(), Λμ::Bool=false )::Float64
  # TODO: can be solved as the null space basis, but requires proper scaling
  res = zeros(1)
  # @show round.(mus, digits=3)
  gg = (x) -> (solveresid2DLinear(res, x, mus, Lambdas, diffop=diffop))^2
  initr = initrange[2]-initrange[1]
  # TODO -- do we need periodicmanifold here?
  x0 = [initr*rand()+initrange[1]]
  r = Optim.optimize(gg, x0, method)
  # @show r.minimizer[1]
  return periodicmanifold(r.minimizer[1])
end

#