text:
using ADCME
using PyPlot
using ADCMEKit

n = 200
pml = 30
C = 1000.0
NT = 1000
Δt = 1.0/NT
x0 = LinRange(0, 1, n+1)
h = 1/n
xE = Array((0:n+2pml)*h .- pml*h)
xH = (xE[2:end]+xE[1:end-1])/2
N = n + 2pml + 1

σE = zeros(N)
for i = 1:pml
    d = i*h
    σE[pml + n + 1 + i] = C* (d/(pml*h))^3
    σE[pml+1-i] = C* (d/(pml*h))^3
end

σH = zeros(N-1)
for i = 1:pml
    d = (i-1/2)*h
    σH[pml + n + i] = C* (d/(pml*h))^3
    σH[pml+1-i] = C* (d/(pml*h))^3
end

function ricker(dt = 0.002, f0 = 5.0)
    nw = 2/(f0*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    b = (π*f0*t).^2
    w = @. (1 - 2b)*exp(-b)
end
R = ricker()
if length(R)<NT+1
    R = [R;zeros(NT+1-length(R))]
end
R = R[1:NT+1]
# tn = ( 0:NT ) *Δt
# R = @. exp( -20(tn-0.3)^2)
# error()
R_ = constant(R)

cH = ones(length(xH)) * 2.
cH[pml + Int(round(0.5÷h)): pml + Int(round(0.75÷h))] .= 1.0
cH[pml + Int(round(0.75÷h)): end] .= 1.
cE = (cH[1:end-1]+cH[2:end])/2
Z = zeros(N)
Z[pml + pml÷2] = 1.0
Z = Z[2:end-1]

function condition(i, E_arr, H_arr)
    i<=NT+1
end

function body(i, E_arr, H_arr)
    E = read(E_arr, i-1)
    H = read(H_arr, i-1)
    ΔH = cH * (E[2:end]-E[1:end-1])/h - σH*H
    H += ΔH * Δt
    ΔE = cE * (H[2:end]-H[1:end-1])/h - σE[2:end-1]*E[2:end-1] + R_[i] * Z
    E = scatter_add(E, 2:N-1, ΔE * Δt)
    # E = scatter_update(E, N÷2, R_[i])
    i+1, write(E_arr, i, E), write(H_arr, i, H)
end

E_arr = TensorArray(NT+1)
H_arr = TensorArray(NT+1)

E_arr = write(E_arr, 1, zeros(N))
H_arr = write(H_arr, 1, zeros(N-1))

i = constant(2, dtype = Int32)

_, E, H = while_loop(condition, body, [i, E_arr, H_arr])

E = stack(E)
H = stack(H)

sess = Session(); init(sess)
E_, H_ = run(sess, [E, H])

pl, = plot([], [], ".-")
xlim(-0.5,1.5)
ylim(minimum(E_), maximum(E_))
xlabel("x")
ylabel("y")
t = title("time = 0.0000")
function update(i)
    t.set_text("time = $(round(i*Δt, digits=4))")
    pl.set_data([xE E_[i,:]]'|>Array)
end
p = animate(update, 1:10:NT+1)
# saveanim(p, "fdtd.gif")


text:
################################################################################
# Dynamics Residual
################################################################################

function dynamics_residual(model::AbstractGameModel, pdtraj::PrimalDualTraj{KN}, k::Int) where {KN}
    zk1 = pdtraj.pr[k+1]::KN
    xk1 = state(zk1)
    zk = pdtraj.pr[k]::KN
    return dynamics_residual(model, zk, xk1)
end

function dynamics_residual(model::AbstractGameModel, zk::AbstractKnotPoint, xk1::SVx) where {SVx}
    return discrete_dynamics(RK2, model, zk) - xk1
end

################################################################################
# Dynamics Jacobian
################################################################################

function ∇dynamics!(∇f::SM, model::AbstractGameModel, pdtraj::PrimalDualTraj{KN}, k::Int) where {KN,SM}
    zk = pdtraj.pr[k]::KN
    return ∇dynamics!(∇f, model, zk)
end

function ∇dynamics!(∇f::SM, model::AbstractGameModel, zk::AbstractKnotPoint) where {SM}
    return RobotDynamics.discrete_jacobian!(RK2, ∇f, model, zk)
end

ext: jl


text:
module TestPadding
	using Base.Test
	using DataArrays

	dv = @data ones(3)
	@assert isequal(dv, padNA(dv, 0, 0))
	@assert length(padNA(dv, 2, 0)) == length(dv) + 2
	@assert length(padNA(dv, 0, 2)) == length(dv) + 2
	@assert length(padNA(dv, 2, 2)) == length(dv) + 4
end


text:
export CGPFunctions, ImgType

module CGPFunctions

using Images
using ImageMorphology
using ImageSegmentation
using LinearAlgebra

global arity = Dict()

SorX = Union{Symbol, Expr}
ImgType = Union{Array{UInt8,2}, Matrix{UInt8}}

function scaled(x::Float64)
    if isnan(x)
        return 0.0
    end
    min(max(x, -1.0), 1.0)
end

function fgen(name::Symbol, ar::Int, s1::SorX, iotype::U; safe::Bool=false) where {U <: Union{DataType, Union}}
    if safe
        @eval function $name(x::T, y::T, p::Array{Float64}=Float64[])::T where {T <: $(iotype)}
            try
                return $s1
            catch
                return x
            end
        end
    else
        @eval function $name(x::T, y::T, p::Array{Float64}=Float64[])::T where {T <: $(iotype)}
            $s1
        end
    end
    arity[String(name)] = ar
end

# Image processing

"""
    rescale_uint_img(x::AbstractArray)::Array{UInt8}

Rescale input array in [0, 255] and convert data to UInt8.
"""
function rescale_uint_img(x::AbstractArray)::Array{UInt8}
    mini, maxi = minimum(x), maximum(x)
    if mini == maxi
        return floor.(UInt8, mod.(x, 255))
    end
    m = (convert(Array{Float64}, x) .- mini) .* (255 / (maxi - mini))
    floor.(UInt8, m)
end

function remove_details(x::ImgType, p::Float64)::ImgType
    n_passes = ceil(Int64, 5 * p)
    remove_details(x, n_passes)
end

function remove_details(x::ImgType, n_passes::Int64)::ImgType
    for i in 1:n_passes
        x = ImageMorphology.erode(x)
    end
    for i in 1:n_passes
        x = ImageMorphology.dilate(x)
    end
    x
end

function felzenszwalb_segmentation(x::ImgType, p::Float64)::ImgType
    min_size = floor(Int64, 10 * p)
    segments = ImageSegmentation.felzenszwalb(x, 50, min_size)
    rescale_uint_img(segments.image_indexmap)
end

function components_segmentation(x::ImgType)::ImgType
    label = label_components(x)
    m = rescale_uint_img(label)
    remove_details(m, 1) # more passes?
end

function make_boxes(x::AbstractArray)::ImgType
    labels = label_components(x)
    boxes = component_boxes(labels)
    # Filter error boxes reaching typemax(Int64) or typemin(Int64)
    #=
    to_remove = Int64[]
    maxi = maximum(size(x))
    for i in eachindex(boxes)
        if any(boxes[i][1] .< 0) || any(boxes[i][2] .< 0) || any(boxes[i][1] .> maxi) || any(boxes[i][2] .> maxi)
            push!(to_remove, i)
        end
    end
    deleteat!(boxes, to_remove)
    =#
    # Compute areas and sort boxes
    areas = zeros(length(boxes))
    for i in eachindex(boxes)
        areas[i] = abs(boxes[i][1][1] - boxes[i][2][1]) * abs(boxes[i][1][2] - boxes[i][2][2])
    end
    boxes = reverse(boxes[sortperm(areas)])
    boxes_img = zeros(Int64, size(x))
    incr = 1
    for i in 2:length(boxes)
        boxes_img[
            boxes[i][1][1]:boxes[i][2][1],
            boxes[i][1][2]:boxes[i][2][2]
        ] .+= incr
        incr += 1
    end
    rescale_uint_img(boxes_img)
end

function box_segmentation(x::ImgType)::ImgType
    m = remove_details(x, 1) # more passes?
    label = label_components(m)
    make_boxes(label)
end

#=
function binary(x::ImgType, t::UInt8)::ImgType
    out = copy(x)
    out[x .< t] .= 0x00
    out[x .>= t] .= 0xff
    out
end

function binary(x::ImgType, t::Int64)::ImgType
    binary(x, convert(UInt8, t))
end

function binary(x::ImgType, p::Float64)::ImgType
    t = floor(UInt8, 255 * p)
    binary(x, t)
end
=#

function binary(x::ImgType, p::Float64)::ImgType
    reshape([x[i] < floor(UInt8, 255 * p) ? 0x00 : 0xff for i in eachindex(x)], size(x))
end

#=
function threshold(x::ImgType, t::UInt8, reverse::Bool=false)::ImgType
    out = copy(x)
    if reverse
        out[x .> t] .= 0x00
    else
        out[x .< t] .= 0x00
    end
    out
end

function threshold(x::ImgType, t::Int64)::ImgType
    threshold(x, convert(UInt8, t))
end

function threshold(x::ImgType, p::Float64)::ImgType
    if p < 0.5
        t = floor(UInt8, 255 * 2 * p)
        threshold(x, t, false)
    else
        t = floor(UInt8, 255 * 2 * (p - 0.5))
        threshold(x, t, true)
    end
end
=#

function threshold(x::ImgType, t::UInt8, reverse::Bool=false)::ImgType
    reshape([reverse ? (x[i]>t ? 0x00 : x[i]) : (x[i]<t ? 0x00 : x[i]) for i in eachindex(x)], size(x))
end

function threshold(x::ImgType, p::Float64)::ImgType
    if p < 0.5
        t = floor(UInt8, 255 * 2 * p)
        threshold(x, t, false)
    else
        t = floor(UInt8, 255 * 2 * (p - 0.5))
        threshold(x, t, true)
    end
end

function index1d_to_index2d(x::AbstractArray, index::Int64)::Tuple{Int64,Int64}
    ci = CartesianIndices(size(x))
    return ci[index][1], ci[index][2]
end

function motion_capture!(x::ImgType, p::Array{Float64})::ImgType
    if length(p) == length(x)
        out = x .- convert(Array{UInt8}, reshape(p, size(x)))
    else
        out = x
        resize!(p, length(x))
    end
    p[:] = x[:]
    return out
end

function motion_distances!(x::ImgType, p::Array{Float64})::ImgType
    if length(p) >= length(x)
        diff = x .- convert(Array{UInt8}, reshape(p[end-length(x)+1:end], size(x)))
        dilated = dilate(diff)
        binarized = binary(dilated, 3)
        labels = label_components(binarized)
        centroids = component_centroids(labels)
        popfirst!(centroids) # remove background centroid
        # boxes = component_boxes(labels)
        # Find reference point
        index1d = max(convert(Int64, ceil(p[1] * length(x))), 1)
        r, c = index1d_to_index2d(x, index1d)
        # Compute distances
        n_centroids = length(centroids)
        distances = zeros(n_centroids)
        for i in 1:n_centroids
            distances[i] = norm(centroids[i] .- (r, c))
        end
        distances_img = zeros(size(x))
        for i in eachindex(distances_img)
            if labels[i] > 0
                distances_img[i] = distances[labels[i]]
            end
        end
        out = rescale_uint_img(distances_img)
    else
        out = x
        resize!(p, length(p) + length(x))
    end
    p[end-length(x)+1:end] = x[:] # save first parameters
    return out
end

fgen(:f_dilate, 1, :(ImageMorphology.dilate(x)), ImgType)
fgen(:f_erode, 1, :(ImageMorphology.erode(x)), ImgType)
fgen(:f_subtract, 2, :(x .- y), ImgType)
fgen(:f_remove_details, 1, :(remove_details(x, p[1])), ImgType)
fgen(:f_make_boxes, 1, :(make_boxes(x)), ImgType)
# Segmentation
fgen(:f_felzenszwalb_segmentation, 1, :(felzenszwalb_segmentation(x, p[1])), ImgType)
fgen(:f_components_segmentation, 1, :(components_segmentation(x)), ImgType)
fgen(:f_box_segmentation, 1, :(box_segmentation(x)), ImgType)
# Thresholds and values
# fgen(:f_negative, 1, :(0xff .- x), ImgType) # Same as f_bitwise_not
fgen(:f_threshold, 1, :(threshold(x, p[1])), ImgType)
fgen(:f_threshold_bis, 1, :(threshold_bis(x, p[1])), ImgType)
fgen(:f_binary, 1, :(binary(x, p[1])), ImgType)
fgen(:f_motion_capture, 1, :(motion_capture!(x, p)), ImgType)
fgen(:f_motion_distances, 1, :(motion_distances!(x, p)), ImgType)
# Filtering
fgen(:f_corners, 1, :(rescale_uint_img(Images.fastcorners(x))), ImgType)
fgen(:f_gaussian, 1, :(rescale_uint_img(ImageFiltering.imfilter(x, Images.Kernel.gaussian(ceil(Int64, 5 * p[1]))))), ImgType)
fgen(:f_laplacian, 1, :(rescale_uint_img(ImageFiltering.imfilter(x, Images.Kernel.Laplacian()))), ImgType)
fgen(:f_sobel_x, 1, :(rescale_uint_img(ImageFiltering.imfilter(x, Images.Kernel.sobel()[2]))), ImgType)
fgen(:f_sobel_y, 1, :(rescale_uint_img(ImageFiltering.imfilter(x, Images.Kernel.sobel()[1]))), ImgType)
fgen(:f_canny, 1, :(rescale_uint_img(Images.canny(x, (Images.Percentile(80), Images.Percentile(20))))), ImgType)
fgen(:f_edges, 1, :(rescale_uint_img(Images.imedge(x)[3])), ImgType)
# Morphological functions
fgen(:f_opening, 1, :(ImageMorphology.opening(x)), ImgType)
fgen(:f_closing, 1, :(ImageMorphology.closing(x)), ImgType)
fgen(:f_tophat, 1, :(ImageMorphology.tophat(x)), ImgType)
fgen(:f_bothat, 1, :(ImageMorphology.bothat(x)), ImgType)
fgen(:f_morphogradient, 1, :(ImageMorphology.morphogradient(x)), ImgType)
fgen(:f_morpholaplace, 1, :(rescale_uint_img(ImageMorphology.morpholaplace(x))), ImgType)
# Bitwise operations
fgen(:f_bitwise_not, 1, :(.~x), ImgType)
fgen(:f_bitwise_and, 2, :(0xff .* (x .& y)), ImgType)
fgen(:f_bitwise_or, 2, :(0xff .* (x .| y)), ImgType)
fgen(:f_bitwise_xor, 2, :(0xff .* (x .⊻ y)), ImgType)
# fgen(:f_compare_eq_img, 2, :(OpenCV.compare(x, y, OpenCV.CMP_EQ)), OpenCV.InputArray)
# fgen(:f_compare_ge_img, 2, :(OpenCV.compare(x, y, OpenCV.CMP_GE)), OpenCV.InputArray)
# fgen(:f_compare_le_img, 2, :(OpenCV.compare(x, y, OpenCV.CMP_LE)), OpenCV.InputArray)

# Mathematical

fgen(:f_add, 2, :((x + y) / 2.0), Float64)
fgen(:f_subtract, 2, :(abs(x - y) / 2.0), Float64)
fgen(:f_mult, 2, :(x * y), Float64)
fgen(:f_div, 2, :(scaled(x / y)), Float64)
fgen(:f_abs, 1, :(abs(x)), Float64)
fgen(:f_sqrt, 1, :(sqrt(abs(x))), Float64)
fgen(:f_pow, 2, :(abs(x) ^ abs(y)), Float64)
fgen(:f_exp, 1, :((2 * (exp(x+1)-1.0))/(exp(2.0)-1.0) -1), Float64)
fgen(:f_sin, 1, :(sin(x)), Float64)
fgen(:f_cos, 1, :(cos(x)), Float64)
fgen(:f_tanh, 1, :(tanh(x)), Float64)
fgen(:f_sqrt_xy, 2, :(sqrt(x^2 + y^2) / sqrt(2.0)), Float64)
fgen(:f_lt, 2, :(Float64(x < y)), Float64)
fgen(:f_gt, 2, :(Float64(x > y)), Float64)

# Logical
fgen(:f_and, 2, :(Float64((&)(Int(round(x)), Int(round(y))))), Float64)
fgen(:f_or, 2, :(Float64((|)(Int(round(x)), Int(round(y))))), Float64)
fgen(:f_xor, 2, :(Float64(xor(Int(abs(round(x))), Int(abs(round(y)))))), Float64)
fgen(:f_not, 1, :(1 - abs(round(x))), Float64)

end


text:
@doc raw"""
    AcousticPerturbationEquations2D(v_mean, c_mean, rho_mean)

!!! warning "Experimental code"
    This system of equations is experimental and may change in any future release.

Acoustic perturbation equations in two space dimensions with constant mean flow. `v_mean` is a tuple
of size two holding the mean velocities in each direction. `rho_mean` and `c_mean` are scalars
representing the mean density and mean speed of sound, respectively.

The equations are based on the APE-4 system introduced in the following paper:

R. Ewert, W. Schröder
"Acoustic perturbation equations based on flow decomposition via source filtering",
Journal of Computational Physics,
Volume 188, Issue 2,
2003,
[DOI: 10.1016/S0021-9991(03)00168-2](https://doi.org/10.1016/S0021-9991(03)00168-2)
"""
struct AcousticPerturbationEquations2D{RealT<:Real} <: AbstractAcousticPerturbationEquations{2, 3}
  v_mean::SVector{2, RealT}
  c_mean::RealT
  rho_mean::RealT
end

function AcousticPerturbationEquations2D(v_mean::NTuple{2,<:Real}, c_mean::Real, rho_mean::Real)
  return AcousticPerturbationEquations2D(SVector(v_mean), c_mean, rho_mean)
end

function AcousticPerturbationEquations2D(; v_mean::NTuple{2,<:Real}, c_mean::Real, rho_mean::Real)
  return AcousticPerturbationEquations2D(SVector(v_mean), c_mean, rho_mean)
end

get_name(::AcousticPerturbationEquations2D) = "AcousticPerturbationEquations2D"
varnames(::typeof(cons2cons), ::AcousticPerturbationEquations2D) = ("v1_prime", "v2_prime", "p_prime")
varnames(::typeof(cons2prim), ::AcousticPerturbationEquations2D) = ("v1_prime", "v2_prime", "p_prime")

"""
    initial_condition_constant(x, t, equations::AcousticPerturbationEquations2D)

A constant initial condition.
"""
function initial_condition_constant(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = 0.0

  return SVector(v1_prime, v2_prime, p_prime)
end

"""
initial_condition_gauss(x, t, equation::AcousticPerturbationEquations2D)

A Gaussian pulse.
"""
function initial_condition_gauss(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = exp(-4*(x[1]^2 + x[2]^2))

  return SVector(v1_prime, v2_prime, p_prime)
end

"""
    initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref).
"""
function initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2D)
  c = 2.0
  A = 0.2
  L = 2.0
  f = 2.0 / L
  a = 1.0
  omega = 2 * pi * f
  init = c + A * sin(omega * (x[1] + x[2] - a*t))

  v1_prime = init
  v2_prime = init
  p = init^2

  return SVector(v1_prime, v2_prime, p)
end

"""
  source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref).
"""
function source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D)
  @unpack v_mean, rho_mean, c_mean = equations

  c = 2.0
  A = 0.2
  L = 2.0
  f = 2.0 / L
  a = 1.0
  omega = 2 * pi * f

  si, co = sincos(omega * (x[1] + x[2] - a * t))
  tmp = v_mean[1] + v_mean[2] - a

  du1 = du2 = A * omega * co * (2 * c + tmp + 2/rho_mean * A * si)
  du3 = A * omega * co * (2 * c_mean^2 * rho_mean + 2 * c * tmp + 2 * A * tmp * si)

  return SVector(du1, du2, du3)
end

# Calculate 1D flux for a single point
@inline function flux(u, orientation, equations::AcousticPerturbationEquations2D)
  v1_prime, v2_prime, p_prime = u
  @unpack v_mean, rho_mean, c_mean = equations

  if orientation == 1
    f1 = v_mean[1] * v1_prime + v_mean[2] * v2_prime + p_prime / rho_mean
    f2 = zero(eltype(u))
    f3 = c_mean^2 * rho_mean * v1_prime + v_mean[orientation] * p_prime
  else
    f1 = zero(eltype(u))
    f2 = v_mean[1] * v1_prime + v_mean[2] * v2_prime + p_prime / rho_mean
    f3 = c_mean^2 * rho_mean * v2_prime + v_mean[orientation] * p_prime
  end

  return SVector(f1, f2, f3)
end

function flux_lax_friedrichs(u_ll, u_rr, orientation, equations::AcousticPerturbationEquations2D)
  f_ll = flux(u_ll, orientation, equations)
  f_rr = flux(u_rr, orientation, equations)

  v_ll = u_ll[orientation] + equations.v_mean[orientation]
  v_rr = u_rr[orientation] + equations.v_mean[orientation]
  c0 = equations.c_mean
  speed = max(abs(v_ll), abs(v_rr)) + c0

  return 0.5 * ( (f_ll + f_rr) - speed * (u_rr - u_ll) )
end

@inline have_constant_speed(::AcousticPerturbationEquations2D) = Val(true)

@inline function max_abs_speeds(equations::AcousticPerturbationEquations2D)
  return abs.(equations.v_mean) .+ equations.c_mean
end

# Convert conservative variables to primitive
@inline cons2prim(u, equations::AcousticPerturbationEquations2D) = u

# Convert conservative variables to entropy variables
@inline cons2entropy(u, equations::AcousticPerturbationEquations2D) = u


text:
# -----------------------------------------------------------------------------------------
# This file includes the single argument definition of `scitype` method and corresponding
# convenience functions. It also includes definitions for `scitype/Scitype` of different
# objects with respect to the `DefaultConvention``.
# -----------------------------------------------------------------------------------------
"""
The scientific type (interpretation) of `X`, distinct from its
machine type.

### Examples
```
julia> scitype(3.14)
Continuous

julia> scitype([1, 2, missing])
AbstractVector{Union{Missing, Count}}

julia> scitype((5, "beige"))
Tuple{Count, Textual}

julia> using CategoricalArrays

julia> X = (gender = categorical(['M', 'M', 'F', 'M', 'F']),
     ndevices = [1, 3, 2, 3, 2])

julia> scitype(X)
Table{Union{AbstractVector{Count}, AbstractVector{Multiclass{2}}}}
```
"""
scitype(X) = ST.scitype(X, CONV)

function ST.scitype(@nospecialize(X), C::DefaultConvention)
    return _scitype(X, C, vtrait(X))
end

function _scitype(X, C, ::Val{:other})
    return ST.fallback_scitype(X, C)
end

# -----------------------------------------------------------------------------------------
# Basic (scalar) scitypes
ST.scitype(::Integer, ::DefaultConvention) = Count
ST.scitype(::AbstractFloat, ::DefaultConvention) = Continuous
ST.scitype(::AbstractString, ::DefaultConvention) = Textual
ST.scitype(::TimeType, ::DefaultConvention) = ScientificTimeType
ST.scitype(::Time, ::DefaultConvention) = ScientificTime
ST.scitype(::Date, ::DefaultConvention) = ScientificDate
ST.scitype(::DateTime, ::DefaultConvention) = ScientificDateTime

# -----------------------------------------------------------------------------------------
# Image scitypes
ST.scitype(img::Arr{<:Gray,2}, ::DefaultConvention) = GrayImage{size(img)...}
ST.scitype(img::Arr{<:AbstractRGB,2}, ::DefaultConvention) =
ColorImage{size(img)...}

# -----------------------------------------------------------------------------------------
# CategoricalArray scitype
function ST.scitype(c::Cat, ::DefaultConvention)
    nc = length(levels(c.pool))
    return ifelse(c.pool.ordered, OrderedFactor{nc}, Multiclass{nc})
end

const CatArrOrSub{T, N} =
    Union{CategoricalArray{T, N}, SubArray{T, N, <:CategoricalArray}}

function ST.scitype(A::CatArrOrSub{T, N}, ::DefaultConvention) where {T, N}
    nlevels = length(levels(A))
    S = ifelse(isordered(A), OrderedFactor{nlevels}, Multiclass{nlevels})
    T >: Missing && (S = Union{S, Missing})
    return AbstractArray{S, N}
end

# -----------------------------------------------------------------------------------------
## Table scitypes

function _scitype(X, C, ::Val{:table})
    if Tables.columnaccess(X)
        return cols_scitype(X)
    else
        return rows_scitype(X)
    end
end

function cols_scitype(X)
    cols = Tables.columns(X)
    sch = Tables.schema(cols)
    return _cols_scitype(cols, sch)
end

function _cols_scitype(cols, ::Union{Nothing, Tables.Schema{nothing, nothing}})
    col_names = Tables.columnnames(cols)
    scitypes = map(col_names) do name
        scitype(Tables.getcolumn(cols, name))
    end
    return Table{Union{scitypes...}}
end

function _cols_scitype(cols, sch::Tables.Schema{names, types}) where {names, types}
    N = length(names)
    if N <= COLS_SPECIALIZATION_THRESHOLD
        return __cols_scitype(cols, sch)
    else
        scitypes = if types === nothing
            Type[scitype(Tables.getcolumn(cols, names[i])) for i in Base.OneTo(N)]
        else
            Type[
                scitype(
                    Tables.getcolumn(
                    cols, fieldtype(types, i), i, names[i])
                ) for i in Base.OneTo(N)
            ]

        end
        return Table{Union{scitypes...}}
    end
end

@inline function __cols_scitype(
    cols,
    sch::Tables.Schema{names, types}
) where {names, types}
    N = length(names)
    if @generated
        stypes = if types === nothing
            (
                :(scitype(Tables.getcolumn(cols, $(Meta.QuoteNode(names[i])))))
                for i in Base.OneTo(N)
            )
        else
            (
                quote
                    scitype(
                        Tables.getcolumn(
                            cols, $(fieldtype(types, i)), $i,  $(Meta.QuoteNode(names[i]))
                        );
                    )
                end
                for i in Base.OneTo(N)
            )
        end
        return :(Table{Union{$(stypes...)}})
    else
        stypes = if types === nothing
            (scitype(Tables.getcolumn(cols, names[i])) for i in Base.OneTo(N))
        else
            (
                scitype(Tables.getcolumn(cols, fieldtype(types, i), i, names[i]))
                for i in Base.OneTo(N)
            )
        end
        return Table{Union{stypes...}}
    end
end

function rows_scitype(X)
    rows = Tables.rows(X)
    sch = Tables.schema(rows)
    return _rows_scitype(rows, sch)
end

function _rows_scitype(rows, sch::Tables.Schema{names, types}) where {names, types}
    rows_ = RowsWrapper(rows)
    cols = if length(names) <= ROWS_SPECIALIZATION_THRESHOLD
        Tables.columns(rows_)
    else
        Tables.dictcolumntable(rows_)
    end
    return _cols_scitype(cols, sch)
end

function _rows_scitype(rows, sch::Union{Nothing, Tables.Schema{nothing, nothing}})
    cols = Tables.dictcolumntable(RowsWrapper(rows))
    return _cols_scitype(cols, sch)
end

# -----------------------------------------------------------------------------------------
# Distributions scitype

const Dist = Distributions

scalar_scitype(::Type) = Unknown
scalar_scitype(::Type{<:Dist.Discrete}) = Count
scalar_scitype(::Type{<:Dist.Continuous}) = Continuous

function space_scitype(
    variate_form::Type{<:Dist.ArrayLikeVariate{0}},
    value_support
)
    return scalar_scitype(value_support)
end

function space_scitype(
    variate_form::Type{<:Dist.ArrayLikeVariate{N}},
    value_support
) where N
    return AbstractArray{scalar_scitype(value_support), N}
end

function ST.scitype(::Distributions.Sampleable{F, S}, ::DefaultConvention) where {F, S}
    return Sampleable{space_scitype(F, S)}
end

function ST.scitype(::Distributions.Distribution{F,S}, ::DefaultConvention) where {F, S}
    return Density{space_scitype(F, S)}
end

# -----------------------------------------------------------------------------------------
# `Scitype` definitions for fast array broadcasting

ST.Scitype(::Type{<:Integer}, ::DefaultConvention) = Count
ST.Scitype(::Type{<:AbstractFloat}, ::DefaultConvention) = Continuous
ST.Scitype(::Type{<:AbstractString}, ::DefaultConvention) = Textual
ST.Scitype(::Type{<:TimeType}, ::DefaultConvention) = ScientificTimeType
ST.Scitype(::Type{<:Date}, ::DefaultConvention) = ScientificDate
ST.Scitype(::Type{<:Time}, ::DefaultConvention) = ScientificTime
ST.Scitype(::Type{<:DateTime}, ::DefaultConvention) = ScientificDateTime
# ST.Scitype(::Type{<:PersistenceDiagram}, ::DefaultConvention) = PersistenceDiagram

# -----------------------------------------------------------------------------------------
## convenience methods
"""
    elscitype(A)
Return the element scientific type of an abstract array `A`. By definition, if
`scitype(A) = AbstractArray{S,N}`, then `elscitype(A) = S`.
"""
elscitype(X) = elscitype(collect(X))
elscitype(X::Arr) = eltype(scitype(X))

"""
    scitype_union(A)
Return the type union, over all elements `x` generated by the iterable `A`,
of `scitype(x)`. See also [`scitype`](@ref).
"""
scitype_union(X) = ST.scitype_union(X, DefaultConvention())


text:
"""
    download_fhir_json_schema(url::AbstractString)

Download and unzip the file at `url`.

## Example
```julia
julia> url = "https://www.hl7.org/fhir/R4/definitions.json.zip"
julia> Generate.download_fhir_json_schema(url)
```
"""
@inline function download_fhir_json_schema(url::AbstractString)
    io = IOBuffer()
    Downloads.download(url, io)
    r = ZipFile.Reader(io)
    return _process_zip_reader(r)
end

@inline function _process_zip_reader(r::ZipFile.Reader)
    result = Dict{String, Union{String, Dict}}()
    files = r.files
    for i in 1:length(files)
        file = files[i]
        filename = file.name
        result[filename] = _process_zip_file(filename, file)
    end
    return result
end

@inline function _process_zip_file(filename::AbstractString, file::ZipFile.ReadableFile)
    filecontents = read(file, String)::String
    if endswith(lowercase(strip(filename)), ".zip")
        io = IOBuffer(filecontents)
        r = ZipFile.Reader(io)
        return _process_zip_reader(r)
    end
    return filecontents
end


text:
@safetestset "correct max phase grad" begin
    using JultraDark:max_phase_grad

    n = 16

    @testset "purely real" begin
        a = reshape(Vector{Complex}(1:n^3), (n, n, n))
        @test max_phase_grad(a) == 0
    end

    @testset "purely imaginary" begin
        a = reshape(Vector{Complex}(1:n^3), (n, n, n)) * im
        @test max_phase_grad(a) == 0
    end

    @testset "real with 1 row negated" begin
        a = reshape(Vector{Complex}(1:n^3), (n, n, n))
        a[1, :, :] *= -1
        @test max_phase_grad(a) ≈ π
    end

    @testset "real with 1 row negated" begin
        a = reshape(Vector{Complex}(1:n^3), (n, n, n))
        a[:, 3, :] *= -1
        @test max_phase_grad(a) ≈ π
    end

    @testset "real with 1 row negated" begin
        a = reshape(Vector{Complex}(1:n^3), (n, n, n))
        a[:, :, 5] *= -1
        @test max_phase_grad(a) ≈ π
    end

    @testset "real with 1 row imaginary" begin
        a = reshape(Vector{Complex}(1:n^3), (n, n, n))
        a[1, :, :] *= im
        @test max_phase_grad(a) ≈ π / 2
    end

end

@safetestset "terminate if grad too large" begin
    using JultraDark
    using JultraDark: evolve_to!, max_phase_grad

    n = 16

    # Set up grid with discontinuous phase
    grids = Grids(1.0, n)
    grids.ψx .= 1
    grids.ψx[1, 1, 1] = -1

    grids_orig = deepcopy(grids)

    options = JultraDark.Config.SimulationConfig(10, t->1)
    output_config = OutputConfig(mktempdir(), [1, 2])

    @test_throws "Phase gradient is too large to start" simulate(grids, options, output_config)
    @test grids.ψx == grids_orig.ψx

    @test_throws "Phase gradient is too large to continue" evolve_to!(0, 1, grids, output_config, options)
    @test grids.ψx == grids_orig.ψx

end


text:
using MechanicalSketch
import MechanicalSketch: Point, rect, arrow, line
import MechanicalSketch: EM, TXTOFFV, O, PT, Point, WI, HE
import MechanicalSketch: RGB, HSL, color_with_lumin, background, sethue, text, fontsize
import MechanicalSketch: settext, setfont, setopacity, FS, color_with_lumin, color_with_lumin2
import MechanicalSketch: empty_figure, PALETTE, finish, paint, ∙
import MechanicalSketch: transform, @layer, do_action, getmatrix, grestore, gsave, scale
import MechanicalSketch: box_fill_outline, label, row_height, pixelwidth, @import_expand
import MechanicalSketch: scale_pt_to_unit, lumin, label_boxed
if !@isdefined °
    @import_expand(~m, s, °, rad)
end
#let
    typcol = PALETTE[9];
    bckcol = color_with_lumin(typcol, 0.5 * lumin(typcol))
    empty_figure(filename = joinpath(@__DIR__, "test_3.png");
        backgroundcolor = bckcol, hue = "white")

    # Straight arrow width
    pt = O + ( -6EM, 2EM)
    arrow(pt, pt + (0, -5EM) )
    arrow(pt, pt + (10EM, -2EM), arrowheadlength=EM, arrowheadangle = pi/12, linewidth = PT)
    arrow(pt, pt + (8EM, 0.5EM), arrowheadlength=EM, arrowheadangle = pi/12, linewidth = PT)

    # Palette
    ncol = length(PALETTE)
    dw =  WI / ncol
    for i = 1:ncol
        posy = -HE/2
        curcol = PALETTE[i]
        sethue(curcol)
        xs = -WI/2 + (i-1) * dw
        ps = Point(xs, posy)
        rect( ps, dw, EM, :fill)
        sethue("white")
        fontsize(FS)
        text("$i L= $(round(HSL(curcol).l, digits = 2))",
                    ps + (0, EM - TXTOFFV))
        fontsize(FS/2)
        setfont("Calibri", FS / 2)
        posy = HE / 2 - 22 * 0.5 * EM
        text("lumin", Point(xs, posy))
        text("lumin2", Point(xs + dw /2, posy))
        for (j, lum) in enumerate(2:-0.1:0)
            curcol = PALETTE[i]
            # Left rectancle, simple luminance
            sethue(color_with_lumin(curcol, lum * 100))
            posy += 0.5 * EM
            ps = Point(xs, posy)
            rect( ps , dw / 2 , EM / 2 , :fill)
            str = string(round(lum, digits = 2))
            sethue("black")
            text(str,  ps + (0, -TXTOFFV + 0.5 *EM) )
            sethue("white")
            text(str,  ps + (dw / 4, -TXTOFFV + 0.5 *EM) )
            # Right rectangle, lumin2
            sethue(color_with_lumin2(curcol, lum * 2))
            rect(ps + (dw/2, 0), dw / 2, EM /2, :fill)
            sethue("black")
            settext("" * str * "",
                    ps + (dw/2, -TXTOFFV + 0.6 *EM) , markup=true)
            sethue("white")
            settext("" * str * "",
                        ps + (dw/2 + dw/ 4, -TXTOFFV + 0.6 *EM) , markup=true)
        end
    end
    # Arc arrows, left screen
    sethue(PALETTE[3])
    y = -HE / 4 + 0EM
    x = - WI / 3
    cpl = O + (x, y)
    text("Unitless angles: Always clockwise", cpl + (0, EM), halign = :center)
    startangle = 0.0
    iter  = enumerate(range(-π , π, step = π / 4))
    for (i, endangle) in iter
        # center, radius, start (cw), end(cw)
        radius = (3 + i / 4 )EM
        arrow(cpl, radius, startangle, endangle;
               arrowheadlength = EM, arrowheadangle = pi / 12, linewidth = 2PT)
    end
    for (i, endangle) in iter
        radius = (3 + i / 4 )EM
        str = "$(round(startangle / π; digits = 2)) to $(round(endangle / π; digits = 2))∙π"
        endpoint = cpl + radius .* (cos(endangle), sin(endangle))
        label_boxed(endpoint, str)
    end

    # Arc arrows with units, mid-screen
    cpm = O + (0, y)
    text("Uniful angles, specified clockwise", cpm + (0, EM), halign = :center)
    startangle = 0.0∙rad
    iter  = enumerate(range(-π, π, step = π / 4)∙rad)
    rdeg(a) = round(°, a; digits = 0) |> typeof(1°)
    for (i, endangle) in iter
        radius = (3 + i / 4 )EM * scale_pt_to_unit(m)
        arrow(cpm, radius, startangle, endangle;
              arrowheadlength = EM, arrowheadangle = pi / 12, linewidth = 2PT, clockwise = true)
    end
    for (i, endangle) in iter
        radius = (3 + i / 4 )EM
        str = "$(rdeg(startangle)) to $(rdeg(endangle))"
        endpoint = cpm + radius .* (cos(endangle), -sin(endangle))
        label_boxed(endpoint, str)
    end

    # Arc arrows with units, right-screen
    cpr = O + (-x, y)
    text("Uniful angles, default counterclockwise", cpr + (0, EM), halign = :center)
    startangle = 0.0∙rad
    iter  = enumerate(range(-π, π, step = π / 4)∙rad)
    rdeg(a) = round(°, a; digits = 0) |> typeof(1°)
    for (i, endangle) in iter
        radius = (3 + i / 4 )EM * scale_pt_to_unit(m)
        arrow(cpr, radius, startangle, endangle;
              arrowheadlength = EM, arrowheadangle = pi / 12, linewidth = 2PT)
    end
    for (i, endangle) in iter
        radius = (3 + i / 4 )EM
        str = "$(rdeg(startangle)) to $(rdeg(endangle))"
        endpoint = cpr + radius .* (cos(endangle), -sin(endangle))
        label_boxed(endpoint, str)
    end

    finish()

#end


text:
# Copied from NEP-PACK

using IterativeSolvers
using LinearAlgebra
using Random

"""
   (Z,a,H,hist)=tiar(Float64,prob,btilde,m,Afact)

#from: https://github.com/nep-pack/NonlinearEigenproblems.jl/blob/master/src/method_tiar.jl

tiar(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS,][proj_solve=false,][inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

Run the tensor infinite Arnoldi method on the nonlinear eigenvalue problem stored in `nep`. This is equivalent to `iar`, but handles orthogonalization with
a tensor representation.

The target `σ` is the center around which eiganvalues are computed.
The value `γ` corresponds
to scaling and specifying a shift and scaling is effectively the same as the
transformation `λ=γs+σ` where `s`
is now the eigenvalue parameter. If you want eigenvalues in a disk centered,
select `σ` as the center
of the disk and `γ` as the radius.
The vector `v` is the starting vector for constructing the Krylov space. The orthogonalization
method, used in contructing the orthogonal basis of the Krylov space, is specified by
`orthmethod`, see the package `IterativeSolvers.jl`.
The iteration is continued until `neigs` Ritz pairs have converged.
This function throws a `NoConvergenceException` if the wanted eigenpairs are not computed after `maxit` iterations.
However, if `neigs` is set to `Inf` the iteration is continued until `maxit` iterations without an error being thrown.
The parameter `proj_solve` determines if the Ritz paris are extracted using the Hessenberg matrix (false),
or as the solution to a projected problem (true). If true, the method is descided by `inner_solver_method`, and the
logging of the inner solvers are descided by `inner_logger`, which works in the same way as `logger`.

See [`augnewton`](@ref) for other parameters.

### References
Algorithm 2 in Jarlebring, Mele, Runborg, The Waveguide Eigenvalue Problem and the Tensor Infinite Arnoldi Method

"""
function tiar(
    ::Type{T},
    nep,
    b,maxit,precomputed_factorization=nothing;
    orthmethod::Type{T_orth}=DGKS,
    tol=eps(real(T))*10000,
    σ=zero(T),
    γ=one(T),
    check_error_every=1)  where{T,T_orth<:IterativeSolvers.OrthogonalizationMethod}

    start_time=time();
    hist=Dict(:time_count => NaN*zeros(N))

    if (isnothing(precomputed_factorization))
        Msolve=factorize(compute_Mder(nep,σ));
    else
        Msolve=precomputed_factorization;
    end

    # Ensure types σ and v are of type T
    σ=T(σ)
    v=Array{T,1}(b)

    # initialization
    n = size(b,1); m = maxit;

    # initialize variables
    a  = zeros(T,m+1,m+1,m+1);
    Z  = zeros(T,n,m+1);
    t  = zeros(T,m+1);
    tt = zeros(T,m+1);
    g  = zeros(T,m+1,m+1);
    f  = zeros(T,m+1,m+1);
    ff = zeros(T,m+1,m+1);
    H  = zeros(T,m+1,m);
    h  = zeros(T,m+1);
    hh = zeros(T,m+1);
    y  = zeros(T,n,m+1);
    α=Array{T,1}(γ.^(0:m)); α[1]=zero(T);

    # local M0inv::LinSolver=create_linsolver(linsolvercreator,nep,σ)

    err = NaN*ones(m+1,m+1);
    λ=zeros(T,m+1); Q=zeros(T,n,m+1);
    Z[:,1]=v; Z[:,1]=Z[:,1]/norm(Z[:,1]);
    a[1,1,1]=one(T);

    print("TIAR Iteration:");

    k=1; conv_eig=0;
    for k=1:m
        print("$k ")

        # computation of y[:,2], ..., y[:,k+1]
        y[:,2:k+1]=Z[:,1:k]*transpose(a[1:k,k,1:k])
        broadcast!(/,view(y,:,2:k+1),view(y,:,2:k+1),(1:k)')

        # computation of y[:,1]
        YY=y[:,1:k+1]*Diagonal(α[1:k+1])
        y[:,1] = compute_Mlincomb(nep,σ,YY);
        y[:,1] = -(Msolve\y[:,1]);

        # Gram–Schmidt orthogonalization in Z
        Z[:,k+1]=y[:,1];
        t[k+1] = orthogonalize_and_normalize!(view(Z,:,1:k), view(Z,:,k+1), view(t,1:k), orthmethod)

        # compute the matrix G
        for l=1:k+1
            for i=2:k+1
                g[i,l]=a[i-1,k,l]/(i-1);
            end
         g[1,l]=t[l];
        end

        # compute h (orthogonalization with tensors factorization)
        h = zero(h)
        Ag = zero(h[1:k])
        for l=1:k
            mul!(Ag,a[1:k,1:k,l]',g[1:k,l])
            h[1:k] .+= Ag;
        end

        # compute the matrix F
        f=g;
        Ah = zero(f[1:k+1,1])
        for l=1:k
            mul!(Ah,a[1:k+1,1:k,l],h[1:k])
            f[1:k+1,l] .-= Ah;
        end

        # re-orthogonalization
        # compute hh (re-orthogonalization with tensors factorization)
        hh = zero(hh)
        Af = zero(hh[1:k])
        for l=1:k
            mul!(Af,a[1:k,1:k,l]',f[1:k,l])
            hh[1:k] .+= Af;
        end

        # compute the matrix FF
        ff=f;
        Ah=zero(ff[1:k+1,1])
        for l=1:k
            mul!(Ah,a[1:k+1,1:k,l],hh[1:k])
            ff[1:k+1,l] .-= Ah;
        end

        # update the orthogonalization coefficients
        h=h+hh; f=ff;
        β=norm(view(f,1:k+1,1:k+1)); # equivalent to Frobenius norm

        # extend the matrix H
        H[1:k,k]=h[1:k]; H[k+1,k]=β;

        # extend the tensor
        for i=1:k+1
            for l=1:k+1
                a[i,k+1,l]=f[i,l]/β;
            end
        end
        hist[:time_count][k]=time()-start_time;

    end
    return (Z,a,H,hist)
end

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
using Base.Test
using CUDArt
using Knet

for N in 1:5
    for A in (CudaArray, Array)
        for T in (Float64, Float32, Int64, Int32)
            S = tuple(rand(1:20,N)...)
            @show (A,T,S)
            a = convert(A{T}, rand(T,S))
            b = convert(A{T}, rand(T,csize(a,rand(1:20))))
            d = KUdense(copy(a))
            @test atype(d)==A
            @test d.arr == a
            @test pointer(d) == pointer(d.arr) == pointer(d.ptr)
            for fname in (:eltype, :length, :ndims, :size, :strides, :isempty)
                @test (@eval $fname($d)==$fname($a))
            end
            for fname in (:size, :stride)
                for n in 1:N
                    @test (@eval $fname($d,$n)==$fname($a,$n))
                end
            end

            if isa(a,Array)
                r = rand(1:length(d))
                a[r] = d[r] = 42
                @test a[r] == d[r] == 42
                rr = map(x->rand(1:x), size(d))
                a[rr...] = d[rr...] = 53
                @test a[rr...] == d[rr...] == 53
            end

            dd = copy(d)
            copy!(dd, b)
            @test dd.arr == b
            copy!(dd, d)
            @test dd.arr == d.arr
            @test pointer(dd) == pointer(dd.arr) == pointer(dd.ptr)
            @test pointer(dd) != pointer(d)

            q = csize(d, rand(1:ccount(d)))
            qq = map(x->(1:x), q)
            resize!(d, q)
            @test to_host(d.arr) == to_host(a)[qq...]

            copy!(d, a)
            setseed(42); rand!(a)
            setseed(42); rand!(d)
            @test d.arr == a
            if isa(eltype(a), FloatingPoint)
                setseed(42); randn!(a)
                setseed(42); randn!(d)
                @test d.arr == a
            end

            cd = cpucopy(d)
            @test pointer(cd) == pointer(cd.arr) == pointer(cd.ptr)
            @test pointer(cd) != pointer(d)
            gd = gpucopy(d)
            @test pointer(gd) == pointer(gd.arr) == pointer(gd.ptr)
            @test pointer(gd) != pointer(d)
            ca = cpucopy(a)
            @test cd.arr == ca
            ga = gpucopy(a)
            @test gd.arr == ga

            # cslice!
            a1 = ndims(a)>1 ? a : reshape(a,length(a),1)
            r1 = rand(1:ccount(a1))
            r2 = rand(1:ccount(a1))
            r1 <= r2 || ((r1,r2)=(r2,r1))
            cslice!(d, a1, r1:r2)
            rr = UnitRange[map(x->(1:x),size(a1)[1:end-1])..., (r1:r2)]
            @test to_host(d.arr) == to_host(a1)[rr...]

            # ccopy!
            a2 = copy(a1)
            d2 = copy(d)
            rand!(d2)
            ccopy!(a2,r1,d2)
            b2 = to_host(d2.arr)
            rr1 = UnitRange[map(x->(1:x),size(a1)[1:end-1])..., (1:r1-1)]
            !isempty(to_host(a1)[rr1...]) && (b2 = cat(ndims(b2), to_host(a1)[rr1...], b2))
            rr2 = UnitRange[map(x->(1:x),size(a1)[1:end-1])..., (r2+1:ccount(a1))]
            !isempty(to_host(a1)[rr2...]) && (b2 = cat(ndims(b2), b2, to_host(a1)[rr2...]))
            @test to_host(a2) == b2

            # ccat!
            da = KUdense(copy(a1))
            db = KUdense(copy(b))
            ccat!(da, db)
            @test to_host(da.arr) == cat(ndims(a1), to_host(a1), to_host(b))

            # uniq!
            for i=1:5; ccat!(da, db); end
            ua = unique(to_host(da.arr), ndims(da))
            uniq!(da)
            @test to_host(da.arr) == ua
        end
    end
end


text:
function unsafe_sizes(buff::Ptr{UInt8})
    r = Lib.blosc_cbuffer_sizes(buff)
    return (r.nbytes, r.cbytes)
end

"""
    sizes(buff::Vector{UInt8}, offset = 1)

Get information about a compressed buffer at offset `offset` (1-indexed), as Tuple with values:

- the number of uncompressed bytes

- the number of compressed bytes
"""
function sizes(buff::Vector{UInt8}, offset = 1)
    @boundscheck checkbounds(buff, offset)
    r = GC.@preserve buff Lib.blosc_cbuffer_sizes(pointer(buff, offset))
    return (r.nbytes, r.cbytes)
end

"""
    uncompressed_sizof(buff::Vector{UInt8}, offset = 1)

Uncompressed size of data in buffer `buff` in bytes`
"""
Base.@propagate_inbounds uncompressed_sizeof(buff::Vector{UInt8}, offset = 1) = sizes(buff)[1]

"""
    uncompressed_length(::Type{T}, buff::Vector{UInt8}, offset = 1)

Uncompressed length of `Vector{T}` then decompressed from buffer `buff`
"""
Base.@propagate_inbounds uncompressed_length(::Type{T},buff::Vector{UInt8}, offset = 1) where {T} = ceil(Int, sizes(buff)[1] / sizeof(T))

"""
    compressed_sizeof(buff::Vector{UInt8}, offset = 1)

Compressed size of data in buffer `buff` in bytes
"""
Base.@propagate_inbounds compressed_sizeof(buff::Vector{UInt8}, offset = 1) = sizes(buff)[2]


text:
# Uniform rand with range e
randrealuni() = Real(MathConstants.e * rand())  # may Euler's number give us good luck
randrealuni(args...) = map(Real, 2 * rand(args...))

const Transformable = Union{TransformDistribution, SimplexDistribution, PDMatDistribution}

#################################
# Single-sample initialisations #
#################################

init(dist::Transformable) = inittrans(dist)
init(dist::Distribution) = rand(dist)

inittrans(dist::UnivariateDistribution) = invlink(dist, randrealuni())
inittrans(dist::MultivariateDistribution) = invlink(dist, randrealuni(size(dist)[1]))
inittrans(dist::MatrixDistribution) = invlink(dist, randrealuni(size(dist)...))

################################
# Multi-sample initialisations #
################################

init(dist::Transformable, n::Int) = inittrans(dist, n)
init(dist::Distribution, n::Int) = rand(dist, n)

inittrans(dist::UnivariateDistribution, n::Int) = invlink(dist, randrealuni(n))
function inittrans(dist::MultivariateDistribution, n::Int)
  return invlink(dist, randrealuni(size(dist)[1], n))
end
function inittrans(dist::MatrixDistribution, n::Int)
  return invlink(dist, [randrealuni(size(dist)...) for _ in 1:n])
end


text:
#db_access.jl: Write/read data to/from HDF5 "database".
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
"Exception: Invalid element of curricculum"
struct InvalidElement <: Exception
	msg
end

"Exception: Missing path in database"
struct MissingPath <: Exception
	msg
end

#==Constants
===============================================================================#
const DATA_UNKNOWNSRCDOC = "Source inconnue???"
const DATA_NEWDOMAIN = "Nouveau domaine"
const DATA_NEWATTENTE = "Nouvelle attente"
const DATA_NEWATTENTE_SHORT = "NA"
const DATA_NEWCONTENT = "Nouveau contenu"

const DATA_NODOMAIN = (0, "", "Aucun domaine")
const DATA_NOATTENTE = (0, "", "Aucune attente", "Aucune attente")
const DATA_NOCONTENT = (0, "", "Aucun contenu")

const MSG_NODATA = "Aucune donnée trouvée."
const MSG_NOSRCDOCINFO = "Source non-spécifiée."
const MSG_NOSUBJECT = "Aucun sujet sélectionné!"
const MSG_CANNOTADDENTRY = "Incapable d'ajouter un item!"

const ERROR_NOSUBJECT = InvalidElement(MSG_NOSUBJECT)
const ERROR_CANNOTADDENTRY = ErrorException(MSG_CANNOTADDENTRY)

#Attribute identifiers:
const AID_NELEM = "NELEM"
const AID_SRCDOC = "SRCDOC"
const AID_DESCR = "DESCR"
const AID_DESCR_SHORT = "DESCR_SHORT"

#==Path utils
===============================================================================#
getpath_subject(sel::ExploreSelection) = "/$(sel.subject)"
function getpath_grade(sel::ExploreSelection)
	pfx = getpath_subject(sel)
	gstr = getid_grade(sel.grade_idx)
	return "$pfx/$gstr"
end
function getpath_domains(sel::ExploreSelection)
	pfx = getpath_grade(sel)
	return "$pfx/domaines"
end
function getpath_seldomain(sel::ExploreSelection)
	pfx = getpath_domains(sel)
	return "$pfx/$(sel.domain_idx)"
end
function getpath_attentes(sel::ExploreSelection)
	pfx = getpath_seldomain(sel)
	return "$pfx/attentes"
end
function getpath_selattente(sel::ExploreSelection)
	pfx = getpath_attentes(sel)
	return "$pfx/$(sel.attente_idx)"
end
function getpath_content(sel::ExploreSelection)
	pfx = getpath_selattente(sel)
	return "$pfx/contenus"
end

#==Constructors
===============================================================================#
function InvalidElement(msg::String)
	return InvalidElement(msg)
end

function MissingPath(msg::String)
	return MissingPath(msg)
end

#==Base accessors for database
===============================================================================#
#Open group, or return nothing !exist:
function gopen_nothing(o, path::String)
	if !HDF5.exists(o, path);
		return nothing
	else
		return HDF5.g_open(o, path)
	end
end
#Returns subgroup, or nothing
function getsubgrp_nothing(grp::HDF5.HDF5Group, name::String)
	try #Don't know how to tell if subgroup exists (only if SOMETHING exists)
		return HDF5.g_open(grp, name)
	catch
		return nothing
	end
end

#Open group, or create if !exist:
function gopen_create(o, path::String)
	if !HDF5.exists(o, path);
		return HDF5.g_create(o, path)
	else
		return HDF5.g_open(o, path)
	end
end

#==Read in datasets & attributes (deal with missing data)
===============================================================================#
read_attr(::Nothing, name::String, expType::Type, default) = default
function read_attr(grp::HDF5.HDF5Group, name::String, expType::Type, default)
	if !HDF5.exists(HDF5.attrs(grp), name); return default; end
	attr = HDF5.a_read(grp, name)
	if !isa(attr, expType); return default; end
	return attr
end
function reada_nelem(grp)
	nelem = read_attr(grp, AID_NELEM, Int64, 0)
	return Int(max(0, nelem)) #Limit >= 0
end
reada_srcdoc(grp) = read_attr(grp, AID_SRCDOC, String, MSG_NOSRCDOCINFO)
reada_descr(grp) = read_attr(grp, AID_DESCR, String, "")
reada_shortdescr(grp) = read_attr(grp, AID_DESCR_SHORT, String, "")

function readds_safe(grp::HDF5.HDF5Group, name::String, expType::Type, default)
	if !HDF5.exists(grp, name); return default; end
	v = HDF5.d_read(grp, name)
	if !isa(v, expType); return default; end
	return v
end

#==Write datasets & attributes:
===============================================================================#
function write_attr(grp::HDF5.HDF5Group, name::String, v)
	if HDF5.exists(HDF5.attrs(grp), name)
		HDF5.a_delete(grp, name)
	end
	HDF5.a_write(grp, name, v)
end
writea_nelem(grp, v::Vector) = write_attr(grp, AID_NELEM, length(v))
writea_srcdoc(grp, v::String) = write_attr(grp, AID_SRCDOC, v)
writea_descr(grp, v::String) = write_attr(grp, AID_DESCR, v)
writea_shortdescr(grp, v::String) = write_attr(grp, AID_DESCR_SHORT, v)

function writeds_safe(grp::HDF5.HDF5Group, name::String, v)
	if HDF5.exists(grp, name); HDF5.o_delete(grp, name); end
	HDF5.d_write(grp, name, v)
	return
end

#==Validation
===============================================================================#
function validate_subject(sel::ExploreSelection)
	if isempty_subject(sel.subject); throw(ERROR_NOSUBJECT); end
	return
end
function validate_grade(sel::ExploreSelection)
	validate_subject(sel)
	if sel.grade_idx < 1; throw(InvalidElement(getpath_grade(sel))); end
	return
end
function validate_domain(sel::ExploreSelection)
	validate_grade(sel)
	if sel.domain_idx < 1; throw(InvalidElement(getpath_domains(sel))); end
	return
end
function validate_attente(sel::ExploreSelection)
	validate_domain(sel)
	if sel.attente_idx < 1; throw(InvalidElement(getpath_attentes(sel))); end
	return
end
function validate_content(sel::ExploreSelection)
	validate_attente(sel)
	if sel.content_idx < 1; throw(InvalidElement(getpath_content(sel))); end
	return
end

#Must actually exist in db:
#-------------------------------------------------------------------------------
function validateexists_simple(db::HDF5.HDF5File, path::String)
	grp = gopen_nothing(db, path)
	if isnothing(grp); throw(MissingPath(path)); end
	return
end
function validateexists_listitem(db::HDF5.HDF5File, path::String, itemidx::Int)
	grp = gopen_nothing(db, path)
	if isnothing(grp); throw(MissingPath(path)); end
	if itemidx > reada_nelem(grp); throw(MissingPath(path)); end
	return
end
function validateexists_subject(db::HDF5.HDF5File, sel::ExploreSelection)
	validate_subject(sel)
	validateexists_simple(db, getpath_subject(sel))
end
function validateexists_grade(db::HDF5.HDF5File, sel::ExploreSelection)
	validateexists_subject(db, sel)
	validateexists_simple(db, getpath_grade(sel))
end
function validateexists_domain(db::HDF5.HDF5File, sel::ExploreSelection)
	validateexists_grade(db, sel)
	validateexists_listitem(db, getpath_domains(sel), sel.domain_idx)
end
function validateexists_attente(db::HDF5.HDF5File, sel::ExploreSelection)
	validateexists_domain(db, sel)
	validateexists_listitem(db, getpath_attentes(sel), sel.attente_idx)
end
function validateexists_content(db::HDF5.HDF5File, sel::ExploreSelection)
	validateexists_attente(db, sel)
	validateexists_listitem(db, getpath_content(sel), sel.content_idx)
end

#==Read in high-level fields. Returns nothing on error.
===============================================================================#
read_subjects(db::HDF5.HDF5File) = names(gopen_nothing(db, "/"))

function read_sourcedoc(db::HDF5.HDF5File, sel::ExploreSelection)
	ggrp = gopen_nothing(db, getpath_grade(sel))
	if isnothing(ggrp); return nothing; end
	return reada_srcdoc(ggrp)
end

function read_domain_list(db::HDF5.HDF5File, sel::ExploreSelection)
	result = []
	grp = gopen_nothing(db, getpath_domains(sel))
	if isnothing(grp)
		return nothing
		#unreachable test data:
		for i in 1:3
			id = getid_domain(i)
			descr = "Domaine $id"
			push!(result, (i, id, descr))
		end
		return result
	end

	#Read in data:
	nelem = reada_nelem(grp)
	for i in 1:nelem
		id = getid_domain(i)
		dgrp = getsubgrp_nothing(grp, "$i")
		descr = reada_descr(dgrp)
		push!(result, (i, id, descr))
	end
	return result
end

function read_attente_list(db::HDF5.HDF5File, sel::ExploreSelection)
	result = []
	grp = gopen_nothing(db, getpath_attentes(sel))
	if isnothing(grp)
		return nothing
		#unreachable test data:
		for i in 1:3
			id = getid_attente(sel.domain_idx, i)
			shortdescr = "Attente $id"
			descr = "Description $id"
			push!(result, (i, id, shortdescr, descr))
		end
	end

	#Read in data:
	nelem = reada_nelem(grp)
	for i in 1:nelem
		id = getid_attente(sel.domain_idx, i)
		agrp = getsubgrp_nothing(grp, "$i")
		shortdescr = reada_shortdescr(agrp)
		descr = reada_descr(agrp)
		push!(result, (i, id, shortdescr, descr))
	end
	return result
end

function read_content_list(db::HDF5.HDF5File, sel::ExploreSelection)
	result = []
	grp = gopen_nothing(db, getpath_content(sel))
	if isnothing(grp); return nothing; end

	#Read in data:
	nelem = reada_nelem(grp)
	for i in 1:nelem
		id = getid_content(sel.domain_idx, sel.attente_idx, i)
		descr = readds_safe(grp, "$i", String, "")
		push!(result, (i, id, descr))
	end
	return result
end

#==Write out high-level fields. Use createslot_* or validateexists_* first!
===============================================================================#
function write_sourcedoc(db::HDF5.HDF5File, sel::ExploreSelection, v::String)
	grp = gopen_create(db, getpath_grade(sel))
	writea_srcdoc(grp, v)
end
function write_domain_descr(db::HDF5.HDF5File, sel::ExploreSelection, descr::String)
	grp = gopen_create(db, getpath_seldomain(sel))
	writea_descr(grp, descr)
end
function write_attente_descr(db::HDF5.HDF5File, sel::ExploreSelection, shortdescr::String, descr::String)
	grp = gopen_create(db, getpath_selattente(sel))
	writea_shortdescr(grp, shortdescr)
	writea_descr(grp, descr)
end
function write_content_descr(db::HDF5.HDF5File, sel::ExploreSelection, descr::String)
	grp = gopen_create(db, getpath_content(sel))
	writeds_safe(grp, "$(sel.content_idx)", descr)
end

#==Create slot: Add new item to db (throw exception if cannot; return true if success)
===============================================================================#
function createslot_subject(db::HDF5.HDF5File, sel::ExploreSelection)
	validate_subject(sel)
	grp = gopen_create(db, getpath_subject(sel))
	return true
end
#Don't really need to explicitly create subject. Create source doc info instead:
function createslot_srcdoc(db::HDF5.HDF5File, sel::ExploreSelection)
	validate_grade(sel) #Otherwise, will generate garbage
	srcstr = read_sourcedoc(db, sel)
	if !isnothing(srcstr); return false; end #Don't overwrite, but no exception
	grp = gopen_create(db, getpath_subject(sel))
	write_sourcedoc(db, sel, DATA_UNKNOWNSRCDOC)
	return true #Created new slot
end
function createslot_domain(db::HDF5.HDF5File, sel::ExploreSelection)
	createslot_srcdoc(db, sel) #Just in case doesn't exist / also validates
	#Don't really need to validate anything else before creating...
	dlist = read_domain_list(db, sel)
	if isnothing(dlist); dlist = []; end
	push!(dlist, DATA_NEWDOMAIN) #Add something to compute length
	sel.domain_idx = length(dlist)
	grp = gopen_create(db, getpath_domains(sel))
	write_domain_descr(db, sel, DATA_NEWDOMAIN)
	writea_nelem(grp, dlist)
	return
end
function createslot_attente(db::HDF5.HDF5File, sel::ExploreSelection)
	validateexists_domain(db, sel) #Don't create if parent nodes don't exist
	dlist = read_attente_list(db, sel)
	if isnothing(dlist); dlist = []; end
	push!(dlist, DATA_NEWATTENTE) #Add something to compute length
	sel.attente_idx = length(dlist)
	grp = gopen_create(db, getpath_attentes(sel))
	write_attente_descr(db, sel, DATA_NEWATTENTE_SHORT, DATA_NEWATTENTE)
	writea_nelem(grp, dlist)
	return
end
function createslot_content(db::HDF5.HDF5File, sel::ExploreSelection)
	validateexists_attente(db, sel) #Don't create if parent nodes don't exist
	dlist = read_content_list(db, sel)
	if isnothing(dlist); dlist = []; end
	push!(dlist, DATA_NEWCONTENT) #Add something to compute length
	sel.content_idx = length(dlist)
	grp = gopen_create(db, getpath_content(sel))
	write_content_descr(db, sel, DATA_NEWCONTENT)
	writea_nelem(grp, dlist)
	return
end

#==Write functions for list of elements
===============================================================================#
#=NOTE
 - Write whole lists at once & update element count.
 - Not very robust.
 - Meant to manually add data using scripts.
=#
function write_domain_list(db::HDF5.HDF5File, sel::ExploreSelection, list::Vector{String})
	sel = deepcopy(sel)
	for (i, descr) in enumerate(list)
		sel.domain_idx = i
		write_domain_descr(db, sel, descr)
	end

	grp = gopen_create(db, getpath_domains(sel))
	writea_nelem(grp, list)
end
function write_attente_list(db::HDF5.HDF5File, sel::ExploreSelection, list::Vector)
	sel = deepcopy(sel)
	for (i, elem) in enumerate(list)
		sel.attente_idx = i
		(shortdescr, descr) = elem
		write_attente_descr(db, sel, shortdescr, descr)
	end

	grp = gopen_create(db, getpath_attentes(sel))
	writea_nelem(grp, list)
end
function write_content_list(db::HDF5.HDF5File, sel::ExploreSelection, list::Vector)
	sel = deepcopy(sel)
	for (i, elem) in enumerate(list)
		sel.content_idx = i
		write_content_descr(db, sel, elem)
	end

	grp = gopen_create(db, getpath_content(sel))
	writea_nelem(grp, list)
end

#Last line


text:
using HorizonSideRobots
include("Functions.jl")

struct BackPath
    sides::NTuple{2,HorizonSide}
    path::Vector{Int}

    function BackPath(robot, sides::NTuple{2,HorizonSide}=(Sud,West))
        local path = Int[]
        while !isborder(robot,sides[1]) || !isborder(robot,sides[2])
            for s in sides
                push!(path, steps_to_border!(robot, s))
            end
        end
        return new(reverse(reverse_side.(sides)), reverse(path))
    end
end

function back!(robot, backpath::BackPath)
    i=1
    for n in backpath.path
        movements!(robot, backpath.sides[i], n)
        i = i % 2 + 1
    end
end

num_steps(back_path::BackPath) = sum(back_path.path)

#back_path = BackPath(r)
#back!(r, back_path)

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
@testset "Sentiment" begin
    m = SentimentAnalyzer()

    d = StringDocument("a very nice thing that everyone likes")

    @test m(d) > 0.5

    d = StringDocument("a horrible thing that everyone hates")

    @test m(d) < 0.5
end


text:
########################
maxwell_boltzmann_primal(f) = f.* log.(f) - f

function computeEntropyH(α, m, ω,η::Function, η_prime_dual::Function)
    result = sum(η(η_prime_dual.(α' * m)[:] ).* ω )
    return result
end

#############################################

function ComputeSphericalBasisKarth(quadpts::Matrix{Float64}, polyDegree::Int64, spatialDim::Int64)
    monomialBasis = zeros(GetBasisSize(polyDegree,spatialDim),size(quadpts,1))

    if spatialDim == 3
        for idx_quad in 1:(size(quadpts)[1])
            monomialBasis[:,idx_quad]  = ComputeShpericalBasis3DPtKarth(quadpts[idx_quad,1],quadpts[idx_quad,2],quadpts[idx_quad,3],polyDegree)
        end
    elseif spatialDim == 1
        for idx_quad in 1:(size(quadpts)[1])
            monomialBasis[:,idx_quad]= ComputeShpericalBasis1DPtKarth(quadpts[idx_quad,1],polyDegree)
        end
    else
        println("dimension 2 currently not supported")
    end

    return monomialBasis

end

function ComputeShpericalBasis3DPtKarth(pointX::Float64,pointY::Float64,pointZ::Float64,polyDegree::Int64)
    idx_vector = 1
    spatialDim = 3
    basisLen = GetBasisSize(polyDegree,spatialDim)
    basisAtPt = ones(1,basisLen)
    for idx_degree in 0:polyDegree
        for a in 0:idx_degree
            for b in 0:(idx_degree-a)
                c = idx_degree - a - b
                basisAtPt[idx_vector] = Power(pointX,a)*Power(pointY,b)*Power(pointZ,c)
                idx_vector = idx_vector+1
            end
        end
    end

    return basisAtPt
end

function ComputeShpericalBasis1DPtKarth(pointX::Float64,polyDegree::Int64)
    idx_vector = 1
    spatialDim = 1
    basisLen = GetBasisSize(polyDegree,spatialDim)
    basisAtPt = ones(1,basisLen)
    for a in 0:polyDegree
                basisAtPt[idx_vector] = Power(pointX,a)
                idx_vector = idx_vector+1
    end

    return basisAtPt
end

function ComputeSphericalBasisAnalytical(quadpts::Matrix{Float64})
    # Hardcoded solution for L = 1, spatialDim = 3

    monomialBasis = zeros(4,size(quadpts,1))
    for idx_quad in 1:(size(quadpts)[1])
        monomialBasis[1,idx_quad]  = 1
        monomialBasis[2,idx_quad]  = quadpts[idx_quad,3] # z
        monomialBasis[3,idx_quad]  = quadpts[idx_quad,2] # y
        monomialBasis[4,idx_quad]  = quadpts[idx_quad,1] # x
    #    monomialBasis[5,idx_quad]  = quadpts[idx_quad,3] * quadpts[idx_quad,3] # z*z
    #    monomialBasis[6,idx_quad]  = quadpts[idx_quad,2] * quadpts[idx_quad,3] # y*z
    #    monomialBasis[7,idx_quad]  = quadpts[idx_quad,2] * quadpts[idx_quad,2] # y*y
    #    monomialBasis[8,idx_quad]  = quadpts[idx_quad,1] * quadpts[idx_quad,3] # x*z
    #    monomialBasis[9,idx_quad]  = quadpts[idx_quad,1] * quadpts[idx_quad,2] # x*y
    #    monomialBasis[10,idx_quad] = quadpts[idx_quad,1] * quadpts[idx_quad,1] # x*x
    end

    #monomialBasis = zeros(3,size(quadpts,1))
    #for idx_quad in 1:(size(quadpts)[1])
    #    monomialBasis[1,idx_quad]  = 1
    #    monomialBasis[2,idx_quad]  = quadpts[idx_quad,1] # z
    #    monomialBasis[3,idx_quad]  = quadpts[idx_quad,1] * quadpts[idx_quad,1] # z*z
    #end

    return monomialBasis
end

function Power(  basis,  exponent )
    if exponent == 0
        return 1.0
    end
    result = basis;
    for i in 2:exponent
        result = result * basis
    end
    return result
end

function GetBasisSize(LMaxDegree, spatialDim)
    basisLen = 0
    for idx_degree in 0:LMaxDegree
        basisLen = basisLen + GetCurrDegreeSize(idx_degree,spatialDim)
    end
    return convert(Int,basisLen)
end

function GetCurrDegreeSize(currDegree, spatialDim)
    return factorial(currDegree + spatialDim -1)/(factorial(currDegree)*factorial(spatialDim-1 ))
end


text:
##                                                                            ##
##              Problem set #2 for Structural modelling course                ##
##                               Nicolás Urdaneta                             ##
##                                                                            ##

clearconsole()
#=
Pkg.add("Optim")
Pkg.add("HTTP")
Pkg.add("GLM")
=#

using LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables, Distributions, Optim, HTTP, GLM

## 1. Basic optimization in Julia

f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
negf(x) =-copy(f(x))
startval = rand(1)
result = optimize(negf, startval, LBFGS())
println(result)

## 2. Optimization with data

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
	ssr = (y.-X*beta)'*(y.-X*beta)
	return ssr
end

beta_hat_ols = optimize(b-> ols(b,X,y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations = 100_000, show_trace = true))

println(beta_hat_ols.minimizer)

	## Just checking:

	bols = inv(X'X)*X'y
	df.white = df.race.==1
	bols_lm = lm(@formula(married ~ age+white+collgrad), df)

## 3. Optimization logit lokelihood

# Log likelihood
# y(lnP) + (1-y)ln(1-P)

function logit_ll(alpha, X, Y)

	ll_i = (y.==1).*log.(exp.(X*alpha)./(1 .+exp.(X*alpha))).+(y.==0).*log.(1 ./(1 .+exp.(X*alpha)))

	ll = -sum(ll_i)

	return ll

end

alpha_hat_optim = optimize(a -> logit_ll(a, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(alpha_hat_optim.minimizer)

## 4. Check logit using glm

alpha_hat_glm = glm(@formula(married ~ age + white + collgrad), df, Binomial(), LogitLink())
println(alpha_hat_glm)

## 5. Multinomial logit

	 # Join categories
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==8 ,:occupation] .= 7
df[df.occupation.==9 ,:occupation] .= 7
df[df.occupation.==10,:occupation] .= 7
df[df.occupation.==11,:occupation] .= 7
df[df.occupation.==12,:occupation] .= 7
df[df.occupation.==13,:occupation] .= 7
freqtable(df, :occupation) # problem solved

	# Redefine X and Y as we changed their rows
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function mlogit(alpha, X, Y)

	K = size(X,2)
	J = length(unique(y))
	N = size(X,1)

	bigY = zeros(N,J)
	for j=1:J
		bigY[:,j] = y.==j
	end

	bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

	num = zeros(N,J)
	dem = zeros(N)
	for j=1:J
		num[:,j] = exp.(X*bigAlpha[:,j])
		dem .+= num[:,j]
	end

	P = num ./ repeat(dem,1,J)

	ll = -sum( bigY.*log.(P))

end

alpha_zero = zeros(6*size(X,2))
alpha_rand = rand(6*size(X,2))
alpha_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]

alpha_start = alpha_true.*rand(size(alpha_true))
println(size(alpha_true))

alpha_hat_optim = optimize(a -> mlogit(a, X, y), alpha_start, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
alpha_hat_mle = alpha_hat_optim.minimizer
println(alpha_hat_mle)

## Stanrd erorrs

	## Adjust mlogit function
function mlogit_hessian(alpha, X, Y)

	K = size(X,2)
	J = length(unique(y))
	N = size(X,1)

	bigY = zeros(N,J)
	for j=1:J
		bigY[:,j] = y.==j
	end

	bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

	T = promote_type(eltype(X),eltype(alpha)) # this line is new
	num = zeros(T,N,J) ## Added T first
	dem = zeros(T,N) ## Added T first
	for j=1:J
		num[:,j] = exp.(X*bigAlpha[:,j])
		dem .+= num[:,j]
	end

	P = num ./ repeat(dem,1,J)

	ll = -sum( bigY.*log.(P))

end

# declare that the objective function is twice differentiable
td = TwiceDifferentiable(b -> mlogit_for_h(b, X, y), alpha_start; autodiff = :forward)
# run the optimizer
alpha_hat_optim_ad = optimize(td, alpha_zero, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
alpha_hat_mle_ad = alpha_hat_optim_ad.minimizer
# evaluate the Hessian at the estimates
H  = Optim.hessian!(td, alpha_hat_mle_ad)
# standard errors = sqrt(diag(inv(H))) [usually it's -H but we've already multiplied the obj fun by -1]
alpha_hat_mle_ad_se = sqrt.(diag(inv(H)))
println([alpha_hat_mle_ad alpha_hat_mle_ad_se]) # these standard errors match Stata


text:
module Thebes

using Luxor
using StaticArrays
using Rotations

include("Point3D.jl")
include("utils.jl")
include("Projection.jl")
include("Object.jl")
include("pin.jl")
include("text.jl")

export project, Projection, newprojection,
       Point3D,

       helloworld, eyepoint, centerpoint, uppoint,
       perspective,

       pin,

       hiddensurface, wireframe,

       Object, Teapot, Cube, Tetrahedron, Pyramid,

       carpet, drawcube, axes3D,
       between, make, distance, midpoint,

       rotateX, rotateY, rotateZ, rotateby!, rotateby,

       moveby!, moveby,

       objecttopoly,
       scaleby!, sortfaces!,
       sphericaltocartesian, cartesiantospherical,
       dotproduct3D, magnitude, anglebetweenvectors,
       surfacenormal, face, crossproduct3D, pointsperpendicular,

       text3D

const CURRENTPROJECTION = [newprojection(Point3D(100, 100, 100), Point3D(0, 0, 0), Point3D(0, 0, 1))]

end


text:
using CxxROOT
CxxROOT.LoadLibraryMap()
tcxx = icxx"""new TCxx("C++","Cxx.jl C++ Interpreter");"""
try
icxx"""
gROOT->GetClass("TGenPhaseSpace");
gROOT->GetClass("TH2F");
gROOT->GetClass("TCanvas");
gROOT->GetClass("TGHtml");
"""
end
app = icxx"""
new TApplication("ROOT Application", NULL, NULL);
"""
icxx"""
if (!gROOT->GetClass("TGenPhaseSpace")) gSystem->Load("libPhysics");
"""
h2 = icxx"""
    new TH2F("h2","h2", 50,1.1,1.8, 50,1.1,1.8);
"""
icxx"""
TLorentzVector target(0.0, 0.0, 0.0, 0.938);
TLorentzVector beam(0.0, 0.0, .65, .65);
TLorentzVector W = beam + target;

//(Momentum, Energy units are Gev/C, GeV)
Double_t masses[3] = { 0.938, 0.139, 0.139} ;

TGenPhaseSpace event;
event.SetDecay(W, 3, masses);

for (Int_t n=0;n<100000;n++) {
  Double_t weight = event.Generate();

  TLorentzVector *pProton = event.GetDecay(0);

  TLorentzVector *pPip    = event.GetDecay(1);
  TLorentzVector *pPim    = event.GetDecay(2);

  TLorentzVector pPPip = *pProton + *pPip;
  TLorentzVector pPPim = *pProton + *pPim;

  $h2->Fill(pPPip.M2() ,pPPim.M2() ,weight);
}
$h2->Draw();
"""
task = @async while true
sleep(0.1)
icxx"gSystem->ProcessEvents();"
end
c1 = icxx"""(TCanvas *)gROOT->GetListOfCanvases()->FindObject("c1");"""


text:
using Revise
using Plots
using Basins
using DynamicalSystems
using RegionTrees

# C. Grebogi, S. W. McDonald, E. Ott, J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983
ds = Systems.grebogi_map(rand(2))
integ  = integrator(ds)

θg=range(0,2π,length=1000)
xg=range(-1,1,length=1000)

@time bsn=Basins.basins_map2D_tree(θg, xg, integ; r_init=0.2, r_max=0.01)

function evaluate(basin,v)
    return findleaf(basin,v).data
end

function compute_frac_dim(bsn)
    root_n(x,n) = foldl((y,_) -> parent(y), 1:n; init=x)
    N = zeros(Int, 20)
    wd=1e3
    k_mx = 0
    for lf in allcells(bsn.basin)
        k = 1
        while !isnothing(root_n(lf,k))
            k += 1
        end
        N[k] += 1
        wd = min(wd,lf.boundary.widths[1])
        k_mx = max(k_mx,k)
    end
    @show N
    ε = reverse(wd * 2 .^ range(0,k_mx-1,step=1))
    return ε,N[1:k_mx]
end

ε,N = compute_frac_dim(bsn)

elog = log10.(1 ./ ε)
Nlog = log10.(N)
D = linear_region(elog[4:end],Nlog[4:end])
@show D

plt = plot(xlim=(θg[1], θg[end]), ylim=(xg[1], xg[end]), legend=nothing)
heatmap!(plt, θg, xg, (x, y) -> evaluate(bsn.basin, SVector(x, y)), fill=true)

#
#
# # plt = plot(xlim=(-1.5, 1.5), ylim=(-1.5, 1.5), legend=nothing)
#  for leaf in allleaves(bsn.basin)
#      v = hcat(collect(vertices(leaf.boundary))...)
#      #plot!(plt, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]])
#      #@show leaf.data
#      plot!(plt, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]], linecolor=:black, lw=0.4)
#
#  end
#  plt


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI   = point(1,5,3)
const pII  = point(2,5,3)
const pIII = point(7,1,0)

const qI  = point(10,11,12)
const qII  = point(10,11,13)
const qIII = point(11,11,12)

const P = simplex(pI,pII,pIII)
const Q = simplex(qI,qII,qIII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 16
pd_ref = PositiveDistance4D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
			return(exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(u,v)
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, pd_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm =[]
n1 = []
n2 = []
n3 = []
for i in 2:1:14
   Accuracy = i
   pd = PositiveDistance4D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, pd)

   num_pts = length(pd.qps)^4
   push!(n1,num_pts)
   push!(res_tp,int_tp)
end

for i in 2:1:7
   Accuracy = i
   pd_s = PositiveDistance4D_S(sing,SauterSchwab3D._shunnham2D(Accuracy))

   int_sp = sauterschwab_parameterized(INTEGRAND, pd_s)

   num_pts = length(pd_s.qps)^2
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   pd_gm = PositiveDistance4D_S(sing,SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1))

   int_gm = sauterschwab_parameterized(INTEGRAND, pd_gm)

   num_pts = length(pd_gm.qps)^2
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(5,1e4),ylims=(1e-9,1e-1))
plot!(xlabel="Quad. points/Func. evals.", ylabel="Rel. Error.", title="Positive Distance 4D",legend=:bottomleft)

using BenchmarkTools

ref = -0.0024840549158923397 - 0.005653396751321127im

pd = PositiveDistance4D(sing, SauterSchwab3D._legendre(4,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, pd)

num_pts = length(pd.qps)^4
println(num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println(err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, pd) end

pd_s = PositiveDistance4D_S(sing,SauterSchwab3D._shunnham2D(4))

   int_sp = sauterschwab_parameterized(INTEGRAND, pd_s)

   num_pts = length(pd_s.qps)^2
println(num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println(err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, pd_s) end
#=

pd_gm = PositiveDistance4D_S(SauterSchwab3D._grundmannMoeller2D(2*4-1))

int_gm = sauterschwab_parameterized(INTEGRAND, pd_gm)

num_pts = length(pd_gm.qps)^2
println(num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println(err_gm)
=#


text:
struct ZeroBlockInitializer end
const zero_blocks = ZeroBlockInitializer()

import Base.Cartesian: @nloops, @ntuple
import BlockArrays: BlockArray, undef_blocks, blockaxes
import FillArrays: Zeros, Fill

"""
Initialise a BlockArray where each block can have a different type.

    BlockArray{T}(zero_blocks, blocksize...)
"""
@generated function BlockArray{T}(::ZeroBlockInitializer, blocksizes::Vararg{AbstractVector{Int},N}) where {T,N}
    return quote
        B = BlockArray{T,N,Array{AbstractArray{T,N},N}}(undef_blocks, blocksizes...)
        axs = axes(B)
        @nloops $N block dim->blockaxes(B,dim) begin
            block_indices = @ntuple $N block
            indices = getindex.(axs, block_indices)
            block_size = length.(indices)
            B[block_indices...] = Zeros{T}(block_size...)
        end
        return B
    end
end

@generated function BlockArray{T}(::ZeroBlockInitializer, blocksizes::Vararg{AbstractVector{Int},N}) where {T <: Array,N}
    return quote
        B = BlockArray{T,N,Array{AbstractArray{T,N},N}}(undef_blocks, blocksizes...)
        axs = axes(B)
        @nloops $N block dim->blockaxes(B,dim) begin
            block_indices = @ntuple $N block
            indices = getindex.(axs, block_indices)
            block_size = length.(indices)
            B[block_indices...] = Fill{T}(T(), block_size...)
        end
        return B
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.58980225336104 <= q <= 16.332536463331582)

                     add_NL_constraint(m, :(swish(0.1097699546178803 + 0.5720816873107069*swish(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + 0.11808627755560508*swish(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.12730937872071246*swish(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.6166642606025916*swish(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + 0.45996993866072433*swish(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + swish(-0.8288330741680832 + -0.5580861219922526*swish(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.8786358511296881*swish(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.6915444540383735*swish(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.5554450106282389*swish(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.20937059223327026*swish(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + swish(0.13615042538554878 + 0.08502572267265318*swish(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.9980047609363689*swish(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.5678523898140919*swish(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.29618367003382495*swish(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.7256587322261439*swish(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + swish(-0.16139551949919406 + -0.5886735107938796*swish(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.02630736114123211*swish(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.03441572289730743*swish(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.3935924617025486*swish(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + 0.3425056157482871*swish(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + swish(-0.5329713534130147 + -0.3183971019485696*swish(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + 0.5847814320797045*swish(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.7330877632484647*swish(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + -0.37731199331975374*swish(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.08817423127621415*swish(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions