text:
@doc raw"""
    cudnnMultiHeadAttnForward(weights, queries, keys, values; o...)
    cudnnMultiHeadAttnForward!(out, weights, queries, keys, values; o...)
    cudnnMultiHeadAttnForward(weights, queries, keys, values, d::cudnnAttnDescriptor; o...)
    cudnnMultiHeadAttnForward!(out, weights, queries, keys, values, d::cudnnAttnDescriptor; o...)

Return the multi-head attention result with `weights`, `queries`, `keys`, and `values`,
overwriting `out` if provided, according to keyword arguments or the attention descriptor
`d`.  The multi-head attention model can be described by the following equations:

```math
\begin{aligned}
&h_i = (W_{V,i} V) \operatorname{softmax}(\operatorname{smScaler}(K^T W^T_{K,i}) (W_{Q,i} q))
&\operatorname(MultiHeadAttn)(q,K,V,W_Q,W_K,W_V,W_O) = \sum_{i=1}^{\operatorname{nHeads}-1} W_{O,i} h_i
\end{aligned}
```

The input arguments are:
* `out`: Optional output tensor.
* `weights`: A weight buffer that contains ``W_Q, W_K, W_V, W_O``.
* `queries`: A query tensor ``Q`` which may contain a batch of queries (the above equations were for a single query vector ``q`` for simplicity).
* `keys`: The keys tensor ``K``.
* `values`: The values tensor ``V``.

Keyword arguments describing the tensors:
* `axes::Vector{cudnnSeqDataAxis_t} = [CUDNN_SEQDATA_VECT_DIM, CUDNN_SEQDATA_BATCH_DIM, CUDNN_SEQDATA_TIME_DIM, CUDNN_SEQDATA_BEAM_DIM]`: an array of length 4 that specifies the role of (Julia) dimensions. VECT has to be the first dimension, all 6 permutations of the remaining three are supported.
* `seqLengthsQO::Vector{<:Integer}`: sequence lengths in the queries and out containers. By default sequences are assumed to be full length of the TIME dimension.
* `seqLengthsKV::Vector{<:Integer}`: sequence lengths in the keys and values containers. By default sequences are assumed to be full length of the TIME dimension.

Keyword arguments describing the attention operation when `d` is not given:
* `attnMode::Unsigned = CUDNN_ATTN_QUERYMAP_ALL_TO_ONE | CUDNN_ATTN_DISABLE_PROJ_BIASES`: bitwise flags indicating various attention options. See cudnn docs for details.
* `nHeads::Integer = 1`: number of attention heads.
* `smScaler::Real = 1`: softmax smoothing (`1.0 >= smScaler >= 0.0`) or sharpening (`smScaler > 1.0`) coefficient. Negative values are not accepted.
* `mathType::cudnnMathType_t = math_mode()`: NVIDIA Tensor Core settings.
* `qProjSize, kProjSize, vProjSize, oProjSize`: vector lengths after projections, set to 0 by default which disables projections.
* `qoMaxSeqLength::Integer`: largest sequence length expected in queries and out, set to their TIME dim by default.
* `kvMaxSeqLength::Integer`: largest sequence length expected in keys and values, set to their TIME dim by default.
* `maxBatchSize::Integer`: largest batch size expected in any container, set to the BATCH dim of queries by default.
* `maxBeamSize::Integer`: largest beam size expected in any container, set to the BEAM dim of queries by default.

Other keyword arguments:
* `residuals = nothing`: optional tensor with the same size as queries that can be used to implement residual connections (see figure in cudnn docs). When residual connections are enabled, the vector length in `queries` should match the vector length in `out`, so that a vector addition is feasible.
* `currIdx::Integer = -1`: Time-step (0-based) in queries to process. When the `currIdx` argument is negative, all ``Q`` time-steps are processed. When `currIdx` is zero or positive, the forward response is computed for the selected time-step only. The latter input can be used in inference mode only, to process one time-step while updating the next attention window and ``Q``, ``K``, ``V`` inputs in-between calls.
* `loWinIdx, hiWinIdx::Array{Cint}`: Two host integer arrays specifying the start and end (0-based) indices of the attention window for each ``Q`` time-step. The start index in ``K``, ``V`` sets is inclusive, and the end index is exclusive. By default set at 0 and `kvMaxSeqLength` respectively.
"""
cudnnMultiHeadAttnForward, cudnnMultiHeadAttnForward!

# The axes argument in the functions below specifies the role of the axes using Julia order: VECT,BATCH,TIME,BEAM by default. Missing trailing dims assumed 1.
const cudnnSeqDataDefaultAxes = [CUDNN_SEQDATA_VECT_DIM, CUDNN_SEQDATA_BATCH_DIM, CUDNN_SEQDATA_TIME_DIM, CUDNN_SEQDATA_BEAM_DIM]

# Public methods
cudnnMultiHeadAttnForward(w,q,k,v; o...) = cudnnMultiHeadAttnForward(w,q,k,v, cudnnAttnDescriptor(q,k,v;o...); o...)
cudnnMultiHeadAttnForward(w,q,k,v, d::cudnnAttnDescriptor; o...) = cudnnMultiHeadAttnForward!(cudnnAttnOutput(q,k,v,d), w,q,k,v,d; o...)
cudnnMultiHeadAttnForward!(out, w,q,k,v; o...) = cudnnMultiHeadAttnForward!(out, w,q,k,v, cudnnAttnDescriptor(q,k,v;o...); o...)

function cudnnMultiHeadAttnForward!(
    out, weights, queries, keys, values, attnDesc::cudnnAttnDescriptor;

    # Input tensor descriptors
    axes::Vector{cudnnSeqDataAxis_t} = cudnnSeqDataDefaultAxes,
    seqLengthsQO::Vector{<:Integer} = fill(Cint(sdim(queries,axes,CUDNN_SEQDATA_TIME_DIM)), sdim(queries,axes,CUDNN_SEQDATA_BATCH_DIM)*sdim(queries,axes,CUDNN_SEQDATA_BEAM_DIM)),
    seqLengthsKV::Vector{<:Integer} = fill(Cint(sdim(keys,axes,CUDNN_SEQDATA_TIME_DIM)), sdim(keys,axes,CUDNN_SEQDATA_BATCH_DIM)*sdim(keys,axes,CUDNN_SEQDATA_BEAM_DIM)),
    devSeqLengthsQO::CuVector{Cint} = convert(CuVector{Cint}, seqLengthsQO),
    devSeqLengthsKV::CuVector{Cint} = convert(CuVector{Cint}, seqLengthsKV),
    qDesc::cudnnSeqDataDescriptor = cudnnSeqDataDescriptor(queries; axes, seqLengthArray=seqLengthsQO),
    kDesc::cudnnSeqDataDescriptor = cudnnSeqDataDescriptor(keys;    axes, seqLengthArray=seqLengthsKV),
    vDesc::cudnnSeqDataDescriptor = cudnnSeqDataDescriptor(values;  axes, seqLengthArray=seqLengthsKV),
    oDesc::cudnnSeqDataDescriptor = cudnnSeqDataDescriptor(out;     axes, seqLengthArray=seqLengthsQO),

    # forw parameters
    residuals = nothing,
    currIdx::Integer = -1,
    loWinIdx::Union{Array{Cint},Nothing} = nothing,
    hiWinIdx::Union{Array{Cint},Nothing} = nothing,
    workspace::Union{CuArray,Nothing}    = nothing,
    reserveSpace::Union{CuArray,Nothing} = nothing,

    # Buffers for gradients
    dweights::Ref = Ref{Any}(),
    dqueries::Ref = Ref{Any}(),
    dkeys::Ref    = Ref{Any}(),
    dvalues::Ref  = Ref{Any}(),
    o...
)
    d = cudnnGetAttnDescriptor(attnDesc)
    dt = juliaDataType(d.dataType)
    @assert dt == eltype(out) == eltype(queries) == eltype(keys) == eltype(values)
    qSize = (d.qProjSize > 0 ? d.qProjSize : size(queries,1))
    kSize = (d.kProjSize > 0 ? d.kProjSize : size(keys,1))
    @assert kSize == qSize  "key size $kSize does not match query size $qSize."
    vSize = (d.vProjSize > 0 ? d.vProjSize : size(values,1))
    @assert size(keys)[2:end] == size(values)[2:end]  "key tensor $(size(keys)) does not match value tensor $(size(values))"
    oSize = (d.oProjSize > 0 ? d.oProjSize : d.nHeads * vSize)
    oDims = (oSize, size(queries)[2:end]...)
    @assert size(out) == oDims  "output size should be $(oDims)"
    @assert residuals === nothing || size(residuals) == oDims  "residual size should be $(oDims)"
    loWinIdx === nothing ? loWinIdx = fill(Cint(0), d.qoMaxSeqLength) : @assert length(loWinIdx) == d.qoMaxSeqLength
    hiWinIdx === nothing ? hiWinIdx = fill(typemax(Cint), d.qoMaxSeqLength) : @assert length(hiWinIdx) == d.qoMaxSeqLength

    @assert axes[1] == CUDNN_SEQDATA_VECT_DIM  "The most inner dimension of the containers should be the vector dimension"
    @assert d.smScaler >= 0  "smScaler should be non-negative"
    @assert d.qoMaxSeqLength >= sdim(queries, axes, CUDNN_SEQDATA_TIME_DIM)
    @assert d.kvMaxSeqLength >= sdim(keys,    axes, CUDNN_SEQDATA_TIME_DIM)
    @assert d.maxBatchSize   >= sdim(queries, axes, CUDNN_SEQDATA_BATCH_DIM)
    @assert d.maxBeamSize    >= sdim(queries, axes, CUDNN_SEQDATA_BEAM_DIM)
    @assert sdim(keys, axes, CUDNN_SEQDATA_BATCH_DIM) == sdim(queries, axes, CUDNN_SEQDATA_BATCH_DIM)  "keys/values and queries have different batch sizes"
    if d.attnMode & CUDNN_ATTN_QUERYMAP_ONE_TO_ONE > 0
        @assert sdim(keys, axes, CUDNN_SEQDATA_BEAM_DIM) == sdim(queries, axes, CUDNN_SEQDATA_BEAM_DIM)  "keys/values and queries have different beam sizes when attnMode is CUDNN_ATTN_QUERYMAP_ONE_TO_ONE"
    else
        @assert sdim(keys, axes, CUDNN_SEQDATA_BEAM_DIM) == 1  "keys/values should have beam=1 when attnMode is CUDNN_ATTN_QUERYMAP_ALL_TO_ONE"
    end

    # Backward called separately on each variable. We will calculate all gradients on first call. Use `dready` to avoid subsequent calls.
    dready = Ref{Bool}(false)   # this will be turned to `true` by the first backward call.

    cudnnMultiHeadAttnForwardAD(
        weights, queries, keys, values, residuals;
        dready, dweights, dqueries, dkeys, dvalues, # dresiduals is equal to dout
        attnDesc, currIdx, loWinIdx, hiWinIdx,
        devSeqLengthsQO, devSeqLengthsKV,
        qDesc, kDesc, vDesc, oDesc,
        out, workspace, reserveSpace)
end

# AD method
function cudnnMultiHeadAttnForwardAD(
    weights, queries, keys, values, residuals;
    dready, dweights, dqueries, dkeys, dvalues,
    attnDesc, currIdx, loWinIdx, hiWinIdx,
    devSeqLengthsQO, devSeqLengthsKV,
    qDesc, kDesc, vDesc, oDesc,
    out, workspace, reserveSpace
)
    # Cannot use @workspace here because it is shared between forw and back calls
    (weightSize, workspaceSize, reserveSpaceSize) = cudnnMultiHeadAttnBuffers(attnDesc)
    if workspaceSize > 0 && workspace === nothing; workspace = cudnnTempSpace(workspaceSize); end
    if reserveSpaceSize > 0 && reserveSpace === nothing; reserveSpace = cudnnTempSpace(reserveSpaceSize); end
    @assert sizeof(weights) >= weightSize  "weights should be at least $weightSize bytes."
    @assert sizeof(workspace) >= workspaceSize  "worksSpace should be at least $workspaceSize bytes"
    @assert sizeof(reserveSpace) >= reserveSpaceSize  "reserveSpace should be at least $reserveSpaceSize bytes"

    cudnnMultiHeadAttnForward(
        handle(), attnDesc, currIdx,
        loWinIdx, hiWinIdx,
        devSeqLengthsQO, devSeqLengthsKV,
        qDesc, queries, something(residuals, CU_NULL),
        kDesc, keys,
        vDesc, values,
        oDesc, out,
        sizeof(weights), something(weights, CU_NULL),
        sizeof(workspace), something(workspace, CU_NULL),
        sizeof(reserveSpace), something(reserveSpace, CU_NULL)
    )
    return out
end

# Helper methods

function cudnnAttnDescriptor(
    queries, keys, values;
    axes = cudnnSeqDataDefaultAxes,
    attnMode::Unsigned = CUDNN_ATTN_QUERYMAP_ALL_TO_ONE | CUDNN_ATTN_DISABLE_PROJ_BIASES |> Cuint,
    nHeads::Integer = Cint(1),
    smScaler::Real = Cdouble(1),
    # dataType::DataType = eltype(queries),
    # computePrec::DataType = eltype(queries),  ## No other option according to 8.0.2
    mathType::cudnnMathType_t = math_mode(),
    # attnDropout::Real = 0, ## The dropout option is currently not supported by the multi-head attention API
    # postDropout::Real = 0, ## The dropout option is currently not supported by the multi-head attention API
    qProjSize::Integer = 0, # Use zero to disable the corresponding projection
    kProjSize::Integer = 0,
    vProjSize::Integer = 0,
    oProjSize::Integer = 0,
    qoMaxSeqLength::Integer = sdim(queries,axes,CUDNN_SEQDATA_TIME_DIM),
    kvMaxSeqLength::Integer = sdim(keys,axes,CUDNN_SEQDATA_TIME_DIM),
    maxBatchSize::Integer = sdim(queries,axes,CUDNN_SEQDATA_BATCH_DIM),
    maxBeamSize::Integer = sdim(queries,axes,CUDNN_SEQDATA_BEAM_DIM),
    o...
)
    cudnnAttnDescriptor(
        Cuint(attnMode),
        Cint(nHeads),
        Cdouble(smScaler),
        cudnnDataType(eltype(queries)),    # dataType
        cudnnDataType(eltype(queries)),    # computePrec
        mathType,
        C_NULL,  # attnDropout
        C_NULL,  # postDropout
        Cint(sdim(queries,axes,CUDNN_SEQDATA_VECT_DIM)), # qSize
        Cint(sdim(keys,   axes,CUDNN_SEQDATA_VECT_DIM)), # kSize
        Cint(sdim(values, axes,CUDNN_SEQDATA_VECT_DIM)), # vSize
        Cint(qProjSize),
        Cint(kProjSize),
        Cint(vProjSize),
        Cint(oProjSize),
        Cint(qoMaxSeqLength),
        Cint(kvMaxSeqLength),
        Cint(maxBatchSize),
        Cint(maxBeamSize)
    )
end

function cudnnGetAttnDescriptor(d::cudnnAttnDescriptor)
    (attnMode, nHeads, smScaler, dataType, computePrec, mathType, attnDropoutDesc, postDropoutDesc, qSize, kSize, vSize, qProjSize, kProjSize, vProjSize, oProjSize, qoMaxSeqLength, kvMaxSeqLength, maxBatchSize, maxBeamSize) = (Ref{Cuint}(), Ref{Cint}(), Ref{Cdouble}(), Ref{cudnnDataType_t}(), Ref{cudnnDataType_t}(), Ref{cudnnMathType_t}(), Ref{cudnnDropoutDescriptor_t}(), Ref{cudnnDropoutDescriptor_t}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}())
    cudnnGetAttnDescriptor(d, attnMode, nHeads, smScaler, dataType, computePrec, mathType, attnDropoutDesc, postDropoutDesc, qSize, kSize, vSize, qProjSize, kProjSize, vProjSize, oProjSize, qoMaxSeqLength, kvMaxSeqLength, maxBatchSize, maxBeamSize)
    (attnMode, nHeads, smScaler, dataType, computePrec, mathType, attnDropoutDesc, postDropoutDesc, qSize, kSize, vSize, qProjSize, kProjSize, vProjSize, oProjSize, qoMaxSeqLength, kvMaxSeqLength, maxBatchSize, maxBeamSize) = (x->x[]).((attnMode, nHeads, smScaler, dataType, computePrec, mathType, attnDropoutDesc, postDropoutDesc, qSize, kSize, vSize, qProjSize, kProjSize, vProjSize, oProjSize, qoMaxSeqLength, kvMaxSeqLength, maxBatchSize, maxBeamSize))
    return (; attnMode, nHeads, smScaler, dataType, computePrec, mathType, attnDropoutDesc, postDropoutDesc, qSize, kSize, vSize, qProjSize, kProjSize, vProjSize, oProjSize, qoMaxSeqLength, kvMaxSeqLength, maxBatchSize, maxBeamSize)
end

function cudnnAttnOutput(queries, keys, values, attnDesc::cudnnAttnDescriptor)
    d = cudnnGetAttnDescriptor(attnDesc)
    vSize = (d.vProjSize > 0 ? d.vProjSize : size(values,1))
    oSize = (d.oProjSize > 0 ? d.oProjSize : d.nHeads * vSize)
    oDims = (oSize, size(queries)[2:end]...)
    out = similar(values, oDims)
    out .= 0 # currIdx >= 0 only fills part of this, zero the rest for consistency
    return out
end

function cudnnMultiHeadAttnBuffers(attnDesc::cudnnAttnDescriptor; training=false)
    weightSize, workspaceSize = Ref{Csize_t}(0), Ref{Csize_t}(0)
    # Assigning NULL to the reserveSpaceSizeInBytes argument indicates that the user does not plan to invoke multi-head attention gradient functions
    reserveSpaceSize = training ? Ref{Csize_t}(0) : C_NULL
    cudnnGetMultiHeadAttnBuffers(handle(), attnDesc, weightSize, workspaceSize, reserveSpaceSize)
    return (weightSize[], workspaceSize[], reserveSpaceSize === C_NULL ? 0 : reserveSpaceSize[])
end

# If we have less than four dims, assume they are padded by 1s on the right for Julia, on the left for cudnn.
# We keep `axes` in Julia order, e.g. axes[1] refers to the function of the first Julia dimension and is always VECT.
"""
    sdim(x,axes,dim)
    sdim(x,axes)

The first form returns the size of `x` in the dimension specified with
`dim::cudnnSeqDataAxis_t` (e.g. CUDNN_SEQDATA_TIME_DIM), i.e. return `size(x,i)` such that
`axes[i]==dim`.

The second form returns an array of length 4 `dims::Vector{Cint}` such that `dims[1+dim] ==
sdim(x,axes,dim)` where `dim::cudnnSeqDataAxis_t` specifies the role of the dimension
(e.g. dims[CUDNN_SEQDATA_TIME_DIM]==5).

The `axes::Vector{cudnnSeqDataAxis_t}` argument is an array of length 4 that
specifies the role of Julia dimensions, e.g. `axes[3]=CUDNN_SEQDATA_TIME_DIM`.
"""
function sdim(x,axes,dim)
    for i in 1:length(axes)
        if axes[i] === dim      # axes[i] = CUDNN_SEQDATA_XXX_DIM
            return size(x,i)
        end
    end
    error("Cannot find $dim in axes")
end

function sdim(x,axes)
    dims = Array{Cint}(undef, 4)
    for dim in (CUDNN_SEQDATA_VECT_DIM, CUDNN_SEQDATA_BATCH_DIM, CUDNN_SEQDATA_TIME_DIM, CUDNN_SEQDATA_BEAM_DIM)
        dims[1+dim] = sdim(x,axes,dim)
    end
    return dims                 # dims[1+CUDNN_SEQDATA_XXX_DIM] = how many XXX
end

# Alternative cudnnSeqDataDescriptor constructor for array
function cudnnSeqDataDescriptor(
    array;
    axes::Vector{cudnnSeqDataAxis_t} = cudnnSeqDataDefaultAxes,
    dimA::Vector{Cint} = sdim(array,axes),
    seqLengthArray::Vector{<:Integer} = fill(Cint(sdim(array,axes,CUDNN_SEQDATA_TIME_DIM)), sdim(array,axes,CUDNN_SEQDATA_BATCH_DIM)*sdim(array,axes,CUDNN_SEQDATA_BEAM_DIM)), # cudnn-doc: The seqLengthArray[] must specify all sequence lengths in the container so the total size of this array should be dimA[CUDNN_SEQDATA_BATCH_DIM] * dimA[CUDNN_SEQDATA_BEAM_DIM].
    paddingFill::Ptr{Cvoid} = C_NULL, # cudnn-doc: Currently, the only supported value for paddingFill is NULL which means this option should be ignored.
)
    nbDims::Cint = CUDNN_SEQDATA_DIM_COUNT          # Currently, the value of this argument should be four. The actual size of the dimA[] and axes[] arrays should be declared using the CUDNN_SEQDATA_DIM_COUNT macro.
    @assert length(axes) == length(dimA) == CUDNN_SEQDATA_DIM_COUNT # cudnn-doc: The number of active dimensions in the dimA[] and axes[] arrays is defined by the nbDims argument.
    seqLengthArraySize = Csize_t(sdim(array,axes,CUDNN_SEQDATA_BATCH_DIM) * sdim(array,axes,CUDNN_SEQDATA_BEAM_DIM))
    @assert length(seqLengthArray) == seqLengthArraySize
    cudnnSeqDataDescriptor(cudnnDataType(eltype(array)), nbDims, dimA, reverse(axes), # cudnn uses reverse order for dims
                           seqLengthArraySize, convert(Vector{Cint}, seqLengthArray),
                           paddingFill)
end

cudnnSeqDataDescriptor(::Nothing; o...) = nothing


text:
#Simulate a determinantally-thinned Poisson point process on a rectangle
#Author: H. Paul Keeler, 2019.

#Note: Need the .+ for adding a scalar to an array
#Also need . for sqrt, exp, cos, sin etc and assigning scalars to arrays
#Index and Boolean arrrays need to be vectors eg v=zeros(n),  NOT v=zeros(n,1)

#clearconsole(); #for clearing Julia REPL console

using Distributions #for random simulations
using PyPlot #uses Python for plotting
using LinearAlgebra;

#include some helper functions
include("funSimSimpleLDPP.jl"); #for simulations
include("funLtoK.jl"); #for converting a non-singular L matrix to a singular matrix

#set random seed for reproducibility
#rng(1);

numbSim=10^4; #number of simulations

###START -- Parameters -- ###START
#Poisson point process parameters
lambda=10; #intensity (ie mean density) of the Poisson process

#choose kernel
choiceKernel=1; #1 for Gaussian (ie squared exponetial );2 for Cauchy
sigma=1;# parameter for Gaussian and Cauchy kernel
alpha=1;# parameter for Cauchy kernel

#Simulation window parameters
xMin=0;xMax=1;yMin=0;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; #rectangle dimensions
areaTotal=xDelta*yDelta; #area of rectangle
###END -- Parameters -- ###END

### START -- Simulate a Poisson point process on a rectangle ###START
#Simulate Poisson point process
numbPoints=rand(Poisson(areaTotal*lambda)); #Poisson number of points
xx=xDelta*(rand(numbPoints)).+xMin;#x coordinates of Poisson points
yy=xDelta*(rand(numbPoints)).+yMin;#y coordinates of Poisson points
### END -- Simulate a Poisson point process on a rectangle --END ###

#numbPoints=5; xx=-2:1:2; yy=3:7 #TEMP

### START -- CREATE L matrix -- START ###
#all squared distances of x/y difference pairs
xxDiff=kron(xx,ones(1,numbPoints))-kron(ones(numbPoints,1),transpose(xx));
yyDiff=kron(yy,ones(1,numbPoints))-kron(ones(numbPoints,1),transpose(yy));
rrDiffSquared=(xxDiff.^2+yyDiff.^2);
if choiceKernel==1
    ##Gaussian/squared exponential kernel
    L=lambda*exp.(-(rrDiffSquared)/sigma^2);
elseif choiceKernel==2
    ##Cauchy kernel
    L=lambda./(1+rrDiffSquared/sigma^2).^(alpha+1/2);
else
    println("choiceKernel has to be equal to 1 or 2.");
end

### END-- CREATE L matrix -- ### END
L=Symmetric(L); #convert to Symmetric marix

### START Testing DPP simulation START ###
#Retrieve eigenvalues and eigenvectors
eigenVectL=eigvecs(L); #eigen decomposition -- vectors
eigenValL=eigvals(L); #eigen decomposition -- values

#run simulations with tests
global probX_i_Emp=zeros(numbPoints); #initialize variables
indexTest=2:3; #choose a subset of [1 numbPoints]
global probTestEmp=0; #initialize variables
#loop through for each simulation
for ss=1:numbSim
    #run determinantal simuation
    indexDPP=funSimSimpleLDPP(eigenVectL,eigenValL); #returns index
    global probX_i_Emp[indexDPP]=probX_i_Emp[indexDPP].+1;

    countTemp=0; #initialize count
    for ii=1:length(indexTest)
        #check that each point of test subset appears
        countTemp=countTemp+any(indexDPP.==indexTest[ii]);
    end
    global probTestEmp=probTestEmp.+(countTemp.==length(indexTest));
end

#empirically estimate the probabilities of each point appearing
probX_i_Emp=probX_i_Emp./numbSim
println("probX_i_Emp = ", probX_i_Emp);

#calculate exactly the probabilities of each point appearing
K=funLtoK(L);
probX_i_Exact=diag(K)
println("probX_i_Exact = ", probX_i_Exact);

#empirically estimate the probabilities of test subset appearing
probTestEmp=probTestEmp./numbSim
println("probTestEmp = ", probTestEmp);

#calculate exactly the probabilities of test subset appearing
probTestExact=det(K[indexTest,indexTest])
println("probTestExact = ", probTestExact);

###END Testing DPP simulation END###

ext: jl


text:
using PyPlot, Revise, transD_GP, DelimitedFiles, Random, BenchmarkTools
## set geometry
rMin = 100 #m
rMax = 20000  #m
nRx = 150
freqs   = [0.1, 0.3, 0.7] #Hz
# receiver geometry
rRx   = collect(LinRange(rMin,rMax,nRx))  # Ranges to receivers   (m)
zRx = 975.     # Depth of Receivers, if this is an array it calls normal FHT without lagged convolution
zTx = 1000.
RxAzim=40.
TxDip=20.
# Note that the receiver depth needs to be in same model layer as transmitter.
# model
zfixed   = [-1e6,    0,      1000,   2000, 2100]
rhofixed = [1e13,    0.3,    1,      100,  1]
## set operator with lagged convolution
F = transD_GP.CSEM1DEr.RadialErLagged(zTx    = [zTx],
                      rRx    = rRx,
                      freqs  = freqs,
                      zRx    = [zRx],
                      RxAzim = RxAzim,
                      TxDip  = TxDip)

## model one wavenumber
krho = 2.0593e-10
rxno = 1
txno = 1
Ek = transD_GP.CSEM1DEr.getCSEM1DKernelsEr!(F, krho, freqs[1], zfixed, rhofixed, rxno, txno)
## timing
@info "Lagged convolution time"
transD_GP.CSEM1DEr.getfield!(F, zfixed, rhofixed)
@btime transD_GP.CSEM1DEr.getfield!($F, $zfixed, $rhofixed)
# Discrete Hankel Tx
G = transD_GP.CSEM1DEr.RadialErDHT(rRx  = F.rRx,
                          freqs  = F.freqs,
                          zRx    = F.zRx[1]*ones(length(F.rRx)),
                          RxAzim = RxAzim*ones(length(F.rRx)),
                          zTx    = F.zTx[1]*ones(length(F.rRx)),
                          TxDip  = TxDip*ones(length(F.rRx))
                          )
@info "Vanilla DHT time"
transD_GP.CSEM1DEr.getfield!(G, zfixed, rhofixed)
@btime transD_GP.CSEM1DEr.getfield!($G, $z, $ρ)
## plot
figure()
subplot(121)
step(log10.(rhofixed[2:end]), zfixed[2:end])
xlim(-1, 2)
ylim(1000, 3700)
grid()
gca().invert_yaxis()
subplot(122)
semilogy(rRx, abs.(F.Er), label="julia lagged HT")
semilogy(rRx, abs.(G.Er), label="julia vanilla DHT")
include("csem_Er_response.jl")
semilogy(rRx, abs.(Erm), "--", label="matlab 2012 code")
grid()
legend()
## random, many layers
Random.seed!(435)
nlayers = 100
z = [-1e6, 0, 1000, 1000 .+ cumsum(50*rand(nlayers-3))...]
ρ = [1e13, 0.3, 1, 10 .^(1 .+ 1.5*rand(nlayers-3))...]
figure()
subplot(121)
step(log10.(ρ[2:end]), z[2:end])
xlim(1,2.5 )
grid(); gca().invert_yaxis()
@info "timing for $nlayers layers"
transD_GP.CSEM1DEr.getfield!(F, z, ρ)
@btime transD_GP.CSEM1DEr.getfield!($F, $z, $ρ)
subplot(122)
semilogy(rRx, abs.(F.Er), label="julia")


text:
# ---
# title: 899. Orderly Queue
# id: problem899
# author: Indigo
# date: 2021-02-01
# difficulty: Hard
# categories: Math, String
# link: <https://leetcode.com/problems/orderly-queue/description/>
# hidden: true
# ---
#
# A string `S` of lowercase letters is given.  Then, we may make any number of
# _moves_.
#
# In each move, we choose one of the first `K` letters (starting from the left),
# remove it, and place it at the end of the string.
#
# Return the lexicographically smallest string we could have after any number of
# moves.
#
#
#
# **Example 1:**
#
#
#
#     Input: S = "cba", K = 1
#     Output: "acb"
#     Explanation:
#     In the first move, we move the 1st character ("c") to the end, obtaining the string "bac".
#     In the second move, we move the 1st character ("b") to the end, obtaining the final result "acb".
#
#
# **Example 2:**
#
#
#
#     Input: S = "baaca", K = 3
#     Output: "aaabc"
#     Explanation:
#     In the first move, we move the 1st character ("b") to the end, obtaining the string "aacab".
#     In the second move, we move the 3rd character ("c") to the end, obtaining the final result "aaabc".
#
#
#
#
# **Note:**
#
#   1. `1 <= K <= S.length <= 1000`
#   2. `S` consists of lowercase letters only.
#
#
## @lc code=start
using LeetCode

function orderly_queue(S::String, K::Int)
    if K == 1
        res = S
        for i in 2:length(S)
            res = min(res, S[i:end] * S[1:i-1])
        end
        return res
    end
    c_units = codeunits(S)[:]
    return c_units |> sort! |> String
end
## @lc code=end


text:
module Exchanges
    # Подключаем библиотеку, связанную с временем
    include("../../utils/Utils.jl")
    import .Utils

    using AMQPClient

    export declare, delete

    """Функция удаления обменника.

    channel: канал связи с сервисом сообщений,
    exchange_id: идентификатор обменника
    """
    function delete(channel::AMQPClient.MessageChannel, exchange_id::String)::Bool
        # Выводим информационное сообщение
        @info Utils.Chronometer.message_with_time("Запущена функция удаления обменника...")
        # Логе для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ" channel
        @debug "> идентификатор обменника" exchange_id
        # Пытаемся удалить обменник
        try
            # Вызываем операцию удаления обменника по идентификатору
            result::Bool = AMQPClient.exchange_delete(channel, exchange_id)
            # Дополняем логирование для режима отладки
            @debug ">> результат операции удаления обменника" result
            # В зависимости от результат присваиваем заначение сообщени
            # для информационного уровня логирования и выводим его
            message::String = result ? "Обменник успешно удалён" : "Не удалось удалить обменник"
            @info string("=>", message)
            # Возвращаем результат выполения операции
            return result
        catch error
            #
            @error Utils.Chronometer.message_with_time("Удаление обменника завершилось ошибкой ") error
            throw(error)
        end
    end

    """Функция регистрации обменника.

    channel: сканал связи с сервисом сообщений,
    exc_id: идентификатор обменника,
    exc_type: тип обменника
    """
    function declare(channel::AMQPClient.MessageChannel, exc_id::String, exc_type::String)::Dict{String, Any}
        # Выводим информационное сообщение
        @info Utils.Chronometer.message_with_time("Запущена функция регистрации обменника...")
        # Логерд для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с сервером RabbitMQ" channel
        @debug "> идентификатор обменника" exc_id
        @debug "> тип обменника" exc_type
        # Инициализируем переменную типа обменника с пустым значением
        exchange_type = nothing
        # Задаём тип обменника в зависимости от переданного наименования типа
        if exc_type == "direct"
            exchange_type = AMQPClient.EXCHANGE_TYPE_DIRECT
        elseif exc_type == "fanout"
            exchange_type = AMQPClient.EXCHANGE_TYPE_FANOUT
        elseif exc_type == "topic"
            exchange_type = AMQPClient.EXCHANGE_TYPE_TOPIC
        elseif exc_type == "headers"
            exchange_type = AMQPClient.EXCHANGE_TYPE_HEADERS
        else
            exchange_type = AMQPClient.EXCHANGE_TYPE_DIRECT
        end
        # Дополняем логирвание для режима отладки
        @debug ">> вычисленный тип обменника" exchange_type
        # Пытаемся объявить обменник
        try
            # Вызываем операцию объявления (регистрации) обменника на сервер RabbitMQ
            result::Bool = AMQPClient.exchange_declare(channel, exc_id, exchange_type)
            # Дополняем логирование для режима отладки
            @debug ">> результат операции регистрации обменника" result
            # В зависимости от результата операции задаём значение сообщения для
            # информационног уровня логирования и выводим его
            message::String = result ? "Регистрация обменника прошла успешно" : "Не удалось ззарегистрировать обменник"
            @info string("=> ", message)
            # В качестве результата возвращаем словарь, содержащий ключи:
            # NAME - наименование обменника
            # TYPE - тип обменника
            return Dict{String, Any}(
                "NAME" => exc_id,
                "TYPE" => exchange_type
            )
        catch error
            # Обрабатываем ошибки: выводим в терминал и кидаем исключение
            @error Utils.Chronometer.message_with_time("Удаление обменника завершилось ошибкой") error
            throw(error)
        end
    end
end


text:
@testset "Launch Configuration" begin
    kernel() = nothing

    device = AMDGPU.default_device()
    queue = AMDGPU.default_queue(device)

    # Group/grid size selection and aliases
    for (groupsize,gridsize) in (
        (1,1),
        (2,4),
        (1024,1024),

        ((1,1),(2,2)),
        ((1024,1),(1024,2)),

        ((1,1,1),(2,2,2)),
        ((1024,1,1),(1024,2,2)),

        ((1,1,1),2),
        (1,(1024,1,1)),
    )
        eval(:(@roc groupsize=$groupsize $kernel()))
        eval(:(@roc groupsize=$groupsize gridsize=$gridsize $kernel()))
        eval(:(@roc gridsize=$gridsize $kernel()))

        threads = groupsize
        blocks = gridsize .÷ groupsize
        eval(:(@roc threads=$threads $kernel()))
        eval(:(@roc blocks=$blocks $kernel()))
        eval(:(@roc threads=$threads blocks=$blocks $kernel()))
    end

    # Device/queue selection and aliases
    # FIXME: Test that device/queue are used!
    eval(:(@roc device=$device $kernel()))
    eval(:(@roc device=$device queue=$queue $kernel()))
    eval(:(@roc queue=$queue $kernel()))
    eval(:(@roc agent=$device $kernel()))
    eval(:(@roc stream=$queue $kernel()))
    eval(:(@roc agent=$device queue=$queue $kernel()))

    # Group size validity
    @test_throws ArgumentError eval(:(@roc groupsize=0 $kernel()))
    eval(:(@roc groupsize=1024 $kernel()))
    @test_throws ArgumentError eval(:(@roc groupsize=1025 $kernel()))

    # No-launch
    @test eval(:(@roc launch=true $kernel())) isa AMDGPU.RuntimeEvent
    @test eval(:(@roc launch=false $kernel())) isa AMDGPU.HostKernel
    @test_throws Exception eval(:(@roc launch=1 $kernel())) # TODO: ArgumentError
end

@testset "No-argument kernel" begin
    kernel() = nothing

    wait(@roc kernel())
end

@testset "Kernel argument alignment" begin
    function kernel(x, y)
        if Int64(x) != y
            error("Fail!")
        end
        nothing
    end
    x = rand(UInt32)
    y = Int64(x)
    wait(@roc kernel(x, y))
end

@testset "Signal waiting" begin
    sig = @roc identity(nothing)
    wait(sig)
    wait(sig.event)
    wait(sig.event.signal)
    wait(sig.event.signal.signal[])
end

@testset "Custom signal" begin
    sig = HSASignal()
    sig2 = @roc signal=sig identity(nothing)
    @test sig2.event.signal == sig
    wait(sig)
    wait(sig2)
end

if length(AMDGPU.get_agents(:gpu)) > 1
    @testset "Multi-GPU" begin
        # HSA will throw if the compiler and launch use different agents
        a1, a2 = AMDGPU.get_agents(:gpu)
        wait(@roc device=AMDGPU.RuntimeDevice(a1) identity(nothing))
        wait(@roc device=AMDGPU.RuntimeDevice(a2) identity(nothing))
    end
else
    @warn "Only 1 GPU detected; skipping multi-GPU tests"
    @test_broken "Multi-GPU"
end


text:
import BinDeps
using DataDeps

function with_accept(f, manual_overwrite)
    auto_accept = if manual_overwrite == nothing
        get(ENV, "DATADEPS_ALWAYS_ACCEPT", false)
    else
        manual_overwrite
    end
    withenv(f, "DATADEPS_ALWAYS_ACCEPT" => string(auto_accept))
end

function datadir(depname, dir = nothing; i_accept_the_terms_of_use = nothing)
    with_accept(i_accept_the_terms_of_use) do
        if dir == nothing
            # use DataDeps defaults
            @datadep_str depname
        else
            # use user-provided dir
            if isdir(dir)
                dir
            else
                DataDeps.env_bool("DATADEPS_DISABLE_DOWNLOAD") && error("DATADEPS_DISABLE_DOWNLOAD enviroment variable set. Can not trigger download.")
                DataDeps.download(DataDeps.registry[depname], dir)
                dir
            end
        end
    end::String
end

function datafile(depname, filename, dir = nothing; recurse = true, kw...)
    path = joinpath(datadir(depname, dir; kw...), filename)
    if !isfile(path)
        @warn "The file \"$path\" does not exist, even though the dataset-specific folder does. This is an unusual situation that may have been caused by a manual creation of an empty folder, or manual deletion of the given file \"$filename\"."
        if dir == nothing
            @info "Retriggering DataDeps.jl for \"$depname\""
            download_dep(depname; kw...)
        else
            @info "Retriggering DataDeps.jl for \"$depname\" to \"$dir\"."
            download_dep(depname, dir; kw...)
        end
        if recurse
            datafile(depname, filename, dir; recurse = false, kw...)
        else
            error("The file \"$path\" still does not exist. One possible explaination could be a spelling error in the name of the requested file.")
        end
    else
        path
    end::String
end

function download_dep(depname, dir = DataDeps.determine_save_path(depname); kw...)
    DataDeps.download(DataDeps.registry[depname], dir; kw...)
end

function download_docstring(modname, depname)
    """
    The corresponding resource file(s) of the dataset is/are
    expected to be located in the specified directory `dir`. If
    `dir` is omitted the directories in
    `DataDeps.default_loadpath` will be searched for an existing
    `$(depname)` subfolder. In case no such subfolder is found,
    `dir` will default to `~/.julia/datadeps/$(depname)`. In the
    case that `dir` does not yet exist, a download prompt will be
    triggered. You can also use `$(modname).download([dir])`
    explicitly for pre-downloading (or re-downloading) the
    dataset. Please take a look at the documentation of the
    package DataDeps.jl for more detail and configuration
    options.
    """
end


text:
"""
List of variables derived in this notebook:

- climatologies
- global means
- zonal means
- geographic maps
- transect transports
- MOC, MHT

Sample workflow:

```
## Setup Computation Parameters
@everywhere sol0="r2"
@everywhere nam="THETA"
@everywhere calc="clim"
@everywhere kk=1

## Preliminary Steps
@everywhere include("ECCO_pkg_grid_etc.jl")
@everywhere pth_in,pth_out,pth_tmp,sol,nt,list_steps=ECCO_path_etc(sol0,calc,nam)
!isdir(pth_out) ? mkdir(pth_out) : nothing
!isdir(pth_tmp) ? mkdir(pth_tmp) : nothing

## Main Computation
include("ECCO_standard_analysis.jl")
```
"""

## climatological mean

@everywhere function comp_clim(tmp_m,tmp_s1,tmp_s2,kk,m)
    nm=length(m:12:nt)
    tmp_m[:,:,m].=0.0
    tmp_s1[:,:,m].=0.0
    tmp_s2[:,:,m].=0.0
    for t in m:12:nt
        tmp=read_monthly(sol,nam,t,list_steps)
        ndims(tmp)>1 ? tmp=tmp[:,kk] : nothing
        tmp_m[:,:,m]=tmp_m[:,:,m]+1.0/nm*γ.write(tmp)
        tmp_s1[:,:,m]=tmp_s1[:,:,m]+γ.write(tmp)
        tmp_s2[:,:,m]=tmp_s2[:,:,m]+γ.write(tmp).^2
    end
end

function main_clim(sol,nam,kk=1)

    tmp_s1 = SharedArray{Float64}(γ.ioSize...,12)
    tmp_s2 = SharedArray{Float64}(γ.ioSize...,12)
    tmp_m = SharedArray{Float64}(γ.ioSize...,12)

    tmp=read_monthly(sol,nam,1,list_steps)
    ndims(tmp)>1 ? nz=size(tmp,2) : nz=1
    nz==1 ? kk=1 : nothing
    nz>1 ? suff=Printf.@sprintf("_k%02d",kk) : suff=""

    @sync @distributed for m in 1:12
        comp_clim(tmp_m,tmp_s1,tmp_s2,kk,m)
    end

    tmp0=read(tmp_m[:],γ)

    tmp=1.0/nt*sum(tmp_s1,dims=3)
    tmp1=read(tmp[:],Γ.XC)

    tmp=1/nt*sum(tmp_s2,dims=3)-tmp.^2
    tmp[findall(tmp.<0.0)].=0.0
    tmp=sqrt.(nt/(nt-1)*tmp)
    tmp2=read(tmp[:],Γ.XC)

    fil_out=joinpath(pth_tmp,nam*suff*".jld2")
    save(fil_out,"mean",tmp1,"std",tmp2,"mon",tmp0)

    return true
end

## global mean

@everywhere function comp_glo(glo,t)
    tmp=read_monthly(sol,nam,t,list_steps)
    if calc=="glo2d"
        tmp=[nansum(tmp[i,j].*Γ.RAC[i]) for j in 1:nr, i in eachindex(Γ.RAC)]
    else
        tmp=[nansum(tmp[i,j].*Γ.hFacC[i,j].*Γ.RAC[i]*Γ.DRF[j]) for j in 1:nr, i in eachindex(Γ.RAC)]
    end
    glo[:,t]=nansum(tmp,2)
end

function main_glo(sol,nam)

    glo = SharedArray{Float64}(nr,nt)
    @sync @distributed for t in 1:nt
        comp_glo(glo,t)
    end

    if calc=="glo2d"
        tmp=[glo[r,t]/Γ.tot_RAC[r] for t in 1:nt, r in 1:nr]
    else
        tmp=[nansum(glo[:,t])/nansum(Γ.tot_VOL) for t in 1:nt]
    end
    save_object(joinpath(pth_tmp,calc*".jld2"),tmp)
end

##

@everywhere function comp_msk0(msk0,zm0,l)
    lats=load(joinpath(pth_tmp,calc*"_lats.jld2"),"single_stored_object")
    dlat=lats[2]-lats[1]
    la0=lats[l]-dlat/2
    la1=lats[l]+dlat/2
    if la1<0.0
        msk=1.0*(Γ.YC.>=la0)*(Γ.YC.<la1)
    elseif la0>0.0
        msk=1.0*(Γ.YC.>la0)*(Γ.YC.<=la1)
    else
        msk=1.0*(Γ.YC.>=la0)*(Γ.YC.<=la1)
    end
    msk[findall(msk.==0.0)].=NaN;
    msk0[:,:,l]=write(msk*Γ.RAC)

    tmp2=[nansum(Γ.mskC[i,j].*msk[i].*Γ.RAC[i]) for j in 1:nr, i in eachindex(Γ.RAC)]
    zm0[l,:]=1.0 ./nansum(tmp2,2)
end

@everywhere function comp_zonmean(zm,t)
    lats=load(joinpath(pth_tmp,calc*"_lats.jld2"),"single_stored_object")
    nl=length(lats)
    tmp=read_monthly(sol,nam,t,list_steps)
    for l in 1:nl
        mskrac=read(msk0[:,:,l],γ)
        tmp1=[nansum(tmp[i,j].*mskrac[i]) for j in 1:nr, i in eachindex(Γ.RAC)]
        zm[l,:,t]=nansum(tmp1,2).*zm0[l,:]
    end
end

@everywhere function comp_zonmean2d(zm,t)
    lats=load(joinpath(pth_tmp,calc*"_lats.jld2"),"single_stored_object")
    nl=length(lats)
    tmp=read_monthly(sol,nam,t,list_steps)
    for l in 1:nl
        mskrac=read(msk0[:,:,l],γ)
        tmp1=[nansum(tmp[i].*mskrac[i]) for i in eachindex(Γ.RAC)]
        zm[l,t]=nansum(tmp1)*zm0[l,1]
    end
end

function main_zonmean(sol,nam)
    dlat=2.0
    lats=(-90+dlat/2:dlat:90-dlat/2)
    save_object(joinpath(pth_tmp,calc*"_lats.jld2"),lats)
    nl=length(lats)

    msk0 = SharedArray{Float64}(γ.ioSize...,nl)
    zm0 = SharedArray{Float64}(nl,nr)
    @sync @distributed for l in 1:nl
        comp_msk0(msk0,zm0,l)
    end
    save_object(joinpath(pth_tmp,calc*"_zm0.jld2"),zm0)
    save_object(joinpath(pth_tmp,calc*"_msk0.jld2"),msk0)

    #to speed up main loop, reuse:
    #- precomputed msk*RAC once and for all
    #- precomputed 1.0./nansum(tmp2,2)

    @everywhere msk0=load(joinpath(pth_tmp,calc*"_msk0.jld2"),"single_stored_object")
    @everywhere zm0=load(joinpath(pth_tmp,calc*"_zm0.jld2"),"single_stored_object")

    if (calc=="zonmean")
        zm = SharedArray{Float64}(nl,nr,nt)
        @sync @distributed for t in 1:nt
            comp_zonmean(zm,t)
        end
    else
        zm = SharedArray{Float64}(nl,nt)
        @sync @distributed for t in 1:nt
            comp_zonmean2d(zm,t)
        end
    end
    save_object(joinpath(pth_tmp,calc*".jld2"),zm)

    return true
end

##

@everywhere function comp_overturn(ov,t)
    U=read_monthly(sol,"UVELMASS",t,list_steps)
    V=read_monthly(sol,"VVELMASS",t,list_steps)
    (Utr,Vtr)=UVtoTransport(U,V,Γ)
    #integrate across latitude circles
    for z=1:nr
        UV=Dict("U"=>Utr[:,z],"V"=>Vtr[:,z],"dimensions"=>["x","y"])
        [ov[l,z,t]=ThroughFlow(UV,LC[l],Γ) for l=1:nl]
    end
    #integrate from bottom
    ov[:,:,t]=reverse(cumsum(reverse(ov[:,:,t],dims=2),dims=2),dims=2)
    #
    true
end

function main_overturn(sol,nam)
    ov = SharedArray{Float64}(nl,nr,nt)
    @sync @distributed for t in 1:nt
        comp_overturn(ov,t)
    end

    save_object(joinpath(pth_tmp,calc*".jld2"),ov)
	"Done with overturning"
end

##

@everywhere function comp_MHT(MHT,t)
    U=read_monthly(sol,"ADVx_TH",t,list_steps)
    V=read_monthly(sol,"ADVy_TH",t,list_steps)
    U=U+read_monthly(sol,"DFxE_TH",t,list_steps)
    V=V+read_monthly(sol,"DFyE_TH",t,list_steps)

    [U[i][findall(isnan.(U[i]))].=0.0 for i in eachindex(U)]
    [V[i][findall(isnan.(V[i]))].=0.0 for i in eachindex(V)]

    Tx=0.0*U[:,1]
    Ty=0.0*V[:,1]
    [Tx=Tx+U[:,z] for z=1:nr]
    [Ty=Ty+V[:,z] for z=1:nr]

    UV=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"])
    [MHT[l,t]=1e-15*4e6*ThroughFlow(UV,LC[l],Γ) for l=1:nl]
end

function main_MHT(sol,nam)
    MHT = SharedArray{Float64}(nl,nt)
    @sync @distributed for t in 1:nt
        comp_MHT(MHT,t)
    end
    save_object(joinpath(pth_tmp,calc*".jld2"),MHT)
	"Done with MHT"
end

##

@everywhere function comp_trsp(trsp,t)
    U=read_monthly(sol,"UVELMASS",t,list_steps)
    V=read_monthly(sol,"VVELMASS",t,list_steps)
    (Utr,Vtr)=UVtoTransport(U,V,Γ)
    #integrate across transport lines
    for z=1:nr
        UV=Dict("U"=>Utr[:,z],"V"=>Vtr[:,z],"dimensions"=>["x","y"])
        [trsp[itr,z,t]=ThroughFlow(UV,msk_trsp[itr],Γ) for itr=1:ntr]
    end
end

function main_trsp(sol,nam)
    trsp = SharedArray{Float64}(ntr,nr,nt)
    @sync @distributed for t in 1:nt
        comp_trsp(trsp,t)
    end

    trsp=[(nam=list_trsp[itr],val=trsp[itr,:,:]) for itr=1:ntr]
    save_object(joinpath(pth_tmp,calc*".jld2"),trsp)
	"Done with transports"
end

function main_function(calc,sol,nam,kk)
    if calc=="clim"
        main_clim(sol,nam,kk)
    elseif (calc=="glo2d")||(calc=="glo3d")
        main_glo(sol,nam)
    elseif (calc=="zonmean")||(calc=="zonmean2d")
        main_zonmean(sol,nam)
    elseif (calc=="overturn")
        main_overturn(sol,nam)
    elseif (calc=="MHT")
        main_MHT(sol,nam)
    elseif (calc=="trsp")
        main_trsp(sol,nam)
    else
        println("unknown calc")
    end
end


text:
"""
The `Romberg` module provides the functions `romberg` and `romberg!` that
perform a Romberg integration of discrete 1-dimensional data `y` sampled at
equally spaced points over `x`.

Romberg integration combines the trapezoid method of integration with Richardson
extrapolation. Significantly lower error estimates of definite integrals can be
obtained using Romberg integration over the trapezoid method alone, at only
slightly higher computational cost.
"""
module Romberg

using Trapz

export romberg, romberg!

maxsteps(N) = trunc(Int, log2(prevpow(2, N)))

"""
    romberg(x::AbstractRange, y::AbstractVector)

Integrate `y` sampled at `x` by Romberg integration applying the maximum possible
number of Richardson extrapolation steps.

Romberg integration requires equally spaced points. This is enforced by
requiring `x` to be an `AbstractRange`, rather than a dense vector.

The algorithm is most efficient if `length(x) == 2ⁿ + 1` for any positive integer
`n`. However, the function will still work for any `length(x)`.

# Examples

```jldoctest
julia> x = range(0, π, length=2^8+1);

julia> romberg(x, sin.(x))
1.9999999999999996
```
"""
function romberg(x::AbstractRange{Tx}, y::AbstractVector{Ty}) where {Tx,Ty}
    N = length(x)

    # Integral over nothing
    N <= 1 && return zero(promote_type(Tx,Ty))

    @boundscheck begin
        N == length(y) || throw(DimensionMismatch("length of `y` not equal to length of `x`"))

        # NOTE: by requiring x::AbstractRange, a fixed step size is guaranteed
    end

    pN = prevpow(2, N-1) + 1
    if pN == N
        max_steps = maxsteps(N)
        return integrate(x, y, max_steps)
    else
        return recursive_integration(x, y)
    end
end

"""
    romberg(x::AbstractRange, y::AbstractVector, max_steps::Integer)

Integrate `y` sampled at `x` by Romberg integration applying `max_steps` of
Richardson extrapolation.

`max_steps` is the number of Richardson extrapolation steps applied (also equal
to the number of trapezoid integrations - 1). In general, the larger it is, the
lower error in the integration. The largest `max_steps` can be is
`log2(prevpow(2, length(x)))`.

When specifying `max_steps`, `length(x) - 1` must be a power of 2.

# Examples

```jldoctest
julia> x = range(0, π, length=2^8+1);

julia> romberg(x, sin.(x), 8)
1.9999999999999996
```

```jldoctest
julia> romberg(x, sin.(x), 9)
ERROR: DomainError with 9:
`max_steps` cannot exceed `log2(prevpow(2, length(x)))` = 8
[...]
```
"""
function romberg(x::AbstractRange{Tx}, y::AbstractVector{Ty}, max_steps::Integer) where {Tx,Ty}
    N = length(x)

    # Integral over nothing
    N <= 1 && return zero(promote_type(Tx,Ty))

    @boundscheck begin
        max_steps <= maxsteps(N) || throw(DomainError(max_steps, "`max_steps` cannot exceed `log2(prevpow(2, length(x)))` = $(maxsteps(length(x)))"))
        ispow2(N-1) || throw(DomainError(length(x), "`length(x) - 1` must be a power of 2"))
        N == length(y) || throw(DimensionMismatch("length of `y` not equal to length of `x`"))

        # NOTE: by requiring x::AbstractRange, a fixed step size is guaranteed
    end

    return integrate(x, y, max_steps)
end

"""
    romberg!(R::AbstractMatrix, x::AbstractRange, y::AbstractVector)

Integrate `y` sampled at `x` by Romberg integration, filling in the square
matrix `R` in place.

The size of `R` determines the number of Richardson extrapolation steps. If `R`
has size 6×6, 5 extrapolation steps will be performed. In general, for a matrix
with dimensions `L×L`, `max_steps` corresponds to `L - 1`.

This function is useful to see how the Romberg integration progressed.

# Examples

```jldoctest
julia> x = range(0, π, length=2^8+1);

julia> R = zeros(4, 4);

julia> romberg!(R, x, sin.(x))
4×4 Array{Float64,2}:
 1.92367e-16  0.0      0.0      0.0
 1.5708       2.0944   0.0      0.0
 1.89612      2.00456  1.99857  0.0
 1.97423      2.00027  1.99998  2.00001
```
"""
function romberg!(R::AbstractMatrix, x::AbstractRange, y::AbstractVector)
    N = length(x)

    N <= 1 && return fill!(R, 0)

    @boundscheck begin
        # Assume `size(R, 1)` is L (`max_steps + 1`)
        size(R,1) <= maxsteps(N)+1 || throw(DomainError(size(R,1), "dimenensions of `R` cannot be greater than `log2(prevpow(2, length(x))) + 1`"))
        isequal(size(R)...) || throw(DimensionMismatch("`R` must be square"))
        N == length(y) || throw(DimensionMismatch("length of `y` not equal to length of `x`"))

        # NOTE: by requiring x::AbstractRange, a fixed step size is guaranteed
    end

    integrate!(R, x, y)
end

function recursive_integration(x, y)
    # Needed if `length(x) - 1` is not a power of 2

    N = length(x)

    N <= 1 && return zero(promote_type(x,y))

    pN = prevpow(2, N-1) + 1

    if pN == N
        # `ispow2(N - 1) == true`
        return integrate(x, y, maxsteps(N))
    elseif pN < N
        # `ispow2(N - 1) == false`
        return integrate(x[1:pN], view(y,1:pN), maxsteps(pN)) + recursive_integration(x[pN:N], view(y,pN:N))
    else
        return trapz(x, y)
    end
end

@inline function integrate(x::AbstractRange{Tx}, y::AbstractVector{Ty}, max_steps::Integer) where {Tx,Ty}
    # `max_steps` are extrapolation steps, size of `R` is `max_steps + 1`
    L = max_steps + 1

    # `Rp` is "Rprevious", `Rc` is "Rcurrent"
    Rp = Array{promote_type(Tx,Ty)}(undef, L)
    Rc = similar(Rp)

    trapezoid!(Rp, x, y)
    extrapolate!(Rp, Rc)

    # best estimate
    return Rc[end]
end

@inline function integrate!(R::AbstractMatrix, x::AbstractRange, y::AbstractVector)
    L = size(R, 1)

    trapezoid!(view(R,:,1), x, y)
    extrapolate!(R)

    return R
end

@inline function trapezoid!(R, x, y)
    N = length(x)
    @inbounds for i in eachindex(R)
        step_size = div(N-1, 2^(i-1))

        if step_size == 0
            # only triggered when N == 1?
            idxs = 1:N-1:N
        else
            idxs = 1:step_size:N
        end

        R[i] = trapz(x[idxs], view(y,idxs))
    end
end

@inline function extrapolate!(Rp, Rc)
    L = length(Rp)
    @inbounds for j = 2:L
        # Precompute reused values
        pow_prev = 4^(j-1)
        den = inv(pow_prev - 1)

        @inbounds for i = j:L
            Rc[i] = (pow_prev*Rp[i] - Rp[i-1]) * den
        end

        # Only copying the relevant values! (not necessary on last loop)
        if j < L
            copyto!(Rp, j, Rc, j, L-j+1)
        end
    end
end

@inline function extrapolate!(R)
    L = size(R,1)
    @inbounds for j = 2:L
        # Precompute reused values
        pow_prev = 4^(j-1)
        den = inv(pow_prev - 1)

        @inbounds for i = j:L
            R[i,j] = (pow_prev*R[i,j-1] - R[i-1,j-1]) * den
        end
    end
end

function error_test(R, max_steps)
    # See http://www.sci.utah.edu/~beiwang/teaching/cs6210-fall-2016/lecture24.pdf

    # The ratio of the difference between successive entries in column `j`
    # should be approximately 4^j.

    ratio = zeros(max_steps-2, max_steps-2)
    for j = 1:max_steps-2
        ratio[j:end,j] .= (R[1+j:end-1,j] - R[j:end-2,j]) ./ (R[2+j:end,j] - R[1+j:end-1,j])
    end

    return ratio
end

end # module


text:
__precompile__()

using PyCall
using ImageView, GtkReactive, Colors, Images

using BackgroundSegmenter

@pyimport skvideo as skv
@pyimport numpy as np
@pyimport skvideo.io as skvio

function load(infile, limit=0)
    V = np.load(infile)
    (t, n, m) = size(V)
    if limit == 0
        return V
    end
    left = convert(Int, round((n - limit) / 2))
    top = convert(Int, round((m - limit) / 2))
    return V[:, left:(left+limit), top:(top+limit)]
end

function time_results(V, fgbg, M, t)
    @time for i in 2:t
        apply!(M, (view(V, :, :, i)), (view(fgbg, :, :, i)));
    end
end

# V = permutedims(load("data/lamprey2.npy"), [2, 3, 1]);
# V = V[300:450, 150:400, :];
V = permutedims(skvio.vread("data/lamprey4.avi")[:, :, :, 3], [2, 3, 1]);
f1 = BackgroundSegmenter.Filter(1, 0.2);
f2 = BackgroundSegmenter.Filter(3, 0.08);
f3 = BackgroundSegmenter.Filter(5, 0.04);
fgbg = zeros(V);
(n, m, t) = size(V)
M = MixtureModel(n, m);
apply!(M, (view(V, :, :, 1)), (view(fgbg, :, :, 1)));
time_results(V, fgbg, M, t);

f = view(fgbg, 1, :, :);
f = filter_components(f, 12);
@time for i in 2:t
    f = view(fgbg, i, :, :);
    f = filter_components(f, 12);
end
@time fgbg = filter_components(fgbg, 40);
# W = @view V[:, :, 2:end];
# colors = distinguishable_colors(blobs+1);
# components = Array{RGB{N0f8}}(size(fgbg));
# for idx in eachindex(fgbg)
#     components[idx] = colors[fgbg[idx]+1] * (W[idx] * 1.0/255.0);
# end
# imshow(components)
# @time fgbg1 = apply(f1, fgbg);
# A0 = V;
# A1 = scale(V, fgbg1);
# guidict = imshow(hcat(A0, A1))
# @time fgbg2 = apply(f2, fgbg);
# @time fgbg3 = apply(f3, fgbg);
# A2 = scale(V, fgbg2);
# A3 = scale(V, fgbg3);

# C = zeros(UInt8, 2n, 2m, t);
# C[1:n, 1:m, :] = A0;
# C[n+1:end, 1:m, :] = A2;
# C[1:n, m+1:end, :] = A1;
# C[n+1:end, m+1:end, :] = A3;

# idx = annotate!(guidict, AnnotationBox((20, 30), (40, 20), linewidth=2, color=RGB(0, 1, 0)))
# mm =  [MixtureModel(5) for i in 1:n, j in 1:m];
# cut = MinCut(n, m);
# for i in 1:5
#     time_results(V, fgbg, cut,  mm, i)
# end
# @time  for i in 6:t-1
#     fgbg[i, :, :] = label_components(apply_mrf!(cut, mm, (@view V[i, :, :]), 2.0), 12)
# end
# time_results(V, fgbg, cut,  mm, t)

# A = permutedims(V .* fgbg, [2, 3, 1])


text:
using Thebes

function convert_off_file(fname)
    infile = readlines(open(fname))
    filter!(l -> !isempty(l), infile)
    @show infile[1]
    if ! startswith(infile[1], "OFF")
        error("not an OFF file")
    end
    nverts, nfaces, nedges = map(Meta.parse, split(infile[2]))
    vertices = Point3D[]
    for i in 3:nverts+2
        a, b, c = map(Meta.parse, split(strip(infile[i])))
        push!(vertices, Point3D(a, -b, c))
    end

    faces = Array{Int64, 1}[]
    for i in nverts + 2 + 1: nverts + 2 + nfaces
        vals = map(Meta.parse, split(strip(infile[i]), r" |\t"))
        ta = Int64[]
        for j in 2:vals[1] + 1
            push!(ta, vals[j] + 1)
        end
        push!(faces, ta)
    end
    return (vertices, faces)
end

# cd(expanduser("/tmp/geometry"))
# iobuffer = IOBuffer()
#
# for f in ["boxcube.off", "boxtorus.off", "concave.off", "cone.off", "cross.off", "cube.off",
# "cuboctahedron.off", "dodecahedron.off", "geodesic.off", "helix2.off", "icosahedron.off",
# "icosidodecahedron.off", "octahedron.off", "octtorus.off", "rhombicosidodecahedron.off",
# "rhombicuboctahedron.off", "rhombitruncated_cubeoctahedron.off",
# "rhombitruncated_icosidodecahedron.off", "snub_cube.off", "snub_dodecahedron.off",
# "sphere2.off", "tet3d.off", "tetrahedron.off", "triangle.off", "truncated_cube.off",
# "truncated_dodecahedron.off", "truncated_icosahedron.off", "truncated_octahedron.off",
# "truncated_tetrahedron.off"]
#     fname, _ = splitext(f)
#     object = convert_off_file(f)
#     println(iobuffer, fname, " ", object)
# end
#
#=open("/tmp/objects.txt", "w") do f
    write(f, String(take!(iobuffer)))
end
for f in readdir()
object = convert_off_file(f)
println(iobuffer, f, " ", object)
end
=#

sphere1  = convert_off_file("/tmp/Untitled.obj.off")


text:
using LinearMaps
using ProgressMeter
using Printf

using Pkg
const cuda = if "CUDArt" ∈ keys(Pkg.installed())
    using CUDArt
    using CUBLAS
    using CUSPARSE
    true
else
    false
end

function SI(v, base = 1000)
    if v > 1 && v < 1e24
        prefixes = ["k", "M", "G", "T", "P", "E"]
        p = min(max(floor(Int, log(v)/log(base)), 0), length(prefixes))
        v /= base^p
        @sprintf("%.2f %s", v, p != 0 ? prefixes[p] : "")
    else
        @sprintf("%.3g ", v)
    end
end

function integrate(observe::Function,
                   v₀::AbstractVector,
                   tmax::Number, steps::Integer,
                   propagator::MagnusPropagator;
                   save_intermediate::Bool = false,
                   verbose::Bool = false)
    τ = tmax/steps
    V = similar(v₀, eltype(v₀), length(v₀), save_intermediate ? steps + 1 : 1)
    copy!(view(V, :, 1), v₀)
    cur = 1
    next = save_intermediate ? 2 : 1
    verbose && println("$(typeof(propagator))")
    prog = Progress(steps, 0.1, "Integrating ")
    tic()
    for i = 1:steps
        propagator((i-1)*τ, τ, view(V, :, cur), view(V, :, next))

        cur = next
        observe(view(V, :, next), i, τ)
        save_intermediate && (next += 1)
        verbose && ProgressMeter.update!(prog, i)
    end
    ms = toq()*1000
    verbose && println("Grid points/ms: ", SI(length(v₀)*steps/ms))
    Dict(:V => V,
         :milliseconds => ms,
         :performance => length(v₀)*steps/ms)
end

integrate(v₀::AbstractVector,
          tmax::Number, steps::Integer,
          propagator::MagnusPropagator;
          save_intermediate::Bool = false,
          verbose::Bool = false) =
              integrate((V,i,τ) -> (),
                        v₀,
                        tmax, steps, propagator;
                        save_intermediate = save_intermediate,
                        verbose = verbose)

function integrate(observe::Function,
                   v₀::AbstractVector,
                   tmax::Number, steps::Integer,
                   B::AbstractMatrix,
                   f::Function,
                   C::AbstractMatrix,
                   a::Number = 1;
                   propagator = :cfet,
                   exponentiator = :lanczos,
                   lanczos_m = 30,
                   mode = :cpu,
                   save_intermediate::Bool = false,
                   verbose::Bool = false)
    if mode == :gpu
        !cuda && error("Cuda not available, gpu integration impossible")
        verbose && println("Active device: $(CUDArt.name(CUDArt.device_properties(device())))")
        v₀ = CudaArray(v₀)
        B = upload(B)
        C = upload(C)
    end

    B = LinearMap(B)
    C = LinearMap(C)

    Exp = if exponentiator == :lanczos
        LanczosExponentiator(lanczos_m, v₀)
    elseif exponentiator == :scalar
        ScalarExponentiator(v₀)
    else
        error("Unknown exponentiator, $(string(exponentiator))")
    end

    propagator = if propagator == :cfet
        CFET4BfCPropagator(B, f, C, a, Exp)
    elseif propagator == :midpoint
        MidpointPropagator(t -> B + f(t)*C, a, Exp)
    else
        error("Unknown propagator, $(string(propagator))")
    end

    results = integrate(observe,
                        v₀,
                        tmax, steps, propagator;
                        save_intermediate = save_intermediate,
                        verbose = verbose)
    if mode == :gpu
        Dict(results..., :V => to_host(results[:V]))
    else
        results
    end
end

integrate(v₀::AbstractVector,
          tmax::Number, steps::Integer,
          B::AbstractMatrix,
          f::Function,
          C::AbstractMatrix,
          a::Number = 1;
          propagator = :cfet,
          exponentiator = :lanczos,
          lanczos_m = 30,
          mode = :cpu,
          save_intermediate::Bool = false,
          verbose::Bool = false) =
              integrate((V,i,τ) -> (),
                        v₀,
                        tmax, steps,
                        B, f, C, a;
                        propagator = propagator,
                        exponentiator = exponentiator,
                        lanczos_m = lanczos_m,
                        mode = mode,
                        save_intermediate = save_intermediate,
                        verbose = verbose)

export integrate


text:
# using GeometryPrimitives: orthoaxes
using FFTW
export MaxwellGrid, MaxwellData, t2c, c2t, kx_c2t, kx_t2c, zx_t2c, kxinv_c2t, kxinv_t2c, ε⁻¹_dot, ε_dot_approx, M, M̂, P, P̂, Mₖ, M̂ₖ, t2c!, c2t!, kcross_c2t!, kcross_t2c!, zcross_t2c!, kcrossinv_c2t!, kcrossinv_t2c!, ε⁻¹_dot!, ε_dot_approx!, M!, M̂!, P!, P̂!, calc_kpg, H_Mₖ_H

struct MaxwellGrid
    Δx::Float64
    Δy::Float64
    Δz::Float64
    Nx::Int64
    Ny::Int64
    Nz::Int64
    δx::Float64
    δy::Float64
    δz::Float64
    x::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    y::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    z::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    g⃗::Array{Array{Float64,1},3}
	𝓕::FFTW.cFFTWPlan
	𝓕⁻¹::AbstractFFTs.ScaledPlan
	𝓕!::FFTW.cFFTWPlan
	𝓕⁻¹!::AbstractFFTs.ScaledPlan
end

MaxwellGrid(Δx::Float64,Δy::Float64,Δz::Float64,Nx::Int,Ny::Int,Nz::Int) = MaxwellGrid(
    Δx,
    Δy,
    Δz,
    Nx,
    Ny,
    Nz,
    Δx / Nx,    # δx
    Δy / Ny,    # δy
    Δz / Nz,    # δz
    ( ( Δx / Nx ) .* (0:(Nx-1))) .- Δx/2.,  # x
    ( ( Δy / Ny ) .* (0:(Ny-1))) .- Δy/2.,  # y
    ( ( Δz / Nz ) .* (0:(Nz-1))) .- Δz/2.,  # z
    [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)], # g⃗
    # (𝓕 = plan_fft(randn(ComplexF64, (3,Nx,Ny,Nz))); inv(𝓕); 𝓕),  # planned FFT operator 𝓕
    # (𝓕! = plan_fft!(randn(ComplexF64, (3,Nx,Ny,Nz))); inv(𝓕!); 𝓕!), # planned in-place FFT operator 𝓕!
	plan_fft(randn(ComplexF64, (3,Nx,Ny,Nz)),(2:4)),  # planned FFT operator 𝓕
	plan_ifft(randn(ComplexF64, (3,Nx,Ny,Nz)),(2:4)),
	plan_fft!(randn(ComplexF64, (3,Nx,Ny,Nz)),(2:4)),
	plan_ifft!(randn(ComplexF64, (3,Nx,Ny,Nz)),(2:4)), # planned in-place FFT operator 𝓕!
)

MaxwellGrid(Δx::Float64,Δy::Float64,Nx::Int,Ny::Int) = MaxwellGrid(
    Δx,
    Δy,
    1.,
    Nx,
    Ny,
    1,
    Δx / Nx,    # δx
    Δy / Ny,    # δy
    1.,    # δz
    ( ( Δx / Nx ) .* (0:(Nx-1))) .- Δx/2.,  # x
    ( ( Δy / Ny ) .* (0:(Ny-1))) .- Δy/2.,  # y
    0.0:1.0:0.0,  # z
    [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(1,1.0)], # g⃗
    # (𝓕 = plan_fft(randn(ComplexF64, (3,Nx,Ny,1))); inv(𝓕); 𝓕),  # planned FFT operator 𝓕
    # (𝓕! = plan_fft!(randn(ComplexF64, (3,Nx,Ny,1))); inv(𝓕!); 𝓕!), # planned in-place FFT operator 𝓕!
	plan_fft(randn(ComplexF64, (3,Nx,Ny,1)),(2:4)),  # planned FFT operator 𝓕
	plan_ifft(randn(ComplexF64, (3,Nx,Ny,1)),(2:4)),
	plan_fft!(randn(ComplexF64, (3,Nx,Ny,1)),(2:4)),
	plan_ifft!(randn(ComplexF64, (3,Nx,Ny,1)),(2:4)), # planned in-place FFT operator 𝓕!
)

mutable struct MaxwellData
    k::Float64
    ω²::Float64
    ω²ₖ::Float64
    ω::Float64
    ωₖ::Float64
    H⃗::Array{ComplexF64,2}
    H::Array{ComplexF64,4}
    e::Array{ComplexF64,4}
    d::Array{ComplexF64,4}
    grid::MaxwellGrid
	Δx::Float64
    Δy::Float64
    Δz::Float64
    Nx::Int64
    Ny::Int64
    Nz::Int64
	Neigs::Int64
    δx::Float64
    δy::Float64
    δz::Float64
    x::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    y::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    z::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    g⃗::Array{Array{Float64,1},3}
    mn::Array{Float64,5}
	kpg_mag::Array{Float64,3}
    𝓕::FFTW.cFFTWPlan
	𝓕⁻¹::AbstractFFTs.ScaledPlan
    𝓕!::FFTW.cFFTWPlan
	𝓕⁻¹!::AbstractFFTs.ScaledPlan
end

MaxwellData(k::Float64,g::MaxwellGrid,Neigs::Int64) = MaxwellData(
    k,
    0.0,
    0.0,
    0.0,
    0.0,
    randn(ComplexF64,(2*g.Nx*g.Ny*g.Nz,Neigs)),
    randn(ComplexF64,(2,g.Nx,g.Ny,g.Nz)),
    randn(ComplexF64,(3,g.Nx,g.Ny,g.Nz)),
    randn(ComplexF64,(3,g.Nx,g.Ny,g.Nz)),
    g,
    g.Δx,
    g.Δy,
    g.Δz,
    g.Nx,
    g.Ny,
    g.Nz,
	Neigs,
    g.δx,       # δx
    g.δy,       # δy
    g.δz,       # δz
    g.x,        # x
    g.y,        # y
    g.z,        # z
    g.g⃗,
    calc_kpg(k,g.g⃗)[2], # ( (kpg_mag, kpg_mn) = calc_kpg(k,g.g⃗); kpg_mn), #( (kpg_mag, kpg_mn) = calc_kpg(k,g.Δx,g.Δy,g.Δz,g.Nx,g.Ny,g.Nz); kpg_mn),  # mn
	calc_kpg(k,g.g⃗)[1], # kpg_mag,
    g.𝓕,
	g.𝓕⁻¹,
    g.𝓕!,
	g.𝓕⁻¹!,
)

MaxwellData(k::Float64,g::MaxwellGrid) = MaxwellData(k,g,1)
MaxwellData(k::Float64,Δx::Float64,Δy::Float64,Δz::Float64,Nx::Int,Ny::Int,Nz::Int) = MaxwellData(k,MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz))
MaxwellData(k::Float64,Δx::Float64,Δy::Float64,Nx::Int,Ny::Int) = MaxwellData(k,MaxwellGrid(Δx,Δy,Nx,Ny))

# non-Mutating Operators

function calc_kpg(kz::T,g⃗::Array{Array{T,1},3})::Tuple{Array{T,3},Array{T,5}} where T <: Real
	g⃗ₜ_zero_mask = Zygote.@ignore( [ sum(abs2.(gg[1:2])) for gg in g⃗ ] .> 0. );
	g⃗ₜ_zero_mask! = Zygote.@ignore( .!(g⃗ₜ_zero_mask) );

	ŷ = [0.; 1. ;0.]
	k⃗ = [0.;0.;kz]
	# @tullio kpg[a,i,j,k] := k⃗[a] - g⃗[i,j,k][a] nograd=g⃗ fastmath=false
	@tullio kpg[a,i,j,k] := k⃗[a] - g⃗[i,j,k][a] fastmath=false
	@tullio kpg_mag[i,j,k] := sqrt <| kpg[a,i,j,k]^2 fastmath=false
	zxinds = [2; 1; 3]
	zxscales = [-1; 1. ;0.] #[[0. -1. 0.]; [-1. 0. 0.]; [0. 0. 0.]]
	@tullio kpg_nt[a,i,j,k] := zxscales[a] * kpg[zxinds[a],i,j,k] * g⃗ₜ_zero_mask[i,j,k] + ŷ[a] * g⃗ₜ_zero_mask![i,j,k]  nograd=(zxscales,zxinds,ŷ,g⃗ₜ_zero_mask,g⃗ₜ_zero_mask!) fastmath=false
	@tullio kpg_nmag[i,j,k] := sqrt <| kpg_nt[a,i,j,k]^2 fastmath=false
	@tullio kpg_n[a,i,j,k] := kpg_nt[a,i,j,k] / kpg_nmag[i,j,k] fastmath=false
	xinds1 = [2; 3; 1]
	xinds2 = [3; 1; 2]
	@tullio kpg_mt[a,i,j,k] := kpg_n[xinds1[a],i,j,k] * kpg[xinds2[a],i,j,k] - kpg[xinds1[a],i,j,k] * kpg_n[xinds2[a],i,j,k] nograd=(xinds1,xinds2) fastmath=false
	@tullio kpg_mmag[i,j,k] := sqrt <| kpg_mt[a,i,j,k]^2 fastmath=false
	@tullio kpg_m[a,i,j,k] := kpg_mt[a,i,j,k] / kpg_mmag[i,j,k] fastmath=false
	kpg_mn_basis = [[1. 0.] ; [0. 1.]]
	@tullio kpg_mn[a,b,i,j,k] := kpg_mn_basis[b,1] * kpg_m[a,i,j,k] + kpg_mn_basis[b,2] * kpg_n[a,i,j,k] nograd=kpg_mn_basis fastmath=false
	return kpg_mag, kpg_mn
end

function calc_kpg(kz::T,Δx::T,Δy::T,Δz::T,Nx::Int64,Ny::Int64,Nz::Int64)::Tuple{Array{T,3},Array{T,5}} where T <: Real
	g⃗ = Zygote.@ignore( [ [gx;gy;gz] for gx in collect(fftfreq(Nx,Nx/Δx)), gy in collect(fftfreq(Ny,Ny/Δy)), gz in collect(fftfreq(Nz,Nz/Δz))] )
	# g⃗ = [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)]
	calc_kpg(kz,Zygote.dropgrad(g⃗))
end

# Non-Mutating Operators

"""
    t2c: v⃗ (transverse vector) → a⃗ (cartesian vector)
"""
function t2c(H,mn)
    @tullio h[a,i,j,k] := H[b,i,j,k] * mn[a,b,i,j,k] fastmath=false
end

"""
    c2t: a⃗ (cartesian vector) → v⃗ (transverse vector)
"""
function c2t(h,mn)
    @tullio H[a,i,j,k] := h[b,i,j,k] * mn[b,a,i,j,k] fastmath=false
end

"""
    kx_t2c: a⃗ (cartesian vector) = k⃗ × v⃗ (transverse vector)
"""
function kx_t2c(H,mn,kpg_mag)
	kxscales = [-1.; 1.]
	kxinds = [2; 1]
    @tullio d[a,i,j,k] := kxscales[b] * H[kxinds[b],i,j,k] * mn[a,b,i,j,k] * kpg_mag[i,j,k] nograd=(kxscales,kxinds) fastmath=false
end

"""
    kx_c2t: v⃗ (transverse vector) = k⃗ × a⃗ (cartesian vector)
"""
function kx_c2t(e⃗,mn,kpg_mag)
	kxscales = [-1.; 1.]
    kxinds = [2; 1]
    @tullio H[b,i,j,k] := kxscales[b] * e⃗[a,i,j,k] * mn[a,kxinds[b],i,j,k] * kpg_mag[i,j,k] nograd=(kxinds,kxscales) fastmath=false
end

"""
    kxinv_t2c: compute a⃗ (cartestion vector) st. v⃗ (cartesian vector from two trans. vector components) ≈ k⃗ × a⃗
    This neglects the component of a⃗ parallel to k⃗ (not available by inverting this cross product)
"""
function kxinv_t2c(H,mn,kpg_mag)
	kxinvscales = [1.; -1.]
	kxinds = [2; 1]
    @tullio e⃗[a,i,j,k] := kxscales[b] * H[kxinds[b],i,j,k] * mn[a,b,i,j,k] / kpg_mag[i,j,k] nograd=(kxscales,kxinds) fastmath=false
end

"""
    kxinv_c2t: compute  v⃗ (transverse 2-vector) st. a⃗ (cartestion 3-vector) = k⃗ × v⃗
    This cross product inversion is exact because v⃗ is transverse (perp.) to k⃗
"""
function kxinv_c2t(d⃗,mn,kpg_mag)
	kxscales = [1.; -1.]
    kxinds = [2; 1]
    @tullio H[b,i,j,k] := kxscales[b] * d⃗[a,i,j,k] * mn[a,kxinds[b],i,j,k] / kpg_mag[i,j,k] nograd=(kxinds,kxscales) fastmath=false
end

"""
    zx_t2c: a⃗ (cartesian vector) = ẑ × v⃗ (transverse vector)
"""
function zx_t2c(H,mn)
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxH[a,i,j,k] := zxscales[a] * H[b,i,j,k] * mn[zxinds[a],b,i,j,k] nograd=(zxscales,zxinds) fastmath=false
end

"""
    ε⁻¹_dot_t: e⃗  = ε⁻¹ ⋅ d⃗ (transverse vectors)
"""
function ε⁻¹_dot_t(d⃗,ε⁻¹)
	@tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * fft(d⃗,(2:4))[b,i,j,k] fastmath=false
	return ifft(e⃗,(2:4))
end

"""
    ε⁻¹_dot: e⃗  = ε⁻¹ ⋅ d⃗ (cartesian vectors)
"""
function ε⁻¹_dot(d⃗,ε⁻¹)
	@tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * d⃗[b,i,j,k] fastmath=false
	# @tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * d⃗[b,i,j,k] / 2 + ε⁻¹[b,a,i,j,k] * d⃗[b,i,j,k] / 2 fastmath=false
end

"""
    ε_dot_approx: approximate     d⃗  = ε ⋅ e⃗
                    using         d⃗  ≈  e⃗ * ( 3 / Tr(ε⁻¹) )
    (all cartesian vectors)
"""
function ε_dot_approx(e⃗,ε⁻¹)
    @tullio d⃗[b,i,j,k] := e⃗[b,i,j,k] * 3 / ε⁻¹[a,a,i,j,k] fastmath=false
end

function M(H,ε⁻¹,mn,kpg_mag)
    -kx_c2t(ε⁻¹_dot_t(kx_t2c(H,mn,kpg_mag),ε⁻¹),mn,kpg_mag)
end

function M(H,ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)
    -kx_c2t( 𝓕⁻¹ * ε⁻¹_dot( 𝓕 * kx_t2c(H,mn,kpg_mag), ε⁻¹), mn,kpg_mag)
end

function M(Hin::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag)::Array{ComplexF64,1}
    HinA = reshape(Hin,(2,size(ε⁻¹)[end-2:end]...))
    return vec(M(HinA,ε⁻¹,mn,kpg_mag))
end

function M(Hin::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::Array{ComplexF64,1}
    HinA = reshape(Hin,(2,size(ε⁻¹)[end-2:end]...))
    return vec(M(HinA,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹))
end

M̂(ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> M(H,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹)::AbstractArray{ComplexF64,1},*(2,size(ε⁻¹)[end-2:end]...),ishermitian=true,ismutating=false)

function P(Hin::AbstractArray{ComplexF64,4},ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::AbstractArray{ComplexF64,4}
    kxinv_c2t( 𝓕 * ε_dot_approx( 𝓕⁻¹ * kxinv_t2c(H,mn,kpg_mag),ε⁻¹),mn,kpg_mag)
end

function P(Hin::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹)::AbstractArray{ComplexF64,1}
    HinA = reshape(Hin,(2,size(ε⁻¹)[end-2:end]...))
    return vec(P(HinA,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹))
end

P̂(ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> P(H,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹)::AbstractArray{ComplexF64,1},*(2,size(ε⁻¹)[end-2:end]...),ishermitian=true,ismutating=false)

# function Mₖ(H,ε⁻¹,mn,kpg_mag)
function Mₖ(H::T,ε⁻¹,mn,kpg_mag)::T where T<:AbstractArray{ComplexF64,4}
    kx_c2t(ε⁻¹_dot_t(zx_t2c(H,mn),ε⁻¹),mn,kpg_mag)
end

function Mₖ(H,ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)
    kx_c2t( 𝓕⁻¹ * ε⁻¹_dot( 𝓕 * zx_t2c(H,mn), ε⁻¹), mn,kpg_mag)
end

# function Mₖ(H::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag)::AbstractArray{ComplexF64,1}
function Mₖ(H::T,ε⁻¹,mn,kpg_mag)::T where T<:AbstractArray{ComplexF64,1}
    Ha = reshape(H,(2,size(ε⁻¹)[end-2:end]...))
    return vec(Mₖ(Ha,ε⁻¹,mn,kpg_mag))
end

function Mₖ(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{T,5},mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::Array{ComplexF64,1} where T
    Ha = reshape(H,(2,size(ε⁻¹)[end-2:end]...))
    return vec(Mₖ(Ha,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹))
end

M̂ₖ(ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> Mₖ(H,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹)::AbstractArray{ComplexF64,1},*(2,size(ε⁻¹)[end-2:end]...),ishermitian=true,ismutating=false)

function H_Mₖ_H(H::AbstractArray{ComplexF64,4},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn)
	kxinds = [2; 1]
	kxscales = [-1.; 1.]
	# temp1 = zx_t2c(H,mn)
	# temp2 = ε⁻¹_dot_t(temp1,ε⁻¹)
	# @tullio out[_] := conj.(H)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * temp[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) nograd=(kxscales,kxinds) fastmath=false
	@tullio out := conj.(H)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * ε⁻¹_dot_t(zx_t2c(H,mn),ε⁻¹)[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) nograd=(kxscales,kxinds) fastmath=false
	return abs(out[1])
	# -real( dot(H, kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(H,mn), (2:4) ), ε⁻¹), (2:4)), mn,kpg_mag) ) )
end

# function H_Mₖ_H(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{ComplexF64,5},kpg_mag,mn,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::Float64
# 	kxinds = [2; 1]
# 	kxscales = [-1.; 1.]
# 	@tullio out[_] := conj.(H)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * ( 𝓕⁻¹ * ε⁻¹_dot( 𝓕 * zx_t2c(H,mn), ε⁻¹) )[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) fastmath=false
# 	return abs(out[1])
# end

# H_Mₖ_H(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{ComplexF64,5},kpg_mag,mn,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)  = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn,𝓕,𝓕⁻¹)
H_Mₖ_H(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn) = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn)
# H_Mₖ_H(H::AbstractArray{ComplexF64,2},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)  = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn,𝓕,𝓕⁻¹)
H_Mₖ_H(H::AbstractArray{ComplexF64,2},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn) = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn)

###### Mutating Operators #######

function t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ds.e[1,i,j,k] = ( Hin[1,i,j,k] * ds.mn[1,1,i,j,k] + Hin[2,i,j,k] * ds.mn[1,2,i,j,k] )
        @fastmath @inbounds ds.e[2,i,j,k] = ( Hin[1,i,j,k] * ds.mn[2,1,i,j,k] + Hin[2,i,j,k] * ds.mn[2,2,i,j,k] )
    	@fastmath @inbounds ds.e[3,i,j,k] = ( Hin[1,i,j,k] * ds.mn[3,1,i,j,k] + Hin[2,i,j,k] * ds.mn[3,2,i,j,k] )
	end
    return ds.e
end

function c2t!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ds.e[1,i,j,k] =  Hin[1,i,j,k] * ds.mn[1,1,i,j,k] + Hin[2,i,j,k] * ds.mn[2,1,i,j,k] + Hin[3,i,j,k] * ds.mn[3,1,i,j,k]
        @fastmath @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.mn[1,2,i,j,k] + Hin[2,i,j,k] * ds.mn[2,2,i,j,k] + Hin[3,i,j,k] * ds.mn[3,2,i,j,k]
    end
    return ds.e
end

function zcross_t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ds.e[1,i,j,k] = -Hin[1,i,j,k] * ds.mn[2,1,i,j,k] - Hin[2,i,j,k] * ds.mn[2,2,i,j,k]
        @fastmath @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.mn[1,1,i,j,k] + Hin[2,i,j,k] * ds.mn[1,2,i,j,k]
    end
    return ds.e
end

function kcross_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ds.d[1,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[1,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[1,1,i,j,k] ) * -ds.kpg_mag[i,j,k]
        @fastmath @inbounds ds.d[2,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[2,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[2,1,i,j,k] ) * -ds.kpg_mag[i,j,k]
        @fastmath @inbounds ds.d[3,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[3,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[3,1,i,j,k] ) * -ds.kpg_mag[i,j,k]
    end
    return ds.d
end

function kcross_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
		@fastmath @inbounds  ds.H[1,i,j,k] =  (	ds.e[1,i,j,k] * ds.mn[1,2,i,j,k] + ds.e[2,i,j,k] * ds.mn[2,2,i,j,k] + ds.e[3,i,j,k] * ds.mn[3,2,i,j,k]	) * -ds.kpg_mag[i,j,k]
		@fastmath @inbounds  ds.H[2,i,j,k] =  (	ds.e[1,i,j,k] * ds.mn[1,1,i,j,k] + ds.e[2,i,j,k] * ds.mn[2,1,i,j,k] + ds.e[3,i,j,k] * ds.mn[3,1,i,j,k]	) * ds.kpg_mag[i,j,k]
    end
    return ds.H
end

function ε⁻¹_dot!(ε⁻¹::Array{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ds.e[1,i,j,k] =  ε⁻¹[1,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,1,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,1,i,j,k]*ds.d[3,i,j,k]
        @fastmath @inbounds ds.e[2,i,j,k] =  ε⁻¹[1,2,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,2,i,j,k]*ds.d[3,i,j,k]
        @fastmath @inbounds ds.e[3,i,j,k] =  ε⁻¹[1,3,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,3,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,3,i,j,k]*ds.d[3,i,j,k]
        # ds.e[1,i,j,k] =  ε⁻¹[1,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[1,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[1,3,i,j,k]*ds.d[3,i,j,k]
        # ds.e[2,i,j,k] =  ε⁻¹[2,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[2,3,i,j,k]*ds.d[3,i,j,k]
        # ds.e[3,i,j,k] =  ε⁻¹[3,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[3,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,3,i,j,k]*ds.d[3,i,j,k]
    end
    return ds.e
end

function M!(ε⁻¹::Array{Float64,5},ds::MaxwellData)::Array{ComplexF64,4}
    kcross_t2c!(ds);
    mul!(ds.d,ds.𝓕!,ds.d);
    ε⁻¹_dot!(ε⁻¹,ds);
	mul!(ds.e,ds.𝓕⁻¹!::AbstractFFTs.ScaledPlan,ds.e);
    kcross_c2t!(ds)
end

function M!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{Float64,5},ds::MaxwellData)::Array{ComplexF64,1}
    @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
    M!(ε⁻¹,ds);
    @inbounds Hout .= vec(ds.H)
end

M̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData) = LinearMap{ComplexF64}((2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true) do y::AbstractVector{ComplexF64},x::AbstractVector{ComplexF64}
    M!(y,x,ε⁻¹,ds)::AbstractArray{ComplexF64,1}
    end

# function M̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData)
#     function f!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         M!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(f!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end

### Preconditioner P̂ & Component Operators (approximate inverse operations of M̂) ###

function kcrossinv_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds scale::Float64 = inv(ds.kpg_mag[i,j,k])
        @fastmath @inbounds ds.e[1,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[1,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[1,1,i,j,k] ) * scale
        @fastmath @inbounds ds.e[2,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[2,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[2,1,i,j,k] ) * scale
        @fastmath @inbounds ds.e[3,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[3,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[3,1,i,j,k] ) * scale
    end
    return ds.e
end

function kcrossinv_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds scale = -inv(ds.kpg_mag[i,j,k])
        @fastmath @inbounds ds.H[1,i,j,k] =  (	ds.d[1,i,j,k] * ds.mn[1,2,i,j,k] + ds.d[2,i,j,k] * ds.mn[2,2,i,j,k] + ds.d[3,i,j,k] * ds.mn[3,2,i,j,k]	) * -scale
        @fastmath @inbounds ds.H[2,i,j,k] =  (	ds.d[1,i,j,k] * ds.mn[1,1,i,j,k] + ds.d[2,i,j,k] * ds.mn[2,1,i,j,k] + ds.d[3,i,j,k] * ds.mn[3,1,i,j,k]	) * scale
    end
    return ds.H
end

function ε_dot_approx!(ε⁻¹::AbstractArray{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ε_ave = 3. * inv( ε⁻¹[1,1,i,j,k] + ε⁻¹[2,2,i,j,k] + ε⁻¹[3,3,i,j,k] ) # tr(ε⁻¹[:,:,i,j,k])
        @fastmath @inbounds ds.d[1,i,j,k] =  ε_ave * ds.e[1,i,j,k]
        @fastmath @inbounds ds.d[2,i,j,k] =  ε_ave * ds.e[2,i,j,k]
        @fastmath @inbounds ds.d[3,i,j,k] =  ε_ave * ds.e[3,i,j,k]
    end
    return ds.d
end

function P!(ε⁻¹::AbstractArray{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    kcrossinv_t2c!(ds);
    # ds.𝓕⁻¹! * ds.e;
    # ldiv!(ds.e,ds.𝓕!,ds.e)
	mul!(ds.e,ds.𝓕⁻¹!,ds.e);
    ε_dot_approx!(ε⁻¹,ds);
    # ds.𝓕! * ds.d;
    mul!(ds.d,ds.𝓕!,ds.d);
    kcrossinv_c2t!(ds)
end

function P!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,1}
    @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
    P!(ε⁻¹,ds);
    @inbounds Hout .= vec(ds.H)
end

P̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData) = LinearMap{ComplexF64}((2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true) do y::AbstractVector{ComplexF64},x::AbstractVector{ComplexF64}
	P!(y,x,ε⁻¹,ds)::AbstractArray{ComplexF64,1}
    end


text:
### Extract Terms

"""
extract_patterns: Extract Patterns from Text
```julia
extract_patterns(s,:puncts)

```
"""
function extract_patterns(text::AbstractString, pattern)
        return collect(x.match for x in eachmatch(REGEX_DICT[pattern], text))
end

"""
extract_numbers: Extract all numbers from a given text
"""
function extract_numbers(text::AbstractString)
        res = [m.match for m in eachmatch(NUMBERS_REGEX,text)]
        return res
    end

"""
extract_phonenumbers: Extract all phone numbers from a given text
"""
function extract_phonenumbers(text::AbstractString)
        res = [m.match for m in eachmatch(PHONENUMBERS_REGEX,text)]
        return res
    end

"""
extract_currencies: Extract all currencies from a given text
"""
function extract_currencies(text::AbstractString)
        res = [m.match for m in eachmatch(CURRENCY_REGEX,text)]
        return res
    end

"""
extract_btc: Extract all Bitcoin BTC Address/Values from a given text
"""
function extract_btc(text::AbstractString)
        res = [m.match for m in eachmatch(BTC_ADDRESS_REGEX,text)]
        return res
    end

"""
extract_mastercards: Extract all MasterCard Address/Values from a given text
"""
function extract_mastercards(text::AbstractString)
        res = [m.match for m in eachmatch(MASTERCard_REGEX,text)]
        return res
    end

"""
extract_visacards: Extract all VisaCards Address/Values from a given text
"""
function extract_visacards(text::AbstractString)
        res = [m.match for m in eachmatch(VISACard_REGEX,text)]
        return res
    end

"""
extract_emails: Extract all emails from a given text
"""
function extract_emails(text::AbstractString)
        res = [m.match for m in eachmatch(EMAIL_REGEX,text)]
        return res
    end

"""
extract_urls: Extract all urls from a given text
"""
function extract_urls(text::AbstractString)
        res = [m.match for m in eachmatch(URL_REGEX,text)]
        return res
    end

"""
extract_puncts: Extract all punctuations from a given text
"""
function extract_puncts(text::AbstractString)
        res = [m.match for m in eachmatch(PUNCTS_REGEX,text)]
        return res
    end

"""
extract_special_characters: Extract all special characters from a given text
"""
function extract_special_characters(text::AbstractString)
        res = [m.match for m in eachmatch(SPECIAL_CHARACTERS_REGEX,text)]
        return res
    end

"""
extract_emojis: Extract all emojis from a given text
"""
function extract_emojis(text::AbstractString)
        res = [m.match for m in eachmatch(EMOJI_REGEX,text)]
        return res
    end

"""
extract_userhandles: Extract all userhandles and mentions from a given text
"""
function extract_userhandles(text::AbstractString)
        res = [m.match for m in eachmatch(USER_HANDLES_REGEX,text)]
        return res
    end

"""
extract_hashtags: Extract all hashtags from a given text
"""
function extract_hashtags(text::AbstractString)
        res = [m.match for m in eachmatch(HASHTAG_REGEX,text)]
        return res
    end

"""
extract_userhandles: Extract all dates from a given text
"""
function extract_dates(text::AbstractString)
        res = [m.match for m in eachmatch(DATE_REGEX,text)]
        return res
    end

"""
extract_pobox: Extract all Post Office Box from a given text
"""
function extract_pobox(text::AbstractString)
        res = [m.match for m in eachmatch(PoBOX_REGEX,text)]
        return res
    end

"""
extract_streetaddress: Extract all Street Address from a given text
"""
function extract_streetaddress(text::AbstractString)
        res = [m.match for m in eachmatch(STREET_ADDRESS_REGEX,text)]
        return res
    end

"""
extract_md5sha: Extract all MD5 or SHA Hashes from a given text
"""
function extract_md5sha(text::AbstractString)
        res = [m.match for m in eachmatch(MD5_SHA_REGEX,text)]
        return res
    end

# returns RegexMatch Type
# function extract_numbers(s::AbstractString)
#            res = eachmatch(r"\d+",s)
#            return collect(res)
# end

function extract_terms_in_bracket(text::AbstractString;bracket_form="curly")
        CURLY_BRACKET_REGEX = r"\{(.*?)\}"
        SQUARE_BRACKET_REGEX = r"\[(.*?)\]"
        NORMAL_BRACKET_REGEX = r"\((.*?)\)"
        if (bracket_form == "curly") | (bracket_form == "{}")
            result = [m.match for m in eachmatch(CURLY_BRACKET_REGEX,text)]
        elseif (bracket_form == "square") | (bracket_form == "[]")
            result = [m.match for m in eachmatch(SQUARE_BRACKET_REGEX,text)]
        elseif (bracket_form == "normal") | (bracket_form == "()")
            result = [m.match for m in eachmatch(NORMAL_BRACKET_REGEX,text)]
           end
        return result
end

"""
extract_stopwords: returns a list of stopwords in the given text
param:
-----
 - lang: specify the language for the stopwords : en,es,fr,ru,yo,de
"""
function extract_stopwords(text::AbstractString,lang="en")
        if lang == "en"
                stopwords_in_use = STOPWORDS_en
    elseif lang == "es"
        stopwords_in_use = STOPWORDS_es
    elseif lang == "fr"
        stopwords_in_use = STOPWORDS_fr
    elseif lang == "ru"
        stopwords_in_use = STOPWORDS_ru
    elseif lang == "yo"
        stopwords_in_use = STOPWORDS_yo
    elseif lang == "de"
        stopwords_in_use = STOPWORDS_de
    else
        stopwords_in_use = STOPWORDS_en
    end

        alltokens = []
        for w in split(lowercase(text)," ")
                if w in stopwords_in_use
                        push!(alltokens,w)
                end
        end
        return Set(alltokens)

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RunCommandInput <: SwaggerModel
    commandId::Any # spec type: Union{ Nothing, String } # spec name: commandId
    script::Any # spec type: Union{ Nothing, Vector{String} } # spec name: script
    parameters::Any # spec type: Union{ Nothing, Vector{RunCommandInputParameter} } # spec name: parameters

    function RunCommandInput(;commandId=nothing, script=nothing, parameters=nothing)
        o = new()
        validate_property(RunCommandInput, Symbol("commandId"), commandId)
        setfield!(o, Symbol("commandId"), commandId)
        validate_property(RunCommandInput, Symbol("script"), script)
        setfield!(o, Symbol("script"), script)
        validate_property(RunCommandInput, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        o
    end
end # type RunCommandInput

const _property_map_RunCommandInput = Dict{Symbol,Symbol}(Symbol("commandId")=>Symbol("commandId"), Symbol("script")=>Symbol("script"), Symbol("parameters")=>Symbol("parameters"))
const _property_types_RunCommandInput = Dict{Symbol,String}(Symbol("commandId")=>"String", Symbol("script")=>"Vector{String}", Symbol("parameters")=>"Vector{RunCommandInputParameter}")
Base.propertynames(::Type{ RunCommandInput }) = collect(keys(_property_map_RunCommandInput))
Swagger.property_type(::Type{ RunCommandInput }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RunCommandInput[name]))}
Swagger.field_name(::Type{ RunCommandInput }, property_name::Symbol) =  _property_map_RunCommandInput[property_name]

function check_required(o::RunCommandInput)
    (getproperty(o, Symbol("commandId")) === nothing) && (return false)
    true
end

function validate_property(::Type{ RunCommandInput }, name::Symbol, val)
end


text:
# categorical distribution sample

module CategoricalDistribution

  function main()
    # dice
    pi_ = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6]
    c = [0, 1, 0, 0, 0, 0] # get 2, 1-of-K coding
    println(cat(c, pi_))
  end

  function cat(c, pi_)
    elems = [pik_^ck for (ck, pik_) in zip(c, pi_)]
    return prod(elems)
  end

end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
    AbstractConcatBasis

is an abstract type representing a block diagonal basis but
with modified axes.
"""
abstract type AbstractConcatBasis{T} <: Basis{T} end

copy(S::AbstractConcatBasis) = S

@simplify function *(D::Derivative, S::AbstractConcatBasis)
    axes(D,2) == axes(S,1) || throw(DimensionMismatch())
    args = arguments.(Ref(ApplyLayout{typeof(*)}()), Derivative.(axes.(S.args,1)) .* S.args)
    all(length.(args) .== 2) || error("Not implemented")
    concatbasis(S, map(first, args)...) * mortar(Diagonal([map(last, args)...]))
end

"""
    PiecewiseBasis(args...)

is an analogue of `Basis` that takes the union of the first axis,
and the second axis is a blocked concatenatation of args.
If there is overlap, it uses the first in order.
"""
struct PiecewiseBasis{T, Args} <: AbstractConcatBasis{T}
    args::Args
end

PiecewiseBasis{T}(args...) where T = PiecewiseBasis{T,typeof(args)}(args)
PiecewiseBasis(args...) = PiecewiseBasis{mapreduce(eltype,promote_type,args)}(args...)
PiecewiseBasis{T}(args::AbstractVector) where T = PiecewiseBasis{T,typeof(args)}(args)
PiecewiseBasis(args::AbstractVector) = PiecewiseBasis{eltype(eltype(args))}(args)

concatbasis(::PiecewiseBasis, args...) = PiecewiseBasis(args...)

axes(A::PiecewiseBasis) = (union(axes.(A.args,1)...), blockedrange(size.(A.args,2)))

==(A::PiecewiseBasis, B::PiecewiseBasis) = all(A.args .== B.args)

function QuasiArrays._getindex(::Type{IND}, A::PiecewiseBasis{T}, (x,j)::IND) where {IND,T}
    Jj = findblockindex(axes(A,2), j)
    J = Int(block(Jj))
    @boundscheck x in axes(A,1) || throw(BoundsError(A, (x,j)))
    x in axes(A.args[J],1) && return A.args[J][x, blockindex(Jj)]
    zero(T)
end

"""
    VcatBasis

is an analogue of `Basis` that vcats the values.
"""
struct VcatBasis{T, Args} <: AbstractConcatBasis{T}
    args::Args
    function VcatBasis{T, Args}(args::Args) where {T,Args}
        ax = axes(args[1],1)
        for a in Base.tail(args)
            ax == axes(a,1) || throw(ArgumentError("Must be defined on same"))
        end
        new{T,Args}(args)
    end
end

VcatBasis{T}(args...) where T = VcatBasis{T,typeof(args)}(args)
VcatBasis(args::Vararg{Any,N}) where N = VcatBasis{SVector{N,mapreduce(eltype,promote_type,args)}}(args...)

concatbasis(::VcatBasis, args...) = VcatBasis(args...)

axes(A::VcatBasis{<:Any,<:Tuple}) = (axes(A.args[1],1), blockedrange(size.(A.args,2)))

==(A::VcatBasis, B::VcatBasis) = all(A.args .== B.args)

function QuasiArrays._getindex(::Type{IND}, A::VcatBasis{T}, (x,j)::IND) where {IND,T<:SVector{N}} where N
    Jj = findblockindex(axes(A,2), j)
    J = Int(block(Jj))
    T(zeros(J-1)..., A.args[J][x, blockindex(Jj)], zeros(N-J)...) # TODO: type stable via @generated
end

"""
    VcatBasis

is an analogue of `Basis` that hvcats the values, so they are matrix valued.
"""
struct HvcatBasis{T, Args} <: AbstractConcatBasis{T}
    n::Int
    args::Args
    function HvcatBasis{T, Args}(n::Int, args::Args) where {T,Args}
        ax = axes(args[1],1)
        for a in Base.tail(args)
            ax == axes(a,1) || throw(ArgumentError("Must be defined on same"))
        end
        new{T,Args}(n, args)
    end
end

HvcatBasis{T}(n::Int, args...) where T = HvcatBasis{T,typeof(args)}(n, args)
HvcatBasis(n::Int, args::Vararg{Any,N}) where N = HvcatBasis{Matrix{mapreduce(eltype,promote_type,args)}}(n, args...)

concatbasis(S::HvcatBasis, args...) = HvcatBasis(S.n, args...)

axes(A::HvcatBasis{<:Any,<:Tuple}) = (axes(A.args[1],1), blockedrange(size.(A.args,2)))

==(A::HvcatBasis, B::HvcatBasis) = all(A.args .== B.args)

function QuasiArrays._getindex(::Type{IND}, A::HvcatBasis{T}, (x,j)::IND) where {T,IND}
    Jj = findblockindex(axes(A,2), j)
    J = Int(block(Jj))
    hvcat(A.n, zeros(J-1)..., A.args[J][x, blockindex(Jj)], zeros(length(A.args)-J)...)::T
end

@simplify *(D::Derivative, H::ApplyQuasiMatrix{<:Any,typeof(hcat)}) = hcat((Ref(D) .* H.args)...)


text:
@inline splitview(a::Tuple...) = splitview(Float64,a)
@inline splitview(::Type{T},a::Tuple...) where {T} = splitview(T,a)
@inline splitview(x::AbstractVector,a::Tuple...) = splitview(x,a)

@generated function splitview(::Type{T},arr::NTuple{N,Tuple}) where {T,N}
    quote
      len = 0
      @nexprs $N (n)->(s_n = len+1; e_n = len+prod(arr[n]); len = e_n)
      x = Array{T}(undef, len)
      X = @ntuple $N (n)->(reshape(view(x,s_n:e_n),arr[n]))
      start = @ntuple $N (n)->(s_n)
      stop = @ntuple $N (n)->(e_n)
      return x,X,start,stop
    end
end

@generated function splitview(x::AbstractVector{T},arr::NTuple{N,Tuple}) where {T,N}
    quote
      len = 0
      @nexprs $N (n)->(s_n = len+1; e_n = len+prod(arr[n]); len = e_n)
      X = @ntuple $N (n)->(reshape(view(x,s_n:e_n),arr[n]))
      start = @ntuple $N (n)->(s_n)
      stop = @ntuple $N (n)->(e_n)
      return X,start,stop
    end
end

# unexported from Base
const FastContiguousSubArray{T,N,P,I<:Tuple{Union{Colon, UnitRange}, Vararg{Any}}} = SubArray{T,N,P,I,true}

parentindex(A::FastContiguousSubArray, i::Int) = A.offset1 + i
parentindex(A::FastContiguousSubArray, i::Int...) = A.offset1 + prod(i)
parentindex(A::FastContiguousSubArray, i::NTuple{N,Int}) where {N} = A.offset1 + prod(i)

@inline parentindex(A::ReshapedArray, i::Int...) = LinearIndices(size(A))[i...]
parentindex(A::ReshapedArray, i::NTuple{N,Int}) where {N} = LinearIndices(size(A))[i...]

vecidx(A::Array, i::Int) = i

vecidx(A::FastContiguousSubArray, i::Int) = vecidx(A.parent, parentindex(A, i))
@inline vecidx(A::FastContiguousSubArray, i::Int...) = vecidx(A.parent, parentindex(A, i...))
vecidx(A::FastContiguousSubArray, i::NTuple{N,Int}) where {N} = vecidx(A.parent, parentindex(A, i))

@inline vecidx(A::ReshapedArray, i::Int...) = vecidx(A.parent, parentindex(A, i...))
vecidx(A::ReshapedArray, i::NTuple{N,Int}) where {N} = vecidx(A.parent, parentindex(A, i))


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_opt_rho1_9(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_opt_rho1_9!(A_copy)
end

@inline function exp_sastre_opt_rho1_9!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=-4.479293476259415
    coeff2=-1.965591461532149
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-4.479293476259415
    coeff2=-1.965591461532149
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.5431852019042509
    coeff2=2.808044477000553
    coeff3=3.072468368540587
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=1.8758020055089961
    coeff2=-4.077994868395093
    coeff3=-0.5726777535092064
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=3.513744518233991
    coeff2=3.480716229698247
    coeff3=-7.395064443323127
    coeff4=3.9135175595994687
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.02417939999948854
    coeff2=0.013482072269100465
    coeff3=-0.0015522395764164978
    coeff4=-2.310576048441552e-6
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-2.7514538554861674
    coeff2=28.533085265354913
    coeff3=-32.70207976001747
    coeff4=15.117957339910209
    coeff5=-31.69337657641421
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=1.4649998087002614
    coeff2=-1.1631539168725962
    coeff3=18.10426755193074
    coeff4=-13.187053051683513
    coeff5=27.426929433347695
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-2.57802584497714
    coeff2=-1.9681270306651504
    coeff3=-0.5172562598353871
    coeff4=0.2838459715474883
    coeff5=-0.0930303538906508
    coeff6=-0.17367745193891393
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=2.7927443516334343
    coeff2=-2.5111253025765686e-6
    coeff3=1.1720013130685409e-7
    coeff4=2.1593312631548876e-10
    coeff5=2.7906362162645103e-10
    coeff6=-5.055543522880249e-19
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=9.679532976243642
    coeff2=3.9359833423253106
    coeff3=6.090302916482898
    coeff4=-2.7031124428719058
    coeff5=-2.542481119421314
    coeff6=1.6620101377768384
    coeff7=3.4265640842763014
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
"""
Attitude objects contains Euler angles (rad) and quaternions.

# Constructors
    Attitude(psi, theta, phi)
    Attitude(q0, q1, q2, q3)
"""
struct Attitude
    # Euler angles
    psi::Number
    theta::Number
    phi::Number
    # Quaternions
    q0::Number
    q1::Number
    q2::Number
    q3::Number
end

show(io::IO, att::Attitude) = print(
    io,
    "Attitude: ψ=$(rad2deg(att.psi))º, θ=$(rad2deg(att.theta))º, ϕ=$(rad2deg(att.phi))º",
    )

# Constructor from Euler angles
Attitude(psi, theta, phi) = Attitude(psi, theta, phi,
                                     euler2quaternion(psi, theta, phi)...)

# Constructor from quaternions
Attitude(q0, q1, q2, q3) = Attitude(quaternion2euler(q0, q1, q2, q3)...,
                                    q0, q1, q2, q3)

# Getters
"""
    get_euler_angles(att::Attitude)

Return Euler angles ψ, θ, ϕ (yaw, pith, roll).
"""
get_euler_angles(att::Attitude) = [att.psi, att.theta, att.phi]

"""
    get_quaternions(att::Attitude)

Return quaternions.
"""
get_quaternions(att::Attitude) = [att.q0, att.q1, att.q2, att.q3]

function isapprox(x::Attitude, y::Attitude; rtol=1e-8, atol=0.0, nans=false)

    return isapprox(
        [x.psi, x.theta, x.phi, x.q0, x.q1, x.q2, x.q3],
        [y.psi, y.theta, y.phi, y.q0, y.q1, y.q2, y.q3],
        rtol=rtol, atol=atol, nans=nans,
        )
end


text:
using FactCheck, MarketData
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("last methods takes the first observations") do

    context("defaults to n=1") do
        @fact last(cl).values[1] --> 21.9
    end

    context("takes other n values") do
        @fact last(cl, 2).values[1] --> 22.43
    end

end