text:
using Requires
using ..Mesh.Grids: EveryDirection, VerticalDirection, HorizontalDirection
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
end
using KernelAbstractions.Extras: @unroll

const _M = Grids._M

@doc """
    knl_apply_filter!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{direction},
                      Q, ::Val{states}, filtermatrix,
                      elems) where {dim, N, nstate, states, direction}

Computational kernel: Applies the `filtermatrix` to the `states` of `Q`.

The `direction` argument is used to control if the filter is applied in the
horizontal and/or vertical reference directions.
""" knl_apply_filter!
@kernel function knl_apply_filter!(
    ::Val{dim},
    ::Val{N},
    ::Val{nstate},
    ::Val{direction},
    Q,
    ::Val{states},
    filtermatrix,
    elems,
) where {dim, N, nstate, direction, states}
    @uniform begin
        FT = eltype(Q)

        Nq = N + 1
        Nqk = dim == 2 ? 1 : Nq

        if direction isa EveryDirection
            filterinξ1 = filterinξ2 = true
            filterinξ3 = dim == 2 ? false : true
        elseif direction isa HorizontalDirection
            filterinξ1 = true
            filterinξ2 = dim == 2 ? false : true
            filterinξ3 = false
        elseif direction isa VerticalDirection
            filterinξ1 = false
            filterinξ2 = dim == 2 ? true : false
            filterinξ3 = dim == 2 ? false : true
        end

        nfilterstates = length(states)
    end

    s_filter = @localmem FT (Nq, Nq)
    s_Q = @localmem FT (Nq, Nq, Nqk, nfilterstates)
    l_Qfiltered = @private FT (nfilterstates,)

    e = @index(Group, Linear)
    i, j, k = @index(Local, NTuple)

    @inbounds begin
        s_filter[i, j] = filtermatrix[i, j]

        @unroll for fs in 1:nfilterstates
            l_Qfiltered[fs] = zero(FT)
        end

        ijk = i + Nq * ((j - 1) + Nq * (k - 1))

        @unroll for fs in 1:nfilterstates
            s_Q[i, j, k, fs] = Q[ijk, states[fs], e]
        end

        if filterinξ1
            @synchronize
            @unroll for n in 1:Nq
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += s_filter[i, n] * s_Q[n, j, k, fs]
                end
            end

            if filterinξ2 || filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    s_Q[i, j, k, fs] = l_Qfiltered[fs]
                    l_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ2
            @synchronize
            @unroll for n in 1:Nq
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += s_filter[j, n] * s_Q[i, n, k, fs]
                end
            end

            if filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    s_Q[i, j, k, fs] = l_Qfiltered[fs]
                    l_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ3
            @synchronize
            @unroll for n in 1:Nqk
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += s_filter[k, n] * s_Q[i, j, n, fs]
                end
            end
        end

        # Store result
        ijk = i + Nq * ((j - 1) + Nq * (k - 1))
        @unroll for fs in 1:nfilterstates
            Q[ijk, states[fs], e] = l_Qfiltered[fs]
        end

        @synchronize
    end
end

@kernel function knl_apply_TMAR_filter!(
    ::Val{nreduce},
    ::Val{dim},
    ::Val{N},
    Q,
    ::Val{filterstates},
    vgeo,
    elems,
) where {nreduce, dim, N, filterstates}
    @uniform begin
        FT = eltype(Q)

        Nq = N + 1
        Nqj = dim == 2 ? 1 : Nq

        nfilterstates = length(filterstates)
        nelemperblock = 1
    end

    l_Q = @private FT (nfilterstates, Nq)
    l_MJ = @private FT (Nq,)

    s_MJQ = @localmem FT (Nq * Nqj, nfilterstates)
    s_MJQclipped = @localmem FT (Nq * Nqj, nfilterstates)

    e = @index(Group, Linear)
    i, j = @index(Local, NTuple)

    @inbounds begin
        # loop up the pencil and load Q and MJ
        @unroll for k in 1:Nq
            ijk = i + Nq * ((j - 1) + Nqj * (k - 1))

            @unroll for sf in 1:nfilterstates
                s = filterstates[sf]
                l_Q[sf, k] = Q[ijk, s, e]
            end

            l_MJ[k] = vgeo[ijk, _M, e]
        end

        @unroll for sf in 1:nfilterstates
            MJQ, MJQclipped = zero(FT), zero(FT)

            @unroll for k in 1:Nq
                MJ = l_MJ[k]
                Qs = l_Q[sf, k]
                Qsclipped = Qs ≥ 0 ? Qs : zero(Qs)

                MJQ += MJ * Qs
                MJQclipped += MJ * Qsclipped
            end

            ij = i + Nq * (j - 1)

            s_MJQ[ij, sf] = MJQ
            s_MJQclipped[ij, sf] = MJQclipped
        end
        @synchronize

        @unroll for n in 11:-1:1
            if nreduce ≥ 2^n
                ij = i + Nq * (j - 1)
                ijshift = ij + 2^(n - 1)
                if ij ≤ 2^(n - 1) && ijshift ≤ Nq * Nqj
                    @unroll for sf in 1:nfilterstates
                        s_MJQ[ij, sf] += s_MJQ[ijshift, sf]
                        s_MJQclipped[ij, sf] += s_MJQclipped[ijshift, sf]
                    end
                end
                @synchronize
            end
        end

        @unroll for sf in 1:nfilterstates
            qs_average = s_MJQ[1, sf]
            qs_clipped_average = s_MJQclipped[1, sf]

            r = qs_average > 0 ? qs_average / qs_clipped_average : zero(FT)

            s = filterstates[sf]
            @unroll for k in 1:Nq
                ijk = i + Nq * ((j - 1) + Nqj * (k - 1))

                Qs = l_Q[sf, k]
                Q[ijk, s, e] = Qs ≥ 0 ? r * Qs : zero(Qs)
            end
        end
    end
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
import Base: @deprecate, @deprecate_binding, depwarn

###
### v0.15 deprecations
###

### Add empty exists method for JLD,MAT to extend to smooth over deprecation process PR#790
export exists
function exists end

### Changed in PR#776
@deprecate create_dataset(parent::Union{File,Group}, path::AbstractString, dtype::Datatype, dspace::Dataspace,
    lcpl::Properties, dcpl::Properties, dapl::Properties, dxpl::Properties) HDF5.Dataset(HDF5.API.h5d_create(parent, path, dtype, dspace, lcpl, dcpl, dapl), HDF5.file(parent), dxpl) false

### Changed in PR#798
@deprecate get_dims(dspace::Union{Dataspace,Dataset,Attribute}) get_extent_dims(dspace) false
@deprecate set_dims!(dset::Dataspace) set_extent_dims(dset) false

### Changed in PR #844
@deprecate silence_errors(f::Function) f()

for name in names(API; all=true)
    if name ∉ names(HDF5; all=true) && startswith(uppercase(String(name)), "H")
        depmsg = ", use HDF5.API.$name instead."
        @eval Base.@deprecate_binding $name API.$name false $depmsg
    end
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ b2d786ec-7f73-11ea-1a0c-f38d7b6bbc1e
md"""
# The Basel problem

_Leonard Euler_ proved in 1741 that the series

```math
\frac{1}{1} + \frac{1}{4} + \frac{1}{9} + \cdots
```

converges to

```math
\frac{\pi^2}{6}.
```
"""

# ╔═╡ b2d79330-7f73-11ea-0d1c-a9aad1efaae1
n = 1:100000

# ╔═╡ b2d79376-7f73-11ea-2dce-cb9c449eece6
seq = n .^ -2

# ╔═╡ b2d792c2-7f73-11ea-0c65-a5042701e9f3
sqrt(sum(seq) * 6.0)

# ╔═╡ Cell order:
# ╟─b2d786ec-7f73-11ea-1a0c-f38d7b6bbc1e
# ╠═b2d792c2-7f73-11ea-0c65-a5042701e9f3
# ╠═b2d79330-7f73-11ea-0d1c-a9aad1efaae1
# ╠═b2d79376-7f73-11ea-2dce-cb9c449eece6


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

using Markdown
using InteractiveUtils

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

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, StatsPlots, Plots
	Random.seed!(0)
 end

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

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	mu = 20
	d1, d2 = Normal(mu,mu), Exponential(mu)

	n = 100
	data1 = rand(d1,n)
	data2 = rand(d2,n)
end

# ╔═╡ bdabe54e-250a-11eb-1714-25afa772347e
begin
	qqnorm(data1, c=:blue, ms=3, msw=0, label="Normal Data")
	qqnorm!(data2, c=:red, ms=3, msw=0, label="Exponential Data",
			xlabel="Normal Theoretical Quantiles",
			ylabel="Data Quantiles", legend=true)
end

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

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═bdabe54e-250a-11eb-1714-25afa772347e
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using Endpoints, Base.Test, HTTP

include("json2.jl")

Endpoints.@GET (@__MODULE__) "test" ()->"testing GET resource"
Endpoints.@GET @__MODULE__ "test/{var::String}" var->"testing $var resource"
Endpoints.@GET @__MODULE__ "test/int/{var::Int}" var->"testing int $var resource"
Endpoints.@GET @__MODULE__ "test/type/{::Type{T}}" T->parse(T, "1")
# Issue #2
path = "/test/dynamic/path"
Endpoints.@GET path ()->"testing path from variable"

type Bod
    id::Int
    nm::String
end

Endpoints.@POST "test/post" body body->"do you like my $body"
Endpoints.@POST "test/post/andarg/{arg}" body (arg, body)->"do you like my $arg and $body"
Endpoints.@POST "test/post/typedbod" body::Bod bod->"nice bod: $(bod.id), $(bod.nm)"

HTTP.serve(handler = Endpoints.handler)


text:
facts("Utilities") do

    context("Base.abs returns all positive numbers") do
        @fact sum((abs(cl .- op) .>= 0).values) --> length(cl)
    end

    context("pad correctly pads NaN values for non-existent values") do
        @fact sum((abs(cl .- op) .>= 0).values) --> length(cl)
    end

    context("pad correctly pads NaN values for missing values") do
        @fact sum((abs(cl .- op) .>= 0).values) --> length(cl)
    end
end


text:
export Isaac64, rand, sysRandomSeed
type IsaacState
    a::Uint64
    b::Uint64
    aa::Uint64
    bb::Uint64
    cc::Uint64
    RANDSIZL::Int
    RANDSIZ::Int
end
type Isaac64
    state::IsaacState
    randrsl::Vector{Uint64}
    mm ::Vector{Uint64}
    randcnt::Int
end
function Isaac64(seed::Vector{Uint64}=Uint64[])
    state,randrsl, mm = main(seed)
    return Isaac64(state, randrsl, mm, state.RANDSIZ+1)
end
function Isaac64(state::IsaacState, randrsl::Vector{Uint64}, mm::Vector{Uint64})
    return Isaac64(state, randrsl, mm, state.RANDSIZ+1)
end

function rand{T<:Union(Uint8,Uint16, Uint32, Uint64, Int8, Int16, Int32, Int64)}(rng::Isaac64, t::Type{T})
    return convert(t, rand(rng))
end
function rand(rng::Isaac64, t::Type{Float64})
    return float64(string("0.", string(int128(rand(rng)))))
end
function rand(rng::Isaac64)
    if rng.randcnt != 1
        return rng.randrsl[rng.randcnt -=1]
    else
        rng.randcnt=rng.state.RANDSIZ
        isaac64rng(rng.state, rng.randrsl, rng.mm)
        return rng.randrsl[rng.randcnt]
    end
end
function ind(s::IsaacState, mm, x)
    return mm[1 + (x & ((s.RANDSIZ-1)<<3))>>3]
end

function rngstep(s::IsaacState, randrsl, mm, a_mix, m1, m2, r, i)
    x = mm[m1 + i]
    s.a = a_mix + mm[m2 + i]
    mm[m1 + i] = y = ind(s, mm, x) + s.a + s.b
    randrsl[r + i] = s.b = ind(s, mm, y >> s.RANDSIZL) + x
end
function isaac64rng(rng::Isaac64)
    isaac64rng(rng.state, rng.randrsl, rng.mm)
end
function isaac64rng(s::IsaacState, randrsl, mm)
    s.a = s.aa
    s.cc += 1
    s.b = s.bb + s.cc
    m1 = r = 0
    mid = m2 = s.RANDSIZ >> 1
    for i = 1:4:mid-2
        rngstep(s, randrsl, mm, ~(s.a $ (s.a << 21)), m1, m2, r, i)
        rngstep(s, randrsl, mm,   s.a $ (s.a >> 5),   m1, m2, r, i + 1)
        rngstep(s, randrsl, mm,   s.a $ (s.a << 12),  m1, m2, r, i + 2)
        rngstep(s, randrsl, mm,   s.a $ (s.a >> 33),  m1, m2, r, i + 3)
    end
    m2 = 0
    r = m1 = s.RANDSIZ >> 1
    for i = 1:4:mid-2
        rngstep(s, randrsl, mm, ~(s.a $ (s.a << 21)), m1, m2, r, i)
        rngstep(s, randrsl, mm,   s.a $ (s.a >> 5),   m1, m2, r, i + 1)
        rngstep(s, randrsl, mm,   s.a $ (s.a << 12),  m1, m2, r, i + 2)
        rngstep(s, randrsl, mm,   s.a $ (s.a >> 33),  m1, m2, r, i + 3)
    end
    s.bb = s.b
    s.aa = s.a;
end

macro mix()
    quote
        a -= e; f $= h >> 9;  h += a
        b -= f; g $= a << 9;  a += b
        c -= g; h $= b >> 23; b += c
        d -= h; a $= c << 15; c += d
        e -= a; b $= d >> 14; d += e
        f -= b; c $= e << 20; e += f
        g -= c; d $= f >> 17; f += g
        h -= d; e $= g << 14; g += h
    end
end
function randinit(randrsl, mm, RANDSIZL, RANDSIZ, flag=false)
    a = b = c = d = e = f = g = h = 0x9e3779b97f4a7c13
    for i = 1:4
        @mix
    end
    for i = 1:8:RANDSIZ
        if flag
            a += mm[i  ]; b += mm[i+1]; c += mm[i+2]; d += mm[i+3]
            e += mm[i+4]; f += mm[i+5]; g += mm[i+6]; h += mm[i+7]
        end
        @mix
            mm[i  ] = a; mm[i+1] = b; mm[i+2] = c; mm[i+3] = d
            mm[i+4] = e; mm[i+5] = f; mm[i+6] = g; mm[i+7] = h
    end
    if flag
        for i = 1:8:RANDSIZ
              a += mm[i  ]; b += mm[i+1]; c += mm[i+2]; d += mm[i+3]
              e += mm[i+4]; f += mm[i+5]; g += mm[i+6]; h += mm[i+7]
            @mix
              mm[i  ] = a; mm[i+1] = b; mm[i+2] = c; mm[i+3] = d
              mm[i+4] = e; mm[i+5] = f; mm[i+6] = g; mm[i+7] = h
        end
    end
    z = 0x0000000000000000
    s2 = IsaacState(z, z, z, z, z, RANDSIZL, RANDSIZ)
    isaac64rng(s2, randrsl, mm)
    return s2
end
function main(randrsl::Vector{Uint64}=Uint64[])
    RANDSIZL = 8
    RANDSIZ = length(randrsl)
    if RANDSIZ == 0
        RANDSIZ = 256
        resize!(randrsl, RANDSIZ)
        fill!(randrsl, uint64(0))
    elseif RANDSIZ % 8 != 0
        error("dimension of seeding array must be a factor of 8")
    end
    mm = zeros(Uint64, RANDSIZ)
    s = randinit(randrsl, mm, RANDSIZL, RANDSIZ, true)
    for i = 1:2
        isaac64rng(s, randrsl, mm)
    end
    return s, randrsl, mm
end

function sysRandomSeed{T<:Number}(size=256, t::Type{T}=Uint64, sysrand_path="/dev/urandom")
    seed_array = resize!(t[], size)
    f=open(sysrand_path)
    read!(f, seed_array)
    close(f)
    return seed_array
end

max_stars_repo_name: pik/isaac-jl

max_issues_repo_name: pik/isaac-jl

max_forks_repo_name: pik/isaac-jl


text:
module Bumper

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Bumper" => Ahorn.EntityPlacement(
        Maple.Bumper
    )
)

Ahorn.nodeLimits(entity::Maple.Bumper) = 0, 1

sprite = "objects/Bumper/Idle22.png"

function Ahorn.selection(entity::Maple.Bumper)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        return [Ahorn.getSpriteRectangle(sprite, x, y), Ahorn.getSpriteRectangle(sprite, nx, ny)]
    end

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bumper)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        theta = atan(y - ny, x - nx)
        Ahorn.drawArrow(ctx, x, y, nx + cos(theta) * 8, ny + sin(theta) * 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bumper, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end


text:
"""
    ks89RecursiveResidual(setting; indices, k)
Calculate recursive residual for the given regression setting and observation.
# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `indices::ArrayInt,1`: Indices of observations used in the linear model.
- `k::Int`: Observation indice the recursive residual is calculated for.
# Notes
    This is a helper function for the ks89 function and it is not directly used.
# References
Kianifard, Farid, and William H. Swallow. "Using recursive residuals, calculated on
adaptively-ordered observations, to identify outliers in linear regression."
Biometrics (1989): 571-585.
"""
function ks89RecursiveResidual(setting::RegressionSetting, indices::Array{Int,1}, k::Int)
    X, y = @extractRegressionSetting setting
    return ks89RecursiveResidual(X, y, indices, k)
end

function ks89RecursiveResidual(X::Array{Float64,2}, y::Array{Float64,1}, indices::Array{Int,1}, k::Int)
    useX = X[indices, :]
    useY = y[indices]
    olsreg = ols(useX, useY)
    betas = coef(olsreg)
    XX = inv(useX'useX)
    w = (y[k] - X[k,:]' * betas) / sqrt(1 + X[k,:]' * XX * X[k,:])
    return w
end

"""
    ks89(setting; alpha = 0.05)
Perform the Kianifard & Swallow (1989) algorithm for the given regression setting.
# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `alpha::Float64`: Optional argument of the probability of rejecting the null hypothesis.
# Examples
```julia-repl
julia> reg0001 = createRegressionSetting(@formula(stackloss ~ airflow + watertemp + acidcond), stackloss)
julia> ks89(reg0001)
2-element Array{Int64,1}:
  4
 21
```
# References
Kianifard, Farid, and William H. Swallow. "Using recursive residuals, calculated on
adaptively-ordered observations, to identify outliers in linear regression."
Biometrics (1989): 571-585.
"""
function ks89(setting::RegressionSetting; alpha=0.05)
    X = designMatrix(setting)
    y = responseVector(setting)
    return ks89(X, y, alpha=alpha)
end

function ks89(X::Array{Float64,2}, y::Array{Float64,1}; alpha=0.05)
    stdres = studentizedResiduals(X, y)
    orderingindices = sortperm(abs.(stdres))
    n, p = size(X)
    basisindices = orderingindices[1:p]
    w = zeros(Float64, n)
    s = zeros(Float64, n)
    ws = zeros(Float64, n)
    @inbounds for i in (p + 1):n
        index = orderingindices[i]
        w[index] = ks89RecursiveResidual(X, y, basisindices, index)
        s[index] = jacknifedS(X, y, index)
        ws[index] = w[index] / s[index]
        basisindices = orderingindices[1:i]
    end
    td = TDist(n - p - 1)
    q = quantile(td, alpha)
    result = filter(i -> abs.(ws[i]) > abs(q), 1:n)
    return result
end

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
"""
    bet(
        mag::AbstractArray{<:Real, 3},
        vsz::NTuple{3, Real},
        betargs::AbstractString = "-m -n -f 0.5"
    ) -> Array{Bool, 3}

Interface to FSL's bet.

### Arguments
- `mag::AbstractArray{<:Real, 3}`: magnitude image
- `vsz::NTuple{3, Real}`: voxel size
- `betargs::AbstractString = "-m -n -f 0.5"`: bet options

### Returns
- `Array{Bool, 3}`: binary brain mask
"""
function bet(
    mag::AbstractArray{<:Real, 3},
    vsz::NTuple{3, Real},
    betargs::AbstractString = "-m -n -f 0.5"
)
    if Sys.which("bet") === nothing
        error("bet: command not found")
    end

    magfile = tempname(cleanup=false) * ".nii"
    maskfile = replace(magfile, ".nii" => "_mask.nii.gz")

    nii = NIVolume(mag, voxel_size=vsz)
    niwrite(magfile, nii)

    try
        # bet appends '_mask' to output filename regardless of user choice
        run(`bet $magfile $magfile $(split(betargs))`)
        nii = niread(maskfile)
    finally
        rm(magfile, force=true)
        rm(maskfile, force=true)
    end

    return Array{Bool, 3}(nii.raw)
end


text:
# # Statistics
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/07_Stats.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/07_Stats.ipynb)
#
# ## Missing values
#
# Julia, like R, has a value to represent
# [missing values](https://docs.julialang.org/en/v1/manual/missing/index.html):

data = [ 1.0, 2.0, missing, 4.0 ]

# This value implements three-valued logic:

false & missing

#-

true & missing

# You can use `ismissing` or `skipmissing` when necessary:

sum(data)

#-

ismissing.(data)

#-

sum(data[.!(ismissing.(data))])

#-

sum(skipmissing(data))

# ## DataFrames
#
# It is very useful to work with tabular data. One of the most simplest Julia
# packages for that is [DataFrames](http://juliadata.github.io/DataFrames.jl/stable/).

using DataFrames

# To read this kind of files, you can use the
# [CSV package](https://juliadata.github.io/CSV.jl/stable/).

using CSV

# For example, the `pdb_chain_taxonomy.tsv.gz` file that has a summary of the
# NCBI tax_id(s), scientific_name(s) and chain type for each PDB chain that
# has been processed in the SIFTS database. This table should be downloaded
# from the [SIFTS site](https://www.ebi.ac.uk/pdbe/docs/sifts/quick.html).
#
# ```julia
# table_path = download(
#     "ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/flatfiles/tsv/pdb_chain_taxonomy.tsv.gz",
#     "pdb_chain_taxonomy.tsv.gz")
#
# run(`gunzip $table_path`)
# ```
#
# If that fails, you can use the first lines stored in the data folder:
# ```julia
# using  JuliaForBioinformatics
# data_path = abspath(pathof(JuliaForBioinformatics), "..", "..", "data")
# table_path = joinpath(data_path, "pdb_chain_taxonomy_head.tsv")
# ```

table_path = "pdb_chain_taxonomy.tsv"

df = CSV.read(table_path,
    header = 2,  ## the header is in the second line
    delim = '\t',  ## delimiter is TAB instead of ','
    quotechar='`'  ## file don't use "" to quote, e.g.: "Bacillus coli" Migula 1895
    )

# ### Examples
#
# Select human PDB chains:

df[:TAX_ID] .== 9606

#-

df[df[:TAX_ID] .== 9606, [:PDB, :CHAIN]] |> unique

# You can use `|>` for easy function chaining.
#
# What are the species with more PDB chains?

count_df = by(df, :TAX_ID, Count = :TAX_ID => length)
sort!(count_df, :Count, rev=true)

# #### Exercise 1
#
# What are the species with more PDBs (not PDB chains)?
#
# **Hint:** You can use
# [anonymous functions](https://docs.julialang.org/en/v1/manual/functions/index.html#man-anonymous-functions-1):

f(x) = 2x + 1
g(x) = sin(π*x)

#-

x -> f(g(x))

# or function composition (using `∘`, `\circ<TAB>`):

f ∘ g

#-

## ...your solution...

# ## Plots
#
# They are multiple plotting packages in Julia. Here I will show
# [StatsPlots](https://github.com/JuliaPlots/StatsPlots.jl), an extension of
# [Plots](docs.juliaplots.org/latest/) for statistical plotting.
# However, if you love the grammar of graphics, you will be more comfortable
# with [Gadfly](https://github.com/GiovineItalia/Gadfly.jl).
#

using StatsPlots

#-

@df count_df bar(:TAX_ID, :Count)

#-

@df count_df marginalhist(:TAX_ID, :Count)

#-

@df count_df violin(:Count)

#-

@df count_df boxplot!([1.0], :Count, bar_width=0.1)

# #### Exercise 2
#
# Do a histogram and a density plot of the variable `:Count`.
# Hint: Use `normalize=true`

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
macro newton_tr_trace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(gr)
                dt["h(x)"] = copy(H)
                dt["delta"] = copy(delta)
                dt["interior"] = interior
                dt["hard case"] = hard_case
                dt["reached_subproblem_solution"] = reached_subproblem_solution
                dt["lambda"] = lambda
            end
            grnorm = norm(gr, Inf)
            update!(tr,
                    iteration,
                    f_x,
                    grnorm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

# Check whether we are in the "hard case".
#
# Args:
#  H_eigv: The eigenvalues of H, low to high
#  qg: The inner product of the eigenvalues and the gradient in the same order
#
# Returns:
#  hard_case: Whether it is a candidate for the hard case
#  lambda_1_multiplicity: The number of times the lowest eigenvalue is repeated,
#                         which is only correct if hard_case is true.
function check_hard_case_candidate(H_eigv, qg)
    @assert length(H_eigv) == length(qg)
    if H_eigv[1] >= 0
        # The hard case is only when the smallest eigenvalue is negative.
        return false, 1
    end
    hard_case = true
    lambda_index = 1
    hard_case_check_done = false
    while !hard_case_check_done
        if lambda_index > length(H_eigv)
            hard_case_check_done = true
        elseif abs(H_eigv[1] - H_eigv[lambda_index]) > 1e-10
            # The eigenvalues are reported in order.
            hard_case_check_done = true
        else
            if abs(qg[lambda_index]) > 1e-10
                hard_case_check_done = true
                hard_case = false
            end
            lambda_index += 1
        end
    end

    hard_case, lambda_index - 1
end

# Choose a point in the trust region for the next step using
# the interative (nearly exact) method of section 4.3 of Nocedal and Wright.
# This is appropriate for Hessians that you factorize quickly.
#
# Args:
#  gr: The gradient
#  H:  The Hessian
#  delta:  The trust region size, ||s|| <= delta
#  s: Memory allocated for the step size, updated in place
#  tolerance: The convergence tolerance for root finding
#  max_iters: The maximum number of root finding iterations
#
# Returns:
#  m - The numeric value of the quadratic minimization.
#  interior - A boolean indicating whether the solution was interior
#  lambda - The chosen regularizing quantity
#  hard_case - Whether or not it was a "hard case" as described by N&W
#  reached_solution - Whether or not a solution was reached (as opposed to
#      terminating early due to max_iters)
function solve_tr_subproblem!{T}(gr::Vector{T},
                                 H::Matrix{T},
                                 delta::T,
                                 s::Vector{T};
                                 tolerance=1e-10,
                                 max_iters=5)
    n = length(gr)
    delta_sq = delta^2

    @assert n == length(s)
    @assert (n, n) == size(H)
    @assert max_iters >= 1

    # Note that currently the eigenvalues are only sorted if H is perfectly
    # symmetric.  (Julia issue #17093)
    H_eig = eigfact(Symmetric(H))
    min_H_ev, max_H_ev = H_eig[:values][1], H_eig[:values][n]
    H_ridged = copy(H)

    # Cache the inner products between the eigenvectors and the gradient.
    qg = Array(T, n)
    for i=1:n
        qg[i] = vecdot(H_eig[:vectors][:, i], gr)
    end

    # Function 4.39 in N&W
    function p_sq_norm(lambda, min_i)
        p_sum = 0.
        for i = min_i:n
            p_sum += qg[i]^2 / (lambda + H_eig[:values][i])^2
        end
        p_sum
    end

    # These values describe the outcome of the subproblem.  They will be
    # set below and returned at the end.
    interior = true
    hard_case = false
    reached_solution = true

    if min_H_ev >= 1e-8 && p_sq_norm(0.0, 1) <= delta_sq
        # No shrinkage is necessary: -(H \ gr) is the minimizer
        interior = true
        reached_solution = true
        s[:] = -(H_eig[:vectors] ./ H_eig[:values]') * H_eig[:vectors]' * gr
        lambda = 0.0
    else
        interior = false

        # The hard case is when the gradient is orthogonal to all
        # eigenvectors associated with the lowest eigenvalue.
        hard_case_candidate, min_H_ev_multiplicity =
            check_hard_case_candidate(H_eig[:values], qg)

        # Solutions smaller than this lower bound on lambda are not allowed:
        # they don't ridge H enough to make H_ridge PSD.
        lambda_lb = -min_H_ev + max(1e-4, abs(min_H_ev) * 1e-4)
        lambda = lambda_lb

        hard_case = false
        if hard_case_candidate
            # The "hard case". lambda is taken to be -min_H_ev and we only need
            # to find a multiple of an orthogonal eigenvector that lands the
            # iterate on the boundary.

            # Formula 4.45 in N&W
            p_lambda2 = p_sq_norm(lambda, min_H_ev_multiplicity + 1)
            if p_lambda2 > delta_sq
                # Then we can simply solve using root finding.
                # Set a starting point between the minimum and largest eigenvalues.
                lambda = lambda_lb + 0.01 * (max_H_ev - lambda_lb)
            else
                hard_case = true
                reached_solution = true

                tau = sqrt(delta_sq - p_lambda2)

                # I don't think it matters which eigenvector we pick so take
                # the first.
                for i=1:n
                    s[i] = tau * H_eig[:vectors][i, 1]
                    for k=(min_H_ev_multiplicity + 1):n
                        s[i] = s[i] +
                               qg[k] * H_eig[:vectors][i, k] / (H_eig[:values][k] + lambda)
                    end
                end
            end
        end

        if !hard_case
            # Algorithim 4.3 of N&W, with s insted of p_l for consistency with
            # Optim.jl

            for i=1:n
                H_ridged[i, i] = H[i, i] + lambda
            end

            reached_solution = false
            for iter in 1:max_iters
                lambda_previous = lambda

                # Version 0.5 requires an exactly symmetric matrix, but
                # version 0.4 does not have this function signature for chol().
                R = VERSION < v"0.5-" ? chol(H_ridged): chol(Hermitian(H_ridged))
                s[:] = -R \ (R' \ gr)
                q_l = R' \ s
                norm2_s = vecdot(s, s)
                lambda_update = norm2_s * (sqrt(norm2_s) - delta) / (delta * vecdot(q_l, q_l))
                lambda += lambda_update

                # Check that lambda is not less than lambda_lb, and if so, go
                # half the way to lambda_lb.
                if lambda < (lambda_lb + 1e-8)
                    lambda = 0.5 * (lambda_previous - lambda_lb) + lambda_lb
                end

                for i=1:n
                    H_ridged[i, i] = H[i, i] + lambda
                end

                if abs(lambda - lambda_previous) < tolerance
                    reached_solution = true
                    break
                end
            end
        end
    end

    m = vecdot(gr, s) + 0.5 * vecdot(s, H_ridged * s)

    return m, interior, lambda, hard_case, reached_solution
end

immutable NewtonTrustRegion{T <: Real} <: Optimizer
    initial_delta::T
    delta_hat::T
    eta::T
    rho_lower::T
    rho_upper::T
end

NewtonTrustRegion(; initial_delta::Real = 1.0,
                    delta_hat::Real = 100.0,
                    eta::Real = 0.1,
                    rho_lower::Real = 0.25,
                    rho_upper::Real = 0.75) =
  NewtonTrustRegion(initial_delta, delta_hat, eta, rho_lower, rho_upper)

function optimize{T}(d::TwiceDifferentiableFunction,
                     initial_x::Vector{T},
                     mo::NewtonTrustRegion,
                     o::OptimizationOptions)

    @assert(mo.delta_hat > 0, "delta_hat must be strictly positive")
    @assert(0 < mo.initial_delta < mo.delta_hat, "delta must be in (0, delta_hat)")
    @assert(0 <= mo.eta < mo.rho_lower, "eta must be in [0, rho_lower)")
    @assert(mo.rho_lower < mo.rho_upper, "must have rho_lower < rho_upper")
    @assert(mo.rho_lower >= 0.)

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(initial_x), copy(initial_x)

    # Count the total number of iterations
    iteration = 0

    # Track calls to function and gradient
    f_calls, g_calls = 0, 0

    # Count number of parameters
    n = length(x)

    # Maintain current gradient in gr
    gr = Array(T, n)

    # The current search direction
    s = Array(T, n)

    # Store f(x), the function value, in f_x
    f_x_previous, f_x = NaN, d.fg!(x, gr)

    # We need to store the previous gradient in case we reject a step.
    gr_previous = copy(gr)

    f_calls, g_calls = f_calls + 1, g_calls + 1

    # Store the hessian in H
    H = Array(T, n, n)
    d.h!(x, H)

    # Keep track of trust region sizes
    delta = copy(mo.initial_delta)

    # Record attributes of the subproblem in the trace.
    hard_case = false
    reached_subproblem_solution = true
    interior = true
    lambda = NaN

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @newton_tr_trace

    # Assess multiple types of convergence
    x_converged, f_converged = false, false
    g_converged = vecnorm(gr, Inf) < o.g_tol

    # Iterate until convergence
    converged = g_converged
    while !converged && iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        # Find the next step direction.
        m, interior, lambda, hard_case, reached_subproblem_solution =
            solve_tr_subproblem!(gr, H, delta, s)

        # Maintain a record of previous position
        copy!(x_previous, x)

        # Update current position
        for i in 1:n
            @inbounds x[i] = x[i] + s[i]
        end

        # Update the function value and gradient
        copy!(gr_previous, gr)
        f_x_previous, f_x = f_x, d.fg!(x, gr)
        f_calls, g_calls = f_calls + 1, g_calls + 1

        # Update the trust region size based on the discrepancy between
        # the predicted and actual function values.  (Algorithm 4.1 in N&W)
        f_x_diff = f_x_previous - f_x
        if abs(m) <= eps(T)
            # This should only happen when the step is very small, in which case
            # we should accept the step and assess_convergence().
            rho = 1.0
        elseif m > 0
            # This can happen if the trust region radius is too large and the
            # Hessian is not positive definite.  We should shrink the trust
            # region.
            rho = mo.rho_lower - 1.0
        else
            rho = f_x_diff / (0 - m)
        end

        if rho < mo.rho_lower
            delta *= 0.25
        elseif (rho > mo.rho_upper) && (!interior)
            delta = min(2 * delta, mo.delta_hat)
        else
            # else leave delta unchanged.
        end

        if rho > mo.eta
            # Accept the point and check convergence
            x_converged,
            f_converged,
            g_converged,
            converged = assess_convergence(x,
                                           x_previous,
                                           f_x,
                                           f_x_previous,
                                           gr,
                                           o.x_tol,
                                           o.f_tol,
                                           o.g_tol)
            if !converged
                # Only compute the next Hessian if we haven't converged
                d.h!(x, H)
            end
        else
            # The improvement is too small and we won't take it.

            # If you reject an interior solution, make sure that the next
            # delta is smaller than the current step.  Otherwise you waste
            # steps reducing delta by constant factors while each solution
            # will be the same.
            x_diff = x - x_previous
            delta = 0.25 * sqrt(vecdot(x_diff, x_diff))

            f_x = f_x_previous
            copy!(x, x_previous)
            copy!(gr, gr_previous)

        end

        @newton_tr_trace
    end

    return MultivariateOptimizationResults("Newton's Method with Trust Region",
                                           initial_x,
                                           x,
                                           Float64(f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           o.x_tol,
                                           f_converged,
                                           o.f_tol,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           g_calls)
end


text:
"""
    row_mul(v::AbstractArray, A::AbstractMatrix)
Multiply all the rows of a matrix `A` by some vector `v` element by element.

"""
function row_mul(v::AbstractVector, A::AbstractMatrix)

    #Error
    if length(v) != size(A,2)
        return ErrorException("The number of colums of the matrix ($(size(A,2))) must be equal to the length of the vector ($(length(v)))" )
    end

    B = similar(A)

    @simd for i in 1:size(A,1)
        @simd for j in 1:size(A,2)
            @inbounds B[i,j] = v[j]*A[i,j]
        end
    end

    return B
end

"""
    rand_exclusive(v::AbstractArray,n::Integer)
Returns an array with `n` different random elements of the array `v`.
"""
function rand_exclusive(v::AbstractVector,n::Integer)

    #Error checks
    if n>length(v)
        return ErrorException("Length of the array `v` ($(length(v))) must be less or equal than `n` ($n)")
    end

    a = Array{Number}(undef,n)
    a[1] = rand(v)

    for i in 2:n
        tmp = rand(v)
        while tmp in a[1:i-1]
            tmp = rand(v)
        end
        a[i] = tmp
    end
    return a
end

function best_individue(fob::Function,A::AbstractMatrix)
    vals = findmin(@inbounds [fob(A[i,:]...) for i in 1:size(A,1)])
    return (vals[1],A[vals[2],:])
end

"""
    rand_cauchy_trunc(μ::Real, c::Real)
Returns a random number with a *Cauchy's distribution* with **mean `μ`** and truncated between (0,1]. If the value calculated exceds 1, the returned value is 1.
If the value is less or equal to zero, the value es calculated again.
"""
function rand_cauchy_trunc(μ::Real, c::Real)
    r = c*tan(pi*(rand()-0.5)) + μ
    if r > 1
        return 1
    elseif r <= 0
        return rand_cauchy_trunc(μ::Real, c::Real)
    end
    return r
end

"""
    rand_normal_trunc(μ::Real, σ::Real)
Returns a random number with a *normal* distribution with **mean `μ`** and **standard deviation `σ`** truncated between [0,1].
"""
function rand_normal_trunc(μ::Real, σ::Real)
    r = σ*randn() + μ
    if r > 1
        return 1
    elseif r < 0
        return 0
    end
    return r
end

"""
    Lehmer_mean(X::AbstractArray,p::Real)
Return the **Lehmer mean** of degree `p` of an array `X`.
"""
function Lehmer_mean(X::AbstractVector, p::Real)
    return sum(X.^p)/sum(X.^(p-1))
end

"""
    pbest(fob::Function, A::AbstractMatrix, p::Real)
Return the top 100*p% best individues of a population evaluated in a function `fob`.
"""
function pbest(fob::Function, A::AbstractMatrix, p::Real)
    len = size(A,1)
    plen = floor(Integer,p*len)
    ppop = sortperm(@inbounds [fob(A[i,:]...) for i in 1:len])[1:plen]
    return A[ppop,:]
end


text:
#=

This file is the entry-point for the vSmartMOM module.

It includes this module's source files and exports the relevant keywords.

=#

module CoreRT

using Interpolations               # For interpolating the vmr's
using LinearAlgebra                # For linear algebra routines
using ProgressMeter                # Showing progress in for loops
using Distributions                # Distributions of aerosols
using ..Scattering                 # Use scattering module
using ..Absorption                 # Use absorption module
using ...vSmartMOM         # Use parent RadiativeTransfer module
using ...Architectures             # Use Architectures module

using CUDA                         # GPU CuArrays and functions
using KernelAbstractions           # Abstracting code for CPU/GPU
using KernelAbstractions.Extras
using CUDAKernels

using FastGaussQuadrature          # Computes quadrature points (Gauss, legendre, Radau,...)
using TimerOutputs                 # For timing sections of the code
using Parameters                   # For keyword arguments in structs
using DocStringExtensions          # For documenting
using YAML                         # For reading properties files
using ForwardDiff                  # Automatic Differentiation
using NNlib                        # For batched multiplications
import NNlib.batched_mul           # Required to overwrite batched_mul for Duals
using NCDatasets                   # For loading absco lookup tables

import Base.show                   # For overloading show for custom types

# More threads in LA wasn't really helpful, can be turned off here:
# LinearAlgebra.BLAS.set_num_threads(1)

# Constants and Types
include("constants.jl")                     # Scientific constants
include("types.jl")                         # All custom types for this module

# Solvers
include("CoreKernel/elemental.jl")          # Elemental
include("CoreKernel/doubling.jl")           # Doubling
include("CoreKernel/interaction.jl")        # Interaction
include("CoreKernel/rt_kernel.jl")          # Handle Core RT (Elemental/Doubling/Interaction)

include("postprocessing_vza.jl")            # Postprocess (Azimuthal Weighting)
include("rt_run.jl")                        # Starting point for RT

# GPU
include("gpu_batched.jl")                   # Batched operations

# Utilities / Helper Functions
include("atmo_prof.jl")                     # Helper Functions for Handling Atmospheric Profiles
include("rt_helper_functions.jl")           # Miscellaneous Utility Functions
include("rt_set_streams.jl")                # Set streams before RT
include("parameters_from_yaml.jl")          # Loading in parameters from YAML file
include("model_from_parameters.jl")         # Converting parameters to derived model attributes
include("show_utils.jl")                    # Pretty-printing objects

# Surfaces
include("lambertian_surface.jl")            # Lambertian Surface

# Functions to export
export parameters_from_yaml,                # Getting parameters from a file
       model_from_parameters,               # Converting the parameters to model
       rt_run,                              # Run the RT code
       default_parameters                   # Set of default parameters

# Export types to show easily
export GaussQuadFullSphere, LambertianSurfaceScalar, LambertianSurfaceSpectrum

end


text:
function render_gui(window, clear_color; fps)
    CImGui.Render()
    GLFW.MakeContextCurrent(window)
    display_w, display_h = GLFW.GetFramebufferSize(window)
    glViewport(0, 0, display_w, display_h)
    glClearColor(clear_color...)
    glClear(GL_COLOR_BUFFER_BIT)
    ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())
    GLFW.MakeContextCurrent(window)
    GLFW.SwapBuffers(window)
    sleep(1/fps) # hacky update rate
end

function destroy_gui(ctx, window)
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end

function init_gui()
    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1300, 800, "Lab Instruments")
    @assert window != C_NULL

    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()
    CImGui.StyleColorsClassic()
    fonts_dir = joinpath(@__DIR__, "..", "fonts")
    fonts = CImGui.GetIO().Fonts
    CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)
    # CImGui.SetWindowFontScale(Cfloat(2.0)) # crashes julia
    return window, ctx
end

function ShowMenuWindow(psu_conf, dmm_conf, fgen_conf, scope_conf, relays_conf)
	CImGui.Begin("Menu")
		@c CImGui.Checkbox("psu", &psu_conf.active)
		CImGui.SameLine(), @c CImGui.Checkbox("dmm", &dmm_conf.active)
		CImGui.SameLine(), @c CImGui.Checkbox("fgen", &fgen_conf.active)
		CImGui.SameLine(), @c CImGui.Checkbox("scope", &scope_conf.active)
		CImGui.SameLine(), @c CImGui.Checkbox("relays", &relays_conf.active)
	CImGui.End()
	return psu_conf, dmm_conf, fgen_conf, scope_conf, relays_conf
end

#function ShowMenuWindow(psu_conf)
#	CImGui.Begin("Menu")
#		@c CImGui.Checkbox("psu", &psu_conf.active)
#	CImGui.End()
#	return psu_conf
#end

function draw_toggle_button(button)
	# hack for toggle button
	if button.state == false
		CImGui.PushID(1)
		CImGui.PushStyleColor(CImGui.ImGuiCol_Button, CImGui.HSV(button.off_color...))
		CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonHovered, CImGui.HSV(button.hover_color...))
		CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonActive, CImGui.HSV(button.off_color...))
	else
		CImGui.PushID(1)
		CImGui.PushStyleColor(CImGui.ImGuiCol_Button, CImGui.HSV(button.on_color...))
		CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonHovered, CImGui.HSV(button.hover_color...))
		CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonActive, CImGui.HSV(button.on_color...))
	end
	#pressed = CImGui.Button(button.name, CImGui.ImVec2(30,20))  # 30px width, 20px height
	pressed = CImGui.Button(button.name)
	pressed && toggle!(button)

	CImGui.PopStyleColor(3)
	CImGui.PopID()
	return button, pressed
end

# PSU
function draw_psu_info(args...)
	CImGui.Text(args[1])
	for arg in args[2:end]
		CImGui.SameLine(), CImGui.Text(arg)
	end
	return nothing
end

# DMM
function draw_dmm_info(args...)
	CImGui.Text(args[1])
	for arg in args[2:end]
		CImGui.SameLine(), CImGui.Text(arg)
	end
	return nothing
end

# Scope
function draw_scope_info(args...)
	CImGui.Text(args[1])
	for arg in args[2:end]
		CImGui.SameLine(), CImGui.Text(arg)
	end
	return nothing
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
#=
This script shows how a simple MLP net may be used
for regression. It shows how data in memory may be
used for training and evaluation, and how to obtain
the predictions from the trained net.
=#
using MXNet
using Distributions
using PyPlot

# data generating process
generate_inputs(mean, var, size) = rand(MvNormal(mean, var), size)
output(data) = sin(data[1,:]).*sin(data[2,:])./(data[1,:].*data[2,:])

# create training and evaluation data sets
mean=[0.0;0.0]
var=[1.0 0.0;0.0 1.0]
samplesize  = 5000
TrainInput = generate_inputs(mean, var, samplesize)
TrainOutput = output(TrainInput)
ValidationInput = generate_inputs(mean, var, samplesize)
ValidationOutput = output(ValidationInput)

# how to set up data providers using data in memory
batchsize = 100 # can adjust this later, but must be defined now for next line
trainprovider = mx.ArrayDataProvider(:data => TrainInput, batch_size=batchsize, shuffle=true, :label => TrainOutput)
evalprovider = mx.ArrayDataProvider(:data => ValidationInput, batch_size=batchsize, shuffle=true, :label => ValidationOutput)

# create a two hidden layer MPL: try varying num_hidden, and change tanh to relu,
# or add/remove a layer
data = mx.Variable(:data)
label = mx.Variable(:label)
net  = @mx.chain    mx.FullyConnected(data = data, num_hidden=10) =>
                    mx.Activation(act_type=:tanh) =>
                    mx.FullyConnected(num_hidden=3) =>
                    mx.Activation(act_type=:tanh) =>
                    mx.FullyConnected(num_hidden=1)

# squared error loss is appropriate for regression, don't change
cost = mx.LinearRegressionOutput(data = net, label=label)

# final model definition, don't change, except if using gpu
model = mx.FeedForward(cost, context=mx.cpu())

# set up the optimizer: select one, explore parameters, if desired
#optimizer = mx.SGD(lr=0.01, momentum=0.9, weight_decay=0.00001)
optimizer = mx.ADAM()

# train, reporting loss for training and evaluation sets
# initial training with small batch size, to get to a good neighborhood
batchsize = 100
mx.fit(model, optimizer, initializer=mx.NormalInitializer(0.0,0.1), eval_metric=mx.MSE(), trainprovider, eval_data=evalprovider, n_epoch = 500)
# more training with the full sample
batchsize = samplesize
mx.fit(model, optimizer, eval_metric=mx.MSE(), trainprovider, eval_data=evalprovider, n_epoch = 500)

# obtain predictions
plotprovider = mx.ArrayDataProvider(:data => ValidationInput, :label => ValidationOutput)
fit = mx.predict(model, plotprovider)
plot(ValidationOutput,fit',".")
xlabel("true")
ylabel("predicted")
title("outputs: true versus predicted. 45º line is what we hope for")


text:
# ---
# title: 1281. Subtract the Product and Sum of Digits of an Integer
# id: problem1281
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/subtract-the-product-and-sum-of-digits-of-an-integer/description/>
# hidden: true
# ---
#
# Given an integer number `n`, return the difference between the product of its
# digits and the sum of its digits.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 234
#     Output: 15
#     Explanation:
#     Product of digits = 2 * 3 * 4 = 24
#     Sum of digits = 2 + 3 + 4 = 9
#     Result = 24 - 9 = 15
#
#
# **Example 2:**
#
#
#
#     Input: n = 4421
#     Output: 21
#     Explanation: Product of digits = 4 * 4 * 2 * 1 = 32
#     Sum of digits = 4 + 4 + 2 + 1 = 11
#     Result = 32 - 11 = 21
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
const U = UInt8

"""
    encode(x) -> Vector{$U}

Encode some [ETF](http://erlang.org/doc/apps/erts/erl_ext_dist.html) data.
"""
encode(x) = U[VERSION_MAGIC; _encode(x)]

_encode(x) = error("Unknown ETF type $(typeof(x))")

# Integers

# Note: Encoding with precision does not include the tag byte.
_encode(x::Integer, p::Integer) = reverse(reinterpret(U, [x]))[end-p+1:end]
function _encode(x::Integer)
     return if 0 <= x <= 255
         U[97, x]
     elseif typemin(Int32) <= x <= typemax(Int32)
         U[98; _encode(x, 4)]
     else
         n = Int(x < 0)
         digs = digits(U, x; base=256)
         if -(BigInt(256)^255) < x < BigInt(256)^255
             U[110; round(U, log(256, abs(x))) + 1; n; digs]
         else
             U[111; _encode(round(UInt32, log(256, abs(x))) + 1, 4); n; digs]
         end
     end
end

# Floats

_encode(x::AbstractFloat) = U[V, 70, reverse(reinterpret(U, [x]))...]

# Tuples

function _encode(x::Tuple)
    len = length(x)
    tag = len <= typemax(U) ? 104 : 105
    return U[tag; _encode(len)[2:end]; vcat(_encode.(x)...)]
end

# Atoms

function _encode(x::Symbol)
    # This one is a weird, ATOM_EXT is deprecated but latest OTP still uses it.
    # Also the length limit for atoms is supposedly 255, which would fit in SMALL_ATOM_EXT.
    # However atoms of any length up to 255 are encoded as ATOM_EXT.
    # We're just going to do what the spec actually says, except for the atom length limit.
    s = Vector{U}(string(x))
    len = length(s)
    return len <= typemax(U) ? U[119; len; s] : U[118; _encode(len, 2); s]
end
_encode(x::Bool) = _encode(Symbol(x))
_encode(::Nothing) = _encode(:nil)

# Lists

function _encode(x::AbstractArray)
    return isempty(x) ? U[106] : U[108; _encode(length(x), 4); _encode.(x)...; 106]
end
_encode(x::Vector{U}) = U[107; _encode(length(x), 2); x]

# Maps

function _encode(x::AbstractDict)
    return U[
        116;
        _encode(length(x), 4);
        vcat(_encode.(vcat([[p.first, p.second] for p in x]...))...);
    ]
end

# Strings

_encode(x::AbstractString) = U[109; _encode(length(x), 4); Vector{U}(x)]


text:
using NBodyIPs, StaticArrays, JuLIP, BenchmarkTools

# using NBodyIPs.Polys: invariants, invariants_d, invariants_ed
# using JuLIP.Potentials: evaluate, evaluate_d
const BLI = NBodyIPs.BLInvariants

println("------------------------------------------")
println("  [1] Bond-Length Invariants Evaluation:")
println("------------------------------------------")
for r in [ (@SVector [1.234]),
           (@SVector rand(3)),
           (@SVector rand(6)),
           (@SVector rand(10)) ]
   println("dim = $(length(r))")
   print("     invariants: ")
   @btime BLI.invariants($r)
   print("   invariants_d: ")
   @btime BLI.invariants_d($r)
   print("  invariants_ed: ")
   @btime BLI.invariants_ed($r)
end

const BAI = NBodyIPs.BAInvariants

println("------------------------------------------")
println("  [2] Bond-Angle Invariants Evaluation:")
println("------------------------------------------")
for x in [ (@SVector [1.234]),
           (@SVector rand(3)),
           (@SVector rand(6)),
           (@SVector rand(10)) ]
   println("dim = $(length(x))")
   print("     invariants: ")
   @btime BAI.invariants($x)
   print("   invariants_d: ")
   @btime BAI.invariants_d($x)
   print("  invariants_ed: ")
   @btime BAI.invariants_ed($x)
end

# println("----------------------------------------")
# println("   [2] `NBody` Evaluation")
# println("----------------------------------------")
#
# r0 = rnn(:Cu)
# at = rattle!(bulk(:Cu, cubic=true) * 2, 0.02)
# @show length(at)
#
# TRANSFORM = let r0 = r0
#    (@analytic r -> (r0/r)^3)
# end
# rcut3 = 3.1 * r0
# D3 = Dictionary(TRANSFORM, (:cos, 0.66*rcut3, rcut3) )
# rcut4 = 2.1 * r0
# D4 = Dictionary(TRANSFORM, (:cos, 0.66*rcut4, rcut4) )
# rcut5 = 1.5 * r0
# D5 = Dictionary(TRANSFORM, (:cos, 0.66*rcut5, rcut5) )
#
# n = 10
# println("[1] Quick profiling for a N-body with $n basis functions")
# r = 1.0 + rand(SVector{3, Float64})
# V3 = NBody( [tuple([rand(0:4, 3);0]...) for n = 1:n], 1.0+rand(n), D3 )
# print("     V3: "); @btime evaluate($V3, $r)
# print("  @D V3: "); @btime evaluate_d($V3, $r)
# # print("  nlist: "); @btime neighbourlist($at, $rcut3)
# # print(" energy: "); @btime energy($V3, $at)
# # print(" forces: "); @btime forces($V3, $at)
#
# function rand_tuple_4b()
#    t = rand(0:5, 7)
#    t[rand(1:6, 5)] = 0
#    return tuple(t...)
# end
#
# r = 1.0 + rand(SVector{6, Float64})
# V4 = NBody( [rand_tuple_4b() for _ = 1:n], 1.0+rand(n), D4 )
# print("     V4: "); @btime evaluate($V4, $r)
# print("  @D V4: "); @btime evaluate_d($V4, $r)
# # print("  nlist: "); @btime neighbourlist(at, rcut4)
# # print(" energy: "); @btime energy(V4, at)
# # print(" forces: "); @btime forces(V4, at)
#
# function rand_tuple_5b()
#    t = rand(0:6, 11)
#    t[rand(1:10, 9)] = 0
#    return tuple(t...)
# end
#
# r = 1.0 + rand(SVector{10, Float64})
# V5 = NBody( [rand_tuple_5b() for _ = 1:n], 1.0+rand(n), D5 )
# print("     V5: "); @btime evaluate($V5, $r)
# print("  @D V5: "); @btime evaluate_d($V5, $r)
# # print("  nlist: "); @btime neighbourlist(at, rcut4)
# # print(" energy: "); @btime energy(V4, at)
# # print(" forces: "); @btime forces(V4, at)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using Dubins
using Memento

using Test

setlevel!(getlogger(Dubins), "error")

@testset "Dubins" begin
    include("test_api.jl")
    include("test_path.jl")
end


text:
function solve(connections, use_twice=false)

    paths = Vector()
    visited = Set()
    function dfs(curr, can_still_use_twice=false, path=Vector())
        key = curr * (curr in visited ? "2" : "")
        push!(path, curr)
        if lowercase(curr) == curr
            push!(visited, key)
        end
        if curr == "end"
            push!(paths, join(path, ','))
        else
            for connection in connections[curr]
                if connection ∉ visited
                    dfs(connection, can_still_use_twice, copy(path))
                elseif can_still_use_twice && connection ∉ ["start", "end"]
                    dfs(connection, false, copy(path))
                end
            end
        end
        delete!(visited, key)
    end

    dfs("start", use_twice)
    return length(paths)
end

function main()
    connections = Dict()
    for line in readlines()
        from, to = split(line, '-')
        push!(get!(connections, from, Vector()), to)
        push!(get!(connections, to, Vector()), from)
    end
    println(solve(connections))
    println(solve(connections, true))
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
@doc raw"""
Constructs constraint from rating data and related variable tuple.

# Constraint

``` var1[r[1], t]^2 + var2[r[1], t]^2 <= r[2]^2 ```

where r in rating data and t in time steps.

# LaTeX

`` x_1^2 + x_2^2 \leq r^2 ``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* rating_data::Vector{Tuple{String, Float64}} : rating data name (1) and value (2)
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol} : the names of the variables
- : var_names[1] : var1
- : var_names[2] : var2
"""
function rating_constraint!(
    optimization_container::OptimizationContainer,
    rating_data::Vector{Tuple{String, Float64}},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol},
)
    time_steps = model_time_steps(optimization_container)
    var1 = get_variable(optimization_container, var_names[1])
    var2 = get_variable(optimization_container, var_names[2])
    add_cons_container!(
        optimization_container,
        cons_name,
        [r[1] for r in rating_data],
        time_steps,
    )
    constraint = get_constraint(optimization_container, cons_name)

    for r in rating_data
        for t in time_steps
            constraint[r[1], t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                var1[r[1], t]^2 + var2[r[1], t]^2 <= r[2]^2
            )
        end
    end

    return
end


text:
function next_chebyshev(previous::Sum_Of_Functions, two_previous::Sum_Of_Functions, dim_name::Symbol)
    return PE_Function(2.0,Dict{Symbol,PE_Unit}(dim_name => PE_Unit(0.0,0.0,1))) * previous - two_previous
end

function get_first_kind_Chebyshevs(num::Int, dim_name::Symbol)
    if num < 1
        error("Select a positive number of Chebyshevs to return")
    else
        first_kind_chebyshevs = Array{Sum_Of_Functions}(undef, num)
        first_kind_chebyshevs[1] = Sum_Of_Functions([PE_Function(1.0,Dict{Symbol,PE_Unit}(dim_name => PE_Unit(0.0,0.0,0)))])
        if num > 1
            first_kind_chebyshevs[2] = Sum_Of_Functions([PE_Function(1.0,Dict{Symbol,PE_Unit}(dim_name => PE_Unit(0.0,0.0,1)))])
        end
        for n in 3:num
            first_kind_chebyshevs[n] = next_chebyshev(first_kind_chebyshevs[n-1], first_kind_chebyshevs[n-2], dim_name)
        end
        return first_kind_chebyshevs
    end
end
function get_second_kind_Chebyshevs(num::Int, dim_name::Symbol)
    if num < 1
        error("Select a positive number of Chebyshevs to return")
    else
        second_kind_chebyshevs = Array{Sum_Of_Functions}(undef, num)
        second_kind_chebyshevs[1] = Sum_Of_Functions([PE_Function(1.0,Dict{Symbol,PE_Unit}(dim_name => PE_Unit(0.0,0.0,0)))])
        if num > 1
            second_kind_chebyshevs[2] = Sum_Of_Functions([PE_Function(2.0,Dict{Symbol,PE_Unit}(dim_name => PE_Unit(0.0,0.0,1)))])
        end
        for n in 3:num
            second_kind_chebyshevs[n] = next_chebyshev(second_kind_chebyshevs[n-1], second_kind_chebyshevs[n-2], dim_name)
        end
        return second_kind_chebyshevs
    end
end

"""
    get_chevyshevs_up_to(num::Int, first_kind::Bool = true; dim_name::Symbol = default_symbol)
    Output all chebyshev polynomials up to degree num.
"""
function get_chevyshevs_up_to(num::Int, first_kind::Bool = true; dim_name::Symbol = default_symbol)
    if first_kind
        return get_first_kind_Chebyshevs(num,dim_name)
    else
        return get_second_kind_Chebyshevs(num,dim_name)
    end
end


text:
module Containers

include("module_common.jl")

const _libcloud_container_types = PyCall.PyNULL()
const _libcloud_container_providers = PyCall.PyNULL()
const _libcloud_container_base = PyCall.PyNULL()

immutable ContainerDriver
    o::PyObject

    function ContainerDriver(provider::Compat.String, args...; kwargs...)
        cls = _libcloud_container_providers[:get_driver](provider)
        ContainerDriver(cls(args...; kwargs...))
    end
    ContainerDriver(driver::PyObject) = new(driver)
end
show(io::IO, c::ContainerDriver) = print(io, c.o[:__str__]())

immutable ContainerImage
    o::PyObject

    id::Nullable{Compat.String}
    name::Compat.String
    path::Compat.String
    version::Compat.String
    extra::Dict
    driver::ContainerDriver

    function ContainerImage(o::PyObject)
        new(o, o[:id], o[:name], o[:path], o[:version], o[:extra], ContainerDriver(o[:driver]))
    end
end
PyObject(o::ContainerImage) = o.o
convert(::Type{ContainerImage}, o::PyObject) =  ContainerImage(o)
show(io::IO, o::ContainerImage) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_container_base, "ContainerImage") ContainerImage

immutable ContainerCluster
    o::PyObject

    id::Nullable{Compat.String}
    name::Compat.String
    extra::Dict
    driver::ContainerDriver

    function ContainerCluster(o::PyObject)
        new(o, o[:id], o[:name], o[:extra], ContainerDriver(o[:driver]))
    end
end
PyObject(o::ContainerCluster) = o.o
convert(::Type{ContainerCluster}, o::PyObject) =  ContainerCluster(o)
show(io::IO, o::ContainerCluster) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_container_base, "ContainerCluster") ContainerCluster

immutable ClusterLocation
    o::PyObject

    id::Compat.String
    name::Compat.String
    country::Compat.String
    driver::ContainerDriver

    function ClusterLocation(o::PyObject)
        new(o, o[:id], o[:name], o[:country], ContainerDriver(o[:driver]))
    end
end
PyObject(o::ClusterLocation) = o.o
convert(::Type{ClusterLocation}, o::PyObject) =  ClusterLocation(o)
show(io::IO, o::ClusterLocation) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_container_base, "ClusterLocation") ClusterLocation

immutable Container
    o::PyObject

    id::Nullable{Compat.String}
    name::Compat.String
    image::ContainerImage
    state::Compat.String
    ip_addresses::Array
    extra::Dict
    driver::ContainerDriver

    function Container(o::PyObject)
        new(o, o[:id], o[:name], o[:image], o[:state], o[:ip_addresses], o[:extra], ContainerDriver(o[:driver]))
    end
end
PyObject(o::Container) = o.o
convert(::Type{Container}, o::PyObject) =  Container(o)
show(io::IO, o::Container) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_container_base, "Container") Container

const _containerdriver_fns = [
    :install_image, :list_images, :list_containers, :deploy_container, :get_container,
    :start_container, :stop_container, :restart_container, :destroy_container,
    :list_locations, :create_cluster, :destroy_cluster, :list_clusters, :get_cluster
]

const _container_cluster_fns = [ :list_containers, :destroy ]
const _container_image_fns = [ :deploy ]
const _container_fns = [ :start, :stop, :restart, :destroy ]

for f in union(Set(_containerdriver_fns), Set(_base_driver_fns))
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_container_base, "ContainerDriver", $sf) $(f)(d::ContainerDriver, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
end

for f in _container_cluster_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_container_base, "ContainerCluster", $sf) $(f)(c::ContainerCluster, args...; kwargs...) = c.o[$(sf)](args...; kwargs...)
end

for f in _container_image_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_container_base, "ContainerImage", $sf) $(f)(c::ContainerImage, args...; kwargs...) = c.o[$(sf)](args...; kwargs...)
end

for f in _container_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_container_base, "Container", $sf) $(f)(c::Container, args...; kwargs...) = c.o[$(sf)](args...; kwargs...)
end

# Initialize container service
function __init__()
    copy!(_libcloud_container_types, pyimport("libcloud.container.types"))
    copy!(_libcloud_container_providers, pyimport("libcloud.container.providers"))
    copy!(_libcloud_container_base, pyimport("libcloud.container.base"))

    global const ContainerProvider = pywrap(_libcloud_container_types[:Provider])
    global const ContainerState = pywrap(_libcloud_container_types[:ContainerState])

    _map_types(_libcloud_container_base, (ContainerImage, ContainerCluster, ClusterLocation, Container))
end

# types
export ContainerProvider, ContainerDriver, ContainerState, ContainerImage, ContainerCluster, ClusterLocation, Container
# container driver functions
export install_image, list_images, list_containers, deploy_container, get_container,
       start_container, stop_container, restart_container, destroy_container,
       list_locations, create_cluster, destroy_cluster, list_clusters, get_cluster
# container cluster functions
export list_containers, destroy
# container image functions
export deploy
# container functions
export start, stop, restart, destroy

end # module Containers


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

using Markdown
using InteractiveUtils

# ╔═╡ 9661e304-39b3-11eb-1373-c3959d5b37ff
using Plots

# ╔═╡ 838d642e-399f-11eb-36f4-6fd348a4f84d
begin
	import DarkMode
	DarkMode.enable()
end

# ╔═╡ 4a5a6590-399e-11eb-0931-c737ec2e34e4
struct rainsimulator
	l::Float64 # path length
	h::Float64 # modelregion height
	vr::Float64 # rain velocity
 	vm::Float64 # agent velocity
	θ::Float64 # angle with respect to vertical
	ρR::Float64 # rain density
	a::Float64 # head diameter
	b::Float64 # agent height
	Δt::Float64 #timestep
end

# ╔═╡ 206e2cca-399f-11eb-2ee0-d9f9456ebbdb
mutable struct rain
	x_pos::Array
	y_pos::Array
end

# ╔═╡ 3e015876-39a3-11eb-21a1-d90460d77835
mutable struct agent
	x_pos::Float64
	y_pos::Float64
end

# ╔═╡ 2f4c2286-39ab-11eb-1c23-e998a265526a
mutable struct response
	time::Array
	cummulative_hits::Array
end

# ╔═╡ dcc75c1e-39ad-11eb-0b0b-13add6a0ec03

# ╔═╡ 6c5ea9ba-39a2-11eb-0fa2-71931a4498ee

# ╔═╡ ae4cc5ea-39a1-11eb-2079-69a20a3b9a45
function init_field(config::rainsimulator)
	density = config.ρR
	dim_x = config.l
	dim_y = config.h
	numofraindrops = convert(Int64,floor(dim_x*dim_y*density))
	rain(dim_x*rand(numofraindrops),dim_y*rand(numofraindrops))
end

# ╔═╡ 442d4544-39a7-11eb-385c-035da6732479
function update_rain!(raininstanz::rain,config::rainsimulator)
	θ = config.θ
	vr = config.vr
	Δt = config.Δt
	Δx = -cosd(θ)
	Δy = -sind(θ)
	raininstanz.x_pos = raininstanz.x_pos .+ (Δx * vr* Δt)
	raininstanz.y_pos = raininstanz.y_pos .+ (Δy * vr* Δt)
	for i in 1:length(raininstanz.x_pos)
		if raininstanz.x_pos[i] < 0
			raininstanz.x_pos[i] += config.l
		end
		if raininstanz.x_pos[i] > config.l
			raininstanz.x_pos[i] -= config.l
		end
		if raininstanz.y_pos[i] < 0
			raininstanz.y_pos[i] += config.h
		end
		if raininstanz.y_pos[i] > config.h
			raininstanz.y_pos[i] -= config.h
		end
	end
end

# ╔═╡ d7a666de-39a7-11eb-0537-ebd5786d023a
function update_agent!(agentinstanz::agent,config::rainsimulator)
	vm = config.vm
	Δt = config.Δt
	agentinstanz.x_pos += vm * Δt
end

# ╔═╡ 6f2bde14-39ab-11eb-2197-d1af40ec3dd8
function check_hit_drops!(
		responseinstanz::response,
		raininstanz::rain,
		agentinstanz::agent,
		config::rainsimulator)
	x_a = agentinstanz.x_pos
	y_a = 0.0
	x_b = agentinstanz.x_pos + config.a
	y_b = config.b
	y = raininstanz.y_pos
	push!(responseinstanz.cummulative_hits,0)
	for i in 1:length(raininstanz.x_pos)
		if (raininstanz.x_pos[i] > x_a) && (raininstanz.x_pos[i] < x_b)
			if (raininstanz.y_pos[i] > y_a) && (raininstanz.y_pos[i] < y_b)
				responseinstanz.cummulative_hits[end] += 1
				raininstanz.x_pos[i] = config.l*rand()
				raininstanz.y_pos[i] = config.h*rand()
			end
		end
	end
end

# ╔═╡ 68a7ee40-39ac-11eb-3b89-b9a85874f469
function maketimestep!(responseinstanz::response,
		raininstanz::rain,
		agentinstanz::agent,
		config::rainsimulator)
	update_rain!(raininstanz, config)
	update_agent!(agentinstanz, config)
	check_hit_drops!(responseinstanz,raininstanz,agentinstanz,config)
	push!(responseinstanz.time,responseinstanz.time[end] + config.Δt)
end

# ╔═╡ 266f8b74-39af-11eb-28fb-d7deda1a1a36
function run(config::rainsimulator)
	responseinstanz = response([0.0],[0])
	agentinstanz = agent(0,0)
	rainfield = init_field(config)
	n = ceil(config.l/(config.vm*config.Δt))
	print(n)
	for i in 1:n
		maketimestep!(responseinstanz,rainfield,agentinstanz, config)
	end
	responseinstanz
end

# ╔═╡ f5c27d76-39af-11eb-3e30-d37e2db01542
begin
	c1=rainsimulator(100, 100, 5.1, 1, 89, 1, 0.5, 2, 0.01)
	c2=rainsimulator(100, 100, 5.1, 5.1, 89, 1, 0.5, 2, 0.01)
	c3=rainsimulator(100, 100, 5.1, 10, 89, 1, 0.5, 2, 0.01)
	c4=rainsimulator(100, 100, 5.1, 1, -89, 1, 0.5, 2, 0.01)
	c5=rainsimulator(100, 100, 5.1, 5.1, -89, 1, 0.5, 2, 0.01)
	c6=rainsimulator(100, 100, 5.1, 10, -89, 1, 0.5, 2, 0.01)
end

# ╔═╡ 7cb9b56a-39b0-11eb-31a0-b98dcd35ad08
begin
	output1=run(c1)
	output2=run(c2)
	output3=run(c3)
	output4=run(c4)
	output5=run(c5)
	output6=run(c6)
end

# ╔═╡ 86fef7d2-68a0-11eb-1974-4d2695221bdf
output1.cummulative_hits

# ╔═╡ b41fdc84-39b3-11eb-0107-2fff3564a831

# ╔═╡ e814282e-39b3-11eb-0698-5b00781a15c0
begin
	pl = plot()
	plot!(output1.time,output1.cummulative_hits)
	plot!(output2.time,output2.cummulative_hits)
	plot!(output3.time,output3.cummulative_hits)
	pl
end

# ╔═╡ b9893c4a-39b5-11eb-042a-99e76c8f37f9
begin
	pl2 = plot()
	plot!(output1.time,cumsum(output1.cummulative_hits))
	plot!(output2.time,cumsum(output2.cummulative_hits))
	plot!(output3.time,cumsum(output3.cummulative_hits))
	plot!(output4.time,cumsum(output4.cummulative_hits))
	plot!(output5.time,cumsum(output5.cummulative_hits))
	plot!(output6.time,cumsum(output6.cummulative_hits))
	pl2
end

# ╔═╡ Cell order:
# ╠═838d642e-399f-11eb-36f4-6fd348a4f84d
# ╠═4a5a6590-399e-11eb-0931-c737ec2e34e4
# ╠═206e2cca-399f-11eb-2ee0-d9f9456ebbdb
# ╠═3e015876-39a3-11eb-21a1-d90460d77835
# ╠═2f4c2286-39ab-11eb-1c23-e998a265526a
# ╠═dcc75c1e-39ad-11eb-0b0b-13add6a0ec03
# ╠═6c5ea9ba-39a2-11eb-0fa2-71931a4498ee
# ╠═ae4cc5ea-39a1-11eb-2079-69a20a3b9a45
# ╠═442d4544-39a7-11eb-385c-035da6732479
# ╠═d7a666de-39a7-11eb-0537-ebd5786d023a
# ╠═6f2bde14-39ab-11eb-2197-d1af40ec3dd8
# ╠═68a7ee40-39ac-11eb-3b89-b9a85874f469
# ╠═266f8b74-39af-11eb-28fb-d7deda1a1a36
# ╠═f5c27d76-39af-11eb-3e30-d37e2db01542
# ╠═7cb9b56a-39b0-11eb-31a0-b98dcd35ad08
# ╠═86fef7d2-68a0-11eb-1974-4d2695221bdf
# ╠═9661e304-39b3-11eb-1373-c3959d5b37ff
# ╠═b41fdc84-39b3-11eb-0107-2fff3564a831
# ╠═e814282e-39b3-11eb-0698-5b00781a15c0
# ╠═b9893c4a-39b5-11eb-042a-99e76c8f37f9


text:
#===============================================================================

Using all of the digits 1 through 9 and concatenating them freely to form
decimal integers, different sets can be formed. Interestingly with the set
{2,5,47,89,631}, all of the elements belonging to it are prime.

How many distinct sets containing each of the digits one through nine exactly
once contain only prime elements?

--------------------------------------------------------------------------------

The algorithm proposed is based on the association of a number to its digits.

Let n=d1...dn be a number made of n digits. Then

    encode(n) = 0 if n contains 0 or if a digit appears multiple times.
              = (Σ 2^di)/2

The division by 2 results from the fact the 0 is not allowed, so Σ 2^di is
always odd. Thus the range of encode(n) is 1 to 511 (9 bits set to 1).

Each prime up to 987654321 are encoded as explained above, and each code is
associated with the number of primes it corresponds to, excepted for primes
encoded as 0, which are discarded.

A similar code is also associated to a set of primes, by computing it over all
the digits of all the primes belonging to the set.

With this, a prime p can be added to the set S if encode(p) and encode(S) have
no digits in common, i.e. if encode(p) & encode(S) is null. The code of the new
set is encode(p) | encode(S).

Of course, this applies to all primes whose code is the same as encode(p).

Then the solution is built by exploring recursively the tree of the
possibilities.

===============================================================================#

import Primes

function encode(n::Int)
    result::UInt16 = 0

    while (n != 0)
        n, digit = divrem(n, 10)

        bit = 1 << digit

        if (digit == 0 || result & bit != 0)
            return 0
        end

        result |= bit
    end

    result >> 1
end

primes = Primes.primes(1, 987654321);

function compute_codes(primes)
    result = zeros(Int64, 511)

    for p in primes
        code = encode(p)
        if code != 0
            result[code] += 1
        end
    end

    result;
end

const allcodes = compute_codes(primes)

const fullcode = encode(123456789)

function compute_nbr_sets(code_idx, code_of_set, nbr_of_prime_sets)

    if code_of_set == fullcode
        return nbr_of_prime_sets
    end

    if code_idx == 512
        return 0
    end

    result = compute_nbr_sets(code_idx+1, code_of_set, nbr_of_prime_sets);

    nbr_assoc_primes = allcodes[code_idx]

    if nbr_assoc_primes != 0 && code_of_set & code_idx == 0
        result += compute_nbr_sets(code_idx+1,
                                   code_of_set|code_idx,
                                   nbr_assoc_primes * nbr_of_prime_sets)
    end

    result
end

println("Euler 118: ", compute_nbr_sets(1, 0, 1));

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
function state1(m,c,k,f,dt,u0,v0)

	#--------------------------------------------------------------------------
	# Esta funcion calcula la respuesta de un sistema de 1 dgl aplicando el
	# metodo del espacio de los estados
	#
	# Datos de entrada
	# m          masa (kg)
	# k          rigidez (N/m)
	# c          amortiguamiento
	# f          vector de fuerzas (N)
	# dt         incremento de tiempo
	# u0         posicon inicial
	# v0         velocidad inicial
	#
	# Datos de salida
	# u          desplazamiento de un sistema de 1 gdl sometido a x
	# v          velocidad del sistema
	# a          aceleracion del sistema
	#
	# javier.cara@upm.es, january 2016
	#--------------------------------------------------------------------------

	# numero de instantes de tiempo
	N=length(f)

	# matrices Ac y Bc
	Ac = [0 1;-k/m -c/m]
	Bc = [0;1/m]

	# matrices A y B
	Ad = expm(Ac*dt)
	Bd = (Ad-eye(2))*inv(Ac)*Bc

	# calculo de los estados
	X = zeros(2,N)

	# posicion y velocidad iniciales:
	# en el instante cero tenemos las condiciones iniciales y el primer valor de f
	X[:,1] = [u0;v0]
	for r in 1:N-1
		 X[:,r+1] = Ad*X[:,r] + Bd*f[r]
	end

	u = X[1,:]
	v = X[2,:]
	a = - k/m*u - c/m*v + 1/m*f

	return u,v,a

end

# =============================================
function state1_test()

	m = 1.
	w=2*pi # frecuencia propia del sistema
	z=0.02 # tasa de amortiguamiento
	k = m*w^2
	c = 2*z*m*w
	dt = 0.01
	f=randn(1,100)
	u0 = 0.
	v0 = 0.

	u,v,a = state1(m,c,k,f,dt,u0,v0)

	return u,v,a

end


text:
function codegen(::Type{Val{:iterate}}, ::Type{SMMALA}, job::BasicMCJob)
  result::Expr
  update = []
  noupdate = []
  burninbody = []
  body = []

  vform = variate_form(job.pstate)
  if vform != Univariate && vform != Multivariate
    error("Only univariate or multivariate parameter states allowed in SMMALA code generation")
  end

  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    push!(body, :(_job.sstate.tune.proposed += 1))
  end

  push!(body, :(_job.sstate.μ = _job.pstate.value+0.5*_job.sstate.tune.step*_job.sstate.oldfirstterm))

  if vform == Univariate
    push!(body, :(_job.sstate.pstate.value = _job.sstate.μ+_job.sstate.sqrttunestep*_job.sstate.cholinvtensor*randn()))
  elseif vform == Multivariate
    push!(
      body,
      :(_job.sstate.pstate.value = _job.sstate.μ+_job.sstate.sqrttunestep*_job.sstate.cholinvtensor*randn(_job.pstate.size))
    )
  end

  push!(body, :(_job.parameter.uptotensorlogtarget!(_job.sstate.pstate)))

  if vform == Multivariate && job.sampler.transform != nothing
    push!(body, :(_job.sstate.pstate.tensorlogtarget = _job.sampler.transform(_job.sstate.pstate.tensorlogtarget)))
  end

  push!(body, :(_job.sstate.ratio = _job.sstate.pstate.logtarget-_job.pstate.logtarget))

  if vform == Univariate
    push!(
      body,
      :(
        _job.sstate.ratio += (
          0.5*(
            log(_job.sstate.tune.step*_job.sstate.oldinvtensor)
            +abs2(_job.sstate.pstate.value-_job.sstate.μ)*_job.pstate.tensorlogtarget/_job.sstate.tune.step
          )
        )
      )
    )
  elseif vform == Multivariate
    push!(
      body,
      :(
        _job.sstate.ratio += (
          0.5*(
            logdet(_job.sstate.tune.step*_job.sstate.oldinvtensor)
            +dot(
              _job.sstate.pstate.value-_job.sstate.μ,
              _job.pstate.tensorlogtarget*(_job.sstate.pstate.value-_job.sstate.μ)
            )/_job.sstate.tune.step
          )
        )
      )
    )
  end

  push!(body, :(_job.sstate.newinvtensor = inv(_job.sstate.pstate.tensorlogtarget)))

  push!(body, :(_job.sstate.newfirstterm = _job.sstate.newinvtensor*_job.sstate.pstate.gradlogtarget))

  push!(body, :(_job.sstate.μ = _job.sstate.pstate.value+0.5*_job.sstate.tune.step*_job.sstate.newfirstterm))

  if vform == Univariate
    push!(
      body,
      :(
        _job.sstate.ratio -= (
          0.5*(
            log(_job.sstate.tune.step*_job.sstate.newinvtensor)
            +abs2(_job.pstate.value-_job.sstate.μ)*_job.sstate.pstate.tensorlogtarget/_job.sstate.tune.step
          )
        )
      )
    )
  elseif vform == Multivariate
    push!(
      body,
      :(
        _job.sstate.ratio -= (
          0.5*(
            logdet(_job.sstate.tune.step*_job.sstate.newinvtensor)
            +dot(
              _job.pstate.value-_job.sstate.μ,
              _job.sstate.pstate.tensorlogtarget*(_job.pstate.value-_job.sstate.μ)
            )/_job.sstate.tune.step
          )
        )
      )
    )
  end

  if vform == Univariate
    push!(update, :(_job.pstate.value = _job.sstate.pstate.value))
    push!(update, :(_job.pstate.gradlogtarget = _job.sstate.pstate.gradlogtarget))
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      push!(update, :(_job.pstate.gradloglikelihood = _job.sstate.pstate.gradloglikelihood))
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      push!(update, :(_job.pstate.gradlogprior = _job.sstate.pstate.gradlogprior))
    end
    push!(update, :(_job.pstate.tensorlogtarget = _job.sstate.pstate.tensorlogtarget))
    if in(:tensorloglikelihood, job.outopts[:monitor]) && job.parameter.tensorloglikelihood! != nothing
      push!(update, :(_job.pstate.tensorloglikelihood = _job.sstate.pstate.tensorloglikelihood))
    end
    if in(:tensorlogprior, job.outopts[:monitor]) && job.parameter.tensorlogprior! != nothing
      push!(update, :(_job.pstate.tensorlogprior = _job.sstate.pstate.tensorlogprior))
    end
    push!(update, :(_job.sstate.oldinvtensor = _job.sstate.newinvtensor))
    push!(update, :(_job.sstate.cholinvtensor = chol(_job.sstate.newinvtensor, :L)))
    push!(update, :(_job.sstate.oldfirstterm = _job.sstate.newfirstterm))
  elseif vform == Multivariate
    push!(update, :(_job.pstate.value = copy(_job.sstate.pstate.value)))
    push!(update, :(_job.pstate.gradlogtarget = copy(_job.sstate.pstate.gradlogtarget)))
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      push!(update, :(_job.pstate.gradloglikelihood = copy(_job.sstate.pstate.gradloglikelihood)))
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      push!(update, :(_job.pstate.gradlogprior = copy(_job.sstate.pstate.gradlogprior)))
    end
    push!(update, :(_job.pstate.tensorlogtarget = copy(_job.sstate.pstate.tensorlogtarget)))
    if in(:tensorloglikelihood, job.outopts[:monitor]) && job.parameter.tensorloglikelihood! != nothing
      push!(update, :(_job.pstate.tensorloglikelihood = copy(_job.sstate.pstate.tensorloglikelihood)))
    end
    if in(:tensorlogprior, job.outopts[:monitor]) && job.parameter.tensorlogprior! != nothing
      push!(update, :(_job.pstate.tensorlogprior = copy(_job.sstate.pstate.tensorlogprior)))
    end
    push!(update, :(_job.sstate.oldinvtensor = copy(_job.sstate.newinvtensor)))
    push!(update, :(_job.sstate.cholinvtensor = chol(_job.sstate.newinvtensor, Val{:L})))
    push!(update, :(_job.sstate.oldfirstterm = copy(_job.sstate.newfirstterm)))
  end
  push!(update, :(_job.pstate.logtarget = _job.sstate.pstate.logtarget))
  if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
    push!(update, :(_job.pstate.loglikelihood = _job.sstate.pstate.loglikelihood))
  end
  if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
    push!(update, :(_job.pstate.logprior = _job.sstate.pstate.logprior))
  end
  if in(:accept, job.outopts[:diagnostics])
    push!(update, :(_job.pstate.diagnosticvalues[1] = true))
    push!(noupdate, :(_job.pstate.diagnosticvalues[1] = false))
  end
  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    push!(update, :(_job.sstate.tune.accepted += 1))
  end

  push!(body, Expr(:if, :(_job.sstate.ratio > 0 || (_job.sstate.ratio > log(rand()))), Expr(:block, update...), noupdate...))

  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    push!(burninbody, :(rate!(_job.sstate.tune)))

    if isa(job.tuner, AcceptanceRateMCTuner)
      push!(burninbody, :(tune!(_job.sstate.tune, _job.tuner)))
      push!(burninbody, :(_job.sstate.sqrttunestep = sqrt(_job.sstate.tune.step)))
    end

    if job.tuner.verbose
      fmt_iter = format_iteration(ndigits(job.range.burnin))
      fmt_perc = format_percentage()

      push!(burninbody, :(println(
        "Burnin iteration ",
        $(fmt_iter)(_job.sstate.tune.totproposed),
        " of ",
        _job.range.burnin,
        ": ",
        $(fmt_perc)(100*_job.sstate.tune.rate),
        " % acceptance rate"
      )))
    end

    push!(burninbody, :(reset_burnin!(_job.sstate.tune)))

    push!(body, Expr(
      :if,
      :(_job.sstate.tune.totproposed <= _job.range.burnin && mod(_job.sstate.tune.proposed, _job.tuner.period) == 0),
      Expr(:block, burninbody...)
    ))
  end

  if !job.plain
    push!(body, :(produce()))
  end

  @gensym _iterate

  result = quote
    function $_iterate(_job::BasicMCJob)
      $(body...)
    end
  end

  result
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule battery_cli_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("battery_cli")
JLLWrappers.@generate_main_file("battery_cli", UUID("dd074c01-bd9e-510d-b91d-361c690844c9"))
end  # module battery_cli_jll


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const GLMClassifier_PostEvaluationTransform = (;[
    Symbol("Logit") => Int32(0),
    Symbol("Probit") => Int32(1),
]...)

const GLMClassifier_ClassEncoding = (;[
    Symbol("ReferenceClass") => Int32(0),
    Symbol("OneVsRest") => Int32(1),
]...)

mutable struct GLMClassifier_DoubleArray <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GLMClassifier_DoubleArray(; kwargs...)
        obj = new(meta(GLMClassifier_DoubleArray), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct GLMClassifier_DoubleArray
const __meta_GLMClassifier_DoubleArray = Ref{ProtoMeta}()
function meta(::Type{GLMClassifier_DoubleArray})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GLMClassifier_DoubleArray)
            __meta_GLMClassifier_DoubleArray[] = target = ProtoMeta(GLMClassifier_DoubleArray)
            pack = Symbol[:value]
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{Float64}]
            meta(target, GLMClassifier_DoubleArray, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GLMClassifier_DoubleArray[]
    end
end
function Base.getproperty(obj::GLMClassifier_DoubleArray, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    else
        getfield(obj, name)
    end
end

mutable struct GLMClassifier <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GLMClassifier(; kwargs...)
        obj = new(meta(GLMClassifier), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct GLMClassifier
const __meta_GLMClassifier = Ref{ProtoMeta}()
function meta(::Type{GLMClassifier})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GLMClassifier)
            __meta_GLMClassifier[] = target = ProtoMeta(GLMClassifier)
            fnum = Int[1,2,3,4,100,101]
            pack = Symbol[:offset]
            allflds = Pair{Symbol,Union{Type,String}}[:weights => Base.Vector{GLMClassifier_DoubleArray}, :offset => Base.Vector{Float64}, :postEvaluationTransform => Int32, :classEncoding => Int32, :stringClassLabels => StringVector, :int64ClassLabels => Int64Vector]
            oneofs = Int[0,0,0,0,1,1]
            oneof_names = Symbol[Symbol("ClassLabels")]
            meta(target, GLMClassifier, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_GLMClassifier[]
    end
end
function Base.getproperty(obj::GLMClassifier, name::Symbol)
    if name === :weights
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{GLMClassifier_DoubleArray}
    elseif name === :offset
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :postEvaluationTransform
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :classEncoding
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :stringClassLabels
        return (obj.__protobuf_jl_internal_values[name])::StringVector
    elseif name === :int64ClassLabels
        return (obj.__protobuf_jl_internal_values[name])::Int64Vector
    else
        getfield(obj, name)
    end
end

export GLMClassifier_PostEvaluationTransform, GLMClassifier_ClassEncoding, GLMClassifier_DoubleArray, GLMClassifier


text:
abstract type AbstractTerm{𝒯} end
abstract type AbstractGravity{𝒯} <: AbstractTerm{𝒯} end

struct Gravity{𝒯} <: AbstractGravity{𝒯} end
struct Buoyancy{𝒯} <: AbstractGravity{𝒯} end

@inline calc_component!(source, ::Nothing, state, _...) = nothing

@inline function calc_component!(source, ::Gravity, state, aux, physics)
    ρ = state.ρ
    ∇Φ = aux.∇Φ

    source.ρu -= ρ * ∇Φ

    nothing
end

# really ρe should be ρeᵢₙₜ
@inline function calc_component!(source, ::Gravity{(:ρ, :ρu, :ρe)}, state, aux, physics)
    ρ = state.ρ
    ρu = state.ρu
    ∇Φ = aux.∇Φ

    source.ρu -= ρ * ∇Φ
    source.ρe -= ρu' * ∇Φ

    nothing
end

@inline function calc_component!(source, ::Buoyancy, state, aux, physics)
    ρθ = state.ρθ
    k = vertical_unit_vector(aux.orientation, aux)
    α = physics.params.α
    g = physics.params.g

    source.ρu -= -α * g * k * ρθ

    nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
@test first(evaluate5(A♣, K♣, Q♣, J♣, T♣)) == 1  #  Royal Flush
@test first(evaluate5(K♣, Q♣, J♣, T♣, 9♣)) == 2  #  King-High Straight Flush
@test first(evaluate5(Q♣, J♣, T♣, 9♣, 8♣)) == 3  #  Queen-High Straight Flush
@test first(evaluate5(J♣, T♣, 9♣, 8♣, 7♣)) == 4  #  Jack-High Straight Flush
@test first(evaluate5(T♣, 9♣, 8♣, 7♣, 6♣)) == 5  #  Ten-High Straight Flush
@test first(evaluate5(9♣, 8♣, 7♣, 6♣, 5♣)) == 6  #  Nine-High Straight Flush
@test first(evaluate5(8♣, 7♣, 6♣, 5♣, 4♣)) == 7  #  Eight-High Straight Flush
@test first(evaluate5(7♣, 6♣, 5♣, 4♣, 3♣)) == 8  #  Seven-High Straight Flush
@test first(evaluate5(6♣, 5♣, 4♣, 3♣, 2♣)) == 9  #  Six-High Straight Flush
@test first(evaluate5(5♣, 4♣, 3♣, 2♣, A♣)) == 10 #  Five-High Straight Flush


text:
using Test
using CollectiveSpins
import LinearAlgebra

#all functions should have immediate rejection if k exceeds corner of BZ

@testset "collective_modes" begin

cs = CollectiveSpins
la = LinearAlgebra
k0 = 2pi
tol = 1e-05
tol_1D = 1e-02
tol_2D = 1e-01
N = 600

a_list = collect(1:5)*0.1 .- 0.001
for a in a_list
    #assert that lighcone exists a < lambda0/2
    @test CollectiveSpins.collective_modes.Gamma_k_chain(pi/a, a, [1,0,0]) == 0 #arbitrary, unnormalized initializations throughout
    @test CollectiveSpins.collective_modes.Gamma_k_chain(pi/a, a, [0,1,0.5]) == 0
end

a_list = collect(1:5)*1e+08
for a in a_list
    #assert lim a--> Inf, Omega_k --> 0
    @test abs(cs.collective_modes.Omega_k_chain(0.2*a, a, [1,0,0])) < tol
    @test abs(cs.collective_modes.Omega_k_chain(0.2*a, a, [0,1,0.5])) < tol
    #assert lim a--> Inf, Delta_k --> gamma0 (1)
    @test abs(cs.collective_modes.Gamma_k_chain(0, a, [1,0,0]) - 1) < tol
    @test abs(cs.collective_modes.Gamma_k_chain(0, a, [0,1,0.5]) - 1) < tol
end

a_list = [0.07, 0.11, 0.2, 0.51, 1.7]
k_list = [0.4pi, 0.1pi, 0pi, 0.9pi, 0.3pi]
k_list = k_list ./ a_list
Omega_x_calculated = [-6.24945, -12.2009, -2.18103, -0.11888, 0.0165684]
Omega_yz_calculated = [3.42368, 4.99508, 1.28359, -0.494113, 0.0327376]
Gamma_x_calculated = [0, 5.40947, 3.75, 0.325667, 1.00791]
Gamma_yz_calculated = [0, 4.11345, 1.875, 1.30775, 0.819573]
for index = 1:length(a_list)
    #assert that selection of arbitrary a and k values match simplified, semi-manual calculation
    @test abs(cs.collective_modes.Omega_k_chain(k_list[index], a_list[index], [1,0,0]) - Omega_x_calculated[index]) < tol
    @test abs(cs.collective_modes.Omega_k_chain(k_list[index], a_list[index], [0,5,22]) - Omega_yz_calculated[index]) < tol
    @test abs(cs.collective_modes.Gamma_k_chain(k_list[index], a_list[index], [55,0,0]) - Gamma_x_calculated[index]) < tol
    @test abs(cs.collective_modes.Gamma_k_chain(k_list[index], a_list[index], [0,87,2.3+2im]) - Gamma_yz_calculated[index]) < tol
end

#real-space Dyadic Green's function
function Gr(r_vec::Array{<:Number, 1})
    r = la.norm(r_vec)
    return exp(im*k0*r)/(4pi*k0^2*r^3) * ( la.UniformScaling(k0^2*r^2 + im*k0*r - 1) + (-k0^2*r^2 - 3*im*k0*r + 3)*r_vec*la.transpose(r_vec)/r^2 )
end

#1D k-space Green's function by summation of real-space
function Gk_sumR_1D(k, a, polarization::Array{<:Number, 1})
    polarization = la.normalize(polarization)
    Omegak = 0
    Gammak = 0
    for n = -N:N
        if !(n == 0)
            r_vec = zeros(3)
            r_vec[1] = n*a
            Omegak += la.conj(la.transpose(polarization))*real(Gr(r_vec) * exp(Complex(-im*k*r_vec[1])))*polarization
            Gammak += la.conj(la.transpose(polarization))*imag(Gr(r_vec) * exp(Complex(-im*k*r_vec[1])))*polarization
        end
    end
    if  abs(k) > k0
        return -3pi*real(Omegak)/k0, 0
    end
    return -3pi*real(Omegak)/k0, 6pi*real(Gammak)/k0 + 1
end

a_list = [0.8, 0.25, 1.2, 0.3, 0.5, 1.7, 0.6]
k_list = [0.1pi/0.8, -0.1pi/0.25, 0.3pi/1.2, -0.1pi/0.3, -0.2pi/0.5, 0.3pi/1.7, 0.2pi/0.6]
polarization_list = [[1, 1, -1], [-1, -im, -1], [1, 0, im], [0, -im, 5], [1, 0, 0], [-1, im, 1], [0, 1, -2], [-4, -1, 2]]
for index = 1:length(a_list)
    #assert that 1D k-space calculations with G(k) match G_k(r) definition. Convergence of real-space sum version is unstable by nature
    #thus only moderate a values are taken and tolerance is quite large (1e-02). This is especially true around the edge of the lightcone.
    Omegak_real, Gammak_real = Gk_sumR_1D(k_list[index], a_list[index], polarization_list[index])
    @test abs(cs.collective_modes.Omega_k_chain(k_list[index], a_list[index], polarization_list[index]) -
    Omegak_real) < tol_1D
    @test abs(cs.collective_modes.Gamma_k_chain(k_list[index], a_list[index], polarization_list[index]) -
    Gammak_real) < tol_1D
end

#2D k-space Green's function by summation of real-space
function Gk_sumR(k_vec::Array{<:Number, 1}, a_vec1::Array{<:Number, 1}, a_vec2::Array{<:Number, 1}, polarization::Array{<:Number, 1})
    polarization = la.normalize(polarization)
    Omegak = 0
    Gammak = 0
    for n1 = -N:N
        for n2 = -N:N
            if !(n1 == 0) || !(n2 == 0)
                r_vec = zeros(3)
                r_vec[2:3] = n1*a_vec1 + n2*a_vec2
                Omegak += la.conj(la.transpose(polarization))*real(Gr(r_vec) * exp(Complex(-im*la.dot(k_vec, r_vec[2:3]))))*polarization
                Gammak += la.conj(la.transpose(polarization))*imag(Gr(r_vec) * exp(Complex(-im*la.dot(k_vec, r_vec[2:3]))))*polarization
            end
        end
    end
    if la.norm(k_vec) > k0
        return -3pi*real(Omegak)/k0, 0
    end
    return -3pi*real(Omegak)/k0, 6pi*real(Gammak)/k0 + 1
end

a_vec1_list = [[0.8, 0], [0.25, 0], [1.2, 0], [0.3, 0], [0.25/2, -0.25*sqrt(3)/2], [0.5/2, -0.5/2*sqrt(3)/2], [1.7, 0.25], [0.2, 0.6]]
a_vec2_list = [[0, 0.8], [0, 0.25], [0, 0.4], [0, 0.7], [0.25/2, 0.25*sqrt(3)/2], [0.5/2, 0.5*sqrt(3)/2], [0.25, 0.9], [0.7, 0.3]]
k_vec_list = [[-0.1pi/0.8, 0.3pi/0.8], [-0.1pi/0.25, -0.1pi/0.25], [0.6pi/1.2, 0.99pi/0.4], [0.1pi/0.3, -0.2pi/0.7], [0.2pi/0.25, 0.5pi/0.25], [1pi/0.5, 1pi/0.5], [0, 0.2pi/0.8], [-0.1pi/0.5, -0.2pi/0.4]]
polarization_list = [[1, 1, 1], [-1, im, 1], [-1, 0, -im], [0, im, 5], [1, 0, 0], [-1, im, 1], [0, 1, 2], [4, -1, -2]]
for index = 1:length(a_vec1_list)
    #assert that 2D k-space calculations with G(k) match G_k(r) definition. Convergence of real-space sum version is unstable by nature
    #thus only moderate a values are taken and tolerance is quite large (1e-01). This is especially true around the edge of the lightcone.
    Omegak_real, Gammak_real = Gk_sumR(k_vec_list[index], a_vec1_list[index], a_vec2_list[index], polarization_list[index])
    @test abs(cs.collective_modes.Omega_k_2D(k_vec_list[index], a_vec1_list[index], a_vec2_list[index], polarization_list[index]) -
    Omegak_real) < tol_2D
    @test abs(cs.collective_modes.Gamma_k_2D(k_vec_list[index], a_vec1_list[index], a_vec2_list[index], polarization_list[index]) -
    Gammak_real) < tol_2D
end

end #module


text:
"""
[Napoleon' theorem](https://en.wikipedia.org/wiki/Napoleon%27s_theorem)
states that if equilateral triangles are constructed on the sides
of any triangle, either all outward or all inward, the lines connecting the centres of those
equilateral triangles themselves form an equilateral triangle.
"""
module Napoleon

using Plots
using SymPy
using PlaneGeometry
using PlaneGeometry.GeoPlots

export
    napoleon_proof, napoleon_draw, napoleon_rand, napoleon_tri

"""
    napoleon_draw(📐::Triangle, labels='A':'C')

Verify Napoleon's theorem for the triangle `📐`.  The function works the same ways as
`napoleon_tri` but return `(🖼, hold)` where `🖼` is a plot and `hold` indicates where the theorem
holds.
"""
function napoleon_draw(tri::Triangle, labels='A':'C')
    trishape=shape(tri)

    plot(trishape, leg=false, fill=(0, :green), aspect_ratio=:equal, fillalpha= 0.2)

    scatter!(trishape.x, trishape.y, color=:red, series_annotations = text.(labels, :bottom))

    outer_tri = outer_equitriangles(tri)
    for t in outer_tri
        plot!(shape(t), leg=false, fill=(0, :pink), aspect_ratio=:equal, fillalpha=0.7)
    end

    outer_circles = circumcircle.(vertices.(outer_tri))

    outer_centers = PlaneGeometry.center.(outer_circles);

    nap_tri = Triangle(outer_centers...);

    plt = plot!(shape(nap_tri), leg=false, fill=(0, :blue), fillalpha=0.7, aspect_ratio=:equal)

    hold = isequilateral(nap_tri)

    plt, hold
end

"""
    napoleon_draw(A, B, C)

Verify Napoleon's theorem.
"""
function napoleon_draw(A, B, C)
    tri = Triangle(A, B, C)

    napoleon_draw(tri)
end

"""
    napoleon_draw(xA, yA, xB, yB, xC, yC)

Verify Napoleon's theorem.
"""
function napoleon_draw(xA, yA, xB, yB, xC, yC)
    A = Point(xA, yA);
    B = Point(xB, yB);
    C = Point(xC, yC);

    napoleon_draw(A, B, C)
end

"""
    napoleon_rand()

Verify Napoleon's theorem for a random triangle.
"""
function napoleon_rand()
    pts = rand(0:1//10:1,6);
    # Todo: sometimes this fails. Fix it.
    plt, hold = napoleon_draw(pts...)
end

"""
    napoleon_proof()

Prove Napoleon's theorem.
"""
function napoleon_proof()
    @vars by cx positive=true;
    @vars cy;

    A = Point(0, 0); B = Point(0, by); C = Point(cx, cy);
    tri = Triangle(A, B, C)

    npt = napoleon_tri(tri)
    isequilateral(npt)
end

"""
    napoleon_tri(📐)

Find the triangle that connects the center of three outer equilateral triangles of `📐`.
"""
function napoleon_tri(tri)
    outer_tri = outer_equitriangles(tri);
    outer_circles = circumcircle.(vertices.(outer_tri));
    outer_centers = PlaneGeometry.center.(outer_circles);
    Triangle(outer_centers...)
end

end