text:
@testset "Power Models SOCWR tests" begin

@testset "case 5" begin
    println("==================================")
    println("SOCWRPowerModel case5.m")
    println("==================================")

    pm = build_generic_model("data/pglib_opf_case5_pjm.m", SOCWRPowerModel, PowerModels.post_ots)
    m = pm.model
    @variable(m, 0 <= aeiou <= 1)
    @NLconstraint(m, aeiou^2== 1)

    setsolver(m, juniper_strong_no_restart)
    status = solve(m)

    @test status == :Optimal || status == :LocalOptimal

    juniper_val = getobjectivevalue(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 14999.7, atol=1e0)
end

@testset "case 5 w inc constr" begin
    println("==============================================")
    println("SOCWRPowerModel case5.m no incumbent constr")
    println("==============================================")

    pm = build_generic_model("data/pglib_opf_case5_pjm.m", SOCWRPowerModel, PowerModels.post_ots)
    m = pm.model
    @variable(m, 0 <= aeiou <= 1)
    @NLconstraint(m, aeiou^2== 1)

    solver = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        strong_restart = false,
        incumbent_constr = true,
        traverse_strategy = :DFS
    )
    setsolver(m, solver)
    status = solve(m)

    @test status == :Optimal || status == :LocalOptimal

    juniper_val = getobjectivevalue(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test m.internalModel.options.strong_restart == false
    @test m.internalModel.options.incumbent_constr == true
    @test isapprox(juniper_val, 14999.7, atol=1e0)
end

@testset "case 5 epsilon obj constr" begin
    println("==============================================")
    println("SOCWRPowerModel case5.m epsilon constr")
    println("==============================================")

    pm = build_generic_model("data/pglib_opf_case5_pjm.m", SOCWRPowerModel, PowerModels.post_ots)
    m = pm.model
    @variable(m, 0 <= aeiou <= 1)
    @NLconstraint(m, aeiou^2== 1)

    solver = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        strong_branching_perc = 0, # this defaults to 2 variables
        strong_restart = false,
        obj_epsilon = 0.5,
        traverse_strategy = :DFS
    )
    setsolver(m, solver)
    status = solve(m)

    @test status == :Optimal || status == :LocalOptimal

    juniper_val = getobjectivevalue(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 14999.7, atol=1e0)
end

end


text:
"""
Implementation of the 1-Rule-algorithm, which finds classification rules from a set of instances.
The result is a one-level deicsion tree (here of type `OneTree`).
For more information about the 1-Rule-algorithm see e.g.: https://datacadamia.com/data_mining/one_rule

The implementation of this algorithm is intended for teaching purposes; it's not for production use.
"""

module OneRule

using DataFrames

export  OneTree,
        OneNode,
        get_best_tree,
        get_nodes

include("nodes.jl")
include("trees.jl")

end # module

ext: jl

max_stars_repo_name: roland-KA/OneRule

max_issues_repo_name: roland-KA/OneRule

max_forks_repo_name: roland-KA/OneRule


text:
"""
```
hessian!(m::AbstractDSGEModel, x::Vector{T}, data::AbstractArray;
         verbose::Symbol = :none) where {T<:AbstractFloat}
```

Compute Hessian of DSGE posterior function evaluated at x.
"""
function hessian!(m::AbstractDSGEModel,
                  x::Vector{T},
                  data::AbstractArray;
                  verbose::Symbol = :none) where T<:AbstractFloat
    DSGE.update!(m, x)

    # Index of free parameters
    para_free      = [!θ.fixed for θ in m.parameters]
    para_free_inds = findall(para_free)

    # Compute hessian only for freem parameters with indices less than max. Useful for
    # testing purposes.
    max_free_ind = n_hessian_test_params(m)
    if max_free_ind < maximum(para_free_inds)
        para_free_inds = para_free_inds[1:max_free_ind]
    end

    n_para = length(x)
    hessian  = zeros(n_para, n_para)

    # x_hessian is the vector of free params
    # x_model is the vector of all params
    x_model = copy(x)
    x_hessian = x_model[para_free_inds]

    function f_hessian(x_hessian)
        x_model[para_free_inds] = x_hessian
        return -posterior!(m, x_model, data)
    end

    distr=use_parallel_workers(m)
    hessian_free, has_errors = hessizero(f_hessian, x_hessian;
        check_neg_diag=true, verbose=verbose, distr=distr)

    # Fill in rows/cols of zeros corresponding to location of fixed parameters
    # For each row corresponding to a free parameter, fill in columns corresponding to free
    # parameters. Everything else is 0.
    for (row_free, row_full) in enumerate(para_free_inds)
        hessian[row_full,para_free_inds] = hessian_free[row_free,:]
    end

    return hessian, has_errors
end


text:
# Arnoldi/Lanczos iteration algorithms (with custom IOP)

#######################################
# Output type/cache
"""
    KrylovSubspace{T}(n,[maxiter=30]) -> Ks

Constructs an uninitialized Krylov subspace, which can be filled by `arnoldi!`.

The dimension of the subspace, `Ks.m`, can be dynamically altered but should
be smaller than `maxiter`, the maximum allowed arnoldi iterations.

    getV(Ks) -> V
    getH(Ks) -> H

Access methods for the (extended) orthonormal basis `V` and the (extended)
Gram-Schmidt coefficients `H`. Both methods return a view into the storage
arrays and has the correct dimensions as indicated by `Ks.m`.

    resize!(Ks, maxiter) -> Ks

Resize `Ks` to a different `maxiter`, destroying its contents.

This is an expensive operation and should be used scarcely.
"""
mutable struct KrylovSubspace{B, T, U}
    m::Int        # subspace dimension
    maxiter::Int  # maximum allowed subspace size
    beta::B       # norm(b,2)
    V::Matrix{T}  # orthonormal bases
    H::Matrix{U}  # Gram-Schmidt coefficients (real for Hermitian matrices)
    KrylovSubspace{T,U}(n::Integer, maxiter::Integer=30) where {T,U} = new{real(T), T, U}(
        maxiter, maxiter, zero(real(T)), Matrix{T}(undef, n, maxiter + 1),
        fill(zero(U), maxiter + 1, maxiter))
    KrylovSubspace{T}(n::Integer, maxiter::Integer=30) where {T} = KrylovSubspace{T,T}(n, maxiter)
end
getH(Ks::KrylovSubspace) = @view(Ks.H[1:Ks.m + 1, 1:Ks.m])
getV(Ks::KrylovSubspace) = @view(Ks.V[:, 1:Ks.m + 1])
function Base.resize!(Ks::KrylovSubspace{B,T,U}, maxiter::Integer) where {B,T,U}
    V = Matrix{T}(undef, size(Ks.V, 1), maxiter + 1)
    H = fill(zero(U), maxiter + 1, maxiter)
    Ks.V = V; Ks.H = H
    Ks.m = Ks.maxiter = maxiter
    return Ks
end
function Base.show(io::IO, Ks::KrylovSubspace)
    println(io, "$(Ks.m)-dimensional Krylov subspace with fields")
    println(io, "beta: $(Ks.beta)")
    print(io, "V: ")
    println(IOContext(io, :limit => true), getV(Ks))
    print(io, "H: ")
    println(IOContext(io, :limit => true), getH(Ks))
end

# Helper functions that returns the real part if that is what is
# required (for Hermitian matrices), otherwise returns the value
# as-is.
coeff(::Type{T},α::T) where {T} = α
coeff(::Type{U},α::T) where {U<:Real,T<:Complex} = real(α)

#######################################
# Arnoldi/Lanczos with custom IOP
"""
    arnoldi(A,b[;m,tol,opnorm,iop]) -> Ks

Performs `m` anoldi iterations to obtain the Krylov subspace K_m(A,b).

The n x (m + 1) basis vectors `getV(Ks)` and the (m + 1) x m upper Hessenberg
matrix `getH(Ks)` are related by the recurrence formula

```
v_1=b,\\quad Av_j = \\sum_{i=1}^{j+1}h_{ij}v_i\\quad(j = 1,2,\\ldots,m)
```

`iop` determines the length of the incomplete orthogonalization procedure [^1].
The default value of 0 indicates full Arnoldi. For symmetric/Hermitian `A`,
`iop` will be ignored and the Lanczos algorithm will be used instead.

Refer to `KrylovSubspace` for more information regarding the output.

Happy-breakdown occurs whenver `norm(v_j) < tol * opnorm(A, Inf)`, in this case
the dimension of `Ks` is smaller than `m`.

[^1]: Koskela, A. (2015). Approximating the matrix exponential of an
advection-diffusion operator using the incomplete orthogonalization method. In
Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353).
Springer, Cham.
"""
function arnoldi(A, b; m=min(30, size(A, 1)), tol=1e-7, opnorm=LinearAlgebra.opnorm, iop=0, cache=nothing )
    Ks = KrylovSubspace{eltype(b)}(length(b), m)
    arnoldi!(Ks, A, b; m=m, tol=tol, opnorm=opnorm, iop=iop)
end

"""
    arnoldi_step!(j, iop, n, A, V, H)

Take the `j`:th step of the Lanczos iteration.
"""
function arnoldi_step!(j::Integer, iop::Integer, A,
                       V::AbstractMatrix{T}, H::AbstractMatrix{U}) where {T,U}
    x,y = @view(V[:, j]),@view(V[:, j+1])
    mul!(y, A, x)
    @inbounds for i = max(1, j - iop + 1):j
        alpha = coeff(U, dot(@view(V[:, i]), y))
        H[i, j] = alpha
        axpy!(-alpha, @view(V[:, i]), y)
    end
    beta = norm(y)
    H[j+1, j] = beta
    @. y /= beta
    beta
end
"""
    arnoldi!(Ks,A,b[;tol,m,opnorm,iop]) -> Ks

Non-allocating version of `arnoldi`.
"""
function arnoldi!(Ks::KrylovSubspace{B, T, U}, A, b::AbstractVector{T};
                  tol::Real=1e-7, m::Int=min(Ks.maxiter, size(A, 1)),
                  opnorm=LinearAlgebra.opnorm, iop::Int=0, cache=nothing) where {B, T <: Number, U <: Number}
    if ishermitian(A)
        return lanczos!(Ks, A, b; tol=tol, m=m, opnorm=opnorm)
    end
    if m > Ks.maxiter
        resize!(Ks, m)
    else
        Ks.m = m # might change if happy-breakdown occurs
    end
    V, H = getV(Ks), getH(Ks)
    vtol = tol * opnorm(A, Inf)
    if iop == 0
        iop = m
    end
    # Safe checks
    n = size(V, 1)
    @assert length(b) == size(A,1) == size(A,2) == n "Dimension mismatch"
    # Arnoldi iterations (with IOP)
    fill!(H, zero(U))
    Ks.beta = norm(b)
    @. V[:, 1] = b / Ks.beta
    @inbounds for j = 1:m
        beta = arnoldi_step!(j, iop, A, V, H)
        if beta < vtol # happy-breakdown
            Ks.m = j
            break
        end
    end
    return Ks
end

"""
    lanczos_step!(j, m, n, A, V, H)

Take the `j`:th step of the Lanczos iteration.
"""
function lanczos_step!(j::Integer, A,
                       V::AbstractMatrix{T},
                       α::AbstractVector{U},
                       β::AbstractVector{B}) where {B,T,U}
    x,y = @view(V[:, j]),@view(V[:, j+1])
    mul!(y, A, x)
    α[j] = coeff(U, dot(x, y))
    axpy!(-α[j], x, y)
    j > 1 && axpy!(-β[j-1], @view(V[:, j-1]), y)
    β[j] = norm(y)
    @. y /= β[j]
    β[j]
end

macro diagview(A,d::Integer=0)
    s = d<=0 ? 1+abs(d) : :(m+$d)
    quote
        m = size($(esc(A)),1)
        @view($(esc(A))[($s):m+1:end])
    end
end

"""
    realview(R, V)

Returns a view of the real components of the complex vector `V`.
"""
realview(::Type{R}, V::AbstractVector{C}) where {R,C<:Complex} =
    @view(reinterpret(R, V)[1:2:end])
"""
    realview(R, V)

Returns a view of the real components of the real vector `V`.
"""
realview(::Type{R}, V::AbstractVector{R}) where {R} = V

"""
    lanczos!(Ks,A,b[;tol,m,opnorm]) -> Ks

A variation of `arnoldi!` that uses the Lanczos algorithm for
Hermitian matrices.
"""
function lanczos!(Ks::KrylovSubspace{B, T, U}, A, b::AbstractVector{T};
                  tol=1e-7, m=min(Ks.maxiter, size(A, 1)),
                  opnorm=LinearAlgebra.opnorm, cache=nothing) where {B, T <: Number, U <: Number}
    if m > Ks.maxiter
        resize!(Ks, m)
    else
        Ks.m = m # might change if happy-breakdown occurs
    end
    V, H = getV(Ks), getH(Ks)
    vtol = tol * opnorm(A, Inf)
    # Safe checks
    n = size(V, 1)
    @assert length(b) == size(A,1) == size(A,2) == n "Dimension mismatch"
    # Lanczos iterations
    fill!(H, zero(T))
    Ks.beta = norm(b)
    @. V[:, 1] = b / Ks.beta
    α = @diagview(H)
    # β is always real, even though α may (in general) be complex.
    β = realview(B, @diagview(H,-1))
    @inbounds for j = 1:m
        if vtol > lanczos_step!(j, A, V, α, β)
            # happy-breakdown
            Ks.m = j
            break
        end
    end
    copyto!(@diagview(H,1), β[1:end-1])
    return Ks
end


text:
@testset "caches" begin

function f!(fvec, x)
    fvec[1] = (x[1]+3)*(x[2]^3-7)+18
    fvec[2] = sin(x[2]*exp(x[1])-1)
end

df = OnceDifferentiable(f!, rand(2), rand(2))

# Test that the cache is actually passed all the way by verifying that it's
# modified during the solve calls.
nc = NLsolve.NewtonCache(df)
nc.p .= 22.0
ncp = copy(nc.p)
r = NLsolve.newton(df, [ 1.; 1.], 0.1, 0.1, 100, false, false, false, LineSearches.Static(), nc)
@test !(ncp == nc.p)

ac = NLsolve.AndersonCache(df, NLsolve.Anderson{10}())
ac.x .= 22.0
acx = copy(ac.x)
r = NLsolve.anderson(df, [ 1.; 1.], 0.1, 0.1, 100, false, false, false, 0.9, 1, 0, ac)
@test !(acx == ac.x)

ntc = NLsolve.NewtonTrustRegionCache(df)
ntc.r_predict .= 22.0
ntcr_predict = copy(ntc.r_predict)
r = NLsolve.trust_region(df, [ 1.; 1.], 0.1, 0.1, 100, false, false, false, 1.0, true, ntc)
@test !(ntcr_predict == ntc.r_predict)

end


text:
using MCPhylo
using Test

@testset "insert_node!" begin
    tree = MCPhylo.parsing_newick_string("(A,B,(C,D,E)F)G;")
    insert_tree = MCPhylo.parsing_newick_string("(A,B,((C,D,E)no_name)F)G;")
    insert_tree_newick = newick(insert_tree)
    insert_tree2 = MCPhylo.parsing_newick_string("(A,B,(((C,D)no_name,E)no_name)F)G;")
    insert_tree_newick2 = newick(insert_tree2)
    MCPhylo.number_nodes!(tree)

    children = [find_by_name(tree, "C"), find_by_name(tree, "D"),
                find_by_name(tree, "E")]
    mother = find_by_name(tree, "F")
    inserted_node = insert_node!(mother, children)
    @test newick(tree) == insert_tree_newick

    pop!(children)
    insert_node!(inserted_node, children)
    @test newick(tree) == insert_tree_newick2

    new_node = Node("G")
    append!(children, new_node)
    @test_throws AssertionError insert_node!(mother, children)
 end

@testset "delete_node!" begin
    tree = MCPhylo.parsing_newick_string("(A,B,(C,D)E)F;")
    MCPhylo.number_nodes!(tree)
    remove_tree = MCPhylo.parsing_newick_string("(A,B,C,D)F;")
    remove_tree_newick = newick(remove_tree)
    remove_tree2 = MCPhylo.parsing_newick_string("(B,C,D)F;")
    remove_tree2_newick = newick(remove_tree2)

    delete_node!(find_by_name(tree, "E"))
    @test newick(tree) == remove_tree_newick
    delete_node!(find_by_name(tree, "A"))
    @test newick(tree) == remove_tree2_newick
    @test_throws ArgumentError delete_node!(find_by_name(tree, "F"))
end

@testset "find_lca" begin
    tree = MCPhylo.parsing_newick_string("(A,B,(C,(D,E)F)G)H;")
    MCPhylo.set_binary!(tree)

    @testset "find_lca by name" begin
        @test find_lca(tree, ["A","C"]) == tree
        @test find_lca(tree, ["D", "E"]) == find_by_name(tree, "F")
        @test find_lca(tree, ["D", "F"]) == find_by_name(tree, "F")
        @test find_lca(tree, ["C", "D"]) == find_by_name(tree, "G")
        @test find_lca(tree, ["A","B","C"]) == tree
        @test find_lca(tree, ["D", "E", "F", "G"]) == find_by_name(tree, "G")
        @test find_lca(tree, ["D", "E", "C"]) == find_by_name(tree, "G")
    end # find_lca by name

    @testset "find_lca by nodes" begin
        nodes = [find_by_name(tree, "A"), find_by_name(tree, "B"),
                 find_by_name(tree, "C")]
        @test find_lca(tree, nodes) == tree

        nodes = [find_by_name(tree, "D"), find_by_name(tree, "E"),
                 find_by_name(tree, "F"), find_by_name(tree, "G")]
        @test find_lca(tree, nodes) == find_by_name(tree, "G")
    end # find_lca by nodes
end

@testset "check_leafsets" begin
    tree1 = MCPhylo.parsing_newick_string("(((A,B),C),(D,E))")
    tree2 = MCPhylo.parsing_newick_string("((A,C),(B,D,E))")
    tree3 = MCPhylo.parsing_newick_string("(((G,C),A),D,F)")
    tree4 = MCPhylo.parsing_newick_string("(((B,C),A),D,E)")
    tree5 = MCPhylo.parsing_newick_string("(((B,C),A),D,F)")

    trees = [tree1, tree2, tree3, tree4, tree5]
    MCPhylo.number_nodes!.(trees)
    MCPhylo.set_binary!.(trees)

    @test_throws ArgumentError MCPhylo.check_leafsets(trees)
end


text:
module KmerAnalysis

export
    CountMode,
    Canonical,
    CANONICAL,
    NonCanonical,
    NONCANONICAL,

    # MerCount
    MerCount,
    mer,
    freq,
    merge,
    merge_into!,
    collapse_into_counts,
    collapse_into_counts!,
    collect_mers,
    collect_mers!,

    # Counters
    serial_mem,
    count!,
    dist_mem,

    # K-mer frequency spectra
    KmerFrequencySpectra,
    spectra,
    spectra!,
    find_next_peak,
    find_spectra_peaks,

    # Indexed k-mer counts
    IndexedCounts,
    add_count!,
    project_count

using BioSequences, ReadDatastores, Distributed

abstract type CountMode end

struct Canonical    <: CountMode end
struct NonCanonical <: CountMode end

@inline (::Canonical)(x) = canonical(x)
@inline (::NonCanonical)(x) = fwmer(x)

const CANONICAL = Canonical()
const NONCANONICAL = NonCanonical()

abstract type AbstractKmerCounter{M<:AbstractMer} <: Function end

include("MerCount.jl")
include("counters/serial_mem.jl")
include("counters/dist_mem.jl")
include("IndexedCounts.jl")
include("KmerFrequencySpectra.jl")
include("kci.jl")

end # module


text:
function smooth_pi0{T<:FloatingPoint}(pValues::Vector{T}, lambda = 0.05:0.05:0.95)
    validPvalues(pValues)
    #validPvalues(lambda) ## TODO check bounds
    n = length(pValues)
    ## CHCK: Only for smoothing (due to cubic spline)?
    if length(lambda) < 4
        throw(ArgumentError())
    end
    if !issorted(lambda)
        throw(DomainError())
    end
    pi0 = Float64[mean(pValues .>= l) / (1-l) for l in lambda]
    ##CoordInterpGrid(lambda, pi0, BCnil, InterpCubic) ## broken
end


text:
type Foo
    x::Int
end

function plot_err(x::Foo)
    check_defined_err(:PyPlot)
    d = 1:x.x
    plot(d, d.^2)
end

function plot_eval(x::Foo)
    if !isdefined(:PyPlot)
        eval(Expr(:using, :PyPlot))
    end
    d = 1:x.x
    PyPlot.plot(d, d.^2)
end

a = Foo(10)
try
    plot_err(a)
catch err
    println(err)
end

using PyPlot

plot_eval(a)

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
# Methods for absolute and relative abundances

abundancetable(df::DataFrame) = ComMatrix(
    convert(Matrix{Float64}, df[2:end]),
        string.(df[1]),
        string.(names(df[2:end]))
    )

abundancetable(table::AbstractArray{T,2},
    site = ["sample_$x" for x in axes(table, 2)],
    species = ["feature_$x" for x in axes(table, 1)]
    ) where T<:Real = ComMatrix(Float64.(table), species, site)

"""
Filter an abundance table to the top `n` features accross all samples

This function also adds a row for "other", which sums the abundances of the
remaining features.
"""
function filterabund(abun::AbstractComMatrix, n::Int=minimum(10, nfeatures(abun)))
    # TODO: add prevalence filter

    totals = featuretotals(abun)

    srt = sortperm(totals, rev=true)

    newabun = getfeature(abun, srt[1:n])

    remainder = [sum(occurrences(abun)[srt[n+1:end], i]) for i in 1:size(abun, 2)]'
    newabun = vcat(newabun, remainder)
    newrows = cat(featurenames(abun)[srt[1:n]], ["other"], dims=1)

    return abundancetable(newabun, samplenames(abun), newrows)
end

filterabund(df::DataFrame, n::Int=10) = filterabund(abundancetable(df), n)

function rownormalize!(abt::AbstractComMatrix)
    for i in 1:size(abt, 1)
        rowmax = maximum(occurrences(abt)[i,:])
        for j in 1:size(abt, 2)
            occurrences(abt)[i,j] /= rowmax
        end
    end
end

function rownormalize(abt::AbstractComMatrix)
    abt = deepcopy(abt)
    rownormalize!(abt)
    return abt
end

function colnormalize!(abt::AbstractComMatrix)
    for j in 1:size(abt, 2)
        colmax = maximum(occurrences(abt)[:,j])
        for i in 1:size(abt, 1)
            occurrences(abt)[i,j] /= colmax
        end
    end
end

function colnormalize(abt::AbstractComMatrix)
    abt = deepcopy(abt)
    rownormalize!(abt)
    return abt
end

function relativeabundance!(a::AbstractComMatrix; kind::Symbol=:fraction)
    in(kind, [:percent, :fraction]) || error("Invalid kind: $kind")
    if eltype(a.occurrences) != Float64; a.occurrences = Float64.(a.occurrences) end
    for i in 1:nsamples(a)
        s = sum(getsample(a,i))
        s == 0 && continue
        for x in 1:nfeatures(a)
            kind == :fraction ? occurrences(a)[x,i] /= s : occurrences(a)[x,i] /= (s / 100.)
        end
    end
end

function relativeabundance(a::AbstractComMatrix; kind::Symbol=:fraction)
    relab = deepcopy(a)
    relativeabundance!(relab, kind=kind)
    return relab
end


text:
using StableRNGs, Random, Test
#include("./functions/SmoothFunctions/smoothFunction.jl")
#include("functions/SmoothFunctions/smoothFunction.jl")

#TODO: needs resolution for windows PC
#=
if Sys.iswindows()
slash =  "\"
else slash = "/";
end

=#
#checking path difference:
windowspath = "src\functions\SmoothFunctions\smoothFunction.jl"
macospath = "/functions/SmoothFunctions/smoothFunction.jl"
#include("src\functions\SmoothFunctions\smoothFunction.jl") "./functions/SmoothFunctions/smoothFunction.jl")

slash = Sys.iswindows() ? slash =  "\" : slash = "/"
#windowspath = "src\functions\SmoothFunctions\smoothFunction.jl"
#macospath = "/functions/SmoothFunctions/smoothFunction.jl"
"""=#

include("src/functions/SmoothFunctions/smoothFunction.jl")
#"/functions/SmoothFunctions/smoothFunction.jl")

function ziggurat(smooth) end

#--- testing area

lobound = 1
upbound = 128
stepSize = 1 # TODO: which steepsize
# 1 is the least Positive integer in Julia
# thus its the best in this contest

#--- RNG

"""
Reproducibility

By using an RNG parameter initialized with a given seed,
you can reproduce the same pseudorandom number sequence

when running your program multiple times. However,
a minor release of Julia (e.g. 1.3 to 1.4)

may change the sequence of pseudorandom numbers
generated from a specific seed,

 in particular if MersenneTwister is used.
 (Even if the sequence produced by a low-level function

 like rand does not change, the output of higher-level
 functions like randsubseq may change due to algorithm updates.)

Rationale: guaranteeing that pseudorandom streams never change
"Prohibits many algorithmic improvements"

If you need to guarantee exact reproducibility of random data, it is advisable to simply save the data (e.g. as a supplementary attachment in a scientific publication). (You can also, of course, specify a particular Julia version and package manifest, especially if you require bit reproducibility.)

Software tests that rely on specific "random" data should also generally either save the data, embed it into the test code, or
use third-party packages like StableRNGs.jl.
 On the other hand, tests that should pass for most random data (e.g. testing A \ (A*x) ≈ x for a random matrix A = randn(n,n))
 can use an RNG, with a fixed seed, to ensure that simply running the test many times does not encounter a failure due to very improbable data (e.g. an extremely ill-conditioned matrix).

The statistical distribution from which random samples are drawn is guaranteed to be the same across any minor Julia releases.
"""

"""
Results from the 'MersenneTwister'
were unsatisfactory, I'm afraid
TODO: delete block once an alternative been found

"""
#NOT RECOMMENDED USING THE MersenneTwister (very bad in number generation)

using Random
seed = 1234; #1= rand(2)# 2 random variables N/a in julia 1.5
rng()
rng = MersenneTwister(seed);
rand(rng, 2) == x1;
reseed = Random.seed!(1234);
Random.seed!(MersenneTwister(seed), seed) # (should be) a reproducible condition
rand(rng, 2) == x1

seed1 = rand(Random.seed!(rng), Bool) # not reproducible

seed2 = rand(Random.seed!(rng), Bool) # (not) reproducible

seed3 = rand(MersenneTwister(), Bool) # not rng has no seed. Infer reproducible either

#---StableRNG - for stability
"""from the docs:
StableRNG is currently an alias for LehmerRNG,
and implements a well understood
linear congruential generator (LCG);
an LCG is not state of the art,
but is fast and is believed to have
reasonably good Statistical Properties [1],
suitable at least for tests of a wide range of packages.
The choice of this particular RNG is based on its simplicity,
 which limits the chances for bugs. Note that only
StableRNG is exported from the package,
& should be the only type used in client code;
LehmerRNG might be renamed, or might be made a
distinct type from StableRNG in any upcoming minor
(i.e. non-breaking) release.
"""
rng = StableRNG(123)
A = randn(rng, 10, 10)
num1 = randn(rng, 1)
num2 = randn(rng, 1) # instead of randn(10, 10)
@test inv(inv(A)) ≈ A
#---range
range1(a = lobound, b = upbound) = lobound:stepSize:length(upbound)
Random.Sampler
rand(range1)

range2(a = lobound, b = upbound) = (lobound in Base.range(lobound, upbound; (upbound - lobound), stepSize))

r = range1(2, 8)

res = rand([rng = GLOBAL_RNG], [r])

#--- Base.://
Function
//(num = 3, den = 5)
"""
Divide two integers or rational numbers, giving a Rational result.
"""
#---range

function range_step_stop_length(a, step, len::Integer)
    start = a - step * (len - oneunit(len))
    if start isa Signed
        # overflow in recomputing length from stop is okay
        return StepRange{typeof(start),typeof(step)}(start, step, convert(typeof(start), a))
    end
    return StepRangeLen{typeof(start),typeof(start),typeof(step)}(start, step, len)
end

for i in range(1000)
    @test(range_step_stop_length)

end
# Stop and length as the only argument
#working
function range_stop_length(a, len::Integer)
    # overflow in recomputing length from stop is okay
    len = length(b - a)

    # overflow in recomputing length from stop is okay
    #return UnitRange(start, oftype(start, a))
    #end

    return StepRangeLen{typeof(start),typeof(start),typeof(step)}(start, step, len)
end
# Stop and length as the only argument

function range_stop_length(a=0,b=100)
# Start and length as the only argument
        # overflow in recomputing length from stop is okay
		 len = length(b - a)

	type_b = typeof(b)

# double return

#len = length(b - a);

#working
function range_stop_length(a = 0, b = 100)
    # Start and length as the only argument
     len = length(b - a)
    type_b = typeof(b)

    return StepRangeLen{type_b,typeof(a)}(a, step, len)
end

using Distributions

function algorithm2()
    """Rejection sampling method : 3 step process
    does PRN means RNG? (maybenot, more likely a whole PDF)
    idea: probability distribution P(x) is subsumed by a set of boxes, resembling a ziggurat

   1. Generate 2 PRNs x, y #from a uniform distribution , together, form a point

   2. If y < P(x) return x # accepted retrun true too if possible

   3. else point rejected
   then go back & draw another point (select & test )

   --------------
In a ziggurat algorithm,
the desired probability distribution P(x) lies beneath a
stack of rectangular layers. Layers are designed such that
they all contain the exact same area, so that each box can be randomly
chosen with equal probability using a uniform random integer i
to ensure uniform coverage of P(x). The height fi and length Xi of
each ziggurat layer are pre-calculated in lookup tables.
 Because of these lookup tables, ziggurat algorithms are most
efficiently implemented on systems with large caches (e.g. modern CPUs,
but not current GPUs)

X

  """
using Distributions
  i = Uniform()

function algorithm1()
    """
    1. Generate U1, i #from a uniform distribution
    2. If i < Lmax, return U1Xi
    3. Generate j from A
    4. Ifj=0, returna value from the tail 5. Generate U2
    6. x ← Xj + U1(Xj−1 − Xj) # Q1. what is U1(Xj−1 − Xj), again?
    7. If (fi−1 − fi)U2 < P(x), return x
    8. Go to 4.
    """
    i = Distributions.Uniform(0, 1)
    U1 = Uniform(0, 1)
    #3. Generate j from A # Q.what is lmax?
    rng = StableRNG(123)
    Lmax = [ rand(rng, 1)
    if (i < Lmax)
        return U1Xi #returns a max  value
    end #otherwise  create a random array A
    A = smooth()  # randn(rng, 0, 1) #generate A matrix , via a stable RNG from 0 to 1
    j = A[rand(1:end)] #randn(rng, 1) #pick 1 value for R  j from A  #done

for i in enumerate ([i in rand(0,1,1])
#DEFINE X
j = A[rand(1:end)] #pick j randomly (from a Uniform distribution)
#https://discourse.julialang.org/t/generate-random-value-from-a-given-function-out-of-box/5793/4
#using StatsBase, Gadfly

x = linspace(0,π,100)
P = smooth(x)     # pdf
P = P/sum(P)   #cdf

#r = StatsBase.sample(1:100, Weights(P),10000)
plot(x=x[r],Geom.histogram)

if j isa 0
        return returnValfromTrail()

        If (fi−1 − fi) U2 < P(x) return x end

    end
end

"""
total number of layers in the modified algorihtm Lmax cannot be determined,
until the last layer is calculated,
"""
#testing area

x = X[j] + U1(Xj−1 − Xj)

ε = max(x ,[fi +(x−X[i])*(fi−1 −fi)−P(x)] )

ε=f[i−1] +Log(f[i] −f[i−1])+Xi)*(fi −fi−1))

# end

n = 1000000000
for j in numerate(n)

end

end

function returnValfromTrail() end

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using BitBasis
using Documenter
using DocThemeIndigo

indigo = DocThemeIndigo.install(BitBasis)

makedocs(
    modules = [BitBasis],
    repo = "https://github.com/QuantumBFS/BitBasis.jl/blob/{commit}{path}#{line}",
    format = Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
        canonical = "https://quantumbfs.github.io/BitBasis.jl/dev/",
        assets = String[indigo, "assets/favicon.ico"],
    ),
    sitename = "BitBasis.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = ["Home" => "index.md", "Tutorial" => "tutorial.md", "Manual" => "man.md"],
)

deploydocs(repo = "github.com/QuantumBFS/BitBasis.jl.git")


text:
using Distributed

addprocs(Sys.CPU_THREADS ÷ 2)

@everywhere using ElectrostaticPIC1D, FFTW, JLD2, Plots, Random, SpecialFunctions, Test, ThreadsX
@everywhere Random.seed!(0)

@everywhere function go(;NG=256, NPPCPS=16, L=1.0, nvortexpergrid=1,
    cleandatadir=true, dryrunabort=false, datadir="data")
#  NG = 64 # number of grid points
#  NPPCPS = 16 # number of particles per cell per species
#  nvortexpergrid = 1 # number of vortices per grid, should really be an Int
#  L = 1.0 # length of domain

  Δ = L / NG
  NS = 2 # number of species
  NP = NPPCPS * NG * NS # total number of particles
  q0 = me = Te = 1 # electron charge, mass and temperature are all unity
  mi = 1836 * 2 * me
  Ti = Te
  # unit temperature, unit electron mass, so Debye length is 2π / Πₑ
  # Set Debye length = Δ, so Πₑ = 1 / Δ
  debyelength = Δ
  electronplasmafrequency = sqrt(Te / me) / debyelength
  density = electronplasmafrequency^2 / q0^2 * me # plasma number density for electrons
  weight = density * L / (NPPCPS * NG) # weight of each particle

  electronnuclide = Nuclide(-q0, me) # electron
  ionnuclide = Nuclide(q0, mi) # ion

  cleandatadir && run(`rm -rf $datadir`)

  particleshapes = ((DeltaFunctionShape(), "Delta"),
                    (GaussianShape(Δ), "Gaussian"),
                    #(BSpline{0}(Δ), "BSpline0"),
                    #(BSpline{1}(Δ), "BSpline1"),
                    #(BSpline{2}(Δ), "BSpline2"),
                    )

  fieldsolvers = (#(FourierField(NG,L), "Fourier"),
                  (LSFEMField(NG,L,BSpline{1}(Δ)), "LSFEM_BSpline1"),
                  (LSFEMField(NG,L,BSpline{2}(Δ)), "LSFEM_BSpline2"),
                  (LSFEMField(NG,L,GaussianShape(Δ * √2)), "LSFEM_Gaussian"),
                  #(GalerkinFEMField(NG,L,BSpline{1}(Δ), BSpline{2}(Δ)), "Galerkin_BSpline1_BSpline2"),
                  #(FiniteDifferenceField(NG,L,order=1), "FiniteDifference1"),
                  #(FiniteDifferenceField(NG,L,order=2), "FiniteDifference2"),
                  #(FiniteDifferenceField(NG,L,order=4), "FiniteDifference4"),
                  )

  particleics = Dict()

  x1 = (bitreverse.(0:NPPCPS * NG-1) .+ 2.0^63) / 2.0^64
  xe = mod.(x1, L)
  xi = mod.(x1, L)
  ve = erfinv.(2 .* mod.(x1 .+ π, 1) .- 1) .* √(Te / me)
  vi = erfinv.(2 .* mod.(x1 .+ π, 1) .- 1) .* √(Ti / mi)
  @assert all(isfinite, xe)
  @assert all(isfinite, xi)
  @assert all(isfinite, ve)
  @assert all(isfinite, vi)

  particleics[:Quiet] = deepcopy((xe, ve, xi, vi))

  function innerloop(inputs)

    particleshape, particletypename = inputs[1]

    field, fieldtypename = inputs[2]

    icname, particleics = inputs[3]
    xe, ve, xi, vi = particleics

    run(`mkdir -p $datadir/$fieldtypename/$particletypename/$icname`)
    stub = "$(datadir)/$(fieldtypename)/$(particletypename)/$(icname)/"

    electrons = Species([Particle(electronnuclide, particleshape; x=xe[i], v=ve[i], w=weight)
                    for i in eachindex(xe)])
    ions = Species([Particle(ionnuclide, particleshape; x=xi[i], v=vi[i], w=weight)
                    for i in eachindex(xi)])

    plasma = Plasma([electrons, ions])

    ti = LeapFrogTimeIntegrator(plasma, field; cflmultiplier=1/2)

    expectedenergy = weight * NP / 2
    expectedmomentum = 0.0
    chargedensitynormalisation = weight * NP * (abs(charge(electronnuclide)) + abs(charge(ionnuclide)))
    electricfieldnormalisation = chargedensitynormalisation * L / NG

    x = cellcentres(field)
    dryrunabort && return nothing

    try
      sim = Simulation(plasma, field, ti, diagnosticdumpevery=2,
                       endtime=10.0, filenamestub=stub)

      init!(sim) # set up to start
      @time run!(sim)

      fileindices = 0:ElectrostaticPIC1D.diagnosticdumpcounter(sim)
      Ex = zeros(length(x), length(fileindices))
      Rho = zeros(length(x), length(fileindices))
      fieldenergies = []
      particleenergies = []
      particlemomenta = []
      times = []
      timesdr = []

      for (i, fileindex) in enumerate(fileindices)
        sim_i = ElectrostaticPIC1D.load(stub, fileindex)

        Ex[:, i] .= sim_i.field.electricfield.(x) ./ electricfieldnormalisation
        Rho[:, i] .= sim_i.field.chargedensity.(x) ./ chargedensitynormalisation

        push!(timesdr, sim_i.time[])

        i % 10 == 0 || continue

        push!(times, sim_i.time[])
        particleenergy = energy(sim_i.plasma)
        fieldenergy = energydensity(sim_i.field) * L
        particlemomentum = momentum(sim_i.plasma)

#        plot!(x, Ex[:, i], label="E")
#        plot!(x, Rho[:, i], label="ρ")
        push!(fieldenergies, fieldenergy)
        push!(particleenergies, particleenergy)
        push!(particlemomenta, particlemomentum)
      end
      @save stub * "energydata.jld2" fieldenergies particleenergies expectedenergy times
      @save stub * "momentumdata.jld2" particlemomenta times
      totalenergies = fieldenergies + particleenergies

      y = fieldenergies / totalenergies[1]
      plot(times, log10.(y), label="field energy")

      y = particleenergies / totalenergies[1]
      plot!(times, log10.(abs.(y)), label="particle energy")

      y = (fieldenergies .+ particleenergies) / expectedenergy
      plot!(times, log10.(abs.(y .- y[1])), label="energy error")

      y = particlemomenta
      plot!(times, log10.(abs.(y .- y[1])), label="momentum")

      xlabel!("Time", legend=:bottomright)
      xlims!(0, maximum(times))
      ylims!(-17, 1.0)
      savefig(stub * "plot.pdf")

      for (Z, fieldstring) ∈ ((Ex, "Ex"), (Rho, "Rho"))
        LogAbsFFTEx = log10.(abs.(fft(Z')));
        LogAbsFFTEx[:, 1] *= NaN; # zeros on the log10 colorbar destroy the contrast
        LogAbsFFTEx[:, (end÷2) + 1 ] *= NaN;
        ks = (-NG÷2 : NG÷2 - 1)
        NT = length(timesdr)
        NT_2 = Int(floor(NT/2))

        ws = (0:NT_2-1) ./ (timesdr[end] - timesdr[1]) / (electronplasmafrequency / 2π)
        Z1 = fftshift(LogAbsFFTEx)[end-length(ws)+1:end, :]
        heatmap(ks, ws, Z1)
        xlabel!("Wavenumber [Waves per box]")
        ylabel!("Frequency [Π]")
        ylims!(0, min(3, maximum(ws)))
        xlims!(ks[1], ks[end])
        savefig(stub * "$(fieldstring)_dispersionrelation.pdf")
      end
      @save stub * "dispersionrelation.jld2" Ex Rho timesdr x L electronplasmafrequency

    catch e
      @show fieldtypename, particletypename
      @show e
      rethrow(e)
    end

    return nothing
  end

  inputs = Iterators.product(particleshapes, fieldsolvers, particleics)

  map(innerloop, inputs)

end

@testset "Simulations" begin
  go(; datadir="data_maxwellian")
end


text:
# ------------------------------------
# Tests adapted from Git2Go Library
# ------------------------------------
context("test create / lookup ref") do
    test_path = joinpath(pwd(), "testrepo")
    repo = create_test_repo(test_path)
    try
        cid, tid = seed_test_repo(repo)

        _ = create_ref(repo, "refs/tags/tree", tid, force=true)
        tag = lookup_ref(repo, "refs/tags/tree")
        @test isa(tag, GitReference{Oid})

        ref = lookup_ref(repo, "HEAD")
        @test isa(ref, GitReference{Sym})

        @test target(ref) == nothing
        ref = resolve(ref)
        @test isa(ref, GitReference{Oid})
        @test isa(target(ref), Oid)
        @test symbolic_target(ref) == ""
        @test hex(cid) == hex(target(ref))

        _ = rename(tag, "refs/tags/renamed", force=false)
        tag = lookup_ref(repo, "refs/tags/renamed")
        @test isa(tag, GitReference{Oid})
    finally
        close(repo)
        LibGit2.free!(repo)
        Base.gc()
        rm(test_path, recursive=true)
    end
end

#=
context() do
    test_path = joinpath(pwd(), "testrepo")
    repo = create_test_repo(test_path)
    try
        cid, tid = seed_test_repo(repo)

        sig = Signature("test", "test@test.com")
        idx = GitIndex(repo)
        add!(idx, "README")
        tid = write_tree!(idx)

        message = "This is a commit\n"
        tree = lookup_tree(repo, tid)

        cid = commit(repo, "HEAD", sig, sig, message, tree)

        _ = create_ref(repo, "refs/heads/one",   cid, force=true)
        _ = create_ref(repo, "refs/heads/two",   cid, force=true)
        _ = create_ref(repo, "refs/heads/three", cid, force=true)

        expected = [join(["refs/heads", x], "/")
                    for x in ["master","one","two","three"]]
        test_names = AbstractString[]
        for r in iter_refs(repo)
            push!(test_names, name(r))
        end
        sort!(expected)
        sort!(test_names)
        for (exp, tst) in zip(expected, test_names)
            @test exp == tst
        end
        # test glob
        expected = ["refs/heads/two", "refs/heads/three"]
        test_names = AbstractString[]
        for r in iter_refs(repo, "refs/heads/t*")
            push!(test_names, name(r))
        end
        sort!(expected)
        sort!(test_names)
        for (exp, tst) in zip(expected, test_names)
            @test exp == tst
        end
    finally
        close(repo)
        LibGit2.free!(repo)
        Base.gc()
        rm(test_path, recursive=true)
    end
end
=#

# -----------------------------------------
# Tests adapted from Ruby's Rugged Library
# -----------------------------------------
with_repo_access() do test_repo, path
    context("test reference validity") do
        valid = "refs/foobar"
        invalid = "refs/nope^*"
        @test is_valid_ref(valid) == true
        @test is_valid_ref(invalid) == false
    end

    # test can handle exceptions

    context("test list references") do
        tmp = map((r) -> replace(name(r), "refs/", ""), foreach(GitReference, test_repo))
        @test join(sort(tmp), ":") == "heads/master:heads/packed:notes/commits:tags/v0.9:tags/v1.0"
    end

    context("test can filter refs with regex") do
        tmp = map((r) -> replace(name(r), "refs/", ""), foreach(GitReference, test_repo, "refs/tags/*"))
        refs = join(sort(tmp), ":")
        @test refs == "tags/v0.9:tags/v1.0"
    end

    context("test can filter refs with string") do
        tmp = map((r) -> replace(name(r), "refs/", ""), foreach(GitReference, test_repo, "*0.9*"))
        refs = join(sort(tmp), ":")
        @test refs == "tags/v0.9"
    end
end

with_repo_access() do test_repo, path

  context("test can open reference") do
    ref = lookup_ref(test_repo, "refs/heads/master")
    @test target(ref) == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test isa(ref, GitReference{Oid})
    @test name(ref) ==  "refs/heads/master"
    @test peel(ref) == nothing
  end

  context("test can open a symbolic reference") do
    ref = lookup_ref(test_repo, "HEAD")
    @test symbolic_target(ref) == "refs/heads/master"
    @test isa(ref, GitReference{Sym})

    resolved = resolve(ref)
    @test isa(resolved, GitReference{Oid})
    @test target(resolved) == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test target(resolved) == peel(ref)
  end

  context("test looking up missing ref returns nil") do
    ref = lookup_ref(test_repo, "lol/wut")
    @test ref == nothing
  end

  context("test load reflog") do
      ref = lookup_ref(test_repo, "refs/heads/master")
      #@test has_reflog(ref) == true
      rlog = GitReflog(ref)
      entry = rlog[2]
      @test isa(entry, GitReflogEntry)
      @test entry.id_old == Oid("8496071c1b46c854b31185ea97743be6a8774479")
      @test entry.id_new == Oid("5b5b025afb0b4c913b4c338a42934a3863bf3644")
      @test entry.message == "commit: another commit"
      @test isa(entry.committer, Signature)
      @test email(entry.committer) == "schacon@gmail.com"
  end

  context("test reference exists") do
    @test exists(test_repo, "refs/heads/master") == true
    @test exists(test_repo, "lol/wut") == false
  end

  context("test load packed ref") do
    ref = lookup_ref(test_repo, "refs/heads/packed")
    @test target(ref) == Oid("41bc8c69075bbdb46c5c6f0566cc8cc5b46e8bd9")
    @test isa(ref, GitReference{Oid})
    @test name(ref) == "refs/heads/packed"
  end

  context("test resolve head") do
    ref = lookup_ref(test_repo, "HEAD")
    #TODO: make target work for both direct and symboloic
    @test symbolic_target(ref) == "refs/heads/master"
    @test isa(ref, GitReference{Sym})

    head = resolve(ref)
    @test target(head) == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test isa(head, GitReference{Oid})
  end

  context("test reference to tag") do
    ref = lookup_ref(test_repo, "refs/tags/v1.0")
    @test target(ref) == Oid("0c37a5391bbff43c37f0d0371823a5509eed5b1d")
    @test peel(ref) == Oid("5b5b025afb0b4c913b4c338a42934a3863bf3644")
  end
end

with_tmp_repo_access("test create ref") do test_repo, path
    create_ref(test_repo,
               "refs/heads/unit_test",
               "refs/heads/master")

    create_ref(test_repo,
               "refs/heads/unit_test",
               "refs/heads/master",
               force=true)
end

with_tmp_repo_access("test create ref with unicode") do test_repo, path
    UNICODE_REF_NAME = "A\314\212ngstro\314\210m"
    create_ref(test_repo,
               "refs/heads/$UNICODE_REF_NAME",
               "refs/heads/master")
    refs = map(r -> normalize_string(replace(name(r), "refs/", ""), :NFC), foreach(GitReference, test_repo))
    @test normalize_string("heads/$UNICODE_REF_NAME", :NFC) in refs
end

with_tmp_repo_access("test lookup non-existant ref") do test_repo, path
    ref = lookup_ref(test_repo, "lol/wut")
    @test ref == nothing
end

with_tmp_repo_access("test name ref") do test_repo, path
    @test realpath(workdir(test_repo)) == realpath(path)

    o = Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    ref = create_ref(test_repo, "refs/heads/unit_test", o)

    @test o == target(ref)
    @test isa(ref, GitReference{Oid})
    @test name(ref) == "refs/heads/unit_test"
    #delete!(test_repo, ref)
end

with_tmp_repo_access("test rename ref") do test_repo, path
    ref = create_ref(test_repo,
      "refs/heads/unit_test",
      Oid("36060c58702ed4c2a40832c51758d5344201d89a"))

    @test target(ref) == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test isa(ref, GitReference{Oid})
    @test name(ref) == "refs/heads/unit_test"

    new_ref = rename(ref, "refs/heads/rug_new_name")
    @test name(new_ref) == "refs/heads/rug_new_name"
    #delete!(test_repo, new_ref)
end

with_tmp_repo_access("test set target ref") do test_repo, path
    ref = create_ref(test_repo,
                     "refs/heads/unit_test",
                     Oid("36060c58702ed4c2a40832c51758d5344201d89a"))

    @test target(ref) == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test isa(ref, GitReference{Oid})
    @test name(ref) == "refs/heads/unit_test"
    new_ref = set_target(ref, Oid("5b5b025afb0b4c913b4c338a42934a3863bf3644"))
    t = target(new_ref)
    @test t == Oid("5b5b025afb0b4c913b4c338a42934a3863bf3644")
    #delete!(test_repo, new_ref)
end

with_tmp_repo_access("test create symbolic target") do test_repo, path
    ref1 = create_ref(test_repo, "refs/heads/Ångström", "refs/heads/master")
    ref2 = create_ref(test_repo, "refs/heads/foobar", "refs/heads/Ångström")
    @test name(ref1) == "refs/heads/Ångström"
    @test symbolic_target(ref2) ==  "refs/heads/Ångström"
end

with_tmp_repo_access("test reflog") do test_repo, path
    ref = create_ref(test_repo,
                     "refs/heads/test-reflog",
                     Oid("36060c58702ed4c2a40832c51758d5344201d89a"))
    log!(ref, "", Signature("foo", "foo@bar"))
    log!(ref, "commit: bla bla", Signature("foo", "foo@bar"))
    rlog = GitReflog(ref)
    # TODO: this fails for travis (cannot reproduce in local tests)
    # travis states that length(rlog) == 3 with fist reflog entry
    # containing no information (unknown user/email)
    @test length(rlog) == 3

    #@test rlog[2].id_old == Oid("0000000000000000000000000000000000000000")
    @test rlog[2].id_old == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[2].id_new == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[2].message == ""
    @test name(rlog[2].committer) == "foo"

    @test rlog[end].id_old == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[end].id_new == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[end].message == "commit: bla bla"
    @test name(rlog[end].committer) == "foo"
    @test email(rlog[end].committer) == "foo@bar"
end

with_tmp_repo_access("test ref log with config") do test_repo, path
    ref = create_ref(test_repo,
                     "refs/heads/test-reflog",
                     Oid("36060c58702ed4c2a40832c51758d5344201d89a"))
    testname = "Julia User"
    testemail = "user@julia.com"

    cfg = GitConfig(test_repo)
    cfg["user.name"] = testname
    cfg["user.email"] = testemail

    log!(ref)
    log!(ref, "commit: bla bla")

    rlog = GitReflog(ref)
    # TODO: this fails for travis (cannot reproduce in local tests)
    # travis states that length(rlog) == 3 with fist reflog entry
    # containing no information (unknown user/email)
    @test length(rlog) == 3

    #@test rlog[end-1].id_old == Oid("0000000000000000000000000000000000000000")
    @test rlog[2].id_old == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[2].id_new == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[2].message == ""
    @test name(rlog[2].committer) == "Julia User"
    @test email(rlog[2].committer) == "user@julia.com"

    @test rlog[end].id_old == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[end].id_new == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test rlog[end].message == "commit: bla bla"
    @test name(rlog[end].committer) == "Julia User"
    @test email(rlog[end].committer) == "user@julia.com"
end


text:
include("required.jl")

## Conditional arguments
# outcome = "raw"
# outcome = "bart"
# save_additional_figures = true

# Make sure "outcome" is defined
if !(@isdefined outcome) || !(outcome in ["raw", "bart"])
    @error "'outcome' must be either 'raw' or 'bart'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Species occupancy

# Prep values
include("04_full-extent.jl")
include("05_subareas.jl")

# Count sites with occurrences per species
function getoccupancy(Yobs; count_sites=false)
    occupancy = vec(sum(Yobs; dims=1))
    if !count_sites
        occupancy = occupancy ./ size(Yobs, 1)
    end
    return occupancy
end
occupancy_count = getoccupancy(Yobs; count_sites=true)
occupancy = getoccupancy(Yobs)
occupancy_NE = getoccupancy(Yobs_NE)
occupancy_SW = getoccupancy(Yobs_SW)

# Visualize
scatter(
    occupancy;
    formatter=:plain,
    legend=:none,
    xlabel="Species ID number",
    ylabel="Species occupancy",
)

# Classify as rare species
# Based on local occupancy
function getrarespecies(occupancy::Vector{<:AbstractFloat}, threshold::AbstractFloat)
    @assert all(x -> 0.0 <= x <= 1.0, occupancy) "occupancy must be between 0 and 1"
    rarespecies = Vector{Union{Missing,Int64}}(undef, length(occupancy))
    for (i, op) in enumerate(occupancy)
        if iszero(op)
            rarespecies[i] = missing
        else
            rarespecies[i] = op < threshold ? 1 : 0
        end
    end
    return rarespecies
end
# Based on total occupancy
function getrarespecies(
    occupancy::Vector{<:T}, total_occupancy::Vector{<:T}, threshold::T
) where {T<:AbstractFloat}
    @assert all(x -> 0.0 <= x <= 1.0, occupancy) "occupancy must be between 0 and 1"
    @assert all(x -> 0.0 <= x <= 1.0, total_occupancy) "occupancy must be between 0 and 1"
    rarespecies = Vector{Union{Missing,Int64}}(undef, length(occupancy))
    for (i, op) in enumerate(occupancy)
        if iszero(op)
            rarespecies[i] = missing
        else
            rarespecies[i] = total_occupancy[i] < threshold ? 1 : 0
        end
    end
    return rarespecies
end

# Local occupancy
threshold = 0.4
rarespecies = getrarespecies(occupancy, threshold)
rarespecies_NE = getrarespecies(occupancy_NE, 0.4)
rarespecies_SW = getrarespecies(occupancy_SW, 0.4)

# Total occupancy
rarespecies_NE_total = getrarespecies(occupancy_NE, occupancy, 0.4)
rarespecies_SW_total = getrarespecies(occupancy_SW, occupancy, 0.4)

# Percentage of rarespecies
rarespecies_p = mean(skipmissing(rarespecies))
rarespecies_p_NE = mean(skipmissing(rarespecies_NE))
rarespecies_p_SW = mean(skipmissing(rarespecies_SW))
rarespecies_p_NE_total = mean(skipmissing(rarespecies_NE_total))
rarespecies_p_SW_total = mean(skipmissing(rarespecies_SW_total))

# Wrap all steps as single function
function get_rarespecies_p(Yobs, threshold)
    occupancy = getoccupancy(Yobs)
    rarespecies = getrarespecies(occupancy, threshold)
    rarespecies_p = mean(skipmissing(rarespecies))
    return rarespecies_p
end
get_rarespecies_p(Yobs, threshold)

# LCBD-richness correlation
get_eusrr(richness, lcbd) = corspearman(collect(richness), collect(lcbd))
eusrr = get_eusrr(richness_full, lcbd_full)
eusrr_NE = get_eusrr(richness_NE, lcbd_NE)
eusrr_SW = get_eusrr(richness_SW, lcbd_SW)

## Plot EUSRR ~ rare species percentage
# Similar to Fig. 3 of Yao et al. 2021

# Subareas & rarity based on local occupancy
scatter(
    [rarespecies_p rarespecies_p_NE rarespecies_p_SW],
    [eusrr eusrr_NE eusrr_SW];
    xlabel="Percentage of rare species (%)",
    ylabel="EUSRR",
    xlim=(0.40, 1.0),
    ylim=(-1.0, 1.0),
    label=["Total" "NE" "SW"],
    legendtitle="Regions",
    legendtitlefontsize=9,
) |> x -> hline!(x, [0.0]; style=:dash, c=:grey, label=:none)
# EUSRR negative whatever subarea & percentage of rare species?

# Subareas & rarity based on total occupancy
scatter(
    [rarespecies_p rarespecies_p_NE_total rarespecies_p_SW_total],
    [eusrr eusrr_NE eusrr_SW];
    xlabel="Percentage of rare species (%)",
    ylabel="EUSRR",
    xlim=(0.40, 1.0),
    ylim=(-1.0, 1.0),
    label=["Total" "NE" "SW"],
    legendtitle="Regions",
    legendtitlefontsize=9,
) |> x -> hline!(x, [0.0]; style=:dash, c=:grey, label=:none)
# EUSRR negative whatever subarea & percentage of rare species?

# Effect of thresholds
# Thresholds variation only
thresholds = [0.1, 0.2, 0.3, 0.4, 0.5]
[get_rarespecies_p(Yobs, t) for t in thresholds]
# Thresholds and region
Ys = [Yobs, Yobs_NE, Yobs_SW]
rarespecies_matrix = [get_rarespecies_p(y, t) for t in thresholds, y in Ys]

# Subareas & thresholds
p_eusrr =
    scatter(
        rarespecies_matrix,
        [eusrr eusrr_NE eusrr_SW];
        xlabel="Percentage of rare species (%)",
        ylabel="EUSRR",
        xlim=(0.0, 1.0),
        ylim=(-1.0, 1.0),
        label=["Total" "NE" "SW"],
        legendtitle="Regions",
        legendtitlefontsize=9,
    ) |> x -> hline!(x, [0.0]; style=:dash, c=:grey, label=:none)
# EUSRR negative whatever subarea, percentage of rare species & threshold
# Whatever the threshold, EUSRR negative
if (@isdefined save_additional_figures) && save_additional_figures == true
    savefig(
        p_eusrr, joinpath("fig", outcome, "06_$(outcome)_eusrr.png")
    )
end

p_thresholds = plot(
    thresholds,
    rarespecies_matrix;
    xlabel="Rare species threshold",
    ylabel="Rare species percentage",
    labels=["Total" "NE" "SW"],
    legend=:bottomright,
)

if (@isdefined save_additional_figures) && save_additional_figures == true
    savefig(
        p_thresholds, joinpath("fig", outcome, "06_$(outcome)_thresholds.png")
    )
end

## Effect of scaling on rarespecies

# Scaling EUSRR
rarespecies_scaling = Vector{Float64}()
eusrr_scaling = Vector{Float64}()
for sc in subarea_coords
    distributions_sc = [d[sc] for d in distributions]
    Y_sc = Ymatrix(distributions_sc)
    Yobs_sc = Ymatrix(distributions_sc; observed=true)
    richness_sc = richness(Y_sc, distributions_sc[1])
    lcbd_sc = lcbd(Y_sc, distributions_sc[1]; relative=false)

    rarespecies_sc = get_rarespecies_p(Yobs_sc, threshold)
    eusrr_sc = get_eusrr(richness_sc, lcbd_sc)

    push!(rarespecies_scaling, rarespecies_sc)
    push!(eusrr_scaling, eusrr_sc)
end
rarespecies_scaling
eusrr_scaling

# Plot scaling EUSRR
scatter(
    rarespecies_scaling,
    eusrr_scaling;
    xlabel="Percentage of rare species (%)",
    ylabel="EUSRR",
    xlim=(0.0, 1.0),
    ylim=(-1.0, 1.0),
    label=["Total" "NE" "SW"],
    legendtitle="Regions",
    legendtitlefontsize=9,
) |> x -> hline!(x, [0.0]; style=:dash, c=:grey, label=:none)
# Huh ok...
show(stdout, "text/plain", rarespecies_scaling)
show(stdout, "text/plain", eusrr_scaling)
# Whatever

## Spatial distribution of rare species proportion

# Rarespecies proportion per site
function get_site_rarespecies(Yobs, rarespecies, richness_layer)
    if any(ismissing, rarespecies)
        rarespecies = replace(rarespecies, missing => 0)
    end
    Yrare = Yobs + repeat(permutedims(rarespecies), size(Yobs, 1)) .- 1.0
    Yrare_p = sum(isone, Yrare; dims=2) ./ sum(Yobs; dims=2)

    rarespecies_layer = copy(richness_layer)
    inds_obs = findall(!isnothing, rarespecies_layer.grid)
    rarespecies_layer.grid[inds_obs] .= vec(Yrare_p)
    return rarespecies_layer
end

# Full extent
rarespecies_layer = get_site_rarespecies(Yobs, rarespecies, richness_full)
plot_layer(rarespecies_layer; c=:viridis)

# NE local rarity
rarespecies_layer_NE = get_site_rarespecies(Yobs_NE, rarespecies_NE, richness_NE)
p1 = plot_layer(rarespecies_layer_NE; c=:viridis)

# NE total rarity
rarespecies_layer_NE_total = get_site_rarespecies(Yobs_NE, rarespecies_NE_total, richness_NE)
p2 = plot_layer(rarespecies_layer_NE_total; c=:viridis)

# SW local rarity
rarespecies_layer_SW = get_site_rarespecies(Yobs_SW, rarespecies_SW, richness_SW)
p3 = plot_layer(rarespecies_layer_SW; c=:viridis)

# SW total rarity
rarespecies_layer_SW_total = get_site_rarespecies(Yobs_SW, rarespecies_SW_total, richness_SW)
p4 = plot_layer(rarespecies_layer_SW_total; c=:viridis)

p_subareas = plot(p1, p2, p3, p4; dpi=200)

## Ascending & descending parts

# Check relationship plots
combined_plot

# Separate ascending & descending portions in plots
function ascending_plots(richness, lcbd, rarespecies)
    # Choose threshold
    ascending_threshold = minimum(lcbd)
    min_indx = findall(x -> x == ascending_threshold, lcbd.grid)
    abs_min = median(richness.grid[min_indx])

    # Apply threshold
    binlayer = replace(richness, Pair.(unique(richness), unique(richness) .> abs_min)...)

    # Thresholded rarity map subpanel
    lplot = plot_layer(binlayer; c=cgrad(:PuOr; rev=true), colorbar=:none)

    # Create two group labels
    bin_names = map(x -> isone(x) ? "Ascending" : "Descending", collect(binlayer))
    isascending = isequal.(bin_names, "Ascending")

    # Select occurrence values only
    rare_values = collect(rarespecies)

    # Density subpanel
    dplot = density(rare_values[isascending]; c=:PuOr, label="Ascending")
    density!(rare_values[.!isascending]; c=cgrad(:PuOr; rev=true), label="Descending")
    plot!(; xlabel="Rare species percentage", ylabel="Probability Density", bottommargin=4.0mm)

    l = @layout [a{0.5w} b{0.45w}]
    return plot(lplot, dplot; size=(900, 300), layout=l)
end

# Create subarea plots
p_asc1 = ascending_plots(richness_NE, lcbd_NE, rarespecies_layer_NE)
p_asc2 = ascending_plots(richness_SW, lcbd_SW, rarespecies_layer_SW)

# Fix axes & dimensions
yticks!(p_asc1[1], 40:2:50)
yticks!(p_asc2[1], 30:2:40)
xticks!(p_asc1[1], -80:4:-60)
xticks!(p_asc2[1], -120:4:-100)

# Combine plots
asc_plots = plot(
    p_asc1,
    p_asc2;
    layout=(2, 1),
    size=(900, 660),
    title=["a) Northeast subregion" "" "b) Southwest subregion" ""],
    titleloc=:left,
    dpi=600
)

# Export plot
# save_additional_figures = true
if (@isdefined save_additional_figures) && save_additional_figures == true
    savefig(
        asc_plots,
        joinpath("fig", outcome, "06_$(outcome)_rare-species.png"),
    )
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
__precompile__(true)
module OMRemote

  # (C) Christian Kral, 2017-2019, BSD-3 license
  # OMRemote remote access to run OpenModelica simulations from Julia
  # https://gitlab.com/christiankral/RemoteOM.git

  using OMJulia
  using Base

  export envOM, simulateOM

  function envOM(env...)
    # Set environment variables
    if length(env) == 0
      println("Setting environment for OpenModelica: OPENMODELICAHOME=/usr/")
      ENV["OPENMODELICAHOME"] = "/usr/"
    else
      println("Setting environment for OpenModelica: OPENMODELICAHOME="*env)
      ENV["OPENMODELICAHOME"] = env[1]
    end
  end

  """
  # Function call

  resFile = `simulateOM(model;sysLibs="Modelica",
  sysVers = "", workFiles = "", files = "",
  workDir = "/work", simDir = "/tmp/OpenModelica",
  StartTime = 0, StopTime = 0, Tolerance = 1E-6, copy = false)`

  # Description

  This function simulates and OpenModelica model and returns the file name
  of the result file

  # Variables

  `model` Modelica experiment to be simulated (incl. dot notation)

  `sysLibs` Modelica system libraries included by OpenModelica; default
            is the Modelica Standard Library (MLS) "Modelica"; several
            libraries may be separated by `:`

  `sysVers` Versions of system libaries (::String); the number of versions seperated
            by `:` has to be equal to the number of system libraries seperated
            by `:`

  `workFiles` Modelica files to be loaded, which, e.g., contain `model`; the
              Modelica files have to be located in in `workDir`; several
              files may be separated by `:`

  `files` Additional Modelica files to be loaded; the absolute file path has
          to be provided in order to read the file properly; several
          files may be separated by `:`

  `workDir` Working directory; default value = "/work"

  `simDir` Simulation directory, default value = "/tmp/OpenModelica"

  `StartTime` Start time of simulation; default value = 0

  `StopTime` Stop time of simulation; if the StopTime <= StartTime,
             the original setting from the Experiment
             annotation is used; default value = 0

  `Tolerance` Simulation tolerance; default value = 1E-6

  `NumberOfIntervals` Number of output intervals; default value = 500

  `copy` Copy simulation result created by OpenModelica to directory `workDir`,
         if `true`; default value = false; in case of `true` the copied file
         name is assigned to the return variable `resFile`, otherwise the
         original file location created by OpenModelica is assigned to `resFile`
  """
  function simulateOM(model;sysLibs="Modelica",sysVers="",workFiles="",files="",
                      workDir="/work",simDir="/tmp/OpenModelica",
                      StartTime=0,StopTime=0,Tolerance=1E-6,
                      NumberOfIntervals=500, copy=true)

    localDir=pwd()
    # Load MSL and libraries (separated by :) and simulate model
    # Result determined by DyMat is return argument

    # Set environment variable for system libraries
    envOM()

    # Exstension of result file, created by OpenModelica
    resExt = "_res.mat"
    rm(simDir,force=true,recursive=true)
    mkdir(simDir)
    cd(simDir)
    omc = OMJulia.OMCSession()
    # Change directory in OM
    OMJulia.sendExpression(omc, "cd(\""*simDir*"\")")

    # Set system directory
    # sysDir = ENV["OPENMODELICAHOME"]*"/"*sysdir*"/"

    # Load system libraries
    sysLibsList = split(sysLibs,":")
    sysVersList = split(sysVers,":")
    for k in collect(1:length(sysLibsList))
      # Determine name of system library
      sysLib = sysLibsList[k]
      if length(sysLibsList)>length(sysVersList)
        # Not every sysLib entry has a version number
        # In this case the version number is not treated at all
        if sysLib!=""
          print("  loadModel("*sysLib*")")
          OMJulia.sendExpression(omc, "loadModel("*sysLib*")")
          if status==true
            println(" successful")
          else
            println(" failed")
          end
        end
      elseif length(sysLibsList)==length(sysVersList)
        # Every sysLib entry has a version number
        # In this case the version number is treated
        # Determine version of system library
        sysVer = sysVersList[k]
        if sysLib!=""
          if sysVer==""
            print("  loadModel("*sysLib*")")
            status = OMJulia.sendExpression(omc, "loadModel("*sysLib*")")
          else
            print("  loadModel("*sysLib*",{\""*sysVer*"\"})")
            status = OMJulia.sendExpression(omc, "loadModel("*sysLib*",{\""*sysVer*"\"})")
          end
          if status==true
            println(" successful")
          else
            println(" failed")
          end
        end
      else
        error("sysLibs and sysVers string entries, separated by : are not equal")
      end
    end

    # Load Modelica work files
    workFilesList = split(workFiles,":")
    for workFile in workFilesList
      if workFile != ""
        print("  loadFile("*workDir*"/"*workFile*")")
        OMJulia.sendExpression(omc, "loadFile(\""*workDir*"/"*workFile*"\")")
        if status==true
          println(" successful")
        else
          println(" failed")
        end
      end
    end

    # Load additional Modelica files
    filesList = split(files,":")
    for file in filesList
      if file!=""
        print("  loadFile("*file*")")
        status = OMJulia.sendExpression(omc, "loadfile(\"*file*\")")
        if status==true
        else
          println(" failed")
          println(" successful")
        end
      end
    end

    print("  instantiateModel("*model*")")
    status = OMJulia.sendExpression(omc, "instantiateModel("*model*")")
    if status!=""
      println(" successful")
    else
      println(" failed")
    end

    # Assemble options
    simulateOpts=",tolerance="*string(Tolerance)
    if StartTime!=0
      simulateOpts=", startTime="*string(StartTime)
    end
    if StopTime>StartTime
      simulateOpts=", stopTime="*string(StopTime)
    end
    if NumberOfIntervals!=500
      simulateOpts=", numberOfIntervals="*string(NumberOfIntervals)
    end
    # Run simulaton
    print("  simulate("*model*simulateOpts*")")
    status = OMJulia.sendExpression(omc, "simulate("*model*simulateOpts*")")
    if status!=""
      println(" successful")
    else
      println(" failed")
    end

    simFile = simDir*"/"*model*resExt
    if copy
      resFile = workDir*"/"*model*resExt
      println("  copy result file from: "*simDir*"/"*model*resExt)
      println("                     to: "*resFile)
      cp(simFile,resFile,force=true)
      # Return handle and result file
      return resFile
    else
      return simFile
    end
    cd(localDir)
  end

end


text:
basic(x,p::AbstractVector) = @. (1-(exp(-p[1]*x)))*p[2]
inverse(y,p::AbstractVector) = @. -(log(1-(y/p[2])))/p[1]
##
function basic(xaxis,case::String)
    println(case)
    if case == "low"
        p = [0.075,32.5]
    elseif case ==  "high"
        p = [0.075,57.5]
    else
        error("case not recognized")
    end
    basic(xaxis,p)
end

function inverse(xaxis,case::String)
    if case == :low
        p = [0.075,32.5]
    elseif case ==  :high
        p = [0.075,57.5]
    end
    inverse(xaxis,p)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 87f59dc7-5149-4eb6-9d81-440ee8cecd72
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using MLCourse, Plots, MLJ, DataFrames, Random, CSV, Flux, Distributions,
          StatsPlots, MLJFlux, OpenML, Random
    import MLCourse: heaviside
end

# ╔═╡ 83e2c454-042f-11ec-32f7-d9b38eeb7769
begin
    using PlutoUI
    PlutoUI.TableOfContents()
end

# ╔═╡ c95ad296-5191-4f72-a7d7-87ddebc43a65
md"# Solving the XOR Problem with Learned Feature Vectors

In the following example we define a loss function that learns the features (weights of the hidden neurons) and the \"regression coefficients\" (output weights) with gradient descent.
"

# ╔═╡ fde9639d-5f41-4037-ab7b-d3dbb09e8d3d
function xor_generator(; n = 200)
	x = 2 * rand(n, 2) .- 1
	DataFrame(X1 = x[:, 1], X2 = x[:, 2],
		      y = (x[:, 1] .> 0) .⊻ (x[:, 2] .> 0))
end

# ╔═╡ ea7fc5b6-79cc-4cc5-820e-5c55da83207e
begin
    xor_data = xor_generator()
    xor_input = Array(select(xor_data, Not(:y)))
    xor_target = xor_data.y
end;

# ╔═╡ a2539e77-943e-464c-9ad3-c8542eab36ab
begin
    f(x, w₁, w₂, w₃, w₄, β) = β[1] * relu.(x * w₁) .+
                              β[2] * relu.(x * w₂) .+
                              β[3] * relu.(x * w₃) .+
                              β[4] * relu.(x * w₄) .+
                              β[5]
    f(x, θ) = f(x, θ[1:2], θ[3:4], θ[5:6], θ[7:8], θ[9:13])
    function log_reg_loss(θ)
        p = σ.(f(xor_input, θ))
        -mean(@. xor_target * log(p) + (1 - xor_target) * log(1 - p))
    end
end

# ╔═╡ 258a3459-3fef-4655-830c-3bdf11eb282d
md"seed = $(@bind(seed, Select([string(i) for i in 1:20])))

t = $(@bind t Slider(1:10^2))
"

# ╔═╡ 49c2b3a1-50c0-4fb9-a06b-2de7be216702
begin
    Random.seed!(Meta.parse(seed))
    θ₀ = .1 * randn(13)
    path = [copy(θ₀)]
    MLCourse.gradient_descent(log_reg_loss, θ₀, .05, 10^4,
                              callback = x -> push!(path, copy(x)))
end;

# ╔═╡ 9aa58d2b-783a-4b74-ae36-f2ff0fb0be3f
md"The green arrows in the figure above are the weights w₁, …, w₄ of the four (hidden) neurons. The coloring of the dots is based on classification at decision threshold 0.5, i.e. if the activation of the output neuron is above 0.5 for a given input point, this point is classified as red (green, otherwise). Note that the initial configuration of the green vectors depends on the initialization of gradient descent and therefore on the seed. For some seeds the optimal solution to the xor-problem is not found by gradient descent."

# ╔═╡ b650ddb4-5854-4862-ba26-0bf92f77c3df
begin
    idxs = min.(max.(1:100, floor.(Int, 1.097.^(1:100))), 10^4)
    losses = log_reg_loss.(path[idxs])
end;

# ╔═╡ 16d0808e-4094-46ff-8f92-1ed21aa6191b
let idx = idxs[t]
    θ = path[idx]
    prediction = σ.(f(xor_input, θ)) .> .5
    p1 = scatter(xor_data.X1, xor_data.X2, c = prediction .+ 1,
                 xlim = (-1.5, 1.5), ylim = (-1.5, 1.5),
                 xlabel = "X1", ylabel = "X2")
    for i in 1:4
        plot!([0, θ[i*2 - 1]], [0, θ[i*2]], c = :green, w = 3, arrow = true)
    end
    hline!([0], c = :black)
    vline!([0], c = :black)
    p2 = plot(idxs, losses, ylim = (0, 1), xscale = :log10, xlabel = "gradient descent step t",
              ylabel = "loss", title = "learning curve")
    scatter!([idxs[t]], [losses[t]], c = :red)
    plot(p1, p2, layout = (1, 2), size = (700, 400), legend = false)
end

# ╔═╡ 8c54bfb2-54c5-4777-ad53-97926da97f7e
Markdown.parse("
## Avoiding Local Minima with Overparametrization
The following code fits 10 times with different initial guesses for the
parameters an MLP with 4 hidden neurons and 10 times an MLP with 10 hidden
neurons to our xor data. We see that training and test accuracy is better for
the larger network. The reason is that gradient descent for the the larger
network does not get stuck in suboptimal solutions.

Because it takes a bit of time to run the code, it is included here in text form.
Copy-paste it to a cell to run it.
```julia
begin
    using MLJFlux
    function fit_xor(n_hidden)
        x = select(xor_data, Not(:y))
        y = coerce(xor_data.y, Binary)
        builder = MLJFlux.Short(n_hidden = n_hidden,
                                dropout = 0,
                                σ = relu)
        mach = machine(NeuralNetworkClassifier(builder = builder,
                                               batch_size = 32,
                                               epochs = 10^4),
                       x, y) |> fit!
        xor_test = xor_generator(n = 10^4)
        x_test = select(xor_test, Not(:y))
        y_test = coerce(xor_test.y, Binary)
        DataFrame(n_hidden = n_hidden,
                  training_accuracy = mean(predict_mode(mach, x) .== y),
                  test_accuracy = mean(predict_mode(mach, x_test) .== y_test))
    end
    xor_results = vcat([fit_xor(n_hidden) for n_hidden in [4, 10], _ in 1:10]...)
    @df xor_results dotplot(:n_hidden, :training_accuracy,
                            xlabel = \"number of hidden units\",
                            ylabel = \"accuracy\",
                            label = \"training\",
                            yrange = (.5, 1.01),
                            aspect_ratio = 20,
                            legend = :bottomright,
                            xticks = [4, 10])
    @df xor_results dotplot!(:n_hidden, :test_accuracy, label = \"test\")
end
```
$(MLCourse.embed_figure("overparametrized.png"))
")

# ╔═╡ 5c2154ef-2a35-4d91-959d-f72100049894
md"# Multilayer Perceptrons

Below you see three equivalent ways to implement a neural network:
1. In the `expanded` form all parameters and activation functions are explicitely written out in \"one line\". With more neurons or hidden layers, this form quickly becomes annoying to write.
2. In the `neuron` form we see clearly that the neural network is composed of multiple simple neurons with different input weights and biases.
3. The `matrix` form is the most succinct form and most commonly used in modern implementations. In this form the biases of each layer are organized in a vector and the input weights of all neurons in the same layer are organized in a matrix.
"

# ╔═╡ 71dabb22-816b-4aa2-9a29-f2894989b2ea
begin
    function neural_net_expanded(x,
                                 w₁₀¹, w₁₁¹, w₂₀¹, w₂₁¹, w₃₀¹, w₃₁¹, w₄₀¹, w₄₁¹, g¹,
                                 w₁₀², w₁₁², w₁₂², w₁₃², w₁₄², g²)
        g²(w₁₀² +
           w₁₁² * g¹(w₁₀¹ + w₁₁¹ * x) +
           w₁₂² * g¹(w₂₀¹ + w₂₁¹ * x) +
           w₁₃² * g¹(w₃₀¹ + w₃₁¹ * x) +
           w₁₄² * g¹(w₄₀¹ + w₄₁¹ * x))
    end
    neuron(x, b, w, g) = g.(b .+ w' * x)
    function neural_net_neurons(x,
                                w₁₀¹, w₁₁¹, w₂₀¹, w₂₁¹, w₃₀¹, w₃₁¹, w₄₀¹, w₄₁¹, g¹,
                                w₁₀², w₁₁², w₁₂², w₁₃², w₁₄², g²)
        second_layer_input = [neuron(x, w₁₀¹, w₁₁¹, g¹),
                              neuron(x, w₂₀¹, w₂₁¹, g¹),
                              neuron(x, w₃₀¹, w₃₁¹, g¹),
                              neuron(x, w₄₀¹, w₄₁¹, g¹)]
        neuron(second_layer_input, w₁₀², [w₁₁², w₁₂², w₁₃², w₁₄²], g²)
    end
    function neural_net_matrix(x,
                               w₁₀¹, w₁₁¹, w₂₀¹, w₂₁¹, w₃₀¹, w₃₁¹, w₄₀¹, w₄₁¹, g¹,
                               w₁₀², w₁₁², w₁₂², w₁₃², w₁₄², g²)
        b¹ = [w₁₀¹, w₂₀¹, w₃₀¹, w₄₀¹]
        w¹ = [w₁₁¹
              w₂₁¹
              w₃₁¹
              w₄₁¹]
        b² = [w₁₀²]
        w² = [w₁₁² w₁₂² w₁₃² w₁₄²]
        neural_net_matrix(x, b¹, w¹, g¹, b², w², g²)
    end
    function neural_net_matrix(x, b¹, w¹, g¹, b², w², g²)
        g².(b² .+ w² * g¹.(b¹ .+ w¹ * x))
    end
end;

# ╔═╡ 653ea4e8-47b2-4d9f-bc5c-583468d6e1ae
md"## 1D to 1D
**First Layer** activation function ``g^{(1)}`` = $(@bind g1 Select(string.([relu, selu, sigmoid, tanh, identity, heaviside])))

``w_{10}^{(1)}`` = $(@bind w101 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{11}^{(1)}`` = $(@bind w111 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{20}^{(1)}`` = $(@bind w201 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{21}^{(1)}`` = $(@bind w211 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{30}^{(1)}`` = $(@bind w301 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{31}^{(1)}`` = $(@bind w311 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{40}^{(1)}`` = $(@bind w401 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{41}^{(1)}`` = $(@bind w411 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

**Second Layer** activation function ``g^{(2)}`` = $(@bind g2 Select(string.([identity, relu, selu, sigmoid, tanh, heaviside])))

``w_{10}^{(2)}`` = $(@bind w102 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{11}^{(2)}`` = $(@bind w112 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{12}^{(2)}`` = $(@bind w122 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{13}^{(2)}`` = $(@bind w132 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{14}^{(2)}`` = $(@bind w142 Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
"

# ╔═╡ 1171bcb8-a2af-49f3-9a7e-2e7cac34f9f4
let
    x = -5:.01:5
    _g1 = getproperty(@__MODULE__, Symbol(g1))
    _g2 = getproperty(@__MODULE__, Symbol(g2))
    y = neural_net_expanded.(x, w101, w111, w201, w211, w301, w311, w401, w411, _g1, w102, w112, w122, w132, w142, _g2)
    py = plot(x, y, xlabel = "x", ylabel = "y")
    p1 = plot(x, neuron.(x, w101, w111, _g1), xlabel = "x", ylabel = "hidden neuron 1")
    p2 = plot(x, neuron.(x, w201, w211, _g1), xlabel = "x", ylabel = "hidden neuron 2")
    p3 = plot(x, neuron.(x, w301, w311, _g1), xlabel = "x", ylabel = "hidden neuron 3")
    p4 = plot(x, neuron.(x, w401, w411, _g1), xlabel = "x", ylabel = "hidden neuron 4")
    plot(plot(p1, p2, p3, p4, layout = (1, 4)), py, layout = (2, 1), legend = false)
end

# ╔═╡ 1c70c3c4-bded-44ec-ac40-c71268b7307e
md"## 2D to 1D
**First Layer** activation function ``g^{(1)}`` = $(@bind g1_2d Select(string.([relu, selu, sigmoid, tanh, identity, heaviside]), default = string(:heaviside)))

``w_{10}^{(1)}`` = $(@bind w101_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{11}^{(1)}`` = $(@bind w111_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{12}^{(1)}`` = $(@bind w121_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{20}^{(1)}`` = $(@bind w201_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{21}^{(1)}`` = $(@bind w211_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{22}^{(1)}`` = $(@bind w221_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{30}^{(1)}`` = $(@bind w301_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{31}^{(1)}`` = $(@bind w311_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{32}^{(1)}`` = $(@bind w321_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{40}^{(1)}`` = $(@bind w401_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{41}^{(1)}`` = $(@bind w411_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{42}^{(1)}`` = $(@bind w421_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

**Second Layer** activation function ``g^{(2)}`` = $(@bind g2_2d Select(string.([identity, relu, selu, sigmoid, tanh, heaviside])))

``w_{10}^{(2)}`` = $(@bind w102_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{11}^{(2)}`` = $(@bind w112_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{12}^{(2)}`` = $(@bind w122_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))

``w_{13}^{(2)}`` = $(@bind w132_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
``w_{14}^{(2)}`` = $(@bind w142_2d Slider(-3.:.1:3., default = floor(60*(rand()-.5))/10, show_value = true))
"

# ╔═╡ 15cd82d3-c255-4eb5-9f0d-1e662c488027
md"If you want to tinker more with simple neural networks there is [this nice web-application](http://playground.tensorflow.org/)."

# ╔═╡ 9e6ad99f-8cd4-4372-a4a2-f280e49e21a3
begin
    function tocol(y)
        ma, mi = maximum(y), minimum(y)
        palette(:jet1, 101)[min.(101, max.(1, reshape(floor.(Int, (y .+ 6) * 5) .+ 1, :)))]
    end
end;

# ╔═╡ 66bc40f3-7710-4892-ae9f-79b9c850b0ea
let x = -5:.1:5
    grid = MLCourse.grid(x, x, output_format = DataFrame)
    X = Array(Array(grid)')
    _g1 = getproperty(@__MODULE__, Symbol(g1_2d))
    _g2 = getproperty(@__MODULE__, Symbol(g2_2d))
    b1 = [w101_2d, w201_2d, w301_2d, w401_2d]
    w1 = [w111_2d w121_2d
          w211_2d w221_2d
          w311_2d w321_2d
          w411_2d w421_2d]
    b2 = [w102_2d]
    w2 = [w112_2d w122_2d w132_2d w142_2d]
    y = neural_net_matrix(X, b1, w1, _g1, b2, w2, _g2)
    p1 = scatter(grid.X, grid.Y, c = tocol(neuron(X, b1[1], w1[1, :], _g1)), title = "hidden 1")
    p2 = scatter(grid.X, grid.Y, c = tocol(neuron(X, b1[2], w1[2, :], _g1)), title = "hidden 2")
    p3 = scatter(grid.X, grid.Y, c = tocol(neuron(X, b1[3], w1[3, :], _g1)), title = "hidden 3")
    p4 = scatter(grid.X, grid.Y, c = tocol(neuron(X, b1[4], w1[4, :], _g1)), title = "hidden 4")
    py = scatter(grid.X, grid.Y, c = tocol(y), title = "y",)
    plot(plot(p1, p2, p3, p4, layout = (1, 4)),
         py, layout = (2, 1), markersize = 2, markerstrokewidth = 0,
                 xlabel = "x₁", ylabel = "x₂", legend = false, aspect_ratio = 1)
end

# ╔═╡ 9d3e7643-34fd-40ee-b442-9d2a434f30e0
Markdown.parse("# Regression with MLPs

## Fitting the Weather Data

Let us find out if with an MLP we can get better predictions than with multiple
linear regression on the weather data set.

The following code takes roughly 2 minutes to run.

```julia
weather = CSV.read(joinpath(@__DIR__, \"..\", \"data\", \"weather2015-2018.csv\"),
                   DataFrame)
weather_test = CSV.read(joinpath(@__DIR__, \"..\", \"data\", \"weather2019-2020.csv\"),
                        DataFrame)
weather_x = float.(select(weather, Not([:LUZ_wind_peak, :time]))[1:end-5, :])
weather_y = weather.LUZ_wind_peak[6:end]
weather_test_x = float.(select(weather_test, Not([:LUZ_wind_peak, :time]))[1:end-5, :])
weather_test_y = weather_test.LUZ_wind_peak[6:end]

Random.seed!(31)
nn = NeuralNetworkRegressor(builder = MLJFlux.Short(n_hidden = 128,
                                                    dropout = .5,
                                                    σ = relu),
                            optimiser = ADAMW(),
                            batch_size = 128,
                            epochs = 150)
mach = machine(Standardizer() |> nn,
               weather_x, weather_y)
fit!(mach, verbosity = 2)
```

!!! note

    We use here a `Pipeline` (`|>`) that standardizes with a `Standardizer`
    the input first to mean zero and standard deviation one.
    We do this, because the usual initializations
    of neural networks work best with standardized input.

    To learn more about the usual initializations you can have a look at this [blog post](https://towardsdatascience.com/weight-initialization-in-neural-networks-a-journey-from-the-basics-to-kaiming-954fb9b47c79) (optional).

!!! note

    Our builder constructs a neural network with one hidden layer of 128 relu neurons.
    In the construction we specify that we would like to train with 50% dropout
    to reduce the risk of overfitting.

Let us evalute the learned machine.

```julia
(training_error = rmse(predict(mach, weather_x), weather_y),
 test_error = rmse(predict(mach, weather_test_x), weather_test_y))
```
We find a training error of 7.38 and a test error of 8.65.
Both errors are lower than what we found with multiple linear regression
(training error ≈ 8.09, test_error ≈ 8.91). However, comparing the predictions
to the ground truth we see that a lot of the variability is still uncaptured
by the model. This means most likely that the irreducible error is pretty high in this dataset.
```julia
scatter(predict(mach, weather_test_x), weather_test_y,
        xlabel = \"predicted\", ylabel = \"data\")
plot!(identity, legend = false)
```
$(MLCourse.embed_figure("weather_mlp.png"))
")

# ╔═╡ 0ba71a2a-f2f9-4fc7-aa81-416799e79e57
Markdown.parse("## Parametrizing Variances of Log-Normal Distributions with MLPs

For this example we do not use the abstractions provided by MLJ. Instead we use
directly Flux. You only need to do this if you want to implement some non-standard neural network.

Let us load the data.
```julia
weather = CSV.read(joinpath(@__DIR__, \"..\", \"data\", \"weather2015-2018.csv\"),
                   DataFrame)
standardize(x) = (x .- mean(x)) ./ std(x)
weather_input = Float32.(standardize(weather.LUZ_pressure[1:end-5])')
weather_output = Float32.(standardize(log.(weather.LUZ_wind_peak[6:end])))
```

We define a simple neural network with one hidden layer of 50 tanh neurons.
```julia
nn = Chain(Dense(1, 50, tanh), Dense(50, 2))
```

Let us now define the negative log-likelihood loss function.
```julia
function loss(x, y)
    output = nn(x)
    m = output[1, :]
    s = softplus.(output[2, :])
    mean((m .- y) .^ 2 ./ (2 * s) .+ log.(s))
end
```
This function computes the two output dimensions of the neural network.
The first output dimension is our mean. The second output dimension is transformed
to positive values with the `softplus` function ``x \\mapsto \\log(\\exp(x) + 1)``
to enode the variance of the normal distribution. In the log-likelihood loss we
drop the constant term ``\\frac12\\log2\\pi``.

Now we train the neural network for 50 epochs.
```julia
opt = ADAMW()
p = Flux.params(nn) # these are the parameters to be adapted.
data = Flux.DataLoader((weather_input, weather_output), batchsize = 32)
for _ in 1:50
    Flux.Optimise.train!(loss, p, data, opt)
end
```

Let us now plot the result.
```julia
grid = collect(-5:.1:3.5)'
pred = nn(grid)
m = pred[1, :]
s = sqrt.(softplus.(pred[2, :]))
histogram2d(weather_input', weather_output, bins = (250, 200),
            markersize = 3, xlabel = \"LUZ_pressure [hPa]\",
            ylabel = \"LUZ_wind_peak [km/h]\", label = nothing, cbar = false)
plot!(grid', m, c = :red, w = 3, label = \"mean\")
plot!(grid', m .+ s, c = :red, linestyle = :dash, w = 3, label = \"± standard deviation\")
plot!(grid', m .- s, c = :red, linestyle = :dash, w = 3, label = nothing)
```

$(MLCourse.embed_figure("weather_mlp_normal.png"))

We see that the neural network found a slightly non-linear relationship between
input and average output and the variance of the noise is smaller for high
pressure values than for low pressure values.
")

# ╔═╡ 4bd63efb-0aef-41cd-a9b7-9fa78db107a0
Markdown.parse("
## Fitting a Polynomial in 1000 Input Dimensions
The following artificial data set has ``n=10^5`` points of ``p = 1000`` dimensions.
The response ``y = x_7^2 - 1`` depends only on the seventh input dimension.
Although neural networks are not specifically designed to fit polynomials,
this example demonstrates that with only 50 hidden relu-neurons an accurate fit
can be found.

WARNING: it takes more than 10 minutes to run this code.
```julia
let
    x = DataFrame(randn(10^5, 1000), :auto)
    y = x.x7 .^2 .- 1
    mach = fit!(machine(NeuralNetworkRegressor(builder = MLJFlux.Short(n_hidden = 50,
                                                                       dropout = 0,
                                                                       σ = relu),
                                               batch_size = 128,
                                               epochs = 10^2),
                        x, y), verbosity = 2)
    grid = -3:.1:3
    x_test = DataFrame(randn(length(grid), 1000), :auto)
    x_test.x7 = grid
    plot(grid, grid.^2 .- 1, label = \"ground truth\", w = 3, legend = :bottomright)
    plot!(grid, predict(mach, x_test),
          label = \"fit\", xlabel = \"x₇\", ylabel = \"y\", w = 3)
end
```

$(MLCourse.embed_figure("poly.png"))
")

# ╔═╡ 14a8eacb-352e-4c3b-8908-2a6f77ffc6fe
md"# Classification with MLPs

## Fitting MNIST

The following code to fit the MNIST data with a multilayer perceptron takes a
few minutes to run.

```julia
mnist_x, mnist_y = let df = OpenML.load(554) |> DataFrame
    coerce!(df, :class => Multiclass)
    coerce!(df, Count => MLJ.Continuous)
    Float32.(df[:, 1:end-1] ./ 255),
    df.class
end
mnist_mach = machine(NeuralNetworkClassifier(
                         builder = MLJFlux.@builder(Chain(Dense(n_in, 100, relu),
                                                          Dense(100, n_out))),
                         batch_size = 32,
                         epochs = 20),
                      mnist_x[1:60000, :],
                      mnist_y[1:60000])
fit!(mnist_mach, verbosity = 2)
mean(predict_mode(mnist_mach, mnist_x[60001:70000, :]) .!= mnist_y[60001:70000])
```
The result is a misclassification rate of approximately 2.2%, which is better
than our first nearest-neigbor result. Additionally, prediction of the class
for new images is much faster than with nearest neighbor approaches, because the
information about the training data is stored in the weights of the neural network
and no explicit comparison with all training images is required.
"

# ╔═╡ 2f034d19-1c00-4a10-a880-99436ab00957
md"# Exercises

## Conceptual
1. To get a feeling for the kind of functions of one predictor can be fitted with neural networks, we will draw ``y`` as a function of ``x`` for some values of the weights. It may be helpful to sketch this neural network with the input neuron, the hidden neurons and the output neuron and label the connections with the weights.
    * Draw in the same figure ``a_1^{(1)} = g(w_{10}^{(1)} + w_{11}^{(1)} x)``, ``a_2^{(1)}=g(w_{20}^{(1)} + w_{21}^{(1)} x)`` and ``\bar y = w_0^{(2)} + w_1^{(2)}a_1^{(1)} +  w_2^{(2)}a_2^{(1)}`` as a function of ``x``. Use ``w_{10}^{(1)} = 0``, ``w_{11}^{(1)} = 1``, ``w_{20}^{(1)} = - 2``, ``w_{21}^{(1)} = 2``,  ``w_0^{(2)} = 1``, ``w_1^{(2)} = 2``, ``w_2^{(2)} = 1`` and use the rectified linear activation function ``g = \mbox{relu}``. At which ``x``-values does the slope change? Give the answer in terms of the weights.
    * Draw a similar graph for ``w_{11}^{(1)} < 0`` and ``w_{10}^{(1)}/w_{11}^{(1)}<w_{20}^{(1)}/w_{21}^{(1)}``, e.g. with ``w_{10}^{(1)} = 0``, ``w_{11}^{(1)} = -1``, ``w_{20}^{(1)} = 2``, ``w_{21}^{(1)} = 2``,  ``w_0^{(2)} = 1``, ``w_1^{(2)} = 2``, ``w_2^{(2)} = 1``.
    * How does the graph above change, if ``w_1^{(2)} = -2`` and ``w_2^{(2)} = -1``?
    * Let us assume we add more neurons to the same hidden layer, i.e. we have ``a_1^{(1)}, \ldots, a_{d^{(1)}}^{(1)}`` activations in the first layer. How would the graph look differently in this case? Draw a sketch and describe in one sentence the qualitative difference.
    * Let us assume we add instead more hidden layers with relu-activations. Would the graph of this neural network look qualitatively different from the ones we have drawn so far in this exercise?
    * (optional) Show that a neural network with one hidden layer of 3 relu-neurons can perfectly fit any continuous piece-wise linear function of the form ``y = \left\{\begin{array}{ll} a_1 + b_1 x & x < c_1 \\ a_2 + b_2 x & c_1 \leq x < c_2 \\ a_3 + b_3 x & c_2 \leq x \end{array}\right. \hspace{1cm}  ``
                    with ``c_1 =  \frac{a_1 - a_2}{b_2 - b_1} < c_2 = \frac{a_2 - a_3}{b_3 - b_2}``. Express ``a_1, a_2, a_3`` and ``b_1, b_2, b_3`` in terms of the network weights. There are multiple solutions; find one of them.
1. Consider a neural network with two hidden layers: ``p = 4`` input units, ``2`` units in the first hidden layer, ``3`` units in the second hidden layer, and a single output.
    - Draw a picture of the network.
    - How many parameters are there?
    - Assume the output of this network is the mean of a conditional normal distribution. Write the negative log-likelihood loss using matrix notation.
1. Let us assume you want to predict the proportion ``y`` of citizens voting for one of two parties based on 113 different features that depend e.g. on the voting history, income distribution or demographics (see e.g. [here](https://www.ozy.com/news-and-politics/the-forecast-the-methodology-behind-our-2020-election-model/379778/) for an actual example). Because the proportion is a number between 0 and 1, it is reasonable to model it as a sample from a [beta distribution](https://en.wikipedia.org/wiki/Beta_distribution) ``p(y) = \frac{y^{\alpha - 1}(1-y)^{\beta - 1}}{B(\alpha, \beta)}``, where ``\alpha>0`` and ``\beta>0`` are shape parameters and ``B`` is the beta function. You want parametrize the conditional distribution of proportions ``y`` given predictors ``x`` with a neural network.
    * How many input units should this neural network have?
    * How many output units should this neural network have?
    * What kind of activation function would you choose for the output layer?
    * Which loss function would you use to estimate the network parameters?
1. (Optional) Grant Sanderson has some beautiful [videos about neural networks](https://www.3blue1brown.com/topics/neural-networks). Have a look at them, if you are interested.

## Applied
1. In this exercise our goal is to find a good machine learning model to predict the fat content of a meat sample on the basis of its near infrared absorbance spectrum. We use the Tecator data set `OpenML.describe_dataset(505)`. The first 100 columns of this data set contain measurements of near infrared absorbance at different frequencies for different pieces of meat. *Hint:* you can select all these columns based on name with `select(data, r\"absorbance\")` (The \"r\" in this command stands for Regex and it means that all columns with name containing the word \"absorbance\" should be selected). The column `:fat` contains the fat content of each piece of meat. We will use the first 172 data points for training and validation and the rest of the data points as a test set. Take our recipe for supervised learning (last slide of the presentation on \"Model Assessment and Hyperparameter Tuning\") as a guideline.
    - Have a look at the raw data by e.g. checking if there are missing values and looking at a correlation plot.
    - Fit some multiple linear regression models (with e.g. regularization constants tuned with cross-valdiation), compute the `rmse` of your best model on the test set and create a scatter plot that shows the actual fat content of the test data point versus the predicted fat content.
    - Fit some neural network model (with 2 hyper-paramters of your choice tuned by cross-validation; warning: it may take quite some time to fit if you use a high number for `nfolds`), compute the `rmse` of your best model on the test set and create a scatter plot that shows the actual fat content of the test data point versus the predicted fat content. *Hint 1:* standardization of input and output matters. *Hint 2:* If you want to tune some neural network parameters in a pipeline you can usually access them in the `range` function as `:(neural_network_regressor.builder.dropout)`, for example.
"

# ╔═╡ 7c6951e8-726a-4181-b15f-b629a2835d03
MLCourse.list_notebooks(@__FILE__)

# ╔═╡ 8c72c4b5-452d-4ba3-903b-866cac1c799d
MLCourse.footer()

# ╔═╡ Cell order:
# ╠═87f59dc7-5149-4eb6-9d81-440ee8cecd72
# ╟─c95ad296-5191-4f72-a7d7-87ddebc43a65
# ╠═fde9639d-5f41-4037-ab7b-d3dbb09e8d3d
# ╠═ea7fc5b6-79cc-4cc5-820e-5c55da83207e
# ╠═a2539e77-943e-464c-9ad3-c8542eab36ab
# ╠═49c2b3a1-50c0-4fb9-a06b-2de7be216702
# ╟─258a3459-3fef-4655-830c-3bdf11eb282d
# ╟─16d0808e-4094-46ff-8f92-1ed21aa6191b
# ╟─9aa58d2b-783a-4b74-ae36-f2ff0fb0be3f
# ╟─b650ddb4-5854-4862-ba26-0bf92f77c3df
# ╟─8c54bfb2-54c5-4777-ad53-97926da97f7e
# ╟─5c2154ef-2a35-4d91-959d-f72100049894
# ╠═71dabb22-816b-4aa2-9a29-f2894989b2ea
# ╟─653ea4e8-47b2-4d9f-bc5c-583468d6e1ae
# ╟─1171bcb8-a2af-49f3-9a7e-2e7cac34f9f4
# ╟─1c70c3c4-bded-44ec-ac40-c71268b7307e
# ╟─66bc40f3-7710-4892-ae9f-79b9c850b0ea
# ╟─15cd82d3-c255-4eb5-9f0d-1e662c488027
# ╟─9e6ad99f-8cd4-4372-a4a2-f280e49e21a3
# ╟─9d3e7643-34fd-40ee-b442-9d2a434f30e0
# ╟─0ba71a2a-f2f9-4fc7-aa81-416799e79e57
# ╟─4bd63efb-0aef-41cd-a9b7-9fa78db107a0
# ╟─14a8eacb-352e-4c3b-8908-2a6f77ffc6fe
# ╟─2f034d19-1c00-4a10-a880-99436ab00957
# ╟─7c6951e8-726a-4181-b15f-b629a2835d03
# ╟─83e2c454-042f-11ec-32f7-d9b38eeb7769
# ╟─8c72c4b5-452d-4ba3-903b-866cac1c799d

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
export ppsimple

function ppsimple(pold::PolynomialProbability)
    pnew = PolynomialProbability()

    for pi = 1:length(pold)
        piece = pold[pi]
        if (!isempty(piece.polys))
            order = sortperm(piece.exp, rev=true)
            allexp = piece.exp[order]
            allpoly = piece.polys[order, :]

            # Combine like exponents
            npoly = allpoly[1, :]
            nexp = [allexp[1]]
            for ii = 2:length(allexp)
                if allexp[ii] == nexp[end]
                    npoly[end, :] = npoly[end, :] + allpoly[ii, :]
                else
                    push!(nexp, allexp[ii])
                    npoly = vcat(npoly, allpoly[ii, :])
                end
            end

            # Eliminate zero rows and columns
            nonzero = npoly .== 0 # just initialize
            nonzero[:, end:-1:1] = cumsum(float(npoly[:, end:-1:1] .!= 0), 2) .!= 0
            colskept = sum(nonzero, 1) .!= 0
            rowskept = sum(nonzero, 2) .!= 0
            npoly = npoly[rowskept, colskept]
            nexp = nexp[collect(rowskept)]
        else
            npoly = zeros(1, 1)
            nexp = [0.]
        end

        append!(pnew, piece.t0, piece.t1, npoly, nexp)
    end

    pnew
end


text:
"""
function lineup
"""
function lineup(event_data::PublicWyscoutLoader, game_id::Int)
    tmp = filter("match_id" => ==(game_id), event_data.match_index)
    data = get_matchs(event_data, tmp[1,:db_matches])
    data = filter("wyId"=> ==(game_id), data)[:,"teamsData"]
    data = DataFrame(data)
    data = Dict(i => DataFrame(data[:,i]) for i in names(data))

    return data
end

"""
    function get_df()
"""
function get_df(path::String)

    data = DataFrame(jsontable(JSON3.read(read(path))))

    return data
end


text:
randc(m,n) = randn(m,n)+im*randn(m,n)

function grf(m,n,C::StationaryCovarianceFunction)
    tx = collect(0:m-1)
    ty = collect(0:n-1)
    Rjs = zeros(n,m)
    for i in 1:m, j in 1:n
        Rjs[j,i] = C([tx[1],ty[1]],[tx[i],ty[j]])
    end
    Sjs = [Rjs Rjs[:,end] Rjs[:,end:-1:2]]
    Sjs2 = [Sjs ; Sjs[end,:] ; Sjs[end:-1:2,1] Sjs[end:-1:2,end:-1:2]]
    Γ = fft(Sjs2)
    Z = randc(size(Sjs2)...)
    fft(sqrt(Γ./(4*m*n)).*Z)[1:m,1:n]
end


text:
using Documenter, MdSlicing

makedocs(sitename="My Documentation")


text:
using Pkg

Pkg.develop([
    PackageSpec(
        name = "FoldsChainRules",
        path = dirname(@__DIR__),
        # url = ...,
    ),
    PackageSpec(
        name = "FoldsChainRulesBenchmarks",
        path = dirname(@__DIR__),
        # url = ...,
        subdir = "benchmark/FoldsChainRulesBenchmarks",
    ),
])

Pkg.instantiate()