text:
"""
    function _rippa(A,a)

Estimate the Leave-One-Out (LOO) errors using rippas method. Complexity of 𝛰(3)
compared to calculating the exact LOO at a cost of 𝛰(4). `A` is the RBF matrix
and `a` is the weights of the RBF.
"""
function _rippa(A,a)
    N = size(A, 1)           #Number of Leave-One-Out (LOO) errors
    E = Array{Float64}(undef,N)

    _rippa!(E,A,a)
end

"""
    function _rippa!(E,A,a)

Same as _rippa but inplace.
"""
function _rippa!(E,A,a)
    N = size(A, 1)           #Number of Leave-One-Out (LOO) errors
    xₖ = A\I
    @inbounds for k = 1:N
        E[k] = a[k] / xₖ[k,k]   #Estimated error for k-th subset
    end

    return E
end

"""
    function _rmse_rippa!(A,a)

Same as _rippa! but calculates the RMSE value directly.
"""
function _rmse_rippa!(A,a)
    N = size(A, 1)           #Number of Leave-One-Out (LOO) errors
    xₖ = A\I
    RMSE = 0.0
    @inbounds for k = 1:N
        RMSE += (a[k] / xₖ[k,k])^2   #Accumulate square of estimated error for k-th subset
    end

    return sqrt(RMSE/N)
end

"""
    function RMSErrorLOO(interp,samples::Array{Float64,2},plan;
    cond_max = 1/eps(Float64)/10000, rippa = false)

Calculate the Leave-One-Out RMS error for a interpolation method.
`rippa` can be used to calculate the approximation of the LOO error
for Radial Basis Functions at a cost of 𝛰(3)
(compared to 𝛰(4)). `cond_max` sets the maximum allowed condition number for
matrix `A` used in the RBF calculation.
"""
function RMSErrorLOO(interp, samples, plan, smooth;
    cond_max::Float64=1e6, cond_check::Bool=false, rippa::Bool=false, rbf_dist_metric = Euclidean())

    #initiate arrays
    N = length(samples)
    E = Array{Float64}(undef,N)       #Error
    A = zeros(Float64, N, N)    #RBF matrix A
    ests = Array{Float64}(undef,N)    #Function estimate based on LOO
    LOOinds = Array{Int}(undef,N - 1, N)
    for i = 1:N
        LOOinds[:,i] = filter(x -> x != i, 1:N) # Get the leave one out sub indices
    end

    RMSErrorLOO!(A, ests, LOOinds, interp, samples, plan, smooth;
    cond_max=cond_max, rippa=rippa, rbf_dist_metric = rbf_dist_metric)
end

function RMSErrorLOO!(A, ests, LOOinds, interp::U,
 samples, plan::T, smooth; cond_max::Float64=1e6, cond_check::Bool=false,
  rippa::Bool=false, rbf_dist_metric = Euclidean()) where T <: AbstractArray where U <: AbstractArray

    N = length(samples)
    RMSE = Inf

    #perform Leave-One-Out estimation and calculate error
    if rippa
        @assert typeof(interp) <: Vector{ScatteredInterpolation.RadialBasisFunction} "Rippas
              algorithm only available for Radial Basis Functions"

        #interpolation object trained on the entire dataset
        itp, A = interpolate(interp, plan,
            samples, returnRBFmatrix=true, metric = rbf_dist_metric, smooth=smooth)

        #check the conditioning of matrix A for RBFs
        if cond_check && (cond(A) > cond_max)
            #RMSE is already Inf
        else
            #RBF error estimation based on Rippas algorithm
            RMSE = _rmse_rippa!(A,itp.w)
        end
    else
        for i = 1:N

            if typeof(smooth) <: AbstractVector
                loo_smooth = smooth[LOOinds[:,i]]
            else
                loo_smooth = smooth
            end

            #interpolation object trained on the LOO information
            itp, A = interpolate(interp[LOOinds[:,i]], plan[:,LOOinds[:,i]],
                samples[LOOinds[:,i]], returnRBFmatrix=true, metric = rbf_dist_metric, smooth=loo_smooth)

            #check the conditioning of matrix A for RBFs
            if cond_check && (cond(A) > cond_max)
                #RMSE is already Inf
            else
                #evaluate the interpolation object in the LOO position
                ests[i] = ScatteredInterpolation.evaluate(itp, plan[:,i])[1]
            end
        end
        RMSE = sqrt(mean((samples .- ests).^2))
    end

    return RMSE
end

function RMSErrorLOO!(A, ests, LOOinds, interp, samples, plan::T, smooth;
 cond_max::Float64=1e6, cond_check::Bool=false, rippa::Bool=false, rbf_dist_metric = Euclidean()) where T <: AbstractArray

    N = length(samples)
    RMSE = Inf

    #perform Leave-One-Out estimation and calculate error
    if rippa
        @assert typeof(interp) <: ScatteredInterpolation.RadialBasisFunction "Rippas algorithm only available for Radial Basis Functions"

        #interpolation object trained on the entire dataset
        itp, A = interpolate(interp, plan,
         samples, returnRBFmatrix=true, metric = rbf_dist_metric, smooth=smooth)

        #check the conditioning of matrix A for RBFs
        if cond_check && (cond(A) > cond_max)
            #RMSE is already Inf
        else
            #RBF error estimation based on Rippas algorithm
            RMSE = _rmse_rippa!(A,itp.w)
        end
    else
        for i = 1:N
            if typeof(smooth) <: AbstractVector
                loo_smooth = smooth[LOOinds[:,i]]
            else
                loo_smooth = smooth
            end

            #interpolation object trained on the LOO information
            itp, A = interpolate(interp, plan[:,LOOinds[:,i]],
                samples[LOOinds[:,i]], returnRBFmatrix=true, metric = rbf_dist_metric, smooth=loo_smooth)

            #check the conditioning of matrix A for RBFs
            if cond_check && (cond(A) > cond_max)
                #RMSE is already Inf
            else
                #evaluate the interpolation object in the LOO position
                ests[i] = ScatteredInterpolation.evaluate(itp, plan[:,i])[1]
            end
        end
        RMSE = sqrt(mean((samples .- ests).^2))
    end

    return RMSE
end

"""
    function scale(old_x::T,new_min,new_max;old_min=minimum(old_x),
    old_max=maximum(old_x)) where T <: Real

Scale a scalar to match a new range.
"""
function scale(old_x::T,new_min,new_max;old_min=minimum(old_x),old_max=maximum(old_x)) where T <: Real

    old_x = (((old_x - old_min) * (new_max - new_min)) / (old_max - old_min)) + new_min
end

"""
    function scale(old_x::Array{T,1},new_min,new_max;old_min=minimum(old_x),
    old_max=maximum(old_x)) where T <: Real

Scale a vector to match a new range.
"""
function scale(old_x::Array{T,1},new_min,new_max;old_min=minimum(old_x),old_max=maximum(old_x)) where T <: Real

    newX = (((old_x .- old_min) .* (new_max .- new_min)) ./ (old_max .- old_min)) .+ new_min
end

"""
    function scale(old_x::Array{T,2},direction::Int,new_min,new_max;
    old_min=minimum(old_x,direction),old_max=maximum(old_x,direction)) where T <: Real

Scale a 2D matrix to match a new range along the specified `direction`. Specify
the old extrema as scalars to use the same factor per row/column, or as a
vector to scale each row/column by a unique value.
"""
function scale(old_x::Array{T,2},direction::Int,new_min,new_max;old_min=minimum(old_x,dims=direction),old_max=maximum(old_x,dims=direction)) where T <: Real

    if length(old_min) == length(old_max) == 1
        newX = mapslices(x->scale(x,new_min,new_max;old_min=old_min,old_max=old_max), old_x, dims=direction)
    else
        old_min_itr = (x for x in old_min)
        old_max_itr = (x for x in old_max)

        mapfun = function (x)
            (omin, old_min_itr) = Iterators.peel(old_min_itr)
            (omax, old_max_itr) = Iterators.peel(old_max_itr)
            scale(x,new_min,new_max;old_min=omin,old_max=omax)
        end
        newX = mapslices(mapfun, old_x, dims=direction)
    end
    return newX
end

scale(x::Missing,min_val,max_val;old_min,old_max) = missing

"""
    preprocess_point(points,scaling;base_scale::Array{Float64,2})

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point(points,scaling;base_scale::Array{Float64,2})

    min_max = extrema(base_scale, dims=2)

    preprocessed_point = similar(points)
    for i = 1:size(preprocessed_point,1)
        preprocessed_point[i,:] = scale(points[i,:],-1.0*scaling[i],1.0*scaling[i],
        old_min = min_max[i][1], old_max = min_max[i][2])
    end
    return preprocessed_point
end

"""
    preprocess_point!(preprocessed_point,points,scaling;old_min::Array{Float64,2},old_max::Array{Float64,2})

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point!(preprocessed_point,points,scaling;old_min::Array{Float64,2},old_max::Array{Float64,2})

    for i = 1:size(preprocessed_point,1)
        preprocessed_point[i,:] = scale(points[i,:],-1.0*scaling[i],1.0*scaling[i],
        old_min = old_min[i], old_max = old_max[i])
    end
    return preprocessed_point
end

"""
    preprocess_point!(preprocessed_point,points,scaling::Bool;old_min::Array{Float64,2},old_max::Array{Float64,2}) where U

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point!(preprocessed_point,points,scaling::Bool;old_min::Array{Float64,2},old_max::Array{Float64,2}) where U
    preprocessed_point = points
end

"""
    preprocess_point(points,scaling::Bool;base_scale::Array{Float64,2}) where U

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point(points,scaling::Bool;base_scale::Array{Float64,2}) where U
    points
end

"""
    surrogate_evaluate(preprocessed_est_point,estimation_point,itp,optres,old_min,old_max)

Processes the input to scale correctly based on the options used.
"""
function surrogate_evaluate(preprocessed_est_point,estimation_point,itp,optres,old_min,old_max)
    preprocessed_est_point = preprocess_point!( preprocessed_est_point,
                                                estimation_point, optres.scaling;
                                                old_min=old_min,
                                                old_max=old_max,
                                                )

    # Evaluate the interpolation object
    res = ScatteredInterpolation.evaluate(itp, preprocessed_est_point)[1]
end

"""
    construct_search_range(plan::Array{Float64,2}, variable_kernel_width,
    min_rbf_width, max_rbf_width, variable_dim_scaling,
    min_scale, max_scale, smooth, max_smooth)

Creates array of tuples for the optimization of the RBF hyperparameters.
"""
function construct_search_range(plan::Array{Float64,2}, variable_kernel_width,
    min_rbf_width, max_rbf_width, variable_dim_scaling,
    min_scale, max_scale, smooth, max_smooth)

    n_dims, n_samples = size(plan)
    sr = Array{Tuple{Float64,Float64},1}()

    # Add the kernel width and type search range
    variable_kernel_width ? n_kerns = n_samples : n_kerns = 1
    push!(sr,   create_sr(  (min_range=min_rbf_width, max_range=max_rbf_width, n_times=n_kerns),
    (min_range=0.0, max_range=1.0, n_times=n_kerns))...)

    # Add the dimensional scaling search range
    variable_dim_scaling ? n_dim_scales = n_dims : n_dim_scales = 0
    push!(sr,   create_sr(  (min_range=min_scale, max_range=max_scale, n_times=n_dim_scales))...)

    # Add the ridge regression smoothing search range
    (smooth == :variable) && (n_smooth = n_samples)
    (smooth == :single)   && (n_smooth = 1)
    (smooth == :single_user)   && (n_smooth = 0)
    (smooth == false)     && (n_smooth = 0)
    push!(sr,   create_sr(  (min_range=0.0, max_range=max_smooth, n_times=n_smooth))...)

    return sr
end

"""
    create_sr(vargs::NamedTuple{(:min_range, :max_range, :n_times),Tuple{Float64,Float64,Int64}}...)

Facilitate the creation of array of tuples for the optimization of the RBF hyperparameters.
"""
function create_sr(vargs::NamedTuple{(:min_range, :max_range, :n_times),Tuple{Float64,Float64,Int64}}...)
    sr = Array{Tuple{Float64,Float64},1}()
    for varg in vargs
        for i = 1:varg[3]
            push!(sr,(varg[1],varg[2]))
        end
    end
    return sr
end

"""
    extract_vector_range(vargs::Int64...)

Supply length of subsets contained in a vector and receive a tuple containing
all the ranges needed to extract each subset in the vector.
"""
function extract_vector_range(vargs::Int...)
    output = Array{Union{UnitRange{Int}, Int, Bool},1}()

    count = 0
    for (i,varg) in enumerate(vargs)
        if varg == 0
            push!(output, false)
        else
            push!(output, (count+1):(count+varg))
            count += varg
            (length(output[end]) == 1) && (output[end] = output[end][1])
        end
    end
    return Tuple(output)
end

"""
    extract_bboptim_hypers(bboptim_fcall_vector,plan,kerns,
    variable_kernel_width,variable_dim_scaling,
    smooth,smooth_user)

Get the optimized RBF hyperparameters in useful format for further use.
"""
function extract_bboptim_hypers(bboptim_fcall_vector,plan,kerns,
    variable_kernel_width,variable_dim_scaling,
    smooth,smooth_user)

    n_dims, n_samples = size(plan)

    # Kernel width and type length
    variable_kernel_width ? n_kerns = n_samples : n_kerns = 1

    # Dimensional scaling length
    variable_dim_scaling ? n_dim_scales = n_dims : n_dim_scales = 0

    # Ridge regression smoothing length
    (smooth == :variable) && (n_smooth = n_samples)
    (smooth == :single)   && (n_smooth = 1)
    (smooth == :single_user)   && (n_smooth = 0)
    (smooth == false)     && (n_smooth = 0)

    width_inds, kernel_float_inds, scaling_inds, smooth_inds = extract_vector_range(n_kerns,
                                                            n_kerns,
                                                            n_dim_scales,
                                                            n_smooth
                                                            )

    # Arrange the width and smoothing results
    width = bboptim_fcall_vector[width_inds]
    kernel_float = bboptim_fcall_vector[kernel_float_inds]
    !(scaling_inds == false) ? scaling = bboptim_fcall_vector[scaling_inds] : scaling = scaling_inds
    (smooth == false) &&            (smooth = smooth_inds)
    (smooth == :single) &&          (smooth = bboptim_fcall_vector[smooth_inds])
    (smooth == :variable) &&        (smooth = bboptim_fcall_vector[smooth_inds])
    (smooth == :single_user) &&     (smooth = smooth_user)

    # Arrange the RBF kernel result
    kern_ind = round.(Int,scale(kernel_float,1,length(kerns),old_min=0,old_max=1))
    if variable_kernel_width
        kern = Vector{ScatteredInterpolation.RadialBasisFunction}(undef,size(plan,2))
        for i = 1:size(plan,2)
            kern[i] = kerns[kern_ind[i]](width[i])
        end
    elseif !variable_kernel_width
        kern = kerns[kern_ind](width)
    end

    return kern, scaling, smooth
end

"""
    interp_obj(inpt::Vector{Float64}, kerns, samples,
    plan::Array{Float64,2}; rippa::Bool = false,
    variable_kernel_width::Bool = true, variable_dim_scaling::Bool = true,
    smooth = false, cond_max=cond_max, rbf_dist_metric = Distances.Euclidean(),
    smooth_user::Float64 = 0.0)

Objective function for optimization of interpolation kernel function and width.
"""
function interp_obj(inpt::Vector{Float64}, A, ests, LOOinds, kerns, samples,
        plan::Array{Float64,2}; rippa::Bool = false,
        variable_kernel_width::Bool = true, variable_dim_scaling::Bool = true,
        smooth::Union{Bool, Symbol} = false, cond_max::Float64=Inf,  cond_check::Bool=false, rbf_dist_metric = Distances.Euclidean(),
        smooth_user::Float64 = 0.0)

    kern, scaling, smooth = extract_bboptim_hypers( inpt,plan,kerns,variable_kernel_width,
                                                    variable_dim_scaling,smooth,
                                                    smooth_user)

    #Preprocess the plan based on the settings used.
    preprocessed_plan = preprocess_point(plan,scaling,base_scale=plan)

    #Wrapped in try catch-block to catch failure to solve linear eq. system
    E = try
        E = RMSErrorLOO!(A, ests, LOOinds, kern, samples, preprocessed_plan, smooth; rippa = rippa,
        cond_max = cond_max, rbf_dist_metric = rbf_dist_metric, cond_check=cond_check)
    catch
        E = Inf
    end

    return E
end

"""
    rbf_hypers_opt(samples_org::Array{Float64,2}, plan::Array{Float64,2}, options::Options)

Optimization function of Radial Basis Function kernel and width.
"""
function rbf_hypers_opt(samples_org::Array{Float64,2}, plan::Array{Float64,2}, options::Options)

    @unpack rippa, variable_kernel_width, variable_dim_scaling, rbf_opt_method,
            min_rbf_width, max_rbf_width, min_scale, max_scale, cond_max,
            rbf_dist_metric, rbf_opt_gens, kerns, rbf_opt_pop,
            smooth, max_smooth, smooth_user, cond_check, constrained_seed_gens = options

    # Sample order for ScatteredInterpolation
    samples = vec(samples_org)

    # Pre-allocate for inplace functions
    N = length(samples) #Number of Leave-One-Out (LOO) errors
    A = zeros(Float64, N, N)    #RBF matrix A
    ests = Array{Float64}(undef,N)    #Function estimate based on LOO
    LOOinds = Array{Int}(undef,N - 1, N)
    for i = 1:N
        LOOinds[:,i] = filter(x -> x != i, 1:N) # Get the leave one out sub indices
    end

    # Create the hyperparameter search range based on the input options for a constrained space
    sr = construct_search_range(plan, false,
    min_rbf_width, max_rbf_width, variable_dim_scaling,
    min_scale, max_scale, smooth, max_smooth)

    #Improve the interpolant by first optimising the hyperparameters in a constrained search space
    if constrained_seed_gens > 0

        # RBF hyperparameter objective function
        itp_obj = function (x)
            interp_obj(x,A,ests,LOOinds,kerns,samples,plan;
                    rippa=rippa, variable_kernel_width=false,
                    variable_dim_scaling=variable_dim_scaling, smooth=smooth,
                    cond_max=cond_max,rbf_dist_metric=rbf_dist_metric,cond_check=cond_check)
        end

        @timeit_debug "hyper opt constrained" begin
        # Optimize the interpolant hyperparameters
        res = bboptimize(itp_obj;
                Method=rbf_opt_method,SearchRange=sr, MaxFuncEvals=rbf_opt_gens,
                TraceMode=:silent, rbf_dist_metric=rbf_dist_metric,
                TargetFitness = 1e-5, FitnessTolerance = 1e-6,
                PopulationSize = rbf_opt_pop,
                MaxStepsWithoutProgress=rbf_opt_gens,
                MaxNumStepsWithoutFuncEvals=rbf_opt_gens,
                );
        end
        #@show best_fitness(res)
        kern, scaling, smoothing = extract_bboptim_hypers( res.archive_output.best_candidate,
                                                        plan,kerns,false,
                                                        variable_dim_scaling,smooth,smooth_user)

        x0 = res.archive_output.best_candidate
    else
        x0 = BlackBoxOptim.Utils.latin_hypercube_sampling(first.(sr),last.(sr),1)
    end

    # Create the hyperparameter search range based on the input options
    sr = construct_search_range(plan, variable_kernel_width,
                                min_rbf_width, max_rbf_width, variable_dim_scaling,
                                min_scale, max_scale, smooth, max_smooth)

    # RBF hyperparameter objective function
    itp_obj = function (x)
        interp_obj(x,A,ests,LOOinds,kerns,samples,plan;
                rippa=rippa, variable_kernel_width=variable_kernel_width,
                variable_dim_scaling=variable_dim_scaling, smooth=smooth,
                cond_max=cond_max,rbf_dist_metric=rbf_dist_metric,cond_check=cond_check)
    end

    #First two values are the kern and kern_width, x0[3:end] is scaling and smooth
    #kern    kern_width    scaling    smooth
    if variable_kernel_width
        x1=[[x0[1] for _ in 1:size(plan,2)]; [x0[2] for _ in 1:size(plan,2)]; x0[3:end]]
    else
        x1 = x0
    end
    @assert length(x1)==length(sr)

    #Create population from LHC using the results from the constrained optimisation
    #of hyperparameters (if available), else supply random population.
    nvals=length(sr)
    if constrained_seed_gens > 0
        pop_constrained = [x1[i] for i in 1:nvals, _ in 1:1]
        pop = Array{Float64,2}(undef,nvals,rbf_opt_pop)
        pop[:,1] = pop_constrained
        for i in size(pop_constrained,2)+1:size(pop,2)
            pop[:,i] = BlackBoxOptim.Utils.latin_hypercube_sampling(first.(sr),last.(sr),1)
        end
    else
        pop = BlackBoxOptim.Utils.latin_hypercube_sampling(first.(sr),last.(sr),rbf_opt_pop)
    end

    @timeit_debug "hyper opt full" begin
    # Optimize the interpolant hyperparameters
    res = bboptimize(itp_obj;
            Method=rbf_opt_method,SearchRange=sr, MaxFuncEvals=rbf_opt_gens,
            TraceMode=:silent, rbf_dist_metric=rbf_dist_metric,
            TargetFitness = 1e-5, FitnessTolerance = 1e-6,
            PopulationSize = rbf_opt_pop,
            MaxStepsWithoutProgress=rbf_opt_gens,
            MaxNumStepsWithoutFuncEvals=rbf_opt_gens,
            Population=pop
            );
    end
    #@show best_fitness(res)
    kern, scaling, smoothing = extract_bboptim_hypers( res.archive_output.best_candidate,
                                                    plan,kerns,variable_kernel_width,
                                                    variable_dim_scaling,smooth,smooth_user)

    # Return the optimized hyperparameters in the correct type
    return RBFHypers(kern, scaling, smoothing), best_fitness(res)
end


text:
#=
MurmurHash3 was written by Austin Appleby, and is placed in the public
domain.
=#
const c1 = 0xcc9e2d51
const c2 = 0x1b873593

function _murmur32(data::Union{AbstractString,AbstractVector{UInt8}}, seed::UInt32)::UInt32
    # head
    len = UInt32(length(data))
    nblocks = div(len, 4)
    last = nblocks*4
    h = seed
    blocks = Ptr{UInt32}(pointer(data))

    # body
    @inbounds for block = 1:nblocks
        k = unsafe_load(blocks, block)
        k *= c1; k = rotl32(k, 15); k *= c2; h = h ⊻ k
        h = rotl32(h, 13); h *= UInt32(5); h += 0xe6546b64
    end

    # tail
    k = UInt32(0)
    remainder = len & 3
    tail = Ptr{UInt8}(pointer(data))

    if remainder == 3 k = k ⊻ UInt32(unsafe_load(tail, last+3)) << 16 end
    if remainder >= 2 k = k ⊻ UInt32(unsafe_load(tail, last+2)) <<  8 end
    if remainder >= 1
        k = k ⊻ unsafe_load(tail, last+1)
        k *= c1; k = rotl32(k, 15); k *= c2; h = h ⊻ k
    end

    # finalization
    h = h ⊻ len
    h = fmix32(h)
end

function _check_seed(k::Union{AbstractString,AbstractVector{UInt8}}, s::Integer)
    try
        _murmur32(k, UInt32(s))
    catch InexactError
        throw(DomainError(s, "`seed` must be in [0..$(typemax(UInt32))]"))
    end
end

murmur32(k::AbstractVector{UInt8}, s::Integer)    = _check_seed(k, s)
murmur32(k::AbstractVector{UInt8})                = _murmur32(k, UInt32(0))
murmur32(k::AbstractString, s::Integer)           = _check_seed(k, s)
murmur32(k::AbstractString)                       = _murmur32(k, UInt32(0))
murmur32(k::AbstractVector{<:Number}, s::Integer) = _check_seed(reinterpret(UInt8, k), s)
murmur32(k::AbstractVector{<:Number})             = murmur32(reinterpret(UInt8, k))
murmur32(k::AbstractArray{<:Number}, s::Integer)  = murmur32(vec(k), s)
murmur32(k::AbstractArray{<:Number})              = murmur32(vec(k))

"""
    murmur32(key::AbstractArray{<:Number}[, seed])
    murmur32(key::AbstractString[, seed])

Hash `key` to a `UInt32` using the 32 bit MurmurHash3 hashing function, optionally
seeded with `seed`.
"""
murmur32

@inline rotl32(x::UInt32, r::Integer)::UInt32 = (x << r) | (x >> (32 - r))

@inline function fmix32(h::UInt32)::UInt32
    h = h ⊻ h >> 16
    h *= 0x85ebca6b
    h = h ⊻ h >> 13
    h *= 0xc2b2ae35
    h = h ⊻ h >> 16
end

ext: jl


text:
@doc raw"""
    RationalQuadraticKernel([α=1 [,β=1]])

The rational-quadratic kernel is given by:

```math
\kappa(\mathbf{x},\mathbf{y})
= \left(1 +\alpha ||\mathbf{x},\mathbf{y}||^2\right)^{-\beta}
\qquad \alpha > 0, \; \beta > 0
```

where ``\alpha`` is a scaling parameter and ``\beta`` is a shape parameter. This kernel can
be seen as an infinite sum of Gaussian kernels. If one sets ``\alpha = \alpha_0 / \beta``,
then taking the limit ``\beta \rightarrow \infty`` results in the Gaussian kernel with
scaling parameter ``\alpha_0``.
"""
struct RationalQuadraticKernel{T<:AbstractFloat} <: MercerKernel{T}
    alpha::HyperParameter{T}
    beta::HyperParameter{T}
    RationalQuadraticKernel{T}(α::Real, β::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
        HyperParameter(convert(T,β), interval(OpenBound(zero(T)), nothing))
    )
end
function RationalQuadraticKernel(α::T1 = 1.0, β::T2 = one(T1)) where {T1<:Real,T2<:Real}
    RationalQuadraticKernel{floattype(T1,T2)}(α, β)
end

@inline rationalquadratickernel(z::T, α::T, β::T) where {T<:AbstractFloat} = (1 + α*z)^(-β)

@inline basefunction(::RationalQuadraticKernel) = SquaredEuclidean()
@inline function kappa(κ::RationalQuadraticKernel{T}, z::T) where {T}
    rationalquadratickernel(z, getvalue(κ.alpha), getvalue(κ.beta))
end


text:
# List of analysis for performance evaluation of search processes
# Ideally short functions

"""primal_integral(v, t; )

Compute the primal integral
"""
function primal_integral(v, t; BKS = 0, T = 1000)
    return 100*(sum(v[i-1](t[i] - t[i-1]) + v[end]*T-t[end], 1:length(v)) / (T * BKS) - 1)
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
abstract type CNFLayer <: Function end
Flux.trainable(m::CNFLayer) = (m.p,)

"""
Constructs a continuous-time recurrent neural network, also known as a neural
ordinary differential equation (neural ODE), with fast gradient calculation
via adjoints [1] and specialized for density estimation based on continuous
normalizing flows (CNF) [2] with a direct computation of the trace
of the dynamics' jacobian. At a high level this corresponds to the following steps:

1. Parameterize the variable of interest x(t) as a function f(z,θ,t) of a base variable z(t) with known density p_z;
2. Use the transformation of variables formula to predict the density p_x as a function of the density p_z and the trace of the Jacobian of f;
3. Choose the parameter θ to minimize a loss function of p_x (usually the negative likelihood of the data);

After these steps one may use the NN model and the learned θ to predict the density p_x for new values of x.

```julia
DeterministicCNF(model,tspan,basedist=nothing,monte_carlo=false,args...;kwargs...)
```
Arguments:
- `model`: A Chain neural network that defines the dynamics of the model.
- `basedist`: Distribution of the base variable. Set to the unit normal by default.
- `tspan`: The timespan to be solved on.
- `kwargs`: Additional arguments splatted to the ODE solver. See the
  [Common Solver Arguments](https://diffeq.sciml.ai/dev/basics/common_solver_opts/)
  documentation for more details.
Ref
[1]L. S. Pontryagin, Mathematical Theory of Optimal Processes. CRC Press, 1987.
[2]R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud. Neural Ordinary Differential Equations. arXiv preprint at arXiv1806.07366, 2019.
[3]W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, D. Duvenaud. FFJORD: Free-Form Continuous Dynamic For Scalable Reversible Generative Models. arXiv preprint at arXiv1810.01367, 2018.

"""
struct DeterministicCNF{M,P,RE,Distribution,T,A,K} <: CNFLayer
    model::M
    p::P
    re::RE
    basedist::Distribution
    tspan::T
    args::A
    kwargs::K

    function DeterministicCNF(model,tspan,args...;p = nothing,basedist=nothing,kwargs...)
        _p,re = Flux.destructure(model)
        if p === nothing
            p = _p
        end
        if basedist === nothing
            size_input = size(model[1].W)[2]
            basedist = MvNormal(zeros(size_input), I + zeros(size_input,size_input))
        end
        new{typeof(model),typeof(p),typeof(re),typeof(basedist),typeof(tspan),typeof(args),typeof(kwargs)}(
            model,p,re,basedist,tspan,args,kwargs)
    end
end

"""
Constructs a continuous-time recurrent neural network, also known as a neural
ordinary differential equation (neural ODE), with fast gradient calculation
via adjoints [1] and specialized for density estimation based on continuous
normalizing flows (CNF) [2] with a stochastic approach [2] for the computation of the trace
of the dynamics' jacobian. At a high level this corresponds to the following steps:

1. Parameterize the variable of interest x(t) as a function f(z,θ,t) of a base variable z(t) with known density p_z;
2. Use the transformation of variables formula to predict the density p_x as a function of the density p_z and the trace of the Jacobian of f;
3. Choose the parameter θ to minimize a loss function of p_x (usually the negative likelihood of the data);

After these steps one may use the NN model and the learned θ to predict the density p_x for new values of x.

```julia
FFJORD(model,basedist=nothing,monte_carlo=false,tspan,args...;kwargs...)
```
Arguments:
- `model`: A Chain neural network that defines the dynamics of the model.
- `basedist`: Distribution of the base variable. Set to the unit normal by default.
- `tspan`: The timespan to be solved on.
- `kwargs`: Additional arguments splatted to the ODE solver. See the
  [Common Solver Arguments](https://diffeq.sciml.ai/dev/basics/common_solver_opts/)
  documentation for more details.
Ref
[1]L. S. Pontryagin, Mathematical Theory of Optimal Processes. CRC Press, 1987.
[2]R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud. Neural Ordinary Differential Equations. arXiv preprint at arXiv1806.07366, 2019.
[3]W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, D. Duvenaud. FFJORD: Free-Form Continuous Dynamic For Scalable Reversible Generative Models. arXiv preprint at arXiv1810.01367, 2018.

"""
struct FFJORD{M,P,RE,Distribution,T,A,K} <: CNFLayer
    model::M
    p::P
    re::RE
    basedist::Distribution
    tspan::T
    args::A
    kwargs::K

    function FFJORD(model,tspan,args...;p = nothing,basedist=nothing,kwargs...)
        _p,re = Flux.destructure(model)
        if p === nothing
            p = _p
        end
        if basedist === nothing
            size_input = size(model[1].W)[2]
            basedist = MvNormal(zeros(size_input), I + zeros(size_input,size_input))
        end
        new{typeof(model),typeof(p),typeof(re),typeof(basedist),typeof(tspan),typeof(args),typeof(kwargs)}(
            model,p,re,basedist,tspan,args,kwargs)
    end
end

function jacobian_fn(f, x::AbstractVector)
   y::AbstractVector, back = Zygote.pullback(f, x)
   ȳ(i) = [i == j for j = 1:length(y)]
   vcat([transpose(back(ȳ(i))[1]) for i = 1:length(y)]...)
 end

function cnf(du,u,p,t,re)
    z = @view u[1:end-1]
    m = re(p)
    J = jacobian_fn(m, z)
    trace_jac = length(z) == 1 ? sum(J) : tr(J)
    du[1:end-1] = m(z)
    du[end] = -trace_jac
end

function ffjord(du,u,p,t,re,e,monte_carlo,regularize)
    m = re(p)
    if regularize
        z = @view u[1:end-3]
        _, back = Zygote.pullback(m,z)
        eJ = back(e)[1]
        if monte_carlo
            trace_jac = (eJ.*e)[1]
        else
            J = jacobian_fn(m, z)
            trace_jac = length(z) == 1 ? sum(J) : tr(J)
        end
        du[1:end-3] = m(z)
        du[end-2] = -trace_jac
        du[end-1] = sum(abs2, m(z))
        du[end] = norm(eJ)^2
    else
        z = @view u[1:end-1]
        _, back = Zygote.pullback(m,z)
        eJ = back(e)[1]
        if monte_carlo
            trace_jac = (eJ.*e)[1]
        else
            J = jacobian_fn(m, z)
            trace_jac = length(z) == 1 ? sum(J) : tr(J)
        end
        du[1:end-1] = m(z)
        du[end] = -trace_jac
    end
end

function (n::DeterministicCNF)(x,p=n.p)
    cnf_ = (du,u,p,t)->cnf(du,u,p,t,n.re)
    prob = ODEProblem{true}(cnf_,vcat(x,0f0),n.tspan,p)
    sense = InterpolatingAdjoint(autojacvec = false)
    pred = solve(prob,n.args...;sensealg=sense,n.kwargs...)[:,end]
    pz = n.basedist
    z = pred[1:end-1]
    delta_logp = pred[end]
    logpz = logpdf(pz, z)
    logpx = logpz .- delta_logp
    return logpx[1]
end

function (n::FFJORD)(x,p=n.p,regularize=false,monte_carlo=true)
    e = randn(Float32,length(x))
    pz = n.basedist
    sense = InterpolatingAdjoint(autojacvec = false)
    if regularize
        ffjord_ = (du,u,p,t)->ffjord(du,u,p,t,n.re,e,monte_carlo,regularize)
        prob = ODEProblem{true}(ffjord_,vcat(x,0f0,0f0,0f0),n.tspan,p)
        pred = solve(prob,n.args...;sensealg=sense,n.kwargs...)[:,end]
        z = pred[1:end-3]
        delta_logp = pred[end-2]
        λ₁ = pred[end-1]
        λ₂ =  pred[end]
        logpz = logpdf(pz, z)
        logpx = logpz .- delta_logp
        return logpx[1], λ₁, λ₂
    else
        ffjord_ = (du,u,p,t)->ffjord(du,u,p,t,n.re,e,monte_carlo,regularize)
        prob = ODEProblem{true}(ffjord_,vcat(x,0f0),n.tspan,p)
        pred = solve(prob,n.args...;sensealg=sense,n.kwargs...)[:,end]
        z = pred[1:end-1]
        delta_logp = pred[end]
        logpz = logpdf(pz, z)
        logpx = logpz .- delta_logp
        return logpx[1]
    end
end


text:
"""
    StateVecs

Provides a state vector that is includes information
of the domain decomposition into a grid mean,
environment and updrafts.
"""
module StateVecs

using ..FiniteDifferenceGrids

# Forwarded methods
import ..DomainDecompositions:
    var_names,
    DomainIdx,
    allcombinations,
    gridmean,
    environment,
    updraft,
    over_sub_domains,
    alldomains,
    get_i_var,
    var_string,
    var_suffix
using ..DomainDecompositions

export StateVec, over_sub_domains, Cut, Dual, allcombinations
export var_names, var_string, var_suffix
export assign!, assign_real!, assign_ghost!, extrap!, extrap_0th_order!
export compare

struct FieldsPerElement
    vals::Vector
end

"""
    StateVec{VM}

A state vector containing the number of subdomains,
`n_subdomains`, a `NamedTuple` variable mapper, a tuple
of the variable names, and a vector of vectors, containing
the values for all of the variables.
"""
struct StateVec{VM}
    """
    A `DomainDecomposition`, containing the map between subdomains variables and the state vector index
    """
    domain_decomp::DomainDecomposition{VM}
    """
    A vector of vectors containing all values for all variables.
    """
    fields::Vector{FieldsPerElement}
end

"""
    StateVec(vars::Tuple{Vararg{Tuple{Symbol, DD}}}, grid::Grid{T}) where {DD, T}

Return a state vector, given a tuple of tuples of variable
names and the number of their subdomains.
"""
function StateVec(vars, grid::Grid{T}, domain_set::DomainSet) where {T}
    domain_decomp = DomainDecomposition(vars, domain_set)
    n_vars = sum([sum(domain_set, dss) for (v, dss) in vars])
    all_vars = Any[eltype(grid.Δz)(0) for v in 1:n_vars]
    fields = [FieldsPerElement(deepcopy(all_vars)) for k in over_elems(grid)]
    return StateVec{typeof(domain_decomp.var_mapper)}(domain_decomp, fields)
end

function Base.show(io::IO, sv::StateVec)
    println(io, "\n----------------------- State Vector")
    println(io, sv.domain_decomp)
    vn = var_names(sv)
    headers = [
        length(over_sub_domains(sv, ϕ)) == 1 ? string(ϕ) :
            string(ϕ) * "_" * string(i) for ϕ in vn
        for i in over_sub_domains(sv, ϕ)
    ]
    n_vars = length(headers)
    n_elem = length(sv.fields)
    data = reshape(
        [
            sv[ϕ, k, i] for ϕ in vn for i in over_sub_domains(sv, ϕ)
            for k in 1:n_elem
        ],
        n_elem,
        n_vars,
    )
    Base.print_matrix(io, [reshape(Text.(headers), 1, n_vars); data])
    println(io, "\n-----------------------")
end

# Overload
var_names(sv::StateVec) = var_names(sv.domain_decomp)
alldomains(sv::StateVec) = alldomains(sv.domain_decomp)
DomainIdx(sv::StateVec) = DomainIdx(sv.domain_decomp)
allcombinations(sv::StateVec) = allcombinations(sv.domain_decomp)
gridmean(sv::StateVec) = gridmean(sv.domain_decomp)
environment(sv::StateVec) = environment(sv.domain_decomp)
updraft(sv::StateVec) = updraft(sv.domain_decomp)
get_i_var(sv::StateVec, ϕ::Symbol, i) = get_i_var(sv.domain_decomp, ϕ, i)

"""
    over_sub_domains(sv::StateVec, ϕ::Symbol)

Get list of indexes over all subdomains for variable `ϕ`.
"""
over_sub_domains(sv::StateVec, ϕ::Symbol) =
    over_sub_domains(sv.domain_decomp, ϕ)

function Base.getindex(sv::StateVec, ϕ::Symbol, k, i = 0)
    i == 0 && (i = gridmean(sv))
    i_sv = DomainDecompositions.get_i_state_vec(sv.domain_decomp, ϕ, i)
    @inbounds sv.fields[k].vals[i_sv]
end
function Base.setindex!(sv::StateVec, val, ϕ::Symbol, k, i = 0)
    i == 0 && (i = gridmean(sv))
    @inbounds i_sv =
        DomainDecompositions.get_i_state_vec(sv.domain_decomp, ϕ, i)
    @inbounds sv.fields[k].vals[i_sv] = val
end

abstract type AbstractCut{I} end

"""
    Cut{I} <: AbstractCut{I}

A Cut struct used to slice the state
vector along the grid-element dimension.
This is used as an API to pass Cuts
into local derivative/interpolation routines.
"""
struct Cut{I} <: AbstractCut{I}
    e::I
end
function Base.getindex(sv::StateVec, ϕ::Symbol, cut::Cut, i = 0)
    @inbounds [sv[ϕ, k, i] for k in (cut.e - 1):(cut.e + 1)]
end

"""
    Dual{I} <: AbstractCut{I}

A Dual struct used to slice the state
vector along the grid-element dimension.
In addition, this interpolates fields
[(f[k-1]+f[k])/2, (f[k]+f[k+1])/2]
"""
struct Dual{I} <: AbstractCut{I}
    e::I
end
function Base.getindex(sv::StateVec, ϕ::Symbol, dual::Dual, i = 0)
    @inbounds [
        (sv[ϕ, k, i] + sv[ϕ, k + 1, i]) / 2 for k in (dual.e - 1):(dual.e)
    ]
end

Base.isnan(sv::StateVec) = any([any(isnan.(fpe.vals)) for fpe in sv.fields])
Base.isinf(sv::StateVec) = any([any(isinf.(fpe.vals)) for fpe in sv.fields])

"""
    var_suffix(sv::StateVec, ϕ::Symbol, i=0)

A suffix string for variable `ϕ` indicating its sub-domain.
"""
var_suffix(sv::StateVec, ϕ::Symbol, i = 0) = var_suffix(sv.domain_decomp, ϕ, i)

"""
    var_string(sv::StateVec, ϕ::Symbol, i=0)

A string of the variable `ϕ` with the appropriate suffix indicating its sub-domain.
"""
var_string(sv::StateVec, ϕ::Symbol, i = 0) = var_string(sv.domain_decomp, ϕ, i)

"""
    assign!(sv::StateVec, var_names, grid::Grid, val, i=1)

Assign value `val` to variable `ϕ` for all ghost points.
"""
function assign!(sv::StateVec, var_names, grid::Grid{FT}, val::FT) where {FT}
    gm, en, ud, sd, al = allcombinations(sv)
    !(var_names isa Tuple) && (var_names = (var_names,))
    @inbounds for k in over_elems(grid),
        ϕ in var_names,
        i in over_sub_domains(sv, ϕ)

        sv[ϕ, k, i] = val
    end
end

"""
    assign!(sv::StateVec, grid::Grid, val)

Assign value `val` to all variables in state vector.
"""
function assign!(sv::StateVec, grid::Grid{FT}, val::FT) where {FT}
    gm, en, ud, sd, al = allcombinations(sv)
    @inbounds for k in over_elems(grid),
        ϕ in var_names(sv),
        i in over_sub_domains(sv, ϕ)

        sv[ϕ, k, i] = val
    end
end

"""
    assign_real!(sv::StateVec, ϕ::Symbol, grid::Grid, vec_real, i=0)

Assign the real elements to the given vector `vec_real`.
"""
function assign_real!(sv::StateVec, ϕ::Symbol, grid::Grid, vec_real, i = 0)
    i == 0 && (i = gridmean(sv))
    k_real = 1
    @assert length(vec_real) == length(over_elems_real(grid))
    @inbounds for k in over_elems_real(grid)
        sv[ϕ, k, i] = vec_real[k_real]
        k_real += 1
    end
end

"""
    assign_ghost!(sv::StateVec, ϕ::Symbol, grid::Grid, val, i=0)

Assign value `val` to variable `ϕ` for all ghost points.
"""
function assign_ghost!(sv::StateVec, ϕ::Symbol, grid::Grid, val, i = 0)
    i == 0 && (i = gridmean(sv))
    @inbounds for k in over_elems_ghost(grid)
        sv[ϕ, k, i] = val
    end
end

"""
    assign_ghost!(sv::StateVec, ϕ::Symbol, grid::Grid, val, b::ZBoundary, i=1)

Assign value `val` to variable `ϕ` for all ghost points.
"""
function assign_ghost!(
    sv::StateVec,
    ϕ::Symbol,
    grid::Grid,
    val,
    b::ZBoundary,
    i = 0,
)
    i == 0 && (i = gridmean(sv))
    @inbounds for k in over_elems_ghost(grid, b)
        sv[ϕ, k, i] = val
    end
end

"""
    extrap!(sv::StateVec, ϕ::Symbol, grid::Grid, i=1)

Extrapolate variable `ϕ` to the first ghost point.
"""
function extrap!(sv::StateVec, ϕ::Symbol, grid::Grid, i = 0)
    i == 0 && (i = gridmean(sv))
    @inbounds for b in (Zmin(), Zmax())
        kg, ki, kii = boundary_points(grid, b)
        sv[ϕ, kg, i] = 2 * sv[ϕ, ki, i] - sv[ϕ, kii, i]
    end
end

"""
    extrap!(sv::StateVec, ϕ::Symbol, grid::Grid, ::ZBoundary, i=1)

Extrapolate variable `ϕ` to the first ghost point.
"""
function extrap!(sv::StateVec, ϕ::Symbol, grid::Grid, b::ZBoundary, i = 0)
    i == 0 && (i = gridmean(sv))
    kg, ki, kii = boundary_points(grid, b)
    @inbounds sv[ϕ, kg, i] = 2 * sv[ϕ, ki, i] - sv[ϕ, kii, i]
end

"""
    extrap!(sv::StateVec, ϕ::Symbol, grid::Grid, dual_val, b::ZBoundary, i=1)

Extrapolate variable `ϕ` to the first ghost point.
"""
function extrap!(
    sv::StateVec,
    ϕ::Symbol,
    grid::Grid,
    dual_val,
    b::ZBoundary,
    i = 0,
)
    i == 0 && (i = gridmean(sv))
    ki = first_interior(grid, b)
    kg = first_ghost(grid, b)
    @inbounds sv[ϕ, kg, i] = 2 * dual_val - sv[ϕ, ki, i]
end

"""
    extrap_0th_order!(sv::StateVec, ϕ::Symbol, grid::Grid, i=1)

Extrapolate variable `ϕ` to the first ghost point using zeroth order approximation.
"""
function extrap_0th_order!(sv::StateVec, var_names, grid::Grid, i = 0)
    !(var_names isa Tuple) && (var_names = (var_names,))
    i == 0 && (i = gridmean(sv))
    @inbounds for ϕ in var_names
        @inbounds for b in (Zmin(), Zmax())
            kg, ki, kii = boundary_points(grid, b)
            sv[ϕ, kg, i] = sv[ϕ, ki, i]
        end
    end
end

include("boundary_conditions.jl")
include("export_funcs.jl")
include("domain_decomposition_funcs.jl")

end


text:
using Test

using Mill
using Mill: nobs, mapdata
using Mill: BagConv, convsum, bagconv, legacy_bagconv, _convshift, ∇wbagconv, ∇xbagconv, ∇convsum
using Mill: ngrams, countngrams
using Mill: p_map, inv_p_map, r_map, inv_r_map, _bagnorm
using Mill: Maybe
using Mill: gradtest, gradf

using Base.Iterators: partition, product
using Base: CodeUnits
using ChainRulesCore
using ChainRulesCore: NotImplementedException
using Combinatorics
using DataFrames
using Documenter
using Flux
using Flux: onehot, onehotbatch
using HierarchicalUtils
using PooledArrays
using Random
using SparseArrays

using BenchmarkTools: @btime

areequal(x) = true
areequal(x, y, zs...) = isequal(x, y) && areequal(y, zs...)

Random.seed!(25)

const BAGS = [
      length2bags([1, 2, 3]),
      length2bags([4, 1, 1]),
      length2bags([2, 2, 2])
]

const BAGS2 = [
        length2bags(ones(Int, 10)),
        length2bags(2 .* ones(Int, 5)),
        length2bags([5, 5]),
        length2bags([10]),
        length2bags([3, 4, 3]),
        AlignedBags([1:3, 0:-1, 0:-1, 4:7, 0:-1, 8:10]),
        AlignedBags([0:-1, 1:5, 0:-1, 0:-1, 0:-1, 6:10]),
        ScatteredBags([collect(1:3), collect(7:10), collect(4:6)]),
        ScatteredBags([collect(7:10), Int[], collect(1:3), Int[], collect(4:6), Int[]]),
        ScatteredBags([Int[], collect(1:10), Int[]]),
]

const BAGS3 = [
         (AlignedBags([1:2, 3:4, 0:-1]),
          ScatteredBags([[2, 3, 4], [1], Int[]]),
          AlignedBags([1:4, 0:-1, 5:8, 0:-1])),
         (AlignedBags([0:-1, 1:2, 3:4]),
          ScatteredBags([[1], [2], [3, 4]]),
          AlignedBags([0:-1, 1:7, 0:-1, 8:8])),
         (AlignedBags([0:-1, 0:-1, 1:2, 3:4]),
          ScatteredBags([[2, 4], Int[], [3, 1], Int[]]),
          AlignedBags([1:1, 2:2, 0:-1, 3:8])),
         (AlignedBags([0:-1, 1:2, 3:4, 0:-1]),
          ScatteredBags([Int[], [1, 3], [2, 4], Int[]]),
          AlignedBags([0:-1, 1:2, 3:6, 7:8]))
        ]

function Mill.unpack2mill(ds::LazyNode{:Sentence})
    s = split.(ds.data, " ")
    x = NGramMatrix(reduce(vcat, s))
    BagNode(ArrayNode(x), Mill.length2bags(length.(s)))
end

_init_agg(t::Type{<:Number}, d) = randn(t, d)
_init_agg(t::Type{<:Integer}, d) = rand(t, d)

# initialize to randn to test a lot of values even though we then initialize to zero
nonparam_aggregations(t::Type, d) = vcat(
                                 SegmentedSum(_init_agg(t, d)),
                                 SegmentedMean(_init_agg(t, d)),
                                 SegmentedMax(_init_agg(t, d)))

param_aggregations(t::Type, d) = vcat(
                                      SegmentedPNorm(_init_agg(t, d), _init_agg(t, d), _init_agg(t, d)),
                                      SegmentedLSE(_init_agg(t, d), _init_agg(t, d)))

all_aggregations(t::Type, d) = vcat(nonparam_aggregations(t, d), param_aggregations(t, d))

@testset "Doctests" begin
    DocMeta.setdocmeta!(Mill, :DocTestSetup, quote
        using Mill, Flux, Random, SparseArrays, Setfield, HierarchicalUtils
    end; recursive=true)
    doctest(Mill)
end

for test_f in readdir(".")
    (endswith(test_f, ".jl") && test_f != "runtests.jl") || continue
    @testset verbose = true "$test_f" begin
        include(test_f)
    end
end


text:
using BenchmarkTools

# see this section of docs https://docs.julialang.org/en/v1/base/arrays/#Broadcast-and-vectorization-1

function G_single(t::T) where T<:Real

   G = 1. + 2. *exp(-t*(2. /3.)) + 4. *exp(-t*(4. /5.))

end

function G_vec(t::Vector{T}) where T<:Real

   G = 1. .+ 2. .* exp.(-t .* (2. /3.)) .+ 4. .* exp.(-t .* (4. /5.))

end

function G_arraycomp(t::Vector{T}) where T<:Real
    [G_single(ti) for ti in t]
end

function G_loop(t::Vector{T}) where T<:Real

    out = Array{T}(undef, length(t))

    for (i,ti) in enumerate(t)
        out[i] = G_single(ti)
    end

    out
end

t = collect(0.0:0.1:10000)

println("test vectorised version of function")
@btime G_vec($t)
println("\n")

println("test broadcasting using . notation")
@btime G_single.($t)
println("\n")

println("test broadcast function explicitly (should be the same as . notation)")
@btime broadcast($G_single, $t)
println("\n")

println("test using array comprehension")
@btime G_arraycomp($t)
println("\n")

println("test using for loop")
@btime G_loop($t)
println("\n")

println("Benchmarks complete")

max_stars_repo_path: 20190321_alex follow up 1.jl

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_path: 20190321_alex follow up 1.jl

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_path: 20190321_alex follow up 1.jl

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using ITensors
using ITensors.Strided
using LinearAlgebra
using Printf
using Random

Random.seed!(1234)
BLAS.set_num_threads(1)
Strided.set_num_threads(1)
ITensors.enable_threaded_blocksparse()
#ITensors.disable_threaded_blocksparse()

let
  N = 100

  sites = siteinds("S=1",N,conserve_qns=true)

  ampo = AutoMPO()
  for j in 1:N-1
    ampo .+= 0.5, "S+", j, "S-", j+1
    ampo .+= 0.5, "S-", j, "S+", j+1
    ampo .+=      "Sz", j, "Sz", j+1
  end
  H = MPO(ampo, sites)

  # This step makes the MPO more sparse.
  # It generally improves DMRG performance
  # at large bond dimensions but makes DMRG slower at
  # small bond dimensions.
  H = splitblocks(linkinds, H)

  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  psi0 = randomMPS(sites,state,10)

  # Plan to do 5 DMRG sweeps:
  sweeps = Sweeps(5)
  # Set maximum MPS bond dimensions for each sweep
  setmaxdim!(sweeps, 10,20,100,100,200)
  # Set maximum truncation error allowed when adapting bond dimensions
  setcutoff!(sweeps, 1E-10)
  @show sweeps

  # Run the DMRG algorithm, returning energy and optimized MPS
  energy, psi = dmrg(H,psi0, sweeps)
  @printf("Final energy = %.12f\n",energy)
end


text:
using LinearAlgebra

mutable struct ActorCritic{T,TP,TV,OP,OV} <: AbstractPolicyAgent where {T<:Real,TP<:AbstractPolicy,TV<:AbstractFuncApprox,OP<:AbstractOptimizer,OV<:AbstractOptimizer}
    π::TP
    vf::TV
    popt::OP
    vopt::OV
    γ::T
    gp::Array{T}
    gv::Array{T}
    function ActorCritic(π::AbstractPolicy, vf::AbstractFuncApprox, popt::AbstractOptimizer, vopt::AbstractOptimizer, γ::T) where {T<:Real}
        tp = eltype(get_params(π))
        sp = size(get_params(π))
        sv = size(get_params(vf))
        π2 = clone(π)
        vf2 = clone(vf)
        popt2 = clone(popt, π2)
        vopt2 = clone(vopt, vf2)
        new{tp,typeof(π2),typeof(vf2),typeof(popt2),typeof(vopt2)}(π2, vf2, popt2, vopt2, tp(γ), zeros(tp, sp), zeros(tp, sv))
    end
end

function act!(agent::ActorCritic, env::AbstractEnvironment, rng::AbstractRNG)
    v = call_gradient!(agent.gv, agent.vf, env.state, 1)
    logp = get_action_gradient_logp!(agent.gp, agent.π, env.state, rng)

    r = step!(env, agent.π.action, rng)
    target = 0.
    vnext = 0.
    if is_terminal(env)
        target = r
    else
        vnext = call(agent.vf, env.state, 1)
        target = r + agent.γ * vnext
    end

    delta = target - v

    update!(agent.vopt, delta, agent.gv)
    update!(agent.popt, delta, agent.gp)

    if is_terminal(env)
        new_episode!(agent, rng)
    end
end

function act!(agent::ActorCritic{T,TP,TV,OP,Parl2AccumulatingTraceOptimizer{T,TV}}, env::AbstractEnvironment, rng::AbstractRNG) where {T,TP,TV,OP}
    v = call_gradient!(agent.gv, agent.vf, env.state, 1)
    logp = get_action_gradient_logp!(agent.gp, agent.π, env.state, rng)

    r = step!(env, agent.π.action, rng)
    target = 0.
    vnext = 0.
    gvnext = zeros(T, size(agent.gv))
    if is_terminal(env)
        target = r
    else
        vnext = call_gradient!(gvnext, agent.vf, env.state, 1)
        target = r + agent.γ * vnext
    end

    delta = target - v - 0.01 * logp

    update!(agent.vopt, delta, agent.gv, gvnext)
    update!(agent.popt, delta, agent.gp)

    if is_terminal(env)
        new_episode!(agent, rng)
    end
end

function new_episode!(agent::ActorCritic, rng::AbstractRNG)
    new_episode!(agent.vopt)
    new_episode!(agent.popt)
end

function clone(agent::ActorCritic)::ActorCritic
    a = ActorCritic(agent.π, agent.vf, agent.popt, agent.vopt, agent.γ)
    return a
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
module Zarr

import JSON
import Blosc

include("metadata.jl")
include("Storage/Storage.jl")
include("Compressors.jl")
include("ZArray.jl")
include("ZGroup.jl")

export ZArray, ZGroup, zopen, zzeros, zcreate, storagesize, storageratio,
  zinfo, DirectoryStore, S3Store, zgroup

end # module


text:
struct ViolinGeometry <: Gadfly.GeometryElement
    order::Int
    tag::Symbol
end
ViolinGeometry(; order=1, tag=empty_tag) = ViolinGeometry(order, tag)

const violin = ViolinGeometry

element_aesthetics(::ViolinGeometry) = [:x, :y, :color]

default_statistic(::ViolinGeometry) = Gadfly.Stat.violin()

function render(geom::ViolinGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    # TODO: What should we do with the color aesthetic?

    Gadfly.assert_aesthetics_defined("Geom.violin", aes, :y, :width)
    Gadfly.assert_aesthetics_equal_length("Geom.violin", aes, :y, :width)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = fill(theme.default_color, length(aes.y))
    aes = Gadfly.inherit(aes, default_aes)

    # Group y, width and color by x
    ux = unique(aes.x)
    grouped_color = Dict(x => first(aes.color[aes.x.==x]) for x in ux)
    grouped_y = Dict(x => aes.y[aes.x.==x] for x in ux)
    grouped_width = Dict(x => aes.width[aes.x.==x] for x in ux)

    kgy = keys(grouped_y)
    violins = [vcat([(x - w/2, y) for (y, w) in zip(grouped_y[x], grouped_width[x])],
                reverse!([(x + w/2, y) for (y, w) in zip(grouped_y[x], grouped_width[x])]))
                for x in kgy]
    colors = [grouped_color[x] for x in kgy]

    ctx = context(order=geom.order)
    compose!(ctx, Compose.polygon(violins, geom.tag), fill(colors))

    compose!(ctx, svgclass("geometry"))

end


text:
using ReachabilityAnalysis, SparseArrays, JLD2

LazySets.set_ztol(Float64, 1e-14)

const x25 = [zeros(24); 1.0; zeros(23)]
const x25e = vcat(x25, 0.0)

building_path = joinpath(@__DIR__, "building.jld2")

function building_BLDF01()
    @load building_path A B
    n = size(A, 1)
    U = Interval(0.8, 1.0)
    S = @system(x' = Ax + Bu, u ∈ U, x ∈ Universe(n))

    # initial states
    center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
    radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
    X0 = Hyperrectangle(center_X0, radius_X0)

    prob_BLDF01 = InitialValueProblem(S, X0)
end

using ReachabilityAnalysis: add_dimension

function building_BLDC01()
    @load building_path A B
    n = size(A, 1)
    U = Interval(0.8, 1.0)

    # initial states
    center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
    radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
    X0 = Hyperrectangle(center_X0, radius_X0)

    Ae = add_dimension(A)
    Ae[1:n, end] = B
    prob_BLDC01 = @ivp(x' = Ae * x, x(0) ∈ X0 × U)
end


text:
include("system_rhs.jl")
include("kkt_system_solver.jl")
include("schur.jl")
include("schur_direct.jl")
include("symmetric.jl")
include("clever_symmetric.jl")


text:
#=
Taken from https://github.com/JuliaPlots/StatsMakie.jl/blob/master/src/typerecipes/distribution.jl
The StatMakie.jl package is licensed under the MIT "Expat" License:
    Copyright (c) 2018: Pietro Vertechi. =#

# pick a nice default x range given a distribution
function default_range(dist::Distribution, alpha=0.0001)
    minval = isfinite(minimum(dist)) ? minimum(dist) : quantile(dist, alpha)
    maxval = isfinite(maximum(dist)) ? maximum(dist) : quantile(dist, 1 - alpha)
    minval..maxval
end

isdiscrete(::Distribution) = false
isdiscrete(::Distribution{<:VariateForm,<:Discrete}) = true

support(dist::Distribution) = default_range(dist)
support(dist::Distribution{<:VariateForm,<:Discrete}) = UnitRange(endpoints(default_range(dist))...)

convert_arguments(P::PlotFunc, dist::Distribution) = convert_arguments(P, support(dist), dist)

function convert_arguments(P::PlotFunc, x::Union{Interval,AbstractVector}, dist::Distribution)
    default_ptype = isdiscrete(dist) ? ScatterLines : Lines
    ptype = plottype(P, default_ptype)
    to_plotspec(ptype, convert_arguments(ptype, x, x -> pdf(dist, x)))
end
# -----------------------------------------------------------------------------
# qqplots (M. K. Borregaard implementation from StatPlots)

@recipe(QQNorm) do scene
    default_theme(scene, Scatter)
end

@recipe(QQPlot) do scene
    default_theme(scene, Scatter)
end

convert_arguments(::Type{<:QQNorm}, args...; kwargs...) =
    convert_arguments(QQPlot, Distributions.Normal, args...; kwargs...)

convert_arguments(::Type{<:QQPlot}, args...; kwargs...) =
    convert_arguments(Scatter, qqbuild(loc(args...)...); kwargs...)

function convert_arguments(P::PlotFunc, h::QQPair; qqline=:identity)
    line = if qqline in (:fit, :quantile, :identity, :R)
        xs = [extrema(h.qx)...]
        if qqline == :identity
            ys = xs
        elseif qqline == :fit
            itc, slp = hcat(fill!(similar(h.qx), 1), h.qx) \ h.qy
            ys = slp .* xs .+ itc
        else # if qqline == :quantile || qqline == :R
            quantx, quanty = quantile(h.qx, [0.25, 0.75]), quantile(h.qy, [0.25, 0.75])
            slp = diff(quanty) ./ diff(quantx)
            ys = quanty .+ slp .* (xs .- quantx)
        end
        Point{2,Float32}.(xs, ys)
    else
        nothing
    end
    ptype = plottype(Scatter, P)
    PlotSpec{ptype}(h, line)
end

used_attributes(::Type{<:QQNorm}, args...) = (:qqline,)
used_attributes(::Type{<:QQPlot}, args...) = (:qqline,)
used_attributes(::PlotFunc, ::QQPair, args...) = (:qqline,)

function plot!(p::Combined{T,<:Tuple{QQPair,L}}) where {T,L}
    plot!(p, Scatter, attributes_from(Scatter, p), lift(h -> Point{2,Float32}.(h.qx, h.qy), p[1]))
    L !== Nothing && plot!(p, LineSegments, attributes_from(Lines, p), p[2])
end

loc(D::Type{T}, x) where T <: Distribution = Distributions.fit(D, x), x
loc(D, x) = D, x


text:
import PredictMD
import Test

# Test.@test(typeof(PredictMD.registry_url_list()) <: Vector{String})
# Test.@test(length(PredictMD.registry_url_list()) > 0)


text:
include("./ADMM.jl")
# include("./ChainedADMM.jl")


text:
# ---
# title: 836. Rectangle Overlap
# id: problem836
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/rectangle-overlap/description/>
# hidden: true
# ---
#
# An axis-aligned rectangle is represented as a list `[x1, y1, x2, y2]`, where
# `(x1, y1)` is the coordinate of its bottom-left corner, and `(x2, y2)` is the
# coordinate of its top-right corner. Its top and bottom edges are parallel to
# the X-axis, and its left and right edges are parallel to the Y-axis.
#
# Two rectangles overlap if the area of their intersection is **positive**. To
# be clear, two rectangles that only touch at the corner or edges do not
# overlap.
#
# Given two axis-aligned rectangles `rec1` and `rec2`, return `true` _if they
# overlap, otherwise return_`false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: rec1 = [0,0,2,2], rec2 = [1,1,3,3]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: rec1 = [0,0,1,1], rec2 = [1,0,2,1]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: rec1 = [0,0,1,1], rec2 = [2,2,3,3]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `rect1.length == 4`
#   * `rect2.length == 4`
#   * `-109 <= rec1[i], rec2[i] <= 109`
#   * `rec1[0] <= rec1[2]` and `rec1[1] <= rec1[3]`
#   * `rec2[0] <= rec2[2]` and `rec2[1] <= rec2[3]`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using BenchmarkTools
using DitherPunk
using TestImages
using ImageCore

on_CI = haskey(ENV, "GITHUB_ACTIONS")

img_gray = testimage("fabio_gray_256")
img_color = testimage("fabio_color_256")

## Define color scheme
white = RGB{Float32}(1, 1, 1)
yellow = RGB{Float32}(1, 1, 0)
green = RGB{Float32}(0, 0.5, 0)
orange = RGB{Float32}(1, 0.5, 0)
red = RGB{Float32}(1, 0, 0)
blue = RGB{Float32}(0, 0, 1)
cs = [white, yellow, green, orange, red, blue]

# Use one representative algorithm of each type
algs = Dict(
    "error diffusion" => FloydSteinberg(),
    "ordered dithering" => Bayer(; level=3),
    "closest color" => ClosestColor(),
    "threshold dithering" => WhiteNoiseThreshold(),
)

# Tag what is tested on each algorithm type
SUITE = BenchmarkGroup()
SUITE["error diffusion"] = BenchmarkGroup(["binary", "color"])
SUITE["ordered dithering"] = BenchmarkGroup(["binary", "color"])
SUITE["threshold dithering"] = BenchmarkGroup(["binary"])
SUITE["closest color"] = BenchmarkGroup(["binary", "color"])

println(SUITE)

for (algname, alg) in algs
    SUITE[algname]["binary new"] = @benchmarkable dither($(img_gray), $(alg))
    SUITE[algname]["binary inplace"] = @benchmarkable dither!($(copy(img_gray)), $(alg))

    SUITE[algname]["per-channel new"] = @benchmarkable dither($(img_color), $(alg))
    SUITE[algname]["per-channel inplace"] = @benchmarkable dither!(
        $(copy(img_color)), $(alg)
    )

    if "color" in SUITE[algname].tags
        SUITE[algname]["color new"] = @benchmarkable dither($(img_color), $(alg), $(cs))
        SUITE[algname]["color inplace"] = @benchmarkable dither!(
            $(copy(img_color)), $(alg), $(cs)
        )
    end
end


text:
import .ExpressionExplorer: SymbolsState, FunctionName

"Information container about the cells to run in a reactive call and any cells that will err."
Base.@kwdef struct TopologicalOrder
	input_topology::NotebookTopology
	"Cells that form a directed acyclic graph, in topological order."
	runnable::Vector{Cell}
	"Cells that are in a directed cycle, with corresponding `ReactivityError`s."
	errable::Dict{Cell,ReactivityError}
end


text:
using patterns
using Test

tests = ["chain_of_responsibility",
         "composite",
         "decorator",
         "iterator",
         "observer",
         "template",
         "state",
         "singleton"]

@info "Running tests:"

for t in tests
    tfile = string(t, ".jl")
    println()
    println()
    include(tfile)
end


text:
###Measurement data of a Predator-Prey model
### Column 1: population size of prey
### Column 2: population size of predators
function predatorpreydataset1()
    #create the data
    timepoints = linspace(0,50,6)
    measurements =
        [51.7002527839933	3.62886163423551;
         88.4059566464213	11.8592208786134;
         50.5971643551162	33.1299499686354;
         25.2455746232269	30.1443180597855;
         42.2660932030265	9.00691042344652;
         51.4088102060744	29.8816967299902]

    #return the data
    timepoints,measurements
end


text:
################################################################################
#
#  NfAbs/NonSimple.jl : non-simple absolute number fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017, 2018: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2018 Tommy Hofmann, Claus Fieker
#
################################################################################

export NfAbsNS, NfAbsNSElem

export issimple, simple_extension

@inline base_ring(K::NfAbsNS) = FlintQQ

@inline degree(K::NfAbsNS) = K.degree

@inline degrees(K::NfAbsNS) = K.degrees

@inline ngens(K::NfAbsNS) = length(K.pol)

################################################################################
#
#  Copy
#
################################################################################

function Base.deepcopy_internal(a::NfAbsNSElem, dict::ObjectIdDict)
  # TODO: Fix this once deepcopy is fixed for fmpq_mpoly
  # z = NfAbsNSElem(Base.deepcopy_internal(data(a), dict))
  z = NfAbsNSElem(Base.deepcopy(data(a)))
  z.parent = parent(a)
  return z
end

#julia's a^i needs copy
function Base.copy(a::NfAbsNSElem)
  return parent(a)(a.data)
end

################################################################################
#
#  Comply with Nemo ring interface
#
################################################################################

elem_type(::Type{NfAbsNS}) = NfAbsNSElem

elem_type(::NfAbsNS) = NfAbsNSElem

parent_type(::Type{NfAbsNSElem}) = NfAbsNS

needs_parentheses(::NfAbsNSElem) = true

isnegative(x::NfAbsNSElem) = Nemo.isnegative(data(x))

show_minus_one(::Type{NfAbsNSElem}) = true

function iszero(a::NfAbsNSElem)
  reduce!(a)
  return iszero(data(a))
end

function isone(a::NfAbsNSElem)
  reduce!(a)
  return isone(data(a))
end

Nemo.zero(K::NfAbsNS) = K(Nemo.zero(parent(K.pol[1])))

Nemo.one(K::NfAbsNS) = K(Nemo.one(parent(K.pol[1])))

Nemo.one(a::NfAbsNSElem) = one(a.parent)

dot(a::NfAbsNSElem, b::Union{Integer, fmpz}) = a * b

dot(a::Union{Integer, fmpz}, b::NfAbsNSElem) = b * a

################################################################################
#
#  Random
#
################################################################################

function rand(K::NfAbsNS, r::UnitRange)
  # TODO: This is super slow
  b = basis(K)
  z = K()
  for i in 1:degree(K)
    z += rand(r) * b[i]
  end
  return z
end

################################################################################
#
#  Promotion
#
################################################################################

if isdefined(Nemo, :promote_rule1)
  Nemo.promote_rule{T <: Integer}(::Type{NfAbsNSElem}, ::Type{T}) = NfAbsNSElem

  Nemo.promote_rule(::Type{NfAbsNSElem}, ::Type{fmpz}) = NfAbsNSElem

  Nemo.promote_rule(::Type{NfAbsNSElem}, ::Type{fmpq}) = NfAbsNSElem

  function Nemo.promote_rule(::Type{NfAbsNSElem}, ::Type{U}) where {U}
    Nemo.promote_rule(fmpq, U) == fmpq ? NfAbsNSElem : Nemo.promote_rule1(U, NfAbsNSElem)
  end
else
  Nemo.promote_rule{T <: Integer}(::Type{NfAbsNSElem}, ::Type{T}) = NfAbsNSElem

  Nemo.promote_rule(::Type{NfAbsNSElem}, ::Type{fmpz}) = NfAbsNSElem

  Nemo.promote_rule(::Type{NfAbsNSElem}, ::Type{fmpq}) = NfAbsNSElem

  function Nemo.promote_rule(::Type{NfAbsNSElem}, ::Type{U}) where {U <: Nemo.RingElement}
    Nemo.promote_rule(fmpq, U) == fmpq ? NfAbsNSElem : Union{}
  end
end

################################################################################
#
#  Field access
#
################################################################################

@inline Nemo.data(a::NfAbsNSElem) = a.data

@inline Nemo.parent(a::NfAbsNSElem) = a.parent::NfAbsNS

issimple(a::NfAbsNS) = false

issimple(::Type{NfAbsNS}) = false

function basis(K::NfAbsNS)
  if isdefined(K, :basis)
    return copy(K.basis)::Vector{NfAbsNSElem}
  else
    g = gens(K)
    b = NfAbsNSElem[]
    for i in CartesianRange(Tuple(1:degrees(K)[i] for i in 1:ngens(K)))
      push!(b, prod(g[j]^(i[j] - 1) for j=1:length(i)))
    end
    K.basis = b
    return copy(b)::Vector{NfAbsNSElem}
  end
end

# Given an exponent vector b, the following function returns the index of
# the basis element corresponding to b.
function monomial_to_index(K::NfAbsNS, b::Vector{T}) where {T}
  n = ngens(K)
  idx = b[n]
  for j in n-1:-1:1
    idx *= degrees(K)[j]
    idx += b[j]
  end
  return Int(idx + 1)
end

################################################################################
#
#  Reduction
#
################################################################################

function reduce!(a::NfAbsNSElem)
  q, a.data = divrem(a.data, parent(a).pol)
  return a
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, a::NfAbsNS)
  print(io, "Non-simple number field with defining polynomials ", a.pol)
end

#TODO: this is a terrible show func.
function Base.show(io::IO, a::NfAbsNSElem)
  f = data(a)
  show(io, f)
end

################################################################################
#
#  Unary operators
#
################################################################################

function Base.:(-)(a::NfAbsNSElem)
  return parent(a)(-data(a))
end

################################################################################
#
#  Binary operators
#
################################################################################

function Base.:(+)(a::NfAbsNSElem, b::NfAbsNSElem)
  return parent(a)(data(a) + data(b))
end

function Base.:(-)(a::NfAbsNSElem, b::NfAbsNSElem)
  return parent(a)(data(a) - data(b))
end

function Base.:(*)(a::NfAbsNSElem, b::NfAbsNSElem)
  return parent(a)(data(a) * data(b))
end

function Base.:(//)(a::NfAbsNSElem, b::NfAbsNSElem)
  return div(a, b)
end

function Nemo.div(a::NfAbsNSElem, b::NfAbsNSElem)
  return a * inv(b)
end

Nemo.divexact(a::NfAbsNSElem, b::NfAbsNSElem) = div(a, b)

################################################################################
#
#  Powering
#
################################################################################

function Base.:(^)(a::NfAbsNSElem, b::Integer)
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  else#if mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

function Base.:(^)(a::NfAbsNSElem, b::fmpz)
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  else# mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

################################################################################
#
#  Comparison
#
################################################################################

function Base.:(==)(a::NfAbsNSElem, b::NfAbsNSElem)
  reduce!(a)
  reduce!(b)
  return data(a) == data(b)
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function Nemo.mul!(c::NfAbsNSElem, a::NfAbsNSElem, b::NfAbsNSElem)
  mul!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

function Nemo.add!(c::NfAbsNSElem, a::NfAbsNSElem, b::NfAbsNSElem)
  add!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

function Nemo.addeq!(b::NfAbsNSElem, a::NfAbsNSElem)
  addeq!(b.data, a.data)
  b = reduce!(b)
  return b
end

#

function Nemo.mul!(c::NfAbsNSElem, a::NfAbsNSElem, b::fmpz)
  mul!(c.data, a.data, parent(c.data)(b))
  c = reduce!(c)
  return c
end

################################################################################
#
#  Conversion to matrix
#
################################################################################

function elem_to_mat_row!(M::fmpz_mat, i::Int, d::fmpz, a::NfAbsNSElem)
  K = parent(a)
  # TODO: This is super bad
  # Proper implementation needs access to the content of the underlying
  # fmpq_mpoly

  for j in 1:cols(M)
    M[i, j] = zero(FlintZZ)
  end

  one!(d)

  if length(data(a)) == 0
    return nothing
  end

  z = zero_matrix(FlintQQ, 1, cols(M))
  elem_to_mat_row!(z, 1, a)
  z_q = FakeFmpqMat(z)

  for j in 1:cols(M)
    M[i, j] = z_q.num[1, j]
  end

  ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), d, z_q.den)

  return nothing
end

function elem_to_mat_row!(M::fmpq_mat, i::Int, a::NfAbsNSElem)
  K = parent(a)
  for j in 1:cols(M)
    M[i, j] = zero(FlintQQ)
  end
  adata = data(a)
  for j in 1:length(adata)
    exps = Nemo._get_termexp_ui(adata, j)
    k = monomial_to_index(K, exps)
    M[i, k] = coeff(adata, j)
  end
  return M
end

function elem_from_mat_row(K::NfAbsNS, M::fmpq_mat, i::Int)
  a = K()
  b = basis(K)
  for c = 1:cols(M)
    a += M[i, c]*b[c]
  end
  return a
end

function elem_from_mat_row(K::NfAbsNS, M::fmpz_mat, i::Int, d::fmpz)
  a = K()
  b = basis(K)
  for c = 1:cols(M)
    a += M[i, c]*b[c]
  end
  return divexact(a, d)
end

function SRow(a::NfAbsNSElem)
  sr = SRow(FlintQQ)
  adata = data(a)
  for i=1:length(adata)
    # TODO: Do this inplace with preallocated exps array
    exps = Nemo._get_termexp_ui(adata, i)
    push!(sr.pos, monomial_to_index(parent(a), exps))
    push!(sr.values, coeff(adata, i))
  end
  p = sortperm(sr.pos)
  sr.pos = sr.pos[p]
  sr.values = sr.values[p]
  return sr
end

################################################################################
#
#  Minimal polynomial
#
################################################################################

function minpoly_dense(a::NfAbsNSElem)
  K = parent(a)
  n = degree(K)
  M = zero_matrix(FlintQQ, degree(K)+1, degree(K))
  z = a^0
  elem_to_mat_row!(M, 1, z)
  z *= a
  elem_to_mat_row!(M, 2, z)
  i = 2
  Qt, _ = PolynomialRing(FlintQQ,"t", cached=false)
  while true
    if n % (i-1) == 0 && rank(M) < i
      N = nullspace(sub(M, 1:i, 1:cols(M))')
      @assert N[1] == 1
      v = Vector{fmpq}(i)
      for j in 1:i
        v[j] = N[2][j, 1]
      end
      #f = Qt([N[2][j, 1] for j=1:i])
      f = Qt(v)
      return f*inv(lead(f))
    end
    z *= a
    elem_to_mat_row!(M, i+1, z)
    i += 1
  end
end

function minpoly_sparse(a::NfAbsNSElem)
  K = parent(a)
  n = degree(K)
  M = SMat(FlintQQ)
  z = a^0
  sz = SRow(z)
  i = 0
  push!(sz.values, FlintQQ(1))
  push!(sz.pos, n+i+1)
  push!(M, sz)
  z *= a
  sz = SRow(z)
  i = 1
  Qt, t = PolynomialRing(FlintQQ, "t", cached = false)
  while true
    if n % i == 0
      echelon!(M)
      fl, so = cansolve_ut(sub(M, 1:i, 1:n), sz)
      if fl
        so = mul(so, sub(M, 1:i, n+1:cols(M)))
        # TH: If so is the zero vector, we cannot use the iteration,
        # so we do it by hand.
        if length(so.pos) == 0
          f = t^i
        else
          f = t^i - sum(c*t^(k-1) for (k,c) = so)
        end
        return f
      end
    end
    push!(sz.values, FlintQQ(1))
    push!(sz.pos, n+i+1)
    push!(M, sz)
    z *= a
    sz = SRow(z)
    i += 1
    if i > degree(parent(a))
      error("too large")
    end
  end
end

function minpoly(a::NfAbsNSElem)
  return minpoly_sparse(a)
end

################################################################################
#
#  Characteristic polynomial
#
################################################################################

function charpoly(a::NfAbsNSElem)
  f = minpoly(a)
  return f^div(degree(parent(a)), degree(f))
end

################################################################################
#
#  Inverse
#
################################################################################

function inv(a::NfAbsNSElem)
  if iszero(a)
    error("division by zero")
  end
  f = minpoly(a)
  z = parent(a)(coeff(f, degree(f)))
  for i=degree(f)-1:-1:1
    z = z*a + coeff(f, i)
  end
  return -z*inv(coeff(f, 0))
end

################################################################################
#
#  Norm
#
################################################################################

function norm(a::NfAbsNSElem)
  f = minpoly(a)
  return (-1)^degree(parent(a)) * coeff(f, 0)^div(degree(parent(a)), degree(f))
end

################################################################################
#
#  Trace
#
################################################################################

function trace(a::NfAbsNSElem)
  f = minpoly(a)
  return -coeff(f, degree(f)-1)*div(degree(parent(a)), degree(f))
end

################################################################################
#
#  Representation matrix
#
################################################################################

function representation_matrix(a::NfAbsNSElem)
  K = parent(a)
  b = basis(K)
  M = zero_matrix(FlintQQ, degree(K), degree(K))
  for i=1:degree(K)
    elem_to_mat_row!(M, i, a*b[i])
  end
  return M
end

function representation_matrix_q(a::NfAbsNSElem)
  M = representation_matrix(a)
  # TODO: This is suboptimal.
  Mf = FakeFmpqMat(M)
  return Mf.num, Mf.den
end

################################################################################
#
#  Substitution
#
################################################################################

# TODO: Preallocate the exps array
function msubst(f::fmpq_mpoly, v::Array{T, 1}) where {T}
  n = length(v)
  @assert n == nvars(parent(f))
  r = FlintQQ()
  for i=1:length(f)
    exps = Nemo._get_termexp_ui(f, i)
    r += coeff(f, i) * prod(v[j]^exps[j] for j=1:n)
  end
  return r
end

################################################################################
#
#  Morphisms
#
################################################################################

mutable struct NfAbsToNfAbsNS <: Map{AnticNumberField, NfAbsNS, HeckeMap, NfAbsToNfAbsNS}
  header::MapHeader{AnticNumberField, NfAbsNS}
  prim_img::NfAbsNSElem
  emb::Array{nf_elem, 1}

  function NfAbsToNfAbsNS(K::AnticNumberField, L::NfAbsNS, a::NfAbsNSElem, emb::Array{nf_elem, 1})
    function image(x::nf_elem)
      # x is an element of K
      f = x.parent.pol.parent(x)
      return f(a)
    end

    function preimage(x::NfAbsNSElem)
      return msubst(data(x), emb)
    end

    z = new()
    z.prim_img = a
    z.emb = emb
    z.header = MapHeader(K, L, image, preimage)
    return z
  end

  function NfAbsToNfAbsNS(K::AnticNumberField, L::NfAbsNS, a::NfAbsNSElem)
    function image(x::nf_elem)
      # x is an element of K
      f = x.parent.pol.parent(x)
      return f(a)
    end

    z = new()
    z.prim_img = a
    z.header = MapHeader(K, L, image)
    return z
  end
end

mutable struct NfAbsNSToNfAbsNS <: Map{NfAbsNS, NfAbsNS, HeckeMap, NfAbsNSToNfAbsNS}
  header::MapHeader{NfAbsNS, NfAbsNS}
  emb::Array{NfAbsNSElem, 1}

  function NfAbsNSToNfAbsNS(K::NfAbsNS, L::NfAbsNS, emb::Array{NfAbsNSElem, 1})
    function image(x::NfAbsNSElem)
      # x is an element of K
      return msubst(data(x), emb)
    end

    z = new()
    z.emb = emb
    z.header = MapHeader(K, L, image)
    return z
  end
end

# TODO: The following is opposite to our new convention
function Base.:(*)(f::NfAbsNSToNfAbsNS, g::NfAbsNSToNfAbsNS)
  domain(f) == codomain(g) || throw("Maps not compatible")
  a = gens(domain(g))
  return NfAbsNSToNfAbsNS(domain(g), codomain(f), [ f(g(x)) for x in a])
end

function Base.:(==)(f::NfAbsNSToNfAbsNS, g::NfAbsNSToNfAbsNS)
  if domain(f) != domain(g) || codomain(f) != codomain(g)
    return false
  end

  L = domain(f)

  for a in gens(L)
    if f(a) != g(a)
      return false
    end
  end

  return true
end

################################################################################
#
#  Simple extensions
#
################################################################################

function simple_extension(K::NfAbsNS)
  n = ngens(K)
  g = gens(K)

  pe = g[1]
  i = 1
  ind = Int[1]
  f = minpoly(pe)
  #TODO: use resultants rather than minpoly??
  while i < n
    i += 1
    j = 1
    f = minpoly(pe + j * g[i])
    while degree(f) < prod(total_degree(K.pol[k]) for k in 1:i)
      j += 1
      f = minpoly(pe + j * g[i])
    end
    push!(ind, j)
    pe += j * g[i]
  end
  Ka, a = number_field(f)
  k = base_ring(K)
  M = zero_matrix(k, degree(K), degree(K))
  z = one(K)
  elem_to_mat_row!(M, 1, z)
  elem_to_mat_row!(M, 2, pe)
  z *= pe
  for i=3:degree(K)
    z *= pe
    elem_to_mat_row!(M, i, z)
  end
  N = zero_matrix(k, 1, degree(K))
  b = basis(Ka)
  emb = typeof(b)()
  for i=1:n
    elem_to_mat_row!(N, 1, g[i])
    s = solve(M', N')
    push!(emb, sum(b[j] * s[j,1] for j in 1:degree(Ka)))
  end

  return Ka, NfAbsToNfAbsNS(Ka, K, pe, emb)
end

################################################################################
#
#  Constructors and parent object overloading
#
################################################################################

function number_field(f::Array{fmpq_poly, 1}, s::String="_\$")
  S = Symbol(s)
  n = length(f)
  Qx, x = PolynomialRing(FlintQQ, n, s)
  K = NfAbsNS(fmpq_mpoly[f[i](x[i]) for i=1:n],
              Symbol[Symbol("$s$i") for i=1:n])
  K.degrees = [degree(f[i]) for i in 1:n]
  K.degree = prod(K.degrees)
  return K, gens(K)
end

gens(K::NfAbsNS) = [K(x) for x = gens(parent(K.pol[1]))]

function (K::NfAbsNS)(a::fmpq_mpoly)
  q, w = divrem(a, K.pol)
  z = NfAbsNSElem(w)
  z.parent = K
  return z
end

(K::NfAbsNS)(a::Integer) = K(parent(K.pol[1])(a))

(K::NfAbsNS)(a::Rational{T}) where {T <: Integer} = K(parent(K.pol[1])(a))

(K::NfAbsNS)(a::fmpz) = K(parent(K.pol[1])(a))

(K::NfAbsNS)(a::fmpq) = K(parent(K.pol[1])(a))

(K::NfAbsNS)() = zero(K)


text:
module Volcanito

export
    @select,
    @where,
    @group_by,
    @aggregate_vector,
    @order_by,
    @limit,
    @inner_join,
    @left_join,
    @right_join,
    @outer_join,
    materialize

import DataFrames
import MacroTools: inexpr, postwalk
import Printf: @printf

include("query/expression_operations/escape_symbols.jl")
include("query/expression_operations/column_name.jl")
include("query/expression_operations/remove_backticks.jl")
include("query/expression_operations/validate.jl")
include("query/expression_operations/get_alias.jl")
include("query/expression_operations/find_column_names.jl")
include("query/expression_operations/index_column_names.jl")
include("query/expression_operations/gensym_index.jl")
include("query/expression_operations/rewrite_column_names.jl")
include("query/expression_operations/safe_tuple_name.jl")
include("query/expression_operations/predicates.jl")
include("query/expression_operations/fusion.jl")
include("query/expression_operations/passes/locals.jl")
include("query/expression_operations/passes/lift.jl")
include("query/expression_operations/passes/tvl.jl")
include("query/expression_operations/forms/tuple_form.jl")
include("query/expression_operations/forms/broadcast_form.jl")
include("query/expression_operations/forms/vector_form.jl")
include("query/expression_operations/expression/column_names_tuple_expr.jl")
include("query/expression_operations/expression/expression.jl")
include("query/expression_operations/join_utils.jl")

include("query/macros/select.jl")
include("query/macros/where.jl")
include("query/macros/group_by.jl")
include("query/macros/aggregate_vector.jl")
include("query/macros/order_by.jl")
include("query/macros/limit.jl")
include("query/macros/inner_join.jl")
include("query/macros/left_join.jl")
include("query/macros/right_join.jl")
include("query/macros/outer_join.jl")

include("logical_nodes/logical_node.jl")
include("logical_nodes/projection.jl")
include("logical_nodes/selection.jl")
include("logical_nodes/group_by.jl")
include("logical_nodes/aggregate_vector.jl")
include("logical_nodes/order_by.jl")
include("logical_nodes/limit.jl")
include("logical_nodes/join.jl")

include("physical_operations/dataframes/iterators.jl")
include("physical_operations/dataframes/materialize.jl")
include("physical_operations/dataframes/generate_column.jl")
include("physical_operations/dataframes/projection.jl")
include("physical_operations/dataframes/selection.jl")
include("physical_operations/dataframes/group_by.jl")
include("physical_operations/dataframes/aggregate_vector.jl")
include("physical_operations/dataframes/order_by.jl")
include("physical_operations/dataframes/limit.jl")
include("physical_operations/dataframes/join.jl")

end


text:
"""
    WSymbol

A Wolfram language symbol. The `W""` string macro can be used as a short form.
"""
struct WSymbol
    name::String
end
WSymbol(sym::Symbol) = WSymbol(String(sym))
Base.print(io::IO, s::WSymbol) = print(io, s.name)
Base.show(io::IO, s::WSymbol) = print(io, 'W', '"', s.name, '"')

macro W_str(str)
    WSymbol(str)
end

struct WReal
    value::String
end
Base.show(io::IO, x::WReal) = print(io, x.value)

struct WInteger
    value::String
end
Base.show(io::IO, x::WInteger) = print(io, x.value)

"""
    WExpr

A Wolfram language expression.
"""
struct WExpr
    head
    args
end

function Base.print(io::IO, w::WExpr)
    print(io, w.head)
    print(io, '[')
    join(io, w.args, ", ")
    print(io, ']')
end
function Base.show(io::IO, w::WExpr)
    show(io, w.head)
    print(io, '(')
    isfirst = true
    for arg in w.args
        if !isfirst
            print(io, ", ")
        else
            isfirst = false
        end
        show(io, arg)
    end
    print(io, ')')
end

(w::WSymbol)(args...) = WExpr(w, args)
(w::WExpr)(args...) = WExpr(w, args)


text:
"""
This is a mess at the moment. :)
"""

include("viaCheb.jl")
include("classical.jl")
using .bycheb
using LinearAlgebra, Plots
using .Classical
using Plots

method="Explicit"
#method="Implicit"
#method="ImplicitPinv"
#method="Mixed"
μ = 0.5

print("Method: $method\n")

normby = "sum" #true

FACTOR = 0.1

α= 1.0    * FACTOR
ν = 1.0/3   * FACTOR * 9
s = 1.0   * FACTOR * 0.3
Lrel = 5

#N=40
N = floor(Int, 73) #64 / √2) # floor(Int, 73 * √2)

M = 500.0
nMx=0

print("M: $M\n")

TINY = 0.1#1.0#0.1 #0.1

β = 1.0
print("β: $β\n")

#################
rc = Classical.ring(32;xf=4.0)
fig = Plots.plot()

#####

function eye(n::Integer)
    diagm(ones(n))
end

nn = N÷2 + N%2
r = bycheb.ring(N,(α,ν,s,Lrel))
Lcrit = r.Ls["Lcrit"]
Lvisc = r.Ls["Lvisc"]
#Plots.vline!([L0-Lcrit],label="crit",legend=:bottomleft)
Plots.vline!([Lcrit/Lvisc],label="crit",legend=:bottom)#:topleft)
#Plots.vline!([L0-Lvisc],label="visc")
Plots.vline!([Lvisc/Lvisc],label="visc")
#display(fig)
ws = sqrt(3*ν/s) # wave speed
v=[ones(nn);zeros(nn)]
#v[1:nn÷2].=0.0
#Σ = view(v,    1:nn )
#Φ = view(v, nn+1:e)
minAd1 = minimum( abs.(diff(r.y)./(abs.(r.u[1:end-1]).+ws)))
minAd2 = minimum( abs.(diff(r.y)./(abs.(r.u[2:end]).+ws)))
minDiff = minimum( diff(r.y).^2 / (3*ν) )
dt = TINY * 0.5 * min(minAd1, minAd2, minDiff)
print("dt: $dt \n")

IM = eye(2*nn) - dt*(-r.AA+r.B)
MM = eye(2*nn) - μ*dt*(-r.AA+r.B)
#PinvIM = pinv(IM)
Nm = r.Ls["L0"]^5/(5α) / dt# advective time is L0^5/(5α)/dt
NNN = M * Nm
Np = 1002003 #trunc(Int,Nm*0.001)#0

print("Nm: $Nm   NNN: $NNN \n")

rc.β = β
for j in 1:1000000
    Classical.iterate!(rc)
end
Plots.plot!(rc.x_[1:24],rc.σ_[1:24];legend=false)
display(fig)

function normΣ(w,nn,flag)
    if flag=="L1 norm"
        w./=(sum(w[1:nn]) / nn) #not quite right but close?
    elseif flag=="maxnorm"
        w./=maximum(w[1:nn])
    else
        @warn "Error in normΣ.\n"
    end
end
function myplot(r,v;nM="final",first=true)
    ynn=r.y[1:nn]
    Σ = v[1:nn]
    Φ = v[nn+1:end]
    L0 = r.Ls["L0"]
    local Lcrit = r.Ls["Lcrit"]
    local Lvisc = r.Ls["Lvisc"]
    NX = r.parmd["X"] # key dimensionless number

    XXXv = (L0.-ynn)./Lvisc
    XXXc = (L0.-ynn)./Lcrit
    if first
        f = Plots.plot
    else
        f = Plots.plot!
    end
    f(XXXv,Σ,title="X: $NX  α: $α  ν: $ν\ns: $s  N: $N  M: $M  Lrel: $Lrel",label=:none,xlims=(0.0,(L0/Lvisc)))

    #Plots.vline!([L0],label="moon")
end
function myplot!(r,v;nM="final")
    myplot(r,v;nM=nM,first=false)
end

for i in 1:trunc(Int,1000*Np)#NNN)
    global v
    local w
    if method=="Explicit"
        w = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; # w./=sum(w[1:nn]);
        if normby=="max"
            w./=maximum(w[1:nn])
        elseif normby=="sum"
            w./=sum(w[1:nn]/nn);
        elseif normby=="bdy"
            if w[1] <= 0
                w./=sum(w[1:nn]/nn)/nn
            else
                w./= w[1]
                #print("Got here!")
            end
        else
            @warn "Failure\n"
        end#if
        bycheb.setκ!(r,w[1:nn],β)
        v = w + dt*(-r.AA+r.B)*w
    elseif method=="Implicit"
        w = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; # w[1]=0.0; w[nn+1]=0.0;
        if normby=="max"
            w./=maximum(w[1:nn])
        elseif normby=="sum"
            w./=sum(w[1:nn]/nn);
        elseif normby=="bdy"
            if w[1] <= 0
                w./=sum(w[1:nn]/nn)/nn
            else
                w./= w[1]
                #print("Got here!")
            end
        else
            @warn "Failure\n"
        end#if
        bycheb.setκ!(r,w[1:nn],β)
        local IM = eye(2*nn) - dt*(-r.AA+r.B)
        v = IM\w
    elseif method=="ImplicitPinv"
        w = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; w./=sum(w[1:nn]);
        if normbymax
            w./=maximum(w[1:nn])
        else
            w./=sum(w[1:nn]);
        end#if
        bycheb.setκ!(r,w[1:nn],β)
        IM = eye(2*nn) - dt*(-r.AA+r.B)
        local PinvIM = pinv(IM)
        for k in 1:100
            w = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; w./=sum(w[1:nn]);
            if normbymax
                w./=maximum(w[1:nn])
            else
                w./=sum(w[1:nn]);
            end#if
            v = PinvIM * w
        end
    elseif method=="Mixed"
        w  = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; w./=sum(w[1:nn]);
        vv = w + μ*dt*(-r.AA+r.B)*w
        v = MM\vv
    else
        print("WARNING: failure.\n")
    end
    if i%Np==0
        global nMx+=1
        myplot!(r,v;nM="$nMx")
        display(fig)
    end
end
myplot!(r,v)
display(fig)

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
module Nebula

    include("Clock.jl")
    include("SDLwrapper.jl")

end #module Nebula


text:
################
# AbstNSCode
################

abstract type AbstNSCode <: Function end

################
# NSCodedict0
################

const _NSCodedict0 = Dict{Symbol, Function}()

_NSCodedict0[:cls] = nsc -> nsc.__cls
_NSCodedict0[:init] = nsc -> nsc.__init[1]

_NSCodedict0[:cst] = nsc -> NSCodecst(nsc)

_NSCodedict0[:dfn] = nsc -> NSCodeTag{:dfn, false}(nsc)
_NSCodedict0[:req] = nsc -> NSCodeTag{:req, false}(nsc)
_NSCodedict0[:prp] = nsc -> NSCodeTag{:prp, false}(nsc)
_NSCodedict0[:fnc] = nsc -> NSCodeTag{:fnc, false}(nsc)
_NSCodedict0[:mth] = nsc -> NSCodeTag{:mth, false}(nsc)

_NSCodedict0[:_instances] = nsc -> [i for (a, k, i) ∈ nsc.__instances]
_NSCodedict0[:_clr_instances] =
    nsc -> (deleteat!(nsc.__instances, 1:length(nsc.__instances)); return)


text:
@testset "anomalies.jl" begin

    data_surv = CSO.dataset("cso_raw")
    data_trudeau = CSO.dataset("Montreal-Trudeau_15min_precipitations")
    data = innerjoin(data_surv, data_trudeau, on=:Date)

    dropmissing!(data, :Duration)
    addColOverflow!(data, "Duration", "Surverse", 3)

    transformedDF = deleteAnomalies!(data, "Surverse", "d24h", 2.0)
    testDF = filter([:Surverse, :d24h] => (x, y) -> (x == 1 && y <= 2), transformedDF)

    #There is no day with overflow and without precipitation in the DF
    @test nrow(testDF) == 0

    #Ancient DF has been changed
    @test nrow(data) == nrow(transformedDF)

end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
using Makie
using Distributions

const d = Normal(0, 1)

𝓝(x) = pdf(d, x)
prior(x, a, b) = (1 - a) * 𝓝(x) + a * 𝓝(x - b)
function main(a, b)
    xs = -5:0.01:5
    ys = [prior(x,a,b) for x in xs]
    lines(xs, ys)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots, StaticArrays, Parameters
default(colorbar=false)

rotate(u, v, w) = let c = cosd, s = sind
    SMatrix{3, 3}(
        c(u)*c(v), c(u)*s(v)*s(w)-s(u)*c(w), c(u)*s(v)*c(w)+s(u)*s(w),
        s(u)*c(v), s(u)*s(v)*s(w)+c(u)*c(w), c(u)*s(v)*s(w)-c(u)*s(w),
            -s(v),      c(v)*s(w),                c(v)*c(w))'
end

"""Assume `surffunc(u, v)` returns `(x, y, z, h)`."""
function plot_surface(u, v, surffunc, param; R = rotate(0, 0, 0), kwargs...)
    xyzh = surffunc.(u', v, Ref(param))
    xyz = (((x, y, z, h),) -> R * SVector(x, y, z)).(xyzh)
    x, y, z = ((a -> a[i]).(xyz) for i in 1:3)
    h = (((x, y, z, h),) -> h).(xyzh)
    surface(x, y, z; fill_z = h, kwargs...)
end

# %%
pyplot(fmt=:png)

function torusfunc(u, v, param)
    @unpack a, b, d = param
    x = (b + a * cos(u)) * cos(v)
    y = (b + a * cos(u)) * sin(v)
    z =      a * sin(u)
    h = d' * SVector(x, y, z)
    x, y, z, h
end

torusparam = (a = 5, b = 10, d = SVector(1, 0, 1))
n = 50
u_torus = v_torus = range(0, 2π; length = n + 1)
plot_surface(u_torus, v_torus, torusfunc, torusparam; size=(500, 400),
    lims=(-20, 20), color=:CMRmap, R=rotate(30, 30, 30))

# %%
pyplot(fmt=:png)
@time anim = @animate for t in range(0, 360; length=41)[1:end-1]
    plot_surface(u_torus, v_torus, torusfunc, torusparam; size=(400, 300),
        lims=(-20, 20), color=:CMRmap, R=rotate(30, 30, 30),
        camera=(t, 20), ticks=false)
end
PyPlot.clf()
gif(anim, "torus.gif")

# %%
pyplot(fmt=:png)

function spherefunc(u, v, param)
    @unpack r, d, f = param
    x = r * cos(u) * cos(v)
    y = r * cos(u) * sin(v)
    z = r * sin(u)
    h = f(d' * SVector(x, y, z))
    x, y, z, h
end

sphereparam = (
    r = 20,
    d = SVector(0, 0, 1),
    f = h -> (1 + abs(h - 1)) * sin(h)
)
n = 50
u_sphere = range(-π/2, π/2; length = n + 1)
v_sphere = range(0, 2π; length = 2n + 1)
plot_surface(u_sphere, v_sphere, spherefunc, sphereparam; size=(500, 400),
    lims=(-20, 20), color=:gist_earth, R=rotate(30, 30, 30),
    camera=(30, 20))

# %%
pyplot(fmt=:png)
@time anim = @animate for t in range(0, 360; length=41)[1:end-1]
    plot_surface(u_sphere, v_sphere, spherefunc, sphereparam; size=(400, 300),
        lims=(-20, 20), color=:gist_earth, R=rotate(30, 30, 30),
        camera=(t, 20), ticks=false)
end
PyPlot.clf()
gif(anim, "sphere.gif")

# %%
plotly(fmt=:auto)
plot_surface(u_torus, v_torus, torusfunc, torusparam; size=(500, 500),
    lims=(-20, 20), color=:CMRmap, R=rotate(30, 30, 30))

# %%
plotly(fmt=:auto)
plot_surface(u_sphere, v_sphere, spherefunc, sphereparam; size=(500, 500),
    lims=(-20, 20), color=:gist_earth, R=rotate(30, 30, 30))

# %%
using Pkg
println("Julia v", VERSION)
Pkg.status("IJulia")
Pkg.status("PyCall"; mode = PKGMODE_MANIFEST)
Pkg.status("Plots")
Pkg.status("PyPlot")
Pkg.status("Plotly")
Pkg.status("StaticArrays")
Pkg.status("Parameters")

# %%

max_stars_repo_path: 0016/3d colored surface plots.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/3d colored surface plots.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/3d colored surface plots.jl

max_forks_repo_name: genkuroki/public


text:
using NetworkLayout, Roots
export MultigraphPlot, relabel, get_curves

function relabel(edgelist::Vector{Tuple{Int,Int}})
    edgelist_relab = Tuple{Int,Int}[]
    vmap = Dict{Int,Int}()
    i = 1
    for edge in edgelist
        for v in edge
            if v ∉ keys(vmap)
                vmap[v] = i
                i+=1
            end
        end
        push!(edgelist_relab, map(x->vmap[x], edge))
    end
    rev_vmap = Dict(v=>k for (k,v) in vmap)
    return edgelist_relab, rev_vmap
end

function relabel(edgelist::Vector{Tuple{String,String}})
    edgelist_relab = Tuple{Int,Int}[]
    vmap = Dict{String,Int}()
    i = 1
    for edge in edgelist
        for v in edge
            if v ∉ keys(vmap)
                vmap[v] = i
                i+=1
            end
        end
        push!(edgelist_relab, map(x->vmap[x], edge))
    end
    rev_vmap = Dict(v=>k for (k,v) in vmap)
    return edgelist_relab, rev_vmap
end

node_edge_intercept(x,w,h,r) = (4*h/(w^2) * x*(w-x))^2 + x^2 - r^2

function get_arc(
    src, dst;
    bendprop=0.1
    )
    # f(x) = 4h/w² x(w-x) (equation for arc from (0,0)->(w,0))
    # We then then translate and rotate this
    w = euclidean(src,dst)
    h = bendprop*euclidean(src,dst)
    θ = angle(complex(dst...)-complex(src...))
    cosθ, sinθ = (cos(θ), sin(θ))
    C = 4*h/w^2
    rawx = 0:0.01:w
    srcx,srcy = src
    xvals = [x*cosθ - C*x*(w-x)*sinθ + srcx for x in rawx]
    yvals = [x*sinθ + C*x*(w-x)*cosθ + srcy for x in rawx]
    return xvals, yvals
end

function get_arc_shorten(
    src, dst;
    bendprop=0.1,
    shorten=0.2
    )

    w = euclidean(src,dst)
    h = bendprop*euclidean(src,dst)
    x_shorten = find_zero(x->node_edge_intercept(x,w,h,shorten), (0.0,shorten))
    θ = angle(complex(dst...)-complex(src...))
    cosθ, sinθ = (cos(θ), sin(θ))
    C = 4*h/w^2
    rawx = x_shorten:0.01:(w-x_shorten)
    srcx,srcy = src
    xvals = [x*cosθ - C*x*(w-x)*sinθ + srcx for x in rawx]
    yvals = [x*sinθ + C*x*(w-x)*cosθ + srcy for x in rawx]
    return xvals, yvals
end

function get_curves(
    edges::AbstractVector,
    locs::AbstractVector;
    bendprop=0.1,
    shorten=0.2
    )

    xout, yout = (Vector{Float64}[],Vector{Float64}[])

    for (i,j) in edges
        xtmp, ytmp = get_arc_shorten(locs[i],locs[j],bendprop=bendprop,shorten=shorten)
        push!(xout,xtmp)
        push!(yout,ytmp)
    end

    return xout, yout

end

@userplot MultigraphPlot
@recipe function f(
    h::MultigraphPlot;
    relabelnodes=true,
    shorten=0.1,
    bendprop=0.1,
    seed=1
    )
    edgelist = h.args[1]

    do_relabel = (relabelnodes) || (eltype(edgelist)==Tuple{String,String})
    if do_relabel
        edgelist, rev_vmap = relabel(edgelist)
        V = length(keys(rev_vmap))
        adjmat = zeros(Int, V, V)
        for (i,j) in edgelist
            adjmat[i,j] += 1
        end
    end
    aspect_ratio --> 1
    locs = NetworkLayout.spring(adjmat,seed=seed)

    @series begin
        # seriestype := :line
        linecolor := :blue
        arrow --> true
        get_curves(edgelist, locs, shorten=shorten, bendprop=bendprop)
    end

    @series begin
        seriestype := :scatter
        markercolor := :blue
        locs
    end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

struct Stage{P}
    parameters::P

    function Stage()
        return new{Nothing}(nothing)
    end

    function Stage(parameters::P) where P
        return new{P}(parameters)
    end
end
parameter_type(stage::Stage{P}) where P = P

struct StageParameters{NT <: NamedTuple}
    names::Vector{Symbol}
    defaults::NT

    function StageParameters(; kw...)
        defaults = values(kw)
        names = collect(keys(defaults))
        NT = typeof(defaults)
        return new{NT}(names, defaults)
    end

    function StageParameters(names::Vector{Symbol}; kw...)
        defaults = values(kw)
        NT = typeof(defaults)
        return new{NT}(names, defaults)
    end
end

function parameters(stage_params::StageParameters; kw...)
    d = Dict(kw)
    params = if isempty(d)
        stage_params.defaults
    else
        filter!(p -> p.first ∈ stage_params.names, d)
        merge(stage_params.defaults, d)
    end
    if length(params) != length(stage_params.names)
        missing = filter(n -> !(n ∈ keys(params)), stage_params.names)
        isare = length(missing) == 1 ? "is" : "are"
        error("Not enough parameters specified. $(join(missing, ',')) $isare missing.")
    end
    return params
end