text:
# NucleicAcids
# ============
#
# DNA and RNA nucleotide types.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# NucleicAcid Encoding
# -------------------
#
# Unambiguous nucleotides are represented in one-hot encoding as follows:
#
#   | NucleicAcid | Bits |
#   | ----------- | ---- |
#   |     A       | 0001 |
#   |     C       | 0010 |
#   |     G       | 0100 |
#   |    T/U      | 1000 |
#
# Ambiguous nucleotides are bitwise OR of these four nucleotides. For example, R
# , A or G, is represented as 0101 (= A: 0001 | G: 0100). The gap symbol is
# always 0000.  The meaningful four bits are stored in the least significant
# bits of a byte.

abstract NucleicAcid
bitstype 8 DNA <: NucleicAcid
bitstype 8 RNA <: NucleicAcid

# Conversion from/to integers
# ---------------------------

Base.convert(::Type{DNA}, nt::UInt8) = reinterpret(DNA, nt)
Base.convert(::Type{RNA}, nt::UInt8) = reinterpret(RNA, nt)
Base.convert(::Type{UInt8}, nt::DNA) = reinterpret(UInt8, nt)
Base.convert(::Type{UInt8}, nt::RNA) = reinterpret(UInt8, nt)
Base.convert{T<:Number,S<:NucleicAcid}(::Type{T}, nt::S) = convert(T, UInt8(nt))
Base.convert{T<:Number,S<:NucleicAcid}(::Type{S}, nt::T) = convert(S, UInt8(nt))

# Bit Operations
# --------------

function Base.:~{N<:NucleicAcid}(x::N)
    return reinterpret(N, ~reinterpret(UInt8, x) & 0b1111)
end

function Base.:|{N<:NucleicAcid}(x::N, y::N)
    return reinterpret(N, reinterpret(UInt8, x) | reinterpret(UInt8, y))
end

function Base.:&{N<:NucleicAcid}(x::N, y::N)
    return reinterpret(N, reinterpret(UInt8, x) & reinterpret(UInt8, y))
end

function Base.:-{N<:NucleicAcid}(x::N, y::N)
    return Int(x) - Int(y)
end

function Base.:-{N<:NucleicAcid}(x::N, y::Integer)
    return x + (-y)
end

function Base.:+{N<:NucleicAcid}(x::N, y::Integer)
    return reinterpret(N, (UInt8(x) + y % UInt8) & 0b1111)
end

function Base.isless{N<:NucleicAcid}(x::N, y::N)
    return isless(reinterpret(UInt8, x), reinterpret(UInt8, y))
end

@inline function Base.count_ones(nt::NucleicAcid)
    return count_ones(reinterpret(UInt8, nt))
end

function Base.trailing_zeros(nt::NucleicAcid)
    return trailing_zeros(reinterpret(UInt8, nt))
end

function gap{N<:NucleicAcid}(::Type{N})
    return reinterpret(N, 0b0000)
end

"""
    isGC(nt::NucleicAcid)

Test if `nt` is surely either guanine or cytosine.
"""
function isGC(nt::NucleicAcid)
    bits = reinterpret(UInt8, nt)
    return bits != 0 && (bits & 0b1001) == 0
end

"""
    ispurine(nt::NucleicAcid)

Test if nucleotide is surely a purine.
"""
@inline function ispurine(nt::NucleicAcid)
    bits = reinterpret(UInt8, nt)
    return bits != 0 && (bits & 0b1010) == 0
end

"""
    ispyrimidine(nt::NucleicAcid)

Test if nucleotide is surely a pyrimidine.
"""
@inline function ispyrimidine(nt::NucleicAcid)
    bits = reinterpret(UInt8, nt)
    return bits != 0 && (bits & 0b0101) == 0
end

"""
    isambiguous(nt::NucleicAcid)

Test if `nt` is ambiguous nucleotide.
"""
@inline function isambiguous(nt::NucleicAcid)
    return count_ones(nt) > 1
end

"""
    iscertain(nt::NucleicAcid)

Test if `nt` is a non-ambiguous nucleotide e.g. ATCG.
"""
@inline function iscertain(nt::NucleicAcid)
    return count_ones(nt) == 1
end

"""
    isgap(nt::NucleicAcid)
Test if `nt` is a gap.
"""
@inline function isgap(nt::NucleicAcid)
    return count_ones(nt) == 0
end

"""
    complement(nt::NucleicAcid)

Return the complementary nucleotide of `nt`.
"""
function complement(nt::NucleicAcid)
    bits = compatbits(nt)
    return reinterpret(
        typeof(nt),
        (bits & 0x01) << 3 | (bits & 0x08) >> 3 |
        (bits & 0x02) << 1 | (bits & 0x04) >> 1)
end

function Base.isvalid{T<:NucleicAcid}(::Type{T}, x::Integer)
    return 0 ≤ x < 16
end

function Base.isvalid(nt::NucleicAcid)
    return reinterpret(UInt8, nt) ≤ 0b1111
end

"""
    iscompatible(x, y)

Return `true` if and only if `x` and `y` are compatible with each other (i.e.
`x` and `y` can be the same symbol).

`x` and `y` must be the same type (`DNA`, `RNA` or `AminoAcid`).

# Examples

```julia
julia> iscompatible(DNA_A, DNA_A)
true

julia> iscompatible(DNA_C, DNA_N)  # DNA_N can be DNA_C
true

julia> iscompatible(DNA_C, DNA_R)  # DNA_R (A or G) cannot be DNA_C
false

julia> iscompatible(AA_A, AA_X)    # AA_X can be AA_A
true

```
"""
@inline function iscompatible{T<:NucleicAcid}(x::T, y::T)
    return compatbits(x) & compatbits(y) != 0
end

# Return the compatibility bits of `nt`.
@inline function compatbits(nt::NucleicAcid)
    return reinterpret(UInt8, nt)
end

# NucleicAcid encoding definition
# ------------------------------

# DNA

# lookup table for characters
const char_to_dna = [0x80 for _ in 0x00:0xff]
const dna_to_char = Vector{Char}(16)

# derived from "The DDBJ/ENA/GenBank Feature Table Definition"
# §7.4.1 Nucleotide base code (IUPAC)
# http://www.insdc.org/documents/feature_table.html#7.4.1
for (char, doc, bits) in [
        ('-', "DNA Gap",                                   0b0000),
        ('A', "DNA Adenine",                               0b0001),
        ('C', "DNA Cytosine",                              0b0010),
        ('G', "DNA Guanine",                               0b0100),
        ('T', "DNA Thymine",                               0b1000),
        ('M', "DNA Adenine or Cytosine",                   0b0011),
        ('R', "DNA Adenine or Guanine",                    0b0101),
        ('W', "DNA Adenine or Thymine",                    0b1001),
        ('S', "DNA Cytosine or Guanine",                   0b0110),
        ('Y', "DNA Cytosine or Thymine",                   0b1010),
        ('K', "DNA Guanine or Thymine",                    0b1100),
        ('V', "DNA Adenine, Cytosine or Guanine",          0b0111),
        ('H', "DNA Adenine, Cytosine or Thymine",          0b1011),
        ('D', "DNA Adenine, Guanine or Thymine",           0b1101),
        ('B', "DNA Cytosine, Guanine or Thymine",          0b1110),
        ('N', "DNA Adenine, Cytosine, Guanine or Thymine", 0b1111)]
    var = Symbol("DNA_", char != '-' ? char : "Gap")
    @eval begin
        @doc $(doc) const $(var) = reinterpret(DNA, $(bits))
        char_to_dna[$(Int(char)+1)] = char_to_dna[$(Int(lowercase(char))+1)] = $(bits)
        dna_to_char[$(Int(bits)+1)] = $(char)
    end
end

@eval alphabet(::Type{DNA}) = $(tuple([reinterpret(DNA, x)
                                                 for x in 0b0000:0b1111]...))

const ACGT = (DNA_A, DNA_C, DNA_G, DNA_T)
const ACGTN = (DNA_A, DNA_C, DNA_G, DNA_T, DNA_N)

# RNA

# lookup table for characters
const char_to_rna = [0x80 for _ in 0x00:0xff]
const rna_to_char = Vector{Char}(16)

for (char, doc, dna) in [
        ('-', "RNA Gap",                                  DNA_Gap),
        ('A', "RNA Adenine",                              DNA_A  ),
        ('C', "RNA Cytosine",                             DNA_C  ),
        ('G', "RNA Guanine",                              DNA_G  ),
        ('U', "RNA Uracil",                               DNA_T  ),
        ('M', "RNA Adenine or Cytosine",                  DNA_M  ),
        ('R', "RNA Adenine or Guanine",                   DNA_R  ),
        ('W', "RNA Adenine or Uracil",                    DNA_W  ),
        ('S', "RNA Cytosine or Guanine",                  DNA_S  ),
        ('Y', "RNA Cytosine or Uracil",                   DNA_Y  ),
        ('K', "RNA Guanine or Uracil",                    DNA_K  ),
        ('V', "RNA Adenine, Cytosine or Guanine",         DNA_V  ),
        ('H', "RNA Adenine, Cytosine or Uracil",          DNA_H  ),
        ('D', "RNA Adenine, Guanine or Uracil",           DNA_D  ),
        ('B', "RNA Cytosine, Guanine or Uracil",          DNA_B  ),
        ('N', "RNA Adenine, Cytosine, Guanine or Uracil", DNA_N  )]
    var = Symbol("RNA_", char != '-' ? char : "Gap")
    @eval begin
        @doc $(doc) const $(var) = reinterpret(RNA, $(dna))
        char_to_rna[$(Int(char)+1)] = char_to_rna[$(Int(lowercase(char)+1))] = reinterpret(UInt8, $(dna))
        rna_to_char[$(Int(dna)+1)] = $(char)
    end
end

@eval alphabet(::Type{RNA}) = $(tuple([reinterpret(RNA, x)
                                                 for x in 0b0000:0b1111]...))

const ACGU = (RNA_A, RNA_C, RNA_G, RNA_U)
const ACGUN = (RNA_A, RNA_C, RNA_G, RNA_U, RNA_N)

# Print functions
# ---------------

function Base.convert(::Type{DNA}, c::Char)
    if c > '\xff'
        throw(InexactError())
    end
    @inbounds dna = char_to_dna[Int(c) + 1]
    if !isvalid(DNA, dna)
        throw(InexactError())
    end
    return reinterpret(DNA, dna)
end

function Base.convert(::Type{RNA}, c::Char)
    if c > '\xff'
        throw(InexactError())
    end
    @inbounds rna = char_to_rna[Int(c) + 1]
    if !isvalid(RNA, rna)
        throw(InexactError())
    end
    return reinterpret(RNA, rna)
end

function Base.convert(::Type{Char}, nt::DNA)
    return dna_to_char[convert(UInt8, nt) + 1]
end

function Base.convert(::Type{Char}, nt::RNA)
    return rna_to_char[convert(UInt8, nt) + 1]
end

function Base.show(io::IO, nt::DNA)
    if isvalid(nt)
        if nt == DNA_Gap
            write(io, "DNA_Gap")
        else
            write(io, "DNA_", Char(nt))
        end
    else
        write(io, "Invalid DNA")
    end
    return
end

function Base.show(io::IO, nt::RNA)
    if isvalid(nt)
        if nt == RNA_Gap
            write(io, "RNA_Gap")
        else
            write(io, "RNA_", Char(nt))
        end
    else
        write(io, "Invalid RNA")
    end
    return
end

function Base.print(io::IO, nt::NucleicAcid)
    if !isvalid(nt)
        throw(ArgumentError("nucleic acid is invalid"))
    end
    write(io, Char(nt))
    return
end


text:
using Test
using ExactDiagonalization.FED
using QuantumLattices
import QuantumLattices.Mathematics.AlgebraOverFields: idtype
using Arpack: eigs

@testset "GBasis" begin
    basis=GBasis{UInt}(2)
    @test basis|>kind==basis|>typeof|>kind==:g
    @test basis|>idtype==basis|>typeof|>idtype==Int
    @test dimension(basis)==2^2
    @test id(basis)==2
    @test basis[4]==3
    @test searchsortedfirst(basis,3)==4
    @test string(basis)=="GBasis(2):\n  0\n  1\n  10\n  11\n"
    @test repr(basis)=="GBasis(2)"
end

@testset "matrix" begin
    lattice=Lattice("L2P",[Point(PID(1),(0.0,)),Point(PID(2),(1.0,))])
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=1,nspin=2,nnambu=2),lattice.pids)
    table=Table(config,usualfockindextotuple)
    t=Hopping{'F'}(:t,1.0,1)
    U=Hubbard{'F'}(:U,1.0)
    μ=Onsite{'F'}(:μ,-0.5)
    g=Generator((t,U,μ),Bonds(lattice),config,table,false)
    m=matrix(expand(g),GBasis{UInt}(4))
    @test eigs(m,which=:SR,ritzvec=false)[1][1]≈-2.5615528128088303
end

ext: jl


text:
abstract type Data end
abstract type Query end
abstract type Queries end

const QueryIndex = Int
const Interval = Tuple{Int64, Int64}

struct MWState
    real::Data
    synthetic::Data
    queries::Queries
    real_answers::Array{Float64, 1}
    measurements::Dict{Int, Float64}
    scale::Float64
end

struct MWParameters
    epsilon::Float64
    iterations::Int64
    repetitions::Int64
    noisy_init::Bool
    verbose::Bool
    init_budget::Float64
    noisy_max_budget::Float64
end

struct Tabular <: Data
    data::Array{Float64, 2}
end

"""
    get(qs::Queries, qindex::QueryIndex)

Returns Query corresponding to given QueryIndex.
"""
function get(qs::Queries, qindex::QueryIndex)
    throw("`get` not implemented for `$(typeof(qs))`.")
end

"""
    evaluate(q::Query, d::Data)

Returns float value of Query evaluated on Data.
"""
function evaluate(q::Query, d::Data)
    throw("`evaluate` not implemented for `$(typeof(q))`, `$(typeof(d))`.")
end

"""
    evaluate(qs::Queries, d::Data)

Returns vector of floats corresponding to Queries evaluated on Data.
"""
function evaluate(qs::Queries, d::Data)
    throw("`evaluate` not implemented for `$(typeof(qs))`, `$(typeof(d))`.")
end

"""
    initialize(qs::Queries, d::Data, ps::MWParameters)

Returns MWState initialization for given queries, data, and parameters.
"""
function initialize(qs::Queries, d::Data, ps::MWParameters)
    throw("`initialize` not implemented for `$(typeof(qs))`, `$(typeof(d))`.")
end

"""
    update!(q::Query, d::Data, error::Float64)

Performs multiplicative weights update of synthetic data `d` for given query `q`
and error value `error`.
"""
function update!(q::Query, d::Data, error::Float64)
    throw("`update!` not implemented for `$(typeof(q))`, `$(typeof(d))`.")
end

"""
   normalize!(d::Data)

Normalizes data.
"""
function normalize!(d::Data)
    throw("`normalize!` not implemented for `$(typeof(d))`.")
end


text:
using .RCall

function RCall.sexp(f::RLEVector)
    RCall.R"library(GenomicRanges)"
    v = values(f)
    w = widths(f)
    RCall.R"Rle($v, $w)"
end

function rcopy(::Type{RLEVector}, s::Ptr{S4Sxp})
    RLEVector(Rcall.rcopy(s[:values]), cumsum(RCall.rcopy(s[:lengths])))
end

rcopytype(::Type{RCall.RClass{:Rle}}, s::Ptr{RCall.S4Sxp}) = RLEVector


text:
# Implementation of the Genetic algorithm BRKGA in Julia
module BRKGA

export Parameters, BRKGAProblem
export start

abstract BRKGAProblem

param(prob::BRKGAProblem) = prob.param
nalleles(prob::BRKGAProblem)  = error("Please implement nalleles function for $(typeof(x)) subtype of BRKGAProblem")
decoder(prob::BRKGAProblem)  = error("Please implement decoder function for $(typeof(x)) subtype of BRKGAProblem")

type Individual
    chromosome::Vector{Float64}
    fitness::Float64
end

function Individual(prob::BRKGAProblem)
    chr = rand(nalleles(prob))
    fit = decoder(prob, chr)
    Individual(chr, fit)
end

function Individual(prob::BRKGAProblem, chr::Vector{Float64})
    fit = decoder(prob, chr)
    Individual(chr, fit)
end

fitness(ind::Individual)          = ind.fitness
chromosome(ind::Individual)       = ind.chromosome
allele(ind::Individual, idx::Int) = ind.chromosome[idx]

type Parameters
    proballeles::Float64 # Crossover probability
    maxpop::Int64        # Max number of the population (has to be bigger than elitepop+ elitepop+ nonelitepop)
    maxgen::Int64        # Max number of generations
    elitepop::Int64      # Number of elite individuals
    randind::Int64       # Number of new random individuals inserted in each generation
    crosstype::Int64     # Crossover type
end

Parameters(proballeles::Float64, maxpop::Int64, maxgen::Int64, felite::Float64, frand::Float64, crosstype::Int64) =
 Parameters(proballeles, maxpop, maxgen, floor(Int64, maxpop*felite), floor(Int64, maxpop*frand), crosstype)

# Util function for parameters
probal(para::Parameters)    = para.proballeles
maxpop(para::Parameters)    = para.maxpop
maxgen(para::Parameters)    = para.maxgen
neli(para::Parameters)      = para.elitepop
crossind(para::Parameters)  = maxpop(para) -(neli(para) + randind(para))
randind(para::Parameters)   = para.randind
crosstype(para::Parameters) = para.crosstype

"""

Make new population crossing-over an elite individual with non-elite to construct a new individual
"""
function crosspop(prob::BRKGAProblem, pop::Vector{Individual})
    para = param(prob)
    newpop = Array(Individual, crossind(para))
    for i = 1:crossind(para)
        # Choose parent
        indpeli = floor(Int64, rand()*neli(para)) +1

        # Choose the other parent
        if crosstype(para) == 1
            # Join elite with non-elite or random
            indpneli = neli(para) + floor(Int64, rand()*(crossind(para)+randind(para))) +1
        elseif crosstype(para) == 2
            # Join elite with elite, non-elite or random
            indpneli = floor(Int64, rand()*maxpop(para)) +1
            # Test if they are the same
            if indpneli == peli
                if indpeli < neli(para)
                    indpeli += 1
                else
                    indpeli -= 1
                end
            end
        end

        peli = pop[indpeli]   # Elite parent
        pneli = pop[indpneli] # Non-elite parent

        # The first parent is always the one with less fitness
        if fitness(peli) > fitness(pneli)
            aux = peli
            peli = pneli
            pneli = aux
        end

        # Joint the parents
        chr = Array(Float64, nalleles(prob))
        for j = 1:nalleles(prob)
            r = rand()

            # Allele form the non-elite parent
            if r > probal(para)
                chr[j] = allele(pneli, j)
            else # Allele form the elite parent
                chr[j] = allele(peli, j)
            end
        end
        newpop[i] = Individual(prob, chr)
    end
    return newpop
end

"""

Create a random population
"""
function randpop(prob::BRKGAProblem, n::Int)
    para = param(prob)
    newpop = Array(Individual, n)
    # Generates a random solution
    for i = 1:n
        newpop[i] = Individual(prob)
    end
    return newpop
end

"""

Start the GA algorithm
"""
function start(prob::BRKGAProblem; targetcost=0, maxstableit=Inf)
    tic()
    startt = time()
    para = param(prob)
    pop = randpop(prob, maxpop(para))

    lastsol = Inf
    g = 0
    stableit = 0 # Iterations without improvement
    cursol = 0
    bestsoltime = 0
    bestsolution = Inf
    for g = 1:maxgen(para)
        sort!(pop, by=fitness)
        cropop = crosspop(prob, pop)
        rndpop = randpop(prob, randind(para))
        pop = vcat(pop[1:neli(para)], cropop, rndpop)
        cursol = fitness(pop[1])

        info("Generation: $g current solution: $(cursol)")

        if cursol < targetcost
            break
        end
        if isapprox(lastsol, cursol)
            info("Iteration without improvement $(stableit)")
            stableit += 1
            if stableit >= maxstableit
                break
            end
        else
            stableit = 0
        end
        if cursol < bestsolution
            bestsoltime = time() -startt
            bestsolution = cursol
        end
        lastsol = cursol
    end
    totaltime = toq()
    info("Final time $(totaltime)")
    return cursol, g, totaltime, bestsoltime
end

end # BRKGA


text:
module MultiDimDictionaries
using Reexport
@reexport using Dictionaries

import Base:
  convert,
  copy,
  keys,
  getindex,
  get,
  isassigned,
  setindex!,
  insert!,
  delete!,
  length,
  show,
  Tuple,
  eltype,
  ==,
  similar,
  vcat,
  hcat,
  hvncat
import Dictionaries:
  deletetoken!,
  issettable,
  istokenizable,
  isinsertable,
  gettokenvalue,
  gettoken!,
  merge,
  haskey,
  settokenvalue!,
  set!

include("tuple_convert.jl")
include("linearindex.jl")

_isless(n1::Integer, n2::Integer) = (n1 < n2)
_isless(n1, n2) = false

function expand_dims!(dims, key)
  N = length(key)
  if length(dims) < N
    append!(dims, zeros(N - length(dims)))
  end
  @assert length(dims) ≥ N
  for n in 1:N
    if _isless(dims[n], key[n])
      dims[n] = key[n]
    end
  end
  return dims
end

max_key_length(dictionary::Dictionary) = maximum(length.(keys(dictionary)); init=0)

function default_dims(dictionary::Dictionary)
  isempty(dictionary) && return Int[]
  N = max_key_length(dictionary)
  # Determine the dimensions automatically from the keys
  dims = fill(0, N)
  for key in keys(dictionary)
    expand_dims!(dims, key)
  end
  return dims
end

struct Private end

# TODO: implement slicing with a `to_keys` generalization of `to_indices`.
struct MultiDimDictionary{I<:Tuple,T} <: AbstractDictionary{I,T}
  dictionary::Dictionary{I,T}
  dims::Vector{Int}
  function MultiDimDictionary{I,T}(
    ::Private, dictionary::Dictionary, dims::Vector{Int}
  ) where {I<:Tuple,T}
    # TODO: `dims` may be longer than `default_dims` if
    # a key was removed. Should `dims` be adjusted
    # when keys are removed?
    # @assert all(dims .≥ default_dims(dictionary))
    return new{I,T}(dictionary, dims)
  end
end

function copy(dictionary::MultiDimDictionary{I,T}) where {I,T}
  return MultiDimDictionary{I,T}(
    Private(), deepcopy(dictionary.dictionary), copy(dictionary.dims)
  )
end

function MultiDimDictionary{I,T}(::Private, dictionary::Dictionary, dims) where {I<:Tuple,T}
  return MultiDimDictionary{I,T}(Private(), dictionary, collect(dims))
end

function MultiDimDictionary{I,T}(
  dictionary::Dictionary; dims=default_dims(dictionary)
) where {I<:Tuple,T}
  return MultiDimDictionary{I,T}(Private(), dictionary, dims)
end

function MultiDimDictionary{I}(dictionary::Dictionary; kwargs...) where {I<:Tuple}
  return MultiDimDictionary{I,eltype(dictionary)}(dictionary; kwargs...)
end

function MultiDimDictionary(dictionary::Dictionary{I,T}; kwargs...) where {I<:Tuple,T}
  return MultiDimDictionary{I,T}(dictionary; kwargs...)
end

function MultiDimDictionary(dictionary::Dictionary{I,T}; kwargs...) where {I,T}
  d = Dictionary{Tuple,T}(tuple_convert.(keys(dictionary)), dictionary)
  return MultiDimDictionary(d; kwargs...)
end

function MultiDimDictionary{I,T}(; kwargs...) where {I<:Tuple,T}
  return MultiDimDictionary{I,T}(Dictionary{I,T}(); kwargs...)
end

MultiDimDictionary(indexable) = MultiDimDictionary(Dictionary(indexable))

function MultiDimDictionary{I}(indexable) where {I<:Tuple}
  return MultiDimDictionary{I}(Dictionary(indexable))
end

function MultiDimDictionary{I,T}(indexable) where {I<:Tuple,T}
  return MultiDimDictionary{I,T}(Dictionary{I,T}(indexable))
end

MultiDimDictionary(inds, values) = MultiDimDictionary(Dictionary(inds, values))

function MultiDimDictionary{I}(inds, values) where {I<:Tuple}
  return MultiDimDictionary{I}(Dictionary{I}(inds, values))
end

function MultiDimDictionary{I,T}(inds, values) where {I<:Tuple,T}
  return MultiDimDictionary{I,T}(Dictionary{I,T}(inds, values))
end

MultiDimDictionary(; kwargs...) = MultiDimDictionary(Dictionary{Tuple}(); kwargs...)

multidimdictionary(iter) = MultiDimDictionary(dictionary(iter))

keys(dictionary::MultiDimDictionary) = keys(dictionary.dictionary)

function gettokenvalue(dictionary::MultiDimDictionary, token)
  return gettokenvalue(dictionary.dictionary, token)
end
function gettoken!(dictionary::MultiDimDictionary{I}, token::I) where {I}
  return gettoken!(dictionary.dictionary, token)
end

# Trait distinguishing if an index is a slice or just an element
abstract type IndexType end

struct SliceIndex <: IndexType end
struct ElementIndex <: IndexType end

# In general assume just an index
IndexType(::Any) = ElementIndex()

IndexType(::Colon) = SliceIndex()
IndexType(::Vector) = SliceIndex()
IndexType(::AbstractRange) = SliceIndex()

# For multi-index, check if any are SliceIndex
function IndexType(i1, i...)
  return any(x -> IndexType(x) isa SliceIndex, (i1, i...)) ? SliceIndex() : ElementIndex()
end

function getindex(dictionary::MultiDimDictionary, index::Tuple)
  return getindex(IndexType(index...), dictionary, index)
end

function getindex(dictionary::MultiDimDictionary, i...)
  return getindex(dictionary, tuple(i...))
end

function getindex(::ElementIndex, dictionary::MultiDimDictionary, index::Tuple)
  return getindex(dictionary.dictionary, index)
end

function getindex(::ElementIndex, dictionary::MultiDimDictionary, i...)
  return getindex(ElementIndex(), dictionary, tuple(i...))
end

function getindex(dictionary::MultiDimDictionary, index::LinearIndex)
  return getindex(dictionary.dictionary.values, index.I)
end

function isassigned(dictionary::MultiDimDictionary, index::Tuple)
  return isassigned(dictionary.dictionary, index)
end
isassigned(dictionary::MultiDimDictionary, i...) = isassigned(dictionary, (i...,))

haskey(dictionary::MultiDimDictionary, index::Tuple) = haskey(dictionary.dictionary, index)
haskey(dictionary::MultiDimDictionary, i...) = haskey(dictionary, (i...,))

issettable(dictionary::MultiDimDictionary) = issettable(dictionary.dictionary)

function setindex!(
  dictionary::MultiDimDictionary{I,T}, value::T, index::I
) where {I<:Tuple,T}
  expand_dims!(dictionary.dims, index)
  set!(dictionary.dictionary, index, value)
  return dictionary
end

function setindex!(dictionary::MultiDimDictionary{I,T}, value, index::Tuple) where {I,T}
  return setindex!(dictionary, convert(T, value), convert(I, index))
end

function setindex!(dictionary::MultiDimDictionary, value, i...)
  return setindex!(dictionary, value, tuple(i...))
end

function set!(dictionary::MultiDimDictionary{I,T}, index::I, value::T) where {I<:Tuple,T}
  expand_dims!(dictionary.dims, index)
  set!(dictionary.dictionary, index, value)
  return dictionary
end

function set!(dictionary::MultiDimDictionary{I,T}, index::Tuple, value) where {I,T}
  return set!(dictionary, convert(I, index), convert(T, value))
end

function setindex!(dictionary::MultiDimDictionary, value, index::LinearIndex)
  # XXX: Expand dimensions
  setindex!(dictionary.dictionary.values, value, index.I)
  return dictionary
end

isinsertable(dictionary::MultiDimDictionary) = isinsertable(dictionary.dictionary)
istokenizable(dictionary::MultiDimDictionary) = istokenizable(dictionary.dictionary)
function settokenvalue!(dictionary::MultiDimDictionary{<:Any,T}, i, value::T) where {T}
  return settokenvalue!(dictionary.dictionary, i, value)
end

function insert!(dictionary::MultiDimDictionary{I,T}, index::I, value::T) where {I<:Tuple,T}
  insert!(dictionary.dictionary, index, value)
  return dictionary
end

function insert!(dictionary::MultiDimDictionary{I,T}, index::Tuple, value) where {I,T}
  return insert!(dictionary, convert(I, index), convert(T, value))
end

function delete!(dictionary::MultiDimDictionary, index::Tuple)
  delete!(dictionary.dictionary, index)
  return dictionary
end

function delete!(dictionary::MultiDimDictionary, i...)
  return delete!(dictionary, tuple(i...))
end

function deletetoken!(dictionary::MultiDimDictionary, index::Tuple)
  deletetoken!(dictionary.dictionary, index)
  return dictionary
end

function deletetoken!(dictionary::MultiDimDictionary, i...)
  return deletetoken!(dictionary, tuple(i...))
end

#
# Slicing
#

index_in_slice(index, slice::Colon) = true
index_in_slice(index, slice::AbstractRange) = (index ∈ slice)
index_in_slice(index, slice::Vector) = (index ∈ slice)

function index_in_slice(index, slice)
  return index == slice
end

function index_in_slice(index::Tuple, index_slice::Tuple)
  if length(index) > length(index_slice)
    # This makes keys like `("X", 1, 2)` get
    # included in the indexing `dictionary["X", :]`.
    if !isa(last(index_slice), Colon)
      return false
    end
  elseif length(index) < length(index_slice)
    # This makes keys like `("X", 1)` get
    # included in the indexing `dictionary["X", 1, :]`.
    if any(x -> !isa(x, Colon), index_slice[(length(index) + 1):end])
      return false
    end
  end
  for dim in 1:min(length(index), length(index_slice))
    if !index_in_slice(index[dim], index_slice[dim])
      return false
    end
  end
  return true
end

# Drop singleton dimensions when slicing
function slice_index(index::Tuple, index_slice::Tuple)
  keep_dims = Int[]
  for dim in eachindex(index)
    if has_trailing_colon(dim, index_slice) || IndexType(index_slice[dim]) isa SliceIndex
      push!(keep_dims, dim)
    end
  end
  return index[keep_dims]
end

# Drop singleton dimensions when slicing
function slice_indices(
  ::SliceIndex, dictionary::MultiDimDictionary{I}, index_slice::Tuple
) where {I<:Tuple}
  sliced_indices = Indices{Tuple}()
  for index in keys(dictionary)
    if index_in_slice(index, index_slice)
      sliced_index = slice_index(index, index_slice)
      insert!(sliced_indices, sliced_index)
    end
  end
  return sliced_indices
end

function getindex_dropdims(
  ::SliceIndex, dictionary::MultiDimDictionary{I}, index_slice::Tuple
) where {I<:Tuple}
  indices = slice_indices_no_dropdims(SliceIndex(), dictionary, index_slice)
  sliced_indices = slice_indices(SliceIndex(), dictionary, index_slice)
  return MultiDimDictionary(sliced_indices, getindices(dictionary.dictionary, indices))
end

# Don't drop singleton dimensions when slicing
function slice_indices_no_dropdims(
  ::SliceIndex, dictionary::MultiDimDictionary{I}, index_slice::Tuple
) where {I<:Tuple}
  indices = Indices{I}()
  for index in keys(dictionary)
    if index_in_slice(index, index_slice)
      insert!(indices, index)
    end
  end
  return indices
end

function getindex(
  ::SliceIndex, dictionary::MultiDimDictionary{I}, index_slice::Tuple
) where {I<:Tuple}
  indices = slice_indices_no_dropdims(SliceIndex(), dictionary, index_slice)
  # Need convert because of:
  # https://github.com/andyferris/Dictionaries.jl/issues/97
  subdictionary = convert(
    typeof(dictionary.dictionary), getindices(dictionary.dictionary, indices)
  )
  return MultiDimDictionary(indices, subdictionary)
end

# Special version for `dictionary[[("X", 1), ("X", 2)]]`
function getindex(dictionary::MultiDimDictionary{I}, index::Vector) where {I<:Tuple}
  indices = Indices{I}()
  for key in keys(dictionary)
    if index_in_slice(key, index)
      insert!(indices, key)
    end
  end
  subdictionary = convert(
    typeof(dictionary.dictionary), getindices(dictionary.dictionary, indices)
  )
  return MultiDimDictionary(subdictionary)
end

function getindex(::SliceIndex, dictionary::MultiDimDictionary, indices...)
  return getindex(SliceIndex(), dictionary, tuple(indices...))
end

#
# Merging/disjoint unions
#

# Disjoint union without any relabelling (assumes the dictionaries already have
# unique keys).
function merge(
  dictionary1::MultiDimDictionary{I1,T1}, dictionary2::MultiDimDictionary{I2,T2}
) where {I1,T1,I2,T2}
  # TODO: Use `promote(dictionary1, dictionary2)` instead.
  I = promote_type(I1, I2)
  T = promote_type(T1, T2)
  return merge(MultiDimDictionary{I,T}(dictionary1), MultiDimDictionary{I,T}(dictionary2))
end

function merge(
  dictionary1::MultiDimDictionary{I,T}, dictionary2::MultiDimDictionary{I,T}
) where {I,T}
  dims = max.(dictionary1.dims, dictionary2.dims)
  return MultiDimDictionary(merge(dictionary1.dictionary, dictionary2.dictionary); dims)
end

_add(n1::Integer, n2::Integer) = n1 + n2
_add(n1, n2::Integer) = n2
_add(n1::Integer, n2) = n1

"""
    _shift_key(key::Tuple, dims::Tuple, dim::Int, dim_key=1)

Either shift the key in dimension `dim` by `dims[d]`, or if
`dim > length(key)` then extend the key to length `dim` by padding
with 1s and insert `dim_key` at dimension `dim`.

For example:

f = x -> x + 3

_shift_key(Tuple(2, 2, 2), [3, 3, 3], 2) -> Tuple(2, 5, 2)

_shift_key(Tuple(2, 2, 2), [3, 3, 3], 3) -> Tuple(2, 2, 5)

_shift_key(Tuple(2, 2, 2), [3, 3, 3], 4) -> Tuple(2, 2, 2, 1)

_shift_key(Tuple(2, 2, 2), [3, 3, 3], 4, "X") -> Tuple(2, 2, 2, "X")

_shift_key(Tuple(2, 2, 2), [3, 3, 3], 0) -> Tuple(1, 2, 5, 2)

_shift_key(Tuple(2, 2, 2), [3, 3, 3], 0, "X") -> Tuple("X", 2, 5, 2)

This is an "injection": https://en.wikipedia.org/wiki/Disjoint_union
"""
function _shift_key(key::Tuple, dims::Vector{Int}, dim::Int, dim_key=1)
  if dim < 0
    error("Not implemented dim=$dim")
  elseif dim == 0
    return (dim_key, key...)
  elseif dim > length(key)
    return ntuple(j -> j == dim ? dim_key : (j > length(key) ? 1 : key[j]), dim)
  end
  return ntuple(j -> j == dim ? _add(dims[dim], key[dim]) : key[j], length(key))
end

# This is a "disjoint union": https://en.wikipedia.org/wiki/Disjoint_union
# "⊔" can be typed by \sqcup<tab>
function hvncat(
  dim::Int,
  dictionary1::MultiDimDictionary,
  dictionary2::MultiDimDictionary;
  new_dim_keys=(1, 2),
)
  shifted_keys1 = map(
    key -> _shift_key(key, zero(dictionary1.dims), dim, new_dim_keys[1]), keys(dictionary1)
  )
  shifted_keys2 = map(
    key -> _shift_key(key, dictionary1.dims, dim, new_dim_keys[2]), keys(dictionary2)
  )

  ## # map doesn't narrow the types properly, so narrow them by hand
  ## narrowed_keytype1 = mapreduce(typeof, promote_type, shifted_keys1)
  ## narrowed_keytype2 = mapreduce(typeof, promote_type, shifted_keys2)

  narrowed_keytype1 = Tuple
  narrowed_keytype2 = Tuple

  shifted_dictionary1 = MultiDimDictionary{narrowed_keytype1}(
    shifted_keys1, values(dictionary1)
  )
  shifted_dictionary2 = MultiDimDictionary{narrowed_keytype2}(
    shifted_keys2, values(dictionary2)
  )
  return merge(shifted_dictionary1, shifted_dictionary2)
end

function vcat(dictionary1::MultiDimDictionary, dictionary2::MultiDimDictionary; kwargs...)
  return hvncat(1, dictionary1, dictionary2; kwargs...)
end

function hcat(dictionary1::MultiDimDictionary, dictionary2::MultiDimDictionary; kwargs...)
  return hvncat(2, dictionary1, dictionary2; kwargs...)
end

# TODO: define `disjoint_union(dictionaries...; dim::Int, new_dim_keys)` to do a disjoint union
# of a number of dictionaries.
function disjoint_union(
  dictionary1::MultiDimDictionary, dictionary2::MultiDimDictionary; dim::Int=0, kwargs...
)
  return hvncat(dim, dictionary1, dictionary2; kwargs...)
end

function ⊔(dictionary1::MultiDimDictionary, dictionary2::MultiDimDictionary; kwargs...)
  return disjoint_union(dictionary1, dictionary2; kwargs...)
end

#
# exports
#

export MultiDimDictionary, LinearIndex, disjoint_union, ⊔

end


text:
import .Plots.gui
import .Plots.plot
import .Plots.plot!
import .Plots.annotate!

function plot(env::CartPoleEnv; kwargs...)
    a, d = env.action, env.done
    x, xdot, theta, thetadot = env.state
    l = 2 * env.params.halflength
    xthreshold = env.params.xthreshold
    # set the frame
    plot(
        xlims=(-xthreshold, xthreshold),
        ylims=(-.1, l + 0.1),
        legend=false,
        border=:none,
    )
    # plot the cart
    plot!([x - 0.5, x - 0.5, x + 0.5, x + 0.5], [-.05, 0, 0, -.05];
        seriestype=:shape,
    )
    # plot the pole
    plot!([x, x + l * sin(theta)], [0, l * cos(theta)];
        linewidth=3,
    )
    # plot the arrow
    if isa(action_space(env), Base.OneTo)
        a = a == 2 ? 1 : -1
    end
    plot!([x + sign(a)*0.5, x + sign(a)*0.7], [ -.025, -.025];
        linewidth=3,
        arrow=true,
        color=:black,
    )
    # if done plot pink circle in top right (green in t > max steps)
    if d
        color = :pink
        if env.t > env.params.max_steps
            color = :green
        end
        plot!([xthreshold - 0.2], [l];
            marker=:circle,
            markersize=20,
            markerstrokewidth=0.,
            color=color,
        )
    end

    p = plot!(;kwargs...)
    gui(p)
    p
end

# adapted from https://github.com/JuliaML/Reinforce.jl/blob/master/src/envs/mountain_car.jl
height(xs) = sin(3 * xs) * 0.45 + 0.55
rotate(xs, ys, θ) = xs * cos(θ) - ys * sin(θ), ys * cos(θ) + xs * sin(θ)
translate(xs, ys, t) = xs .+ t[1], ys .+ t[2]

function plot(env::MountainCarEnv; kwargs...)
    x, v = env.state
    d = env.done

    plot(
        xlims=(env.params.min_pos - 0.1, env.params.max_pos + 0.2),
        ylims=(-.1, height(env.params.max_pos) + 0.1),
        legend=false,
        border=:none,
    )
    # plot the terrain
    xs = LinRange(env.params.min_pos, env.params.max_pos, 100)
    ys = height.(xs)
    plot!(xs, ys)

    # plot the car
    θ = cos(3 * x)
    carwidth = 0.05
    carheight = carwidth / 2
    clearance = 0.2 * carheight
    xs = [-carwidth / 2, -carwidth / 2, carwidth / 2, carwidth / 2]
    ys = [0, carheight, carheight, 0]
    ys .+= clearance
    xs, ys = rotate(xs, ys, θ)
    xs, ys = translate(xs, ys, [x, height(x)])
    plot!(xs, ys; seriestype=:shape)

     # if done plot pink circle in top right (green if reached the goal)
     if d
        color = :pink
        if x >= env.params.goal_pos && v >= env.params.goal_velocity
            color = :green
        end

        plot!([env.params.max_pos+0.1], [height(env.params.max_pos)+0.1];
            marker=:circle,
            markersize=20,
            markerstrokewidth=0.,
            color=color,
        )
    end

    # Add annotation of action on the bottom
    pos_x = (env.params.max_pos + env.params.min_pos) / 2
    pos_y = -0.05
    a = round(env.action, digits=2)
    annotate!(pos_x, pos_y, "Action: $a", :black)
    p = plot!(;kwargs...)
    gui(p)
    p
end

function plot(env::PendulumEnv; kwargs...)
    size = 1.0
    width = 0.01
    height = 1.3 * size
    s = env.state
    d = env.done
    cθ, sθ, ω = pendulum_observation(s)

    plot(
        xlims=(-size, size),
        ylims=(-size, size),
        legend=false,
        border=:none,
    )

    # Plot pendulum
    xs = [-width, -width, width, width]
    ys = [0, height, height, 0]

    plot!(xs * cθ - ys * sθ, ys * cθ + xs * sθ;
          seriestype=:shape,
          color=6,
          )

    # Plot anchor
    plot!([0.0], [0.0];
          marker=:circle,
          markersize=20,
          markerstrokewidth=2.,
          color=4,
          )

    # Plot arrow
    arrow_radius = 0.2
    torque = env.action
    xs = LinRange(2*arrow_radius/3, -2*arrow_radius/3, 100)
    if torque > 0
        xs = -xs
    end
    ys = sqrt.(arrow_radius^2 .- xs.^2)
    plot!(xs, ys;
          arrow=true,
          linewidth=2,
          color=2,
          fillalpha = abs(torque) / env.params.max_torque,
          )

    if d
      plot!([size - 0.2], [size + 0.1];
            marker=:circle,
            markersize=20,
            markerstrokewidth=0.,
            color=:pink,
        )
    end

    p = plot!(;kwargs...)
    gui(p)
    p
end


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

function naive_inputs(store::AbstractStore, get_state=Redux.get_state)
    if ReduxImGui.InputText(store, s->get_state(s).input_texts["basic_input_text"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        @info "Current string: "*ReduxImGui.get_string(get_state(store).input_texts["basic_input_text"])
    end
    CImGui.SameLine()
    ShowHelpMarker("USER:\nHold SHIFT or use mouse to select text.\n"*"CTRL+Left/Right to word jump.\n"*"CTRL+A or double-click to select all.\n"*"CTRL+X,CTRL+C,CTRL+V clipboard.\n"*"CTRL+Z,CTRL+Y undo/redo.\n"*"ESCAPE to revert.")

    if ReduxImGui.InputTextWithHint(store, s->get_state(s).input_text_with_hints["basic_input_text_with_hint"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
    end

    if ReduxImGui.InputInt(store, s->get_state(s).input_ints["basic_input_int"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_values(get_state(store).input_ints["basic_input_int"]) |> first
        @info "Current value: $v"
    end
    CImGui.SameLine()
    ShowHelpMarker("You can apply arithmetic operators +,*,/ on numerical values.\n  e.g. [ 100 ], input \'*2\', result becomes [ 200 ]\nUse +- to subtract.\n")

    if ReduxImGui.InputFloat(store, s->get_state(s).input_floats["basic_input_float"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).input_floats["basic_input_float"])
        @info "Current value: $v"
    end

    if ReduxImGui.InputDouble(store, s->get_state(s).input_doubles["basic_input_double"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).input_doubles["basic_input_double"])
        @info "Current value: $v"
    end

    if ReduxImGui.InputFloat(store, s->get_state(s).input_floats["basic_input_scientific"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).input_floats["basic_input_scientific"])
        @info "Current value: $v"
    end
    CImGui.SameLine()
    ShowHelpMarker("You can input value using the scientific notation,\n  e.g. \"1e+8\" becomes \"100000000\".\n")
end


text:
T=3
eps_exp=(1.)/(3.)
tsname="joCoreBlock - basic"
@testset "$tsname" begin
for t=1:T # start test loop

    s=rand(3:9,3)
    u=rand(7:11,2)
mo=vec([0 rand(0:u[1]) u[1]]) #rand(0:9,3)
no=vec([0 rand(0:u[2]) u[2]]) #rand(0:9,3)
m=max((mo+s)...)
n=max((no+circshift(s,-1))...)
    a=rand(ComplexF64,s[1],s[2])
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64)
    b=rand(ComplexF64,s[2],s[3])
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64)
    c=rand(ComplexF64,s[3],s[1])
    C=joMatrix(c;DDT=ComplexF32,RDT=ComplexF64)
bm=zeros(ComplexF64,m,n)
bm[(mo[1]+1):(mo[1]+s[1]),(no[1]+1):(no[1]+s[2])]+=a
bm[(mo[2]+1):(mo[2]+s[2]),(no[2]+1):(no[2]+s[3])]+=b
bm[(mo[3]+1):(mo[3]+s[3]),(no[3]+1):(no[3]+s[1])]+=c
    CB=joCoreBlock(A,B,C;moffsets=mo,noffsets=no)
    CBm=CB.m
    CBn=CB.n

    verbose && println("$tsname ($CBm,$CBn)")
    @testset "$CBm x $CBn" begin
        vn=rand(ComplexF32,CBn);
        vm=rand(ComplexF64,CBm);
        mvn=rand(ComplexF32,CBn,2);
        mvm=rand(ComplexF64,CBm,2);
        @test norm(bm*vn-CB*vn)<eps(norm(bm*vn))^eps_exp
        @test norm(bm*mvn-CB*mvn)<eps(norm(bm*mvn))^eps_exp
        @test norm(transpose(bm)*vm-transpose(CB)*vm)<eps(norm(transpose(bm)*vm))^eps_exp
        @test norm(transpose(bm)*mvm-transpose(CB)*mvm)<eps(norm(transpose(bm)*mvm))^eps_exp
        @test norm(adjoint(bm)*vm-adjoint(CB)*vm)<eps(norm(adjoint(bm)*vm))^eps_exp
        @test norm(adjoint(bm)*mvm-adjoint(CB)*mvm)<eps(norm(adjoint(bm)*mvm))^eps_exp
        @test norm(conj(bm)*vn-conj(CB)*vn)<eps(norm(conj(bm)*vn))^eps_exp
        @test norm(conj(bm)*mvn-conj(CB)*mvn)<eps(norm(conj(bm)*mvn))^eps_exp
        @test norm((-bm)*vn-(-CB)*vn)<eps(norm((-bm)*vn))^eps_exp
        @test norm((-bm)*mvn-(-CB)*mvn)<eps(norm((-bm)*mvn))^eps_exp
    end

end # end test loop
end
tsname="joCoreBlock - weighted"
@testset "$tsname" begin
for t=1:T # start test loop

    s=rand(3:9,3)
    u=rand(7:11,2)
    w=rand(ComplexF64,3)
mo=vec([0 rand(0:u[1]) u[1]]) #rand(0:9,3)
no=vec([0 rand(0:u[2]) u[2]]) #rand(0:9,3)
m=max((mo+s)...)
n=max((no+circshift(s,-1))...)
    a=rand(ComplexF64,s[1],s[2])
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64)
    b=rand(ComplexF64,s[2],s[3])
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64)
    c=rand(ComplexF64,s[3],s[1])
    C=joMatrix(c;DDT=ComplexF32,RDT=ComplexF64)
bm=zeros(ComplexF64,m,n)
bm[(mo[1]+1):(mo[1]+s[1]),(no[1]+1):(no[1]+s[2])]+=w[1]*a
bm[(mo[2]+1):(mo[2]+s[2]),(no[2]+1):(no[2]+s[3])]+=w[2]*b
bm[(mo[3]+1):(mo[3]+s[3]),(no[3]+1):(no[3]+s[1])]+=w[3]*c
    CB=joCoreBlock(A,B,C;moffsets=mo,noffsets=no,weights=w)
    CBm=CB.m
    CBn=CB.n

    verbose && println("$tsname ($CBm,$CBn)")
    @testset "$CBm x $CBn" begin
        vn=rand(ComplexF32,CBn);
        vm=rand(ComplexF64,CBm);
        mvn=rand(ComplexF32,CBn,2);
        mvm=rand(ComplexF64,CBm,2);
        @test norm(bm*vn-CB*vn)<eps(norm(bm*vn))^eps_exp
        @test norm(bm*mvn-CB*mvn)<eps(norm(bm*mvn))^eps_exp
        @test norm(transpose(bm)*vm-transpose(CB)*vm)<eps(norm(transpose(bm)*vm))^eps_exp
        @test norm(transpose(bm)*mvm-transpose(CB)*mvm)<eps(norm(transpose(bm)*mvm))^eps_exp
        @test norm(adjoint(bm)*vm-adjoint(CB)*vm)<eps(norm(adjoint(bm)*vm))^eps_exp
        @test norm(adjoint(bm)*mvm-adjoint(CB)*mvm)<eps(norm(adjoint(bm)*mvm))^eps_exp
        @test norm(conj(bm)*vn-conj(CB)*vn)<eps(norm(conj(bm)*vn))^eps_exp
        @test norm(conj(bm)*mvn-conj(CB)*mvn)<eps(norm(conj(bm)*mvn))^eps_exp
        @test norm((-bm)*vn-(-CB)*vn)<eps(norm((-bm)*vn))^eps_exp
        @test norm((-bm)*mvn-(-CB)*mvn)<eps(norm((-bm)*mvn))^eps_exp
    end

end # end test loop
end


text:
export irreps
"""
    irreps(g::FiniteGroup; χ::AbstractVector)

Get the irreducible representations.
"""
function irreps(g::AbstractFiniteGroup, χ::AbstractVector)
    isone(χ[1]) && return [χ[inclass(g, i)] * ones(1,1) for i = 1:order(g)]
    reg = regular_rep(g)
    project_out_rep(g, χ, reg)
end

function irreps(g::AbstractFiniteGroup, χ::Union{AbstractMatrix, CharacterTable})
    reg = regular_rep(g)
    [project_out_rep(g, χ[i, :], reg) for i = 1:size(χ, 1)]
end

function irreps(g::AbstractFiniteGroup, c::Characters)
    irreps(g, c.χ)
end

function irreps(g::AbstractFiniteGroup; R::Bool=false)
    R ? real_irreps(g) : irreps(g, charactertable(g))
end

#-------------------------------------------------------------------------------
# Project out irreps from regular representation
#-------------------------------------------------------------------------------
export regular_rep, proj_operator
"""
    regular_rep(g::AbstractFiniteGroup)

Return regular representation.
"""
function regular_rep(g::AbstractFiniteGroup)
    n = order(g)
    rep = Vector{SparseMatrixCSC{Int64, Int64}}(undef, n)
    all_ones = fill(1, n)
    @threads for i = 1:n
        rep[i] = sparse([g[i,j] for j=1:n], 1:n, all_ones, n, n)
    end
    rep
end
#-------------------------------------------------------------------------------
function proj_operator(r::AbstractVector{<:AbstractMatrix}, χ::AbstractVector{<:Number}, g::AbstractFiniteGroup)
    dtype = promote_type(eltype(χ), eltype.(r)...)
    m = zeros(dtype, size(r[1]))
    if length(χ) == length(class(g))
        for i = 1:order(g)
            m += conj(χ[inclass(g, i)]) * r[i]
        end
    elseif length(χ) == order(g)
        for i = 1:order(g)
            m += χ[i] * r[i]
        end
    else
        error("Invalid length $(length(χ)) for χ.")
    end
    m
end

function proj_operator(r::AbstractVector{<:AbstractMatrix}, χ::AbstractVector{<:Number})
    dtype = promote_type(eltype(χ), eltype.(r)...)
    m = zeros(dtype, size(r[1]))
    for i = 1:length(r)
        m += χ[i] * r[i]
    end
    m
end
#-------------------------------------------------------------------------------
function project_out_rep(
    g::AbstractFiniteGroup,
    χ::AbstractVector{<:Number},
    reg::AbstractVector{<:AbstractMatrix};
    tol::Real=1e-7
)
    D = Int(χ[1])
    isone(D) && return [χ[inclass(g, i)] * ones(1,1) for i = 1:order(g)]
    vs = begin
        e, v = proj_operator(reg, χ, g) |> Hermitian |> eigen
        v[:, e .> 1e-2]
    end
    vsi = vs'
    classes = class(g)
    preg = Vector{Matrix{eltype(vs)}}(undef, length(classes))
    @threads for i = 1:length(classes)
        preg[i] = vsi * reg[classes[i][1]] * vs
    end
    v0 = find_least_degen(preg, D, tol=tol)
    P = krylov_space(reg[2:order(g)], vs*v0, D)
    rep = transform_rep(P', reg, P)
    isone(check_real_rep(g, χ)) ? real_rep(rep) : rep
end
#-------------------------------------------------------------------------------
function find_least_degen(rep, D::Integer; dim::Integer=D^2, tol::Real=1e-7)
    min_degen = dim
    min_vs = nothing
    for mat in rep
        e, v = eigen(mat)
        e_split = spectrum_split(e, tol=tol)
        for ei in e_split
            if length(ei) == D
                return v[:,ei[1]]
            elseif length(ei) < min_degen
                min_degen = length(ei)
                min_vs = v[:, ei]
            end
        end
    end
    # If single representation matrix is not enough to get rid of ther degeneracy,
    # we then calculate the rep. mats. restriceted to the degenerate space.
    min_vs = svd(min_vs).U
    rep2 = (min_vs' * mat * min_vs for mat in rep)
    min_vs * find_least_degen(rep2, D, dim=min_degen, tol=tol)
end
#-------------------------------------------------------------------------------
"""
    krylov_space(mats::AbstractVector{<:AbstractMatrix}, v0::AbstractVector, n::Integer; tol)

Find the krylov space of the group action from initial vector `v0`.
"""
function krylov_space(
    mats::AbstractVector{<:AbstractMatrix},
    v0::AbstractVector,
    n::Integer;
    tol::Real=1e-3
)
    isone(n) && return reshape(v0, :, 1)
    r, i, nmat = 1, 0, length(mats)
    vs = v0
    while r < n
        i = mod(i, nmat) + 1
        vs = hcat(vs, mats[i] * vs)
        U, S = svd(vs)
        pos = S .> tol
        vs = U[:, pos]
        r = size(vs, 2)
    end
    @assert r == n "Dimension of Krylov space not match: expect d = $n, get d = $r"
    vs
end

#-------------------------------------------------------------------------------
# Realify
#-------------------------------------------------------------------------------
export check_real_rep, real_rep, real_irreps
"""
    check_real_rep(g::AbstractFiniteGroup, χ)

Check whether a representation is real/complex/pseudo-real:
+1 : Real
 0 : Complex
-1 : Pseudo-real
"""
function check_real_rep(g::AbstractFiniteGroup, χ)
    n = sum(χ[inclass(g, g[i,i])] for i=1:order(g)) / order(g)
    if abs(n-1) < 1e-7          # Real
        1
    elseif abs(n+1) < 1e-7      # Complex
        -1
    elseif abs(n) < 1e-7        # Pseudo-real
        0
    else
        error("Sum of χ(g²) = $n ≠ ±g,0")
    end
end
#-------------------------------------------------------------------------------
"""
    real_rep(r::AbstractVector{<:AbstractMatrix})

Return real representation for a complex-value real representation.

For a representation D(g), it is a real representation means that
    S⋅D̄(g)⋅S⁺ = D(g) => Sᵀ = S,
Find the sqaure root of S, i.e.,
    S = W⋅W => W⋅D̄(t)⋅W⁺ = W⁺⋅D(g)⋅W
The real representation is W⁺⋅D(g)⋅W.
"""
function real_rep(r::AbstractVector{<:AbstractMatrix})
    R = sum(kron(m, m) for m in r)
    v = eigen(R).vectors
    U = reshape(v[:, end], size(r[1])) * sqrt(size(r[1], 1))
    W = unitarysqrt(U)
    Wi = W'
    rep = Vector{Matrix{Float64}}(undef, length(r))
    Threads.@threads for i = 1:length(r)
        rep[i] = real.(Wi * r[i] * W)
    end
    rep
end
#-------------------------------------------------------------------------------
"""
    real_irreps(g::AbstractFiniteGroup)

Get all irreducible real representations of g.
"""
function real_irreps(g::AbstractFiniteGroup)
    ct = charactertable(g)
    rows = single_complex_row(ct)
    [real_irreps(g, ct[r, :]) for r in rows]
end

function real_irreps(g::AbstractFiniteGroup, χ::AbstractVector)
    rep = irreps(g, χ)
    if promote_type(typeof.(χ)...) <: Real
        rep
    else
        [[real(m) imag(m); -imag(m) real(m)] for m in rep]
    end
end
#-------------------------------------------------------------------------------
function single_complex_row(ct::CharacterTable)
    nrow = length(ct)
    row = fill(true, nrow)
    for i = 1:nrow
        row[i] || continue
        if !(eltype(ct[i]) <: Real)
            χc = conj(ct[i])
            for j = i+1:nrow
                ct[j] == χc && (row[j] = false)
            end
        end
    end
    (1:nrow)[row]
end

#-------------------------------------------------------------------------------
# Helper
#-------------------------------------------------------------------------------
export transform_rep
"""
    transform_rep(U, D1, V)

Transform the basis of the representationc `D1`.

The inputs is `U`, `D1`, and `V`, the output is:
    D₂(g) = U⋅D₁(g)⋅V
"""
function transform_rep(U::AbstractMatrix{Tu}, D1::AbstractVector{<:AbstractMatrix{Tr}}, V::AbstractMatrix{Tv}) where {Tu, Tr, Tv}
    D2 = Vector{Matrix{promote_type(Tu, Tr, Tv)}}(undef, length(D1))
    @threads for i = 1:length(D1)
        D2[i] = U * D1[i] * V
    end
    D2
end

function transform_rep(rep::AbstractVector{<:AbstractMatrix}, v::AbstractMatrix)
    if norm(v' * v - I) < 1e-10
        transform_rep(v', rep, v)
    else
        transform_rep(inv(v), rep, v)
    end
end
#-------------------------------------------------------------------------------
"""
Eigen decomposition of unitary matrix.

The eigen vectors are unitary.
"""
function unitaryeigen(U::AbstractMatrix)
    val, vec = eigen(U)
    # Require the eigen vector of a unitary matrix to be unitary.
    # Here we enforce this by orthogonalized the vectors in the degenerate space.
    for s in spectrum_split(val)
        isone(length(s)) && continue
        vec[:, s] .= svd(vec[:, s]).U
    end
    val, vec
end
#-------------------------------------------------------------------------------
function unitarysqrt(U::AbstractMatrix)
    val, vec = unitaryeigen(U)
    sval = @. sqrt(Complex(val))
    vec * Diagonal(sval) * vec'
end

#-------------------------------------------------------------------------------
# Misc
#-------------------------------------------------------------------------------
export oplus, unitary_rep, equivalent_transform
"""
Direct sum of matrices or representations.
"""
function oplus(mats::AbstractMatrix...)
    dn = size.(mats, 1)
    D = sum(dn)
    dtype = promote_type(eltype.(mats)...)
    m = zeros(dtype, D, D)
    s = 1
    for i in 1:length(mats)
        d = dn[i]
        m[s:s+d-1, s:s+d-1] .= mats[i]
        s += d
    end
    m
end

function oplus(reps::AbstractVector{<:AbstractMatrix{<:Number}}...)
    [oplus((rep[i] for rep in reps)...) for i=1:length(reps[1])]
end
#-------------------------------------------------------------------------------
"""
Convert a representation to unitary representation.

Construct the matrix H = ∑D⁺(g)D(g) = v⋅e⋅v⁺ = (X⋅X⁺)⁻¹ = (X⁻¹)⁺⋅X⁻¹
    => X = v⋅(√e)⁻¹⋅v⁺, X⁻¹ = v⋅(√e)⋅v⁺
The new representation D'(g) := X⁻¹⋅D(g)⋅X is unitary:
    D'⁺(g)⋅D'(g) = X⁺⋅D(g)⋅(X⁻¹)⁺⋅X⁻¹⋅D(g)⋅X = X⁺⋅H⋅X = I.
"""
function unitary_rep(rep::AbstractVector{<:AbstractMatrix})
    H = sum(m' * m for m in rep)
    e, v = eigen(Hermitian(H))
    es, vd = sqrt.(e), v'
    X = v * Diagonal(1 ./ es) * vd
    Xi = v * Diagonal(es) * vd
    transform_rep(Xi, rep, X)
end
#-------------------------------------------------------------------------------
"""
Return the transformation matrix between two equivalent irreducible representation
`rep1`, `rep2`, satisfying:
    rep2[i] = U⁻¹ * rep1[i] * U.
"""
function equivalent_transform(rep1::AbstractVector{<:AbstractMatrix}, rep2::AbstractVector{<:AbstractMatrix})
    D = size(rep1[1], 1)
    H = sum(kron(conj(rep2[i]), rep1[i]) for i = 1:length(rep1))
    e, v = eigen(H)
    @assert abs(e[end]) - abs(e[end-1]) > 1e-3 "Encounter possible degeneracy, implying reducibility."
    reshape(v[:, end], D, D) * sqrt(D)
end
#-------------------------------------------------------------------------------
"""
Check whether a group is legit
"""
function check_rep(g::AbstractFiniteGroup, r::AbstractVector{<:AbstractMatrix}; tol::Real=1e-7)
    Q = Matrix{Bool}(undef, order(g), order(g))
    for k = 1:order(g)
        for l = 1:order(g)
            m = g[k, l]
            Q[k, l] = norm(r[m] - r[k] * r[l]) < tol
        end
    end
    all(Q)
end
#-------------------------------------------------------------------------------
function check_unitary(r::AbstractVector{<:AbstractMatrix}; tol::Real=1e-7)
    Q = Vector{Bool}(undef, length(r))
    for i = 1:length(r)
        Q[i] = norm(r[i]' * r[i] - I) < tol
    end
    all(Q)
end


text:
@testset "Label propagation" begin
    n=10
    g10 = CompleteGraph(n)
    for g in testgraphs(g10)
      z = copy(g)
      for k=2:5
          z = blkdiag(z, g)
          add_edge!(z, (k-1)*n, k*n)
          c, ch = @inferred(label_propagation(z))
          a = collect(n:n:k*n)
          a = Int[div(i-1,n)+1 for i=1:k*n]
          # check the number of communities
          @test length(unique(a)) == length(unique(c))
          # check the partition
          @test a == c
      end
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Pkg
Pkg.activate("examples")
Pkg.instantiate()

import GPLinearODEMaker; GLOM = GPLinearODEMaker

kernel, n_kern_hyper = GLOM.include_kernel("se")

n = 100
xs = 20 .* sort(rand(n))
noise1 = 0.1 .* ones(n)
noise2 = 0.2 .* ones(n)
y1 = sin.(xs) .+ (noise1 .* randn(n))
y2 = cos.(xs) .+ (noise2 .* randn(n))

ys = collect(Iterators.flatten(zip(y1, y2)))
noise = collect(Iterators.flatten(zip(noise1, noise2)))

glo = GLOM.GLO(kernel, n_kern_hyper, 2, 2, xs, ys; noise = noise, a=[[1. 0.1];[0.1 1]])
total_hyperparameters = append!(collect(Iterators.flatten(glo.a)), [10])
workspace = GLOM.nlogL_matrix_workspace(glo, total_hyperparameters)

using Optim

initial_x = GLOM.remove_zeros(total_hyperparameters)
f(non_zero_hyper::Vector{T} where T<:Real) = GLOM.nlogL_GLOM!(workspace, glo, non_zero_hyper)
function g!(G::Vector{T}, non_zero_hyper::Vector{T}) where T<:Real
    G[:] = GLOM.∇nlogL_GLOM!(workspace, glo, non_zero_hyper)
end
function h!(H::Matrix{T}, non_zero_hyper::Vector{T}) where T<:Real
    H[:, :] = GLOM.∇∇nlogL_GLOM!(workspace, glo, non_zero_hyper)
end

function optim_cb(x::OptimizationState)
    println()
    if x.iteration > 0
        println("Iteration:     ", x.iteration)
        println("Time so far:   ", x.metadata["time"], " s")
        println("Currnet score: ", x.value)
        println("Gradient norm: ", x.g_norm)
        println()
    end
    return false
end

# @time result = optimize(f, initial_x, NelderMead()) # 26s
# @time result = optimize(f, g!, initial_x, LBFGS()) # 40s
@time result = optimize(f, g!, h!, initial_x, NewtonTrustRegion(), Optim.Options(callback=optim_cb, g_tol=1e-6, iterations=50)) # 15s

fit_total_hyperparameters = GLOM.reconstruct_total_hyperparameters(glo, result.minimizer)

n_samp_points = convert(Int64, max(500, round(2 * sqrt(2) * length(glo.x_obs))))
x_samp = collect(range(minimum(glo.x_obs); stop=maximum(glo.x_obs), length=n_samp_points))
n_total_samp_points = n_samp_points * glo.n_out
n_meas = length(glo.x_obs)

mean_GP, σ, mean_GP_obs, Σ = GLOM.GP_posteriors(glo, x_samp, fit_total_hyperparameters; return_mean_obs=true)

#=
n_show = 5
show_curves = zeros(n_show, n_total_samp_points)
L = GLOM.ridge_chol(Σ).L
for i in 1:n_show
    show_curves[i, :] = L * randn(n_total_samp_points) + mean_GP
end

using Plots

function make_plot(output::Integer, label::String; show_draws::Bool=false)
    sample_output_indices = output:glo.n_out:n_total_samp_points
    obs_output_indices = output:glo.n_out:length(ys)
    p = scatter(xs, ys[obs_output_indices], yerror=noise1, label=label)
    if show_draws
        for i in 1:n_show
            plot!(x_samp, show_curves[i, sample_output_indices], leg=false)
        end
    end
    plot!(x_samp, mean_GP[sample_output_indices]; ribbon=σ[sample_output_indices], alpha=0.3, label="GP")
    return p
end

plot(make_plot(1, "Sin"), make_plot(2, "Cos"), layout=(2,1), size=(960,540))
savefig("example/simple_ode.png")
=#

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
module SatisfiabilityInterface

using ModelingToolkit

export BoundedInteger, Model, solve, make_vector

import Base: push!

include("modeling_toolkit_interface.jl")
include("bounded_integer.jl")
include("constraints.jl")
include("model.jl")
include("sat_problem.jl")
include("solver.jl")
include("symbolic_problem.jl")

end


text:
using DICOM
function read_dicom_dir(dicomdir, filename_base="image_")
    datap = Dict()
    filename_base = "image_"
    dicomdir="C:\\Users\\Jirik/data/medical/orig/3Dircadb1.1/MASKS_DICOM/portalvein/"
    files = readdir(dicomdir)
    files_filtered = [fn for fn in files if findfirst(".", fn) == nothing ]

    function get_int(fn)
        return parse(Int, fn[length(filename_base)+1:end])
    end
    files_sorted = sort(files_filtered, by=get_int)

    dcmData = dcm_parse(dicomdir * files_sorted[1])
    image = dcmData[tag"Pixel Data"]
    pxsp = dcmData[tag"Pixel Spacing"]
    slth = dcmData[tag"Slice Thickness"]
    datap["voxelsize_mm"] = [slth, pxsp[1], pxsp[2]]
    # slloc1 = dcmData[tag"Slice Location"]
    # if length(files_sorted) > 1
    #     dcmData2 = dcm_parse(dicomdir * files_sorted[2])
    #     # slloc2 = dcmData2[tag"Slice Location"]
    #     datap["voxelsize_mm"] = [abs(slloc2 - slloc1), pxsp[1], pxsp[2]]
    #
    # else
    #     datap["voxelsize_mm"] = [pxsp[1], pxsp[2]]
    # end

    newdata = Array{eltype(image), 3}(undef, length(files_sorted), size(image,1), size(image,2))
    for i=1:length(files_sorted)
        dcmData = dcm_parse(dicomdir * files_sorted[i])
        image = dcmData[tag"Pixel Data"]
        newdata[i, :, :] = image
    end
    datap["data3d"] = newdata
    return datap
end


text:
using Base.LinAlg
import Base.LinAlg.givensAlgorithm
# Below code stolen from newer Base

"""
    lowrankupdate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

Update a Cholesky factorization `C` with the vector `v`. If `A = C[:U]'C[:U]` then `CC = cholfact(C[:U]'C[:U] + v*v')` but the computation of `CC` only uses `O(n^2)` operations. The input factorization `C` is updated in place such that on exit `C == CC`. The vector `v` is destroyed during the computation.
"""
function lowrankupdate!(C::Cholesky, v::StridedVector)
    A = C.factors
    n = length(v)
    if size(C, 1) != n
        throw(DimensionMismatch("updating vector must fit size of factorization"))
    end
    if C.uplo == 'U'
        conj!(v)
    end

    for i = 1:n

        # Compute Givens rotation
        c, s, r = givensAlgorithm(A[i,i], v[i])

        # Store new diagonal element
        A[i,i] = r

        # Update remaining elements in row/column
        if C.uplo == 'U'
            for j = i + 1:n
                Aij = A[i,j]
                vj  = v[j]
                A[i,j]  =   c*Aij + s*vj
                v[j]    = -s'*Aij + c*vj
            end
        else
            for j = i + 1:n
                Aji = A[j,i]
                vj  = v[j]
                A[j,i]  =   c*Aji + s*vj
                v[j]    = -s'*Aji + c*vj
            end
        end
    end
    return C
end

"""
    lowrankdowndate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

Downdate a Cholesky factorization `C` with the vector `v`. If `A = C[:U]'C[:U]` then `CC = cholfact(C[:U]'C[:U] - v*v')` but the computation of `CC` only uses `O(n^2)` operations. The input factorization `C` is updated in place such that on exit `C == CC`. The vector `v` is destroyed during the computation.
"""
function lowrankdowndate!(C::Cholesky, v::StridedVector)
    A = C.factors
    n = length(v)
    if size(C, 1) != n
        throw(DimensionMismatch("updating vector must fit size of factorization"))
    end
    if C.uplo == 'U'
        conj!(v)
    end

    for i = 1:n

        Aii = A[i,i]

        # Compute Givens rotation
        s = conj(v[i]/Aii)
        s2 = abs2(s)
        if s2 > 1
            throw(LinAlg.PosDefException(i))
        end
        c = sqrt(1 - abs2(s))

        # Store new diagonal element
        A[i,i] = c*Aii

        # Update remaining elements in row/column
        if C.uplo == 'U'
            for j = i + 1:n
                vj = v[j]
                Aij = (A[i,j] - s*vj)/c
                A[i,j] = Aij
                v[j] = -s'*Aij + c*vj
            end
        else
            for j = i + 1:n
                vj = v[j]
                Aji = (A[j,i] - s*vj)/c
                A[j,i] = Aji
                v[j] = -s'*Aji + c*vj
            end
        end
    end
    return C
end


text:
"""
# Strength of connection metrics for AMG

This module implements various strength of connection metrics. The strength of
connection matrix measures how closely related two vertices are. Large
values indicate that vertices are close together. Low values indicate they
are far appart. Zero values indicate that vertices are not connected.

This modules borrows heavily from pyamg. See [https://github.com/pyamg/pyamg/blob/master/pyamg/strength.py](https://github.com/pyamg/pyamg/blob/master/pyamg/strength.py).
"""
module Strength

using SparseRelaxation

"""
    algebraic_distance{T}( A :: AbstractSparseMatrix{T} # Matrix to compute distance on
                         , k :: Int # number of smoothing iterations
                         , R :: Int # number of test vectors
                         , p # norm to measure in, use `Inf` for inf norm
                         , theta :: T = 0.1 # drop tolerance for small distances
                         , smoother :: Function = weighted_jacobi! # smoother to use
                         )

Computes a strength of connection matrix using algebraic distance. This metric
takes `R` test vectors and smooths them `k` times using `smoother`. The
resulting strength of connection between `i` and `j` is the `norm` of the
difference between row `i` and `j` of the test vectors. The sparsity pattern of
`A` is mantained. `smoother` must be inplace. `theta` determines the tolerance
for dropping small distances. `theta` must be in (0, 1).
"""
function algebraic_distance{T}( A :: AbstractSparseMatrix{T}
                              , k :: Int
                              , R :: Int
                              , p :: Float64
                              , theta :: T = 0.1
                              , smoother :: Function = weighted_jacobi!
                              )
    m,n = size(A)
    tvs = rand(n, R) - 0.5 # generate test vectors
    b = zeros(n) # zero rhs

    # smooth test vectors
    for _ in 1:k
        for i in R
            smoother(A, tvs[:,i], b)
        end
    end

    xs,ys = findn(A)
    dists = if p == Inf
                maximum(tvs[xs] - tvs[ys], 2)
            else
                (sum(abs(tvs[xs] - tvs[ys]).^p, 2)./R).^(1/p)
            end

    M = sparse(xs, ys, dists) # removes zero entries were i == j

    M.nzval = 1 ./ M.nzval # invert values; large tv differences indicate large distance

    filter_distances!(M, theta) # remove small values
    M = M + speye(M)            # add identity back
    normalize_cols!(M)          # scale columns

    M
end

"""
    filter_distances!{T}(A :: AbstractSparseMatrix{T}, theta :: T)

Removes entry `e` from strength of connection if `e < theta * max(colvalues)`.
"""
function filter_distances!{T}(A :: AbstractSparseMatrix{T}, theta :: T)
    rows = rowvals(A)
    vals = nonzeros(A)
    m, n = size(A)
    @inbounds for col in 1:n
        maxval :: T = maximum(vals[nzrange(A, col)])

        for j in nzrange(A, col)
            row = rows[j]
            val = vals[j]
            if val < theta * maxval
                vals[j] = 0
            end
        end
    end

    SparseMatrix.dropzeros!(A)
end

"""
    normalize_cols!{T}(A :: AbstractSparseMatrix{T})

Scale each column in `A` by its largest value.
"""
function normalize_cols!{T}(A :: AbstractSparseMatrix{T})
    rows = rowvals(A)
    vals = nonzeros(A)
    m, n = size(A)
    @inbounds for col in 1:n
        r = nzrange(A, col)
        maxval :: T = maximum(vals[r])
        vals[r] ./ maxval
    end
end

end


text:
#PGPE nonlinearity
#out of place
function nlin(c)
    ψ = c2x(c,tinfo)
    return x2c(abs2.(ψ).*ψ,tinfo)
end

#in place
function nlin!(dc,c,tinfo)
    c2x!(ψ,c,tinfo)
    x2c!(dc,abs2.(ψ).*ψ,tinfo)
end

#dPGPE in reservoir frame
#out of place
function Lgp(c,p,t)
     return P.*(-im*(1-im*γ)*((en - μ).*c .+ g*nlin(c)))
end

#in place
function Lgp!(dc,c,p,t)
    nlin!(dc,c,tinfo)
    dc .= P.*(-im*(1-im*γ)*((en - μ).*c .+ g*dc))
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDA_Runtime_jll
using Base
using Base: UUID
using LazyArtifacts
Base.include(@__MODULE__, joinpath("..", ".pkg", "platform_augmentation.jl"))
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDA_Runtime")
JLLWrappers.@generate_main_file("CUDA_Runtime", UUID("76a88914-d11a-5bdc-97e0-2f5a05c973a2"))
end  # module CUDA_Runtime_jll


text:
@doc doc"""
The function returns a function that maps each index value to the
domain of the edge values. Note that index 0 will be mapped to the
first entry of the edge values in order to allow smoother graphics
with Gadfly.
"""->
function getInd2RangeFunc(edges::Array{Float64, 1})
    ## calculate slope
    m = (edges[end] - edges[1])./(length(edges) - 1)

    ## intercept
    c = edges[1]

    mapFunc = x -> c + m*x
    return mapFunc
end

@doc doc"""
Group observations by two-dimensional grid and apply function to
grouped observations. Internally, this makes use of function `by` from
the `DataFrames` package.
Observations in the first interval will get a location index of 0.
This will allow better plotting with Gadfly.
Edges must be monotonically increasing, and no point may lie outside
of edge range.
Observations with `NA` for either x or y value will be removed.
"""->
function localAppl(df::DataFrame, xCol::Symbol, yCol::Symbol,
                   xE::Array{Float64, 1}, yE::Array{Float64, 1},
                   aggrFunc::Function)

    ## remove NAs
    validObsInds = !isna(df[xCol]) & !isna(df[yCol])
    dfShort = df[validObsInds, :]

    if any(dfShort[xCol] .< xE[1])
        error("observations lower than smallest edge: adapt x range")
    end
    if any(dfShort[xCol] .> xE[end])
        error("observations larger than highest edge: adapt x range")
    end
    if any(dfShort[yCol] .< yE[1])
        error("observations lower than smallest edge: adapt y range")
    end
    if any(dfShort[yCol] .> yE[end])
        error("observations larger than highest edge: adapt y range")
    end

    nObs = size(dfShort, 1)

    ## preallocation
    xInds = zeros(Int64, nObs)
    yInds = zeros(Int64, nObs)

    # find respective grid interval
    for ii = 1:nObs
        ## values in first interval get mapped to index 0
        x = searchsortedfirst(xE, dfShort[ii, xCol]) - 2
        y = searchsortedfirst(yE, dfShort[ii, yCol]) - 2

        xInds[ii] = x
        yInds[ii] = y
    end
    extDf = [dfShort DataFrame(xInds = xInds, yInds = yInds)]

    # group by grid interval, apply function
    res = by(extDf, [:xInds, :yInds], aggrFunc)

    names!(res, [xCol, yCol, :value])
    return res
end

function plotLocalProperties(df::DataFrame, xCol::Symbol,
                             yCol::Symbol,
                             xEdges::Array{Float64, 1},
                             yEdges::Array{Float64, 1},
                             aggrFunc::Function;
                             xlab = string(xCol),
                             ylab = string(yCol),
                             ttl = "")

    res = localAppl(df, xCol, yCol, xEdges, yEdges, aggrFunc)

    ## get maximum indices
    xMax = maximum(res[xCol])
    yMax = maximum(res[yCol])

    ## get functions to transform inds to values
    xInd2Range = getInd2RangeFunc(xEdges)
    yInd2Range = getInd2RangeFunc(yEdges)

    currPlot = Gadfly.plot(res, x=xCol, y=yCol, color=:value, Geom.rectbin,
                           Scale.x_continuous(labels=x ->
                           @sprintf("%2.2f", xInd2Range(x)),
                                              minvalue=0, maxvalue=xMax),
                           Scale.y_continuous(labels=x ->
                           @sprintf("%2.2f", yInd2Range(x)),
                                              minvalue=0, maxvalue=yMax),
                           Guide.xlabel(xlab),
                           Guide.ylabel(ylab),
                           Guide.title(ttl))
    return currPlot
end


text:
module io

    #supported drivers
    __drivers__ = [
        "mysql",
        "mariadb",
        "sqlite"
    ]

    function driver_support(driver_name)
        return driver_name in __drivers__
    end

    function driver_families()
        result = []
        for driver_name in __drivers__
            append!(result, [driver_class(driver_name)])
        end

        return unique(result)
    end

    function driver_class(driver_name)
        result = nothing
        if driver_name in ["mysql", "mariadb"]
            result = "mysql-family"
        elseif driver_name in ["sqlite"]
            result = "sqlite-family"
        end

        return result
    end

    function create_jdbc_url(con_info::Dict)
        jdbc_url = nothing
        if con_info["driver"] == "mysql-family"
            jdbc_url = string("jdbc:mysql://", con_info["host"], ":", con_info["port"], "/", con_info["schema"], "?user=", con_info["uid"], "&password=", con_info["psswd"], con_info["args"])
        elseif con_info["driver"] == "sqlite-family"
            jdbc_url = string("jdbc:sqlite:", con_info["schema"])
        end

        return jdbc_url
    end

end


text:
#
# This example assumes below packages are available in your environment
#

using Turing, MCMCChains, Distributions, ParetoSmooth, Random
using StatsFuns

Random.seed!(5574)

data = rand(Normal(0, 1), 30)

@model function model(y)
    μ ~ Normal(0, 1)
    σ ~ truncated(Cauchy(0, 1), 0.0, Inf)

    y .~ Normal(μ, σ)
end

chains = sample(model(data), NUTS(1000, .65), MCMCThreads(), 1000, 3)

# method for MCMCChains
function pointwise_loglikes(chain::Chains, data, ll_fun)
    samples = Array(Chains(chain, :parameters).value)
    pointwise_loglikes(samples, data, ll_fun)
end

# generic method for arrays
function pointwise_loglikes(samples::Array{Float64,3}, data, ll_fun)
    n_data = length(data)
    n_samples, n_chains = size(samples)[[1,3]]
    pointwise_lls = fill(0.0, n_data, n_samples, n_chains)
    for c in 1:n_chains
        for s in 1:n_samples
            for d in 1:n_data
                pointwise_lls[d,s,c] = ll_fun(samples[s,:,c], data[d])
            end
        end
    end
    return pointwise_lls
end

function compute_loo(psis_output, pointwise_lls)
    dims = size(pointwise_lls)
    lwp = deepcopy(pointwise_lls)
    lwp += psis_output.weights;
    lwpt = reshape(lwp, dims[1], dims[2] * dims[3])';
    loos = reshape(logsumexp(lwpt; dims=1), size(lwpt, 2));
    return sum(loos)
end

# compute the pointwise log likelihoods where indices correspond to [data, sample, chain]
pointwise_lls = pointwise_loglikes(chains, data, (p,d)->logpdf(Normal(p...), d))

# compute the psis object
psis_output = psis(pointwise_lls)

# return loo based on Rob's example
loo = compute_loo(psis_output, pointwise_lls)


text:
ClockCycle = Int

abstract type DispatchCostModel end

Base.@kwdef struct DefaultDispatchCostModel <: DispatchCostModel
    skip::ClockCycle
    static_dispatch::ClockCycle
    dynamic_dispatch::ClockCycle
end
const basemodel = DefaultDispatchCostModel(
    skip                = 3,
    static_dispatch     = 8,
    dynamic_dispatch    = 2000,
)

function costof(
            typefreqs, # iterable over Frequency-s
            fixtypes,
            costmodel::DispatchCostModel = basemodel
        )
    total = 0
    for f in typefreqs
        skipcost, isstatic = calc_skipcost(f.type, f.freq, fixtypes, costmodel)
        dispatchcost = f.freq * (isstatic ? costmodel.static_dispatch
                                          : costmodel.dynamic_dispatch)
        total += skipcost + dispatchcost
    end
    return total
end

# Return (skipcost, isstatic)
function calc_skipcost(type, freq, fixtypes, costmodel)
    idx = findfirst(type, fixtypes)
    perskip_cost = costmodel.skip * freq
    if isnothing(idx)
        return perskip_cost * length(fixtypes), false
    end
    return perskip_cost * idx, true
end


text:
#
###
# Matrix Operator
###

""" 1D mat-vecs """
struct MatrixOp{T,Tm<:AbstractMatrix{T}} <: AbstractOperator{T,1}
    mat::Tm
end

@forward MatrixOp.mat (
                       Base.size, Base.adjoint, #Base.:*, Base.:\
                       issquare,
                       SciMLBase.has_ldiv, SciMLBase.has_ldiv!
                      )

function Base.:*(A::MatrixOp, u::AbstractField{<:Number,1})
    A.mat * _vec(u)
end

function Base.:\(A::MatrixOp, u::AbstractField{<:Number,1})
    A.mat \ _vec(u)
end

function LinearAlgebra.mul!(v::AbstractField{<:Number,1}, A::MatrixOp, u::AbstractField{<:Number,1})
    mul!(_vec(v), A.mat, _vec(u))
    return v
end

function LinearAlgebra.ldiv!(v::AbstractField{<:Number,1}, A::MatrixOp, u::AbstractField{<:Number,1})
    ldiv!(_vec(v), A.mat, _vec(u))
    return v
end

###
# Diagonal Operator
###

""" Diagonal Scaling Operator """
struct DiagonalOp{T,D,Tdiag<:AbstractField{T,D}} <: AbstractOperator{T,D}
    diag::Tdiag # diagonal vector
end

Base.size(A::DiagonalOp) = size(Diagonal(A.diag))
Base.adjoint(A::DiagonalOp) = A
Base.inv(A::DiagonalOp) = DiagonalOp(1 ./ A.diag)

SciMLBase.has_ldiv(::DiagonalOp) = true
SciMLBase.has_ldiv!(::DiagonalOp) = true
issquare(::DiagonalOp) = true

function Base.:*(A::DiagonalOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    Diagonal(A.diag) * _vec(u)
end

function Base.:\(A::DiagonalOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    Diagonal(A.diag) \ _vec(u)
end

function LinearAlgebra.mul!(v::AbstractField{<:Number,D}, A::DiagonalOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    mul!(_vec(v), Diagonal(A.diag), _vec(u))
    return v
end

function LinearAlgebra.ldiv!(v::AbstractField{<:Number,D}, A::DiagonalOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    ldiv!(_vec(v), Diagonal(A.diag), _vec(u))
    return v
end

function LinearAlgebra.ldiv!(A::DiagonalOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    ldiv!(Diagonal(A.diag), _vec(u))
    return u
end

# fusion
for op in (
           :+, :-, :*,
          )
    @eval function Base.$op(A::DiagonalOp, B::DiagonalOp)
        Diag = $op(Diagonal(A.diag), Diagonal(B.diag))
        DiagonalOp(Diag.diag)
    end

    @eval function Base.$op(λ::Number, A::DiagonalOp)
        diag = $op(λ, A.diag)
        DiagonalOp(diag)
    end

    @eval function Base.$op(A::DiagonalOp, λ::Number)
        diag = $op(A.diag, λ)
        DiagonalOp(diag)
    end
end

function Base.:/(A::DiagonalOp, λ::Number)
    diag = A.diag ./ λ
    DiagonalOp(diag)
end

function Base.:/(A::DiagonalOp, B::DiagonalOp)
    Diag = Diagonal(A.diag) / Diagonal(B.diag)
    DiagonalOp(Diag.diag)
end

function Base.:\(A::DiagonalOp, B::DiagonalOp)
    Diag = Diagonal(A.diag) \ Diagonal(B.diag)
    DiagonalOp(Diag.diag)
end

LinearAlgebra.:lmul!(a::Number,B::DiagonalOp) = lmul!(a,B.diag)
LinearAlgebra.:rmul!(A::DiagonalOp,b::Number) = rmul!(A.diag,b)

###
# Tensor Product Operators
###

"""
Tensor product operator

(B ⊗ A) * u
"""
function tensor_product!(V,U,A,B,cache) # 2D
    """ V .= A * U * B' """
    mul!(cache, A, U)
    mul!(V, cache, B')
end

"""
Tensor product operator

(C ⊗ B ⊗ A) * u
"""
function tensor_product!(V,U,A,B,Ct,cache1,cache2) # 3D
    szU = size(U)
    U_re = _reshape(U, (szU[1], szU[2]*szU[3]))
    mul!(cache1, A, U_re)

    # B op - write to cache2. use views

    szC = size(cache2)
    C_re = _reshape(cache2, (szC[1]*szC[2], szC[3]))
    mul!(V, C_re, Ct')

    V
end

###
# TensorProductOp2D
###

""" 2D Tensor Product Operator """
struct TensorProductOp2D{T,Ta,Tb,Tc} <: AbstractTensorProductOperator{T,2}
    A::Ta
    B::Tb

    cache::Tc
    isunset::Bool

    function TensorProductOp2D(A, B, cache = nothing, isunset = cache === nothing)
        T = promote_type(eltype(A), eltype(B))
        new{T,typeof(A),typeof(B),typeof(cache)}(A, B, cache, isunset)
    end
end

Base.size(A::TensorProductOp2D) = size(A.A) .* size(A.B)
issquare(A::TensorProductOp2D) = issquare(A.A) & issquare(A.B)

function Base.adjoint(A::TensorProductOp2D)
    if issquare(A)
        TensorProdOp(A.A', A.B', A.cache)
    else
        TensorProdOp(A.A', A.B')
    end
end

function Base.:*(A::TensorProductOp2D{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    v = copy(u)
    @set! v.array = A.A * u.array * A.B' # get type information from u
end

function init_cache(A::TensorProductOp2D, U)
    cache = A.A * U
end

function LinearAlgebra.mul!(v::AbstractField{<:Number,2}, A::TensorProductOp2D, u::AbstractField{<:Number,2})
    U = u.array
    V = v.array

    if A.isunset
        cache = init_cache(A, U)
        A = set_cache(A, cache)
        mul!(V, cache, A.B')
        return v
    end

    tensor_product!(V, U, A.A, A.B, A.cache)
    return v
end

# fusion
for op in (
           :+ , :- , :* , :/, :\,
          )
    @eval function Base.$op(A::TensorProductOp2D, B::TensorProductOp2D)
        A = $op(A.A, B.A)
        B = $op(A.B, B.B)
        TensorProductOp2D(A, B)
    end
end
#


text:
#####
## Static definitions for use with TintiNet.jl
#

# 1. Functions to pass a single sample through the network

"""
    network_pass_single_GPU(model, x::Vector{String})

    This function receives a single sequence `x` as a `Vector{String}`,
    applies a preprocessing function and sends it to the GPU, passing it through the network `model`,
    bringing the result back to CPU and returning it as a CPU object.
"""
function network_pass_single_GPU(model, x::Vector{String})

    x_gpu = fastaVocab(x)
    x_gpu = reshape(x_gpu, :, 1) |> gpu
    network_outputs = model(x_gpu) |> cpu
    CUDA.unsafe_free!(x_gpu)

    return(network_outputs)

end

"""
    network_pass_single_CPU(model, x::Vector{String})

    This function receives a single sequence `x` as a `Vector{String}`,
    applies a preprocessing function and passes through the network `model` in the CPU,
    returning the result.
"""
function network_pass_single_CPU(model, x::Vector{String})

    x = fastaVocab(x)
    x = reshape(x_gpu, :, 1)
    network_outputs = model(x_gpu)

    return(network_outputs)

end

# 2. Functions to pass a batch of samples through the network

"""
    network_pass_batched_GPU(model, x::Vector{String})

    This function receives a batch of sequences `x` as a `Vector{Vector{String}}`,
    applies a preprocessing function and sends it to the GPU, passing it through the network `model`,
    bringing the result back to CPU and returning it as a CPU object.
"""
function network_pass_batched_GPU(model, x::Vector{Vector{String}})

    x_gpu = fastaVocab(preprocess(x, 128)) |> gpu
    network_outputs = model(x_gpu) |> cpu

    return(network_outputs)

end

"""
    network_pass_single_CPU(model, x::Vector{String})

    This function receives a batch of sequences `x` as a `Vector{Vector{String}}`,
    applies a preprocessing function and passesng it through the network `model`,
    returning the network outputs.
"""
function network_pass_batched_CPU(model, x::Vector{Vector{String}})

    x_cpu = fastaVocab(preprocess(x, 128))
    network_outputs = model(x_cpu)

    return(network_outputs)

end

# 3. Function to compute classifier predictions using the network_pass functions

"""
    compute_classifier_predictions(model, sequences; batched=true, batch_size=64, use_gpu=true, sleep_time_seconds=0.02)

    This function receives a `Vector{String}` `v` and returns a same-length `Vector{String}`
    with every element of `v` padded to length `n` with symbol `p`
"""
function compute_classifier_predictions(
            model,
            sequences::Vector{Vector{String}};
            batched=false,
            eval_batch_size=64,
            use_gpu=true,
            sleep_time_seconds=0.02)

    # Pre-allocate the vector containing the prediction strings
    predictions_vector = Vector{Vector{String}}(undef, length(sequences))

    if batched # Utilize a batched network pass

        prediction_position = 1
        sequence_batches = Flux.Data.DataLoader(sequences; batchsize=eval_batch_size, shuffle=false)

        for batch in sequence_batches

            seq_lens = length.(batch)
            batch_outputs = use_gpu ? network_pass_batched_GPU(model, batch) : network_pass_batched_CPU(model, batch)
            batch_prediction = Flux.onecold(batch_outputs)

            for i in 1:length(batch)
                predictions_vector[prediction_position] = ss3Vocab[vec(batch_prediction[1:seq_lens[i], i])]
                prediction_position += 1
            end

#            sleep(sleep_time_seconds)
        end

    else # Utilize a single-sequence network pass

        for i in 1:length(sequences)
            single_prediction = use_gpu ? network_pass_single_GPU(model, sequences[i]) : network_pass_single_CPU(model, sequences[i])
            predictions_vector[i] = ss3Vocab[vec(Flux.onecold(single_prediction))]
            sleep(sleep_time_seconds)
        end

    end # end if batched

    return(predictions_vector)

end

# 3. Function to compute regressor predictions using the network_pass functions

"""
    compute_regressor_predictions(model, sequences; batched=true, batch_size=64, use_gpu=true, sleep_time_seconds=0.02)

    This function receives a `Vector{String}` `v` and returns a same-length `Vector{String}`
    with every element of `v` padded to length `n` with symbol `p`
"""
function compute_regressor_predictions(
            model,
            sequences::Vector{Vector{String}};
            batched=false,
            eval_batch_size=64,
            use_gpu=true,
            sleep_time_seconds=0.02)

    # Pre-allocate the vector containing the prediction strings
    predicted_phis_vector = Vector{Vector{Float64}}(undef, length(sequences))
    predicted_psis_vector = Vector{Vector{Float64}}(undef, length(sequences))
    predicted_accs_vector = Vector{Vector{Float64}}(undef, length(sequences))

    if batched # Utilize a batched network pass

        prediction_position = 1
        sequence_batches = Flux.Data.DataLoader(sequences; batchsize=eval_batch_size, shuffle=false)

        for batch in sequence_batches

            seq_lens = length.(batch)
            batch_outputs = use_gpu ? network_pass_batched_GPU(model, batch) : network_pass_batched_CPU(model, batch)

            for i in 1:length(batch)

                predicted_phis_vector[prediction_position] = atand.(single_prediction[1, j, 1:seq_lens[prediction_position]], single_prediction[2, j, 1:seq_lens[prediction_position]])
                predicted_psis_vector[prediction_position] = atand.(single_prediction[3, j, 1:seq_lens[prediction_position]], single_prediction[4, j, 1:seq_lens[prediction_position]])
                predicted_accs_vector[prediction_position] = (single_prediction[5, j, 1:seq_lens[prediction_position]] .* 100) .+ 100
                prediction_position += 1

            end

#            sleep(sleep_time_seconds)
        end

    else # Utilize a single-sequence network pass

        for i in 1:length(sequences)
            single_prediction = use_gpu ? network_pass_single_GPU(model, sequences[i]) : network_pass_single_CPU(model, sequences[i])

#        @show single_prediction[1,1] # For debugging

            predicted_phis_vector[i] = [ atand(single_prediction[1, j], single_prediction[2, j]) for j in 1:length(sequences[i]) ]
            predicted_psis_vector[i] = [ atand(single_prediction[3, j], single_prediction[4, j]) for j in 1:length(sequences[i]) ]
            predicted_accs_vector[i] = [ (single_prediction[5, j]*100)+100 for j in 1:length(sequences[i]) ]
            sleep(sleep_time_seconds)
        end

    end # end if batched

    return(predicted_phis_vector, predicted_psis_vector, predicted_accs_vector)

end


text:
include("Structures/TS.jl")
include("Structures/SEPP.jl")


text:
"""
    temperature_date()

Returns a vector of length 774 = 31*24 with the hourly temperature at LAX
in May 2016 (discussed on pages 259 and 266).
"""
function temperature_data()
# Hourly temperature at LAX in May 2016 in a vector of length 31*24 = 744.  Measurements are hourly at 12:53 AM to 11:53PM.
t = [
59.0, 57.9, 57.0, 57.0, 55.9, 57.0, 57.0, 60.1, 62.1, 64.0, 64.9, 66.0, 69.1, 68.0, 69.1, 66.0, 64.0, 62.1, 60.1, 59.0, 57.9, 57.9, 59.0, 57.9,  # May 1
57.9, 57.9, 57.0, 57.0, 57.0, 57.9, 57.9, 59.0, 61.0, 62.1, 64.0, 64.0, 64.0, 64.9, 64.9, 64.9, 64.0, 62.1, 61.0, 59.0, 59.0, 59.0, 59.0, 59.0,
59.0, 59.0, 57.9, 57.9, 57.9, 57.0, 59.0, 62.1, 64.0, 64.9, 66.0, 66.9, 66.9, 66.0, 68.0, 66.9, 64.0, 61.0, 60.1, 59.0, 57.9, 57.0, 57.9, 57.9,
57.0, 57.0, 57.0, 57.0, 57.9, 57.9, 57.9, 59.0, 59.0, 57.9, 59.0, 62.1, 64.0, 64.0, 64.0, 64.0, 63.0, 61.0, 61.0, 60.1, 60.1, 60.1, 60.1, 60.1,
60.1, 59.0, 59.0, 59.0, 60.1, 60.1, 60.1, 60.1, 61.0, 63.0, 63.0, 63.0, 64.9, 66.0, 68.0, 66.9, 66.9, 64.0, 63.0, 62.1, 62.1, 60.1, 60.1, 59.0,
57.9, 57.9, 59.0, 59.0, 57.9, 57.9, 59.0, 60.1, 60.1, 61.0, 63.0, 64.9, 64.9, 66.0, 66.9, 66.9, 64.9, 64.0, 62.1, 60.1, 60.1, 60.1, 60.1, 59.0,
59.0, 60.1, 59.0, 59.0, 59.0, 59.0, 57.0, 57.9, 59.0, 60.1, 59.0, 64.0, 66.0, 66.0, 66.0, 64.9, 64.9, 64.0, 62.1, 61.0, 61.0, 60.1, 60.1, 60.1,
60.1, 59.0, 59.0, 57.9, 57.9, 57.9, 57.9, 59.0, 61.0, 63.0, 64.0, 64.9, 66.0, 66.0, 64.9, 64.9, 64.0, 63.0, 62.1, 61.0, 61.0, 60.1, 60.1, 60.1,
59.0, 60.1, 59.0, 59.0, 57.9, 57.9, 59.0, 60.1, 61.0, 61.0, 62.1, 63.0, 64.0, 64.9, 64.0, 64.9, 64.0, 63.0, 61.0, 60.1, 59.0, 59.0, 59.0, 59.0,
59.0, 59.0, 57.9, 57.9, 59.0, 59.0, 60.1, 60.1, 61.0, 62.1, 63.0, 64.0, 64.9, 66.9, 66.9, 66.9, 64.9, 63.0, 62.1, 61.0, 60.1, 60.1, 60.1, 60.1,
60.1, 60.1, 59.0, 60.1, 59.0, 60.1, 60.1, 61.0, 61.0, 62.1, 63.0, 64.9, 66.9, 66.9, 66.0, 64.9, 63.0, 62.1, 61.0, 61.0, 60.1, 59.0, 59.0, 57.9,
57.9, 57.9, 57.0, 57.9, 57.9, 57.9, 57.9, 59.0, 61.0, 64.0, 64.0, 64.0, 64.9, 64.9, 64.9, 64.0, 62.1, 61.0, 61.0, 60.1, 60.1, 60.1, 60.1, 60.1,
59.0, 59.0, 60.1, 60.1, 60.1, 60.1, 60.1, 61.0, 62.1, 62.1, 63.0, 64.0, 69.1, 68.0, 68.0, 66.0, 64.9, 63.0, 63.0, 62.1, 62.1, 62.1, 62.1, 61.0,
63.0, 62.1, 61.0, 61.0, 61.0, 61.0, 61.0, 62.1, 62.1, 64.9, 64.9, 64.9, 66.9, 66.9, 66.0, 66.0, 64.9, 64.0, 62.1, 61.0, 61.0, 61.0, 61.0, 62.1,
63.0, 62.1, 62.1, 62.1, 62.1, 62.1, 62.1, 62.1, 63.0, 64.0, 66.0, 68.0, 70.0, 69.1, 68.0, 66.9, 66.9, 64.9, 64.0, 62.1, 61.0, 61.0, 61.0, 61.0,
61.0, 62.1, 62.1, 62.1, 61.0, 62.1, 62.1, 63.0, 63.0, 63.0, 64.9, 66.9, 68.0, 68.0, 68.0, 68.0, 66.9, 64.9, 64.0, 62.1, 63.0, 62.1, 62.1, 62.1,
62.1, 61.0, 61.0, 61.0, 61.0, 61.0, 62.1, 62.1, 64.0, 64.0, 66.0, 66.9, 66.9, 68.0, 68.0, 66.9, 66.0, 64.9, 63.0, 62.1, 62.1, 62.1, 61.0, 61.0,
62.1, 61.0, 61.0, 61.0, 61.0, 61.0, 61.0, 61.0, 61.0, 62.1, 62.1, 63.0, 63.0, 62.1, 63.0, 62.1, 62.1, 61.0, 61.0, 61.0, 61.0, 61.0, 61.0, 61.0,
61.0, 62.1, 61.0, 61.0, 61.0, 61.0, 61.0, 62.1, 62.1, 63.0, 64.0, 64.9, 64.9, 66.9, 66.9, 66.9, 64.9, 63.0, 62.1, 61.0, 61.0, 61.0, 61.0, 60.1,
60.1, 60.1, 60.1, 59.0, 57.9, 59.0, 59.0, 61.0, 62.1, 66.0, 64.0, 64.0, 64.0, 64.9, 64.9, 62.1, 63.0, 61.0, 61.0, 59.0, 59.0, 57.9, 57.9, 57.9,
57.0, 57.0, 55.9, 57.0, 57.0, 55.9, 57.0, 60.1, 62.1, 62.1, 64.0, 64.0, 64.9, 64.0, 64.0, 63.0, 63.0, 62.1, 61.0, 60.1, 60.1, 60.1, 59.0, 59.0,
57.9, 57.9, 57.0, 54.0, 55.0, 54.0, 57.0, 59.0, 61.0, 64.9, 66.0, 66.9, 66.9, 66.0, 66.0, 66.0, 64.9, 64.0, 64.0, 61.0, 61.0, 60.1, 60.1, 59.0,
57.9, 57.0, 57.9, 57.0, 55.9, 57.0, 57.9, 61.0, 63.0, 64.0, 66.9, 69.1, 70.0, 70.0, 70.0, 66.0, 66.0, 63.0, 62.1, 60.1, 59.0, 59.0, 57.9, 57.9,
57.9, 57.9, 57.9, 57.9, 57.9, 57.9, 59.0, 62.1, 62.1, 63.0, 64.0, 66.0, 66.0, 64.9, 64.9, 64.9, 64.9, 63.0, 61.0, 61.0, 60.1, 60.1, 60.1, 59.0,
57.9, 57.9, 57.9, 57.0, 57.0, 57.0, 57.9, 60.1, 61.0, 62.1, 63.0, 64.0, 64.0, 64.0, 64.0, 66.0, 64.9, 64.0, 62.1, 60.1, 60.1, 59.0, 57.9, 57.9,
57.9, 57.9, 57.9, 55.9, 55.9, 57.0, 59.0, 61.0, 63.0, 64.0, 66.0, 66.0, 68.0, 68.0, 66.9, 64.9, 64.0, 63.0, 63.0, 61.0, 60.1, 60.1, 60.1, 60.1,
59.0, 59.0, 59.0, 59.0, 57.9, 57.9, 60.1, 61.0, 63.0, 62.1, 64.9, 66.0, 66.9, 68.0, 68.0, 66.9, 66.0, 64.0, 61.0, 61.0, 60.1, 60.1, 60.1, 60.1,
60.1, 61.0, 61.0, 60.1, 60.1, 60.1, 60.1, 61.0, 62.1, 63.0, 63.0, 64.0, 66.0, 66.9, 66.0, 64.0, 63.0, 62.1, 61.0, 60.1, 59.0, 59.0, 59.0, 60.1,
60.1, 60.1, 60.1, 60.1, 60.1, 59.0, 60.1, 61.0, 62.1, 64.9, 64.9, 66.0, 66.9, 66.9, 66.0, 66.0, 66.0, 64.0, 62.1, 61.0, 60.1, 60.1, 60.1, 59.0,
59.0, 59.0, 57.9, 57.9, 59.0, 60.1, 60.1, 61.0, 63.0, 63.0, 64.0, 64.9, 66.9, 68.0, 66.0, 66.0, 64.9, 63.0, 62.1, 61.0, 60.1, 61.0, 60.1, 59.0,
59.0, 59.0, 59.0, 59.0, 60.1, 60.1, 61.0, 61.0, 62.1, 63.0, 63.0, 64.9, 68.0, 68.0, 68.0, 66.0, 64.9, 63.0, 62.1, 62.1, 62.1, 62.1, 62.1, 63.0 ];

return t
end

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
module testMoments
using Test
using LinearAlgebra

srcDir = dirname(dirname(@__FILE__))*"\\src\\Moments\\"
include(srcDir*"Moments.jl")
import .Moments

const la = LinearAlgebra
const mom = Moments

@testset "eᵢ" begin
    n =7
    for k in 1:n
        e = mom.eᵢ(n,k)
        @test length(e) == n
        @test maximum(e) == 1
        @test minimum(e) == 0
        @test sum(e) == 1
    end
end

@testset "var_kron" begin
    A = repeat([[1,0,1]],2,2)
    B = repeat([[0,1,0]],2,2)
    A_tens_B =  mom.var_kron(A,B)
    @test A_tens_B == repeat([[1,1,1]],4,4)

    n = rand(2:4)
    m = rand(2:5)
    k = n*m
    A = [rand(0:2,5) for i in 1:n, j in 1:n]
    B = [rand(0:2,5) for i in 1:m, j in 1:m]
    A_tens_B =  mom.var_kron(A,B)
    @test size(A_tens_B) == (k,k)

    A =  mom.make_mon_expo(4,(1,1))
    B =  mom.make_mon_expo(4,(1,1);isle = false)
    C = mom.var_kron(A,B)
    for i in 1:4, j in 1:4
        rc(k) = (1+(k-1)*4):(4+(k-1)*4)
        @assert C[rc(i),rc(j)] == [A[i,j]] .+ B
    end

    # Need more
end

@testset "make_mom_expo_mat" begin
    mon_expo_mat =  mom.make_mon_expo(3,(1,1);isle = true)
    @test la.diag(mon_expo_mat) == [ [0, 0, 0],
                                  [2, 0, 0],
                                  [0, 2, 0],
                                  [0, 0, 2]]

    n = 4
    t = 2

    M_vec =  mom.make_mon_expo(n,t)
    M_mat =  mom.make_mon_expo(n,(t,t))

    ## Complex moments

    t = 2
    d = (2,2)
    n = sum(2 .* d)
    M_vec =  mom.make_mon_expo(d,t)
    @assert length(M_vec) == binomial(n+t,t)
    M_mat =  mom.make_mon_expo(d,(t,t))
    @assert M_vec == M_mat[1,:]

end

@testset "expo_conj" begin
    @test [mom.expo_conj((2,2), 1:8)] == [[3, 4, 1, 2, 7, 8, 5, 6]]
    @test [mom.expo_conj((3,3), 1:12)] == [[4, 5, 6, 1, 2, 3, 10, 11, 12, 7, 8, 9]]
end

@testset "make_mom_expo_keys" begin
    n = 3
    t = (1,1)
    mek = mom.make_mon_expo(n,2)
    @test mek == [[0, 0, 0],
                            [1, 0, 0],
                            [0, 1, 0],
                            [0, 0, 1],
                            [2, 0, 0],
                            [1, 1, 0],
                            [1, 0, 1],
                            [0, 2, 0],
                            [0, 1, 1],
                            [0, 0, 2]]

    for i in 1:3
        n = rand(3:6,1)[1]
        t = rand(2:4,1)[1]
        mek = mom.make_mon_expo(n,t*2)
        @test size(unique(mek))[1]  == size(mek)[1]
    end
end

## Real

@testset "make_xxᵀ_tens_yyᵀ" begin
    d = (3,3)
    xxᵀtensyyᵀ = mom.make_xxᵀ_tens_yyᵀ(d)
    @test size(xxᵀtensyyᵀ) == (prod(d),prod(d))

    @test diag(xxᵀtensyyᵀ) ==  [[2, 0, 0, 2, 0, 0],
                                [2, 0, 0, 0, 2, 0],
                                [2, 0, 0, 0, 0, 2],
                                [0, 2, 0, 2, 0, 0],
                                [0, 2, 0, 0, 2, 0],
                                [0, 2, 0, 0, 0, 2],
                                [0, 0, 2, 2, 0, 0],
                                [0, 0, 2, 0, 2, 0],
                                [0, 0, 2, 0, 0, 2]]

    @test sum.(xxᵀtensyyᵀ) == repeat([4],prod(d),prod(d))

    for i in 1:5
    d = rand(2:5,2)
    xxᵀtensyyᵀ = mom.make_xxᵀ_tens_yyᵀ(d)
    halfsum(arr) = [sum(arr[1:d[1]]),sum(arr[d[1]+1:end])]
    @test halfsum.(xxᵀtensyyᵀ) == repeat([[2,2]],prod(d),prod(d))
    end
end

@testset "get_ℝ_block_diagᵀ" begin
    d = (3,3)
    t = (4,4)
    ℝ_block_diag = mom.get_ℝ_block_diag(d,t)

    @test sum([size(ℝ_block_diag[b])[1] for b ∈ keys(ℝ_block_diag)]) == binomial(sum(d)+t[1],t[1])
end

## Complex

@testset "make_xx̄ᵀ_tens_yȳᵀ" begin
    d = rand(2:5,2)
    xx̄ᵀ_tens_yȳᵀ = mom.make_xx̄ᵀ_tens_yȳᵀ(d)
end

@testset "get_ℂ_block_diagᵀ" begin

for b in keys(xx̄yȳMM_blocks)
    B = xx̄yȳMM_blocks[b]

    split_expo(ααᶥββᶥ,d) =     (ααᶥββᶥ[1:d[1]],
                                ααᶥββᶥ[d[1]+1:2*d[1]],
                                ααᶥββᶥ[1+2*d[1]:2*d[1]+d[2]],
                                ααᶥββᶥ[1+2*d[1]+d[2]:end])

    Iᵗ_temp = map(x ->  sum.(split_expo(x,d)),B)
    map(x ->  split_expo(x,d),B)
    r_list  = map(v -> v[1]-v[2],Iᵗ_temp) ; s_list  = map(v -> v[3]-v[4],Iᵗ_temp)
    @test (unique([r_list...])...,unique([s_list...])...) == b .* 2
end

pf(x) = prod(factorial.(x))
@test split_expo([1,1,1,2,2,2,3,3,3,4,4,4],(3,3)) == ([1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4])
γ,γᶥ,δ,δᶥ = split_expo([2,0,0,0,0,0,1,0,0,0,0,0],(3,3))
η,ηᶥ,ζ,ζᶥ = split_expo([0,0,0,0,0,0,2,0,0,0,0,0],(3,3))
a = ((-1)^sum(δᶥ + ζᶥ) * (im)^sum(δ+δᶥ+ζ+ζᶥ))* pf(γ + δ) * pf(γᶥ + δᶥ) * pf(η + ζ) * pf(ηᶥ + ζᶥ)
b = pf(γ)*pf(γᶥ)*pf(δ)*pf(δᶥ)*pf(η)*pf(ηᶥ)*pf(ζ)*pf(ζᶥ)

t = (2,2)
d = (2,2)
xx̄yȳMM_blocks = Moments.get_xx̄yȳMM_blocks(d,t)
MM = Moments.make_mon_expo(d,t)
size(MM)[1] == sum([size(xx̄yȳMM_blocks[b])[1] for b in [keys(xx̄yȳMM_blocks)...]])

γγᶥδδᶥζζᶥηηᶥ = Moments.get_γγᶥδδᶥζζᶥηηᶥ(d,t)
@assert size([keys(γγᶥδδᶥζζᶥηηᶥ)...])[1] == size(Moments.make_mon_expo(d,2*t[1]))[1]
for b in [keys(γγᶥδδᶥζζᶥηηᶥ)...]
    @assert [b] == unique(γγᶥδδᶥζζᶥηηᶥ[b][1] .+ γγᶥδδᶥζζᶥηηᶥ[b][2])
end

end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
using Test
using FlightMechanics
using FlightMechanics.Models
using FlightMechanics.Aircrafts

qinf = 10.
Sw = 5.
qSw = qinf*Sw
b = 100.
c = 10.
α, β = 0.0, 0.5

point = [1., 1., 1.]
fcoeff_body = [1., 2., 3.]
mcoeff_body = [-1., -2., -3.]

cfm_body = PointForcesMoments(point, fcoeff_body, mcoeff_body)
pfm_body = PointForcesMoments(point, qSw*fcoeff_body, qSw*[b, c, b].*mcoeff_body)

cfm_wind = PointForcesMoments(point,
                              body2wind(fcoeff_body..., α, β),
                              body2wind([b, c, b].*mcoeff_body..., α, β)./[b, c, b]
                              )
pfm_wind = PointForcesMoments(point,
                              body2wind(qSw*fcoeff_body..., α, β),
                              body2wind(qSw*[b, c, b].*mcoeff_body..., α, β)
                              )

aero_base = F16Aerodynamics(pfm_wind, cfm_wind, pfm_body, cfm_body)

aero1 = aerodynamics_from_body_total(pfm_body, α, β, qinf, Sw, b, c, aero_base)
aero2 = aerodynamics_from_body_coeff(cfm_body, α, β, qinf, Sw, b, c, aero_base)
aero3 = aerodynamics_from_wind_total(pfm_wind, α, β, qinf, Sw, b, c, aero_base)
aero4 = aerodynamics_from_wind_coeff(cfm_wind, α, β, qinf, Sw, b, c, aero_base)

for (ii, aero) in enumerate([aero1, aero2, aero3, aero4])
    @test isapprox(get_wind_pfm(aero), get_wind_pfm(aero_base), rtol=1e-5)
    @test isapprox(get_wind_adim_pfm(aero), get_wind_adim_pfm(aero_base), rtol=1e-5)
    @test isapprox(get_body_pfm(aero), get_body_pfm(aero_base), rtol=1e-5)
    @test isapprox(get_body_adim_pfm(aero), get_body_adim_pfm(aero_base), rtol=1e-5)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 8ff88bf0-14cb-11ec-2f1c-8f7745e8a6b2
begin
    using Pkg
    Pkg.activate("..")
    using Revise, Javis, JUtils, LaTeXStrings
end

# ╔═╡ b3c04a1d-41d3-458f-8624-3cb6c114f210
function ground(args...)
    background("black")
    sethue("white")
end

# ╔═╡ 7ac8f5b3-fa39-4e15-841b-5093437afbce
L"\huge{\theta(t) = \pi e^{-\frac{t}{\tau}} \cos(t\omega)}"

# ╔═╡ 4bd25c64-0913-4720-81d9-bd197def92ed
begin
    pendulum_video = Video(800, 800)
    Background(1:500, ground)
    JUtils.pendulum(first_frame = 1, last_frame = 500)
    v = render(pendulum_video, pathname = "../output/pendulum.gif", framerate = 30)
end

# ╔═╡ Cell order:
# ╠═8ff88bf0-14cb-11ec-2f1c-8f7745e8a6b2
# ╠═b3c04a1d-41d3-458f-8624-3cb6c114f210
# ╠═7ac8f5b3-fa39-4e15-841b-5093437afbce
# ╠═4bd25c64-0913-4720-81d9-bd197def92ed