text:
using StatsBase
using UnPack

function rand_sample(P, q)
    L = Vector{Int64}(undef, q)
    for i in 1:q
        L[i] = sample(1:length(P), Weights(P))
        # L[i] = sample(1:sizeof(P), weights(P), 1, replace=true)
    end
    return L
end

function not_in_search_space(position, lb, ub)
    return any(position .<= lb) || any(position .>= ub)
end

function CCE(complex_pops, cce_params)
    # Step 1. Initialize
    @unpack q, alpha, beta, lb, ub, obj_func = cce_params
    n_pop = length(complex_pops)

    # Step 2. Assign weights
    P = [2*(n_pop+1-i) / (n_pop*(n_pop+1)) for i in 1:n_pop]

    # Calculate Population Range (Smallest Hypercube)
    new_lb = complex_pops[1].position
    new_ub = complex_pops[1].position
    for i in 2:n_pop
        new_lb = min.(new_lb, complex_pops[i].position)
        new_ub = max.(new_ub, complex_pops[i].position)
    end

    # CCE main loop
    for it in 1:beta
        # Step 3. Select parents
        L = rand_sample(P, q)
        B = complex_pops[L]

        # Step 4. Generate Offspring
        for k in 1:alpha
            # a) Sort population
            sorted_indexs = sortperm(B)
            sort!(B)
            L[:] = L[sorted_indexs]

            # Calculate the centroid
            g = zeros(length(lb))
            for i in 1:q-1
                g .= g .+ B[i].position
            end
            g .= g ./ (q-1)

            # b) Reflection step
            reflection = deepcopy(B[end])
            reflection.position = 2 .* g .- B[end].position # newly generated point using reflection
            if not_in_search_space(reflection.position, lb, ub)
                reflection.position = uniform_rand(new_lb, new_ub)
            end
            reflection.cost = obj_func(reflection.position)

            if reflection.cost < B[end].cost
                B[end] = deepcopy(reflection)
            else # Contraction
                contraction = deepcopy(B[end])
                contraction.position = (g .+ B[end].position) ./ 2
                contraction.cost = obj_func(contraction.position)

                if contraction.cost < B[end].cost
                    B[end] = deepcopy(contraction)
                else
                    B[end].position = uniform_rand(new_lb, new_ub)
                    B[end].cost = obj_func(B[end].position)
                end
            end

        end

        complex_pops[L] = B
    end
    return complex_pops
end


text:
# todo: potentially add caching. only execute if input filename has been modified after the splitdir contents
function split_raw_data(j)
  for filename in j.inputs
    wrkrs = workers()
    input_filename = joinpath(datadir, filename)
    !isfile(input_filename) && return true # not in data/ but already in split/
    info("Splitting $input_filename using $(length(wrkrs)) workers")
    output_filenames = Array(String, length(wrkrs))
    for i in 1:length(wrkrs)
      output_dir = joinpath(splitdir, string(wrkrs[i]))
      !isdir(output_dir) && mkdir(output_dir)
      output_filenames[i] = joinpath(output_dir, filename)
    end
    # open the files for writing (1 file per core)
    handles = [open(f, "w") for f in output_filenames]
    handle = 1
    for ln in readlines(input_filename)
      current_output_file = output_filenames[handle]
      # println("current handle: $(handle) and current output filename $(output_filenames[handle])")
      current_handle = handles[handle]
      write(current_handle, ln) # write line into current handle
      handle += 1 # get next handle
      if handle > length(handles)
        handle = 1 # cycle back to first handle once handles have been exhausted
      end
    end
    map(close, handles) # close all output handles
  end
  return true
end

job_output_dir(j) = joinpath(outputdir, string(j.jobid))

function create_job_area(j)
  info("Creating job area...")
  joboutputdir = job_output_dir(j)
  isdir(joboutputdir) && rm(joboutputdir, force=true, recursive=true)
  mkdir(joboutputdir)
end

function read_sink(j, phase)
  inputfile = joinpath(job_output_dir(j), string(myid()) * "." * phase)
  open(inputfile)
end

function write_sink(j)
  outputfile = joinpath(splitdir, string(myid()), j.outputfilename)
  isfile(outputfile) && rm(outputfile)
  open(outputfile, "w")
end

ext: jl

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
export MonteCarloLearner,
       FIRST_VISIT,
       EVERY_VISIT,
       NO_SAMPLING,
       ORDINARY_IMPORTANCE_SAMPLING,
       WEIGHTED_IMPORTANCE_SAMPLING

using .Utils: CachedSampleAvg
using StatsBase: countmap

abstract type AbstractVisitType end
struct FirstVisit <: AbstractVisitType end
struct EveryVisit <: AbstractVisitType end
const FIRST_VISIT = FirstVisit()
const EVERY_VISIT = EveryVisit()

abstract type AbstractSamplingType end
struct NoSampling <: AbstractSamplingType end
struct OrdinaryImportanceSampling <: AbstractSamplingType end
struct WeightedImportanceSampling <: AbstractSamplingType end
const NO_SAMPLING = NoSampling()
const ORDINARY_IMPORTANCE_SAMPLING = OrdinaryImportanceSampling()
const WEIGHTED_IMPORTANCE_SAMPLING = WeightedImportanceSampling()

"""
    MonteCarloLearner(; approximator::A, γ = 1.0, α = 1.0, kind = FIRST_VISIT, sampling = NO_SAMPLING, returns = CachedSampleAvg())
"""
struct MonteCarloLearner{T,A,R,S} <: AbstractLearner
    approximator::A
    γ::Float64
    α::Float64
    returns::R

    MonteCarloLearner(
        ;
        approximator::A,
        γ = 1.0,
        α = 1.0,
        kind = FIRST_VISIT,
        sampling = NO_SAMPLING,
        returns = CachedSampleAvg(),
    ) where {A} =
        new{typeof(kind),A,typeof(returns),typeof(sampling)}(approximator, γ, α, returns)
end

(learner::MonteCarloLearner)(s) = learner.approximator(s)
(learner::MonteCarloLearner)(obs::Observation) = learner.approximator(get_state(obs))
(learner::MonteCarloLearner)(obs::Observation, a) = learner.approximator(get_state(s), a)

function update!(
    learner::MonteCarloLearner{<:FirstVisit,<:AbstractVApproximator,<:Any,<:NoSampling},
    transitions,
)
    states, rewards = transitions.states, transitions.rewards
    V, γ, α, Returns, G, T = learner.approximator,
        learner.γ,
        learner.α,
        learner.returns,
        0.0,
        length(states)
    seen_states = countmap(states)

    for t = T:-1:1
        S, R = states[t], rewards[t]
        G = γ * G + R
        if seen_states[S] == 1  # first visit
            update!(V, S => α * (Returns(S, G) - V(S)))
            delete!(seen_states, S)
        else
            seen_states[S] -= 1
        end
    end
end

function update!(
    learner::MonteCarloLearner{
        <:FirstVisit,
        <:AbstractVApproximator,
        <:CachedSampleAvg,
        <:OrdinaryImportanceSampling,
    },
    transitions,
    weights,
)
    states, rewards = transitions.states, transitions.rewards
    V, γ, α, Returns, G, ρ, T = learner.approximator,
        learner.γ,
        learner.α,
        learner.returns,
        0.0,
        1.0,
        length(states)
    seen_states = countmap(states)

    for t = T:-1:1
        S, R = states[t], rewards[t]
        G = γ * G + R
        ρ *= weights[t]
        if seen_states[S] == 1  # first visit
            update!(V, S => α * (Returns(S, ρ * G) - V(S)))
            delete!(seen_states, S)
        else
            seen_states[S] -= 1
        end
    end
end

function update!(
    learner::MonteCarloLearner{
        <:FirstVisit,
        <:AbstractVApproximator,
        <:Tuple{CachedSum,CachedSum},
        <:WeightedImportanceSampling,
    },
    transitions,
    weights,
)
    states, rewards = transitions.states, transitions.rewards
    V, γ, α, (G_cached, ρ_cached), G, ρ, T = learner.approximator,
        learner.γ,
        learner.α,
        learner.returns,
        0.0,
        1.0,
        length(states)
    seen_states = countmap(states)

    for t = T:-1:1
        S, R = states[t], rewards[t]
        G = γ * G + R
        ρ *= weights[t]
        if seen_states[S] == 1  # first visit
            numerator = G_cached(S, ρ * G)
            denominator = ρ_cached(S, ρ)
            val = denominator == 0 ? 0 : numerator / denominator
            update!(V, S => α * (val - V(S)))
            delete!(seen_states, S)
        else
            seen_states[S] -= 1
        end
    end
end

function update!(
    learner::MonteCarloLearner{<:EveryVisit,<:AbstractVApproximator},
    transitions,
)
    states, rewards = transitions.states, transitions.rewards
    α, γ, V, Returns, G = learner.α, learner.γ, learner.approximator, learner.returns, 0.0
    for (s, r) in Iterators.reverse(zip(states, rewards))
        G = γ * G + r
        update!(V, s => α * (Returns(s, G) - V(s)))
    end
end

function update!(
    learner::MonteCarloLearner{<:FirstVisit,<:AbstractQApproximator},
    transitions,
)
    states, actions, rewards = transitions.states, transitions.actions, transitions.rewards
    α, γ, Q, Returns, G, T = learner.α,
        learner.γ,
        learner.approximator,
        learner.returns,
        0.0,
        length(states)
    seen_pairs = countmap(zip(states, actions))

    for t = T:-1:1
        S, A, R = states[t], actions[t], rewards[t]
        pair = (S, A)
        G = γ * G + R
        if seen_pairs[pair] == 1  # first visit
            update!(Q, pair => α * (Returns(pair, G) - Q(S, A)))
            delete!(seen_pairs, pair)
        else
            seen_pairs[pair] -= 1
        end
    end
end

function update!(
    learner::MonteCarloLearner{<:EveryVisit,<:AbstractQApproximator},
    transitions,
)
    states, actions, rewards = transitions.states, transitions.actions, transitions.rewards
    α, γ, Q, Returns, G = learner.α, learner.γ, learner.approximator, learner.returns, 0.0
    for (s, a, r) in Iterators.reverse(zip(states, actions, rewards))
        G = γ * G + r
        update!(Q, (s, a) => α * (Returns((s, a), G) - Q(s, a)))
    end
end

function extract_transitions(
    buffer::EpisodeTurnBuffer,
    ::MonteCarloLearner{T,A},
) where {T,A<:AbstractQApproximator}
    if isfull(buffer)
        @views (
            states = state(buffer)[1:end-1],
            actions = action(buffer)[1:end-1],
            rewards = reward(buffer)[2:end],
        )
    else
        nothing
    end
end


text:
facts("Reward types") do
    @fact AbstractArrayReward <: AbstractReward => true
    @fact ArrayReward <: AbstractArrayReward => true
    @fact SparseReward <: AbstractArrayReward => true
end

facts("ArrayReward construction") do
    @fact typeof(Reward(rand(5))) => ArrayReward{Float64,1}
    @fact typeof(Reward(rand(5, 3))) => ArrayReward{Float64,2}
    @fact typeof(Reward(rand(5, 5, 3))) => ArrayReward{Float64,3}
    @fact_throws ErrorException Reward(rand(5, 4, 3))
    @fact typeof(Reward(sprand(5, 3, 1/3))) => SparseReward{Float64,Int}
end

facts("1d ArrayReward methods") do
    R = ArrayReward([1, 3, 1, 1, 2])
    @fact R[2] => 3
    @fact reward(R, 2) => 3
    @fact reward(R, 2, 3) => 3
    @fact reward(R, 2, 1, 3) => 3
    @fact_throws BoundsError R[2, 3]
    @fact_throws BoundsError R[2, 1, 3]
    @fact num_states(R) => 5
    @fact_throws ErrorException num_actions(R)
end

facts("2d ArrayReward methods") do
    R = ArrayReward([1 2 3; 2 3 1; 3 1 2; 2 1 3; 1 3 2])
    @fact R[2, 3] => 1
    @fact reward(R, 2, 3) => 1
    @fact reward(R, 2, 1, 3) => 1
    @fact_throws BoundsError R[2, 1, 3]
    @fact_throws MethodError reward(R, 2)
    @fact num_states(R) => 5
    @fact num_actions(R) => 3
end

facts("3d ArrayReward methods") do
    R = ArrayReward(
        cat(
            3,
            [1 2 3; 2 3 1; 3 1 2],
            [2 3 1; 3 1 2; 2 1 3],
        )
    )
    @fact R[2, 3, 2] => 2
    @fact reward(R, 2, 3, 2) => 2
    @fact num_states(R) => 3
    @fact num_actions(R) => 2
    @fact_throws MethodError reward(R, 2, 3)
    @fact_throws MethodError reward(R, 2)
end

facts("SparseReward methods") do
    R = SparseReward(sparse([1, 2, 3, 4, 5], [1, 2, 3, 2, 1], [1, 3, 1, 1, 2]))
    @fact R[2, 2] => 3
    @fact R[2, 1] => 0
    @fact reward(R, 2, 2) => 3
    @fact reward(R, 2, 1) => 0
    @fact reward(R, 2, 1, 2) => 3
    @fact num_states(R) => 5
    @fact num_actions(R) => 3
    @fact_throws MethodError reward(R, 2)
    @fact_throws BoundsError R[2, 1, 2]
end


text:
"""

The exported names are
$(EXPORTS)
"""
module Visualization

using DocStringExtensions
using WriteVTK
using WriteVTK.VTKCellTypes: VTK_VERTEX
using WriteVTK.VTKCellTypes: VTK_LINE
using WriteVTK.VTKCellTypes: VTK_TRIANGLE
using WriteVTK.VTKCellTypes: VTK_QUAD
using WriteVTK.VTKCellTypes: VTK_TETRA
using WriteVTK.VTKCellTypes: VTK_HEXAHEDRON
using WriteVTK.VTKCellTypes: VTK_QUADRATIC_TRIANGLE
using WriteVTK.VTKCellTypes: VTK_QUADRATIC_TETRA
using WriteVTK.VTKCellTypes: VTK_QUADRATIC_EDGE
using WriteVTK.VTKCellTypes: VTK_QUADRATIC_HEXAHEDRON
#using WriteVTK.VTKCellTypes: VTK_BIQUADRATIC_HEXAHEDRON
using WriteVTK.VTKCellTypes: VTK_QUADRATIC_QUAD
using WriteVTK.VTKCellTypes: VTK_BIQUADRATIC_QUAD
using WriteVTK.VTKCellTypes: VTK_WEDGE
using WriteVTK.VTKCellTypes: VTK_PYRAMID
using WriteVTK.VTKCellTypes: VTK_LAGRANGE_QUADRILATERAL
using WriteVTK.VTKCellTypes: VTK_LAGRANGE_TRIANGLE
using WriteVTK.VTKCellTypes: VTK_LAGRANGE_HEXAHEDRON
using WriteVTK.VTKCellTypes: VTK_LAGRANGE_TETRAHEDRON

using Gridap.Helpers
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Arrays
using Gridap.Polynomials
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.ReferenceFEs: _find_unique_with_indices
using Gridap.ReferenceFEs: SerendipityPolytope
using Gridap.FESpaces
using FillArrays
using Gridap.CellData

import Gridap.Geometry: get_reffes
import Gridap.Geometry: get_cell_type
import Gridap.Geometry: get_node_coordinates
import Gridap.Geometry: get_cell_nodes

import AbstractTrees

export writevtk
export createvtk
export write_vtk_file
export print_op_tree
export visualization_data

include("VisualizationData.jl")
include("Vtk.jl")

include("PrintOpTrees.jl")

end # module


text:
"Information container about the cells to run in a reactive call and any cells that will err."
struct CellTopology
	"Cells that form a directed acyclic graph, in topological order."
	runnable::Array{Cell, 1}
	"Cells that are in a directed cycle, with corresponding `ReactivityError`s."
	errable::Dict{Cell, ReactivityError}
end

"Return a `CellTopology` that lists the cells to be evaluated in a single reactive run, in topological order. Includes the given roots."
function dependent_cells(notebook::Notebook, roots::Array{Cell, 1})::CellTopology
	entries = Cell[]
	exits = Cell[]
	errable = Dict{Cell, ReactivityError}()

	function dfs(cell::Cell)
		if cell in keys(errable)
			return
		elseif cell in exits
			return
		elseif length(entries) > 0 && entries[end] == cell
			return # a cell referencing itself is legal
		elseif cell in entries
			currently_in = setdiff(entries, exits)
			cycle = currently_in[findfirst(currently_in .== [cell]):end]
			for cell in cycle
				errable[cell] = CyclicReferenceError(cycle)
			end
			return
		end

		push!(entries, cell)
		assigners = where_assigned(notebook, cell.symstate.assignments)
		if length(assigners) > 1
			for cell in assigners
				errable[cell] = MultipleDefinitionsError(cell, assigners)
			end
		end
		referencers = where_referenced(notebook, cell.symstate.assignments)
		dfs.(setdiff(union(assigners, referencers), [cell]))
		push!(exits, cell)
	end

	# reversing roots for a minor issue: mutually unrelated cells would otherwise be run in reverse order
	dfs.(reverse(roots))
	ordered = reverse(exits)
	CellTopology(setdiff(ordered, keys(errable)), errable)
end

function disjoint(a::Set, b::Set)
	!any(x in a for x in b)
end

"Return the cells that reference any of the given symbols. Recurses down functions calls, but not down cells."
function where_referenced(notebook::Notebook, symbols::Set{Symbol})::Array{Cell, 1}
	return filter(notebook.cells) do cell
		if !disjoint(symbols, cell.symstate.references)
			return true
		end
        for func in cell.symstate.funccalls
            if haskey(notebook.combined_funcdefs, func)
                if !disjoint(symbols, notebook.combined_funcdefs[func].references)
                    return true
                end
            end
		end
		return false
	end
end

"Return the cells that assign to any of the given symbols. Recurses down functions calls, but not down cells."
function where_assigned(notebook::Notebook, symbols::Set{Symbol})::Array{Cell, 1}
	return filter(notebook.cells) do cell
		if !disjoint(symbols, cell.symstate.assignments)
			return true
		end
        for func in cell.symstate.funccalls
            if haskey(notebook.combined_funcdefs, func)
                if !disjoint(symbols, notebook.combined_funcdefs[func].assignments)
                    return true
                end
            end
		end
		return false
	end
end

"Update the combined collection of function definitions, where multiple specialisations of a function are combined into a single `SymbolsState`."
function update_funcdefs!(notebook::Notebook)
	# TODO: optimise
	combined = notebook.combined_funcdefs = Dict{Symbol, SymbolsState}()

	for cell in notebook.cells
		for (func, symstate) in cell.symstate.funcdefs
			if haskey(combined, func)
				combined[func] = symstate ∪ combined[func]
			else
				combined[func] = symstate
			end
		end
	end
end

"Return all variables that a cell references, including those referenced through function calls."
function all_references(notebook::Notebook, cell::Cell)::Set{Symbol}
	calls = all_recursed_calls!(notebook, cell.symstate)
	filter!(calls) do call
		call in keys(notebook.combined_funcdefs)
	end
	return union(cell.symstate.references, (notebook.combined_funcdefs[func].references for func in calls)...)
end

"Return all variables that a cell assigns to, including mutable globals assigned through function calls."
function all_assignments(notebook::Notebook, cell::Cell)::Set{Symbol}
	calls = all_recursed_calls!(notebook, cell.symstate)
	filter!(calls) do call
		call in keys(notebook.combined_funcdefs)
	end
	return union(cell.symstate.assignments, (notebook.combined_funcdefs[func].assignments for func in calls)...)
end

"Return all functions called by a cell, and all functions called by those functions, et cetera."
function all_recursed_calls!(notebook::Notebook, symstate::SymbolsState, found::Set{Symbol}=Set{Symbol}())::Set{Symbol}
	for func in symstate.funccalls
		if func in found
			# done
		else
            push!(found, func)
            if haskey(notebook.combined_funcdefs, func)
                inner_symstate = notebook.combined_funcdefs[func]
                all_recursed_calls!(notebook, inner_symstate, found)
            end
		end
	end

	return found
end


text:
module TestInference

include("preamble.jl")

f1() = plan(
    ND(1) => (Index(:i),),
    ND(2) => (Index(:i), Index(:j)),
)

f2() = plan(
    ND(1) => (:i,),
    ND(2) => (:i, :j),
)

@testset "plan" begin
    if VERSION < v"1.2-"
        @test_broken_inferred f1()
        @test_broken_inferred f2()
    else
        @test_inferred f1()
        @test_inferred f2()
    end
end

@testset "axisfor" begin
    patterns = AccessPattern.((
        ND(1) => (:i,),
        ND(2) => (:i, :j),
    ))
    @test_inferred _axisfor(patterns, Index(:i))
    @test_inferred _axisfor(patterns, Index(:j))
end

end  # module


text:
#! format: off

abstract type AuxVariableType end

struct AuxVarKey{T <: AuxVariableType, U <: PSY.Component} <: OptimizationContainerKey
    aux_var_type::Type{T}
    device_type::Type{U}
end

function encode_key(::AuxVarKey{T, U}) where {T <: AuxVariableType, U <: PSY.Component}
    return Symbol("$(IS.strip_module_name(string(T)))_$(IS.strip_module_name(string(U)))")
end

struct TimeDurationOn <: AuxVariableType end
struct TimeDurationOff <: AuxVariableType end

""" Auxiliary Variable for Thermal Generation Models that solve for power above min"""
struct PowerOutput <: AuxVariableType end


text:
export GencodeGene,
    GencodeTranscript,
    GencodeExon,
    gencode_dir,
    gencode_load,
    gencode_load_file,
    gencode_locate,
    gencode_genes,
    gencode_transcript

include("config.jl")
include("index.jl")
include("search.jl")


text:
import Base.convert
import Base.==

const NODE = UInt8(0)
const LEAF = UInt8(1)

abstract type Tree end

struct Leaf <: Tree
    x::Vector{UInt8} # Bytes
end

==(a::Leaf, b::Leaf) = a.x == b.x

struct Node <: Tree
    x::Vector{Tree} # Leaf or Node
end

==(a::Node, b::Node) = a.x == b.x

Node() = Node([])

Base.push!(n::Node, y) = push!(n.x, y)

toint(x) = reinterpret(UInt32, x[4:-1:1])[1] ### TOREMOVE

function parseb(x)

    if x[1] == LEAF

        L = interpret(UInt32, x[2:5])

        bytes = x[6:5+L]
        leaf = Leaf(bytes)

        if length(x) == L + 5
            rest = []
        else
            rest = x[L+6:end]
        end

        return leaf, rest

    elseif x[2] == NODE

        N = interpret(UInt32, x[2:5])

        rest = x[6:end]

        node = Node()

        for i in 1:N
            head, tail = parseb(rest)
            push!(node, head)
            rest = tail
        end

        return node, rest
    end
end

decode(x::Vector{UInt8}) = parseb(x)[1]
decode(x::AbstractString) = decode(hex2bytes(replace(x, " "=>""))) # I could have optional arguments here as well

function tobin(leaf::Leaf)

    N = UInt32(length(leaf.x))

    Nbin = interpret(Vector{UInt8}, N)
    bin = UInt8[LEAF, Nbin..., leaf.x...]

    return bin
end

function tobin(node::Node)

    N = UInt32(length(node.x))
    Nbin = interpret(Vector{UInt8}, N)

    data = UInt8[]

    for n in node.x
        b = tobin(n)
        append!(data, b)
    end

    bin = UInt8[NODE, Nbin..., data...]

    return bin
end

encode(x::Tree) = tobin(x)

convert(::Type{T}, x::Leaf) where T <: Integer = interpret(T, x.x)

function convert(::Type{String}, x::Leaf)
    return String(copy(x.x))
end

function convert(::Type{Vector{T}}, x::Node) where T #<: Integer
    return T[convert(T, i) for i in x.x]
end

function convert(cfact::Type{T}, x::Node) where T <: Tuple
     return Tuple((convert(ci, xi) for (xi, ci) in zip(x.x, cfact.types)))
end

function Leaf(x::Signed)

    bytes = interpret(Vector{UInt8}, x)

    # Adding a redundant byte to ensure that the number is positive.
    if bytes[1] > 127
        return Leaf(UInt8[0, bytes...])
    else
        return Leaf(bytes)
    end
end

Leaf(x::Unsigned) = Leaf(interpret(Vector{UInt8}, x))

function Leaf(x::Integer, k::Integer)
    leaf = Leaf(x)
    N = findfirst(x -> x != UInt8(0), leaf.x)
    bytes = leaf.x[N:end]
    pad = k - length(bytes)

    newleaf = Leaf(UInt8[zeros(UInt8, pad)...,bytes...])
    return newleaf
end

function Leaf(x::AbstractString)
    bytes = Vector{UInt8}(x)
    return Leaf(bytes)
end

Tree(x::Any) = Leaf(x)
Tree(x::Node) = x
Tree(x::Leaf) = x
Tree(x::Tuple) = Node(x)

function Node(x::Tuple)
    node = Node()
    for i in x
        r = Tree(i)
        push!(node, r)
    end
    return node
end

############################ COMPOSITE TYPE PARSING ############################

using CryptoGroups: Hash, PGroup, ECGroup, Group, ElGamal, value, specialize, a, b, spec, generator, <|, name, ECPoint, field, gx, gy

(h::Hash)(t::Tree) = h(convert(Vector{UInt8}, t))  ### need to relocate

# Need something smarter in the end
bitlength(::Type{G}) where G <: PGroup = bitlength(modulus(G))
bitlength(x::PGroup) = bitlength(modulus(x))

bitlength(::Type{ECGroup{P}}) where P <: ECPoint = bitlength(modulus(field(P)))
bitlength(g::G) where G <: ECGroup = bitlength(G)

Tree(x::PGroup; L = bitlength(x)) = Leaf(value(x), div(L + 1, 8, RoundUp))

# Probably I will need to replace
convert(::Type{G}, x::Leaf) where G <: PGroup = convert(G, convert(BigInt, x))

### Note that only PrimeCurves are supported.
convert(::Type{G}, x::Node) where G <: ECGroup = G <| convert(Tuple{BigInt, BigInt}, x)

function Tree(g::G; L = bitlength(G)) where G <: ECGroup

    gxleaf = Leaf(value(gx(g)), div(L + 1, 8, RoundUp))
    gyleaf = Leaf(value(gy(g)), div(L + 1, 8, RoundUp))

    gtree = Tree((gxleaf, gyleaf))

    return gtree
end

function Tree(x::Vector{<:Group})
    L = bitlength(x[1])
    s = Tree[Tree(i, L = L) for i in x]
    return Node(s)
end

function marshal(x::PGroup)

    java_name = "com.verificatum.arithm.ModPGroup"
    p = modulus(x)
    q = order(x)
    g = Tree(x)
    e = UInt32(1)

    msg = (java_name, (p, q, g, e))

    tree = Tree(msg)

    return tree
end

normalize_ecgroup_name(x::String) = replace(x, "_"=>"-")
normalize_ecgroup_name(x::Symbol) = normalize_ecgroup_name(String(x))

function marshal(g::ECGroup)

    java_name = "com.verificatum.arithm.ECqPGroup"

    @assert spec(g) == spec(name(g)) "wrong group name"

    v_name = normalize_ecgroup_name(name(g))

    msg = (java_name, v_name)

    tree = Tree(msg)

    return tree
end

function unmarshal(tree::Tree)

    group_type = convert(String, tree.x[1])

    if group_type == "com.verificatum.arithm.ModPGroup"
        _unmarshal_pgroup(tree.x[2])
    elseif group_type == "com.verificatum.arithm.ECqPGroup"
        _unmarshal_ecgroup(tree.x[2])
    else
        error("Unrecognized group type: $group_type")
    end
end

function _unmarshal_pgroup(x::Node)

    (p, q, g, e) = convert(Tuple{BigInt, BigInt, BigInt, UInt32}, x)

    G = specialize(PGroup, p, q)

    x = G <| g

    return x
end

spec_name(x::String) = Symbol(replace(x, "-"=>"_"))

function _unmarshal_ecgroup(x::Leaf)

    group_spec_str = convert(String, x)
    name = spec_name(group_spec_str)

    group_spec = spec(name)
    G = specialize(ECGroup, group_spec; name)
    g = G <| generator(group_spec)

    return g
end

#function convert(::Type{ElGamal{G}}, tree::Tree) where G <: PGroup
#    𝐚, 𝐛 = convert(Tuple{Vector{BigInt}, Vector{BigInt}}, tree)

function convert(::Type{ElGamal{G}}, tree::Tree) where G <: Group
    𝐚, 𝐛 = convert(Tuple{Vector{G}, Vector{G}}, tree)
    𝐞 = ElGamal{G}(𝐚, 𝐛)
    return 𝐞
end

function Tree(𝐞::ElGamal{<:Group})
    𝐚 = a(𝐞)
    𝐛 = b(𝐞)
    tree = Tree((𝐚, 𝐛))
end


text:
using DatasetManager
using Documenter

DocMeta.setdocmeta!(DatasetManager, :DocTestSetup, :(using DatasetManager); recursive=true)
makedocs(;
    modules=[DatasetManager],
    authors="Allen Hill <allenofthehills@gmail.com> and contributors",
    repo="https://github.com/NantelBiomechLab/DatasetManager.jl/blob/{commit}{path}#L{line}",
    sitename="DatasetManager.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://NantelBiomechLab.github.io/DatasetManager.jl",
        assets=String[
            "assets/css/custom.css"
        ],
        highlights=["matlab"],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => [
            "Describing datasets" => "examples/datasets-examples.md",
            "Defining `AbstractSource` subtypes" => "examples/abstractsource-subtypes.md"
        ],
        "Julia Reference" => "julia-reference.md",
        "MATLAB Reference" => "matlab-reference.md"
    ],
)

deploydocs(;
    repo="github.com/NantelBiomechLab/DatasetManager.jl",
    push_preview=true
)


text:
module Epochs

using EarthOrientation
using ERFA

import Base: +, -, ==, isapprox, isless

using ..TimeScales, ..Periods
import ..TimeScales: acronyms
import ..LeapSeconds: leapseconds, LSK_DATA

export Epoch, julian, julian1, julian2, +, -, ==, isapprox, isless,
    leapseconds, jd2000, jd1950, mjd, timescale

const date_fmt = dateformat"yyyy-mm-ddTHH:MM:SS.sss"

struct Epoch{S,T<:Number}
    jd1::T
    jd2::T
    """
        Epoch{T}(jd1, jd2=0.0) where {T}

    Construct an `Epoch` with timescale `T` from a two-part Julian date.

    # Example

    ```jldoctest
    julia> Epoch{TT}(2.4578265e6, 0.30440190993249416)
    2017-03-14T07:18:20.325 TT
    ```
    """
    function Epoch{S}(jd1::T, jd2::T=zero(T)) where {S, T<:Number}
        new{S::TimeScale,T}(jd1, jd2)
    end
end

function Base.show(io::IO, ep::Epoch)
    print(io, "$(Dates.format(DateTime(ep), date_fmt)) $(timescale(ep))")
end

timescale(ep::Epoch{S}) where {S} = S

"""
    Epoch{T}(year, month, day,
        hour=0, minute=0, seconds=0, milliseconds=0) where {T}

Construct an `Epoch` with timescale `T` at the given date and time.

# Example

```jldoctest
julia> Epoch{TT}(2017, 3, 14, 7, 18, 20, 325)
2017-03-14T07:18:20.325 TT
```
"""
function Epoch{T}(year, month, day,
    hour=0, minute=0, seconds=0, milliseconds=0) where {T}
    jd, jd1 = ERFA.dtf2d(string(T),
    year, month, day, hour, minute, seconds + milliseconds/1000)
    Epoch{T}(jd, jd1)
end

"""
    Epoch{T}(dt::DateTime) where {T}

Convert a `DateTime` object to an `Epoch` with timescale `T`.

# Example

```jldoctest
julia> Epoch{TT}(DateTime(2017, 3, 14, 7, 18, 20, 325))
2017-03-14T07:18:20.325 TT
```
"""
function Epoch{T}(dt::DateTime) where {T}
    Epoch{T}(Dates.year(dt), Dates.month(dt), Dates.day(dt),
        Dates.hour(dt), Dates.minute(dt),
        Dates.second(dt) + Dates.millisecond(dt)/1000)
end

"""
    DateTime(ep::Epoch{T}) where T

Convert an `Epoch` with timescale `T` to a `DateTime` object.

# Example

```jldoctest
julia> DateTime(Epoch{TT}(2017, 3, 14, 7, 18, 20, 325))
2017-03-14T07:18:20.325
```
"""
function Base.DateTime(ep::Epoch{T}) where {T}
    dt = ERFA.d2dtf(string(T), 3, julian1(ep), julian2(ep))
    DateTime(dt...)
end

"""
    Epoch{T}(timestamp::AbstractString,
        fmt::DateFormat=dateformat"yyyy-mm-ddTHH:MM:SS.sss") where {T}

Construct an `Epoch` with timescale `T` from a timestamp. Optionally a `DateFormat`
object can be passed which improves performance if many date strings need to be
parsed and the format is known in advance.

# Example

```jldoctest
julia> Epoch{TT}("2017-03-14T07:18:20.325")
2017-03-14T07:18:20.325 TT
```
"""
Epoch{T}(str::AbstractString, fmt=date_fmt) where {T} = Epoch{T}(DateTime(str, fmt))

for scale in acronyms
    epoch = Symbol(scale, "Epoch")
    @eval begin
        const $epoch = Epoch{$scale}
        export $epoch
    end
end

julian1(ep) = ep.jd1
julian2(ep) = ep.jd2
julian(ep) = julian1(ep) + julian2(ep)
mjd(ep) = julian(ep) - MJD
jd2000(ep) = julian(ep) - J2000
jd1950(ep) = julian(ep) - J1950

(::Second)(ep::Epoch, base=0.0) = seconds((julian1(ep) - base + julian2(ep)) * days)
(::Minute)(ep::Epoch, base=0.0) = minutes((julian1(ep) - base + julian2(ep)) * days)
(::Hour)(ep::Epoch, base=0.0) = hours((julian1(ep) - base + julian2(ep)) * days)
(::Day)(ep::Epoch, base=0.0) = days((julian1(ep) - base + julian2(ep)) * days)
(::Year)(ep::Epoch, base=0.0) = years((julian1(ep) - base + julian2(ep)) * days)
(::Century)(ep::Epoch, base=0.0) = centuries((julian1(ep) - base + julian2(ep)) * days)

dut1(ep::Epoch) = getΔUT1(julian(ep))
leapseconds(ep::Epoch) = leapseconds(julian(ep))

function isapprox(a::Epoch{T}, b::Epoch{T}) where {T}
    return julian(a) ≈ julian(b)
end

function (==)(a::Epoch{T}, b::Epoch{T}) where {T}
    return DateTime(a) == DateTime(b)
end

isless(ep1::Epoch{T}, ep2::Epoch{T}) where {T} = julian(ep1) < julian(ep2)

function (+)(ep::Epoch{S,T1}, p::Period{U,T2}) where {T1,T2,S,U<:TimeUnit}
    delta = get(days(p))
    if delta >= oneunit(T2)
        ep1 = Epoch{S}(julian1(ep) + delta, julian2(ep))
    else
        ep1 = Epoch{S}(julian1(ep), julian2(ep) + delta)
    end
end

function (-)(ep::Epoch{S,T1}, p::Period{U,T2}) where {T1,T2,S,U<:TimeUnit}
    delta = get(days(p))
    if delta >= oneunit(T2)
        ep1 = Epoch{S}(julian1(ep) - delta, julian2(ep))
    else
        ep1 = Epoch{S}(julian1(ep), julian2(ep) - delta)
    end
end

function (-)(ep1::Epoch{T}, ep2::Epoch{T}) where {T}
    ((julian1(ep1) - julian1(ep2)) + (julian2(ep1) - julian2(ep2))) * days
end

include("conversions.jl")

end


text:
module reducedspin

using QuantumOptics, Base.Cartesian
using Combinatorics: combinations

export ReducedSpinBasis, reducedspintransition, reducedspinstate, reducedsigmap,
    reducedsigmam, reducedsigmax, reducedsigmay, reducedsigmaz, reducedsigmapsigmam

import Base: ==

using ..interaction, ..CollectiveSpins

"""
    ReducedSpinBasis(N, M)

Basis for a system of N spin 1/2 systems, up to the M'th excitation.
"""
struct ReducedSpinBasis{N, M, MS, T<:Int} <: Basis
    shape::Vector{T}
    indexMapper::Vector{Pair{Vector{T},T}}

    function ReducedSpinBasis(N::T, M::T, MS::T) where T<:Int
        if N < 1
            throw(DimensionMismatch())
        end
        if M < 1
            throw(DimensionMissmatch())
        end
        if M < MS
            throw(DimensionMissmatch())
        end

        dim = sum(binomial(N, k) for k=MS:M)
        inds = Vector{T}[]
        for k=MS:M
            append!(inds, collect(combinations(1:N,k)))
        end
        @assert length(inds)==dim
        indexMapper = (inds .=> [1:dim;])
        new{N, M, MS, T}([dim], indexMapper)
    end
end
ReducedSpinBasis(N::Int, M::Int) = ReducedSpinBasis(N, M, 0)
function Base.show(stream::IO, b::ReducedSpinBasis)
    write(stream, "ReducedSpin(N=$(b.N), M=$(b.M), MS=$(b.MS))")
end
function Base.getproperty(b::ReducedSpinBasis{N,M,MS}, s::Symbol) where {N,M,MS}
    if s === :N
        N
    elseif s === :M
        M
    elseif s === :MS
        MS
    else
        getfield(b, s)
    end
end

==(b1::T, b2::T) where T<:ReducedSpinBasis = true

"""
    index(b::ReducedSpinBasis, x:Vector{Int})

    Get the state index given excitation's positions.
"""
function index(b::ReducedSpinBasis, x::Vector{Int})
    @assert length(x) <= b.M
    @assert length(x) >= b.MS
    x_ = sort(x)
    i = findfirst(y->y[1]==x_, b.indexMapper)
    isa(i, Nothing) && throw(BoundsError(b, x_))
    return b.indexMapper[i][2]
end
index(b::ReducedSpinBasis, x::Int) = index(b, [x])
index(b::ReducedSpinBasis, x) = index(b, convert(Vector{Int}, x))

"""
    reducedspintransition(b::ReducedSpinBasis, to::Vector{Int}, from::Vector{Int})

Transition operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|``, where to and from are given as vectors denoting the excitations.
"""
function reducedspintransition(b::ReducedSpinBasis, to::Vector{Int}, from::Vector{Int})
    op = SparseOperator(b)
    op.data[index(b, to), index(b, from)] = 1.
    op
end
reducedspintransition(b::ReducedSpinBasis, to, from) = reducedspintransition(b, convert(Vector{Int}, to), convert(Vector{Int}, from))
reducedspintransition(b::ReducedSpinBasis, to::Int, from::Int) = reducedspintransition(b, [to], [from])

"""
    reducedsigmap(b::ReducedSpinBasis, j::Int)

Sigma Plus Operator for the j-th particle.
"""
function reducedsigmap(b::ReducedSpinBasis, j::Int)
    N = b.N
    M = b.M
    MS = b.MS
    @assert MS!=M

    op = SparseOperator(b)
    for m=MS+1:M
        to, from = transition_idx(b, j, m)
        for i=1:length(to)
            op.data[to[i],from[i]] = 1.0
        end
    end
    return op
end

"""
    transition_idx(b, j, m)

Find the indices of all basis states where a transition from the `m-1` excitation
manifold into the `m` excitation manifold can occur by raising atom `j`
"""
function transition_idx(b, j, m)
    # Find all states with correct number of excitations
    to_ = filter(x->length(x[1])==m, b.indexMapper)
    from_ = filter(x->length(x[1])==m-1, b.indexMapper)

    # to all states that involve j but only from states that do not involve j
    filter!(x->j∈x[1],to_)
    filter!(x->!(j∈x[1]),from_)

    # Get and return the corresponding indices
    to = [t[2] for t=to_]
    from = [f[2] for f=from_]
    return to, from
end

"""
    projector_idx(b, j, m)

Find the indices of all basis states in the `m` excitation manifold involving
the atom `j` in order to build a projector onto the upper state of this atom.
"""
function projector_idx(b,j,m)
    i_ = filter(x->length(x[1])==m, b.indexMapper)
    filter!(x->j∈x[1],i_)
    return getindex.(i_,2)
end

"""
    reducedsigmam(b::ReducedSpinBasis, j::Int)

Sigma Minus Operator for the j-th particle.
"""
function reducedsigmam(b::ReducedSpinBasis, j::Int)
    return dagger(reducedsigmap(b, j))
end

"""
    reducedsigmax(b::ReducedSpinBasis, j::Int)

Sigma-X Operator for the j-th particle.
"""
function reducedsigmax(b::ReducedSpinBasis, j::Int)
    return reducedsigmap(b, j) + reducedsigmam(b, j)
end

"""
    reducedsigmay(b::ReducedSpinBasis, j::Int)

Sigma-Y Operator for the j-th particle.
"""
function reducedsigmay(b::ReducedSpinBasis, j::Int)
    return im*(-reducedsigmap(b, j) + reducedsigmam(b, j))
end

"""
    reducedsigmaz(b::ReducedSpinBasis, j::Int)

Sigma-Z Operator for the j-th particle.
"""
function reducedsigmaz(b::ReducedSpinBasis, j::Int)
    proj = SparseOperator(b)
    for m=b.MS:b.M
        i = projector_idx(b,j,m)
        for i_=i
            proj.data[i_,i_] += 1.0
        end
    end
    return 2*proj - one(b)
end

"""
    reducedspinstate(b::ReducedSpinBasis, inds::Vector{Int})

State where the excitations are placed in the atoms given by `inds`. Note, that
`b.MS <= length(inds) <= b.M` must be satisfied.

# Examples
```julia-repl
julia> b = CollectiveSpins.ReducedSpinBasis(4,2)
ReducedSpin(N=4, M=2, MS=0)

julia> GS = CollectiveSpins.reducedspinstate(b,[]) # get the ground state
Ket(dim=11)
  basis: ReducedSpin(N=4, M=2, MS=0)
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

julia> ψ2 = CollectiveSpins.reducedspinstate(b,[1,2]) # First and second atom excited
Ket(dim=11)
  basis: ReducedSpin(N=4, M=2, MS=0)
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
```
"""
function reducedspinstate(b::ReducedSpinBasis, inds::Vector{Int})
    basisstate(b, index(b, inds))
end
reducedspinstate(b::ReducedSpinBasis, n) = reducedspinstate(b, convert(Vector{Int}, n))
reducedspinstate(b::ReducedSpinBasis, n::Int) = reducedspinstate(b, [n])

"""
    reducedsigmapsigmam(b::ReducedSpinBasis, i, j)

Create the operator product σᵢ⁺σⱼ⁻ directly on a [`ReducedSpinBasis`](@ref).
Useful especially for a basis where only one excitation is included, since then the
single operators are zero (do not conserve excitation number), but the product is not.
"""
function reducedsigmapsigmam(b::ReducedSpinBasis, i, j)
    op = SparseOperator(b)
    for idx1 in b.indexMapper
        if j in idx1[1]
            for idx2 in b.indexMapper
                i in idx2[1] && (op.data[idx2[2], idx1[2]] = 1)
            end
        end
    end
    return op
end

"""
    Hamiltonian(S::SpinCollection, M::Int=1)

Builds the dipole-dipole Hamiltonian.

# Arguments
* S: SpinCollection
* M: Number of excitations.
"""
function Hamiltonian(S::SpinCollection, M::Int=1)
    N = length(S.spins)
    b = ReducedSpinBasis(N, M)
    OmegaM = interaction.OmegaMatrix(S)

    H = SparseOperator(b)
    for i=1:N
        if S.spins[i].delta != 0.
            H += 0.5*S.spins[i].delta*reducedsigmaz(b,i)
        end
    end
    for i=1:N, j=1:N
        if i != j
            H += OmegaM[i, j]*reducedsigmapsigmam(b,i,j)
        end
    end

    return H
end

"""
    JumpOperators(S::SpinCollectino, M::Int=1)

Gamma Matix and Jump Operators for dipole-dipole interactino.

# Arguments
* S: Spin collection.
* M: Number of excitations.
"""
function JumpOperators(S::SpinCollection, M::Int=1)
        N = length(S.spins)
        b = ReducedSpinBasis(N, M)

        sm(j) = reducedsigmam(b, j)
        Jumps = [ sm(j) for j=1:N]
        GammaM = interaction.GammaMatrix(S)

        return GammaM, Jumps
    end

"""
    timeevolution(T, S::SpinCollection, psi0; kwargs...)
"""
function timeevolution(T, S::SpinCollection, psi0::Union{Ket{B}, DenseOpType{B, B}}; kwargs...) where B <: ReducedSpinBasis
    M = isa(psi0, Ket) ? psi0.basis.M : psi0.basis_l.M
    H = Hamiltonian(S, M)
    GammaM, J = JumpOperators(S, M)
    return QuantumOptics.timeevolution.master(T, psi0, H, J; rates=GammaM, kwargs...)
end

"""
    Hamiltonian_nh(S::SpinCollection,  M::Int=1)

Non-Hermitian Hamiltonian

# Argument
* S: SpinCollection
* M: Number of excitations.
"""
function Hamiltonian_nh(S::SpinCollection, M::Int=1, MS::Int=1)
    N = length(S.spins)
    b = ReducedSpinBasis(N, M, MS)
    OmegaM = interaction.OmegaMatrix(S)
    GammaM = interaction.GammaMatrix(S)

    H = SparseOperator(b)
    for i=1:N
        if S.spins[i].delta != 0.
            H += 0.5*S.spins[i].delta * reducedsigmaz(b,i)
        end
    end
    for i=1:N, j=1:N
            H += (OmegaM[i, j] - 0.5im*GammaM[i, j])*reducedsigmapsigmam(b,i,j)
    end

    return H
end

"""
    schroedinger_nh(T, S::SpinCollection, psi0; kwargs...)

    Time evolution of the non-Hermitian Hamiltonian in the single excitation manifold (M=1) via the Schrödinger equation.
"""
function schroedinger_nh(T, S::SpinCollection, psi0::Ket{B}; kwargs...) where B <: ReducedSpinBasis
    @assert psi0.basis.M == psi0.basis.MS == 1
    H = Hamiltonian_nh(S, 1, 1)

    return QuantumOptics.timeevolution.schroedinger(T, psi0, H; kwargs...)
end

end # module