text:
using NeRCA
using Plots
using LinearAlgebra
using Optim
using JuMP
using GLPK
using Ipopt

const n = 1.35
const c = 2.99792458e8 / 1e9

struct SingleDUTrack
    distance
    t_closest
    z_closest
    dir_z
    t₀
end

function single_du_params(track::NeRCA.Track)
    pos = track.pos
    dir = Direction(normalize(track.dir))
    t₀ = track.time
    proj = pos ⋅ dir
    z_closest = (pos.z - dir.z*proj) / (1 - dir.z^2)
    t_closest = t₀ + (z_closest * dir.z - proj)/c
    p_t_closest = pos + c * (t_closest - t₀) * dir
    d_closest = √(p_t_closest[1]^2 + p_t_closest[2]^2)

    #SingleDUTrack(d_closest, t_closest, z_closest, dir.z, t₀)
    d_closest, t_closest, z_closest, dir.z, t₀
end

function make_cherenkov_calc(d_closest, t_closest, z_closest, dir_z, t₀)
    d_γ(z) = n/√(n^2 - 1) * √(d_closest^2 + (z-z_closest)^2 * (1 - dir_z^2))
    z -> begin
        (t_closest - t₀) + 1/c * ((z - z_closest)*dir_z + (n^2 - 1)/n * d_γ(z))
    end
end

dom_z_positions = range(120, length=18, stop=800)

track = NeRCA.Track([0, 1, 0.5], [0, 500, 500], 0)
sdp = single_du_params(track)
ccalc = make_cherenkov_calc(sdp...)
times = ccalc.(dom_z_positions) * 1e9
plot(times, dom_z_positions, seriestype = :scatter)

function make_quality_function(positions, times)
    function quality_function(d_closest, t_closest, z_closest, dir_z, t₀)
        ccalc = make_cherenkov_calc(d_closest, t_closest, z_closest, dir_z, t₀)
        expected_times = ccalc.(positions)
        return sum((times - expected_times).^2)
    end
    return quality_function
end

qfunc = make_quality_function(dom_z_positions, ccalc.(dom_z_positions))

model = Model(with_optimizer(Ipopt.Optimizer))

register(model, :qfunc, 5, qfunc, autodiff=true)

@variable(model, -1000 <= d_closest <= 1000)
@variable(model, -10000 <= t_closest <= 10000)
@variable(model, -1000 <= z_closest <= 1000)
@variable(model, -1 <= dir_z <= 1)
@variable(model, -1000 <= t₀ <= 1000)

@NLobjective(model, Min, qfunc(d_closest, t_closest, z_closest, dir_z, t₀))

optimize!(model)
termination_status(model)
value(d_closest)
value(t_closest)
value(z_closest)
value(dir_z)
value(t₀)

ccalc = make_cherenkov_calc(value(d_closest), value(t_closest), value(z_closest), value(dir_z), value(t₀))
times = ccalc.(dom_z_positions)
plot!(times, dom_z_positions, seriestype = :scatter)


text:
function templateTestTreduce1(arr::AbstractArray)::Bool
	init = randInit();
	return reduce(+, arr; init = init) == treduce(+, arr; init = init);
end

function testTreduce1()::Bool
	return templateTestTreduce1(randArray());
end

function testTreduce1Offset()
	return templateTestTreduce1(offsetArray());
end

function testTreduce1Empty()::Bool
	return templateTestTreduce1([]);
end

#####

function templateTestTreduce2(arr::AbstractArray)::Bool
	return reduce(+, arr) == treduce(+, arr);
end

function testTreduce2()::Bool
	templateTestTreduce2(randArray());
end

function testTreduce2Offset()
	templateTestTreduce2(offsetArray());
end

function testTreduce2Empty()
	treduce(identity, []);
end

ext: jl


text:
# This file is part of Intersections.jl

export binarysearch, doublingsearch, doublingsearchrev, seqsearch, seqsearchrev

"""
	binarysearch(A, x, sp=1, ep=length(A))

Finds the insertion position of `x` in `A` in the range `sp:ep`
"""
function binarysearch(A, x, sp=1, ep=length(A))
	while sp < ep
		mid = div(sp + ep, 2)
		if x <= _get_key(A, mid)
			ep = mid
		else
			sp = mid + 1
		end
	end

	x <= _get_key(A, sp) ? sp : sp + 1
end

"""
	doublingsearch(A, x, sp=1, ep=length(A))

Finds the insertion position of `x` in `A`, starting at `sp`
"""
function doublingsearch(A, x, sp=1, ep=length(A))
	p = 0
    i = 1

    while sp+i <= ep && _get_key(A, sp+i) < x
		p = i
		i += i
    end

    binarysearch(A, x, sp + p, min(ep, sp+i))
end

"""
	doublingsearchrev(A, x, sp=1, ep=length(A))

Finds the insertion position of `x` in `A`, starting at the end
"""
function doublingsearchrev(A, x, sp=1, ep=length(A))
	x > _get_key(A, ep) && return ep + 1

	i = 1
	p = ep
    while p >= sp && x <= _get_key(A, p)
		ep = p
		p -= i
		i += i
    end

	binarysearch(A, x, max(p, sp), ep)
end

"""
	seqsearchrev(A, x, sp=1, ep=length(A))

Reverse sequential search, i.e., it starts from `ep` to `sp`
"""
function seqsearchrev(A, x, sp=1, ep=length(A))
	pos = ep + 1
	while pos > sp && x <= _get_key(A, pos-1)
        pos -= 1
    end

	pos
end

"""
	searchrev(A, x, sp=1, ep=length(A))

Sequential search, i.e., it starts from `sp` to `ep`
"""
function seqsearch(A, x, sp=1, ep=length(A))
	x > _get_key(A, ep) && return ep + 1
	while sp <= ep && _get_key(A, sp) < x
        sp += 1
    end

	sp > ep ? ep : sp
end


text:
#=
Implements the Kalman filter for a linear Gaussian state space model.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-07-29

References
----------

http://quant-econ.net/jl/kalman.html

TODO: Do docstrings here after implementing LinerStateSpace
=#

type Kalman
    A
    G
    Q
    R
    k
    n
    cur_x_hat
    cur_sigma
end

# Initializes current mean and cov to zeros
function Kalman(A, G, Q, R)
    k = size(G, 1)
    n = size(G, 2)
    xhat = n == 1 ? zero(eltype(A)) : zeros(n)
    Sigma = n == 1 ? zero(eltype(A)) : zeros(n, n)
    return Kalman(A, G, Q, R, k, n, xhat, Sigma)
end

function set_state!(k::Kalman, x_hat, Sigma)
    k.cur_x_hat = x_hat
    k.cur_sigma = Sigma
    Void
end

"""
Updates the moments (cur_x_hat, cur_sigma) of the time t prior to the
time t filtering distribution, using current measurement y_t.
The updates are according to
    x_{hat}^F = x_{hat} + Sigma G' (G Sigma G' + R)^{-1}
                    (y - G x_{hat})
    Sigma^F = Sigma - Sigma G' (G Sigma G' + R)^{-1} G
                Sigma

#### Arguments

- `k::Kalman` An instance of the Kalman filter
- `y` The current measurement

"""
function prior_to_filtered!(k::Kalman, y)
    # simplify notation
    G, R = k.G, k.R
    x_hat, Sigma = k.cur_x_hat, k.cur_sigma

    # and then update
    if k.k > 1
        reshape(y, k.k, 1)
    end
    A = Sigma * G'
    B = G * Sigma' * G' + R
    M = A * inv(B)
    k.cur_x_hat = x_hat + M * (y - G * x_hat)
    k.cur_sigma = Sigma - M * G * Sigma
    Void
end

"""
Updates the moments of the time t filtering distribution to the
moments of the predictive distribution, which becomes the time
t+1 prior

#### Arguments
- `k::Kalman` An instance of the Kalman filter

"""
function filtered_to_forecast!(k::Kalman)
    # simplify notation
    A, Q = k.A, k.Q
    x_hat, Sigma = k.cur_x_hat, k.cur_sigma

    # and then update
    k.cur_x_hat = A * x_hat
    k.cur_sigma = A * Sigma * A' + Q
    Void
end

"""
Updates cur_x_hat and cur_sigma given array `y` of length `k`.  The full
update, from one period to the next

#### Arguments

- `k::Kalman` An instance of the Kalman filter
- `y` An array representing the current measurement

"""
function update!(k::Kalman, y)
    prior_to_filtered!(k, y)
    filtered_to_forecast!(k)
    Void
end

function stationary_values(k::Kalman)
    # simplify notation
    A, Q, G, R = k.A, k.Q, k.G, k.R

    # solve Riccati equation, obtain Kalman gain
    Sigma_inf = solve_discrete_riccati(A', G', Q, R)
    K_inf = A * Sigma_inf * G' * inv(G * Sigma_inf * G' + R)
    return Sigma_inf, K_inf
end


text:
mutable struct Position
    x::Int
    y::Int
end

struct InventoryItem
    id::Int
    quantity::Int
end

mutable struct Player
    pos::Position
    char::Int
    emp::Int
    inv::Vector{InventoryItem}
    rank::Int
    id::String
end

const PlayerCache = Dict{String, Player}


text:
"""
    approach_wall!(initxy, dx, dy, stepsize) -> initxy::Array{Float64,1}

Shorten proposed vector into a parylene wall face with no catalase.
"""
function approach_wall!(initxy::Array{Float64,1},
                        dx::Float64, dy::Float64,
                        stepsize::Float64)
    newx = initxy[1] + stepsize*dx
    newy = initxy[2] + stepsize*dy
    # Recursively shorten the proposed vector
    if inwalls(newx, newy)
        return approach_wall!(initxy, dx, dy, stepsize * 0.5)
    else
        initxy[1] = newx
        initxy[2] = newy
        return initxy
    end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
struct RoleSelection <: AbstractSelection
    r::Symbol
end
if_role(r::Symbol) = RoleSelection(r)
if_role(r::Vector{Symbol}) = RoleSelection.(r)
resolve(df, s::RoleSelection) = resolve(df, r.roles[s.r])

resolveroles(df, d) = Dict(k => first.(getselectionpairs(df, v)) for (k,v) in d)
add_roles(r, roles...) = merge!(r.roles, Dict(roles))


text:
## This example comes from https://diffeqflux.sciml.ai/dev/Flux/

using ComponentArrays
using OrdinaryDiffEq
using Plots
using UnPack

using DiffEqFlux: sciml_train
using Flux: glorot_uniform, ADAM, σ, relu
using Optim: LBFGS

struct MeasurementNoise{T}
    sigma::T
end
Base.:(+)(array, mn::MeasurementNoise) = array .+ randn(size(array)...).*mn.sigma

# Problem setup
u0 = Float32[2.; 0.]
datasize = 50
tspan = (0.0f0, 7f0)
tspan2 = (0.0f0, 25f0)

# Make truth data
function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end

t = range(tspan[1], tspan[2], length = datasize)
# t = Float32.(vcat(range(0.0, 0.9, length=10), 10 .^ range(log10(tspan[1]+1), log10(tspan[2]), length=datasize-10)))
t = Float32.([0; 10 .^ range(log10(tspan[1] + 0.01), log10(tspan[2]), length=datasize-1)])
prob = ODEProblem(trueODEfunc, u0, tspan2)
ode_sol = solve(prob, Tsit5())
ode_data = Array(ode_sol(t)) + MeasurementNoise(0.1)

# Function for creating neural layer components
neural_layer(in, out) = ComponentArray{Float32}(W=glorot_uniform(out, in), b=zeros(out))

# Dense neural layer function
dense(layer, activation=identity) = u -> activation.(layer.W * u + layer.b)

# Neural ODE function
dudt(u, p, t) = u.^3 |> dense(p.L1, σ) |> dense(p.L2)

prob = ODEProblem(dudt, u0, tspan2)

# Create optimization parameter vector
n_hidden = 50
layers = (L1=neural_layer(2, n_hidden), L2=neural_layer(n_hidden, 2))
θ = ComponentArray(u=u0, p=layers)

# Prediction and loss functions
predict_n_ode(θ) = Array(solve(prob, Tsit5(), u0=θ.u, p=θ.p, saveat=t))
full_sol(θ) = solve(prob, Tsit5(), u0=θ.u, p=θ.p)

function loss_n_ode(θ)
    pred = predict_n_ode(θ)

    loss = sum(abs2, ode_data .- pred)/datasize + 0.1*(sum(abs, θ.p)/length(θ.p))
    return loss, pred
end
loss_n_ode(θ)

# anim = Animation()

 # Callback function to observe training
cb = function (θ, loss, pred; doplot=false)
    display(loss)
    # plot current prediction against data
    train_sol = full_sol(θ)

    pl_1 = plot(train_sol, vars=1)
    plot!(pl_1, ode_sol, vars=1)
    scatter!(pl_1, t, pred[1,:])
    scatter!(pl_1, t, ode_data[1,:], legend=false)

    pl_2 = plot(train_sol, vars=2)
    plot!(pl_2, ode_sol, vars=2)
    scatter!(pl_2, t, pred[2,:])
    scatter!(pl_2, t, ode_data[2,:], legend=false)

    pl_3 = plot(train_sol, vars=(1,2), label = "prediction")
    plot!(pl_3, ode_sol, vars=(1,2), label = "truth")
    scatter!(pl_3, pred[1,:], pred[2,:], label = "predicted data")
    scatter!(pl_3, ode_data[1,:], ode_data[2,:], label = "measured data")
    plot!(pl_3, hcat(pred[1,:], ode_data[1,:])', hcat(pred[2,:], ode_data[2,:])', label = false, color=:lightgray, legend=:bottomright)

    display(plot(plot(pl_1, pl_2, layout=(2,1), size=(400,500)), pl_3, layout=(1,2), size=(950,500)))
    # frame(anim)
    return false
end

# Display the ODE with the initial parameter values.
cb(θ, loss_n_ode(θ)...)

data = Iterators.repeated((), 1000)

res1 = sciml_train(loss_n_ode, θ, ADAM(0.05); maxiters=100, save_best=true, cb=cb)
cb(res1.minimizer, loss_n_ode(res1.minimizer)...; doplot=true)

res2 = sciml_train(loss_n_ode, res1.minimizer, LBFGS(), cb=cb)
cb(res2.minimizer, loss_n_ode(res2.minimizer)...; doplot=true)

# gif(anim, "DiffEqFlux.gif", fps=15)


text:
@testset "map/" begin

    include("map_translate.jl")
    include("map_affine.jl")
    include("map_quadratic.jl")
    include("transform_src_insert_dest_map.jl")

end


text:
using LinearAlgebra, DataFrames,Optim, ForwardDiff, BenchmarkTools,Distributions,
    Expectations, QuantEcon, Statistics, GLM
using Distributions: invsqrt2π, log2π, sqrt2, invsqrt2
using DCDC
using Test
using Distributed, Suppressor
using Plots

begin "Dynamic Decision Process"
    σ₀ = 1;
    β = 0.8;
    nM = 50;
    nT = 5;
    ddc = DynamicDecisionProcess(σ₀,0.8);
    plot(ddc.ValueFn.xdata,ddc.ValueFn.y);
    plot!(ddc.PolicyFn.xdata,ddc.PolicyFn.y);
    data = simulate_ddc(nM,nT,ddc);
end

begin "GLM"
    x = randn(300,3);
    y = x * [1,3,4] + randn(300);
    testdata = DataFrame(x);
    testdata.y = y;
    ols = lm(@formula(y ~ x1+x2+x3), testdata);
    @show stderror(ols);
    # @show ols.model.pp.beta0;
end

using ForwardDiff
begin "ForwardDiff"
    h(x) = sin(x[1]) + x[1] * x[2] + sinh(x[1] * x[2]) # multivariate.
    x = [1.4 2.2]
    @show ForwardDiff.gradient(h,x) # use AD, seeds from x
    #Or, can use complicated functions of many variables
    f(x) = sum(sin, x) + prod(tan, x) * sum(sqrt, x)
    g = (x) -> ForwardDiff.gradient(f, x); # g() is now the gradient
    @show g(rand(20)); # gradient at a random point

    function squareroot(x) #pretending we don't know sqrt()
        z = copy(x) # Initial starting point for Newton’s method
        while abs(z*z - x) > 1e-13
            z = z - (z*z-x)/(2z)
        end
        return z
    end
    sqrt(2.0)
    dsqrt(x) = ForwardDiff.derivative(squareroot, x)
    dsqrt(2.0)
end

begin "Multi-variate function: Can we forwarddiff?"
    func = (x,y) -> x^3 + y^2;
    kk  = y->(x->func(x,y));
    a,b = randn(2);
    ans1 = kk(b)(a);
    ans2 = func(a,b);
    @show ans1 == ans2;
    dfunx = (x,y) -> ForwardDiff.derivative(z->kk(y)(z),x);
    dfunx2 = (x,y) -> ForwardDiff.derivative(z->func(z,y),x);
    deriv1 =dfunx(a,b);
    deriv2 = dfunx2(a,b);
    @show deriv1 == deriv2;
end

begin "Test transition derivative"
    dtrans_x = (x,y) -> ForwardDiff.derivative( z-> ddc.trans(z,y),x) #Derivative w.r.t first component
    dtrans_c = (x,y) -> ForwardDiff.derivative( z-> ddc.trans(y,z),x) #Derivative w.r.t second component
    dtrans_x(1,3) == 1.05
end

using Flux
using Flux.Tracker
using Flux.Tracker: update!
begin "Test flux: graident and jacobian"
    f(x) = 3x^2 + 2x + 1
    # df/dx = 6x + 2
    df(x) = Tracker.gradient(f, x)[1]
    df(2);
    A = rand(2,2);
    f(x) = A * x
    x0 = [0.1, 2.0]
    f(x0)
    Flux.jacobian(f, x0)
end

using TimerOutputs
begin
    # Create the timer object
    to = TimerOutput()
    # Time something with an assigned label
    @timeit to "sleep" sleep(0.3)
    # Data is accumulated for multiple calls
    for i in 1:100
        @timeit to "loop" 1+1
    end
    # Nested sections are possible
    @timeit to "nest 1" begin
        @timeit to "nest 2" begin
            @timeit to "nest 3.1" rand(10^3)
            @timeit to "nest 3.2" rand(10^4)
            @timeit to "nest 3.3" rand(10^5)
        end
        rand(10^6)
    end
end

# Expectations
begin "Test Expectatoin"
    dist = Normal();
    E = expectation(dist, Gaussian; n = 301)
    f = x -> x^2
    expectation(f, dist)
end

using Optim
using Optim: converged, maximum, maximizer, minimizer, iterations #some extra functions
begin "Test Optim"
      result = optimize(x-> x^2, -2.0, 1.0)
      @show converged(result) || error("Failed to converge in $(iterations(result)) iterations")
      @show xmin = result.minimizer
      @show result.minimum
end


text:
# TODO (reactivate)

# n = 1
# T = Float64
#
# function F(x::Vector, b::Vector)
#     b[:] = x.^2
# end
#
# function J(x::Vector, A::Matrix)
#     A[:,:] = 2x
# end
#
#
# x = ones(T, n)
# nl = NewtonSolver(x, F)
# solve!(nl)
# # println(nl.status.i, ", ", nl.status.rₐ,", ",  nl.status.rᵣ,", ",  nl.status.rₛ)
# @test_approx_eq_eps(nl.x, zeros(T, n), 1E-7)
#
#
# x = ones(T, n)
# nl = NewtonSolver(x, F, J=J)
# solve!(nl)
# # println(nl.status.i, ", ", nl.status.rₐ,", ",  nl.status.rᵣ,", ",  nl.status.rₛ)
# @test_approx_eq_eps(nl.x, zeros(T, n), 1E-7)


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 12/15/2014

module AbstractEncounterDBNImpl

export AbstractEncounterDBN

abstract AbstractEncounterDBN

end

module AbstractEncounterDBNInterfaces

export
       # common
       get,
       getInitialState

# common

function get()
end

function getInitialState()
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
module TestUtilities

using Test
using MLJBase

export @testset_accelerated, include_everywhere, @test_mach_sequence,
    @test_model_sequence

using ComputationalResources
using ComputationalResources: CPUProcesses

macro testset_accelerated(name::String, var, ex)
    testset_accelerated(name, var, ex)
end
macro testset_accelerated(name::String, var, opts::Expr, ex)
    testset_accelerated(name, var, ex; eval(opts)...)
end
function testset_accelerated(name::String, var, ex; exclude=[])
    final_ex = quote
       local $var = CPU1()
        @testset $name $ex
    end

    resources = AbstractResource[CPUProcesses(), CPUThreads()]

    for res in resources
        if any(x->typeof(res)<:x, exclude)
            push!(final_ex.args, quote
               local $var = $res
               @testset $(name*" ($(typeof(res).name))") begin
                   @test_broken false
               end
            end)
        else
            push!(final_ex.args, quote
               local $var = $res
               @testset $(name*" ($(typeof(res).name))") $ex
            end)
        end
    end
    # preserve outer location if possible
    if ex isa Expr && ex.head === :block && !isempty(ex.args) && ex.args[1] isa LineNumberNode
        final_ex = Expr(:block, ex.args[1], final_ex)
    end
    return esc(final_ex)
end

function sedate!(fit_ex)
    kwarg_exs = filter(fit_ex.args) do arg
        arg isa Expr && arg.head == :kw
    end
    keys = map(kwarg_exs) do arg
        arg.args[1]
    end
    :verbosity in keys &&
        error("You cannot specify `verbosity` in @test_mach_sequence "*
              "or @test_model_sequence. ")
    push!(fit_ex.args, Expr(:kw, :verbosity, -5000))
    return fit_ex
end

macro test_mach_sequence(fit_ex, sequence_exs...)
    sedate!(fit_ex)
    seq = gensym(:sequence)
    esc(quote
        MLJBase.flush!(MLJBase.MACHINE_CHANNEL)
        $fit_ex
        local $seq = MLJBase.flush!(MLJBase.MACHINE_CHANNEL)
        # for s in $seq
        #     println(s)
        # end
        @test $seq in [$(sequence_exs...)]
    end)
end

# function weakly_in(object:Tuple{Symbol,Model}, itr)
#     for tup in itr
#         tup[1] === object[1] && tup[2] == tup

macro test_model_sequence(fit_ex, sequence_exs...)
    sedate!(fit_ex)
    seq = gensym(:sequence)
    esc(quote
        MLJBase.flush!(MLJBase.MACHINE_CHANNEL)
        $fit_ex
        local $seq = map(MLJBase.flush!(MLJBase.MACHINE_CHANNEL)) do tup
            (tup[1], tup[2].model)
        end
        # for s in $seq
        #     println(s)
        # end
        @test $seq in [$(sequence_exs...)]
    end)
end

end


text:
using Reactive, SignalView, GeometryTypes, Colors, GLWindow, GLVisualize, GLAbstraction
import Images

window = glscreen()

# loadasset loads data from the GLVisualize asset folder and is defined as
# FileIO.load(assetpath(name))
doge = Images.restrict(Images.restrict(loadasset("doge.png"))).data

# Convert to RGBA{Float32}. Float for filtering and 32 because it fits the GPU better
img = map(RGB{Float32}, doge)
# create a slider that goes from 1-20 in 0.1 steps
sigma = Signal(1f0)

# performant conversion to RGBAU8, implemted with a functor
# in 0.5 anonymous functions offer the same speed, so this wouldn't be needed
"""
Applies a gaussian filter to `img` and converts it to RGBA{U8}
"""
function gauss(img, sigma)
	Images.imfilter_gaussian(img, [sigma, sigma])
end

function laplacian(img)
    map(RGB{Float32}, Images.imfilter(img, Images.imlaplacian()))
end

function postprocess(img, scale)
    map(img) do color
        scaled = color .* scale
        RGB{Float32}(1f0-scaled.r, 1f0-scaled.g, 1f0-scaled.b)
    end
end
gaussed = map(gauss, Signal(img), sigma)
println("gaussed ", typeof(gaussed))

edges  = map(laplacian, gaussed)
println("edges ", typeof(edges))

dilated = map(postprocess, edges, Signal(4.0f0))
println("dilated ", typeof(dilated))

#pass the signal root to view_signal
view_signal(sigma, window)

renderloop(window)


text:
"""
    getindex(g, label)

Alias for [`get_vertex(g, label)`](@ref), intended for dictionary-like use.
"""
Base.getindex(g::AbstractDataGraph{T,VL}, label::VL) where {T,VL} = get_vertex(g, label)


text:
struct SpinOrbitInteraction <: OneBodyOperator end

atomic_hamiltonian(::Atom{T,B,O,TC,CV,P}) where {T,B,O,TC,CV,P<:Union{RelativisticPseudoPotential,PseudoPotentials.ScalarSORelativisticPseudoPotential}} =
    FieldFreeOneBodyHamiltonian() + CoulombInteraction() + SpinOrbitInteraction()

Base.show(io::IO, ::SpinOrbitInteraction) = write(io, "V̂ₛₒ")

"""
    mⱼ(o::SpinOrbital)

Return the total mⱼ that `o` couples to in LS coupling, i.e. `mℓ+mₛ`.
"""
mⱼ(o::SpinOrbital) = sum(o.m)

function Base.iszero(me::OrbitalMatrixElement{1,A,SpinOrbitInteraction,B}) where {A<:SpinOrbital,B<:SpinOrbital}
    a = me.a[1]
    b = me.b[1]

    ℓa = first(jmⱼ(a))
    ℓb = first(jmⱼ(b))
    # This is correct for both nonrelativistic as well as relativistic
    # spin-orbitals, since for the former, first(jmⱼ(o)),mⱼ(o) returns
    # ℓ_o,m_j_o, and for the latter j_o,m_j_o. The pseudo-potentials
    # are formed using 𝒫(ℓjmⱼ) projectors, which are diagonal in the
    # basis of relativistic spin-orbitals (except for the principal
    # quantum number n), but block-diagonal in the basis of
    # non-relativistic spin-orbitals, with the entries of the blocks
    # given by the Clebsch–Gordan coefficients ⟨ℓm_ℓsmₛ|ℓsjmⱼ⟩.
    # Additionally, for relativistic spin-orbitals, we must check ℓ.
    a.orb.ℓ ≠ b.orb.ℓ || ℓa ≠ ℓb || iszero(ℓa) || mⱼ(a) != mⱼ(b)
end

function get_operator(::SpinOrbitInteraction, atom::Atom, a::SpinOrbital, b::SpinOrbital; kwargs...)
    if a == b && atom.potential isa RelativisticPseudoPotential
        # In case of two-component relativistic pseudopotentials, the
        # orbital-diagonal spin–orbit contribution is already
        # accounted for.
        return 0
    end
    R = radial_basis(atom)
    r = locs(R)
    Vₛₒ = spin_orbit_potential(atom.potential, r, a, b)
    operator(Diagonal(Vₛₒ), R)
end


text:
using Test
using Dash

@testset "default args" begin
    app = dash()
    @test app.root_path == pwd()
    @test isempty(app.config.external_stylesheets)
    @test isempty(app.config.external_scripts)
    @test app.config.requests_pathname_prefix == "/"
    @test app.config.routes_pathname_prefix == "/"
    @test app.config.assets_folder == "assets"
    @test app.config.assets_url_path == "assets"
    @test app.config.assets_ignore == ""

    @test app.config.serve_locally == true
    @test app.config.suppress_callback_exceptions == false
    @test app.config.eager_loading == false

    @test isempty(app.config.meta_tags)
    @test app.index_string == Dash.default_index
    @test app.config.assets_external_path == nothing

    @test app.config.include_assets_files == true
    @test app.config.show_undo_redo == false
end

@testset "setted args" begin
    app = dash(; external_stylesheets=["https://test.css"])
    @test app.config.external_stylesheets == ["https://test.css"]

    app = dash(; external_stylesheets=[Dict("url" => "https://test.css", "integrity" => "integrity")])
    @test app.config.external_stylesheets == [
        Dict("url" => "https://test.css", "integrity" => "integrity")
        ]

    app = dash(; external_scripts = ["http://test.js"])
    @test app.config.external_scripts == ["http://test.js"]

    app = dash(; external_scripts=[Dict("url" => "https://test.js", "integrity" => "integrity")])
    @test app.config.external_scripts == [
        Dict("url" => "https://test.js", "integrity" => "integrity")
        ]

    app = dash(; url_base_pathname = "/base/")
    @test app.config.url_base_pathname == "/base/"
    @test app.config.requests_pathname_prefix == "/base/"
    @test app.config.routes_pathname_prefix == "/base/"

    app = dash(; requests_pathname_prefix = "/prefix/")
    @test app.config.requests_pathname_prefix == "/prefix/"
    @test app.config.routes_pathname_prefix == "/"

    app = dash(; routes_pathname_prefix = "/prefix/")
    @test app.config.requests_pathname_prefix == "/prefix/"
    @test app.config.routes_pathname_prefix == "/prefix/"

    app = dash(; requests_pathname_prefix = "/reg/prefix/", routes_pathname_prefix = "/prefix/")

    @test app.config.requests_pathname_prefix == "/reg/prefix/"
    @test app.config.routes_pathname_prefix == "/prefix/"

    @test_throws ErrorException app = dash(; url_base_pathname = "/base")
    @test_throws ErrorException app = dash(; url_base_pathname = "base/")
    @test_throws ErrorException app = dash(; url_base_pathname = "/", routes_pathname_prefix = "/prefix/")
    @test_throws ErrorException app = dash(; requests_pathname_prefix = "/prefix")
    @test_throws ErrorException app = dash(; routes_pathname_prefix = "/prefix")
    @test_throws ErrorException app = dash(; requests_pathname_prefix = "prefix/")
    @test_throws ErrorException app = dash(; routes_pathname_prefix = "prefix/")
    @test_throws ErrorException app = dash(; requests_pathname_prefix = "/reg/prefix/", routes_pathname_prefix = "/ix/")

    app = dash(; assets_folder = "images")
    @test app.config.assets_folder == "images"

    app = dash(; assets_url_path = "/images")
    @test app.config.assets_url_path == "images"

    app = dash(; assets_ignore = "ignore")
    @test app.config.assets_ignore == "ignore"

    app = dash(; serve_locally = false)
    @test app.config.serve_locally == false

    app = dash(; suppress_callback_exceptions = true)
    @test app.config.suppress_callback_exceptions == true

    app = dash(; eager_loading = false)
    @test app.config.eager_loading == false

    app = dash(; meta_tags = [Dict(["name"=>"test", "content" => "content"])])
    @test app.config.meta_tags == [Dict(["name"=>"test", "content" => "content"])]

    @test_throws ErrorException app = dash(; index_string = "")
    app = dash(; index_string = "{%app_entry%}{%config%}{%scripts%}")
    @test app.index_string == "{%app_entry%}{%config%}{%scripts%}"

    app = dash(; assets_external_path = "external")
    @test app.config.assets_external_path == "external"

    app = dash(; include_assets_files = true)
    @test app.config.include_assets_files == true

    app = dash(; show_undo_redo = false)
    @test app.config.show_undo_redo == false

end


text:
module Divergences

using StatsFuns
using LoopVectorization
using Parameters
import LoopVectorization: vifelse
import VectorizationBase: andmask
import Distances: PreMetric

abstract type AbstractDivergence end
abstract type Divergence <: AbstractDivergence end
abstract type AbstractModifiedDivergence <: AbstractDivergence end

struct CressieRead{T} <: Divergence
    α::T
    function CressieRead(α::T) where T<:Real
        #@assert isempty(findall((in)([-1, 0]), α)) "CressieRead is defined for all α != {-1,0}"
        new{T}(α)
    end
end

struct ChiSquared  <: Divergence end
struct KullbackLeibler  <: Divergence end
struct ReverseKullbackLeibler <: Divergence end
struct Hellinger <: Divergence end

struct ModifiedDivergence{D, T} <: AbstractModifiedDivergence
    d::D
    m::NamedTuple{(:γ₀, :γ₁, :γ₂, :ρ), Tuple{T, T, T, T}}
end

struct FullyModifiedDivergence{D, T} <: AbstractModifiedDivergence
    d::D
    m::NamedTuple{(:γ₀, :γ₁, :γ₂, :ρ, :g₀, :g₁, :g₂, :φ), Tuple{T, T, T, T, T, T, T, T}}
end

function ModifiedDivergence(D::Divergence, ρ::Real)
    @assert ρ > 1 "A ModifiedDivergence requires ρ > 1"
    γ₀ = eval(D, [ρ])[1]
    γ₁ = gradient(D, [ρ])[1]
    γ₂ = hessian(D, [ρ])[1]
    ModifiedDivergence(D, (γ₀=γ₀, γ₁=γ₁, γ₂=γ₂, ρ=ρ))
end

function FullyModifiedDivergence(D::Divergence, φ::Real, ρ::Real)
    @assert ρ > 1 "A ModifiedDivergence requires ρ > 1"
    @assert φ < 1 && φ > 0 "A ModifiedDivergence requires  φ ∈ (0,1)"
    γ₀ = eval(D, [ρ])[1]
    γ₁ = gradient(D, [ρ])[1]
    γ₂ = hessian(D, [ρ])[1]
    g₀ = eval(D, [φ])[1]
    g₁ = gradient(D, [φ])[1]
    g₂ = hessian(D, [φ])[1]
    FullyModifiedDivergence(D, (γ₀=γ₀, γ₁=γ₁, γ₂=γ₂, ρ=ρ, g₀=g₀, g₁=g₁, g₂=g₂, φ=φ))
end

const 𝒦ℒ=KullbackLeibler
const ℬ𝓊𝓇ℊ=ReverseKullbackLeibler
const 𝒞ℛ=CressieRead
const ℋ𝒟=Hellinger
const χ²=ChiSquared
include("divs.jl")

export
    # KL
    KullbackLeibler,
    # RKL
    ReverseKullbackLeibler,
    # HD
    Hellinger,
    # CR
    CressieRead,
    #
    ChiSquared,
    # Modified
    ModifiedDivergence,
    # FullyModified
    FullyModifiedDivergence,
    # Abbr.
    𝒦ℒ,
    ℬ𝓊𝓇ℊ,
    𝒞ℛ,
    ℋ𝒟,
    χ²
end


text:
using Requests

# TODO: Should probably move this outside of this repo as it forces us to install Requests
# which isn't usually needed.

function travis_benchmark(slug, token, pkg, baseline, candidate)
    slug = replace(slug, "/", "%2F")

    # Using the Travis API v3 beta. Documentation currently unreleased.
    # https://docs.travis-ci.com/user/triggering-builds
    body = Dict(
        "request" => Dict(
            "message" => "Performance benchmark test: \"$baseline\" vs. \"$candidate\"",
            "branch" => candidate,
            "config" => Dict(
                "os" => "linux",  # Guess we don't have to restrict the OS
                "julia" => "nightly",  # Need LibGit2 support (in 0.5+)
                "env" => Dict(
                    "matrix" => ["TEST=performance"],  # For future use
                ),
                "install" => [
                    "git remote set-branches origin $baseline $candidate",
                    "if [[ -a .git/shallow ]]; then git fetch --unshallow; fi",
                    "git branch $baseline origin/$baseline",  # Now both candidate and baseline are local branches
                    "julia -e 'Pkg.clone(pwd()); Pkg.build(\"$pkg\");'",
                    "julia -e 'Pkg.clone(\"https://github.com/invenia/PkgBenchmarks.jl\");'",  # Temporary
                ],
                "script" => [
                    # Note: Due to using `Pkg.clone` the origin actually refers to local branches in the `pwd()` repo.
                    "julia -e 'using PkgBenchmarks, BenchmarkTools; a, b = benchmark(\"$pkg\", \"origin/$baseline\", \"origin/$candidate\"); j = judge(minimum(a), minimum(b)); println(leaves(j)); showall(a); showall(b); exit(!isempty(regressions(j)))'",
                ],
                "after_success" => "echo hi",  # TODO: find a better way to override "after_success"
            ),
        )
    )

    post(
        "https://api.travis-ci.org/repo/$slug/requests",
        json = body,
        headers = Dict(
            "Travis-API-Version" => "3",
            "Authorization" => "token \"$token\"",
        ),
    )
end


text:
"""
    PODBasis{T}(coefficients::Matrix{T}, modes::Matrix{T})

Datastructure to store a Proper Orthogonal Decomposition basis.

The original data which the POD basis is representing can be reconstructed by right-
multiplying the modes with the coefficients, i.e. `A = P.modes*P.coefficients`.
"""
struct PODBasis{T}
    coefficients::Matrix{T}
    modes::Matrix{T}
end


text:
include("../oracles/portfolio_oracle.jl")

using Test, Random, LinearAlgebra

@testset "Portfolio Oracle" begin
    Random.seed!(342)

    d = 10
    factors = 4
    F = 0.05*rand(d, factors)

    mu = rand(d)
    r = mu + F*randn(factors) + 0.2*randn(d)

    println("R is : $r")

    Sigma = F*F' + Matrix(0.2I, d, d)

    (z1, w1) = portfolio_simplex_oracle_jump_basic(-r, Sigma, 0.1)

    oracle = portfolio_simplex_jump_setup(Sigma, 0.1)
    (z2, w2) = oracle(-r)

    @test isapprox(z1, z2; atol = 0.0001)
    @test isapprox(w1, w2; atol = 0.0001)
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
using Test
using CartesianMesh
#using Revise
using CutCell
include("useful_routines.jl")

cm = CutCell.CellMap([0.,2.],[3.,3.])
@test CutCell.dimension(cm) == 2
@test CutCell.jacobian(cm) == [1.5,0.5]
@test allapprox(cm([0.,0.]),[1.5,2.5])
@test allapprox(cm([1.,0.]),[3.,2.5])

mesh = UniformMesh([0.,0.],[1.,1.],[2,2])
@test allequal(CutCell.elements_per_mesh_side(mesh),[2,2])
cellmaps = CutCell.cell_maps(mesh)
cm1 = CutCell.CellMap([0.,0.],[0.5,0.5])
cm2 = CutCell.CellMap([0.,0.5],[0.5,1.])
cm3 = CutCell.CellMap([0.5,0.],[1.,0.5])
cm4 = CutCell.CellMap([0.5,0.5],[1.,1.])
@test allequal(cellmaps,[cm1,cm2,cm3,cm4])

cm = CutCell.CellMap([0.,0.],[1.,1.])
@test CutCell.determinant_jacobian(cm) ≈ 0.25

cm = CutCell.CellMap([1.,1.],[5.,8.])
@test CutCell.determinant_jacobian(cm) ≈ 7.0

cm = CutCell.CellMap([2.,1.],[4.,2.])
@test allapprox(CutCell.inverse(cm,[3,1.5]),[0.,0.])
@test allapprox(CutCell.inverse(cm,[4,2.]),[1.,1.])
@test allapprox(CutCell.inverse(cm,[2,1]),[-1.,-1.])


text:
# function to run a Julia script outside of the current environment
function script(code; wrapper=``, args=``)
    if Base.JLOptions().project != C_NULL
        args = `$args --project=$(unsafe_string(Base.JLOptions().project))`
    end
    mktemp() do path, io
        write(io, code)
        flush(io)
        cmd = `$wrapper $(Base.julia_cmd()) $args $path`
        # redirect stderr to stdout to have it picked up by Weave.jl
        run(pipeline(ignorestatus(cmd), stderr=stdout))
    end
    nothing
end


text:
@testset "TX hash" begin

    chain = make_chain(2)
    blk1 = Block(chain[0])
    tx1 = blk1.transactions[1]
    hash1 = double_sha256(tx1)
    @test string(hash1, base = 16) == "4a5e1e4baab89f3a32518a88c31bc87f618f76673e2cc77ab2127b7afdeda33b"

end


text:
# ArbComplex structs hold complex balls given as ArbReal pairs

mutable struct ArbComplex{P}  <: Number  # P is the precision in bits
                         #      real midpoint
    real_mid_exp::Int    # fmpz         exponent of 2 (2^exp)
    real_mid_size::UInt  # mp_size_t    nwords and sign (lsb holds sign of significand)
    real_mid_d1::UInt    # significand  unsigned, immediate value or the initial span
    real_mid_d2::UInt    #   (d1, d2)   the final part indicating the significand, or 0
                         #      real radius
    real_rad_exp::Int    # fmpz       exponent of 2 (2^exp)
    real_rad_man::UInt   # mp_limb_t  radius, unsigned by definition
                         #      imaginary midpoint
    imag_mid_exp::Int    # fmpz         exponent of 2 (2^exp)
    imag_mid_size::UInt  # mp_size_t    nwords and sign (lsb holds sign of significand)
    imag_mid_d1::UInt    # significand  unsigned, immediate value or the initial span
    imag_mid_d2::UInt    #   (d1, d2)   the final part indicating the significand, or 0
                         #      imaginary radius
    imag_rad_exp::Int    # fmpz       exponent of 2 (2^exp)
    imag_rad_man::UInt   # mp_limb_t  radius, unsigned by definition

    function ArbComplex{P}() where {P}
        z = new{P}()
        ccall(@libarb(acb_init), Cvoid, (Ref{ArbComplex},), z)
        finalizer(clear_acb, z)
        return z
    end
end

clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x)

ArbComplex{P}(x::ArbComplex{P}) where {P} = x
ArbComplex(x::ArbComplex{P}) where {P} = x

ArbComplex{P}(x::Missing) where {P} = missing
ArbComplex(x::Missing) = missing

ArbComplex(x, y, prec::Int) = prec>=MINIMUM_PRECISION ? ArbComplex{workingbits(prec)}(x, y) : throw(DomainError("bit precision $prec < $MINIMUM_PRECISION"))

function ArbComplex(x::T, y::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number}
    if base === 10
        digits = digits > 0 ? bits4digits(digits) : (bits > 0 ? bits : DEFAULT_PRECISION.x)
    elseif base === 2
        digits = bits > 0 ? bits : (digits > 0 ? digits : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbComplex(x, y, digits)
end

function ArbComplex(x::T1, y::T2; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T1<:Number, T2<:Number}
    ArbComplex(promote(x, y)...,)
end

const Analytic = Cint(0) # prefer the non-analytic versions

@inline sign_bit(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size)
@inline sign_bits(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size), isodd(x.imag_mid_size)
@inline sign_bit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = isodd(x.real_mid_size)
@inline sign_bit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = isodd(x.imag_mid_size)

# ArbComplex{P}(sqrt(2.0))
function ArbComplex{P}(rea::Float64) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_d), Cvoid, (Ref{ArbComplex}, Float64), z, rea)
    return z
end

function ArbComplex{P}(rea::Cint) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_si), Cvoid, (Ref{ArbComplex}, Cint), z, rea)
    return z
end

function ArbComplex{P}(rea::ArbFloat{P}) where {P}
    x = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, rea)
    y = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), y, x)
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_arb), Cvoid, (Ref{ArbComplex}, Ref{ArbReal}), z, y)
    return z
end

function ArbComplex{P}(rea::ArbReal{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_arb), Cvoid, (Ref{ArbComplex}, Ref{ArbReal}), z, rea)
    return z
end

# ArbComplex{P}(sqrt(2.0), inv(17.0))
function ArbComplex{P}(re::Float64, im::Float64) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_d_d), Cvoid, (Ref{ArbComplex}, Float64, Float64), z, re, im)
    return z
end

function ArbComplex{P}(x::Int64, y::Int64) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_si_si), Cvoid, (Ref{ArbComplex}, Cint, Cint, Clong), z, x, y, P)
    return z
end

function copy(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}), z, x)
    return z
end

deepcopy(x::ArbComplex{P}) where {P} = copy(x)

function ArbComplex{P}(x::ArbComplex{Q}, roundingmode::RoundingMode) where {P,Q}
    z = ArbComplex{P}()
    rounding = match_rounding_mode(roundingmode)
    res = ccall(@libarb(acb_set_round), Cint,
                (Ref{ArbComplex}, Ref{ArbComplex}, Clong, Cint), z, x, P, rounding)
    return z
end

ArbComplex(x::ArbFloat{P}) where {P} = ArbComplex{P}(x)
ArbComplex(x::ArbReal{P}) where {P} = ArbComplex{P}(x)
ArbComplex(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ArbComplex{P}(x,y)
ArbComplex(x::ArbReal{P}, y::ArbReal{P}) where {P} = ArbComplex{P}(x,y)
ArbComplex(x::ArbFloat{P}, y::ArbReal{P}) where {P} = ArbComplex{P}(ArbReal{P}(x),y)
ArbComplex(x::ArbReal{P}, y::ArbFloat{P}) where {P} = ArbComplex{P}(x,ArbReal{P}(y))

function ArbComplex{P}(x::T) where {P,T<:Integer}
    y = ArbReal{P}(x)
    z = ArbComplex{P}(y)
    return z
end

function ArbComplex{P}(x::T1, y::T2) where {P,T1<:Union{Integer,AbstractFloat},T2<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z  = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::Irrational{S}) where {P,S}
    y = ArbReal{P}(x)
    z = ArbComplex{P}(y)
    return z
end

function ArbComplex{P}(x::Irrational{S1}, y::Irrational{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z  = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::Irrational{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::T, y::Irrational{S}) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::T, y::ArbFloat{P}) where {P,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::T, y::ArbReal{P}) where {P,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    z = ArbComplex{P}(x1, y)
    return z
end
function ArbComplex{P}(x::ArbReal{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}}
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x, y1)
    return z
end

function ArbComplex{P}(x::ArbFloat{P}, y::ArbReal{P}) where {P}
    x1 = ArbReal{P}(x)
    z = ArbComplex{P}(x1, y)
    return z
end
function ArbComplex{P}(x::ArbReal{P}, y::ArbFloat{P}) where {P}
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x, y1)
    return z
end

function ArbComplex{P}(x::T, y::ArbFloat{S}) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::T, y::ArbReal{S}) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbReal{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{S1}, y::ArbFloat{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbReal{S1}, y::ArbReal{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{S1}, y::ArbReal{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbReal{S1}, y::ArbFloat{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function midpoint(x::ArbComplex{P}) where {P}
    xreal = real(x)
    ximag = imag(x)
    zreal = ArbReal{P}()
    zimag = ArbReal{P}()
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zreal, xreal)
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zimag, ximag)
    return ArbComplex{P}(zreal, zimag)
end

function radius(x::ArbComplex{P}) where {P}
    xreal = real(x)
    ximag = imag(x)
    zreal = ArbReal{P}()
    zimag = ArbReal{P}()
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zreal, xreal)
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zimag, ximag)
    return ArbComplex{P}(zreal, zimag)
end

trunc(x::ArbComplex{P}) where {P} = ArbComplex{P}(trunc(real(x)), trunc(imag(x)))
trunc(::Type{T}, x::ArbComplex{P}) where {P, T} = T(trunc(real(x))), T(trunc(imag(x)))

floor(x::ArbComplex{P}) where {P} = ArbComplex{P}(florr(real(x)), floor(imag(x)))
floor(::Type{T}, x::ArbComplex{P}) where {P, T} = T(floor(real(x))), T(floor(imag(x)))

ceil(x::ArbComplex{P}) where {P} = ArbComplex{P}(ceil(real(x)), ceil(imag(x)))
ceil(::Type{T}, x::ArbComplex{P}) where {P, T} = T(ceil(real(x))), T(ceil(imag(x)))


text:
# This file defines [`testlearner`](#) that constructs a [`Learner`](#)
# with a very simple optimisation problem. This should be used in tests
# for callbacks.

struct TestModel
    coeff
end
TestModel(coeff::Number) = TestModel([coeff])
Flux.trainable(m::TestModel) = (m.coeff,)
(m::TestModel)(x) = x .* m.coeff
Flux.@functor TestModel

function testbatch(batchsize, coeff)
    xs = rand(1.:100., batchsize)
    return (xs, xs .* coeff)
end

function testbatches(n::Int, coeff, batchsize = 8)
    (testbatch(batchsize, coeff) for _ ∈ 1:n)
end

"""
    testlearner(callbacks...[; opt, nbatches, coeff, batchsize, kwargs...])

Construct a [`Learner`](#) with a simple optimization problem. This
learner should be used in tests that require training a model, e.g.
for callbacks.
"""
function testlearner(
        args...;
        opt = Descent(0.001),
        nbatches = 16,
        coeff = 3,
        batchsize = 8,
        usedefaultcallbacks = false,
        kwargs...)
    model = TestModel(rand())
    data = collect(testbatches(nbatches, coeff, batchsize))
    Learner(
        model,
        (data, data),
        opt,
        Flux.mae,
        args...;
        usedefaultcallbacks = usedefaultcallbacks,
        kwargs...)
end


text:
# randint
n = 1000
x = zeros(Int, n)
for i=1:n
    k = randint([0.1,0.1,0.8])
    x[i] = k
    @test 1 <= k <= 3
end
@test count(k->k==3, x) > n/2

for i=1:n
    k = randint([0.2,0.0,0.8])
    @test 1 <= k <= 3
    @test k != 2
end

for i=1:n
    k = randint([1.0,0.0,0.])
    @test k == 1
end

@test 0 == randint([0.0,0.0,0.])

# DiscreteDistribution
n = 1000
x = zeros(Int, n)
d = DiscreteDistribution([0.1,0.1,0.8])
for i=1:n
    k = rand(d)
    x[i] = k
    @test 1 <= k <= 3
end
@test count(k->k==3, x) > n/2

d = DiscreteDistribution([0.2,0.,0.8])
for i=1:n
    k = rand(d)
    @test 1 <= k <= 3
    @test k != 2
end

d = DiscreteDistribution([1.0,0.0,0.])
for i=1:n
    k = rand(d)
    @test k == 1
end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
using Documenter
using DocThemePC2
using ThreadPinning
using Literate
using LinearAlgebra
# using MKL # optional
# BLAS.set_num_threads(1)

const ci = get(ENV, "CI", "") == "true"

@info "Preparing DocThemePC2"
DocThemePC2.install(@__DIR__)

@info "Generating Documenter.jl site"
makedocs(;
    sitename="ThreadPinning.jl",
    authors="Carsten Bauer",
    modules=[ThreadPinning],
    checkdocs=:exports,
    # doctest = ci,
    pages=[
        "ThreadPinning" => "index.md",
        "Examples" => [
            "Pinning Julia Threads" => "examples/ex_pinning_julia_threads.md",
            "Autochecking BLAS Thread Settings" => "examples/ex_blas.md",
            "Measuring Core-to-Core Latency" => "examples/ex_core2core_latency.md",
        ],
        "Explanations" => [
            "Why Pin Julia Threads?" => "explanations/why.md",
            "Julia Threads + BLAS Threads" => "explanations/blas.md",
            "How It Works" => "explanations/how.md",
        ],
        "References" => [
            "API" => "refs/api.md",
            "Latency" => "refs/latency.md",
            "LibX" => "refs/libX.md",
            "Utility" => "refs/utility.md",
        ],
    ],
    # assets = ["assets/custom.css", "assets/custom.js"]
    repo="https://github.com/carstenbauer/ThreadPinning.jl/blob/{commit}{path}#{line}",
    format=Documenter.HTML(; collapselevel=1),#, assets = ["assets/favicon.ico"])
)

if ci
    @info "Deploying documentation to GitHub"
    deploydocs(;
        repo="github.com/carstenbauer/ThreadPinning.jl.git",
        devbranch="main",
        # push_preview=true,
        # target = "site",
    )
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_fielddata

"""
    vtk_fielddata(;kwargs...)
    vtk_fielddata(children::Any;kwargs...)
    vtk_fielddata(children_maker::Function;kwargs...)

A FieldData component.
FieldData is exposing a FieldData to a downstream element
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
"""
function vtk_fielddata(; kwargs...)
        available_props = Symbol[:children]
        wild_props = Symbol[]
        return Component("vtk_fielddata", "FieldData", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_fielddata(children::Any; kwargs...) = vtk_fielddata(;kwargs..., children = children)
vtk_fielddata(children_maker::Function; kwargs...) = vtk_fielddata(children_maker(); kwargs...)

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
mutable struct Node
    position :: Vector{Float64} # [x, y, z], z is optional
    DOFS :: Vector{Bool} # [dof1, dof2, dof3, ...] type of analysis deduced from length
    x :: Float64 # individual cartesian coordinates of Node
    y :: Float64 #
    z :: Float64 #
    load :: Vector{Float64} # external load, if applicable
    reaction :: Vector{Float64} # external reaction, if applicable
    disp :: Vector{Float64} # displacements of dofs (after analysis)
    elements :: Vector{Tuple{Int64,Int64}} # [(element index, 0 = start / 1 = end), ...]
    globalIndex :: Vector{Int64} #positions in global DOF order

    # create node with assigned DOFs
    function Node(position, DOFS::Vector{Bool})
        node = new(position, DOFS)
        if length(position) == 3
            node.x, node.y, node.z = node.position
        else
            node.x, node.y = node.position
        end
        node.elements = Vector{Tuple{Int64,Int64}}([])
        return node
    end

    # node with predefined DOFs for typical scenarios
    function Node(position, type::Symbol, fixity::Symbol)
        dims = length(position)
        node = new(position, dofMaker(dims, type, fixity))
        if length(position) == 3
            node.x, node.y, node.z = node.position
        else
            node.x, node.y = node.position
        end
        node.elements = Vector{Tuple{Int64,Int64}}([])
        return node
    end

    # basic node generation (position only)
    function Node(position)
        node = new(position)
        if length(position) == 3
            node.x, node.y, node.z = node.position
        else
            node.x, node.y = node.position
        end
        node.elements = Vector{Tuple{Int64,Int64}}([])
        return node
    end
end

mutable struct Element
    nodeIndex :: Vector{Int64} # index of start/end nodes
    type :: Symbol # :truss or :dims
    posStart :: Vector{Float64} # coordinates of start node
    posEnd :: Vector{Float64} # coordinates of end node
    length :: Float64 # element length
    k :: Symmetric{Float64, Matrix{Float64}} # stiffness matrix
    dofIndex ::Vector{Int64} # DOF indices in global expanded coordinates
    localForce :: Vector{Float64}  #list of internal loads in local coordinates
    axialForce :: Float64 #-ve compression, +ve tension
    globalForce :: Vector{Float64}
    A :: Float64 # section area
    E :: Float64 # modulus of elasticity
    G :: Float64 # shear modulus, OPTIONAL
    Iz :: Float64 # strong axis moment of inertia, OPTIONAL
    Iy :: Float64 # weak axis moment of inertia, OPTIONAL
    J :: Float64 # torsional constant, OPTIONAL
    Ψ :: Float64 # angle of roll w/r/t local x axis, for 3d frames
    R :: Matrix{Float64} # rotation matrix
    LCS :: Vector{Vec3{Float64}} #local coordinate system

    # Basic generator without material properties
    function Element(nodes::Vector{Node}, nodeIndex::Vector{Int64}, type::Symbol)
        element = new(nodeIndex)
        element.posStart = nodes[nodeIndex[1]].position
        element.posEnd = nodes[nodeIndex[2]].position
        element.length = norm(element.posEnd .- element.posStart)
        element.type = type
        return element
    end

    # truss generator, 2 and 3 dimensions
    function Element(nodes::Vector{Node}, nodeIndex::Vector{Int64}, E::Float64, A::Float64)
        element = new(nodeIndex)
        element.posStart = nodes[nodeIndex[1]].position
        element.posEnd = nodes[nodeIndex[2]].position
        element.length = norm(element.posEnd .- element.posStart)
        element.E = E
        element.A = A
        element.type = :truss
        return element
    end

    # 2D frame generator
    function Element(nodes::Vector{Node}, nodeIndex::Vector{Int64}, E::Float64, A::Float64, I::Float64)
        element = new(nodeIndex)
        element.posStart = nodes[nodeIndex[1]].position
        element.posEnd = nodes[nodeIndex[2]].position
        element.length = norm(element.posEnd .- element.posStart)
        element.E = E
        element.A = A
        element.Iz = I
        element.type = :frame
        return element
    end

    # 3D frame generator
    function Element(nodes::Vector{Node}, nodeIndex::Vector{Int64}, E::Float64, A::Float64, G::Float64, Iz::Float64, Iy::Float64, J::Float64)
        element = new(nodeIndex)
        element.posStart = nodes[nodeIndex[1]].position
        element.posEnd = nodes[nodeIndex[2]].position
        element.length = norm(element.posEnd .- element.posStart)
        element.E = E
        element.A = A
        element.G = G
        element.Iz = Iz
        element.Iy = Iy
        element.J = J
        element.Ψ = pi/2 #default
        element.type = :frame
        element.LCS = lcs(element, element.Ψ)
        return element
    end

end

# creates a Load type with proper DOF indexing
mutable struct Load
    position :: Vector{Float64} #approximate position of load
    load :: Vector{Float64} #value of load [p1, p2, p3, ..., pdof]
    index :: Int64 #index of node assigned to load

    function Load(nodes::Vector{Node}, position::Vector{Float64}, load::Vector{Float64}; tol = 0.5)
        load = new(position, load)
        load.index = findmin([norm(n.position .- load.position) for n in nodes])[2]
        return load
    end
end

mutable struct Structure
    nodes :: Vector{Node} # array of Nodes
    elements :: Vector{Element} # array of Elements
    loads :: Vector{Load} # array of Loads
    dims :: Int64 # 2, 3
    displacedNodes :: Vector{Node} # array of displaced nodes
    displacedElements :: Vector{Element} # array of displaced elements
    DOFS :: Vector{Bool} # degrees of freedom
    K :: SparseMatrixCSC{Float64,Int64} # global stiffness matrix
    F :: Vector{Float64} # external loads
    U :: Vector{Float64} # Displacements after analysis
    reactions :: Vector{Float64} # external reactions
    compliance :: Float64 # compliance after analysis
    nNodes :: Int64 # number of nodes
    nElements :: Int64 #number of elements
    nDOFS :: Int64 #number of DOFs (total)
    freeDOFS :: Vector{Int64} #array of free DOFs

    #general structure definition with all necesary components
    function Structure(nodes::Vector{Node}, elements::Vector{Element}, loads::Vector{Load})

        # initialize
        structure = new(nodes, elements, loads)

        # deduce dimensions
        structure.dims = length(nodes[1].position)

        # global DOF vector
        structure.DOFS = vcat([node.DOFS for node in nodes]...)

        # general information
        structure.nNodes = length(nodes)
        structure.nElements = length(elements)
        structure.nDOFS = length(structure.DOFS)
        structure.freeDOFS = findall(structure.DOFS)

        # create global load vector
        structure.F = zeros(structure.nDOFS)
        nodalDOFLength = length(nodes[1].DOFS)

        for load in loads
            idx = load.index * nodalDOFLength - (nodalDOFLength - 1) .+ collect(0:nodalDOFLength-1)
            structure.F[idx] .= load.load
        end

        return structure
    end

    # Basic structure definition without loads
    function Structure(nodes::Vector{Node}, elements::Vector{Element})
        structure = new(nodes, elements)
        structure.dims = length(nodes[1].position)
        structure.DOFS = vcat([node.DOFS for node in nodes]...)

        structure.nNodes = length(nodes)
        structure.nElements = length(elements)
        structure.nDOFS = length(structure.DOFS)
        structure.freeDOFS = findall(structure.DOFS)

        return structure
    end

end

mutable struct Geometry
    structure :: Structure
    nodes :: Vector #{Point3{Float64}}
    elements :: Vector{Vector} #Vector{Vector{Point3{Float64}}}
    axialForce :: Vector{Float64}
    areas :: Vector{Float64}
    displacedNodes :: Vector #Vector{Point3{Float64}}
    displacedElements :: Vector{Vector} #Vector{Vector{Point3{Float64}}}
    midpoints :: Vector #Vector{Point3{Float64}}
    loads :: Vector{Vector{Float64}}
    reactions :: Vector #Vector{Vec3{Float64}}
    nodeLabels :: Vector{String}
    elementLabels :: Vector{String}
    pinDOFS :: NTuple #{6, Vector{Float64}}
    fixDOFS :: NTuple #{6, Vector{Float64}}

    function Geometry(structure::Structure; nodeAnalysis = false, SF = :auto)
        geometry = new(structure)

        if structure.dims == 3
            geometry.nodes = Point3.([node.position for node in structure.nodes])
            geometry.elements = [[Point3(element.posStart), Point3(element.posEnd)] for element in structure.elements]
        else
            geometry.nodes = Point3.([[node.position; 0.0] for node in structure.nodes])
            geometry.elements = [[Point3([element.posStart; 0.0]), Point3([element.posEnd; 0.0])] for element in structure.elements]
        end

        areas = [element.A for element in structure.elements]
        geometry.areas = areas ./ maximum(areas)
        # geometry.midpoints = [(e.posStart + e.posEnd) / 2 for e in structure.elements]

        geometry.midpoints = sum.(geometry.elements) ./ 2

        geometry.nodeLabels = ["N" * string(i) for i = 1:structure.nNodes]
        geometry.elementLabels = ["E" * string(i) for i = 1:structure.nElements]

        geometry.pinDOFS = pinDofConverter(structure)
        geometry.fixDOFS = fixDofConverter(structure)

        if isdefined(structure, :displacedNodes)
            if structure.dims == 3
                geometry.displacedNodes = Point3.([node.position for node in structure.displacedNodes])
                geometry.displacedElements = [[Point3(element.posStart), Point3(element.posEnd)] for element in structure.displacedElements]
            else
                geometry.displacedNodes = Point3.([[node.position; 0.0] for node in structure.displacedNodes])
                geometry.displacedElements = [[Point3([element.posStart; 0.0]), Point3([element.posEnd; 0.0])] for element in structure.displacedElements]
            end
        end

        if isdefined(structure, :loads)
            geometry.loads = loadConverter(structure; scaleFactor = SF)
        end

        if isdefined(structure.elements[1], :localForce)
            geometry.axialForce = [element.axialForce for element in structure.elements]
            geometry.reactions = reactionConverter(structure)
        end

        return geometry
    end

end

max_stars_repo_name: keithjlee/Asap

max_issues_repo_name: keithjlee/Asap

max_forks_repo_name: keithjlee/Asap


text:
using Gadfly
using Main.DistMesh2D

function fh(x, y)::Real
    return 1.0
end

function dc(p)
    return -1 + abs(4 - sqrt(sum(p .^ 2)))
end

function dca(p)
    return sqrt(sum(p .^ 2)) - 2
end

function dcrcc(p)
    return max(drectangle(p, -100.0, 100.0, -100.0, 100.0), - dcircle(p, 20.0, 10.0, 15.0))
end

function eplot(x,y)
    Gadfly.plot(x = x, y = y, Geom.path, Coord.cartesian(fixed = true))
end

#x, y = distmesh2d(dcrcc, fh, [-100.0 -100.0; 100.0 100.0], 7.0, [-100.0 -100.0; -100.0 100.0; 100.0 -100.0; 100.0 100.0])
function fdd(p)
    return dcircle(p, 0.0, 0.0, 1.0)
end
x, y = distmesh2d(fdd, fh, [-1.0 -1.0; 1.0 1.0], 0.2)
eplot(x,y)

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
module MyFunctions

export fact

function fact(n::Integer)
    if n < 0
        throw(DomainError(n, "fact argument must be positive"))
    else
        result = 1
        for i in 2:n
            result *= i
        end
        return result
    end
end

end


text:
"""
    equiSpacedPosition(x, knotsIn)
Returns the bin in which x is contained for the bins defined by the
limits in knotsIn. By default if the value is in the limit it is pPosition
at the right bin except for the las limit which is positioned to the left.
#Examples
```julia julia-repl
julia> equiSpacedPosition(2, [1,2,3,4,5,6,7,8])
2
```
"""
function equiSpacedPosition(x, knotsIn)

    @assert all(diff(knotsIn).>=0)    # bin limits are monotonically increaseing

    nbins = length(knotsIn) - 1
    lsubd = knotsIn[2] - knotsIn[1]

    proximityTol = 1e-4
    xinDist = abs(knotsIn[1]-x)
    xfinDist = abs(x - knotsIn[end])

    if x <= knotsIn[1]

        xinDist/lsubd < proximityTol ? s = 1 : throw(DomainError(x, "The value is out of the domain (left)"))

    elseif x >= knotsIn[end]

        xfinDist/lsubd < proximityTol ? s = nbins : throw(DomainError(x, "The value is out of the domain (right)"))

    else
        s = floor(Int,(x - knotsIn[1])/lsubd) + 1

        s > nbins ? s = nbins : nothing
    end
    return s
end

"""
    evalBSplineBasis(x, knots)
Builds the evaluation of the cubic BSpline basis corresponding to
the `knots` at `x`.
#Examples
```julia julia-repl
julia> evalBSplineBasis(10.,buildKnots(1.,10., 12))
12-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.16666666666666663
 0.6666666666666666
 0.16666666666666666
```
"""
function evalBSplineBasis(x, knots)
    transfMatrix = (1. / 6.) * [ -1. 3. -3. 1.;
                                  3. -6. 3. 0.;
                                 -3. 0. 3. 0.;
                                  1. 4. 1. 0.]
    nvertices = length(knots)
    lsubd = knots[2] - knots[1]
    knotsIn = knots[2:end-1] # the internal vertices are all except extrema

    bsplineBasis = zeros(promote_type(typeof(x), eltype(knots)), nvertices)

    s = equiSpacedPosition(x, knotsIn)

    ξ = (x - knotsIn[s])/lsubd
    ξvlocal = [ξ^3, ξ^2, ξ, one(promote_type(typeof(x), eltype(knots)))]

    bsplineBasis[s:s+3] = ξvlocal'*transfMatrix

    return bsplineBasis
end

"""
    evalBSplineBasis2dim(x, y, knotsx, knotsy)
Builds the evaluation of the dyadic product of 2 cubic BSplines
basis, the first corresponding to `knotsx` at `x` and the second
corresponding to `knotsy`at `y`. The output is already vectorized
instead of being a matrix.
#Examples
```julia julia-repl
julia> evalBSplineBasis2dim(3.,2.,buildKnots(1.,10., 12),buildKnots(0.,5.,14))
168-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
```
"""
function evalBSplineBasis2dim(x, y, knotsx, knotsy)

    bsplineBasisx = evalBSplineBasis(x, knotsx)
    bsplineBasisy = evalBSplineBasis(y, knotsy)
    return vec(bsplineBasisx*bsplineBasisy')
end

"""
    evalBSplineBasisD1(x, knots)
Builds the evaluation of the derivative of the cubic BSpline basis corresponding to
the `knots` at `x`.
#Examples
```julia julia-repl
julia> evalBSplineBasisD1(10.,buildKnots(1.,10., 12))
12-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  ⋮
 -0.5
  0.0
  0.5
```
"""
function evalBSplineBasisD1(x, knots)

    transfMatrix = (1. / 6.) * [ -1. 3. -3. 1.;
                                   3. -6. 3. 0.;
                                   -3. 0. 3. 0.;
                                    1. 4. 1. 0.]
    nvertices = length(knots)
    lsubd = knots[2] - knots[1]
    knotsIn = knots[2:end-1] # the internal vertices are all except extrema
    bsplineBasisD1 = zeros(typeof(x), nvertices)

    s = equiSpacedPosition(x, knotsIn)

    ξ = (x - knotsIn[s])/lsubd
    ξvlocalD1 = [3. * ξ^2, 2. * ξ, 1., 0.]

    bsplineBasisD1[s:s+3] = ξvlocalD1'*transfMatrix*(1/lsubd)
    return bsplineBasisD1
end

function evalBSplineBasis2dimD1x(x, y, knotsx, knotsy)
    dbsplineBasisx = evalBSplineBasisD1(x, knotsx)
    bsplineBasisy = evalBSplineBasis(y, knotsy)
    return vec(dbsplineBasisx*bsplineBasisy')
end

function evalBSplineBasis2dimD1y(x, y, knotsx, knotsy)
    bsplineBasisx = evalBSplineBasis(x, knotsx)
    dbsplineBasisy = evalBSplineBasisD1(y, knotsy)
    return vec(bsplineBasisx*dbsplineBasisy')
end

"""
    evalBSpline(x, knots, vertices)
Evaluation of the cubic BSpline at `x` for vertices values `vertices`
and the BSpline basis defined by `evalBSplineBasis(x, knots)`.
#Examples
```julia julia-repl
julia> evalBSpline(9.,buildKnots(1.,10., 12), ones(12))
0.9999999999999999
```
"""
function evalBSpline(x, knots, vertices)
    return evalBSplineBasis(x, knots)'*vertices
end

evalBSpline(x, bs::BSpline) = evalBSpline(x, bs.knots, bs.vertices)

"""
    evalBSpline2dim(x, y, knotsx, knotsy, vertices)
Evaluation of the cubic bidimensional BSpline with knots `knotsx`, `knotsy`
at `(x,y)` for vertices `vertices`.

#Examples
```julia julia-repl
julia> knotsx = knotsy = 0:0.1:10;
julia> vertices = rand(length(knotsx)*length(knotsy ));
julia> evalBSpline2dim(2., 3., knotsx, knotsy, vertices)
0.3800038047730809
```
"""
function evalBSpline2dim(x, y, knotsx, knotsy, vertices)
    return evalBSplineBasis2dim(x, y, knotsx, knotsy)'*vertices
end

evalBSpline2dim(x, y, bs2dim::BSpline2dim) = evalBSpline2dim(x, y, bs2dim.knotsx, bs2dim.knotsy, bs2dim.vertices)

"""
    evalBSplineD1(x, knots, vertices)
Evaluating the derivative ofthe cubic BSpline at `x` for vertices
values `vertices` and the BSpline basis derivative defined by
`evalBSplineBasisD1(x, knots)`.
#Examples
```julia julia-repl
julia> evalBSplineD1(9.,buildKnots(1.,10., 12), ones(12))
0.0
```
"""
function evalBSplineD1(x, knots, vertices)
    return evalBSplineBasisD1(x, knots)'*vertices
end

evalBSplineD1(x, bs::BSpline) = evalBSplineD1(x, bs.knots, bs.vertices)

"""
    evalBSpline2dimD1x(x, y, knotsx, knotsy, vertices)
Evaluation of the derivative with respect to `x`of the
cubic bidimensional BSpline with knots `knotsx`, `knotsy`
at `(x,y)` for vertices `vertices`.

#Examples
```julia julia-repl
julia> knotsx = knotsy = 0:0.1:10;
julia> vertices = rand(length(knotsx)*length(knotsy ));
julia> evalBSpline2dimD1x(2., 3., knotsx, knotsy, vertices)
-2.0863387254122316
```
"""
function evalBSpline2dimD1x(x, y, knotsx, knotsy, vertices)
    return evalBSplineBasis2dimD1x(x, y, knotsx, knotsy)'*vertices
end

evalBSpline2dimD1x(x, y, bs2dim::BSpline2dim) = evalBSpline2dimD1x(x, y, bs2dim.knotsx, bs2dim.knotsy, bs2dim.vertices)

"""
    evalBSpline2dimD1y(x, y, knotsx, knotsy, vertices)
Evaluation of the derivative with respect to `y`of the
cubic bidimensional BSpline with knots `knotsx`, `knotsy`
at `(x,y)` for vertices `vertices`.

#Examples
```julia julia-repl
julia> knotsx = knotsy = 0:0.1:10;
julia> vertices = rand(length(knotsx)*length(knotsy ));
julia> evalBSpline2dimD1y(2., 3., knotsx, knotsy, vertices)
2.345106484778582
```
"""
function evalBSpline2dimD1y(x, y, knotsx, knotsy, vertices)
    return evalBSplineBasis2dimD1y(x, y, knotsx, knotsy)'*vertices
end

evalBSpline2dimD1y(x, y, bs2dim::BSpline2dim) = evalBSpline2dimD1y(x, y, bs2dim.knotsx, bs2dim.knotsy, bs2dim.vertices)

"""
    nVerticesNeededLeftRight(knots, xin, xend)
Given the knots vector `knots` and the limits `xin` and `xend`,
determines the number of vertices that knots needs  on the left
and on the right needed to reach the values `xin` and `xend`.

#Examples
```julia julia-repl
julia> knots = 0.:1.:10.;
julia> nVerticesNeededLeftRight(knots, -2, 14)
(3, 5)
```
"""
function nVerticesNeededLeftRight(knots, xin, xend)
    knotsInmin = knots[2]
    knotsInmax = knots[end-1]

    lsubd = knots[2] - knots[1]

    intervLeft = (knotsInmin - xin)/lsubd
    intervRight = (xend - knotsInmax)/lsubd

    intervLeft > 0 ? nKnotsLeft = Int64(ceil(intervLeft)) : nKnotsLeft = 0
    intervRight > 0 ? nKnotsRight = Int64(ceil(intervRight)) : nKnotsRight = 0
    return nKnotsLeft, nKnotsRight
end

"""
    extendedKnots(knots, nKnotsLeft, nKnotsRight)
Given the knots vector `knots` and `nKnotsLeft` and `nKnotsRight`
extends the `knots` on the left `nKnotsLeft` and on the right `nKnotsRight`.

#Examples
```julia julia-repl
julia> knots = 0.:1.:10.;
julia> extendedKnots(knots, 2, 2)
15-element Vector{Float64}:
 -2.0
 -1.0
  ⋮
 11.0
 12.0
```
"""
function extendedKnots(knots, nKnotsLeft, nKnotsRight)

    lsubd = knots[2] - knots[1]

    addKnotsLeft = range(knots[1]-lsubd*nKnotsLeft, knots[1]-lsubd, length = nKnotsLeft)
    addKnotsRight = range(knots[end]+lsubd, knots[end]+lsubd*nKnotsRight, length = nKnotsRight)

    knotsNew = vcat(addKnotsLeft, knots, addKnotsRight)

    return knotsNew
end

"""
    expandBSpline(knots, vertices, xin, xend)
Given the knots and the vertices of a given BSpline with a certain range of validity,
it returns a new BSpline which extends the range of validity of the previous BSpline
betwenn `xin` and `xend`.
#Examples
```julia julia-repl
julia> bs2 = expandBSpline(bs1, 0., 5.);
julia> bs2.knots
20-element Vector{Float64}:
 -0.6666666666666663
 -0.33333333333333304
  2.220446049250313e-16
  ⋮
  5.333333333333333
  5.666666666666666
julia> bs2.vertices
20-element Vector{Float64}:
  -1.8150147193435557
  -1.2594147415453065
  -0.7038147637472744
   ⋮
  27.07391859179092
  29.851651947327774
```
"""
function expandBSpline(knots, vertices,
                       xin, xend)

    nvertices = length(knots)

    nKnotsLeft, nKnotsRight = nVerticesNeededLeftRight(knots, xin, xend)
    knotsNew = extendedKnots(knots, nKnotsLeft, nKnotsRight)
    nverticesNew = length(knotsNew)

    A = zeros(nverticesNew, nverticesNew)

    D2 = D2matrix(nverticesNew)
    Ω2 = Ω2matrix(nverticesNew)
    W = Wmatrix(nverticesNew)

    A[nKnotsLeft+1 : (nKnotsLeft + nvertices), nKnotsLeft+1 : (nKnotsLeft + nvertices)] =Float64.(I(nvertices))

    for i = nKnotsLeft+1:(nKnotsLeft + nvertices)
        W[i,i] = 1e4
    end

    Ω2[diagind(Ω2)] .= 1.

    Asys = A'*W*A + D2'*Ω2*D2
    bsys = A'*W*vcat(zeros(nKnotsLeft), vertices, zeros(nKnotsRight))

    verticesNew = Asys\bsys

    return BSpline(knotsNew, verticesNew)
end

expandBSpline(bs::BSpline, xin, xend) = expandBSpline(bs.knots, bs.vertices, xin, xend)

"""
    expandBSpline2dim(knotsx, knotsy, vertices, xin, xend, yin, yend)
Given a 2D BSpline by `knotsx` , `knotsy` and `vetices`
compute the expanded 2D BSpline to the limits `xin, `xend`, `yin`, `yend`.

#Examples
```julia julia-repl
julia> knotsx = knotsy = 0:0.5:10;
julia> vertices = rand(length(knotsx)*length(knotsy ));
julia> expandBSpline2dim(knotsx, knotsy, vertices, -1., 11., -1., 11.)
BSpline2dim{Float64}([-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0  …  7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5],
[-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0  …  7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5],
[-1.0029300060952653, -0.4790706367303448, 0.03153695697931029, 0.48370119176194715, 0.8342290770083228, 1.1030443114892081, 1.2623356760827873, 1.2585422976688359, 1.131589871095645, 0.9860138921208538
…  0.8984551404840291, 0.7049426036058777, 0.4552726974231055, 0.2849091663041328, 0.22446630431980158, 0.22686278460500844, 0.3058408589841374, 0.4651356516733504, 0.6732304620674567, 0.8966217568026448])
```
"""
function expandBSpline2dim(knotsx, knotsy, vertices,
                           xin, xend, yin, yend)
    nverticesx = length(knotsx)
    nverticesy = length(knotsy)

    nKnotsLeftx, nKnotsRightx = nVerticesNeededLeftRight(knotsx, xin, xend)
    knotsNewx = extendedKnots(knotsx, nKnotsLeftx, nKnotsRightx)

    nverticesxNew = length(knotsNewx)

    nKnotsLefty, nKnotsRighty = nVerticesNeededLeftRight(knotsy, yin, yend)
    knotsNewy = extendedKnots(knotsy, nKnotsLefty, nKnotsRighty)

    nverticesyNew = length(knotsNewy)

    indicesVector = reshape(1:nverticesxNew*nverticesyNew, nverticesxNew, nverticesyNew)

    yindVerticesNew = zeros(nverticesxNew*nverticesyNew)
    NNew = zeros(nverticesxNew*nverticesyNew, nverticesxNew*nverticesyNew)

    # Identify a priori known vertices

    verticesMatrix = reshape(vertices, nverticesx, nverticesy)

    w = 1e10

    for i = (nKnotsLeftx+1):(nverticesx+nKnotsLeftx)
        for j = (nKnotsLefty+1):(nverticesy+nKnotsLefty)
            locali = i - nKnotsLeftx  ;localj = j - nKnotsLefty
            yindVerticesNew[indicesVector[i,j]] = verticesMatrix[locali,localj]
            NNew[indicesVector[i,j],indicesVector[i,j]] = 1
        end
    end

    Wmat = Diagonal(fill(w, size(NNew,1)))

    Px, Py = PD2matrices2dim(nverticesxNew, nverticesyNew)
    λx = 1.
    λy = 1.

    Asys = NNew'*Wmat*NNew + λx*Px + λy*Py
    bsys = NNew'*Wmat*yindVerticesNew

    verticesNew =  Asys\bsys

    return BSpline2dim(knotsNewx, knotsNewy, verticesNew)
end

expandBSpline2dim(bs2dim::BSpline2dim, xin, xend, yin, yend) = expandBSpline2dim(bs2dim.knotsx, bs2dim.knotsy, bs2dim.vertices, xin, xend, yin, yend)

## Document
function buildSystemMatrix(N;
                            W = Wmatrix(size(N,1)),
                           Ω1 = Ω1matrix(size(N,2)),
                           Ω2 = Ω2matrix(size(N,2)),
                           Ω3 = Ω3matrix(size(N,2)))
    nvertices = size(N,2)
    D1 = D1matrix(nvertices)
    D2 = D2matrix(nvertices)
    D3 = D3matrix(nvertices)

    Asys = N'*W*N + D1'*Ω1*D1 + D2'*Ω2*D2 + D3'*Ω3*D3

    return Asys
end
## Document
function buildSystemIndVect(N, y;
                            W = Wmatrix(size(N,1)))
    bsys = N'*W*y
    return bsys
end