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

export QuEST32, QuEST64

module _QuEST_jl_Build_Setup

    include(joinpath("..","deps","build_setup.jl"))

end #^ module _QuEST_jl_Build_Setup

using Libdl: dlopen, dlclose, RTLD_LAZY, RTLD_DEEPBIND, RTLD_GLOBAL

const _QUEST_LIB = Ref{Ptr{Cvoid}}(0)

"""
Function `_dlopen_QuEST()` — loads QuEST dll, returns precision code.
"""
function _dlopen_QuEST(num::Int) ::Cint
    @assert num ∈ [32,64]
    let numstr = string(num)

        if _QuEST_jl_Build_Setup.EXPERT_BUILD
            global _QUEST_LIB[] =
                dlopen("libQuEST_"*numstr,
                       RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
        else
            global _QUEST_LIB[] =
                dlopen(joinpath(@__DIR__,"..","deps","QuEST","build"*numstr,"QuEST","libQuEST"),
                       RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
        end

    end #^ let

    return ccall(:getQuEST_PREC , Cint, (),)
end #^ fn _dlopen_QuEST

function _dlclose_QuEST() ::Nothing
    dlclose(_QUEST_LIB[])
    _QUEST_LIB[] = 0
    nothing;
end

"""
Functions `_quest_mtx_𝑛()` for 𝑛 ∈ {2,4} — convert Julia matrix into QuEST matrix
"""
function _quest_mtx_2(U ::Matrix{Complex{FLT}}) ::QuEST_Types.ComplexMatrix2 where FLT <: Real
    @assert size(U) == (2,2)
    u = QuEST_Types.ComplexMatrix2(
        ( (real(U[1,1]), real(U[1,2])), (real(U[2,1]), real(U[2,2])) ),
        ( (imag(U[1,1]), imag(U[1,2])), (imag(U[2,1]), imag(U[2,2])) )  )
    return u
end

"""
Functions `_quest_mtx_𝑛()` for 𝑛 ∈ {2,4} — convert Julia matrix into QuEST matrix
"""
function _quest_mtx_4(U ::Matrix{Complex{FLT}}) ::QuEST_Types.ComplexMatrix4 where FLT <: Real
    @assert size(U) = (4,4)
    u = QuEST_Types.ComplexMatrix4(
        ( ( real(U[1,1]), real(U[1,2]), real(U[1,3]), real(U[1,4]) ),
          ( real(U[2,1]), real(U[2,2]), real(U[2,3]), real(U[2,4]) ),
          ( real(U[3,1]), real(U[3,2]), real(U[3,3]), real(U[3,4]) ),
          ( real(U[4,1]), real(U[4,2]), real(U[4,3]), real(U[4,4]) ) ),
        ( ( imag(U[1,1]), imag(U[1,2]), imag(U[1,3]), imag(U[1,4]) ),
          ( imag(U[2,1]), imag(U[2,2]), imag(U[2,3]), imag(U[2,4]) ),
          ( imag(U[3,1]), imag(U[3,2]), imag(U[3,3]), imag(U[3,4]) ),
          ( imag(U[4,1]), imag(U[4,2]), imag(U[4,3]), imag(U[4,4]) ) )  )
end

################################################################################
# 32 bit QuEST
################################################################################

"""
Module `QuEST`𝑥𝑦 — Julia wrapper for QuEST with 𝑥𝑦-bit floating point precision.
"""
module QuEST32
export QuEST_Types
export QubitIdx

import .._dlopen_QuEST

function QuEST_init() ::Nothing
    prec = _dlopen_QuEST(32)
    @assert prec == 1 "Wrong precision. Please rebuild the package"
end

module QuEST_Types  # 32 bit qreal
    using CEnum

    const MPI_FLOAT  = Cfloat   # This is needed, because ...
    const MPI_DOUBLE = Cdouble  # ...

    include(joinpath("..","deps","questclang_common_32.jl"))

end #^ module QuEST_Types

const QubitIdx = Cint
const Qreal    = Float32
@assert sizeof(Qreal)==sizeof(QuEST_Types.qreal)

#
# Include C-wrappers, which will use the correct data types
#

using .QuEST_Types

include("data_structure_functions.jl")
include("QASM_logging.jl")
include("debugging.jl")
include("operators.jl")
include("decoherence.jl")
include("state_init.jl")
include("unitaries.jl")
include("calculations.jl")
include("gates.jl")

end #^ module QuEST32

################################################################################
# 64 bit QuEST
################################################################################

"""
Module `QuEST`𝑥𝑦 — Julia wrapper for QuEST with 𝑥𝑦-bit floating point precision.
"""
module QuEST64
export QuEST_Types
export QubitIdx

import .._dlopen_QuEST

function QuEST_init() ::Nothing
    prec = _dlopen_QuEST(64)
    @assert prec == 2 "Wrong precision. Please rebuild the package"
end

module QuEST_Types  # 64 bit qreal
    using CEnum

    const MPI_FLOAT  = Cfloat   # This is needed, because ...
    const MPI_DOUBLE = Cdouble  # ...

    include(joinpath("..","deps","questclang_common_64.jl"))

end #^ module QuEST_Types

const QubitIdx = Cint
const Qreal    = Float64
@assert sizeof(Qreal)==sizeof(QuEST_Types.qreal)

#
# Include C-wrappers, which will use the correct data types
#

using .QuEST_Types

include("data_structure_functions.jl")
include("QASM_logging.jl")
include("debugging.jl")
include("operators.jl")
include("decoherence.jl")
include("state_init.jl")
include("unitaries.jl")
include("calculations.jl")
include("gates.jl")

end #^ module QuEST64

end # module QuEST_jl
# EOF


text:
#=============================#
#== put needed modules here ==#
#=============================#
import JuliaChem

#================================#
#== JuliaChem execution script ==#
#================================#
function minimal_rhf(input_file)
  try
    #== read in input file ==#
    molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
      output=0)

    #== generate basis set ==#
    mol, basis = JuliaChem.JCBasis.run(molecule, model;
      output=2)

    #JuliaChem.JCMolecule.run(mol)

    #== perform scf calculation ==#
    rhf_energy = Dict()
    if haskey(keywords, "scf")
      rhf_energy = JuliaChem.JCRHF.Energy.run(mol, basis, keywords["scf"];
        output=0)
    else
      rhf_energy = JuliaChem.JCRHF.Energy.run(mol, basis;
        output=0)
    end

    #display(rhf_energy["Density"]); println()
    #display(rhf_energy["Energy-Weighted Density"]); println()

    #== perform gradient ==#
    #rhf_gradient = JuliaChem.JCGrad.run(mol, basis; output=2)

    return rhf_energy
  catch e
    bt = catch_backtrace()
    msg = sprint(showerror, e, bt)
    println(msg)

    JuliaChem.finalize()
    exit()
  end
end

ext: jl


text:
#=
using Pkg

tmp_packages = ["Kronecker", "Documenter"]

push!(LOAD_PATH,"../src/")

Pkg.activate(".")

Pkg.add.(tmp_packages) # IMPORTANT
=#

using Documenter, Kronecker, LinearAlgebra

makedocs(sitename="Kronecker.jl",
        authors = "Michiel Stock",
        format = :html,
        modules = [Kronecker])

deploydocs(
        repo = "github.com/MichielStock/Kronecker.jl.git",
        )


text:
using Assignment03
using Test

bases = ["ATGC"...; string.(["ATGC"...])]
cbases = ["TACG"...]
notbases = ["BDXZ"...; string.(["BDXZ"...])]
pur = ('A','G')
pyr = ('T','C')

@testset "Assignment03" begin

@testset "Question 1" begin
    @test all(x-> complement(x) isa Char, bases)

    for i in 1:4
        @test complement(bases[i]) == cbases[i]
        @test complement(bases[i+4]) == cbases[i]
    end
    for b in notbases
        @test_throws ErrorException complement(b)
    end
end

@testset "Question 2" begin
    @test all(ispurine, pur)
    @test all(ispyrimidine, pyr)
    @test all(p-> ispurine(p) == !ispyrimidine(p), bases)

    for b in notbases
        @test_throws ErrorException ispurine(b)
        @test_throws ErrorException ispyrimidine(b)
    end
end

@testset "Question 3" begin
    @test all(p-> base_type(p) == "purine", pur)
    @test all(p-> base_type(p) == "pyrimidine", pyr)
    for i in 1:4
        @test base_type(bases[i]) == base_type(bases[i+4])
        @test base_type(bases[i]) != base_type(cbases[i])
    end
    for b in notbases
        @test_throws ErrorException base_type(b)
    end
end

@testset "Question 4" begin
    @test gc_content(join(bases)) == 0.5
    @test gc_content(join(bases[1:5])) == 0.4
    @test gc_content(lowercase(join(bases))) == 0.5
    @test gc_content(lowercase(join(bases[1:5]))) == 0.4
end

@testset "Question 5" begin
    for b in bases, f in [complement, ispurine, ispyrimidine, base_type]
        @test f(b) == f(lowercase(b))
    end
end

end # Assignment03


text:
"""
## beemat!

This subroutine forms the strain-displacement matrix for
axisymmetric solids subjected to non-axisymmetric loading.

### Method
```julia
beemat!(bee, deriv)
```

### Arguments
```julia
* bee::Matrix{Float64}         : Bee matrix (Updated)
* deriv::Matrix{Float64}       : Derivative
```
"""
function beemat!(bee::Matrix{Float64},deriv::Matrix{Float64})
  #
  # This subroutine forms the bee matrix in 2-d (ih=3 or 4) or 3-d (ih=6).
  #
  ih = size(bee,1)
  nod = size(deriv,2)
  if ih == 3 || ih == 4
   for m in 1:nod
     k = 2m
     l = k-1
     x = deriv[1,m]
     y = deriv[2,m]
     bee[1,l] = x
     bee[3,k] = x
     bee[2,k] = y
     bee[3,l] = y
   end
  elseif ih == 6
   for m in 1:nod
     n = 3m
     k = n-1
     l = k-1
     x = deriv[1,m]
     y = deriv[2,m]
     z = deriv[3,m]
     bee[1,l] = x
     bee[4,k] = x
     bee[6,n] = x
     bee[2,k] = y
     bee[4,l] = y
     bee[5,n] = y
     bee[3,n] = z
     bee[5,k] = z
     bee[6,l] = z
   end
  else
   println("Wrong dimension for nst in bee matrix.")
  end
end


text:
# src/conversion.jl
# Convert Rationals into other Numeral Types

import Base: convert, Bool, float, AbstractFloat, big, trunc, floor, ceil, round, promote_rule

"""```
convert(R::Type{<:AbstractRatio}, x::AbstractRatio)
    === build_ratio(R, x)
```
Convert an abstract ratio into another one through the generic construction.
"""
function convert(R::Type{<:AbstractRatio}, x::AbstractRatio)
    return build_ratio(R, x)
end

"""```
convert(R::Type{<:AbstractRational}, x::AbstractRatio)
    === resolve(R, x)
```
Convert a ratio `x` to a rational through canonical resolution.
"""
function convert(R::Type{<:AbstractRational}, x::AbstractRatio)
    return resolve(R, x)
end

"""```
convert(::Type{BigFloat}, x::AbstractRatio)
    === BigFloat(numerator(x)) / denominator(x)
```
Convert a ratio `x` to a high precision `BigFloat`.
"""
function convert(::Type{BigFloat}, x::AbstractRatio)
    return BigFloat(numerator(x)) / denominator(x)
end

"""```
convert(::Type{T <: AbstractFloat}, x::AbstractRatio{S})
```
Convert a ratio `x` to an floating point type.
"""
function convert(::Type{T}, x::AbstractRatio{S}) where {T<:AbstractFloat,S}
    P = promote_type(T, S)
    return convert(T, convert(P, numerator(x)) / convert(P, denominator(x)))
end

"""```
convert(R::Type{<:AbstractRatio{T}}, n::Integer)
    === build_ratio(R, convert(T, n), one(T))
```
Convert an integer `n` into a ratio.
"""
function convert(R::Type{<:AbstractRatio{T}}, n::Integer) where {T}
    return build_ratio(R, convert(T, n), one(T))
end

"""```
convert(R::Type{<:AbstractRatio{T}}, x::Rational)
    === build_ratio(R, convert(T, numerator(x)), convert(T, denominator(x)))
```
Convert a standard rational `x` to a ratio.
"""
function convert(R::Type{<:AbstractRatio{T}}, x::Rational) where {T}
    return build_ratio(R, convert(T, numerator(x)), convert(T, denominator(x)))
end

"""```
convert(::Type{Rational{T}}, x::AbstractRatio)
    === convert(T, numerator(x)) // convert(T, denominator(x))
```
Convert a ratio `x` to a standard rational.
"""
function convert(::Type{Rational{T}}, x::AbstractRatio) where {T}
    return convert(T, numerator(x)) // convert(T, denominator(x))
end

"""```
Bool(x::AbstractRatio)
```
Convert a ratio `x` to a boolean value.
"""
function Bool(x::AbstractRatio)
    return iszero(x) ? false : isone(x) ? true : throw(InexactError(:Bool, Bool, x))
end

"""```
float(::Type{<:AbstractRatio{T}})
    === float(T)
```
Return the floating point type for the underlying integer type of a ratio.
"""
float(::Type{<:AbstractRatio{T}}) where {T} = float(T)

"""```
AbstractFloat(x::AbstractRatio)
    === float(numerator(x)) / float(denominator(x))
```
Convert a ratio `x` to an abstract floating point number.
"""
AbstractFloat(x::AbstractRatio) = float(numerator(x)) / float(denominator(x))

"""```
big(x::AbstractRatio)
    === build_ratio(typeof(x), big.(numden(x))...)
```
Convert a ratio `x` into a high-precision `big` ratio.
"""
big(x::AbstractRatio) = build_ratio(typeof(x), big.(numden(x))...)

"""```
trunc(::Type{T}, x::AbstractRatio)
    === convert(T, div(numden(x)...))
```
Truncate a ratio `x` to the integer type `T`.
"""
trunc(::Type{T}, x::AbstractRatio) where {T} = convert(T, div(numden(x)...))

"""```
trunc(x::AbstractRatio{T})
    === build_ratio(typeof(x), trunc(T, x))
```
Truncate a ratio `x` at the nearest integer.
"""
trunc(x::AbstractRatio{T}) where {T} = build_ratio(typeof(x), trunc(T, x))

"""```
floor(::Type{T}, x::AbstractRatio)
    === convert(T, fld(numden(x)...))
```
Round a ratio `x` to the highest integer of type `T` less than or equal to `x`.
"""
floor(::Type{T}, x::AbstractRatio) where {T} = convert(T, fld(numden(x)...))

"""```
floor(x::AbstractRatio{T})
    === build_ratio(typeof(x), floor(T, x))
```
Round a ratio `x` to the highest integer less than or equal to `x`.
"""
floor(x::AbstractRatio{T}) where {T} = build_ratio(typeof(x), floor(T, x))

"""```
ceil(::Type{T}, x::AbstractRatio)
    === convert(T, cld(numden(x)...))
```
Round a ratio `x` to the least integer of type `T` greater than or equal to `x`.
"""
ceil(::Type{T}, x::AbstractRatio) where {T} = convert(T, cld(numden(x)...))

"""```
ceil(x::AbstractRatio{T})
    === build_ratio(typeof(x), ceil(T, x))
```
Round a ratio `x` to the least integer greater than or equal to `x`.
"""
ceil(x::AbstractRatio{T}) where {T} = build_ratio(typeof(x), ceil(T, x))

"""```
round(::Type{T}, x::AbstractRational, r::RoundingMode = RoundNearest)
```
Round a rational `x` according to the rounding mode `r`.
"""
function round(::Type{T}, x::AbstractRational, r::RoundingMode = RoundNearest) where {T}
    return _round_rational(T, x, r)
end

"""```
round(::Type{T}, x::AbstractRatio{Bool})
```
Round a boolean ratio `x` to an integer of type `T`.
"""
function round(::Type{T}, x::AbstractRatio{Bool}, ::RoundingMode = RoundNearest) where {T}
    if denominator(x) == false && T <: Integer
        throw(DivideError())
    end
    return convert(T, x)
end

"""```
round(x::AbstractRatio, r::RoundingMode)
    === round(typeof(x), x, r)
```
Round the ratio `x` to the integer determined by the rounding mode `r`.
"""
round(x::AbstractRatio, r::RoundingMode) = round(typeof(x), x, r)

"""```
round(x::AbstractRatio{T})
    === build_ratio(typeof(x), round(T, x))
```
Round the ratio `x` to the nearest integer.
"""
round(x::AbstractRatio{T}) where {T} = build_ratio(typeof(x), round(T, x))

"""```
promote_rule(::Type{<:AbstractRatio{T}}, ::Type{S<:AbstractFloat})
    === promote_type(T, S)
```
Return promotion rule for ratios to abstract floats.
"""
function promote_rule(::Type{<:AbstractRatio{T}}, ::Type{S}) where {T,S<:AbstractFloat}
    return promote_type(T, S)
end

"""```
promote_rule(::Type{<:AbstractRatio{T}}, ::Type{Rational{S}})
    === Rational{promote_type(T, S)}
```
Return promotion rule for ratios to standard rationals.
"""
function promote_rule(::Type{<:AbstractRatio{T}}, ::Type{Rational{S}}) where {T,S}
    return Rational{promote_type(T, S)}
end


text:
# interval.jl

iscompact(d::TypedEndpointsInterval{:closed,:closed}) = true
iscompact(d::TypedEndpointsInterval) = false

##################
### An interval
##################

approx_indomain(x, d::AbstractInterval, tolerance) =
    (x <= rightendpoint(d)+tolerance) && (x >= leftendpoint(d)-tolerance)

function point_in_domain(d::AbstractInterval)
    isempty(d) && throw(BoundsError())
    mean(d)
end

isapprox(d1::AbstractInterval, d2::AbstractInterval; kwds...) =
    isapprox(leftendpoint(d1), leftendpoint(d2); kwds...) &&
    isapprox(rightendpoint(d1), rightendpoint(d2); kwds...)

boundary(d::TypedEndpointsInterval{:closed,:closed}) = Point(leftendpoint(d)) ∪ Point(rightendpoint(d))
boundary(d::TypedEndpointsInterval{:closed,:open}) = Point(leftendpoint(d))
boundary(d::TypedEndpointsInterval{:open,:closed}) = Point(rightendpoint(d))
boundary(d::TypedEndpointsInterval{:open,:open,T}) where T = EmptySpace{T}()
boundary(d::AbstractInterval) = boundary(Interval(d))

# We extend some functionality of intervals to mapped intervals
const MappedInterval{D <: AbstractInterval,T} = MappedDomain{D,T}

endpoints(d::MappedInterval) = forward_map(d) .* endpoints(source(d))
for op in (:leftendpoint, :rightendpoint)
    @eval $op(d::MappedInterval) = forward_map(d) * $op(source(d))
end

for op in (:isleftopen, :isrightopen, :isleftclosed, :isrightclosed)
    @eval $op(d::MappedInterval) = $op(source(d))
end

## Some special intervals follow, e.g.:
# - the unit interval [0,1]
# - the 'Chebyshev' interval [-1,1]
# - ...
# Unlike a generic interval, these specific intervals have no data, and only one
# type parameter (T).

"""
The abstract type `FixedInterval` is the supertype of intervals with endpoints
determined by the type, rather than field values. Examples include `UnitInterval`
and `ChebyshevInterval`.
"""
abstract type FixedInterval{L,R,T} <: TypedEndpointsInterval{L,R,T} end
const ClosedFixedInterval{T} = FixedInterval{:closed,:closed,T}

"""
Return an interval that is similar to the given interval, but with endpoints
`a` and `b` instead.
"""# Assume a closed interval by default
similar_interval(d::ClosedFixedInterval{T}, a, b) where {T} = ClosedInterval{float(T)}(a, b)

"The closed unit interval [0,1]."
struct UnitInterval{T} <: ClosedFixedInterval{T} end

UnitInterval() = UnitInterval{Float64}()

endpoints(d::UnitInterval{T}) where {T} = (zero(T), one(T))

"The closed interval [-1,1]."
struct ChebyshevInterval{T} <: ClosedFixedInterval{T}
end

ChebyshevInterval() = ChebyshevInterval{Float64}()

endpoints(d::ChebyshevInterval{T}) where {T} = (-one(T),one(T))

"The half-open positive halfline `[0,∞)`."
struct HalfLine{T} <: FixedInterval{:closed,:open,T} end
HalfLine() = HalfLine{Float64}()

endpoints(d::HalfLine{T}) where {T} = (zero(T), T(Inf))

indomain(x, d::HalfLine) = x >= 0

function similar_interval(d::HalfLine, a, b)
    @assert a == 0
    @assert isinf(b) && b > 0
    d
end

point_in_domain(d::HalfLine) = zero(eltype(d))

"The open negative halfline `(-∞,0)`."
struct NegativeHalfLine{T} <: FixedInterval{:open,:open,T} end
NegativeHalfLine() = NegativeHalfLine{Float64}()

endpoints(d::NegativeHalfLine{T}) where {T} = (-T(Inf), zero(T))

# Open at both endpoints

indomain(x, d::NegativeHalfLine) = x < 0

function similar_interval(d::NegativeHalfLine, a, b)
    @assert isinf(a) && a < 0
    @assert b == 0
    d
end

point_in_domain(d::NegativeHalfLine) = -one(eltype(d))

similar_interval(d::Interval{L,R,T}, a, b) where {L,R,T} =
    Interval{L,R,float(T)}(a, b)

#########################################
# A few set operations with known result
#########################################

# We define an exhaustive list of combinations of the four fixed intervals
# combined above in the routines 'intersect', 'union' and 'setdiff' where
# the output is known explicitly.

# Since fixed intervals are fully determined by their type,
# the result of intersect, union or setdiff is always known for two
# domains of the same type.
intersect(d1::D, d2::D) where {D <: FixedInterval} = d1
union(d1::D, d2::D) where {D <: FixedInterval} = d1
setdiff(d1::D, d2::D) where {D <: FixedInterval} = EmptySpace{eltype(D)}()

# [0,1] ∩ [-1,1] = [0,1]
intersect(d1::UnitInterval{T}, d2::ChebyshevInterval{T}) where {T} = UnitInterval{T}()
intersect(d1::ChebyshevInterval{T}, d2::UnitInterval{T}) where {T} = UnitInterval{T}()
# [0,1] ∩ [0,∞) = [0,1]
intersect(d1::UnitInterval{T}, d2::HalfLine{T}) where {T} = UnitInterval{T}()
intersect(d1::HalfLine{T}, d2::UnitInterval{T}) where {T} = UnitInterval{T}()
# [0,1] ∩ (-∞,0) = {}
intersect(d1::UnitInterval{T}, d2::NegativeHalfLine{T}) where {T} = EmptySpace{T}()
intersect(d1::NegativeHalfLine{T}, d2::UnitInterval{T}) where {T} = EmptySpace{T}()
# [-1,1] ∩ [0,∞) = [0,1]
intersect(d1::ChebyshevInterval{T}, d2::HalfLine{T}) where {T} = UnitInterval{T}()
intersect(d1::HalfLine{T}, d2::ChebyshevInterval{T}) where {T} = UnitInterval{T}()
# [0,∞) ∩ (-∞,0) = {}
intersect(d1::HalfLine{T}, d2::NegativeHalfLine{T}) where {T} = EmptySpace{T}()
intersect(d1::NegativeHalfLine{T}, d2::HalfLine{T}) where {T} = EmptySpace{T}()

# [0,1] ∪ [-1,1] = [-1,1]
union(d1::UnitInterval{T}, d2::ChebyshevInterval{T}) where {T} = ChebyshevInterval{T}()
union(d1::ChebyshevInterval{T}, d2::UnitInterval{T}) where {T} = ChebyshevInterval{T}()
# [0,1] ∪ [0,∞) = [0,∞)
union(d1::UnitInterval{T}, d2::HalfLine{T}) where {T} = HalfLine{T}()
union(d1::HalfLine{T}, d2::UnitInterval{T}) where {T} = HalfLine{T}()

# (-∞,0) ∪ [0,∞) = (-∞,∞)
# Note: T<:real to ensure that FullSpace{T} is not larger than intended.
union(d1::NegativeHalfLine{T}, d2::HalfLine{T}) where {T<:Real} = FullSpace{T}()
union(d1::HalfLine{T}, d2::NegativeHalfLine{T}) where {T<:Real} = FullSpace{T}()

# [0,1] ∖ [-1,1] = {}
setdiff(d1::UnitInterval{T}, d2::ChebyshevInterval{T}) where {T} = EmptySpace{T}()
# [0,1] ∖ [0,∞) = {}
setdiff(d1::UnitInterval{T}, d2::HalfLine{T}) where {T} = EmptySpace{T}()
# [0,1] ∖ (-∞,0) = [0,1]
setdiff(d1::UnitInterval{T}, d2::NegativeHalfLine{T}) where {T} = UnitInterval{T}()
# [-1,1] ∖ (-∞,0) = [0,1]
setdiff(d1::ChebyshevInterval{T}, d2::NegativeHalfLine{T}) where {T} = UnitInterval{T}()
# [0,∞) ∖ (-∞,0) = [0,∞)
setdiff(d1::HalfLine{T}, d2::NegativeHalfLine{T}) where {T} = HalfLine{T}()
# (-∞,0) ∖ [0,1] = (-∞,0)
setdiff(d1::NegativeHalfLine{T}, d2::UnitInterval{T}) where {T} = NegativeHalfLine{T}()
# (-∞,0) ∖ [0,∞) = (-∞,0)
setdiff(d1::NegativeHalfLine{T}, d2::HalfLine{T}) where {T} = NegativeHalfLine{T}()

#################################
# Conversions between intervals
#################################
# Avoid depcrecated warning: Warning: Constructors no longer fall back to `convert`.
# example: A constructor `AbstractInterval{Float64}(::IntervalSets.ClosedInterval{Float64})` should be defined instead.

function convert(::Type{UnitInterval{T}}, d::AbstractInterval) where {T}
    endpoints(d) == (0,1) || throw(InexactError(:convert,UnitInterval,d))
    UnitInterval{T}()
end

function convert(::Type{ChebyshevInterval{T}}, d::AbstractInterval) where {T}
    endpoints(d) == (-1,1)|| throw(InexactError(:convert,ChebyshevInterval,d))
    ChebyshevInterval{T}()
end

UnitInterval{T}(d::AbstractInterval) where T = convert(UnitInterval{T}, d)
ChebyshevInterval{T}(d::AbstractInterval) where T = convert(ChebyshevInterval{T}, d)

########################
# Arithmetic operations
########################

# Some computations with intervals simplify without having to use a mapped domain.
# This is only the case for Interval{L,R,T}, and not for any of the FixedIntervals
# because the endpoints of the latter are, well, fixed.

-(d::ChebyshevInterval) = d
-(d::AbstractInterval) = similar_interval(d, -rightendpoint(d), -leftendpoint(d))

for op in (:+, :-), Inter in (:AbstractInterval, :ClosedInterval)
    @eval $op(d::$Inter, x::Real) = similar_interval(d, $op(leftendpoint(d),x), $op(rightendpoint(d),x))
end

for Inter in (:AbstractInterval, :ClosedInterval)
    @eval begin
        +(x::Real, d::$Inter) = similar_interval(d, x+leftendpoint(d), x+rightendpoint(d))
        -(x::Real, d::$Inter) = similar_interval(d, x-rightendpoint(d), x-leftendpoint(d))
    end
end

for op in (:*, :/)
    @eval function $op(d::AbstractInterval, x::Real)
        if x ≥ 0 # -{x : 0 ≤ x ≤ 1} should be {x : -1 ≤ x ≤ 0}, not empty set {x : 0 ≤ x ≤ -1}
            similar_interval(d, $op(leftendpoint(d),x), $op(rightendpoint(d),x))
        else
            similar_interval(d, $op(rightendpoint(d),x), $op(leftendpoint(d),x))
        end
    end
end

for op in (:*, :\), Inter in (:AbstractInterval, :ClosedInterval)
    @eval function $op(x::Real, d::$Inter)
        if x ≥ 0 # -{x : 0 ≤ x ≤ 1} should be {x : -1 ≤ x ≤ 0}, not empty set {x : 0 ≤ x ≤ -1}
            similar_interval(d, $op(x,leftendpoint(d)), $op(x,rightendpoint(d)))
        else
            similar_interval(d, $op(x,rightendpoint(d)), $op(x,leftendpoint(d)))
        end
    end
end

function show(io::IO, d::ChebyshevInterval)
    print(io, Interval(d))
    print(io, " (Chebyshev)")
end
function show(io::IO, d::UnitInterval)
    print(io, Interval(d))
    print(io, " (Unit)")
end
function show(io::IO, d::HalfLine)
    print(io, Interval(d))
    print(io, " (HalfLine)")
end

function setdiff(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where T
    a1 = leftendpoint(d1)
    b1 = rightendpoint(d1)
    a2 = leftendpoint(d2)
    b2 = rightendpoint(d2)

    isempty(d1) && return d1
    isempty(d2) && return d1
    b1 < a2 && return d1
    a1 < a2 ≤ b1 ≤ b2 && return (a1 .. a2)
    a1 < a2 ≤ b2 < b1 && return UnionDomain(a1 .. a2) ∪ UnionDomain(b2 .. b1)
    a2 ≤ a1 < b2 < b1 && return (b2 .. b1)
    a2 ≤ a1 ≤ b1 ≤ b2 && return EmptySpace{T}()

    @assert b2 ≤ a1
    d1
end

setdiff(d1::AbstractInterval, d2::AbstractInterval) = setdiff(promote(d1,d2)...)

function *(map::AffineMap, domain::AbstractInterval)
    le = map*leftendpoint(domain)
    re = map*rightendpoint(domain)
    if le<re
        similar_interval(domain,le,re)
    else
        similar_interval(domain,re,le)
    end
end


text:
# Variable Number of Arguments

funvar(a, b, x...) = (a, b, x)

println(funvar(1, 2))
println(funvar(1, 2, 3))
println(funvar(1, 2, 3, 4, 5))

# Optional Arguments

incrementor(x, inc = 1) = x + inc
println(incrementor(1))
println(incrementor(1, 2))

# Keyword Arguments
funcalc(a, b; operator = +) = operator(a, b)

println(funcalc(1, 2))
println(funcalc(1, 2, operator = *))

# Do block syntax for function Arguments

v = map(x -> begin
  if (x < 0)
    return 0
  else
    return 1
  end
end, [-1, 0, 1])

v1 = map([-1, 0, 1]) do x
  if (x < 0)
    return 0
  else
    return 1
  end
end

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
#
# @description Defines the Cryptography.Ciphers module.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export BlockSizeError

# @description Defines cipher block size error exception type.
struct BlockSizeError <: Exception end

# Helpers
include("Ciphers/conversion.jl")

# Abstract types.
include("Ciphers/stream_cipher.jl")
include("Ciphers/block_cipher.jl")
include("Ciphers/multiblock_cipher.jl")

# Block ciphers.
include("Ciphers/null_block_cipher.jl")
include("Ciphers/feistel_cipher.jl")
include("Ciphers/des_cipher.jl")
include("Ciphers/triple_des_cipher.jl")
include("Ciphers/blowfish_cipher.jl")
include("Ciphers/rc4_cipher.jl")

# Multi-block ciphers.
include("Ciphers/cbc_mode_cipher.jl")

# Stream ciphers.


text:
using Benchmark
using VirtualArrays

function print_seperator(file::IOStream, name::AbstractString)
    seperator = "################################################################################"
    println(file)
    println(file, seperator)
    println(file, "# ", name)
    println(file, seperator)
    println(file)
end

open(joinpath(dirname(@__FILE__), "benchmark_result"), "w") do benchmark_file
    num_parents = 10
    len = 100
    num_test = 100
    last = len * num_parents
    half = div(last, 2)
    half_plus_1 = half + 1

    parents = []
    for i in 1:num_parents
        push!(parents, rand(Int64, len))
    end

    print_seperator(benchmark_file, "constructor")

    virtual_cat_bench() = virtual_cat(1, parents...)
    virtual_vcat_bench() = virtual_vcat(parents...)
    virtual_hcat_bench() = virtual_hcat(parents...)
    constructor_1() = VirtualArray{Int64, 1}(parents...)
    constructor_2() = VirtualArray{Int64, 1}(1, parents...)
    constructor_3() = VirtualArray{Int64, 2}(1, parents...)
    constructor_4() = VirtualArray{Int64, 2}(2, parents...)
    cat_bench() = cat(1, parents...)

    println(benchmark_file, benchmark(virtual_cat_bench, "virtual_cat", "virtual_cat(1, parents...)", num_test))
    println(benchmark_file, benchmark(virtual_vcat_bench, "virtual_vcat", "virtual_vcat(parents...)", num_test))
    println(benchmark_file, benchmark(virtual_hcat_bench, "virtual_hcat", "virtual_hcat(parents...)", num_test))
    println(benchmark_file, benchmark(constructor_1, "constructor_1", "VirtualArray{Int64, 1}(parents...)", num_test))
    println(benchmark_file, benchmark(constructor_2, "constructor_2", "VirtualArray{Int64, 1}(1, parents...)", num_test))
    println(benchmark_file, benchmark(constructor_3, "constructor_3", "VirtualArray{Int64, 1}(1, parents...)", num_test))
    println(benchmark_file, benchmark(constructor_4, "constructor_4", "VirtualArray{Int64, 2}(2, parents...)", num_test))
    println(benchmark_file, compare([virtual_cat_bench, virtual_vcat_bench, virtual_hcat_bench, constructor_1, constructor_2, constructor_3, constructor_4, cat_bench], num_test))
    v = VirtualArray{Int64, 1}(1, parents...)
    n = cat(1, parents...)

    print_seperator(benchmark_file, "getindex")

    get_index_virtual_bench1() = v[1]
    get_index_virtual_bench2() = v[last]
    get_index_virtual_bench3() = v[end]
    get_index_normal_bench1() = n[1]
    get_index_normal_bench2() = n[last]
    get_index_normal_bench3() = n[end]

    println(benchmark_file, benchmark(get_index_virtual_bench1, "getting first index", "v[1]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_bench2, "getting last index", "v[last]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_bench3, "getting last index with end", "v[end]", num_test))
    println(benchmark_file, benchmark(get_index_normal_bench1, "getting from normal array first index", "n[1]", num_test))
    println(benchmark_file, benchmark(get_index_normal_bench2, "getting from normal array last index", "n[last]", num_test))
    println(benchmark_file, benchmark(get_index_normal_bench3, "getting from normal array last index with end", "n[end]", num_test))
    println(benchmark_file, compare([get_index_virtual_bench1, get_index_virtual_bench2, get_index_virtual_bench3, get_index_normal_bench1, get_index_normal_bench2, get_index_normal_bench3], num_test))

    print_seperator(benchmark_file, "ranged getindex")

    get_index_virtual_small_start_bench() = v[1:1]
    get_index_virtual_small_end_bench() =   v[last:last]
    get_index_virtual_front_half_bench() =  v[1:half]
    get_index_virtual_end_half_bench() =    v[half_plus_1:last]
    get_index_virtual_entire_bench1() =     v[1:last]
    get_index_virtual_entire_bench2() =     v[:]
    get_index_normal_small_start_bench() =  n[1:1]
    get_index_normal_small_end_bench() =    n[last:last]
    get_index_normal_front_half_bench() =   n[1:half]
    get_index_normal_end_half_bench() =     n[half_plus_1:last]
    get_index_normal_entire_bench1() =      n[1:last]
    get_index_normal_entire_bench2() =      n[:]

    println(benchmark_file, benchmark(get_index_virtual_small_start_bench, "getting range from virtual first index", "v[1:1]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_small_end_bench, "getting range from virtual last index", "v[last:last]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_front_half_bench, "getting range from virtual first half", "v[1:half]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_end_half_bench, "getting range from virtual last half", "v[half_plus_1:last]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_entire_bench1, "getting range from virtual whole with number", "v[1:last]", num_test))
    println(benchmark_file, benchmark(get_index_virtual_entire_bench2, "getting range from virtual whole with colon", "v[:]", num_test))
    println(benchmark_file, benchmark(get_index_normal_small_start_bench, "getting range from normal first index", "n[1:1]", num_test))
    println(benchmark_file, benchmark(get_index_normal_small_end_bench, "getting range from normal last index", "n[last:last]", num_test))
    println(benchmark_file, benchmark(get_index_normal_front_half_bench, "getting range from normal first half", "n[1:half]", num_test))
    println(benchmark_file, benchmark(get_index_normal_end_half_bench, "getting range from normal last half", "n[half_plus_1:last]", num_test))
    println(benchmark_file, benchmark(get_index_normal_entire_bench1, "getting range from normal whole with number", "n[1:last]", num_test))
    println(benchmark_file, benchmark(get_index_normal_entire_bench2, "getting range from normal whole with colon", "n[:]", num_test))
    println(benchmark_file, compare(
        [
            get_index_virtual_small_start_bench, get_index_virtual_small_end_bench,
            get_index_virtual_front_half_bench, get_index_virtual_end_half_bench,
            get_index_virtual_entire_bench1, get_index_virtual_entire_bench2,
            get_index_normal_small_start_bench, get_index_normal_small_end_bench,
            get_index_normal_front_half_bench, get_index_normal_end_half_bench,
            get_index_normal_entire_bench1, get_index_normal_entire_bench2
        ], num_test))

    print_seperator(benchmark_file, "setindex")

    set_index_virtual_bench1() = v[1] = 1
    set_index_virtual_bench2() = v[last] = 1
    set_index_virtual_bench3() = v[end] = 1
    set_index_normal_bench1() = n[1] = 1
    set_index_normal_bench2() = n[last] = 1
    set_index_normal_bench3() = n[end] = 1

    println(benchmark_file, benchmark(set_index_virtual_bench1, "setting first index", "v[1] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_bench2, "setting last index", "v[last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_bench3, "setting last index with end", "v[end] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_bench1, "setting from normal array first index", "n[1] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_bench2, "setting from normal array last index", "n[last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_bench3, "setting from normal array last index with end", "n[end] = 1", num_test))
    println(benchmark_file, compare([set_index_virtual_bench1, set_index_virtual_bench2, set_index_virtual_bench3, set_index_normal_bench1, set_index_normal_bench2, set_index_normal_bench3], num_test))

    print_seperator(benchmark_file, "ranged setindex")

    set_index_virtual_small_start_bench() = v[1:1] = 1
    set_index_virtual_small_end_bench() =   v[last:last] = 1
    set_index_virtual_front_half_bench() =  v[1:half] = 1
    set_index_virtual_end_half_bench() =    v[half_plus_1:last] = 1
    set_index_virtual_entire_bench1() =     v[1:last] = 1
    set_index_virtual_entire_bench2() =     v[:] = 1
    set_index_normal_small_start_bench() =  n[1:1] = 1
    set_index_normal_small_end_bench() =    n[last:last] = 1
    set_index_normal_front_half_bench() =   n[1:half] = 1
    set_index_normal_end_half_bench() =     n[half_plus_1:last] = 1
    set_index_normal_entire_bench1() =      n[1:last] = 1
    set_index_normal_entire_bench2() =      n[:] = 1

    println(benchmark_file, benchmark(set_index_virtual_small_start_bench, "setting range from virtual first index", "v[1:1] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_small_end_bench, "setting range from virtual last index", "v[last:last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_front_half_bench, "setting range from virtual first half", "v[1:half] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_end_half_bench, "setting range from virtual last half", "v[half_plus_1:last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_entire_bench1, "setting range from virtual whole with number", "v[1:last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_virtual_entire_bench2, "setting range from virtual whole with colon", "v[:] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_small_start_bench, "setting range from normal first index", "n[1:1] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_small_end_bench, "setting range from normal last index", "n[last:last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_front_half_bench, "setting range from normal first half", "n[1:half] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_end_half_bench, "setting range from normal last half", "n[half_plus_1:last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_entire_bench1, "setting range from normal whole with number", "n[1:last] = 1", num_test))
    println(benchmark_file, benchmark(set_index_normal_entire_bench2, "setting range from normal whole with colon", "n[:] = 1", num_test))
    println(benchmark_file, compare(
        [
            set_index_virtual_small_start_bench, set_index_virtual_small_end_bench,
            set_index_virtual_front_half_bench, set_index_virtual_end_half_bench,
            set_index_virtual_entire_bench1, set_index_virtual_entire_bench2,
            set_index_normal_small_start_bench, set_index_normal_small_end_bench,
            set_index_normal_front_half_bench, set_index_normal_end_half_bench,
            set_index_normal_entire_bench1, set_index_normal_entire_bench2
        ], num_test))

    print_seperator(benchmark_file, "addition")

    addition_virtual_bench() = v + 1
    addition_normal_bench() = n + 1

    println(benchmark_file, benchmark(addition_virtual_bench, "addition virtual array", "v + 1", num_test))
    println(benchmark_file, benchmark(addition_normal_bench, "addition normal array", "n + 1", num_test))
    println(benchmark_file, compare([addition_virtual_bench, addition_normal_bench], num_test))

    print_seperator(benchmark_file, "subtraction")

    subtraction_virtual_bench() = v - 1
    subtraction_normal_bench() = n - 1

    println(benchmark_file, benchmark(subtraction_virtual_bench, "subtraction virtual array", "v - 1", num_test))
    println(benchmark_file, benchmark(subtraction_normal_bench, "subtraction normal array", "n - 1", num_test))
    println(benchmark_file, compare([subtraction_virtual_bench, subtraction_normal_bench], num_test))

    print_seperator(benchmark_file, "multiplication")

    multiplication_virtual_bench() = v * 2
    multiplication_normal_bench() = n * 2

    println(benchmark_file, benchmark(multiplication_virtual_bench, "multiplication virtual array", "v * 1", num_test))
    println(benchmark_file, benchmark(multiplication_normal_bench, "multiplication normal array", "n * 1", num_test))
    println(benchmark_file, compare([multiplication_virtual_bench, multiplication_normal_bench], num_test))

    print_seperator(benchmark_file, "division")

    division_virtual_bench() = v / 2
    division_normal_bench() = n / 2

    println(benchmark_file, benchmark(division_virtual_bench, "division virtual array", "v / 1", num_test))
    println(benchmark_file, benchmark(division_normal_bench, "division normal array", "n / 1", num_test))
    println(benchmark_file, compare([division_virtual_bench, division_normal_bench], num_test))

    print_seperator(benchmark_file, "size")

    size_virtual_bench() = size(v)
    size_normal_bench() = size(n)

    println(benchmark_file, benchmark(size_virtual_bench, "size of a virtual array", "size(v)", num_test))
    println(benchmark_file, benchmark(size_normal_bench, "size of a normal array", "size(n)", num_test))
    println(benchmark_file, compare([size_virtual_bench, size_normal_bench], num_test))

    print_seperator(benchmark_file, "length")

    length_virtual_bench() = length(v)
    length_normal_bench() = length(n)

    println(benchmark_file, benchmark(length_virtual_bench, "length of a virtual array", "length(v)", num_test))
    println(benchmark_file, benchmark(length_normal_bench, "length of a normal array", "length(n)", num_test))
    println(benchmark_file, compare([length_virtual_bench, length_normal_bench], num_test))

    print_seperator(benchmark_file, "eachindex")

    eachindex_virtual_bench() = eachindex(v)
    eachindex_normal_bench() = eachindex(n)

    println(benchmark_file, benchmark(eachindex_virtual_bench, "eachindex of a virtual array", "eachindex(v)", num_test))
    println(benchmark_file, benchmark(eachindex_normal_bench, "eachindex of a normal array", "eachindex(n)", num_test))
    println(benchmark_file, compare([eachindex_virtual_bench, eachindex_normal_bench], num_test))

    standard() = 0

    print_seperator(benchmark_file, "virtual array")

    println(benchmark_file, compare(
    [
        standard,
        virtual_cat_bench, virtual_vcat_bench, virtual_hcat_bench, constructor_1, constructor_2, constructor_3, constructor_4,
        get_index_virtual_bench1, get_index_virtual_bench2, get_index_virtual_bench3,
        get_index_virtual_small_start_bench, get_index_virtual_small_end_bench,
        get_index_virtual_front_half_bench, get_index_virtual_end_half_bench,
        get_index_virtual_entire_bench1, get_index_virtual_entire_bench2,
        set_index_virtual_bench1, set_index_virtual_bench2, set_index_virtual_bench3,
        set_index_virtual_small_start_bench, set_index_virtual_small_end_bench,
        set_index_virtual_front_half_bench, set_index_virtual_end_half_bench,
        set_index_virtual_entire_bench1, set_index_virtual_entire_bench2,
        addition_virtual_bench,
        subtraction_virtual_bench,
        multiplication_virtual_bench,
        division_virtual_bench,
        size_virtual_bench,
        length_virtual_bench,
        eachindex_virtual_bench
    ],
    num_test), "\n")

    print_seperator(benchmark_file, "normal array")

    println(benchmark_file, compare(
    [
        standard,
        cat_bench,
        get_index_normal_bench1, get_index_normal_bench2, get_index_normal_bench3,
        get_index_normal_small_start_bench, get_index_normal_small_end_bench,
        get_index_normal_front_half_bench, get_index_normal_end_half_bench,
        get_index_normal_entire_bench1, get_index_normal_entire_bench2,
        set_index_normal_bench1, set_index_normal_bench2, set_index_normal_bench3,
        set_index_normal_small_start_bench, set_index_normal_small_end_bench,
        set_index_normal_front_half_bench, set_index_normal_end_half_bench,
        set_index_normal_entire_bench1, set_index_normal_entire_bench2,
        addition_normal_bench,
        subtraction_normal_bench,
        multiplication_normal_bench,
        division_normal_bench,
        size_normal_bench,
        length_normal_bench,
        eachindex_normal_bench
    ],
    num_test), "\n")

    print_seperator(benchmark_file, "overall")

    println(benchmark_file, compare(
    [
        standard,
        virtual_cat_bench, virtual_vcat_bench, virtual_hcat_bench, constructor_1, constructor_2, constructor_3, constructor_4, cat_bench,
        get_index_virtual_bench1, get_index_virtual_bench2, get_index_virtual_bench3, get_index_normal_bench1, get_index_normal_bench2, get_index_normal_bench3,
        get_index_virtual_small_start_bench, get_index_virtual_small_end_bench,
        get_index_virtual_front_half_bench, get_index_virtual_end_half_bench,
        get_index_virtual_entire_bench1, get_index_virtual_entire_bench2,
        get_index_normal_small_start_bench, get_index_normal_small_end_bench,
        get_index_normal_front_half_bench, get_index_normal_end_half_bench,
        get_index_normal_entire_bench1, get_index_normal_entire_bench2,
        set_index_virtual_bench1, set_index_virtual_bench2, set_index_virtual_bench3, set_index_normal_bench1, set_index_normal_bench2, set_index_normal_bench3,
        set_index_virtual_small_start_bench, set_index_virtual_small_end_bench,
        set_index_virtual_front_half_bench, set_index_virtual_end_half_bench,
        set_index_virtual_entire_bench1, set_index_virtual_entire_bench2,
        set_index_normal_small_start_bench, set_index_normal_small_end_bench,
        set_index_normal_front_half_bench, set_index_normal_end_half_bench,
        set_index_normal_entire_bench1, set_index_normal_entire_bench2,
        addition_virtual_bench, addition_normal_bench,
        subtraction_virtual_bench, subtraction_normal_bench,
        multiplication_virtual_bench, multiplication_normal_bench,
        division_virtual_bench, division_normal_bench,
        size_virtual_bench, size_normal_bench,
        length_virtual_bench, length_normal_bench,
        eachindex_virtual_bench, eachindex_normal_bench
    ],
    num_test), "\n")
end


text:
println(true && false)
println(true && true)
println(false && true)
println(false && false)

println(true || false)
println(true || true)
println(false || true)
println(false || false)

x = true
println(!x)

x = !x
println(!x)

# priorités
println(!false && false)
println(!true || true)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
"""
    envar_penalty(estim_settings::EstimSettings, Σ::AbstractArray{Float64}, Ψ::SubArray{Float64}, Φ::FloatArray)

Compute the value of the loglikelihood penalty.
"""
envar_penalty(estim_settings::EstimSettings, Σ::AbstractArray{Float64}, Ψ::SubArray{Float64}, Φ::FloatArray) = tr(inv(Σ)*(envar_penalty_ridge(estim_settings, Ψ) + envar_penalty_lasso(estim_settings, Ψ, Φ)))::Float64;

# Ridge and LASSO components
envar_penalty_ridge(estim_settings::EstimSettings, Ψ::SubArray{Float64}) = ((1-estim_settings.α)/2) .* Symmetric(Ψ*estim_settings.Γ*Ψ')::SymMatrix;
envar_penalty_lasso(estim_settings::EstimSettings, Ψ::SubArray{Float64}, Φ::FloatArray) = (estim_settings.α/2) .* Symmetric((Ψ .* sqrt.(Φ))*estim_settings.Γ*(Ψ .* sqrt.(Φ))')::SymMatrix;

"""
    ksmoother_ecm!(estim_settings::EstimSettings, kalman_settings::MutableKalmanSettings, status::KalmanStatus)

Kalman smoother: RTS smoother from the last evaluated time period in status to t==0.

This instance of the smoother returns the ECM statistics and updates the initial conditions in KalmanSettings.

# Model
The state space model used below is,

``Y_{t} = B*X_{t} + e_{t}``

``X_{t} = C*X_{t-1} + v_{t}``

Where ``e_{t} ~ N(0, R)`` and ``v_{t} ~ N(0, V)``.

# Arguments
- `estim_settings`: EstimSettings struct
- `kalman_settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function ksmoother_ecm!(estim_settings::EstimSettings, kalman_settings::MutableKalmanSettings, status::KalmanStatus)

    # Memory pre-allocation
    E = zeros(estim_settings.n, estim_settings.n);
    F = zeros(estim_settings.n, estim_settings.np);
    G = zeros(estim_settings.np, estim_settings.np);
    Xs = status.X_post;
    Ps = status.P_post;

    # Loop over t
    for t=status.t:-1:2

        # Pointers
        Xp = status.history_X_prior[t];
        Pp = status.history_P_prior[t];
        Xf_lagged = status.history_X_post[t-1];
        Pf_lagged = status.history_P_post[t-1];

        # Smoothed estimates for t-1
        J1 = compute_J1(Pf_lagged, Pp, kalman_settings);
        Xs_old = backwards_pass(Xf_lagged, J1, Xs, Xp);
        Ps_old = backwards_pass(Pf_lagged, J1, Ps, Pp);

        # Update ECM statistics
        update_ecm_stats!(estim_settings, Xs, Xs_old, Ps, Ps_old, E, F, G);

        # Update Xs and Ps
        Xs = copy(Xs_old);
        Ps = copy(Ps_old);
    end

    # Pointers
    Xp = status.history_X_prior[1];
    Pp = status.history_P_prior[1];

    # Compute smoothed estimates for t==0
    J1 = compute_J1(kalman_settings.P0, Pp, kalman_settings);
    kalman_settings.X0 = backwards_pass(kalman_settings.X0, J1, Xs, Xp);
    kalman_settings.P0 = backwards_pass(kalman_settings.P0, J1, Ps, Pp);

    # Update ECM statistics
    update_ecm_stats!(estim_settings, Xs, kalman_settings.X0, Ps, kalman_settings.P0, E, F, G);

    # Use Symmetric for E and G
    E_sym = Symmetric(E)::SymMatrix;
    G_sym = Symmetric(G)::SymMatrix;

    # Return ECM statistics
    return E_sym, F, G_sym;
end

"""
    update_ecm_stats!(estim_settings::EstimSettings, Xs::FloatVector, Xs_old::FloatVector, Ps::SymMatrix, Ps_old::SymMatrix, E::FloatArray, F::FloatArray, G::FloatArray)

Update the ECM statistics.
"""
function update_ecm_stats!(estim_settings::EstimSettings, Xs::FloatVector, Xs_old::FloatVector, Ps::SymMatrix, Ps_old::SymMatrix, E::FloatArray, F::FloatArray, G::FloatArray)

    # Views
    Xs_view = @view Xs[1:estim_settings.n];
    Ps_view = @view Ps[1:estim_settings.n,1:estim_settings.n];
    Xs_old_view = @view Xs_old[1:estim_settings.np];
    Ps_old_view = @view Ps_old[1:estim_settings.np,1:estim_settings.np];
    PPs_view = @view Ps[1:estim_settings.n,estim_settings.n+1:end];

    # Update ECM statistics
    E .+= Xs_view*Xs_view' + Ps_view;
    F .+= Xs_view*Xs_old_view' + PPs_view;
    G .+= Xs_old_view*Xs_old_view' + Ps_old_view;
end

"""
    ecm(estim_settings::EstimSettings)

Estimate an elastic-net VAR(p) using the ECM algorithm in Pellegrino (2020).

# Arguments
- `estim_settings`: settings used for the estimation

# References
Pellegrino (2020)
"""
function ecm(estim_settings::EstimSettings)

    # Check inputs
    check_bounds(estim_settings.p, 1);
    check_bounds(estim_settings.λ, 0);
    check_bounds(estim_settings.α, 0, 1);
    check_bounds(estim_settings.β, 1);
    check_bounds(estim_settings.max_iter, 3);
    check_bounds(estim_settings.max_iter, estim_settings.prerun);
    check_bounds(estim_settings.n, 2); # It supports only multivariate models (for now ...)

    # Initialise using the coordinate descent algorithm
    verb_message(estim_settings.verb, "ecm > initialisation");
    Ψ_init, Σ_init = coordinate_descent(estim_settings);

    #=
    The state vector includes additional n terms with respect to the standard VAR companion form representation.
    This is to estimate the lag-one covariance smoother as in Watson and Engle (1983).
    =#

    # State-space parameters
    kalman_settings = MutableKalmanSettings(estim_settings.Y,
                                            [Matrix{Float64}(I, estim_settings.n, estim_settings.n) zeros(estim_settings.n, estim_settings.np)],
                                            Symmetric(Matrix{Float64}(I, estim_settings.n, estim_settings.n).*estim_settings.ε)::SymMatrix,
                                            ext_companion_form(Ψ_init, Σ_init)...);

    # Initialise additional variables
    Ψ = @view kalman_settings.C[1:estim_settings.n, 1:estim_settings.np];
    Φ = @. 1 / (abs(Ψ) + estim_settings.ε);

    # ECM controls
    pen_loglik_old = -Inf;
    pen_loglik_new = -Inf;

    # Run ECM
    for iter=1:estim_settings.max_iter

        # Run Kalman filter
        status = KalmanStatus();
        for t=1:kalman_settings.T
            kfilter!(kalman_settings, status);
        end

        if iter > estim_settings.prerun

            # New penalised loglikelihood
            Σ = Symmetric(@view parent(kalman_settings.V)[1:estim_settings.n, 1:estim_settings.n]);
            pen_loglik_new = status.loglik - envar_penalty(estim_settings, Σ, Ψ, Φ);
            verb_message(estim_settings.verb, "ecm > iter=$(iter-estim_settings.prerun), penalised loglik=$(round(pen_loglik_new, digits=5))");

            # Stop when the ECM algorithm converges
            if iter > estim_settings.prerun+1
                if isconverged(pen_loglik_new, pen_loglik_old, estim_settings.tol, estim_settings.ε, true)
                    verb_message(estim_settings.verb, "ecm > converged!\n");
                    break;
                end
            end

            # Store current run information
            pen_loglik_old = pen_loglik_new;

        else
            verb_message(estim_settings.verb, "ecm > prerun $iter (out of $(estim_settings.prerun))");
        end

        E, F, G = ksmoother_ecm!(estim_settings, kalman_settings, status);

        # VAR(p) coefficients
        Φ = @. 1 / (abs(Ψ) + estim_settings.ε);
        for i=1:estim_settings.n
            Φ_i = @view Φ[i, :];
            F_i = @view F[i,:];
            XX_i = Symmetric(G + estim_settings.Γ.*((1-estim_settings.α)*I + estim_settings.α.*Diagonal(Φ_i)))::SymMatrix;
            kalman_settings.C[i, 1:estim_settings.np] = inv(XX_i)*F_i;
        end

        # Covariance matrix of the VAR(p) residuals
        parent(kalman_settings.V)[1:estim_settings.n, 1:estim_settings.n] =
            Symmetric(E-F*Ψ'-Ψ*F'+Ψ*G*Ψ' + (1-estim_settings.α).*(Ψ*estim_settings.Γ*Ψ') + estim_settings.α.*((Ψ.*sqrt.(Φ))*estim_settings.Γ*(Ψ.*sqrt.(Φ))'))::SymMatrix ./ estim_settings.T;
    end

    # Return output
    out_kalman_settings = ImmutableKalmanSettings(estim_settings.Y_output,
                                                  kalman_settings.B, kalman_settings.R,
                                                  kalman_settings.C, kalman_settings.V,
                                                  kalman_settings.X0, kalman_settings.P0);
    return out_kalman_settings;
end


text:
import Base.close, Base.start, Base.next, Base.done

export tabix, query

type Tabix
    ptr :: Ptr{Uint8}
    is_open :: Bool
end

function Base.close(tbi :: Tabix)
    if tbi.is_open
        ccall((:ti_close, "libtabix"), Void, (Ptr{Uint8},), tbi.ptr)
        tbi.is_open = false
    end
    tbi
end

function tabix(path :: String)
    ptr = ccall((:ti_open, "libtabix"),
                Ptr{Uint8},
                (Ptr{Uint8}, Int),
                bytestring(path),
                0)
    tbi = Tabix(ptr, true)
    finalizer(tbi, Base.close)
    tbi
end

# Iterator stuff

type TabixIterator
    ptr :: Ptr{Uint8}
    owner :: Tabix
    contig :: String
    range :: Range1{Int}

    next_line :: String
    done :: Bool
end

function finalize(tbii :: TabixIterator)
    ccall((:ti_iter_destroy, "libtabix"),
          Void,
          (Ptr{Uint8},),
          tbii.ptr)
    tbii.done = true
end

function TabixIterator(tbi :: Tabix, contig :: String, range :: Range1{Int})
    ptr = ccall((:ti_query, "libtabix"),
                Ptr{Uint8},
                (Ptr{Uint8}, Ptr{Uint8}, Int, Int),
                tbi.ptr, bytestring(contig), range[1], range[end])
    if (ptr == C_NULL) # Zero search results
        tbii = TabixIterator(ptr, tbi, contig, range, "", true)
    else
        tbii = TabixIterator(ptr, tbi, contig, range, "", false)
        finalizer(tbii, finalize)
        advance!(tbii)
    end
    tbii
end

function advance!(tbii :: TabixIterator)
    line = ccall((:ti_read, "libtabix"),
                 Ptr{Uint8},
                 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}),
                 tbii.owner.ptr, tbii.ptr, 0)
    if (line == C_NULL)
        tbii.done = true
    else
        tbii.next_line = bytestring(line)
    end
end

# FIXME: make a new TabixIterator
Base.start(tbii :: TabixIterator) = nothing
function Base.next(tbii :: TabixIterator, nil)
    data = tbii.next_line
    advance!(tbii)
    (data, nothing)
end
Base.done(tbii :: TabixIterator, nil) = tbii.done

function query(tbi :: Tabix, contig :: String, range :: Range1{Int})
    TabixIterator(tbi, contig, range)
end


text:
import MechanicalSketch
import MechanicalSketch: mm, kg, m, set_scale_sketch, HE, WI, EM, PT, circle, text, N, kN, PALETTE
import MechanicalSketch: color_with_lumin, empty_figure, background, sethue, O, °, finish
import MechanicalSketch: line, setline, setopacity, scale_pt_to_unit, empty_figure
let
# Rope data from https://www.hendrikvedergroup.com/_asset/_public/Hendrik-Veder-Group/Downloads/8896-03-Dyneema-folder-offset_LR-v4.pdf
# Diameters below 12 mm are excluded here, because they vary too much for reasonable generalization. Properties are generally better for
# smaller rope
ROPE_D = [12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44,
    46, 48, 50, 52, 56, 60, 64, 68, 72, 80, 88, 96]mm
ROPE_WEIGHT = [9.50, 12.80, 16, 20.80, 25.50, 30.50, 35.80, 41, 46.50, 52, 57, 62.50, 68, 74, 84, 93,
    102, 111, 121, 131, 141, 163, 175, 200, 226, 254, 313, 379, 451]kg/100m
ROPE_MINIMUM_BREAK = [161.90, 215.80, 269.80, 343.30, 407.10, 490.50, 569, 647.40, 725.90, 799.50, 868.20, 941.70, 1020.20,
    1098.70, 1245.80, 1373.40, 1491.10, 1618.60, 1755.90, 1893.30, 2020.80, 2315.10, 2472, 2766.30, 3099.90, 3413.80, 4139.70,
    4934.30, 5768.10]kN
ρ_SK75 = 0.975 * 0.001kg/(10mm)^3
σ_TS_SK75 = 3600N/mm^2

BACKCOLOR = color_with_lumin(PALETTE[8], 90);
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_17.png"))
    background(BACKCOLOR)
    sethue(PALETTE[5])
    set_scale_sketch(sum(ROPE_D) * 1.3, HE)
end
restart()

# Horizontal positions per rope
pts = let
    posx = [0.0]mm
    pop!(posx)
    for (i,d) in enumerate(ROPE_D)
        p = if i == 1
                -0.45 * WI * scale_pt_to_unit(m)
        else
            posx[i-1] + max(60.0mm, d*1.1)
        end
        push!(posx, p)
    end
    pts = map(x-> O +(x, 700.0mm), posx)
end
let
    circle.(pts, ROPE_D ./ 2,:fill)
    labs = string.(ROPE_D)
    text.(labs, O + (0, EM) .+pts, -60.0°)
end

circlearea(d) = π / 4 * d^2
circleradius(A) = sqrt(A/π)
# Plot packed area
let
    A_rope = ROPE_WEIGHT ./ ρ_SK75 .|> mm^2
    sethue(color_with_lumin(PALETTE[5], 80))
    circle.(pts , circleradius.(A_rope),:fill)
end

"Relative filled areas, actual to full circle"
area_rel =  (ROPE_WEIGHT ./ ρ_SK75 .|> mm^2) ./ circlearea.(ROPE_D)
# Plot relative filled areas
setline(3PT)
bottoms = map(pt -> pt + (0.0mm, -600.0mm), pts)
tops = map(pt -> pt + (0.0mm, 300.0mm), bottoms)
line.(bottoms, tops, :stroke)
stri = "Fill factors and approximation (green line)"
text(stri, tops[1] + (0.0, -0.5EM))

sethue(PALETTE[4])
map(tops, bottoms, area_rel) do t, b, rel
    line(b, b + (t - b) * rel, :stroke)
end
# lsq method applied
fill_factor(d) = (-0.003*d/mm + 0.853)
area_filled(d) = fill_factor(d)* circlearea(d)

# Plot appproximated area
sethue(PALETTE[3])
setopacity(0.5)
map(tops, bottoms, ROPE_D) do t, b, d
    rel = fill_factor(d)
    top = b + (t - b) * rel
    line(top + (-0.5EM, 0.0), top + (0.5EM, 0.0), :stroke)
end
setopacity(1.0)

ideal_spin_strength(d) = area_filled(d) * σ_TS_SK75 |> kN
spinfactors = ROPE_MINIMUM_BREAK./ ideal_spin_strength.(ROPE_D)
#lsq method applied
spin_factor(d) = 5*10^-5 * (d/mm)^2 + -0.0059*(d/mm) + 0.5504
breaking_strength(d) = spinfactor(d) * ideal_spin_strength(d)

# Plot spin factors
setline(3PT)
sethue(color_with_lumin(PALETTE[5], 80))
bottoms = map(pt -> pt + (0.0mm, -1200.0mm), pts)
tops = map(pt -> pt + (0.0mm, 300.0mm), bottoms)
line.(bottoms, tops, :stroke)
stri = "Spin factors and approximation (green line)"
text(stri, tops[1] + (0.0, -0.5EM))

sethue(PALETTE[4])
map(tops, bottoms, spinfactors) do t, b, rel
    line(b, b + (t - b) * rel, :stroke)
end

# Plot appproximated spin factor
sethue(PALETTE[3])
setopacity(0.5)
map(tops, bottoms, ROPE_D) do t, b, d
    rel = spin_factor(d)
    top = b + (t - b) * rel
    line(top + (-0.5EM, 0.0), top + (0.5EM, 0.0), :stroke)
end
setopacity(1.0)

finish()
# Reset default
set_scale_sketch()

nothing
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(6, 6, 2, 1, 1, Line(2, 3)),
  :properties => [
    5.0e9 6.0e4;
    1.0e9  2.0e4],
  :etype => [1, 1, 1, 2, 2, 2],
  :x_coords => [0.0, 6.0, 6.0, 12.0, 12.0, 14.0],
  :y_coords => [0.0, 0.0, -4.0, 0.0, -5.0, 0.0],
  :g_num => [
    1 2 4 3 3 5;
    2 4 6 2 4 4],
  :support => [
    (1, [0 0 1]),
    (3, [0 0 0]),
    (5, [0 0 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -60.0 -60.0]),
    (2, [0.0 -180.0 -80.0]),
    (4, [0.0 -140.0 133.33]),
    (6, [0.0 -20.0 6.67])
    ],
  :penalty => 1e19,
  :eq_nodal_forces_and_moments => [
    (1, [0.0 -60.0 -60.0 0.0 -60.0 60.0]),
    (2, [0.0 -120.0 -140.0 0.0 -120.0 140.0]),
    (3, [0.0 -20.0 -6.67 0.0 -20.0 6.67])
  ]
)

data |> display
println()

@time fem, dis_df, fm_df = p44(data)
println()

display(dis_df)
println()
display(fm_df)
println()

if VERSION.minor < 6
  using Plots
  gr(size=(400,600))

  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3)
  titles = ["p44.1 rotations", "p44.1 y shear force", "p44.1 z moment"]
  moms = vcat(
    convert(Array, fm_df[:, :z1_Moment]),
    convert(Array, fm_df[:, :z2_Moment])[end]
  )
  fors = vcat(
    convert(Array, fm_df[:, :y1_Force]),
    convert(Array, fm_df[:, :y2_Force])[end]
  )
  x_coords = data[:x_coords]

  p[1] = plot(fem.displacements[3,:], ylim=(-0.002, 0.002),
    xlabel="node", ylabel="rotation [radians]", color=:red,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[1], leg=false)
  p[2] = plot(fors, lab="y Shear force", ylim=(-150.0, 250.0),
    xlabel="node", ylabel="shear force [N]", color=:blue,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[2], leg=false)
  p[3] = plot(moms, lab="z Moment", ylim=(-10.0, 150.0),
    xlabel="node", ylabel="z moment [Nm]", color=:green,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[3], leg=false)

  plot(p..., layout=(3, 1))
  savefig(ProjDir*"/Ex44.1.png")

end


text:
export Particles

module Particles

using ..Ranger

const MAX_PARTICLE_SPEED = 20.0
const MAX_PARTICLE_LIFETIME = 1.0
const MAX_PARTICLES = 50
const MAX_PARTICE_SIZE = 10.0

include("abstracts.jl")
include("activator_360.jl")
include("particle.jl")
include("particle_system.jl")

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


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

mutable struct ManagementPolicyDefinition <: SwaggerModel
    actions::Any # spec type: Union{ Nothing, ManagementPolicyAction } # spec name: actions
    filters::Any # spec type: Union{ Nothing, ManagementPolicyFilter } # spec name: filters

    function ManagementPolicyDefinition(;actions=nothing, filters=nothing)
        o = new()
        validate_property(ManagementPolicyDefinition, Symbol("actions"), actions)
        setfield!(o, Symbol("actions"), actions)
        validate_property(ManagementPolicyDefinition, Symbol("filters"), filters)
        setfield!(o, Symbol("filters"), filters)
        o
    end
end # type ManagementPolicyDefinition

const _property_map_ManagementPolicyDefinition = Dict{Symbol,Symbol}(Symbol("actions")=>Symbol("actions"), Symbol("filters")=>Symbol("filters"))
const _property_types_ManagementPolicyDefinition = Dict{Symbol,String}(Symbol("actions")=>"ManagementPolicyAction", Symbol("filters")=>"ManagementPolicyFilter")
Base.propertynames(::Type{ ManagementPolicyDefinition }) = collect(keys(_property_map_ManagementPolicyDefinition))
Swagger.property_type(::Type{ ManagementPolicyDefinition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagementPolicyDefinition[name]))}
Swagger.field_name(::Type{ ManagementPolicyDefinition }, property_name::Symbol) =  _property_map_ManagementPolicyDefinition[property_name]

function check_required(o::ManagementPolicyDefinition)
    (getproperty(o, Symbol("actions")) === nothing) && (return false)
    true
end

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


text:
using Expectations, Distributions

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test, Compat.Random
else
    using Test, Random, Statistics, LinearAlgebra
end

@testset "Iterable distributions" begin include("iterable.jl") end


text:
abstract FiniteHorizonAgent
nS(m::FiniteHorizonAgent) = size(m.policy, 1)
horizon(m::FiniteHorizonAgent) = size(m.policy, 2)
nA(m::FiniteHorizonAgent) = size(m.rewards, 2)

sample_action(m::FiniteHorizonAgent, s, t) = m.policy[s, t]
function observe!(m::FiniteHorizonAgent, s, a, t, r, sn)
    m.ns[sn, s, a, t] += 1
    m.n[s, a, t] +=  1
    if !m.rewards_known
        m.rewards[s, a, t] += r
    end
end

include("psrl.jl")
include("ubev.jl")
include("mbie.jl")
include("mormax.jl")
include("ucrl2.jl")
include("medianpac.jl")
include("delayedql.jl")
include("oim.jl")
include("ucfh.jl")


text:
###################################################################################################
#                                        Production Conditions
###################################################################################################
function can_wait(actr)
    c1(actr) = isempty(actr.visual_location.buffer)
    c2(actr) = isempty(actr.visual.buffer)
    c3(actr) = !actr.visual.state.busy
    c4(actr) = !actr.motor.state.busy
    return c1, c2, c3, c4
end

function can_find(actr)
    c1(actr) = isempty(actr.visual_location.buffer)
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end

function can_attend(actr)
    c1(actr) = !isempty(actr.visual_location.buffer)
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end
###################################################################################################
#                                        Production Actions
###################################################################################################
function find_action(actr, args...)
    vo = actr.visicon[1]
    chunk = vo_to_chunk(actr, vo)
    add_to_buffer!(actr.visual_location, chunk)
    return nothing
end

function attend_action(actr, task, args...)
    buffer = actr.visual_location.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual_location)
    attending!(actr, chunk, task)
    return nothing
end

function respond_action(actr, task)
    clear_buffer!(actr.visual)
    key = "sb"
    responding!(actr, task, key)
    return nothing
end