text:
# Constructors.
"""
    LinearOperator(Mv, Mtu, Maw, M; symmetric=false, hermitian=false)
Construct a linear operator from a dense or sparse matrix, using `Mv` as storage space
for `M * v` and `Mtu` as storage space for `transpose(M) * u` and `Maw` to store
`adjoint(M) * w`. Use the optional keyword arguments to indicate whether the operator
is symmetric and/or hermitian.
"""
function LinearOperator(
  Mv::AbstractVector{T},
  Mtu::AbstractVector{T},
  Maw::AbstractVector{T},
  M::AbstractMatrix{T};
  symmetric = false,
  hermitian = false,
) where {T}
  nrow, ncol = size(M)
  @assert length(Mv) == nrow
  @assert length(Mtu) == ncol
  @assert length(Maw) == ncol
  prod! = @closure (res, v, α, β) -> mul!(res, M, v, α, β)
  tprod! = @closure (res, u, α, β) -> mul!(res, transpose(M), u, α, β)
  ctprod! = @closure (res, w, α, β) -> mul!(res, adjoint(M), w, α, β)
  LinearOperator{T}(nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, Mv, Mtu, Maw)
end

  """
  LinearOperator(Mv, M :: Symmetric{<:Real})

Construct a linear operator from a symmetric real square matrix `M` with preallocation
using `Mv` as storage space.
"""
LinearOperator(
  Mv::AbstractVector{T},
  M::Union{SymTridiagonal{T}, Symmetric{T}},
  ) where {T <: Real} = LinearOperator(Mv, Mv, Mv, M, symmetric = true, hermitian = true)

  """
    LinearOperator(M)

Constructs a linear operator from a matrix.
"""
function LinearOperator(M::AbstractMatrix{T}; storagetype = Vector{T}, symmetric = false, hermitian = false) where {T}
  nrow, ncol = size(M)
  local Mv, Mtu, Maw
  if T <: Real
    if symmetric
      Maw = Mtu = Mv = storagetype(undef, nrow)
    else
      Mv = storagetype(undef, nrow)
      Maw = Mtu = storagetype(undef, ncol)
    end
  else
    if symmetric && hermitian # Actually real, but T is not
      Maw = Mtu = Mv = storagetype(undef, nrow)
    elseif symmetric
      Mtu = Mv = storagetype(undef, nrow)
      Maw = storagetype(undef, ncol)
    elseif hermitian
      Mv = storagetype(undef, nrow)
      Maw = Mtu = storagetype(undef, ncol)
    else
      Mv = storagetype(undef, nrow)
      Mtu = storagetype(undef, ncol)
      Maw = storagetype(undef, ncol)
    end
  end
  LinearOperator(Mv, Mtu, Maw, M, symmetric = symmetric, hermitian = hermitian)
end

"""
    LinearOperator(M)
Constructs a linear operator from a symmetric tridiagonal matrix. If
its elements are real, it is also Hermitian, otherwise complex
symmetric.
"""
function LinearOperator(M::SymTridiagonal{T}; storagetype = Vector{T}, kwargs...) where {T}
  nrow, ncol = size(M)
  Mv = storagetype(undef, nrow)
  hermitian = eltype(M) <: Real
  Maw = hermitian ? Mv : storagetype(undef, ncol)
  LinearOperator(Mv, Mv, Maw, M, symmetric = true, hermitian = hermitian)
end

"""
    LinearOperator(M)
Constructs a linear operator from a symmetric matrix. If
its elements are real, it is also Hermitian, otherwise complex
symmetric.
"""
function LinearOperator(M::Symmetric{T}; storagetype = Vector{T}, kwargs...) where {T}
  nrow, ncol = size(M)
  Mv = storagetype(undef, nrow)
  hermitian = eltype(M) <: Real
  Maw = hermitian ? Mv : storagetype(undef, ncol)
  LinearOperator(Mv, Mv, Maw, M, symmetric = true, hermitian = hermitian)
end

"""
    LinearOperator(M)
Constructs a linear operator from a Hermitian matrix. If
its elements are real, it is also symmetric.
"""
function LinearOperator(M::Hermitian{T}; storagetype = Vector{T}, kwargs...) where {T}
  nrow, ncol = size(M)
  Mv = storagetype(undef, nrow)
  symmetric = eltype(M) <: Real
  Mtu = symmetric ? Mv : storagetype(undef, ncol)
  LinearOperator(Mv, Mtu, Mv, M, symmetric = symmetric, hermitian = true)
end

# the only advantage of this constructor is optional args
# use LinearOperator{Float64} if you mean real instead of complex
"""
    LinearOperator(nrow, ncol, symmetric, hermitian, prod,
                    [tprod=nothing,
                    ctprod=nothing])
Construct a linear operator from functions.
"""
function LinearOperator(
  nrow::I,
  ncol::I,
  symmetric::Bool,
  hermitian::Bool,
  prod,
  tprod = nothing,
  ctprod = nothing,
) where {I<:Integer}
  T = hermitian ? (symmetric ? Float64 : ComplexF64) : ComplexF64
  Mv = Vector{T}(undef, nrow)
  Mtu = symmetric ? Mv : Vector{T}(undef, ncol)
  Maw = hermitian ? Mv : Vector{T}(undef, ncol)
  LinearOperator{T}(nrow, ncol, symmetric, hermitian, prod, tprod, ctprod, Mv, Mtu, Maw)
end

"""
    LinearOperator(type, nrow, ncol, symmetric, hermitian, prod,
                    [tprod=nothing,
                    ctprod=nothing])
Construct a linear operator from functions where the type is specified as the first argument.
Notice that the linear operator does not enforce the type, so using a wrong type can
result in errors. For instance,
```
A = [im 1.0; 0.0 1.0] # Complex matrix
op = LinearOperator(Float64, 2, 2, false, false, v->A*v, u->transpose(A)*u, w->A'*w)
Matrix(op) # InexactError
```
The error is caused because `Matrix(op)` tries to create a Float64 matrix with the
contents of the complex matrix `A`.
"""
function LinearOperator(
  ::Type{T},
  nrow::I,
  ncol::I,
  symmetric::Bool,
  hermitian::Bool,
  prod,
  tprod = nothing,
  ctprod = nothing,
) where {T, I<:Integer}
  Mv = zeros(T, nrow)
  Mtu = symmetric ? Mv : Vector{T}(undef, ncol)
  Maw = hermitian ? Mv : Vector{T}(undef, ncol)
  LinearOperator{T}(nrow, ncol, symmetric, hermitian, prod, tprod, ctprod, Mv, Mtu, Maw)
end


text:
using RandomNumbers
import Random: seed!
using Test: @test, @test_throws
using Printf: @printf

function compare_dirs(dir1::AbstractString, dir2::AbstractString)
    files1 = readdir(dir1)
    files2 = readdir(dir2)
    @test files1 == files2

    for file in files1
        file1 = joinpath(dir1, file)
        file2 = joinpath(dir2, file)
        lines1 = readlines(file1)
        lines2 = readlines(file2)
        @test lines1 == lines2
    end
end

strip_cr(line::String) = replace(line, r"\r\n$" => "\n")

ext: jl


text:
# updates.jl
#
# Defines the functions that perform the diagonal update, and also
# that build the linked list and operator cluster update

#  (-2,i) is an off-diagonal site operator h(sigma^+_i + sigma^-_i)
#  (-1,i) is a diagonal site operator h
#  (0,0) is the identity operator I - NOT USED IN THE PROJECTOR CASE
#  (i,j) is a diagonal bond operator J(sigma^z_i sigma^z_j)
@inline isdiagonal(op::NTuple{2,Int}) = (op[1] != -2)
@inline isidentity(op::NTuple{2,Int}) = (op[1] == 0)
@inline issiteoperator(op::NTuple{2,Int}) = (op[1] < 0)
@inline isbondoperator(op::NTuple{2,Int}) = (op[1] > 0)

function mc_step!(f::Function, qmc_state::BinaryGroundState, H::TFIM)
    diagonal_update!(qmc_state, H)

    cluster_data = linked_list_update(qmc_state, H)

    f(cluster_data, qmc_state, H)

    cluster_update!(cluster_data, qmc_state, H)
end

mc_step!(qmc_state, H) = mc_step!((args...) -> nothing, qmc_state, H)

########################## finite-beta #######################################

function mc_step_beta!(f::Function, qmc_state::BinaryThermalState, H::TFIM, beta::Real; eq::Bool = false)
    num_ops = diagonal_update_beta!(qmc_state, H, beta; eq = eq)

    cluster_data = linked_list_update_beta(qmc_state, H)

    f(cluster_data, qmc_state, H)

    cluster_update_beta!(cluster_data, qmc_state, H)

    return num_ops
end

mc_step_beta!(qmc_state, H, beta; eq = false) = mc_step_beta!((args...) -> nothing, qmc_state, H, beta; eq = eq)

# returns true if operator insertion succeeded
function insert_diagonal_operator!(qmc_state::BinaryQMCState{N}, H::TFIM{N}, spin_prop::BitArray{N}, n::Int) where N
    site1, site2 = op = rand(H.op_sampler)
    @inbounds if issiteoperator(op) || spin_prop[site1] == spin_prop[site2]
        qmc_state.operator_list[n] = op
        return true
    else
        return false
    end
end

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

function diagonal_update!(qmc_state::BinaryGroundState{N}, H::TFIM{N}) where N
    spin_prop = copyto!(qmc_state.propagated_config, qmc_state.left_config)  # the propagated spin state

    for (n, op) in enumerate(qmc_state.operator_list)
        if !isdiagonal(op)
            @inbounds spin_prop[op[2]] ⊻= 1  # spinflip
        else
            success = false
            while !success
                success = insert_diagonal_operator!(qmc_state, H, spin_prop, n)
            end
        end
    end

    # DEBUG
    # if spin_prop != qmc_state.right_config  # check the spin propagation for error
    #     error("Basis state propagation error in diagonal update!")
    # end
end

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

function linked_list_update(qmc_state::BinaryGroundState{N}, H::TFIM{N}) where N
    Ns = nspins(H)
    spin_left = qmc_state.left_config

    len = 2 * Ns
    @simd for op in qmc_state.operator_list
        len += ifelse(issiteoperator(op), 2, 4)
    end

    # retrieve linked list data structures
    LinkList = qmc_state.linked_list  # needed for cluster update
    LegType = qmc_state.leg_types

    # A diagonal bond operator has non trivial associates for cluster building
    Associates = qmc_state.associates

    First = qmc_state.first

    # The first N elements of the linked list are the spins of the LHS basis state
    @inbounds @simd for i in 1:Ns
        LegType[i] = spin_left[i]
        First[i] = i
        Associates[i] = (0, 0, 0)
    end

    spin_prop = copyto!(qmc_state.propagated_config, spin_left)  # the propagated spin state
    idx = Ns

    # Now, add the 2M operators to the linked list. Each has either 2 or 4 legs
    @inbounds for op in qmc_state.operator_list
        if issiteoperator(op)
            site = op[2]
            # lower or left leg
            idx += 1
            LinkList[idx] = First[site]
            LegType[idx] = spin_prop[site]
            current_link = idx

            if !isdiagonal(op)  # off-diagonal site operator
                spin_prop[site] ⊻= 1  # spinflip
            end

            LinkList[First[site]] = current_link  # completes backwards link
            First[site] = current_link + 1
            Associates[idx] = (0, 0, 0)

            # upper or right leg
            idx += 1
            LegType[idx] = spin_prop[site]
            Associates[idx] = (0, 0, 0)
        else  # diagonal bond operator
            site1, site2 = op

            # lower left
            idx += 1
            LinkList[idx] = First[site1]
            LegType[idx] = spin_prop[site1]
            current_link = idx

            LinkList[First[site1]] = current_link  # completes backwards link
            First[site1] = current_link + 2
            vertex1 = current_link
            Associates[idx] = (vertex1 + 1, vertex1 + 2, vertex1 + 3)

            # lower right
            idx += 1
            LinkList[idx] = First[site2]
            LegType[idx] = spin_prop[site2]
            current_link = idx

            LinkList[First[site2]] = current_link  # completes backwards link
            First[site2] = current_link + 2
            Associates[idx] = (vertex1, vertex1 + 2, vertex1 + 3)

            # upper left
            idx += 1
            LegType[idx] = spin_prop[site1]
            Associates[idx] = (vertex1, vertex1 + 1, vertex1 + 3)

            # upper right
            idx += 1
            LegType[idx] = spin_prop[site2]
            Associates[idx] = (vertex1, vertex1 + 1, vertex1 + 2)
        end
    end

    # The last N elements of the linked list are the final spin state
    @inbounds for i in 1:Ns
        idx += 1
        LinkList[idx] = First[i]
        LegType[idx] = spin_prop[i]
        LinkList[First[i]] = idx
        Associates[idx] = (0, 0, 0)
    end

    # DEBUG
    # if spin_prop != spin_right
    #     @debug "Basis state propagation error: LINKED LIST"
    # end

    return len
end

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

function cluster_update!(lsize::Int, qmc_state::BinaryGroundState{N}, H::TFIM{N}) where N
    Ns = nspins(H)
    spin_left, spin_right = qmc_state.left_config, qmc_state.right_config
    operator_list = qmc_state.operator_list

    LinkList = qmc_state.linked_list
    LegType = qmc_state.leg_types
    Associates = qmc_state.associates

    in_cluster = zeros(Int, lsize)
    cstack = Stack{Int}()  # This is the stack of vertices in a cluster
    ccount = 0  # cluster number counter

    @inbounds for i in 1:lsize
        # Add a new leg onto the cluster
        if (in_cluster[i] == 0 && Associates[i] === (0, 0, 0))
            ccount += 1
            push!(cstack, i)
            in_cluster[i] = ccount

            flip = rand(Bool)  # flip a coin for the SW cluster flip
            if flip
                LegType[i] ⊻= 1  # spinflip
            end

            while !isempty(cstack)
                leg = LinkList[pop!(cstack)]

                if in_cluster[leg] == 0
                    in_cluster[leg] = ccount  # add the new leg and flip it
                    if flip
                        LegType[leg] ⊻= 1
                    end

                    # now check all associates and add to cluster
                    assoc = Associates[leg]  # a 3-tuple
                    if assoc !== (0, 0, 0)
                        for a in assoc
                            push!(cstack, a)
                            in_cluster[a] = ccount
                            if flip
                                LegType[a] ⊻= 1
                            end
                        end
                    end
                end

            end
        end
    end

    # map back basis states and operator list
    @inbounds for i in 1:Ns
        spin_left[i] = LegType[i]  # left basis state
        spin_right[i] = LegType[lsize-Ns+i]  # right basis state
    end

    ocount = Ns + 1  # next on is leg Ns + 1
    @inbounds for (n, op) in enumerate(operator_list)
        if isbondoperator(op)
            ocount += 4
        else
            if LegType[ocount] == LegType[ocount+1]  # diagonal
                operator_list[n] = (-1, op[2])
            else  # off-diagonal
                operator_list[n] = (-2, op[2])
            end
            ocount += 2
        end
    end

end

#############################################################################
#############  FINITE BETA FUNCTIONS BELOW ##################################
#############################################################################

function resize_op_list!(qmc_state::BinaryThermalState, new_size::Int)
    operator_list = filter!(!isidentity, qmc_state.operator_list)
    len = length(operator_list)

    if len < new_size
        tail = init_op_list(new_size - len)
        append!(operator_list, tail)
    end

    len = 4*length(operator_list)
    # these are going to be overwritten by linked_list_update_beta
    # which will be called right after the diagonal update that called this function
    resize!(qmc_state.linked_list, len)
    resize!(qmc_state.leg_types, len)
    resize!(qmc_state.associates, len)
end

function diagonal_update_beta!(qmc_state::BinaryThermalState, H::TFIM, beta::Real; eq::Bool = false)

    # define the Metropolis probability as a constant
    # https://pitp.phas.ubc.ca/confs/sherbrooke2012/archives/Melko_SSEQMC.pdf
    # equation 1.42
    P_norm = beta * H.P_normalization

    num_ids = count(isidentity, qmc_state.operator_list)
    P_remove = (num_ids + 1) / P_norm
    P_accept = P_norm / num_ids

    spin_prop = copyto!(qmc_state.propagated_config, qmc_state.left_config)  # the propagated spin state

    @inbounds for (n, op) in enumerate(qmc_state.operator_list)
        if !isdiagonal(op)
            spin_prop[op[2]] ⊻= 1  # spinflip
        elseif !isidentity(op)
            if rand() < P_remove
                qmc_state.operator_list[n] = (0, 0)
                num_ids += 1
                P_remove = (num_ids + 1) / P_norm
                P_accept = P_norm / num_ids
            end
        else
            if rand() < P_accept
                success = insert_diagonal_operator!(qmc_state, H, spin_prop, n)

                if success
                    # save one operation lol
                    P_remove = num_ids / P_norm
                    num_ids -= 1
                    P_accept = P_norm / num_ids
                end
            end
        end
    end

    # DEBUG
    # if spin_prop != qmc_state.right_config  # check the spin propagation for error
    #     error("Basis state propagation error in diagonal update!")
    # end

    total_list_size = length(qmc_state.operator_list)
    num_ops = total_list_size - num_ids

    if eq && 1.2*num_ops > total_list_size
        resize_op_list!(qmc_state, round(Int, 1.5*num_ops, RoundUp))
    end

    return num_ops
end

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

function linked_list_update_beta(qmc_state::BinaryThermalState, H::TFIM)
    Ns = nspins(H)
    spin_left, spin_right = qmc_state.left_config, qmc_state.right_config

    len = 0
    @simd for op in qmc_state.operator_list
        if issiteoperator(op)
            len += 2
        elseif isbondoperator(op)
            len += 4
        end
    end

    # initialize linked list data structures
    LinkList = qmc_state.linked_list  # needed for cluster update
    LegType = qmc_state.leg_types

    # A diagonal bond operator has non trivial associates for cluster building
    Associates = qmc_state.associates

    First = fill!(qmc_state.first, 0)  #initialize the First list
    Last = fill!(qmc_state.last, 0)   #initialize the Last list
    idx = 0

    spin_prop = copyto!(qmc_state.propagated_config, spin_left)  # the propagated spin state

    # Now, add the 2M operators to the linked list. Each has either 2 or 4 legs
    @inbounds for op in qmc_state.operator_list
        if issiteoperator(op)
            site = op[2]
            # lower or left leg
            idx += 1
            LinkList[idx] = First[site]
            LegType[idx] = spin_prop[site]
            current_link = idx

            if !isdiagonal(op)  # off-diagonal site operator
                spin_prop[site] ⊻= 1  # spinflip
            end

            if First[site] != 0
                LinkList[First[site]] = current_link  # completes backwards link
            else
                Last[site] = current_link
            end
            First[site] = current_link + 1
            Associates[idx] = (0, 0, 0)

            # upper or right leg
            idx += 1
            LegType[idx] = spin_prop[site]
            Associates[idx] = (0, 0, 0)
        elseif isbondoperator(op)  # diagonal bond operator
            site1, site2 = op

            # lower left
            idx += 1
            LinkList[idx] = First[site1]
            LegType[idx] = spin_prop[site1]
            current_link = idx

            if First[site1] != 0
                LinkList[First[site1]] = current_link  # completes backwards link
            else
                Last[site1] = current_link
            end

            First[site1] = current_link + 2
            vertex1 = current_link
            Associates[idx] = (vertex1 + 1, vertex1 + 2, vertex1 + 3)

            # lower right
            idx += 1
            LinkList[idx] = First[site2]
            LegType[idx] = spin_prop[site2]
            current_link = idx

            if First[site2] != 0
                LinkList[First[site2]] = current_link  # completes backwards link
            else
                Last[site2] = current_link
            end

            First[site2] = current_link + 2
            Associates[idx] = (vertex1, vertex1 + 2, vertex1 + 3)

            # upper left
            idx += 1
            LegType[idx] = spin_prop[site1]
            Associates[idx] = (vertex1, vertex1 + 1, vertex1 + 3)

            # upper right
            idx += 1
            LegType[idx] = spin_prop[site2]
            Associates[idx] = (vertex1, vertex1 + 1, vertex1 + 2)
        end
    end

    #Periodic boundary conditions for finite-beta
    @inbounds for i in 1:Ns
        if First[i] != 0  #This might be encountered at high temperatures
            LinkList[First[i]] = Last[i]
            LinkList[Last[i]] = First[i]
        end
    end

    # DEBUG
    # if spin_prop != spin_right
    #     @debug "Basis state propagation error: LINKED LIST"
    # end

    return len
end

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

function cluster_update_beta!(lsize::Int, qmc_state::BinaryThermalState, H::TFIM)
    Ns = nspins(H)
    spin_left, spin_right = qmc_state.left_config, qmc_state.right_config
    operator_list = qmc_state.operator_list

    LinkList = qmc_state.linked_list
    LegType = qmc_state.leg_types
    Associates = qmc_state.associates

    in_cluster = zeros(Int, lsize)
    cstack = Stack{Int}()  # This is the stack of vertices in a cluster
    ccount = 0  # cluster number counter

    @inbounds for i in 1:lsize
        # Add a new leg onto the cluster
        if (in_cluster[i] == 0 && Associates[i] === (0, 0, 0))
            ccount += 1
            push!(cstack, i)
            in_cluster[i] = ccount

            flip = rand(Bool)  # flip a coin for the SW cluster flip
            if flip
                LegType[i] ⊻= 1  # spinflip
            end

            while !isempty(cstack)
                leg = LinkList[pop!(cstack)]

                if in_cluster[leg] == 0
                    in_cluster[leg] = ccount  # add the new leg and flip it
                    if flip
                        LegType[leg] ⊻= 1
                    end

                    # now check all associates and add to cluster
                    assoc = Associates[leg]  # a 3-tuple
                    if assoc !== (0, 0, 0)
                        for a in assoc
                            push!(cstack, a)
                            in_cluster[a] = ccount
                            if flip
                                LegType[a] ⊻= 1
                            end
                        end
                    end
                end

            end
        end
    end

    #println(in_cluster)

    # map back basis states and operator list
    First = qmc_state.first
    Last = qmc_state.last
    @inbounds for i in 1:Ns
        if First[i] != 0
            spin_left[i] = LegType[Last[i]]  # left basis state
            spin_right[i] = LegType[First[i]]  # right basis state
        else
            #randomly flip spins not connected to operators
            spin_left[i] = spin_right[i] = rand(Bool)
        end

    end

    ocount = 1  # first leg
    @inbounds for (n, op) in enumerate(operator_list)
        if isbondoperator(op)
            ocount += 4
        elseif !isidentity(op)
            if LegType[ocount] == LegType[ocount+1]  # diagonal
                operator_list[n] = (-1, op[2])
            else  # off-diagonal
                operator_list[n] = (-2, op[2])
            end
            ocount += 2
        end
    end

end


text:
using LaTeXStrings

using Plots
using DFWannier

T=Float64
cell_param = 4.3392*T[ 0.8424444 -0.4863855 -3.7e-8;-1.5e-8 0.972771 6.0e-9;-5.7e-8 6.0e-9 1.5226926]
atom_pos   = [cell_param*[1/3 ,2/3, 0.3194],cell_param*[2/3,1/3,0.7363],cell_param*[0.0,0.0,0.0]]
# atoms = [[PhysAtom(Point3D(atom_pos[1]),0.210752) for i=1:4];[PhysAtom(Point3D(atom_pos[2]),0.296619) for i=1:4];[PhysAtom(Point3D(atom_pos[3]), -0.770879) for i=1:4]]
atoms = [[PhysAtom(Point3D(atom_pos[1]),-0.0462138) for i=1:4];[PhysAtom(Point3D(atom_pos[2]),0.355636) for i=1:4];[PhysAtom(Point3D(atom_pos[3]), -0.768408) for i=1:4]]
BiTeI = WannierModel{T}("/home/ponet/Documents/PhD/BiTeI/NSOC/","/home/ponet/Documents/PhD/BiTeI/NSOC/2BiTeI_bands.out",atoms);
# x = WannierModel{T}("/Users/ponet/Documents/Fysica/PhD/BiTeI/NSOC/","/Users/ponet/Documents/Fysica/PhD/BiTeI/NSOC/2BiTeI_bands.out",atoms);

calculate_angmom(BiTeI.wfcs[6],BiTeI.wfcs[7])
calculate_angmom(BiTeI.wfcs[6],BiTeI.wfcs[8])

tbbands_soc_tot = calculate_eig_soc(x)
df_bands_soc = read_qe_bands_file("/home/ponet/Documents/PhD/BiTeI/SOC/2BiTeI_bands.out",T)
plot(df_bands_soc,tbbands_soc_tot,ylim=[5,12])
plot(df_bands_soc,ylim=[-5,5],fermi=8.1176)
minimum(df_bands_soc[39].eigvals)-maximum(df_bands_soc[38].eigvals)
for (i,wfc) in enumerate(x.wfcs)
  wfc.atom = atoms[i]
end

tbbands_tot = calculate_eig_cm(x);
df_bands = read_qe_bands_file("/home/ponet/Documents/PhD/BiTeI/NSOC/2BiTeI_bands.out",T)
tbbands    = calculate_eig_cm_angmom(x,90:0.2:110);
tbbandssoc = calculate_eig_cm_angmom_soc(x,test_poitns);
tbbandssoc = calculate_eig_cm_angmom_soc(x,90:0.2:110);

plot(tbbandssoc)
)plot(df_bands_soc[21:44]
#----------all the plots--------
plot(plot([tbbandssoc[9:10]...,tbbands[5]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(2.5:0.2:4)),plot([tbbandssoc[9:10]...,tbbands[5]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot([tbbandssoc[9:10]...,tbbands[5]],:angmom1_x,label = ["SOC" "SOC" "No SOC"],leg=false,title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(-1:0.2:0.2)),plot([tbbandssoc[9:10]...,tbbands[5]],:spin1_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[9:10]...,tbbands[5]],:angmom1_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title="",yticks=collect(-0.2:0.2:0.6)),plot([tbbandssoc[9:10]...,tbbands[5]],:spin1_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
plot(plot([tbbandssoc[5:6]...,tbbands[3]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(3.4:0.2:5)),plot([tbbandssoc[5:6]...,tbbands[3]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)",yticks=collect(-3:0.2:-1)), plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_x,leg=false,label = ["SOC" "SOC" "No SOC"],title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(-0.5:0.2:1)),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
plot(plot([tbbandssoc[7:8]...,tbbands[4]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[7:8]...,tbbands[4]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot([tbbandssoc[7:8]...,tbbands[4]],:angmom2_x,leg=false,label = ["SOC" "SOC" "No SOC"],title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[7:8]...,tbbands[4]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[7:8]...,tbbands[4]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),plot([tbbandssoc[7:8]...,tbbands[4]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
#-----------tbbands vs dfbands-----------------#

dfbands = read_bands_file("/Users/ponet/Documents/Fysica/PhD/BiTeI/NSOC/02BiTeI_bands.out",Float64);
plot(dfbands[8:20],tbbands,legend=false)
outputs = pull_outputs(BiTeI_soc)

BiTeI_nsoc = load_server_job("ponet@10.255.9.115","BiTeI/NSOC","/home/ponet/Documents/PhD/BiTeI/NSOC")
BiTeI_soc = load_server_job("ponet@10.255.9.115","BiTeI/SOC","/home/ponet/Documents/PhD/BiTeI/SOC")
 print_flow(BiTeI_nsoc)
add_calculation!(BiTeI_nsoc,default_inputs[:projwfc],4)
#  BiTeI_nsoc = load_job(phd_dir*"BiTeI/NSOC",server = default_server,server_dir = "BiTeI/NSOC")
 BiTeI_soc = load_server_job("ponet@10.255.9.115","BiTeI/SOC","/home/ponet/Documents/PhD/BiTeI/SOC")
atomic_positions = get_data(BiTeI_soc,"bands",:atomic_positions)
Te_pos = atomic_positions[:Te]
I_pos  = atomic_positions[]
atomic_positions[:Te] = I_pos
atomic_positions[:I] = Te_pos
atomic_positions[:I] = [Point3D{Float32}(1/3,2/3,0.3076)]
atomic_positions[:Te] = [Point3D{Float32}(2/3,1/3,0.7482)]
change_data!(BiTeI_soc, ["bands","scf"],:atomic_positions, atomic_positions)
cell_param = Float32[0.842444445  -0.486385500  -0.000000037;-0.000000015   0.972770975   0.000000006;-0.000000057   0.000000006   1.522692506]
change_flags!(BiTeI_nsoc,Dict(:dis_win_min => 0.f0,:dis_win_max => 10.f0, :dis_froz_min => 0.f0,:dis_froz_max => 10.f0))
outputs = pull_outputs(BiTeI_nsoc,extras=["*.xsf","*r.dat"])

bands = read_qe_bands_file(BiTeI_nsoc.local_dir*"2BiTeI_bands.out")
using Plots

change_flow!(BiTeI_soc,[("vc",false),("nscf",false),("bands",true)])
fermi = read_fermi_from_qe_file(BiTeI_nsoc.local_dir*"1BiTeI_scf.out")
set_flags!(job_nsoc,:projwfc,Dict(:kresolveddos=>true))
set_flags!(job_nsoc,:system,Dict(:occupations => "'smearing'",:smearing => "'cold'",:degauss => 0.01))
change_flow!(job_nsoc,[(5,true),(6,false),(7,false)])
print_flag(job_nsoc,:A)
print_flow(job
plot(bands,ylim=[-10,15])
submit_job(BiTeI_nsoc)

zmat_Bi,ticks = read_qe_kpdos(BiTeI_nsoc.local_dir*"pwo.pdos_atm#3(Bi)_wfc#2(p)")
zmat_Te,ticks = read_qe_kpdos(BiTeI_nsoc.local_dir*"pwo.pdos_atm#2(Te)_wfc#2(p)")
zmat_I,ticks = read_qe_kpdos(BiTeI_nsoc.local_dir*"pwo.pdos_atm#1(I)_wfc#2(p)")
plot(heatmap(zmat_Bi,title = "Bi_p"), heatmap(zmat_Te, title="Te_p"), heatmap(zmat_I,title = "I_p"),ticks=ticks)


text:
module EuclideanDistanceMatrices

using SCS
using Convex
using Random
using Distances
using LinearAlgebra
using Distributions
using TripletEmbeddings

import Base: size, getindex, *, ndims, show
import LinearAlgebra: rank

export EDM, MaskedEDM, NoisyEDM, NoisyMaskedEDM,
       J,

       # EDMs.jl
       masked, isedm, nitems,

       # utilities.jl
       relativeerror, randommask, idxmask,

       # completion.jl
       complete,
       AlternatingDescent, SDP # losses

include("EDMs.jl")
include("masks.jl")
include("triplets.jl")
include("utilities.jl")
include("completion.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 63ba08cc-59a8-11eb-0a0f-27efac60d779
using Pkg, DrWatson

# ╔═╡ 6db218c6-59a8-11eb-2a8b-7107354cf590
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ 51fc19b8-59a8-11eb-2214-15aca59b807b
md" ## Clip-08-01-06s.jl"

# ╔═╡ 8aaa4bcc-59a8-11eb-2003-f1213b116565
begin
	df = CSV.read(sr_datadir("rugged.csv"), DataFrame)
	dropmissing!(df, :rgdppc_2000)
	dropmissing!(df, :rugged)
	df.log_gdp = log.(df[:, :rgdppc_2000])
	df.log_gdp_s = df.log_gdp / mean(df.log_gdp)
	df.rugged_s = df.rugged / maximum(df.rugged)
	data = (N = size(df, 1), G = df.log_gdp_s,
		R = df.rugged_s)
	PRECIS(df[:, [:log_gdp, :log_gdp_s, :rugged, :rugged_s]])
end

# ╔═╡ d7d3e626-6a45-11eb-1820-a3ec98e556b1
stan8_1_1 = "
parameters {
	real a;
	real b;
}

model {
	a ~ normal(1, 1);
	b ~ normal(0, 1);
}
";

# ╔═╡ a3e7c070-6a46-11eb-072b-3943db854020
begin
	m8_1_1s = SampleModel("m8.1.1s", stan8_1_1)
	rc8_1_1s = stan_sample(m8_1_1s)
	if success(rc8_1_1s)
		post8_1_1s_df = read_samples(m8_1_1s, :dataframe)
		PRECIS(post8_1_1s_df[:, [:a, :b]])
	end
end

# ╔═╡ ac6a235c-6a67-11eb-0d6c-c38eb2162598
begin
	x = 0:0.01:1
	p1 = plot(;ylim=(0.6, 1.4), leg=false)
	for i in 1:20:4000
		y = post8_1_1s_df.a[i] .+ post8_1_1s_df.b[i] .* x
		plot!(x, y, color=:lightgrey)
	end
	plot!(x, mean(post8_1_1s_df.a) .+ mean(post8_1_1s_df.b) .* x, color=:darkblue)
end;

# ╔═╡ 7f74d204-6a6a-11eb-3931-4de1dfb8d527
stan8_1_2 = "
parameters {
	real a;
	real b;
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.1);
}
";

# ╔═╡ 8f79ecfc-6a6a-11eb-1173-2fa6e4609d72
begin
	m8_1_2s = SampleModel("m8.1.2s", stan8_1_2)
	rc8_1_2s = stan_sample(m8_1_2s)
	if success(rc8_1_2s)
		post8_1_2s_df = read_samples(m8_1_2s, :dataframe)
		PRECIS(post8_1_2s_df[:, [:a, :b]])
	end
end

# ╔═╡ baaabbf4-6a6a-11eb-091d-cd8a0e80004f
begin
	p2 = plot(;ylim=(0.6, 1.4), leg=false)
	for i in 1:80:4000
		y = post8_1_2s_df.a[i] .+ post8_1_2s_df.b[i] .* x
		plot!(x, y, color=:lightgrey)
	end
	plot!(x, mean(post8_1_2s_df.a) .+ mean(post8_1_2s_df.b) .* x, color=:darkblue)

	plot(p1, p2, layout=(1, 2))
end

# ╔═╡ 9a39b1c0-6a6c-11eb-05ce-ab82a8d58dfe
stan8_1_3 = "
parameters {
	real a;
	real b;
real sigma;
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
sigma ~ exponential(1);
}
";

# ╔═╡ 6269a87e-6a6c-11eb-2e14-39df73ddd6d8
begin
	m8_1_3s = SampleModel("m8.1.3s", stan8_1_3)
	rc8_1_3s = stan_sample(m8_1_3s; data)
	if success(rc8_1_3s)
		post8_1_3s_df = read_samples(m8_1_3s, :dataframe)
		PRECIS(post8_1_3s_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ 6b462290-6a6c-11eb-2d9b-97241415b825
begin
	plot(;ylim=(0.6, 1.4), leg=false)
	for i in 1:80:4000
		y = post8_1_3s_df.a[i] .+ post8_1_3s_df.b[i] .* x
		plot!(x, y, color=:lightgrey)
	end
	plot!(x, mean(post8_1_3s_df.a) .+ mean(post8_1_3s_df.b) .* x, color=:darkblue)
end

# ╔═╡ e890c59c-6a68-11eb-0d33-4d52a21d0ccf
stan8_1_4 = "
data {
	int N;
	vector[N] G;
	vector[N] R;
}

parameters {
	real a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	mu = a + b * (R - 0.125);
}

model {
	a ~ normal(1, 1);
	b ~ normal(0, 1);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
";

# ╔═╡ 8e815d5c-6a66-11eb-21dc-99734c31f4e1
begin
	m8_1_4s = SampleModel("m8.1.4s", stan8_1_4)
	rc8_1_4s = stan_sample(m8_1_4s; data)
	if success(rc8_1_4s)
		post8_1_4s_df = read_samples(m8_1_4s, :dataframe)
		PRECIS(post8_1_4s_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ 04a39568-6a6c-11eb-19ea-f55be3914a51
begin
	plot(;ylim=(0.6, 1.4), leg=false)
	for i in 1:80:4000
		y = post8_1_4s_df.a[i] .+ post8_1_4s_df.b[i] .* x
		plot!(x, y, color=:lightgrey)
	end
	plot!(x, mean(post8_1_4s_df.a) .+ mean(post8_1_4s_df.b) .* x, color=:darkblue)
end

# ╔═╡ 45767e2e-6a63-11eb-3e12-354a7e32a374
md" ## End of clip-08-01-06s.jl"

# ╔═╡ Cell order:
# ╟─51fc19b8-59a8-11eb-2214-15aca59b807b
# ╠═63ba08cc-59a8-11eb-0a0f-27efac60d779
# ╠═6db218c6-59a8-11eb-2a8b-7107354cf590
# ╠═8aaa4bcc-59a8-11eb-2003-f1213b116565
# ╠═d7d3e626-6a45-11eb-1820-a3ec98e556b1
# ╠═a3e7c070-6a46-11eb-072b-3943db854020
# ╠═ac6a235c-6a67-11eb-0d6c-c38eb2162598
# ╠═7f74d204-6a6a-11eb-3931-4de1dfb8d527
# ╠═8f79ecfc-6a6a-11eb-1173-2fa6e4609d72
# ╠═baaabbf4-6a6a-11eb-091d-cd8a0e80004f
# ╠═9a39b1c0-6a6c-11eb-05ce-ab82a8d58dfe
# ╠═6269a87e-6a6c-11eb-2e14-39df73ddd6d8
# ╠═6b462290-6a6c-11eb-2d9b-97241415b825
# ╠═e890c59c-6a68-11eb-0d33-4d52a21d0ccf
# ╠═8e815d5c-6a66-11eb-21dc-99734c31f4e1
# ╠═04a39568-6a6c-11eb-19ea-f55be3914a51
# ╟─45767e2e-6a63-11eb-3e12-354a7e32a374


text:
### Sampleability

abstract type Sampleability end

mutable struct Deterministic <: Sampleability end

mutable struct Random <: Sampleability end

### Abstract variable

abstract type Variable{S<:Sampleability} end

VariableVector{V<:Variable} = Vector{V}

eltype{S<:Sampleability}(::Type{Variable{S}}) = S

vertex_key(v::Variable) = v.key
vertex_index(v::Variable) = v.index
is_indexed(v::Variable) = v.index > 0 ? true : false

keys(variables::VariableVector) = Symbol[v.key for v in variables]
indexes(variables::VariableVector) = Integer[v.index for v in variables]

sort_by_index(vs::VariableVector) = vs[[v.index for v in vs]]

function default_state(v::Variable, v0, outopts::Dict=Dict())
  local vstate::VariableState

  if isa(v0, VariableState)
    vstate = v0
  elseif isa(v0, Number) ||
    (isa(v0, Vector) && issubtype(eltype(v0), Number)) ||
    (isa(v0, Matrix) && issubtype(eltype(v0), Number))
    if isa(v, Parameter)
      vstate = default_state(v, v0, outopts)
    else
      vstate = default_state(v, v0)
    end
  else
    error("Variable state or state value of type $(typeof(v0)) not valid")
  end

  vstate
end

default_state(v::VariableVector, v0::Vector, outopts::Vector) =
  VariableState[default_state(v[i], v0[i], outopts[i]) for i in 1:length(v0)]

function default_state(v::VariableVector, v0::Vector, outopts::Vector, dpindex::IntegerVector)
  opts = fill(Dict(), length(v0))
  for i in 1:length(dpindex)
    opts[dpindex[i]] = outopts[i]
  end
  default_state(v, v0, opts)
end

show(io::IO, v::Variable) = print(io, "Variable [$(v.index)]: $(v.key) ($(typeof(v)))")

### Deterministic Variable subtypes

## Constant

mutable struct Constant <: Variable{Deterministic}
  key::Symbol
  index::Integer
end

Constant(key::Symbol) = Constant(key, 0)

default_state(variable::Constant, value::N) where {N<:Number} = BasicUnvVariableState(value)
default_state(variable::Constant, value::Vector{N}) where {N<:Number} = BasicMuvVariableState(value)
default_state(variable::Constant, value::Matrix{N}) where {N<:Number} = BasicMavVariableState(value)

show(io::IO, ::Type{Constant}) = print(io, "Constant")

dotshape(variable::Constant) = "trapezium"

## Hyperparameter

const Hyperparameter = Constant

## Data

mutable struct Data <: Variable{Deterministic}
  key::Symbol
  index::Integer
  update::Union{Function, Void}
end

Data(key::Symbol, index::Integer) = Data(key, index, nothing)
Data(key::Symbol, update::Union{Function, Void}) = Data(key, 0, update)
Data(key::Symbol) = Data(key, 0, nothing)

default_state(variable::Data, value::N) where {N<:Number} = BasicUnvVariableState(value)
default_state(variable::Data, value::Vector{N}) where {N<:Number} = BasicMuvVariableState(value)
default_state(variable::Data, value::Matrix{N}) where {N<:Number} = BasicMavVariableState(value)

show(io::IO, ::Type{Data}) = print(io, "Data")

dotshape(variable::Data) = "box"

## Transformation

mutable struct Transformation <: Variable{Deterministic}
  key::Symbol
  index::Integer
  transform::Function
  states::VariableStateVector
end

Transformation(key::Symbol, index::Integer, transform::Function=()->()) =
  Transformation(key, index, transform, VariableState[])

Transformation(key::Symbol, transform::Function=()->(), states::VariableStateVector=VariableState[]) =
  Transformation(key, 0, transform, states)

default_state(variable::Transformation, value::N) where {N<:Number} = BasicUnvVariableState(value)
default_state(variable::Transformation, value::Vector{N}) where {N<:Number} = BasicMuvVariableState(value)
default_state(variable::Transformation, value::Matrix{N}) where {N<:Number} = BasicMavVariableState(value)

show(io::IO, ::Type{Transformation}) = print(io, "Transformation")

dotshape(variable::Transformation) = "polygon"


text:
@testset "Distributions" begin

# Testing the second derivative of distributions provided by Distributions.jl
@test isapprox(_PMDSE.heslogpdf(_DST.Exponential(100.0),100.0),  0.0,     atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(_DST.Weibull(3.0,100.0),2.0),   -24.75,    atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(_DST.Normal(100.0,5.0),45.0),   -0.04,    atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(_DST.LogNormal(1.0, 1.0),3.0), 0.07275451 , atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(_DST.Beta(2.0,4.0),0.5),        -16.00,   atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(_DST.Gamma(2.0, 3.0),10.0),  -0.01,  atol=1e-8)

# An extended Beta distribution with support [0,1] should reduce to the 'normal'
# Beta distribution
alt = _DST.Beta(2.0,4.0)
dst = ExtendedBeta(2.0,4.0,0.0,1.0)
@test isapprox(_DST.mean(alt), _PMDSE.mean(dst),                          atol=1e-8)
@test isapprox(_DST.mode(alt), _PMDSE.mode(dst),                          atol=1e-8)
@test isapprox(_DST.skewness(alt), _PMDSE.skewness(dst),                  atol=1e-8)
@test isapprox(_DST.pdf(alt,0.5), _PMDSE.pdf(dst,0.5),                    atol=1e-8)
@test isapprox(_DST.logpdf(alt,0.5), _PMDSE.logpdf(dst,0.5),              atol=1e-8)
@test isapprox(_DST.gradlogpdf(alt,0.5), _PMDSE.gradlogpdf(dst,0.5),      atol=1e-8)

# Testing logpdf/gradlogpdf/heslogpdf of the extended beta distribution
dst = ExtendedBeta(2.0,4.0,10.0,100.0)
@test isapprox(_PMDSE.logpdf(dst,50.0), log(_PMDSE.pdf(dst,50.0)),          atol=1e-8)
@test isapprox(_PMDSE.gradlogpdf(dst,50.0),               -0.035,         atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(dst,50.0),                -0.001825,      atol=1e-8)

# Testing logpdf/gradlogpdf/heslogpdf of gaussian mixture models
beta = _PMDSE.ExtendedBeta(1.6339, 20.9022, 0.001, 2.0)
gmm = _DST.MixtureModel([_DST.Normal(0.14585265, sqrt(0.01144241))], [1.0])
g1 = _DST.Normal(gmm.components[1].μ, gmm.components[1].σ)
@test isapprox(_DST.logpdf(g1, beta.min),      _PMDSE.logpdf(gmm, beta.min),     atol = 1e-8)
@test isapprox(_DST.logpdf(g1, beta.max),      _PMDSE.logpdf(gmm, beta.max),     atol = 1e-8)
@test isapprox(_DST.gradlogpdf(g1, beta.min),  _PMDSE.gradlogpdf(gmm, beta.min), atol = 1e-8)
@test isapprox(_DST.gradlogpdf(g1, beta.max),  _PMDSE.gradlogpdf(gmm, beta.max), atol = 1e-8)
@test isapprox(_PMDSE.heslogpdf(g1, beta.min), _PMDSE.heslogpdf(gmm, beta.min),  atol = 1e-2)
@test isapprox(_PMDSE.heslogpdf(g1, beta.max), _PMDSE.heslogpdf(gmm, beta.max),  atol = 1e-2)

p = _Poly.Polynomial([0,1,1])
@test isapprox(_PMDSE.logpdf(p, 0.0), 0.0, atol=1e-8)
@test isapprox(_PMDSE.logpdf(p, 1.0), 2.0, atol=1e-8)
@test isapprox(_PMDSE.gradlogpdf(p, 1.0), 3.0, atol=1e-8)
@test isapprox(_PMDSE.gradlogpdf(p, 0.0), 1.0, atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(p, 0.0), 2.0, atol=1e-8)
@test isapprox(_PMDSE.heslogpdf(p, 1.0), 2.0, atol=1e-8)

end


text:
export pause, @assert_raises, @timeit, @p_str, @p_mstr, @dq_str, @sq_str, @s_mstr

macro p_str(s) s end
macro p_mstr(s) s end
macro dq_str(s) "\"" * s * "\"" end
macro sq_str(s) "'" * s * "'" end
macro s_mstr(s) strip(lstrip(s))  end

function pause(msg="PAUSED")
    println(msg)
    readline()
end

macro assert_raises(exc, code)
    quote
        try
            $(code)
        catch e
            isa(e, $(exc)) && return true
        end
        return false
    end
end

macro timeit(run, setup, n)
    quote
        eval($(setup))
        f = eval($(run))
        [f(i) for i in 1:100] # compile
        # benchmark(f, string($(run)), 500000)["AverageWall"][1]
        # n = 500000
        td = Array(Uint64, $n)
        q = 0.
        for (i, x) in zip(1:$n, rand(Uint64, $n))
            t_start = time_ns()
            q += f(x)
            t_end = time_ns()
            td[i] = t_end - t_start
        end
        println(q)  # so it can't optimize out no matter how smart
        mean(td)
    end
end


text:
"""
    lu_decompose(mat)
Decomposes a `n x n` non singular matrix into a lower triangular matrix (L) and an upper triangular matrix (U)
"""
function lu_decompose(mat)
	n = mat |> size |> first
	L = zeros(n, n)
	U = zeros(n, n)

	for i in 1:n
		for j in i:n
			s = 0
			for k in 1:i
				s += L[i,k] * U[k,j]
			end
			U[i,j] = mat[i,j] - s
		end

		for k in i:n
			if i == k
				L[i,i] = 1
    			else
				s = 0
				for j in 1:i
					s += L[k,j] * U[j,i]
				end
				L[k,i] = (mat[k,i] - s) / U[i,i]
			end
		end
	end

	return L, U
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
################################################
const default_umfpack_pivot_tolerance=SuiteSparse.UMFPACK.umf_ctrl[3+1]

"""
$(TYPEDEF)

Control parameters for Newton method.

Newton's method solves ``F(u)=0`` by the iterative procedure ``u_{i+1}=u_{i} - d_i F'(u_i)^{-1}F(u_i)``
starting with some inital value ``u_0``, where ``d_i`` is the damping.

$(TYPEDFIELDS)
"""
mutable struct NewtonControl

    """
    Tolerance (in terms of norm of Newton update):
    terminate if ``\\Delta_i=||u_{i+1}-u_i||_\\infty <`` `tol_absolute`.

    Default value: `1.0e-10`.
    """
    tol_absolute::Float64

    """
    Tolerance (relative to the first residual):
    terminate if ``\\Delta_i/\\Delta_0<`` `tol_relative`.

    Default value: `1.0e-10`.
    """
    tol_relative::Float64

    """
    Tolerance for roundoff error detection:
    terminate if   ``|\\;||u_{i+1}||_1 - ||u_{i}||_1\\;|/ ||u_{i}||_1<`` `tol_round` occured `max_round` times in a row.

    Default value: `1.0e-10`.
    """
    tol_round::Float64

    """
    Tolerance for monotonicity test:
    terminat with error if ``\\Delta_i/\\Delta_{i-1}>`` `1/tol_mono`.

    Default value: 1.0e-3
    """
    tol_mono::Float64

    """
    Initial damping parameter ``d_0``.

    Default value: `1.0`.

    To handle convergence problems, set this to a value less than 1.
    """
    damp_initial::Float64

    """
    Damping parameter growth factor: ``d_{i+1}=\\min(d_i\\cdot`` `max_growth` ``,1)``

    Default value: `1.2`

    Generally it should be set to a value between 1 and 2.
    """
    damp_growth::Float64

    """
    Maximum number of iterations.

    Default value: 100
    """
    max_iterations::Int32

    """
    Maximum number of reuses of lu factorization.
    It this value is 0, linear systems are solved by a sparse direct solver, and it's LU factorization is called in every Newton step.

    Otherwise, a BICGstab iterative method is used for linear system solution with an LU factorization as preconditioner which is updated only every `max_lureuse` Newton step.

    Default value: 0.
    """
    max_lureuse::Int32

    """
    Maximum number of consecutive iterations within roundoff error tolerance

    Default value: 1000 (effectively switching of this criterion).
    """
    max_round::Int32

    """
    Tolerance of iterative linear solver.

    Default value: `1.0e-4`.
    """
    tol_linear::Float64

    """
    Factorization kind for linear sytems (see ExtendableSparse.jl)
    """
    factorization::AbstractFactorization

    """
    Verbosity flag

    Default value: `false`.
    """
    verbose::Bool

    """
    Handle exceptions in [`embed!`](@ref) and  [`evolve!`](@ref) methods. If `true`, exceptions in Newton solves are catched, the embedding resp. time step is lowered, and solution is retried.

    Default value: `false`
    """
    handle_exceptions::Bool

    """
    Initial parameter step for [`embed!`](@ref) method.

    Default value: `1.0`
    """
    Δp::Float64

    """
    Maximal parameter step for [`embed!`](@ref) method

    Default value: `1.0`
    """
    Δp_max::Float64

    """
    Minimal parameter step for [`embed!`](@ref) method.

    Default value: `1.0e-3`
    """
    Δp_min::Float64

    """
    Time step for [`evolve!`](@ref) method.

    Default value: 0.1
    """
    Δt::Float64

    """
    Maximal time step for [`evolve!`](@ref) method.

    Default value: 1
    """
    Δt_max::Float64

    """
    Minimal time step for [`evolve!`](@ref) method.

    Default value: `1.0e-3`
    """
    Δt_min::Float64

    """
    Maximal step size growth for  [`evolve!`](@ref) method.

    Default: 1.2
    """
    Δt_grow::Float64

    """
    Optimal size of update for  [`evolve!`](@ref) method.
    The algorithm keeps this value approximately constant.

    Default: 0.1
    """
    Δu_opt::Float64

    """
    force first timestep
    """
    force_first_step::Bool

    """
    Edge cutoff for rectangular triangles.

    Default value: `0.0`.
    """
    edge_cutoff::Float64

    """
    Pivot tolerance for umfpack

    Default value: $(default_umfpack_pivot_tolerance).
    """
    umfpack_pivot_tolerance::Float64

    NewtonControl()=NewtonControl(new())
end

################################################
"""
$(TYPEDSIGNATURES)

Default constructor
"""
function NewtonControl(this)
    this.tol_absolute=1.0e-10
    this.tol_relative=1.0e-10
    this.tol_round=1.0e-10
    this.tol_mono=1.0e-3
    this.damp_initial=1.0
    this.damp_growth=1.2
    this.max_lureuse=0
    this.tol_linear=1.0e-4
    this.max_round=1000
    this.verbose=false
    this.max_iterations=100
    this.Δp=1
    this.Δp_max=1
    this.Δp_min=1.0e-3
    this.Δt=0.1
    this.Δt_max=1
    this.Δt_min=1.0e-3
    this.Δt_grow=1.2
    this.Δu_opt=0.1
    this.force_first_step=false
    this.handle_exceptions=false
    this.edge_cutoff=0.0
    this.umfpack_pivot_tolerance=default_umfpack_pivot_tolerance
    this.factorization=LUFactorization()
    return this
end

"""
````
timesteps!(control,Δt; grow=1.0)
````

Modify control data such that the time steps are fixed to a
geometric sequence such that Δt_new=Δt_old*grow
"""
function fixed_timesteps!(control,Δt; grow=1.0)
    control.Δt=Δt
    control.Δt_max=Δt
    control.Δt_min=Δt
    control.Δt_grow=grow
    control.Δu_opt=floatmax()
    control
end

function Base.show(io::IO, this::NewtonControl)
    for name in fieldnames(typeof(this))
        @printf("%16s = ",name)
        println(io,getfield(this,name))
    end
end
9


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""NodeDaemonEndpoints lists ports opened by daemons running on the Node.

    IoK8sApiCoreV1NodeDaemonEndpoints(;
        kubeletEndpoint=nothing,
    )

    - kubeletEndpoint::IoK8sApiCoreV1DaemonEndpoint : Endpoint on which Kubelet is listening.
"""
mutable struct IoK8sApiCoreV1NodeDaemonEndpoints <: SwaggerModel
    kubeletEndpoint::Any # spec type: Union{ Nothing, IoK8sApiCoreV1DaemonEndpoint } # spec name: kubeletEndpoint

    function IoK8sApiCoreV1NodeDaemonEndpoints(;kubeletEndpoint=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeDaemonEndpoints, Symbol("kubeletEndpoint"), kubeletEndpoint)
        setfield!(o, Symbol("kubeletEndpoint"), kubeletEndpoint)
        o
    end
end # type IoK8sApiCoreV1NodeDaemonEndpoints

const _property_map_IoK8sApiCoreV1NodeDaemonEndpoints = Dict{Symbol,Symbol}(Symbol("kubeletEndpoint")=>Symbol("kubeletEndpoint"))
const _property_types_IoK8sApiCoreV1NodeDaemonEndpoints = Dict{Symbol,String}(Symbol("kubeletEndpoint")=>"IoK8sApiCoreV1DaemonEndpoint")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeDaemonEndpoints }) = collect(keys(_property_map_IoK8sApiCoreV1NodeDaemonEndpoints))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeDaemonEndpoints }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeDaemonEndpoints[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeDaemonEndpoints }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeDaemonEndpoints[property_name]

function check_required(o::IoK8sApiCoreV1NodeDaemonEndpoints)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NodeDaemonEndpoints }, name::Symbol, val)
end


text:
using Seahawk
using Test

@testset "Seahawk.jl" begin
@testset "Cosine" begin
    results = Seahawk.search(stdout, "cosine"; silent=true)
    # check if Base.cos was found
    found = false
    for result in results
        if result[2].name == "cos" && result[2].mod == "Base"
            found = true
            break
        end
    end
    Seahawk.print_results(results[1:2])
    @test found
end

@testset "Pkg.add" begin
    results = Seahawk.search(stdout, "add package"; silent=true)
    # check if Base.cos was found
    found = false
    for result in results
        if result[2].name == "add" && result[2].mod == "Pkg"
            found = true
            break
        end
    end
    @test found
end

@testset "Seahawk.search" begin
    results = Seahawk.search(stdout, "Search docstrings"; silent=true)
    # check if Base.cos was found
    found = false
    for result in results
        if result[2].name == "search" && result[2].mod == "Seahawk"
            found = true
            break
        end
    end
    @test found
end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.concept.type.entitytype;
#
# public class EntityTypeSteps {}

# Scenario: Entity types can be created
@then("entity(person) is null: false") do context
    attr_person = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    @expect attr_person !== nothing
end

@then("entity(person) get supertype: entity") do context
    attr_person = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    attr_super = g.get_supertype(g.as_remote(attr_person, context[:transaction]))
    @expect attr_super.label.name == "entity"
end

# Scenario: Entity types can be deleted
@then("entity(company) is null: false") do context
    attr_comp = g.get(ConceptManager(context[:transaction]), EntityType, "company")
    @expect attr_comp !== nothing
end

@when("delete entity type: company") do context
    attr_comp = g.get(ConceptManager(context[:transaction]), EntityType, "company")
    g.delete(g.as_remote(attr_comp, context[:transaction]))
end

@then("entity(company) is null: true") do context
    attr_comp = g.get(ConceptManager(context[:transaction]), EntityType, "company")
    @expect attr_comp === nothing
end

@then("entity(entity) get subtypes do not contain:") do context
    types_inp = [db[1] for db in context.datatable]
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "entity")
    sub_types = g.get_subtypes(g.as_remote(attr, context[:transaction]))
    for i in 1:length(types_inp)
        inp_type = attr = g.get(ConceptManager(context[:transaction]), EntityType, types_inp[i])
        @expect !in(inp_type, sub_types)
    end
end

@when("delete entity type: person") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    g.delete(g.as_remote(attr, context[:transaction]))
end

@then("entity(person) is null: true") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    @expect attr === nothing
end

# Scenario: Entity types that have instances cannot be deleted
@when("\$x = entity(person) create new instance") do context
    ent_type = g.get(context[:cm], EntityType, "person")
    res = g.create(g.as_remote(ent_type, context[:transaction]))
    context[:x] = res
end

@then("delete entity type: person; throws exception") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    try
        g.delete(g.as_remote(attr_comp, context[:transaction]))
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can change labels
@then("entity(person) get label: person") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    @expect attr.label.name == "person"
end

@when("entity(person) set label: horse") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    g.set_label(g.as_remote(attr, context[:transaction]), "horse")
end

@then("entity(horse) is null: false") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "horse")
    @expect attr !== nothing
end

@then("entity(horse) get label: horse") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "horse")
    @expect attr.label.name == "horse"
end

@when("entity(horse) set label: animal") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "horse")
    g.set_label(g.as_remote(attr, context[:transaction]), "animal")
end

@then("entity(horse) is null: true") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "horse")
    @expect attr === nothing
end

@then("entity(animal) is null: false") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "animal")
    @expect attr !== nothing
end

@then("entity(animal) get label: animal") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "animal")
    @expect attr.label.name == "animal"
end

# Scenario: Entity types can be set to abstract
@then("entity(person) is abstract: true") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    res_abstract = g.is_abstract(g.as_remote(attr, context[:transaction]))
    @expect res_abstract
end

@then("entity(person) create new instance; throws exception") do context
    try
        ent_type = g.get(context[:cm], EntityType, "person")
        res = g.create(g.as_remote(ent_type, context[:transaction]))
    catch ex
        @expect ex !== nothing
    end
end

@then("entity(company) is abstract: false") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "company")
    res_abstract = g.is_abstract(g.as_remote(attr, context[:transaction]))
    @expect res_abstract === false
end

@when("entity(company) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), EntityType, "company")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

@then("entity(company) is abstract: true") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "company")
    res_abstract = g.is_abstract(g.as_remote(attr, context[:transaction]))
    @expect res_abstract
end

@then("entity(company) create new instance; throws exception") do context
    try
        ent_type = g.get(context[:cm], EntityType, "company")
        res = g.create(g.as_remote(ent_type, context[:transaction]))
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can be subtypes of other entity types
@when("put entity type: man") do context
    g.put(ConceptManager(context[:transaction]),
                            EntityType,
                            "man")
end

@when("put entity type: woman") do context
    g.put(ConceptManager(context[:transaction]),
                            EntityType,
                            "woman")
end

@when("put entity type: cat") do context
    g.put(ConceptManager(context[:transaction]),
                            EntityType,
                            "cat")
end

@when("put entity type: animal") do context
    g.put(ConceptManager(context[:transaction]),
                            EntityType,
                            "animal")
end

@when("entity(man) set supertype: person") do context
    _set_supertype_for_type(context, EntityType, "man", "person")
end

@when("entity(woman) set supertype: person") do context
    _set_supertype_for_type(context, EntityType, "woman", "person")
end

@when("entity(person) set supertype: animal") do context
    _set_supertype_for_type(context, EntityType, "person", "animal")
end

@when("entity(cat) set supertype: animal") do context
    _set_supertype_for_type(context, EntityType, "cat", "animal")
end

@then("entity(man) get supertype: person") do context
    res = _is_supertype_for_type(context, EntityType, "man", "person")
    @expect res
end

@then("entity(woman) get supertype: person") do context
    res = _is_supertype_for_type(context, EntityType, "woman", "person")
    @expect res
end

@then("entity(person) get supertype: animal") do context
    res = _is_supertype_for_type(context, EntityType, "person", "animal")
    @expect res
end

@then("entity(cat) get supertype: animal") do context
    res = _is_supertype_for_type(context, EntityType, "cat", "animal")
    @expect res
end

@then("entity(man) get supertypes contain:") do context
    res_contain = _supertypes_contain(context, EntityType, "man")
    @expect all(res_contain)
end

@then("entity(woman) get supertypes contain:") do context
    res_contain = _supertypes_contain(context, EntityType, "woman")
    @expect all(res_contain)
end

@then("entity(person) get supertypes contain:") do context
    res_contain = _supertypes_contain(context, EntityType, "person")
    @expect all(res_contain)
end

@then("entity(cat) get supertypes contain:") do context
    res_contain = _supertypes_contain(context, EntityType, "cat")
    @expect all(res_contain)
end

@then("entity(man) get subtypes contain:") do context
    res_contain = _subtypes_contain(context, EntityType, "man")
    @expect all(res_contain)
end

@then("entity(woman) get subtypes contain:") do context
    res_contain = _subtypes_contain(context, EntityType, "woman")
    @expect all(res_contain)
end

@then("entity(person) get subtypes contain:") do context
    res_contain = _subtypes_contain(context, EntityType, "person")
    @expect all(res_contain)
end

@then("entity(cat) get subtypes contain:") do context
    res_contain = _subtypes_contain(context, EntityType, "cat")
    @expect all(res_contain)
end

@then("entity(animal) get subtypes contain:") do context
    res_contain = _subtypes_contain(context, EntityType, "animal")
    @expect all(res_contain)
end

@then("entity(entity) get subtypes contain:") do context
    res_contain = _subtypes_contain(context, EntityType, "entity")
    @expect all(res_contain)
end

# Scenario: Entity types cannot subtype itself
@when("entity(person) set supertype: person; throws exception") do context
    try
        _set_supertype_for_type(context, EntityType, "person", "person")
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can have keys
@then("entity(person) get owns key types contain:") do context
    res = _get_owns_contain(context, EntityType, "person", true)
    @expect all(res)
end

# Scenario: Entity types can only commit keys if every instance owns a distinct key
@when("entity(person) set owns key type: email; throws exception") do context
    attr = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    owns_attr = g.get(ConceptManager(context[:transaction]), AttributeType, "email")
    try
        g.set_owns(g.as_remote(attr, context[:transaction]),
                owns_attr,
                true)
    catch ex
        @expect ex !== nothing
    end
end

@when("\$alice = attribute(email) as(string) put: alice@typedb.ai") do context
    attr = g.get(context[:cm], AttributeType, "email")
    context[:alice] = g.put(g.as_remote(attr, context[:transaction]), "alice@typedb.ai")
end

@when("\$bob = attribute(email) as(string) put: bob@typedb.ai") do context
    attr = g.get(context[:cm], AttributeType, "email")
    context[:bob] = g.put(g.as_remote(attr, context[:transaction]), "bob@typedb.ai")
end

@when("entity \$b set has: \$bob") do context
    set_has(context[:transaction], context[:b], context[:bob])
end

# Scenario: Entity types can unset keys
@when("entity(person) unset owns key type: email") do context
    attr_owns = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    attr_owned = g.get(ConceptManager(context[:transaction]), AttributeType, "email")
    g.unset_owns(g.as_remote(attr_owns, context[:transaction]),
                    attr_owned)
end

@when("entity(person) unset owns key type: username") do context
    attr_owns = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    attr_owned = g.get(ConceptManager(context[:transaction]), AttributeType, "username")
    g.unset_owns(g.as_remote(attr_owns, context[:transaction]),
                    attr_owned)
end

@then("entity(person) get owns key types do not contain:") do context
    res = _get_owns_contain(context, EntityType, "person", true)
    @expect !all(res)
end

# Scenario: Entity types cannot have keys of attributes that are not keyable
@given("entity(person) set owns key type: age") do context
    g.con.entity_set_owns("person", "age", context[:transaction], true)
end

@given("entity(person) set owns key type: name") do context
    g.con.entity_set_owns("person", "name", context[:transaction], true)
end

@given("entity(person) set owns key type: timestamp") do context
    g.con.entity_set_owns("person", "timestamp", context[:transaction], true)
end

@given("entity(person) set owns key type: is-open; throws exception") do context
    try
        g.con.entity_set_owns("person", "timestamp", context[:transaction])
    catch ex
        @expect ex !== nothing
    end
end

@given("entity(person) set owns key type: rating; throws exception") do context
    try
        g.con.entity_set_owns("person", "rating", context[:transaction], true)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can have attributes
@then("entity(person) get owns attribute types contain:") do context
    res = _get_owns_contain(context, EntityType, "person", false)
    @expect all(res)
end

# Scenario: Entity types can unset owning attributes
@when("entity(person) unset owns attribute type: age") do context
    attr_owns = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    attr_owned = g.get(ConceptManager(context[:transaction]), AttributeType, "age")
    g.unset_owns(g.as_remote(attr_owns, context[:transaction]),
                    attr_owned)
end

@then("entity(person) get owns attribute types do not contain:") do context
    res = _get_owns_contain(context, EntityType, "person", false)
    @expect all(res) === false
end

@when("entity(person) unset owns attribute type: name") do context
    attr_owns = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    attr_owned = g.get(ConceptManager(context[:transaction]), AttributeType, "name")
    g.unset_owns(g.as_remote(attr_owns, context[:transaction]),
                    attr_owned)
end

# Scenario: Entity types cannot unset owning attributes that are owned by existing instances
@when("\$alice = attribute(name) as(string) put: alice") do context
    attr = g.get(context[:cm], AttributeType, "name")
    context[:alice] = g.put(g.as_remote(attr, context[:transaction]), "alice")
end

@when("entity(person) unset owns attribute type: name; throws exception") do context
    attr_owns = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    attr_owned = g.get(ConceptManager(context[:transaction]), AttributeType, "name")

    try
        g.unset_owns(g.as_remote(attr_owns, context[:transaction]),
                        attr_owned)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can inherit keys and attributes
@when("put attribute type: reference, with value type: string") do context
    g.put(ConceptManager(context[:transaction]),
            AttributeType,
            "reference",
            VALUE_TYPE.STRING)
end

@when("put entity type: customer") do context
    g.put(ConceptManager(context[:transaction]),
            EntityType,
            "customer")
end

@when("entity(customer) set supertype: person") do context
    _set_supertype_for_type(context, EntityType, "customer", "person")
end

@given("entity(customer) set owns key type: reference") do context
    g.con.entity_set_owns("customer", "reference", context[:transaction], true)
end

@given("entity(customer) set owns attribute type: rating") do context
    g.con.entity_set_owns("customer", "rating", context[:transaction], false)
end

@then("entity(customer) get owns key types contain:") do context
    res = _get_owns_contain(context, EntityType, "customer", true)
    @expect all(res)
end

@when("put attribute type: points, with value type: double") do context
    g.put(ConceptManager(context[:transaction]),
                            AttributeType,
                            "points",
                            VALUE_TYPE.DOUBLE)
end

@when("put entity type: subscriber") do context
    g.put(ConceptManager(context[:transaction]),
                            EntityType,
                            "subscriber")
end

@when("entity(subscriber) set supertype: customer") do context
    _set_supertype_for_type(context, EntityType, "subscriber", "customer")
end

@when("entity(subscriber) set owns key type: license") do context
    g.con.entity_set_owns("subscriber", "license", context[:transaction], true)
end

@when("entity(subscriber) set owns attribute type: points") do context
    g.con.entity_set_owns("subscriber", "points", context[:transaction], false)
end

@then("entity(subscriber) get owns key types contain:") do context
    res = _get_owns_contain(context, EntityType, "subscriber", true)
    @expect all(res)
end

@then("entity(customer) get owns attribute types contain:") do context
    res = _get_owns_contain(context, EntityType, "customer", false)
    @expect all(res)
end

@then("entity(subscriber) get owns attribute types contain:") do context
    res = _get_owns_contain(context, EntityType, "subscriber", false)
    @expect all(res)
end

# Scenario: Entity types can inherit keys and attributes that are subtypes of each other
@when("attribute(username) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), AttributeType, "username")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

@when("attribute(score) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), AttributeType, "score")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

@when("attribute(reference) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), AttributeType, "reference")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

@when("attribute(reference) set supertype: username") do context
    _set_supertype_for_type(context, AttributeType, "reference", "username")
end

@when("attribute(rating) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), AttributeType, "rating")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

@when("attribute(rating) set supertype: score") do context
    _set_supertype_for_type(context, AttributeType, "rating", "score")
end

@when("entity(customer) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

@when("attribute(license) set supertype: reference") do context
    _set_supertype_for_type(context, AttributeType, "license", "reference")
end

@when("attribute(points) set supertype: rating") do context
    _set_supertype_for_type(context, AttributeType, "points", "rating")
end

@when("entity(subscriber) set abstract: true") do context
    attr_name = g.get(ConceptManager(context[:transaction]), EntityType, "subscriber")
    g.set_abstract(g.as_remote(attr_name, context[:transaction]))
end

# Scenario: Entity types can override inherited keys and attributes

@when("put attribute type: work-email, with value type: string") do context
    g.put(ConceptManager(context[:transaction]),
                            AttributeType,
                            "work-email",
                            VALUE_TYPE.STRING)
end

@when("attribute(work-email) set supertype: email") do context
    _set_supertype_for_type(context, AttributeType, "work-email", "email")
end

@when("put attribute type: nick-name, with value type: string") do context
    g.put(ConceptManager(context[:transaction]),
                            AttributeType,
                            "nick-name",
                            VALUE_TYPE.STRING)
end

@when("attribute(nick-name) set supertype: name") do context
    _set_supertype_for_type(context, AttributeType, "nick-name", "name")
end

@when("entity(customer) set owns key type: work-email as email") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "work-email")
    attr_overriden = g.get(ConceptManager(context[:transaction]), AttributeType, "email")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                true,
                attr_overriden)
end

@when("entity(customer) set owns attribute type: nick-name as name") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "nick-name")
    attr_overriden = g.get(ConceptManager(context[:transaction]), AttributeType, "name")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                false,
                attr_overriden)
end

@then("entity(customer) get owns key types do not contain:") do context
    res = _get_owns_contain(context, EntityType, "customer", true)
    @expect all(res) === false
end

@then("entity(customer) get owns attribute types do not contain:") do context
    res = _get_owns_contain(context, EntityType, "customer", false)
    @expect all(res) === false
end

@when("entity(subscriber) set owns key type: license as reference") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "subscriber")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "license")
    attr_overriden = g.get(ConceptManager(context[:transaction]), AttributeType, "reference")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                true,
                attr_overriden)
end

@when("entity(subscriber) set owns attribute type: points as rating") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "subscriber")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "points")
    attr_overriden = g.get(ConceptManager(context[:transaction]), AttributeType, "rating")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                false,
                attr_overriden)
end

@then("entity(subscriber) get owns key types do not contain:") do context
    res = _get_owns_contain(context, EntityType, "subscriber", true)
    @expect all(res) === false
end

@then("entity(subscriber) get owns attribute types do not contain:") do context
    res = _get_owns_contain(context, EntityType, "subscriber", false)
    @expect all(res) === false
end

# Scenario: Entity types can override inherited attributes as keys
@when("entity(customer) set owns key type: username as name") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "username")
    attr_overriden = g.get(ConceptManager(context[:transaction]), AttributeType, "name")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                true,
                attr_overriden)
end

# Scenario: Entity types can redeclare inherited attributes as keys (which will override)
@when("entity(customer) set owns key type: email") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "email")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                true)
end

@when("entity(subscriber) set supertype: person") do context
    _set_supertype_for_type(context, EntityType, "subscriber", "person")
end

@when("entity(subscriber) set owns key type: email") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "subscriber")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "email")

    g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                attr_set_owns,
                true)
end

# Scenario: Entity types cannot redeclare inherited attributes as attributes
@when("entity(customer) set owns attribute type: name; throws exception") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "name")

    try
        g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                    attr_set_owns,
                    false)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot redeclare inherited keys as keys or attributes
@when("entity(customer) set owns key type: email; throws exception") do context
    entity_owns = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    attr_set_owns = g.get(ConceptManager(context[:transaction]), AttributeType, "email")

    try
        g.set_owns(g.as_remote(entity_owns, context[:transaction]),
                    attr_set_owns,
                    true)
    catch ex
        @expect ex !== nothing
    end
end

@when("entity(customer) set owns attribute type: email; throws exception") do context
    try
        g.con.entity_set_owns("customer", "email", context[:transaction], true)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot redeclare inherited key attribute types
@when("put attribute type: customer-email, with value type: string") do context
    g.put(ConceptManager(context[:transaction]),
            AttributeType,
            "customer-email",
            VALUE_TYPE.STRING)
end

@when("attribute(customer-email) set supertype: email") do context
    _set_supertype_for_type(context, AttributeType, "customer-email", "email")
end

@when("entity(customer) set owns key type: customer-email") do context
    g.con.entity_set_owns("customer", "customer-email", context[:transaction], true)
end

@given("entity(subscriber) set owns key type: email; throws exception") do context
    try
        g.con.entity_set_owns("subscriber", "email", context[:transaction], true)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot redeclare overridden key attribute types
@given("entity(subscriber) set owns key type: customer-email; throws exception") do context
    try
        g.con.entity_set_owns("subscriber", "customer-email", context[:transaction], true)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot redeclare inherited owns attribute types
@when("put attribute type: customer-name, with value type: string") do context
    g.put(ConceptManager(context[:transaction]),
            AttributeType,
            "customer-name",
            VALUE_TYPE.STRING)
end

@when("attribute(customer-name) set supertype: name") do context
    _set_supertype_for_type(context, AttributeType, "customer-name", "name")
end

@when("entity(customer) set owns attribute type: customer-name") do context
    g.con.entity_set_owns("customer", "customer-name", context[:transaction], false)
end

@given("entity(subscriber) set owns attribute type: name; throws exception") do context
    try
        g.con.entity_set_owns("subscriber", "name", context[:transaction], true)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot redeclare overridden owns attribute types
@given("entity(subscriber) set owns attribute type: customer-name; throws exception") do context
    try
        g.con.entity_set_owns("subscriber", "customer-name", context[:transaction], false)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot redeclare overridden owns attribute types
@when("attribute(email) set supertype: username") do context
    _set_supertype_for_type(context, AttributeType, "email", "username")
end

@when("attribute(first-name) set supertype: name") do context
    _set_supertype_for_type(context, AttributeType, "first-name", "name")
end

@when("entity(person) set owns key type: email as username; throws exception") do context
    try
        g.con.entity_set_owns("person", "email", context[:transaction], false, "username")
    catch ex
        @expect ex !== nothing
    end
end

@when("entity(person) set owns attribute type: first-name as name; throws exception") do context
    try
        g.con.entity_set_owns("person", "first-name", context[:transaction], false, "name")
    catch ex
        @expect ex !== nothing
    end
end

# entity(customer) set owns attribute type: email as username; throws exception
@when("entity(customer) set owns attribute type: email as username; throws exception") do context
    try
        g.con.entity_set_owns("customer", "email", context[:transaction], false, "username")
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot override inherited keys and attributes other than with their subtypes
@when("entity(customer) set owns key type: reference as username; throws exception") do context
    try
        g.con.entity_set_owns("customer", "reference", context[:transaction], true, "username")
    catch ex
        @expect ex !== nothing
    end
end

@then("entity(customer) set owns attribute type: rating as name; throws exception") do context
    try
        g.con.entity_set_owns("customer", "rating", context[:transaction], false, "name")
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can play role types
@then("entity(person) get playing roles contain:") do context
    res = _get_playing_roles_contain("person", EntityType, context)
    @expect all(res)
end

@then("relation(marriage) get role(husband) get players contain:") do context
    res = _get_players_contain("marriage","husband", context)
    @expect all(res)
end

@then("relation(marriage) get role(wife) get players contain:") do context
    res = _get_players_contain("marriage","wife", context)
    @expect all(res)
end

# Scenario: Entity types can unset playing role types
@then("entity(person) unset plays role: marriage:husband") do context
    role = RoleType(g.Label("marriage","husband"), false)
    attr = g.get(ConceptManager(context[:transaction]),
                    EntityType,
                    "person")

    g.unset_plays(g.as_remote(attr, context[:transaction]),role)
end

@then("entity(person) get playing roles do not contain:") do context
    res = _get_playing_roles_contain("person", EntityType, context)
    @expect all(res) === false
end

@then("relation(marriage) get role(husband) get players do not contain:") do context
    res = _get_players_contain("marriage","husband", context, should_contain = false)
    @expect all(res) === false
end

@then("entity(person) unset plays role: marriage:wife") do context
    role = RoleType(g.Label("marriage","wife"), false)
    attr = g.get(ConceptManager(context[:transaction]),
                    EntityType,
                    "person")

    g.unset_plays(g.as_remote(attr, context[:transaction]),role)
end

@then("relation(marriage) get role(wife) get players do not contain:") do context
    res = _get_players_contain("marriage","wife", context)
    @expect all(res) === false
end

# Scenario: Attempting to unset playing a role type that an entity type cannot actually play throws
@then("entity(person) unset plays role: marriage:husband; throws exception") do context
    role = RoleType(g.Label("marriage","husband"), false)
    attr = g.get(ConceptManager(context[:transaction]),
                    EntityType,
                    "person")
    try
        g.unset_plays(g.as_remote(attr, context[:transaction]),role)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot unset playing role types that are currently played by existing instances
@when("\$m = relation(marriage) create new instance") do context
    mar_type = g.get(ConceptManager(context[:transaction]), RelationType, "marriage")
    context[:m] = g.create(g.as_remote(mar_type, context[:transaction]))
end

@then("entity(person) unset plays role: marriage:wife; throws exception") do context
    role = RoleType(g.Label("marriage","wife"), false)
    attr = g.get(ConceptManager(context[:transaction]),
                    EntityType,
                    "person")
    try
        g.unset_plays(g.as_remote(attr, context[:transaction]),role)
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types can inherit playing role types
@when("put relation type: parentship") do context
    g.put(ConceptManager(context[:transaction]),
        RelationType,
        "parentship")
end

@when("relation(parentship) set relates role: parent") do context
    g.set_relates(context[:transaction],
        g.Label("","parentship"),
                "parent")
end

@when("relation(parentship) set relates role: child") do context
    g.set_relates(context[:transaction],
        g.Label("","parentship"),
                "child")
end

@when("entity(animal) set plays role: parentship:parent") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "animal")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("parentship", "parent"),false))
end

@when("entity(animal) set plays role: parentship:child") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "animal")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("parentship", "child"),false))
end

@when("put relation type: sales") do context
    g.put(ConceptManager(context[:transaction]),
        RelationType,
        "sales")
end

@when("relation(sales) set relates role: buyer") do context
    g.set_relates(context[:transaction],
        g.Label("","sales"),
                "buyer")
end

@when("entity(customer) set plays role: sales:buyer") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "customer")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("sales", "buyer"),false))
end

@then("entity(customer) get playing roles contain:") do context
    res = _get_playing_roles_contain("customer", EntityType, context)
    @expect all(res)
end

@then("entity(animal) get playing roles contain:") do context
    res = _get_playing_roles_contain("animal", EntityType, context)
    @expect all(res)
end

# Scenario: Entity types can inherit playing role types that are subtypes of each other
@when("put relation type: fathership") do context
    g.put(ConceptManager(context[:transaction]),
        RelationType,
        "fathership")
end

@when("relation(fathership) set supertype: parentship") do context
    _set_supertype_for_type(context, RelationType, "fathership", "parentship")
end

@when("relation(fathership) set relates role: father as parent") do context
    g.set_relates(context[:transaction],
        g.Label("","fathership"),
                "father",
                "parent")
end

@when("entity(person) set plays role: parentship:parent") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("parentship", "parent"),false))
end

@when("entity(person) set plays role: parentship:child") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("parentship", "child"),false))
end

@when("entity(man) set plays role: fathership:father") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "man")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("fathership", "father"),false))
end

@then("entity(man) get playing roles contain:") do context
    res = _get_playing_roles_contain("man", EntityType, context)
    @expect all(res)
end

@when("put relation type: mothership") do context
    g.put(ConceptManager(context[:transaction]),
        RelationType,
        "mothership")
end

@when("relation(mothership) set supertype: parentship") do context
    _set_supertype_for_type(context, RelationType, "mothership", "parentship")
end

@when("relation(mothership) set relates role: mother as parent") do context
    g.set_relates(context[:transaction],
        g.Label("","mothership"),
                "mother",
                "parent")
end

@when("entity(woman) set plays role: mothership:mother") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "woman")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("mothership", "mother"),false))
end

@then("entity(woman) get playing roles contain:") do context
    res = _get_playing_roles_contain("woman", EntityType, context)
    @expect all(res)
end

# Scenario: Entity types can override inherited playing role types

@when("entity(man) set plays role: fathership:father as parentship:parent") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "man")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]),
        g.RoleType(g.Label("fathership", "father"),false),
        g.RoleType(g.Label("parentship", "parent"),false))
end

@then("entity(man) get playing roles do not contain:") do context
    res = _get_playing_roles_contain("man", EntityType, context)
    @expect all(res) === false
end

@when("entity(woman) set plays role: mothership:mother as parentship:parent") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "woman")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]),
        g.RoleType(g.Label("mothership", "mother"),false),
        g.RoleType(g.Label("parentship", "parent"),false))
end

@then("entity(woman) get playing roles do not contain:") do context
    res = _get_playing_roles_contain("woman", EntityType, context)
    @expect all(res) === false
end

# Scenario: Entity types cannot redeclare inherited/overridden playing role types
@when("entity(boy) set supertype: man") do context
    _set_supertype_for_type(context, EntityType, "boy", "man")
end

@when("entity(boy) set plays role: parentship:parent; throws exception") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "boy")
    try
        g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("parentship", "parent"),false))
    catch ex
        @expect ex !== nothing
    end
end

@then("entity(boy) set plays role: fathership:father; throws exception") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "boy")
    try
        g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("parentship", "parent"),false))
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot override declared playing role types
@then("entity(person) set plays role: fathership:father as parentship:parent; throws exception") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    try
        g.set_plays(g.as_remote(ent_pers,
            context[:transaction]),
            g.RoleType(g.Label("fathership", "father"),false),
            g.RoleType(g.Label("parentship", "parent"),false))
    catch ex
        @expect ex !== nothing
    end
end

# Scenario: Entity types cannot override inherited playing role types other than with their subtypes
@when("relation(fathership) set relates role: father") do context
    g.set_relates(context[:transaction],
        g.Label("","fathership"),
                "father")
end

@then("entity(man) set plays role: fathership:father as parentship:parent; throws exception") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "man")
    try
        g.set_plays(g.as_remote(ent_pers,
            context[:transaction]),
            g.RoleType(g.Label("fathership", "father"),false),
            g.RoleType(g.Label("parentship", "parent"),false))
    catch ex
        @expect ex !== nothing
    end
end

@when("relation(marriage) set plays role: locates:located") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), RelationType, "marriage")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("locates", "located"),false))
end

@then("relation(marriage) get playing roles contain:") do context
    res = _get_playing_roles_contain("marriage", RelationType, context)
    @expect all(res)
end

@when("put relation type: organises") do context
    g.put(context[:cm], RelationType, "organises")
end

@when("relation(organises) set relates role: organiser") do context
    g.set_relates(context[:transaction],
        g.Label("","organises"),
                "organiser")
end

@when("relation(organises) set relates role: organised") do context
    g.set_relates(context[:transaction],
        g.Label("","organises"),
                "organised")
end

@when("relation(marriage) set plays role: organises:organised") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), RelationType, "marriage")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("organises", "organised"),false))
end


text:
module Parameters

end # module


text:
module STFT
"""
Compute the short-time Fourier transform of a signal
Parameters:
    -- data: The signal to find the STFT of
    -- windowFunc: An array sampled from the window function to be applied before taking the FFTs. The length of this function defines the window length
    -- windowStep: The step ( in samples) to move the window between FFTs
    -- low: The lowest harmonic to include in the output
    -- high: The highest harmonic to include in the output
    -- transformFunction: A function to apply to the windowed, FFT'd pieces of _data_
    -- processingFunction: A function to do the processing. This function is passed
       the window index and the windowed, FFT'd piece and must return a value of type
       outType. For instance, to compute the amplitude along some path, processingFunction should be (dat, i)->dat[path[i]]
"""
function stft_processed{dataType,windowType,outType}(data::AbstractArray{dataType, 1}, windowFunc::AbstractArray{windowType, 1}, windowStep::Integer, low::Integer, high::Integer, transformFunction::Function, processingFunction::Function, outT::Type{outType})
    dataLength = size(data, 1)
    windowLength = size(windowFunc)[1]
    nHarmonics = 1 + high - low
    nSnippets = Int(floor((dataLength-windowLength)/windowStep))
    endSample = nSnippets*windowStep
    if low > windowLength
	  throw(ArgumentError("The low harmonic (parameter 'low') must be less than the length of the window"))
    end
    if low < 1
	    throw(ArgumentError("The low harmonic (parameter 'low') must be at least 1"))
    end
    STFT_processed = Array{outType, 1}(nSnippets)
    temp = Array{Float64, 1}(windowLength)
    tempFFT = Array{Float64, 1}(windowLength)
    tempFFTProcessed = Array{Float64, 1}(nHarmonics)
    windowPos = 1
    harmonic = 1
    plan = plan_rfft(temp, flags = FFTW.MEASURE)

    for sample in 1:windowStep:endSample
        temp = windowFunc.*view(data, sample:(sample+windowLength-1))
        tempFFT = plan*temp

        harmonic = 1
        for i in low:high
            tempFFTProcessed[harmonic] = transformFunction(tempFFT[i].re, tempFFT[i].im)
	    STFT_processed[windowPos] = processingFunction(tempFFTProcessed, windowPos)
            harmonic += 1
        end

        windowPos += 1
    end

    return STFT_processed
end

"""
Compute the short-time Fourier transform of a signal
Parameters:
    -- data: The signal to find the STFT of
    -- windowFunc: An array sampled from the window function to be applied before taking the FFTs. The length of this function defines the window length
    -- windowStep: The step ( in samples) to move the window between FFTs
    -- low: The lowest harmonic to include in the output
    -- high: The highest harmonic to include in the output
    -- transformFunction: A function to apply to the windowed, FFT'd pieces of _data_
"""
function stft{dataType,windowType}(data::AbstractArray{dataType, 1}, windowFunc::AbstractArray{windowType, 1}, windowStep::Integer, low::Integer, high::Integer, transformFunction::Function)
    dataLength = size(data, 1)
    windowLength = size(windowFunc)[1]
    nHarmonics = 1 + high - low
    nSnippets = Int(floor((dataLength-windowLength)/windowStep))
    endSample = nSnippets*windowStep
    if low > windowLength
	  throw(ArgumentError("The low harmonic (parameter 'low') must be less than the length of the window"))
    end
    if low < 1
	    throw(ArgumentError("The low harmonic (parameter 'low') must be at least 1"))
    end
    STFT = Array{Float64, 2}(nHarmonics, nSnippets)
    temp = Array{Float64, 1}(windowLength)
    tempFFT = Array{Float64, 1}(windowLength)
    windowPos = 1
    harmonic = 1
    plan = plan_rfft(temp, flags = FFTW.MEASURE)

    for sample in 1:windowStep:endSample
        temp = windowFunc.*view(data, sample:(sample+windowLength-1))
        tempFFT = plan*temp

        harmonic = 1
        for i in low:high
            STFT[harmonic, windowPos] = transformFunction(tempFFT[i].re, tempFFT[i].im)
            harmonic += 1
        end

        windowPos += 1
    end

    return STFT
end

"""Compute the short-time Fourier transform, using a linear detrend on each windowed snippet"""
function detrended_stft{dataType<:Any,windowType<:Any}(data::AbstractArray{dataType, 1}, startSample::Integer, endSample::Integer, windowFunc::AbstractArray{windowType, 1}, windowStep::Integer, low::Integer, high::Integer, transformFunction::Function)
  if endSample == -1
    endSample = size(data, 1)
  end

    dataLength = size(data, 1)
    usedDataLength = endSample - startSample
    windowLength = size(windowFunc)[1]
    nHarmonics = 1 + high - low
    nSnippets = Int(floor((usedDataLength-windowLength)/windowStep))
    endSample = nSnippets*windowStep
    if low > windowLength
	  throw(ArgumentError("The low harmonic (parameter 'low') must be less than the length of the window"))
    end
    if low < 1
	    throw(ArgumentError("The low harmonic (parameter 'low') must be at least 1"))
    end
    STFT = Array{Float64, 2}(nHarmonics, nSnippets)
    temp = Array{Float64, 1}(windowLength)
    tempFFT = Array{Float64, 1}(windowLength)
    windowPos = 1
    harmonic = 1
    plan = plan_rfft(temp, flags = FFTW.MEASURE)

    for sample in startSample:windowStep:endSample
	temp = collect(view(data, sample:(sample+windowLength-1)))
	(a, b) = linreg(collect(1.0:length(temp)), temp)
	for i in eachindex(temp)
	    temp[i] -= (a+b*i) #FIXME: Do we really want to subtract out the intercept here?
	    temp[i] *= windowFunc[i]
	end
        tempFFT = plan*temp

        harmonic = 1
        for i in low:high
            STFT[harmonic, windowPos] = transformFunction(tempFFT[i].re, tempFFT[i].im)
            harmonic += 1
        end

        windowPos += 1
    end

    return STFT
end

include("centroid.jl")
include("envelope.jl")
include("windows.jl")
include("transformFunctions.jl")

export stft
export detrended_stft
end


text:
module FactorioUtils

    export
        # Blueprint
        decode_blueprint,
        encode_blueprint

    include("blueprint.jl")
end # FactorioUtils Module


text:
using ThresholdStability
using Documenter

# DocMeta.setdocmeta!(ThresholdStability, :DocTestSetup, :(using ThresholdStability); recursive=true)

makedocs(
    authors="Sam Wycherley",
    sitename="ThresholdStability.jl",
    format=Documenter.HTML(prettyurls=get(ENV, "CI", "false") == "true"
    ),
    pages=[
        "Home" => "index.md",
        "Stability" => "stability.md",
        "CKSVAR" => "cksvar.md"
    ],
)

deploydocs(;
    repo="github.com/samwycherley/ThresholdStability.jl",
)