text:
# Solution of Hamilton-Jacobi-Bellman Equations
# Michele Fornino

# # Dependencies
# using LinearAlgebra
# using SparseArrays
# using MarkovChains
# using ProgressMeter

# import Base:
# 	show

# Temporary Modules
# using Distributed
# using BenchmarkTools
# using Debugger
# using PyPlot

#---------------------------------------------------------------------
# 						HJB EQUATION TYPE
# --------------------------------------------------------------------

"""
HJBEquation

Type for storing a HJB Equation problem, i.e., a continuous time DP problem.

Can be constructed directly by supplying:

* flow_returns function mapping states and controls to flow returns.
* state_equations function mapping states and controls to drifts of endogenous state variables.
* optimal_controls function returning optimized controls as a function of the states and of ∇V, the gradient of the value function.
* state_constraints function, returning the gradient of the value function for a choice of controls that sets drifts to zero.
* state_grids_endog a tuple with an array of equally spaced gridpoints for each endogenous state variable.
* exog_shock a ContinuousTimeMarkovChain object summarizing the exogenous forcing term.
* discount_rate a.k.a. ρ

See also <https://en.wikipedia.org/wiki/Hamilton%E2%80%93Jacobi%E2%80%93Bellman_equation>.

"""
struct HJBEquation

	# Basic objects
	discount_rate::Real									# Discount factor
	flow_returns::Function										# Returns function
	state_equations::Function 		 							# State Equations
	optimal_controls::Function 									# Optimal choice of controls from FOC(s)
	state_constraints::Function 								# Function to impose state constraints
	exog_infinitesimal_generator::AbstractMatrix{<:Real}	# Jump-Diffusion infinitesimal generator
	state_grids::Tuple{Vararg{AbstractVector{<:Real}}}					# Grids state variables

	# Useful objects
	N::Integer									# Number of endogenou states
	M::Integer									# Number of control variables
	nstates::Integer
	dx::Tuple{Vararg{<:Real}}	# diff(state_grids)
	state_grids_lengths::Tuple{Vararg{<:Integer}}
	state_grids_strides::Tuple{Vararg{<:Integer}}
	mc_state_grid::AbstractArray{<:AbstractVector{<:Real}}

# KEEP GOING

	# Constructor
	function HJBEquation(flow_returns::Function,
						 state_equations::Function,
						 optimal_controls::Function,
						 state_constraints::Function,
						 state_grids_endog::Tuple{Vararg{AbstractVector{<:Real}}},
						 exog_shock::ContinuousTimeMarkovChain,
						 discount_rate::Real)

		# Implement checks as needed

		state_grids = (state_grids_endog..., exog_shock.states)
		state_grids_lengths = length.(state_grids)
		state_grids_strides = Int.(cumprod(state_grids_lengths) ./ state_grids_lengths)
		N = length(state_grids_endog)
		M = length(optimal_controls([state_i[1] for state_i in state_grids],
							 		zeros(Float64, N)))
		nstates = prod([n for n in state_grids_lengths])
		nstates_endog = Int.(nstates / state_grids_lengths[end])
		exog_infinitesimal_generator = kron(exog_shock.infinitesimal_generator,
											sparse(I, nstates_endog, nstates_endog))
		dx = map(x -> x[2] - x[1], state_grids[1:N])
		mc_state_grid = collect.(Iterators.product(state_grids...))

		new(discount_rate,
			flow_returns,
			state_equations,
			optimal_controls,
			state_constraints,
			exog_infinitesimal_generator,
			state_grids,
			N,
			M,
			nstates,
			dx,
			state_grids_lengths,
			state_grids_strides,
			mc_state_grid)
	end
end

# Overloading the show operator to get pretty-printing behavior.
function Base.show(io::IO, ::MIME"text/plain", hjb::HJBEquation)
	print(io, "Hamilton-Jacobi-Bellman Equation.\n")
	print(io, "\nUse .field notation to get access to members.\n")
end

#---------------------------------------------------------------------
# 						HJB SOLUTION TYPE
# --------------------------------------------------------------------
struct HJBSolution
	infinitesimal_generator::AbstractMatrix{<:Real}
	mc_state_grid::AbstractArray{<:AbstractVector{<:Real}}
	value_function::AbstractArray{<:Real}
	policy_functions::AbstractArray{<:Real}
end

# Overloading the show operator to get pretty-printing behavior.
function Base.show(io::IO, ::MIME"text/plain", hjb::HJBSolution)
	print(io, "Solution of Hamilton-Jacobi-Bellman Equation.\n")
	print(io, "\nUse .field notation to get access to members.\n")
end

# Wrapper for stationary_distribution
function MarkovChains.stationary_distribution(sol::HJBSolution)
	return stationary_distribution(ContinuousTimeMarkovChain(sol.infinitesimal_generator, sol.mc_state_grid))
end

# Wrapper for random_sample with random guess according to statdist if not provided.
function MarkovChains.random_sample(sol::HJBSolution)
	idx_start = rand(stationary_distribution(sol))
	return random_sample(sol, idx_start)
end

function MarkovChains.random_sample(sol::HJBSolution, idx_start::Integer)
	return random_sample(ContinuousTimeMarkovChain(sol.infinitesimal_generator, sol.mc_state_grid), idx_start)
end

#---------------------------------------------------------------------
# 						HJB SOLVE ROUTINE
# --------------------------------------------------------------------

# Wrapper for case in which guess is not provided. Use naive guess of controller who sets drifts to zero.
function solve(hjb::HJBEquation)

	# Construct standard guess based on "stay put" policy.
	∇V₀ = hjb.state_constraints.(hjb.mc_state_grid[:])
	c₀ = hjb.optimal_controls.(hjb.mc_state_grid[:], ∇V₀)
	guess = hjb.flow_returns.(hjb.mc_state_grid[:], c₀) ./ hjb.discount_rate

	return solve(hjb, guess)
end

function solve(hjb::HJBEquation,
			   guess::AbstractArray;  # THIS NEEDS TO BE OPTIONAL PARAMETER
			   step_size::Tr=1e3,
			   tol::Tr=1e-8,
			   maxit::Ti where Ti<:Integer=1000) where Tr<:Real

	# Parameters for loop
	Δ = step_size
	ρ = hjb.discount_rate

	# GUESS OF THE VALUE FUNCTION
	V_next = guess

	# PREALLOCATE OBJECTS
	infinitesimal_generator = sparse([], [], [], hjb.nstates, hjb.nstates)
	policy_functions = Matrix{Float64}(undef, hjb.nstates, hjb.M)

	c₊ = Matrix{Float64}(undef, hjb.nstates, hjb.M)
	c₋ = similar(c₊)

	ẋ₊ = Matrix{Float64}(undef, hjb.nstates, hjb.N)
	ẋ₋ = similar(ẋ₊)

	up = Matrix{Bool}(undef, hjb.nstates, hjb.N)
	dwn = similar(up)
	stay = similar(up)

	R = Vector{Float64}(undef, hjb.nstates)

	lin_idx = Vector{Int64}(1:hjb.nstates)

	function controls(states::Vector{<:Int64}, ∇V::Matrix{Float64})::Matrix{Float64}
		c = Matrix{Float64}(undef, length(states), hjb.M)
		for state in states
			c[state, :] .= hjb.optimal_controls(hjb.mc_state_grid[state], ∇V[state,:])
		end
		return c
	end

	function drifts(states::Vector{<:Int64}, c::Matrix{Float64})::Matrix{Float64}
		ẋ = Matrix{Float64}(undef, length(states), hjb.N)
		dx = [dx for dx in hjb.dx[1:hjb.N]]
		for state in states
			ẋ[state, :] .= hjb.state_equations(hjb.mc_state_grid[state], c[state, :]) ./ dx
		end
		# ẋ[abs.(ẋ) .< 1e-16] .= 0.0
		return ẋ
	end

	function returns(states::Vector{<:Int64}, c::Matrix{Float64})::Vector{Float64}
		R = Vector{Float64}(undef, length(states))
		for state in states
			R[state] = hjb.flow_returns(hjb.mc_state_grid[state], c[state, :])
		end
		return R
	end

	function constraints(states::Vector{<:Int64})
		∇R = Matrix{Float64}(undef, length(states), hjb.N)
		for state in states
			∇R[state, :] .= hjb.state_constraints(hjb.mc_state_grid[state])
		end
		return ∇R
	end

	# Enforce state constraints here!
	∇V₀ = constraints(lin_idx)
	∇V₀_cart = reshape(∇V₀, (hjb.state_grids_lengths..., hjb.N)...)
	∇V₊_cart = copy(∇V₀_cart)
	∇V₋_cart = copy(∇V₀_cart)
	c₀ = controls(lin_idx, ∇V₀)

	# Create sparse matrix rows, cols, vals triplet. Start with shock
	rows = Vector{Int64}(undef, 3 * (hjb.N+1) * hjb.nstates)
	cols = similar(rows)
	vals = Vector{Float64}(undef, 3 * (hjb.N+1) * hjb.nstates)
	r_mc , c_mc, v_mc = findnz(hjb.exog_infinitesimal_generator)
	len_mc = length(r_mc)
	ran = 1:len_mc
	rows[ran] = view(r_mc, :)
	cols[ran] = view(c_mc, :)
	vals[ran] = view(v_mc, :)

	# Initialize iterators
	p = ProgressThresh(log10(tol), desc="Iterating until convergence. Showing log₁₀(supnorm(Vₙ₊₁, Vₙ)): ")
	sup = 1e8
	it = 0
	while it < maxit

		# Update
		it += 1
		V = reshape(V_next, hjb.state_grids_lengths...)

		# Compute marginal value for each endogenous state, imposing the state constraints
		for dim in 1:hjb.N

			view_∇V₊_cart = selectdim(∇V₊_cart[axes(V)..., dim], dim, 1:hjb.state_grids_lengths[dim]-1)
			view_∇V₋_cart = selectdim(∇V₋_cart[axes(V)..., dim], dim, 2:hjb.state_grids_lengths[dim])

			∇V₊_cart[parentindices(view_∇V₊_cart)..., dim] = (diff(V, dims = dim) ./ hjb.dx[dim])[:]
			∇V₋_cart[parentindices(view_∇V₋_cart)..., dim] = (diff(V, dims = dim) ./ hjb.dx[dim])[:]

			view_∇V₊_cart = selectdim(∇V₊_cart[axes(V)..., dim], dim, hjb.state_grids_lengths[dim])
			view_∇V₋_cart = selectdim(∇V₋_cart[axes(V)..., dim], dim, 1)
		end

		c₊ = controls(lin_idx, reshape(∇V₊_cart, (hjb.nstates, hjb.N)...))
		c₋ = controls(lin_idx, reshape(∇V₋_cart, (hjb.nstates, hjb.N)...))

		ẋ₊ = drifts(lin_idx, c₊)
		ẋ₋ = drifts(lin_idx, c₋)

		# Upwinding
		up .= ẋ₊ .> 0.0
		dwn .= ẋ₋ .< 0.0
		stay .= (.~ up) .& (.~ dwn)

		# Construct upwind scheme
		# ∇V = reshape(∇V₊_cart, (hjb.nstates, hjb.N)...) .* up +
		# 	 reshape(∇V₋_cart, (hjb.nstates, hjb.N)...) .* dwn +
		# 	 ∇V₀ .* stay

		# Optimal choice of controls
		# c = controls(lin_idx, ∇V)
		c = c₊ .* up .+ c₋ .* dwn .+ c₀ .* stay

		# Maximized returns function
		R = returns(lin_idx, c)

		accumul = len_mc
		for dim in 1:hjb.N

			# Going UP
			up_dim = view(up, :, dim)
			len = sum(up_dim)
			ran = accumul .+ (1:len)
			rows[ran] = view(lin_idx, up_dim)
			cols[ran] = view(lin_idx, up_dim) .+ hjb.state_grids_strides[dim]
			vals[ran] = view(ẋ₊, up_dim, dim)
			accumul += len
			ran = accumul .+ (1:len)
			rows[ran] = view(lin_idx, up_dim)
			cols[ran] = view(lin_idx, up_dim)
			vals[ran] = .- view(ẋ₊, up_dim, dim)
			accumul += len

			# Going DOWN
			dwn_dim = view(dwn, :, dim)
			len = sum(dwn_dim)
			ran = accumul .+ (1:len)
			rows[ran] = view(lin_idx, dwn_dim)
			cols[ran] = view(lin_idx, dwn_dim) .- hjb.state_grids_strides[dim]
			vals[ran] = .- view(ẋ₋, dwn_dim, dim)
			accumul += len
			ran = accumul .+ (1:len)
			rows[ran] = view(lin_idx, dwn_dim)
			cols[ran] = view(lin_idx, dwn_dim)
			vals[ran] = view(ẋ₋, dwn_dim, dim)
			accumul += len

		end

		# Build Sparse transition matrix A (infinitesimal generator)
		A = sparse(view(rows, 1:accumul),
				   view(cols, 1:accumul),
				   view(vals, 1:accumul),
				   hjb.nstates,
				   hjb.nstates)

		# Compute value function update
		V_next = (sparse( (1/Δ + ρ)I, hjb.nstates, hjb.nstates) .- A) \
				 (R + reshape(V, hjb.nstates) ./ Δ)

		# Compute supnorm
		sup = maximum(abs.(V[:] - V_next[:]))
		update!(p, log10(sup))

		# Decision point for convergence
		if sup <= tol
			infinitesimal_generator = A
			policy_functions = c
			break
		end
	end

	if (it == maxit) & (sup > tol)
		error("Reached iteration $it/$maxit without converging.")
	end

	return HJBSolution(infinitesimal_generator,
					   hjb.mc_state_grid[:],
					   V_next,
					   policy_functions)
end


text:
@symbol_func function K_LP(cur_reactor::AbstractReactor)
  cur_K = cur_reactor.C_ejima

  cur_K += cur_reactor.l_i / 2

  cur_K += _external_inductance_term(cur_reactor)

  cur_K *= ( 4e-7 * cur_reactor.pi )

  cur_K
end

@symbol_func function _external_inductance_term(cur_reactor::AbstractReactor)
  cur_epsilon = cur_reactor.epsilon

  cur_kappa = kappa_x(cur_reactor)

  cur_numerator = _b_HN(cur_epsilon)

  cur_numerator -= _a_HN(cur_epsilon)

  cur_denominator = cur_kappa

  cur_denominator *= _d_HN(cur_epsilon)

  cur_denominator /= ( 1 - cur_epsilon )

  cur_denominator += 1

  cur_term = cur_numerator

  cur_term /= cur_denominator

  cur_term
end

function _a_HN(cur_epsilon::SymEngine.Basic, use_default::Bool=false)
  use_default && return _a_HN_default(cur_epsilon)
  symbols(:a_HN)
end

function _b_HN(cur_epsilon::SymEngine.Basic, use_default::Bool=false)
  use_default && return _b_HN_default(cur_epsilon)
  symbols(:b_HN)
end

function _d_HN(cur_epsilon::SymEngine.Basic, use_default::Bool=false)
  use_default && return _d_HN_default(cur_epsilon)
  symbols(:d_HN)
end

function _a_HN(cur_epsilon::AbstractFloat)
  _a_HN_default(cur_epsilon)
end

function _b_HN(cur_epsilon::AbstractFloat)
  _b_HN_default(cur_epsilon)
end

function _d_HN(cur_epsilon::AbstractFloat)
  _d_HN_default(cur_epsilon)
end

function _a_HN_default(cur_epsilon::AbstractSymbol)
  cur_a = 2.0

  cur_a += 9.25 * sqrt(cur_epsilon)

  cur_a -= 1.21 * cur_epsilon

  cur_a
end

function _b_HN_default(cur_epsilon::AbstractSymbol)
  cur_b = 1.0

  cur_b += 1.81 * sqrt(cur_epsilon)

  cur_b += 2.05 * cur_epsilon

  cur_b *= log( 8.0 / cur_epsilon )

  cur_b
end

function _d_HN_default(cur_epsilon::AbstractSymbol)
  cur_d = 1.0

  cur_d += 2 * cur_epsilon ^ 4

  cur_d -= 6 * cur_epsilon ^ 5

  cur_d += 3.7 * cur_epsilon ^ 6

  cur_d *= 0.73 * sqrt(cur_epsilon)

  cur_d
end

ext: jl


text:
cauchy_contour(a, radius = 1, ::Type{T} = real(typeof(a))) where {T} =
    radius*ComplexUnitCircle{T}() .+ a

export cauchy_integral
"""
    cauchy_integral(f [; a = 0, n  = 1, radius = 1])

Evaluate the integral of `n!/(2πi) f(z)/(z-a)^n` along a circular contour
with radius `radius` around the point `a`.

If the integrand is analytic inside the contour except at the pole around `z=a`
then its value should correspond to the derivative `f^(n)(a)`.
"""
function cauchy_integral(f; a = zero(ComplexF64), radius = 1, n = 1)
    T = real(typeof(a))
    factorial(n)/(2*T(pi)*im) *
        integral( z-> f(z)/(z-a)^(n+1), cauchy_contour(a, radius))
end


text:
import BinaryProvider

const _pkg_assets = joinpath(dirname(@__DIR__), "assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

const tmp = joinpath(_pkg_assets, "tmp")
!isdir(tmp) && mkdir(tmp)

download("https://github.com/KaTeX/KaTeX/releases/download/v0.10.0/katex.tar.gz", joinpath(tmp, "katex.tar.gz"))

BinaryProvider.gen_unpack_cmd(joinpath(tmp, "katex.tar.gz"), tmp) |> run

const katexdir = joinpath(tmp, "katex")
const katexfontdir = joinpath(katexdir, "fonts")

const fontdir = joinpath(_pkg_assets, "fonts")
!isdir(fontdir) && mkdir(fontdir)
for font in readdir(katexfontdir)
    mv(joinpath(katexfontdir, font), joinpath(_pkg_assets, "fonts", font), force = true)
end

for file in ["katex.min.js", "katex.min.css", joinpath("contrib", "auto-render.min.js")]
    mv(joinpath(katexdir, file), joinpath(_pkg_assets, splitdir(file)[2]), force = true)
end

rm(tmp, recursive=true)


text:
# Rohan McLure, Australian National University (2019)
# Distributed implementation of the (Star) Stencil Parallel Research Kernel
# Making use of the ArrayChannels.jl library for Julia
# https://github.com/rohanmclure/ArrayChannels.jl
@everywhere using ArrayChannels

function main()
    if length(ARGS) != 3 || nprocs() == 1
        println("Usage: julia -p <# workers> stencil_array_channels.jl <# iterations> <# order> <# radius>")
        exit(1)
    end

    argv = map(x -> parse(Int, x), ARGS)
    iterations, order, r = argv
    procs_x, procs_y = factor(nworkers())

    # Channel setup
    out_channels = Matrix{Vector{Union{ArrayChannel,Nothing}}}(undef, procs_y, procs_x)
    in_channels = Matrix{Vector{Union{ArrayChannel,Nothing}}}(undef, procs_y, procs_x)
    dims = Matrix{Tuple{Int,Int}}(undef, procs_y, procs_x)

    # Dimension setup
    std_width = Int(ceil(order / procs_x))
    std_height = Int(ceil(order / procs_y))
    last_width = order % std_width == 0 ? std_width : order % std_width
    last_height = order % std_height == 0 ? std_height : order % std_height

    for i in 1:procs_y
        for j in 1:procs_x
            h = i == procs_y ? last_height : std_height
            w = j == procs_x ? last_width : std_width

            v = Vector{Union{ArrayChannel, Nothing}}(undef, 4)
            id = workers()[(i-1)*procs_x+j]
            v[1] = i != 1 ? ArrayChannel(Float64, [id, id-procs_x], r, w) : nothing
            v[2] = j != procs_x ? ArrayChannel(Float64, [id, id+1], h, r) : nothing
            v[3] = i != procs_y ? ArrayChannel(Float64, [id, id+procs_x], r, w) : nothing
            v[4] = j != 1 ? ArrayChannel(Float64, [id, id-1], h, r) : nothing

            dims[i,j] = (h,w)
            out_channels[i,j] = v
        end
    end

    for i in 1:procs_y
        for j in 1:procs_x
            v = Vector{Union{ArrayChannel,Nothing}}(undef, 4)
            v[1] = i != 1 ? out_channels[i-1,j][3] : nothing
            v[2] = j != procs_x ? out_channels[i,j+1][4] : nothing
            v[3] = i != procs_y ? out_channels[i+1,j][1] : nothing
            v[4] = j != 1 ? out_channels[i,j-1][2] : nothing
            in_channels[i,j] = v
        end
    end
    futures = Matrix{Future}(undef, procs_y, procs_x)
    @sync for i in 1 : procs_y
        for j in 1:procs_x
            futures[i,j] = @spawnat workers()[(i-1)*procs_x+j] begin
                do_stencil(order, iterations, dims[i,j], (procs_x, procs_y), r, out_channels[i,j], in_channels[i,j])
            end
        end
    end

    # Concatenate the outputs
    results = map(fetch, futures)
    norm = sum(map(x->x[1], results))
    time = maximum(map(x->x[2], results))

    # Verify the result
    active_points = (order-2*r)^2
    norm /= active_points
    epsilon=1.e-8
    reference_norm = 2*(iterations+1)
    if abs(norm-reference_norm) < epsilon
        println("Solution validates")
        stencil_size = 4*r+1 # Star
        flops = (2*stencil_size+1) * active_points
        avgtime = time/iterations
        println("Rate (MFlops/s): $(1.e-6*flops/avgtime) Avg time (s):  $avgtime")
    else
        println("ERROR: L1 norm = $norm Reference L1 norm = $reference_norm")
        exit(9)
    end
end

@everywhere function do_stencil(
    order,
    iterations::Int,
    local_dims::Tuple{Int,Int},
    proc_dims::Tuple{Int,Int},
    r::Int,
    out_channels::Vector{Union{ArrayChannel,Nothing}},
    in_channels::Vector{Union{ArrayChannel,Nothing}}
)
    # Flags for whether workers exist above, rightward ... of this rank
    ub, rb, db, lb = out_channels .!== nothing
    h, w = local_dims

    # Initialise A
    A = zeros(Float64, (local_dims .+ 2*r)...)
    A_central = view(A, r+1:r+h, r+1:r+w)
    procs_x, procs_y = proc_dims
    id = myid() - 2
    p_i, p_j = (id ÷ procs_x, id % procs_x)
    # Start coordinate in global coordinates
    std_height, std_width = Int(ceil(order / procs_y)), Int(ceil(order / procs_x))
    i_start, j_start = p_i * std_height, p_j * std_width
    i_o_start, j_o_start = (ub ? 1 : r+1), (lb ? 1 : r+1)
    i_o_end, j_o_end = (db ? h : h-r), (rb ? w : w-r)
    for i in 1:h
        for j in 1:w
            A_central[i,j] = (p_i*std_height + i - 1.0) + (p_j*std_width + j - 1.0)
        end
    end

    # Initialise output array
    B = zeros(Float64, local_dims...)
    # Discrete Divergence operator
    W = zeros(Float64, 2*r+1, 2*r+1)
    for i in 1:r
        W[r+1, r+i+1] = 1.0/(2*i*r)
        W[r+i+1, r+1] = 1.0/(2*i*r)
        W[r+1, r-i+1] = -1.0/(2*i*r)
        W[r-i+1, r+1] = -1.0/(2*i*r)
    end

    u_out_v, r_out_v, d_out_v, l_out_v = [
        view(A, r+1:2*r, r+1:r+w),
        view(A, r+1:r+h, w+1:r+w),
        view(A, h+1:r+h, r+1:r+w),
        view(A, r+1:r+h, r+1:2+r)
    ]
    u_in_v, r_in_v, d_in_v, l_in_v = [
        view(A, 1:r, r+1:r+w),
        view(A, r+1:r+h, r+w+1:2*r+w),
        view(A, r+h+1:2*r+h, r+1:r+w),
        view(A, r+1:r+h, 1:r)
    ]
    # Proc ID references
    id = myid()
    u_id, r_id, d_id, l_id = [
        id - procs_x,
        id + 1,
        id + procs_x,
        id - 1
    ]

    # Channels
    u_out_c, r_out_c, d_out_c, l_out_c = out_channels
    u_in_c, r_in_c, d_in_c, l_in_c = in_channels

    local t0, t1

    for k in 0 : iterations
        if k == 1
            t0 = time_ns()
        end
        # Send / Receive
        @sync begin
            if ub
                @async begin
                    copy!(u_out_c, u_out_v)
                    put!(u_out_c, u_id)
                end
                @async begin
                    take!(u_in_c, u_id)
                    copy!(u_in_v, u_in_c)
                end
            end
            if rb
                @async begin
                    copy!(r_out_c, r_out_v)
                    put!(r_out_c, r_id)
                end
                @async begin
                    take!(r_in_c, r_id)
                    copy!(r_in_v, r_in_c)
                end
            end
            if db
                @async begin
                    copy!(d_out_c, d_out_v)
                    put!(d_out_c, d_id)
                end
                @async begin
                    take!(d_in_c, d_id)
                    copy!(d_in_v, d_in_c)
                end
            end
            if lb
                @async begin
                    copy!(l_out_c, l_out_v)
                    put!(l_out_c, l_id)
                end
                @async begin
                    take!(l_in_c, l_id)
                    copy!(l_in_v, l_in_c)
                end
            end
        end

        for j in j_o_start : j_o_end
            for i in i_o_start : i_o_end
                # Star case
                tmp = 0.0
                for jj in -r:r
                    @inbounds tmp += W[r+1, r+1+jj] * A[i+r, j+jj+r]
                end
                for ii in -r:-1
                    @inbounds tmp += W[r+1+ii, r+1] * A[i+ii+r, j+r]
                end
                for ii in 1:r
                    @inbounds tmp += W[r+1+ii, r+1] * A[i+ii+r, j+r]
                end
                @inbounds B[i,j] += tmp
            end
        end

        for j in 1:w
            for i in 1:h
                @inbounds A[i+r, j+r] += 1.0
            end
        end
    end

    t1 = time_ns()

    norm = 0.0
    for i in i_o_start : i_o_end
        for j in j_o_start : j_o_end
            norm += abs(B[i,j])
        end
    end
    return (norm, (t1-t0)*1.e-9)
end

function factor(x::Int)
    s = Int(floor(sqrt(x)))
    for p in s:-1:1
        if x % p == 0
            return (p, x ÷ p)
        end
    end
end

main()


text:
# using Ipaper
# using Test

@testset "missing" begin
  @test drop_missing([1, 2, missing]) == [1, 2, 0]
  @test drop_missing([1, 2, missing], -999) == [1, 2, -999]
  @test drop_missing([1, 2, 3], -999) == [1, 2, 3]

  @test isequal(to_missing([1, 2, 0]), [1, 2, missing])
  @test isequal(to_missing([1, 2, 3], 3), [1, 2, missing])

  @test replace_value!([1, 2], 2, 3) == [1, 3]
end


text:
const STOPWORDS_en = ["a", "about", "above", "across", "after", "afterwards", "again", "against", "ain", "all", "almost", "alone", "along", "already", "also", "although", "always", "am", "among", "amongst", "amount", "an", "and", "another", "any", "anyhow", "anyone", "anything", "anyway", "anywhere", "are", "aren", "aren't", "around", "as", "at", "back", "be", "became", "because", "become", "becomes", "becoming", "been", "before", "beforehand", "behind", "being", "below", "beside", "besides", "between", "beyond", "both", "bottom", "but", "by", "ca", "call", "can", "cannot", "could", "couldn", "couldn't", "d", "did", "didn", "didn't", "do", "does", "doesn", "doesn't", "doing", "don", "don't", "done", "down", "due", "during", "each", "eight", "either", "eleven", "else", "elsewhere", "empty", "enough", "even", "ever", "every", "everyone", "everything", "everywhere", "except", "few", "fifteen", "fifty", "first", "five", "for", "former", "formerly", "forty", "four", "from", "front", "full", "further", "get", "give", "go", "had", "hadn", "hadn't", "has", "hasn", "hasn't", "have", "haven", "haven't", "having", "he", "hence", "her", "here", "hereafter", "hereby", "herein", "hereupon", "hers", "herself", "him", "himself", "his", "how", "however", "hundred", "i", "if", "in", "indeed", "into", "is", "isn", "isn't", "it", "it's", "its", "itself", "just", "keep", "last", "latter", "latterly", "least", "less", "ll", "m", "ma", "made", "make", "many", "may", "me", "meanwhile", "might", "mightn", "mightn't", "mine", "more", "moreover", "most", "mostly", "move", "much", "must", "mustn", "mustn't", "my", "myself", "name", "namely", "needn", "needn't", "neither", "never", "nevertheless", "next", "nine", "no", "nobody", "none", "noone", "nor", "not", "nothing", "now", "nowhere", "o", "of", "off", "often", "on", "once", "one", "only", "onto", "or", "other", "others", "otherwise", "our", "ours", "ourselves", "out", "over", "own", "part", "per", "perhaps", "please", "put", "quite", "rather", "re", "really", "regarding", "s", "same", "say", "see", "seem", "seemed", "seeming", "seems", "serious", "several", "shan", "shan't", "she", "she's", "should", "should've", "shouldn", "shouldn't", "show", "side", "since", "six", "sixty", "so", "some", "somehow", "someone", "something", "sometime", "sometimes", "somewhere", "still", "such", "t", "take", "ten", "than", "that", "that'll", "the", "their", "theirs", "them", "themselves", "then", "thence", "there", "thereafter", "thereby", "therefore", "therein", "thereupon", "these", "they", "third", "this", "those", "though", "three", "through", "throughout", "thru", "thus", "to", "together", "too", "top", "toward", "towards", "twelve", "twenty", "two", "under", "unless", "until", "up", "upon", "us", "used", "using", "various", "ve", "very", "via", "was", "wasn", "wasn't", "we", "well", "were", "weren", "weren't", "what", "whatever", "when", "whence", "whenever", "where", "whereafter", "whereas", "whereby", "wherein", "whereupon", "wherever", "whether", "which", "while", "whither", "who", "whoever", "whole", "whom", "whose", "why", "will", "with", "within", "without", "won", "won't", "would", "wouldn", "wouldn't", "y", "yet", "you", "you'd", "you'll", "you're", "you've", "your", "yours", "yourself", "yourselves"]
const STOPWORDS_es = ["agregó", "poder", "debe", "cerca", "intentais", "está", "todas", "durante", "pesar", "estoy", "ello", "buenas", "hacemos", "adelante", "mismas", "trabajo", "quién", "desde", "cierta", "suyas", "próximos", "ahí", "ningunos", "casi", "dicho", "cual", "ex", "vais", "empleais", "despues", "siendo", "general", "sean", "otras", "deprisa", "dónde", "esta", "aún", "nueva", "cuántos", "intenta", "porque", "bueno", "podría", "vuestras", "debido", "través", "enseguida", "lleva", "sé", "alrededor", "horas", "trabajar", "mía", "hacen", "medio", "algún", "mayor", "habrá", "conseguimos", "mis", "tendrá", "mios", "aproximadamente", "eramos", "buen", "además", "emplean", "quiere", "conocer", "expresó", "ti", "empleo", "buenos", "cuanto", "sus", "no", "tras", "raras", "ante", "cuándo", "posible", "estan", "vaya", "vuestro", "éstas", "arribaabajo", "arriba", "trabajas", "cualquier", "mismo", "salvo", "sino", "este", "contra", "delante", "soyos", "asi", "de", "dentro", "informó", "podriamos", "sabemos", "trabaja", "hacer", "eres", "aquellas", "partir", "otros", "quienes", "cuánta", "mediante", "qué", "empleas", "haces", "aquí", "todo", "parece", "dicen", "ser", "vuestros", "somos", "tuvo", "total", "nos", "va", "quiza", "segundo", "si", "cuánto", "luego", "estar", "ampleamos", "hacia", "mia", "informo", "también", "existe", "igual", "eso", "te", "ella", "nadie", "estas", "ningunas", "largo", "consideró", "llegó", "quizás", "momento", "diferentes", "ademas", "nuestros", "llevar", "detras", "dieron", "propio", "ese", "sois", "cuáles", "ciertos", "suya", "alguna", "manifestó", "esa", "veces", "actualmente", "todavía", "muchas", "tengo", "tan", "serán", "hizo", "nuevo", "gran", "se", "intentamos", "es", "demás", "cuantas", "vosotros", "aquellos", "tiempo", "tuyos", "lugar", "cuatro", "estamos", "toda", "alguno", "como", "hago", "propios", "mal", "ocho", "tres", "diferente", "el", "mí", "adrede", "breve", "hemos", "fueron", "nosotros", "verdad", "nunca", "sobre", "pasada", "después", "quizá", "con", "intento", "podria", "contigo", "mencionó", "quedó", "modo", "misma", "apenas", "ha", "encima", "consigue", "nuevos", "pudo", "poca", "consigues", "otra", "bajo", "aquello", "saben", "hacerlo", "tuya", "fuera", "dia", "cuando", "dio", "mejor", "ellas", "tenga", "al", "tambien", "conseguir", "podemos", "estados", "hicieron", "tenía", "hablan", "despacio", "según", "junto", "así", "lejos", "último", "conmigo", "habia", "pocos", "antaño", "cierto", "vamos", "incluso", "día", "queremos", "temprano", "demasiado", "más", "supuesto", "muchos", "primera", "algunas", "propia", "debajo", "ninguno", "tú", "estado", "mismos", "pero", "proximo", "usas", "esos", "podrán", "muy", "repente", "trabajan", "ya", "usamos", "emplear", "nuevas", "dos", "ésas", "intentas", "cuanta", "soy", "tercera", "quizas", "ver", "hasta", "considera", "las", "cuenta", "alli", "aquéllos", "aquéllas", "decir", "ahi", "usais", "paìs", "sabes", "próximo", "ultimo", "usan", "lado", "indicó", "pronto", "éstos", "podrían", "aquél", "primeros", "cuántas", "sigue", "ahora", "intentan", "ningún", "os", "teneis", "añadió", "cuál", "estos", "primer", "unas", "estaba", "excepto", "bien", "antano", "solamente", "aseguró", "pocas", "buena", "tanto", "cuantos", "afirmó", "podriais", "por", "pues", "aqui", "yo", "acuerdo", "ellos", "él", "dado", "menos", "comentó", "son", "estará", "tal", "tenemos", "tuyas", "le", "aquélla", "unos", "atras", "realizado", "respecto", "tendrán", "nada", "consiguen", "trabajais", "eras", "sera", "explicó", "intentar", "señaló", "cómo", "cinco", "bastante", "esas", "consigo", "segunda", "haber", "solos", "seis", "propias", "ayer", "final", "había", "sabe", "sí", "claro", "tienen", "mío", "usa", "varias", "fue", "verdadero", "ésos", "mucho", "mi", "fui", "realizó", "algo", "era", "antes", "nosotras", "para", "hay", "haceis", "sólo", "sin", "sola", "enfrente", "cuales", "grandes", "principalmente", "tus", "vez", "sea", "qeu", "segun", "peor", "un", "poco", "la", "manera", "mucha", "embargo", "puedo", "encuentra", "fin", "detrás", "gueno", "su", "anterior", "da", "ir", "tampoco", "tuyo", "vosotras", "en", "lo", "todavia", "han", "parte", "solo", "pasado", "dar", "allí", "aunque", "saber", "creo", "podeis", "últimas", "días", "donde", "deben", "tarde", "esto", "estaban", "ése", "sería", "dias", "cosas", "aquella", "mas", "ni", "fuimos", "voy", "pueda", "aun", "mientras", "pueden", "tu", "tenido", "solas", "tener", "cada", "última", "habían", "entonces", "todos", "siguiente", "haciendo", "tiene", "últimos", "mías", "ninguna", "trabajamos", "ejemplo", "sabeis", "ambos", "ciertas", "dijo", "hubo", "usted", "podrá", "los", "mio", "nuestro", "hace", "les", "usar", "éste", "menudo", "aquel", "realizar", "primero", "hecho", "dijeron", "haya", "habla", "dan", "siempre", "quiénes", "eran", "trata", "uso", "nuestra", "podrias", "ustedes", "quien", "he", "dejó", "varios", "míos", "pais", "mias", "del", "poner", "estuvo", "nuestras", "están", "que", "existen", "vuestra", "hoy", "me", "suyo", "ésta", "valor", "será", "entre", "algunos", "podrian", "puede", "siete", "van", "dice", "ésa", "sido", "otro", "uno", "verdadera", "estais", "una"]
const STOPWORDS_fr = ["allo", "derriere", "onzième", "oust", "vôtres", "durant", "troisièmement", "rendre", "serait", "olé", "sept", "comme", "souvent", "suivants", "soit", "desormais", "hou", "ore", "quand", "lesquelles", "boum", "toi", "prealable", "dix-huit", "sixième", "trois", "ceux-ci", "derrière", "différentes", "hi", "ouf", "houp", "fait", "o", "etant", "dont", "suffisante", "aucune", "suivre", "minimale", "pendant", "sauf", "specifiques", "seul", "vais", "suivantes", "beau", "d’", "possibles", "avait", "leur", "plouf", "peux", "quels", "dixième", "différent", "huitième", "puis", "vlan", "relativement", "anterieures", "était", "bravo", "tiens", "ayant", "seulement", "quelle", "pfft", "suivant", "particulièrement", "quatorze", "pouah", "deuxièmement", "seize", "ohé", "quinze", "autrui", "alors", "tien", "probable", "divers", "avons", "parfois", "pan", "dessus", "ça", "chez", "allô", "celles-là", "notre", "aux", "quel", "etc", "pouvait", "ô", "troisième", "seraient", "apres", "restant", "pure", "soi", "hors", "memes", "de", "auraient", "rare", "soixante", "sacrebleu", "celui-ci", "mon", "ho", "subtiles", "s’", "hep", "qu'", "reste", "on", "celui-là", "pas", "quiconque", "toutes", "par", "deja", "naturelle", "specifique", "nos", "va", "si", "malgre", "doit", "lesquels", "quatre-vingt", "tel", "ces", "depuis", "onze", "brrr", "ce", "tous", "aie", "lui", "te", "celui", "moi", "douze", "cependant", "zut", "mes", "au", "laquelle", "gens", "différente", "voici", "là", "sien", "parle", "nôtres", "c’", "toujours", "quatre", "restrictif", "basee", "trop", "très", "peu", "tes", "huit", "chaque", "procedant", "superpose", "tsoin", "avec", "speculatif", "vont", "parce", "pense", "voilà", "telle", "outre", "se", "vous", "eux", "es", "désormais", "etais", "personne", "elle", "ceci", "proche", "passé", "quoique", "sont", "siennes", "derniere", "aura", "clac", "m'", "néanmoins", "pire", "quant-à-soi", "rares", "flac", "vu", "tres", "meme", "uns", "vif", "laisser", "l'", "t'", "ont", "hue", "ès", "auquel", "suivante", "assez", "quelques", "egales", "fi", "vas", "nombreuses", "etre", "sent", "crac", "ses", "clic", "precisement", "ha", "necessairement", "psitt", "importe", "premièrement", "debout", "encore", "malgré", "font", "ouverts", "vous-mêmes", "remarquable", "auront", "ai", "tend", "possessifs", "couic", "quant", "exterieur", "tels", "plus", "toi-même", "elles-mêmes", "avaient", "moins", "pu", "sait", "tente", "suffit", "sienne", "trente", "suis", "toc", "tenir", "feront", "ouverte", "ailleurs", "enfin", "nous-mêmes", "tellement", "allaient", "bah", "possible", "abord", "pour", "eux-mêmes", "elle-même", "neuvième", "floc", "dring", "stop", "celles", "moindres", "c'", "sa", "vé", "hui", "uniques", "moi-meme", "afin", "celles-ci", "exactement", "faisant", "lès", "neanmoins", "semblable", "et", "treize", "allons", "probante", "las", "tardive", "aucun", "uniformement", "nous", "ou", "retour", "vives", "celle-là", "differentes", "moyennant", "tout", "dix", "devant", "vive", "cette", "cet", "dite", "bien", "lorsque", "absolument", "quelconque", "juste", "hum", "hop", "directe", "ait", "different", "rend", "je", "chacune", "â", "ceux", "nombreux", "na", "après", "cinquante", "egale", "peuvent", "surtout", "hé", "comment", "son", "chère", "nul", "avais", "ta", "auxquels", "effet", "il", "votre", "le", "sein", "avant", "certaines", "ma", "est", "t’", "ceux-là", "notamment", "sera", "delà", "près", "sans", "diverse", "pourrais", "mien", "vifs", "siens", "vivat", "lui-même", "mais", "elles", "première", "qui", "plutôt", "autres", "tsouin", "touchant", "dernier", "plein", "hormis", "suffisant", "suit", "premier", "ils", "n'", "vingt", "avoir", "merci", "cinq", "eu", "contre", "combien", "parlent", "quanta", "étais", "restent", "dans", "celle", "tienne", "dessous", "aussi", "nouveau", "unique", "sapristi", "anterieure", "duquel", "excepté", "fais", "bas", "longtemps", "neuf", "holà", "semblaient", "telles", "hurrah", "plusieurs", "a", "quatrièmement", "donc", "differents", "particulier", "strictement", "chacun", "quarante", "as", "pourrait", "certaine", "semble", "etaient", "mille", "pourquoi", "un", "bat", "tic", "maint", "différents", "la", "nôtre", "oh", "attendu", "lequel", "eh", "aujourd", "leurs", "paf", "diverses", "dix-sept", "douzième", "autrement", "maintenant", "s'", "da", "etait", "n’", "moi-même", "en", "cinquantaine", "desquels", "parler", "vers", "certes", "unes", "lors", "dit", "dehors", "lui-meme", "pfut", "d'", "façon", "j’", "étaient", "tant", "peut", "partant", "desquelles", "hein", "rien", "possessif", "aupres", "auxquelles", "envers", "quelles", "pur", "via", "toutefois", "ni", "ainsi", "revoilà", "chers", "pres", "parseme", "quelque", "tu", "ton", "maximale", "toute", "cher", "autre", "deux", "étant", "cela", "hélas", "celle-ci", "chiche", "faisaient", "environ", "pff", "autrefois", "té", "des", "devra", "hem", "i", "miennes", "chut", "cent", "du", "necessaire", "soi-même", "aujourd'hui", "dits", "ouste", "relative", "les", "certain", "revoici", "concernant", "miens", "doivent", "euh", "non", "tac", "vôtre", "anterieur", "seule", "egalement", "m’", "dire", "sur", "dès", "selon", "aurait", "six", "rarement", "pif", "car", "été", "cinquième", "certains", "seront", "dix-neuf", "même", "quelqu'un", "être", "bigre", "ouvert", "cinquantième", "semblent", "naturelles", "que", "dedans", "beaucoup", "directement", "vos", "ah", "jusque", "sous", "ci", "tenant", "une", "l’", "ouias", "mêmes", "parmi", "me", "septième", "quoi", "naturel", "devers", "où", "tiennes", "qu’", "deuxième", "particulière", "extenso", "à", "chères", "comparables", "entre", "mince", "ne", "permet", "quatrième", "mienne", "puisque", "jusqu", "ollé", "j'", "sinon", "compris", "comparable"]
const STOPWORDS_ru = ["моего", "все", "нашей", "такая", "нею", "его", "ем", "от", "моя", "свою", "своих", "бы", "будучи", "нет", "нему", "сами", "ела", "это", "когда", "одним", "кого", "которою", "ней", "одному", "свое", "этою", "нашею", "моими", "мое", "всей", "одних", "мочь", "ест", "они", "который", "наш", "тою", "он", "кто", "всю", "чтобы", "наше", "такого", "нашим", "всея", "из", "которого", "эти", "по", "те", "них", "едим", "мой", "буду", "может", "на", "ей", "мы", "есть", "одни", "самим", "свой", "которым", "таком", "этих", "мне", "могут", "всеми", "нее", "она", "такой", "для", "такую", "теми", "вот", "которой", "сама", "саму", "её", "наши", "которыми", "что", "оно", "с", "одну", "этот", "будете", "могу", "себе", "такие", "ними", "которую", "тех", "могла", "самом", "нам", "одно", "тот", "ему", "всем", "этими", "но", "своей", "во", "будем", "всём", "вами", "своем", "сам", "эту", "их", "тобою", "том", "ею", "еще", "одна", "своего", "и", "был", "будь", "всех", "комья", "как", "такому", "самими", "можешь", "можете", "одними", "только", "чего", "моё", "моей", "свои", "своя", "тобой", "той", "будешь", "один", "моим", "мог", "своею", "же", "ты", "или", "этому", "могли", "могло", "своим", "ел", "моею", "собой", "наших", "которые", "вся", "моему", "мою", "не", "мои", "я", "ее", "уже", "да", "одного", "этой", "была", "нашем", "можем", "а", "этим", "которому", "само", "эта", "за", "себя", "этом", "кому", "нашему", "будьте", "чему", "наса", "неё", "одном", "которое", "томах", "него", "того", "вам", "емъ", "тому", "нами", "которая", "таким", "тем", "ко", "у", "ещё", "самих", "в", "так", "такое", "своём", "меня", "самого", "всему", "будет", "ими", "моём", "при", "собою", "одной", "ешь", "та", "тебя", "чем", "ним", "оне", "вы", "всего", "к", "об", "будут", "ком", "если", "наша", "ту", "были", "было", "моем", "едят", "могите", "всею", "своими", "своему", "до", "то", "нашими", "нашего", "такою", "нем", "им", "которых", "нас", "имъ", "мною", "одною", "самому", "вас", "кем", "мной", "нём", "всё", "тебе", "котором", "моги", "такими", "чём", "моих", "весь", "нашу", "быть", "о", "своё", "таких", "этого"]
const STOPWORDS_yo = ["fún", "è", "ní", "rẹ̀", "k", "sínú", "àwọn", "n", "j", "á", "ti", "ọ", "mo", "mọ̀", "ù", "ì", "kì", "tí", "ṣùgbọ́n", "p", "máa", "nítorí", "nígbà", "jẹ", "u", "m", "púpọ̀", "í", "padà", "jẹ́", "náà", "bí", "bá", "wọ́n", "ẹ", "wọn", "ni", "àti", "b", "gbogbo", "yìí", "ò", "ṣé", "́", "t", "f", "láti", "̣", "ẹmọ́", "bẹ̀rẹ̀", "ọ̀pọ̀lọpọ̀", "w", "̀", "mi", "ǹ", "sí", "d", "r", "pẹ̀lú", "wà", "inú", "ó", "òun", "wá", "à", "y", "kò", "fẹ́", "i", "e", "ńlá", "sì", "nǹkan", "lè", "kan", "jù", "l", "pé", "é", "s", "o", "ọjọ́", "g", "ń", "ṣe", "an", "a", "lọ", "ú", "kí", "ṣ"]
const STOPWORDS_de = ["dieselbe", "nahm", "drei", "eigen", "vielen", "guter", "nun", "würden", "zu", "sollten", "derselben", "vergangenen", "dasselbe", "grosse", "könnte", "á", "wegen", "bereits", "demgegenüber", "acht", "mögen", "meinem", "rechter", "jedermanns", "weiteres", "dir", "trotzdem", "wohl", "wird", "des", "durchaus", "her", "je", "außer", "machte", "infolgedessen", "erst", "später", "während", "derselbe", "unsere", "bisher", "wir", "großer", "allen", "dürft", "ihre", "hätten", "in", "jahr", "wenn", "dermassen", "ihrem", "gibt", "recht", "erstes", "schon", "überhaupt", "sei", "eigenen", "achte", "hat", "jemand", "eines", "ausser", "sonst", "hoch", "siebentes", "immer", "unter", "gekonnt", "denn", "habe", "hast", "jenes", "einmal", "alles", "sollte", "lange", "jeder", "weniges", "gar", "wem", "seine", "fünfte", "währenddessen", "elf", "gegen", "dahin", "wenigstens", "besser", "zwar", "weiteren", "allein", "seitdem", "kein", "sechstes", "vier", "damals", "zehnten", "endlich", "der", "früher", "wurde", "war", "mögt", "los", "rechte", "zunächst", "dieses", "dasein", "ganzen", "siebte", "nur", "gemusst", "ihm", "jedermann", "große", "drin", "jetzt", "eine", "danach", "solchem", "mussten", "doch", "siebtes", "niemanden", "ausserdem", "etwas", "davon", "dementsprechend", "da", "hatten", "dem", "am", "jemanden", "jahren", "hier", "er", "sollen", "sowie", "kann", "tag", "nie", "deren", "auf", "magst", "hatte", "zusammen", "konnten", "zehn", "ganz", "diesem", "desselben", "statt", "nachdem", "gab", "meines", "seien", "einigen", "neunter", "zehntes", "weitere", "im", "manches", "würde", "hin", "bald", "als", "daneben", "lang", "kleiner", "willst", "sieben", "mag", "ach", "gleich", "daß", "deinem", "selbst", "gesagt", "siebente", "zwei", "darfst", "möchte", "beide", "dies", "denselben", "worden", "ihn", "allgemeinen", "das", "wie", "zwanzig", "sind", "zeit", "demzufolge", "demgemäss", "groß", "kleine", "uhr", "fünftes", "gut", "nicht", "meinen", "so", "nein", "ihres", "welcher", "ihren", "zweites", "weniger", "manchen", "darf", "dann", "einige", "allem", "ab", "ist", "anderen", "bekannt", "möglich", "werde", "meiner", "dabei", "sein", "mir", "dank", "dein", "diejenige", "soll", "jeden", "manchem", "kam", "jener", "dafür", "heißt", "dieselben", "beim", "dich", "leicht", "weit", "solang", "aller", "zur", "über", "demselben", "alle", "ihr", "ersten", "viele", "jene", "en", "einen", "beispiel", "durfte", "noch", "rechten", "dessen", "keinen", "tat", "ganzes", "dagegen", "wessen", "einander", "seinem", "kommen", "richtig", "neuntes", "eigener", "seinen", "großen", "ihrer", "a", "darunter", "darauf", "ehrlich", "mittel", "bin", "waren", "was", "übrigens", "dort", "vom", "wollte", "darum", "vierten", "dritter", "drittes", "macht", "wenige", "unserer", "genug", "einiges", "neunten", "bist", "mochte", "dahinter", "ein", "siebenten", "unser", "allerdings", "daran", "dritte", "gerade", "sehr", "muß", "dazu", "den", "dieser", "teil", "muss", "erster", "damit", "neben", "neue", "mich", "oben", "sie", "kannst", "gewesen", "keiner", "könnt", "keinem", "nach", "mehr", "will", "musste", "einiger", "deiner", "jedem", "natürlich", "eigenes", "ende", "schlecht", "achten", "werden", "zweiter", "jedoch", "sechste", "sechster", "jahre", "ich", "tel", "wann", "zweite", "vergangene", "gutes", "neuen", "seines", "die", "gekannt", "wollt", "offen", "heisst", "ins", "rechtes", "dadurch", "satt", "ihnen", "solches", "zehnter", "kurz", "grosses", "von", "können", "seiner", "na", "niemandem", "zehnte", "uns", "anderem", "einem", "solchen", "hinter", "währenddem", "geworden", "rund", "jemandem", "manche", "einmaleins", "grossen", "durften", "siebter", "welche", "siebten", "diejenigen", "denen", "achter", "müsst", "weil", "diesen", "deswegen", "kleines", "sechsten", "solche", "ganze", "müssen", "nichts", "ag", "derjenigen", "fünfter", "mochten", "hätte", "konnte", "sondern", "wahr", "vielleicht", "auch", "sich", "wollten", "indem", "dürfen", "grosser", "lieber", "haben", "etwa", "daher", "daraus", "darüber", "eben", "davor", "ohne", "zweiten", "fünf", "mancher", "zuerst", "wo", "außerdem", "jenem", "sagte", "man", "andere", "neun", "zurück", "gemocht", "welches", "oder", "gehabt", "gewollt", "heute", "um", "geht", "geschweige", "niemand", "dazwischen", "dermaßen", "also", "musst", "wurden", "jenen", "demgemäß", "keine", "wer", "durch", "tage", "daselbst", "deshalb", "wirklich", "mein", "darin", "vielem", "viel", "bis", "euch", "vierte", "entweder", "vor", "gedurft", "ja", "morgen", "dritten", "werdet", "siebenter", "du", "leider", "irgend", "seit", "machen", "wollen", "besonders", "zum", "ob", "zwischen", "warum", "wart", "ganzer", "gegenüber", "kaum", "tun", "andern", "wäre", "achtes", "aber", "ging", "einer", "besten", "wenig", "kommt", "gute", "weiter", "und", "großes", "vierter", "erste", "gemacht", "mit", "fünften", "für", "wirst", "kleinen", "viertes", "jede", "meine", "aus", "solcher", "gern", "diese", "welchem", "oft", "beiden", "habt", "sechs", "zugleich", "sah", "bei", "gross", "es", "tagen", "dass", "sagt", "seid", "wen", "eigene", "derjenige", "neunte", "ebenso", "an", "deine", "anders", "wieder", "welchen", "gehen"]


text:
depth = divedivedive(0)
println("A clean dive reaches a depth of ", depth, ".")

depth = 0
try
    divedivedive()
end
println("A dirty dive reaches a depth of ", depth, ".")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BioMedQuery.Processes
using BioMedQuery.Entrez
using BioMedQuery.UMLS
using JLD
using MySQL
using BCBIVizUtils
using Colors
using JLD

function load_occurrance_file(; results_dir = "results/",)
    file  = jldopen(results_dir* "/occur_sp.jdl", "r")
    occur = read(file, "occur")
    # display(occur)
    file  = jldopen(results_dir* "/mesh2ind.jdl", "r")
    mesh2ind = read(file, "mesh2ind")

    # display(mesh2ind)
    #invert map
    ind2label = Dict()
    for key  in keys(mesh2ind)
        val = mesh2ind[key]
        ind2label[val] = key
    end
    labels = Array{ASCIIString}(length(ind2label))
    for i in range(1,length(ind2label))
        labels[i] = ind2label[i]
    end

    coo = occur* occur.'

    graph = BCBIVizUtils.coo2graph(coo)
    colors = distinguishable_colors(length(ind2label), RGB(1,0,0))

    return graph, labels, colors
end

function plot3D()
    graph, labels, colors = init()
    BCBIVizUtils.graph2plotlyjs3D(graph, labels, colors)
end

"""
    epilepsy_comorbidities(;    results_dir = "results/",
                                max_articles = typemax(Int64),
                                run_pubmed_search_and_save = true,
                                run_mesh2umls_map = true,
                                filter_semantic_occurrences  = true,
                                comorbidity_concept = "Disease or Syndrome")

PuMed/Medline comorbidities for Epilepsy

Code by: Isabel Restrepo
BCBI - Brown University
Version: Julia 0.6

Before running, configure your MySQL onnection by adding to ~/.juliarc.jl
the following variables:

ENV["PUBMEDMINER_DB_HOST"]="db_host"
ENV["PUBMEDMINER_DB_USER"]="your_user"
ENV["PUBMEDMINER_DB_PSSWD"]="your_password"

You alseo need to configure the email addess associated with you NCBI account,
and USER and PASSWD associated with your UMLS/UTS account.
"""
function epilepsy_comorbidities(; results_dir = "results/",
                                  max_articles = typemax(Int64),
                                  run_pubmed_search_and_save = true,
                                  run_mesh2umls_map = false,
                                  filter_semantic_occurrences  = false,
                                  comorbidity_concept = "Disease or Syndrome",
                                  verbose = false, overwrite=true)

    # 1. Database configuration
    db_host = ENV["PUBMEDMINER_DB_HOST"]
    mysql_usr=ENV["PUBMEDMINER_DB_USER"]
    mysql_pswd=ENV["PUBMEDMINER_DB_PSSWD"]
    dbname="pubmed_comorbidities"
    db = nothing

    db_config = Dict(:host=>db_host,
                  :dbname=>dbname,
                  :username=>mysql_usr,
                  :pswd=>mysql_pswd,
                  :overwrite=>overwrite,
                  :tablename=>"article_epilepsy")

    # 2. PubMed Search
    info("----------------Start: pubmed_search_and_save")
    @time if run_pubmed_search_and_save
        email= ENV["NCBI_EMAIL"]
        search_term="epilepsy[mh]"

        db = pubmed_pmid_search_and_save(email, search_term, max_articles,
        save_pmid_mysql, db_config, verbose)
    end
    info("----------------Done: pubmed_search_and_save")

    # #3. MAP all MeSH to UMLS concepts
    # info("----------------Start: run_mesh2umls_map")
    # @time if run_mesh2umls_map

    #     if db == nothing
    #         db = mysql_connect(db_host, mysql_usr, mysql_pswd, dbname)
    #     end
    #     user = ENV["UMLS_USER"]
    #     psswd = ENV["UMLS_PSSWD"]
    #     credentials = Credentials(user, psswd)
    #     append = false

    #     map_mesh_to_umls_async!(db, credentials; append_results=append)
    # end
    # info("----------------Done: run_mesh2umls_map")

    #4. Filter those that are only Psychiatric Mental Disorders
    info("----------------Start: filter_semantic_occurrences")
    @time if filter_semantic_occurrences

        if db == nothing
            db = mysql_connect(db_host, mysql_usr, mysql_pswd, dbname)
        end

        if !isdir(results_dir)
            mkdir(results_dir)
        end

        occur_path = results_dir*"/occur_sp.jdl"
        labels2ind_path = results_dir*"/labels2ind.jdl"

        labels2ind, occur = umls_semantic_occurrences(db, comorbidity_concept)

        println("-------------------------------------------------------------")
        println("Output Descritor to Index Dictionary")
        println(labels2ind)
        println("-------------------------------------------------------------")

        println("-------------------------------------------------------------")
        println("Output Data Matrix")
        println(occur)
        println("-------------------------------------------------------------")

        # save(occur_path, "occur", occur)
        jldopen(occur_path, "w") do file
         write(file, "occur", occur)
        end
        jldopen(labels2ind_path, "w") do file
         write(file, "labels2ind", labels2ind)
        end
    end
    info("----------------Done: filter_semantic_occurrences")

    #5. Plot comorbidities graph

end

# epilepsy_comorbidities(results_dir = "/Users/isa/dropbox_brown/results/epilepsy_pubmed_comorbidities",
# max_articles = 5, run_pubmed_search_and_save = false)
epilepsy_comorbidities(results_dir = "/Users/isa/dropbox_brown/results/epilepsy_pubmed_comorbidities")

#---------------times:

# INFO: ----------------Done: pubmed_search_and_save
# 52.191503 seconds (1.56 M allocations: 98.617 MiB, 0.06% gc time)
#Maybe this time can be improved by saving all ids at once.

# INFO: ----------------Done: run_mesh2umls_map

# INFO: ----------------Done: filter_semantic_occurrences


text:
@testset "Test default NCU operation" begin
    # 5 qubit test
    num_qubits = 5
    cct_s = QuantExQASM.Algorithms.create_ncu_circuit(num_qubits)

    @test begin
        result = get_statevector_using_qiskit( cct_s, big_endian=true )
        result[2^num_qubits] ≈ 1.0 + 0im
    end
    @test begin
        result = get_statevector_using_picoquant( cct_s, big_endian=true )
        result[2^num_qubits] ≈ 1.0 + 0im
    end
end

@testset "Test NCU+auxiliary (5+3+1) operation" begin
    # 5+3+1 qubit test
    num_qubits = 9
    cct_s = QuantExQASM.Algorithms.create_ncu_circuit(num_qubits, true)

    @test begin
        result = get_statevector_using_qiskit( cct_s, big_endian=true )
        result[2^Int(ceil(num_qubits/2)) + 2^(num_qubits-1)] ≈ 1.0 + 0im
    end
    @test begin
        result = get_statevector_using_picoquant( cct_s, big_endian=true )
        result[2^Int(ceil(num_qubits/2)) + 2^(num_qubits-1)] ≈ 1.0 + 0im
    end
end


text:
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H, D},
    ::Type{S},
) where {
    H <: PSY.HydroGen,
    D <: AbstractHydroDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, H, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H, HydroDispatchReservoirFlow},
    ::Type{S},
) where {H <: PSY.HydroGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_limit_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, H, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, HydroDispatchReservoirStorage, S)

    return
end

#=
# All Hydro UC formulations are currently not supported
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H,D},
    ::Type{S};

) where {H<:PSY.HydroGen,D<:AbstractHydroUnitCommitment,S<:PM.AbstractPowerModel}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S,get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S,get_feedforward(model))
    commitment_constraints!(psi_container, devices, model, S,get_feedforward(model))
    feedforward!(psi_container, H,get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S)

    return
end
=#

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H, D},
    ::Type{S},
) where {
    H <: PSY.HydroGen,
    D <: AbstractHydroDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, H, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H, HydroDispatchReservoirFlow},
    ::Type{S},
) where {H <: PSY.HydroGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_limit_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, H, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, HydroDispatchReservoirFlow, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H, HydroDispatchReservoirStorage},
    ::Type{S},
) where {H <: PSY.HydroGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    energy_variables!(psi_container, devices)
    spillage_variables!(psi_container, devices)

    #Initial Conditions
    storage_energy_init(psi_container.initial_conditions, devices)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    energy_balance_constraint!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, H, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, HydroDispatchReservoirStorage, S)

    return
end

#=
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H,D},
    ::Type{S};

) where {H<:PSY.HydroGen,D<:AbstractHydroUnitCommitment,S<:PM.AbstractActivePowerModel}

    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    commitment_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S,get_feedforward(model))
    commitment_constraints!(psi_container, devices, model, S,get_feedforward(model))
    feedforward!(psi_container, H,get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S)

    return
end
=#

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{H, FixedOutput},
    ::Type{S},
) where {H <: PSY.HydroGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(H, sys)

    if !validate_available_devices(H, devices)
        return
    end

    nodal_expression!(psi_container, devices, S)

    return
end

# Currently no Hydro device supports a Unit commiment formulation
#=
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{PSY.HydroDispatch, D},
    ::Type{S};

) where {D <: AbstractHydroUnitCommitment, S <: PM.AbstractPowerModel}
    @warn("The Formulation $(D) only applies to Dispatchable Hydro, *
               Consider Changing the Device Formulation to FixedOutput")

    construct_device!(
        psi_container,
        DeviceModel(PSY.HydroDispatch, FixedOutput),
        S

    )
end
=#

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{PSY.HydroDispatch, FixedOutput},
    ::Type{S},
) where {S <: PM.AbstractPowerModel}
    devices = get_available_components(PSY.HydroDispatch, sys)

    if !validate_available_devices(PSY.HydroDispatch, devices)
        return
    end

    nodal_expression!(psi_container, devices, S)

    return
end


text:
abstract type sCPAAlphaModel <: AlphaModel end

struct sCPAAlphaParam <: EoSParam
    c1::SingleParam{Float64}
end

@newmodelsimple sCPAAlpha sCPAAlphaModel sCPAAlphaParam

export sCPAAlpha
function sCPAAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["SAFT/CPA/sCPA/sCPA_like.csv"]; userlocations=userlocations, ignore_missing_singleparams=["Mw"], verbose=verbose)
    c1 = params["c1"]
    packagedparams = sCPAAlphaParam(c1)
    model = sCPAAlpha(packagedparams, verbose=verbose)
    return model
end

function α_function(model::CubicModel,V,T,z,alpha_model::sCPAAlphaModel)
    Tc = model.params.Tc.values
    Tr = @. T/Tc
    c1  = alpha_model.params.c1.values
    α  = @. (1+c1*(1-√(Tr)))^2
    return α
end


text:
using AOM
using Test

@testset "AOM.jl" begin
    # Write your own tests here.
end


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

import Statistics: mean

"""
    overlay(t1, t2, ...)

Construct a TSeries in which each observation is taken from the first
non-missing observation in the list of arguments. A missing observation is one
for which [`istypenan`](@ref) returns `true`.

All TSeries in the argument list must be of the same frequency. The data type of
the resulting TSeries is computed by the standard promotion of numerical types
in Julia. Its range is the union of the ranges of the arguments.
"""
@inline overlay(ts::Vararg{TSeries{F},N}) where {F <: Frequency,N} = overlay(mapreduce(rangeof, union, ts), ts...)

"""
    overlay(rng, t1, t2, ...)

If the first argument is a range (must be of the same frequency), that becomes
the range of the resulting TSeries.
"""
function overlay(rng::AbstractRange{MIT{F}}, ts::Vararg{TSeries{F},N}) where {F <: Frequency,N}
    T = mapreduce(eltype, promote_type, ts)
    ret = TSeries(rng, typenan(T))
    # na = collection of periods where the entry of ret is missing (typenan(T))
    na = collect(rng)
    for t in ts
        if isempty(na)
            # if na is empty, then we've assigned all slots
            break
        end
        # keep = periods that are not yet assigned and t has valid values in them
        keep = intersect(na, rangeof(t)[values(@. !istypenan(t))])
        # assign
        ret[keep] = t[keep]
        # update na by removing the periods we just assigned
        na = setdiff(na, keep)
    end
    return ret
end

function _valid_range(t::TSeries)
    fd = firstdate(t)
    ld = lastdate(t)
    while fd <= ld && istypenan(t[ld])
        ld -= 1
    end
    while fd <= ld && istypenan(t[fd])
        fd += 1
    end
    return fd:ld
end

Base.strip(t::TSeries) = getindex(t, _valid_range(t))
strip!(t::TSeries) = resize!(t, _valid_range(t))

"""
    fconvert(F, t)

Convert the time series `t` to the desired frequency `F`.
"""
fconvert(F::Type{<:Frequency}, t::TSeries; args...) = error("""
Conversion from $(frequencyof(t)) to $F not implemented.
""")

# do nothing when the source and target frequencies are the same.
fconvert(::Type{F}, t::TSeries{F}) where {F <: Frequency} = t

"""
    fconvert(F1, t::TSeries{F2}; method) where {F1 <: YPFrequency, F2 <: YPFrequency}

Convert between frequencies of the [`YPFrequency`](@ref) variety.

TODO: describe `method` when converting to a higher frequency (interpolation)
TODO: describe `method` when converting to a lower frequency (aggregation)

"""
function fconvert(F::Type{<:YPFrequency{N1}}, t::TSeries{<:YPFrequency{N2}}; method=nothing) where {N1,N2}
    args = Dict()
    if method !== nothing
        args[:method] = method
    end
    N1 > N2 ? _to_higher(F, t; args...) : _to_lower(F, t; args...)
end

function _to_higher(F::Type{<:YPFrequency{N1}}, t::TSeries{<:YPFrequency{N2}}; method=:const) where {N1,N2}
    (np, r) = divrem(N1, N2)
    if r != 0
        throw(ArgumentError("Cannot convert to higher frequency with $N1 ppy from $N2 ppy - not an exact multiple."))
    end
    # np = number of periods of the destination frequency for each period of the source frequency
    (y1, p1) = mit2yp(firstindex(t))
    # (y2, p2) = yp(lastindex(t))
    fi = MIT{F}(y1, (p1 - 1) * np + 1)
    # lastindex_s = pp(y2, p2*np; N=N1))
    if method == :const
        return TSeries(fi, repeat(t.values, inner=np))
    else
        throw(ArgumentError("Conversion method not available: $(method)."))
    end
end

function _to_lower(F::Type{<:YPFrequency{N1}}, t::TSeries{<:YPFrequency{N2}}; method=:mean) where {N1,N2}
    (np, r) = divrem(N2, N1)
    # println("np = $np, r = $r")
    if r != 0
        throw(ArgumentError("Cannot convert to lower frequency with $N1 from $N2 - not an exact multiple."))
    end
    (y1, p1) = mit2yp(firstindex(t))
    (d1, r1) = divrem(p1 - 1, np)
    fi = MIT{F}(y1, d1 + 1) + (r1 > 0)
    # println("y1 = $y1, p1 = $p1, d1 = $d1, r1 = $r1, fi = $fi")
    (y2, p2) = mit2yp(lastindex(t))
    (d2, r2) = divrem(p2 - 1, np)
    li = MIT{F}(y2, d2 + 1) - (r2 < np - 1)
    # println("y2 = $y2, p2 = $p2, d2 = $d2, r2 = $r2, li = $li")
    ret = TSeries(eltype(t), fi:li)
    vals = t[begin + (r1 > 0) * (np - r1):end - (r2 < np-1)*(1+r2)].values
    # println("vals = $vals")
    if method == :mean
        copyto!(ret, mean(reshape(vals, np, :); dims=1))
    elseif method == :sum
        copyto!(ret, sum(reshape(vals, np, :); dims=1))
    elseif method == :begin
        copyto!(ret, reshape(vals, np, :)[begin, :])
    elseif method == :end
        copyto!(ret, reshape(vals, np, :)[end, :])
    else
        throw(ArgumentError("Conversion method not available: $(method)."))
    end
    return ret
end


text:
function __init__svhn2()
    DEPNAME = "SVHN2"
    TRAINDATA = "train_32x32.mat"
    TESTDATA  = "test_32x32.mat"
    EXTRADATA = "extra_32x32.mat"

    register(DataDep(
        DEPNAME,
        """
        Dataset: The Street View House Numbers (SVHN) Dataset
        Authors: Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, Andrew Y. Ng
        Website: http://ufldl.stanford.edu/housenumbers
        Format: Cropped Digits (Format 2 on the website)
        Note: for non-commercial use only

        [Netzer et al., 2011]
            Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, Andrew Y. Ng
            "Reading Digits in Natural Images with Unsupervised Feature Learning"
            NIPS Workshop on Deep Learning and Unsupervised Feature Learning 2011

        The dataset is split up into three subsets: 73257
        digits for training, 26032 digits for testing, and
        531131 additional to use as extra training data.

        The files are available for download at the official
        website linked above. Note that using the data
        responsibly and respecting copyright remains your
        responsibility. For example the website mentions that
        the data is for non-commercial use only. Please read
        the website to make sure you want to download the
        dataset.
        """,
        "http://ufldl.stanford.edu/housenumbers/" .* [TRAINDATA, TESTDATA, EXTRADATA],
        "2fa3b0b79baf39de36ed7579e6947760e6241f4c52b6b406cabc44d654c13a50"
    ))
end

"""
    SVHN2(; Tx=Float32, split=:train, dir=nothing)
    SVHN2([Tx, split])

The Street View House Numbers (SVHN) Dataset.

- Authors: Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, Andrew Y. Ng
- Website: http://ufldl.stanford.edu/housenumbers

SVHN was obtained from house numbers in Google Street View
images. As such they are quite diverse in terms of orientation
and image background. Similar to MNIST, SVHN has 10 classes (the
digits 0-9), but unlike MNIST there is more data and the images
are a little bigger (32x32 instead of 28x28) with an additional
RGB color channel. The dataset is split up into three subsets:
73257 digits for training, 26032 digits for testing, and 531131
additional to use as extra training data.

# Arguments

$ARGUMENTS_SUPERVISED_ARRAY
- `split`: selects the data partition. Can take the values `:train:`, `:test` or `:extra`.

# Fields

$FIELDS_SUPERVISED_ARRAY
- `split`.

# Methods

$METHODS_SUPERVISED_ARRAY
- [`convert2image`](@ref) converts features to `RGB` images.

# Examples

```julia-repl
julia> using MLDatasets: SVHN2

julia> using MLDatasets: SVHN2

julia> dataset = SVHN2()
SVHN2:
  metadata    =>    Dict{String, Any} with 2 entries
  split       =>    :train
  features    =>    32×32×3×73257 Array{Float32, 4}
  targets     =>    73257-element Vector{Int64}

julia> dataset[1:5].targets
5-element Vector{Int64}:
 1
 9
 2
 3
 2

julia> dataset.metadata
Dict{String, Any} with 2 entries:
  "n_observations" => 73257
  "class_names"    => ["1", "2", "3", "4", "5", "6", "7", "8", "9", "0"]
```
"""
struct SVHN2 <: SupervisedDataset
    metadata::Dict{String, Any}
    split::Symbol
    features::Array{<:Any, 4}
    targets::Vector{Int}
end

SVHN2(; split=:train, Tx=Float32, dir=nothing) = SVHN2(Tx, split; dir)

function SVHN2(Tx::Type, split::Symbol=:train; dir=nothing)
    DEPNAME = "SVHN2"
    TRAINDATA = "train_32x32.mat"
    TESTDATA  = "test_32x32.mat"
    EXTRADATA = "extra_32x32.mat"
    CLASSES = [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]
    @assert split ∈ [:train, :test, :extra]
    if split == :train
        PATH = TRAINDATA
    elseif split == :test
        PATH = TESTDATA
    else
        PATH = EXTRADATA
    end

    path = datafile(DEPNAME, PATH, dir)
    vars = matread(path)
    images = vars["X"]::Array{UInt8,4}
    labels = vars["y"]
    images = permutedims(images, (2, 1, 3, 4))
    features = bytes_to_type(Tx, images)
    targets = Vector{Int}(vec(labels))

    metadata = Dict{String, Any}()
    metadata["n_observations"] = size(features)[end]
    metadata["class_names"] = string.(CLASSES)

    return SVHN2(metadata, split, features, targets)
end

convert2image(::Type{<:SVHN2}, x) = convert2image(CIFAR10, x)

# DEPRECATED INTERFACE, REMOVE IN v0.7 (or 0.6.x)
function Base.getproperty(::Type{SVHN2}, s::Symbol)
    if s == :traintensor
        @warn "SVHN2.traintensor() is deprecated, use `SVHN2(split=:train).features` instead." maxlog=2
        traintensor(T::Type=N0f8; kws...) = traintensor(T, :; kws...)
        traintensor(i; kws...) = traintensor(N0f8, i; kws...)
        function traintensor(T::Type, i; dir=nothing)
            SVHN2(; split=:train, Tx=T, dir)[i][1]
        end
        return traintensor
    elseif s == :testtensor
        @warn "SVHN2.testtensor() is deprecated, use `SVHN2(split=:test).features` instead."  maxlog=2
        testtensor(T::Type=N0f8; kws...) = testtensor(T, :; kws...)
        testtensor(i; kws...) = testtensor(N0f8, i; kws...)
        function testtensor(T::Type, i; dir=nothing)
            SVHN2(; split=:test, Tx=T, dir)[i][1]
        end
        return testtensor
    elseif s == :trainlabels
        @warn "SVHN2.trainlabels() is deprecated, use `SVHN2(split=:train).targets` instead."  maxlog=2
        trainlabels(; kws...) = trainlabels(:; kws...)
        function trainlabels(i; dir=nothing)
            SVHN2(; split=:train, dir)[i][2]
        end
        return trainlabels
    elseif s == :testlabels
        @warn "SVHN2.testlabels() is deprecated, use `SVHN2(split=:test).targets` instead." maxlog=2
        testlabels(; kws...) = testlabels(:; kws...)
        function testlabels(i; dir=nothing)
            SVHN2(; split=:test, dir)[i][2]
        end
        return testlabels
    elseif s == :traindata
        @warn "SVHN2.traindata() is deprecated, use `SVHN2(split=:train)[]` instead." maxlog=2
        traindata(T::Type=N0f8; kws...) = traindata(T, :; kws...)
        traindata(i; kws...) = traindata(N0f8, i; kws...)
        function traindata(T::Type, i; dir=nothing)
            SVHN2(; split=:train, Tx=T, dir)[i]
        end
        return traindata
    elseif s == :testdata
        @warn "SVHN2.testdata() is deprecated, use `SVHN2(split=:test)[]` instead."  maxlog=2
        testdata(T::Type=N0f8; kws...) = testdata(T, :; kws...)
        testdata(i; kws...) = testdata(N0f8, i; kws...)
        function testdata(T::Type, i; dir=nothing)
            SVHN2(; split=:test, Tx=T, dir)[i]
        end
        return testdata
    elseif s == :convert2image
        @warn "SVHN2.convert2image(x) is deprecated, use `convert2image(SVHN2, x)` instead"
        return x -> convert2image(SVHN2, x)
    else
        return getfield(SVHN2, s)
    end
end


text:
"""
    readpng(pathname)

Read a PNG file.

This returns a image object suitable for placing on the current drawing with `placeimage()`.
You can access its `width` and `height` fields:

    image = readpng("/tmp/test-image.png")
    w = image.width
    h = image.height
"""
function readpng(pathname)
    im = Cairo.read_from_png(pathname)
    if iszero(im.width) || iszero(im.height)
        throw(error("can't read this PNG image from $(pathname). Either it's not a valid PNG file, or the format is different from what Cairo is expecting."))
    end
    return im
end

function paint_with_alpha(ctx::Cairo.CairoContext, a = 0.5)
    Cairo.paint_with_alpha(get_current_cr(), a)
end

"""
    placeimage(img, xpos, ypos; centered=false)

Place a PNG image on the drawing at (`xpos`/`ypos`). The image `img` has been previously
loaded using `readpng()`.

Use keyword `centered=true` to place the center of the image at the position.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos; centered=false)
    if centered == true
        w, h = img.width, img.height
        xpos, ypos = xpos - (w/2), ypos - (h/2)
    end
    Cairo.set_source_surface(get_current_cr(), img, xpos, ypos)
    # no alpha
    Cairo.paint(get_current_cr())
end

"""
    placeimage(img, pos; centered=false)

Place the top left corner of the PNG image on the drawing at `pos`.

Use keyword `centered=true` to place the center of the image at the position.
"""
placeimage(img::Cairo.CairoSurface, pt::Point; kwargs...) = placeimage(img, pt.x, pt.y; kwargs...)

"""
    placeimage(img, xpos, ypos, a; centered=false)

Place a PNG image on the drawing at (`xpos`/`ypos`) with transparency `a`.

Use keyword `centered=true` to place the center of the image at the position.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos, alpha; centered=false)
    if centered == true
        w, h = img.width, img.height
        xpos, ypos = xpos - (w/2), ypos - (h/2)
    end
    Cairo.set_source_surface(get_current_cr(), img, xpos, ypos)
    paint_with_alpha(get_current_cr(), alpha)
end

"""
    placeimage(img, pos, a; centered=false)

Place a PNG image on the drawing at `pos` with transparency `a`.

Use keyword `centered=true` to place the center of the image at the position.
"""
placeimage(img::Cairo.CairoSurface, pt::Point, alpha; kwargs...) =
  placeimage(img::Cairo.CairoSurface, pt.x, pt.y, alpha; kwargs...)


text:
@inline signbit(a::DoubleFloat{T}) where {T} = signbit(HI(a))
@inline sign(a::DoubleFloat{T}) where {T} = sign(HI(a))

@inline function (-)(a::DoubleFloat{T}) where {T}
    if iszero(LO(a))
        DoubleFloat(-HI(a), LO(a))
    else
        DoubleFloat(-HI(a), -LO(a))
    end
end

@inline function abs(a::DoubleFloat{T}) where {T}
    if HI(a) >= 0.0
        a
    else # HI(a) < 0.0
        -a
    end
end

@inline function flipsign(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    signbit(y) ? -x : x
end
@inline function flipsign(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    signbit(y) ? -x : x
end

@inline function copysign(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    signbit(y) ? -(abs(x)) : abs(x)
end
@inline function copysign(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    signbit(y) ? -(abs(x)) : abs(x)
end

function frexp(x::DoubleFloat{T}) where {T<:AbstractFloat}
    frhi, exphi = frexp(HI(x))
    frlo, explo = frexp(LO(x))
    return (frhi, exphi), (frlo, explo)
end

function ldexp(x::DoubleFloat{T}, exponent::I) where {T, I<:Integer}
    return DoubleFloat(ldexp(HI(x), exponent), ldexp(LO(x), exponent))
end

function ldexp(dhi::Tuple{T,T}, dlo::Tuple{T,T}) where {T<:AbstractFloat}
    return DoubleFloat(ldexp(dhi[1], dhi[2]), ldexp(dlo[1], dlo[2]))
end

function exponent(x::DoubleFloat{T}) where {T<:AbstractFloat}
    ehi = exponent(HI(x))
    elo = exponent(LO(x))
    return ehi, elo
end

function significand(x::DoubleFloat{T}) where {T<:AbstractFloat}
    shi = significand(HI(x))
    slo = significand(LO(x))
    return shi, slo
end

function signs(x::DoubleFloat{T}) where {T<:AbstractFloat}
    shi = sign(HI(x))
    slo = sign(LO(x))
    return shi, slo
end

function nextfloat(x::DoubleFloat{T}) where {T<:AbstractFloat}
    !isfinite(x) && return(x)
    signbit(x) && return -prevfloat(-x)
    iszero(LO(x)) && return DoubleFloat(HI(x), eps(HI(x)))
    return DoubleFloat{T}(HI(x), nextfloat(LO(x)))
end

function prevfloat(x::DoubleFloat{T}) where {T<:AbstractFloat}
    !isfinite(x) && return(x)
    signbit(x) && return -nextfloat(-x)
    iszero(LO(x)) && return DoubleFloat(HI(x), -eps(HI(x)))
    return DoubleFloat{T}(HI(x), prevfloat(LO(x)))
end

function prevfloat(x::DoubleFloat{T}, n::Int) where {T<:AbstractFloat}
     return nextfloat(x, -n)
end

function nextfloat(x::DoubleFloat{T}, n::Int) where {T<:AbstractFloat}
   !isfinite(x) && return(x)
   signbit(x) && return -prevfloat(-x, n)
   iszero(LO(x)) && return DoubleFloat(HI(x), sign(x)*T(n)*eps(HI(x)))
   signbit(LO(x)) && return DoubleFloat(HI(x), -prevfloat(-LO(x),n))
   return DoubleFloat(HI(x), nextfloat(LO(x),n))
end

eps(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(eps(eps( one(T) ))/2)
eps(x::DoubleFloat{T}) where {T<:AbstractFloat} = abs(x) * eps(typeof(x))

ulp(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(eps(eps( one(T) ))/4)
ulp(x::DoubleFloat{T}) where {T<:AbstractFloat} = abs(x) * ulp(typeof(x))

@inline function intpart(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    ihi = modf(hi)[2]
    ilo = modf(lo)[2]
    return add_2(ihi, ilo)
end

@inline function fracpart(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    fhi = modf(hi)[1]
    flo = modf(lo)[1]
    return add_2(fhi, flo)
end

function intpart(x::DoubleFloat{T}) where {T<:AbstractFloat}
    ihi, ilo = intpart(HILO(x))
    return DoubleFloat(ihi, ilo)
end

function fracpart(x::DoubleFloat{T}) where {T<:AbstractFloat}
    fhi, flo = fracpart(HILO(x))
    return DoubleFloat(fhi, flo)
end

function modf(x::DoubleFloat{T}) where {T<:AbstractFloat}
    hi, lo = HILO(x)
    fhi, ihi = modf(hi)
    flo, ilo = modf(lo)
    ihi, ilo = add_2(ihi, ilo)
    fhi, flo = add_2(fhi, flo)
    i = DoubleFloat(ihi, ilo)
    f = DoubleFloat(fhi, flo)
    return f, i
end

function fmod(fpart::DoubleFloat{T}, ipart::DoubleFloat{T}) where {T<:AbstractFloat}
   return ipart + fpart
end

function fmod(parts::Tuple{DoubleFloat{T}, DoubleFloat{T}}) where {T<:AbstractFloat}
   return parts[1] + parts[2]
end


text:
using Intervals
using TimeZones
using Test
using Dates

@testset "Intervals" begin
    include("inclusivity.jl")
    include("endpoint.jl")
    include("interval.jl")
    include("anchoredinterval.jl")
    include("comparisons.jl")
end


text:
"""
    legendre_point(p::T) where {T<:Integer}

Calculate Legendre points of polynomial degree p
"""
legendre_point(p::T) where {T<:Integer} = gausslegendre(p + 1)[1]

"""
    ∂legendre(p::T, x) where {T<:Integer}

Calculate derivatives of Legendre polynomials of degree p at location x
"""
∂legendre(p::T, x) where {T<:Integer} = last(sf_legendre_Pl_deriv_array(p, x)[2])

function ∂legendre(p::I, x::T) where {I<:Integer,T<:AbstractArray{<:Real,1}}
    Δ = similar(x)
    for i in eachindex(Δ)
        Δ[i] = ∂legendre(p, x[i])
    end

    return Δ
end

"""
    ∂radau(p::TI, x::TU) where {TI<:Integer,TU<:Union{Real,AbstractArray{<:Real,1}}}

Calculate derivatives of Radau polynomials of degree p at location x
"""
function ∂radau(p::TI, x::TU) where {TI<:Integer,TU<:Union{Real,AbstractArray{<:Real,1}}}
    Δ = ∂legendre(p, x)
    Δ_plus = ∂legendre(p + 1, x)

    dgl = @. (-1.0)^p * 0.5 * (Δ - Δ_plus)
    dgr = @. 0.5 * (Δ + Δ_plus)

    return dgl, dgr
end


text:
##
using Plots
using FFTW
using DSP
sr = 48000
baseFreq = 20
overSamp = 2

maxHarms = floor(Int, sr / (3.0+baseFreq) + 0.5)
v = maxHarms
v -= 1
v |= v >> 1
v |= v >> 2
v |= v >> 4
v |= v >> 8
v |= v >> 16
v += 1
tableLen = v * 2 * overSamp

topFreq = baseFreq * 2 / sr

## sawtooth
nHarms = 32
a = zeros(Float64,tableLen)
for i=1:nHarms
    temp = -1.0/i
    a[i+1]=-temp
    a[tableLen-i+1]=temp
end
# plot(a)

A = fft(a)
plot([real.(A),imag.(A)])

##
shift = 1280
x = imag.(A)
xs = circshift(x,shift)
p1 = plot([x,xs])
X = ifft(xs)
p2 = plot([a[1:32],abs.(X[1:32])])
p3 = plot(unwrap(angle.(X[1:32])))
#p2 = plot([a[1:32],real.(X[1:32]),imag.(X[1:32])],ylims=(-1,1))
plot(p1,p2,p3,layout=(3,1))

##
anim = @animate for shift=range(1,tableLen,step=64)
    x = imag.(A)
    xs = circshift(x,shift)
    p1 = plot([x,xs])
    X = ifft(xs)
    p2 = plot([a[1:nHarms],abs.(X[1:nHarms])])
    p3 = plot(unwrap(angle.(X[1:nHarms])),ylim=(-nHarms*pi,nHarms*pi))
    #p2 = plot([a[1:32],real.(X[1:32]),imag.(X[1:32])],ylims=(-1,1))
    plot(p1,p2,p3,layout=(3,1),size=(600,800))
end

gif(anim,"analysis/test.gif",fps=30)


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

using Markdown
using InteractiveUtils

# ╔═╡ 96694b30-f034-11ea-1e07-f10a958db49f
md""" So far we have discussed limits in Calculus of various types, continuity, mean value theorem, etc..

In this unit we will discuss derivatives and their applications.

Given two points $(x_1,y_1), (x_2,y_2)$ we can calculate the slope $m$ between them:

$$m = \frac{y_2-y_1}{x_2-x_1}$$
"""

# ╔═╡ 6271a790-f035-11ea-116b-11e6dcef81b1
md"""We can calculate slope of linear functions, ever since we took our first Algebra course.

But what about "non-linear" fuctions?

Let's consider our friendly quadradic fuction $f(x) = x^2$.

Here's a graph:"""

# ╔═╡ 3283a680-f050-11ea-13da-6f5150506375
md"What about slope between two points on the graph?"

# ╔═╡ db616880-f052-11ea-25e3-1bc2914021fb
md"""Hmm. That doesn't look like that'll convey useful information.

What about if we look at two closer points?"""

# ╔═╡ e76835ee-f039-11ea-002c-c5eea484adf4
using Gadfly

# ╔═╡ 621a3160-f038-11ea-007d-c7b61f010b89
begin
	p = plot((z)->z^2, -1.0,1.0)
	#push!(p, layer(x=[slider], y =[slider^2], size=[.05],color=[colorant"gold"]))
	push!(p, Coord.cartesian(fixed=true))
end

# ╔═╡ 8d3b5ff0-f050-11ea-1d4f-173996efe68b
begin
	p1 = plot((z)->z^2, -1.0,1.0)
	push!(p1,layer((z)->(.416*z+.4583),-.5, .92))
	push!(p1, Coord.cartesian(fixed=true))
	push!(p1,Guide.title("Slope between (-0.5,.25) and (.92,.85)"))
end

# ╔═╡ 333c9510-f036-11ea-3e76-d1005102754d
begin
	p2 = plot((z)->z^2, -1.5,1.5)
	push!(p2,layer((z)->(1.6911*z-.7058),.5, 1.5))
	push!(p2, Coord.cartesian(fixed=true))
	push!(p2,Guide.title("Slope between (.75,.5625) and (.92,.85)"))
end

# ╔═╡ Cell order:
# ╟─e76835ee-f039-11ea-002c-c5eea484adf4
# ╟─96694b30-f034-11ea-1e07-f10a958db49f
# ╟─6271a790-f035-11ea-116b-11e6dcef81b1
# ╟─621a3160-f038-11ea-007d-c7b61f010b89
# ╟─3283a680-f050-11ea-13da-6f5150506375
# ╟─8d3b5ff0-f050-11ea-1d4f-173996efe68b
# ╟─db616880-f052-11ea-25e3-1bc2914021fb
# ╟─333c9510-f036-11ea-3e76-d1005102754d


text:
"""
    NDIndex(i, j, k...)   -> I
    NDIndex((i, j, k...)) -> I

A multidimensional index that refers to a single element. Each dimension is represented by
a single `Int` or `StaticInt`.

```julia
julia> using ArrayInterface: NDIndex

julia> using Static

julia> i = NDIndex(static(1), 2, static(3))
NDIndex(static(1), 2, static(3))

julia> i[static(1)]
static(1)

julia> i[1]
1

```
"""
struct NDIndex{N,I<:Tuple{Vararg{Any,N}}} <: AbstractCartesianIndex{N}
    index::I

    global _NDIndex(index::Tuple{Vararg{Any,N}}) where {N} = new{N,typeof(index)}(index)

    function NDIndex{N,I}(index::I) where {N,I<:Tuple{Vararg{Integer,N}}}
        for i in index
            (i isa Int) || i isa StaticInt || throw(ArgumentError("NDIndex does not support values of type $(typeof(i))"))
        end
        return new{N,I}(index)
    end

    NDIndex{N}(index::Tuple) where {N} = _ndindex(static(N), _flatten(index...))
    NDIndex{N}(index...) where {N} = _ndindex(static(N), _flatten(index...))

    NDIndex{0}(::Tuple{}) = new{0,Tuple{}}(())
    NDIndex{0}() = NDIndex{0}(())

    NDIndex(index::Tuple) = _NDIndex(_flatten(index...))
    NDIndex(index...) = _NDIndex(_flatten(index...))
end

_ndindex(n::StaticInt{N}, i::Tuple{Vararg{Union{Int,StaticInt},N}}) where {N} = _NDIndex(i)
function _ndindex(n::StaticInt{N}, i::Tuple{Vararg{Any,M}}) where {N,M}
    M > N && throw(ArgumentError("input tuple of length $M, requested $N"))
    return _NDIndex(_fill_to_length(i, n))
end
_fill_to_length(x::Tuple{Vararg{Any,N}}, n::StaticInt{N}) where {N} = x
@inline function _fill_to_length(x::Tuple{Vararg{Any,M}}, n::StaticInt{N}) where {M,N}
    return _fill_to_length((x..., static(1)), n)
end

_flatten(i::StaticInt{N}) where {N} = (i,)
_flatten(i::Integer) = (Int(i),)
_flatten(i::Base.AbstractCartesianIndex) = _flatten(Tuple(i)...)
@inline _flatten(i::Integer, I...) = (canonicalize(i), _flatten(I...)...)
@inline function _flatten(i::Base.AbstractCartesianIndex, I...)
    return (_flatten(Tuple(i)...)..., _flatten(I...)...)
 end
Base.Tuple(index::NDIndex) = index.index

Static.dynamic(x::NDIndex) = CartesianIndex(dynamic(Tuple(x)))
Static.static(x::CartesianIndex) = _NDIndex(static(Tuple(x)))
Static.known(::Type{NDIndex{N,I}}) where {N,I} = known(I)

Base.show(io::IO, i::NDIndex) = (print(io, "NDIndex"); show(io, Tuple(i)))

# length
Base.length(::NDIndex{N}) where {N} = N
Base.length(::Type{NDIndex{N,I}}) where {N,I} = N
known_length(::Type{NDIndex{N,I}}) where {N,I} = N

# indexing
@propagate_inbounds function getindex(x::NDIndex{N,T}, i::Int)::Int where {N,T}
    return Int(getfield(Tuple(x), i))
end
@propagate_inbounds function getindex(x::NDIndex{N,T}, i::StaticInt{I}) where {N,T,I}
    return getfield(Tuple(x), I)
end
@propagate_inbounds Base.getindex(x::NDIndex, i::Integer) = ArrayInterface.getindex(x, i)

# Base.get(A::AbstractArray, I::CartesianIndex, default) = get(A, I.I, default)
# eltype(::Type{T}) where {T<:CartesianIndex} = eltype(fieldtype(T, :I))

Base.setindex(x::NDIndex, i, j) = NDIndex(Base.setindex(Tuple(x), i, j))

# equality
Base.:(==)(x::NDIndex{N}, y::NDIndex{N}) where N = Tuple(x) == Tuple(y)

# zeros and ones
Base.zero(::NDIndex{N}) where {N} = zero(NDIndex{N})
Base.zero(::Type{NDIndex{N}}) where {N} = _NDIndex(ntuple(_ -> static(0), Val(N)))
Base.oneunit(::NDIndex{N}) where {N} = oneunit(NDIndex{N})
Base.oneunit(::Type{NDIndex{N}}) where {N} = _NDIndex(ntuple(_ -> static(1), Val(N)))

@inline function Base.split(i::NDIndex, V::Val)
    i, j = split(Tuple(i), V)
    return NDIndex(i), NDIndex(j)
end

# arithmetic, min/max
@inline Base.:(-)(i::NDIndex{N}) where {N} = NDIndex{N}(map(-, Tuple(i)))
@inline function Base.:(+)(i1::NDIndex{N}, i2::NDIndex{N}) where {N}
    return _NDIndex(map(+, Tuple(i1), Tuple(i2)))
end
@inline function Base.:(-)(i1::NDIndex{N}, i2::NDIndex{N}) where {N}
    return _NDIndex(map(-, Tuple(i1), Tuple(i2)))
end
@inline function Base.min(i1::NDIndex{N}, i2::NDIndex{N}) where {N}
    return _NDIndex(map(min, Tuple(i1), Tuple(i2)))
end
@inline function Base.max(i1::NDIndex{N}, i2::NDIndex{N}) where {N}
    return _NDIndex(map(max, Tuple(i1), Tuple(i2)))
end
@inline Base.:(*)(a::Integer, i::NDIndex{N}) where {N} = _NDIndex(map(x->a*x, Tuple(i)))
@inline Base.:(*)(i::NDIndex, a::Integer) = *(a, i)

Base.CartesianIndex(x::NDIndex) = CartesianIndex(Tuple(x))

# comparison
@inline function Base.isless(x::NDIndex{N}, y::NDIndex{N}) where {N}
    return Bool(_isless(static(0), Tuple(x), Tuple(y)))
end

Static.lt(x::NDIndex{N}, y::NDIndex{N}) where {N} = _isless(static(0), Tuple(x), Tuple(y))

_final_isless(c::Int) = c === 1
_final_isless(::StaticInt{N}) where {N} = static(false)
_final_isless(::StaticInt{1}) = static(true)
_isless(c::C, x::Tuple{}, y::Tuple{}) where {C} = _final_isless(c)
function _isless(c::C, x::Tuple, y::Tuple) where {C}
    return _isless(icmp(c, x, y), Base.front(x), Base.front(y))
end
icmp(::StaticInt{0}, x::Tuple, y::Tuple) = icmp(last(x), last(y))
icmp(::StaticInt{N}, x::Tuple, y::Tuple) where {N} = static(N)
function icmp(cmp::Int, x::Tuple, y::Tuple)
    if cmp === 0
        return icmp(Int(last(x)), Int(last(y)))
    else
        return cmp
    end
end
icmp(a, b) = _icmp(lt(a, b), a, b)
_icmp(::True, a, b) = static(1)
_icmp(::False, a, b) = __icmp(Static.eq(a, b))
function _icmp(x::Bool, a, b)
    if x
        return 1
    else
        return __icmp(a == b)
    end
end
__icmp(::True) = static(0)
__icmp(::False) = static(-1)
function __icmp(x::Bool)
    if x
        return 0
    else
        return -1
    end
end

#  Necessary for compatibility with Base
# In simple cases, we know that we don't need to use axes(A). Optimize those
# until Julia gets smart enough to elide the call on its own:
@inline function Base.to_indices(A, inds, I::Tuple{NDIndex, Vararg{Any}})
    return Base.to_indices(A, inds, (Tuple(I[1])..., tail(I)...))
end
# But for arrays of CartesianIndex, we just skip the appropriate number of inds
@inline function Base.to_indices(A, inds, I::Tuple{AbstractArray{NDIndex{N}}, Vararg{Any}}) where N
    _, indstail = IteratorsMD.split(inds, Val(N))
    return (Base.to_index(A, I[1]), Base.to_indices(A, indstail, tail(I))...)
end


text:
module MyExample

greet() = print("Hello World!")
include("extra_file.jl")

display(my_f(2,1))
end # module


text:
mutable struct NSFVecCoefs
    A::Matrix{Float64}
    B::Matrix{Float64}
    C::Matrix{Float64}
    T::Union{Matrix{Float64},Nothing}
    NSFVecCoefs() = new()
    NSFVecCoefs(A,B,C,T) = new(A,B,C,T)
end

mutable struct NSFVector
    size::Int
    moments::Int
    scales::Int
    coefs::Vector{NSFVecCoefs}
end

function NSFVector_decomp(size::Int,
        moments::Int,
        scales::Int,
        X::Matrix{Float64}
        )::NSFVector

    if size%(1<<scales) != 0
        throw("size not divisible by 2^scales")
    end

    H,G = get_qmf(moments)
    Y = NSFVector(size,moments,scales,Vector{NSFVecCoefs}(undef,scales))
    T = X
    for j in 1:scales
        A,B,C,T = dwtmat(size,T,moments,H,G)
        Y.coefs[j] = NSFVecCoefs(A,B,C,nothing)
        size ÷= 2
    end
    Y.coefs[scales].T = T

    return Y
end

function NSFVector_reconst(X::NSFVector)::Matrix{Float64}
    n = X.size÷(2^X.scales)
    H,G = get_qmf(X.moments)
    T = X.coefs[X.scales].T
    for j in X.scales:-1:1
        T = idwtmat(n,X.moments,H,G,X.coefs[j].A,X.coefs[j].B,X.coefs[j].C,T)
        n *= 2
    end
    return T
end

max_stars_repo_name: xinshuoyang/FactNSF

max_issues_repo_name: xinshuoyang/FactNSF

max_forks_repo_name: xinshuoyang/FactNSF


text:
"""
# GoldenSection
## Constructor
```julia
    GoldenSection(;)
```

## Description
The `GoldenSection` method seeks to minimize a univariate function on an interval
`[a, b]`. At all times the algorithm maintains a tuple of three minimizer candidates
`(c, d, e)` where ``c<d<e`` such that the ratio of the largest to the smallest interval
is the Golden Ratio.

## References
https://en.wikipedia.org/wiki/Golden-section_search
"""
struct GoldenSection <: UnivariateOptimizer end

Base.summary(::GoldenSection) = "Golden Section Search"

function optimize(f, x_lower::T, x_upper::T,
     mo::GoldenSection;
     rel_tol::T = sqrt(eps(T)),
     abs_tol::T = eps(T),
     iterations::Integer = 1_000,
     store_trace::Bool = false,
     show_trace::Bool = false,
     callback = nothing,
     show_every = 1,
     extended_trace::Bool = false,
     nargs...) where T <: AbstractFloat
    if x_lower > x_upper
        error("x_lower must be less than x_upper")
    end
    t0 = time()
    options = (store_trace=store_trace, show_trace=show_trace, show_every=show_every, callback=callback)
    # Save for later
    initial_lower = x_lower
    initial_upper = x_upper

    golden_ratio::T = 0.5 * (3.0 - sqrt(5.0))

    new_minimizer = x_lower + golden_ratio*(x_upper-x_lower)
    new_minimum = f(new_minimizer)
    best_bound = "initial"
    f_calls = 1 # Number of calls to f

    iteration = 0
    converged = false

    # Trace the history of states visited
    tr = OptimizationTrace{T, typeof(mo)}()
    tracing = store_trace || show_trace || extended_trace || callback !== nothing
    stopped_by_callback = false
    if tracing
        # update trace; callbacks can stop routine early by returning true
        state = (new_minimizer=new_minimizer,
                 x_lower=x_lower,
                 x_upper=x_upper,
                 best_bound=best_bound,
                 new_minimum=new_minimum)
        stopped_by_callback = trace!(tr, nothing, state, iteration, mo, options, time()-t0)
    end

    while iteration < iterations && !stopped_by_callback

        x_tol = rel_tol * abs(new_minimizer) + abs_tol

        x_midpoint = (x_upper+x_lower)/2

        if abs(new_minimizer - x_midpoint) <= 2*x_tol - (x_upper-x_lower)/2
            converged = true
            break
        end

        iteration += 1

        if x_upper - new_minimizer > new_minimizer - x_lower
            new_x = new_minimizer + golden_ratio*(x_upper - new_minimizer)
            new_f = f(new_x)
            f_calls += 1
            if new_f < new_minimum
                x_lower = new_minimizer
                best_bound = "lower"
                new_minimizer = new_x
                new_minimum = new_f
            else
                x_upper = new_x
                best_bound = "upper"
            end
        else
            new_x = new_minimizer - golden_ratio*(new_minimizer - x_lower)
            new_f = f(new_x)
            f_calls += 1
            if new_f < new_minimum
                x_upper = new_minimizer
                best_bound = "upper"
                new_minimizer = new_x
                new_minimum = new_f
            else
                x_lower = new_x
                best_bound = "lower"
            end
        end

        if tracing
            # update trace; callbacks can stop routine early by returning true
            state = (new_minimizer=new_minimizer,
                     x_lower=x_lower,
                     x_upper=x_upper,
                     best_bound=best_bound,
                     new_minimum=new_minimum)
            stopped_by_callback = trace!(tr, nothing, state, iteration, mo, options, time()-t0)
        end
    end

    return UnivariateOptimizationResults(mo,
                                         initial_lower,
                                         initial_upper,
                                         new_minimizer,
                                         new_minimum,
                                         iteration,
                                         iteration == iterations,
                                         converged,
                                         rel_tol,
                                         abs_tol,
                                         tr,
                                         f_calls)
end

function trace!(tr, d, state, iteration, method::GoldenSection, options, curr_time=time())
    dt = Dict()
    dt["time"] = curr_time
    dt["minimizer"] = state.new_minimizer
    dt["x_lower"] = state.x_lower
    dt["x_upper"] = state.x_upper
    T = eltype(state.new_minimum)

    update!(tr,
            iteration,
            state.new_minimum,
            T(NaN),
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end


text:
abstract type RunMode end

function setup end
function num_workers end
function array_initializer end
function start_computation end

_array_constructor(type, dims...) = fill(convert(type, NaN), dims)
_sharedarray_constructor(type, dims...) = fill!(SharedArray{type}(dims), NaN)

struct Serial <: RunMode end

setup(::Serial) = nothing
num_workers(::Serial) = 1
array_initializer(::Serial) = _array_constructor

start_computation(f, ::Serial, parameter_chunks) = f(parameter_chunks[1])

struct Threaded <: RunMode
    Nthreads::Int
    function Threaded(Nthreads=Threads.nthreads())
        if Threads.nthreads() < Nthreads
            logmsg("Runmode Threaded: $Nthreads requested but Julia has only $(Threads.nthreads())! Continuing with lower thread count.")
            Nthreads = Threads.nthreads()
        end
        new(Nthreads)
    end
end

setup(::Threaded) = nothing
num_workers(t::Threaded) = t.Nthreads
array_initializer(::Threaded) = _array_constructor

function start_computation(f, ::Threaded, parameter_chunks)
    @sync for chunk in parameter_chunks
        Threads.@spawn begin
            logmsg("Range: $(chunk) on thread #$(Threads.threadid())")
            f(chunk)
        end
    end
end

struct Parallel <: RunMode
    Nprocs::Int
    dosetup::Bool
    NBLASthreads::Int
    useMKL::Bool
    function Parallel(Nprocs, dosetup=true, NBLASthreads=1, useMKL=false)
        new(Nprocs, dosetup, NBLASthreads)
    end
end
function Parallel(; Nprocs=length(Sys.cpu_info()), dosetup=true, NBLASthreads=1, useMKL=false)
    Parallel(Nprocs, dosetup, NBLASthreads, useMKL)
end

setup(p::Parallel) = p.dosetup && _initialize_procs(p.Nprocs, p.NBLASthreads, p.useMKL)
num_workers(p::Parallel) = p.Nprocs
array_initializer(::Parallel) = _sharedarray_constructor

function _initialize_procs(total, num_BLAS_threads, useMKL)
    logmsg("Initialzing worker processes")
    to_add = total
    if workers() != [1]
        to_add = total-length(workers())
    end
    if to_add > 0
        logmsg("Adding $to_add new processes")
        addprocs(to_add; topology=:master_worker)
    end
    @everywhere @eval begin
        import Pkg
        Pkg.activate(".")
        using SimLib, LinearAlgebra
        if $useMKL
            using MKL
            logmsg("BLAS config: $(LinearAlgebra.BLAS.get_config())")
        end
        let blas_threads = $num_BLAS_threads
            logmsg("BLAS threads = $blas_threads")
            LinearAlgebra.BLAS.set_num_threads(blas_threads)
        end
    end
end

function start_computation(f, ::Parallel,  parameter_chunks)
    @sync for (p, chunk) in zip(workers(), parameter_chunks)
        @async remotecall_wait(p) do
            logmsg("Range: $(chunk) on process #$(myid())")
            f(chunk)
        end
    end
end


text:
"""
DimensionalData.jl provides types and methods for indexing with named dimensions,
using named dimensions in Base and Statistics methods instead of Integer dims,
and selecting data from dimension values instead of using indices directly.

## Dimensions

Dimensions are simply types that wrap values. They both store dimension values
and are used for dimension lookup or indices, ranges or dimension number.
`X`, `Y`, `Z` and `Time` are the unexported defaults, add this line to use them:
```julia
using DimensionalData: X, Y, Z, Time
```

A generalised [`Dim`](@ref) type is available to use arbitrary symbols to name dimensions.
Custom dimensions can be defined using the [`@dim`](@ref) macro.

We can use dim wrappers for indexing, so that the dimension order in the underlying array
does not need to be known:

```
a[X(1:10), Y(1:4)]
```

## Selectors

Selectors find indices in the dimension based on values `At`, `Near`, or `Between`
the index value(s).

We can use selectors in conjuction with dim wrappers:

```julia
a[X(1:10), Y<|At(25.7)]
```

Without dim wrappers selectors must be in the right order:

```julia
usin Unitful
a[Near(23u"s"), Between(10.5u"m", 50.5u"m")]
```

Dim types or objects can be used instead of a dimension number in many
Base and Statistics methods:

```julia
mean(a; dims=X)
std(a; dims=Y())
```

"""
module DimensionalData

using ConstructionBase, LinearAlgebra, RecipesBase, Statistics, Dates

using Base: tail, OneTo

export AbstractDimension, Dim

export Selector, Near, Between, At

export Locus, Center, Start, End, UnknownLocus

export Sampling, SingleSample, MultiSample, UnknownSampling

export Order, Ordered, Unordered

export Grid, UnknownGrid

export AbstractCategoricalGrid, CategoricalGrid

export IndependentGrid, AbstractAlignedGrid, AlignedGrid, BoundedGrid, RegularGrid

export DependentGrid, TransformedGrid, LookupGrid

export AbstractDimensionalArray, DimensionalArray

export Metadata, AbstractArrayMetadata, ArrayMetadata, AbstractDimMetadata, DimMetadata

export dims, refdims, metadata, name, shortname,
       val, label, units, order, bounds, locus, grid, <|

include("interface.jl")
include("metadata.jl")
include("order.jl")
include("grid.jl")
include("dimension.jl")
include("selector.jl")
include("array.jl")
include("methods.jl")
include("primitives.jl")
include("utils.jl")
include("plotrecipes.jl")

end