text:
#####################
# auxiliary functions
#####################

"""
Calculate new value of pixel using mask.

#Arguments
- `copied::Array{Float64,2}`: copy of an color matrice.
- `mask::Array{Float64,2}`: mask used in transition.
- `x::Int64`: first coordinate of left-upper corner of copied.
- `y::Int64`: second coordinate of left-upper corner of copied.
"""
function layingMask(copied, mask, x, y)
    dimM = size(mask)
    return sum(copied[x:x + dimM[2] - 1, y:y + dimM[1] - 1] .* mask) / sum(mask)
end

"""
Assign to the pixel minimum value from surrounding ones, determined by mask.

#Arguments
- `copied::Array{Float64,2}`: copy of an color matrice.
- `mask::Array{Float64,2}`: mask used in transition.
- `x::Int64`: first coordinate of left-upper corner of copied.
- `y::Int64`: second coordinate of left-upper corner of copied.
"""
function minimumValue(copied, mask, x, y)
    dimM = size(mask)
    return minimum(copied[x:x + dimM[2] - 1, y:y + dimM[1] - 1])
end

"""
Assign to the pixel maximum value from the surrounding ones, determined by mask.

#Arguments
- `copied::Array{Float64,2}`: copy of an color matrice .
- `mask::Array{Float64,2}`: mask used in transition.
- `x::Int64`: first coordinate of left-upper corner of copied.
- `y::Int64`: second coordinate of left-upper corner of copied.
"""
function maximumValue(copied, mask, x, y)
    dimM = size(mask)
    return maximum(copied[x:x + dimM[2] - 1, y:y + dimM[1] - 1])
end

"""
Convert given picture with specified type of mask.
Copy every RGB matrix and reproduce extreme pixels.
Create new matrix with values calculated by layingMask function and return list of new R, G, B matrices.

#Arguments
- `picture::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `typeofmask::Array{Float64,2}`: mask, determined by typeofmask function.
"""
function converting(picture, typeofmask)
    mask = typeofmask
    dimR = size(picture[1])
    dimM = size(mask)
    R = copy(picture[1])
    G = copy(picture[2])
    B = copy(picture[3])
    newColorList = Array([])

    @simd for matrix in [R, G, B]
        copied = matrix

        helpHorizontalTop = copied[1,:]    # first row of matrix
        for n in 1:((dimM[1] - 1) / 2)
            copied = [helpHorizontalTop';copied] # adding upper bufor ((dimM[1]-1)/2) times
        end

        helpHorizontalBottom = copied[end,:]    # last row of matrix
        for n in 1:((dimM[1] - 1) / 2)
            copied = [copied;helpHorizontalBottom'] # adding bottom bufor ((dimM[1]-1)/2) times
        end

        helpVerticalLeft = copied[:,1]  # first column of matrix
        for i in 1:((dimM[2] - 1) / 2)
            copied = [helpVerticalLeft  copied]  # adding left bufor ((dimM[1]-1)/2) times
        end

        helpVerticalRight = copied[:,end]  # last column of matrix
        for i in 1:((dimM[2] - 1) / 2)
            copied = [copied  helpVerticalRight]  # adding right bufor ((dimM[2]-1)/2) times
        end

        dimC = size(copied)
        newPixelValues = Array{Float64}(undef, dimR[1], dimR[2])  # reserving memory
        matrix = collect(layingMask(copied, mask, x, y) for x in 1:Int64(dimC[1] - dimM[1] + 1), y in 1:Int64(dimC[2] - dimM[2] + 1))

        push!(newColorList, matrix)
    end # loop for R G B matrices
    return newColorList
end # converting

######################
# Low-pass filters
######################

"""
Generate mask to blurr the picture.

#Arguments
- `dim::Int64`: length of side of a mask.
"""
function maskAverage(dim)
    mask = ones(Float64, dim, dim)
    return mask
end

"""
Generate mask to blurr the picture.

#Arguments
- `dim::Int64`: length of side of a mask.
"""
function circle(dim)
    mask = ones(Float64, dim, dim)
    mask[1,1] = mask[1,dim] = mask[dim,1] = mask[dim,dim] = 0
    return mask
end

"""
Generate mask to blurr the picture.

#Arguments
- `dim::Int64`: length of side of a mask.
- `param k::Float64`: value of central pixel.
"""
function LP3(dim, k)
    mask = ones(Float64, dim, dim)
    mask[Int64((dim - 1) / 2 + 1),Int64((dim - 1) / 2 + 1)] = k
    return mask
end

######################
# High-pass filters
######################

"""
Generate mask to sharpen a picture.

#Arguments
- `dim::Int64`: length of side of a mask.
- `param k:Float64`: value of central pixel.
"""
function meanRemoval(dim, k)
    mask = -1(ones(Float64, dim, dim))
    mask[Int64((dim - 1) / 2 + 1),Int64((dim - 1) / 2 + 1)] = k
    return mask
end

"""
Generate mask to sharpen a picture.

#Arguments
- `dim::Int64`: length of side of a mask.
- `param k:Float64`: value of central pixel.
"""
function HP1(dim, k)
    mask = -1(ones(Float64, dim, dim))
    mask[1,1] = mask[1,dim] = mask[dim,1] = mask[dim,dim] = 0
    mask[Int64((dim - 1) / 2 + 1),Int64((dim - 1) / 2 + 1)] = k
    return mask
end

"""
Generate mask to sharpen a picture.

#Arguments
- `dim::Int64`: length of side of a mask.
- `param k:Float64`: value of central pixel.
"""
function HP2(dim, k)
    mask = -2(ones(Float64, dim, dim))
    mask[1,1] = mask[1,dim] = mask[dim,1] = mask[dim,dim] = 1
    mask[Int64((dim - 1) / 2 + 1),Int64((dim - 1) / 2 + 1)] = k
    return mask
end

#####
# Edge detection
#####
"""
Generate mask to detect edges on the picture in given direction.

#Arguments
- `dim::Int64`: length of side of a mask.
- `direction::String`: abbreviation of the name of direction:
                     "h" horizontal edges,
                     "v" vertical edges,
                     "d" diagonal edges.
"""
function detectingEdges(dim, direction)
    mask = (zeros(Float64, dim, dim))
    mask[Int64((dim - 1) / 2 + 1),Int64((dim - 1) / 2 + 1)] = 1

    if direction == "h"
        mask[1,Int64((dim - 1) / 2 + 1)] = -1
    elseif direction == "v"
       mask[Int64((dim - 1) / 2 + 1),1] = -1
    elseif direction == "d"
       mask[1,1] = -1
    end

    return mask
end

"""
Convert a given picture with a mask. Copy every RGB matrix and reproduce extreme pixels.
Create new matrix filled with values calculated by minimumValue function
and return list of of new R G and B values.

#Arguments
- `picture::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `dimension::Int64`: length of side of a mask.
"""
function mini(picture, dimension)
    dimR = size(picture[1])
    dimM = dimension
    mask = zeros(dimM, dimM)
    R = copy(picture[1])
    G = copy(picture[2])
    B = copy(picture[3])
    newColorList = Array([])

    @simd for matrix in [R, G, B]
        copied = matrix

        helpHorizontalTop = copied[1,:]    # first row of matrix
        for n in 1:((dimM - 1) / 2)
            copied = [helpHorizontalTop';copied] # adding upper bufor ((dimM[1]-1)/2) times
        end

        helpHorizontalBottom = copied[end,:]    # last row of matrix
        for n in 1:((dimM - 1) / 2)
            copied = [copied;helpHorizontalBottom'] # adding bottom bufor ((dimM[1]-1)/2) times
        end

        helpVerticalLeft = copied[:,1]  # first column of matrix
        for i in 1:((dimM - 1) / 2)
            copied = [helpVerticalLeft  copied]  # adding left bufor ((dimM[1]-1)/2) times
        end

        helpVerticalRight = copied[:,end]  # last column of matrix
        for i in 1:((dimM - 1) / 2)
            copied = [copied  helpVerticalRight]  # adding right bufor ((dimM[2]-1)/2) times
        end

        dimC = size(copied)
        newPixelValues = Array{Float64}(undef, dimR[1], dimR[2])  # reserving memory
        matrix = collect(minimumValue(copied, mask, x, y) for x in 1:Int64(dimC[1] - dimM + 1), y in 1:Int64(dimC[2] - dimM + 1))

        push!(newColorList, matrix)
    end # loop for R G B matrices
    return newColorList
end # mini

"""
Convert a given picture with mask. Copy every RGB matrix and reproduce extreme pixels.
Then create new matrix with values calculated by maximumValue function
and return list of new R G and B values.

#Arguments
- `picture::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `dimension::Int64`: length of a side of the mask.
"""
function maxi(picture, dimension)
    dimR = size(picture[1])
    dimM = dimension
    mask = zeros(dimM, dimM)
    R = copy(picture[1])
    G = copy(picture[2])
    B = copy(picture[3])
    newColorList = Array([])

    @simd for matrix in [R, G, B]
        copied = matrix

        helpHorizontalTop = copied[1,:]    # first row of matrix
        for n in 1:((dimM - 1) / 2)
            copied = [helpHorizontalTop';copied] # adding upper bufor ((dimM[1]-1)/2) times
        end

        helpHorizontalBottom = copied[end,:]    # last row of matrix
        for n in 1:((dimM - 1) / 2)
            copied = [copied;helpHorizontalBottom'] # adding bottom bufor ((dimM[1]-1)/2) times
        end

        helpVerticalLeft = copied[:,1]  # first column of matrix
        for i in 1:((dimM - 1) / 2)
            copied = [helpVerticalLeft  copied]  # adding left bufor ((dimM[1]-1)/2) times
        end

        helpVerticalRight = copied[:,end]  # last column of matrix
        for i in 1:((dimM - 1) / 2)
            copied = [copied  helpVerticalRight] # adding right bufor ((dimM[2]-1)/2) times
        end

        dimC = size(copied)
        newPixelValues = Array{Float64}(undef, dimR[1], dimR[2])  # reserving memory
        matrix = collect(maximumValue(copied, mask, x, y) for x in 1:Int64(dimC[1] - dimM + 1), y in 1:Int64(dimC[2] - dimM + 1))

        push!(newColorList, matrix)
    end # loop for R G B matrices
    return newColorList
end # maxi

########
# WHATIF
########
"""
Convert blue matrice of a given picture with mask. Copy every of RGB matrices
and reproduce extreme pixels of B matrix.
Then create new B matrix with values calculated by layingMask function
and return list of R G and new B values.

#Arguments
- `picture::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `typeofmask::Array{Float64,2}`: mask which is a result of typeofmask function.
"""
function whatIfJustBlue(picture, typeofmask)
    mask = typeofmask
    dimR = size(picture[1])
    dimM = size(mask)
    R = copy(picture[1])
    G = copy(picture[2])
    B = copy(picture[3])
    newColorList = [R,G]

    copied = B

    helpHorizontalTop = copied[1,:]
    for n in 1:((dimM[1] - 1) / 2)
        copied = [helpHorizontalTop';copied]
    end

    helpHorizontalBottom = copied[end,:]
    for n in 1:((dimM[1] - 1) / 2)
        copied = [copied;helpHorizontalBottom']
    end

    helpVerticalLeft = copied[:,1]
    for i in 1:((dimM[2] - 1) / 2)
        copied = [helpVerticalLeft  copied]
    end

    helpVerticalRight = copied[:,end]
    for i in 1:((dimM[2] - 1) / 2)
        copied = [copied  helpVerticalRight]
    end

    dimC = size(copied)
    newPixelValues = Array{Float64}(undef, dimR[1], dimR[2]) # reserving memory
    B = collect(layingMask(copied, mask, x, y) for x in 1:Int64(dimC[1] - dimM[1] + 1), y in 1:Int64(dimC[2] - dimM[2] + 1))
    push!(newColorList, B)

    return newColorList
end # whatIfJustBlue

max_stars_repo_name: o-Mateo-o/ImageEditorApp

max_issues_repo_name: o-Mateo-o/ImageEditorApp

max_forks_repo_name: o-Mateo-o/ImageEditorApp


text:
"""
biconnected component decomposition
----------------------------------

Any connected graph decomposes into a tree of biconnected components
called the block-cut tree of the graph. The blocks are attached to each
other at shared vertices called cut vertices or articulation points.

This implementation is based on the algorithm provided by Tarjan
in "Depth-First Search and Linear Graph Algorithms".
"""
mutable struct Biconnected_components_output
    map::Vector{Int64} #biconnected_component_number
    articulation_points::Vector{Bool}
    number::Int64
    A::MatrixNetwork #MatrixNetwork
end

"""
`biconnected_components!`
---
This function returns the number of biconnected components in the
underlying graph. It expects an undirected graph as its input.

Functions
---------
- `number = biconnected_components!(A::MatrixNetwork, articulation::Vector{Bool}, map::Vector{Int64})`

Inputs
------
- `A`: the adjacency matrix.
- `articulation`: A boolean array, where each element is initialized to false.
- `map`: Vector of size equal to the number of edges.

Returns
-------
- `cn`: The number of biconnected components in the graph

Example
-------
A = load_matrix_network("biconnected_example")
B = MatrixNetwork(A)
number_of_components = biconnected_components!(B, zeros(Bool,0), zeros(Int64,0))
"""
function biconnected_components!(A::MatrixNetwork, articulation::Vector{Bool}, map::Vector{Int64})
    n=length(A.rp)-1
    rp=A.rp
    ci=A.ci
    components = length(map) >= length(ci)
    cn=1
    low=zeros(Int64,n)
    dt=-1*ones(Int64,n)
    pred=zeros(Int64,n)
    t=0
    largest_bcc=0
    edge_count_bcc=0
    rs=Tuple{Int,Int,Int}[]
    cs=Tuple{Int,Int,Int}[]
    root_children= 0
    art = length(articulation) >= n

    #start dfs at 1.
    for sv=1:n
        v=sv
        if dt[v]>0
            continue
        end
        root_children = 0
        low[v]=t
        dt[v]=t
        t=t+1
        w = ci[rp[v]]
        if (v == w) && (rp[v]+1 == rp[v+1])
            if components
                map[rp[v]]=cn
            end
            cn = cn+1
            if art
                articulation[v]=1
            end
        else
            push!(rs, (v,v,rp[v])) #add v,rp[v] to stack
        end

        while size(rs, 1) > 0
            (g_parent, parent, children_index) = rs[end]
            child_index = children_index
            if child_index != rp[parent+1]
               child = ci[child_index]
               children_index = children_index+1
               (g_parent, parent, temp) = pop!(rs)
               push!(rs, (g_parent, parent, children_index))
               if (g_parent == child) || (parent == child)
                   continue   # skip tree edges and self loops
               end
               if dt[child] >= 0
                    if dt[child] <= dt[parent]
                        if dt[child] < low[parent] #Update Low
                            low[parent] = dt[child]
                        end
                        push!(cs, (parent,child,child_index))
                    end
                else
                    low[child] = dt[child] = t
                    t+=1
                    push!(rs, (parent, child, rp[child]))
                    push!(cs, (parent,child,child_index))
                end
            else
                pop!(rs)
                if size(rs,1) > 1
                    if low[parent] >= dt[g_parent]  # g_parent is an Articulation point
                        if art
                            articulation[g_parent]=1
                        end
                        while true
                            (u,v,child_index)=pop!(cs)
                            if components
                                map[child_index]=cn #assign the component number of the corresponding edge
                            end
                            if u == g_parent && v == parent
                               break
                            end
                        end
                        cn+=1
                    end
                    if low[parent] < low[g_parent]
                       low[g_parent] = low[parent]
                    end
                elseif size(rs,1) == 1
                    root_children +=1
                    while true
                        (u,v,child_index) = pop!(cs)
                        if components
                            map[child_index]=cn
                        end
                        if u == g_parent && v == parent
                            break
                        end
                    end
                    cn+=1
                end
            end
        end
    end
    cn=cn-1

    return cn
end

"""
`biconnected_components`
-----------------------
This function requires a symmetric matrix as input. Depending on the user input this
function returns either the biconnected component number associated with each edge
or articlation points or both.

Inputs
---------
- `A`: the adjacency matrix
- Optional Keyword inputs
  - `art=true`: returns the articulation points of the graph.
  - `components=true`:returns the biconnected component labels associated with each
  edge.

Returns
-------
- Returns a `Biconnected_components_output` type which includes
`map` : biconnected component labels associated with each edge,
`articulation_points`: boolean array that signifies whether a vertex is an articulation point and
`number`: Number of biconnected components in the graph.

Example
-------
A = load_matrix_network("biconnected_example")
B = MatrixNetwork(A)
bcc = biconnected_components(B)
map = bcc.map
articulation_vector = bcc.articulation_points
number_of_components = bcc.number
"""
function biconnected_components(A::MatrixNetwork; art::Bool = true, components::Bool = true)
    map = components ? zeros(Int64, length(A.ci)) : zeros(Int64, 0)
    articulation = art ? zeros(Bool, A.n) : zeros(Bool, 0)
    cn = biconnected_components!(A, articulation, map)
    return Biconnected_components_output(map,articulation,cn,A)
end

###############################
##    Conversion Functions    #
###############################

#CSC
biconnected_components(A::SparseMatrixCSC;kwargs...) = biconnected_components(MatrixNetwork(A);kwargs...)

#Triplet
biconnected_components(ei::Vector{Int64},ej::Vector{Int64};kwargs...) = biconnected_components(MatrixNetwork(ei,ej);kwargs...)

#CSR
biconnected_components(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64; kwargs...) where {T} = biconnected_components(MatrixNetwork(n,rp,ci,vals);kwargs...)

ext: jl


text:
module DetectionIoTools

using WAV
using ThreadTools
using Base.Threads
using Serialization
using Statistics

export mapsoundfiles, thing2file, save_interesting, save_interesting_concat

"""
    mapsoundfiles(f::F, files, [windowlength]) where F

map a Function over a list of wav files.

#Arguments:
- `f`: a Function
- `files`: DESCRIPTION
- `windowlength`: Optional. If supplied, `f` will be applied over non-overlapping windows for the sound, and a `Vector{Vector{RT}}` is returned, where `RT` is the return type of `f`.
"""
function mapsoundfiles(f::F,files,windowlength) where F
    GC.gc()
    # extension === nothing || (files = filter(f->splitext(f)[end] == extension, files))
    @time embeddings = tmap1(max(nthreads()-1,1),files) do file
        @info "Reading file $file"
        sound = wavread(file)[1]
        map(f,Iterators.partition(sound, windowlength))
    end
end

function mapsoundfiles(f::F,files) where F
    GC.gc()
    # extension === nothing || (files = filter(f->splitext(f)[end] == extension, files))
    @time embeddings = tmap1(max(nthreads()-1,1),files) do file
        @info "Reading file $file"
        sound = wavread(file)[1]
        f(sound)
    end
end

"""
    thing2file(thing, things [, files])

Finds the file corresponding to the `thing` amongst `things`. `things` can be a vector or a vector of vectors.
"""
function thing2file(thing, things::AbstractVector{<:AbstractVector}, files::AbstractVector = 1:length(things))
    findres = findfirst.(==(thing), things)
    fileno = findfirst(!=(nothing), findres)
    files[fileno]
end

function thing2file(thing, things, files::AbstractVector = 1:length(things))
    fileno = findfirst(==(thing), things)
    fileno === nothing && error("Did not find thing in things")
    files[fileno]
end

function WAV.wavplay(x::AbstractArray,fs)
    path = "/tmp/a.wav"
    wavwrite(x,path,Fs=fs)
    @info "Wrote to file $path"
    Threads.@spawn run(`totem $path`)
end

function WAV.wavplay(path::String)
    Threads.@spawn run(`totem $path`)
end

save_interesting(sounds, inds, args...) = save_interesting(sounds, findall(inds), args...)

"""
    save_interesting(sounds::Vector, inds::Vector{Int}, contextwindow=1)

Saves interesting wav files to disk, together with one concatenated file which contains all the itneresting sounds. The paths will be printed in the terminal.

# Arguments:
- `inds`: A list of indices to interesting sound clips
- `contextwindow`: Saves this many clips before and after
"""
function save_interesting(sounds::AbstractVector{T}, inds::Vector{Int}, contextwindow=1, fs=41000) where T
    tempdir = mktempdir()
    for ind ∈ inds
        extended_inds = max(1, ind-contextwindow):min(length(sounds), ind+contextwindow)
        sound = map(extended_inds) do i
            getsound(sounds[i])
        end
        sound = reduce(vcat, sound)[:]
        filename = T <: AbstractString ? splitpath(files[ind])[end] : string(ind)
        tempfile = joinpath(tempdir, filename*".wav")
        sound .-= mean(Float32.(sound))
        sound .*= 1/maximum(abs.(sound))
        wavwrite(sound, tempfile, Fs=fs)
        println(tempfile)
    end
    save_interesting_concat(sounds, inds, tempdir, fs=fs)
end

getsound(sound) = sound
function getsound(sound::AbstractString)
    ext = splitext(sound)[2]
    if ext == "wav"
        wavread(sound)
    else
        deserialize(sound)
    end
end

function save_interesting_concat(sounds, inds::Vector{Int}, tempdir=mktempdir(); fs=41000)
    sound = map(inds) do i
        sound = getsound(sounds[i])
        sound .-= mean(Float32.(sound))
        sound .*= 1/maximum(abs.(sound))
    end
    sound = reduce(vcat, sound)[:]
    tempfile = joinpath(tempdir, "concatenated.wav")
    wavwrite(sound, tempfile, Fs=fs)
    println(tempfile)
    tempdir
end

end


text:
#RSDeltaSigmaPort: Base plot generation tools
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#
const linlin = cons(:a, xyaxes=set(xscale=:lin, yscale=:lin))
const linlog = cons(:a, xyaxes=set(xscale=:lin, yscale=:log))
const loglin = cons(:a, xyaxes=set(xscale=:log, yscale=:lin))
const loglog = cons(:a, xyaxes=set(xscale=:log, yscale=:log))

const dfltline = cons(:a, line=set(style=:solid, color=:blue, width=2))

#==Waveform builders
===============================================================================#
function _wfrm__(x::Vector, y::Vector, sweepid::String)
	validatexy_vec(x, y)
	return DataF1(x,y)
end
_wfrm__(x::AbstractVector, y::Vector, sweepid::String) = _wfrm__(collect(x), y, sweepid)
function _wfrm__(x::Array{Float64,2}, y::Array{Float64,2}, sweepid::String)
	validatexy_1param(x, y)
	N = size(x, 1)
	wfrm = fill(DataRS{DataF1}, PSweep(sweepid, collect(1:N))) do i
		return DataF1(x[i,:], y[i,:]) #Returns from implicit "do" function
	end
	return wfrm
end

"""`waveform(x, y; sweepid="i")`

Create a waveform object with a given string to identify the sweep.
"""
waveform(x, y; sweepid::String="i") = _wfrm__(x, y, sweepid)

#==logsmooth
===============================================================================#
"""`(f, p) = logsmooth(X,inBin,nbin=8,n=3)`

Smooth the fft, X, and convert it to dB.

 - Use nbin bins from 0 to 3*inBin, thereafter increase bin sizes by a factor
   of 1.1, staying less than 2^10.
 - For the n sets of bins inBin+[0:2], 2*inBin+[0:2], ... n*inBin+[0:2], don't
   do averaging. This way, the noise BW and the scaling of the tone and its
   harmonics are unchanged.

Unfortunately, harmonics above the nth appear smaller than they
really are because their energy is averaged over many bins.
"""
function logsmooth(X, inBin; nbin::Int=8, n::Int=3)
#= Comment (MALaforge): VERIFYME
 - Goes infinite @inBin because startbin[] values jump back on 1st iteration of
   `for i in 1:n`.
 - Was this ententional to help user locate inBin?
=#
	N=length(X); N2=div(N,2)
	f1 = mod(inBin-1, nbin) + 1
	startbin = vcat(f1:nbin:inBin-1, inBin:inBin+2)
	for i in 1:n
		startbin = vcat(startbin, startbin[end]+1:nbin:i*inBin-1, i*inBin .+ (0:2))
	end
	m = startbin[length(startbin)]+nbin
	while m < N2
		startbin = vcat(startbin, m)
		nbin = min(nbin*1.1, 2^10)
		m = round(Int, m+nbin)
	end
	stopbin = vcat(startbin[2:length(startbin)] .- 1, N2)
	f = ((startbin .+ stopbin)/2 .- 1)/N
	p = zeros(size(f))

	for i in 1:length(f)
		p[i] = dbp(norm(X[startbin[i]:stopbin[i]])^2 / (stopbin[i]-startbin[i]+1))
	end

	return (f, p)
end

#==plotUsage()
===============================================================================#
"""`plotUsage(sv,colors=[:blue, :green, :red, :yellow])`
Plot the element usage for a multi-elemet DAC.
The colors are for sv = 1,-1,i,-i.
"""
function plotUsage(sv,colors=[:blue, :green, :red, :yellow])
	(M, T) = size(sv)
	T2 = ceil(Int, (T+1)/2)

	plot = cons(:plot, linlin, title = "", legend=false,
		xyaxes=set(xmin=0, xmax=T, ymin=0, ymax=M),
		labels=set(xaxis="Time", yaxis="Element Number"),
	)

	#Plot the grid
	x = [(0:T)'; (0:T)']
		x = x[:]
	y = [zeros(1,T2); M*ones(2,T2); zeros(1,T2)]
		y = y[1:2*(T+1)]
	yw = waveform(x, y) #Convert to waveform
	push!(plot,
		cons(:wfrm, yw, line=set(style=:solid, color=:black, width=2), label=""),
	)

	M2 = ceil(Int, (M+1)/2)
	x = [zeros(1,M2); T*ones(2,M2); zeros(1,M2)]
		x = x[1:2*(M+1)]
	y = [0:M; 0:M]
		y = y[:]
	yw = waveform(x, y) #Convert to waveform
	push!(plot,
		cons(:wfrm, yw, line=set(style=:solid, color=:black, width=2), label=""),
	)

	@warn("TODO (plotUsage): draw polygons.")

#=Create polygons
	#Not yet supported by EasyPlot (would have to manipulate InspectDR manually)
	for t in 1:T
		for i in 1:M
			if sv[i,t] == 1
				fill([t-1 t-1 t t],[i-1 i i i-1],colors(1));
			elseif sv(i,t) == -1
				fill([t-1 t-1 t t],[i-1 i i i-1],colors(2));
			elseif sv(i,t) == 1i
				fill([t-1 t-1 t t],[i-1 i i i-1],colors(3));
			elseif sv(i,t) == -1i
				fill([t-1 t-1 t t],[i-1 i i i-1],colors(4));
			end
		end
	end
=#
	return plot
end

#Last line


text:
@testset "uniform sampler" begin

	cc = dummyChoiceContext()

	@testset "getparams, numparams, paramranges, sample" begin

		@testset "using default construction" begin

			s = DataGenerators.UniformSampler()

			@testset "numparams and paramranges" begin
				@test DataGenerators.numparams(s) == 2
				prs = DataGenerators.paramranges(s)
				@test typeof(prs) <: Vector{Tuple{Float64,Float64}}
				@test prs == fill((-realmax(Float64), realmax(Float64)), 2)
			end

			@testset "default params" begin
				@test DataGenerators.getparams(s) == [-1.0, 1.0]
			end

			@mtestset "default sampling" reps=Main.REPS alpha=Main.ALPHA begin
				x, trace = DataGenerators.sample(s, (0,1), cc)
				@test typeof(x) <: Float64
				@test -1.0 <= x <= 1.0
		        @mtest_distributed_as Uniform(-1.0,1.0) x
			end

		end

		@testset "non-default construction" begin

			s = DataGenerators.UniformSampler([-100.7,129.3762])

			@testset "constructor with params" begin

				@test DataGenerators.getparams(s) == [-100.7,129.3762]

				@mtestset "is consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
					x, trace = DataGenerators.sample(s, (0,1), cc)
			        @mtest_distributed_as Uniform(-100.7,129.3762) x
				end

			end

		end

	end

	@testset "setparams" begin

	    s = DataGenerators.UniformSampler()
	    prs = DataGenerators.paramranges(s)
	    midparams = map(pr->robustmidpoint(pr[1],pr[2]), prs)

	    @testset "valid parameters $params" for params in [[8.2,50.78], [-3113.10, -228.99],  [-6.198, 6.23],]

	        DataGenerators.setparams!(s, params)

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Uniform(params[1],params[2]) x
			end

	    end

		@testset "setparams adjusts order of parameters" begin

			DataGenerators.setparams!(s, [50.77,8.3])
			@test DataGenerators.getparams(s) == [8.3,50.77]

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Uniform(8.3,50.77) x
			end

		end

		@testset "handles full Float64 range" begin
			params = [-realmax(Float64), realmax(Float64)]
			DataGenerators.setparams!(s, params)

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        # @mtest_distributed_as Uniform(-realmax(Float64),realmax(Float64)) x
				# Uniform returns Inf for this range, so instead:
				@test -realmax(Float64) <= x <= realmax(Float64)
				@mtest_that_sometimes x < -realmax(Float64)/2
				@mtest_that_sometimes -realmax(Float64)/2 <= x < 0
				@mtest_that_sometimes  0 <= x < realmax(Float64)/2
				@mtest_that_sometimes  realmax(Float64)/2 <= x
			end

		end

	    @testset "boundary values for parameter index $pidx" for pidx in 1:length(prs)

	        pr = prs[pidx]
	        params = copy(midparams)

			@testset "bound index $bidx" for bidx in 1:2

	            params[pidx] = pr[bidx]
		        DataGenerators.setparams!(s, params)

				@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
		        	x, trace = DataGenerators.sample(s, (0,1), cc)
					sortedparams = sort(params)
			        @mtest_distributed_as Uniform(sortedparams[1],sortedparams[2]) x
				end

				@testset "range check exception" begin
		            params[pidx] = bidx == 1 ? prevfloat(pr[bidx]) : nextfloat(pr[bidx])
		            @test_throws ErrorException DataGenerators.setparams!(s, params)
				end

			end

	    end

		@testset "lower bound equals upper bound" begin

			DataGenerators.setparams!(s, [0.712, 0.712])

			@test DataGenerators.getparams(s) == [0.712, 0.712]

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
				@test x == 0.712
			end

		end

		@testset "samples boundaries of support" begin

			DataGenerators.setparams!(s, [87.23, nextfloat(nextfloat(87.23))])

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
				@mtest_values_are [87.23,nextfloat(87.23),nextfloat(nextfloat(87.23))] x
			end

		end

	    @testset "wrong number of parameters" begin
	        @test_throws ErrorException DataGenerators.setparams!(s, midparams[1:end-1])
	        @test_throws ErrorException DataGenerators.setparams!(s, [midparams; 0.5])
	    end

	end

	@testset "estimateparams" begin

	   	otherparams = [-42.9, 42.2]

	    @testset "from parameters $params" for params in [[20.55, 29.12], [50.22, 50.22],]
		# equal bounds case is motivated because special handling of equal bound sets type of s.distribution to a number rather than an instance of Uniform

	        s1 = DataGenerators.UniformSampler(params)
	        s2 = DataGenerators.UniformSampler(otherparams)
	        traces = map(1:100) do i
	            x, trace = DataGenerators.sample(s1, (0,1), cc)
	            trace
	        end
	        estimateparams!(s2, traces)

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
				if params[1] != params[2]
					@mtest_distributed_as Uniform(params[1],params[2]) x
				else
					@test x == params[1]
				end
			end

	    end

	    @testset "too few traces" begin

	        params = [101.0, 106.0]
	        s1 = DataGenerators.UniformSampler(params)
	        s2 = DataGenerators.UniformSampler(otherparams)
	        traces = map(1:1) do i
	            x, trace = DataGenerators.sample(s1, (0,1), cc)
	            trace
	        end
	        estimateparams!(s2, traces)

			@mtestset "consistent with uniform" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as Uniform(-42.9,42.2) x
			end

	    end

	end

end


text:
t0 = time_ns()

# test
import CSV
import DataFrames
import DataFrames: DataFrame
import Random
import OnlineStats
import StatsBase

using Distributed # only for Tests?
using Statistics
using Pkg # only for Tests?
using Test # only for Tests?
# using Pkg; Pkg.activate(normpath(joinpath(pwd(),".."));Pkg.instantiate();#Pkg.test()
using DecisionTrees

println("Time for using & import: ",round((time_ns() - t0) / 1e9, digits=3),"s")

# set directories
global pkgdir = ""
global testdir = joinpath(dirname(@__FILE__), "") # joinpath(@__DIR__,"..","test",filename)
# the following if clause is for "direct execution" (e.g. in VScode where @__FILE__ does not seem to work)
if testdir == ""
    try
        global pkgdir = normpath(joinpath(dirname(pathof(DecisionTrees)), ".."))
        global testdir = joinpath(pkgdir, "test")
    catch
    end
else
    global pkgdir = normpath(joinpath(testdir, ".."))
end
datadir = joinpath(pkgdir, "data")

global graphvizexe
global graphvizsetup_is_working

Random.seed!(1239946)
DecisionTreesTests = @testset verbose = true "DecisionTrees" begin
    t0runtests = time_ns()
    println("DTM: Running tests...") # Time = $(t0runtests)")

    global eltypesData1 = [Int,	Float64,	Float64,	Float64,	Float64,	Float64,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	String,	String,	String,	String,	String,	String,	Int,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	Int,	String,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	Int,	Int,	Int,	String,	String,	String,	String,	String,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Float64,	Int,	Float64,	Float64,	Float64,	Float64,	Int,	Int,	Int,	Int,	Int,	Int]

    @info("DTM: Note that the tests produce a relatively long log and write a number of files in temporary directories. You can disregard the log, if all tests pass.")
    @warn("DTM: So far only very rudimentary tests have been implemented!")

    include("graphvizsetup.jl")

    include("smoketests.jl") #running a number of functions to see if anything 'smokes' (runs into an unexpected error)
    include("errors_and_warnings.jl") #testing a few warnings and errors
    include("functionalityTests.jl") #testing resulting values of models
    include("helperfunctions.jl")

    @testset verbose = true "Other Tests" begin
        @test 1 == 1
    end
    println("DTM: Testing finished. Time needed:")
    TimeForTesting = (time_ns() - t0runtests) / 1e9
    println((round(TimeForTesting, digits=0)), "s")
    println(round(TimeForTesting / 60, digits=1), "m")
    println(round(TimeForTesting / 3600, digits=3), "h")
end
# dump(DecisionTreesTests)


text:
using Test, TestPackage3

@static if VERSION > v"1.3"
  @test hello3("Julia") == "Hello, Julia"
elseif VERSION > v"1.2"
  @test domath3(2.0) ≈ 7.0
else
  multiply3(2.0) == 8.0
end


text:
using Distributions
using UnPack

"""
  Option

Any contract representing the right, but not the
obligation to buy or sell a derivative at some point in time.
"""
abstract type Option end
abstract type AmericanOption <: Option end
abstract type EuropeanOption <: Option end

# NOTE: American option doesn't have volatility (σ). Can it be Calculated?

"""
  AmericanCall Option

An option which can be exercised at any point up to maturity.
"""
struct AmericanCall <: AmericanOption
    S::Float64  # stock price
    K::Float64  # strike price
    u::Float64  # up-trend
    d::Float64  # down-trend
    r::Float64  # annual rate
    R::Float64  # risk-free rate => 1 + r
    q::Float64  # (R - d) / (u - d)
    T::Float64  # maturity in months
    t::Float64  # point in time
    function AmericanCall(
        S::Float64, K::Float64, u::Float64, d::Float64, r::Float64, T::Float64
    )
        return new(S, K, u, d, r / 12, 1 + r / 12, ((1 + r / 12) - d) / (u - d), T, 0.0)
    end
    function AmericanCall(
        S::Float64, K::Float64, u::Float64, d::Float64, r::Float64, T::Float64, t::Int64
    )
        return new(
            S, K, u, d, r / 12, 1 + r / 12, ((1 + r / 12) - d) / (u - d), T, Float64(t)
        )
    end
end

"""
  AmericanPut Option

An option which can be exercised at any point up to maturity.
"""
struct AmericanPut <: AmericanOption
    S::Float64  # stock price
    K::Float64  # strike price
    u::Float64  # up-trend
    d::Float64  # down-trend
    r::Float64  # annual rate
    R::Float64  # risk-free rate => 1 + r
    q::Float64  # (R - d) / (u - d)
    T::Float64  # maturity in months
    t::Float64  # point in time
    function AmericanPut(
        S::Float64, K::Float64, u::Float64, d::Float64, r::Float64, T::Float64
    )
        return new(S, K, u, d, r / 12, 1 + r / 12, ((1 + r / 12) - d) / (u - d), T, 0.0)
    end
    function AmericanPut(
        S::Float64, K::Float64, u::Float64, d::Float64, r::Float64, T::Float64, t::Int64
    )
        return new(
            S, K, u, d, r / 12, 1 + r / 12, ((1 + r / 12) - d) / (u - d), T, Float64(t)
        )
    end
end

"""
  EuropeanCall Option

An option which can only be exercised at maturity.
"""
struct EuropeanCall <: EuropeanOption
    S::Float64  # stock price
    K::Float64  # strike price
    r::Float64  # rate
    σ::Float64  # volatility
    T::Float64  # maturity period
    t::Float64  # point in time
    function EuropeanCall(S::Float64, K::Float64, r::Float64, σ::Float64, T::Float64)
        return new(S, K, r, σ, T, 0.0)
    end
    function EuropeanCall(
        S::Float64, K::Float64, r::Float64, σ::Float64, T::Float64, t::Int64
    )
        return new(S, K, r, σ, T, Float64(t))
    end
end

"""
  EuropeanPut Option

An option which can only be exercised at maturity.
"""
struct EuropeanPut <: EuropeanOption
    S::Float64  # stock price
    K::Float64  # strike price
    r::Float64  # rate
    σ::Float64  # volatility
    T::Float64  # maturity period
    t::Float64  # point in time
    function EuropeanPut(S::Float64, K::Float64, r::Float64, σ::Float64, T::Float64)
        return new(S, K, r, σ, T, 0.0)
    end
    function EuropeanPut(
        S::Float64, K::Float64, r::Float64, σ::Float64, T::Float64, t::Int64
    )
        return new(S, K, r, σ, T, Float64(t))
    end
end

# #######################
# Option Pricing formulas
# #######################

"""
    price(e::European)

Returns a tuple of the call & put premiums of a European Option.
"""
function price(e::EuropeanOption)
    @unpack S, K, r, σ, T, t = e
    d₁ = (log(S / K) + (r + 0.5σ^2) * (T - t)) / (σ * √(T - t))
    d₂ = d₁ - σ * √(T - t)
    C = S * cdf(Normal(), d₁) - K * exp(-r * T) * cdf(Normal(), d₂)
    P = C + K * exp(-r * T) - S  # put-call parity
    return (C, P)
end

# #######################
# Option Greeks
# #######################

# TODO: add functions for checking 'in' and 'out' of money

"""
    isat(o::Option)

Returns true if an option is "at the money".
"""
function isat(o::Option)
    return o.S == o.K
end

"""
    Δ(o::Option, t=0.0)

Returns Δ of a European call & put option at some point in time `t`.

Δ is the sensitivity of the option premium with respect to the underlying derivative's
price.

# Delta Hedging

- A portfolio of derivatives has a ``Δ`` equal to the sum of each component's ``Δ``
- A portfolio is delta neutral when its ``Δ = 0``

## Dynamic Hedging

A strategy of rebalancing a portfolio's components such that the portfolio delta is zero.
"""
function Δ(o::Option)
    @unpack S, K, r, σ, T, t = o
    d₁ = (log(S / K) + (r + 0.5σ^2) * (T - t)) / (σ * √(T - t))
    C = cdf(Normal(), d₁)
    P = -cdf(Normal(), -d₁)
    return (C, P)
end

"""
    Γ(o::Option)

The Gamma measures the rate of change of a derivative's ``Δ`` with respect to its
underlying stock price.
"""
function Γ(o::Option)
    @unpack S, K, r, σ, T, t = o
    d₁ = (log(S / K) + (r + 0.5σ^2) * (T - t)) / (σ * √(T - t))
    return pdf(Normal(), d₁) / (S * σ * √(T - t))
end

"""
    Θ(o::Option)

The Theta measures the rate of change of a derivative's price with respect to time (daily).
"""
function Θ(o::Option)
    @unpack S, K, r, σ, T, t = o
    d₁ = (log(S / K) + (r + 0.5σ^2) * (T - t)) / (σ * √(T - t))
    d₂ = d₁ - σ * √(T - t)
    first = (S * pdf(Normal(), d₁) * σ) / (2 * √(T - t))
    second = r * K * exp(-r * (T - t)) * cdf(Normal(), d₂)
    return (-first - second) / 365
end

# NOTE: Might make sense to re-work method to calculate derivative value.

"""
    dprice(o::Option, δS, δt)

Return the value of a derivative using its option Greeks.

`dprice` returns a tuple of values for call & put respectively.

We know that:

'''math
\begin{aligned}
  Δ &= \frac{δP}{δS} \\
  Γ &= \frac{δ^{2}P}{δS^{2}} \\
  Θ &= \frac{δP}{δt} \\
\end{aligned}

∴ f ≈ Δ * δS + \frac{1}{2}Γ * (δS)^{2} + Θ * δt
'''
"""
function dprice(o::Option, δS, δt)
    return (cprice(o, δS, δt), pprice(o, δS, δt))
end

"""
    cprice(o::Option, δS, δt)

Return underlying derivative value of call option.
"""
function cprice(o::Option, δS, δt)
    C, _ = price(o)
    return C + Δ(o)[1] * δS + 0.5Γ(o) * δS^2 + Θ(o) * δt
end

# WARNING: Need to verifty `pprice` formula

"""
    pprice(o::Option, δS, δt)

Return underlying derivative value of put option.
"""
function pprice(o::Option, δS, δt)
    _, P = price(o)
    return P + Δ(o)[2] * δS + 0.5Γ(o) * δS^2 + Θ(o) * δt
end


text:
module MIPSolver2017_05

using JuMP

export solveMIP

include("branch_and_bound.jl")

end


text:
using LibCEED, Printf, StaticArrays

include("common.jl")

function transform_mesh_coords!(dim, mesh_size, mesh_coords)
    @witharray coords = mesh_coords begin
        @inbounds @simd for i = 1:mesh_size
            # map [0,1] to [0,1] varying the mesh density
            coords[i] = 0.5 + 1.0/sqrt(3.0)*sin((2.0/3.0)*pi*(coords[i] - 0.5))
        end
    end
    exact_sa = (dim == 1 ? 2 : dim == 2 ? 4 : 6)
end

function run_ex2(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, gallery)
    ncompx = dim
    prob_size < 0 && (prob_size = 256*1024)

    mesh_order = max(mesh_order, sol_order)
    sol_order = mesh_order

    ceed = Ceed(ceed_spec)
    mesh_basis =
        create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS)
    sol_basis =
        create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS)

    nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size)
    println("Mesh size: ", nxyz)

    # Build CeedElemRestriction objects describing the mesh and solution discrete
    # representations.
    mesh_size, mesh_restr, _ = build_cartesian_restriction(
        ceed,
        dim,
        nxyz,
        mesh_order,
        ncompx,
        num_qpts,
        mode=RestrictionOnly,
    )
    sol_size, _, qdata_restr_i = build_cartesian_restriction(
        ceed,
        dim,
        nxyz,
        sol_order,
        div(dim*(dim + 1), 2),
        num_qpts,
        mode=StridedOnly,
    )
    sol_size, sol_restr, sol_restr_i = build_cartesian_restriction(
        ceed,
        dim,
        nxyz,
        sol_order,
        1,
        num_qpts,
        mode=RestrictionAndStrided,
    )
    println("Number of mesh nodes     : ", div(mesh_size, dim))
    println("Number of solution nodes : ", sol_size)

    # Create a CeedVector with the mesh coordinates.
    mesh_coords = CeedVector(ceed, mesh_size)
    set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)

    # Apply a transformation to the mesh.
    exact_sa = transform_mesh_coords!(dim, mesh_size, mesh_coords)

    # Create the Q-function that builds the diffusion operator (i.e. computes its
    # quadrature data) and set its context data.
    if !gallery
        @interior_qf build_qfunc = (
            ceed,
            dim=dim,
            (J, :in, EVAL_GRAD, dim, dim),
            (w, :in, EVAL_WEIGHT),
            (qdata, :out, EVAL_NONE, dim*(dim + 1)÷2),
            begin
                Jinv = inv(J)
                qdata .= setvoigt(w*det(J)*Jinv*Jinv')
            end,
        )
    else
        build_qfunc = create_interior_qfunction(ceed, "Poisson$(dim)DBuild")
    end

    # Create the operator that builds the quadrature data for the diffusion
    # operator.
    build_oper = Operator(
        ceed,
        qf=build_qfunc,
        fields=[
            (gallery ? :dx : :J, mesh_restr, mesh_basis, CeedVectorActive()),
            (gallery ? :weights : :w, ElemRestrictionNone(), mesh_basis, CeedVectorNone()),
            (:qdata, qdata_restr_i, BasisCollocated(), CeedVectorActive()),
        ],
    )

    # Compute the quadrature data for the diffusion operator.
    elem_qpts = num_qpts^dim
    num_elem = prod(nxyz)
    qdata = CeedVector(ceed, num_elem*elem_qpts*div(dim*(dim + 1), 2))
    print("Computing the quadrature data for the diffusion operator ...")
    flush(stdout)
    apply!(build_oper, mesh_coords, qdata)
    println(" done.")

    # Create the Q-function that defines the action of the diffusion operator.
    if !gallery
        @interior_qf apply_qfunc = (
            ceed,
            dim=dim,
            (du, :in, EVAL_GRAD, dim),
            (qdata, :in, EVAL_NONE, dim*(dim + 1)÷2),
            (dv, :out, EVAL_GRAD, dim),
            begin
                dXdxdXdxT = getvoigt(qdata)
                dv .= dXdxdXdxT*du
            end,
        )
    else
        apply_qfunc = create_interior_qfunction(ceed, "Poisson$(dim)DApply")
    end

    # Create the diffusion operator.
    oper = Operator(
        ceed,
        qf=apply_qfunc,
        fields=[
            (:du, sol_restr, sol_basis, CeedVectorActive()),
            (:qdata, qdata_restr_i, BasisCollocated(), qdata),
            (:dv, sol_restr, sol_basis, CeedVectorActive()),
        ],
    )

    # Compute the mesh surface area using the diff operator:
    #                                             sa = 1^T \cdot abs( K \cdot x).
    print("Computing the mesh surface area using the formula: sa = 1^T.|K.x| ...")
    flush(stdout)

    # Create auxiliary solution-size vectors.
    u = CeedVector(ceed, sol_size)
    v = CeedVector(ceed, sol_size)
    # Initialize 'u' with sum of coordinates, x+y+z.
    u[] = 0.0
    @witharray_read(
        x_host = mesh_coords,
        size = (mesh_size÷dim, dim),
        @witharray(u_host = u, size = (sol_size, 1), sum!(u_host, x_host))
    )

    # Apply the diffusion operator: 'u' -> 'v'.
    apply!(oper, u, v)
    sa = witharray_read(x -> sum(abs, x), v, MEM_HOST)

    println(" done.")
    @printf("Exact mesh surface area    : % .14g\n", exact_sa)
    @printf("Computed mesh surface area : % .14g\n", sa)
    @printf("Surface area error         : % .14g\n", sa - exact_sa)
end

run_ex2(
    ceed_spec="/cpu/self",
    dim=3,
    mesh_order=4,
    sol_order=4,
    num_qpts=6,
    prob_size=-1,
    gallery=false,
)

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
# Julia wrapper for header: /usr/include/scip/cons_conjunction.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrConjunction(scip)
    ccall((:SCIPincludeConshdlrConjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsConjunction(scip, cons, name, nconss, conss, enforce, check, _local, modifiable, dynamic)
    ccall((:SCIPcreateConsConjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_CONS}}, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nconss, conss, enforce, check, _local, modifiable, dynamic)
end

function SCIPcreateConsBasicConjunction(scip, cons, name, nconss, conss)
    ccall((:SCIPcreateConsBasicConjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_CONS}}), scip, cons, name, nconss, conss)
end

function SCIPaddConsElemConjunction(scip, cons, addcons)
    ccall((:SCIPaddConsElemConjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_CONS}), scip, cons, addcons)
end


text:
include("AgilentDSOX4024A.jl")
include("AgilentDSOX4034A.jl")
include("scope_common.jl")

#const SignalGenAndScope = Union{AgilentDSOX4024A, Keysight33612A}

# Priority 1 when defined
#instrument_reset(obj::Instr{AgilentDSOX4024A})    = write(obj, "*RST for 402")
# Priority 2
#instrument_reset(obj::Instr{T}) where (T <: SignalGenAndScope)   = write(obj, "*RST for SpecialGroup")
# Priority 3
#instrument_reset(obj::Instr{T}) where (T <: Oscilloscope)   = write(obj, "*RST for Oscilloscope")

#instrument_reset(obj::Instr{AgilentDSOX4034A})    = write(obj, "*RST for 403")

# Include model specific files here


text:
################
# Example data
################
N = 50
x_uncertain = [UncertainValue(Normal, x, rand(Uniform(0.1, 0.8))) for x in rand(N)]
y_uncertain = [UncertainValue(Normal, y, rand(Uniform(0.1, 0.8))) for y in rand(N)]
x = UncertainValueDataset(x_uncertain)
y = UncertainValueDataset(y_uncertain)

time_uncertain = [UncertainValue(Normal, i, 1) for i = 1:length(x)];
time_certain = [CertainValue(i) for i = 1:length(x)];
timeinds_x = UncertainIndexDataset(time_uncertain)
timeinds_y = UncertainIndexDataset(time_certain)

X = UncertainIndexValueDataset(timeinds_x, x)
Y = UncertainIndexValueDataset(timeinds_y, y);

###########################################################
# SequentialResampling
###########################################################
seq = SequentialResampling(StrictlyIncreasing())

@test resample(X, seq) isa NTuple{2, Vector{<:Real}}
@test resample(X, seq)[1] |> length == N
@test resample(X, seq)[2] |> length == N
@test resample(Y, seq) isa NTuple{2, Vector{<:Real}}
@test resample(Y, seq)[1] |> length == N
@test resample(Y, seq)[2] |> length == N

###########################################################
# SequentialInterpolatedResampling
###########################################################
seq_intp = SequentialInterpolatedResampling(StrictlyIncreasing(), RegularGrid(0, N, 2))

@test resample(X, seq_intp) isa Tuple{AbstractRange, Vector{<:Real}}
@test resample(X, seq_intp)[1] |> length == length(0:2:N)
@test resample(X, seq_intp)[2] |> length == length(0:2:N)
@test resample(Y, seq_intp) isa Tuple{AbstractRange, Vector{<:Real}}
@test resample(Y, seq_intp)[1] |> length == length(0:2:N)
@test resample(Y, seq_intp)[2] |> length == length(0:2:N)

#################################
# ConstrainedValueResampling
#################################

c = ConstrainedValueResampling(10, TruncateStd(1))

r = resample(timeinds_y, c)
@test length(r) == c.n
@test length(r[1]) == N

r = resample(timeinds_x, c)
@test length(r) == c.n
@test length(r[1]) == N

r = resample(x, c)
@test length(r) == c.n
@test length(r[1]) == N

r = resample(y, c)
@test length(r) == c.n
@test length(r[1]) == N

#################################
# ConstrainedIndexValueResampling
#################################

c = ConstrainedIndexValueResampling(10, (TruncateStd(1), TruncateStd(1)))
rx = resample(X, c)
ry = resample(Y, c)

@test rx isa Vector{Tuple{Vector{Float64},Vector{Float64}}}
@test ry isa Vector{Tuple{Vector{Float64},Vector{Float64}}}

@test length(rx) == c.n
@test length(ry) == c.n
@test length(rx[1]) == 2
@test length(ry[1]) == 2
@test length(rx[1][1]) == N
@test length(ry[1][2]) == N
@test length(rx[2][1]) == N
@test length(ry[2][2]) == N


text:
module TestRemove

using ZulipReminderBot
using ZulipReminderBot: Adapter, ZulipRequest, process_remove, process
using ZulipReminderBot: heartbeat!, TimedMessage, Message, Opts
using Dates

using Test

import DBInterface

struct MockDB end
struct MockDBStmt
    s::String
end

function ZulipReminderBot.prepare(db::Adapter{MockDB, MockDB}, query)
    MockDBStmt(query)
end
function DBInterface.execute(stmt::MockDBStmt, x)
    if stmt.s == "DELETE FROM messages WHERE id = ? AND msg_sender_id = ?" && x == (1, 1)
        return []
    elseif stmt.s == "DELETE FROM messages WHERE id = ? AND msg_sender_id = ?" && x == (2, 1)
        return []
    end
end

function DBInterface.close!(stmt::MockDBStmt) end

function mockmessage(content; sender_id = 1, stream = "Stream1", topic = "Topic1", type = "stream")
    return Message(stream, topic, type, sender_id, content)
end

@testset "test process_remove command from same sender" begin
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1"))
    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 50, cin, cout)

    zreq = ZulipRequest()
    zreq.message.sender_id = 1
    zreq.data = "remove 1, 2"
    db = Adapter(MockDB(), MockDB)

    res = process_remove(zreq, db, cin, Dates.now(), nothing)
    @test res == "Messages were removed from the schedule"
    heartbeat!(sched, 51, cin, cout)
    @test length(sched) == 0
    @test !isready(cin)
    @test !isready(cout)
end

@testset "test process_remove command from different senders" begin
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1"))
    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 50, cin, cout)

    zreq = ZulipRequest()
    zreq.message.sender_id = 2
    zreq.data = "remove 1, 2"
    db = Adapter(MockDB(), MockDB)

    res = process_remove(zreq, db, cin, Dates.now(), nothing)
    @test res == "Messages were removed from the schedule"
    heartbeat!(sched, 51, cin, cout)
    @test length(sched) == 1
    @test sched[1].id == 1 && sched[1].msg.sender_id == 1
    @test !isready(cin)
    @test !isready(cout)
end

@testset "test process from same sender" begin
    opts = Opts("abc123", "", 0)
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1"; type = "private"))
    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 50, cin, cout)

    zreq = ZulipRequest()
    zreq.message.sender_id = 1
    zreq.data = "remove 1, 2"
    zreq.token = "abc123"
    db = Adapter(MockDB(), MockDB)

    res = process(zreq, db, cin, Dates.now(), opts)
    @test res == "{\"content\":\"Messages were removed from the schedule\"}"
    heartbeat!(sched, 51, cin, cout)
    @test length(sched) == 0
    @test !isready(cin)
    @test !isready(cout)
end

@testset "test process from different senders" begin
    opts = Opts("abc123", "", 0)
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1"; type = "private"))
    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 50, cin, cout)

    zreq = ZulipRequest()
    zreq.message.sender_id = 2
    zreq.data = "remove 1, 2"
    zreq.token = "abc123"
    db = Adapter(MockDB(), MockDB)

    res = process(zreq, db, cin, Dates.now(), opts)
    @test res == "{\"content\":\"Messages were removed from the schedule\"}"
    heartbeat!(sched, 51, cin, cout)
    @test length(sched) == 1
    @test sched[1].id == 1 && sched[1].msg.sender_id == 1
    @test !isready(cin)
    @test !isready(cout)
end

end # module


text:
#VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ReTerms

using DataArrays, DataFrames, NLopt, StatsBase

export LMM,ReMat,VectorReMat,ColMajorLowerTriangular,DiagonalLowerTriangular

export AIC, BIC, fixef, lowerbd, objective, pwrss

using Base.LinAlg.BlasInt

import Base: ==

include("utils.jl")
include("blockmats.jl")
include("remat.jl")
include("paramlowertriangular.jl")
include("cfactor.jl")
include("inject.jl")
include("pls.jl")

end # module


text:
using Test
using NNlib
using Zygote
using NNlibCUDA
using ForwardDiff: Dual
using Statistics: mean
using CUDA
CUDA.allowscalar(false)

include("test_utils.jl")

if CUDA.has_cuda()
    include("activations.jl")
    include("batchedmul.jl")
    include("upsample.jl")
    include("conv.jl")
    include("pooling.jl")
    include("softmax.jl")
    include("batchnorm.jl")
    include("scatter.jl")
    include("gather.jl")
end


text:
abstract type ElementoDetalhe end
struct Moeda <: ElementoDetalhe
    codigo::Symbol
    function Moeda(codigo::Symbol)
        global CONST_MOEDAS
        @assert codigo in CONST_MOEDAS "moeda invalida: $codigo"
        new(codigo)
    end
    Moeda(codigo::String) = Moeda(Symbol(codigo))
end
struct Posicao <: ElementoDetalhe
    codigo::Symbol
    function Posicao(codigo::Symbol)
        global CONST_POSICOES
        @assert codigo in CONST_POSICOES "posicao invalida: $codigo. should be :onshore or :offshore"
        new(codigo)
    end
    "1 = :onshore, 2 = :offshore"
    function Posicao(codigo::String)
        if codigo == "1"
            return Posicao(:onshore)
        elseif codigo == "2"
            return Posicao(:offshore)
        else
            error("codigo deve ser \"1\" ou \"2\"")
        end
    end
end
struct Pais <: ElementoDetalhe
    codigo::Symbol
    function Pais(codigo::Symbol)
        global CONST_PAISES
        @assert codigo in CONST_PAISES "pais invalido: $codigo"
        new(codigo)
    end
    Pais(codigo::String) = Pais(Symbol(codigo))
end

struct DetalheConta
    elementos::Vector{ElementoDetalhe}
    valor::Float64
    function DetalheConta(elementos::Vector{ElementoDetalhe}, valor::Float64)
        @assert length(elementos) > 0 "DetalheConta deve conter ao menos um elemento"
        new(elementos, valor)
    end
end

struct Conta
    codigo::String
    valor::Float64
    detalhes::Vector{DetalheConta}
    function Conta(codigo::String, valor::Float64, detalhes::Vector{DetalheConta})
        if length(detalhes) > 0
            sum_detalhes = sum([d.valor for d in detalhes])
            @assert valor ≈ sum_detalhes "Conta XML2011: soma dos detalhes ($sum_detalhes) deve ser igual ao valor ($valor)"
        end
        return new(codigo, valor, detalhes)
    end
    function Conta(codigo::String, valor::Float64)
        return new(codigo, valor, Vector{DetalheConta}())
    end
    Conta(codigo::String, valor::Int64) = Conta(codigo, Float64(valor))
    function Conta(codigo::String, detalhes::Vector{DetalheConta})
        valor = sum([d.valor for d in detalhes])
        return Conta(codigo, valor, detalhes)
    end
end

get_valor(conta::Conta) = conta.valor
function get_valor(conta::Conta, query::Vector{ElementoDetalhe})
    for d in conta.detalhes
        if all([(q in d.elementos) for q in query])
            return d.valor
        end
    end
    error("valor nao encontrado para os criterios $query")
end

const EMAIL_PATTERN = r"^[a-z0-9!#$%&'*+/=?^_`{|}~-]+(?:\.[a-z0-9!#$%&'*+/=?^_`{|}~-]+)*@(?:[a-z0-9](?:[a-z0-9-]*[a-z0-9])?\.)+[a-z0-9](?:[a-z0-9-]*[a-z0-9])?$"i

struct Responsavel
    nome::String
    telefone::String
    email::String
    function Responsavel(nome::String, telefone::String, email::String)
        @assert occursin(EMAIL_PATTERN, email) "email invalido"
        return new(nome, telefone, email)
    end
end

abstract type TipoEnvio end
struct Inclusao <: TipoEnvio end
encode(tipo::Inclusao)::String = "I"
struct Substituicao <: TipoEnvio end
encode(tipo::Substituicao)::String = "S"
function decode_tipo_envio(tipo::String)
    if tipo == "I"
        return Inclusao()
    elseif tipo == "S"
        return Substituicao()
    else
        error("tipo envio invalido: $tipo")
    end
end

struct Doc2011
    data::Date
    cnpj::String
    tipo::TipoEnvio
    responsavel::Responsavel
    contas::Vector{Conta}
    function Doc2011(data::Date, cnpj::String, tipo::TipoEnvio, responsavel::Responsavel, contas::Vector{Conta})
        doc = new(data, cnpj, tipo, responsavel, contas)
        validar(doc)
        return doc
    end
end

function get_conta(doc::Doc2011, codigo::String)
    for conta in doc.contas
        if conta.codigo == codigo
            return conta
        end
    end
    error("conta $codigo inexistente")
end

function has_conta(doc::Doc2011, codigo::String)
    for conta in doc.contas
        if conta.codigo == codigo
            return true
        end
    end
    return false
end

get_valor(doc::Doc2011, codigo::String) = get_valor(get_conta(doc, codigo))

function get_valor(doc::Doc2011, codigo::String, default::Float64)
    if has_conta(doc, codigo)
        return get_valor(doc, codigo)
    else
        return default
    end
end


text:
nn = Chain(Dense(21, 21, tanh), Dense(21, 21))
nn1 = FastChain(FastDense(21, 21, tanh), FastDense(21, 21))

X = randn(Float32, 21, 10)
Y = rand(Float32, 21, 10)
KitML.sci_train!(nn, (X, Y), ADAM())
KitML.sci_train!(nn, Flux.Data.DataLoader((X, Y)), ADAM(); device = cpu, epoch = 1)
KitML.sci_train(nn1, (X, Y))

loss(p) = nn1(X, p) |> sum
p1 = initial_params(nn1)
KitML.sci_train(loss, p1)

cd(@__DIR__)
model = KitML.load_model("tfmodel.h5"; mode = :tf)
KitML.sci_train!(model, (randn(Float32, 1, 4), randn(Float32, 1, 1)))


text:
module Denom
    export convert_fracts

    function convert_fracts(lst::Array{Array{Int,1},1})::Array{Array{BigInt,1},1}
        rats = map( x -> x[1]//x[2], lst )
        lcmt = lcm( map( BigInt∘denominator, rats ) )
        collect( map( x -> [lcmt*x,lcmt], rats ) )
    end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
function calculateSplitValue(a::DifferenceSplit, fname::Symbol, number_of_num_features::Int, labellist::Vector{T}, sumnumerator::Array{Float64,1}, sumdenominator::Array{Float64,1}, sumweight::Array{Float64,1}, countlistfloat::Array{Float64,1}, minweight::Float64, subs::DTSubsets) where T <: Unsigned
# here randomweight==0
# for subsets, exhaustive search with flipping members (gray code) or "increasing" subset search ({1}, {1,2}, {1,2,3}, .... {1,2,3, ....., n-1,2})
# all input lists (labellist,sumnumerator,sumdenominator,sumweight,countlistfloat) need to be sorted in the same manner
# convention, in the beginning everything is on the right side
    elementsInLeftChildBV = BitVector(undef, length(labellist));fill!(elementsInLeftChildBV, false) # indicates which classes belong to right child
    chosen_subset_bitarray = similar(elementsInLeftChildBV) # BitVector(undef,0)
    numtot = sum(sumnumerator)
    denomtot = sum(sumdenominator)
    weighttot = sum(sumweight)
    weighttot_minw = weighttot - minweight
    sumnl = sumdl = sumwl = 0.0
    val = -Inf
    chosen_sumwl = NaN
    for i in subs
  # @show "vv=$(i)" #should be the element which switches
  # i is an index, it indicates which element of labellist will flip sides for the next calculation
        @inbounds if elementsInLeftChildBV[i]
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true) # updating needs to occur here before we copy the array (in case it is a better split than what we have seen so far)
    # move class from left to right side
            @inbounds sumnl -= sumnumerator[i]
	@inbounds sumdl -= sumdenominator[i]
	@inbounds sumwl -= sumweight[i]
        else
      # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true)
      # move class from right to left side
            @inbounds sumnl += sumnumerator[i]
            @inbounds sumdl += sumdenominator[i]
            @inbounds sumwl += sumweight[i]
        end
        if (sumwl > minweight) && (weighttot_minw > sumwl) # do we have enough exposure? is the split valid?
            valnew = abs(sumnl / sumdl - (numtot - sumnl) / (denomtot - sumdl))
	  if valnew > val
                val = valnew
                chosen_subset_bitarray .= elementsInLeftChildBV
                chosen_sumwl = sumwl
            end
        end
    end
    if isfinite(val)
  # warning: if the labellist is not sorted 1:n we need to be careful here!
        chosen_subset = labellist[chosen_subset_bitarray]
    else
        chosen_subset = Array{T}(undef, 0)
    end
    return val, chosen_subset, chosen_sumwl, weighttot - chosen_sumwl
end

# todo/tbd assert eltype of lablelist is a UInt
function calculateSplitValue(a::DifferenceSplit, fname::Symbol, number_of_num_features::Int, labellist::Vector{T}, sumnumerator::Array{Float64,1}, sumdenominator::Array{Float64,1}, sumweight::Array{Float64,1}, countlistfloat::Array{Float64,1}, minweight::Float64, subs::DTSubsets, feature_column_id::Int) where T <: Unsigned
# here randomweight>0
# for subsets, exhaustive search with flipping members (gray code) or "increasing" subset search ({1}, {1,2}, {1,2,3}, .... {1,2,3, ....., n-1,2})
# all input lists (labellist,sumnumerator,sumdenominator,sumweight,countlistfloat) need to be sorted in the same manner
# convention, in the beginning everything is on the right side
    elementsInLeftChildBV = BitVector(undef, length(labellist));fill!(elementsInLeftChildBV, false) # indicates which classes belong to right child
    numtot = sum(sumnumerator)
    denomtot = sum(sumdenominator)
    weighttot = sum(sumweight)
    weighttot_minw = weighttot - minweight
    sumnl = sumdl = sumwl = 0.0
    this_splitlist = Array{Splitdef{T}}(undef, 0)

    for i in subs
    # i is an index, it indicates which element of labellist will flip sides for the next calculation
        @inbounds if elementsInLeftChildBV[i]
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true) # updating needs to occur here before we copy the array (in case it is a better split than what we have seen so far)
    # move class from left to right side
            @inbounds sumnl -= sumnumerator[i]
	@inbounds sumdl -= sumdenominator[i]
	@inbounds sumwl -= sumweight[i]
        else
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            @inbounds elementsInLeftChildBV[i] = xor(elementsInLeftChildBV[i], true)
      # move class from right to left side
            @inbounds sumnl += sumnumerator[i]
            @inbounds sumdl += sumdenominator[i]
            @inbounds sumwl += sumweight[i]
        end
        if (sumwl > minweight) && (weighttot_minw > sumwl) # do we have enough exposure? is the split valid?
            valnew = abs(sumnl / sumdl - (numtot - sumnl) / (denomtot - sumdl))
            feature_column_id2 = feature_column_id < 0 ? abs(feature_column_id) + number_of_num_features : feature_column_id
	  push!(this_splitlist, Splitdef{T}(feature_column_id, feature_column_id2, fname, labellist[elementsInLeftChildBV], valnew, sumwl, weighttot - sumwl))
        end
    end
    return this_splitlist
end


text:
function canFgBg(effect::Maple.Effect)
    return true, true
end

editingOptions(entity::Maple.Effect) = Dict{String, Any}()
editingIgnored(entity::Maple.Effect) = String[]
editingOrder(entity::Maple.Effect) = String[
    "name", "only", "exclude", "tag",
    "flag", "notflag"
]

function registerPlacements!(placements::Array{Type{Maple.Effect{T}} where T, 1}, loaded::Array{String, 1})
    empty!(placements)

    for modul in loaded
        if hasModuleField(modul, "placements")
            res = getModuleField(modul, "placements")

            if isa(res, Array)
                append!(placements, res)

            else
                push!(placements, res)
            end
        end
    end
end

const loadedEffects = joinpath.(abs"effects", readdir(abs"effects"))
const effectPlacements = Array{Type{Maple.Effect{T}} where T, 1}()

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
function rm_glb_pca(x::Neuro1DRealSignal; remove_ratio = 0.9)::Neuro1DRealSignal
    local _x = x.signal'
    local M = fit(PCA, _x, pratio=1.0)
    local p_ratio = M.prinvars ./ M.tprinvar

    local cum_ratio = 0
    local pc_idx = 0

    for i = 1:length(p_ratio)
        local _cum_ratio = cum_ratio + p_ratio[i]

        if (_cum_ratio >= remove_ratio)
            pc_idx = i
            break
        end

        cum_ratio = _cum_ratio
    end

    local new_mod = PCA(
        M.mean,
        M.proj[:, pc_idx:end],
        M.prinvars[pc_idx:end],
        sum(M.prinvars[pc_idx:end]),
        M.tvar)

    Neuro1DRealSignal(
        signal = convert(
            Matrix{Float64},
            MultivariateStats.reconstruct(new_mod, transform(new_mod, _x)
        )'),
        sample_rate = x.sample_rate
    )
end

function rm_glb_mean(x::Neuro1DRealSignal)::Neuro1DRealSignal
    glb_μ = mean(x.signal, dims = 2)

    βs = cov(x.signal, glb_μ) ./ var(glb_μ)

    Neuro1DRealSignal(
        signal = x.signal - glb_μ * βs',
        sample_rate = x.sample_rate
    )
end

function rm_glb(
    x::Neuro1DRealSignal,
    method = :glb_mean;
    kwargs...)::Neuro1DRealSignal

    if method == :glb_mean
        rm_glb_mean(x)
    elseif method == :pca
        rm_glb_pca(x; kwargs...)
    else
        throw(ArgumentError("invalid method"))
    end
end


text:
using CSV
using Plots

include("mlp.jl")

using .MLP

data = CSV.File("data/breast-cancer-wisconsin.data", header=false)
data = broadcast(collect.(data)) do x
    if x[7] == "?"
        x[7] = 0
    else
        x[7] = parse(Int, x[7])
    end
    x
end
data = collect(hcat(convert(Array{Array{Int64, 1}, 1}, data)...)')

d_tra = Dict(
    2 => [0, 1],
    4 => [1, 0],
)

data_X = data[:, 2:end-1] / 10
data_Y = collect(hcat(broadcast(data[:, end]) do x
    d_tra[x]
end...)')

data_X = Float64.(data_X)
data_Y = Float64.(data_Y)

let quz = "question-2", n_hidden_layers=5, n_epochs=800, η=0.9
    mlp = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=n_epochs, η=η)
    fit!(mlp)
    save_figures(mlp, "figs/$(quz)_")
    println("$(quz) done.")
end

let quz = "question-3", n_hidden_layers=5, η=0.9
    mlp = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=800, η=η)
    fit!(mlp)
    mlp_10 = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=10, η=η)
    fit!(mlp_10)
    mlp_100_000 = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=100_000, η=η)
    fit!(mlp_100_000)

    plt_1 = plot(
        mlp.loss_history,
        label="n_epochs=800, last=$(mlp.loss_history[end])",
        title="loss",
        legend=:outertopright,
        size=(1024, 512),
    )
    plot!(plt_1, mlp_10.loss_history, label="n_epochs=10, last=$(mlp_10.loss_history[end])")
    plot!(plt_1, mlp_100_000.loss_history, label="n_epochs=10,000, last=$(mlp_100_000.loss_history[end])")
    savefig(plt_1, "figs/$(quz)_loss_fig.png")

    plt_2 = plot(
        [mlp.train_accuracy_history mlp.test_accuracy_history],
        label=["trainset, n_epochs=800, last=$(mlp.train_accuracy_history[end])" "testset, n_epochs=800, last=$(mlp.test_accuracy_history[end])"],
        title="accuracy",
        legend=:outertopright,
        size=(1024, 512),
    )
    plot!(
        plt_2,
        [mlp_10.train_accuracy_history mlp_10.test_accuracy_history],
        label=["trainset, n_epochs=10, last=$(mlp_10.train_accuracy_history[end])" "testset, n_epochs=10, last=$(mlp_10.test_accuracy_history[end])"],
    )
    plot!(
        plt_2,
        [mlp_100_000.train_accuracy_history mlp_100_000.test_accuracy_history],
        label=["trainset, n_epochs=100,000, last=$(mlp_100_000.train_accuracy_history[end])" "testset, n_epochs=100,000, last=$(mlp_100_000.test_accuracy_history[end])"],
    )
    savefig(plt_2, "figs/$(quz)_accuracy_fig.png")
    println("$(quz) done.")
end

let quz = "question-4", n_epochs=800, η=0.9
    mlp = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=5, n_epochs=n_epochs, η=η)
    fit!(mlp)
    mlp_2 = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=2, n_epochs=n_epochs, η=η)
    fit!(mlp_2)
    mlp_30 = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=30, n_epochs=n_epochs, η=η)
    fit!(mlp_30)

    plt_1 = plot(
        mlp.loss_history,
        label="n_hidden_layers=5, last=$(mlp.loss_history[end])",
        title="loss",
        legend=:outertopright,
        size=(1024, 512),
    )
    plot!(plt_1, mlp_2.loss_history, label="n_hidden_layers=2, last=$(mlp_2.loss_history[end])")
    plot!(plt_1, mlp_30.loss_history, label="n_hidden_layers=30, last=$(mlp_30.loss_history[end])")
    savefig(plt_1, "figs/$(quz)_loss_fig.png")

    plt_2 = plot(
        [mlp.train_accuracy_history mlp.test_accuracy_history],
        label=["trainset, n_hidden_layers=5, last=$(mlp.train_accuracy_history[end])" "testset, n_hidden_layers=5, last=$(mlp.test_accuracy_history[end])"],
        title="accuracy",
        legend=:outertopright,
        size=(1024, 512),
    )
    plot!(
        plt_2,
        [mlp_2.train_accuracy_history mlp_2.test_accuracy_history],
        label=["trainset, n_hidden_layers=2, last=$(mlp_2.train_accuracy_history[end])" "testset, n_hidden_layers=2, last=$(mlp_2.test_accuracy_history[end])"],
    )
    plot!(
        plt_2,
        [mlp_30.train_accuracy_history mlp_30.test_accuracy_history],
        label=["trainset, n_hidden_layers=30, last=$(mlp_30.train_accuracy_history[end])" "testset, n_hidden_layers=30, last=$(mlp_30.test_accuracy_history[end])"],
    )
    savefig(plt_2, "figs/$(quz)_accuracy_fig.png")
    println("$(quz) done.")
end

let quz = "question-5", n_tries=100, n_hidden_layers=5, n_epochs=800, η=0.9
    mlp_l = SimpleMLP{Float64}[]
    for t ∈ 1:n_tries
        mlp = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=n_epochs, η=η)
        fit!(mlp)
        push!(mlp_l, mlp)
    end
    w_mlp = mlp_l[argmax(map(x -> x.loss_history[end], mlp_l))]
    save_figures(w_mlp, "figs/$(quz)_")
    println("$(quz) done.")
end

let quz = "question-6", n_hidden_layers=5, n_epochs=800, η=0.9
    mlp = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=n_epochs, η=η)
    fit!(mlp)
    mlp_a = SimpleMLP{Float64}(data_X, data_Y, n_hidden_layers=n_hidden_layers, n_epochs=n_epochs, η=η)
    adaptive_fit!(mlp_a)

    plt_1 = plot(
        mlp.loss_history,
        label="Gradient Descent, last=$(mlp.loss_history[end])",
        title="loss",
        legend=:outertopright,
        size=(1024, 512),
    )
    plot!(plt_1, mlp_a.loss_history, label="AdaGrad, last=$(mlp_a.loss_history[end])")
    savefig(plt_1, "figs/$(quz)_loss_fig.png")

    plt_2 = plot(
        [mlp.train_accuracy_history mlp.test_accuracy_history],
        label=["trainset, Gradient Descent, last=$(mlp.train_accuracy_history[end])" "testset, Gradient Descent, last=$(mlp.test_accuracy_history[end])"],
        title="accuracy",
        legend=:outertopright,
        size=(1024, 512),
    )
    plot!(
        plt_2,
        [mlp_a.train_accuracy_history mlp_a.test_accuracy_history],
        label=["trainset, AdaGrad, last=$(mlp_a.train_accuracy_history[end])" "testset, AdaGrad, last=$(mlp_a.test_accuracy_history[end])"],
    )
    savefig(plt_2, "figs/$(quz)_accuracy_fig.png")
    println("$(quz) done.")
end


text:
# utility functions

function get_first_nz_idx(x::Array{T}, mask = nothing) where {T<:Real}
    first_nz_idx = 0
    for i in 1:length(x)
        mask != nothing && !mask[i] && continue
        first_nz_idx = i
        abs(x[i]) > eps() && break
    end
    first_nz_idx
end

"""Compute elementary symmetric polynomials for given Λ and k.
"""
function elem_symm_poly(Λ::Array{Float64}, k::Int)
    N = length(Λ)
    poly = zeros(k + 1, N + 1)
    poly[1, :] .= 1
    for l in (1:k) .+ 1
        for n in (1:N) .+ 1
            poly[l, n] = poly[l, n - 1] + Λ[n - 1] * poly[l - 1, n - 1]
        end
    end
    poly
end


text:
using DiffEqBase: set_t!, set_u!, set_ut!
using OrdinaryDiffEq

println("First")
# set_X!(integrator, integrator.X) should not change the result.
@testset "Trivial $setter ($alg, inplace=$iip)" for alg in [RK4, Trapezoid],
    setter in [set_t!, set_u!, set_ut!],
    iip in [false, true]
  if iip
    f = (du,u,p,t) -> (du .= 2u)
    prob = ODEProblem{iip}(f,[1/2],(0.0,2.0))
  else
    f = (u,p,t) -> (2u)
    prob = ODEProblem{iip}(f,1/2,(0.0,2.0))
  end
  t_half = 1.0

  integrator1 = init(prob, alg(); save_everystep=false)
  integrator2 = init(prob, alg(); save_everystep=false)

  step!(integrator1, t_half)  # "inexact" stepping w/o tstops
  if setter === set_t!
    set_t!(integrator1, integrator1.t)
  elseif setter === set_u!
    set_u!(integrator1, integrator1.u)
  else
    set_ut!(integrator1, integrator1.u, integrator1.t)
  end
  solve!(integrator1)

  solve!(integrator2)

  @test integrator1.t == integrator2.t
  if alg === Trapezoid
    rtol = integrator1.opts.reltol * 100
    atol = integrator1.opts.abstol
    @test integrator1.u ≈ integrator2.u  rtol=rtol atol=atol
  else
    @test integrator1.u == integrator2.u
  end
end

println("Second")

@testset "Resolve with $setter ($alg, inplace=$iip)" for alg in [RK4, Trapezoid],
    setter in [set_t!, set_ut!],
    iip in [false, true]
  if iip
    f = (du,u,p,t) -> (du .= 2u * cos(2π * t))
    prob1 = ODEProblem{iip}(f,[1/2],(0.0,1.0))
  else
    f = (u,p,t) -> (2u * cos(2π * t))
    prob1 = ODEProblem{iip}(f,1/2,(0.0,1.0))
  end
  prob2 = remake(prob1; tspan=(0.0,2.0))
  integrator1 = init(prob1, alg(); save_everystep=false)
  integrator2 = init(prob2, alg(); save_everystep=false)

  solve!(integrator1)
  if setter === set_t!
    set_t!(integrator1, 0)
  else
    set_ut!(integrator1, integrator1.u, 0)
  end
  solve!(integrator1)

  solve!(integrator2)

  rtol = integrator1.opts.reltol
  atol = integrator1.opts.abstol
  if alg === Trapezoid
    rtol *= 100
  end
  @test integrator1.u ≈ integrator2.u  rtol=rtol atol=atol
end


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

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

# ╔═╡ 150d58e7-0e73-4b36-836c-d81eef531a9c
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using MLJ, MLJLinearModels, Plots, LinearAlgebra, Random, DataFrames, CSV, MLCourse
    import MLCourse: Polynomial
end

# ╔═╡ 2e9ce2a9-217e-4910-b6ce-d174f2f2668e
using PlutoUI; PlutoUI.TableOfContents()

# ╔═╡ 78bdd11d-b6f9-4ba6-8b2e-6189c4005bf1
md"# Ridge Regression (L2 Regularization)
"

# ╔═╡ 64b9cfa0-99f7-439b-b70e-f9266754ff74
md"In the following cell there is some custom code to run ridge regression and the lasso for the simple example of 1-dimensional input. In this example we penalize also the intercept β₀. For ridge regression the solution is
```math
\begin{eqnarray*}
\beta_1 &= \frac{\langle x y \rangle - \frac{\langle x\rangle \langle y\rangle}{1 + \lambda}}{\langle x^2\rangle - \frac{\langle x\rangle^2}{1 + \lambda} + \lambda}\\
\beta_0 &= \frac{\langle y \rangle - \beta_1 \langle x \rangle}{1 + \lambda}
\end{eqnarray*}
```
where ``\langle . \rangle`` denotes the average.

For the lasso, there we run a fixed point iteration that starts at the unregularized solution of linear regression and shrinks β₁ and β₀ towards zero until there is not change anymore.
"

# ╔═╡ 8bd483cc-f490-11eb-38a1-b342dd2551fd
begin
    function ridge_regression(x, y, λ)
       β₁ = (mean(x .* y) - mean(x) * mean(y)/(1 + λ))/
		    (mean(x.^2) - mean(x)^2/(1 + λ) + λ)
       β₀ = (mean(y) - β₁ * mean(x))/(1 + λ)
       (β₀ = β₀, β₁ = β₁)
    end
    function updateβ₀(x̄, ȳ, β₁, l)
        tmp = ȳ - β₁ * x̄
        abs(tmp) > l ? tmp - sign(tmp) * l : 0.
    end
    function updateβ₁(x̄, ȳ, x2, xy, β₀, l)
        tmp = (x̄ * ȳ - xy - x̄ * sign(β₀) * l)
        abs(tmp) > l ? (tmp - sign(tmp) * l)/(x̄^2 - x2) : 0
    end
    function lasso(x, y, λ)
        x̄ = mean(x)
        ȳ = mean(y)
        x2 = mean(x.^2)
        xy = mean(x .* y)
        β₁ = (x̄ * ȳ - xy)/(x̄^2 - x2)
        β₀ = ȳ - β₁ * x̄
        β₀old, β₁old = zero(β₀), zero(β₁)
        while β₀old != β₀ || β₁old != β₁
            β₀old, β₁old = β₀, β₁
            β₁ = updateβ₁(x̄, ȳ, x2, xy, β₀, λ)
            β₀ = updateβ₀(x̄, ȳ, β₁, λ)
        end
       (β₀ = β₀, β₁ = β₁)
    end
end;

# ╔═╡ 3a71618b-0d74-4480-b87b-c6ba87363b5d
md"We generate now some data to observe the effect of regularization when fitting this data."

# ╔═╡ 9e1e8284-a8c1-47a9-83d0-2d8fbd8ce005
n = 30; x = rand(n); y = 2.2x .+ .3 .+ .2randn(n);

# ╔═╡ 1009e251-59af-4f1a-9d0a-e96f4b696cad
md"λ₂ = $(@bind λ₂ Slider(0:.01:5, show_value = true))"

# ╔═╡ 50ac0b07-ffee-40c3-843e-984b3c628282
l2coefs = ridge_regression(x, y, λ₂)

# ╔═╡ 58746554-ca5a-4e8e-97e5-587a9c2aa44c
let r = λ₂ == 0 ? 6 : norm([l2coefs...]),
    ccol = plot_color(:blue, .3),
    path = hcat([[ridge_regression(x, y, l)...] for l in 0:.01:5]...)
    p1 = scatter(x, y, label = "data", xlabel = "x", ylabel = "y",
	             legend = :topleft)
    plot!(x -> l2coefs.β₀ + x * l2coefs.β₁, w = 3, label = "ridge regression")
    p2 = contour(-1:.1:3, -1:.1:3, (β₀, β₁) -> mean((β₀ .+ β₁*x .- y).^2),
                 label = "loss", title = "loss with constraints",
		         legend = :bottomright,
                 levels = 100, aspect_ratio = 1, ylims = (-1, 3), xlims = (-1, 3))
    plot!(t -> r * sin(t), t -> r * cos(t), 0:.001:2π,
          fill = (0, ccol), label = "constraint", color = ccol)
    plot!(path[1, :], path[2, :], label = "path", color = :blue, w = 3)
    scatter!([l2coefs.β₀], [l2coefs.β₁], label = "current fit", markersize = 6, color = :red)
    p3 = plot(0:.01:5, path[1, :], label = "β₀", xlabel = "λ₂", ylabel = "")
    plot!(0:.01:5, path[2, :], label = "β₁", ylims = (0, 2.4))
    scatter!([λ₂], [l2coefs.β₀], label = nothing, markersize = 6, color = :red)
    scatter!([λ₂], [l2coefs.β₁], label = nothing, markersize = 6, color = :red)
    p4 = contour(-1:.1:3, -1:.1:3, (β₀, β₁) -> mean((β₀ .+ β₁*x .- y).^2) + λ₂ * (β₀^2 + β₁^2),
                 label = "loss", title = "regularized loss",
                 levels = 100, aspect_ratio = 1, ylims = (-1, 3), xlims = (-1, 3))
    scatter!([l2coefs.β₀], [l2coefs.β₁], markersize = 6, label = nothing, color = :red)
    plot(p1, p4, p3, p2,
         layout = (2, 2), size = (700, 600), cbar = false)
end

# ╔═╡ f43a82e2-1145-426d-8e0e-5363d1c38ccf
md"# Lasso (L1 Regression)"

# ╔═╡ ff7cc2bf-2a38-46d2-8d11-529159b08c82
md"λ₁ = $(@bind λ₁ Slider(0:.01:1, show_value = true))"

# ╔═╡ 4841f9ba-f3d2-4c65-9225-bc8d0c0a9478
l1coefs = lasso(x, y, λ₁)

# ╔═╡ ed2b7969-79cd-43c8-bcdb-34dab89c2cb0
let r = λ₁ == 0 ? 10 : norm([l1coefs...], 1),
    ccol = plot_color(:blue, .3),
    path = hcat([[lasso(x, y, l)...] for l in 0:.01:1]...)
    p1 = scatter(x, y, label = "data", xlabel = "x", ylabel = "y", legend = :topleft)
    plot!(x -> l1coefs.β₀ + x * l1coefs.β₁, w = 3, label = "lasso")
    p2 = contour(-1:.1:3, -1:.1:3, (β₀, β₁) -> mean((β₀ .+ β₁*x .- y).^2),
                 label = "loss", title = "loss with constraints",
		         legend = :topright,
                 levels = 100, aspect_ratio = 1, ylims = (-1, 3), xlims = (-1, 3))
    plot!([0, r, 0, -r, 0], [r, 0, -r, 0, r],
          fill = (0, ccol), label = "constraint", color = ccol)
    plot!(path[1, :], path[2, :], label = "path", color = :blue, w = 3)
    scatter!([l1coefs.β₀], [l1coefs.β₁], label = "current fit", markersize = 6,
		     color = :red)
    p3 = plot(0:.01:1, path[1, :], label = "β₀", xlabel = "λ₁", ylabel = "")
    plot!(0:.01:1, path[2, :], label = "β₁", ylims = (0, 2.4))
    scatter!([λ₁], [l1coefs.β₀], label = nothing, markersize = 6, color = :red)
    scatter!([λ₁], [l1coefs.β₁], label = nothing, markersize = 6, color = :red)
    p4 = contour(-1:.1:3, -1:.1:3, (β₀, β₁) -> mean((β₀ .+ β₁*x .- y).^2)/2 + λ₁ * (abs(β₀) + abs(β₁)),
                 label = "loss", title = "regularized loss",
                 levels = 100, aspect_ratio = 1, ylims = (-1, 3), xlims = (-1, 3))
    scatter!([l1coefs.β₀], [l1coefs.β₁], markersize = 6, label = nothing, color = :red)
    plot(p1, p4, p3, p2,
         layout = (2, 2), size = (700, 600), cbar = false)
end

# ╔═╡ 4c3c816c-e901-4931-a27c-632b60291ad7
md"Instead of using the custom code to compute the ridge regression and the lasso we could have used some MLJ functions."

# ╔═╡ c1033416-334e-4b0e-b81e-6f9137402730
let mach = machine(RidgeRegressor(lambda = 3.82, penalize_intercept = true),
	               DataFrame(x = x), y)
	fit!(mach, verbosity = 0)
	fitted_params(mach)
end

# ╔═╡ 1dae5378-f3eb-4598-a060-445bfd8afe5e
md"Let us check that we get indeed the same result with our custom method."

# ╔═╡ 15a85810-ccbd-4aa3-98a5-fdcf68c97adb
ridge_regression(x, y, 3.82)

# ╔═╡ 0429acfe-d31e-427a-96d9-deddfa2c30f8
let mach = machine(LassoRegressor(lambda = .1,
	                              # usually the intercept is not penalized,
	                              # but here we do penalize it.
	                              penalize_intercept = true,
	                              # usually the default optimizer is quite good,
	                              # but here we decrease the tolerance to get
	                              # higher precision.
                                  solver = ISTA(tol = 1e-8)),
	            DataFrame(x = x), y)
	fit!(mach, verbosity = 0)
	fitted_params(mach)
end

# ╔═╡ 2ba5b965-1dce-4773-b6eb-fd838876674d
lasso(x, y, .1)

# ╔═╡ 6c87eb35-ddb3-44a3-b4ae-77a371e28960
md"There is also the `ElasticNetRegressor` that allows to fit with L1 and L2 penalties of different strengths. Look up the documentation to learn more about it."

# ╔═╡ ca394b88-06dc-4188-884a-50d7c180aa33
md"# Regularization Examples

## Polynomial Ridge Regression

In `MLJ` can apply ridge regression or the lasso to polynomial regression simply by replacing in the pipeline the `LinearRegressor` with a `RidgeRegressor` or a `LassoRegressor`, for example `mach = Polynomial(degree = 3) |> RidgeRegressor(lambda = 1e-3)`."

# ╔═╡ b45a9739-0f81-4c0c-a93b-434a5af91490
begin
    f(x) = .3 * sin(10x) + .7x
    function regression_data_generator(; n, seed = 3, rng = MersenneTwister(seed))
        x = rand(rng, n)
        DataFrame(x = x, y = f.(x) .+ .1*randn(rng, n))
    end
    regression_data = regression_data_generator(n = 50)
end;

# ╔═╡ 3d50111b-3a08-4a41-96ce-d77a8e37275d
md"degree = $(@bind degree Slider(0:20, show_value = true, default = 20))

$(@bind lambda Slider(-14:.1:.5, default = -4))
"

# ╔═╡ a54e3439-69b8-41c8-bfe0-4575795fb9b8
md"λ = $(lambda == -14 ? 0 : 10.0^lambda)"

# ╔═╡ bdbf0dfd-8da5-4e54-89c4-ef4d6b3796ce
let X = select(regression_data, Not(:y)), y = regression_data.y
    mach = fit!(machine(Polynomial(; degree) |> RidgeRegressor(lambda = lambda == -14 ? 0 : 10.0^lambda),
                        X, y), verbosity = 0)
    p1 = scatter(regression_data.x, y, label = "training data", ylims = (-.1, 1.1))
    plot!(f, label = "generator", c = :green, w = 2)
    grid = 0:.01:1
    pred = predict(mach, (x = grid,))
    plot!(grid, pred,
          label = "fit", w = 3, c = :red, legend = :topleft)
    annotate!([(.28, .6, "reducible error ≈ $(round(mean((pred .- f.(grid)).^2), sigdigits = 3))")])
end

# ╔═╡ fe2fe54f-0163-4f5d-9fd1-3d1aa3580875
begin
    model = Polynomial() |> RidgeRegressor()
    self_tuning_model = TunedModel(model = model,
                                   tuning =  Grid(goal = 500),
                                   resampling = CV(nfolds = 5),
                                   range = [range(model, :(polynomial.degree),
                                                  lower = 1, upper = 20),
                                            range(model, :(ridge_regressor.lambda),
                                                  lower = 1e-12, upper = 1e-3,
                                                  scale = :log)],
                                   measure = rmse)
    self_tuning_mach = machine(self_tuning_model,
                               select(regression_data, :x),
                               regression_data.y) |> fit!
end;

# ╔═╡ f5057d4a-1103-4728-becc-287d93d682ba
plot(self_tuning_mach)

# ╔═╡ bd54bfcd-f682-4b74-8a44-35463d421491
report(self_tuning_mach)

# ╔═╡ 596fd0f2-eee0-46ca-a203-e7cbac6f9788
let
    p1 = scatter(regression_data.x, regression_data.y,
                 label = "training data", ylims = (-.1, 1.1))
    plot!(f, label = "generator", c = :green, w = 2)
    grid = 0:.01:1
    pred = predict(self_tuning_mach, (x = grid,))
    plot!(grid, pred,
          label = "fit", w = 3, c = :red, legend = :topleft)
    annotate!([(.28, .6, "reducible error ≈ $(round(mean((pred .- f.(grid)).^2), sigdigits = 3))")])
end

# ╔═╡ 8e170a5a-9c46-413e-895d-796e178b69df
md"## Multiple Logistic Ridge Regression on the Spam Data

We load here the preprocessed spam data.
"

# ╔═╡ 8e542a48-ed28-4297-b2e8-d6a755a5fdf9
begin
    spam_train = CSV.read(joinpath(dirname(pathof(MLCourse)), "..", "data",
                                   "spam_preprocessed.csv"), DataFrame)
	spam_train.spam_or_ham = String.(spam_train.spam_or_ham)
    coerce!(spam_train, :spam_or_ham => OrderedFactor)
    spam_test = CSV.read(joinpath(dirname(pathof(MLCourse)), "..", "data",
                                  "spam_preprocessed_test.csv"), DataFrame)
	spam_test.spam_or_ham = String.(spam_test.spam_or_ham)
    coerce!(spam_test, :spam_or_ham => OrderedFactor)
end;

# ╔═╡ c5ef5d4e-200d-46d9-86fa-50af1896a6c3
md"The `LogisticClassifier` and the `MultinomialClassifier` have a `penalty` argument that can be used to enforce an L1 or L2 penalty. Look up the documentation to learn more about it."

# ╔═╡ d956613e-db32-488c-8ebb-fd61dfa31e59
spam_fit = fit!(machine(LogisticClassifier(penalty = :l2, lambda = 1e-5),
                        select(spam_train, Not(:spam_or_ham)),
                        spam_train.spam_or_ham));

# ╔═╡ 552f14fc-06c7-4c2a-9515-e64f28828b70
confusion_matrix(predict_mode(spam_fit, select(spam_train, Not(:spam_or_ham))),
                 spam_train.spam_or_ham)

# ╔═╡ ef701511-db7e-4dc0-8d31-ea14471943ab
confusion_matrix(predict_mode(spam_fit, select(spam_test, Not(:spam_or_ham))),
                 spam_test.spam_or_ham)

# ╔═╡ 19bdd76c-4131-422d-983e-1b29cd9edd30
md"We see that the test misclassification rate with regularization
is lower than in our original fit without regularization
(notebook \"Generalized Linear Regression\"; 48 false negatives and 48 false
positives). The misclassification rate on the training set is higher. This
indicates that unregularized logistic regression is too flexible for our spam
data set.
"

# ╔═╡ 13655a50-fbbb-46c7-bdf7-ed5644646966
md"## The Lasso Path for the Weather Data

For the Lasso it is often interesting to see the fitted parameter values for different regularization values (the Lasso path). In the following we use the package `GLMNet` to do so."

# ╔═╡ 1fa932c1-ce29-40ca-a8dc-e636aa2ecf66
weather = CSV.read(joinpath(@__DIR__, "..", "data", "weather2015-2018.csv"), DataFrame);

# ╔═╡ 470dc7f4-04a9-4253-8125-9112778021eb
import GLMNet: glmnet

# ╔═╡ ecf80b6a-1946-46fd-b1b4-bcbe91848e3c
begin
    weather_input = select(weather, Not(:LUZ_wind_peak))[1:end-5, :]
    weather_output = weather.LUZ_wind_peak[6:end]
    weather_fits = glmnet(Array(weather_input), weather_output)
end

# ╔═╡ 4652a904-5edb-463c-a046-5c5d378f7cca
let lambda = log.(weather_fits.lambda),
    col_names = names(weather_input)
    plotly()
    p = plot()
    for i in 1:size(weather_fits.betas, 1)
        plot!(lambda, weather_fits.betas[i, :], label = col_names[i])
    end
    plot!(legend = :outertopright, xlabel = "log(λ)", size = (700, 400))
    gr()
    p
end

# ╔═╡ 40bb385f-1cbd-4555-a8ab-544a67f33595
let lambda = log.(weather_fits.lambda)
    p1 = plot(lambda, 100 * weather_fits.dev_ratio, ylabel = "% variance explained")
    p2 = plot(lambda, reshape(sum(weather_fits.betas .!= 0, dims = 1), :),
              ylabel = "non-zero parameters",
              xlabel = "log(λ)")
    plot(p1, p2, layout = (2, 1), legend = false)
end

# ╔═╡ c9ed011c-8d36-4926-9ec4-84be3b4878d7
scatter(weather_input.BER_wind_peak, weather_output)

# ╔═╡ 8262b948-6d54-4348-87d1-4c762c74db30
md"# Exercises

## Conceptual

#### Exercise 1.
We review here the two formulations of regularization.
A standard approach to solve constraint optimization
problems makes use of Karush-Kuhn-Tucker (KKT) multipliers. For example, to find
the minimum of function ``f(x)`` under the constraint ``g(x) \leq s`` one can define
the loss function ``L(x, \lambda) = f(x) + \lambda(g(x) - s)`` where
``\lambda\geq0`` is a KKT multiplier.  Minimizing the loss both in ``x`` and
``\lambda`` amounts to solving the equations ``\frac{\partial L}{\partial x} = f'(x) + \lambda g'(x) = 0`` and ``\frac{\partial L}{\partial \lambda} = g(x) - s = 0``,
if the solution is on the boundary of the area defined by the inequality
constraint; otherwise one can find the solution by simply solving the
unconstrained problem, i.e. with ``\lambda = 0``. In this formulation one choses
the size ``s`` of the allowed area and finds ``\lambda`` by solving the equations.

Interestingly, the loss function ``L(x) = f(x) + \lambda
g(x)`` has exactly the same partial derivative in ``x`` as ``L(x, \lambda)`` and
therefore all critical points of ``L(x)`` have corresponding critical
points of ``L(x,\lambda)``. Because of this, regularization is often formulated
as \"adding a regularization term to the cost function\". For example, given the loss
    function of linear regression ``L(\beta) = \frac1n\sum_{i=1}^n(y_i - \beta_0 + \beta_1x_{i1} + \cdots + \beta_p x_{ip})^2`` one can define the
        L1-regularized loss function ``L_\mathrm{L1}(\beta) = L(\beta) + \lambda
\|\beta\|_1`` and choose a value for ``\lambda`` instead of the size ``s`` of the
allowed area.

1. Derive how ``\lambda`` in the second formulation depends on ``s`` in the first formulation for ridge regression with standardized one-dimensional input. *Hint:* use the analytical solution at the top of this notebook, the fact that ``\langle x \rangle = 0`` and ``\langle x^2\rangle = 0`` and note that ``\beta_0^2 + \beta_1^2 = s``, if the solution lies on the boundary.
1. Argue, why choosing ``\lambda = 0`` in the second formulation is equivalent to choosing a sufficiently large ``s`` in the first formulation.
2. Argue, why choosing ``\lambda = \infty`` in the second formulation is equivalent to choosing ``s = 0`` in the first formulation.

#### Exercise 2.
Consider a data set with as many data points as predictors ``n = p``.
Assume ``x_{ii} = 1`` and ``x_{ij} = 0`` for all ``i\neq j`` and arbitrary values
``y_i``. To simplify the problem further we perform regression without an
intercept. We would like to study L1- and L2-regularized multiple linear regression.

1. Write the mean squared error loss once with L1 regularization and once with L2 regularization for this setting and the fomulation of regularization with regularization constant ``\lambda``.
2. Show that in the case of L2 regularization the estimated coefficients take the form ``\hat \beta_j = y_j/(1 + \lambda)``.
3. Show that in the case of L1 regularization the estimated coefficients take the form ``\hat \beta_j = y_j - \lambda/2``, if ``y_j > \lambda/2``, ``\hat \beta_j = y_j + \lambda/2``, if ``y_j < -\lambda/2`` and ``\hat \beta_j = 0`` otherwise.
4. Write a brief summary on how the estimated coefficients ``\hat \beta_j`` are changed relative to the unregularized solution for both kinds of regularization.

## Applied

#### Exercise 1.
Create an artificial dataset with 20 points, 4 predictors ``X_1, X_2, X_3, X_4``
and ``Y = X_1 + \epsilon`` with ``\mathrm{Var}(\epsilon) = 0.1^2``.

1. Find with cross-validation and the lasso the best model.
2. Find with cross-validation and the ridge regression the best model.
3. Which of the two best models has the lowest reducible error?
4. Repeat the above 3 steps on an artificial data set with 20 points and 4 predictors with ``Y = 10X_1 + X_2 + .1 * X_3 + .01 * X_4 + \epsilon``.
"

# ╔═╡ e04c5e8a-15f8-44a8-845d-60acaf795813
MLCourse.list_notebooks(@__FILE__)

# ╔═╡ c48dff95-8028-4e97-8ec6-705ea2b9c72e
MLCourse.footer()

# ╔═╡ Cell order:
# ╟─2e9ce2a9-217e-4910-b6ce-d174f2f2668e
# ╠═150d58e7-0e73-4b36-836c-d81eef531a9c
# ╟─78bdd11d-b6f9-4ba6-8b2e-6189c4005bf1
# ╟─64b9cfa0-99f7-439b-b70e-f9266754ff74
# ╠═8bd483cc-f490-11eb-38a1-b342dd2551fd
# ╟─3a71618b-0d74-4480-b87b-c6ba87363b5d
# ╠═9e1e8284-a8c1-47a9-83d0-2d8fbd8ce005
# ╟─1009e251-59af-4f1a-9d0a-e96f4b696cad
# ╠═50ac0b07-ffee-40c3-843e-984b3c628282
# ╟─58746554-ca5a-4e8e-97e5-587a9c2aa44c
# ╟─f43a82e2-1145-426d-8e0e-5363d1c38ccf
# ╟─ff7cc2bf-2a38-46d2-8d11-529159b08c82
# ╠═4841f9ba-f3d2-4c65-9225-bc8d0c0a9478
# ╟─ed2b7969-79cd-43c8-bcdb-34dab89c2cb0
# ╟─4c3c816c-e901-4931-a27c-632b60291ad7
# ╠═c1033416-334e-4b0e-b81e-6f9137402730
# ╟─1dae5378-f3eb-4598-a060-445bfd8afe5e
# ╠═15a85810-ccbd-4aa3-98a5-fdcf68c97adb
# ╠═0429acfe-d31e-427a-96d9-deddfa2c30f8
# ╠═2ba5b965-1dce-4773-b6eb-fd838876674d
# ╟─6c87eb35-ddb3-44a3-b4ae-77a371e28960
# ╟─ca394b88-06dc-4188-884a-50d7c180aa33
# ╠═b45a9739-0f81-4c0c-a93b-434a5af91490
# ╟─3d50111b-3a08-4a41-96ce-d77a8e37275d
# ╟─a54e3439-69b8-41c8-bfe0-4575795fb9b8
# ╟─bdbf0dfd-8da5-4e54-89c4-ef4d6b3796ce
# ╠═fe2fe54f-0163-4f5d-9fd1-3d1aa3580875
# ╠═f5057d4a-1103-4728-becc-287d93d682ba
# ╠═bd54bfcd-f682-4b74-8a44-35463d421491
# ╟─596fd0f2-eee0-46ca-a203-e7cbac6f9788
# ╟─8e170a5a-9c46-413e-895d-796e178b69df
# ╠═8e542a48-ed28-4297-b2e8-d6a755a5fdf9
# ╟─c5ef5d4e-200d-46d9-86fa-50af1896a6c3
# ╠═d956613e-db32-488c-8ebb-fd61dfa31e59
# ╠═552f14fc-06c7-4c2a-9515-e64f28828b70
# ╠═ef701511-db7e-4dc0-8d31-ea14471943ab
# ╟─19bdd76c-4131-422d-983e-1b29cd9edd30
# ╟─13655a50-fbbb-46c7-bdf7-ed5644646966
# ╠═1fa932c1-ce29-40ca-a8dc-e636aa2ecf66
# ╠═470dc7f4-04a9-4253-8125-9112778021eb
# ╠═ecf80b6a-1946-46fd-b1b4-bcbe91848e3c
# ╟─4652a904-5edb-463c-a046-5c5d378f7cca
# ╟─40bb385f-1cbd-4555-a8ab-544a67f33595
# ╠═c9ed011c-8d36-4926-9ec4-84be3b4878d7
# ╟─8262b948-6d54-4348-87d1-4c762c74db30
# ╟─e04c5e8a-15f8-44a8-845d-60acaf795813
# ╟─c48dff95-8028-4e97-8ec6-705ea2b9c72e

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse