text:
module TestUtils

using Base.Test
using DataFrames
using WorldBankDataTd
using Dates

## last day of month
##------------------

expOut = Date(2014, 9, 30)
actOut = WorldBankDataTd.getLastDayOfMonth(Date(2014, 9))
@test expOut == actOut

expOut = Date(2014, 3, 31)
actOut = WorldBankDataTd.getLastDayOfMonth(Date(2014, 3))
@test expOut == actOut

## last day of quarter
##--------------------

expOut = Date(2014, 3, 31)
actOut = WorldBankDataTd.getLastDayOfQuarter(Date(2014))
@test expOut == actOut

expOut = Date(2014, 9, 30)
actOut = WorldBankDataTd.getLastDayOfQuarter(Date(2014, 7))
@test expOut == actOut

expOut = Date(2014, 9, 30)
actOut = WorldBankDataTd.getLastDayOfQuarter(Date(2014, 9, 30))
@test expOut == actOut

## convert dates to Date type
##---------------------------

@test Date(2014, 12, 31) == WorldBankDataTd.formatDate("2014")
@test Date(2014, 3, 31) == WorldBankDataTd.formatDate("2014Q1")
@test Date(2014, 9, 30) == WorldBankDataTd.formatDate("2014Q3")
@test Date(2014, 3, 31) == WorldBankDataTd.formatDate("2014M03")
@test Date(2014, 9, 30) == WorldBankDataTd.formatDate("2014M09")

## conversion to float
##--------------------

@test isequal(WorldBankDataTd.tofloat("NA"), NA)
@test WorldBankDataTd.tofloat("1.23") == 1.23
@test isequal(WorldBankDataTd.tofloat(Void), NA)
@test isequal(WorldBankDataTd.tofloat(Nothing), NA)

## appending to DataArray
##-----------------------

### push Void

da = @data(["one", NA, "two"])
expDa = @data(["one", NA, "two", NA])

WorldBankDataTd.clean_and_push!(da, Void)
@test isequal(da, expDa)

### push ""

da = @data(["one", NA, "two"])
expDa = @data(["one", NA, "two", NA])

WorldBankDataTd.clean_and_push!(da, "")
@test isequal(da, expDa)

### push "NA"

da = @data(["one", NA, "two"])
expDa = @data(["one", NA, "two", NA])

WorldBankDataTd.clean_and_push!(da, "NA")
@test isequal(da, expDa)

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    MHProposalDistTuning

Abstract super-type for Metropolis-Hastings tuning strategies for
proposal distributions.
"""
abstract type MHProposalDistTuning <: MCMCTuningAlgorithm end
export MHProposalDistTuning

"""
    MetropolisHastings

Metropolis-Hastings MCMC sampling algorithm.

Constructors:

```julia
MetropolisHastings(
    proposal::ProposalDistSpec
    weighting::AbstractMCMCWeightingScheme
    tuning::TN = MHProposalDistTuning
)
"""
@with_kw struct MetropolisHastings{
    Q<:ProposalDistSpec,
    WS<:AbstractMCMCWeightingScheme,
    TN<:MHProposalDistTuning,
} <: MCMCAlgorithm
    proposal::Q = MvTDistProposal()
    weighting::WS = RepetitionWeighting()
    tuning::TN = AdaptiveMHTuning()
end

export MetropolisHastings

bat_default(::Type{MCMCSampling}, ::Val{:trafo}, sampler::MetropolisHastings) = PriorToGaussian()

bat_default(::Type{MCMCSampling}, ::Val{:nsteps}, sampler::MetropolisHastings, trafo::AbstractDensityTransformTarget, nchains::Integer) = 10^5

get_mcmc_tuning(algorithm::MetropolisHastings) = algorithm.tuning

mcmc_compatible(::MetropolisHastings, ::AbstractProposalDist, ::NoVarBounds) = true

mcmc_compatible(::MetropolisHastings, proposaldist::AbstractProposalDist, bounds::HyperRectBounds) =
    issymmetric(proposaldist) || all(x -> x == hard_bounds, bounds.bt)

mcmc_compatible(::MetropolisHastings, proposaldist::AbstractProposalDist, bounds::HierarchicalDensityBounds) =
    issymmetric(proposaldist)

mutable struct MHIterator{
    AL<:MetropolisHastings,
    D<:AbstractDensity,
    R<:AbstractRNG,
    PR<:RNGPartition,
    Q<:AbstractProposalDist,
    SV<:DensitySampleVector
} <: MCMCIterator
    algorithm::AL
    density::D
    rng::R
    rngpart_cycle::PR
    info::MCMCIteratorInfo
    proposaldist::Q
    samples::SV
    nsamples::Int64
    stepno::Int64
end

function MHIterator(
    rng::AbstractRNG,
    algorithm::MCMCAlgorithm,
    density::AbstractDensity,
    info::MCMCIteratorInfo,
    x_init::AbstractVector{P},
) where {P<:Real}
    stepno::Int64 = 0

    npar = totalndof(density)

    params_vec = Vector{P}(undef, npar)
    params_vec .= x_init
    !(params_vec in var_bounds(density)) && throw(ArgumentError("Parameter(s) out of bounds"))

    proposaldist = algorithm.proposal(P, npar)

    log_posterior_value = eval_logval(density, params_vec, strict = true)

    T = typeof(log_posterior_value)
    W = sample_weight_type(typeof(algorithm.weighting))

    sample_info = MCMCSampleID(info.id, info.cycle, 1, CURRENT_SAMPLE)
    current_sample = DensitySample(params_vec, log_posterior_value, one(W), sample_info, nothing)
    samples = DensitySampleVector{Vector{P},T,W,MCMCSampleID,Nothing}(undef, 0, npar)
    push!(samples, current_sample)

    nsamples::Int64 = 0

    rngpart_cycle = RNGPartition(rng, 0:(typemax(Int16) - 2))

    chain = MHIterator(
        algorithm,
        density,
        rng,
        rngpart_cycle,
        info,
        proposaldist,
        samples,
        nsamples,
        stepno
    )

    reset_rng_counters!(chain)

    chain
end

function MCMCIterator(
    rng::AbstractRNG,
    algorithm::MetropolisHastings,
    density::AbstractDensity,
    chainid::Integer,
    startpos::AbstractVector{<:Real}
)
    cycle = 0
    tuned = false
    converged = false
    info = MCMCIteratorInfo(chainid, cycle, tuned, converged)
    MHIterator(rng, algorithm, density, info, startpos)
end

@inline _current_sample_idx(chain::MHIterator) = firstindex(chain.samples)
@inline _proposed_sample_idx(chain::MHIterator) = lastindex(chain.samples)

getalgorithm(chain::MHIterator) = chain.algorithm

getdensity(chain::MHIterator) = chain.density

getrng(chain::MHIterator) = chain.rng

mcmc_info(chain::MHIterator) = chain.info

nsteps(chain::MHIterator) = chain.stepno

nsamples(chain::MHIterator) = chain.nsamples

current_sample(chain::MHIterator) = chain.samples[_current_sample_idx(chain)]

sample_type(chain::MHIterator) = eltype(chain.samples)

function reset_rng_counters!(chain::MHIterator)
    set_rng!(chain.rng, chain.rngpart_cycle, chain.info.cycle)
    rngpart_step = RNGPartition(chain.rng, 0:(typemax(Int32) - 2))
    set_rng!(chain.rng, rngpart_step, chain.stepno)
    nothing
end

function samples_available(chain::MHIterator)
    i = _current_sample_idx(chain::MHIterator)
    chain.samples.info.sampletype[i] == ACCEPTED_SAMPLE
end

function get_samples!(appendable, chain::MHIterator, nonzero_weights::Bool)::typeof(appendable)
    if samples_available(chain)
        samples = chain.samples

        for i in eachindex(samples)
            st = samples.info.sampletype[i]
            if (
                (st == ACCEPTED_SAMPLE || st == REJECTED_SAMPLE) &&
                (samples.weight[i] > 0 || !nonzero_weights)
            )
                push!(appendable, samples[i])
            end
        end
    end
    appendable
end

function next_cycle!(chain::MHIterator)
    _cleanup_samples(chain)

    chain.info = MCMCIteratorInfo(chain.info, cycle = chain.info.cycle + 1)
    chain.nsamples = 0
    chain.stepno = 0

    reset_rng_counters!(chain)

    resize!(chain.samples, 1)

    i = _proposed_sample_idx(chain)
    @assert chain.samples.info[i].sampletype == CURRENT_SAMPLE
    chain.samples.weight[i] = 1

    chain.samples.info[i] = MCMCSampleID(chain.info.id, chain.info.cycle, chain.stepno, CURRENT_SAMPLE)

    chain
end

function _cleanup_samples(chain::MHIterator)
    samples = chain.samples
    current = _current_sample_idx(chain)
    proposed = _proposed_sample_idx(chain)
    if (current != proposed) && samples.info.sampletype[proposed] == CURRENT_SAMPLE
        # Proposal was accepted in the last step
        @assert samples.info.sampletype[current] == ACCEPTED_SAMPLE
        samples.v[current] .= samples.v[proposed]
        samples.logd[current] = samples.logd[proposed]
        samples.weight[current] = samples.weight[proposed]
        samples.info[current] = samples.info[proposed]

        resize!(samples, 1)
    end
end

function mcmc_step!(chain::MHIterator)
    _cleanup_samples(chain)

    samples = chain.samples
    algorithm = getalgorithm(chain)

    if !mcmc_compatible(algorithm, chain.proposaldist, var_bounds(getdensity(chain)))
        error("Implementation of algorithm $algorithm does not support current parameter bounds with current proposal distribution")
    end

    chain.stepno += 1
    reset_rng_counters!(chain)

    rng = getrng(chain)
    density = getdensity(chain)

    proposaldist = chain.proposaldist

    # Grow samples vector by one:
    resize!(samples, size(samples, 1) + 1)
    samples.info[lastindex(samples)] = MCMCSampleID(chain.info.id, chain.info.cycle, chain.stepno, PROPOSED_SAMPLE)

    current = _current_sample_idx(chain)
    proposed = _proposed_sample_idx(chain)
    @assert current != proposed

    current_params = samples.v[current]
    proposed_params = samples.v[proposed]

    # Propose new variate:
    samples.weight[proposed] = 0
    proposal_rand!(rng, proposaldist, proposed_params, current_params)
    renormalize_variate!(proposed_params, density, proposed_params)

    current_log_posterior = samples.logd[current]
    T = typeof(current_log_posterior)

    # Evaluate prior and likelihood with proposed variate:
    proposed_log_posterior = eval_logval(density, proposed_params, strict = false)

    samples.logd[proposed] = proposed_log_posterior

    p_accept = if proposed_log_posterior > -Inf
        # log of ratio of forward/reverse transition probability
        log_tpr = if issymmetric(proposaldist)
            T(0)
        else
            log_tp_fwd = distribution_logpdf(proposaldist, proposed_params, current_params)
            log_tp_rev = distribution_logpdf(proposaldist, current_params, proposed_params)
            T(log_tp_fwd - log_tp_rev)
        end

        p_accept_unclamped = exp(proposed_log_posterior - current_log_posterior - log_tpr)
        T(clamp(p_accept_unclamped, 0, 1))
    else
        zero(T)
    end

    @assert p_accept >= 0
    accepted = rand(chain.rng, float(typeof(p_accept))) < p_accept

    if accepted
        samples.info.sampletype[current] = ACCEPTED_SAMPLE
        samples.info.sampletype[proposed] = CURRENT_SAMPLE
        chain.nsamples += 1
    else
        samples.info.sampletype[proposed] = REJECTED_SAMPLE
    end

    delta_w_current, w_proposed = _mh_weights(algorithm, p_accept, accepted)
    samples.weight[current] += delta_w_current
    samples.weight[proposed] = w_proposed

    nothing
end

function _mh_weights(
    algorithm::MetropolisHastings{Q,<:RepetitionWeighting},
    p_accept::Real,
    accepted::Bool
) where {Q,W}
    if accepted
        (0, 1)
    else
        (1, 0)
    end
end

function _mh_weights(
    algorithm::MetropolisHastings{Q,<:ARPWeighting},
    p_accept::Real,
    accepted::Bool
) where {Q,W}
    T = typeof(p_accept)
    if p_accept ≈ 1
        (zero(T), one(T))
    elseif p_accept ≈ 0
        (one(T), zero(T))
    else
        (T(1 - p_accept), p_accept)
    end
end

eff_acceptance_ratio(chain::MHIterator) = nsamples(chain) / nsteps(chain)

ext: jl


text:
using Oceananigans.Architectures
using Oceananigans.Architectures: architecture, arch_array, unified_array, device
using KernelAbstractions: @kernel, @index
using IterativeSolvers, SparseArrays, LinearAlgebra
using Oceananigans.Solvers:
                matrix_from_coefficients,
                update_diag!,
                arch_sparse_matrix,
                constructors,
                validate_settings,
                iterating,
                update_diag!

import Oceananigans.Grids: architecture
import Oceananigans.Solvers: solve!, iterate!

mutable struct UnifiedDiagonalIterativeSolver{G, Dev, R, L, D, M, P, F, T, B}
                       grid :: G
                    devices :: Dev
               problem_size :: R
        matrix_constructors :: L
                   diagonal :: D
                     matrix :: M
                       ρⁱ⁻¹ :: P
                   solution :: F
             matrix_product :: F
           search_direction :: F
                   residual :: F
                  tolerance :: T
                previous_Δt :: B
         maximum_iterations :: Int
                  iteration :: Int
                          n :: Int
end

function UnifiedDiagonalIterativeSolver(coeffs;
                                        grid,
                                        mrg,
                                        maximum_iterations = prod(size(grid)),
                                        tolerance = 1e-13,
                                        reduced_dim = (false, false, false),
                                        placeholder_timestep = -1.0)

    arch = architecture(mrg)
    temp_constructors, temp_diagonal, problem_size = matrix_from_coefficients(CPU(), grid, coeffs, reduced_dim)
    temp_matrix = arch_sparse_matrix(CPU(), temp_constructors)
    N = prod(problem_size)
    M = length(mrg)
    n = div(N, M)

    solution       = similar(unified_array(architecture(mrg), temp_diagonal))
    matrix_product = similar(unified_array(architecture(mrg), temp_diagonal))
    residual       = similar(unified_array(architecture(mrg), temp_diagonal))
    search         = similar(unified_array(architecture(mrg), temp_diagonal))

    placeholder_matrix  = []
    matrix_constructors = []
    diagonal            = []

    devices = mrg.devices

    for (idx, dev) in enumerate(devices)
        switch_device!(dev)
        push!(placeholder_matrix, temp_matrix[n*(idx-1)+1:n*idx, :])
        push!(matrix_constructors, constructors(arch, placeholder_matrix[idx]))
        push!(diagonal, arch_array(arch, temp_diagonal[n*(idx-1)+1:n*idx]))
        placeholder_matrix[idx] = arch_sparse_matrix(arch, placeholder_matrix[idx])
    end
    sync_all_devices!(devices)

    return UnifiedDiagonalIterativeSolver(mrg,
                                          devices,
                                          problem_size,
                                          matrix_constructors,
                                          diagonal,
                                          placeholder_matrix,
                                          zero(eltype(grid)),
                                          solution,
                                          matrix_product,
                                          search,
                                          residual,
                                          tolerance,
                                          placeholder_timestep,
                                          maximum_iterations,
                                          0,
                                          n)
end

architecture(solver::UnifiedDiagonalIterativeSolver) = architecture(solver.grid)

@inline function update_solver!(solver, Δt)
    arch = architecture(solver.grid)
    for (idx, dev) in enumerate(solver.devices)
        switch_device!(dev)
        constr = deepcopy(solver.matrix_constructors[idx])
        update_diag!(constr, arch, solver.n, prod(solver.problem_size), solver.diagonal[idx], Δt, solver.n * (idx - 1))
        solver.matrix[idx] = arch_sparse_matrix(arch, constr)
    end
    sync_all_devices!(solver.devices)
    solver.previous_Δt = Δt
end

function solve!(x, solver::UnifiedDiagonalIterativeSolver, b, Δt, args...)

    # arch = architecture(solver.grid)
    # prefetch_solver!(solver, arch)

    # Initialize
    solver.iteration = 0

    # q = A*x
    q = solver.matrix_product
    unified_mul!(q, solver, solver.solution)

    # r = b - A*x
    solver.residual .= b .- q

    while iterating(solver)
        iterate!(solver.solution, solver, args...)
    end

    return solver.solution
end

function iterate!(x, solver::UnifiedDiagonalIterativeSolver, args...)
    q = solver.matrix_product
    r = solver.residual
    p = solver.search_direction

    z = r
    ρ = dot(z, r)

    if solver.iteration == 0
        p .= z
    else
        β = ρ / solver.ρⁱ⁻¹
        p .= z .+ β .* p
    end

    # q = A*x
    unified_mul!(q, solver, p)
    α = ρ / dot(p, q)

    x .+= α .* p
    r .-= α .* q

    solver.iteration += 1
    solver.ρⁱ⁻¹ = ρ
    @show solver.iteration, solver.ρⁱ⁻¹
    return
end

@inline function unified_mul!(q, solver, x)
    for (idx, dev) in enumerate(solver.devices)
        switch_device!(dev)
        @views q[solver.n*(idx-1)+1:solver.n*idx] .= solver.matrix[idx] * x
    end
    sync_all_devices!(solver.devices)
end

@inline function unified_mul_and_dot!(q, solver, x, ρ)
    α = []
    for (idx, dev) in enumerate(solver.devices)
        switch_device!(dev)
        @views q[solver.n*(idx-1)+1:solver.n*idx] .= solver.matrix[idx] * x
        push!(α, dot(q[solver.n*(idx-1)+1:solver.n*idx], p[solver.n*(idx-1)+1:solver.n*idx]))
    end
    sync_all_devices!(solver.devices)
    return ρ / sum(α)
end

@inline prefetch!(array, bytes, dev::CuDevice) = Mem.prefetch(array.storage.buffer, bytes; device = dev)
@inline prefetch!(array, bytes, dev::CPU)      = nothing

@inline prefetch_solver!(solver, ::CPU) = nothing
@inline function prefetch_solver!(solver, ::GPU)
    bytes = sizeof(eltype(solver.grid)) * prod(solver.problem_size)
    for (idx, dev) in enumerate(solver.devices)
        switch_device!(dev)
        Mem.prefetch(solver.solution.storage.buffer,         bytes; device = dev)
        Mem.prefetch(solver.matrix_product.storage.buffer,   bytes; device = dev)
        Mem.prefetch(solver.residual.storage.buffer,         bytes; device = dev)
        Mem.prefetch(solver.search_direction.storage.buffer, bytes; device = dev)
    end
    for (idx, dev) in enumerate(solver.devices)
        sync_device!(dev)
    end
end

function Base.show(io::IO, solver::UnifiedDiagonalIterativeSolver)
    print(io, "Oceananigans-compatible preconditioned conjugate gradient solver.\n")
    print(io, " Problem size = "  , solver.problem_size, '\n')
    print(io, " Grid = "  , solver.grid, "\n")
    print(io, " Devices = ", solver.devices, "\n")
    print(io, " Divided length = ", solver.n)
    return nothing
end


text:
function SQ._collect(src::AbstractTable, g::SQ.SummarizeNode)
    res = default(src)
    for h in g.helpers
        apply!(res, h, src)
    end
    return res
end

function apply!(tbl, h::SQ.SummarizeHelper, src)::Void
    res_field, f, g, arg_fields = SQ.parts(h)
    T, row_itr = _preprocess(f, src, arg_fields)
    temp = Vector{T}()
    grow_nonnull_output!(temp, f, row_itr)
    tbl[res_field] = NullableArray([g(temp)])
    return
end

@noinline function grow_nonnull_output!(output, f, tpl_itr)
    for (i, tpl) in enumerate(tpl_itr)
        # Automatically lift the function f here.
        if !hasnulls(tpl)
            push!(output, f(map(unsafe_get, tpl)))
        end
    end
    return
end


text:
"""
    convert_to_struct(tree_dict::Dict{String,Any}, tree_obj::Node)

Takes a tree loaded from json and convert it back to a proper internal representation.

# Arguments
- `tree_dict::Dict{String,Any}`: tree as a dictionary after being read from json.
- `tree_obj::Node`: the tree (Node).
"""
function convert_to_struct(tree_dict::Dict{String,Any}, tree_obj::Node)

    # Iterate through each child of a node
    for (child_idx, child) in enumerate(tree_dict["children"])

        # Initialize child
        push!(tree_obj.children, Node([],label=child["taxa_label"]))

        # Add the CA's for the child
        for CA in child["CAs"]
            push!(tree_obj.children[child_idx].CAs, Rule((CA["idxs"]...,), ([ca[1] for ca in CA["char_attr"]]...,), CA["is_pure"], CA["num_group"], CA["num_non_group"], CA["occurances"]))
        end

        # Iterate through all of the child's descendents
        for (desc_idx, descendent) in enumerate(child["children"])

            # Initialize child's descendent
            push!(tree_obj.children[child_idx].children, Node([],label=descendent["taxa_label"]))

            # Add the CA's for the child's descendent
            for CA in descendent["CAs"]
                push!(tree_obj.children[child_idx].children[desc_idx].CAs, Rule((CA["idxs"]...,), ([ca[1] for ca in CA["char_attr"]]...,), CA["is_pure"], CA["num_group"], CA["num_non_group"], CA["occurances"]))
            end

            # Recur on the child's descendent to add the next layer to the tree
            convert_to_struct(descendent, tree_obj.children[child_idx].children[desc_idx])
        end
    end
    return tree_obj
end

"""
    get_next_hit(hitnames::Vector{String}, hit_idx::Int64)

Gets the next best hit returned from a blastn search.

# Arguments
- `hitnames::Vector{String}`: a list of all blastn hitnames.
- `hit_idx::Int64`: index of the current hit.
"""
function get_next_hit(hitnames::Vector{String}, hit_idx::Int64)

    # Initialize varibales
    new_hit = ""
    used_hits = Array{String,1}()
    new_idx = 0

    # Get all currently used hits
    for hit in hitnames[1:hit_idx]
        if !(hit in used_hits)
            push!(used_hits, hit)
        end
    end

    # Iterate through remaining hits to find the next new one
    for hit in hitnames[hit_idx+1:end]
        new_idx += 1
        if !(hit in used_hits)
            new_hit = hit
            break
        end
    end

    # Return the new hit and its index in the hitlist
    return new_hit,new_idx+hit_idx
end

"""
    add_blanks_to_back(subject::String, query::String, new_seq::String,
    subj_len::Int64, query_len::Int64, subj_non_blanks::Int64,
    hitnames::Vector{String}, hit_idx::Int64, character_labels::Dict{String,String})

Adds blanks to the back of a sequence from a blast match.

# Arguments
- `subject::String`: the subject the query is being matched to.
- `query::String`: the query that is having blanks added to it.
- `new_seq::String`: the new sequence (query with added blanks).
- `subj_len::Int64`: length of the subject.
- `query_len::Int64`: length of the query.
- `subj_non_blanks::Int64`: number of non blanks in the subject.
- `hitnames::Vector{String}`: list of blast hits.
- `hit_idx::Int64`: index of the current blast hit.
- `character_labels::Dict{String,String}`: a mapping of the character labels to the corresponding sequences.
"""
function add_blanks_to_back(subject::String, query::String, new_seq::String, subj_len::Int64, query_len::Int64, subj_non_blanks::Int64, hitnames::Vector{String}, hit_idx::Int64, character_labels::Dict{String,String})

    # If empty, return the sequence, and blanks to keep the same length as the subject
    if query == ""
        return new_seq * repeat("-", subj_len)

    # If empty, return the sequence
    elseif subj_len == 0
        return new_seq

    # If subject length equals query length, return the sequence with the rest of the query on the back
    elseif subj_len == query_len
        return new_seq * query

    # ******* We can probably remove this *******
    # If subject contains only blanks
    elseif subj_non_blanks == 0

        # Get the next best hit
        next_hit,hit_idx = get_next_hit(hitnames, hit_idx)

        # If no other hits, add rest of the subject (blanks) onto the sequence
        if next_hit == ""
            return new_seq * subject

        # Found a new hit
        else
            # Get the rest of the new subject
            new_subj = character_labels[next_hit][end-subj_len+1:end]

            # Finish the recursion on the new subject
            return add_blanks_to_back(new_subj, query, new_seq, subj_len, query_len, length(replace(new_subj, "-", "")), hitnames, hit_idx, character_labels)
        end
    end

    # Check if the current subject character is a blank
    if subject[1] == '-'

        # Recur on the rest of the subject and the entire query
        return add_blanks_to_back(subject[2:end], query, new_seq * '-', subj_len-1, query_len, subj_non_blanks, hitnames, hit_idx, character_labels)

    # Else recur on the rest of the subject and the rest of the query
    else
        return add_blanks_to_back(subject[2:end], query[2:end], new_seq * query[1], subj_len-1, query_len-1, subj_non_blanks-1, hitnames, hit_idx, character_labels)
    end
end

"""
    add_blanks_to_front(subject::String, query::String, new_seq::String,
    subj_len::Int64, query_len::Int64, subj_non_blanks::Int64,
    hitnames::Vector{String}, hit_idx::Int64, character_labels::Dict{String,String})

Adds blanks to the front of a sequence from a blast match.

# Arguments
- `subject::String`: the subject the query is being matched to.
- `query::String`: the query that is having blanks added to it.
- `new_seq::String`: the new sequence (query with added blanks).
- `subj_len::Int64`: length of the subject.
- `query_len::Int64`: length of the query.
- `subj_non_blanks::Int64`: number of non blanks in the subject.
- `hitnames::Vector{String}`: list of blast hits.
- `hit_idx::Int64`: index of the current blast hit.
- `character_labels::Dict{String,String}`: a mapping of the character labels to the corresponding sequences.
"""
function add_blanks_to_front(subject::String, query::String, new_seq::String, subj_len::Int64, query_len::Int64, subj_non_blanks::Int64, hitnames::Vector{String}, hit_idx::Int64, character_labels::Dict{String,String})

    # If empty, return the sequence, and blanks to keep the same length as the subject
    if query == ""
        return repeat("-", subj_len) * new_seq

    # If empty, return the sequence
    elseif subj_len == 0
        return new_seq

    # If subject length equals query length, return the sequence with the rest of the query on the front
    elseif subj_len == query_len
        return query * new_seq

    # ******* We can probably remove this *******
    # If subject contains only blanks
    elseif subj_non_blanks == 0

        # Get the next best hit
        next_hit,hit_idx = get_next_hit(hitnames, hit_idx)

        # If no other hits, add rest of the subject (blanks) onto the sequence
        if next_hit == ""
            return subject * new_seq

        # Found a new hit
        else
            # Get the rest of the new subject
            new_subj = character_labels[next_hit][1:subj_len]

            # Finish the recursion on the new subject
            return add_blanks_to_front(new_subj, query, new_seq, subj_len, query_len, length(replace(new_subj, "-", "")), hitnames, hit_idx, character_labels)
        end
    end

    # Check if the current subject character is a blank
    if subject[end] == '-'

        # Recur on the rest of the subject and the entire query
        return add_blanks_to_front(subject[1:end-1], query, '-' * new_seq, subj_len-1, query_len, subj_non_blanks, hitnames, hit_idx, character_labels)

    # Recur on the rest of the subject and the rest of the query
    else
        return add_blanks_to_front(subject[1:end-1], query[1:end-1], query[end] * new_seq, subj_len-1, query_len-1, subj_non_blanks-1, hitnames, hit_idx, character_labels)
    end
end

"""
    get_best_hit(results::Array{BioTools.BLAST.BLASTResult,1}, query::String,
    character_labels::Dict{String,String},  character_labels_no_gaps::Dict{String,String})

Gets the hit from blastn that has the most sequence coverage with no gaps compared to the query sequence.

# Arguments
- `results::Array{BioTools.BLAST.BLASTResult,1}`: blastn results.
- `query::String`: the query that is having blanks added to it.
- `character_labels::Dict{String,String}`: a mapping of the character labels to the corresponding sequences.
- `character_labels_no_gaps::Dict{String,String}`: character labels with gaps removed from sequences.
"""
function get_best_hit(results::Array{BioTools.BLAST.BLASTResult,1}, query::String, character_labels::Dict{String,String},  character_labels_no_gaps::Dict{String,String})

    # Initialize varibales
    error = Array{Int,1}(undef, 0)

    # Iterate through each result
    for result in results

        # Extract the subject from the result
        subject = character_labels[result.hitname]

        subject_no_gaps = character_labels_no_gaps[result.hitname]

        # Get the number of characters that would be lost if we used this results alignment
        try
            orig_subj_start = match(Regex(convert(String, result.hit)), subject_no_gaps).offset
            subj_start = get_adjusted_start(orig_subj_start, subject)
            query_start = match(Regex(convert(String, result.alignment.seq)), query).offset
            subject_front = subject[1:subj_start-1]
            query_front = query[1:query_start-1]
            subject_back = subject[subj_start:end]
            query_back = query[query_start:end]
            lost_characters = abs(length(replace(subject_front, "-" => "")) - length(query_front)) + abs(length(replace(subject_back, "-" => "")) - length(query_back))
            push!(error, lost_characters)
        catch
            push!(error, 1000000)
        end
    end

    # Return the reuslt which minimizes the number of characters lost from the new alignment
    return results[argmin(error)]
end

"""
    add_blanks(query_path::String, db_path::String, character_labels::Dict{String,String},
    character_labels_no_gaps::Dict{String,String} ; return_blast::Bool=false)

Adds blanks to an input sequence given a database.

# Arguments
- `query_path::String`: path to the query file.
- `db_path::String`: path to the blast database.
- `character_labels::Dict{String,String}`: a mapping of the character labels to the corresponding sequences.
- `character_labels_no_gaps::Dict{String,String}`: character labels with gaps removed from sequences.
- `return_blast::Bool=false`: whether to return blast results.
- `protein::Bool=false`: if protein sequence.
"""
function add_blanks(query_path::String, db_path::String, character_labels::Dict{String,String},
                    character_labels_no_gaps::Dict{String,String} ; return_blast::Bool=false,
                    protein::Bool=false)

    # Initialize
    best_result = 0
    orig_subj_start = 0

    # Read in the query sequence from file
    query = ""
    for line in readlines(query_path)
        if occursin(">", line)
            continue
        end
        query = query * line
    end

    # # Get the results from blastn
    if protein
         results = blastp(query_path, db_path, ["-task", "blastp", "-max_target_seqs", 10], db=true)
    else
        results = blastn(query_path, db_path, ["-task", "blastn", "-max_target_seqs", 10], db=true)
    end

    # Extract the best result
    try
        best_result = get_best_hit(results, query, character_labels, character_labels_no_gaps)
    catch
        if return_blast
            return "Error",["Error"]
        else
            return "Error"
        end
    end

    # Get all hitnames
    hitnames = [result.hitname for result in results]

    # Extract the subject as the top hit
    subject = character_labels[best_result.hitname]

    subject_no_gaps = character_labels_no_gaps[best_result.hitname]

    # Convert the best hit to a string and remove gaps that blast added
    best_hit = replace(convert(String, best_result.hit), "-" => "")

    # Get the starting position of the matched subject
    try
        orig_subj_start = match(Regex(best_hit), subject_no_gaps).offset
    catch
        if return_blast
            return "Error",["Error"]
        else
            return "Error"
        end
    end

    subj_start = get_adjusted_start(orig_subj_start, subject)

    # Convert the best hit to a string and remove gaps that blast added
    blast_alignment = replace(convert(String, best_result.alignment.seq), "-" => "")

    # Get the starting position of the query
    query_start = match(Regex(blast_alignment), query).offset

    # Get the fronts and backs of the subject and query
    subject_front = subject[1:subj_start-1]
    query_front = query[1:query_start-1]
    subject_back = subject[subj_start:end]
    query_back = query[query_start:end]

    # Return the query sequence with blanks added to the front and back
    new_seq = add_blanks_to_front(subject_front, query_front, "", subj_start-1, query_start-1, length(replace(subject_front, "-" => "")), hitnames, 1, character_labels) *
    add_blanks_to_back(subject_back, query_back, "", length(subject_back), length(query_back), length(replace(subject_back, "-" => "")), hitnames, 1, character_labels)

    if return_blast
        return new_seq, hitnames
    else
        return new_seq
    end
end


text:
# Adapted from Knet's src/data.jl (author: Deniz Yuret)

struct DataLoader{D,R<:AbstractRNG}
    data::D
    batchsize::Int
    nobs::Int
    partial::Bool
    shuffle::Bool
    rng::R
end

"""
    DataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG)

An object that iterates over mini-batches of `data`,
each mini-batch containing `batchsize` observations
(except possibly the last one).

Takes as input a single data tensor, or a tuple (or a named tuple) of tensors.
The last dimension in each tensor is the observation dimension, i.e. the one
divided into mini-batches.

If `shuffle=true`, it shuffles the observations each time iterations are re-started.
If `partial=false` and the number of observations is not divisible by the batchsize,
then the last mini-batch is dropped.

The original data is preserved in the `data` field of the DataLoader.

# Examples

```jldoctest
julia> Xtrain = rand(10, 100);

julia> array_loader = DataLoader(Xtrain, batchsize=2);

julia> for x in array_loader
         @assert size(x) == (10, 2)
         # do something with x, 50 times
       end

julia> array_loader.data === Xtrain
true

julia> tuple_loader = DataLoader((Xtrain,), batchsize=2);  # similar, but yielding 1-element tuples

julia> for x in tuple_loader
         @assert x isa Tuple{Matrix}
         @assert size(x[1]) == (10, 2)
       end

julia> Ytrain = rand('a':'z', 100);  # now make a DataLoader yielding 2-element named tuples

julia> train_loader = DataLoader((data=Xtrain, label=Ytrain), batchsize=5, shuffle=true);

julia> for epoch in 1:100
         for (x, y) in train_loader  # access via tuple destructuring
           @assert size(x) == (10, 5)
           @assert size(y) == (5,)
           # loss += f(x, y) # etc, runs 100 * 20 times
         end
       end

julia> first(train_loader).label isa Vector{Char}  # access via property name
true

julia> first(train_loader).label == Ytrain[1:5]  # because of shuffle=true
false

julia> foreach(println∘summary, DataLoader(rand(Int8, 10, 64), batchsize=30))  # partial=false would omit last
10×30 Matrix{Int8}
10×30 Matrix{Int8}
10×4 Matrix{Int8}
```
"""
function DataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG)
    batchsize > 0 || throw(ArgumentError("Need positive batchsize"))
    nobs = numobs(data)
    if nobs < batchsize
        @warn "Number of observations less than batchsize, decreasing the batchsize to $nobs"
        batchsize = nobs
    end
    DataLoader(data, batchsize, nobs, partial, shuffle, rng)
end

@propagate_inbounds function Base.iterate(d::DataLoader)
    data = d.data
    if d.shuffle
        data = shuffleobs(d.rng, data)
    end
    gen = eachobs(data; d.batchsize, d.partial)
    res = iterate(gen)
    res === nothing && return
    res[1], (gen, res[2])
end

@propagate_inbounds function Base.iterate(d::DataLoader, state)
    gen, i = state
    res = iterate(gen, i)
    res === nothing && return
    res[1], (gen, res[2])
end

function Base.length(d::DataLoader)
    n = d.nobs / d.batchsize
    d.partial ? ceil(Int, n) : floor(Int, n)
end

Base.IteratorEltype(d::DataLoader) = Base.EltypeUnknown()


text:
mutable struct XWindowManager <: AbstractWindowManager
    conn::Connection
    windows::Dict{xcb_window_t,XCBWindow}
    keymap::Keymap
    callbacks::Dict{XCBWindow, WindowCallbacks}
    function XWindowManager(conn::Connection, windows::Dict{xcb_window_t,XCBWindow}, keymap::Keymap, callbacks::Dict{XCBWindow, WindowCallbacks})
        wm = new(conn, windows, keymap, Dict())
        for (win, cb) ∈ callbacks
            set_callbacks!(wm, win, cb)
        end
        wm
    end
    XWindowManager(conn, windows, keymap, callbacks) = XWindowManager(convert(Connection, conn), convert(Dict{xcb_window_t,XCBWindow}, windows), convert(Keymap, keymap), convert(Dict{XCBWindow,WindowCallbacks}, callbacks))
end

XWindowManager(conn::Connection, windows::AbstractDict) = XWindowManager(conn, windows, Keymap(conn), Dict())
XWindowManager(conn::Connection, windows::Vector{XCBWindow}) = XWindowManager(conn, Dict(win.id => win for win in windows))

function set_callbacks!(wm::XWindowManager, win::XCBWindow, callbacks::WindowCallbacks)
    wm.callbacks[win] = callbacks
    set_event_mask(win, callbacks)
end

function poll_for_event(wm::XWindowManager)
    while true
        event = xcb_poll_for_event(wm.conn)
        event ≠ C_NULL && return event
        yield()
    end
end

function wait_for_event(wm::XWindowManager)
    event = xcb_wait_for_event(wm.conn)
    event == C_NULL ? nothing : event
end

function terminate_window!(wm::XWindowManager, win::XCBWindow)
    delete!(wm.windows, win.id)
    finalize(win)
end

get_window(wm::XWindowManager, id::Integer) = wm.windows[id]
get_window(wm::XWindowManager, event::xcb_xkb_state_notify_event_t) = nothing
get_window(wm::XWindowManager, event::xcb_keymap_notify_event_t) = nothing
get_window(wm::XWindowManager, event) = get_window(wm, window_id(event))

callbacks(wm::XWindowManager, win::XCBWindow) = get(wm.callbacks, win, WindowCallbacks())

"""
Extract a keycode from an key event.
"""
keycode(wm::XWindowManager, details::EventDetails{<:KeyEvent}) = keycode(wm.keymap, details.data.key_name)

"""
Return condensed information regarding a keystroke as a `String`. Includes key name, keycode, character input and key symbol.
"""
keystroke_info(wm::XWindowManager, event_details::EventDetails) = keystroke_info(wm.keymap, event_details)


text:
# Define and register smooth functions
_cos_wave(t, f, A, st, ϕ) = A*cos(2*π*f*(t - st) + ϕ)
_sin_wave(t, f, A, st, ϕ) = A*sin(2*π*f*(t - st) + ϕ)
_damped_sine_wave(t, f, A, st, ϕ, d) = exp((st-t)*d)*A*sin(2*π*f*(t-st) + ϕ)
_ramp(t, δ, st, et, h) = h/(et-st)*(_xH(t, δ, st) - _xH(t, δ, et))
_square_wave(t, δ, f, A, st) = A*2atan(sin(2π*(t-st)*f)/δ)/π
_step(t, δ, h, a) = h*(atan((t-a)/δ)/π + 0.5)
_triangular_wave(t, δ, f, A, st) = A*(1-2acos((1 - δ)sin(2π*(t-st)*f))/π)
_xH(t, δ, tₒ) = (t-tₒ)*(1+((t-tₒ)/sqrt((t-tₒ)^2+δ^2)))/2

@register_symbolic _cos_wave(t, f, A, st, ϕ)
@register_symbolic _sin_wave(t, f, A, st, ϕ)
@register_symbolic _damped_sine_wave(t, f, A, st, ϕ, damping)
@register_symbolic _ramp(t, δ, st, et, h)
@register_symbolic _square_wave(t, δ, f, A, st)
@register_symbolic _step(t, δ, h, a)
@register_symbolic _triangular_wave(t, δ, f, A, st)

"""
```julia
function ConstantVoltage(; name, V)
```

The source for an ideal constant voltage.

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v` and is always constant

# Connectors
- `p`
  Positive pin
- `n`
  Negative pin

# Parameters:
- `V`: [`V`]
  The constant voltage across the terminals of this source
"""
function ConstantVoltage(;name, V)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters V=V
    eqs = [
        v ~ V
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function CosineVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0)
```

A source in which the voltage across its terminals is a cosine function of time.

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Observables
- `offset`: [`V`]
  A constant offset added to the voltage output
- `amplitude`: [`V`]
  The amplitude of the cosine function
- `frequency`: [`Hz`]
  The frequency of the cosine function
- `starttime`: [`s`]
  The time at which the source starts functioning. Before this time, the voltage across
  its terminals is 0.
- `phase`: [`rad`]
  The phase offset of the cosine function
"""
function CosineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0)
    δ = 0.00001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
        phase=phase
    end
    eqs = [
        v ~ _cos_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function ExpSineVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0)
```

A source in which the voltage across its terminals is a damped sine function of time.

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`V`]
  A constant offset added to the voltage output
- `amplitude`: [`V`]
  The amplitude of the damped sine function
- `frequency`: [`Hz`]
  The frequency of the damped sine function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this time, the voltage across
  its terminals is `offset`.
- `phase`: [`rad`]
  The phase offset of the damped sine function
- `damping_coef`: [`1/s`]
  Damping coefficient of the damped sine function
"""
function ExpSineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0)
    δ = 0.00001
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
        phase=phase
        damping=damping
    end
    eqs = [
        v ~ _damped_sine_wave(t, frequency, amplitude, start_time, phase, damping) * _step(t, δ, 1.0, start_time)
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function RampVoltage(; name, offset=0.0, start_time=0.0, duration=1.0, height=1.0)
```

A source in which the voltage across grows linearly from `offset` to `offset+height` over
the time interval `duration` starting at `start_time`

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`V`]
  A constant offset added to the voltage output
- `start_time`: [`s`]
  The time at which the voltage starts growing
- `duration`: [`s`]
  The duration of the ramp (`0.0` gives a step)
- `height`: [`V`]
  The amount that the voltage grows in the time interval
"""
function RampVoltage(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0)
    δ = 0.00001
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        height=height
        start_time=start_time
        duration=duration
    end
    eqs = [
        v ~ _ramp(t, δ, start_time, start_time + duration, height) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function SineVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0)
```

A source in which the voltage across its terminals is a sine function of time.

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`V`]
  A constant offset added to the voltage output
- `amplitude`: [`V`]
  The amplitude of the sine function
- `frequency`: [`Hz`]
  The frequency of the sine function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this time, the voltage across
  its terminals is `offset`.
- `phase`: [`rad`]
  The phase offset of the sine function
"""
function SineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0)
    δ = 0.00001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
        phase=phase
    end
    eqs = [
        v ~ _sin_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function SquareVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
```

A source in which the voltage across its terminals is a square function of time.

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`V`]
  A constant offset added to the voltage output
- `amplitude`: [`V`]
  The amplitude of the square wave function
- `frequency`: [`Hz`]
  The frequency of the square wave function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this time, the voltage across
  its terminals is `offset`.
"""
function SquareVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
    δ = 0.0001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        start_time=start_time
    end
    eqs = [
        v ~ _square_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function StepVoltage(; name, offset=0.0, start_time=0.0, height=1.0)
```

A source in which the voltage across its terminals increases from `offset` to `offset+height` at
`start_time`

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Observables
- `offset`: [`V`]
  A constant offset added to the voltage output
- `start_time`: [`s`]
  The time at which the source starts functioning, and the voltage jumps
- `height`: [`V`]
  Magnitude of increase in voltage
"""
function StepVoltage(;name, offset=0.0, start_time=0.0, height=1.0)
    δ = 0.0001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        height=height
        start_time=start_time
    end
    eqs = [
        v ~ _step(t, δ, height, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function TriangularVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
```

A source in which the voltage across its terminals is a triangular function of time.

# States
- `v(t)`: [`V`]
  The voltage across this source, given by `p.v - n.v`

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Observables
- `offset`: [`V`]
  A constant offset added to the voltage output
- `amplitude`: [`V`]
  Amplitude of the triangular wave function
- `frequency`: [`Hz`]
  Frequency of the triangular wave function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this, the output of the source is
  `offset`
"""
function TriangularVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
    δ = 0.00001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
    end
    eqs = [
        v ~ _triangular_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

# Current Sources ######################################################################################################
"""
```julia
function ConstantCurrent(; name, I = 1.0)
```

The source for an ideal constant current.

# States
- `i(t)`: [`A`]
  The current through this source, which is always constant

# Connectors
- `p`
  Positive pin
- `n`
  Negative pin

# Parameters:
- `I`: [`A`]
  The constant current through the terminals of this source
"""
function ConstantCurrent(;name, I = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters I=I
    eqs = [
        i ~ I
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function CosineCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0)
```

A source in which the current through its terminals is a cosine function of time.

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Observables
- `offset`: [`A`]
  A constant offset added to the current output
- `amplitude`: [`A`]
  The amplitude of the cosine function
- `frequency`: [`Hz`]
  The frequency of the cosine function
- `starttime`: [`s`]
  The time at which the source starts functioning. Before this time, the current through
  its terminals is 0.
- `phase`: [`rad`]
  The phase offset of the cosine function
"""
function CosineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0)
    δ = 0.00001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
        phase=phase
    end
    eqs = [
        i ~ _cos_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function ExpSineCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0)
```

A source in which the current through its terminals is a damped sine function of time.

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`A`]
  A constant offset added to the current output
- `amplitude`: [`A`]
  The amplitude of the damped sine function
- `frequency`: [`Hz`]
  The frequency of the damped sine function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this time, the current through
  its terminals is `offset`.
- `phase`: [`rad`]
  The phase offset of the damped sine function
- `damping_coef`: [`1/s`]
  Damping coefficient of the damped sine function
"""
function ExpSineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0)
    δ = 0.00001
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
        phase=phase
        damping=damping
    end
    eqs = [
        i ~ _damped_sine_wave(t, frequency, amplitude, start_time, phase, damping) * _step(t, δ, 1.0, start_time)
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function RampCurrent(; name, offset=0.0, start_time=0.0, duration=1.0, height=1.0)
```

A source in which the current grows linearly from `offset` to `offset+height` over
the time interval `duration` starting at `start_time`

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`A`]
  A constant offset added to the current output
- `start_time`: [`s`]
  The time at which the current starts growing
- `duration`: [`s`]
  The duration of the ramp (`0.0` gives a step)
- `height`: [`A`]
  The amount that the current grows in the time interval
"""
function RampCurrent(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0)
    δ = 0.00001
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        height=height
        start_time=start_time
        duration=duration
    end
    eqs = [
        i ~ _ramp(t, δ, start_time, start_time + duration, height) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function SineCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0)
```

A source in which the current through its terminals is a sine function of time.

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`A`]
  A constant offset added to the current output
- `amplitude`: [`V`]
  The amplitude of the sine function
- `frequency`: [`Hz`]
  The frequency of the sine function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this time, the current through
  its terminals is `offset`.
- `phase`: [`rad`]
  The phase offset of the sine function
"""
function SineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0)
    δ = 0.00001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
        phase=phase
    end
    eqs = [
        i ~ _sin_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function SquareCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
```

A source in which the current through its terminals is a square function of time.

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Parameters
- `offset`: [`A`]
  A constant offset added to the current output
- `amplitude`: [`A`]
  The amplitude of the square wave function
- `frequency`: [`Hz`]
  The frequency of the square wave function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this time, the current through
  its terminals is `offset`.
"""
function SquareCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
    δ = 0.0001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        start_time=start_time
    end
    eqs = [
        i ~ _square_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function StepCurrent(; name, offset=0.0, start_time=0.0, height=1.0)
```

A source in which the current through its terminals increases from `offset` to `offset+height` at
`start_time`

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Observables
- `offset`: [`A`]
  A constant offset added to the current output
- `start_time`: [`s`]
  The time at which the source starts functioning, and the current jumps
- `height`: [`A`]
  Magnitude of increase in current
"""
function StepCurrent(;name, offset=0.0, start_time=0.0, height=1.0)
    δ = 0.0001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        height=height
        start_time=start_time
    end
    eqs = [
        i ~ _step(t, δ, height, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
```julia
function TriangularCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
```

A source in which the current through its terminals is a triangular function of time.

# States
- `i(t)`: [`A`]
  The current through this source

# Connectors
- `p`
  Positive port
- `n`
  Negative port

# Observables
- `offset`: [`A`]
  A constant offset added to the current output
- `amplitude`: [`A`]
  Amplitude of the triangular wave function
- `frequency`: [`Hz`]
  Frequency of the triangular wave function
- `start_time`: [`s`]
  The time at which the source starts functioning. Before this, the output of the source is
  `offset`
"""
function TriangularCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0)
    δ = 0.00001

    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters begin
        offset=offset
        amplitude=amplitude
        frequency=frequency
        start_time=start_time
    end
    eqs = [
        i ~ _triangular_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset
    ]

    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end


text:
"""
BDD Module
"""

module BDD

import Base
export bdd, header!, var!, node!, not, and, or, xor, imp, ite, todot

"""
AbstractNode{Ts}

An abstract node for BDD. The parameter Ts represents a type for a symbol of variable.
Ts is usually Symbol.
"""

abstract type AbstractNode{Ts} end

function Base.show(io::IO, n::AbstractNode{Ts}) where Ts
    Base.show(io, "node$(n.id)")
end

"""
type alias
"""

const HeaderID = UInt
const NodeID = UInt

"""
struct
"""

mutable struct NodeHeader{Ts}
    id::HeaderID
    level::Int
    label::Ts
end

mutable struct Node{Ts} <: AbstractNode{Ts}
    id::NodeID
    header::NodeHeader{Ts}
    low::AbstractNode{Ts}
    high::AbstractNode{Ts}
end

struct Terminal{Ts} <: AbstractNode{Ts}
    id::NodeID
end

"""
Operators
"""

abstract type AbstractOperator end

struct BDDNot <: AbstractOperator end
struct BDDAnd <: AbstractOperator end
struct BDDOr <: AbstractOperator end
struct BDDXor <: AbstractOperator end

struct UniOperator{Ts}
    op::AbstractOperator
    cache::Dict{NodeID,AbstractNode{Ts}}
end

struct BinOperator{Ts}
    op::AbstractOperator
    cache::Dict{Tuple{NodeID,NodeID},AbstractNode{Ts}}
end

"""
BDDForest

A structure to store the information on BDD.

Todo: Get rid of `mutable` (export `totalnodeid` to BDDNodeManager?)
"""

mutable struct BDDForest{Ts}
    totalnodeid::NodeID
    vars::Vector{NodeHeader{Ts}}
    headers::Dict{Ts,NodeHeader{Ts}}
    utable::Dict{Tuple{HeaderID,NodeID,NodeID},AbstractNode{Ts}}
    zero::Terminal{Ts}
    one::Terminal{Ts}
    notop::UniOperator{Ts}
    andop::BinOperator{Ts}
    orop::BinOperator{Ts}
    xorop::BinOperator{Ts}
end

function bdd(::Type{Ts} = Symbol) where Ts
    zero = Terminal{Ts}(0)
    one = Terminal{Ts}(1)
    start_node_id::NodeID = 2
    v = Vector{NodeHeader{Ts}}[]
    h = Dict{Ts,NodeHeader{Ts}}()
    ut = Dict{Tuple{HeaderID,NodeID,NodeID},AbstractNode{Ts}}()
    notop = UniOperator{Ts}(BDDNot(), Dict{NodeID,AbstractNode{Ts}}())
    andop = BinOperator{Ts}(BDDAnd(), Dict{Tuple{NodeID,NodeID},AbstractNode{Ts}}())
    orop = BinOperator{Ts}(BDDOr(), Dict{Tuple{NodeID,NodeID},AbstractNode{Ts}}())
    xorop = BinOperator{Ts}(BDDXor(), Dict{Tuple{NodeID,NodeID},AbstractNode{Ts}}())
    b = BDDForest{Ts}(2, v, h, ut, zero, one, notop, andop, orop, xorop)
end

function get_next_nodeid!(b::BDDForest{Ts}) where Ts
    id = b.totalnodeid
    b.totalnodeid += 1
    id
end

"""
header!

Get a header for a variable on BDD with a label
"""

function header!(b::BDDForest{Ts}, label::Ts)::NodeHeader{Ts} where Ts
    h = get(b.headers, label) do
        h = NodeHeader{Ts}(HeaderID(length(b.vars)), length(b.vars), label)
        push!(b.vars, h)
        b.headers[label] = h
    end
    h
end

"""
var!

Get a node for a variable
"""

function var!(b::BDDForest{Ts}, label::Ts)::AbstractNode{Ts} where Ts
    h = header!(b, label)
    node!(b, h, b.zero, b.one)
end

"""
node!

Create a BDD node
"""

function node!(b::BDDForest{Ts}, h::NodeHeader{Ts}, low::AbstractNode{Ts}, high::AbstractNode{Ts})::AbstractNode{Ts} where Ts
    if low.id == high.id
        return low
    end
    key = (h.id, low.id, high.id)
    get(b.utable, key) do
        nextid = get_next_nodeid!(b)
        b.utable[key] = Node(nextid, h, low, high)
    end
end

"""
not

Return a node of negation for a given node.
"""

function not(b::BDDForest{Ts}, f::AbstractNode{Ts}) where Ts
    uniapply!(b.notop, b, f)
end

"""
and

Return a node of AND for given nodes
"""

function and(b::BDDForest{Ts}, f::Vararg{AbstractNode{Ts}}) where Ts
    ans = b.one
    for i = 1:length(f)
        ans = binapply!(b.andop, b, ans, f[i])
    end
    ans
end

"""
or

Return a node of OR for given nodes
"""

function or(b::BDDForest{Ts}, f::Vararg{AbstractNode{Ts}}) where Ts
    ans = b.zero
    for i = 1:length(f)
        ans = binapply!(b.orop, b, ans, f[i])
    end
    ans
end

"""
xor

Return a node of XOR for given two nodes
"""

function xor(b::BDDForest{Ts}, f::AbstractNode{Ts}, g::AbstractNode{Ts}) where Ts
    binapply!(b.xorop, b, f, g)
end

"""
imp

Return a node of IMP for given two nodes
"""

function imp(b::BDDForest{Ts}, f::AbstractNode{Ts}, g::AbstractNode{Ts}) where Ts
    or(b, not(b, f), g)
end

"""
imp

Return a node of if-then-else for given three nodes
"""

function ite(b::BDDForest{Ts}, f::AbstractNode{Ts}, g::AbstractNode{Ts}, h::AbstractNode{Ts}) where Ts
    or(b, and(b, f, g), and(b, not(b, f), h))
end

"""
uniapply

Apply operation for unioperator
"""

function uniapply!(op::UniOperator{Ts}, b::BDDForest{Ts}, f::AbstractNode{Ts})::AbstractNode{Ts} where Ts
    return _uniapply!(op.op, op, b, f)
end

function _uniapply!(::AbstractOperator, op::UniOperator{Ts}, b::BDDForest{Ts}, f::Node{Ts})::AbstractNode{Ts} where Ts
    get(op.cache, f.id) do
        n0 = _uniapply!(op.op, op, b, f.low)
        n1 = _uniapply!(op.op, op, b, f.high)
        ans = node!(b, f.header, n0, n1)
        op.cache[f.id] = ans
    end
end

function _uniapply!(::BDDNot, op::UniOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts})::AbstractNode{Ts} where Ts
    if f == b.one
        b.zero
    else
        b.one
    end
end

"""
binapply

Apply operation for binoperator
"""

function binapply!(op::BinOperator{Ts}, b::BDDForest{Ts}, f::AbstractNode{Ts}, g::AbstractNode{Ts})::AbstractNode{Ts} where Ts
    return _binapply!(op.op, op, b, f, g)
end

function _binapply!(::AbstractOperator, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Node{Ts}, g::Node{Ts})::AbstractNode{Ts} where Ts
    key = (f.id, g.id)
    get(op.cache, key) do
        if f.header.level > g.header.level
            n0 = _binapply!(op.op, op, b, f.low, g)
            n1 = _binapply!(op.op, op, b, f.high, g)
            ans = node!(b, f.header, n0, n1)
        elseif f.header.level < g.header.level
            n0 = _binapply!(op.op, op, b, f, g.low)
            n1 = _binapply!(op.op, op, b, f, g.high)
            ans = node!(b, g.header, n0, n1)
        else
            n0 = _binapply!(op.op, op, b, f.low, g.low)
            n1 = _binapply!(op.op, op, b, f.high, g.high)
            ans = node!(b, f.header, n0, n1)
        end
        op.cache[key] = ans
    end
end

function _binapply!(::AbstractOperator, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts}, g::Node{Ts})::AbstractNode{Ts} where Ts
    key = (f.id, g.id)
    get(op.cache, key) do
        n0 = _binapply!(op.op, op, b, f, g.low)
        n1 = _binapply!(op.op, op, b, f, g.high)
        op.cache[key] = node!(b, g.header, n0, n1)
    end
end

function _binapply!(::AbstractOperator, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Node{Ts}, g::Terminal{Ts})::AbstractNode{Ts} where Ts
    key = (f.id, g.id)
    get(op.cache, key) do
        n0 = _binapply!(op.op, op, b, f.low, g)
        n1 = _binapply!(op.op, op, b, f.high, g)
        op.cache[key] = node!(b, f.header, n0, n1)
    end
end

## and

function _binapply!(::BDDAnd, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts}, g::Node{Ts})::AbstractNode{Ts} where Ts
    if f == b.one
        g
    else
        b.zero
    end
end

function _binapply!(::BDDAnd, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Node{Ts}, g::Terminal{Ts})::AbstractNode{Ts} where Ts
    if g == b.one
        f
    else
        b.zero
    end
end

function _binapply!(::BDDAnd, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts}, g::Terminal{Ts})::AbstractNode{Ts} where Ts
    if f == b.one && g == b.one
        b.one
    else
        b.zero
    end
end

## or

function _binapply!(::BDDOr, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts}, g::Node{Ts})::AbstractNode{Ts} where Ts
    if f == b.one
        b.one
    else
        g
    end
end

function _binapply!(::BDDOr, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Node{Ts}, g::Terminal{Ts})::AbstractNode{Ts} where Ts
    if g == b.one
        b.one
    else
        f
    end
end

function _binapply!(::BDDOr, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts}, g::Terminal{Ts})::AbstractNode{Ts} where Ts
    if f == b.zero && g == b.zero
        b.zero
    else
        b.one
    end
end

## xor

function _binapply!(::BDDXor, op::BinOperator{Ts}, b::BDDForest{Ts}, f::Terminal{Ts}, g::Terminal{Ts})::AbstractNode{Ts} where Ts
    if f == g
        b.zero
    else
        b.one
    end
end

"""
todot(forest, f)
Return a string for dot to draw a diagram.
"""

function todot(b::BDDForest{Ts}, f::AbstractNode{Ts}) where Ts
    io = IOBuffer()
    visited = Set{AbstractNode{Ts}}()
    println(io, "digraph { layout=dot; overlap=false; splines=true; node [fontsize=10];")
    _todot!(b, f, visited, io)
    println(io, "}")
    return String(take!(io))
end

function _todot!(b::BDDForest{Ts}, f::Terminal{Ts}, visited::Set{AbstractNode{Ts}}, io::IO)::Nothing where Ts
    if in(f, visited)
        return
    end
    if f == b.zero
        println(io, "\"obj$(f.id)\" [shape = square, label = \"0\"];")
    else
        println(io, "\"obj$(f.id)\" [shape = square, label = \"1\"];")
    end
    push!(visited, f)
    nothing
end

function _todot!(b::BDDForest{Ts}, f::Node{Ts}, visited::Set{AbstractNode{Ts}}, io::IO)::Nothing where Ts
    if in(f, visited)
        return
    end
    println(io, "\"obj$(f.id)\" [shape = circle, label = \"$(f.header.label)\"];")
    _todot!(b, f.low, visited, io)
    _todot!(b, f.high, visited, io)
    println(io, "\"obj$(f.id)\" -> \"obj$(f.low.id)\" [label = \"0\"];")
    println(io, "\"obj$(f.id)\" -> \"obj$(f.high.id)\" [label = \"1\"];")
    push!(visited, f)
    nothing
end

end


text:
using Measurements

import Dates
import ERFA

@testset "Epochs" begin
    @testset "UTC" begin
        before = "2012-06-30T23:59:59.000"
        start = "2012-06-30T23:59:60.000"
        during = "2012-06-30T23:59:60.300"
        after = "2012-07-01T00:00:00.000"

        before_exp = (second=394372833, fraction=0.0)
        start_exp = (second=394372834, fraction=0.0)
        during_exp = (second=394372834, fraction=0.3)
        after_exp = (second=394372835, fraction=0.0)

        before_act = from_utc(before)
        start_act = from_utc(start)
        during_act = from_utc(during)
        after_act = from_utc(after)

        @test before_act.second == before_exp.second
        @test before_act.fraction ≈ before_exp.fraction
        @test start_act.second == start_exp.second
        @test start_act.fraction ≈ start_exp.fraction
        @test during_act.second == during_exp.second
        @test during_act.fraction ≈ during_exp.fraction
        @test after_act.second == after_exp.second
        @test after_act.fraction ≈ after_exp.fraction

        during_dt = from_utc(2012, 6, 30, 23, 59, 60.3)
        during_dt1 = from_utc(2012, 6, 30, 23, 59, 60, 0.3)
        @test during_dt.second == during_exp.second
        @test during_dt.fraction ≈ during_exp.fraction
        @test during_dt1.second == during_exp.second
        @test during_dt1.fraction ≈ during_exp.fraction

        @test to_utc(before_act) == before
        @test to_utc(start_act) == start
        @test to_utc(during_act) == during
        @test to_utc(after_act) == after

        dt = Dates.DateTime(2016, 12, 31, 23, 59, 59, 300)
        @test dt == to_utc(Dates.DateTime, from_utc(dt))

        sixties = AstroTime.DateTime(1961, 3, 5, 23, 4, 12.0)
        sixties_exp = (second=-1225198547, fraction=0.5057117799999999)
        sixties_act = from_utc(sixties)
        sixties_utc = to_utc(AstroTime.DateTime, sixties_act)

        @test sixties_act.second == sixties_exp.second
        @test sixties_act.fraction ≈ sixties_exp.fraction
        @test sixties_utc ≈ sixties
    end
    @testset "Precision" begin
        ep = TAIEpoch(TAIEpoch(2000, 1, 1, 12), 2eps())
        @test ep.second == 0
        @test ep.fraction == 2eps()

        ep += 10000centuries
        @test ep.second == value(seconds(10000centuries))
        @test ep.fraction == 2eps()

        # Issue 44
        elong1 = 0.0
        elong2 = π
        u = 6371.0
        tt = TTEpoch(2000, 1, 1)
        tdb_tt1 = getoffset(TT, TDB, tt.second, tt.fraction, elong1, u, 0.0)
        tdb_tt2 = getoffset(TT, TDB, tt.second, tt.fraction, elong2, u, 0.0)
        Δtdb = tdb_tt2 - tdb_tt1
        tdb1 = TDBEpoch(tdb_tt1, tt)
        tdb2 = TDBEpoch(tdb_tt2, tt)
        @test value(tdb2 - tdb1) ≈ Δtdb
        @test tdb1 != tdb2
        tdb1 = TDBEpoch(tt, elong1, u, 0.0)
        tdb2 = TDBEpoch(tt, elong2, u, 0.0)
        @test value(tdb2 - tdb1) ≈ Δtdb
        @test tdb1 != tdb2

        t0 = TTEpoch(2000, 1, 1, 12, 0, 32.0)
        t1 = TAIEpoch(2000, 1, 1, 12, 0, 32.0)
        t2 = TAIEpoch(2000, 1, 1, 12, 0, 0.0)
        @test_throws MethodError t1 - t0
        @test_throws MethodError t1 < t0
        @test t2 - t1 == -32.0seconds
        @test t2 < t1

        today = TTEpoch(2000, 1, 1, 12, 0, 13.123)
        age_of_the_universe = 13.772e9years
        big_bang = today - age_of_the_universe
        baryons = big_bang + 1e-11seconds
        @test baryons + age_of_the_universe - today == 1e-11seconds

        reception_time = TDBEpoch("2021-07-01T00:00:00.00")
        rtlt_a = seconds(1.5days)
        rtlt_b = rtlt_a + 1e-6seconds
        transmission_time_a = reception_time + rtlt_a
        transmission_time_b = reception_time + rtlt_b
        @test transmission_time_b - transmission_time_a == 1e-6seconds
    end
    @testset "Parsing" begin
        @test AstroTime.TimeScales.tryparse(1.0) === nothing
        @test TAIEpoch("2000-01-01T00:00:00.000") == TAIEpoch(2000, 1, 1)
        @test UT1Epoch("2000-01-01T00:00:00.000") == UT1Epoch(2000, 1, 1)
        @test TTEpoch("2000-01-01T00:00:00.000") == TTEpoch(2000, 1, 1)
        @test TCGEpoch("2000-01-01T00:00:00.000") == TCGEpoch(2000, 1, 1)
        @test TCBEpoch("2000-01-01T00:00:00.000") == TCBEpoch(2000, 1, 1)
        @test TDBEpoch("2000-01-01T00:00:00.000") == TDBEpoch(2000, 1, 1)
        @test Epoch("2000-01-01T00:00:00.000 TAI") == TAIEpoch(2000, 1, 1)
        @test TAIEpoch("2000-001", "yyyy-DDD") == TAIEpoch(2000, 1, 1)
        @test Epoch("2000-001 TAI", "yyyy-DDD ttt") == TAIEpoch(2000, 1, 1)
        @test_throws ArgumentError Epoch("2000-01-01T00:00:00.000")
    end
    @testset "Output" begin
        ep = TAIEpoch(2018, 8, 14, 10, 2, 51.551247436378276)
        @test AstroTime.format(ep, "yyyy-DDDTHH:MM:SS.sss") == "2018-226T10:02:51.551"
        @test AstroTime.format(ep, "HH:MM ttt") == "10:02 TAI"
        @test string(TAI) == "TAI"
        @test string(TT) == "TT"
        @test string(UT1) == "UT1"
        @test string(TCG) == "TCG"
        @test string(TDB) == "TDB"
        @test string(TCB) == "TCB"
    end
    @testset "Arithmetic" begin
        ep = TAIEpoch(2000, 1, 1)
        ep1 = TAIEpoch(2000, 1, 2)
        @test (ep + 1.0seconds) - ep   == 1.0seconds
        @test (ep + 1.0minutes) - ep   == seconds(1.0minutes)
        @test (ep + 1.0hours) - ep     == seconds(1.0hours)
        @test (ep + 1.0days) - ep      == seconds(1.0days)
        @test (ep + 1.0years) - ep     == seconds(1.0years)
        @test (ep + 1.0centuries) - ep == seconds(1.0centuries)
        @test ep < ep1
        @test isless(ep, ep1)
    end
    @testset "Conversion" begin
        include("conversions.jl")
        dt = AstroTime.DateTime(2018, 8, 14, 10, 2, 51.551247436378276)
        ep = TAIEpoch(2018, 8, 14, 10, 2, 51.551247436378276)
        @test TAIEpoch(dt) == ep
        @test TAIEpoch(Dates.DateTime(dt)) == TAIEpoch(2018, 8, 14, 10, 2, 51.551)
        @test TAIEpoch(AstroTime.Date(2018, 8, 14)) == TAIEpoch(2018, 8, 14, 0, 0, 0.0)
        @test now(Epoch) isa TAIEpoch
        @test now(TDBEpoch) isa TDBEpoch

        tt = TTEpoch(2000, 1, 1, 12)
        @test TTEpoch(tt) == tt
        @test Epoch{TerrestrialTime,Float64}(tt) == tt
        @test tt - J2000_EPOCH == 0.0seconds
        tai = TAIEpoch(2000, 1, 1, 12)
        @test tai.second == 0
        @test tai.fraction == 0.0
        @test TTEpoch(tai) ≈ TTEpoch(2000, 1, 1, 12, 0, 32.184)
        @test Epoch(tai, TT) ≈ TTEpoch(2000, 1, 1, 12, 0, 32.184)
        @test TTEpoch(32.184, tai) == TTEpoch(tai)

        ut1 = UT1Epoch(2000, 1, 1)
        ut1_tai = getoffset(ut1, TAI)
        tai = TAIEpoch(ut1)
        tai_tt = getoffset(tai, TT)
        tt = TTEpoch(tai)
        tt_tdb = getoffset(tt, TDB)
        tdb = TDBEpoch(tt)
        tdb_tcb = getoffset(tdb, TCB)
        tcb = TCBEpoch(tdb)
        @test getoffset(ut1, TCB) == ut1_tai + tai_tt + tt_tdb + tdb_tcb
    end
    @testset "TDB" begin
        ep = TTEpoch(2000, 1, 1)
        @test TDBEpoch(ep) ≈ TDBEpoch(ep, 0.0, 0.0, 0.0) rtol=1e-3
        jd1, jd2 = value.(julian_twopart(ep))
        ut = fractionofday(UT1Epoch(ep))
        elong, u, v = abs.(randn(3)) * 1000
        exp = ERFA.dtdb(jd1, jd2, ut, elong, u, v)
        act = getoffset(ep, TDB, elong, u, v)
        @test act ≈ exp
        second, fraction = 394372865, 0.1839999999999975
        offset = getoffset(TT, TDB, second, fraction)
        @test offset ≈ 0.105187547186749e-3
    end
    @testset "Julian Dates" begin
        jd = 0.0days
        ep = TAIEpoch(jd)
        @test ep == TAIEpoch(2000, 1, 1, 12)
        @test julian_period(ep) == 0.0days
        @test julian_period(ep; scale=TT, unit=seconds) == 32.184seconds
        @test julian_period(Float64, ep) == 0.0
        @test j2000(ep) == jd
        jd = 86400.0seconds
        ep = TAIEpoch(jd)
        @test ep == TAIEpoch(2000, 1, 2, 12)
        @test j2000(ep) == days(jd)
        jd = 2.451545e6days
        ep = TAIEpoch(jd, origin=:julian)
        @test ep == TAIEpoch(2000, 1, 1, 12)
        @test julian(ep) == jd
        jd = 51544.5days
        ep = TAIEpoch(jd, origin=:modified_julian)
        @test ep == TAIEpoch(2000, 1, 1, 12)
        @test modified_julian(ep) == jd
        @test_throws ArgumentError TAIEpoch(jd, origin=:julia)
    end
    @testset "Accessors" begin
        @test TAIEpoch(JULIAN_EPOCH - Inf * seconds) == PAST_INFINITY
        @test TAIEpoch(JULIAN_EPOCH + Inf * seconds) == FUTURE_INFINITY
        @test string(PAST_INFINITY) == "-5877490-03-03T00:00:00.000 TAI"
        @test string(FUTURE_INFINITY) == "5881610-07-11T23:59:59.999 TAI"
        y = 2018
        m = 2
        d = 6
        hr = 20
        mn = 45
        sec = 59.371248965
        ep = TAIEpoch(y, m, d, hr, mn, sec)
        @test year(ep) == y
        @test month(ep) == m
        @test day(ep) == d
        @test hour(ep) == hr
        @test minute(ep) == mn
        @test second(Float64, ep) == sec
        @test second(Int, ep) == 59
        @test millisecond(ep) == 371
        @test microsecond(ep) == 248
        @test nanosecond(ep) == 965
        @test subsecond(ep, 9) == nanosecond(ep)
        @test yearmonthday(ep) == (y, m, d)
        @test fractionofsecond(ep) ≈ 0.371248965
        @test AstroTime.Date(ep) == AstroTime.Date(y, m, d)
        @test AstroTime.Time(ep) == AstroTime.Time(hr, mn, sec)
        @test AstroTime.DateTime(ep) == AstroTime.DateTime(y, m, d, hr, mn, sec)
        @test Dates.Date(ep) == Dates.Date(y, m, d)
        @test Dates.Time(ep) == Dates.Time(hr, mn, 59, 371, 248, 965)
        @test Dates.DateTime(ep) == Dates.DateTime(y, m, d, hr, mn, 59, 371)
    end
    @testset "Ranges" begin
        rng = TAIEpoch(2018, 1, 1):TAIEpoch(2018, 1, 2)
        @test step(rng) == 1seconds
        @test length(rng) == 86401
        @test first(rng) == TAIEpoch(2018, 1, 1)
        @test last(rng) == TAIEpoch(2018, 1, 2)
        rng = TAIEpoch(2018, 1, 1):13seconds:TAIEpoch(2018, 1, 1, 0, 1)
        @test step(rng) == 13seconds
        @test last(rng) == TAIEpoch(2018, 1, 1, 0, 0, 52.0)
    end
    @testset "Parametrization" begin
        ep_f64 = TAIEpoch(2000, 1, 1)
        ep_err = TAIEpoch(ep_f64.second, 1.0 ± 1.1)
        Δt = (30 ± 0.1) * seconds
        @test typeof(Δt) == AstroPeriod{AstroTime.Periods.Second,Measurement{Float64}}
        @test typeof(ep_f64) == Epoch{InternationalAtomicTime,Float64}
        @test typeof(ep_err) == Epoch{InternationalAtomicTime,Measurement{Float64}}
        @test typeof(ep_f64 + Δt) == Epoch{InternationalAtomicTime,Measurement{Float64}}
        @test typeof(ep_err + Δt) == Epoch{InternationalAtomicTime,Measurement{Float64}}
        jd1_err = (0.0 ± 0.001) * days
        jd2_err = (0.5 ± 0.001) * days
        ep_jd1 = TAIEpoch(jd1_err)
        @test typeof(ep_jd1) == Epoch{InternationalAtomicTime,Measurement{Float64}}
        ep_jd2 = TAIEpoch(jd1_err, jd2_err)
        @test typeof(ep_jd2) == Epoch{InternationalAtomicTime,Measurement{Float64}}
        ut1_err = UT1Epoch(ep_f64.second, 1.0 ± 1.1)
        tcg_err = TCGEpoch(ut1_err)
        tcb_err = TCBEpoch(ut1_err)
        @test typeof(ut1_err) == Epoch{UniversalTime,Measurement{Float64}}
        @test typeof(tcg_err) == Epoch{GeocentricCoordinateTime,Measurement{Float64}}
        @test typeof(tcb_err) == Epoch{BarycentricCoordinateTime,Measurement{Float64}}
        @test typeof(UT1Epoch(tcg_err)) == Epoch{UniversalTime,Measurement{Float64}}
        @test typeof(UT1Epoch(tcb_err)) == Epoch{UniversalTime,Measurement{Float64}}
    end
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
module Macros

export @def, @load_struct, @load_struct_cat, @save_struct, @save_struct_index, @load_index, @load_index_cat, init!, @load_vec_into_struct

"""
	@def name definition

Defines the macro `name`, which writes the code in `definition`
"""
macro def(name, definition)
  return quote
      macro $(esc(name))()
          esc($(Expr(:quote, definition)))
      end
  end
end

"""
	@load_struct_cat name s v1 v2 ...

Load variables `v1`, `v2`, ... from structure `s` and rename them as suffix `v1_name`, `v2_name`, ...
"""
macro load_struct_cat(name, d, vars...)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :($(Symbol("$(v)_$(name)")) = $d.$v))
	end
	return esc(expr)
end

"""
	@load_index d i v1 v2 ...

Loads variables `v1`, `v2`, ... from structure `d` from the index `i` of the associated fields.
"""
macro load_index(d, i, vars...)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :($v = ($d.$v)[$i]))
	end
	return esc(expr)
end

"""
	@load_index_cat d i cat v1 v2 ...

Loads variables `v1_cat`, `v2_cat`, ... from structure `d` from the index `i` of the associated fields `v1`, `v2`, ... .
"""
macro load_index_cat(d, i, cat, vars...)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :($(Symbol("$(v)_$(cat)")) = ($d.$v)[$i]))
	end
	return esc(expr)
end

"""
	@load_struct d v1 ...

Loads variables `v1`, ... from the structure `d`
"""
macro load_struct(d, vars...)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :($v = $d.$v))
	end
	return esc(expr)
end

"""
	@save_struct d v1 ...

Saves local variables `v1`, ... in the associated slot `d.v1` of structure `d`
"""
macro save_struct(d, vars...)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :($d.$v = $v))
	end
	return esc(expr)
end

"""
	@save_struct_index d i v1 ...

Saves local variables `v1`, ... at the index `i` of the associated slot `(d.v1)[i]` of structure `d`
"""
macro save_struct_index(d, i, vars...)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :(($d.$v)[$i] = $v))
	end
	return esc(expr)
end

"""
	@init_struct s v

Initializes the properties of structure `s` to value `v`
"""
macro init_struct(s, v)
	expr = Expr(:block)
	for v in vars
		push!(expr.args, :($s.$v = $v))
	end
	return expr
end

"""
	init!(s::S)
"""
@generated function init!(s::S, N) where S
	expr = Expr(:block)
	for v in fieldnames(S)
		push!(expr.args, :(s.$v = zeros(N)))
	end
	return expr
end

"""
	@load_vec_into_struct vec s v1 ...

Assign values `vec[1]`, `vec[2]`, ... to  fields `v1`, `v2`, ... of structure `s`
"""
macro load_vec_into_struct(vec, s, vars...)
	expr = Expr(:block)
	for (k,v) in enumerate(vars)
		push!(expr.args, :($s.$v = $vec[$k]))
	end
	return esc(expr)
end

end


text:
#=
is_unique_chars:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-21
=#

module UniqueCharChecker
    # parameters
    const max_len = 128

    # methods
    function check_unique(str)
        println("")
        println("Check string $(str) is unique or not")

        if length(str) > max_len
            println("Too long string")
            return false
        end

        char_set = [false for i in 1:max_len]

        for j in 1:length(str)
            ascii = Int(str[j]) # used as idx of char_set
            if char_set[ascii] == true
                println("Not unique string")
                return false
            end
            char_set[ascii] = true
        end
        println("Unique string")
        return true
    end

    function main()
        @time println(check_unique("abcdefg"))
        @time println(check_unique("aaaaaaaaaaaa"))
        @time println(check_unique("aAbBcCdD"))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .UniqueCharChecker
    UniqueCharChecker.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
struct PaintBucket{Tbnd}
    N::Int64
    x_f_start::Float64
    x_f_stop::Float64
    x_f::LinRange{Float64}
    x_t_start::Float64
    x_t_stop::Float64
    x_t::LinRange{Float64}
    y_f_start::Float64
    y_f_stop::Float64
    y_f::LinRange{Float64}
    y_t_start::Float64
    y_t_stop::Float64
    y_t::LinRange{Float64}
    ε::Array{Complex{Float64},2}
    F::Array{Float64,2}
    ψ_plot::Array{Complex{Float64},3}
    boundary_data::Tbnd
    pml_data::NTuple{4,Tuple{Array{Float64,1},Array{Float64,1}}}

    function PaintBucket(sim::Simulation{Tsys,Tbnd,Discretization{Polar}}, ψ, truncate::Bool) where {Tsys,Tbnd}

        r_f = sim.dis.x[1][1] .+ sim.dis.dx[1]*(-1:sim.dis.N[1])
        θ_f = sim.dis.x[2][1] .+ sim.dis.dx[2]*(-1:sim.dis.N[2])

        x_f_start = minimum(sim.dis.X)-sim.dis.dx[1]/2
        x_f_stop  = maximum(sim.dis.X)+sim.dis.dx[1]/2
        x_f = LinRange(x_f_start, x_f_stop, 2sim.dis.N[1])

        y_f_start = minimum(sim.dis.Y)-sim.dis.dx[1]/2
        y_f_stop  = maximum(sim.dis.Y)+sim.dis.dx[1]/2
        y_f = LinRange(y_f_start, y_f_stop, 2sim.dis.N[1])

        ε_temp = Array{ComplexF64}(undef, sim.dis.N[1]+2, sim.dis.N[2]+2)
        ε_temp[2:end-1,2:end-1] = sim.sys.ε
        @views ε_temp[:,1] = ε_temp[:,end-1]
        @views ε_temp[:,end] = ε_temp[:,2]
        @views ε_temp[1,:] = ε_temp[2,:]
        @views ε_temp[end,:] = ε_temp[end-1,:]
        # ε_itp = extrapolate(scale(interpolate(ε_temp,BSpline(Quadratic(Periodic(OnGrid())))),r_f,θ_f),complex(NaN,NaN))
        ε_itp = extrapolate(scale(interpolate(ε_temp,BSpline(Linear())),r_f,θ_f),complex(NaN,NaN))
        ε = @. ε_itp(hypot(x_f,y_f'),sim.dis.x[2][1]+mod(atan(y_f',x_f)-sim.dis.x[2][1],2π))

        F_temp = Array{Float64}(undef, sim.dis.N[1]+2, sim.dis.N[2]+2)
        F_temp[2:end-1,2:end-1] = sim.sys.F
        @views F_temp[:,1] = F_temp[:,end-1]
        @views F_temp[:,end] = F_temp[:,2]
        @views F_temp[1,:] = F_temp[2,:]
        @views F_temp[end,:] = F_temp[end-1,:]
        # F_itp = extrapolate(scale(interpolate(F_temp,BSpline(Cubic(Periodic(OnGrid())))),r_f,θ_f),NaN)
        F_itp = extrapolate(scale(interpolate(F_temp,BSpline(Linear())),r_f,θ_f),NaN)
        F = @. F_itp(hypot(x_f,y_f'),sim.dis.x[2][1]+mod(atan(y_f',x_f)-sim.dis.x[2][1],2π))

        N, M, X, Y, x_t, y_t, ∂Ω, ψ = apply_truncation(sim, ψ, truncate)

        r_t = x_t[1] .+ sim.dis.dx[1]*(-1:M[1])
        θ_t = y_t[1] .+ sim.dis.dx[2]*(-1:M[2])

        x_t_start = minimum(X)-sim.dis.dx[1]/2
        x_t_stop  = maximum(X)+sim.dis.dx[1]/2
        x_t = LinRange(x_t_start, x_t_stop, 2M[1])

        y_t_start = minimum(Y)-sim.dis.dx[1]/2
        y_t_stop  = maximum(Y)+sim.dis.dx[1]/2
        y_t = LinRange(y_t_start, y_t_stop, 2M[1])

        # root_r = sqrt.(hypot.(X,Y))
        ψ_temp = Array{ComplexF64}(undef, M[1]+2, M[2]+2)
        ψ_plot = Array{ComplexF64}(undef, 2M[1], 2M[1], N)
        for i ∈ 1:N
            @views ψ_temp[2:end-1,2:end-1] = reshape(ψ[:,i]/maximum(abs.(ψ[:,i])), M[1], M[2])#./root_r
            @views ψ_temp[:,1] = ψ_temp[:,end-1]
            @views ψ_temp[:,end] = ψ_temp[:,2]
            @views ψ_temp[1,:] = ψ_temp[2,:]
            @views ψ_temp[end,:] = ψ_temp[end-1,:]
            ψ_itp = extrapolate(scale(interpolate(ψ_temp,BSpline(Cubic(Periodic(OnGrid())))),r_t,θ_t),complex(NaN,NaN))
            ψ_plot[:,:,i] = @. ψ_itp(hypot(x_t,y_t'),sim.dis.x[2][1]+mod(atan(y_t',x_t)-sim.dis.x[2][1],2π))
        end

        boundary_data1_f = (r_f[1]*cos.(θ_f),r_f[1]*sin.(θ_f))
        boundary_data2_f = (r_f[end]*cos.(θ_f),r_f[end]*sin.(θ_f))
        boundary_data3_f = (r_f*cos(θ_f[1]),r_f*sin(θ_f[1]))
        boundary_data4_f = (r_f*cos(θ_f[end]),r_f*sin(θ_f[end]))

        N_pml = 101
        pml_r_1 = vcat(
                    LinRange(sim.dis.x[1][1], sim.dis.x_tr[1][1], N_pml),
                    fill(sim.dis.x_tr[1][1], N_pml),
                    LinRange(sim.dis.x_tr[1][1], sim.dis.x[1][1], N_pml),
                fill(sim.dis.x[1][1],N_pml)
                )
        pml_θ_1 = vcat(
                    fill(sim.dis.x[2][1],N_pml),
                    LinRange(sim.dis.x[2][1], sim.dis.x[2][end], N_pml),
                    fill(sim.dis.x[2][end],N_pml),
                    LinRange(sim.dis.x[2][end], sim.dis.x[2][1], N_pml)
                    )

        pml_r_2 = vcat(
                    LinRange(sim.dis.x_tr[1][end], sim.dis.x[1][end], N_pml),
                    fill(sim.dis.x[1][end], N_pml),
                    LinRange(sim.dis.x[1][end], sim.dis.x_tr[1][end], N_pml),
                    fill(sim.dis.x_tr[1][end],N_pml)
                    )
        pml_θ_2 = vcat(
                    fill(sim.dis.x[2][1],N_pml),
                    LinRange(sim.dis.x[2][1], sim.dis.x[2][end], N_pml),
                    fill(sim.dis.x[2][end],N_pml),
                    LinRange(sim.dis.x[2][end], sim.dis.x[2][1], N_pml)
                    )

        pml_r_3 = vcat(
                    LinRange(sim.dis.x[1][1], sim.dis.x[1][end], N_pml),
                    fill(sim.dis.x[1][end], N_pml),
                    LinRange(sim.dis.x[1][end], sim.dis.x[1][1], N_pml),
                    fill(sim.dis.x[1][1],N_pml)
                    )
        pml_θ_3 = vcat(
                    fill(sim.dis.x[2][1],N_pml),
                    LinRange(sim.dis.x[2][1], sim.dis.x_tr[2][1], N_pml),
                    fill(sim.dis.x_tr[2][1],N_pml),
                    LinRange(sim.dis.x_tr[2][1], sim.dis.x[2][1], N_pml)
                    )

        pml_r_4 = vcat(
                    LinRange(sim.dis.x[1][1], sim.dis.x[1][end], N_pml),
                    fill(sim.dis.x[1][end], N_pml),
                    LinRange(sim.dis.x[1][end], sim.dis.x[1][1], N_pml),
                    fill(sim.dis.x[1][1],N_pml)
                    )
        pml_θ_4 = vcat(
                    fill(sim.dis.x_tr[2][end],N_pml),
                    LinRange(sim.dis.x_tr[2][end], sim.dis.x[2][end], N_pml),
                    fill(sim.dis.x[2][end],N_pml),
                    LinRange(sim.dis.x[2][end], sim.dis.x_tr[2][end], N_pml)
                    )

        pml_data1 = (pml_r_1.*cos.(pml_θ_1),pml_r_1.*sin.(pml_θ_1))
        pml_data2 = (pml_r_2.*cos.(pml_θ_2),pml_r_2.*sin.(pml_θ_2))
        pml_data3 = (pml_r_3.*cos.(pml_θ_3),pml_r_3.*sin.(pml_θ_3))
        pml_data4 = (pml_r_4.*cos.(pml_θ_4),pml_r_4.*sin.(pml_θ_4))

        bnd_data = (boundary_data1_f, boundary_data2_f, boundary_data3_f, boundary_data4_f)

        return new{typeof(bnd_data)}(N, x_f_start, x_f_stop, x_f, x_t_start, x_t_stop, x_t,
            y_f_start, y_f_stop, y_f, y_t_start, y_t_stop, y_t,
            ε, F, ψ_plot,
            bnd_data,
            (pml_data1, pml_data2, pml_data3, pml_data4) )
    end

    function PaintBucket(sim::Simulation{Tsys,Tbnd,Discretization{Cartesian}}, ψ, truncate::Bool) where Tsys<:System where Tbnd<:Boundary

        N, M, X, Y, x_t, y_t, ∂Ω, ψ = apply_truncation(sim, ψ, truncate)

        x_f_start = sim.bnd.∂Ω[1][1]
        x_f_stop  = sim.bnd.∂Ω[1][2]
        x_f = LinRange(sim.dis.x[1][1],sim.dis.x[1][end],length(sim.dis.x[1]))#sim.dis.x[1][:]

        y_f_start = sim.bnd.∂Ω[2][1]
        y_f_stop  = sim.bnd.∂Ω[2][2]
        y_f = LinRange(sim.dis.x[2][1],sim.dis.x[2][end],length(sim.dis.x[2]))#sim.dis.x[2][:]

        x_t_start = ∂Ω[1][1]
        x_t_stop  = ∂Ω[1][2]
        x_t = LinRange(x_t_start, x_t_stop, M[1])

        y_t_start = ∂Ω[2][1]
        y_t_stop  = ∂Ω[2][2]
        y_t = LinRange(y_t_start, y_t_stop, M[2])

        ε = sim.sys.ε
        # ε = (1 .+ 1im*sim.sys.Σe/10).*sim.sys.ε
        F = sim.sys.F

        ψ_plot = Array{ComplexF64}(undef, M[1], M[2], N)
        for i in 1:N
            ψ_plot[:,:,i] = reshape(ψ[:,i]/maximum(abs.(ψ[:,i])), M[1], M[2])
        end

        boundary_data1_f = ([x_f[1],x_f[1]],[y_f[1],y_f[end]])
        boundary_data2_f = ([x_f[end],x_f[end]],[y_f[1],y_f[end]])
        boundary_data3_f = ([x_f[1],x_f[end]],[y_f[1],y_f[1]])
        boundary_data4_f = ([x_f[1],x_f[end]],[y_f[end],y_f[end]])

        N_pml = 2
        pml_x_1 = vcat(
                    LinRange(sim.dis.x[1][1], sim.dis.x_tr[1][1], N_pml),
                    fill(sim.dis.x_tr[1][1], N_pml),
                    LinRange(sim.dis.x_tr[1][1], sim.dis.x[1][1], N_pml),
                    fill(sim.dis.x[1][1],N_pml)
                    )
        pml_y_1 = vcat(
                    fill(sim.dis.x[2][1],N_pml),
                    LinRange(sim.dis.x[2][1], sim.dis.x[2][end], N_pml),
                    fill(sim.dis.x[2][end],N_pml),
                    LinRange(sim.dis.x[2][end], sim.dis.x[2][1], N_pml)
                    )

        pml_x_2 = vcat(
                    LinRange(sim.dis.x_tr[1][end], sim.dis.x[1][end], N_pml),
                    fill(sim.dis.x[1][end], N_pml),
                    LinRange(sim.dis.x[1][end], sim.dis.x_tr[1][end], N_pml),
                    fill(sim.dis.x_tr[1][end],N_pml)
                    )
        pml_y_2 = vcat(
                    fill(sim.dis.x[2][1],N_pml),
                    LinRange(sim.dis.x[2][1], sim.dis.x[2][end], N_pml),
                    fill(sim.dis.x[2][end],N_pml),
                    LinRange(sim.dis.x[2][end], sim.dis.x[2][1], N_pml)
                    )

        pml_x_3 = vcat(
                    LinRange(sim.dis.x[1][1], sim.dis.x[1][end], N_pml),
                    fill(sim.dis.x[1][end], N_pml),
                    LinRange(sim.dis.x[1][end], sim.dis.x[1][1], N_pml),
                    fill(sim.dis.x[1][1],N_pml)
                    )
        pml_y_3 = vcat(
                    fill(sim.dis.x[2][1],N_pml),
                    LinRange(sim.dis.x[2][1], sim.dis.x_tr[2][1], N_pml),
                    fill(sim.dis.x_tr[2][1],N_pml),
                    LinRange(sim.dis.x_tr[2][1], sim.dis.x[2][1], N_pml)
                    )

        pml_x_4 = vcat(
                    LinRange(sim.dis.x[1][1], sim.dis.x[1][end], N_pml),
                    fill(sim.dis.x[1][end], N_pml),
                    LinRange(sim.dis.x[1][end], sim.dis.x[1][1], N_pml),
                    fill(sim.dis.x[1][1],N_pml)
                    )
        pml_y_4 = vcat(
                    fill(sim.dis.x_tr[2][end],N_pml),
                    LinRange(sim.dis.x_tr[2][end], sim.dis.x[2][end], N_pml),
                    fill(sim.dis.x[2][end],N_pml),
                    LinRange(sim.dis.x[2][end], sim.dis.x_tr[2][end], N_pml)
                    )

        pml_data1 = (pml_x_1,pml_y_1)
        pml_data2 = (pml_x_2,pml_y_2)
        pml_data3 = (pml_x_3,pml_y_3)
        pml_data4 = (pml_x_4,pml_y_4)

        M[1]==1 ? x_t = y_t : nothing
        M[1]==1 ? x_f_start = y_f_start : nothing
        M[1]==1 ? x_f_stop = y_f_stop : nothing
        M[1]==1 ? x_f = y_f : nothing
        M[1]==1 ? x_t_start = y_t_start : nothing
        M[1]==1 ? x_t_stop = y_t_stop : nothing

        bnd_data = (boundary_data1_f, boundary_data2_f, boundary_data3_f, boundary_data4_f)
        return new{typeof(bnd_data)}(N, x_f_start, x_f_stop, x_f, x_t_start, x_t_stop, x_t[:],
            y_f_start, y_f_stop, y_f, y_t_start, y_t_stop, y_t[:],
            ε, F, ψ_plot,
            bnd_data,
            (pml_data1, pml_data2, pml_data3, pml_data4) )
    end
end

"""
    p = plot(domain)
"""
@recipe function f(dom::Domain)
    (dom.is_in_domain)
end

"""
    p = plot(sim; by=nothing, truncate=false)
"""
@recipe function f(sim::Simulation; by=nothing, truncate=false)
    pb = PaintBucket(sim, Array{ComplexF64}(undef,prod(sim.dis.N),1), truncate)
    if isnothing(by)
        layout --> (1,3)
        (sim,pb,by,1)
    else
        (sim,pb,by,1,1)
    end
end

# internal use, sim summary
@recipe function f(sim::Simulation, pb::PaintBucket, by, _1)
    bys = [:real, :imag, :F]
    for i ∈ 1:3
        @series begin
            subplot := i
            (sim, pb, bys[i], 1, 1)
        end
    end
end

# interal use only, sim elementary
@recipe function f(sim::Simulation, pb::PaintBucket, by, _1, _2)
    aspect_ratio --> 1
    grid --> false
    xlims --> [pb.x_t_start,pb.x_t_stop]
    ylims --> [pb.y_t_start,pb.y_t_stop]
    colorbar --> false
    overwrite_figure --> false
    levels --> 15
    lw --> 6
    seriestype --> :heatmap
    legend --> false
    framestyle --> :grid

    if by ∈ [:real, real, "real"]
        color --> :sequential
        title --> "Re n(x)"
        p = real(sqrt.(pb.ɛ))
        pm = p[.!isnan.(p)]
        clims := (minimum(pm), maximum(pm))
    elseif by ∈ [:imag, imag, "imag"]
        color --> :diverging
        title --> "Im n(x)"
        p = imag(sqrt.(pb.ɛ))
        pm = p[.!isnan.(p)]
        clims := (-maximum(abs.(pm)), maximum(abs.(pm)))
    elseif by ∈ [:F,"F"]
        color --> :diverging
        title --> "F(x)"
        p = pb.F
        pm = p[.!isnan.(p)]
        clims := (-maximum(abs.(pm)), maximum(abs.(pm)))
    else
        throw(ArgumentError("unrecognized `by` keyword"))
    end
    @series begin
        pb.x_f, pb.y_f, permutedims(p)
    end
    for k ∈ 1:2 # cycle over dims
        for i ∈ eachindex(sim.bnd.bc) # cycle over sides
            if isDirichlet(sim.bnd.bc[k][i])
                linestyle := :solid
            elseif isNeumann(sim.bnd.bc[k][i])
                linestyle := :dash
            end
            if !isnoBL(sim.bnd.bl[k][i])
                if isPML(sim.bnd.bl[k][i])
                    fill_color = :red
                elseif iscPML(sim.bnd.bl[k][i])
                    fill_color = :blue
                end
                @series begin
                    lw := 0
                    seriestype := :path
                    color := fill_color
                    alpha := .25
                    fill --> (0,.5,fill_color)
                    pb.pml_data[2*(k-1)+i]
                end
            end
            if isDirichlet(sim.bnd.bc[k][i]) | isNeumann(sim.bnd.bc[k][i])
                @series begin
                    lw := 2
                    color := :black
                    seriestype := :path
                    pb.boundary_data[i]
                end
            end
        end
    end
end

"""
    p = plot(sim, ψ, [inds; by=nothing, summary=false, truncate=false])
"""
@recipe function f(sim::Simulation, ψ::AbstractArray, inds=1:size(ψ,2); by=nothing, structure=false, truncate=false)
    pb = PaintBucket(sim, ψ[:,inds], truncate)
    n = length(inds)
    if isnothing(by)
        if structure
            layout --> (1+n,3)
        else
            layout --> (n,3)
        end
    else
        if structure
            layout --> (3+n)
        else
            layout --> (n)
        end
    end
    if structure
        @series (sim,pb,by,1)
    end
    if isnothing(by)
        for i ∈ 1:n
            @series begin
                subplot := 3*(structure+i-1)+1
                (sim,pb,i,:real,1,1,1)
            end
            @series begin
                subplot := 3*(structure+i-1)+2
                (sim,pb,i,:imag,1,1,1)
            end
            @series begin
                subplot := 3*(structure+i-1)+3
                (sim,pb,i,:abs2,1,1,1)
            end
        end
    else
        for i ∈ 1:n
            @series begin
                subplot := 3*structure+i
                (sim,pb,i,by,1,1,1)
            end
        end
    end
end

# elementary plot sim of ψ[:,:,i]
@recipe function f(sim::Simulation,pb::PaintBucket,i::Int,by,_1,_2,_3)
    aspect_ratio --> 1
    grid --> false
    xlims --> [pb.x_t_start,pb.x_t_stop]
    ylims --> [pb.y_t_start,pb.y_t_stop]
    colorbar --> false
    overwrite_figure --> false
    levels --> 15
    lw --> 6
    seriestype --> :heatmap
    legend --> false
    framestyle --> :grid

    p = pb.ψ_plot[:,:,i]
    if by ∈ [:real,real,"real"]
        color --> :diverging
        title --> "Re psi(x)"
        p = real(p)
    elseif by ∈ [:imag,imag,"imag"]
        color --> :diverging
        title --> "Im psi(x)"
        p = imag(p)
    elseif by ∈ [:abs,abs,"abs"]
        color --> :sequential
        title --> "|psi(x)|"
        p = abs.(p)
    elseif by ∈ [:abs2,abs2,"abs2"]
        color --> :sequential
        title --> "|psi(x)|²"
        p = abs2.(p)
    else
        throw(ArgumentError("unrecognized `by` keyword"))
    end
    @series begin
        pb.x_t, pb.y_t, permutedims(p)
    end
end

function apply_truncation(sim::Simulation, ψ, truncate::Bool)
    if isempty(ψ)
        N=0
    else
        N = size(ψ,2)
    end
    if !isempty(ψ)
        idx = truncate ? sim.dis.X_idx : 1:prod(sim.dis.N)
    else
        idx = Int[]
    end
    if !truncate || iszero(N)
        M = sim.dis.N
        X, Y = sim.dis.X, sim.dis.Y
        x, y = sim.dis.x[1], sim.dis.x[2]
        ∂Ω = sim.bnd.∂Ω
    else
        M = sim.dis.N_tr
        X = reshape(sim.dis.X[sim.dis.X_idx],M[1],M[2])
        Y = reshape(sim.dis.Y[sim.dis.X_idx],M[1],M[2])
        x, y = sim.dis.x_tr[1], sim.dis.x_tr[2]
        ∂Ω = sim.bnd.∂Ω_tr
        ψ = ψ[idx,:]
    end
    return N, M, X, Y, x, y, ∂Ω, ψ
end

# ################################################################################
# ########## ANIMATION
# ################################################################################
# """
#     iterator = wave(sim, ψ; by=real, n=60, seriestype=:heatmap)
#
# input for Plots.animate
#
# Use cases:
#
# `animate(wave(sim,ψ), file_name)` creates a .gif with filename
#
# `animate(wave(sim,ψ; n=20), file_name, fps=10)` creates a 2 second movie
#
# Note: default `fps`=20, and `n`=60, so default movie is 3 seconds long
# """
# function wave(sim::Simulation, ψ; truncate=true, by=real, n=60)
#     if 1 ∈ sim.dis.N
#         return imap( ϕ->(sim, exp(-1im*ϕ)*ψ, by, truncate, 1), 0:2π/n:2π*(1-1/n))
#     else
#         return imap( ϕ->(sim, exp(-1im*ϕ)*ψ, by, truncate, 1, 1), 0:2π/n:2π*(1-1/n))
#     end
# end