text:
#!/usr/bin/env julia

#------------------------------------------------------------------------------#
#                                 mcmc.jl                                      #
#------------------------------------------------------------------------------#

# First, import some dependencies.

using ArgParse
using FastaIO
using JSON

# Assign the AMBIGS dictionary.

const AMBIGS = Dict(
    "A" => 'A', "C" => 'C', "G" => 'G', "T" => 'T',
    "AC" => 'M', "AG" => 'R', "AT" => 'W', "CG" => 'S',
    "CT" => 'Y', "GT" => 'K', "ACG" => 'V', "ACT" => 'H',
    "AGT" => 'D', "CGT" => 'B', "ACGT" => 'N', "" => '-'
)

#------------------------------------------------------------------------------#
#                          functional definitions                              #
#------------------------------------------------------------------------------#

#----------------------------------index---------------------------------------#

# Define a function that returns all indices at which a string matches a
# certain character.

function index(s::String, c::Char)

    # For each character in the string, check if it matches the input character
    # and if it does, return the index.

    for i in 1:length(s)
        @inbounds s_i = s[i]
        if c == s_i
            return i
        end
    end

    # Return 0 if there are no matches.

    0
end

#---------------------------------countmsa-------------------------------------#

# Define a function to count, for each character of an alphabet, the number
# of times that character appears in each sequence of an MSA file.

function countmsa(msa::String, alphabet::String)

    # Initialize the counters, counts and ncols.

    counts = nothing
    ncols = 0
    tic()

    # For each sequence in the input .msa file,

    for (i, (name, seq)) in enumerate(FastaReader(msa))

        # If the counters are empty, assign to ncols the length of the current
        # sequence and initialize counts as a vector of zeros with length ncols.

        if counts == nothing
            ncols = length(seq)
            counts = zeros(Int64, (ncols, 4))

        # If the counters are not empty, check that ncols == length(seq), and
        # say something if they are not equal.

        elseif length(seq) != ncols
            error("provided file is not an MSA! length($name) = $(length(seq)), not $ncols!")
        end

        # For each index of the current sequence,

        for j in 1:length(seq)

           # Check which character is at the current index and count it.

            @inbounds L = uppercase(seq[j])
            for k in 1:length(alphabet)
                @inbounds c = alphabet[k]
                if L == c
                    counts[j, k] += 1
                end
            end
        end
        progress("loading msa: read $i sequences .. ")
    end
    done()

    # Return the array of counts.

    counts
end

#---------------------------------rategrid-------------------------------------#

# Define a function to generate a rate grid of a given size and error threshold.

function rategrid(n::Integer, error_threshold::Float64)
    update("generating rate grid .. ")

    # Generate a vector of rates increasing on a log scale.

    if n == 10
        rates = [
            0.0001, 0.0005, 0.001, 0.002, 0.004, 0.008 , 0.01, 0.02, 0.05, 0.25
        ]
    else
        rates = logspace(-4,log10(0.1),n)
    end
    #print (rates)

    # Initialize uniques, an empty two-dimensional array. Uniques is initialized
    # as a set to avoid duplicate rows.

    uniques = Set{Vector{Float64}}()

    # Build uniques using the error threshold.

    M = 1. - 3. * error_threshold
	union!(uniques, [[M,error_threshold,error_threshold,error_threshold]])
	union!(uniques, [[error_threshold,M,error_threshold,error_threshold]])
	union!(uniques, [[error_threshold,error_threshold,M,error_threshold]])
	union!(uniques, [[error_threshold,error_threshold,error_threshold,M]])

    # For each value in the rates vector,

    for a in rates

        # Add a row to the uniques array.

        M = 1. - a - 2. * error_threshold

	    for i1 in 1:4
	        for i2 in 1:4
	            if i1 != i2
                    rate_vec = fill(error_threshold, (4))
	                rate_vec[i1] = M
	                rate_vec[i2] = a
	                union!(uniques, Any[rate_vec])
	            end
	        end
	    end

	    # Add four more rows to the uniques array.

        for b in rates
            for c in rates
                M = 1. - (a + b + c)
		        union!(uniques, [[M,a,b,c]])
		        union!(uniques, [[a,M,b,c]])
                union!(uniques, [[a,b,M,c]])
                union!(uniques, [[a,b,c,M]])

            end
        end
    end
    done()

    # Initialize a counter i.

    i = 1

    # Initialize an empty rate grid.

    rg = Array{Float64}(length(uniques),4)

    # For each row in uniques,

    for r in uniques
	    #println (r)

	    # Set the present row of the rate grid equal to the present row of
	    # uniques.

	    rg[i, :] = r

	    # Advance the counter i.

	    i += 1
    end

    # Print the length of uniques and return the rate grid.

    println(length(uniques))
    #print (uniques)
    return rg
end

#----------------------------------lfact---------------------------------------#

# Initialize the lfact_cache

lfact_cache = Float64[0.0]

# Define a function to return the lfact (logarithm of factorial) for a given
# integer.

function lfact(N::Integer)

    # If the input integer is 0, return 0.

    if N == 0
        return 0.0
    end

    # Otherwise, note the current length of the lfact_cache and check if it is
    # less than the input integer.

    n = length(lfact_cache)
    if n < N

        # If so, extend the lfact_cache to be of length equal to the input
        # integer and fill the new space with logarithmically-growing values.

        resize!(lfact_cache, N)
        for i in (n + 1):N
            @inbounds lfact_cache[i] = lfact_cache[i - 1] + log(i)
        end
    end

    # Return the largest value of the lfact_cache.

    @inbounds return lfact_cache[N]
end

#-----------------------------------lmc----------------------------------------#

# Define a function to return the lmc value for a given site.

function lmc(counts::Array{Int64,2}, site::Int64, nchars::Int64)

    # Initialize placeholders.

    c = 0.0
    s = 0

    # For each character,

    for char in 1:nchars
        @inbounds c -= lfact(counts[site, char])
        @inbounds s += counts[site, char]
    end
    return c + lfact(s)
end

#--------------------------------gridscores------------------------------------#

# Approximate negative infinity.

const smin = log(realmin(Float64))

# Define a function to compute grid scores.

function gridscores(counts::Array{Int64,2}, rates::Array{Float64,2}, alphabet::String)
    update("computing grid scores .. ")
    nchars = length(alphabet)
    npoints, _ = size(rates)
    nsites, _ = size(counts)

    # Initialize arrays for conditionals and scalers.

    conditionals = Array{Float64}(npoints, nsites)
    scalers = Array{Float64}(nsites,)

    # Compute log (rates)

    log_rates = log.(rates)

    # For each site of the input .msa,

    for i in 1:nsites

        # Set m to approximately negative infinity

        m = -realmax(Float64)

        # Compute lmc for the current site

        @inbounds c = lmc(counts, i, nchars)

        # For each rate,

        for j in 1:npoints
           s = c

           # for each character,

           for k in 1:nchars
                @inbounds s += counts[i, k] * log_rates[j, k]
            end
            if s > m
                m = s
            end

            # Update the conditionals.

            @inbounds conditionals[j, i] = s
        end

        # Update the scalers.

        @inbounds scalers[i] = m

        # For each rate,

        for j in 1:npoints
            @inbounds s = conditionals[j, i] - m
            if s < smin
                @inbounds conditionals[j, i] = 0
            else
                @inbounds conditionals[j, i] = exp(s)
            end
            #println (conditionals [j, i])
            #if (i == 771) && conditionals[j, i] > 0
            #    println (conditionals[j, i], "\t", join (rates[j , :], "\t"))
            #end
        end

    end
    done()

    # Return the conditionals and scalers.

    return (conditionals, scalers)
end

#---------------------------------ldensity-------------------------------------#

# Define a function to compute ldensity.

function ldensity(weights::Array{Float64,2}, alpha::Float64)

    # Set bounds on the weights array and on alpha to protect against extreme
    # cases.

    if minimum(weights) <= 1e-10
        return 1e-10
    end
    if alpha == 1
        return 0.0
    end

    # Compute the density.

    n = length(weights)
    r = lgamma(n * alpha) - n * lgamma(alpha)
    for i in 1:n
        @inbounds r += log(weights[i]) * (alpha - 1)
    end

    # Return the value of the density given weights and alpha.

    return r
end

#------------------------------------jll---------------------------------------#

# Define a function to compute ll and dir.

function jll(conditionals::Array{Float64,2}, scalers::Array{Float64,1}, weights::Array{Float64,2}, alpha::Float64)
    ll = sum(log.(weights * conditionals) + scalers')
    dir = ldensity(weights, alpha)
    return (ll, dir)
end

#-----------------------------------mcmc---------------------------------------#

# Define a function to carry out the mcmc process. Alpha was once called
# concentration_parameter.

function mcmc(
        conditionals::Array{Float64,2},
        scalers::Array{Float64,1};
        chain_length::Int64=10_000_000,
        burnin_fraction::Float64=0.5,
        expected_nsamples::Int64=100,
        alpha::Float64=0.5)

    # Compute nburnin from the inputs.

    nburnin = round(chain_length * burnin_fraction)

    # Note the shape of conditionals.

    npoints, nsites = size(conditionals)

    # Normalize the conditionals and weights by site.

    normalized_by_site = ones((1, npoints)) * conditionals

    normalized_weights = conditionals * ((1 ./ normalized_by_site) .* eye(nsites))

    # Sum the weights by site and normalize.

    sum_by_site = normalized_weights * ones((nsites, 1))
    weights = transpose((rand((npoints, 1)) .* (1.2 * sum_by_site - 0.8 * sum_by_site)) + 0.8 * sum_by_site)
    weights /= sum(weights)

    # Determine step size.

    stepsize = max(2. / max(npoints, nsites), median(weights))

    # Determine the number of samples and initialize arrays for their weights
    # and likelihoods.

    nsample = div(chain_length - nburnin, expected_nsamples)
    sampled_weights = zeros(Float64, (expected_nsamples, npoints))
    sampled_likelihoods = zeros(Float64, (expected_nsamples,))

    # Compute likelihoods, logsum, and jll for the current site.

    current_site_likelihoods = weights * conditionals
    current_site_logsum = sum(log.(current_site_likelihoods))
    current_ll, current_dir = jll(conditionals, scalers, weights, alpha)

    # Initialize step trackers.

    step = 1
    accepted_steps = 0
    mean_sampled_ll = 0
    sample_idx = 0

    # Initialize some other objects.

    diffvec = Array{Float64}(1, nsites)
    idxs = 1:npoints
    llstr = "(burning in)"

    tic()

    # For each step of the mcmc simulation,

    for step in 1:chain_length

        # Choose two random but separate indices.

        i = rand(idxs)
        j = rand(idxs)
        while i == j
            j = rand(idxs)
        end

        # Choose a random number proportional to stepsize and choose two
        # of the weights.

        change = rand() * stepsize
        @inbounds weights_i = weights[i]
        @inbounds weights_j = weights[j]

        # If weights_i exceeds change,

        if weights_i > change

            # Initialize lldiff

            lldiff = 0

            # For each site,

            for k in 1:nsites

                # Add a corresponding entry to diffvec and update lldiff.

                @inbounds diffvec[k] = (conditionals[j, k] - conditionals[i, k]) * change
                @inbounds lldiff += log(current_site_likelihoods[k] + diffvec[k])
            end

            # Normalize lldiff, generate dirdiff, and infer movecost.

            lldiff -= current_site_logsum
            @inbounds dirdiff = (alpha - 1) * (log((weights_i - change) / weights_i) + log((weights_j + change) / weights_j))
            movecost = exp(lldiff + dirdiff)

            # Use movecost to determine the probability of changing state.

            if rand() <= movecost

                # Update values and weights for the current site.

                current_ll += lldiff
                current_dir += dirdiff

                current_site_likelihoods += diffvec
                current_site_logsum += lldiff

                @inbounds weights[i] -= change
                @inbounds weights[j] += change

                # Advance the accepted_steps counter.

                accepted_steps += 1
            end
        end

        # If this is not a burn-in step,

        if step > nburnin

            # If a sampling interval has elapsed,

            if (step - nburnin) % nsample == 0

                # Advance the sample index counter, and sample the weights and
                # likelihoods from this simulation step.

                sample_idx += 1
                sampled_weights[sample_idx, :] = weights
                sampled_likelihoods[sample_idx] = current_ll
                mean_sampled_ll += current_ll
                llstr = "$(trunc(mean_sampled_ll / sample_idx, 3))"
            end
        end

        # If a sampling interval has elapsed,

        if step % nsample == 0

            # Report progress.

            progress("running MCMC chain: step = $step/$chain_length, mean logL = $llstr, acceptance rate = $(trunc(accepted_steps / step, 3)) .. ")
        end
    end

    done()

    # Return all of the sampled likelihoods and weights.

    return (sampled_likelihoods, sampled_weights)
end

#---------------------------------loadrates------------------------------------#

# Define a function to load the rates file.

function loadrates()
    rates3 = readdlm("data/rates.txt", ' ', Float64)
    nrates, _ = size(rates3)
    rates = Array(Float64, (nrates, 4))
    rates[:, 1:3] = rates3
    for i in 1:nrates
        @inbounds rates[i, 4] = 1.0 - sum(rates[i, 1:3])
    end
    return rates
end

#---------------------------------postproc-------------------------------------#

# Define a function to compute prior and posterior probabilities for the
# conditionals and ws arrays.

function postproc(conditionals::Array{Float64,2}, ws::Array{Float64,2})
    update("computing posterior probabilities .. ")

    # Record the shapes of conditionals and ws.

    nsites, npoints = size(conditionals)
    nsamples, _ = size(ws)

    # Initialize an array to store the priors.

    priors = zeros(Float64, (1, npoints))

    # Deduce the priors from ws and normalize.

    for i in 1:nsamples
        @inbounds priors += reshape(ws[i, :], 1, npoints)
    end
    priors = priors' / nsamples
    priors /= sum(priors)

    # Compute the normalization constant for the posteriors.

    normalization = conditionals * priors

    # Compute the posteriors.

    posteriors = ((1 ./ normalization) .* eye(nsites)) * (conditionals * (priors .* eye(npoints)))
    done()

    # Return the priors and posteriors.

    return (priors, posteriors)
end

#-------------------------------callvariants-----------------------------------#

# Define a function to make variant calls.

function callvariants(
        counts::Array{Int64,2},
        rates::Array{Float64,2},
        priors::Array{Float64,2},
        posteriors::Array{Float64,2},
        alphabet::String;
        target_rate::Float64=0.1,
        posterior_threshold::Float64=0.999,
        json_file=nothing,
        min_coverage::Int64=100)

    update("calling variants .. ")

    # Record the number of posteriors and filter them according to the target
    # rate.

    nsites, npoints = size(posteriors)
    posterior_prob = posteriors * map(x -> x > target_rate, rates)
    majority_prob  = posteriors * map(x -> x >= 0.5, rates)

    # If json_file already exists,

    if json_file != nothing

        # Initialize the overall report and rate dictionaries.

        overall_report = Dict()
        rate_dict = Dict{Int64, Dict{String,Float64}}()

        # For each posterior,

        for i in 1:npoints

            # If the corresponding prior is great enough,

            if priors[i] >= 0.5/nsites

                # Add the rates and prior at this index to the rate dictionary.

                rate_dict[i] = Dict()
                for N in 1:length(alphabet)
                    rate_dict[i][string(alphabet[N])] = rates[i, N]
                end
                rate_dict[i]["weight"] = priors[i]
            end
        end

        # Record some more information.

        posterior_means = posteriors * rates
        posterior_dict = Dict{Int64, Dict{String,Any}}()
        overall_report["priors"] = rate_dict
        consensus = ""
    end

    # Initialize the variants dictionary.

    variants = Dict{Int64,String}()

    # For each site,

    for i in 1:nsites
        s = Set{Char}()

        # Record the majority pair.

        maj_pair = findmax(majority_prob[i, :])

        # If the majority posterior is large enough,

        if maj_pair[1] >= posterior_threshold

            # Add the majority prior to s.

            push!(s, alphabet[maj_pair[2]])

            # For each character in the alphabet,

            for j in 1:length(alphabet)

                # If the posterior at the current character is large enough,

                if posterior_prob[i, j] >= posterior_threshold

                    # Add the current character to s.

                    push!(s, alphabet[j])
                end
            end
        end

        # Organize s, and assign it to the variants at this index.

        s_ = join(sort!(collect(s)))
        #if length(s) > 1
        println("$i: $s_,\t[$(join(counts[i, :], '\t'))]")
        #end
        variants[i] = s_

        # If json_file already exists,

        if json_file != nothing

            # If coverage at this index is above the minimum threshold,

            if sum(counts[i,:]) >= min_coverage

                # Record the index of the maximum posterior probability and
                # update consensus.

                index = indmax(posterior_prob[i,:])
                consensus = *(consensus,string(alphabet[index]))
            else

                # Apply a null update to consensus.

                consensus = *(consensus,"-")
            end

            # Update the posterior dictionary.
            posterior_dict[i] = Dict()
            for N in 1:length(alphabet)
                posterior_dict[i][string(alphabet[N])] = [posterior_prob[i,N], posterior_means[i,N], counts[i,N]]
            end
            #posterior_dict[i] = (String=>Any)[
            #    string(alphabet[N]) => [posterior_prob[i,N], posterior_means[i,N], counts[i,N]] for N = 1:length(alphabet) ]
            posterior_dict[i]["variants"] = s_
        end
    end
    done()

    # If json_file already exists,

    if json_file != nothing

        # update the overall report and write it to json_file.

        overall_report["consensus"] = consensus
        overall_report["posteriors"] = posterior_dict
        json_out = open(json_file, "w")
        JSON.print(json_out, overall_report)
        close(json_out)
    end

    # Return the variants.

    return variants
end

#--------------------------------filterseqs------------------------------------#

# Define a function to filter the sequences.

function filterseqs(msa::String, variants::Dict{Int64,String}, dest::String, alphabet::String)
    tic()

    # Extract a sites list from the variants.

    sites = sort!(collect(keys(variants)))

    # Initialize varmask.

    varmask = zeros(Bool, ((length(variants), 4)))

    # For each site,

    for i in 1:length(sites)
        site = sites[i]

        # For each character in the alphabet,

        for j in 1:length(alphabet)

            # If the current character is a variant at the current site,

            if alphabet[j] in variants[site]

               # Record it in varmask.

                varmask[i, j] = true
            end
        end
    end

    # Open a fasta file.

    FastaWriter(dest == "-" ? STDOUT : dest) do fw

        # For each sequence in the input .msa file,

        for (i, (name, seq)) in enumerate(FastaReader(msa))

            # Initialize the counter v,

            v = 1
            seq_ = Array{Char}(length(seq),)

            # For each index in the current sequence,

            for j in 1:length(seq)
                seq_[j] = seq[j]

                # If the sequence has not ended, and the counter v is equal
                # to the current indes,

                if v < length(sites) && j == sites[v]

                    # Extract the alphabetical index of the character at the
                    # current index.

                    idx = index(alphabet, seq[j])

                    # If the alphabetical index is positive and this character
                    # is not a variant at this site, but the variant at this
                    # site is ambiguous,

                    if idx > 0 && !varmask[v, idx] && variants[sites[v]] in keys(AMBIGS)

                        # Assign the appropriate AMBIGS element to the current
                        # index.

                        seq_[j] = AMBIGS[variants[sites[v]]]
                    end

                    # Advance the counter v.

                    v += 1
                end
            end

            # Report progress.

            writeentry(fw, name, join(seq_, ""))
            progress("filtering unsupported variants: processed $i sequences .. ")
        end
    end
    done()
end

#------------------------------------------------------------------------------#
#                         completion announcement                              #
#------------------------------------------------------------------------------#

# Report

function done()
    print_with_color(:blue, STDERR, "done, took $(toq())s")
    println(STDERR)
end
function progress(args...)
    print_with_color(:blue, STDERR, "\rINFO: ", args...)
end
function update(args...)
    print_with_color(:blue, STDERR, "INFO: ", args...)
    tic()
end

#------------------------------------------------------------------------------#
#                                   main                                       #
#------------------------------------------------------------------------------#

# Define the main function.

function main(args)

    # Parse arguments.

    s = ArgParseSettings()
    s.description = "call variants using a multinomial model sampled by MCMC"

    @add_arg_table s begin
        "--grid-density", "-g"
            arg_type = Int64
            default = 10
        "--chain-length", "-c"
            arg_type = Int64
            default = 2_000_000
        "--burnin-fraction", "-b"
            arg_type = Float64
            default = 0.5
        "--target-rate", "-t"
            arg_type = Float64
            default = 0.01
        "--posterior-threshold", "-p"
            arg_type = Float64
            default = 0.95
        "--filter", "-f"
            arg_type = String
        "--json-report", "-j"
            arg_type = String
        "msa"
            required = true
    end

    parsed_args = parse_args(args, s)

    # Construct the rate grid.

    rates = rategrid(parsed_args["grid-density"], parsed_args["target-rate"])

    # Construct the alphabet.

    const alphabet = "ACGT"

    # Execute countmsa() on the input .msa file and the alphabet.

    counts = countmsa(parsed_args["msa"], alphabet)

    # Construct conditionals and scalers using gridscores().

    conditionals, scalers = gridscores(counts, rates, alphabet)

    # Obtain lls and ws via MCMC simulation.

    lls, ws = mcmc(
        conditionals, scalers,
        chain_length=parsed_args["chain-length"], burnin_fraction=parsed_args["burnin-fraction"])

    # Extract priors and posteriors from conditionals and ws using postproc.

    priors, posteriors = postproc(conditionals', ws)

    # Call variants.

    variants = callvariants(
        counts, rates, priors, posteriors, alphabet,
        target_rate=parsed_args["target-rate"], posterior_threshold=parsed_args["posterior-threshold"], json_file = parsed_args["json-report"])

    # If instructed to, filter the sequences.

    if "filter" in keys(parsed_args)
        filterseqs(parsed_args["msa"], variants, parsed_args["filter"], alphabet)
    end
end

# Execute.

main(ARGS)

max_stars_repo_name: stevenweaver/mmvc

max_issues_repo_name: stevenweaver/mmvc

max_forks_repo_name: stevenweaver/mmvc


text:
@recipe function f(::Type{Val{:tsaplot}}, x, y, z)

    N = length(x)

    grid := false
    legend := false
    linewidth --> 1

    @series begin
        seriestype := :stem
        seriescolor --> :black
        x := x
        y := y
    end

    bound = 1.96 / √N
    linestyle := :dash
    linecolor --> :blue

    @series begin
        seriestype := :hline
        y := -bound
    end

    @series begin
        seriestype := :hline
        y := bound
    end

    @series begin
        seriestype := :hline
        linestyle := :solid
        linecolor := :black
        linewidth := 1
        y := 0.0
    end

end

@userplot Correlogram
@recipe function f(c::Correlogram)
    if length(c.args) != 2 || !(typeof(c.args[1]) <: AbstractVector) ||
        !(typeof(c.args[2]) <: Integer)
        error("Correlograms should be given a vector and an integer.  Got: $(typeof(h.args))")
    end
    z, k = c.args
    seriestype := :tsaplot
    1:k, map(i -> autocorrelation(z, i), 1:k)
end

@userplot Partial_Correlogram
@recipe function f(c::Partial_Correlogram; recursive=true)
    if length(c.args) != 2 || !(typeof(c.args[1]) <: AbstractVector) ||
        !(typeof(c.args[2]) <: Integer)
        error("Correlograms should be given a vector and an integer.  Got: $(typeof(h.args))")
    end
    z, k = c.args
    seriestype := :tsaplot

    if recursive
        ρ = map(i -> autocorrelation(z, i), 1:k)
        T = typeof(zero(eltype(z)) / 1)
        ϕ = zeros(T, k)
        ϕpp = zeros(T, k)
        σ2 = Ref{T}(autocovariance(z, 0))
        for i in 1:k
            levinson_durbin!(ϕ, σ2, ρ, i)
            ϕpp[i] = ϕ[i]
        end
    else
        ϕpp = map(i -> partial_autocorrelation(z, i), 1:k)
    end

    1:k, ϕpp
end

ext: jl


text:
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

"""
The raw transaction format and several peer-to-peer network messages use a type
of variable-length integer to indicate the number of bytes in a following piece
of data.
Bitcoin Core code refers to these variable length integers as compactSize.
Many other documents refer to them as var_int or varInt, but this risks
conflation with other variable-length integer encodings—such as the CVarInt
class used in Bitcoin Core for serializing data to disk. Because it’s used in
the transaction format, the format of compactSize unsigned integers is part of
the consensus rules.
For numbers from 0 to 252, compactSize unsigned integers look like regular
unsigned integers. For other numbers up to 0xffffffffffffffff, a byte is
prefixed to the number to indicate its length—but otherwise the numbers look
like regular unsigned integers in little-endian order.

Value                                    | Bytes Used  | Format
-----------------------------------------|-------------|-----------------------------------------
>= 0 && <= 252                           | 1           | uint8_t
>= 253 && <= 0xffff                      | 3           | 0xfd followed by the number as uint16_t
>= 0x10000 && <= 0xffffffff              | 5           | 0xfe followed by the number as uint32_t
>= 0x100000000 && <= 0xffffffffffffffff  | 9           | 0xff followed by the number as uint64_t

For example, the number 515 is encoded as 0xfd0302.
"""
struct CompactSizeUInt <: Unsigned
    value::Unsigned
    prefix::Union{Nothing, UInt8}
end

CompactSizeUInt(n::Integer) =
    n <  0x00               ? error("Integer is negative!") :
    n <= 0xfc               ? CompactSizeUInt(UInt8(n), nothing) :
    n <= 0xffff             ? CompactSizeUInt(UInt16(n), 0xfd) :
    n <= 0xffffffff         ? CompactSizeUInt(UInt32(n), 0xfe) :
    n <= 0xffffffffffffffff ? CompactSizeUInt(UInt64(n), 0xff) :
    error("Integer is too large!")

function Base.show(io::IO, z::CompactSizeUInt)
    print(io, Int(z.value))
end

String(z::CompactSizeUInt) = string(Int(z.value))

"""
    read(io::IOBuffer)::CompactSizeUInt

Returns a CompactSizeUInt from an IO stream
"""
function CompactSizeUInt(io::IOBuffer)
    p = read(io, 1)[1]
    if p == 0xfd
        n = ltoh(reinterpret(UInt16, read(io, 2))[1])
    elseif p == 0xfe
        n = ltoh(reinterpret(UInt32, read(io, 4))[1])
    elseif p == 0xff
        n = ltoh(reinterpret(UInt64, read(io, 8))[1])
    else
        n = p
        p = nothing
    end
    CompactSizeUInt(n, p)
end

"""
    serialize(n::CompactSizeUInt) -> Vector{UInt8}

Returns the bytes serialization of a CompactSizeUInt
"""
function serialize(n::CompactSizeUInt)
    if n.prefix == nothing
        return [n.value]
    else
        return append!([n.prefix], reinterpret(UInt8, [htol(n.value)]))
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b3e63d70-68c6-11eb-17c2-0f899e780c9e
begin
	using PlutoUI
	using Test
end

# ╔═╡ 2d54bdea-68c6-11eb-3f0f-6978c66b6785
md"""
### Queen Placement problem
"""

# ╔═╡ 2142ac02-68c8-11eb-3c34-35667eb5e7d5
md"""
**Define a type point**.
"""

# ╔═╡ 559b845a-68c6-11eb-37c8-7798fbf254aa
begin
  struct Point{T <: Integer}
  	x::T
  	y::T
  end

  Base.show(io::IO, p::Point) = print(io, "($(p.x), $(p.y))")

  const VP{T} = Vector{Point{T}} where {T <: Integer}
end

# ╔═╡ 6aa8bfde-68c6-11eb-224b-1dd95917caf7
"""
Add a queen at position p to the current board
"""
qadd(board::VP, p::Point) = push!(copy(board), p)

# ╔═╡ 70f930e4-68c6-11eb-2881-5fcbfb3358f1
begin
	Base.:-(p₁::Point, p₂::Point) = (p₁.x - p₂.x, p₁.y - p₂.y)
	Base.abs(p::Point) = (abs(p.x), abs(p.y))

	## Iterator
	Base.collect(p::Point) = [p.x, p.y]

	function Base.iterate(iter::Point, state=(collect(iter), 1))
  		comp, ix = state
  		comp === nothing && return nothing     ## are we done?
  		ix == length(comp) && return (comp[ix], (nothing, ix))
  		(comp[ix], (collect(iter), ix + 1))
	end

	Base.length(iter::Point) where T = length(collect(iter))
	Base.eltype(::Point{T}) where T = T
end

# ╔═╡ cc5b11aa-68c6-11eb-25af-35517b6dd119
p = Point(2, 3)

# ╔═╡ 997cc7b0-68c6-11eb-2a4d-715782d7cc30
with_terminal() do
	for c ∈ p; println(c); end
end

# ╔═╡ e8a4563e-68c6-11eb-2df6-236cdf425e22
begin
	"""
	Check for constriant satisfaction:
	no 2 queens using the same row nor the same columns nor the same diagonals
	"""
	qhit(board::VP, p::Point) = any(map(o_p -> qhit_hlpr(p, o_p), board))

	qhit_hlpr(p₁::Point, p₂::Point) =
    	any(p₁ .== p₂) || abs((p₁ - p₂)[1]) == abs((p₁ - p₂)[2])
	#   check for ≠ row/col or ≠ diags
end

# ╔═╡ f39f8fdc-68c6-11eb-0521-6b474172d441
"""
Try placing the n queens on the n × n board
"""
function queen_solver(n::Int; m::Int=n, board=VP{Int}())
    m < 1 && return board
    for x = 1:n, y = 1:n
      p = Point(x, y)
      if ! qhit(board, p)
          ## then place queen at position Point(x, y)
          s = queen_solver(n; m=m-1, board=qadd(board, p))
          s != nothing && return s
      end
    end

    return nothing
end

# ╔═╡ fd0f3c68-68c6-11eb-1285-41862ff156d1
methods(queen_solver)

# ╔═╡ 0dea7ef8-68c7-11eb-0499-1b1c2e83f31b
with_terminal() do
	n = 8   # for 8 queens
	sol = queen_solver(n)
	println("solution for n: $(n)")
	if sol == nothing
		println("No solution")
	else
		for c ∈ sol; println(c) ; end
	end

end

# ╔═╡ d09cf0ca-68c7-11eb-3bc3-6dadca56a7e0
with_terminal() do
	n = 4   # for 4 queens
	sol = queen_solver(n)
	println("solution for n: $(n)")
	if sol == nothing
		println("No solution")
	else
		for c ∈ sol; println(c) ; end
	end
end

# ╔═╡ be45e20a-68c8-11eb-1281-55f090b4d5bb
with_terminal() do
	n = 3   # for 3 queens
	sol = queen_solver(n)
	println("solution for n: $(n)")
	if sol == nothing
		println("No solution")
	else
		for c ∈ sol; println(c) ; end
	end
end

# ╔═╡ Cell order:
# ╟─2d54bdea-68c6-11eb-3f0f-6978c66b6785
# ╠═b3e63d70-68c6-11eb-17c2-0f899e780c9e
# ╟─2142ac02-68c8-11eb-3c34-35667eb5e7d5
# ╠═559b845a-68c6-11eb-37c8-7798fbf254aa
# ╠═6aa8bfde-68c6-11eb-224b-1dd95917caf7
# ╠═70f930e4-68c6-11eb-2881-5fcbfb3358f1
# ╠═cc5b11aa-68c6-11eb-25af-35517b6dd119
# ╠═997cc7b0-68c6-11eb-2a4d-715782d7cc30
# ╠═e8a4563e-68c6-11eb-2df6-236cdf425e22
# ╠═f39f8fdc-68c6-11eb-0521-6b474172d441
# ╠═fd0f3c68-68c6-11eb-1285-41862ff156d1
# ╠═0dea7ef8-68c7-11eb-0499-1b1c2e83f31b
# ╠═d09cf0ca-68c7-11eb-3bc3-6dadca56a7e0
# ╠═be45e20a-68c8-11eb-1281-55f090b4d5bb

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using Genie.Router
using Genie.Renderer.Html

route("/") do
  html("Hello World")
end

max_stars_repo_name: QuestofIranon/Hello-World-Julia

max_issues_repo_name: QuestofIranon/Hello-World-Julia

max_forks_repo_name: QuestofIranon/Hello-World-Julia


text:
using Test, LinearMaps, LinearAlgebra, BenchmarkTools

@testset "identity/scaling map" begin
    @test_throws ArgumentError LinearMaps.UniformScalingMap(true, -1)
    A = 2 * rand(ComplexF64, (10, 10)) .- 1
    B = rand(size(A)...)
    M = @inferred 1 * LinearMap(A)
    N = @inferred LinearMap(B)
    LC = @inferred M + N
    v = rand(ComplexF64, 10)
    w = similar(v)
    Id = @inferred LinearMap(I, 10)
    @test occursin("10×10 LinearMaps.UniformScalingMap{Bool}", sprint((t, s) -> show(t, "text/plain", s), Id))
    @test_throws ErrorException LinearMaps.UniformScalingMap(1, 10, 20)
    @test_throws ErrorException LinearMaps.UniformScalingMap(1, (10, 20))
    @test size(Id) == (10, 10)
    @test @inferred isreal(Id)
    @test @inferred issymmetric(Id)
    @test @inferred ishermitian(Id)
    @test @inferred isposdef(Id)
    @test Id * v == v
    @test (2 * M' + 3 * I) * v == 2 * A'v + 3v
    @test (3 * I + 2 * M') * v == 2 * A'v + 3v
    @test (2 * M' - 3 * I) * v == 2 * A'v - 3v
    @test (3 * I - 2 * M') * v == -2 * A'v + 3v
    @test (3 * I - 2 * M') * v == -2 * A'v + 3v
    @test transpose(LinearMap(2 * M' + 3 * I)) * v ≈ transpose(2 * A' + 3 * I) * v
    @test LinearMap(2 * M' + 0I)' * v ≈ (2 * A')' * v
    for λ in (0, 1, rand()), α in (0, 1, rand()), β in (0, 1, rand()), sz in (10, (10,5))
        Λ = @inferred LinearMap(λ*I, 10)
        x = rand(Float64, sz)
        y = rand(Float64, sz)
        b = @benchmarkable mul!($y, $Λ, $x, $α, $β)
        @test run(b, samples=3).allocs == 0
        y = deepcopy(x)
        @inferred mul!(y, Λ, x, α, β)
        @test y ≈ λ * x * α + β * x
    end
    for elty in (Float64, ComplexF64), transform in (identity, transpose, adjoint)
        λ = rand(elty)
        x = rand(10)
        J = @inferred LinearMap(LinearMaps.UniformScalingMap(λ, 10))
        @test transform(J) * x == transform(λ) * x
        J = @inferred LinearMap(λ*I, 10)
        @test (λ * J) * x == (J * λ) * x == (λ * λ) * x
    end
    X = rand(10, 10); Y = similar(X)
    @test mul!(Y, Id, X) == X
    @test Id*X isa LinearMap
    @test X*Id isa LinearMap
    @test Matrix(Id*X) == X
    @test Matrix(X*Id) == X
end


text:
@everywhere begin
    check_option(option, value, default) = Dagger.get_options(option, default) == value
    function inc_by_myid(x)
        m = Dagger.@spawn myid()
        fetch(Dagger.@spawn x + m)
    end
    switch_option(option, value, default) = Dagger.with_options(NamedTuple([option=>value])) do
        fetch(Dagger.@spawn Dagger.get_options(option, default))
    end

    struct SpecialFunc
        x::Int
    end
    function (sf::SpecialFunc)(obj)
        if obj isa Dagger.Chunk
            return fetch(obj)+sf.x
        else
            return sf.x
        end
    end
    Dagger.move(from::Union{OSProc,Dagger.ThreadProc}, to::Dagger.ThreadProc, sf::SpecialFunc) =
        (to.owner == 1 && to.tid == 1) ? SpecialFunc(1) : sf
end

@testset "Scope propagation" begin
    first_wid = first(workers())
    last_wid = last(workers())
    for (option, default, value, value2) in [
        # Special handling
        (:scope, AnyScope(), ProcessScope(first_wid), ProcessScope(last_wid)),
        (:processor, OSProc(), Dagger.ThreadProc(first_wid, 1), Dagger.ThreadProc(last_wid, 1)),
        # ThunkOptions field
        (:single, 0, first_wid, last_wid),
        # Thunk field
        (:meta, false, true, false)
    ]
        # Test local and remote default values
        @test Dagger.get_options(option, default) == default
        @test fetch(Dagger.@spawn check_option(option, default, default))

        # Test local propagation
        Dagger.with_options(NamedTuple([option=>value])) do
            fetch(@async @test Dagger.get_options(option, default) == value)
            fetch(Threads.@spawn @test Dagger.get_options(option, default) == value)
        end

        # Test remote option switching
        Dagger.with_options(NamedTuple([option=>value])) do
            @test fetch(Dagger.@spawn switch_option(option, value2, default)) == value2
        end
    end

    # Test scope/single is applied
    for wid in workers()
        for (option, value) in [
            (:scope, ProcessScope(wid)),
            (:single, wid)
        ]
            Dagger.with_options(NamedTuple([option=>value])) do
                @test fetch(Dagger.@spawn inc_by_myid(1)) == 1+wid
                @test fetch(Dagger.@spawn inc_by_myid(2)) == 2+wid
            end
        end
    end

    # Test processor/meta is applied
    sf = SpecialFunc(0)
    obj = Dagger.tochunk(42)
    Dagger.with_options(;scope=Dagger.ExactScope(Dagger.ThreadProc(first_wid,1))) do
        @test fetch(Dagger.@spawn sf(obj)) == 0
        @test fetch(Dagger.@spawn sf(obj)) == 0
    end
    Dagger.with_options(;scope=Dagger.ExactScope(Dagger.ThreadProc(1,1)), processor=Dagger.ThreadProc(1,1), meta=true) do
        @test fetch(Dagger.@spawn sf(obj)) == 43
        @test fetch(Dagger.@spawn sf(obj)) == 43
    end
end


text:
module Assignment06

export count_records,
       fasta_header,
       parse_fasta

include("assignment.jl")

end


text:
"""
    LayerSetting

The current settings of an [`Layer`](@ref) which are saved in `layer.current_setting`.

# Fields
- `opacity::Float64`: the current opacity
- `current_scale::Tuple{Float64, Float64}`: the current scale
- `rotation_angle::Float64`: the angle of rotation of a layer.
- `misc::Dict{Symbol, Any}`: any other misc metadata
"""
mutable struct LayerSetting
    opacity::Float64
    scale::Scale
    rotation_angle::Float64
    misc::Dict{Symbol,Any}

    LayerSetting() = new(1.0, Scale(1.0, 1.0), 0.0, Dict())
end


text:
using Test
using TestSetExtensions
using FluxTraining
using FluxTraining: EpochEnd, LR, getdataloader, getoptimparam, setoptimparam!, protect, Protected, ProtectedException
using Flux: trainable

include("./testdata.jl")


text:
using Test

using MPI

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

# Not possible to PROD a Char (and neither Int8 with OpenMPI)
typs = setdiff([(Base.uniontypes(MPI.MPIDatatype))...], [Char, Int8, UInt8])
for typ in typs
    val = convert(typ,rank + 1)
    B = MPI.Scan(val, MPI.PROD, comm)
    @test B[1] ≈ convert(typ, prod(1:(rank+1)))
end

MPI.Finalize()


text:
const TEST_NAME = "AESNI1x128"

include("common.jl")

using RNG.Random123

r = AESNI1x(123)

test_all(r, 100_000_000)


text:
module Enums
export SectionType,MaterialType,LoadCaseType

@enum SparseSolver begin
    USE_PARDISO
end

@enum SectionType begin
    GENERAL_SECTION=0
    ISECTION=1
    HSECTION=2
    BOX=3
    PIPE=4
    CIRCLE=5
    RECTANGLE=6
end

@enum LoadCaseType begin
    STATIC
    MODAL
    BUCKLING
    TIME_HISTORY
    SPECTRUM
end

@enum MaterialType begin
    GENERAL_MATERIAL
    ISOELASTIC
    UNIAXIAL_METAL
end

@enum MassSource begin
    WEIGHT
    LOAD
end

@enum MassMatrixType begin
    CONCENTRATED
    COORDINATED
end

@enum DampType begin
    CONSTANT
    RAYLEIGH
end

@enum DynamicAlgorithm begin
    CENTRAL_DIFF
    NEWMARK_THETA
    WILSON_BETA
    MODOL_DECOMP
end

@enum ModalType begin
    EIGEN
    RITZ
end

@enum RefCSys begin
    GLOBAL
    LOCAL
end

end


text:
# Basic library loading and API calling

import Base: get

export
    CUDA_VENDOR,
    @cucall

# Early library loading
# TODO: allow overriding the loaded library using ENV
function load_library()
    try
        # TODO: thorough check of vendor?
        return (dlopen("libcuda"), "NVIDIA")
    end

    try
        return (dlopen("libocelot"), "Ocelot")
    end

    error("Could not load CUDA (or any compatible) library")
end
const (libcuda, CUDA_VENDOR) = load_library()

# API call wrapper
macro cucall(f, argtypes, args...)
    quote
        api_function = resolve($f)
        status = ccall(dlsym(libcuda, api_function), Cint, $argtypes, $(args...))
        if status != 0
            err = CuError(Int(status))
            throw(err)
        end
    end
end

api_mapping = Dict{Symbol,Symbol}()
resolve(f::Symbol) = get(api_mapping, f, f)

function initialize_api()
    # Create mapping for versioned API calls
    if haskey(ENV, "CUDA_FORCE_API_VERSION")
        api_version = ENV["CUDA_FORCE_API_VERSION"]
    else
        api_version = driver_version()
    end
    global api_mapping
    if api_version >= 3020
        api_mapping[:cuDeviceTotalMem]   = :cuDeviceTotalMem_v2
        api_mapping[:cuCtxCreate]        = :cuCtxCreate_v2
        api_mapping[:cuMemAlloc]         = :cuMemAlloc_v2
        api_mapping[:cuMemcpyHtoD]       = :cuMemcpyHtoD_v2
        api_mapping[:cuMemcpyDtoH]       = :cuMemcpyDtoH_v2
        api_mapping[:cuMemFree]          = :cuMemFree_v2
        api_mapping[:cuModuleGetGlobal]  = :cuModuleGetGlobal_v2
        api_mapping[:cuMemsetD32]        = :cuMemsetD32_v2
    end
    if api_version >= 4000
        api_mapping[:cuCtxDestroy]       = :cuCtxDestroy_v2
        api_mapping[:cuCtxPushCurrent]   = :cuCtxPushCurrent_v2
        api_mapping[:cuCtxPopCurrent]    = :cuCtxPopCurrent_v2
    end

    # Initialize the driver
    @cucall(:cuInit, (Cint,), 0)
end

function driver_version()
    version_box = ptrbox(Cint)
    @cucall(:cuDriverGetVersion, (Ptr{Cint},), version_box)
    return ptrunbox(version_box)
end

# Box a variable into an array for ccall() passing
ptrbox(T::Type) = Array(T, 1)
ptrbox(T::Type, val) = T[val]
ptrunbox{T}(box::Array{T, 1}) = box[1]
ptrunbox{T}(box::Array{T, 1}, desttype::Type) = convert(desttype, ptrunbox(box))


text:
# Consider creating Source module in source.jl under source/ directory.  Such a module can
# be used by
#    include("source/source.jl")
#    importall (or using?) .Source
# See "importall .LinAlg" in julia-0.6/base/sysimg.jl.

# When we discretize a source, a source located in the interval that begins or ends with a
# ghost point is cumbersome.  In such a case, we have to assign current density at the ghost
# point, which is impossible because there is no entry in the current density array that
# is assigned to ghost points.
#
# Depending on the boundary condition, we have to handle such a case differently.  Suppose
# we are assiging a point source of the primal field type in the direction parallel to (but
# slightly away from) the boundary.  We assign a portion of the source current at the
# nonghost point adjacent to the ghost point.  The question is how to handle the portion of
# the source current to be assigned to the ghost point.
#
# - For Bloch, we have to assign a point source at the nonghost point corresponding to the
# ghost point.  Note that this should be additive to any existing current density at the
# nonghost point.  (The same situation as two point sources in neighboring grid intervals:
# at the shared grid point, the contributions from the two point sources are added.  The
# first and last grid intervals are adjacent through the Bloch boundary condition.)
#
# - For the symmetry boundary, putting current density component parallel to the boundary
# effectively puts another anti-symmetric current density component at the mirror symmetry
# position.  Usually these two current density components are sufficiently well-separated
# and their contributions do not overlap.  However, in the situation considered here, they
# are actually in neighboring cell intervals that share a point on the symmetry boundary.
#
# If a point source is in a cell interval between two grid points A and B, the method to
# assign current densities to A and B is the "sliding" method.  In this method, imagine the
# point source location changes gradually from A to B.  When it is exactly at A, current
# density is assigned only to A, and the magnitude J_A of the current density is determined
# by dividing the current of the point source by the grid cell area occupied by A.  We can
# similarly calculate the current density J_B to assign to B in the case where the point
# source is exactly at B.  (Note that J_A ≠ J_B in general even though the current of the
# point source is the same, because the area occupied by A and B can be different if the
# grid is nonuniform.)  Now, if the point source location changes gradually from A to B, the
# current density to assign to A should decrease linearly from J_A to 0, and the current
# density to assign to B should increase linearly from 0 to J_B.  Therefore, if the point
# source is located between A and B such that the ratio of the distances to A and B is
# x:(1-x), then the current densities to assign to A and B should be J_A * (1-x) and J_B * x,
# respectively.  Note that these becomes J_A and 0 for x = 0 (when the point source is at A)
# and 0 and J_B for x = 1 (when the point source is at B).
#
# The situation of the symmetry boundary condition discussed earlier, where the two point
# sources with opposite directions are placed at the mirror symmetry locations around the
# symmetry boundary, can be thought of as assigning the current densities for the two
# individual point sources using the above described recipe.  Note that as the two point
# sources approach to the symmetry boundary, the portions of the two point sources' currents
# assigned at the shared point on the boundary exactly cancel each other.  Therefore, we
# always assign zero current density on the ghost point, which is at the symmetry boundary.
# This makes us to simply ignore assiging current density to the point on the symmetry
# boundary.  The portion of the total current to assign to the grid point adjacent
# to the symmetry boundary is as usual, and this portion becomes zero as the point source
# approaches the symmetry boundary.  This makes sense, because if the point source is
# exactly on the symmetry boundary, being parallel to the boundary, then there must be no
# current in the system at all.
#
#
#
# Now, consider assignment of a point source corresponding to the dual field parallel to the
# boundary.  For Bloch, we proceed similarly.  For the symmetry boundary, now the image
# source put behind the boundary is in the same direction as the original source in the
# domain.  Also, when the source is sufficiently close to the boundary, the actual and
# image sources are in a single dual grid cell interval, whose  the middle point is the
# boundary point and the end points are dual grid points, which are referred to as A and B
# as before for convenience below, instead of the two adjacent cell intervals.  Suppose A is
# within the domain and B is outside the domain.  We don't have to assign a portion of the
# total current of the real point source to B ourselves, which is behind the symmetry
# boundary: this is something automatically done by imposing symmetry.  However, if the
# image point source existing behind the symmetry boundary assigns some portion of its total
# current to A, which is within the domain, that must be assigned ourselves.
#
# Because the two point sources are at mirror-symmetric locations, the reduction in the
# amount of the current assigned to A by the point source moving away from A towards the
# symmetry boundary is exactly compensated by the increase in the amount of the current
# assigned to A by the image point source to the symmetry boundary from behind the boundary.
# Therefore, the net effect is the constant amount of current assigend to the dual grid
# point, no matter where the point source is located inside this cell interval containing
# the symmetry boundary.

# Note that we need to support both Je and Jm to realize unidirectional sources.

export Source
export add_src!

abstract type Source{K,Kf} end  # K: shape dimension; Kf: field dimension

add_src!(jKd::AbsArrNumber, gt₀::AbsVec{GridType}, bounds::Tuple2{AbsVecReal},
         l::Tuple2{NTuple{K,AbsVecReal}}, ∆l::Tuple2{NTuple{K,AbsVecReal}},
         isbloch::AbsVecBool, srcs::Source...) where {K} =
    add_src!(jKd, SVec{K}(gt₀), (float.(l[nPR]),float.(l[nDL])), (float.(∆l[nPR]), float.(∆l[nDL])), SVec{K}(isbloch), srcs...)

# About the determination of gt_cmp
# In the concrete implementation of add_src!() for each concrete Source type, gt_cmp, which
# is the grid type of the J field component to set up (typically indicated by the
# nw-component) is determined by
#
#     gt_cmp = src.isfield˔shp ? gt₀ : gt_w(nw, gt₀)
#
# gt_cmp is used to obtain the locations of the field by
#
#     lcmp = t_ind(l, gt_cmp)
#
# So, why do we choose gt₀ when the field subspace is orthogonal to the shape subspace?
#
# gt₀ represents the grid type of the corners of Yee's voxel whose edges are composed of the
# field lines.  The w-component of the field is located a half grid point away from these
# voxel corners along the w-direction.  This picture describes how to determine gt_cmp when
# the field subspace is equal to the shape subspace.
#
# However, when the field subspace is orthogonal to the shape subspace, the half grid point
# shift to apply to the w-component occurs in the direction normal to the shape subspace.
# Along that direction, the grid lines are not defined, because l is defined in the shape
# subspace.  Therefore, the half grid point shift does not have any effect in determining
# the locations of the w-component, and we can use gt₀ as the grid type in determining the
# field locations.
#
# For example, In the TE equation, the field subspace for the H-field is the z-direction and
# the shape subspace is the xy-plane.  When the boundary field types are the E-field type
# for both the x- and y-directions, gt₀ for the H-field (Hz) is [DUAL,DUAL].  We can use
# [DUAL,DUAL] in determining the xy-locations of the H-field, as the half grid point shift
# along the z-direction for the H-field does not affect the locations of the H-field.
function add_src!(jKd::AbsArrNumber{K₊₁},  # (K+1)-dimensional array of Je (electric) or Jm (magnetic); first K dimensions specify location; last 1 dimension specify field component
                  gt₀::SVec{K,GridType},  # grid type of voxel corners; generated by ft2gt.(ft, boundft)
                  bounds::Tuple2{SFloat{K}},  # bounds[NEG][k] = boundary of domain at negative end in k-direction
                  l::Tuple2{NTuple{K,VecFloat}},  # l[PRIM][k] = primal vertex locations in k-direction
                  ∆l::Tuple2{NTuple{K,VecFloat}},  # ∆l[PRIM][k] = (∆l at primal vertices in w) == diff(l[DUAL][k] including ghost point)
                  isbloch::SBool{K},  # Bloch boundary conditions
                  src::Source{K,Kf}...  # sources
                  ) where {K,K₊₁,Kf}
    for src = srcs
        add_src!(j3d, ft, boundft, bounds, l, ∆l, isbloch, src)
    end
end

# I could have defined distweights with more abstractly typed arguments, such as c::Real,
# bounds::AbsVecReal, l::AbsVecReal, ∆l::AbsVecReal, and have let the type inference system
# create distweights for each distinct set of types.  However, then I cannot control the
# behavior of the function perfectly, because the two numbers used in division could be
# integer or float.  Therefore, I implement the kernel for a specific type later, and here I
# define a wrapper that converts the arguments into those specific type.
#
# Another advantage of this approach is that even though several different versions of the
# wrapper will be generated, the kernel will be generated only once for the specific type.
# This makes the amount of compiled code minimal.
#
# Note that the situation here is a bit different from create_curl.  There, I wanted to
# create a curl operator of differnt types, such as Float, ComplexF, and even Int.  Therefore,
# I had to allow generation of different kernel code for different argument types.
distweights(c::Real, gt::GridType, bounds::AbsVecReal, l::AbsVecReal, ∆l::AbsVecReal, isbloch::Bool) =
    distweights(float(c), gt, SFloat{2}(bounds), float(l), float(∆l), isbloch)

# In a given Cartesian direction, if the source is located between two grid points,
# calculate the weights to distribute over the two points.  The output is in the form of
# ([ind₁, ind₂], [wt₁, wt₂]), where wtₖ is the weight assigned to the grid point at indₖ.
# Note that the physical dimension of the weight is 1/(length).
#
# If the source is placed exactly at one grid point, weight distribution is unnecessary and
# the output satisfies ind₂ = ind₁ and wt₂ = 0.  Note that the output type remains the same
# as before, which is a useful feature for achieving type stability of the function calling
# this function.
#
# This function works for both longitudinal and transversal discretization with respect to
# polarization.  Luckily, the consideration of the boundary conditions discussed at the
# beginning of this file works for both cases.  For example, suppose we are distributing
# along the x-grid an electric current source.  For Bloch, it is obvious that both
# polarizations are treated equally.  For the symmetry boundary condition, we can handle
# both polarization with the same code by determining whether the boundary condition is
# "zeroing" the source on the boundary or not.  Specifically, the zeroing boundary condition
# is defined as follows (see the definition of `zeroing_bc` in the function definition as
# well):
# - if the boundary condition is zeroing the source, the source put beyond the boundary is
# antisymmetric to the source put within the boundary, so the two sources are superposed
# destructively on the boundary (note the two sources are in adjacent cells), and
# - if the boundary condition is not zeroing the source, the source put beyond the boundary
# is symmetric to the source put within the boundary, so the two sources are superposed
# constructively on the boundary (note the two sources are in the same cell).
#
# Therefore,
# - For Jy and Jz (parallel to the x-boundary), the zeroing boundary condition is the
# symmetric E-field boundary condition.  Because the x-location of Jy and Jz is the E-field
# grid location, we know that a symmetric boundary is the E-field boundary if the x-location
# of Jy and Jz are primal.
# - For Jx (normal to the x-boundary), the zeroing boundary condition is the symmetric
# H-field boundary condition.  Because the x-location of Jx is the H-field grid location, we
# know again that a symmetric boundary is the H-field boundary if the x-location of Jx is
# primal.
#
# Note that regardless of the polarization of the source, the boundary condition is zeroing
# if the location type is primal.
#
# In the implementation, assume that the source "value" to distribute is something that does
# not change with discretization (e.g., Id of a point source and K of a plane source).  The
# returned weights are the factors to multiply to this source value.
#
# When distributing a source between two points on the opposite sides of the Bloch boundary,
# the Bloch phase factor is not multiplied. The source we are trying to distribute by this
# function is a delta-function-like source in the axis of consideration (e.g., axis normal
# to a plane source).  It feels a bit strange to impose phase difference between what were
# originally a single point.  Also, I think it would be much easier to multiply the Bloch
# phase factors to the entire 3D array of source when all the sources are assigned and
# properly distributed with weight factors.  That way, we don't have to implement
# multiplication of the Bloch phase factors in every source type.  However, this decision
# can change later if needed.
function distweights(c::Float,  # location of point (c means center)
                     gt::GridType,  # type of the grid planes
                     bounds::SFloat{2},  # [lneg, lpos]: negative- and positive-end boundaries of domain
                     l::AbsVecFloat,  # location of grid planes without ghost points
                     ∆l::AbsVecFloat,  # size of grid cell centered at l (not lg)
                     isbloch::Bool)  # boundary condition
    bounds[1]≤c≤bounds[2] || throw(ArgumentError("c = $c must be within bounds = $bounds."))

    L = bounds[2] - bounds[1]
    N = length(l)

    # Below, what if N = 1?  I think N = 1 will work fine for Bloch, for which indn = indp
    # = 1 and therefore the contribution from indn and indp will be added in source array
    # construction, and for the symmetry boundary that is incongruent with the field type,
    # for which indn = indp = 1 again but only indn will be considered and the current
    # density to set will remain constant regardless of the source position... Wait, we have
    # contribution from both boundaries in this case, so maybe we should consider one real
    # point source and two imaginary point sources...
    #
    # Anyway, I think the main reason to use N = 1 is for 2D simulation, and in that case
    # we use Bloch in that direction.  So, I think it should be OK to allow N = 1 only for
    # Bloch.
    N > 1 || isbloch || throw(ArgumentError("length(l) = $N must be > 1 for symmetry boundary (= non-Bloch)."))

    # Test if c is in the range where the result is not affected by the boundary.
    # - For l[1:end] being the primal grid locations except for the ghost point,
    #       - for Bloch, l[1] ≤ c < l[end], and
    #       - for symmetry, l[2] ≤ c < l[end].  (Note the minimum boundary beng l[2], not l[1].)
    # - For l[1:end] being the dual grid locations except for the ghost point,
    #       - for Bloch, l[1] ≤ c < l[end], and
    #       - for symmetry, l[1] ≤ c < l[end].
    # Therefore, only the case with the primal field type and symmetry boundary condition is
    # different from other cases.  I will call this a case with a boundary condition
    # "congruent" with the field type, because it is the boundary that is effective for the
    # field type of interest.
    zeroing_bc = gt==PRIM && !isbloch
    indn = zeroing_bc ? 2 : 1
    indp = N

    # Below, note that c < l[indp] is used rather than c ≤ l[indp].  If c = l[indp], then
    # ind₁ = findlast(x->x≤0, l.-c) is indp, and ind₂ = ind₁ + 1 is not well-defined when
    # indp = N.  The use of c < l[indp] avoids such a case.
    bc_!eff = l[indn] ≤ c < l[indp]  # boundary condition does not affect result
    if bc_!eff
        ind₁ = findlast(x->x≤0, l.-c)  # because l is sorted and c ≥ l[1], l[1]-c ≤ 0 and thus ind₁ ≠ nothing
        ind₂ = ind₁ + 1
        ∆c1 = c - l[ind₁]
        ∆lc = l[ind₂] - l[ind₁]
    else  # boundary affects result
        if c < l[indn]  # cannot happen for PRIM && Bloch
            ind₁ = indn
            ind₂ = N
            bnd = bounds[1]
        else # c ≥ l[indp] (not >)
            ind₁ = indp
            ind₂ = 1
            bnd = bounds[2]
        end

        ∆c1 = abs(l[ind₁] - c)
        if gt == PRIM
            ∆lc = abs(l[ind₁] - bnd)
        else  # gt == DUAL
            ∆lc = isbloch ? L-abs(l[ind₁]-l[ind₂]) : 2abs(l[ind₁]-bnd)
        end
    end

    r = ∆c1 / ∆lc

    # Below, scale wt₁ and wt₂ with 1-r and r, respectively,
    # - if c is in the region where the result is not affected by the boundary (bc_!eff), or
    # - even if c is in the region where the result is affected by the boundary,
    #       - if the boundary is Bloch, or
    #       - if the boundary is symmetry it is congruent with field type (zeroing_bc = true).
    wt₁ = 1.0 / ∆l[ind₁]
    (bc_!eff || isbloch || zeroing_bc) && (wt₁ *= 1.0-r)

    # Return only one weight factor
    # - if c is exactly at a grid point, or
    # - even if c is not exactly at a grid point, if it is in the region where the result is
    # affected by the boundary but the boundary is non-Bloch (= symmetry).
    if c==l[ind₁] || (!bc_!eff && !isbloch)
        # We want to return only one weight factor in this case, but in order to keep the
        # output SVec type the same (SVec{2} rather than SVec{1}), introduce a
        # dummy second weight factor.
        ind₂ = ind₁
        wt₂ = 0.0
    else
        # Return two weight factors.  Make sure all the cases handled by this else block
        # have ind₂ defined already.
        wt₂ = 1.0 / ∆l[ind₂]
        (bc_!eff || isbloch) && (wt₂ *= r)
    end

    ind = SVec(ind₁, ind₂)
    wt = SVec(wt₁, wt₂)

    return (ind, wt)
end

# funciton calc_blochphase()
# function j3d2j()

# When defining concrete source types below, try to use the following order of arguments in
# the constructor:
# - geometry (e.g., direction normal to a plane)
# - location (e.g., location for a point, intercept for a plane)
# - polarization
# - strength (e.g., current dipole, current, sheet current density)

# I think I can define a "uniform source" that populates a uniform current inside a given
# geometrical object specified with a GeometryPrimitives type.  This can cover various
# primitive source types, such as a point source, line source, and a plane source.  The
# source assignment function will accept a vector of sources, assign these sources similarly
# as assign_param! in assignment.jl (i.e., iterate over a vector of shapes and assign points
# included in each shape).  A main difference is that I have to find voxels that overlap with
# the shape, instead of points in the shape.  Once such voxels are found, I need to assign
# source current density to the edges of the voxels.  (I have to think about subvoxel
# smoothing in a voxel that is partially included in the shape.)
#
# Even more generally, I think I can create a single source that accepts a geometry and
# source strength function.  The uniform sources are a special case with a constant
# source strength.

include("pointsrc.jl")
include("planesrc.jl")


text:
include("berryphase.jl")

using BerryPhase

const S2 = (length(ARGS) >= 1 ? parse(Int, ARGS[1]) : 1)
const M = (length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 10)

bp = berryphase(dimer(S2,1.0), M)

println("γ/π = ", bp/pi)

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
"""
    mdft(ary, samples[, shift, Q])

Compute and return the 2D DFT via a matrix triple product.
Equivalent to fftshift(fft(ifftshift(ary))) if shift==0, samples=size(ary), Q=1.

samples is the number of samples in the output domain.  If it is an integer, it
is broadcast to both dimensions.  Otherwise it may be a tuple to provide each.

shift is used to move the origin.  If the output grid with shift=0 spanned
-10:1:9 Hz, shift=10 will move the output grid to 0:1:19Hz.  The unit is samples,
not hertz.  The first argument is along x, or the columns of the matrix, not dim 0.

Q controls the oversampling factor.  It is equivalent to zero padding the input
array by the same factor.  E.g. mdft of a 256x256 array with Q=2 is the same
as zero padding it to 512x512 and FFTing that.

See also: [`imdft`](@ref)

# References
"Fast computation of Lyot-style coronagraph propagation" Soummer et al, oe-15-24-15935
"""
function mdft(ary, samples::Integer; shift::Tuple{Real,Real} = (0, 0), Q = 1)
    return mdft(ary, (samples, samples), shift = shift, Q = Q)
end

function mdft(ary, samples::Tuple{Integer,Integer}; shift::Tuple{Real,Real} = (0, 0), Q = 1)
    ξ = fftrange(samples[2])'
    η = fftrange(samples[1])
    if shift[1] != 0
        ξ = ξ .+ shift[1]
    end
    if shift[2] != 0
        η = η .+ shift[2]
    end
    return _mdft(ary, ξ, η, Q)
end

function _mdft(ary, ξ, η, Q)
    n, m = size(ary)
    # X,Y,ξ,η look like -128 : 127, say
    X = fftrange(m)
    Y = fftrange(n)'
    kernel = -1im * 2 * π / Q
    pre = exp.(kernel / n .* (η * Y))
    post = exp.(kernel / m .* (X * ξ))
    return pre * ary * post
end

"""
    imdft(ary, samples[, shift, Q])

Compute and return the 2D inverse DFT via a matrix triple product.
Equivalent to fftshift(ifft(ifftshift(ary))) if shift==0, samples=size(ary), Q=1.

samples is the number of samples in the output domain.  If it is an integer, it
is broadcast to both dimensions.  Otherwise it may be a tuple to provide each.

shift is used to move the origin.  If the output grid with shift=0 spanned
-10:1:9 Hz, shift=10 will move the output grid to 0:1:19Hz.  The unit is samples,
not hertz.  The first argument is along x, or the columns of the matrix, not dim 0.

Q controls the oversampling factor.  It is equivalent to zero padding the input
array by the same factor.  E.g. mdft of a 256x256 array with Q=2 is the same
as zero padding it to 512x512 and FFTing that.

See also: [`mdft`](@ref)

# References
"Fast computation of Lyot-style coronagraph propagation" Soummer et al, oe-15-24-15935
"""
function imdft(ary, samples::Integer; shift::Tuple{Real,Real} = (0, 0), Q = 1)
    return imdft(ary, (samples, samples), shift = shift, Q = Q)
end

function imdft(
    ary,
    samples::Tuple{Integer,Integer};
    shift::Tuple{Real,Real} = (0, 0),
    Q = 1,
)
    ξ = fftrange(samples[2])'
    η = fftrange(samples[1])
    if shift[1] != 0
        ξ = ξ .+ shift[1]
    end
    if shift[2] != 0
        η = η .+ shift[2]
    end
    return _imdft(ary, ξ, η, Q)
end

function _imdft(ary, ξ, η, Q)
    n, m = size(ary)
    # X,Y,ξ,η look like -128 : 127, say
    X = fftrange(m)
    Y = fftrange(n)'
    kernel = 1im * 2 * π / Q
    pre = exp.(kernel / n .* (η * Y))
    post = exp.(kernel / m .* (X * ξ))
    return pre * ary * post
end

"""
    fftrange(n)

Computes a range that always matches FFT expectations.  I.e., for even N ranges
matches the example n=4 => -2:1:1 and for odd N matches the example 5 => -2:1:2.

This can be understood as "the output range shall always contain zero, if there
may be only one value for n÷2, it will be the negative one."

The range is always centered; frequency or sampling bins align to it with
fftshift.  It can be fftshifted to be placed in post-FFT coordinates.
"""
function fftrange(n)
    if iseven(n)
        return -(n ÷ 2):(n÷2-1)
    end
    return -(n ÷ 2):(n÷2)
end

max_stars_repo_name: JuliaOptics/Fourier

max_issues_repo_name: JuliaOptics/Fourier

max_forks_repo_name: JuliaOptics/Fourier


text:
using ITensors, Random, SweepContractor, ITensorNetworkAD
using ITensorNetworkAD.ITensorNetworks: ITensor_networks, line_network, TreeTensor
using ITensorNetworkAD.ITensorAutoHOOT: SubNetwork, batch_tensor_contraction

@testset "test the interface" begin
  LTN = LabelledTensorNetwork{Char}()
  LTN['A'] = Tensor(['D', 'B'], [i^2 - 2j for i in 0:2, j in 0:2], 0, 1)
  LTN['B'] = Tensor(['A', 'D', 'C'], [-3^i * j + k for i in 0:2, j in 0:2, k in 0:2], 0, 0)
  LTN['C'] = Tensor(['B', 'D'], [j for i in 0:2, j in 0:2], 1, 0)
  LTN['D'] = Tensor(['A', 'B', 'C'], [i * j * k for i in 0:2, j in 0:2, k in 0:2], 1, 1)

  sweep = sweep_contract(LTN, 100, 100; fast=true)
  out = ldexp(sweep...)
  @test isapprox(out, contract(ITensor_networks(LTN))[])
end

@testset "test on 2D grid" begin
  Random.seed!(1234)
  ITensors.set_warn_order(100)
  row, column, d = 8, 8, 2
  LTN = LabelledTensorNetwork{Tuple{Int,Int}}()
  for i in 1:row, j in 1:column
    adj = Tuple{Int,Int}[]
    i > 1 && push!(adj, (i - 1, j))
    j > 1 && push!(adj, (i, j - 1))
    i < row && push!(adj, (i + 1, j))
    j < column && push!(adj, (i, j + 1))
    LTN[i, j] = Tensor(adj, randn(d * ones(Int, length(adj))...), i, j)
  end
  tnet = ITensor_networks(LTN)
  element_grouping = line_network(tnet)
  tnet_mat = reshape(tnet, row, column)
  line_grouping = SubNetwork(tnet_mat[:, 1])
  for i in 2:column
    line_grouping = SubNetwork(line_grouping, tnet_mat[:, i]...)
  end

  function get_contracted_peps(rank)
    x = tnet_mat[:, 1]
    for i in 2:(column - 1)
      A = tnet_mat[:, i]
      x = contract(MPO(A), MPS(x); cutoff=1e-15, maxdim=rank)[:]
    end
    out_mps = contract(x..., tnet_mat[:, column]...)
    sweep = sweep_contract(LTN, rank, rank; fast=true)
    out = ldexp(sweep...)
    out2 = batch_tensor_contraction(
      TreeTensor, [element_grouping]; cutoff=1e-15, maxdim=rank
    )
    out3 = batch_tensor_contraction(
      TreeTensor, [line_grouping]; cutoff=1e-15, maxdim=rank, optimize=false
    )
    return out, ITensor(out2[1])[], ITensor(out3[1])[], out_mps[]
  end

  out_true, out_element, out_line, out_mps = get_contracted_peps(d^(Int(row / 2)))
  @test abs((out_true - out_element) / out_true) < 1e-3
  @test abs((out_true - out_line) / out_true) < 1e-3
  @test abs((out_true - out_mps) / out_true) < 1e-3
  for rank in [1, 2, 3, 4, 6, 8, 10, 12, 14, 15, 16]
    out, out_element, out_line, out_mps = get_contracted_peps(rank)
    error_sweepcontractor = abs((out - out_true) / out_true)
    error_element = abs((out_element - out_true) / out_true)
    error_line = abs((out_line - out_true) / out_true)
    error_mps = abs((out_mps - out_true) / out_true)
    print(
      "maxdim, ",
      rank,
      ", error_sweepcontractor, ",
      error_sweepcontractor,
      ", error_element, ",
      error_element,
      ", error_line, ",
      error_line,
      ", error_mps, ",
      error_mps,
      "\n",
    )
  end
end


text:
#ДАНО: Робот - в произвольной клетке поля (без внутренних перегородок и маркеров).
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы.
function mark_frame_perimetr!(r)
    x,y = num_sudwest!(r)
    markframe!(r)
    return_sudwest!(r,x,y)
end
function num_sudwest!(r)
    x=0
    y=0
    while isborder(r,Sud)==0
        move!(r,Sud)
        x+=1
    end
    while isborder(r,West)==0
        move!(r,West)
        y+=1
    end
    return x,y
end
function markframe!(r)
    for side in [Nord,Ost,Sud,West]
        while isborder(r,side)==0
            putmarker!(r)
            move!(r,side)
        end
    end
end
function return_sudwest!(r,x,y)
    for i in 1:x
        move!(r,Nord)
    end
    for i in 1:y
        move!(r,Ost)
    end
end
mark_frame_perimetr!(r)


text:
module Peacock

using PyPlot
using LinearAlgebra, FFTW
using CUDA

include("utils.jl")

include("geometry.jl")
export Geometry

include("plane_wave_basis.jl")
export BrillouinZoneCoordinate

include("modes.jl")
export Mode, HilbertSpace, get_field

include("solver.jl")
export Solver, TE, TM, solve

include("plotting.jl")
export plot

include("band_diagrams.jl")
export plot_band_diagram

include("wilson_loops.jl")
export plot_wilson_loop_winding

# Submodule that defines my commonly used crystals
include("Zoo.jl")

end # module


text:
using EmpiricalMeasures
using ArraysOfArrays
using StatsBase
using LinearAlgebra
using Random
using Distributions
using Test

@testset "EmpiricalMeasures.jl" begin

    @testset "Declaring Empirical Measures" begin

        @testset "MvDiscreteNonParametric from ArraysOfArrays" begin

            Random.seed!(7)
            n = 4
            m = 2
            A = nestedview(rand(n, m)')
            p = normalize!(rand(n), 1)

            # Passing probabilities
            μ = @inferred(MvDiscreteNonParametric{Float64,Float64,typeof(A),typeof(p)}(A, p))
            @test support(μ) == A
            @test length(μ) == m
            @test size(μ) == size(flatview(A)')
            @test probs(μ)   == p

            μ = @inferred(MvDiscreteNonParametric(A, p))
            @test support(μ) == A
            @test length(μ) == m
            @test size(μ) == size(flatview(A)')

            # Without passing probabilities
            μ = @inferred(MvDiscreteNonParametric(A))
            @test support(μ) == A
            @test length(μ) == m
            @test size(μ) == size(flatview(A)')
            @test probs(μ)   == fill(1 / n, n)

            # Array of arrays without ArraysOfArrays.jl
            n, m = 3, 2
            p = ([3 / 5, 1 / 5, 1 / 5])
            A = [[1,0],[1,1],[0,1]]
            μ = @inferred(MvDiscreteNonParametric(A, p))

            @test support(μ) == A
            @test length(μ) == m
            @test size(μ) == (length(A), length(A[1]))
            @test probs(μ) == p

        end

        @testset "MvDiscreteNonParametric from Matrix" begin

            Random.seed!(7)
            n, m = 10, 5
            A = rand(n, m)
            p = normalize!(rand(n), 1)

            # Passing probabilities
            μ = @inferred(MvDiscreteNonParametric(A, p))

            @test flatview(support(μ))' == A
            @test length(μ) == m
            @test size(μ) == size(A)
            @test probs(μ)   == p

            # Without passing probabilities
            μ = @inferred(MvDiscreteNonParametric(A))

            @test flatview(support(μ))' == A
            @test length(μ) == m
            @test size(μ) == size(A)
            @test probs(μ) == fill(1 / n, n)

        end

        @testset "empiricalmeasure" begin

            Random.seed!(7)
            n = 4
            m = 2
            A = nestedview(rand(n, m)')
            p = normalize!(rand(n), 1)

            # Passing probabilities
            μ = @inferred(empiricalmeasure(A, p))
            @test support(μ) == A
            @test length(μ) == m
            @test size(μ) == size(flatview(A)')
            @test probs(μ)   == p

            A = rand(n, m)
            μ = @inferred(empiricalmeasure(A))
            @test flatview(support(μ))' == A
            @test length(μ) == m
            @test size(μ) == size(A)
            @test probs(μ) == fill(1 / n, n)

        end
        @testset "empiricalmeasure 1D dispatch" begin

            n = 4
            A = rand(n)
            p = normalize!(rand(n), 1)
            μ = @inferred(empiricalmeasure(A, p))

            @test typeof(μ) <: DiscreteNonParametric
            @test μ.support == sort(A)
            @test μ.p == p[sortperm(A)]

            μ = @inferred(empiricalmeasure(A))
            @test μ.support == sort(A)
            @test μ.p == fill(1 / n, n)
        end
    end

    @testset "Functionalities" begin

        function variance(d)
            v = zeros(length(d))
            for i in 1:length(d)
                s = flatview(d.support)'[:,i]
                mₛ = mean(d)[i]
                v[i] = sum(abs2.(s .- mₛ), Weights(d.p))
            end
            return v
        end

        function covariance(d)
            n = length(d)
            v = zeros(n, n)
            for i in 1:n, j in 1:n
                s = flatview(d.support)'[:,i]
                mₛ = mean(d)[i]

                u = flatview(d.support)'[:,j]
                mᵤ = mean(d)[j]

                v[i,j] = sum((s .- mₛ) .* (u .- mᵤ), Weights(d.p))
            end
            return v
        end

        Random.seed!(7)
        n, m = 7, 9
        A = rand(n, m)
        p = normalize!(rand(n), 1)
        μ = @inferred(MvDiscreteNonParametric(A, p))

        @test mean(μ) == mean(flatview(μ.support)[:,:], Weights(p), dims=2)[:]
        @test var(μ) ≈ variance(μ)
        @test cov(μ) ≈ covariance(μ)
        @test pdf(μ, flatview(μ.support)) ≈ μ.p
        @test pdf(μ, zeros(m)) == 0.0
        @test entropy(μ) == entropy(μ.p)
        @test entropy(μ, 2) == entropy(μ.p, 2)

    end

    @testset "Sampling" begin
        Random.seed!(7)
        μ = empiricalmeasure(rand(10, 10))
        @test rand(μ) in μ.support

        for sample in eachcol(rand(μ, 10))
           @test sample in μ.support
        end

        A = rand(3, 2)
        μ = empiricalmeasure(A, [0.2,0.5,0.3])

        for i in 1:3
        samples = nestedview(rand(μ, 10000))
            @test abs(mean([A[i,:] == s for s in samples]) - μ.p[i]) < 0.05

        end
    end

end


text:
import Base: setindex!, getindex

export AbstractDAQ, DAQParams, SinkImpedance, SINK_FIFTY_OHM, SINK_HIGH, DAQTxChannelSettings, DAQChannelParams, DAQFeedback, DAQTxChannelParams, DAQRxChannelParams,
       createDAQChannels, createDAQParams, startTx, stopTx, setTxParams, currentFrame, readData,
       numRxChannelsTotal, numTxChannelsTotal, numRxChannelsActive, numTxChannelsActive,
       DAQ, readDataPeriod, currentPeriod, getDAQ, getDAQs,
       channelIdx, limitPeak, sinkImpedance, allowedWaveforms, isWaveformAllowed,
       feedbackChannelID, feedbackCalibration, calibration

abstract type AbstractDAQ <: Device end
abstract type DAQParams <: DeviceParams end

@enum SinkImpedance begin
  SINK_FIFTY_OHM
  SINK_HIGH
end

struct DAQTxChannelSettings
  "Applied channel voltage. Dimensions are (components, channels, periods)."
  amplitudes::Array{typeof(1.0u"V"), 3}
  "Applied channel phase. Dimensions are (components, channels, periods)."
  phases::Array{typeof(1.0u"rad"), 3}
  "Minimum time for changing phase and amplitude to the given settings"
  changeTime::typeof(1.0u"s")
  "Channel mapping from ID to index."
  mapping::Dict{String, Integer}

  function DAQTxChannelSettings(amplitudes, phases, changeTime, mapping)
    if size(amplitudes) != size(phases) || size(amplitudes, 2) != length(mapping)
      error("The sizes of `phases` and `amplitudes` as well as the number of channels have to match.")
    end

    if changeTime < 0.0u"s"
      error("The change time can only be positive.")
    end

    return new(amplitudes, phases, changeTime, mapping)
  end
end

numChannels(settings::DAQTxChannelSettings) = size(settings.amplitudes, 2)
channelIDs(settings::DAQTxChannelSettings) = keys(settings.mapping)
amplitudes(settings::DAQTxChannelSettings) = settings.amplitudes
amplitudes(settings::DAQTxChannelSettings, channelID::String) = settings.amplitudes[:, mapping[channelID], :]
phases(settings::DAQTxChannelSettings) = settings.phases
phases(settings::DAQTxChannelSettings, channelID::String) = settings.phases[:, mapping[channelID], :]
changeTime(settings::DAQTxChannelSettings) = settings.changeTime

abstract type DAQChannelParams end
abstract type TxChannelParams <: DAQChannelParams end
abstract type RxChannelParams <: DAQChannelParams end

Base.@kwdef struct DAQFeedback
  channelID::AbstractString
  calibration::Union{typeof(1.0u"T/V"), Nothing} = nothing
end

Base.@kwdef struct DAQTxChannelParams <: TxChannelParams
  channelIdx::Int64
  limitPeak::typeof(1.0u"V")
  sinkImpedance::SinkImpedance = SINK_HIGH
  allowedWaveforms::Vector{Waveform} = [WAVEFORM_SINE]
  feedback::Union{DAQFeedback, Nothing} = nothing
  calibration::Union{typeof(1.0u"V/T"), Nothing} = nothing
end

Base.@kwdef struct DAQRxChannelParams <: RxChannelParams
  channelIdx::Int64
end

"Create DAQ channel description from device dict part."
function createDAQChannels(dict::Dict{String, Any})
  channels = Dict{String, DAQChannelParams}()
  for (key, value) in dict
    splattingDict = Dict{Symbol, Any}()
    if value["type"] == "tx"
      splattingDict[:channelIdx] = value["channel"]
      splattingDict[:limitPeak] = uparse(value["limitPeak"])

      if haskey(value, "sinkImpedance")
        splattingDict[:sinkImpedance] = value["sinkImpedance"] == "FIFTY_OHM" ? SINK_FIFTY_OHM : SINK_HIGH
      end

      if haskey(value, "allowedWaveforms")
        splattingDict[:allowedWaveforms] = toWaveform.(value["allowedWaveforms"])
      end

      if haskey(value, "feedback")
        channelID=value["feedback"]["channelID"]
        calibration=uparse(value["feedback"]["calibration"])

        splattingDict[:feedback] = DAQFeedback(channelID=channelID, calibration=calibration)
      end

      if haskey(value, "calibration")
        splattingDict[:calibration] = uparse.(value["calibration"])
      end

      channels[key] = DAQTxChannelParams(;splattingDict...)
    elseif value["type"] == "rx"
      channels[key] = DAQRxChannelParams(channelIdx=value["channel"])
    end
  end

  return channels
end

# Generic case
function createDAQChannels(::Type{T}, dict::Dict{String, Any}) where {T <: DAQParams}
  return createDAQChannels(dict)
end

"Create the params struct from a dict. Typically called during scanner instantiation."
function createDAQParams(DAQType::DataType, dict::Dict{String, Any})
  @assert DAQType <: DAQParams "The supplied type `$type` cannot be used for creating DAQ params, since it does not inherit from `DAQParams`."

  # Extract all main section fields which means excluding `channels`
  mainSectionFields = [string(field) for field in fieldnames(DAQType) if field != :channels]

  # Split between main section fields and channels, which are dictionaries
  channelDict = Dict{String, Any}()
  for (key, value) in dict
    if value isa Dict && !(key in mainSectionFields)
      channelDict[key] = value

      # Remove key in order to process the rest with the standard function
      delete!(dict, key)
    end
  end

  splattingDict = dict_to_splatting(dict)
  splattingDict[:channels] = createDAQChannels(DAQType, channelDict)

  try
    return DAQType(;splattingDict...)
  catch e
    if e isa UndefKeywordError
      throw(ScannerConfigurationError("The required field `$(e.var)` is missing in your configuration "*
                                      "for a device with the params type `$DAQType`."))
    elseif e isa MethodError
      @warn e.args e.world e.f
      throw(ScannerConfigurationError("A required field is missing in your configuration for a device "*
                                      "with the params type `$DAQType`. Please check "*
                                      "the causing stacktrace."))
    else
      rethrow()
    end
  end
end

#include("Control.jl")
#include("Plotting.jl")
#include("Parameters.jl")

@mustimplement startTx(daq::AbstractDAQ)
@mustimplement stopTx(daq::AbstractDAQ)
@mustimplement setTxParams(daq::AbstractDAQ, sequence::Sequence)
@mustimplement currentFrame(daq::AbstractDAQ)
@mustimplement readData(daq::AbstractDAQ, startFrame, numFrames)
@mustimplement readDataPeriods(daq::AbstractDAQ, startPeriod, numPeriods)
@mustimplement refToField(daq::AbstractDAQ, d::Int64)

@mustimplement setSequenceParams(daq::AbstractDAQ, sequence::Sequence) # Needs to be able to update seqeuence parameters
@mustimplement prepareSequence(daq::AbstractDAQ, sequence::Sequence) # Sequence can be prepared before started
@mustimplement endSequence(daq::AbstractDAQ) # Sequence can be ended outside of producer
@mustimplement prepareTx(daq::AbstractDAQ, sequence::Sequence; allowControlLoop = true) # Tx can be set outside of producer
# Producer prepares a proper sequence if allowed too, then starts it and writes the resulting chunks to the channel

@mustimplement numTxChannelsTotal(daq::AbstractDAQ)
@mustimplement numRxChannelsTotal(daq::AbstractDAQ)
@mustimplement numTxChannelsActive(daq::AbstractDAQ)
@mustimplement numRxChannelsActive(daq::AbstractDAQ)
@mustimplement numRxChannelsReference(daq::AbstractDAQ)
@mustimplement numRxChannelsMeasurement(daq::AbstractDAQ)
@mustimplement numComponentsMax(daq::AbstractDAQ)
@mustimplement canPostpone(daq::AbstractDAQ)
@mustimplement canConvolute(daq::AbstractDAQ)

getDAQs(scanner::MPIScanner) = getDevices(scanner, AbstractDAQ)
getDAQ(scanner::MPIScanner) = getDevice(scanner, AbstractDAQ)

# function dataConversionFactor(daq::AbstractDAQ) #default
#   factor = zeros(2, numRxChannels(daq))
#   factor[1,:] = 1.0
#   factor[2,:] = 0.0
#   return factor
# end

channel(daq::AbstractDAQ, channelID::AbstractString) = daq.params.channels[channelID]
channelIdx(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).channelIdx
channelIdx(daq::AbstractDAQ, channelIDs::Vector{<:AbstractString}) = [channel(daq, channelID).channelIdx for channelID in channelIDs]
limitPeak(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).limitPeak
sinkImpedance(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).sinkImpedance
allowedWaveforms(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).allowedWaveforms
isWaveformAllowed(daq::AbstractDAQ, channelID::AbstractString, waveform::Waveform) = waveform in allowedWaveforms(daq, channelID)
feedback(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).feedback
feedbackChannelID(daq::AbstractDAQ, channelID::AbstractString) = feedback(daq, channelID).channelID
feedbackCalibration(daq::AbstractDAQ, channelID::AbstractString) = feedback(daq, channelID).calibration
calibration(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).calibration

#### Measurement Related Functions ####
@mustimplement startProducer(channel::Channel, daq::AbstractDAQ, numFrames)
@mustimplement channelType(daq::AbstractDAQ) # What is written to the channel
@mustimplement AsyncBuffer(daq::AbstractDAQ) # Buffer structure that contains channel elements
@mustimplement updateAsyncBuffer!(buffer::AsyncBuffer, chunk) # Adds channel element to buffer
@mustimplement retrieveMeasAndRef!(buffer::AsyncBuffer, daq::AbstractDAQ) # Retrieve all available measurement and reference frames from the buffer

function asyncProducer(channel::Channel, daq::AbstractDAQ, sequence::Sequence; prepTx = true, prepSeq = true, endSeq = true)
  if prepTx
      prepareTx(daq, sequence)
  end
  if prepSeq
      setSequenceParams(daq, sequence)
      prepareSequence(daq, sequence)
  end

  numFrames = acqNumFrames(sequence) * acqNumFrameAverages(sequence)
  endFrame = startProducer(channel, daq, numFrames)

  if endSeq
      endSequence(daq, endFrame)
  end
end

function addFramesToAvg(avgBuffer::FrameAverageBuffer, frames::Array{Float32, 4})
  #setIndex - 1 = how many frames were written to the buffer

  # Compute how many frames there will be
  avgSize = size(avgBuffer.buffer)
  resultFrames = div(avgBuffer.setIndex - 1 + size(frames, 4), avgSize[4])

  result = nothing
  if resultFrames > 0
    result = zeros(Float32, avgSize[1], avgSize[2], avgSize[3], resultFrames)
  end

  setResult = 1
  fr = 1
  while fr <= size(frames, 4)
    # How many left vs How many can fit into avgBuffer
    fit = min(size(frames, 4) - fr, avgSize[4] - avgBuffer.setIndex)

    # Insert into buffer
    toFrames = fr + fit
    toAvg = avgBuffer.setIndex + fit
    avgBuffer.buffer[:, :, :, avgBuffer.setIndex:toAvg] = frames[:, :, :, fr:toFrames]
    avgBuffer.setIndex += length(avgBuffer.setIndex:toAvg)
    fr = toFrames + 1

    # Average and add to result
    if avgBuffer.setIndex - 1 == avgSize[4]
      avgFrame = mean(avgBuffer.buffer, dims=4)[:,:,:,:]
      result[:, :, :, setResult] = avgFrame
      setResult += 1
      avgBuffer.setIndex = 1
    end
  end

  return result
end

function updateFrameBuffer!(measState::SequenceMeasState, daq::AbstractDAQ)
  uMeas, uRef = retrieveMeasAndRef!(measState.asyncBuffer, daq)
  if !isnothing(uMeas)
    #isNewFrameAvailable, fr =
    handleNewFrame(measState.type, measState, uMeas)
    #if isNewFrameAvailable && fr > 0
    #  measState.currFrame = fr
    #  measState.consumed = false
    #end
  end
end

function handleNewFrame(::RegularAsyncMeas, measState::SequenceMeasState, uMeas)
  isNewFrameAvailable = false

  fr = addFramesFrom(measState, uMeas)
  isNewFrameAvailable = true

  return isNewFrameAvailable, fr
end

function handleNewFrame(::FrameAveragedAsyncMeas, measState::SequenceMeasState, uMeas)
  isNewFrameAvailable = false

  fr = 0
  framesAvg = addFramesToAvg(measState.avgBuffer, uMeas)
  if !isnothing(framesAvg)
    fr = addFramesFrom(measState, framesAvg)
    isNewFrameAvailable = true
  end

  return isNewFrameAvailable, fr
end

function addFramesFrom(measState::SequenceMeasState, frames::Array{Float32, 4})
  fr = measState.nextFrame
  to = fr + size(frames, 4) - 1
  limit = size(measState.buffer, 4)
  @info "Add frames $fr to $to to framebuffer with $limit size"
  if to <= limit
    measState.buffer[:,:,:,fr:to] = frames
    measState.nextFrame = to + 1
    return fr
  end
  return -1
end

include("RedPitayaDAQ.jl")
include("DummyDAQ.jl")
include("SimpleSimulatedDAQ.jl")