text:
#!/usr/bin/env julia

# testing isinside(), is point inside polygon
# all the random points are drawn only if they're inside one of the polygons

using Luxor

using Test

using Random
Random.seed!(42)

function ngon_poly(x, y, radius, sides::Int, orientation=0, action=:nothing; close=true)
    [Point(x+cos(orientation + n * 2pi/sides) * radius,
    y+sin(orientation + n * 2pi/sides) * radius) for n in 1:sides]
end

function bounding_b(x, y, o)
    return [Point(x-o, y-o), Point(x-o, y+o), Point(x+o, y+o), Point(x+o, y-o)]
end

function point_inside_polygon(fname)
    Drawing(1200, 1200, fname)
    origin()
    background("grey20")
    polys = [ngon_poly(x, y, 60, 5, pi/2) for x in -500:150:500, y in -500:150:500]
    for pol in polys
        randomhue()
        poly(pol, :fill)
    end
    # test whether random points are inside any of the polygons
    for i in 1:5000
        pt = randompoint(-600, -600, 600, 600)
        for pol in polys
            randomsize = rand(4:12)
            bp = bounding_b(pt.x, pt.y, randomsize)
            # check that all four corners of a 5x5 bounding box are inside this polygon
            # slow but this is Julia so it's still quick...
            if isinside(bp[1], pol) && isinside(bp[2], pol) && isinside(bp[3], pol) && isinside(bp[4], pol)
                randomhue()
                circle(pt.x, pt.y, randomsize, :fill)
            end
        end
    end

    @test finish() == true
    println("...finished test: output in $(fname)")
end

point_inside_polygon("point-inside-polygon.pdf")


text:
using Documenter

sections = ("Dynamics", "Systems", "Controls", "Analysis")
for section in sections
    isdir(joinpath(@__DIR__, "src", section)) || mkdir(joinpath(@__DIR__, "src", section))
end

isdir(joinpath(@__DIR__, "src", "Front Matter")) || mkdir(joinpath(@__DIR__, "src", "Front Matter"))
front = "Front Matter" => [
    "Welcome" => "index.md",
    "Chapter 1: Introduction" => joinpath("Front Matter", "Chapter 1: Introduction.md")
]

chapters = let topic = 1
    pages = [
        begin
            section_string = "Topic $topic: $section"
            topic += 1
            section_string
        end => [
            replace(page, ".md"=>"") => joinpath(section, page)
            for page in readdir(joinpath(@__DIR__, "src", section))
        ]
        for section in sections
    ]
end

makedocs(
    sitename = "Exploring Control Theory",
    format   = Documenter.HTML(),
    authors  = "Joe Carpinelli",
    pages    = vcat(front, chapters)
)

deploydocs(
    repo      = "https://github.com/cadojo/ControlTheoryNotes.jl",
    devbranch = "main",
    versions = ["stable" => "v^", "manual", "v#.#", "v#.#.#"]
)

ext: jl

max_stars_repo_name: cadojo/controls

max_issues_repo_name: cadojo/exploring-control-theory

max_forks_repo_name: cadojo/exploring-control-theory


text:
export
    Coord,
    CoordIm,
    AbstractCoord,

    get_dist,
    get_dist2d,
    get_value

abstract AbstractCoord
type Coord <: AbstractCoord
    x::Cdouble
    y::Cdouble
    z::Cdouble
    h::Cdouble
    p::Cdouble
    r::Cdouble

    Coord()=new()
    Coord(x::Real, y::Real, z::Real, h::Real=0.0, p::Real=0.0, r::Real=0.0) =
        new(x, y, z, h, p, r)
end
immutable CoordIm <: AbstractCoord
    x::Cdouble
    y::Cdouble
    z::Cdouble
    h::Cdouble
    p::Cdouble
    r::Cdouble

    CoordIm(x::Real, y::Real, z::Real, h::Real=0.0, p::Real=0.0, r::Real=0.0) =
        new(x, y, z, h, p, r)
end

copy(c::Coord) = Coord(c.x, c.y, c.z, c.h, c.p, c.r)
copy(c::CoordIm) = CoordIm(c.x, c.y, c.z, c.h, c.p, c.r)
deepcopy(c::AbstractCoord) = copy(c)
function copy!(dest::Coord, src::AbstractCoord)
    dest.x = src.x
    dest.y = src.y
    dest.z = src.z
    dest.h = src.h
    dest.p = src.p
    dest.r = src.r
    dest
end

Base.convert(::Type{CoordIm}, c::Coord) = CoordIm(c.x, c.y, c.z, c.h, c.p, c.r)
Base.convert(::Type{Coord}, c::CoordIm) = Coord(c.x, c.y, c.z, c.h, c.p, c.r)

function ==(a::AbstractCoord, b::AbstractCoord)
    isapprox(a.x, b.x) &&
    isapprox(a.y, b.y) &&
    isapprox(a.z, b.z) &&
    isapprox(a.h, b.h) &&
    isapprox(a.p, b.p) &&
    isapprox(a.r, b.r)
end

show(io::IO, coord::AbstractCoord) = @printf(io, "(%.16e, %.16e, %.16e, %.16e, %.16e, %.16e)", coord.x, coord.y, coord.z, coord.h, coord.p, coord.r)
print(io::IO, coord::AbstractCoord) = @printf(io, "(%.3f, %.3f, %.3f, %.3f, %.3f, %.3f)", coord.x, coord.y, coord.z, coord.h, coord.p, coord.r)

function get_dist(a::AbstractCoord, b::AbstractCoord)
    (ccall( (:coord_getDist, LIB_ODRMGR), Cdouble, (Ptr{Coord},Ptr{Coord}),
        pointer_from_objref(a), pointer_from_objref(b)))
end

function get_dist2d(a::AbstractCoord, b::AbstractCoord)
    (ccall( (:coord_getDist2d, LIB_ODRMGR), Cdouble, (Ptr{Coord},Ptr{Coord}),
        pointer_from_objref(a), pointer_from_objref(b)))
end

function *(a::AbstractCoord, d::Real)
    CoordIm(a.x * d,
              a.y * d,
              a.z * d,
              a.h * d,
              a.p * d,
              a.r * d)
end
function +(a::AbstractCoord, b::AbstractCoord)
    CoordIm(a.x + b.x,
              a.y + b.y,
              a.z + b.z,
              a.h + b.h,
              a.p + b.p,
              a.r + b.r)
end
function -(a::AbstractCoord, b::AbstractCoord)
    CoordIm(a.x - b.x,
              a.y - b.y,
              a.z - b.z,
              a.h - b.h,
              a.p - b.p,
              a.r - b.r)
end

function times!(a::Coord, d::Real)
    a.x *= d
    a.y *= d
    a.z *= d
    a.h *= d
    a.p *= d
    a.r *= d
    a
end
function plus!(a::Coord, b::AbstractCoord)
    a.x += b.x
    a.y += b.y
    a.z += b.z
    a.h += b.h
    a.p += b.p
    a.r += b.r
    a
end
function minus!(a::Coord, b::AbstractCoord)
    a.x -= b.x
    a.y -= b.y
    a.z -= b.z
    a.h -= b.h
    a.p -= b.p
    a.r -= b.r
    a
end
function times!(retval::Coord, a::AbstractCoord, d::Real)
    retval.x = a.x * d
    retval.y = a.y * d
    retval.z = a.z * d
    retval.h = a.h * d
    retval.p = a.p * d
    retval.r = a.r * d
    retval
end
function plus!(retval::Coord, a::AbstractCoord, b::AbstractCoord)
    retval.x = a.x + b.x
    retval.y = a.y + b.y
    retval.z = a.z + b.z
    retval.h = a.h + b.h
    retval.p = a.p + b.p
    retval.r = a.r + b.r
    retval
end
function minus!(retval::Coord, a::AbstractCoord, b::AbstractCoord)
    retval.x = a.x - b.x
    retval.y = a.y - b.y
    retval.z = a.z - b.z
    retval.h = a.h - b.h
    retval.p = a.p - b.p
    retval.r = a.r - b.r
    retval
end

function init!(coord::Coord)
    coord.x = 0.0
    coord.y = 0.0
    coord.z = 0.0
    coord.h = 0.0
    coord.p = 0.0
    coord.r = 0.0
    coord
end

get_value(coord::AbstractCoord) =
    ccall((:coord_getValue, LIB_ODRMGR), Cdouble, (Ptr{Void},), pointer_from_objref(coord))

max_stars_repo_name: sisl/OpenDriveManagerLite

max_issues_repo_name: sisl/OpenDriveManagerLite

max_forks_repo_name: sisl/OpenDriveManagerLite


text:
using Documenter, Behavior

makedocs(
    sitename="Behavior",
    pages = [
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Tutorial" => "tutorial.md",
        "Functions" => "functions.md",
        "Gherkin Experimental" => "gherkin_experimental.md"
    ])


text:
#=
Two LIF E neurons receive a sinusoidal type of input
from am input population
=#
push!(LOAD_PATH, abspath(@__DIR__,".."))

using Test
using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using BenchmarkTools

##

dt = 1E-3
Ttot = 100.0
# 2  LIF neurons
myτ = 0.2
vth = 10.
v_r = -5.0
τrefr= 0.0 # refractoriness
τpcd = 0.2 # post synaptic current decay
ps_e = S.PSLIF(myτ,vth,v_r,τrefr,τpcd,2)

##
# 10 modulated input neurons
function ratefun(t::Float64)
  ω = 20.0
  return 5*(1+sin(2π*t/ω))*0.5
end
τin = 0.233
ps_in = S.PSInputPoissonFt(ratefun,τin,10)

# connection from input to E
# 5 inputs go to one neuron, 5 to the other
w_e_input = sparse(vcat([ones(5)... , zeros(5)...]',
    [zeros(5)... , ones(5)...]'))
lmul!(0.33,w_e_input)
conn_e_in = S.ConnGeneralIF(w_e_input)
## Populations
 # the input evolves in time, so it's a population too
pop_in = S.UnconnectedPopulation(ps_in)
pop_e = S.Population(ps_e,(conn_e_in,ps_in))

myntw = S.RecurrentNetwork(dt,pop_in,pop_e)

## recorders : spikes
rec_spikes_e = S.RecSpikes(ps_e,50.0,Ttot)
rec_spikes_in =  S.RecSpikes(ps_in,50.0,Ttot)
# internal potential
krec = 1
rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot)

## Run

times = (0:myntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_state_e,rec_spikes_e,rec_spikes_in])
S.reset!.([ps_e,ps_in])
# initial conditions
fill!(ps_e.state_now,0.0)

for (k,t) in enumerate(times)
  rec_state_e(t,k,myntw)
  rec_spikes_e(t,k,myntw)
  rec_spikes_in(t,k,myntw)
  S.dynamics_step!(t,myntw)
end
# this is useful for visualization only
S.add_fake_spikes!(1.5vth,rec_state_e,rec_spikes_e)
##

spk_in_dict = S.get_spiketimes_dictionary(rec_spikes_in)

_ = let plt=plot(leg=false)
 for (neu,spkt) in pairs(spk_in_dict)
  scatter!(plt,_->neu,spkt ;
    marker=:vline,color=:white,ylims=(0.5,10.5))
 end
 plt
end

##
S.get_mean_rates(rec_spikes_in,dt,Ttot)

##
_ = let plt=plot(leg=false)
  x=rec_state_e.times
  y = rec_state_e.state_now
  plot!(x,[y[1,:] y[2,:]])
end


text:
using Statistics
using Combinatorics
# using ImageFiltering

"""
    rotate_img_around_center(img, angle = 5pi/6)

Function rotates a single image (or a frame) around the center of the image by
@angle radians.
"""
function rotate_img_around_center(img, angle = 5pi/6)
  θ = angle
  rot = recenter(RotMatrix(θ), [size(img)...] .÷ 2)  # a rotation around the center
  x_translation = 0
  y_translation = 0
  tform = rot ∘ Translation(y_translation, x_translation)
  img2 = warp(img, rot, axes(img))

  return img2
end

"""
    get_local_gradients(video_array, centers, sub_img_size)

Computes the gradients in the subimage, takes the mean of sum of absolute
values of both hotizontal and vertical gradients as a representative of a
subimage.
"""
function get_local_gradients(video_array, centers, sub_img_size)
    @debug "Entering get_local_gradients"
    half_size = ceil(Int,(sub_img_size-1)/2)
    half_range = half_size
    h, w, len = get_video_dimension(video_array)
    extracted_pixels = zeros(sub_img_size, sub_img_size, len)

    @debug "starting frame procesing"
    for frame = 1:len
        img = video_array[frame]
        img_grad = imgradients(img, KernelFactors.ando3, "replicate")
        img_grad_abs = map(abs, img_grad[1]) + map(abs, img_grad[2])
        for index_x = 1:size(centers,2)
            c_x = centers[2, index_x]
            for index_y = 1:size(centers,2)
                c_y = centers[1, index_y]
                sub_img  = img_grad_abs[(c_x-half_size):(c_x+half_size),
                                        (c_y-half_size):(c_y+half_size)]

                extracted_pixels[index_x, index_y, frame] =mean(sub_img)
            end
        end
        # @debug "Next frame" frame
    end
    return extracted_pixels
end

"""
    get_img_gradient(img)

Computes the gradients in the `img`.
"""
function get_img_gradient(img)
    @debug "Entering get_local_gradients"

    img_grad = imgradients(img, KernelFactors.ando3, "replicate")

    grad_1 = img_grad[1] .+ abs(findmin(img_grad[1])[1])
    grad_1 ./= findmax(grad_1)[1]

    grad_2 = img_grad[2] .+ abs(findmin(img_grad[2])[1])
    grad_2 ./= findmax(grad_2)[1]

    grad_sum = grad_1 + grad_2
    grad_sum .+= abs(findmin(grad_sum)[1])
    grad_sum ./= findmax(grad_sum)[1]

    return grad_sum
end

"""
    get_local_img_correlations(img, centers, sub_img_size, shift;
                                                        with_gradient=false)

Computes the correlation between the subimages and subimages shifted by values
from range -`shift`:`shift` and returns array of size
length(`centers`) x length(`centers`).

Each of the subimage is center around values stored in  @centers
"""
function get_local_img_correlations(img, centers, sub_img_size::Int;
                                                    with_gradient=false)
# TODO BUG- function is not workig for even numbers
    half_size = ceil(Int,(sub_img_size-1)/2)
    half_range = half_size#
    h, w = size(img)
    extracted_pixels = zeros(sub_img_size, sub_img_size)
    local_correlation = zeros(size(centers,1))

    if with_gradient
        img = get_img_gradient(img)
    end

    position = 1;
    for index = centers
        c_x = index[1]
        c_y = index[2]
        c_x_range = (c_x-half_range):(c_x+half_range)
        c_y_range = (c_y-half_range):(c_y+half_range)
        subimage = img[c_x_range,c_y_range]
        center = img[c_x_range, c_y_range]

        corelation = center .* subimage
        corelation = sum(corelation)
        local_correlation[position] += corelation

        local_correlation[position] /= 256*(sub_img_size^2)^2
        position += 1;
    end

    return local_correlation
end

"""
    get_local_img_correlations(img,centers, masks)

Takes `img` and computes crosscorrelation with set of `masks` around the
`centers`. Crosscorrelation is computed as convolution of the mask and the area
around coordinates stored in `centres`.
"""
function get_local_img_correlations(img, centers, masks::Vector; with_gradient=false)
    masks_num = length(masks)
    sub_img_size = size(masks[1],1)
	# half_size = ceil(Int,(sub_img_size-1)/2)
    half_size = (sub_img_size)÷2
    half_range = half_size
    h, w = size(img)
    local_correlation = zeros(masks_num, size(centers,1) )
    index = centers[1]
    masks_num = length(masks)

    if with_gradient
        img = get_img_gradient(img)
    end

    # position = 1;
    # for index = centers
    for pos = 1:length(centers)
		# global position
		index = centers[pos]

        c_x = index[1]
        c_y = index[2]
        c_x_range = (c_x-half_range):(c_x+half_range)
        c_y_range = (c_y-half_range):(c_y+half_range)

        center = img[c_x_range, c_y_range]
        # mask_pos = 1
        # for mask in masks
		for mask_pos = 1:length(masks)
			mask = masks[mask_pos]
            corelation = center .* mask

            corelation = sum(corelation)
            local_correlation[mask_pos, pos] += corelation
            local_correlation[mask_pos, pos] /= (sub_img_size^2)
            # local_correlation[position, mask_pos ] =  sum(imfilter(center, mask))/(sub_img_size^2)
            # mask_pos +=1
        end

        # position += 1;
    end

    return local_correlation
end

"""
    extract_pixels_from_img(img, indicies_set, video_dim_tuple)

Takes every frame from @video_array and extracts pixels which indicies are in
@indicies_set, thus creating video with only chosen indicies.
"""
function extract_pixels_from_img(img, indicies_set, video_dim_tuple)
   rows = size(indicies_set,2)
   columns = size(indicies_set,2)
   video_length = video_dim_tuple[3]

   extracted_pixels = zeros(rows, columns, video_length)
   extracted_pixels[:,:,] =
                img[indicies_set[1,:],indicies_set[2,:]]

   return extracted_pixels
end

"""
Returns evenly distributed centers of size `image_size`
"""
function get_local_img_centers(points_per_dim, img_size, shift=0, sub_img_size=0 )
    # TODO Applied teproray solution here, so it works only for local gradients
    # start = 0
    # (points_per_dim>shift) ? start_ind = ceil(Int, points_per_dim/2)+ shift :
    #                         start=shift
    start_ind = ceil(Int, sub_img_size/2)
    min_va,  = findmin(img_size)
    last_ind = min_va - start_ind

    set = broadcast(floor, Int, range(start_ind, step=sub_img_size,  stop=last_ind))
    num_indexes = size(set,1)

    centers = Any[]
    for row = 1:num_indexes
        for column = 1:num_indexes
            push!(centers, CartesianIndex(set[row], set[column]))
        end
    end
    return centers
end

"""
    get_img_local_centers(img_size, sub_img_size=10)

Tiles the image of size @img_size into square subimages of size @sub_img_size
and returns vector with CartesianIndex coordinates of the subimages centre in
original image.

By default, takes smaller value from @img_size and then divides it by
@sub_img_size. Resulting value will be the number of returned subimages per
dimension. If @use_square is set to false, then evry dimension is treated
separately, resulting in rectangular grid.

It is possible to set overlap of the tiles with @overlap parameter. By default
it is set to zero, but can be any pixel value smaller that @sub_img_size. If
@overlap is set to value in range (0,1], a fraction of @sub_img_size is used.
"""
function get_img_local_centers(img_size, sub_img_size=10; use_square=true,
                                overlap=0)

    @assert sub_img_size <= findmin(img_size)[1] "@sub_img_size is bigger than image!"
    @assert sub_img_size > 0 "sub_img_size must be positive number"
    @assert overlap<=sub_img_size "The overlap is biger than subimage size!"
    @assert overlap >= 0 "overlap must be positive"

    centers = CartesianIndex[]

    start_ind = ceil(Int, sub_img_size/2)
    if 2*start_ind == sub_img_size
        start_ind +=1
    end

    if overlap>0 && overlap<1
        overlap = floor(Int, sub_img_size*overlap)
    end

    if use_square
        size_v = findmin(img_size)[1]
        size_h = findmin(img_size)[1]
    else
        size_v = img_size[1]
        size_h = img_size[2]
    end

    last_ind_v = size_v - start_ind # TODO check if it is starting at 1st row, not second
    last_ind_h = size_h - start_ind

    val_range_v = floor.(Int, range(start_ind, step=sub_img_size-overlap,  stop=last_ind_v))
    val_range_h = floor.(Int, range(start_ind, step=sub_img_size-overlap,  stop=last_ind_h))

	if isempty(val_range_v) && size_v <= sub_img_size
		val_range_v = [start_ind]
	end

	if isempty(val_range_h) && size_h <= sub_img_size
		val_range_h = [start_ind]
	end

    num_indexes_v = size(val_range_v,1)
    num_indexes_h = size(val_range_h,1)

    for row = 1:num_indexes_v, column = 1:num_indexes_h
        push!(centers, CartesianIndex(val_range_v[row], val_range_h[column]))
    end
    return centers
end

"""
    vectorize_img(video)

Rearrenges the video so that set of n frames (2D matrix varying in
time) the set of vectors is returned, in which each row is a pixel, and each
column is the value of the pixel in n-th frame.
"""
function vectorize_img(img)
    rows, columns = size(img)
    num_of_elements = rows*columns

    vectorized_img = zeros(num_of_elements)

    index = 1;
    for row=1:rows
        for column=1:columns
            vectorized_img[index] = img[row, column];
            index = index+1;
        end
    end

    return vectorized_img
end

"""
    get_video_mask(points_per_dim, video_dimensions; distribution="uniform", sorted=true, x=1, y=1)

Returns matrix of size @points_per_dim x 2 in which indicies of video frame are
stored. The indicies are chosen based one the @distribution argument. One option
is uniform distribution, the second is random distribution.

Uniform distribution: distance between the points in given dimension is the
 even, but vertical distance may be different from horizontal distance between points. This depends on the size of a frame in a image.

Random distribution: the distance between the points is not constant, because
the points are chosen randomly in the ranges 1:horizontal size of frame,
1:vertical size of frame. The returned values may be sorted in ascending order,
if @sorted=true.
"""
function get_video_mask(points_per_dim, video_dimensions;
                            distribution="uniform", sorted=true, patch_params)
    video_height, video_width,  = video_dimensions
    x=patch_params["x"]
    y=patch_params["y"]
    spread=patch_params["spread"]

    if x == 1
        x = floor(Int,video_width/2)
        @warn "Given x is to close to the border. Seeting the value to " x
    elseif x < ceil(Int,points_per_dim/2)
        x = ceil(Int,points_per_dim/2)
        @warn "Given x is to close to the border. Seeting the value to " x
    elseif x > video_width-ceil(Int,points_per_dim/2)
        x = video_width - ceil(Int,points_per_dim/2)
        @warn "Given x is to close to the border. Seeting the value to " x
    end

    if y == 1
        y = floor(Int,video_height/2)
        @warn "Given y is to close to the border. Seeting the value to " y
    elseif y < ceil(Int,points_per_dim/2)
        y = ceil(Int,points_per_dim/2)
        @warn "Given y is to close to the border. Seeting the value to " y
    elseif y > video_height-ceil(Int,points_per_dim/2)
        y = video_height - ceil(Int,points_per_dim/2)
        @warn "Given y is to close to the border. Seeting the value to " y
    end

    if spread*points_per_dim+x > video_width || spread*points_per_dim+y > video_height
        @warn "Given patch parameters might result in indicies exceeding frame size."
    end

    if distribution == "uniform"
        columns = points_per_dim
        rows = points_per_dim

        # +1 is used so that the number of points returned is as requested
        row_step = floor(Int,video_height/rows)
        column_step = floor(Int,video_width/columns)

        (video_height/row_step != points_per_dim) ? row_step+=1 : row_step
        (video_width/column_step !=
                                points_per_dim) ? column_step+=1 : video_width

        vertical_indicies = collect(1:row_step:video_height)
        horizontal_indicies = collect(1:column_step:video_width)

        vertical_indicies = reshape(vertical_indicies, (1,points_per_dim))
        horizontal_indicies = reshape(horizontal_indicies, (1,points_per_dim))

        indicies_set = [vertical_indicies; horizontal_indicies]
    elseif distribution == "random"
        vertical_indicies = rand(1:video_height,1, points_per_dim)
        horizontal_indicies = rand(1:video_width,1, points_per_dim)

        if sorted
            vertical_indicies = sort(vertical_indicies[1,:])
            horizontal_indicies = sort(horizontal_indicies[1,:])

            vertical_indicies = reshape(vertical_indicies, (1,points_per_dim))
            horizontal_indicies =
                              reshape(horizontal_indicies, (1,points_per_dim))
        end
        indicies_set = [vertical_indicies; horizontal_indicies]
    elseif distribution == "patch"
        indicies_set = [collect(1:spread:(spread*points_per_dim)).+x collect(1:spread:(spread*points_per_dim)).+y]'
    end

   return indicies_set
end

"""
    get_gabor_mask_set(;filt_size=25, σ=[2], theta_rad=[0], λ=[15], γ=[0.2],
                            psi_rad=[0], re_part=true, im_part=false)

Returns set of gabor filters generated with given parameters. Parameters are described
below. Function uses Kernel.gabor() from ImageFiltering.

# Arguments
- `filt_size=30` : controls the patch in which filter is created, not wavelet itself
- `σ=2` : controls the width of the waves and thus number of cycles per unit
- `theta_rad=0` : is the rotation in radians,
- `λ=15` : controls the number of waves within the window- higher values- less waves
- `γ=0.2` : is the aspect ratio; small values give long filters
- `psi_rad=0` : phase in radians
- `re_part::Bool`: determines if real part of the Gabor filter is returned; real
    part is normalized to be in range [-0.5,0.5]
- `im_part::Bool`: determines if imaginary part of the Gabor filter is returned
    imaginary part is normalized to be in range [-0.5,0.5]

if both `re_part` and `im_part` are true, then absolute value of complex number
    of form `re_part + im_part im` is returned (it is also normalized to range
    [-0.5,0.5]).
"""
function get_gabor_mask_set(;filt_size=25, σ=[2], theta_rad=[0], λ=[15], γ=[0.2],
                            psi_rad=[0], re_part=true, im_part=false,
                            do_norm=true, do_negative=true)

    kernels = Any[]
    for sigma = σ
        for angle1 = theta_rad
            θ = angle1; #pi*(angle1/180)
            for lambda in λ
                for gamma in γ
                    for angle2 in psi_rad
                        ψ = angle2; #pi*(angle2/180)
                        kernel = Kernel.gabor(filt_size, filt_size,
                                        sigma,
                                        θ,
                                        lambda,
                                        gamma,
                                        ψ)
                        if re_part && !im_part
                            # @debug "Doing real part"
                            if do_norm
                                kernel[1] .+= abs(findmin(kernel[1])[1])
                                kernel[1] ./= findmax(kernel[1])[1]
                                # @debug "Doing norm"
                                if do_negative
                                    kernel[1] .-= 0.5
                                end
                            end
                            push!(kernels,Gray.((kernel[1])))

                        elseif im_part && !re_part
                            if do_norm
                                kernel[2] .+= abs(findmin(kernel[2])[1])
                                kernel[2] ./= findmax(kernel[2])[1]
                                if do_negative
                                     kernel[2] .-= 0.5;
                                end
                            end
                            push!(kernels,Gray.((kernel[2])))

                        else
                            @debug "Using abs(re(A)+im(A))"
                            result = abs.(kernel[1] + kernel[2]im);
                            if do_norm
                                result .+= abs(findmin(result)[1])
                                result ./= findmax(result)[1]
                            end
                            push!(kernels,Gray.())
                        end
                    end # angle2
                end # gamma
            end # lambda
        end # angle1
    end # sigmas
    return kernels
end

"""
    rearrange_filters_arr(im_filter; showing_number=-1)

Creates image with elements stored in `im_filters`. `showing_number` determines
how many of the element from `im_fiter` are displayed.

`im_filters` is an array with elements of type Matrix{Gray}.
"""
function rearrange_filters_arr(im_filter; showing_number=-1, columns=-1)
    mask_size = size(im_filter[1],1)
    im_filter_num = length(im_filter)
    if showing_number == -1 || showing_number > im_filter_num
        max_indeks = im_filter_num
    else
        max_indeks = showing_number
    end

    if columns == -1
        columns = Int(ceil(sqrt(im_filter_num)))
    end
    rows= Int(ceil(im_filter_num/columns))

    all_filters = zeros(Gray, rows*mask_size, columns*mask_size)

    mask_index = 1
    for row in 1:rows
        start_row = (row-1)*mask_size+1
        row_range = start_row:(start_row+mask_size-1)

        for col = 1:columns
            start_col = (col-1)*mask_size+1
            col_range = start_col:(start_col+mask_size-1)
            if mask_index > max_indeks
                break
            else
                all_filters[row_range, col_range] = im_filter[mask_index]
                mask_index += 1
            end
        end
        if mask_index > max_indeks
            break
        end
    end
    return all_filters
end

# TODO remove img size from arguments
function get_local_correlations(method::String, img, img_size, sub_img_size;
                                                        masks = 0,
                                                        points_per_dim=1,
                                                        shift=0,
                                                        with_grad = true,
                                                        overlap = 0,
                                                        use_square=true)
    if method == "correlation"
        @debug "local correlation"
        centers = get_local_img_centers(points_per_dim, img_size, shift,
                                                            sub_img_size)
        extracted_pixels_matrix = get_local_img_correlations(img, centers,
                                                            sub_img_size, shift)

    elseif  method == "gradient_gabor"
        @info "local gradient gabor comparison"
        centers = get_img_local_centers(img_size, sub_img_size)
        local_correlations = get_local_img_correlations(img, centers, masks;
                                                    with_gradient = with_grad)

    elseif  method == "gabor"
        @debug "local gabor comparison"
        centers = get_img_local_centers(img_size, sub_img_size; overlap = overlap, use_square=use_square)
        local_correlations = get_local_img_correlations(img, centers, masks )

    elseif  method == "gradient"
        @debug "local gradient analysis"
        centers = get_local_img_centers(points_per_dim, img_size, shift,
                                                       sub_img_size)
        local_correlations = get_local_img_correlations(img, centers, sub_img_size;
                                                       with_gradient=with_grad)
    else
        indicies_set = get_video_mask(points_per_dim, img_size,
                            distribution="uniform", patch_params=patch_params)
        local_correlations = extract_pixels_from_img(img, indicies_set,
                                                                    img_size)
    end

    return local_correlations
end


text:
# simulate one generation for given set of parameters

ENV["GKS_ENCODING"] = "utf-8"
using LinearAlgebra, Statistics,  Distributions, Random,  Plots, JLD2
include("sim_functions.jl")
include("show_functions.jl")
include("setup_inital.jl")

#sim_set["shock_type"] = 1  # change to simulate actual shock different from shock used for expectations formation
include("setup_functions.jl")

# reset accouting -- account para not needed
account_econ, account_para = q_set_accounting(sim_set)
Random.seed!(sim_set["r_seed"])
gen_utility, state, account_econ = q_one_generation(state, s_para, econ_set, decide_set, sim_set, account_econ )

# show results
q_show_sim_rep(account_econ)


text:
import Pkg
Pkg.add("POMDPs")
import POMDPs
POMDPs.add_registry()

Pkg.add("Random")
Pkg.add("StatsBase")
Pkg.add("Distributions")
Pkg.add("LinearAlgebra")
Pkg.add("SymPy")

Pkg.add("POMDPModelTools")
Pkg.add("POMDPPolicies")
Pkg.add("POMDPSimulators")
Pkg.add("BeliefUpdaters")
Pkg.add(Pkg.PackageSpec(url = "https://github.com/mexsser/DiscreteValueIteration.jl"))
Pkg.add(Pkg.PackageSpec(url = "https://github.com/JuliaPOMDP/QMDP.jl"))

Pkg.add(Pkg.PackageSpec(url = "https://github.com/mexsser/Geo.jl"))
Pkg.add("Plots")
Pkg.add("Distributed")
Pkg.add("Measures")


text:
using TensorToolbox, Test, LinearAlgebra

println("\n\n****Testing dimtree.jl")

println("\n...Testing creating a dimtree.")
T1=dimtree([8,9,5,6,7],[1,2,3,4])
T2=dimtree([8,9,5,6,7])
@test T1==T2
T3=dimtree([5,6,7,8,9],[1,2,3,4])
T4=dimtree(collect(5:9))
@test T3==T4

println("...Testing function non (number-of-nodes).")
@test non(T1)==non(T3)

println("\n...Created two dimtrees for tensor of order 5, with ",non(T1)," nodes.")

println("\n...Testing functions:")
println(" - height")
h=height(T1)
@test h==4
println(" - children")
c=children(T1,2)
@test c==[4,5]
c=children(T3,3)
@test c==[6,7]
println(" - sibling")
@test sibling(T1,4)==5
@test sibling(T2,5)==4
println(" - parent")
p=parent(T1,4)
@test p==2
p=parent(T3,6)
@test p==3
println(" - is_leaf")
@test is_leaf(T1,5) == true
println(" - is_left")
@test is_left(T1,6) == true
println(" - is_right")
@test is_right(T1,5) == true
println(" - count_leaves")
l=count_leaves(T1,2)
@test l==3
l=count_leaves(T3,3)
@test l==2
println(" - dims")
@test dims(T1,4)==[1,2]
@test dims(T3,3)==[2,3]
println(" - lvl")
@test lvl(T1,9)==3
@test lvl(T2,5)==2
println(" - nodes_on_lvl")
@test nodes_on_lvl(T1,2)==[4,5,6,7]
@test nodes_on_lvl(T3,3)==[8,9]
println(" - node2ind")
@test node2ind(T1,6)==4
@test node2ind(T3,[3,5])==[3,1]
println(" - subnodes")
@test subnodes(T3,4)==[4,8,9]
println(" - left_child_length")
@test left_child_length(T1,2)==2


text:
# ==============================================================================
# Max-Plus Productive (flowshops).
# Note: we use the Julia SimpleWeightedGraphs package for converting the flowshop
# to a grap. We use the Julia TikzPictures package instead of GraphPlot.jl for
# displaying them.
# ==============================================================================

using SimpleWeightedGraphs, TikzPictures

# ==============================================================================
# Save a max-plus matrice to a directed graph.

mpgraph(A::ArrMP) = SimpleWeightedDiGraph(A)
mpgraph(A::SpaMP) = SimpleWeightedDiGraph(A)

# ==============================================================================
# shift

# ==============================================================================
# From the matrix of time of tasks (machines, coins) build a max+ linear system

function flowshop(E::ArrMP{T}) where T
    nmach, npiece = size(E)

    A = mpzeros(T, nmach * npiece, nmach * npiece)
    B = mpzeros(T, nmach * npiece, nmach + npiece)
    C = mpzeros(T, nmach + npiece, nmach * npiece)
    D = mpzeros(T, nmach * npiece, nmach * npiece)
    X0 = mpzeros(T, nmach * npiece, 1)

    l = zeros(Int64, 1, nmach)
    d = zeros(Int64, 1, npiece)
    ij = 0
    for i in 1:npiece
        for j in 1:nmach
            ij = i + (j-1) * npiece
            if iszero(E[j, i])
                !iszero(l[j]) && (l[j] += 1)
                !iszero(d[i]) && (d[i] += 1)
            else
                !iszero(l[j]) && (D[ij, ij - l[j]] = E[j, i - l[j]])
                !iszero(d[i]) && (D[ij, ij - d[i] * npiece] = E[j - d[i],i])
                iszero(d[i])  && (B[ij, i] = mp1);  d[i] = 1
                iszero(l[j])  && (B[ij, j + npiece] = mp1); l[j] = 1
            end
        end
        C[i, ij - (d[i] - 1) * npiece] = E[nmach - d[i] + 1, i]
    end
    for j in 1:nmach
        C[j + npiece, j * npiece - l[j] + 1] = E[j, npiece - l[j] + 1]
    end

    mpsyslin(A, B, C, D, X0)
end

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

#function flowshop_graph(E::arrMP, m::arrMP, p::arrMP)
#    (g, T, N)
#end

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

#function flowshop_simu(E::arrMP, m::arrMP, p::arrMP, u::arrMP)
#    mpsimul()
#end

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

#function critical_graph

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

# SimpleWeightedDiGraph
function LaTeXG(PT::ArrMP, O::Float64)
    # Define colors and shapes for arrows.
    settings="\\tikzset{
              darrow/.style={->,draw=black,line width=1pt},
              barrow/.style={->,draw=blue,line width=1pt},
              yarrow/.style={->,draw=yellow,line width=1pt}}"

    # Define colors for places
    opts="entree/.style={circle,draw=blue,line width=1pt},
          sortie/.style={circle,draw=green,fill=green,line width=1pt},
          vide/.style={circle,draw=black,line width=1pt},
          pleine/.style={circle,draw=green,line width=1pt}"

    nm, np = size(PT)

    code = "%% Piece inputs -- Piece outputs\n"
    for i in 1:np
        code *= "\\node[entree] (pi$(i)) at ($(i),0) {}; \\node[sortie] (po$(i)) at ($(i),$(nm+1+O)) {};\n"
    end

    code *= "%% Machine inputs -- Machine outputs\n"
    for i in 1:nm
        code *= "\\node[entree] (mi$(i)) at (0,$(nm-i+1)) {}; \\node[sortie] (mo$(i)) at ($(np+1+O),$(nm-i+1)) {};\n"
    end

    code *= "%% Places\n"
    for j in 1:np
        code *= "% Row $(j)\n"
        for i in 1:nm
            if !iszero(PT[i,j])
                code *= "\\node[vide] (p$(i)$(j)) at ($(j+O),$(i+O)) {$(i)$(j)};\n"
            end
        end
        for i in 1:nm
            code *= "\\draw[yarrow] (mo$(i)) -- (mi$(i));\n"
        end
        code *= "\\draw[yarrow] (po$(j)) -- (pi$(j));\n"
    end

    code *= "%% Transitions\n"

    TikzPicture(code, options=opts, preamble=settings)
end


text:
using Canvas
using FactCheck

include("macros.jl")


text:
########## Useful functions ##########

suite(a, b) = a <= b ? (a : b) : (a : -1 : b)

########## Point ##########

struct Point
    x::Int
    y::Int
end

x(p::Point) = p.x
y(p::Point) = p.y

Point((x, y)::Tuple{Int, Int}) = Point(x, y)
Point(entry::String) = Point(Tuple(parse.(Int, String.(split(entry, ",")))))

########## Line ##########

struct Line
    points::Vector{Point}
    is_straight::Bool
end

points(l::Line) = l.points
is_straight(l::Line) = l.is_straight

function Line(entry::String, include_diagonals::Bool)
    a, b = Point.(String.(split(entry, " -> ")))
    Line(a, b, include_diagonals)
end

function Line(a::Point, b::Point, include_diagonals)

    is_straight = x(a) == x(b) || y(a) == y(b)
    if include_diagonals
        is_straight = is_straight || (abs(x(a) - x(b)) == abs(y(a) - y(b)))
    end

    is_straight ? points = Point.(suite(x(a), x(b)), suite(y(a), y(b))) : points = [a, b]

    Line(points, is_straight)
end

########## generate_grid and danger_score ##########

function generate_grid(lines::Vector{Line})
    grid = Dict{Point, Int}()
    for line in lines
        if is_straight(line)
            for point in points(line)
                get(grid, point, 0) > 0 ? grid[point] += 1 : grid[point] = 1
            end
        end
    end
    grid
end

danger_score(grid::Dict{Point, Int}) = sum(values(grid) .> 1)

########## Script ##########

function main()

    # open file and store content into an input string
    input = open("inputs//day5.txt", "r") do io
        read(io, String)
    end

    # split input into a data vector
    data = String.(split(input, "\r\n"))

    # part1
    grid1 = generate_grid(Line.(data, false))
    println("part 1: ", danger_score(grid1))

    # part2
    grid2 = generate_grid(Line.(data, true))
    println("part 2: ", danger_score(grid2))

end

main()


text:
#=
Implementation of the 2d unwrapping algorithm described in:

- Herráez, M. A., Burton, D. R., Lalor, M. J., & Gdeisat,
  M. A. (2002). Fast two-dimensional phase-unwrapping algorithm based
  on sorting by reliability following a noncontinuous path. Applied
  Optics, 41(35), 7437. http://dx.doi.org/10.1364/ao.41.007437

=#

using ProgressMeter

export unwrap2d

γ2d(v) = v - trunc(v/2π)*2π
γ2d(a,b) = γ2d(b-a)

# Returns integer number of 2π necessary to bring a to b's branch.
function γi(a,b)
    d = b-a
    if abs(d) > 0.5
        round(Int,d)
    else
        0
    end
end

function diff2(φ)
    D² = zeros(φ)

    m,n = size(φ)

    for i = 1:m
        D²[i,1] = D²[i,end] = Inf
    end
    for j = 1:n
        D²[1,j] = D²[end,j] = Inf
    end

    for i = 2:m-1
        @simd for j = 2:n-1
            @inbounds D²[i,j]  = (γ2d(φ[i-1,j  ], φ[i,j]) - γ2d(φ[i,j], φ[i+1,j  ]))^2 # H
            @inbounds D²[i,j] += (γ2d(φ[i  ,j-1], φ[i,j]) - γ2d(φ[i,j], φ[i  ,j+1]))^2 # V
            @inbounds D²[i,j] += (γ2d(φ[i-1,j-1], φ[i,j]) - γ2d(φ[i,j], φ[i+1,j+1]))^2 # D₁
            @inbounds D²[i,j] += (γ2d(φ[i-1,j+1], φ[i,j]) - γ2d(φ[i,j], φ[i+1,j-1]))^2 # D₂
            @inbounds D²[i,j]  = √D²[i,j]
        end
    end
    D²
end

type Edge
    i
    j
    k
    r
end

function find_edges(R²)
    m,n = size(R²)
    # Find all edges
    V = R²[1:end-1,:] + R²[2:end,:]
    H = R²[:,1:end-1] + R²[:,2:end]
    # Construct edge vector
    edges = Any[nothing,nothing]
    for (t,d,M) in ((1,:v,V), (2,:h, H))
        edges[t] = vec(Edge[Edge(i,j,d,M[i,j]) for i = 1:size(M,1), j = 1:size(M,2)])
    end
    edges = [edges...;]
    # Sort edges according to reliability
    sort!(edges, by = e -> e.r, rev=true)
    # for (i,j) in ((1,1),(1,n),(m,1),(m-1,n-1))
    #     filter!(e -> e.i != i && e.j != j, edges)
    # end
    edges
end

function coords(e::Edge)
    if e.k == :v
        [e.i, e.i+1], e.j
    else
        e.i, [e.j, e.j+1]
    end
end

function coords(c::Tuple)
    if length(c[1]) == 1
        [[c[1];c[1]] c[2]]
    else
        [c[1] [c[2];c[2]]]
    end
end
function coords(c::Tuple,i::Integer)
    if length(c[1]) == 1
        c[1],c[2][i]
    else
        c[1][i],c[2]
    end
end

function form_group!(c, φ, groups, group_map, shifts)
    g_id = length(groups) + 1
    group_map[c...] = g_id
    cc = coords(c)
    gg = γi(φ[c...]...)
    shifts[cc[2,:]] = γi(φ[c...]...)
    append!(groups, Matrix{Int}[cc])
end

function add_to_group!(c, φ, groups, group_map, shifts)
    g = group_map[c...]
    ds = γi(φ[c...]...)
    cc = coords(c)
    if g[1] == 0
        group_map[coords(c,1)...] = g[2]
        groups[g[2]] = [groups[g[2]]; cc[1,:]']
        shifts[cc[1,:]] = shifts[cc[2,:]] + ds
    else
        group_map[coords(c,2)...] = g[1]
        groups[g[1]] = [groups[g[1]]; cc[2,:]']
        shifts[cc[2,:]] = shifts[cc[1,:]] - ds
    end
end

function join_groups!(c, φ, groups, group_map, shifts)
    g = group_map[c...]
    ds = γi(φ[c...]...)
    g1,g2 = groups[g[1]], groups[g[2]]
    l1,l2 = size(g1,1), size(g2,1)
    if l1 > l2
        ds = shifts[g1[1,:]...] - ds
        for t = 1:l2
            group_map[g2[t,:]...] = g[1]
            shifts[g2[t,:]...] += ds
        end
        groups[g[1]] = [groups[g[1]]; g2]
        groups[g[2]] = Matrix{Int}()
    else
        ds = shifts[g1[2,:]...] + ds
        for t = 1:l1
            group_map[g1[t,:]...] = g[2]
            shifts[g1[t,:]...] += ds
        end
        groups[g[2]] = [groups[g[2]]; g1]
        groups[g[1]] = Matrix{Int}(0,0)
    end
end

function unwrap2d(φ)
    φ = copy(φ)
    D² = diff2(φ)
    R² = 1./D²
    φ　/= 2π

    edges = find_edges(R²)

    groups = Matrix{Int}[]
    group_map = zeros(Int, size(φ))

    shifts = zeros(Int, size(φ))

    mg = zeros(length(edges))
    lg = zeros(length(edges))

    @showprogress "Unwrapping " for t = 1:length(edges)
        e = edges[t]
        c = coords(e)
        g = group_map[c...]

        if g[1] == g[2] == 0 # Neither point in any group
            form_group!(c, φ, groups, group_map, shifts)
        elseif g[1] == 0 || g[2] == 0 # One point in a group
            add_to_group!(c, φ, groups, group_map, shifts)
        elseif g[1] != g[2] # Both points in group
            join_groups!(c, φ, groups, group_map, shifts)
        end
        mg[t] = maximum(map(length, groups))
        lg[t] = length(groups)
    end
    2π*(φ+shifts)
end

function unwrap2d_debug(φ)
    φ = copy(φ)
    println("Calculating 2nd difference")
    @time D² = diff2(φ)
    @time R² = 1./D²
    φ　/= 2π

    println("Constructing edges")
    edges = find_edges(R²)

    println("Building groups")
    groups = Matrix{Int}[]
    group_map = zeros(Int, size(φ))

    shifts = zeros(Int, size(φ))

    mg = zeros(length(edges))
    lg = zeros(length(edges))

    @time for t = 1:length(edges)
        e = edges[t]
        c = coords(e)
        g = group_map[c...]

        if g[1] == g[2] == 0 # Neither point in any group
            form_group!(c, φ, groups, group_map, shifts)
        elseif g[1] == 0 || g[2] == 0 # One point in a group
            add_to_group!(c, φ, groups, group_map, shifts)
        elseif g[1] != g[2] # Both points in group
            join_groups!(c, φ, groups, group_map, shifts)
        end
        mg[t] = maximum(map(length, groups))
        lg[t] = length(groups)
    end
    println()
    2π*(φ+shifts),D²,R²,mg,lg,group_map,shifts
end


text:
##############################################################################
##
## Sorting
##
##############################################################################

#########################################
## Permutation & Ordering types/functions
#########################################
# Sorting in Julia works through Orderings, where each ordering is a type
# which defines a comparison function lt(o::Ord, a, b).

# UserColOrdering: user ordering of a column; this is just a convenience container
#                  which allows a user to specify column specific orderings
#                  with "order(column, rev=true,...)"

struct UserColOrdering{T<:ColumnIndex}
    col::T
    kwargs
end

"""
    order(col::ColumnIndex; kwargs...)

Specify sorting order for a column `col` in a data frame.
`kwargs` can be `lt`, `by`, `rev`, and `order` with values
following the rules defined in [`sort!`](@ref).

See also: [`sort!`](@ref), [`sort`](@ref)

# Examples
```jldoctest
julia> df = DataFrame(x = [-3, -1, 0, 2, 4], y = 1:5)
5×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ -3    │ 1     │
│ 2   │ -1    │ 2     │
│ 3   │ 0     │ 3     │
│ 4   │ 2     │ 4     │
│ 5   │ 4     │ 5     │

julia> sort(df, order(:x, rev=true))
5×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 4     │ 5     │
│ 2   │ 2     │ 4     │
│ 3   │ 0     │ 3     │
│ 4   │ -1    │ 2     │
│ 5   │ -3    │ 1     │

julia> sort(df, order(:x, by=abs))
5×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 0     │ 3     │
│ 2   │ -1    │ 2     │
│ 3   │ 2     │ 4     │
│ 4   │ -3    │ 1     │
│ 5   │ 4     │ 5     │
```
"""
order(col::T; kwargs...) where {T<:ColumnIndex} = UserColOrdering{T}(col, kwargs)

# Allow getting the column even if it is not wrapped in a UserColOrdering
_getcol(o::UserColOrdering) = o.col
_getcol(x) = x

###
# Get an Ordering for a single column
###
function ordering(col_ord::UserColOrdering, lt::Function, by::Function,
                  rev::Bool, order::Ordering)
    for (k,v) in pairs(col_ord.kwargs)
        if     k == :lt;    lt    = v
        elseif k == :by;    by    = v
        elseif k == :rev;   rev   = v
        elseif k == :order; order = v
        else
            error("Unknown keyword argument: ", string(k))
        end
    end

    Order.ord(lt,by,rev,order)
end

ordering(col::ColumnIndex, lt::Function, by::Function, rev::Bool, order::Ordering) =
    Order.ord(lt,by,rev,order)

# DFPerm: defines a permutation on a particular DataFrame, using
#         a single ordering (O<:Ordering) or a list of column orderings
#         (O<:AbstractVector{Ordering}), one per DataFrame column
#
#         If a user only specifies a few columns, the DataFrame
#         contained in the DFPerm only contains those columns, and
#         the permutation induced by this ordering is used to
#         sort the original (presumably larger) DataFrame

struct DFPerm{O<:Union{Ordering, AbstractVector}, T<:Tuple} <: Ordering
    ord::O
    cols::T
end

function DFPerm(ords::AbstractVector{O}, cols::T) where {O<:Ordering, T<:Tuple}
    if length(ords) != length(cols)
        error("DFPerm: number of column orderings does not equal the number of columns")
    end
    DFPerm{typeof(ords), T}(ords, cols)
end

DFPerm(o::Union{Ordering, AbstractVector}, df::AbstractDataFrame) =
    DFPerm(o, ntuple(i -> df[!, i], ncol(df)))

# get ordering function for the i-th column used for ordering
col_ordering(o::DFPerm{O}, i::Int) where {O<:Ordering} = o.ord
col_ordering(o::DFPerm{V}, i::Int) where {V<:AbstractVector} = o.ord[i]

function Sort.lt(o::DFPerm, a, b)
    @inbounds for i in 1:length(o.cols)
        ord = col_ordering(o, i)
        col = o.cols[i]
        va = col[a]
        vb = col[b]
        lt(ord, va, vb) && return true
        lt(ord, vb, va) && return false
    end
    false # a and b are equal
end

###
# Get an Ordering for a DataFrame
###

################
## Case 1: no columns requested, so sort by all columns using requested order
################
## Case 1a: single order
######
ordering(df::AbstractDataFrame, lt::Function, by::Function, rev::Bool, order::Ordering) =
    DFPerm(Order.ord(lt, by, rev, order), df)

######
## Case 1b: lt, by, rev, and order are Arrays
######
function ordering(df::AbstractDataFrame,
                  lt::AbstractVector{S}, by::AbstractVector{T},
                  rev::AbstractVector{Bool},
                  order::AbstractVector) where {S<:Function, T<:Function}
    if !(length(lt) == length(by) == length(rev) == length(order) == size(df,2))
        throw(ArgumentError("Orderings must be specified for all DataFrame columns"))
    end
    DFPerm([Order.ord(_lt, _by, _rev, _order) for
            (_lt, _by, _rev, _order) in zip(lt, by, rev, order)], df)
end

################
## Case 2:  Return a regular permutation when there's only one column
################
## Case 2a: The column is given directly
######
ordering(df::AbstractDataFrame, col::ColumnIndex, lt::Function, by::Function,
         rev::Bool, order::Ordering) =
    Perm(Order.ord(lt, by, rev, order), df[!, col])

######
## Case 2b: The column is given as a UserColOrdering
######
ordering(df::AbstractDataFrame, col_ord::UserColOrdering, lt::Function, by::Function,
         rev::Bool, order::Ordering) =
    Perm(ordering(col_ord, lt, by, rev, order), df[!, col_ord.col])

################
## Case 3:  General case: cols is an iterable of a combination of ColumnIndexes and UserColOrderings
################
## Case 3a: None of lt, by, rev, or order is an Array
######
function ordering(df::AbstractDataFrame, cols::AbstractVector, lt::Function,
                  by::Function, rev::Bool, order::Ordering)

    if length(cols) == 0
        return ordering(df, lt, by, rev, order)
    end

    if length(cols) == 1
        return ordering(df, cols[1], lt, by, rev, order)
    end

    # Collect per-column ordering info

    ords = Ordering[]
    newcols = Int[]

    for col in cols
        push!(ords, ordering(col, lt, by, rev, order))
        push!(newcols, index(df)[(_getcol(col))])
    end

    # Simplify ordering when all orderings are the same
    if all([ords[i] == ords[1] for i = 2:length(ords)])
        return DFPerm(ords[1], df[!, newcols])
    end

    return DFPerm(ords, df[!, newcols])
end

######
# Case 3b: cols, lt, by, rev, and order are all arrays
######
function ordering(df::AbstractDataFrame, cols::AbstractVector, lt::AbstractVector{S},
                  by::AbstractVector{T}, rev::AbstractVector{Bool},
                  order::AbstractVector) where {S<:Function, T<:Function}

    if !(length(lt) == length(by) == length(rev) == length(order))
        throw(ArgumentError("All ordering arguments must be 1 or the same length."))
    end

    if length(cols) == 0
        return ordering(df, lt, by, rev, order)
    end

    if length(lt) != length(cols)
        throw(ArgumentError("All ordering arguments must be 1 or the same length " *
                            "as the number of columns requested."))
    end

    if length(cols) == 1
        return ordering(df, cols[1], lt[1], by[1], rev[1], order[1])
    end

    # Collect per-column ordering info

    ords = Ordering[]
    newcols = Int[]

    for i in 1:length(cols)
        push!(ords, ordering(cols[i], lt[i], by[i], rev[i], order[i]))
        push!(newcols, index(df)[(_getcol(cols[i]))])
    end

    # Simplify ordering when all orderings are the same
    if all([ords[i] == ords[1] for i = 2:length(ords)])
        return DFPerm(ords[1], df[!, newcols])
    end

    return DFPerm(ords, df[!, newcols])
end

######
## At least one of lt, by, rev, or order is an array or tuple, so expand all to arrays
######
function ordering(df::AbstractDataFrame, cols::AbstractVector, lt, by, rev, order)
    to_array(src::AbstractVector, dims) = src
    to_array(src::Tuple, dims) = [src...]
    to_array(src, dims) = fill(src, dims)

    dims = length(cols) > 0 ? length(cols) : size(df,2)
    ordering(df, cols,
             to_array(lt, dims),
             to_array(by, dims),
             to_array(rev, dims),
             to_array(order, dims))
end

# an explicit error is thrown as Tuple was supported in the past
ordering(df::AbstractDataFrame, cols::Tuple, args...) =
    throw(ArgumentError("Passing a tuple $cols of column selectors when sorting data " *
                        "frame is not supported. Pass a vector $([cols...]) instead."))

###########################
# Default sorting algorithm
###########################

# TimSort is fast for data with structure, but only if the DataFrame is large enough
# TODO: 8192 is informed but somewhat arbitrary

Sort.defalg(df::AbstractDataFrame) =
    size(df, 1) < 8192 ? Sort.MergeSort : SortingAlgorithms.TimSort

# For DataFrames, we can choose the algorithm based on the column type and requested ordering
function Sort.defalg(df::AbstractDataFrame, ::Type{T}, o::Ordering) where T<:Real
    # If we're sorting a single numerical column in forward or reverse,
    # RadixSort will generally be the fastest stable sort
    if isbitstype(T) && sizeof(T) <= 8 && (o==Order.Forward || o==Order.Reverse)
        SortingAlgorithms.RadixSort
    else
        Sort.defalg(df)
    end
end

Sort.defalg(df::AbstractDataFrame, ::Type, o::Ordering) = Sort.defalg(df)
Sort.defalg(df::AbstractDataFrame, col::ColumnIndex, o::Ordering) =
    Sort.defalg(df, eltype(df[!, col]), o)
Sort.defalg(df::AbstractDataFrame, col_ord::UserColOrdering, o::Ordering) =
    Sort.defalg(df, col_ord.col, o)
Sort.defalg(df::AbstractDataFrame, cols, o::Ordering) = Sort.defalg(df)
Sort.defalg(df::AbstractDataFrame, o::Ordering; alg=nothing, cols=[]) =
    alg != nothing ? alg : Sort.defalg(df, cols, o)

########################
## Actual sort functions
########################

"""
    issorted(df::AbstractDataFrame, cols;
             lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)

Test whether data frame `df` sorted by column(s) `cols`.

`cols` can be any column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).

If `rev` is `true`, reverse sorting is performed. To enable reverse sorting
only for some columns, pass `order(c, rev=true)` in `cols`, with `c` the
corresponding column index (see example below).
See other methods for a description of other keyword arguments.
"""
function Base.issorted(df::AbstractDataFrame, cols=[];
                       lt=isless, by=identity, rev=false, order=Forward)
    # exclude AbstractVector as in that case cols can contain order(...) clauses
    if cols isa MultiColumnIndex && !(cols isa AbstractVector)
        cols = index(df)[cols]
    end
    if cols isa ColumnIndex
        return issorted(df[!, cols], lt=lt, by=by, rev=rev, order=order)
    elseif length(cols) == 1
        return issorted(df[!, cols[1]], lt=lt, by=by, rev=rev, order=order)
    else
        return issorted(1:nrow(df), ordering(df, cols, lt, by, rev, order))
    end
end

# sort and sortperm functions

for s in [:(Base.sort), :(Base.sortperm)]
    @eval begin
        function $s(df::AbstractDataFrame, cols=[];
                    alg=nothing, lt=isless, by=identity, rev=false, order=Forward)
            if !(isa(by, Function) || eltype(by) <: Function)
                msg = "'by' must be a Function or a vector of Functions. " *
                      " Perhaps you wanted 'cols'."
                throw(ArgumentError(msg))
            end
            # exclude AbstractVector as in that case cols can contain order(...) clauses
            if cols isa MultiColumnIndex && !(cols isa AbstractVector)
                cols = index(df)[cols]
            end
            ord = ordering(df, cols, lt, by, rev, order)
            _alg = Sort.defalg(df, ord; alg=alg, cols=cols)
            return $s(df, _alg, ord)
        end
    end
end

"""
    sort(df::AbstractDataFrame, cols;
         alg::Union{Algorithm, Nothing}=nothing, lt=isless, by=identity,
         rev::Bool=false, order::Ordering=Forward)

Return a copy of data frame `df` sorted by column(s) `cols`.

`cols` can be any column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).

If `alg` is `nothing` (the default), the most appropriate algorithm is
chosen automatically among `TimSort`, `MergeSort` and `RadixSort` depending
on the type of the sorting columns and on the number of rows in `df`.
If `rev` is `true`, reverse sorting is performed. To enable reverse sorting
only for some columns, pass `order(c, rev=true)` in `cols`, with `c` the
corresponding column index (see example below).
See [`sort!`](@ref) for a description of other keyword arguments.

# Examples
```jldoctest
julia> df = DataFrame(x = [3, 1, 2, 1], y = ["b", "c", "a", "b"])
4×2 DataFrame
│ Row │ x     │ y      │
│     │ Int64 │ String │
├─────┼───────┼────────┤
│ 1   │ 3     │ b      │
│ 2   │ 1     │ c      │
│ 3   │ 2     │ a      │
│ 4   │ 1     │ b      │

julia> sort(df, :x)
4×2 DataFrame
│ Row │ x     │ y      │
│     │ Int64 │ String │
├─────┼───────┼────────┤
│ 1   │ 1     │ c      │
│ 2   │ 1     │ b      │
│ 3   │ 2     │ a      │
│ 4   │ 3     │ b      │

julia> sort(df, [:x, :y])
4×2 DataFrame
│ Row │ x     │ y      │
│     │ Int64 │ String │
├─────┼───────┼────────┤
│ 1   │ 1     │ b      │
│ 2   │ 1     │ c      │
│ 3   │ 2     │ a      │
│ 4   │ 3     │ b      │

julia> sort(df, [:x, :y], rev=true)
4×2 DataFrame
│ Row │ x     │ y      │
│     │ Int64 │ String │
├─────┼───────┼────────┤
│ 1   │ 3     │ b      │
│ 2   │ 2     │ a      │
│ 3   │ 1     │ c      │
│ 4   │ 1     │ b      │

julia> sort(df, [:x, order(:y, rev=true)])
4×2 DataFrame
│ Row │ x     │ y      │
│     │ Int64 │ String │
├─────┼───────┼────────┤
│ 1   │ 1     │ c      │
│ 2   │ 1     │ b      │
│ 3   │ 2     │ a      │
│ 4   │ 3     │ b      │
```
"""
sort(::AbstractDataFrame, ::Any)

"""
    sortperm(df::AbstractDataFrame, cols;
             alg::Union{Algorithm, Nothing}=nothing, lt=isless, by=identity,
             rev::Bool=false, order::Ordering=Forward)

Return a permutation vector of row indices of data frame `df` that puts them in
sorted order according to column(s) `cols`.

`cols` can be any column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).

If `alg` is `nothing` (the default), the most appropriate algorithm is
chosen automatically among `TimSort`, `MergeSort` and `RadixSort` depending
on the type of the sorting columns and on the number of rows in `df`.
If `rev` is `true`, reverse sorting is performed. To enable reverse sorting
only for some columns, pass `order(c, rev=true)` in `cols`, with `c` the
corresponding column index (see example below).
See other methods for a description of other keyword arguments.

# Examples
```jldoctest
julia> df = DataFrame(x = [3, 1, 2, 1], y = ["b", "c", "a", "b"])
4×2 DataFrame
│ Row │ x     │ y      │
│     │ Int64 │ String │
├─────┼───────┼────────┤
│ 1   │ 3     │ b      │
│ 2   │ 1     │ c      │
│ 3   │ 2     │ a      │
│ 4   │ 1     │ b      │

julia> sortperm(df, :x)
4-element Array{Int64,1}:
 2
 4
 3
 1

julia> sortperm(df, (:x, :y))
4-element Array{Int64,1}:
 4
 2
 3
 1

julia> sortperm(df, (:x, :y), rev=true)
4-element Array{Int64,1}:
 1
 3
 2
 4

 julia> sortperm(df, (:x, order(:y, rev=true)))
 4-element Array{Int64,1}:
  2
  4
  3
  1
```
"""
sortperm(::AbstractDataFrame, ::Any)

Base.sort(df::AbstractDataFrame, a::Algorithm, o::Ordering) =
    df[sortperm(df, a, o),:]
Base.sortperm(df::AbstractDataFrame, a::Algorithm, o::Union{Perm,DFPerm}) =
    sort!([1:size(df, 1);], a, o)
Base.sortperm(df::AbstractDataFrame, a::Algorithm, o::Ordering) =
    sortperm(df, a, DFPerm(o,df))


text:
#---------------------------------------------------------
# # [Estimate fluxes](@id fluxes)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/4_fluxes.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/4_fluxes.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

# This will take you through the process of extracting flux information from a given transport operator.
# It is split into 3 parts
# 1. [Figure out the stencil of the operator](@ref operatorstencil)
# 2. [Partition the operator according to said stencil](@ref operatorpartition)
# 3. [Estimate the flux of a given 3D tracer field](@ref fluxestimate)

# Let's start telling Julia we will be using AIBECS and Plots.

using AIBECS, Plots

#-----------------------------------------------
# ## [1. Operator stencil](@id operatorstencil)
#-----------------------------------------------

# Let's load the OCCA grid information `grd` and transport matrix $\mathbf{T}$ (the variable `T` here).

grd, T = OCCA.load()

# Multiplying $\mathbf{T}$ with a given tracer concentration vector $\boldsymbol{x}$ gives the flux-divergence of that tracer.
# In other words, $\mathbf{T}$ operates on $\boldsymbol{x}$ by exchanging tracers between model-grid boxes.
# But $\mathbf{T}$ only operates on *neighboring* grid boxes.
# I.e., $\mathbf{T}$ usually does not directly exchange tracers between boxes far away from each other.

# AIBECS provides a function, `stencil`, that returns all the relative cartesian indices (AKA the relative sub-indices) of the neighbors used by a given operator.
# In the case of OCCA, the stencil is thus

st = stencil(T, grd)

# There is also a `label` function to, well, label the stencil elements in a more human-readable way.
# (Note `label` is not exported, so you must specify that it comes from the `AIBECS` package, as below.)

[st AIBECS.label.(st)]

# The stencil can be visualized with the `plotstencil` recipe

plotstencil(st)

#md # !!! note
#md #     The continuous equivalent of the matrix $\mathbf{T}$, i.e., the flux-divergence operator $\mathcal{T}$, is a differential operator, which is local by definition.
#md #     It is the spatial discretization of the differential operators that imposes the use of neighboring boxes to estimate spatial gradients.
#nb # > Note:
#nb # > In fact, the continuous equivalent of the matrix $\mathbf{T}$, i.e., the flux-divergence operator $\mathcal{T}$ is a differential operator, which is local by definition.
#nb # > It is the spatial discretization of the differential operators that imposes the use of neighboring boxes to estimate spatial gradients.

# This works for any circulation, so you can swap `grd` and `T` and try again... Here we plot the same figure for a few different circulations available in AIBECS.

plts = Any[]
for Circulation in [OCIM0, OCIM1, OCIM2, OCCA]
    grd, T = Circulation.load()
    push!(plts, plotstencil(stencil(grd, T), title=string(Circulation)))
end
plot(plts..., layout=(2,2))

# These stencils are useful to clarify the origin and destination of tracer fluxes we are going to explore in this guide..

#-----------------------------------------------
# ## [2. Operator Partition](@id operatorpartition)
#-----------------------------------------------

# For a given operator $\mathbf{T}$, to each point (or neighbor) $k$ in its stencil corresponds a unique transport operator $\mathbf{T}_k$, which effectively performs the transfer of tracers from that neighbor.
# In other words, if the stencil has $n$ elements, $\mathbf{T}$ can be partitionned into a family of $n-1$ matrices $\mathbf{T}_k$ that each operate in a single direction of the stencil.

# Let's take the "West" neighbor of the OCCA stencil

dir = st[findfirst(AIBECS.label.(st) .== "West")]

# The AIBECS provides a function `directional_T` that returns the transport operator corresponding to the direction provided from the stencil element.

T_West = directional_transport(T, grd, dir)

# Notice how `T_West` only has nonzeros on the main diagonal and on one off-diagonal, that connect each box index with the index of the box "West" of it.

#-----------------------------------------------
# ## [3. Flux Estimate](@id fluxestimate)
#-----------------------------------------------

# We can look at the depth-integrated mean flux of a fictitious tracer $\boldsymbol{x}$ of concentration 1 mol m⁻³ coming from the west by computing and plotting the vertical integral of $-\mathbf{T}\boldsymbol{x}$:

nwet = count(iswet(grd))
x = ones(nwet)
plotverticalintegral(-T_West * x * u"mol/m^3/s" .|> u"μmol/m^3/s", grd, mask=depthvec(grd) .< 100, color=:seaborn_icefire_gradient, clim=1e2 .* (-1,1))

# To get all the matrices for all directions you can use

directional_transports(T, grd)

#md # !!! warning
#md #     This guide is still a work in progress and likely contains errors


text:
using .BlossomV
using SparseArrays
using LightGraphs: add_edge!
using SimpleWeightedGraphs
using MetaGraphs
using ..Utils

export learn_vtree

#############
# Blossom bottom up method
# This method is only available if the user chooses to import BlossomV manually
#############

"Blossom bottom up method, vars are not used"
function blossom_bottom_up(data::DataFrame;α)
    weight = ones(Float64, num_examples(data))
    (_, mi) = mutual_information(data, weight; α)
    vars = Var.(collect(1:num_features(data)))
    info = round.(Int64, 1000001 .+ to_long_mi(mi, -1, -1000000))

    function f(leaf::Vector{<:Vtree})
        variable_sets = collect.(variables.(leaf))

        # even number of nodes, use blossomv alg
        function blossom_bottom_up_even!(variable_sets)::Tuple{Vector{Tuple{Var, Var}}, Int64}
            # 1. calculate pMI
            pMI = set_mutual_information(info, variable_sets)
            pMI = round.(Int64, pMI)

            # 2. solve by blossomv alg
            len = length(variable_sets)
            m = Matching(len)
            for i in 1 : len, j in i + 1 : len
                add_edge(m, i - 1, j - 1, pMI[i, j]) # blossomv index start from 0
            end

            solve(m)
            all_matches = Set{Tuple{Var, Var}}()
            for v in 1 : len
                push!(all_matches, order_asc(v, get_match(m, v - 1) + 1))
            end

            # 3. calculate scores, map index to var
            all_matches = Vector(collect(all_matches))
            score = 0

            for i in 1 : length(all_matches)
                (x, y) = all_matches[i]
                score += pMI[x, y]
            end

            return (all_matches, score)
        end

        # odd number of nodes, try every 2 combinations
        function blossom_bottom_up_odd!(variable_sets)
            # try all len - 1 conditions, find best score(minimun cost)
            (best_matches, best_score) = (nothing, typemax(Int64))
            len = length(variable_sets)
            for index in 1 : len
                indices = [collect(1:index-1);collect(index+1:len)]
                (matches, score) = blossom_bottom_up_even!(variable_sets[indices])
                if score < best_score
                    (best_matches, best_score) = ([[(indices[l], indices[r]) for (l,r) in matches];[index]], score)
                end
            end
            return (best_matches, best_score)
        end

        if length(variable_sets) % 2 == 0
            (matches, score) = blossom_bottom_up_even!(variable_sets)
        else
            (matches, score) = blossom_bottom_up_odd!(variable_sets)
        end

        pairs = []
        for x in matches
            if x isa Tuple
                push!(pairs, (leaf[x[1]], leaf[x[2]]))
            else
                push!(pairs, leaf[x])
            end
        end
        return pairs
    end
    return f
end


text:
function compute(n::Int)::Int
    sieve = zeros(Int, n)
    for i ∈ 2:n
        sieve[i] += i - 1 - i ÷ 2 - i ÷ 3
        j = 2i
        sieve[j:i:n] -= sieve[i]
    end
    return sum(sieve)
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
#=
 Copyright 2013 - 2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015 Oliver Heimlich (oheim@posteo.de)

 Original author: Marco Nehmeier (unit tests in libieeep1788)
 Converted into portable ITL format by Oliver Heimlich with minor corrections.

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
=#
#Language imports

#Test library imports
using Test
#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal_isCommonInterval_test" begin
    @test iscommon(Interval(-27.0, -27.0))
    @test iscommon(Interval(-27.0, 0.0))
    @test iscommon(Interval(0.0, 0.0))
    @test iscommon(Interval(-0.0, -0.0))
    @test iscommon(Interval(-0.0, 0.0))
    @test iscommon(Interval(0.0, -0.0))
    @test iscommon(Interval(5.0, 12.4))
    @test iscommon(Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023))
    @test iscommon(entireinterval(Float64)) == false
    @test iscommon(∅) == false
    @test iscommon(Interval(-Inf, 0.0)) == false
    @test iscommon(Interval(0.0, Inf)) == false
end

@testset "minimal_isCommonInterval_dec_test" begin
    @test iscommon(DecoratedInterval(Interval(-27.0, -27.0), com))
    @test iscommon(DecoratedInterval(Interval(-27.0, 0.0), com))
    @test iscommon(DecoratedInterval(Interval(0.0, 0.0), com))
    @test iscommon(DecoratedInterval(Interval(-0.0, -0.0), com))
    @test iscommon(DecoratedInterval(Interval(-0.0, 0.0), com))
    @test iscommon(DecoratedInterval(Interval(0.0, -0.0), com))
    @test iscommon(DecoratedInterval(Interval(5.0, 12.4), com))
    @test iscommon(DecoratedInterval(Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023), com))
    @test iscommon(DecoratedInterval(Interval(-27.0, -27.0), trv))
    @test iscommon(DecoratedInterval(Interval(-27.0, 0.0), def))
    @test iscommon(DecoratedInterval(Interval(0.0, 0.0), dac))
    @test iscommon(DecoratedInterval(Interval(-0.0, -0.0), trv))
    @test iscommon(DecoratedInterval(Interval(-0.0, 0.0), def))
    @test iscommon(DecoratedInterval(Interval(0.0, -0.0), dac))
    @test iscommon(DecoratedInterval(Interval(5.0, 12.4), def))
    @test iscommon(DecoratedInterval(Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023), trv))
    @test iscommon(DecoratedInterval(entireinterval(Float64), dac)) == false
    @test iscommon(DecoratedInterval(∅, trv)) == false
    @test iscommon(DecoratedInterval(∅, trv)) == false
    @test iscommon(DecoratedInterval(Interval(-Inf, 0.0), trv)) == false
    @test iscommon(DecoratedInterval(Interval(0.0, Inf), def)) == false
end

@testset "minimal_isSingleton_test" begin
    @test isthin(Interval(-27.0, -27.0))
    @test isthin(Interval(-2.0, -2.0))
    @test isthin(Interval(12.0, 12.0))
    @test isthin(Interval(17.1, 17.1))
    @test isthin(Interval(-0.0, -0.0))
    @test isthin(Interval(0.0, 0.0))
    @test isthin(Interval(-0.0, 0.0))
    @test isthin(Interval(0.0, -0.0))
    @test isthin(∅) == false
    @test isthin(entireinterval(Float64)) == false
    @test isthin(Interval(-1.0, 0.0)) == false
    @test isthin(Interval(-1.0, -0.5)) == false
    @test isthin(Interval(1.0, 2.0)) == false
    @test isthin(Interval(-Inf, -0x1.fffffffffffffp1023)) == false
    @test isthin(Interval(-1.0, Inf)) == false
end

@testset "minimal_isSingleton_dec_test" begin
    @test isthin(DecoratedInterval(Interval(-27.0, -27.0), def))
    @test isthin(DecoratedInterval(Interval(-2.0, -2.0), trv))
    @test isthin(DecoratedInterval(Interval(12.0, 12.0), dac))
    @test isthin(DecoratedInterval(Interval(17.1, 17.1), com))
    @test isthin(DecoratedInterval(Interval(-0.0, -0.0), def))
    @test isthin(DecoratedInterval(Interval(0.0, 0.0), com))
    @test isthin(DecoratedInterval(Interval(-0.0, 0.0), def))
    @test isthin(DecoratedInterval(Interval(0.0, -0.0), dac))
    @test isthin(DecoratedInterval(∅, trv)) == false
    @test isthin(DecoratedInterval(∅, trv)) == false
    @test isthin(DecoratedInterval(entireinterval(Float64), def)) == false
    @test isthin(DecoratedInterval(Interval(-1.0, 0.0), dac)) == false
    @test isthin(DecoratedInterval(Interval(-1.0, -0.5), com)) == false
    @test isthin(DecoratedInterval(Interval(1.0, 2.0), def)) == false
    @test isthin(DecoratedInterval(Interval(-Inf, -0x1.fffffffffffffp1023), dac)) == false
    @test isthin(DecoratedInterval(Interval(-1.0, Inf), trv)) == false
end

@testset "minimal_isMember_test" begin
    @test ∈(-27.0, Interval(-27.0, -27.0))
    @test in(-27.0, Interval(-27.0, -27.0))
    @test ∈(-27.0, Interval(-27.0, 0.0))
    @test in(-27.0, Interval(-27.0, 0.0))
    @test ∈(-7.0, Interval(-27.0, 0.0))
    @test in(-7.0, Interval(-27.0, 0.0))
    @test ∈(0.0, Interval(-27.0, 0.0))
    @test in(0.0, Interval(-27.0, 0.0))
    @test ∈(-0.0, Interval(0.0, 0.0))
    @test in(-0.0, Interval(0.0, 0.0))
    @test ∈(0.0, Interval(0.0, 0.0))
    @test in(0.0, Interval(0.0, 0.0))
    @test ∈(0.0, Interval(-0.0, -0.0))
    @test in(0.0, Interval(-0.0, -0.0))
    @test ∈(0.0, Interval(-0.0, 0.0))
    @test in(0.0, Interval(-0.0, 0.0))
    @test ∈(0.0, Interval(0.0, -0.0))
    @test in(0.0, Interval(0.0, -0.0))
    @test ∈(5.0, Interval(5.0, 12.4))
    @test in(5.0, Interval(5.0, 12.4))
    @test ∈(6.3, Interval(5.0, 12.4))
    @test in(6.3, Interval(5.0, 12.4))
    @test ∈(12.4, Interval(5.0, 12.4))
    @test in(12.4, Interval(5.0, 12.4))
    @test ∈(0.0, entireinterval(Float64))
    @test in(0.0, entireinterval(Float64))
    @test ∈(5.0, entireinterval(Float64))
    @test in(5.0, entireinterval(Float64))
    @test ∈(6.3, entireinterval(Float64))
    @test in(6.3, entireinterval(Float64))
    @test ∈(12.4, entireinterval(Float64))
    @test in(12.4, entireinterval(Float64))
    @test ∈(-71.0, Interval(-27.0, 0.0)) == false
    @test in(-71.0, Interval(-27.0, 0.0)) == false
    @test ∈(0.1, Interval(-27.0, 0.0)) == false
    @test in(0.1, Interval(-27.0, 0.0)) == false
    @test ∈(-0.01, Interval(0.0, 0.0)) == false
    @test in(-0.01, Interval(0.0, 0.0)) == false
    @test ∈(0.000001, Interval(0.0, 0.0)) == false
    @test in(0.000001, Interval(0.0, 0.0)) == false
    @test ∈(111110.0, Interval(-0.0, -0.0)) == false
    @test in(111110.0, Interval(-0.0, -0.0)) == false
    @test ∈(4.9, Interval(5.0, 12.4)) == false
    @test in(4.9, Interval(5.0, 12.4)) == false
    @test ∈(-6.3, Interval(5.0, 12.4)) == false
    @test in(-6.3, Interval(5.0, 12.4)) == false
    @test ∈(0.0, ∅) == false
    @test in(0.0, ∅) == false
    @test ∈(-4535.3, ∅) == false
    @test in(-4535.3, ∅) == false
    @test ∈(-Inf, ∅) == false
    @test in(-Inf, ∅) == false
    @test ∈(Inf, ∅) == false
    @test in(Inf, ∅) == false
    @test ∈(-Inf, entireinterval(Float64)) == false
    @test in(-Inf, entireinterval(Float64)) == false
    @test ∈(Inf, entireinterval(Float64)) == false
    @test in(Inf, entireinterval(Float64)) == false
end

@testset "minimal_isMember_dec_test" begin
    @test ∈(-27.0, DecoratedInterval(Interval(-27.0, -27.0), trv))
    @test in(-27.0, DecoratedInterval(Interval(-27.0, -27.0), trv))
    @test ∈(-27.0, DecoratedInterval(Interval(-27.0, 0.0), def))
    @test in(-27.0, DecoratedInterval(Interval(-27.0, 0.0), def))
    @test ∈(-7.0, DecoratedInterval(Interval(-27.0, 0.0), dac))
    @test in(-7.0, DecoratedInterval(Interval(-27.0, 0.0), dac))
    @test ∈(0.0, DecoratedInterval(Interval(-27.0, 0.0), com))
    @test in(0.0, DecoratedInterval(Interval(-27.0, 0.0), com))
    @test ∈(-0.0, DecoratedInterval(Interval(0.0, 0.0), trv))
    @test in(-0.0, DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ∈(0.0, DecoratedInterval(Interval(0.0, 0.0), def))
    @test in(0.0, DecoratedInterval(Interval(0.0, 0.0), def))
    @test ∈(0.0, DecoratedInterval(Interval(-0.0, -0.0), dac))
    @test in(0.0, DecoratedInterval(Interval(-0.0, -0.0), dac))
    @test ∈(0.0, DecoratedInterval(Interval(-0.0, 0.0), com))
    @test in(0.0, DecoratedInterval(Interval(-0.0, 0.0), com))
    @test ∈(0.0, DecoratedInterval(Interval(0.0, -0.0), trv))
    @test in(0.0, DecoratedInterval(Interval(0.0, -0.0), trv))
    @test ∈(5.0, DecoratedInterval(Interval(5.0, 12.4), def))
    @test in(5.0, DecoratedInterval(Interval(5.0, 12.4), def))
    @test ∈(6.3, DecoratedInterval(Interval(5.0, 12.4), dac))
    @test in(6.3, DecoratedInterval(Interval(5.0, 12.4), dac))
    @test ∈(12.4, DecoratedInterval(Interval(5.0, 12.4), com))
    @test in(12.4, DecoratedInterval(Interval(5.0, 12.4), com))
    @test ∈(0.0, DecoratedInterval(entireinterval(Float64), trv))
    @test in(0.0, DecoratedInterval(entireinterval(Float64), trv))
    @test ∈(5.0, DecoratedInterval(entireinterval(Float64), def))
    @test in(5.0, DecoratedInterval(entireinterval(Float64), def))
    @test ∈(6.3, DecoratedInterval(entireinterval(Float64), dac))
    @test in(6.3, DecoratedInterval(entireinterval(Float64), dac))
    @test ∈(12.4, DecoratedInterval(entireinterval(Float64), trv))
    @test in(12.4, DecoratedInterval(entireinterval(Float64), trv))
    @test ∈(-71.0, DecoratedInterval(Interval(-27.0, 0.0), trv)) == false
    @test in(-71.0, DecoratedInterval(Interval(-27.0, 0.0), trv)) == false
    @test ∈(0.1, DecoratedInterval(Interval(-27.0, 0.0), def)) == false
    @test in(0.1, DecoratedInterval(Interval(-27.0, 0.0), def)) == false
    @test ∈(-0.01, DecoratedInterval(Interval(0.0, 0.0), dac)) == false
    @test in(-0.01, DecoratedInterval(Interval(0.0, 0.0), dac)) == false
    @test ∈(0.000001, DecoratedInterval(Interval(0.0, 0.0), com)) == false
    @test in(0.000001, DecoratedInterval(Interval(0.0, 0.0), com)) == false
    @test ∈(111110.0, DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
    @test in(111110.0, DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
    @test ∈(4.9, DecoratedInterval(Interval(5.0, 12.4), def)) == false
    @test in(4.9, DecoratedInterval(Interval(5.0, 12.4), def)) == false
    @test ∈(-6.3, DecoratedInterval(Interval(5.0, 12.4), dac)) == false
    @test in(-6.3, DecoratedInterval(Interval(5.0, 12.4), dac)) == false
    @test ∈(0.0, DecoratedInterval(∅, trv)) == false
    @test in(0.0, DecoratedInterval(∅, trv)) == false
    @test ∈(0.0, DecoratedInterval(∅, trv)) == false
    @test in(0.0, DecoratedInterval(∅, trv)) == false
    @test ∈(-4535.3, DecoratedInterval(∅, trv)) == false
    @test in(-4535.3, DecoratedInterval(∅, trv)) == false
    @test ∈(-4535.3, DecoratedInterval(∅, trv)) == false
    @test in(-4535.3, DecoratedInterval(∅, trv)) == false
    @test ∈(-Inf, DecoratedInterval(∅, trv)) == false
    @test in(-Inf, DecoratedInterval(∅, trv)) == false
    @test ∈(-Inf, DecoratedInterval(∅, trv)) == false
    @test in(-Inf, DecoratedInterval(∅, trv)) == false
    @test ∈(Inf, DecoratedInterval(∅, trv)) == false
    @test in(Inf, DecoratedInterval(∅, trv)) == false
    @test ∈(Inf, DecoratedInterval(∅, trv)) == false
    @test in(Inf, DecoratedInterval(∅, trv)) == false
    @test ∈(-Inf, DecoratedInterval(entireinterval(Float64), trv)) == false
    @test in(-Inf, DecoratedInterval(entireinterval(Float64), trv)) == false
    @test ∈(Inf, DecoratedInterval(entireinterval(Float64), def)) == false
    @test in(Inf, DecoratedInterval(entireinterval(Float64), def)) == false
end
# FactCheck.exitstatus()


text:
subsolver = GLPK.Optimizer

regularizers = [DontRegularize(),
                RegularizedDecomposition(penaltyterm = Linearized()),
                TrustRegion(),
                LevelSet(penaltyterm = InfNorm())]

aggregators = [DontAggregate(),
               PartialAggregate(2),
               Aggregate(),
               DynamicAggregate(2, SelectUniform(2))]

consolidators = [Consolidate(), DontConsolidate()]

penalizations = [Fixed(),
                 Adaptive()]

executions = [Synchronous(), Asynchronous()]

@testset "Structured Solvers" begin
    @info "Running L-shaped tests..."
    @testset "L-shaped: simple problems" begin
        @testset "L-shaped: simple problems" begin
            for (model,scenarios,res,name) in problems
                tol = 1e-5
                sp = instantiate(model,
                                 scenarios,
                                 optimizer = LShaped.Optimizer)
                @test_throws UnloadableStructure optimize!(sp)
                set_silent(sp)
                for execution in executions, regularizer in regularizers, aggregator in aggregators, consolidator in consolidators
                    set_optimizer_attribute(sp, Execution(), execution)
                    set_optimizer_attribute(sp, Regularizer(), regularizer)
                    set_optimizer_attribute(sp, Aggregator(), aggregator)
                    set_optimizer_attribute(sp, Consolidator(), consolidator)
                    @testset "$(optimizer_name(sp)): $name" begin
                        set_optimizer_attribute(sp, MasterOptimizer(), subsolver)
                        set_optimizer_attribute(sp, SubproblemOptimizer(), subsolver)
                        if name == "Infeasible"
                            with_logger(NullLogger()) do
                                set_optimizer_attribute(sp, FeasibilityCuts(), false)
                                optimize!(sp, crash = Crash.EVP())
                                @test termination_status(sp) == MOI.INFEASIBLE
                                set_optimizer_attribute(sp, FeasibilityCuts(), true)
                            end
                        end
                        optimize!(sp, crash = Crash.EVP())
                        @test termination_status(sp) == MOI.OPTIMAL
                        @test isapprox(objective_value(sp), res.VRP, rtol = tol)
                        @test isapprox(optimal_decision(sp), res.x̄, rtol = sqrt(tol))
                    end
                end
            end
        end
    end
    @info "Running progressive-hedging tests..."
    @testset "Progressive-hedging: simple problems" begin
        for (model,scenarios,res,name) in problems
            tol = 1e-2
            sp = instantiate(model,
                             scenarios,
                             optimizer = ProgressiveHedging.Optimizer)
            @test_throws UnloadableStructure optimize!(sp)
            set_silent(sp)
            for execution in executions, penalizer in penalizations
                set_optimizer_attribute(sp, Execution(), execution)
                set_optimizer_attribute(sp, Penalizer(), penalizer)
                set_optimizer_attribute(sp, SubproblemOptimizer(), subsolver)
                set_optimizer_attribute(sp, PrimalTolerance(), 1e-3)
                set_optimizer_attribute(sp, DualTolerance(), 1e-2)
                set_optimizer_attribute(sp, Penaltyterm(), Linearized(num_breakpoints = 1000, spacing = 0.5))
                @testset "$(optimizer_name(sp)): $name" begin
                    optimize!(sp)
                    @test termination_status(sp) == MOI.OPTIMAL
                    @test isapprox(objective_value(sp), res.VRP, rtol = tol)
                    @test isapprox(optimal_decision(sp), res.x̄, rtol = sqrt(tol))
                end
            end
        end
    end
end


text:
@testset "DemoPage" begin
    root = joinpath("assets", "page")

    # default DemoPage behavior
    page = DemoPage(joinpath(root, "default"))
    @test page.title == "default"
    sec1, sec2 = page.sections
    @test sec1.root == joinpath(page.root, "subsection_1")
    @test sec2.root == joinpath(page.root, "subsection_2")

    # title and order
    page = DemoPage(joinpath(root, "title_and_order"))
    @test page.title == "Custom Title"
    sec1, sec2 = page.sections
    @test sec1.root == joinpath(page.root, "subsection_2")
    @test sec2.root == joinpath(page.root, "subsection_1")

    # template
    page = DemoPage(joinpath(root, "template"))
    @test page.title == "Custom Title"

    # template has higher priority
    @suppress_err page = DemoPage(joinpath(root, "suppressed_title"))
    @test page.title == "Custom Title"
end


text:
export Subject, SubjectFactory

import Base: show, similar

##

struct SubjectListener{I}
    schedulerinstance :: I
    actor
end

Base.show(io::IO, ::SubjectListener) = print(io, "SubjectListener()")

mutable struct SubjectProps
    iscompleted :: Bool
    isfailed    :: Bool
    lasterror   :: Any

    SubjectProps() = new(false, false, nothing)
end

Base.show(io::IO, ::SubjectProps) = print(io, "SubjectProps()")

##

"""
    Subject(::Type{D}; scheduler::H = AsapScheduler())

A Subject is a special type of Observable that allows values to be multicasted to many Observers. Subjects are like EventEmitters.
Every Subject is an Observable and an Actor. You can subscribe to a Subject, and you can call `next!` to feed values as well as `error!` and `complete!`.

Note: By convention, every actor subscribed to a Subject observable is not allowed to throw exceptions during `next!`, `error!` and `complete!` calls.
Doing so would lead to undefined behaviour. Use `safe()` operator to bypass this rule.

See also: [`SubjectFactory`](@ref), [`ReplaySubject`](@ref), [`BehaviorSubject`](@ref), [`safe`](@ref)
"""
struct Subject{D, H, I} <: AbstractSubject{D}
    listeners :: Vector{SubjectListener{I}}
    props     :: SubjectProps
    scheduler :: H

    Subject{D, H, I}(scheduler::H) where { D, H <: AbstractScheduler, I } = new(Vector{SubjectListener{I}}(), SubjectProps(), scheduler)
end

function Subject(::Type{D}; scheduler::H = AsapScheduler()) where { D, H <: AbstractScheduler }
    return Subject{D, H, instancetype(D, H)}(scheduler)
end

Base.show(io::IO, ::Subject{D, H}) where { D, H } = print(io, "Subject($D, $H)")

Base.similar(subject::Subject{D, H}) where { D, H } = Subject(D; scheduler = similar(subject.scheduler))

##

iscompleted(subject::Subject) = subject.props.iscompleted
isfailed(subject::Subject)    = subject.props.isfailed
lasterror(subject::Subject)   = subject.props.lasterror

setcompleted!(subject::Subject)      = subject.props.iscompleted = true
setfailed!(subject::Subject)         = subject.props.isfailed = true
setlasterror!(subject::Subject, err) = subject.props.lasterror = err

##

function on_next!(subject::Subject{D, H, I}, data::D) where { D, H, I }
    for listener in copy(subject.listeners)
        scheduled_next!(listener.actor, data, listener.schedulerinstance)
    end
end

function on_error!(subject::Subject, err)
    if !iscompleted(subject) && !isfailed(subject)
        setfailed!(subject)
        setlasterror!(subject, err)
        for listener in copy(subject.listeners)
            scheduled_error!(listener.actor, err, listener.schedulerinstance)
        end
        unsubscribe_listeners!(subject, subject.listeners)
    end
end

function on_complete!(subject::Subject)
    if !iscompleted(subject) && !isfailed(subject)
        setcompleted!(subject)
        for listener in copy(subject.listeners)
            scheduled_complete!(listener.actor, listener.schedulerinstance)
        end
        unsubscribe_listeners!(subject, subject.listeners)
    end
end

function unsubscribe_listeners!(subject::Subject, listeners)
    foreach((listener) -> unsubscribe!(SubjectSubscription(subject, listener)), copy(listeners))
end

##

function on_subscribe!(subject::Subject{D}, actor) where { D }
    if isfailed(subject)
        error!(actor, lasterror(subject))
        return voidTeardown
    elseif iscompleted(subject)
        complete!(actor)
        return voidTeardown
    else
        instance = makeinstance(D, subject.scheduler)
        return scheduled_subscription!(subject, actor, instance)
    end
end

function on_subscribe!(subject::Subject, actor, instance)
    listener = SubjectListener(instance, actor)
    push!(subject.listeners, listener)
    return SubjectSubscription(subject, listener)
end

##

struct SubjectSubscription{ S, L } <: Teardown
    subject  :: S
    listener :: L
end

as_teardown(::Type{ <: SubjectSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::SubjectSubscription)
    filter!((listener) -> listener !== subscription.listener, subscription.subject.listeners)
    return nothing
end

Base.show(io::IO, ::SubjectSubscription) = print(io, "SubjectSubscription()")

##

"""
    SubjectFactory(scheduler::H) where { H <: AbstractScheduler }

A base subject factory that creates an instance of Subject with specified scheduler.

See also: [`AbstractSubjectFactory`](@ref), [`Subject`](@ref)
"""
struct SubjectFactory{ H <: AbstractScheduler } <: AbstractSubjectFactory
    scheduler :: H
end

create_subject(::Type{L}, factory::SubjectFactory) where L = Subject(L, scheduler = similar(factory.scheduler))

Base.show(io::IO, ::SubjectFactory{H}) where H = print(io, "SubjectFactory($H)")


text:
#==============================================================================#
# RJNN: Forward- and back-propagation
#==============================================================================#

#==============================================================================#
# Notes:
# 1. At the time of writing, @view calls still allocate a small amount of memory
#  on the heap (and add a bit of execution time). These calls are necessary when
#  dealing with the final, partial batch of data, which is handled separately
#  from the other batches in the functions below. If views were 'free', we could
#  simply apply them to all of the batches and make the code a bit neater.
#==============================================================================#

# h(x) = 1.71590471*tanh.(0.66666667x)
# dh(h) = 1.14393647*(1 - 0.33963596h^2) # takes h(x) as an argument.
h(x) = tanh.(x)
dh(h) = (1 - h^2)

"Applies a neural network to a given set of observations."
function predict!(A::AbstractMatrix{Float64}, net,
        data::AbstractMatrix{Float64})
    # A is an auxiliary matrix that will be used to store the results; it should
    # have as many rows as the data matrix and as many columns as the network's
    # output layer.
    k, n = batchsize(net), size(data, 1)
    for i = 1:k:n
        j = min(i+k-1, n) # the final batch might not be full.
        feedforward!(net, @view data[i:j,:])
        if j == i+k-1
            A[i:j,:] = net.activities[end]
        else
            A[i:n,:] = @view net.activities[end][1:n-i+1,:]
        end
    end
end

"Computes and stores activities for all of a network's units."
function feedforward!(net, data_batch::AbstractMatrix{Float64})
    k = size(data_batch, 1)
    partial = k < batchsize(net)
    B = partial ? @view(net.activities[1][1:k,:]) : net.activities[1]
    B .= data_batch
    for l = 1:depth(net)
        # B contains the activities for layer l-1 (stored at index l).
        W, b = net.weights[l], net.biases[l]
        A = partial ? @view(net.activities[l+1][1:k,:]) : net.activities[l+1]
        activities!(A, B, W, b)
        B = A
    end
end

"Computes activities for a single layer of units."
function activities!(A, B, W, b)
    # B contains the activities of the previous layer, and W and b the weights
    # and biases of the current layer. The resulting activities h(B*W + b) will
    # be stored in A.
    mul!(A, B, W)
    A .= A .+ b
    A .= h.(A)
end

"Computes the weight gradient of a network with respect to a given data set."
function backprop!(net, data::AbstractMatrix{Float64},
        targets::AbstractMatrix{Float64})
    # In addition to computing the weight gradient, this function computes the
    # gradient with respect to the activations and inputs as well. One could
    # easily alter it to store these derivatives in an auxiliary matrix (similar
    # to the one used in the prediction function).
    k, n = batchsize(net), size(data, 1)
    for G in net.gradients
        fill!(G, 0.)
    end
    for i = 1:k:n
        j = min(i+k-1, n)
        feedforward!(net, @view data[i:j,:]) # activities.
        feedback!(net, @view targets[i:j,:]) # errors and gradient.
    end
end

"Computes error derivatives with respect to the network's activations and uses
 these to update the weight gradient."
function feedback!(net, target_batch::AbstractMatrix{Float64})
    k = size(target_batch, 1)
    partial = k < batchsize(net)
    A, B = net.activities[end], net.errors[end]
    if partial
        A, B = @views A[1:k,:], B[1:k,:]
    end
    outerrors!(B, A, target_batch) # function barrier for type stability.
    for l = depth(net):-1:1
        # B contains the errors for layer l (stored at index l+1).
        W, G = net.weights[l], net.gradients[l]
        A, E = net.activities[l], net.errors[l]
        if partial
            A, E = @views A[1:k,:], E[1:k,:]
        end
        gradient!(G, B, A) # for the weights connecting layers l-1 and l.
        if l > 1
            errors!(E, B, A, W) # for the units in layer l-1.
        else
            inerrors!(E, B, W) # for the input units.
        end
        B = E
    end
end

outerrors!(E, A, T) = E .= (A .- T) .* dh.(A)
inerrors!(E, B, W) = mul!(E, B, transpose(W))

"Computes error derivatives with respect to the activations of a hidden layer."
function errors!(E, B, A, W)
    # B contains the errors of the next layer, A the activities of the current
    # layer, and W the weights between the two layers. The resulting errors will
    # be stored in E.
    mul!(E, B, transpose(W))
    E .*= dh.(A)
end

"Computes error derivatives with respect to the weights between two layers and
 adds these to the current gradient."
function gradient!(G, B, A)
    # B contains the errors of the next layer, and A the activities of the
    # current one. The gradient with respect to the weights between these layers
    # will be added to G (as in G += A'*B).
    mul!(G, transpose(A), B, 1., 1.)
end

max_stars_repo_name: kevdur/rjnn

max_issues_repo_name: kevdur/rjnn

max_forks_repo_name: kevdur/rjnn