text:
##########
# @d_str #
##########
const ktpat = r"\|.*?\>"
const brpat = r"\<.*?\|"
const inpat = r"\<[^\|]*\|[^\|]*\>"

ktrep(str) = "ket("*str[2:end-1]*")"
brrep(str) = "bra("*str[2:end-1]*")"

function inrep(str)
    i = split(str, '|')
    return "(bra("*i[1][2:end]*")*ket("*i[2][1:end-1]*"))"
end

function prune_dirac(str)
    return replace(replace(replace(str, inpat => inrep), brpat => brrep), ktpat => ktrep)
end

macro d_str(str)
    return esc(Meta.parse(prune_dirac(str)))
end

macro d_mstr(str)
    return esc(quote
        local s;
        for s in filter(x -> ! isempty(x), split(strip($str), '\n'))
            eval(Meta.parse(QuDirac.prune_dirac(s)))
        end
    end)
end

###################################
# Definition String Parsing Utils #
###################################
struct OpDefStr
    op_name::AbstractString
    label_args::AbstractString
    lhs_type::AbstractString
    rhs::AbstractString
end

function OpDefStr(str::AbstractString)
    str = rm_whspace(str)
    left, right = @compat(split(str, '=', limit=2))

    if occursin(ktpat, left)
        lhs_type = "Ket"
    elseif occursin(brpat, left)
        lhs_type = "Bra"
    else
        error("Couldn't detect whether operator is acting on Bra or Ket...")
    end

    if lhs_type == "Ket"
        op_name, label_args = leftside_ket(left)
    else # type_str == "Bra"
        op_name, label_args = leftside_bra(left)
    end

    return OpDefStr(op_name, label_args, lhs_type, prune_dirac(right))
end

struct OpDefExpr
    op_sym::Symbol
    label_args::Union{Symbol, Expr}
    lhs_type::Symbol
    rhs::Union{Symbol, Expr}
end

OpDefExpr(ods::OpDefStr) = OpDefExpr(Symbol(ods.op_name), Meta.parse(ods.label_args), Symbol(ods.lhs_type), Meta.parse(ods.rhs))

rm_whspace(str) = join(split(str, r"\s"))

function leftside_ket(str)
    op_name, label_args = @compat(split(str, '|', limit=2))
    if last(label_args) == '>'
        return op_name, label_args[1:end-1]
    else
        error("Left side of definition string is malformed; couldn't find terminating '>' for the Ket.")
    end
end

function leftside_bra(str)
    label_args, op_name = @compat(split(str, '|', limit=2))
    if first(label_args) == '<'
        return op_name, label_args[2:end]
    else
        error("Left side of definition string is malformed; couldn't find initial '<' for the Bra.")
    end
end

###########
# @def_op #
###########

macro def_op(str)
    result_expr = def_op_expr(OpDefStr(str))
    return esc(result_expr)
end

function def_op_expr(ods::OpDefStr)
    odex = OpDefExpr(ods)

    if isa(odex.label_args, Expr)
        func_on_label_def = quote
            local func_on_args($(odex.label_args.args...)) = $(odex.rhs)
            local func_on_label(label::StateLabel) = func_on_args(label...)
        end
    else # isa(label_expr, Symbol)
        func_on_label_def = quote
            local func_on_args($(odex.label_args)) = $(odex.rhs)
            local func_on_label(label::StateLabel) = func_on_args(first(label))
        end
    end

    coeff = coeff_sym(odex)

    result = quote
        $func_on_label_def

        local func_on_pair;
        function func_on_pair(pair::Pair)
          label, c = pair
          return $(coeff) * func_on_label(label)
        end

        $(odex.op_sym)(state::$(odex.lhs_type)) = sum(func_on_pair, QuDirac.dict(state))
    end

    return result
end

coeff_sym(odex::OpDefExpr) = odex.lhs_type == :Ket ? :c : :(c')

###########
# @rep_op #
###########

macro rep_op(str, bases...)
    result_expr = rep_op_expr(OpDefStr(str), build_prod_basis(bases...))
    return esc(result_expr)
end

build_prod_basis(basis) = basis
build_prod_basis(first, bases...) = :(Iterators.product($first, $(bases...)))

function rep_op_expr(ods::OpDefStr, basis)
    odex = OpDefExpr(ods)
    if odex.lhs_type == :Ket
        return gen_ket_repr_expr(odex, basis)
    else
        return gen_bra_repr_expr(odex, basis)
    end
    return ex
end

function gen_ket_repr_expr(odex::OpDefExpr, basis)
    if isa(odex.label_args, Expr)
        ex = quote
            local func_on_args($(odex.label_args.args...)) = $(odex.rhs) * bra($(odex.label_args.args...))
            $(odex.op_sym) = sum(args->func_on_args(args...), $basis)
        end
    else
        ex = quote
            $(odex.op_sym) = sum($(odex.label_args) -> $(odex.rhs) * bra($(odex.label_args)), $basis)
        end
    end
    return ex
end

function gen_bra_repr_expr(odex::OpDefExpr, basis)
    if isa(odex.label_args, Expr)
        ex = quote
            local func_on_args($(odex.label_args.args...)) = ket($(odex.label_args.args...)) * $(odex.rhs)
            $(odex.op_sym) = sum(args->func_on_args(args...), $basis)
        end
    else
        ex = quote
            $(odex.op_sym) = sum($(odex.label_args) -> ket($(odex.label_args)) * $(odex.rhs), $basis)
        end
    end
    return ex
end

export @d_str,
       @d_mstr,
       @def_op,
       @rep_op


text:
export SpecialEuclidean

"""
$(TYPEDEF)

Rigid transform between two Pose2's, assuming (x,y,theta).

DevNotes
- Maybe with Manifolds.jl, `{T <: IIF.SamplableBelief, S, R, P}`

Related

[`Pose3Pose3`](@ref), [`Point2Point2`](@ref), [`MutablePose2Pose2Gaussian`](@ref), [`DynPose2`](@ref), [`InertialPose3`](ref)
"""
struct Pose2Pose2{T <: IIF.SamplableBelief} <: IIF.AbstractManifoldMinimize
  z::T
end
# convenience and default constructor
Pose2Pose2() = Pose2Pose2(MvNormal(Diagonal([1.0; 1.0; 1.0])))

DFG.getManifold(::Pose2Pose2) = Manifolds.SpecialEuclidean(2)

Pose2Pose2(::UniformScaling) = Pose2Pose2() # MvNormal(zeros(3),LinearAlgebra.diagm([1.0;1.0;1.0])) )

function getSample(cf::CalcFactor{<:Pose2Pose2})

  M = getManifold(cf.factor)
  ϵ = getPointIdentity(Pose2)

  X = sampleTangent(M, cf.factor.z, ϵ)
  return X
end

# Assumes X is a tangent vector
function (cf::CalcFactor{<:Pose2Pose2})(X, p, q)
    @assert X isa ProductRepr "Pose2Pose2 expects measurement sample X to be a Manifolds tangent vector, not coordinate or point representation.  Got X=$X"
    M = getManifold(Pose2)
    q̂ = Manifolds.compose(M, p, exp(M, identity_element(M, p), X)) #for groups
    #TODO allocalte for vee! see Manifolds #412, fix for AD
    Xc = zeros(3)
    vee!(M, Xc, q, log(M, q, q̂))
    return Xc
end

# NOTE, serialization support -- will be reduced to macro in future
# ------------------------------------

"""
$(TYPEDEF)
"""
mutable struct PackedPose2Pose2  <: IIF.PackedInferenceType
  datastr::String
  # PackedPose2Pose2() = new()
  # PackedPose2Pose2(x::String) = new(x)
end
function convert(::Type{Pose2Pose2}, d::PackedPose2Pose2)
  return Pose2Pose2(convert(SamplableBelief, d.datastr))
end
function convert(::Type{PackedPose2Pose2}, d::Pose2Pose2)
  return PackedPose2Pose2(convert(PackedSamplableBelief, d.z))
end

# FIXME, rather have separate compareDensity functions
function compare(a::Pose2Pose2,b::Pose2Pose2; tol::Float64=1e-10)
  return compareDensity(a.z, b.z)
  # TP = true
  # TP = TP && norm(a.z.μ-b.z.μ) < (tol + 1e-5)
  # TP = TP && norm(a.z.Σ.mat-b.z.Σ.mat) < tol
  # return TP
end

## Deprecated
# function Pose2Pose2(mean::Array{Float64,1}, cov::Array{Float64,2}, w::Vector{Float64})
#   @warn "Pose2Pose2(mu,cov,w) is deprecated in favor of Pose2Pose2(T(...)) -- use for example Pose2Pose2(MvNormal(mu, cov))"
#   Pose2Pose2(MvNormal(mean, cov))
# end

ext: jl


text:
# data structures graphical models

export FactorGraph, jsondata

alphabets = [:spin, :boolean, :integer, :integer_pos, :real, :real_pos]

mutable struct FactorGraph{T <: Real}
    order::Int
    varible_count::Int
    alphabet::Symbol
    terms::Dict{Tuple,T} # TODO, would be nice to have a stronger tuple type here
    variable_names::Union{Vector{String}, Nothing}
    #FactorGraph(a,b,c,d,e) = check_model_data(a,b,c,d,e) ? new(a,b,c,d,e) : error("generic init problem")
end

FactorGraph(order::Int, varible_count::Int, alphabet::Symbol, terms::Dict{Tuple,T}) where T <: Real = FactorGraph{T}(order, varible_count, alphabet, terms, nothing)
FactorGraph(matrix::Array{T,2}) where T <: Real = convert(FactorGraph{T}, matrix)
FactorGraph(dict::Dict{Tuple,T}) where T <: Real = convert(FactorGraph{T}, dict)
FactorGraph(list::Array{Any,1}) = convert(FactorGraph, list)

function check_model_data(order::Int, varible_count::Int, alphabet::Symbol, terms::Dict{Tuple,T}, variable_names::Union{Vector{String}, Nothing}) where T <: Real
    if !in(alphabet, alphabets)
        error("alphabet $(alphabet) is not supported")
        return false
    end
    if variable_names != nothing && length(variable_names) != varible_count
        error("expected $(varible_count) but only given $(length(variable_names))")
        return false
    end
    for (k,v) in terms
        if length(k) > order
            error("a term has $(length(k)) indices but should have $(order) indices")
            return false
        end
        for (i,index) in enumerate(k)
            #println(i," ",index)
            if index < 1 || index > varible_count
                error("a term has an index of $(index) but it should be in the range of 1:$(varible_count)")
                return false
            end
            #=
            # TODO see when this should be enforced
            if i > 1
                if k[i-1] > index
                    error("the term $(k) does not have ascending indices")
                end
            end
            =#
        end
    end
    return true
end

function Base.show(io::IO, gm::FactorGraph)
    println(io, "alphabet: ", gm.alphabet)
    println(io, "vars: ", gm.varible_count)
    if gm.variable_names != nothing
        println(io, "variable names: ")
        println(io, "  ", get(gm.variable_names))
    end

    println(io, "terms: $(length(gm.terms))")
    for k in sort(collect(keys(gm.terms)), by=(x)->(length(x),x))
        println("  ", k, " => ", gm.terms[k])
    end
end

function jsondata(gm::FactorGraph{T}) where T <: Real
    data = []
    for k in sort(collect(keys(gm.terms)), by=(x)->(length(x),x))
        push!(data, Dict("term" => k, "weight" => gm.terms[k]))
    end
    return data
end

Base.iterate(gm::FactorGraph, kwargs...) = Base.iterate(gm.terms, kwargs...)

Base.length(gm::FactorGraph) = length(gm.terms)
#Base.size(gm::FactorGraph, a...) = size(gm.terms, a...)

Base.getindex(gm::FactorGraph, i) = gm.terms[i]
Base.keys(gm::FactorGraph) = keys(gm.terms)

function diag_keys(gm::FactorGraph)
    dkeys = Tuple[]
    for i in 1:gm.varible_count
        key = diag_key(gm, i)
        if key in keys(gm.terms)
            push!(dkeys, key)
        end
    end
    return sort(dkeys)
end

diag_key(gm::FactorGraph, i::Int) = tuple(fill(i, gm.order)...)

#Base.diag(gm::FactorGraph{T}) where T <: Real = [ get(gm.terms, diag_key(gm, i), zero(T)) for i in 1:gm.varible_count ]

#Base.DataFmt.writecsv(io, gm::FactorGraph{T}, args...; kwargs...) where T <: Real = writecsv(io, convert(Array{T,2}, gm), args...; kwargs...)

Base.convert(::Type{FactorGraph}, m::Array{T,2}) where T <: Real = convert(FactorGraph{T}, m)
function Base.convert(::Type{FactorGraph{T}}, m::Array{T,2}) where T <: Real
    @assert size(m,1) == size(m,2) #check matrix is square

    @info "assuming spin alphabet"
    alphabet = :spin
    varible_count = size(m,1)

    terms = Dict{Tuple,T}()

    for key in permutations(1:varible_count, 1)
        weight = m[key..., key...]
        if !isapprox(weight, 0.0)
            terms[key] = weight
        end
    end

    for key in permutations(1:varible_count, 2)
        weight = m[key...]
        if !isapprox(weight, 0.0)
            terms[key] = weight
        end

        rev = reverse(key)
        if !isapprox(m[rev...], 0.0) && !isapprox(m[key...], m[rev...])
            delta = abs(m[key...] - m[rev...])
            warn("values at $(key) and $(rev) differ by $(delta), only $(key) will be used")
        end
    end

    return FactorGraph(2, varible_count, alphabet, terms)
end

function Base.convert(::Type{Array{T,2}}, gm::FactorGraph{T}) where T <: Real
    if gm.order != 2
        error("cannot convert a FactorGraph of order $(gm.order) to a matrix")
    end

    matrix = zeros(gm.varible_count, gm.varible_count)
    for (k,v) in gm
        if length(k) == 1
            matrix[k..., k...] = v
        else
            matrix[k...] = v
            r = reverse(k)
            matrix[r...] = v
        end
    end

    return matrix
end

Base.convert(::Type{Dict}, m::Array{T,2}) where T <: Real = convert(Dict{Tuple,T}, m)
function Base.convert(::Type{Dict{Tuple,T}}, m::Array{T,2}) where T <: Real
    @assert size(m,1) == size(m,2) #check matrix is square

    varible_count = size(m,1)

    terms = Dict{Tuple,T}()

    for key in permutations(1:varible_count, 1)
        weight = m[key..., key...]
        if !isapprox(weight, 0.0)
            terms[key] = weight
        end
    end

    for key in permutations(1:varible_count, 2, asymmetric=true)
        if key[1] != key[2]
            weight = m[key...]
            if !isapprox(weight, 0.0)
                terms[key] = weight
            end
        end
    end

    return terms
end

function Base.convert(::Type{FactorGraph}, list::Array{Any,1})
    @info "assuming spin alphabet"
    alphabet = :spin

    max_variable = 0
    max_order = 0
    terms = Dict{Tuple,Float64}()

    for item in list
        term = item["term"]
        weight = item["weight"]
        terms[tuple(term...)] = weight

        @assert minimum(term) > 0
        max_order = max(max_order, length(term))
        max_variable = max(max_variable, maximum(term))
    end

    @info "dectected $(max_variable) variables with order $(max_order)"

    return FactorGraph(max_order, max_variable, alphabet, terms)
end

Base.convert(::Type{FactorGraph}, dict::Dict{Tuple,T}) where T <: Real = convert(FactorGraph{T}, dict)
function Base.convert(::Type{FactorGraph{T}}, dict::Dict{Tuple,T}) where T <: Real
    @info "assuming spin alphabet"
    alphabet = :spin

    max_variable = 0
    max_order = 0
    for (term,weight) in dict
        @assert minimum(term) > 0
        max_order = max(max_order, length(term))
        max_variable = max(max_variable, maximum(term))
    end

    @info "dectected $(max_variable) variables with order $(max_order)"

    return FactorGraph(max_order, max_variable, alphabet, dict)
end

permutations(items, order::Int; asymmetric::Bool = false) = sort(permutations([], items, order, asymmetric))

function permutations(partial_perm::Array{Any,1}, items, order::Int, asymmetric::Bool)
    if order == 0
        return [tuple(partial_perm...)]
    else
        perms = []
        for item in items
            if !asymmetric && length(partial_perm) > 0
                if partial_perm[end] >= item
                    continue
                end
            end
            perm = permutations(vcat(partial_perm, item), items, order-1, asymmetric)
            append!(perms, perm)
        end
        return perms
    end
end


text:
import Dates
using Test
import Humanize

@testset "Humanize" begin

@testset "datasize" begin
    VALUES = [300, 3000, 3000000, 3000000000, 3000000000000, 1e26, 3141592]
    MULTIPLIERS = Dict(:dec => 1.5, :bin => 2.1, :gnu => 1.6)
    RESULTS = Dict(
        :dec => ["450.000 B", "4.500 kB", "4.500 MB", "4.500 GB", "4.500 TB", "150.000 YB", "4.712 MB"],
        :bin => ["630.000 B", "6.152 KiB", "6.008 MiB", "5.867 GiB", "5.730 TiB", "173.708 YiB", "6.292 MiB"],
        :gnu => ["480.000", "4.688K", "4.578M", "4.470G", "4.366T", "132.349Y", "4.794M"]
    )
    @testset "$style" for style in (:dec, :bin, :gnu)
        for (value, result) in zip(VALUES, RESULTS[style])
            size = value * MULTIPLIERS[style]
            @test Humanize.datasize(size, style=style, format="%.3f") == result
        end
    end
end  # testset datasize.

@testset "timedelta" begin
    # Years, months, days, hours, minutes, seconds, result.
    DATA = [((0, 0, 0, 0, 0, 0), "a moment"),
            ((0, 0, 0, 0, 0, 1), "a second"),
            ((0, 0, 0, 0, 0, 30), "30 seconds"),
            ((0, 0, 0, 0, 1, 30), "a minute"),
            ((0, 0, 0, 0, 2, 0), "2 minutes"),
            ((0, 0, 0, 1, 30, 30), "an hour"),
            ((0, 0, 0, 23, 50, 50), "23 hours"),
            ((0, 0, 1, 0, 0, 0), "a day"),
            ((0, 0, 500, 0, 0, 0), "1 year, 4 months"),
            ((0, 0, 365 * 2 + 35, 0, 0, 0), "2 years"),
            ((0, 0, 10000, 0, 0, 0), "27 years"),
            ((0, 0, 365 + 30, 0, 0, 0), "1 year, 1 month"),
            ((0, 0, 365 + 4, 0, 0, 0), "a year"),
            ((0, 0, 35, 0, 0, 0), "a month"),
            ((0, 0, 65, 0, 0, 0), "2 months"),
            ((0, 0, 9, 0, 0, 0), "9 days"),
            ((0, 0, 365, 0, 0, 0), "a year")]

    @testset "datetime diff $output" for (inputs, output) in DATA
            base_datetime = Dates.DateTime(2014, 1, 1, 0, 0, 0)
        new_datetime = Dates.Year(inputs[1]) + base_datetime
        new_datetime += Dates.Month(inputs[2])
        new_datetime += Dates.Day(inputs[3])
        new_datetime += Dates.Hour(inputs[4])
        new_datetime += Dates.Minute(inputs[5])
        new_datetime += Dates.Second(inputs[6])
        @test Humanize.timedelta(new_datetime - base_datetime) == output
    end

    @testset "date diff $output" for (inputs, output) in DATA
        sum(inputs[1:3]) == 0 && continue  # Hour-scale or less.
        base_date = Dates.Date(2014, 1, 1)
        new_date = Dates.Year(inputs[1]) + base_date
        new_date += Dates.Month(inputs[2])
        new_date += Dates.Day(inputs[3])
        @test Humanize.timedelta(new_date - base_date) == output
    end
end  # testset timedelta.

@testset "digitsep" begin
    DATA = ((1, "1"),
            (12, "12"),
            (123, "123"),
            (1234, "1,234"),
            (12345, "12,345"),
            (123456, "123,456"),
            (1234567, "1,234,567"),
            (12345678, "12,345,678"),
            (-1, "-1"),
            (-12, "-12"),
            (-123, "-123"),
            (-1234, "-1,234"),
            (-12345, "-12,345"),
            (-123456, "-123,456"),
            (-1234567, "-1,234,567"),
            (-12345678, "-12,345,678"))
    @testset "digitsep $output" for (input, output) in DATA
        @test Humanize.digitsep(input) == output
    end
end  # testset digitsep.

end  # testset Humanize.


text:
using BosonSampling
using Test
using LinearAlgebra

@testset "BosonSampling.jl" begin

    @testset "scattering matrix" begin

        input_state = [2,0]
	    output_state = [1,1]
	    U = matrix_test(2)
	    @test scattering_matrix(U, input_state, output_state) == [1 2; 1 2]

	    input_state = [1,1]
	    output_state = [0,2]
	    U = matrix_test(2)
	    @test scattering_matrix(U, input_state, output_state) == [2 2; 4 4]

	end

	@testset "process_probability: HOM" begin

		@test process_probability(fourier_matrix(2), [1,1], [0,2]) ≈ 0.5 atol = 1e-8
	    @test process_probability(fourier_matrix(2), [1,1], [2,0]) ≈ 0.5 atol = 1e-8
	    @test process_probability(fourier_matrix(2), [1,1], [1,1]) ≈ 0. atol = 1e-8

		S_bosonic = ones(2,2)
		S_dist = [1 0; 0 1]

		@test process_probability_partial(fourier_matrix(2), S_bosonic, [1,1], [0,2]) ≈ 0.5 atol = 1e-8
		@test process_probability_partial(fourier_matrix(2), S_bosonic, [1,1], [2,0]) ≈ 0.5 atol = 1e-8
		@test process_probability_partial(fourier_matrix(2), S_bosonic, [1,1], [1,1]) ≈ 0. atol = 1e-8
		@test process_probability_partial(fourier_matrix(2), S_dist, [1,1], [0,2]) ≈ 0.25 atol = 1e-8
		@test process_probability_partial(fourier_matrix(2), S_dist, [1,1], [2,0]) ≈ 0.25 atol = 1e-8
		@test process_probability_partial(fourier_matrix(2), S_dist, [1,1], [1,1]) ≈ 0.5 atol = 1e-8

	end

	@testset "probability distribution of photons in modes : HOM" begin

		U = fourier_matrix(2)

		occupancy_vector = [1, 0]
		@test proba_partition_bosonic(U = U, occupancy_vector = occupancy_vector) ≈ [0.5,0,0.5] atol = 1e-5

		occupancy_vector = [0, 1]
		@test proba_partition_bosonic(U = U, occupancy_vector = occupancy_vector) ≈ [0.5,0,0.5] atol = 1e-5

		occupancy_vector = [1, 1]
		@test proba_partition_bosonic(U = U, occupancy_vector = occupancy_vector) ≈ [0,0,1.] atol = 1e-5

	end

	@testset "probability distribution of distinguishable photons in modes : HOM" begin

	    U = fourier_matrix(2)
	    part = [1,0]
	    @test proba_partition_distinguishable(occupancy_vector = part, U = U) ≈ [0.25,0.5,0.25] atol = 1e-8

	end

	@testset "partial distinguishability partitions" begin

		m = 10
		n = 4
		input_state = zeros(Int, m)
		occupancy_vector = zeros(Int, m)
		for i=1:n
			input_state[i] = 1
		end
		occupancy_vector[1] = 1

		U = fourier_matrix(10)

		@test proba_partition_partial(U = U, S = ones(n, n), occupancy_vector = occupancy_vector, input_state = input_state) == proba_partition_bosonic(U = U, occupancy_vector = occupancy_vector, input_state = input_state)
	    @test proba_partition_partial(U = U,  S = Matrix{Float64}(I, n, n), occupancy_vector = occupancy_vector, input_state = input_state) ≈ proba_partition_distinguishable(occupancy_vector = occupancy_vector, U = U, input_state = input_state)

	end

	@testset "theoretical_distribution" begin

		n = 3
		occupation = ModeOccupation(random_occupancy(n,n))
		i = Input{Bosonic}(occupation)
		interf = Fourier(i.r.m)
		p_theo = theoretical_distribution(input=i, interf=interf)

		@test sum(p_theo) ≈ 1 atol=1e-9
		@test all(p -> p>=0, p_theo)

		output_events = output_mode_occupation(n,n)
	    for i = 1:length(output_events)
	        if check_suppression_law(output_events[i])
	            @test p_theo[i] ≈ 0 atol=1e-6
	        end
	    end

	end

	@testset "noisy statistics" begin

		n = 3
		occupation = ModeOccupation(random_occupancy(n,n))
    	i = Input{OneParameterInterpolation}(occupation, 0.5)
	    interf = Fourier(i.r.m)
	    res = noisy_distribution(input=i, reflectivity=0.5, interf=interf)

		p_exact = res[1]
		p_approx = res[2]
		p_samp = res[3]

		@test sum(p_exact) ≈ 1 atol=1e-9
		@test sum(p_approx) ≈ 1 atol=1e-9
		@test sum(p_samp) ≈ 1 atol=1e-9
		@test all(p->p>=0, p_exact)
		@test all(p->p>=0, p_approx)
		@test all(p->p>=0, p_samp)

		i = Input{Bosonic}(ModeOccupation(random_occupancy(n,n)))
		res = noisy_distribution(input=i, reflectivity=0.999, interf=interf, approx=false, samp=false)
		p_exact = res[1]
		output_events = output_mode_occupation(n,n)
		for i = 1:length(output_events)
            if check_suppression_law(output_events[i])
		        @test p_exact[i] ≈ 0 atol=1e-6
		    end
		end

	end

	@testset "suppression law boson samplers" begin

		n = 3
	    interf = Fourier(n)

		input_clifford_sampler = Input{Bosonic}(first_modes(n,n))
	   	input_noisy_sampler = Input{OneParameterInterpolation}(first_modes(n,n), 1.0)

	    out_clifford_sampler = cliffords_sampler(input=input_clifford_sampler, interf=interf)
		out_noisy_sampler = noisy_sampler(input=input_noisy_sampler, reflectivity=1.0, interf=interf)

	    @test !check_suppression_law(out_clifford_sampler)
		@test !check_suppression_law(out_noisy_sampler)

	end

	@testset "examples usage" begin
		@test_nowarn include("example_usage.jl")

	end

end


text:
using NARMAX
using Test

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

max_stars_repo_name: wmkouw/NARMAX

max_issues_repo_name: wmkouw/NARMAX

max_forks_repo_name: wmkouw/NARMAX


text:
using Clang
using LibGit2
using Dates

const C_header = "src/gswteos-10.h"
const C_URI    = "https://github.com/TEOS-10/GSW-C/"
const Rfn      = "../README.md"
const Repodir  = "src/"
const RMnew    = "../README.md.new"

const LIBCLANG_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include", "clang-c") |> normpath

const RMIntros  = [
"# GSW.jl\n\nGSW.jl is a Julia wrapper for ",
"[GSW-C](https://github.com/TEOS-10/GSW-C)", # Link
", which is the C implementation of the Thermodynamic Equation of Seawater 2010 (TEOS-10).\n\n",
]

function get_gitsrc(uri, repodir)
  if isdir(repodir)
    gitrepo = LibGit2.GitRepo(repodir)
  else
    gitrepo = LibGit2.clone(uri, repodir)
  end
  ref   = LibGit2.head(gitrepo)
  hash  = LibGit2.GitHash(ref)
  lhash = string(hash)
  shash = string(LibGit2.GitShortHash(hash, 7))

  return (lhash, shash)
end

function ff(x, y)
  return occursin(r"gsw", y) #ismatch(r"gsw", y)
end

lhash, shash = get_gitsrc(C_URI, Repodir)

uri_hash = C_URI * "tree/" * lhash
tm       = Dates.now(Dates.UTC)
tmstr    = Dates.format(tm, "Y/mm/dd HH:MM")

RMIntros[2] = "[GSW-C#$(shash)]($(uri_hash)) (generated by `utils/generate.jl` at $(tmstr))"

#readme = readstring(Rfn)
readme = read(Rfn, String)
mtchs  = match(r"(\# GSW.jl.*?)(^\#\#.*)"ms, readme)

readme_txt = join(RMIntros, "") * mtchs.captures[2]

#write(RMnew, readme_txt)
#quit()

hdrs = read(C_header, String)

context = Clang.init(
output_file = "../$(Repodir)tgen_gswteos10.jl",
common_file = "../$(Repodir)tgen_gswteos_h.jl",
headers = [C_header],
clang_includes = vcat(LIBCLANG_INCLUDE, CLANG_INCLUDE),
#clang_includes = includes,
clang_diagnostics = true,
header_wrapped=(x,y)->ff(x,y),#(@show x; ismatch(r"gsw", y)),
header_library=x->"libgswteos",
# rewriter = rewriter
)

#context.options = Clang.wrap!() #.InternalOptions(true, true)

run(context)

#@show hdrs


text:
type Encoder
    verbose::Bool
    encoder_functions
    encodes_to_string
    emitter::Emitter

    Encoder(io, verbose=false) =
		new(verbose,
		    Dict{DataType,Function}(),
		    Dict{DataType,Bool}(),
		    make_emitter(io, verbose))
end

function add_encoder(e::Encoder, t::DataType, f::Function, encodes_to_string::Bool)
    e.encoder_functions[t] = f
    e.encodes_to_string[t] = encodes_to_string
end

function encode(e::Encoder, x::Any, askey=false)
    if haskey(e.encoder_functions, typeof(x))
        e.encoder_functions[typeof(x)](e, x, askey)
    else
        encode_value(e, x, askey)
    end
end

function encodes_to_string(e::Encoder, x::Any)
    let t = typeof(x)
        if haskey(e.encoder_functions, t)
            e.encodes_to_string[t](e, x)
        else
            encodes_to_string(e, x)
        end
    end
end

Composite = Union{AbstractArray, Dict}

Simple = Union{AbstractString, Integer, BigInt, AbstractFloat, TSymbol, Symbol, Bool}

function encode_top_level(e::Encoder, x::Simple)
    encode_quoted(e, x)
end

function encode_top_level(e::Encoder, x::Composite)
    encode(e, x)
end

function encode_top_level(e::Encoder, x::Any)
    if encodes_to_string(e, x)
        encode_quoted(e, x)
    else
        encode(e, x, false)
    end
end

function encode_quoted(e::Encoder, x::Any)
    emit_array_start(e.emitter)
    emit(e.emitter, "~#'", true)
    emit_array_sep(e.emitter)
    encode(e, x, false)
    emit_array_end(e.emitter)
end

function encode_value(e::Encoder, s::AbstractString, askey::Bool)
    if startswith(s, "~") || startswith(s, "^")
        emit(e.emitter, "~$s", askey)
    else
        emit(e.emitter, s, askey)
    end
end

function encode_value(e::Encoder, s::Symbol, askey::Bool)
    emit(e.emitter, "~:$s", true)
end

function encode_value(e::Encoder, ts::TSymbol, askey::Bool)
    emit(e.emitter, "~\$$(ts.s)", true)
end

function encode_value(e::Encoder, b::Bool, askey::Bool)
    if askey
        emit(e.emitter, b ? "~?t" : "~?f", askey)
    else
        emit(e.emitter, b)
    end
end

function encode_value(e::Encoder, x::Char, askey::Bool)
    emit(e.emitter, "~c$x", askey)
end

function encode_value(e::Encoder, u::URI, askey::Bool)
    s = string(u)
    emit(e.emitter, "~r$s", askey)
end

function encode_value(e::Encoder, u::TURI, askey::Bool)
    s = u.value
    emit(e.emitter, "~r$s", askey)
end

function encode_value(e::Encoder, u::Base.Random.UUID, askey::Bool)
    s = string(u)
    emit(e.emitter, "~u$s", askey)
end

function encode_value(e::Encoder, x::Void, askey::Bool)
    emit_null(e.emitter, askey)
end

function encode_value(e::Encoder, i::Integer, askey::Bool)
    if askey
        emit(e.emitter, "~i$i", askey)
    elseif (i <= JSON_MAX_INT && i >= JSON_MIN_INT)
        emit(e.emitter, i)
    elseif i > MAX_INT64 || i < MIN_INT64
        emit(e.emitter, "~n$i", askey)
    else
        emit(e.emitter, "~i$i", askey)
    end
end

function encode_value(e::Encoder, x::BigInt, askey::Bool)
    let s = string(x)
      emit(e.emitter, "~n$s", askey)
    end
end

function encode_special_float(emitter::Emitter, x::AbstractFloat, askey::Bool)
    if isnan(x)
        emit(emitter, "~zNaN", askey)
    elseif x == Inf
        emit(emitter, "~zINF", askey)
    elseif x == -Inf
        emit(emitter, "~z-INF", askey)
    else
        return false
    end
    return true
end

function encode_value(e::Encoder, x::Decimal, askey::Bool)
    let s = string(x)
        emit(e.emitter, "~f$s", askey)
    end
end

function encode_value(e::Encoder, x::BigFloat, askey::Bool)
    if !encode_special_float(e.emitter, x, askey)
        let s = string(x)
            emit(e.emitter, "~f$s", askey)
        end
    end
end

function encode_value(e::Encoder, x::AbstractFloat, askey::Bool)
    if !encode_special_float(e.emitter, x, askey)
        if askey
             emit(e.emitter, "~d$x", askey)
         else
             emit_raw(e.emitter, string(x))
         end
    end
end

function encode_iterator(e::Encoder, iter)
    for (i, x) in iter
        emit_array_sep(e.emitter, i)
        encode(e, x, false)
    end
end

function encode_tagged_enumerable(e::Encoder, tag::AbstractString, iter)
    emit_array_start(e.emitter)
    emit_tag(e.emitter, tag)
    emit_array_sep(e.emitter)

    emit_array_start(e.emitter)
    encode_iterator(e, iter)
    emit_array_end(e.emitter)
    emit_array_end(e.emitter)
end

function encode_value(e::Encoder, a::AbstractArray, askey::Bool)
    emit_array_start(e.emitter)
    encode_iterator(e, enumerate(a))
    emit_array_end(e.emitter)
end

function encodes_to_string(e::Encoder, x::AbstractArray)
    false
end

function encode_value(e::Encoder, r::Rational, askey::Bool)
    encode_tagged_enumerable(e, "#ratio", enumerate([num(r), den(r)]))
end

function encodes_to_string(e::Encoder, x::Rational)
    false
end

function encode_value(e::Encoder, x::DateTime, askey::Bool)
    let millis = trunc(Int64, Dates.datetime2unix(x) * 1000)
        emit(e.emitter, "~m$millis", askey)
    end
end

function encode_value(e::Encoder, x::Array{Int8}, askey::Bool)
    encoded = base64encode(x)
    emit(e.emitter, "~b$encoded", askey)
end

function encode_value(e::Encoder, x::Date, askey::Bool)
    encode_value(e, DateTime(x), askey)
end

function encode_value(e::Encoder, x::Tuple, askey::Bool)
    encode_tagged_enumerable(e, "#list", enumerate(x))
end

function encode_value(e::Encoder, x::Cons, askey::Bool)
    encode_tagged_enumerable(e, "#list", enumerate(x))
end

function encodes_to_string(e::Encoder, x::Cons)
    false
end

function encode_value(e::Encoder, x::Nil, askey::Bool)
    encode_tagged_enumerable(e, "#list", [])
end

function encodes_to_string(e::Encoder, x::Nil)
    false
end

function encodes_to_string(e::Encoder, x::Tuple)
    false
end

# Maybe can use AbstractSet but instead but need to figure out version boundary
# for when it's available
function encode_value(e::Encoder, x::Set, askey::Bool)
    encode_tagged_enumerable(e, "#set", enumerate(x))
end

function encode_value(e::Encoder, x::TSet, askey::Bool)
    encode_tagged_enumerable(e, "#set", enumerate(x))
end

function encodes_to_string(e::Encoder, x::Set)
    false
end

function encodes_to_string(e::Encoder, x::TSet)
    false
end

function encode_value(e::Encoder, x::Link, askey::Bool)
    let a = [x.href, x.rel, x.name, x.prompt, x.render]
        encode_tagged_enumerable(e, "#link", enumerate(a))
    end
end

function encodes_to_string(e::Encoder, x::Link)
    false
end

function encode_value(e::Encoder, x::TaggedValue, askey::Bool)
    emit_array_start(e.emitter)
    emit_tag(e.emitter, "#$(x.tag)")
    emit_array_sep(e.emitter)
    encode(e, x.value, false)
    emit_array_end(e.emitter)
end

function encodes_to_string(e::Encoder, x::TaggedValue)
    false
end

function has_stringable_keys(e::Encoder, x::Dict)
    for k in keys(x)
	if ! encodes_to_string(e, k)
            return false
        end
    end
    true
end

function encode_map(e::Encoder, x::Dict)
    emit_array_start(e.emitter)
    emit(e.emitter, "^ ", true)

    for (k, v) in x
        emit_array_sep(e.emitter)
        encode_value(e, k, true)
        emit_array_sep(e.emitter)
        encode_value(e, v, false)
    end

    emit_array_end(e.emitter)
end

function encode_cmap(e::Encoder, x::Dict)
    emit_array_start(e.emitter)
    emit(e.emitter, "~#cmap", true)
    emit_array_sep(e.emitter)

    emit_array_start(e.emitter)
    i = 1
    for (k, v) in x
        emit_array_sep(e.emitter, i)
	i = i + 1
        encode_value(e, k, false)
        emit_array_sep(e.emitter)
        encode_value(e, v, false)
    end
    emit_array_end(e.emitter)
    emit_array_end(e.emitter)
end

function encode_verbose_map(e::Encoder, x::Dict)
    emit_map_start(e.emitter)
    i = 1
    for (k, v) in x
        emit_map_sep(e.emitter, i)
        encode_value(e, k, true)
        emit_key_sep(e.emitter)
        encode_value(e, v, false)
	i = i + i
    end
    emit_map_end(e.emitter)
end

function encode_value(e::Encoder, x::Dict, askey::Bool)
    if !has_stringable_keys(e, x)
        encode_cmap(e, x)
    elseif e.verbose
        encode_verbose_map(e, x)
    else
        encode_map(e, x)
    end
end

function encode_value(e::Encoder, x::Dict{AbstractString}, askey::Bool)
    if e.verbose
        encode_verbose_map(e, x)
    else
        encode_map(e, x)
    end
end

function encode_value(e::Encoder, x::Dict{Symbol}, askey::Bool)
    if e.verbose
        encode_verbose_map(e, x)
    else
        encode_map(e, x)
    end
end

function encodes_to_string(e::Encoder, x::Dict)
    false
end

# Default encoder raises exception.
function encode_value(e::Encoder, x::Any, askey::Bool)
    throw(ArgumentError("Don't know how to encode: $x of type $(typeof(x))."))
end

# Default is to claim we do encode to string.
function encodes_to_string(e::Encoder, x::Any)
    true
end


text:
#REC Market Clearing Module

"""
This function does nothing if the product is not REC
"""
function update_rec_supply_curve!(rec_supply_curve::Vector{Vector{Union{String, Float64}}},
    product::T,
    project::P) where {T <: Product, P <: Project{<: BuildPhase}}

    return rec_supply_curve
end

# This function pushes project bids in the REC supply curve
function update_rec_supply_curve!(rec_supply_curve::Vector{Vector{Union{String, Float64}}},
    product::REC,
    project::P) where P <: Project{<: BuildPhase}

    #Each Element of the REC supply Curve:
    #[1] - Project Name (Symbol)
    #[2] - Project REC output (Float64)
    #[3] - Project REC Bid (Float64)

    push!(rec_supply_curve, [get_name(project),
    get_project_rec_output(project),
    get_project_rec_market_bid(project)])

    return rec_supply_curve
end

"""
This function models the actual REC market clearing.
"""
function rec_market_clearing(rec_requirement::Float64,
    ACP::Float64,
    supply_curve::Vector{Vector{Union{String, Float64}}},
    solver::JuMP.MOI.OptimizerWithAttributes)

    #Each Element of the REC supply Curve:
    #[1] - Project Name (Symbol)
    #[2] - Project REC output (Float64)
    #[3] - Project REC Bid (Float64)

    n_supply_seg = length(supply_curve);

    #----------REC Market Clearing Problem---------------------------------------------------------------------------------
    rec_mkt = JuMP.Model(solver);

    #Define the variables
    JuMP.@variables(rec_mkt, begin
    Q_supply[s=1:n_supply_seg] >= 0 # Quantity of cleared REC supply offers
    v_REC >= 0                      # RPS compliance violation
    end)

    #Functions----------------------------------------------------------------------------------------

    #Cost of procuring REC supply for each segment
    supply_cost(Q_seg, s) = supply_curve[s][3] * Q_seg;

    #Expressions---------------------------------------------------------------------------------------

    #Totoal Cost of procuring REC supply
    JuMP.@expression(rec_mkt, total_supply_cost, sum(supply_cost(Q_supply[s], s) for s = 1:n_supply_seg))

    #Constraints--------------------------------------------------------------------------------------

    #Cleared REC supply limit for each segment
    JuMP.@constraint(rec_mkt, [s=1:n_supply_seg], Q_supply[s] <= supply_curve[s][2])

    #Total cleared REC supply should be greater than or equal to the RPS compliance requirement
    JuMP.@constraint(rec_mkt, mkt_clear, sum(Q_supply[s] for s in 1:n_supply_seg) + v_REC >= rec_requirement)

    #Define Objective Function - Social Welfare Maxmization
    JuMP.@objective(rec_mkt, Min, total_supply_cost + v_REC * ACP)

    println("Actual REC Market Clearing:")
    JuMP.optimize!(rec_mkt)

    println(JuMP.termination_status(rec_mkt))
    println(JuMP.objective_value(rec_mkt))

    #REC Market Clearing Price is the shadow variable of the REC balance constraint
    rec_price = AxisArrays.AxisArray(reshape([JuMP.dual(mkt_clear)], 1,), [1])
    rec_accepted_bid = Dict(supply_curve[s][1] => value.(Q_supply[s]) for s in 1:n_supply_seg)
    #------------------------------------------------------------------------------------------------
    return rec_price, rec_accepted_bid

end


text:
macro monkey_eval_str(code::String)
    quote
        run($(esc(Meta.parse("\"$(escape_string(code))\""))))
    end
end

evaluate(code::String; input = stdin, output = stdout) = begin
    raw_program = parse(code; input, output)
    if !isnothing(raw_program)
        macro_env = Environment(; input, output)
        program = define_macros!(macro_env, raw_program)
        expanded = expand_macros(program, macro_env)

        syntax_check_result = analyze(expanded)
        if isa(syntax_check_result, ErrorObj)
            println(output, syntax_check_result)
            return syntax_check_result
        end

        result = evaluate(expanded, Environment(; input, output))
        if isa(result, ErrorObj)
            println(output, result)
        end

        return result
    end
end

evaluate(::Node, ::Environment) = _NULL
evaluate(node::ExpressionStatement, env::Environment) = evaluate(node.expression, env)
evaluate(node::IntegerLiteral, ::Environment) = IntegerObj(node.value)
evaluate(node::BooleanLiteral, ::Environment) = node.value ? _TRUE : _FALSE
evaluate(node::StringLiteral, ::Environment) = StringObj(node.value)
evaluate(node::FunctionLiteral, env::Environment) =
    FunctionObj(node.parameters, node.body, env)

evaluate(node::ArrayLiteral, env::Environment) = begin
    elements = evaluate_expressions(node.elements, env)
    if length(elements) == 1 && isa(elements[1], ErrorObj)
        return elements[1]
    end
    return ArrayObj(elements)
end

evaluate(node::HashLiteral, env::Environment) = begin
    pairs = Dict{Object,Object}()

    for (key_node, value_node) in node.pairs
        key = evaluate(key_node, env)
        if isa(key, ErrorObj)
            return key
        end

        value = evaluate(value_node, env)
        if isa(value, ErrorObj)
            return value
        end

        pairs[key] = value
    end

    return HashObj(pairs)
end

evaluate(node::Identifier, env::Environment) = begin
    val = get(env, node.value)

    if !isnothing(val)
        return val
    end

    builtin = get_builtin_by_name(node.value)

    if !isnothing(builtin)
        return builtin
    end
end

evaluate(node::PrefixExpression, env::Environment) = begin
    right = evaluate(node.right, env)
    return isa(right, ErrorObj) ? right : evaluate_prefix_expression(node.operator, right)
end

evaluate(node::InfixExpression, env::Environment) = begin
    left = evaluate(node.left, env)
    if isa(left, ErrorObj)
        return left
    end
    right = evaluate(node.right, env)
    return isa(right, ErrorObj) ? right :
           evaluate_infix_expression(node.operator, left, right)
end

evaluate(node::IfExpression, env::Environment) = begin
    condition = evaluate(node.condition, env)

    if isa(condition, ErrorObj)
        return condition
    elseif is_truthy(condition)
        return evaluate(node.consequence, env)
    elseif !isnothing(node.alternative)
        return evaluate(node.alternative, env)
    else
        return _NULL
    end
end

evaluate(node::CallExpression, env::Environment) = begin
    # TODO: Currently, `quote()` only processes its first argument
    if token_literal(node.fn) == "quote"
        if isempty(node.arguments)
            return _NULL
        end
        return QuoteObj(evaluate_unquote_calls(node.arguments[1], env))
    end

    fn = evaluate(node.fn, env)
    if isa(fn, ErrorObj)
        return fn
    end

    args = evaluate_expressions(node.arguments, env)
    if length(args) == 1 && isa(args[1], ErrorObj)
        return args[1]
    end

    return isa(fn, Builtin) ? apply_builtin(fn, args, env) : apply_function(fn, args)
end

evaluate(node::IndexExpression, env::Environment) = begin
    left = evaluate(node.left, env)
    if isa(left, ErrorObj)
        return left
    end

    index = evaluate(node.index, env)
    if isa(index, ErrorObj)
        return index
    end

    return evaluate_index_expression(left, index)
end

evaluate(node::LetStatement, env::Environment) = begin
    val = evaluate(node.value, env)
    if isa(val, ErrorObj)
        return val
    end
    if node.reassign
        return reassign!(env, node.name.value, val)
    else
        set!(env, node.name.value, val)
    end
    return val
end

evaluate(node::ReturnStatement, env::Environment) = begin
    val = evaluate(node.return_value, env)
    return isa(val, ErrorObj) ? val : ReturnValue(val)
end

evaluate(::BreakStatement, env::Environment) = BreakObj()

evaluate(::ContinueStatement, env::Environment) = ContinueObj()

evaluate(node::WhileStatement, env::Environment) = begin
    while true
        condition = evaluate(node.condition, env)
        if isa(condition, ErrorObj)
            return condition
        end

        if is_truthy(condition)
            result = evaluate(node.body, Environment(env))
            if isa(result, ReturnValue) || isa(result, ErrorObj)
                return result
            elseif isa(result, BreakObj)
                break
            elseif isa(result, ContinueObj)
                continue
            end
        else
            break
        end
    end

    return _NULL
end

evaluate(block::BlockStatement, env::Environment) = begin
    result = _NULL

    for statement in block.statements
        result = evaluate(statement, env)
        if isa(result, BreakObj) ||
           isa(result, ContinueObj) ||
           isa(result, ReturnValue) ||
           isa(result, ErrorObj)
            return result
        end
    end

    return result
end

evaluate(program::Program, env::Environment) = begin
    result = _NULL

    for statement in program.statements
        result = evaluate(statement, env)
        if isa(result, ReturnValue)
            return result.value
        elseif isa(result, ErrorObj)
            return result
        end
    end

    return result
end

function evaluate_prefix_expression(operator::String, right::Object)
    if operator == "!"
        return evaluate_bang_operator_expression(right)
    elseif operator == "-"
        return evaluate_minus_prefix_operator_expression(right)
    else
        return ErrorObj("unknown operator: " * operator * type_of(right))
    end
end

function evaluate_infix_expression(operator::String, left::Object, right::Object)
    if type_of(left) != type_of(right)
        return ErrorObj(
            "type mismatch: " * type_of(left) * " " * operator * " " * type_of(right),
        )
    end

    if operator == "=="
        return left == right ? _TRUE : _FALSE
    elseif operator == "!="
        return left != right ? _TRUE : _FALSE
    else
        return ErrorObj(
            "unknown operator: " * type_of(left) * " " * operator * " " * type_of(right),
        )
    end
end

function evaluate_infix_expression(operator::String, left::StringObj, right::StringObj)
    if operator == "+"
        return StringObj(left.value * right.value)
    elseif operator == "=="
        return left.value == right.value ? _TRUE : _FALSE
    elseif operator == "!="
        return left.value != right.value ? _TRUE : _FALSE
    else
        return ErrorObj(
            "unknown operator: " * type_of(left) * " " * operator * " " * type_of(right),
        )
    end
end

function evaluate_infix_expression(operator::String, left::IntegerObj, right::IntegerObj)
    if operator == "+"
        return IntegerObj(left.value + right.value)
    elseif operator == "-"
        return IntegerObj(left.value - right.value)
    elseif operator == "*"
        return IntegerObj(left.value * right.value)
    elseif operator == "/"
        if right.value == 0
            return ErrorObj("divide error: division by zero")
        end
        return IntegerObj(left.value ÷ right.value)
    elseif operator == "<"
        return left.value < right.value ? _TRUE : _FALSE
    elseif operator == ">"
        return left.value > right.value ? _TRUE : _FALSE
    elseif operator == "=="
        return left.value == right.value ? _TRUE : _FALSE
    elseif operator == "!="
        return left.value != right.value ? _TRUE : _FALSE
    else
        return ErrorObj(
            "unknown operator: " * type_of(left) * " " * operator * " " * type_of(right),
        )
    end
end

function evaluate_bang_operator_expression(right::Object)
    if right === _FALSE || right === _NULL
        return _TRUE
    else
        return _FALSE
    end
end

evaluate_minus_prefix_operator_expression(right::Object) =
    ErrorObj("unknown operator: -" * type_of(right))
evaluate_minus_prefix_operator_expression(right::IntegerObj) = IntegerObj(-right.value)

evaluate_index_expression(left::Object, ::Object) =
    ErrorObj("index operator not supported: $(type_of(left))")
evaluate_index_expression(::ArrayObj, index::Object) =
    ErrorObj("unsupported index type: $(type_of(index))")
evaluate_index_expression(hash::HashObj, key::Object) = Base.get(hash.pairs, key, _NULL)
evaluate_index_expression(left::ArrayObj, index::IntegerObj) = begin
    idx = index.value
    max_idx = length(left.elements) - 1
    return 0 <= idx <= max_idx ? left.elements[idx+1] : _NULL
end

function evaluate_expressions(expressions::Vector{Expression}, env::Environment)
    results = Object[]

    for expression in expressions
        evaluated = evaluate(expression, env)
        if isa(evaluated, ErrorObj)
            return [evaluated]
        end
        push!(results, evaluated)
    end

    return results
end

evaluate_unquote_calls(quoted::Node, env::Environment) = modify(
    quoted,
    function modifier(node::Node)
        if isa(node, CallExpression)
            if token_literal(node.fn) == "unquote" && length(node.arguments) >= 1
                Node(evaluate(node.arguments[1], env))
            else
                return CallExpression(
                    node.token,
                    modify(node.fn, modifier),
                    map(expression -> modify(expression, modifier), node.arguments),
                )
            end
        else
            return node
        end
    end,
)

function apply_function(fn::FunctionObj, args::Vector{Object})
    if length(fn.parameters) != length(args)
        return ErrorObj("argument error: wrong number of arguments: got $(length(args))")
    end

    extended_env = extend_function_environment(fn, args)
    evaluated = evaluate(fn.body, extended_env)
    return unwrap_return_value(evaluated)
end

apply_function(fn::Object, ::Vector{Object}) = ErrorObj("not a function: " * type_of(fn))
apply_builtin(fn::Builtin, args::Vector{Object}, env::Environment) =
    fn.fn(args...; env = env)

function extend_function_environment(fn::FunctionObj, args::Vector{Object})
    env = Environment(fn.env)

    for (param, arg) in zip(fn.parameters, args)
        set!(env, param.value, arg)
    end

    return env
end

function unwrap_return_value(obj::Object)
    if isa(obj, ReturnValue)
        return obj.value
    else
        return obj
    end
end

Node(::Object) = NullLiteral(Token(NULL, "null"))
Node(obj::IntegerObj) = IntegerLiteral(Token(INT, string(obj.value)), obj.value)
Node(obj::StringObj) = StringLiteral(Token(STRING, obj.value), obj.value)
Node(obj::BooleanObj) =
    BooleanLiteral(Token(obj.value ? TRUE : FALSE, string(obj.value)), obj.value)
Node(obj::HashObj) = HashLiteral(
    Token(LBRACE, "{"),
    Dict(Node(key) => Node(value) for (key, value) in collect(obj.pairs)),
)
Node(obj::ArrayObj) = ArrayLiteral(Token(LBRACKET, "["), map(Node, obj.elements))
Node(obj::FunctionObj) = FunctionLiteral(Token(FUNCTION, "fn"), obj.parameters, obj.body)
Node(obj::QuoteObj) = obj.node

# TODO: Currently, only top-level macro definitions are allowed. We don’t walk down the Statements and check the child nodes for more.
function define_macros!(env::Environment, program::Program)
    new_statements = []

    for statement in program.statements
        if is_macro_definition(statement)
            add_macro!(env, statement)
        else
            push!(new_statements, statement)
        end
    end

    return Program(new_statements)
end

is_macro_definition(::Node) = false
is_macro_definition(ls::LetStatement) = isa(ls.value, MacroLiteral)

add_macro!(env::Environment, ls::LetStatement) = begin
    name = ls.name.value
    value = ls.value
    set!(env, name, MacroObj(value.parameters, value.body, env))
end

function expand_macros(program::Program, env::Environment)
    return modify(
        program,
        (node::Node) -> begin
            if !isa(node, CallExpression)
                return node
            end

            mc = get_macro_call(node, env)
            if isnothing(mc)
                return node
            end

            args = quote_args(node)
            eval_env = extend_macro_env(mc, args)
            evaluated = evaluate(mc.body, eval_env)

            if !isa(evaluated, QuoteObj)
                error("macro error: we only support returning AST-nodes from macros")
            end

            return evaluated.node
        end,
    )
end

function get_macro_call(node::CallExpression, env::Environment)
    ident = node.fn
    if !isa(ident, Identifier)
        return nothing
    end

    mc = get(env, ident.value)
    if !isa(mc, MacroObj)
        return nothing
    end

    return mc
end

quote_args(node::CallExpression) = map(QuoteObj, node.arguments)

function extend_macro_env(mc::MacroObj, args::Vector{QuoteObj})
    extended_env = Environment(mc.env)
    for (parameter, argument) in zip(mc.parameters, args)
        set!(extended_env, parameter.value, argument)
    end
    return extended_env
end


text:
#-----------------------------------------------------------------------
#       BEGIN   atm.jl
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   STANDARD ATMOSPHERIC PROFILE
#-----------------------------------------------------------------------
"""
    ΔnatmStd( h::Real )

The function `ΔnatmStd` is a profile for the atmospheric index of
of refraction obtained from the U. S. standard atmosphere data. The
argument `h` is height from the Earth's surface in km.

"""
function ΔnatmStd( h::Real ) # h is height from surface in km
    tpfl=typeof(h)
    return  -222.66559677864387/
            (
                99.06212001632495 +
                0.15782316597789148*(7.1540968904797175 + h)^2
            ) +
            253.49880024119224/
            (
                112.75665221238798 +
                0.17966910753975862*(7.156537366807337 + h)^2
            )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   REVISED EDLEN FORMULA
#-----------------------------------------------------------------------
"""
    ΔnREdlen( P::Real , T::Real , λ::Real )

The function `ΔnREdlen` provides the (revised) Edlen formula for
computing the index of refraction from atmospheric parameters `P`, `T`,
and `λ`, corresponding to the respective pressure (in Pascals),
temperature (in Celsius), and optical wavelength (in μm).

"""
function ΔnREdlen( P::Real , T::Real , λ::Real )
    tpfl=typeof(P)

    # dispersion relation
    disp = (tpfl(1e-8))*
           tpfl(8342.54 + 2406147.0/(130-λ^(-2)) + 15998.0/(38.9-λ^(-2)))

    return tpfl(P*disp/96095.43)*
           tpfl(
            ( 1. + (1e-8)*(0.601 - 0.00972*T)*P )
            /
            ( 1. + 0.0036610*T )
           )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#       END     atm.jl
#-----------------------------------------------------------------------


text:
using BacenSGS
using Test

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


text:
"""
@ev_draw msg eq init=false

1) If optional argument init = true, set default layout parameters (see below)
2) Place the text 'msg' on the current drawing at global variable cpt.
3) Latexify expression(s) eq and place it as an svg below the text
4) Evaluate eq in the calling context
5) Update layout parameter cpt - current point, based on layout parameters.

Macro for cutting down on boilerplate code. This is not sufficient for
all cases. When this does not produce expected results,
use @macroexpand to see what occurs, and copy / modify that code.

Also, modifying line spacing (Δcpy) and (scalelatex) and of course
the current point (CP) may increase the
usefulness of this macro.

Note that 'msg' accepts pango-style html-like markup for text formatting.

Additional input: locally defined variables
cpt     Point   Mutating current position, bottom left of displayed first line
Δcpy    Float64  Line spacing after text lines, typically 2 (points)
scalelatex = 3.143

## Example:
```julia
cpt = O + (-WI / 2 + EM, -HE / 2 + 2EM)
Δcpy = 2
scalelatex = 3.143
pt += @ev_draw  "Unit stripping function:" :(ul(x) = x / oneunit(x))
```
"""
macro ev_draw(msg, eq, init = false)
    # Check macro input in this module's context
    @assert typeof(eq) <: Expr   "Second argument not an expression, at \n\t$__source__"
    expr = draw_init(init)
    assertions!(expr, __source__)
    place_msg!(expr, msg)
    place_eqs!(expr, eq)

    # Evaluate the equation / expression / block in the calling context - main scope.
    push!(expr.args, esc(quote
        eval($eq)
    end))

    # If the current drawing is a recording surface, return a snapshot
    push!(expr.args, quote
        if currentdrawing().surfacetype == :rec
            snapshot()
        end
    end)
    expr
end

"""
    eq_draw(msg, eq, init = false)
Skips the evaluation, otherwise same as @ev_draw
"""
macro eq_draw(msg, eq, init = false)
    # Check macro input in this module's context
    @assert typeof(eq) <: Expr   "Second argument not an expression, at \n\t$__source__"
    expr = draw_init(init)
    assertions!(expr, __source__)
    place_msg!(expr, msg)
    place_eqs!(expr, eq)
    # If the current drawing is a recording surface, return a snapshot
    push!(expr.args, quote
        if currentdrawing().surfacetype == :rec
            snapshot()
        end
    end)
    expr
end

"""
    draw_init(init)

Create an expression skeleton for macros `@ev_draw`.
init is an expression which is evaluated in the calling context.
If 'init = true', layout parameters are defined globally
in the calling context:
    cpt Current point
    Δcpy Verttical spacing
    scalelatex

For better control, define these prior to the first call.
"""
function draw_init(init)
    expr = Expr(:block)
    push!(expr.args, esc(quote
        if eval($init) == true
            cpt = O + (-WI / 2 + EM, -HE / 2 + 2EM)
            Δcpy = 2
            scalelatex = 3.143
        end
    end))
    expr
end

"""
    assertions!(expr, __source__)

Check that global variables and functions exist in the calling context, for
macro `@ev_draw`.
"""
function assertions!(expr, __source__)

    booleancheck = Expr(:escape, Expr(:isdefined, :cpt))
    assertion = Expr(:(macrocall), Symbol("@assert"), __source__, booleancheck, "Global 'cpt' (current Point) not defined, at \n\t$__source__")
    push!(expr.args, assertion)

    booleancheck = Expr(:escape, Expr(:isdefined, :Δcpy))
    assertion = Expr(:(macrocall), Symbol("@assert"), __source__, booleancheck, "Global 'Δcpy' (additional line spacing, a unitless number) not defined, at \n\t$__source__")
    push!(expr.args, assertion)

    booleancheck = Expr(:escape, Expr(:isdefined, :scalelatex))
    assertion = Expr(:(macrocall), Symbol("@assert"), __source__, booleancheck, "Global 'scalelatex' (e.g. 3.142, normally based on output from place_img with defined height) not defined, at \n\t$__source__")
    push!(expr.args, assertion)

#    booleancheck = Expr(:escape, Expr(:isdefined, :currentdrawing))
#    assertion = Expr(:(macrocall), Symbol("@assert"), __source__, booleancheck, "Global 'currentdrawing()' not defined, at \n\t$__source__")
#    push!(expr.args, assertion)
    expr
end

"""
    place_msg!(expr, msg)

Place the text (no latex) on drawing, update layout parameters in the calling context
of macro `@ev_draw`.
"""
function place_msg!(expr, msg)
    # Evaluate string(msg) in the calling context to get a string for sure.
    push!(expr.args, esc(quote
        local st = eval(string($msg))
    end))

    # Unfortunately, the 'Text API' does not provide easily interpretable feedback on
    # how much vertical space 'msg' occupies. We rely on constant EM and line counting.
    push!(expr.args, esc(quote
        local textlinebreaks = count(c -> c == '\n' || c == '\r', collect(st))
    end))

    # 'Settext' will displace multi-line output upwards. But 'cpt' keeps track of
    # the next line's baseline. So we move 'cpt' down by 'textlinebreaks' * (EM + Δcpy) .
    push!(expr.args, esc(quote
        local Δy = ($EM + Δcpy) * textlinebreaks
        local Δpt =  Point(0, Δy)
    end))
    push!(expr.args, esc(Expr(:+=, :cpt, :Δpt)))

    # Draw string in the modified position
    push!(expr.args, esc(quote
        settext(st, cpt, markup = true)
    end))
    # Move current point to next base line, as if we were drawing more text.
    push!(expr.args, esc(quote
        Δy = ($EM + Δcpy)
        Δpt =  Point(0, Δy)
    end))
    push!(expr.args, esc(Expr(:+=, :cpt, :Δpt)))
    expr
end

"""
    place_eqs!(expr, eq)

Draw the equation(s). The reference point here is not baseline of text,
but the upper left corner of the svg image. Compensate for that.
Uupdate layout parameters in the calling context of macro `@ev_draw`.
"""
function place_eqs!(expr, eq)
    # Draw the equation(s). The reference point here is not baseline of text,
    # but the upper left corner of the svg image. Compensate for that.
    push!(expr.args, esc(quote
        local ptul, ptbr = draw_expr(cpt + Point(0, -EM), scalelatex, $eq)
        # circle(cpt; r = 40mm) # For debugging of layout parameters.
        # circle(ptul; r = 50mm)
        # circle(ptbr; r = 50mm)
    end))

    # Move to next line. Assume the next line will be text, so add the 'EM' that we skipped before the equation
    push!(expr.args, esc(quote
        Δy = ptbr.y - ptul.y
        Δpt =  Point(0, Δy)
    end))
    push!(expr.args, esc(Expr(:+=, :cpt, :Δpt)))
    expr
end


text:
using ZChop # For chopping small imaginary parts in ρ
using Distributed
using TimerOutputs
using ProgressMeter

function simulate_trajectory(model::Tm,
                             initial_state::Ts,
                             file_channel::Union{Channel,RemoteChannel,Nothing}=nothing,
                             progress_channel::Union{Channel,RemoteChannel,Nothing}=nothing;
                             callbacks=[]) where {Tm <: Model, Ts <:State}

    state = Ts(initial_state) # Copy the initial state

    jx = Array{Float64}(undef, length(get_time(model)))

    # Output variables
    jx = similar(jx)
    jy = similar(jx)
    jz = similar(jx)

    jx2 = similar(jx)
    jy2 = similar(jx)
    jz2 = similar(jx)

    FisherT = similar(jx)
    QFisherT = similar(jx)

    FisherStrong = similar(jx)

    jto = 1 # Counter for the output

    for jt = 1 : model.params.Ntime
        dy = measure_current(state, model)
        tr_ρ, tr_τ = updatestate!(state, model, dy)

        if jt % model.params._outsteps == 0

            jx[jto] = expectation_value!(state, model.Jx)
            jy[jto] = expectation_value!(state, model.Jy)
            jz[jto] = expectation_value!(state, model.Jz)

            jx2[jto] = expectation_value!(state, model.Jx2)
            jy2[jto] = expectation_value!(state, model.Jy2)
            jz2[jto] = expectation_value!(state, model.Jz2)

            # We evaluate the classical FI for the continuous measurement
            FisherT[jto] = real(tr_τ^2)

            # We evaluate the QFI for a final strong measurement done at time t
            QFisherT[jto] = QFI!(state)

            # Fisher of the strong final measurement
            FisherStrong[jto] = FI(state, model.measurement)

            # apply any callbacks
            isnothing(callbacks) || for cb in callbacks
                cb(state, model)
            end

            jto += 1
            isnothing(progress_channel) || put!(progress_channel, true)
        end
    end

    Δjx2 = jx2 - jx.^2
    Δjy2 = jy2 - jy.^2
    Δjz2 = jz2 - jz.^2

    xi2x = squeezing_param(model.params.Nj, Δjx2, jy, jz)
    xi2y = squeezing_param(model.params.Nj, Δjy2, jx, jz)
    xi2z = squeezing_param(model.params.Nj, Δjz2, jx, jy)

    result = (FI=FisherT, QFI=QFisherT,
                    Jx=jx, Jy=jy, Jz=jz,
                    Δjx2=Δjx2, Δjy2=Δjy2, Δjz2=Δjz2,
                    xi2x=xi2x, xi2y=xi2y, xi2z=xi2z, FIstrong=FisherStrong)
    isnothing(file_channel) || put!(file_channel, result)

    GC.gc()
    return result
end

function Eff_QFI_HD_Dicke(Nj::Int64, # Number of spins
    Ntraj::Int64,                    # Number of trajectories
    Tfinal::Real,                    # Final time
    dt::Real;                        # Time step
    κ::Real = 1.,                    # Independent noise strength
    κcoll::Real = 1.,                # Collective noise strength
    ω::Real = 0.0,                   # Frequency of the Hamiltonian
    η::Real = 1.,                    # Measurement efficiency
    outpoints = 0,                   # Number of output points
    to = TimerOutput(), file_channel=nothing)

    @info "Eff_QFI_HD_Dicke starting"
    @info "Parameters" Nj Ntraj Tfinal dt κ κcoll ω η outpoints

    modelparams = ModelParameters(Nj=Nj, kind=κ, Gamma=κcoll, kcoll=0., omega=ω, eta=η, dt=dt, Tfinal=Tfinal, outpoints=outpoints)

    model = LocalDephasingModel(modelparams)
    initial_state = blockdiag_css(Nj)

    # Run evolution for each trajectory, and build up the average
    # for FI and final strong measurement QFI
    @timeit_debug to "trajectories" begin
    #result = @showprogress 1 "Computing..." @distributed (+) for ktraj = 1 : Ntraj
    result = @distributed (+) for ktraj = 1 : Ntraj
        result = simulate_trajectory(model, initial_state)
        hcat(result.FI, result.QFI,
             result.Jx, result.Jy, result.Jz,
             result.Δjx2, result.Δjy2, result.Δjz2,
             result.xi2x, result.xi2y, result.xi2z, result.FIstrong)
    end
    end

    # Evaluate averages
    jx = result[:, 3] / Ntraj
    jy = result[:, 4] / Ntraj
    jz = result[:, 5] / Ntraj

    Δjx2 = result[:, 6] / Ntraj
    Δjy2 = result[:, 7] / Ntraj
    Δjz2 = result[:, 8] / Ntraj

    xi2x = result[:, 9] / Ntraj
    xi2y = result[:, 10] / Ntraj
    xi2z = result[:, 11] / Ntraj

    @info "Time details \n$to"
    return (t=get_time(model),
            FI=result[:,1] / Ntraj,
            QFI=result[:,2] / Ntraj,
            jx=jx, jy=jy, jz=jz,
            Δjx=Δjx2, Δjy=Δjy2, Δjz=Δjz2,
            xi2x=xi2x, xi2y=xi2y, xi2z=xi2z, FIstrong=result[:, 12] / Ntraj)
end


text:
googlewmt_init() = register(DataDep(
    "Google-WMT en-de",
    """
    \"\"\"shows in wmt14 of torchtext
    The WMT 2014 English-German dataset, as preprocessed by Google Brain.

    Though this download contains test sets from 2015 and 2016, the train set
    differs slightly from WMT 2015 and 2016 and significantly from WMT 2017.
    \"\"\"

    contain bpe training set and news testset from 2009~2016 (include origin text,
    tokenized, and bpe versions), and also a bpe.32000 and vocab.32000 (merged vocab)
    """,
    "https://drive.google.com/uc?export=download&id=0B_bZck-ksdkpM25jRUN2X2UxMm8",
    "86f7f6e0bab3259a34712b4c034dc933f406d6735ce03fd6de6b9ccb5191ce2e";
    fetch_method=download_gdrive,
    post_fetch_method=DataDeps.unpack,
))

struct GoogleWMT <: Dataset end

function testfile(::GoogleWMT, year=2014; mode=:bpe)
    !(2009 <= year <= 2016) && error("year shoud be in 2009~2016")
    if mode == :bpe
        en = @datadep_str "Google-WMT en-de/newstest$(year).tok.bpe.32000.en"
        de = @datadep_str "Google-WMT en-de/newstest$(year).tok.bpe.32000.de"
    elseif mode == :tok
        en = @datadep_str "Google-WMT en-de/newstest$(year).tok.en"
        de = @datadep_str "Google-WMT en-de/newstest$(year).tok.de"

    elseif mode == :plain
        en = @datadep_str "Google-WMT en-de/newstest$(year).en"
        de = @datadep_str "Google-WMT en-de/newstest$(year).de"
    else
        error("model should be one of :bpe, :tok, and :plain")
    end
    en, de
end

function trainfile(::GoogleWMT)
    en = datadep"Google-WMT en-de/train.tok.clean.bpe.32000.en"
    de = datadep"Google-WMT en-de/train.tok.clean.bpe.32000.de"
    en, de
end

function get_vocab(::GoogleWMT; mode = :vocab)
    if mode == :vocab
        vf = datadep"Google-WMT en-de/vocab.bpe.32000"
        voc = open(vf) do f
            d = Dict{String, Int}()
            for (i, l) ∈ enumerate(eachline(f))
                d[intern(l)] = i
            end
            d
        end
        return voc
    elseif mode == :bpe
        vf = datadep"Google-WMT en-de/bpe.32000"
        bpe = Bpe(vf; sepsym="@@", have_header=false)
        return bpe
    else
        error("mode should be one of :vocab and :bpe")
    end
end


text:
using FractalDevelopment
using Documenter

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

makedocs(;
    modules=[FractalDevelopment],
    authors="Harry Booth <39346323+harrybooth@users.noreply.github.com> and contributors",
    repo="https://github.com/harrybooth/FractalDevelopment.jl/blob/{commit}{path}#{line}",
    sitename="FractalDevelopment.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://harrybooth.github.io/FractalDevelopment.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/harrybooth/FractalDevelopment.jl",
)

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
using MIDI
using BSON
using Flux
using CUDA
using MusicTransformer
using Coconet

tpq = 960
sixteenth = tpq / 4

BSON.@load "checkpoint/my_model.bson" model

dataset = BSON.load("datasets/Tonnetz_Bach.bson")

model = gpu(model)

test_song = gpu(dataset[:valid][rand(1:length(dataset[:valid]))])
song_mask = Coconet.mask(test_song, 0.9) |> gpu
masked_song = test_song ./ (song_mask[1,1,:,:,:] .+ 1)
pos_data = gpu(get_positional_encoding(test_song))

prediction = model(masked_song, pos_data, mask = song_mask)

loss = test_song .- prediction
loss = sum(loss .* loss)

song = zeros(size(prediction)[2:end])

# for each instrument
for i in 1:size(prediction, 3)
    for j in 1:size(prediction, 2)
        if round(prediction[3,j,i]) > 0
            pitch = tonnetz_to_pitch(prediction[1,j,i], prediction[2,j,i])
            song[j,i] = pitch
        end
    end
end

song

file = MIDIFile(1, tpq, output_tracks)
writeMIDIFile(string("results/output_file_", 1), file)

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet


text:
"""
The type of which all objects in the class are a subtype
"""
abstract type ABasic end

### Terms ###

struct Term{T<:ABasic}
    data::Symbol
    class::Type{T}
    args::NTuple{N,Term} where N
end

# constructors
# Term(data::Symbol, class::Type{<:ABasic}, args::NTuple{N,Term} where N) = Term(data, class, args, Ref{UInt}(0))
Term(data::Symbol, class::Type{<:ABasic}) = Term(data, class, ())

# helper functions for Term
data(self::Term) = self.data
class(self::Term) = self.class
func(self::Term) = func(class(self))
domain(x::Term) = domain(typeof(x))
codomain(x::Term) = codomain(typeof(x))
domain(::Type{Term}) = ∅
codomain(::Type{Term}) = ∅
args(self::Term) = self.args
export domain, codomain, arity

### Symbols ###
"""
function symbols and constant symbols
    """
abstract type AFixed <: ABasic end
"""
any "known" function with arity of one or more
    """
abstract type AFunctionSymbol <: AFixed end
"""
roots of trees can only ever be a subtype of this
can also be called singletons or atoms
"""
abstract type AConstantSymbol <: AFixed end
domain(::Type{<:Term{<:AConstantSymbol}}) = ∅

### Variables ###
"""
unknowns like `x` and `f`
all unknwons have domains (values that they can take in) and
a codomain (values that it can take on)
"""
abstract type AVariable <: ABasic end
struct Variable <: AVariable end

Variable(name::Symbol, codomain::Term) = Term(name, Variable, (codomain,))
Variable(name::String, codomain::Term) = Variable(Symbol(name), codomain)
Base.show(io::IO, self::Term{<:AVariable}) = print(io, data(self))
domain(::Term{<:AVariable}) = ∅
codomain(self::Term{<:AVariable}) = args(self)[1]
export Variable

struct FuncVariable <: AVariable end
FuncVariable(name::Symbol, domain::Term, codomain::Term) = Term(name, FuncVariable, (domain, codomain))
FuncVariable(name::String, domain::Term, codomain::Term) = FuncVariable(Symbol(name), domain, codomain)
(f::Term{FuncVariable})(x...) = Term(data(f), FuncVariable, x)
domain(f::Term{<:FuncVariable}) = args(f)[1]
codomain(f::Term{<:FuncVariable}) = args(f)[2]
export FuncVariable

### Macros for new constants ###
macro new_constant(parent::Symbol, abstract_type::Symbol, name::Symbol)
    x = esc(name)
    return quote
        abstract type $(esc(abstract_type)) <: $(esc(parent)) end
        const $(esc(name)) = Term($(Meta.quot(name)), $(esc(abstract_type)))
        Base.show(io::IO, x::Term{<:$(esc(abstract_type))}) = print(io, $(Meta.quot(name)))
        export $(esc(name))
    end
end

macro new_constant(parent::Symbol, abstract_type::Symbol, name::Symbol, codomain::Symbol)
    x = esc(name)
    return quote
        abstract type $(esc(abstract_type)) <: $(esc(parent)) end
        const $(esc(name)) = Term($(Meta.quot(name)), $(esc(abstract_type)))
        Base.show(io::IO, x::Term{<:$(esc(abstract_type))}) = print(io, $(Meta.quot(name)))
        $(esc(:codomain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(codomain))
        export $(esc(name))
    end
end

macro new_function(parent::Symbol, abstract_type::Symbol, name::Symbol, domain, codomain, arity)
    return quote
        abstract type $(esc(abstract_type)) <: $(esc(parent)) end
        struct $(esc(name)) <: $(esc(abstract_type)) end
        const $(esc(name))(x::Vararg{Term}) = Term($(Meta.quot(name)), $(esc(name)), x)
        $(esc(:domain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(domain))
        $(esc(:codomain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(codomain))
        $(esc(:arity))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(arity))

        export $(esc(name))
    end
end

macro new_function(parent::Symbol, abstract_type::Symbol, name::Symbol, domain, codomain, arity, latex_name, should_export::Bool=true)
    if should_export
        return quote
            abstract type $(esc(abstract_type)) <: $(esc(parent)) end
            struct $(esc(name)) <: $(esc(abstract_type)) end
            const $(esc(name))(x::Vararg{Term}) = Term($(Meta.quot(latex_name)), $(esc(name)), x)
            $(esc(:domain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(domain))
            $(esc(:codomain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(codomain))
            $(esc(:arity))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(arity))

            $(esc(latex_name))(x::Vararg{Term}) = $(esc(name))(x...)

            export $(esc(name)), $(esc(latex_name))
        end
    else
        return quote
            abstract type $(esc(abstract_type)) <: $(esc(parent)) end
            struct $(esc(name)) <: $(esc(abstract_type)) end
            const $(esc(name))(x::Vararg{Term}) = Term($(Meta.quot(latex_name)), $(esc(name)), x)
            $(esc(:domain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(domain))
            $(esc(:codomain))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(codomain))
            $(esc(:arity))(::Type{<:Term{<:$(esc(abstract_type))}}) = $(esc(arity))

            $(esc(latex_name))(x::Vararg{Term}) = $(esc(name))(x...)

            export $(esc(name))
        end
    end
end


text:
# # Estimate Binomial draw probabilility

using DynamicHMCModels, MCMCChains

Random.seed!(1356779)

# Define a structure to hold the data.

Base.@kwdef struct BernoulliProblem
    "Total number of draws in the data."
    n::Int
    "Number of draws ' == 1' "
    obs::Vector{Int}
end;

# Write a function to return properly dimensioned transformation.

make_transformation(model::BernoulliProblem) =
    as((p = as𝕀, ))

# Add data

model = BernoulliProblem(; n = 9, obs = rand(Binomial(9, 2/3), 3))

# Make the type callable with the parameters *as a single argument*.

function (model::BernoulliProblem)(θ)
    @unpack n, obs = model        # extract the data
    @unpack p = θ
    loglikelihood(Binomial(n, p), obs)
end

# Use a flat priors (the default, omitted) for α

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Sample 4 chains

a3d = Array{Float64, 3}(undef, 1000, 1, 4);
for j in 1:4
  global results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000;
    #initialization = (q = (p = 0.2,)),
    #reporter = NoProgressReport()
    )
  global posterior = P.transformation.(results.chain)

  for i in 1:1000
    a3d[i, 1, j] = values(posterior[i].p)
  end
end

# Create MCMCChains object

parameter_names = ["p"]
sections =   Dict(
  :parameters => parameter_names,
)
chns = create_mcmcchains(a3d, parameter_names, sections, start=1)
show(chns)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display

println()
DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics)


text:
# Deprecate UnitQuaternion => QuatRotation
Base.@deprecate_binding UnitQuaternion QuatRotation true


text:
module ShipstabNURBS

using Plots

export BasisFuns,FindSpan,
nrbplot3,nrbplot,
bsplevalc,bsplevals,bsplder,
nrbevalc,nrbevals,nrbder,
nrbcircle,nrbcylinder,
revolX,revolX_half

function FindSpan(u,U)
	#
	# Short description:
	# 	computes the span index and the multiplicity of a a value u using a
	# 	binary search algorithm. Only for clamped b-splines and NURBS.
	#
	# Input parameters:
	# 	u: value between 0 and 1 for which we want to find the span index and
	#   multiplicity
	#
	#   U: the knot vector
	#
	# Output paramters:
	#   SpanIndex: the span index of u. if u belongs to the open interval
	#   [Uk,Uk+1[, we say that it belongs to the k-th knot index. In MatLab,
	#   indexes start from 1 and not 0. So to be mathematically correct,
	#   FindSpan computs the span index + 1.
	#
	#   s: the multiplicity of u. If u is no knot, it has a multiplicity of 0
	#
	# ------------------------------------------------------------------------

	# % first, we delete the first and last elements of the knot vector to avoid
	# % problems while performing de Boor's algorithm. (later, we add 1 to the
	# % found span index because we deleted the first element)
	U = U[2:end-1]
	s=0
	m = length(U) - 1
	low = 1
	high = m+1
	mid = round(Int,(low+high)/2)
	if isempty(findall(isequal(u),U)) == 0 && u!=U[end]
    	s = length(findall(isequal(u),U))
    	SpanIndex = findlast(isequal(u),U) + 1 # (+1, because we set U(1) = [])
    	return Int(SpanIndex),s
	end

	if U[end] == u
    	SpanIndex = findfirst(isequal(u),U)
    	return Int(SpanIndex),s
	end

	SpanIndex = mid
	while u < U[Int(mid)] || u >= U[Int(mid+1)]
    	if u < U[Int(mid)]
        	high = mid
    	else
        	low = mid
		end
    	mid = round((high+low)/2)
    	if mid == SpanIndex
        	mid = mid - 1
    	end
   		SpanIndex = mid
	end
	SpanIndex = SpanIndex + 1
	return Int(SpanIndex),s
end

function BasisFuns(n,p,u,U)
	#
	# Short description:
	#
	# compute the coefficients N0,p(u), N1,p(u), ..., Nn,p(u) for u in [0,1]
	#
	# Input parameters:
	#
	# n: there are n+1 control points
	# p: degree of the curve
	# u: value of the parameter u
	# U: knot vector
	#
	# Output parameters:
	#
	# Array containing the n+1 values of basis functions for a given u
	#
	# (algorithm explained in notes)
	# ------------------------------------------------------------------------
	k = FindSpan(u,U)[1]
	N = zeros(1,n+1)
	N[k] = 1
	for d = 1:p
	    if k-d > 0
    	    N[k-d] = ((U[k+1] - u)/(U[k+1] - U[k-d+1]))*N[k-d+1]
	        if k-d == n+1
	            return N
	        end
	        for i = k-d+1:k-1
	            N[i] = ((u - U[i])/(U[i+d]-U[i]))*N[i] + ((U[i+d+1] - u)/(U[i+d+1] - U[i+1]))*N[i+1]
	        end
	        N[k] = ((u - U[k])/(U[k+d] - U[k]))*N[k]
	    else
	        for i = 1:k-1
	            N[i] = ((u - U[i])/(U[i+d]-U[i]))*N[i] + ((U[i+d+1] - u)/(U[i+d+1] - U[i+1]))*N[i+1]
	        end
	        N[k] = ((u - U[k])/(U[k+d] - U[k]))*N[k]
	    end
	end
	return N
end

function bsplevalc(u,U,p,P)
	# U [1 x m+1]
	#
	# Short description:
	#   evaluates a (clamped) B spline curve for a given value for u using De
	#   Boors's algorithm
	#
	# Input parameters:
	#   u: a value in the knot vector U between 0 and 1
	#
	#   p: the degree of the NURBS curve
	#
	#   P: a [1 X n+1] vector with the cartesian x OR y coördinates
	#   of the n+1 control points
	#
	#   W: matrix containing the weights corresponding to the control points
	#   (see function nrbplot)
	#
	#   U: Knot vector [1 X m+1]
	#
	# Output parameters:
	#   C: the cartesian coördinates of the point on the B spline curve for a
	#   given value for u
	# ------------------------------------------------------------------------

	k,s = FindSpan(u,U)
	h = p-s
	if h <= 0
    	C = P[k-s]
    	return C
	end

	B = zeros(p-s+1,h+1)
	B[:,1] = P[k-p:k-s]

	for r = 1:h
    	for i = r+1:h+1
      		a = (u - U[i+k-p-1])/(U[i+k-r] - U[i+k-p-1])
      		B[i,r+1] = (1-a)*B[i-1,r]+a*B[i,r]
    	end
	end
	C = B[h+1,h+1]
	return C
end

function bsplevals(u,U,v,V,p,q,P)
	#
	# Short description:
	#   evaluates a (clamped) B-spline surface for given values for u and v using
	#   De Boors's algorithm
	#
	# Input parameters:
	#   u,v: a value in the knot vectors U and V between 0 and 1
	#
	#   p,q: the degree in respectively the u and v direction
	#
	#   P: cartesian coordinates of the control points (x,y or z direction, see
	#   input parameters of the function nrbplot3)
	#
	#   W: matrix containing the weights corresponding to the control points
	#   (see function nrbplot3)
	#
	# Output parameters:
	#   C: the cartesian coördinates of the point on the B spline surface for the
	#   given values for u and v
	#
	# ------------------------------------------------------------------------

	c,s = FindSpan(u,U)
	Q = zeros(1,p-s+length(P))
	for i = c-p:c-s
    	Q[i] = bsplevalc(v,V,q,P[i,:])
	end
	C = bsplevalc(u,U,p,Q)
	return C
end

function bsplder(P,U,u)
	Uprime = U[2:end-1]
	Q = zeros(size(P,1),size(P,2)-1)
	n = size(P,2) - 1
	m = length(U) - 1
	p = m-n-1

	for i = 1:size(Q,1)
    	for j = 1:size(Q,2)
	        Q[i,j] = p *(P[i,j+1]-P[i,j])/(U[j+p+1]-U[j+1])
	    end
	end

	T = zeros(1,size(Q,1))
	for i = 1:size(Q,1)
	    T[i] = bsplevalc(u,Uprime,p-1,Q[i,:])
	end
	return T
end

function nrbevalc(u,U,p,P,W)
	#
	# Short description:
	#   evaluates a (clamped) NURBS curve for a given value for u using De
	#   Boors's algorithm
	#
	# Input parameters:
	#   u: a value in the knot vector U between 0 and 1
	#
	#   p: the degree of the NURBS curve
	#
	#   P: a [2 X n+1] matrix with on the first row the cartesian x-coördinates
	#   of the n+1 control points and on the second row the associated
	#   y-coördinates.
	#
	#   W: matrix containing the weights corresponding to the control points
	#   (see function nrbplot)
	#
	#   U: Knot vector [1 X m+1]
	#
	# Output parameters:
	#   C: the cartesian coördinates of the point on the NURBS curve for a
	#   given value for u
	#
	#   Wn: the weight of the calculated point on the NURBS curve (needed for
	#   the calculation of NURBS surfaces)
	# ------------------------------------------------------------------------

	Pw = P.*W
	Pw = vcat(Pw, W)
	Qw = zeros(size(Pw,1),1)
	for i = 1:size(Pw,1)
    	Qw[i] = bsplevalc(u,U,p,Pw[i,:])
	end
	C = Qw[1:end-1]./Qw[end]
	Wn = Qw[end]
	return C,Wn
end

function nrbder(P,U,W,u)
	Pw = [P.*W;W]
	Uprime = U[2:end-1]
	Q = zeros(size(Pw,1),size(Pw,2)-1)
	n = size(Pw,2) - 1
	m = length(U) - 1
	p = m-n-1

	for i = 1:size(Q, 1)
	    for j = 1:size(Q,2)
	        Q[i,j] = p *(Pw[i,j+1]-Pw[i,j])/(U[j+p+1]-U[j+1])
	    end
	end

	Ader = zeros(1,size(Q,1))
	for i = 1:length(Ader)
	    Ader[i] = bsplevalc(u,Uprime,p-1,Q[i,:])
	end

	C,Wn = nrbevalc(u,U,p,P,W)
	T = zeros(1,length(Ader)-1)
	for i = 1:length(Ader)-1
	    T[i] = (Ader[i]-Ader[end]*C[i])/Wn
	end
	return T
end

function nrbevals(u,U,v,V,p,q,P,W)
	#
	# Short description:
	#   evaluates a (clamped) NURBS surface for given values for u and v using
	#   De Boors's algorithm
	#
	# Input parameters:
	#   u,v: a value in the knot vectors U and V between 0 and 1
	#
	#   p,q: the degree in respectively the u and v direction
	#
	#   P: cartesian coordinates of the control points (x,y or z direction, see
	#   input parameters of the function nrbplot3)
	#
	#   W: matrix containing the weights corresponding to the control points
	#   (see function nrbplot3)
	#
	# Output parameters:
	#   C: the cartesian coördinates of the point on the NURBS surface for the
	#   given values for u and v
	#
	# ------------------------------------------------------------------------

	c,s = FindSpan(u,U)
	Q = zeros(2,p-s+length(P))
	for i = c-p:c-s
    	Q[1,i] = nrbevalc(v,V,q,P[i,:]',W[i,:]')[1][1]
		Q[2,i] = nrbevalc(v,V,q,P[i,:]',W[i,:]')[2]
	end
	C, ~ = nrbevalc(u,U,p,Q[1,:]',Q[2,:]')
	return C
end

function nrbplot(P,W,U,npts)
	#
	# Short description:
	#
	#   Constructs a (clamped) NURBS curve
	#
	# Input parameters:
	#
	#   P: a [2 X n+1] matrix with on the first row the cartesian x-coördinates
	#   of the n+1 control points and on the second row the associated
	#   y-coördinates.
	#
	#   W: a [n+1] line vector containing the weights for each control point
	#
	#   U: knot vector
	#
	#   npts: the number of points for which the NURBS curve will be calculated
	#
	# Output parameters:
	#
	#   C: a matrix with on the first row the calculated cartesian
	#   x-coördinates and on the second row the associated y-coördinates
	# ----------------------------------------------------------------------

	n = size(P,2) - 1
	m = length(U) - 1
	p = m-n-1
	u = range(0,1,length=npts)
	C = zeros(2,length(u))

	for i = 1:length(u)
    	C[:,i], ~ = nrbevalc(u[i],U,p,P,W)
	end

	plot(P[1,:],P[2,:], line=(:dash,2))
	plot!(C[1,:],C[2,:], line=(:dot,2))
	return C
end

function nrbcircle(npts)
	#
	# Short description:
	#   construcs a NURBS circle evaluated at npts points
	#
	# Output parameters:
	#   P: control points
	#   W: weight vector
	#   U: knot vector
	#   C: carthesian coördinates of the calculated points
	#
	# (see input parameters of nrbplot for a more detailed description)
	#
	#   --------------------------------------------------------------------

	P = [1 1 0 -1 -1 -1 0 1 1;0 1 1 1 0 -1 -1 -1 0]
	W = [1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1]
	U = [0 0 0 1/4 1/4 1/2 1/2 3/4 3/4 1 1 1]

	C = nrbplot(P,W,U,npts)
	return P,W,U,C
end

function nrbplot3(Px,Py,Pz,W,U,V,nptsu,nptsv,flag_plot)
	# Short description:
	#
	#   Constructs a (clamped) NURBS surface
	#
	# Input parameters:
	#
	#   Px, Py, Pz: [m+1 X n+1] matrices with the cartesian coordinates of the
	#   control points for respectively the x, y and z-axis.
	#
	#          ------------------- > U direction
	#          |p1,1        p1,n+1 |
	#          | .     .           |
	#          | .      .          |
	#          | .       .         |
	#          |pm+1,1     pm+1,n+1|
	#          ----------------
	#          |
	#          v V direction
	#
	#   W: the weight assigned to each control point, [m+1 X n+1] matrix.
	#
	#   U,V: knot vectors in respectively the u- and v-direction
	#
	#   nptsu, nptsv: the number of points that subdivide the parametric space
	#   in respectively the u- and v-direction. It determines the fineness of
	#   the mesh in the physical space.
	#
	# Output parameters:
	#
	#   Cx,Cy,Cz: cartesian coördinates of the NURBS surface in respectively
	#   the x,y and z direction
	#   --------------------------------------------------------------------

	n = size(Px,2) - 1
	m = size(Px,1) - 1
	h = length(U) - 1
	k = length(V) - 1

	q = k-n-1
	p = h-m-1

	u = range(0,1,length=nptsu)
	v = range(0,1,length=nptsv)

	lu = length(u)
	lv = length(v)
	Cx = zeros(lu,lv)
	Cy = zeros(lu,lv)
	Cz = zeros(lu,lv)

	for i = 1:length(u)
    	for j = 1:length(v)
        	Cx[i,j] = nrbevals(u[i],U,v[j],V,p,q,Px,W)[1]
        	Cy[i,j] = nrbevals(u[i],U,v[j],V,p,q,Py,W)[1]
        	Cz[i,j] = nrbevals(u[i],U,v[j],V,p,q,Pz,W)[1]
    	end
	end
	if flag_plot == 1
		plot(Cx[1,:],Cy[1,:],Cz[1,:])
    	for i = 2:length(u)
        	plot!(Cx[i,:],Cy[i,:],Cz[i,:])
    	end
    	for i =1:length(v)
        	plot!(Cx[:,i],Cy[:,i],Cz[:,i])
    	end
     	# for i = 1:size(Px,1)
      #    	plot3(Px(i,:),Py(i,:),Pz(i,:),'ko-','Linewidth',2)
     	# end
     	# for i = 1:size(Px,2)
      #    	plot3(Px(:,i),Py(:,i),Pz(:,i),'ko-','Linewidth',2)
     # 	# end
    	# xaxis!("x-axis")
    	# yaxis!("y-axis")
    	# zaxis!("z-axis")
	end
	return Cx, Cy, Cz
end

function nrbcylinder(nptsu,nptsv)
	#
	# Short description:
	#   constructs a NURBS cylinder evaluated at nptsu x nptsv points
	#
	# Output parameters:
	#
	#   see input parameters of the function nrbplot3 for a detailed
	#   description
	#
	#   --------------------------------------------------------------------

	Pz = [1 1 0 -1 -1 -1 0 1 1;1 1 0 -1 -1 -1 0 1 1]'
	Py = [0 1 1 1 0 -1 -1 -1 0;0 1 1 1 0 -1 -1 -1 0]'
	Px = [-1 -1 -1 -1 -1 -1 -1 -1 -1; 1 1 1 1 1 1 1 1 1]'

	W = [1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1; 1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1]'
	U = [0 0 0 1/4 1/4 1/2 1/2 3/4 3/4 1 1 1]
	V = [0 0 1 1]

	nrbplot3(Px,Py,Pz,W,U,V,nptsu,nptsv,1)
	return Px,Py,Pz,W,U,V
end

function revolX(P,W,U,nptsu,nptsv,flag_plot)
	# Constructs a surface of revolution by means of a b-spline curve in 2D.

	Pj = zeros(3,size(P,2))
	Pj[1,:] = P[1,:]
	Pj[3,:] = P[2,:]
	Wj = W

	Pi = [0 0 0 0 0 0 0 0 0;1 1 0 -1 -1 -1 0 1 1;0 1 1 1 0 -1 -1 -1 0]
	Wi = [1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1]
	V = [0 0 0 1/4 1/4 1/2 1/2 3/4 3/4 1 1 1]

	nPj = size(Pj,2)
	nPi = size(Pi,2)

	xj = Pj[1,:]
	zj = Pj[3,:]

	Px = zeros(nPj,nPi)
	Py = Px
	Pz = Px

	for i = 1:length(xj)
    	Px[i,:] = xj[i].* ones(1,nPi)
    	Py[i,:] = zj[i].* Pi[2,:]
    	Pz[i,:] = zj[i].* Pi[3,:]
	end

	W = zeros(size(Px))
	for i = 1:nPj
    	W[i,:] = Wj[i].*Wi
	end
	nrbplot3(Px,Py,Pz,W,U,V,nptsu,nptsv,flag_plot)
	return Px,Py,Pz,U,V,W
end

function revolX_half(P,W,U,nptsu,nptsv,flag_plot)
	# Constructs a surface of revolution by means of a b-spline curve in 2D.

	Pj = zeros(3,size(P,2))
	Pj[1,:] = P[1,:]
	Pj[3,:] = P[2,:]
	Wj = W

	Pi = [0 0 0 0 0 0 0 0 0;1 1 0 -1 -1 -1 0 1 1;0 1 1 1 0 -1 -1 -1 0]
	Wi = [1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1 sqrt(2)/2 1]
	Pi[:,1:4] = []
	Wi[1:4] = []
	V = [0 0 0 1/2 1/2 1 1 1]

	nPj = size(Pj,2)
	nPi = size(Pi,2)

	xj = Pj[1,:]
	zj = Pj[3,:]

	Px = zeros(nPj,nPi)
	Py = Px
	Pz = Px

	for i = 1:length(xj)
    	Px[i,:] = xj[i].* ones(1,nPi)
    	Py[i,:] = zj[i].* Pi[2,:]
    	Pz[i,:] = zj[i].* Pi[3,:]
	end

	W = zeros(size(Px))
	for i = 1:nPj
    	W[i,:] = Wj[i].*Wi
	end
	nrbplot3(Px,Py,Pz,W,U,V,nptsu,nptsv,flag_plot)
	return Px,Py,Pz,U,V,W
end

end # module

max_stars_repo_name: barche/ShipstabNURBS

max_issues_repo_name: barche/ShipstabNURBS

max_forks_repo_name: barche/ShipstabNURBS


text:
"""
# Object
In access control systems, an object can be any resource (file, record, application, ...) on which operations may be performed (read, modify, hunt, fish, ...).
"""
mutable struct Object <: AbstractObject
    label::String;
end;
function Object(;n::Int = 1)::Vector{Object}
    useless(x) = Object("an object");
    broadcast(useless, repeat([0], n));
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
__precompile__()

module RLEVectors

using Compat
import Compat.String
using Reexport
@reexport using DataFrames

# RLEVector type
export RLEVector, FloatRle, IntegerRle, BoolRle, StringRle, rfirst, rwidth, rlast, rvalue, nrun, similar, collect, similar
import Base: show, length, size, start, next, done, Forward, first, last, step, convert, similar, collect, isequal

# collections
import Base: eltype, unique, minimum, maximum, vcat, pop!, push!, shift!, unshift!, insert!, deleteat!, splice!, resize!, empty!, endof, maxabs, minabs, any, all, in, intersect
export       eltype, unique, minimum, maximum, vcat, pop!, push!, shift!, unshift!, insert!, deleteat!, splice!, resize!, empty!, endof, maxabs, minabs, any, all, in, intersect
export shove!, deleterun!, decrement_run!

# indexing
import Base: getindex, setindex!
export getindex, setindex!, ind2run, setrun!, ind2runcontext, head, tail, RLERangesIterator, each

# describe
import Base: isempty, ndims
export ndims, endtype

# group_generics
import Base: +, -, *, /, ^, .+, .-, .*, ./, .^, div, mod, fld, rem
import Base: ==, >, <, !=, <=, >=, .==, .>, .<, .!=, .<=, .>=, &, |
import Base: abs, sign, sqrt, ceil, floor, trunc, cummax, cummin, cumprod, cumsum, log, log10, log2, log1p, acos, acosh, asin, asinh, atan, atanh
import Base: exp, expm1, cos, cosh, sin, sinh, tan, tanh, gamma, lgamma, digamma, trigamma
import Base: max, min, range, prod, sum, any, all, mean
import Base: in
import Base: indexin, findin, median, findmin, findmax
export .+, .-, .*, ./, .^, div, mod, fld, rem, ==, >, <, !=, <=, >=, .==, .>, .<, .!=, .<=, .>=, &, |
export abs, sign, sqrt, ceil, floor, trunc, cummax, cummin, cumprod, cumsum, log, log10, log2, log1p, acos, acosh, asin, asinh, atan, atanh
export exp, expm1, cos, cosh, sin, sinh, tan, tanh, gamma, lgamma, digamma, trigamma
export max, min, range, prod, sum, any, all, mean
export in
export indexin, findin, median, findmin, findmax
export findin2

# math
import StatsBase: mode, countmap
export mode, countmap

# ranges
export disjoin, disjoin_length, ree, inverse_ree, numruns, rangeMeans

# utils
import Base: searchsortedfirst
export rep, searchsortedfirst

# sorting
import Base.Order: Ordering
import Base.Sort: QuickSortAlg
import Base: sort, sort!, issorted, reverse, reverse!, sortperm, Algorithm
export       sort, sort!, issorted, reverse, reverse!, sortperm, permute_runs

# data frames
import DataFrames: AbstractDataFrame, DataFrame, Index, head, tail
export RLEDataFrame

### Includes
include("utils.jl")
include("runs.jl")
include("RLEVector-type.jl")
include("RLEVectorList-type.jl")
include("RLEDataFrame-type.jl")
include("ranges.jl")
include("describe.jl")
include("indexing.jl")
include("group_generics.jl")
include("collections_api.jl")
include("math.jl")
include("sorting.jl")

if VERSION >= v"0.4.0"
    include("precompile.jl")
    _precompile_()
end

end # Module RLEVectors


text:
# This file includes further code that is added to the module JAC.Basics

using  Dates, Printf, Interact,  LinearAlgebra
using  JAC, ..Basics, ..Defaults,  ..Radial, ..ManyElectron, ..Nuclear   ## , ..RadialIntegrals, ..InteractionStrength, ..AngularMomentum

include("jac-add.jl")
include("jac-analyze.jl")
include("jac-compute.jl")
include("jac-determine.jl")
include("jac-diagonalize.jl")
include("jac-display.jl")
include("jac-estimate.jl")
include("jac-exclude.jl")
include("jac-extract.jl")
include("jac-evaluate.jl")
include("jac-generate.jl")
include("jac-integrate.jl")
include("jac-interpolate.jl")
include("jac-is.jl")
include("jac-merge.jl")
include("jac-modify.jl")
include("jac-perform.jl")
include("jac-provide.jl")
include("jac-plot.jl")
include("jac-read.jl")
include("jac-recast.jl")
include("jac-sort.jl")
include("jac-tabulate.jl")
include("jac-tools.jl")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule wget2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("wget2")
JLLWrappers.@generate_main_file("wget2", UUID("8726a5da-ff56-5c52-916d-726e18caca27"))
end  # module wget2_jll


text:
using Memoize
using Random
using Test
using SpinGlassTensors

D = 10
d = 3
sites = 5
T = Float64

@testset "Results of left_env are correctly cached" begin
    # Cache is global, therefore we need to reset it to obtain meaningful
    # results here.
    empty!(memoize_cache(left_env))

    Random.seed!(69)

    ψ = randn(MPS{T}, sites, D, d)
    ϕ = copy(ψ)

    σ = [1, 2, 3]
    η = [1, 2, 3]

    envs = [left_env(mps, v) for mps in (ψ, ϕ) for v ∈ (σ, η)]

    @testset "Calls made with equal arguments result in a cache hit" begin
        # Should result in 4 calls total (top-level + 3 recursive)
        @test length(memoize_cache(left_env)) == 4
    end

    @testset "Cached results are equal to the ones computed during first call" begin
        @test all(env->env==envs[1], envs)
    end
end

@testset "Results of right_env are correctly cached" begin
    empty!(memoize_cache(right_env))

    Random.seed!(42)

    ψ = randn(MPS{T}, sites, D, d)
    ϕ = copy(ψ)

    W = randn(MPO{T}, sites, D, d)
    V = copy(W)

    σ = [2, 1, 1, 2, 1]
    η = copy(σ)

    envs = [
        right_env(mps, mpo, v) for mps in (ψ, ϕ) for mpo ∈ (W, V) for v ∈ (σ, η)
    ]

    @testset "Calls made with equal arguments result in a cache hit" begin
        # Should result in 6 calls total (top-level + 5 recursive)
        @test length(memoize_cache(right_env)) == 6
    end

    @testset "Cached results are equal to the ones computed during first call" begin
        @test all(env->env==envs[1], envs)
    end
end