text:
"""
    VaR(states,prob,alpha::Float32)

implements the Value-at-Risk at level ``\\alpha`` defined by
```math
VaR_\\alpha (Y) = \\arg \\min_x \\left( x\\in \\mathbb{R} : F_Y(x) \\geq \\alpha \\right),
```
for the random variable ``Y`` defined by `states` and `prob`.
"""
function VaR(x::Vector{Float64}, f::Vector{Float64}, α::Float32)
    ind = sortperm(x[1:length(x)])
    x = x[ind]; f = f[ind];
    i = findfirst(p -> p≥α, cumsum(f))
    if i === nothing
        return x[end]
    else
        return x[i]
    end
end

"""
    CTE(states,prob,alpha::Float32)

implements the Conditional Value-at-Risk at level ``\\alpha`` defined by
```math
CTE_\\alpha (Y) = VaR_\\alpha(Y) + \\frac{1}{1-\\alpha} \\mathbb{E} \\left( Y- VaR_\\alpha (Y) \\right)_+ ,
```
for the random variable ``Y`` defined by `states` and `prob`.
"""
function CTE(x::Vector{Float64}, f::Vector{Float64}, α::Float32)
    x = -x; α = 1-α;
    x_α = VaR(x, f, α)
    if iszero(α)
        return -x_α
    else
        tail = x .≤ x_α
        result = (sum(x[tail] .* f[tail]) - (sum(f[tail]) - α) * x_α) / α
        return -result
    end
end

"""
    EVaR2(states,prob,beta)

Solves the optimization problem associated with the primal formulation of the Entropic Value-at-Risk:

```math
EVaR_\\alpha(Y) = \\min_{x >0} \\frac{1}{x} \\left( \\beta +  \\log\\mathbb{E} e^{xY} \\right),
```
where ``Y`` is the discrete random variable defined by `states` and `prob`.
Here the optimization is done via the goldenSearch optimization routine implemented as part of this package.
"""
function EVaR2(states::Vector{Float64}, prob::Vector{Float64},beta::Float32)
	if sum(prob) == 0
		prob = ones(length(states))./ length(states)
	end
	dualZ = zeros(length(states))
	esssup = maximum(states[vec(prob.>0.0)])
	varm = dot(prob, (states .- dot(prob,states)).^2)
	if beta == 0.0
		Fvalue = dot(prob, states)
		tOpt = 0.0
		dualZ = exp.(tOpt .*states)
	else
		if varm >0.0
			tOpt = sqrt.(2 *beta/varm)
			function objective1(t)
				t=max(t,0.0)
				w = states./t
				offset = maximum(w)
    			we = exp.(w .- offset)
    			s = dot(prob, we)
    			w = log(s) + offset
				tmp = t*beta + t*w
				#tmp = t*beta + t*log(dot(prob,exp.(states./t)))
				return(tmp)
			end
			out = goldenSearch(objective1,tOpt)
			tOpt = out[1]
			Fvalue = out[2]
		else
			tOpt = 0.0
			Fvalue = 0.0
			dualZ = esssup
		end
		#if Fvalue < esssup
			dualZ = exp.(tOpt .*states)
		#else
		#	Fvalue = esssup
		#	dualZ[states .>= esssup] .= 1
		#end
		dualZ = dualZ / dot(prob,dualZ)
	end
	return(Fvalue,tOpt,dualZ) # for numerical stability
end

"""
    EVaR(states,prob,beta::Float32)

Solves the optimization problem associated with the primal formulation of the Entropic Value-at-Risk:

```math
EVaR_\\alpha(Y) = \\min_{x >0} \\frac{1}{x} \\left( \\beta +  \\log\\mathbb{E} e^{xY} \\right),
```
where ``Y`` is the discrete random variable defined by `states` and `prob`. Here, the optimization is done using JuMP and Ipopt.
"""
function EVaR(states::Vector{Float64}, prob::Vector{Float64},beta::Float32)
	if sum(prob) == 0
		prob = ones(length(states))./ length(states)
	end
	dualZ = zeros(length(states))
	esssup = maximum(states[prob.>0.0]) #error
	var = dot(prob, (states .- dot(prob,states)).^2 )
	if beta == 0.0
		EVaR = dot(prob, states)
		tOpt = 0
		dualZ = exp.(tOpt .*states)
	else
		if var >0
			tOpt = sqrt.(2 *beta/var)
			EV = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" =>0))
			function objective1(t)
				return( (beta.+ log(dot(prob,exp.(t.*states))))/t)
			end
			register(EV, :objective1, 1, objective1, autodiff=true)

	    	@variable(EV,t,start = tOpt)
			@constraint(EV, t .>=0)
	    	@NLobjective(EV, Min, objective1(t) )

			JuMP.optimize!(EV)
			EVaR = JuMP.objective_value(EV)
	    	tOpt = JuMP.value(t)
		else
			tOpt = 0
			EVaR = 0
			dualZ = esssup
		end
		if EVaR < esssup
			dualZ = exp.(tOpt .*states)
		else
			EVaR = esssup
			#dualZ[states==esssup] .= 1
		end
		#dualZ = dualZ / dot(prob,dualZ)
	end
	return(EVaR,tOpt,dualZ)
end

"""
    AVaR(states,prob,alpha::Float32)

Solves the optimization problem associated with the primal formulation of the Average Value-at-Risk:

```math
AVaR_\\alpha(Y) = \\min_{x\\in \\mathbb{R}} x + \\frac{1}{1-\\alpha} \\mathbb{E} \\left( Y - x \\right)_+,
```
where ``Y`` is the discrete random variable defined by `states` and `prob`.
"""
function AVaR(states::Vector{Float64},prob::Vector{Float64}, alpha::Float32)
    var2 = dot(prob, (states .- dot(prob,states)).^2);
     #First trivial case
    if alpha == 0.0
        return(dot(prob, states))
  #Now set up optimization problem using goldenSearch
    else
       tOpt = sqrt.(2 *alpha/var2)
       EV = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" =>0))
		function objective1(t)
			return( t + 1/(1-alpha) * dot(prob, max.(states .- t,0)) )
		end
		register(EV, :objective1, 1, objective1, autodiff=true)
	    @variable(EV,t,start = tOpt)
		@NLobjective(EV, Min, objective1(t) )

		JuMP.optimize!(EV)
		AVaR = JuMP.objective_value(EV)
	    tOpt = JuMP.value(t)
    end
    return(AVaR, tOpt)
end


text:
using LAP_julia: load_H, make_flow_from_H, gen_rand_points

## create a birl like structure for with the oxford affine dataset
base_path = "/Users/MrTrololord/Downloads/ox_affine/"
images_path = joinpath(base_path, "images")
folders = filter(x -> !isfile(x), readdir(images_path))

#converted to png
for folder in folders
    my_path = joinpath(base_path, folder)
    for file in filter(x -> occursin("pgm", x), readdir(my_path))
        @info file
        img_path = joinpath(my_path, file)
        file_name = split(file, ".")[1]
        out_path = joinpath(my_path, file_name*".png")
        run(`convert $img_path $out_path`)
    end
end

# all pngs
image_paths = []
for folder in folders
    my_path = joinpath(images_path, folder)
    @info my_path
    for file in filter(x -> occursin("png", x), readdir(my_path))
        push!(image_paths, joinpath(my_path, file))
    end
end

# delete png base path
deleting = "/Users/MrTrololord/Downloads/"
image_paths_relative = map(x -> replace(x, deleting => ""), image_paths)

# get target images --> img1
target_img_paths = []
for folder in folders
    my_path = joinpath(images_path, folder)
    for file in filter(x -> occursin("img1.png", x), readdir(my_path))
        push!(target_img_paths, joinpath(my_path, file))
    end
end

target_img_paths

# paths of target_source_transform triads
target_source_transform_paths = []
for folder in folders
    my_path = joinpath(images_path, folder)
    for target in filter(x -> occursin("img1.png", x), readdir(my_path))
        target_path = joinpath(my_path, target)
        for source in filter(x -> !occursin("img1.png", x) && occursin("png", x), readdir(my_path))
            source_path = joinpath(my_path, source)
            img_number = parse(Int, source[end-4])
            tranform_file = filter(x -> !occursin("img", x) && occursin(string(img_number), x), readdir(my_path))[1]
            tranform_path = joinpath(my_path, tranform_file)
            push!(target_source_transform_paths, [target_path, source_path, tranform_path])
        end
    end
end
target_source_transform_paths = permutedims(hcat(target_source_transform_paths...))

target_source_transform_paths_relative = map(x -> replace(x, deleting => ""), target_source_transform_paths)

using DataFrames, CSV
COUNT = 500

# filter target landmarks not it source img
# create target_landmarks
for folder in folders
    target_img = load(joinpath(base_path, "images", folder, "img1.png"))
    target_lnd_path = joinpath(base_path, "landmarks", folder, "target.csv")
    target_inds = gen_rand_points(target_img, COUNT, "Gridded")
    for k in 2:6
        # get flow
        transform_path = joinpath(base_path, "images", folder, "H1to$(k)p")
        cur_H = LAP_julia.load_H(transform_path)
        cur_flow = LAP_julia.make_flow_from_H(cur_H, size(target_img))
        # shift and filter landmarks
        target_inds = LAP_julia.get_valid_landmarks(cur_flow.*(-1), target_inds)
    end
    @info size(target_inds)
    @info size(target_img)
    LAP_julia.save_landmarks(target_inds, target_lnd_path)
end

# warp target landmarks to source. and create
for (k, (target_path, source_path, transform_path)) in enumerate(eachrow(target_source_transform_paths))
    @info (k, target_path, source_path, transform_path)
    folder = split(target_path, "/")[7]
    target_lnd_path = joinpath(base_path, "landmarks", folder, "target.csv")
    cur_H = LAP_julia.load_H(transform_path)
    target_img = load(target_path)
    cur_flow = LAP_julia.make_flow_from_H(cur_H, size(target_img))
    @info size(cur_flow)
    img_number = transform_path[end-1]
    save_path = joinpath(base_path, "landmarks", folder, "source$img_number.csv")
    @info save_path
    LAP_julia.save_shift_landmarks(target_lnd_path, cur_flow.*(-1), save_path)
end

# create loc table
tmp_df = CSV.read("/Volumes/davidkunz/Documents/bakalarka/birl_small/dataset_small_repaired.csv") |> DataFrame

df = deepcopy(tmp_df)
delete!(df, 1:size(tmp_df,1))

rename!(df, Symbol.(replace.(string.(names(df)), Ref(" "=> "_"))))
rename!(df, Symbol.(names(df)))

eltypes(df)
abc = [:a , :b, :c, :d, :e, :f, :g, :h, :i, :j, :k]
new_df = DataFrame(a = eltypes(df)[1][],
                   b = eltypes(df)[2][],
                   c = eltypes(df)[3][],
                   d = eltypes(df)[4][],
                   e = eltypes(df)[5][],
                   f = eltypes(df)[6][],
                   g = eltypes(df)[7][],
                   h = eltypes(df)[8][],
                   i = eltypes(df)[9][],
                   j = eltypes(df)[10][],
                   k = eltypes(df)[11][])

function rename_dict_keys(from_name, to_name, dict)
    new_dict = Dict()
    rename_dict = Dict(from_name .=> to_name)
    for (key, value) in rename_dict
        new_dict[value] = dict[key]
    end
    return new_dict
end

i = 0
for folder in folders
    for k in 2:6
        global i = i + 1
        row = Dict{String}{Any}()
        row["Column1"] = i
        # images
        row["Target image"] = joinpath("images", folder, "img1.png")
        row["Source image"] = joinpath("images", folder, "img$k.png")
        # lnds
        row["Target landmarks"] = joinpath("landmarks", folder, "target.csv")
        row["Source landmarks"] = joinpath("landmarks", folder, "source$k.csv")
        # sizes
        source_img = load(joinpath(base_path, "images", folder, "img$k.png"))
        row["Image size [pixels]"] = string(size(source_img))
        row["Image diagonal [pixels]"] = round(sqrt(sum(size(source_img).^2)), digits=1)
        row["status"] = "training"
        # empty values
        row["Warped target landmarks"] = missing
        row["Warped source landmarks"] = missing
        row["Execution time [minutes]"] = missing
        global row

        @info rename_dict_keys(names(tmp_df), abc, row)
        push!(new_df, rename_dict_keys(names(tmp_df), abc, row))
    end
end

new_df

# rename back
renamed = rename!(new_df, names(tmp_df))

CSV.write(joinpath(base_path, "loc_table.csv"), renamed)

ext: jl


text:
@mpi_handle Win MPI_Win

const WIN_NULL = _Win(MPI_WIN_NULL)
Win() = Win(WIN_NULL.val)

struct LockType
    val::Cint
end
const LOCK_EXCLUSIVE = LockType(MPI_LOCK_EXCLUSIVE)
const LOCK_SHARED = LockType(MPI_LOCK_SHARED)

"""
    MPI.Win_create(base::Array, comm::Comm; infokws...)

Create a window over the array `base`, returning a `Win` object used by these processes to perform RMA operations

This is a collective call over `comm`.

`infokws` are info keys providing optimization hints.

[`MPI.free`](@ref) should be called on the `Win` object once operations have been completed.
"""
function Win_create(base::AbstractArray{T}, comm::Comm; infokws...) where T
    win = Win()
    # int MPI_Win_create(void *base, MPI_Aint size, int disp_unit, MPI_Info info,
    #                    MPI_Comm comm, MPI_Win *win)
    @mpichk ccall(:MPI_Win_create, Cint,
                  (MPIPtr, Cptrdiff_t, Cint, MPI_Info, MPI_Comm, Ptr{MPI_Win}),
                  base, Cptrdiff_t(length(base)*sizeof(T)), sizeof(T), Info(infokws...), comm, win)
    finalizer(free, win)
    return win
end

"""
    MPI.Win_create_dynamic(comm::Comm; infokws...)

Create a dynamic window returning a `Win` object used by these processes to perform RMA operations

This is a collective call over `comm`.

`infokws` are info keys providing optimization hints.

[`MPI.free`](@ref) should be called on the `Win` object once operations have been completed.
"""
function Win_create_dynamic(comm::Comm; kwargs...)
    win = Win()
    # int MPI_Win_create_dynamic(MPI_Info info, MPI_Comm comm, MPI_Win *win)
    @mpichk ccall(:MPI_Win_create_dynamic, Cint,
                  (MPI_Info, MPI_Comm, Ptr{MPI_Win}),
                  Info(kwargs...), comm, win)
    finalizer(free, win)
    return win
end

"""
    (win, ptr) = MPI.Win_allocate_shared(T, len, comm::Comm; infokws...)

Create and allocate a shared memory window for objects of type `T` of length `len`,
returning a `Win` and a `Ptr{T}` object used by these processes to perform RMA operations

This is a collective call over `comm`.

`infokws` are info keys providing optimization hints.

[`MPI.free`](@ref) should be called on the `Win` object once operations have been completed.
"""
function Win_allocate_shared(::Type{T}, len::Int, comm::Comm; kwargs...) where T
    win = Win()
    out_baseptr = Ref{Ptr{T}}()
    # int MPI_Win_allocate_shared(MPI_Aint size, int disp_unit, MPI_Info info,
    #                             MPI_Comm comm, void *baseptr, MPI_Win *win)
    @mpichk ccall(:MPI_Win_allocate_shared, Cint,
                  (Cptrdiff_t, Cint, MPI_Info, MPI_Comm, Ref{Ptr{T}}, Ptr{MPI_Win}),
                  Cptrdiff_t(len*sizeof(T)), sizeof(T), Info(kwargs...), comm, out_baseptr, win)
    finalizer(free, win)
    return win, out_baseptr[]
end

function free(win::Win)
    if win.val != WIN_NULL.val && !Finalized()
        # int MPI_Win_free(MPI_Win *win)
        @mpichk ccall(:MPI_Win_free, Cint, (Ptr{MPI_Win},), win)
    end
    return nothing
end

function Win_shared_query(win::Win, owner_rank::Int)
    out_len = Ref{Cptrdiff_t}()
    out_sizeT = Ref{Cint}()
    out_baseptr = Ref{Ptr{Cvoid}}()
    # int MPI_Win_shared_query(MPI_Win win, int rank, MPI_Aint *size,
    #                          int *disp_unit, void *baseptr)
    @mpichk ccall(:MPI_Win_shared_query, Cint,
                  (MPI_Win, Cint, Ref{Cptrdiff_t}, Ref{Cint}, Ref{Ptr{Cvoid}}),
                  win, owner_rank, out_len, out_sizeT, out_baseptr)
    out_len[], out_sizeT[], out_baseptr[]
end

function Win_attach(win::Win, base::AbstractArray{T}) where T
    # int MPI_Win_attach(MPI_Win win, void *base, MPI_Aint size)
    @mpichk ccall(:MPI_Win_attach, Cint,
                  (MPI_Win, MPIPtr, Cptrdiff_t),
                  win, base, Cptrdiff_t(sizeof(base)))
end

function Win_detach(win::Win, base::AbstractArray{T}) where T
    # int MPI_Win_detach(MPI_Win win, const void *base)
    @mpichk ccall(:MPI_Win_detach, Cint,
                  (MPI_Win, MPIPtr),
                  win, base)
end

function Win_fence(assert::Integer, win::Win)
    # int MPI_Win_fence(int assert, MPI_Win win)
    @mpichk ccall(:MPI_Win_fence, Cint, (Cint, MPI_Win), assert, win)
end

function Win_flush(rank::Integer, win::Win)
    # int MPI_Win_flush(int rank, MPI_Win win)
    @mpichk ccall(:MPI_Win_flush, Cint,(Cint, MPI_Win), rank, win)
end

function Win_sync(win::Win)
    # int MPI_Win_sync(MPI_Win win)
    @mpichk ccall(:MPI_Win_sync, Cint, (MPI_Win,), win)
end

function Win_lock(lock_type::LockType, rank::Integer, assert::Integer, win::Win)
    # int MPI_Win_lock(int lock_type, int rank, int assert, MPI_Win win)
    @mpichk ccall(:MPI_Win_lock, Cint,
                  (Cint, Cint, Cint, MPI_Win),
                  lock_type.val, rank, assert, win)
end

function Win_unlock(rank::Integer, win::Win)
    # int MPI_Win_unlock(int rank, MPI_Win win)
    @mpichk ccall(:MPI_Win_unlock, Cint, (Cint, MPI_Win), rank, win)
end

function Get(origin_buffer, count::Integer, target_rank::Integer, target_disp::Integer, win::Win)
    T = eltype(origin_buffer)
    # int MPI_Get(void *origin_addr, int origin_count,
    #             MPI_Datatype origin_datatype, int target_rank,
    #             MPI_Aint target_disp, int target_count,
    #             MPI_Datatype target_datatype, MPI_Win win)
    @mpichk ccall(:MPI_Get, Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cptrdiff_t, Cint, MPI_Datatype, MPI_Win),
                  origin_buffer, count, Datatype(T), target_rank, Cptrdiff_t(target_disp), count, Datatype(T), win)
end
function Get(origin_buffer::AbstractArray{T}, target_rank::Integer, win::Win) where T
    count = length(origin_buffer)
    Get(origin_buffer, count, target_rank, 0, win)
end
function Get(origin_value::Ref{T}, target_rank::Integer, win::Win) where T
    Get(origin_value, 1, target_rank, 0, win)
end

function Put(origin_buffer, count::Integer, target_rank::Integer, target_disp::Integer, win::Win)
    # int MPI_Put(const void *origin_addr, int origin_count,
    #             MPI_Datatype origin_datatype, int target_rank,
    #             MPI_Aint target_disp, int target_count,
    #             MPI_Datatype target_datatype, MPI_Win win)
    T = eltype(origin_buffer)
    @mpichk ccall(:MPI_Put, Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cptrdiff_t, Cint, MPI_Datatype, MPI_Win),
                  origin_buffer, count, Datatype(T), target_rank, Cptrdiff_t(target_disp), count, Datatype(T), win)
end
function Put(origin_buffer::AbstractArray{T}, target_rank::Integer, win::Win) where T
    count = length(origin_buffer)
    Put(origin_buffer, count, target_rank, 0, win)
end
function Put(origin_value::Ref{T}, target_rank::Integer, win::Win) where T
    Put(origin_value, 1, target_rank, 0, win)
end

function Fetch_and_op(sourceval, returnval, target_rank::Integer, target_disp::Integer, op::Op, win::Win)
    # int MPI_Fetch_and_op(const void *origin_addr, void *result_addr,
    #                      MPI_Datatype datatype, int target_rank, MPI_Aint target_disp,
    #                      MPI_Op op, MPI_Win win)
    @assert eltype(sourceval) == eltype(returnval)
    T = eltype(sourceval)
    @mpichk ccall(:MPI_Fetch_and_op, Cint,
                  (MPIPtr, MPIPtr, MPI_Datatype, Cint, Cptrdiff_t, MPI_Op, MPI_Win),
                  sourceval, returnval, Datatype(T), target_rank, target_disp, op, win)
end

function Accumulate(origin_buffer, count::Integer, target_rank::Integer, target_disp::Integer, op::Op, win::Win)
    # int MPI_Accumulate(const void *origin_addr, int origin_count,
    #                    MPI_Datatype origin_datatype, int target_rank,
    #                    MPI_Aint target_disp, int target_count,
    #                    MPI_Datatype target_datatype, MPI_Op op, MPI_Win win)
    T = eltype(origin_buffer)
    @mpichk ccall(:MPI_Accumulate, Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cptrdiff_t, Cint, MPI_Datatype, MPI_Op, MPI_Win),
                  origin_buffer, count, Datatype(T), target_rank, Cptrdiff_t(target_disp), count, Datatype(T), op, win)
end

function Get_accumulate(origin_buffer, result_buffer, count::Integer, target_rank::Integer, target_disp::Integer, op::Op, win::Win)
    # int MPI_Get_accumulate(const void *origin_addr, int origin_count,
    #                        MPI_Datatype origin_datatype, void *result_addr,
    #                        int result_count, MPI_Datatype result_datatype,
    #                        int target_rank, MPI_Aint target_disp, int target_count,
    #                        MPI_Datatype target_datatype, MPI_Op op, MPI_Win win)
    @assert eltype(origin_buffer) == eltype(result_buffer)
    T = eltype(origin_buffer)
    @mpichk ccall(:MPI_Get_accumulate, Cint,
                  (MPIPtr, Cint, MPI_Datatype, MPIPtr, Cint, MPI_Datatype, Cint, Cptrdiff_t, Cint, MPI_Datatype, MPI_Op, MPI_Win),
                  origin_buffer, count, Datatype(T), result_buffer, count, Datatype(T), target_rank, Cptrdiff_t(target_disp), count, Datatype(T), op, win)
end


text:
#######################
# ApproximateKNNGraph #
#######################

# Abstract KNN Graph definitions
"""
ApproximateKNNGraph{V, U, D, M} subtypes are weighted, directed graphs where each vertex
has exactly `k` forward edges with weights of type `U`.

`D` is the type of the dataset corresponding to this graph, and `M` is a `Distances.PreMetric`
with result type matching `U`.
"""
abstract type ApproximateKNNGraph{V<:Integer,
                                  U<:Real,
                                  D<:AbstractVector,
                                  M<:PreMetric} <: AbstractGraph{V} end

# interface
Base.eltype(::ApproximateKNNGraph{V}) where V = V
# all knn graphs are directed
LightGraphs.is_directed(::Type{<:ApproximateKNNGraph}) = true
LightGraphs.is_directed(::ApproximateKNNGraph) = true

"""
    knn_diameter(g::ApproximateKNNGraph{V}, v::V)

Compute the diameter of the ball centered on `v` that covers
all of `v`s approximate k-nearest neighbors.
"""
function knn_diameter(g::ApproximateKNNGraph{V}, v::V) where V
    neighbs = outneighbors(g, v)
    return 2 * maximum(weights(g)[neighbs, v])
end


text:
__precompile__()

module Master

using Types
export masternodes, gaussnodes, mkshape, mkmaster;

include("masternodes.jl");
include("gaussnodes.jl");
include("legendrepolynomial.jl");
include("simplexmonomial.jl");
include("tensorproduct.jl");
include("mkshape.jl");

struct MASTERStruct
    nd::IntP; # physical dimension
    porder::IntP; # polynomial degree
    pgauss::IntP; # order of Gauss quadrature
    elemtype::IntP; # type of elements
    nodetype::IntP; # type of nodes
    npe::IntP; # number of polynomials per element
    npf::IntP; # number of polynomials per face
    nge::IntP; # number of Gauss points per element
    ngf::IntP; # number of Gauss points per face
    np1d::IntP; # number of polynomials per 1D element
    ng1d::IntP; # number of Gauss points per 1D element
    perm::Array{IntP,2}; # to get face nodes from element nodes
    shapegt::Array{FloatP,3}; # element shape functions at Gauss points (transpose)
    shapegw::Array{FloatP,3}; # element shape functions at Gauss points multiplied by Gauss weights
    shapfgt::Array{FloatP,3}; # face shape functions at Gauss points (transpose)
    shapfgw::Array{FloatP,3}; # face shape functions at Gauss points multiplied by Gauss weights
    shapent::Array{FloatP,3}; # element shape functions at nodes (transpose)
    shapen::Array{FloatP,3};  # element shape functions at nodes
    shapfnt::Array{FloatP,3}; # face shape functions at nodes (transpose)
    shapfn::Array{FloatP,3};  # face shape functions at nodes
    xpe::Array{FloatP,2}; # nodal points on master element
    gpe::Array{FloatP,2}; # gauss points on master element
    gwe::Array{FloatP,2}; # gauss weighs on master element
    xpf::Array{FloatP,2}; # nodal points on master face
    gpf::Array{FloatP,2}; # gauss points on master face
    gwf::Array{FloatP,2}; # gauss weighs on master face
    shap1dgt::Array{FloatP,3}; # 1D shape functions at Gauss points (transpose)
    shap1dgw::Array{FloatP,3}; # 1D shape functions at Gauss points multiplied by Gauss weights
    shap1dnt::Array{FloatP,3}; # 1D shape functions at nodes (transpose)
    shap1dn::Array{FloatP,3}; # 1D shape functions at nodes
    xp1d::Array{FloatP,2}; # node points on 1D element
    gp1d::Array{FloatP,2}; # gauss points on 1D element
    gw1d::Array{FloatP,2}; # gauss weights on 1D element
    telem::Array{IntP,2}; # element
    tface::Array{IntP,2}; # face
end

function mkmaster(dim::IntP,porder::IntP,pgauss::IntP,elemtype::IntP,nodetype::IntP)

# node positions on the master element and master face
xpe,telem,xpf,tface,perm = masternodes(porder,dim,elemtype);

# Gauss points and weights on the master volume element
gpe, gwe = gaussnodes(pgauss,dim,elemtype);

# shape functions and derivatives on the master volume element
shapeg = mkshape(porder,xpe,gpe,elemtype);
shapen = mkshape(porder,xpe,xpe,elemtype);

npe,nge,nd1 = size(shapeg);
shapegt = zeros(nge,npe,nd1);
shapent = zeros(npe,npe,nd1);
shapegw = zeros(npe,nge,nd1);
Mgwe = zeros(nge,nge);
for i = 1:nge
    Mgwe[i,i] = gwe[i];
end
for d=1:nd1
    shapegt[:,:,d] = shapeg[:,:,d]';
    shapent[:,:,d] = shapen[:,:,d]';
    shapegw[:,:,d] = shapeg[:,:,d]*Mgwe;
end

if dim>1
    # Gauss points and weights on the master face element
    gpf, gwf = gaussnodes(pgauss,dim-1,elemtype);
    # shape functions and derivatives on the master face element
    shapfg = mkshape(porder,xpf,gpf,elemtype);
    shapfn = mkshape(porder,xpf,xpf,elemtype);
else
    gpf   = [0.0];
    gwf   = [1.0];
    shapfg = reshape([1.0],1,1,1);
    shapfn = reshape([1.0],1,1,1);
end

npf = size(shapfg,1);
ngf = size(shapfg,2);
shapfgt = zeros(ngf,npf,dim);
shapfnt = zeros(npf,npf,dim);
shapfgw = zeros(npf,ngf,dim);
Mgwf = zeros(ngf,ngf);
for i = 1:ngf
    Mgwf[i,i] = gwf[i];
end
for d=1:dim
    shapfgt[:,:,d] = shapfg[:,:,d]';
    shapfnt[:,:,d] = shapfn[:,:,d]';
    shapfgw[:,:,d] = shapfg[:,:,d]*Mgwf;
end

xp1d,~,~,~,~ = masternodes(porder,1,elemtype);
gp1d,gw1d = gaussnodes(pgauss,1,elemtype);
shap1dg = mkshape(porder,xp1d,gp1d,elemtype);
shap1dn = mkshape(porder,xp1d,xp1d,elemtype);

np1d = size(shap1dg,1);
ng1d = size(shap1dg,2);
shap1dgt = zeros(ng1d,np1d,2);
shap1dnt = zeros(np1d,np1d,2);
shap1dgw = zeros(np1d,ng1d,2);
Mgw1d = zeros(ng1d,ng1d);
for i = 1:ng1d
    Mgw1d[i,i] = gw1d[i];
end
for d=1:2
    shap1dgt[:,:,d] = shap1dg[:,:,d]';
    shap1dgw[:,:,d] = shap1dg[:,:,d]*Mgw1d;
    shap1dnt[:,:,d] = shap1dn[:,:,d]';
end

master = MASTERStruct(dim, porder, pgauss, elemtype, nodetype, npe, npf,
    nge, ngf, np1d, ng1d, perm, shapegt, shapegw, shapfgt, shapfgw, shapent, shapen,
    shapfnt, shapfn, xpe, gpe, gwe, xpf, gpf, gwf, shap1dgt, shap1dgw, shap1dnt,
    shap1dn, xp1d, gp1d, gw1d, telem, tface);

return master;

end

function mkmaster(app)
    master = mkmaster(app.nd,app.porder,app.pgauss,app.elemtype,app.nodetype);
    return master;
end

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using NNlib, Test

ACTIVATION_FUNCTIONS = [σ, relu, leakyrelu, elu, gelu, swish, selu, softplus, softsign, logcosh];

function test_value_float_precision_preserving(a)
    @testset "$(a): " begin
        for T in [Float32, Float64]
            for val in [-10, -1, 0, 1, 10]
                val = @inferred a(T(val))
                @test typeof(val) == T
            end
        end
    end
end

function test_value_int_input_forces_float64(a)
    @testset "$(a): " begin
        for T in [Int32, Int64]
            for val in [-10, -1, 0, 1, 10]
                val = @inferred a(T(val))
                @test typeof(val) == Float64
            end
        end
    end
end

@testset "Activation Functions" begin
    @test σ(0.0) == 0.5
    @test relu(0.0) == 0.0
    @test leakyrelu(0.0) == 0.0
    @test elu(0.0) == 0.0
    @test gelu(0.0) == 0.0
    @test swish(0.0) == 0.0
    @test softplus(0.0) ≈ log(2.0)
    @test softplus(1e8) ≈ 1e8
    @test softplus(-1e8) ≈ 0.0
    @test softsign(0.0) == 0.0
    @test selu(0.0) == 0.0
    @test logcosh(0.0) == log(cosh(0.0))

    @test σ(1.0) == 1.0 / (1.0 + exp(-1.0))
    @test relu(1.0) == 1.0
    @test leakyrelu(1.0) == 1.0
    @test elu(1.0) == 1.0
    @test gelu(1.0) == 0.8411919906082768
    @test swish(1.0) == 1.0 / (1.0 + exp(-1.0))
    @test softplus(1.0) ≈ log(exp(1.0) + 1.0)
    @test softsign(1.0) == 0.5
    @test selu(1.0) == 1.0507009873554804934193349852946
    @test logcosh(1.0) ≈ log(cosh(1.0))

    @test σ(-1.0) == 1.0 / (1.0 + exp(1.0))
    @test relu(-1.0) == 0.0
    @test leakyrelu(-1.0) == -0.01
    @test elu(-1.0) == exp(-1.0) - 1.0
    @test gelu(-1.0) == -0.15880800939172324
    @test swish(-1.0) == -1.0 / (1.0 + exp(1.0))
    @test softplus(-1.0) ≈ log(exp(-1.0) + 1.0)
    @test softsign(-1.0) == -0.5
    @test selu(-1.0) == 1.0507009873554804934193349852946 * 1.6732632423543772848170429916717 * (exp(-1.0) - 1.0)
    @test log(cosh(-1.0)) ≈ log(cosh(-1.0))

    @testset "Float inference" begin
        test_value_float_precision_preserving.(ACTIVATION_FUNCTIONS)
    end

    @testset "Array input" begin
        x = rand(5)
        for a in ACTIVATION_FUNCTIONS
            @test_throws ErrorException a(x)
        end
    end

    @testset "Test Integer64 and Integer32 inputs will force Float64 outputs" begin
        test_value_int_input_forces_float64.(filter(x -> x != relu, ACTIVATION_FUNCTIONS))

        @testset "relu: " begin
            # relu doesn't have to force floating point outputs
            @test typeof(relu(Int64(1))) == Int64
            @test typeof(relu(Int32(1))) == Int32
        end
    end

    @testset "softmax" begin
        xs = rand(5,5)
        @test all(sum(softmax(xs), dims = 1) .≈ 1)
        @test sum(softmax(vec(xs))) ≈ 1

        xs = [-100_000, -100_000.]
        @test softmax(xs) ≈ [0.5, 0.5]
        @test logsoftmax(xs) ≈ log.([0.5, 0.5])

        xs = rand(5)
        @test softmax(xs) ≈ exp.(xs) ./ sum(exp.(xs))
        @test logsoftmax(xs) ≈ log.(softmax(xs))

        xs = Float32[1, 2, 3000.]
        @test logsoftmax(xs) ≈ [-2999, -2998, 0]

        xs = Float32[1 2 3; 1000 2000 3000]
        @test logsoftmax(xs) ≈ [-999 -1998 -2997; 0 0 0.]

        @test NNlib.∇logsoftmax(ones(size(xs)), xs) ≈ zeros(Float32, size(xs))
        @test NNlib.∇softmax(ones(size(xs)), xs) ≈ zeros(Float32, size(xs))

        # These values precalculated using PyTorch's nn.LogSoftmax
        xs = [
            -0.238639  0.748142 -0.283194 -0.525461 -1.5348   -0.797842;
             0.690384  0.211427  0.254794 -0.213572 -0.314174 -0.372663;
            -1.146370 -0.577988  0.718952  0.919720 -0.620773  0.929977
        ]
        ys = [
            0.237703 -0.621474 0.448193 0.546047 0.564185 0.632273;
            -0.930163 0.0519798 0.0549979 0.3799 -0.477112 0.437428;
            0.69246 0.569494 -0.503191 -0.925947 -0.0870738 -1.0697
        ]
        @test isapprox(NNlib.∇logsoftmax(ones(size(xs)), xs), ys; rtol = 1e-6)
        @test isapprox(NNlib.∇softmax(ones(size(xs)), xs), zeros(size(xs)); atol = 1e-6)
    end

    @testset "elu" begin
        @test elu(42) == 42
        @test elu(42.) == 42.

        @test elu(-4) ≈ (exp(-4) - 1)
    end

    @test leakyrelu( 0.4,0.3) ≈  0.4
    @test leakyrelu(-0.4,0.3) ≈ -0.12

    @testset "logsigmoid" begin
        xs = randn(10,10)
        @test logsigmoid.(xs) ≈ log.(sigmoid.(xs))
        for T in [:Float32, :Float64]
            @eval @test logsigmoid.($T[-100_000, 100_000.]) ≈ $T[-100_000, 0.]
        end
    end

    @test logcosh(1_000.0) + log(2) == 1_000.0
end


text:
module EXPRES-StellarSignals

# Write your package code here.

include("io.jl")
export  make_results_df_results

end


text:
"""
## Indexing
`pa` is an instance of `SeisForwExpt`.
* `pa[:data,i]` : get data for ith supersource (defaults to first), after `update!`
* `pa[:ageom]` : get `ageom` originally input while creating `pa`
* `pa[:srcwav]`: `srcwav` used to create `pa`
"""
function Base.getindex(pa::PBorn, s::Symbol, iss::Int = 1)
    @assert iss ≤ length(pa.ageom)
    @assert s in [:data, :ageom, :srcwav]
    if (s == :data)
        (iss == 0) ? (return pa.data) : (return pa.data[iss])
    elseif (s == :ageom)
        return pa.ageom
    elseif (s == :srcwav)
        return pa.srcwav
    end
end

function Base.show(io::Base.IO, pa::PFdtd)
    println(typeof(pa), "")
    println("pa[:data] : modeled data after running `update!`")
    println("pa[:ageom]")
    println("pa[:srcwav]")
end


text:
module Terrain

include("heightfield.jl")
include("faultline.jl")

end # module


text:
struct GreensMeasurement{OT <: AbstractObservable} <: AbstractMeasurement
    obs::OT
end
function save_measurement(file::JLDFile, m::GreensMeasurement, entryname::String)
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(m))
    write(file, entryname * "/obs", m.obs)
    nothing
end
function _load(data, ::Type{T}) where T <: GreensMeasurement
    if data["VERSION"] == 1
        data["type"](data["obs"])
    else
        _load(data, AbstractMeasurement)
    end
end

struct BosonEnergyMeasurement{OT <: AbstractObservable} <: AbstractMeasurement
    obs::OT
end
function save_measurement(file::JLDFile, m::BosonEnergyMeasurement, entryname::String)
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(m))
    write(file, entryname * "/obs", m.obs)
    nothing
end
function _load(data, ::Type{T}) where T <: BosonEnergyMeasurement
    if data["VERSION"] == 1
        data["type"](data["obs"])
    else
        _load(data, AbstractMeasurement)
    end
end

struct OccupationMeasurement{OT <: AbstractObservable} <: AbstractMeasurement
    obs::OT
end
function OccupationMeasurement(m::GreensMeasurement{<: LightObservable}; capacity=capacity(m.obs.B))
    N = size(m.obs)[1]
    o = LightObservable(
        LogBinner(zeros(Float64, N), capacity=capacity),
        "State resolved Occupation",
        "Observables.jld",
        "Occ"
    )

    N = min(length(o.B.count), length(m.obs.B.count))
    if m.obs.B.count[N] > 1
        throw(OverflowError("OccupationMeasurement does not have enough capacity!"))
    end
    N = findlast(>(1), m.obs.B.count)
    N === nothing && return OccupationMeasurement{typeof(o)}(o)

    for i in 1:N
        o.B.x_sum[i]  .= m.obs.B.count[i] .- diag(m.obs.B.x_sum[i])
        o.B.x2_sum[i] .= m.obs.B.count[i] .- 2diag(m.obs.B.x_sum[i]) .+ diag(m.obs.B.x2_sum[i])
        o.B.count[i]   = m.obs.B.count[i]
        o.B.compressors[i].value .= 1 .- diag(m.obs.B.compressors[i].value)
        o.B.compressors[i].switch = m.obs.B.compressors[i].switch
    end
    OccupationMeasurement{typeof(o)}(o)
end
function save_measurement(file::JLDFile, m::OccupationMeasurement, entryname::String)
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(m))
    write(file, entryname * "/obs", m.obs)
    nothing
end
function _load(data, ::Type{T}) where T <: OccupationMeasurement
    if data["VERSION"] == 1
        data["type"](data["obs"])
    else
        _load(data, AbstractMeasurement)
    end
end

_get_shape(model::Model) = (length(lattice(model)),)
_get_shape(mask::RawMask) = (mask.nsites, mask.nsites)
_get_shape(mask::DistanceMask) = length(mask)
mask(m::AbstractMeasurement) = m.mask

abstract type SpinOneHalfMeasurement <: AbstractMeasurement end

struct ChargeDensityCorrelationMeasurement{
        OT <: AbstractObservable,
        AT <: Array,
        MT <: AbstractMask
    } <: SpinOneHalfMeasurement
    obs::OT
    temp::AT
    mask::MT
end
function save_measurement(file::JLDFile, m::ChargeDensityCorrelationMeasurement, entryname::String)
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(m))
    write(file, entryname * "/obs", m.obs)
    nothing
end
function _load(data, ::Type{T}) where T <: ChargeDensityCorrelationMeasurement
    if data["VERSION"] == 1
        OccupationMeasurement(data["obs"])
    else
        _load(data, AbstractMeasurement)
    end
end

struct MagnetizationMeasurement{
        OTx <: AbstractObservable,
        OTy <: AbstractObservable,
        OTz <: AbstractObservable,
        AT <: AbstractArray,
        ATy <: AbstractArray
    } <: SpinOneHalfMeasurement

    x::OTx
    y::OTy
    z::OTz
    temp::AT
    tempy::ATy
end
function save_measurement(file::JLDFile, m::MagnetizationMeasurement, entryname::String)
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", MagnetizationMeasurement)
    write(file, entryname * "/x", m.x)
    write(file, entryname * "/y", m.y)
    write(file, entryname * "/z", m.z)
    nothing
end
function _load(data, ::Type{T}) where T <: MagnetizationMeasurement
    if data["VERSION"] == 1
        x = data["x"]
        y = data["y"]
        z = data["z"]
        temp = similar(x.B.x_sum[1])
        tempy = similar(y.B.x_sum[1])
        data["type"](x, y, z, temp, tempy)
    else
        _load(data, AbstractMeasurement)
    end
end

struct SpinDensityCorrelationMeasurement{
        OTx <: AbstractObservable,
        OTy <: AbstractObservable,
        OTz <: AbstractObservable,
        AT <: Array,
        MT <: AbstractMask
    } <: SpinOneHalfMeasurement

    x::OTx
    y::OTy
    z::OTz
    temp::AT
    mask::MT
end

struct PairingCorrelationMeasurement{
        OT <: AbstractObservable,
        AT <: Array,
        MT <: DistanceMask
    } <: SpinOneHalfMeasurement
    obs::OT
    temp::AT
    mask::MT
    rsm::RestrictedSourceMask
end

to_logbinner(B::LogBinner) = B
to_logbinner(B::LightObservable) = B.B

function mask2iter(mask, iter::EachSitePairByDistance)
    # A[idxs] should output iter order
    # idxs = [mask_idx@iter1, mask_idx@iter2, ...]
    idxs = Vector{Int64}(undef, ndirections(iter))

    N = sqrt(length(iter))
    @assert N ≈ round(Int64, N)
    N = round(Int64, N)
    mask_idxs = Matrix{Int64}(undef, N, N)
    for (idx, src, trg) in getorder(mask)
        mask_idxs[src, trg] = idx
    end

    for i in 1:ndirections(iter)
        pairs = in_direction(iter, i)
        x = first(pairs)
        idx = mask_idxs[x[1], x[2]]
        for (src, trg) in pairs
            if mask_idxs[src, trg] != idx
                error("Mapping is not unqiue: $(mask_idxs[src, trg]) != $(idx)")
            end
        end
        idxs[i] = idx
    end

    idxs
end

function reorder(input::LogBinner, idxs::Vector)
    B = deepcopy(input)
    for lvl in eachindex(B.compressors)
        permute!(B.compressors[lvl].value, idxs)
        permute!(B.x_sum[lvl], idxs)
        permute!(B.x2_sum[lvl], idxs)
    end
    B
end

function reorder(input::LogBinner, idxs1::Vector, idxs2::Vector, idxs3::Vector)
    B = deepcopy(input)
    for lvl in eachindex(B.compressors)
        B.compressors[lvl].value = B.compressors[lvl].value[idxs1, idxs2, idxs3]
        B.x_sum[lvl] = B.x_sum[lvl][idxs1, idxs2, idxs3]
        B.x2_sum[lvl] = B.x2_sum[lvl][idxs1, idxs2, idxs3]
    end
    B
end

function DQMCMeasurement(mc::DQMC, m::GreensMeasurement)
    obs = to_logbinner(m.obs)
    temp = mean(obs)
    DQMCMeasurement{Greens, Nothing}(greens_kernel, obs, temp)
end

function DQMCMeasurement(mc::DQMC, m::BosonEnergyMeasurement)
    obs = to_logbinner(m.obs)
    DQMCMeasurement{Nothing, Nothing}(energy_boson, obs, 0.0)
end

function DQMCMeasurement(mc::DQMC, m::OccupationMeasurement)
    obs = to_logbinner(m.obs)
    temp = mean(obs)
    DQMCMeasurement{Greens, EachSiteAndFlavor}(occupation_kernel, obs, temp)
end

function DQMCMeasurement(mc::DQMC, m::MagnetizationMeasurement)
    x = to_logbinner(m.x)
    y = to_logbinner(m.y)
    z = to_logbinner(m.z)

    return (
        DQMCMeasurement{Greens, EachSite}(mx_kernel, x, mean(x)),
        DQMCMeasurement{Greens, EachSite}(my_kernel, y, mean(y)),
        DQMCMeasurement{Greens, EachSite}(mz_kernel, z, mean(z))
    )
end

function DQMCMeasurement(mc::DQMC, m::ChargeDensityCorrelationMeasurement)
    iter = EachSitePairByDistance(mc, mc.model)
    idxs = mask2iter(m.mask, iter)
    obs = reorder(to_logbinner(m.obs), idxs)
    temp = mean(obs)
    DQMCMeasurement{Greens, EachSitePairByDistance}(cdc_kernel, obs, temp)
end

function DQMCMeasurement(mc::DQMC, m::SpinDensityCorrelationMeasurement)
    iter = EachSitePairByDistance(mc, mc.model)
    idxs = mask2iter(m.mask, iter)

    x = reorder(to_logbinner(m.x), idxs)
    y = reorder(to_logbinner(m.y), idxs)
    z = reorder(to_logbinner(m.z), idxs)
    temp = mean(x)
    tempy = mean(y)
    return (
        DQMCMeasurement{Greens, EachSitePairByDistance}(sdc_x_kernel, x, temp),
        DQMCMeasurement{Greens, EachSitePairByDistance}(sdc_y_kernel, y, tempy),
        DQMCMeasurement{Greens, EachSitePairByDistance}(sdc_z_kernel, z, temp)
    )
end

function DQMCMeasurement(mc::DQMC, m::PairingCorrelationMeasurement)
    iter = EachSitePairByDistance(mc, mc.model)
    idxs = mask2iter(m.mask, iter)
    pre = to_logbinner(m.obs)
    temp = mean(pre)
    D = size(temp, 2)
    obs = reorder(pre, idxs, idxs[1:D], idxs[1:D])
    DQMCMeasurement{Greens, EachLocalQuadByDistance{D}}(pc_kernel, obs, mean(obs))
end


text:
include("../src/qspice.jl")

using Benchmarks

using QSpice.Gates, QSpice.State

s = fromstates(BELL_STATE, QUBIT0, QUBIT1,
               QUBIT0, QUBIT1, QUBIT1, BELL_STATE,
               BELL_STATE, BELL_STATE)

cn = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]

Benchmarks.@benchmarkable(hadamard!,         nothing, hadamard(s, rand(1:s.bits)),                              nothing)
Benchmarks.@benchmarkable(not!,              nothing, not(s, rand(1:s.bits)),                                   nothing)
Benchmarks.@benchmarkable(cnot!,             nothing, cnot(s, rand(1:s.bits), rand(1:s.bits)),                  nothing)
Benchmarks.@benchmarkable(ccnot!,            nothing, ccnot(s, rand(1:s.bits), rand(1:s.bits), rand(1:s.bits)), nothing)
Benchmarks.@benchmarkable(swap!,             nothing, swap(s, rand(1:s.bits), rand(1:s.bits)),                  nothing)
Benchmarks.@benchmarkable(cswap!,            nothing, cswap(s, rand(1:s.bits), rand(1:s.bits), rand(1:s.bits)), nothing)
Benchmarks.@benchmarkable(sqrtswap!,         nothing, sqrtswap(s, rand(1:s.bits), rand(1:s.bits)),              nothing)
Benchmarks.@benchmarkable(phaseshift!,       nothing, phaseshift(s, rand(1:s.bits), rand(1:2pi)),               nothing)
Benchmarks.@benchmarkable(paulix!,           nothing, paulix(s, rand(1:s.bits)),                                nothing)
Benchmarks.@benchmarkable(pauliy!,           nothing, pauliy(s, rand(1:s.bits)),                                nothing)
Benchmarks.@benchmarkable(pauliz!,           nothing, pauliz(s, rand(1:s.bits)),                                nothing)
Benchmarks.@benchmarkable(measure!,          nothing, measure(s),                                               nothing)
Benchmarks.@benchmarkable(partialmeasure!,   nothing, partialmeasure(s, rand(1:s.bits), rand(1:s.bits),
                                                                        rand(1:s.bits), rand(1:s.bits),
                                                                        rand(1:s.bits), rand(1:s.bits)),        nothing)
Benchmarks.@benchmarkable(unitarycnot!,      nothing, unitary(s, cn, 4, 3),                                     nothing)
Benchmarks.@benchmarkable(qft!,              nothing, qft(s),                                                   nothing)

benchmarks = [hadamard!, not!, cnot!, ccnot!, swap!, cswap!, sqrtswap!, phaseshift!, paulix!,
              pauliy!, pauliz!, measure!, partialmeasure!, unitarycnot!, qft!]

for b in benchmarks
    println("Benchmark: ", b)
    result = Benchmarks.execute(b)
    println(result)
end

max_stars_repo_name: leios/quantum

max_issues_repo_name: leios/quantum

max_forks_repo_name: leios/quantum


text:
using Pkgx
using Test

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


text:
@testset "context_implementations.jl" begin
    # https://github.com/TuringLang/DynamicPPL.jl/issues/129
    @testset "#129" begin
        @model function test(x)
            μ ~ MvNormal(fill(0, 2), 2.0)
            z = Vector{Int}(undef, length(x))
            # `z .~ Categorical.(ps)` cannot be parsed by Julia 1.0
            (.~)(z, Categorical.(fill([0.5, 0.5], length(x))))
            for i in 1:length(x)
                x[i] ~ Normal(μ[z[i]], 0.1)
            end
        end

        test([1, 1, -1])(VarInfo(), SampleFromPrior(), LikelihoodContext())
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

function get_test_model()

    X = Dict{Int, Vector{Float64}}(
         1 => [0.0, 0.0],  2 => [2.0, 0.0],
         3 => [0.0, 1.0],  4 => [2.0, 1.0],
         5 => [0.0, 1.0],  6 => [1.3, 1.0],
         7 => [0.0, 2.0],  8 => [1.3, 2.0],
         9 => [1.3, 1.0], 10 => [2.0, 1.0],
        11 => [1.3, 2.0], 12 => [2.0, 2.0])

    T = Dict{Int, Vector{Float64}}(
         7 => [0.0, 288.0],  8 => [0.0, 288.0],
        11 => [0.0, 288.0], 12 => [0.0, 288.0])

    # volume elements, three bodies
    e1 = Element(Quad4, [1, 2, 4, 3])
    e2 = Element(Quad4, [5, 6, 8, 7])
    e3 = Element(Quad4, [9, 10, 12, 11])
    update!([e1, e2, e3], "geometry", X)
    update!([e1, e2, e3], "youngs modulus", 288.0)
    update!([e1, e2, e3], "poissons ratio", 1/3)
    b1 = Element(Seg2, [7, 8])
    b2 = Element(Seg2, [11, 12])
    update!([b1, b2], "geometry", X)
    update!([b1, b2], "displacement traction force", T)
    body1 = Problem(Elasticity, "body 1", 2)
    body1.properties.formulation = :plane_stress
    push!(body1, e1)
    body2 = Problem(Elasticity, "body 2", 2)
    body2.properties.formulation = :plane_stress
    push!(body2, e2, b1)
    body3 = Problem(Elasticity, "body 3", 2)
    body3.properties.formulation = :plane_stress
    push!(body3, e3, b2)

    # boundary elements for dirichlet dx=0
    dx1 = Element(Seg2, [1, 3])
    dx2 = Element(Seg2, [5, 7])
    update!([dx1, dx2], "geometry", X)
    update!([dx1, dx2], "displacement 1", 0.0)
    bc1 = Problem(Dirichlet, "dx=0", 2, "displacement")
    push!(bc1, dx1, dx2)

    # boundary elements for dirichlet dy=0
    dy1 = Element(Seg2, [1, 2])
    update!(dy1, "geometry", X)
    update!(dy1, "displacement 2", 0.0)
    bc2 = Problem(Dirichlet, "dy=0", 2, "displacement")
    push!(bc2, dy1)

    # mortar boundary between body 1 and body 2
    mel1 = Element(Seg2, [3, 4])
    sel1 = Element(Seg2, [5, 6])
    update!([mel1, sel1], "geometry", X)
    update!(sel1, "master elements", [mel1])
    bc3 = Problem(Mortar, "interface between body 1 and 2", 2, "displacement")
    push!(bc3, mel1, sel1)

    # mortar boundary between body 1 and body 3
    sel2 = Element(Seg2, [9, 10])
    update!(sel2, "geometry", X)
    update!(sel2, "master elements", [mel1])
    bc4 = Problem(Mortar, "interface between body 1 and 3", 2, "displacement")
    push!(bc4, mel1, sel2)

    # mortar boundary between body 2 and body 3
    sel3 = Element(Seg2, [6, 8])
    mel2 = Element(Seg2, [9, 11])
    update!([sel3, mel2], "geometry", X)
    update!(sel3, "master elements", [mel2])
    bc5 = Problem(Mortar, "interface between body 2 and 3", 2, "displacement")
    push!(bc5, sel3, mel2)

    return body1, body2, body3, bc1, bc2, bc3, bc4, bc5
end

#= TODO: Fix test
@testset "test 2d mortar problem with three bodies and shared nodes" begin
    body1, body2, body3, bc1, bc2, bc3, bc4, bc5 = get_test_model()
    solver = Solver(Linear)
    push!(solver, body1, body2, body3, bc1, bc2, bc3, bc4, bc5)
    solver()
    X = e3("geometry", [1.0, 1.0], 0.0)
    u = e3("displacement", [1.0, 1.0], 0.0)
    info("displacement at $X: $u")
    u_expected = [-1/3, 1.0]
    @test isapprox(u, u_expected)
end
=#


text:
# # AutoViz Basics

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/basics.ipynb)

# To run these tutorials, you must be using `julia 1.1+` and have `AutoViz.jl` installed.

# The core function of the `AutoViz.jl` package is `AutoViz.render()`
#
#md # ```@docs
#md # render
#md # ```
#
# In its simplest form, the `render` function only takes one argument:
# an iterable collection of renderable objects, `renderables`.
# In order for an object to be renderable, it needs to implement the
# `add_renderable!()` function.
#
# By default, the `render()` function will do its best to make all the contents of
# the scene fit to the canvas of dimensions `canvas_width` x `canvas_height`.
# More fine-grained control over camera positioning can be achieved by providing
# the `camera` keyword.
# In case a custom camera is used, it is important to call
# `update_camera!(camera, scene)` before rendering.
#
# The Cairo surface to be used for plotting can be specified using the `surface` keyword.
# The `render()` function applies the rendering instructions to the `surface` and
# returns it. In the tutorials, we will denote the returned `surface` as a
# `snapshot` of the scene. Such snapshots can be saved to file using the command
#
# ```julia
# write("snapshot.svg", snapshot)
# ```

# ## Basic rendering

# We start our example by rendering an empty roadway

using AutomotiveDrivingModels
using AutoViz

roadway = gen_straight_roadway(3, 100.0)

snapshot = render([roadway], canvas_height=120)
#md write("empty_roadway.svg", snapshot) # hide

#md # ![empty roadway](empty_roadway.svg)

# We can change the background color and render again

AutoViz.colortheme["background"] = colorant"white"
snapshot = render([roadway], canvas_height=120)
#md write("empty_roadway_whitebg.svg", snapshot) # hide

#md # ![empty roadway](empty_roadway_whitebg.svg)

# Let's add some vehicles to the roadway

car_len = 4.8
car_width = 1.8
def = VehicleDef(AgentClass.CAR, car_len, car_width)
w = DEFAULT_LANE_WIDTH
scene = Scene(Entity{VehicleState,VehicleDef,Int64}, 4)  # allocate a scene for 4 agents

# add three cars
push!.(Ref(scene), [
    Entity(VehicleState(VecSE2(10.0,   w, 0.0), roadway, 4.0 + 2.0randn()), def, 1),
    Entity(VehicleState(VecSE2(40.0, 0.0, 0.0), roadway, 4.0 + 2.0randn()), def, 2),
    Entity(VehicleState(VecSE2(70.0,   w, 0.0), roadway, 4.0 + 2.0randn()), def, 3),
])

# add a pedestrian
push!(scene, Entity(
    VehicleState(VecSE2(50.0, 2w, -π/2), roadway, 2.0),
    VehicleDef(AgentClass.PEDESTRIAN, 1., 1.),
    42
))
snapshot = render([roadway, scene], canvas_height=120)
#md write("roadway_with_cars.svg", snapshot) # hide

#md # ![roadway with cars](roadway_with_cars.svg)

#src TODO: explain different colorthemes

# ## Vehicle Shapes
#
# The `render` function provides some defaults for rendering basic building blocks
# such as entities or roadways.
# If the value of `AutoViz.rendermode` is set to `:basic`, entities are simply
# rendered as rectangles with arrows indicating their velocities.
# This can also be done explicitly via

renderables = [
    roadway,
    (EntityRectangle(entity=x) for x in scene)...,
    (VelocityArrow(entity=x) for x in scene)...,
]
snapshot = render(renderables, canvas_height=120)
#md write("roadway_basic_manual.svg", snapshot) # hide

#md # ![roadway basic manual](roadway_basic_manual.svg)

# The result remains the same. The velocity arrows point to the location at which
# the vehicle would be 1 second in the future.
#
# Setting `AutoViz.rendermode` to `:fancy`, the rectangles are replaced by SVG
# images of cars (or pedestrians).

AutoViz.set_render_mode(:fancy)
snapshot = render([roadway, scene], canvas_height=120)
#md write("roadway_fancy.svg", snapshot) # hide

#md # ![roadway fancy](roadway_fancy.svg)

# Which is shorthand for

renderables = [
    roadway, (FancyCar(car=scene[i]) for i in 1:3)..., FancyPedestrian(ped=scene[4])
]
snapshot = render(renderables, canvas_height=120)
#md write("roadway_fancy_manual.svg", snapshot) # hide

#md # ![roadway fancy manual](roadway_fancy_manual.svg)

# A third visualization mode is available in the form of arrow cars, in which the
# arrow indicates the heading direction of the car but does not scale with speed.

renderables = [
    roadway,
    (ArrowCar(scene[i]) for i in 1:3)...,
    FancyPedestrian(ped=scene[4])
]
snapshot = render(renderables, canvas_height=120)
#md write("roadway_arrow.svg", snapshot) # hide

#md # ![roadway arrow](roadway_arrow.svg)

# It is also possible to render single vehicles

fancy_car = FancyCar(car=Entity(VehicleState(VecSE2(0.,0.), 0.), VehicleDef(), 1))
snapshot = render([fancy_car], camera=StaticCamera(zoom=100.))
#md write("fancy_car.svg", snapshot) # hide

#md # ![fancy car](fancy_car.svg)

# ## Vehicle Colors
#
# By default, the `render` function generates a random color for each entity based
# on its ID using the `id_to_color()` function.
# However, vehicle colors can also be assigned explicitly:

colors = [colorant"red", colorant"blue", colorant"green"]
renderables = [
    roadway,
    (FancyCar(car=scene[i], color=colors[i]) for i in 1:3)...,
    FancyPedestrian(ped=scene[4], color=colorant"yellow")
]
snapshot = render(renderables, canvas_height=120)
#md write("scene_custom_colors.svg", snapshot) # hide

#md # ![cars with custom colors](scene_custom_colors.svg)

# ## Custom Renderable Objects

using Colors

struct MyRenderableCircle
    pos::VecE2
    radius::Float64
    color::Colorant
end

function AutoViz.add_renderable!(rendermodel::RenderModel, circle::MyRenderableCircle)
    ## add the desired render instructions to the rendermodel
    add_instruction!(
        rendermodel, AutoViz.render_circle,
        (circle.pos.x, circle.pos.y, circle.radius, circle.color),
        coordinate_system=:scene
    )
    return rendermodel
end

circles = [MyRenderableCircle(VecE2(4i,3.0*(1+sin(i/4))), .5+rand(), rand(RGB)) for i in 1:20]
snapshot = render([roadway, circles..., scene], canvas_height=120)
#md write("custom_circles.svg", snapshot) # hide

#md # ![custom circles](custom_circles.svg)

# ## Simulation and Animations

# We can simulate the scenario over time and visualize the results using `Reel`
# (based on `ffmpeg`).

using Reel

timestep = 0.1
nticks = 50

models = Dict((i => Tim2DDriver() for i in 1:3))  # car models
models[42] = Tim2DDriver()  # TODO: better pedestrian model

scenes = simulate(scene, roadway, models, nticks, timestep)

animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
    i = Int(floor(t/dt)) + 1
    render([roadway, scenes[i]], canvas_height=120)
end

#md write("roadway_animated.gif", animation) #hide

#md # ![animated roadway](roadway_animated.gif)

# Alternatively, the scene can also be visualized interactively using `Interact`

#src TODO: have this in docs too, not just in notebooks (but should not execute in docs)

#nb using Interact
#nb using Reel
#nb using Blink
#nb
#nb w = Window()
#nb viz = @manipulate for step in 1 : length(scenes)
#nb     render([roadway, scenes[step]])
#nb end
#nb body!(w, viz)

# Another common visualization package is `Reactive`.

max_stars_repo_name: tawheeler/AutoViz

max_issues_repo_name: tawheeler/AutoViz

max_forks_repo_name: tawheeler/AutoViz


text:
# Everything necessary to create a RFF approximation to a prior GP with a stationary kernel
# Currently ony supports SqExponentialKernel with variance and lengthscale
struct RFFBasis
    inner_weights  # lengthscale
    outer_weights  # variance (scaled)
    ω  # Sampled frequencies
    τ  # Sampled phases
end

(ϕ::RFFBasis)(x) = ϕ([x])
function (ϕ::RFFBasis)(x::AbstractVector)
    # size(x): (num_data, input_dims)
    # size(ϕ.ω): (1, num_features)
    x_rescaled = x / ϕ.inner_weights
    ωt_x = x_rescaled * ϕ.ω  # size(ωt_x): (num_data, num_features)
    return ϕ.outer_weights * cos.(ωt_x .+ ϕ.τ')
end

function sample_basis(rng, kernel, num_features=10, input_dims=1)
    p_ω = spectral_distribution(kernel, input_dims)
    ω = rand(rng, p_ω, num_features)
    τ = rand(rng, Uniform(0, 2π), num_features)
    inner, outer = spectral_weights(kernel)
    outer_scaled = outer * √(2/num_features)

    return RFFBasis(inner, outer_scaled, ω, τ)
end

function spectral_distribution(::SqExponentialKernel, input_dims=1)
    return MvNormal(Diagonal(Fill(1., input_dims)))
end

spectral_distribution(k::ScaledKernel, args...) = spectral_distribution(k.kernel, args...)
spectral_distribution(k::TransformedKernel, args...) = spectral_distribution(k.kernel, args...)

spectral_distribution(k::Kernel) =  error("Spectral distribution not implemented for kernel:\n$k")

function spectral_weights(::SqExponentialKernel)
    return 1.0, 1.0
end

function spectral_weights(k::ScaledKernel)
    σ² = k.σ²[1]
    inner, outer = spectral_weights(k.kernel)
    return inner, outer * √σ²
end

function spectral_weights(k::TransformedKernel{<:Any, <:ScaleTransform})
    s = k.transform.s[1]
    inner, outer = spectral_weights(k.kernel)
    return inner / s, outer
end

spectral_weights(k::Kernel) = error("Spectral weights not implemented for kernel:\n$k")

# TODO:
# ARDTransform
# ProductKernel
# SumKernel
# MaternKernel


text:
#FIXME: PyObject not found when pyscf only
import PySCF: scf_rhf, ccsd #, pyscf

# Workaround: repeat the imports.
using PyCall: pyimport
pyscf = pyimport("pyscf")
mp = pyimport("pyscf.mp")   # Had to import mp alone ??!
cc = pyimport("pyscf.cc")   # Had to import mp alone ??!

@testset "CCSD: H₂O" begin

   mol = pyscf.gto.Mole()
   mol.build(verbose = 0,
             basis = "sto-3g",
             atom = """8 0  0.    0
                       1 0 -0.757 0.587
                       1 0  0.757 0.587""",
            )

   eri = mol.intor("cint2e_sph", aosym="s8")
   rhf      = scf_rhf(mol,verbose=false)
   e_rhf_me, e, C = rhf.Emol, rhf.e, rhf.C
   mf = pyscf.scf.RHF(mol) # set the model
   mf.kernel()             # run RHF (populate mf with data)

   e_mp2_pyscf = mp.MP2(mf).kernel()[1]              # pure PySCF MP₂ solution

   e_ccsd_me, e_mp2_me = ccsd(e,C,eri,mol,verbose=false)
   e_ccsd_pyscf     = cc.ccsd.CC(mf).kernel()[1]

   @test isapprox(e_mp2_me, e_mp2_pyscf, atol=1e-7)
   @test isapprox(e_ccsd_me, e_ccsd_pyscf,atol=1e-4)
#   println()
#   println("            E(HF)          Ecorr(MP2)           E(MP2)         Ecorr(CCSD)  ")
#   println("        --------------   --------------   ---------------   ---------------")
#   printfmt(" ME:    {1:13.10f}    {2:13.10f}    {3:13.10f}    {4:13.10f}\n",
#                    e_rhf_me,      e_mp2_me3,     e_rhf_me+e_mp2_me3,    e_ccsd_me)
#   printfmt(" PySCF: {1:13.10f}    {2:13.10f}    {3:13.10f}    {4:13.10f}\n",
#                     e_rhf_pyscf,  e_mp2_pyscf,e_rhf_pyscf+e_mp2_pyscf,e_ccsd_pyscf)
end


text:
module operators

export AbstractOperator, DataOperator, length, basis, dagger, ishermitian, tensor, embed,
        tr, ptrace, normalize, normalize!, expect, variance,
        exp, permutesystems, identityoperator, dense

import Base: ==, +, -, *, /, ^, length, one, exp, conj, conj!, transpose
import LinearAlgebra: tr, ishermitian
import SparseArrays: sparse
import ..bases: basis, tensor, ptrace, permutesystems,
            samebases, check_samebases, multiplicable
import ..states: dagger, normalize, normalize!

using ..sortedindices, ..bases, ..states

"""
Abstract base class for all operators.

All deriving operator classes have to define the fields
`basis_l` and `basis_r` defining the left and right side bases.

For fast time evolution also at least the function
`gemv!(alpha, op::AbstractOperator, x::Ket, beta, result::Ket)` should be
implemented. Many other generic multiplication functions can be defined in
terms of this function and are provided automatically.
"""
abstract type AbstractOperator{BL<:Basis,BR<:Basis} end

"""
Abstract type for operators with a data field.

This is an abstract type for operators that have a direct matrix representation
stored in their `.data` field.
"""
abstract type DataOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} end

# Common error messages
arithmetic_unary_error(funcname, x::AbstractOperator) = throw(ArgumentError("$funcname is not defined for this type of operator: $(typeof(x)).\nTry to convert to another operator type first with e.g. dense() or sparse()."))
arithmetic_binary_error(funcname, a::AbstractOperator, b::AbstractOperator) = throw(ArgumentError("$funcname is not defined for this combination of types of operators: $(typeof(a)), $(typeof(b)).\nTry to convert to a common operator type first with e.g. dense() or sparse()."))
addnumbererror() = throw(ArgumentError("Can't add or subtract a number and an operator. You probably want 'op + identityoperator(op)*x'."))

length(a::AbstractOperator) = length(a.basis_l)::Int*length(a.basis_r)::Int
basis(a::AbstractOperator) = (check_samebases(a); a.basis_l)

# Ensure scalar broadcasting
Base.broadcastable(x::AbstractOperator) = Ref(x)

# Arithmetic operations
+(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Addition", a, b)
+(a::Number, b::AbstractOperator) = addnumbererror()
+(a::AbstractOperator, b::Number) = addnumbererror()

-(a::AbstractOperator) = arithmetic_unary_error("Negation", a)
-(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Subtraction", a, b)
-(a::Number, b::AbstractOperator) = addnumbererror()
-(a::AbstractOperator, b::Number) = addnumbererror()

*(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Multiplication", a, b)
^(a::AbstractOperator, b::Int) = Base.power_by_squaring(a, b)

dagger(a::AbstractOperator) = arithmetic_unary_error("Hermitian conjugate", a)
Base.adjoint(a::AbstractOperator) = dagger(a)

conj(a::AbstractOperator) = arithmetic_unary_error("Complex conjugate", a)
conj!(a::AbstractOperator) = conj(a::AbstractOperator)

dense(a::AbstractOperator) = arithmetic_unary_error("Conversion to dense", a)

transpose(a::AbstractOperator) = arithmetic_unary_error("Transpose", a)

"""
    ishermitian(op::AbstractOperator)

Check if an operator is Hermitian.
"""
ishermitian(op::AbstractOperator) = arithmetic_unary_error(ishermitian, op)

"""
    tensor(x::AbstractOperator, y::AbstractOperator, z::AbstractOperator...)

Tensor product ``\\hat{x}⊗\\hat{y}⊗\\hat{z}⊗…`` of the given operators.
"""
tensor(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Tensor product", a, b)
tensor(op::AbstractOperator) = op
tensor(operators::AbstractOperator...) = reduce(tensor, operators)

"""
    embed(basis1[, basis2], indices::Vector, operators::Vector)

Tensor product of operators where missing indices are filled up with identity operators.
"""
function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
               indices::Vector, operators::Vector{T}) where T<:AbstractOperator

    @assert sortedindices.check_embed_indices(indices)

    N = length(basis_l.bases)
    @assert length(basis_r.bases) == N
    @assert length(indices) == length(operators)

    # Embed all single-subspace operators.
    idxop_sb = [x for x in zip(indices, operators) if typeof(x[1]) <: Int]
    indices_sb = [x[1] for x in idxop_sb]
    ops_sb = [x[2] for x in idxop_sb]

    for (idxsb, opsb) in zip(indices_sb, ops_sb)
        (opsb.basis_l == basis_l.bases[idxsb]) || throw(bases.IncompatibleBases())
        (opsb.basis_r == basis_r.bases[idxsb]) || throw(bases.IncompatibleBases())
    end

    embed_op = tensor([i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...)

    # Embed all joint-subspace operators.
    idxop_comp = [x for x in zip(indices, operators) if typeof(x[1]) <: Array]
    for (idxs, op) in idxop_comp
        embed_op *= embed(basis_l, basis_r, idxs, op)
    end

    return embed_op
end

"""
    embed(basis1[, basis2], indices::Vector, operators::Vector)

Embed operator acting on a joint Hilbert space where missing indices are filled up with identity operators.
"""
function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
               indices::Vector{Int}, op::T) where T<:AbstractOperator
    N = length(basis_l.bases)
    @assert length(basis_r.bases) == N

    reduce(tensor, basis_l.bases[indices]) == op.basis_l || throw(bases.IncompatibleBases())
    reduce(tensor, basis_r.bases[indices]) == op.basis_r || throw(bases.IncompatibleBases())

    index_order = [idx for idx in 1:length(basis_l.bases) if idx ∉ indices]
    all_operators = AbstractOperator[identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in index_order]

    for idx in indices
        pushfirst!(index_order, idx)
    end
    push!(all_operators, op)

    sortedindices.check_indices(N, indices)

    # Create the operator.
    permuted_op = tensor(all_operators...)

    # Create a copy to fill with correctly-ordered objects (basis and data).
    unpermuted_op = copy(permuted_op)

    # Create the correctly ordered basis.
    unpermuted_op.basis_l = basis_l
    unpermuted_op.basis_r = basis_r

    # Reorient the matrix to act in the correctly ordered basis.
    # Get the dimensions necessary for index permuting.
    dims_l = [b.shape[1] for b in basis_l.bases]
    dims_r = [b.shape[1] for b in basis_r.bases]

    # Get the order of indices to use in the first reshape. Julia indices go in
    # reverse order.
    expand_order = index_order[end:-1:1]
    # Get the dimensions associated with those indices.
    expand_dims_l = dims_l[expand_order]
    expand_dims_r = dims_r[expand_order]

    # Prepare the permutation to the correctly ordered basis.
    perm_order_l = [indexin(idx, expand_order)[1] for idx in 1:length(dims_l)]
    perm_order_r = [indexin(idx, expand_order)[1] for idx in 1:length(dims_r)]

    # Perform the index expansion, the permutation, and the index collapse.
    M = (reshape(permuted_op.data, tuple([expand_dims_l; expand_dims_r]...)) |>
         x -> permutedims(x, [perm_order_l; perm_order_r .+ length(dims_l)]) |>
         x -> sparse(reshape(x, (prod(dims_l), prod(dims_r)))))

    unpermuted_op.data = M

    return unpermuted_op
end
embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Int, op::AbstractOperator) = embed(basis_l, basis_r, Int[index], [op])
embed(basis::CompositeBasis, index::Int, op::AbstractOperator) = embed(basis, basis, Int[index], [op])
embed(basis::CompositeBasis, indices::Vector, operators::Vector{T}) where {T<:AbstractOperator} = embed(basis, basis, indices, operators)
embed(basis::CompositeBasis, indices::Vector{Int}, op::AbstractOperator) = embed(basis, basis, indices, op)

"""
    embed(basis1[, basis2], operators::Dict)

`operators` is a dictionary `Dict{Vector{Int}, AbstractOperator}`. The integer vector
specifies in which subsystems the corresponding operator is defined.
"""
function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
               operators::Dict{Vector{Int}, T}) where T<:AbstractOperator
    @assert length(basis_l.bases) == length(basis_r.bases)
    N = length(basis_l.bases)
    if length(operators) == 0
        return identityoperator(T, basis_l, basis_r)
    end
    indices, operator_list = zip(operators...)
    operator_list = [operator_list...;]
    indices_flat = [indices...;]
    start_indices_flat = [i[1] for i in indices]
    complement_indices_flat = Int[i for i=1:N if i ∉ indices_flat]
    operators_flat = T[]
    if all([minimum(I):maximum(I);]==I for I in indices)
        for i in 1:N
            if i in complement_indices_flat
                push!(operators_flat, identityoperator(T, basis_l.bases[i], basis_r.bases[i]))
            elseif i in start_indices_flat
                push!(operators_flat, operator_list[indexin(i, start_indices_flat)[1]])
            end
        end
        return tensor(operators_flat...)
    else
        complement_operators = [identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in complement_indices_flat]
        op = tensor([operator_list; complement_operators]...)
        perm = sortperm([indices_flat; complement_indices_flat])
        return permutesystems(op, perm)
    end
end
embed(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:AbstractOperator} = embed(basis_l, basis_r, Dict([i]=>op_i for (i, op_i) in operators); kwargs...)
embed(basis::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:AbstractOperator} = embed(basis, basis, operators; kwargs...)
embed(basis::CompositeBasis, operators::Dict{Vector{Int}, T}; kwargs...) where {T<:AbstractOperator} = embed(basis, basis, operators; kwargs...)

"""
    tr(x::AbstractOperator)

Trace of the given operator.
"""
tr(x::AbstractOperator) = arithmetic_unary_error("Trace", x)

ptrace(a::AbstractOperator, index::Vector{Int}) = arithmetic_unary_error("Partial trace", a)

"""
    normalize(op)

Return the normalized operator so that its `tr(op)` is one.
"""
normalize(op::AbstractOperator) = op/tr(op)

"""
    normalize!(op)

In-place normalization of the given operator so that its `tr(x)` is one.
"""
normalize!(op::AbstractOperator) = throw(ArgumentError("normalize! is not defined for this type of operator: $(typeof(op)).\n You may have to fall back to the non-inplace version 'normalize()'."))

"""
    expect(op, state)

Expectation value of the given operator `op` for the specified `state`.

`state` can either be a (density) operator or a ket.
"""
expect(op::AbstractOperator{B,B}, state::Ket{B}) where B<:Basis = state.data' * (op * state).data
expect(op::AbstractOperator{B1,B2}, state::AbstractOperator{B2,B2}) where {B1<:Basis,B2<:Basis} = tr(op*state)

"""
    expect(index, op, state)

If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number.
"""
function expect(indices::Vector{Int}, op::AbstractOperator{B1,B2}, state::AbstractOperator{B3,B3}) where {B1<:Basis,B2<:Basis,B3<:CompositeBasis}
    N = length(state.basis_l.shape)
    indices_ = sortedindices.complement(N, indices)
    expect(op, ptrace(state, indices_))
end
function expect(indices::Vector{Int}, op::AbstractOperator{B,B}, state::Ket{B2}) where {B<:Basis,B2<:CompositeBasis}
    N = length(state.basis.shape)
    indices_ = sortedindices.complement(N, indices)
    expect(op, ptrace(state, indices_))
end
expect(index::Int, op::AbstractOperator, state) = expect([index], op, state)
expect(op::AbstractOperator, states::Vector) = [expect(op, state) for state=states]
expect(indices::Vector{Int}, op::AbstractOperator, states::Vector) = [expect(indices, op, state) for state=states]

"""
    variance(op, state)

Variance of the given operator `op` for the specified `state`.

`state` can either be a (density) operator or a ket.
"""
function variance(op::AbstractOperator{B,B}, state::Ket{B}) where B<:Basis
    x = op*state
    state.data'*(op*x).data - (state.data'*x.data)^2
end
function variance(op::AbstractOperator{B,B}, state::AbstractOperator{B,B}) where B<:Basis
    expect(op*op, state) - expect(op, state)^2
end

"""
    variance(index, op, state)

If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number
"""
function variance(indices::Vector{Int}, op::AbstractOperator{B,B}, state::AbstractOperator{BC,BC}) where {B<:Basis,BC<:CompositeBasis}
    N = length(state.basis_l.shape)
    indices_ = sortedindices.complement(N, indices)
    variance(op, ptrace(state, indices_))
end
function variance(indices::Vector{Int}, op::AbstractOperator{B,B}, state::Ket{BC}) where {B<:Basis,BC<:CompositeBasis}
    N = length(state.basis.shape)
    indices_ = sortedindices.complement(N, indices)
    variance(op, ptrace(state, indices_))
end
variance(index::Int, op::AbstractOperator, state) = variance([index], op, state)
variance(op::AbstractOperator, states::Vector) = [variance(op, state) for state=states]
variance(indices::Vector{Int}, op::AbstractOperator, states::Vector) = [variance(indices, op, state) for state=states]

"""
    exp(op::AbstractOperator)

Operator exponential.
"""
exp(op::AbstractOperator) = throw(ArgumentError("exp() is not defined for this type of operator: $(typeof(op)).\nTry to convert to dense operator first with dense()."))

permutesystems(a::AbstractOperator, perm::Vector{Int}) = arithmetic_unary_error("Permutations of subsystems", a)

"""
    identityoperator(a::Basis[, b::Basis])

Return an identityoperator in the given bases.
"""
identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:AbstractOperator} = throw(ArgumentError("Identity operator not defined for operator type $T."))
identityoperator(::Type{T}, b::Basis) where {T<:AbstractOperator} = identityoperator(T, b, b)
identityoperator(op::T) where {T<:AbstractOperator} = identityoperator(T, op.basis_l, op.basis_r)

one(b::Basis) = identityoperator(b)
one(op::AbstractOperator) = identityoperator(op)

# Fast in-place multiplication
"""
    gemv!(alpha, a, b, beta, result)

Fast in-place multiplication of operators with state vectors. It
implements the relation `result = beta*result + alpha*a*b`.
Here, `alpha` and `beta` are complex numbers, while `result` and either `a`
or `b` are state vectors while the other one can be of any operator type.
"""
gemv!() = error("Not Implemented.")

"""
    gemm!(alpha, a, b, beta, result)

Fast in-place multiplication of of operators with DenseOperators. It
implements the relation `result = beta*result + alpha*a*b`.
Here, `alpha` and `beta` are complex numbers, while `result` and either `a`
or `b` are dense operators while the other one can be of any operator type.
"""
gemm!() = error("Not Implemented.")

# Helper functions to check validity of arguments
function check_ptrace_arguments(a::AbstractOperator, indices::Vector{Int})
    if !isa(a.basis_l, CompositeBasis) || !isa(a.basis_r, CompositeBasis)
        throw(ArgumentError("Partial trace can only be applied onto operators with composite bases."))
    end
    rank = length(a.basis_l.shape)
    if rank != length(a.basis_r.shape)
        throw(ArgumentError("Partial trace can only be applied onto operators wich have the same number of subsystems in the left basis and right basis."))
    end
    if rank == length(indices)
        throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use tr() instead."))
    end
    sortedindices.check_indices(length(a.basis_l.shape), indices)
    for i=indices
        if a.basis_l.shape[i] != a.basis_r.shape[i]
            throw(ArgumentError("Partial trace can only be applied onto subsystems that have the same left and right dimension."))
        end
    end
end
function check_ptrace_arguments(a::StateVector, indices::Vector{Int})
    if length(basis(a).shape) == length(indices)
        throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use tr() instead."))
    end
    sortedindices.check_indices(length(basis(a).shape), indices)
end

samebases(a::AbstractOperator) = samebases(a.basis_l, a.basis_r)::Bool
samebases(a::AbstractOperator, b::AbstractOperator) = samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool
check_samebases(a::AbstractOperator) = check_samebases(a.basis_l, a.basis_r)

multiplicable(a::AbstractOperator, b::Ket) = multiplicable(a.basis_r, b.basis)
multiplicable(a::Bra, b::AbstractOperator) = multiplicable(a.basis, b.basis_l)
multiplicable(a::AbstractOperator, b::AbstractOperator) = multiplicable(a.basis_r, b.basis_l)

end # module


text:
module Fun1d

using LinearAlgebra
using QuadGK
using Random

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

export lincom
"""
Calculate a linear combination
"""
function lincom(x1::S, y1::T, x2::S, y2::T, x::S) where {S,T}
    return T(x - x2) / T(x1 - x2) * y1 + T(x - x1) / T(x2 - x1) * y2
end

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

export Domain
"""
A domain (in the continuum)
"""
struct Domain{S}
    xmin::S
    xmax::S
end

function invariant(domain::Domain)
    return domain.xmin < domain.xmax
end

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

export Grid
"""
A discrete domain, with a particular resolution
"""
struct Grid{S}
    domain::Domain{S}
    ncells::Int
end

function invariant(grid::Grid)
    return invariant(grid.domain) && grid.ncells > 0
end

export location
"""
Grid point location
"""
function location(grid::Grid, n::Int)
    return lincom(1, grid.domain.xmin, grid.ncells + 1, grid.domain.xmax, n)
end

export spacing
"""
Grid spacing
"""
spacing(grid::Grid) = (grid.domain.xmax - grid.domain.xmin) / grid.ncells

export basis
"""
Evaluate the basis function for grid point `n` at point `x`
"""
function basis(grid::Grid{S}, n::Int, x::S) where {S}
    @assert 1 ≤ n ≤ grid.ncells + 1
    q = lincom(grid.domain.xmin, S(1), grid.domain.xmax, S(grid.ncells + 1), x)
    y1 = lincom(S(n - 1), S(0), S(n), S(1), q)
    y2 = lincom(S(n + 1), S(0), S(n), S(1), q)
    y = max(zero(S), min(y1, y2))
    return y
end

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

export GridFun
"""
A grid function, i.e. a discretized function living on a grid
"""
struct GridFun{S,T} <: AbstractArray{T,1}
    grid::Grid{S}
    values::Vector{T}

    function GridFun{S,T}(grid::Grid{S}, values::Vector{T}) where {S,T}
        @assert length(values) == grid.ncells + 1
        return new{S,T}(grid, values)
    end
end

function invariant(f::GridFun)
    return invariant(f.grid) && length(f.values) == f.grid.ncells + 1
end

Base.:(==)(f::GridFun, g::GridFun) = f.grid == g.grid && f.values == g.values

function GridFun(grid::Grid{S}, values::Vector{T}) where {S,T}
    return GridFun{S,T}(grid, values)
end

# Iteration
Base.eltype(::Type{GridFun{S,T}}) where {S,T} = T
Base.isempty(::GridFun) = false
Base.iterate(f::GridFun, state...) = iterate(f.values, state...)
# Base.length(f::GridFun) = length(f.values)

# Indexing
Base.firstindex(f::GridFun) = firstindex(f.values)
Base.getindex(f::GridFun, i) = getindex(f.values, i)
Base.lastindex(f::GridFun) = lastindex(f.values)
Base.setindex!(f::GridFun, v, i) = setindex!(f.values, v, i)

# Abstract Array
Base.size(f::GridFun) = size(f.values)
Base.IndexStyle(::GridFun) = IndexLinear()
Base.similar(f::GridFun) = GridFun(f.grid, similar(f.values))
function Base.similar(f::GridFun, ::Type{T}) where {T}
    return GridFun(f.grid, similar(f.values, T))
end
Base.similar(f::GridFun, ::Dims) = similar(f)
Base.similar(f::GridFun, ::Dims, ::Type{T}) where {T} = similar(f, T)

# Broadcasting
Base.BroadcastStyle(::Type{<:GridFun}) = Broadcast.ArrayStyle{GridFun}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFun}},
                      ::Type{T}) where {T}
    f = find_gridfun(bc)
    return similar(f, T)
end
find_gridfun(bc::Base.Broadcast.Broadcasted) = find_gridfun(bc.args)
find_gridfun(args::Tuple) = find_gridfun(find_gridfun(args[1]), Base.tail(args))
find_gridfun(x) = x
find_gridfun(::Tuple{}) = nothing
find_gridfun(a::GridFun, rest) = a
find_gridfun(::Any, rest) = find_gridfun(rest)

# Others
function Base.map(fun, f::GridFun, gs::GridFun...)
    @assert all(g.grid == f.grid for g in gs)
    return GridFun(f.grid, map(fun, f.values, (g.values for g in gs)...))
end

function Base.rand(rng::AbstractRNG, ::Type{GridFun{S,T}},
                   grid::Grid{S}) where {S,T}
    return GridFun(grid, rand(rng, T, grid.ncells + 1))
end
function Base.rand(::Type{GridFun{S,T}}, grid::Grid{S}) where {S,T}
    return rand(Random.GLOBAL_RNG, GridFun{S,T}, grid)
end

function Base.zero(::Type{GridFun{S,T}}, grid::Grid{S}) where {S,T}
    return GridFun(grid, zeros(T, grid.ncells + 1))
end
Base.zero(f::GridFun{S,T}) where {S,T} = zero(GridFun{S,T}, f.grid)

Base.iszero(f::GridFun) = all(iszero, f.values)

Base.:+(f::GridFun) = GridFun(f.grid, +f.values)
Base.:-(f::GridFun) = GridFun(f.grid, -f.values)

function Base.:+(f::GridFun, g::GridFun)
    @assert f.grid == g.grid
    return GridFun(f.grid, f.values + g.values)
end
function Base.:-(f::GridFun, g::GridFun)
    @assert f.grid == g.grid
    return GridFun(f.grid, f.values - g.values)
end

Base.:*(f::GridFun, a::Number) = GridFun(f.grid, f.values * a)
Base.:*(a::Number, f::GridFun) = GridFun(f.grid, a * f.values)
Base.:/(f::GridFun, a::Number) = GridFun(f.grid, f.values / a)
Base.:\(a::Number, f::GridFun) = GridFun(f.grid, a \ f.values)

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

export sample
"""
Sample a function at grid points
"""
function sample(::Type{T}, grid::Grid{S}, fun) where {S,T}
    values = T[fun(location(grid, n)) for n in 1:(grid.ncells + 1)]
    return GridFun(grid, values)
end

export project
"""
Project a function onto the basis functions
"""
function project(::Type{T}, grid::Grid{S}, fun) where {S,T}
    atol = eps(T)^(T(3) / 4)
    n = grid.ncells + 1
    dx = spacing(grid)

    fbvalues = Array{T}(undef, n)
    for i in 1:n
        x0 = location(grid, i)
        xmin = max(grid.domain.xmin, x0 - dx)
        xmax = min(grid.domain.xmax, x0 + dx)
        I, E = quadgk(x -> fun(x) * basis(grid, i, x), xmin, xmax; atol=atol)
        fbvalues[i] = I
    end

    bbdiag = Array{S}(undef, n)
    bbdiag[1] = dx / 3
    bbdiag[2:(end - 1)] .= 2 * dx / 3
    bbdiag[end] = dx / 3
    bbsdiag = Array{S}(undef, n - 1)
    bbsdiag[1:end] .= dx / 6
    bbmatrix = SymTridiagonal(bbdiag, bbsdiag)

    fvalues = bbmatrix \ fbvalues
    return GridFun(grid, fvalues)
end

export evaluate
"""
Evaluate a grid function at an arbitrary point
"""
function evaluate(f::GridFun{S,T}, x::S) where {S,T}
    q = lincom(f.grid.domain.xmin, S(1), f.grid.domain.xmax,
               S(f.grid.ncells + 1), x)
    @assert q ≥ 1 && q ≤ f.grid.ncells + 1
    n = min(f.grid.ncells, floor(Int, q))
    y1 = f.values[n] * basis(f.grid, n, x)
    y2 = f.values[n + 1] * basis(f.grid, n + 1, x)
    y = y1 + y2
    return y::T
end

# TODO: norm2 (both for continuum and grid functions)

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

export integrate
"""
Integrate a grid function over the whole domain
"""
function integrate(f::GridFun{S,T}) where {S,T}
    dx = spacing(f.grid)
    n = f.grid.ncells + 1
    int = zero(T)
    int += f.values[1] / 2
    for i in 2:(n - 1)
        int += f.values[i]
    end
    int += f.values[n] / 2
    int *= dx
    return int
end

export norm2
"""
L2-norm of a grid function
"""
norm2(f::GridFun) = sqrt(integrate(map(abs2, f)))

#function LinearAlgebra.norm(f::GridFun{S,T})where{S,T}

export deriv
"""
Derivative of a grid function
"""
function deriv(f::GridFun{S,T}) where {S,T}
    dx = spacing(f.grid)
    n = f.grid.ncells + 1
    dvalues = Array{T}(undef, n)
    dvalues[1] = (f.values[2] - f.values[1]) / dx
    for i in 2:(n - 1)
        dvalues[i] = (f.values[i + 1] - f.values[i - 1]) / 2dx
    end
    dvalues[n] = (f.values[n] - f.values[n - 1]) / dx
    return GridFun(f.grid, dvalues)
end

export deriv2
"""
Second derivative of a grid function
"""
function deriv2(f::GridFun{S,T}) where {S,T}
    dx = spacing(f.grid)
    n = f.grid.ncells + 1
    d2values = Array{T}(undef, n)
    for i in 2:(n - 1)
        d2values[i] = (f.values[i - 1] - 2 * f.values[i] + f.values[i + 1]) /
                      dx^2
    end
    d2values[1] = d2values[2]
    d2values[n] = d2values[n - 1]
    return GridFun(f.grid, d2values)
end

export deriv2_sphere
"""
Second derivative of a grid function in spherical coordinates,
including the metric term (∂ᵣᵣ + 2/r ∂ᵣ)
"""
function deriv2_sphere(f::GridFun{S,T}) where {S,T}
    dx = spacing(f.grid)
    n = f.grid.ncells + 1
    d2values = Array{T}(undef, n)
    # Use L'Hôpital's rule at r = 0: limit (r→0) f′(r) / r = f′′(0)
    d2values[1] = 6 * (f.values[2] - f.values[1]) / dx^2
    for i in 2:(n - 1)
        r = dx * (i - 1)
        d2values[i] = (f.values[i - 1] - 2 * f.values[i] + f.values[i + 1]) /
                      dx^2 + 2 / r * (f.values[i + 1] - f.values[i - 1]) / 2dx
    end
    r = dx * (n - 1)
    d2values[n] = (f.values[n - 2] - 2 * f.values[n - 1] + f.values[n]) / dx^2 +
                  2 / r * (f.values[n] - f.values[n - 1]) / dx
    return GridFun(f.grid, d2values)
end

function deriv2_sphere_unused(f::GridFun{S,T}) where {S,T}
    dx = spacing(f.grid)
    n = f.grid.ncells + 1
    # 1/r² ∂_r r² ∂ᵣ
    # 1. cell-centred (forward) derivative:
    #    u[i] := (u[i+1] - u[i]) / dx
    # 2. multiply by metric term:
    #    u[i] := ((i+1/2) dx)^2 * u[i]
    # 3. vertex-centred (backward) derivative:
    #    u[i] := (u[i] - u[i-1]) / dx
    # 4. divide by metric term:
    #    u[i] := 1/(i dx)^2 * u[i]
    d2values = Array{T}(undef, n)
    d2values[1] = 6 * (f.values[2] - f.values[1]) / dx^2
    for i in 2:(n - 1)
        j = i - 1
        d2values[i] = ((1 + (1 - 4 * j) / T(4 * j^2)) * f.values[i - 1] -
                       (2 + 2 / T(4 * j^2)) * f.values[i] +
                       (1 + (1 + 4 * j) / T(4 * j^2)) * f.values[i + 1]) / dx^2
    end
    j = n - 1
    d2values[n] = ((1 + (1 - 2 * j) / T(2 * j^2)) * f.values[n - 2] -
                   (2 + 2 / T(2 * j^2)) * f.values[n - 1] +
                   (1 + (1 + 2 * j) / T(2 * j^2)) * f.values[n]) / dx^2
    return GridFun(f.grid, d2values)
end

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

end


text:
"""
A symmetric arrowhead matrix:

A = [a_1           b_1]
    [    a_2       b_2]
    [       a_3    b_3]
    [           ⋱   ⋮ ]
    [b_1 b_2 b_3 ⋯  c ]
"""
struct SymArrow{T} <: AbstractMatrix{T}
    a::Vector{T}
    b::Vector{T}
    c::T
end

size(A::SymArrow) = length(A.a) + 1, length(A.a) + 1

function getindex(A::SymArrow{T}, i::Integer, j::Integer) where T
    n = size(A, 1)
    if i == j
        if i == n
            return A.c
        else
            return A.a[i]
        end
    elseif i == n
        return A.b[j]
    elseif j == n
        return A.b[i]
    else
        return zero(T)
    end
end

issymmetric(A::SymArrow) = true

+(A::SymArrow, B::SymArrow) = SymArrow(A.a+B.a, A.b+B.b, A.c+B.c)
-(A::SymArrow, B::SymArrow) = SymArrow(A.a-B.a, A.b-B.b, A.c-B.c)

function isincreasing(a::Vector)
    ret = true
    @inbounds for i = 2:length(a)
        ret *= a[i] ≥ a[i-1]
    end
    ret
end

function eig(A::SymArrow{T}) where T
    n = size(A, 1)
    λ = eigvals(A)
    Q = zeros(T, n, n)
    A1 = synthesize_symarrow(A, λ)
    a, b = A1.a, A1.b
    nrmb = norm(b)
    for j = 1:n-1
        λj = λ[j]
        if norm(b[j]) < nrmb*eps(T)
            println("EZ out: ", norm(b[j]) < nrmb*eps(T))
            if norm(a[j] - λj) < norm(a[j] - λ[j+1])
                Q[j,j] = one(T)
            else
                Q[j-1,j] = one(T)
            end
        else
            nrm = zero(T)
            @inbounds @simd for i = 1:n-1
                Q[i,j] = b[i]/(λj-a[i])
                !isfinite(b[i]/(λj-a[i])) && println("j = ",j,", and i = ",i," and ",b[i],"  ",λj,"  ",a[i],"  ",b[i]/(λj-a[i]))
                nrm += abs2(Q[i,j])
            end
            Q[n,j] = one(T)
            nrm += abs2(Q[n,j])
            nrm = sqrt(nrm)
            invnrm = inv(nrm)
            @inbounds @simd for i = 1:n
                Q[i,j] *= invnrm
            end
        end
    end
    λn = λ[n]
    if norm(a[n-1] - λn) < eps(T)
        println("EZ out: ", norm(a[n-1] - λn))
        Q[n-1,n] = one(T)
        #else
        #    Q[n,n] = one(T)
        #end
    else
        nrm = zero(T)
        @inbounds @simd for i = 1:n-1
            Q[i,n] = b[i]/(λn-a[i])
            !isfinite(b[i]/(λn-a[i])) && println("j = ",n,", and i = ",i," and ",b[i],"  ",λn,"  ",a[i],"  ",b[i]/(λn-a[i]))
            nrm += abs2(Q[i,n])
        end
        Q[n,n] = one(T)
        nrm += abs2(Q[n,n])
        nrm = sqrt(nrm)
        invnrm = inv(nrm)
        @inbounds @simd for i = 1:n
            Q[i,n] *= invnrm
        end
    end
    λ, Q
end

"""
Synthesize a symmetric arrow matrix by its eigenvalues and its shaft.
"""
function synthesize_symarrow(A::SymArrow{T}, λ::Vector{T}) where T
    n = size(A, 1)
    a = copy(A.a)
    b = Vector{T}(n-1)
    c = zero(T)
    @inbounds for i = 1:n-1
        c += λ[i] - a[i]
    end
    c += λ[n]

    @inbounds for i = 1:n-1
        ai = a[i]
        t = (ai-λ[1])*(λ[n]-ai)
        for j = 1:i-1
            t *= (λ[j+1]-ai)/(a[j]-ai)
        end
        for j = i+1:n-1
            t *= (λ[j]-ai)/(a[j]-ai)
        end
        b[i] = copysign(sqrt(t), A.b[i])
    end

    SymArrow(a, b, c)
end

function eigvals(A::SymArrow{T}; verbose::Bool = false) where T
    a, b, c = A.a, A.b, A.c
    @assert isincreasing(a)
    n = size(A, 1)
    n2 = n^2
    λ = zeros(T, n)

    nrmb2 = abs(dot(b, b))
    absc = abs(c)

    λk = first_initial_guess(a[1], nrmb2, c)
    δλk = λk
    while abs(δλk) > n2*eps(λk)
        δλk = first_pick_zero_update(A, λk)
        λk += isfinite(δλk) ? δλk : zero(T)
        verbose && println(1," This is the update: ", δλk," and this is the condition: ",n2*eps(λk))
    end
    λ[1] = λk

    for j = 2:n-1
        λk = (a[j]+a[j-1])/2
        δλk = λk
        while abs(δλk) > n2*eps(λk)
            δλk = pick_zero_update(A, j, λk)
            λk += isfinite(δλk) ? δλk : zero(T)
            verbose && println(j," This is the update: ", δλk," and this is the condition: ",n2*eps(λk))
        end
        λ[j] = λk
    end

    λk = last_initial_guess(a[n-1], nrmb2, c)
    δλk = λk
    verbose && println(n," This is the update: ", δλk," and this is the condition: ",n2*eps(λk))
    while abs(δλk) > n2*eps(λk)
        δλk = last_pick_zero_update(A, λk)
        λk += isfinite(δλk) ? δλk : zero(T)
        verbose && println(n," This is the update: ", δλk," and this is the condition: ",n2*eps(λk))
    end
    λ[n] = λk

    return λ
end

function eigvals_FMM(A::SymArrow{T}; verbose::Bool = false) where T
    a, b, c = A.a, A.b, A.c
    @assert isincreasing(a)
    n = size(A, 1)
    n2 = n^2
    λ = zeros(T, n)

    nrmb2 = abs(dot(b, b))
    absc = abs(c)

    λk = first_initial_guess(a[1], nrmb2, c)
    δλk = λk
    while abs(δλk) > n2*eps(λk)
        δλk = first_pick_zero_update(A, λk)
        λk += isfinite(δλk) ? δλk : zero(T)
        verbose && println(1," This is the update: ", δλk," and this is the condition: ",n2*eps(λk))
    end
    λ[1] = λk

    λmidk = (a[2:n-1]+a[1:n-2])/2
    δλk = λmidk
    while norm(δλk) > n2*eps(norm(λmidk))
        δλk = pick_zero_update_FMM(A, λmidk)
        δλk = T[isfinite(δλk) ? δλk : zero(T) for δλk in δλk]
        λmidk .+= δλk
        verbose && println("j This is the update: ", norm(δλk)," and this is the condition: ",n2*eps(norm(λk)))
    end
    λ[2:n-1] .= λmidk

    λk = last_initial_guess(a[n-1], nrmb2, c)
    δλk = λk
    while abs(δλk) > n2*eps(λk)
        δλk = last_pick_zero_update(A, λk)
        λk += isfinite(δλk) ? δλk : zero(T)
        verbose && println(n," This is the update: ", δλk," and this is the condition: ",n2*eps(λk))
    end
    λ[n] = λk

    return λ
end


text:
# Code for video 3_4:

Pkg.add("ODBC")
using ODBC

conn = ODBC.connect("irisDSN")
# ODBC Connection Object
# ----------------------
# Connection Data Source: irisDSN
# irisDSN Connection Number: 1
# Contains resultset(s)? No

results = query("select * from iris_species")
3x2 DataFrames.DataFrame
| Row | id | species         |
|-----|----|-----------------|
| 1   | 1  | "I. setosa"     |
| 2   | 2  | "I. versicolor" |
| 3   | 3  | "I. virginica"  |

julia> results = query("select * from iris")
150x5 DataFrames.DataFrame
| Row | sepal_length | sepal_width | petal_length | petal_width | species |
|-----|--------------|-------------|--------------|-------------|---------|
| 1   | 5.1          | 3.5         | 1.4          | 0.2         | 1       |
| 2   | 4.9          | 3.0         | 1.4          | 0.2         | 1       |
| 3   | 4.7          | 3.2         | 1.3          | 0.2         | 1       |
| 4   | 4.6          | 3.1         | 1.5          | 0.2         | 1       |
...

sql1 = "select sepal_length, petal_length, iris_species.species
        from iris, iris_species
        where iris_species.id = iris.species
        and iris_species.id = 2"
results = query(sql1)
50x3 DataFrames.DataFrame
| Row | sepal_length | petal_length | species         |
|-----|--------------|--------------|-----------------|
| 1   | 7.0          | 4.7          | "I. versicolor" |
| 2   | 6.4          | 4.5          | "I. versicolor" |
| 3   | 6.9          | 4.9          | "I. versicolor" |
| 4   | 5.5          | 4.0          | "I. versicolor" |
| 5   | 6.5          | 4.6          | "I. versicolor" |
| 6   | 5.7          | 4.5          | "I. versicolor" |
| 7   | 6.3          | 4.7          | "I. versicolor" |
| 8   | 4.9          | 3.3          | "I. versicolor" |
| 9   | 6.6          | 4.6          | "I. versicolor" |
| 10  | 5.2          | 3.9          | "I. versicolor" |
| 11  | 5.0          | 3.5          | "I. versicolor" |
⋮
| 39  | 5.6          | 4.1          | "I. versicolor" |
| 40  | 5.5          | 4.0          | "I. versicolor" |
| 41  | 5.5          | 4.4          | "I. versicolor" |
| 42  | 6.1          | 4.6          | "I. versicolor" |
| 43  | 5.8          | 4.0          | "I. versicolor" |
| 44  | 5.0          | 3.3          | "I. versicolor" |
| 45  | 5.6          | 4.2          | "I. versicolor" |
| 46  | 5.7          | 4.2          | "I. versicolor" |
| 47  | 5.7          | 4.2          | "I. versicolor" |
| 48  | 6.2          | 4.3          | "I. versicolor" |
| 49  | 5.1          | 3.0          | "I. versicolor" |
| 50  | 5.7          | 4.1          | "I. versicolor" |

sql1 =  "select sepal_length, petal_length, iris_species.species"
sql1 *= " from iris, iris_species"
sql1 *= " where iris_species.id = iris.species"
sql1 *= " and iris_species.id = 2"

results = query(sql1, conn)

conn.resultset
50x3 DataFrames.DataFrame
| Row | sepal_length | petal_length | species         |
|-----|--------------|--------------|-----------------|
| 1   | 7.0          | 4.7          | "I. versicolor" |
| 2   | 6.4          | 4.5          | "I. versicolor" |
| 3   | 6.9          | 4.9          | "I. versicolor" |

by(results, :sepal_length, size)
21x2 DataFrames.DataFrame
| Row | sepal_length | x1    |
|-----|--------------|-------|
| 1   | 4.9          | (1,3) |
| 2   | 5.0          | (2,3) |
| 3   | 5.1          | (1,3) |
| 4   | 5.2          | (1,3) |
| 5   | 5.4          | (1,3) |
| 6   | 5.5          | (5,3) |
| 7   | 5.6          | (5,3) |
| 8   | 5.7          | (5,3) |
| 9   | 5.8          | (3,3) |
| 10  | 5.9          | (2,3) |
| 11  | 6.0          | (4,3) |
| 12  | 6.1          | (4,3) |
| 13  | 6.2          | (2,3) |
| 14  | 6.3          | (3,3) |
| 15  | 6.4          | (2,3) |
| 16  | 6.5          | (1,3) |
| 17  | 6.6          | (2,3) |
| 18  | 6.7          | (3,3) |
| 19  | 6.8          | (1,3) |
| 20  | 6.9          | (1,3) |
| 21  | 7.0          | (1,3) |

updsql1 = "update iris_species
           set species = 'Iris setosa'
           where id = 1"
query(updsql1)
0x0 DataFrames.DataFrame

disconnect(conn)

### SQLite

# creating the table and importing data:
# F:\Julia\Julia for Data Science\Scripts\Section 3\Video 3.4>sqlite3
# SQLite version 3.8.6 2014-08-15 11:46:33
# Enter ".help" for usage hints.
# Connected to a transient in-memory database.
# Use ".open FILENAME" to reopen on a persistent database.
# sqlite> .separator ","
# sqlite> .import iris.csv iris
# sqlite> .save iris.db
# sqlite> select count(*) from iris;
# 150
# sqlite> .ex

# F:\Julia\Julia for Data Science\Scripts\Section 3\Video 3.4>

Pkg.add("SQLite")
using SQLite

db = SQLite.DB("iris.db")
SQLite.DB("iris.db")

res = query(db, "select * from iris")
# Data.Table:
# 150x5 Data.Schema:
#  sepal_length, sepal_width, petal_length, petal_width,    species
#    UTF8String,  UTF8String,   UTF8String,  UTF8String, UTF8String
# NullableArrays.NullableArray{T,1}[NullableArrays.NullableArray{UTF8String,1}["5.1
# ","4.9","4.7","4.6","5.0","5.4","4.6","5.0","4.4","4.9"  …  "6.7","6.9","5.8","6.
# 8","6.7","6.7","6.3","6.5","6.2","5.9"],NullableArrays.NullableArray{UTF8String,1
# }["3.5","3.0","3.2","3.1","3.6","3.9","3.4","3.4","2.9","3.1"  …  "3.1","3.1","2.
# 7","3.2","3.3","3.0","2.5","3.0","3.4","3.0"],NullableArrays.NullableArray{UTF8St
# ring,1}["1.4","1.4","1.3","1.5","1.4","1.7","1.4","1.5","1.4","1.5"  …  "5.6","5.
# 1","5.1","5.9","5.7","5.2","5.0","5.2","5.4","5.1"],NullableArrays.NullableArray{
# UTF8String,1}["0.2","0.2","0.2","0.2","0.2","0.4","0.3","0.2","0.2","0.1"  …  "2.
# 4","2.3","1.9","2.3","2.5","2.3","1.9","2.0","2.3","1.8"],NullableArrays.Nullable
# Array{UTF8String,1}["I. setosa","I. setosa","I. setosa","I. setosa","I. setosa","
# I. setosa","I. setosa","I. setosa","I. setosa","I. setosa"  …  "I. virginica","I.
#  virginica","I. virginica","I. virginica","I. virginica","I. virginica","I. virgi
# nica","I. virginica","I. virginica","I. virginica"]]

typeof(res) # DataStreams.Data.Table{Array{NullableArrays.NullableArray{T,1},1}}
size(res) # (150,5)

# sqlite> create table iris (sepal_length float, sepal_width float, petal_length fl
# oat, petal_width float, species text);

sql = "select species, avg(petal_width) from iris"
sql *= " group by species"
res = query(db, sql)
nrec = size(res)[1] # number of records returned: 3
for i = 1:nrec
    println(res[i, 1], "\n\t Mean(petal_with): ", res[i, 2])
end
Nullable("I. setosa")
         Mean(petal_with): Nullable(0.2439999999999999)
Nullable("I. versicolor")
         Mean(petal_with): Nullable(1.3259999999999998)
Nullable("I. virginica")
         Mean(petal_with): Nullable(2.026)

max_stars_repo_path: Section 3/3_4.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 3/3_4.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 3/3_4.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
"""
clip_point(r, shp; ngrid = 10)
===
- shp: DataFrame with the column of `geom`
"""
function clip_point(r::GeoArray, shp; ngrid = 10)
    lons, lats = get_coords(shp)
    LON, LAT   = coords_xy(r)
    info       = gdalinfo(r)
    cellsize   = info["cellsize"][1]
    delta      = cellsize/2

    I, J = meshgrid(-ngrid:ngrid, -ngrid:ngrid)
    I = I[:]; J = J[:];
    loc = DataFrame(I = I, J = J)

    res = []
    df  = [];
    @inbounds for i = 1:length(lons)
        lon = lons[i]
        lat = lats[i]
        I_lon = findall( ((LON .- delta) .<= lon) .& ((LON .+ delta) .>= lon) )[1]
        I_lat = findall( ((LAT .- delta) .<= lat) .& ((LAT .+ delta) .>= lat) )[1]

        ind_x = I_lon-ngrid:I_lon+ngrid
        ind_y = I_lat-ngrid:I_lat+ngrid
        arr = @view(r.A[ind_x, ind_y])

        r2 = GeoArray(arr, st_bbox(LON[ind_x], LAT[ind_y]))
        d = raster2df(r2)
        d = cbind(d, loc)

        push!(res, r2)
        push!(df, d)
    end
    Dict("raster" => res, "data" => df)
end

"""
cellIJ(r::GeoArray, point::Tuple)

@return
- I_lon, I_lat
"""
function cellIJ(r::GeoArray, point::Tuple)
    info       = gdalinfo(r)
    LON, LAT   = info["coords"]
    # LON, LAT   = coords_xy(r)
    cellsize   = info["cellsize"][1]
    delta      = cellsize/2

    lon, lat = point
    I_lon = findall( ((LON .- delta) .<= lon) .& ((LON .+ delta) .>= lon) )[1]
    I_lat = findall( ((LAT .- delta) .<= lat) .& ((LAT .+ delta) .>= lat) )[1]
    I_lon, I_lat
end

export cellIJ
export clip_point