text:
export PutBlock, put

"""
    PutBlock <: AbstractContainer

Type for putting a block at given locations.
"""
struct PutBlock{N, C, GT <: AbstractBlock} <: AbstractContainer{GT, N}
    content::GT
    locs::NTuple{C, Int}

    function PutBlock{N}(block::GT, locs::NTuple{C, Int}) where {N, M, C, GT <: AbstractBlock{M}}
        @assert_locs_safe N locs
        @assert nqubits(block) == C "number of locations doesn't match the size of block"
        return new{N, C, GT}(block, locs)
    end
end

"""
    put(total::Int, pair)

Create a [`PutBlock`](@ref) with total number of active qubits, and a pair of
location and block to put on.

# Example

```jldoctest
julia> put(4, 1=>X)
nqubits: 4
put on (1)
└─ X gate
```

If you want to put a multi-qubit gate on specific locations, you need to write down all possible locations.

```jldoctest
julia> put(4, (1, 3)=>kron(X, Y))
nqubits: 4
put on (1, 3)
└─ kron
   ├─ 1=>X gate
   └─ 2=>Y gate
```

The outter locations creates a scope which make it seems to be a contiguous two qubits for the block inside `PutBlock`.

!!! tips
    It is better to use [`concentrate`](@ref) instead of `put` for large blocks, since put will use the matrix of its contents
    directly instead of making use of what's in it. `put` is more efficient for small blocks.
"""
put(total::Int, pa::Pair{NTuple{M, Int}, <:AbstractBlock}) where M =
    PutBlock{total}(pa.second, pa.first)
put(total::Int, pa::Pair{Int, <:AbstractBlock}) = PutBlock{total}(pa.second, (pa.first, ))
put(total::Int, pa::Pair{<:Any, <:AbstractBlock}) = PutBlock{total}(pa.second, Tuple(pa.first))

"""
    put(pair) -> f(n)

Lazy curried version of [`put`](@ref).

# Example

```jldoctest
julia> put(1=>X)
(n -> put(n, 1 => X gate))
```
"""
put(pa::Pair) = @λ(n -> put(n, pa))

occupied_locs(x::PutBlock) = x.locs
chsubblocks(x::PutBlock{N}, b::AbstractBlock) where N = PutBlock{N}(b, x.locs)
PreserveStyle(::PutBlock) = PreserveAll()
cache_key(pb::PutBlock) = cache_key(pb.content)

mat(::Type{T}, pb::PutBlock{N, 1}) where {T, N} = u1mat(N, mat(T, pb.content), pb.locs...)
mat(::Type{T}, pb::PutBlock{N, C}) where {T, N, C} = unmat(N, mat(T, pb.content), pb.locs)

function apply!(r::ArrayReg{B, T}, pb::PutBlock{N}) where {B, T, N}
    _check_size(r, pb)
    instruct!(matvec(r.state), mat(T, pb.content), pb.locs)
    return r
end

# specialization
for G in [:X, :Y, :Z, :T, :S, :Sdag, :Tdag]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))
    @eval function apply!(r::ArrayReg, pb::PutBlock{N, C, <:$GT}) where {N, C}
        _check_size(r, pb)
        instruct!(matvec(r.state), Val($(QuoteNode(G))), pb.locs)
        return r
    end
end

Base.adjoint(x::PutBlock{N}) where N = PutBlock{N}(adjoint(content(x)), x.locs)
Base.copy(x::PutBlock{N}) where N = PutBlock{N}(x.content, x.locs)
function Base.:(==)(lhs::PutBlock{N, C, GT}, rhs::PutBlock{N, C, GT}) where {N, C, GT}
    return (lhs.content == rhs.content) && (lhs.locs == rhs.locs)
end

function YaoBase.iscommute(x::PutBlock{N}, y::PutBlock{N}) where N
    if x.locs == y.locs
        return iscommute(x.content, y.content)
    else
        return iscommute_fallback(x, y)
    end
end


text:
# Autogenerated wrapper script for Trilinos_jll for x86_64-w64-mingw32-libgfortran3-cxx03
export libamesos, libaztecoo, libbelos, libbelosepetra, libepetra, libepetraext, libifpack, libisorropia, libloca, liblocaepetra, liblocalapack, libnox, libnoxepetra, libnoxlapack, libsacado, libsimpi, libteuchoscomm, libteuchoscore, libteuchosnumerics, libteuchosparameterlist, libteuchosremainder, libtrilinoscouplings, libtrilinosss, libtriutils, libzoltan

using SuiteSparse_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Trilinos")
JLLWrappers.@declare_library_product(libamesos, "libamesos.dll")
JLLWrappers.@declare_library_product(libaztecoo, "libaztecoo.dll")
JLLWrappers.@declare_library_product(libbelos, "libbelos.dll")
JLLWrappers.@declare_library_product(libbelosepetra, "libbelosepetra.dll")
JLLWrappers.@declare_library_product(libepetra, "libepetra.dll")
JLLWrappers.@declare_library_product(libepetraext, "libepetraext.dll")
JLLWrappers.@declare_library_product(libifpack, "libifpack.dll")
JLLWrappers.@declare_library_product(libisorropia, "libisorropia.dll")
JLLWrappers.@declare_library_product(libloca, "libloca.dll")
JLLWrappers.@declare_library_product(liblocaepetra, "liblocaepetra.dll")
JLLWrappers.@declare_library_product(liblocalapack, "liblocalapack.dll")
JLLWrappers.@declare_library_product(libnox, "libnox.dll")
JLLWrappers.@declare_library_product(libnoxepetra, "libnoxepetra.dll")
JLLWrappers.@declare_library_product(libnoxlapack, "libnoxlapack.dll")
JLLWrappers.@declare_library_product(libsacado, "libsacado.dll")
JLLWrappers.@declare_library_product(libsimpi, "libsimpi.dll")
JLLWrappers.@declare_library_product(libteuchoscomm, "libteuchoscomm.dll")
JLLWrappers.@declare_library_product(libteuchoscore, "libteuchoscore.dll")
JLLWrappers.@declare_library_product(libteuchosnumerics, "libteuchosnumerics.dll")
JLLWrappers.@declare_library_product(libteuchosparameterlist, "libteuchosparameterlist.dll")
JLLWrappers.@declare_library_product(libteuchosremainder, "libteuchosremainder.dll")
JLLWrappers.@declare_library_product(libtrilinoscouplings, "libtrilinoscouplings.dll")
JLLWrappers.@declare_library_product(libtrilinosss, "libtrilinosss.dll")
JLLWrappers.@declare_library_product(libtriutils, "libtriutils.dll")
JLLWrappers.@declare_library_product(libzoltan, "libzoltan.dll")
function __init__()
    JLLWrappers.@generate_init_header(SuiteSparse_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libamesos,
        "bin\\libamesos.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libaztecoo,
        "bin\\libaztecoo.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbelos,
        "bin\\libbelos.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbelosepetra,
        "bin\\libbelosepetra.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libepetra,
        "bin\\libepetra.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libepetraext,
        "bin\\libepetraext.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libifpack,
        "bin\\libifpack.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libisorropia,
        "bin\\libisorropia.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libloca,
        "bin\\libloca.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblocaepetra,
        "bin\\liblocaepetra.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblocalapack,
        "bin\\liblocalapack.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnox,
        "bin\\libnox.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnoxepetra,
        "bin\\libnoxepetra.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnoxlapack,
        "bin\\libnoxlapack.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsacado,
        "bin\\libsacado.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsimpi,
        "bin\\libsimpi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libteuchoscomm,
        "bin\\libteuchoscomm.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libteuchoscore,
        "bin\\libteuchoscore.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libteuchosnumerics,
        "bin\\libteuchosnumerics.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libteuchosparameterlist,
        "bin\\libteuchosparameterlist.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libteuchosremainder,
        "bin\\libteuchosremainder.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtrilinoscouplings,
        "bin\\libtrilinoscouplings.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtrilinosss,
        "bin\\libtrilinosss.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtriutils,
        "bin\\libtriutils.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libzoltan,
        "bin\\libzoltan.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
struct RenderPass
    area::Vk.Rect2D
    samples::Int
    # attachments::Vector{}
end

RenderPass(area; samples = 1) = RenderPass(area, samples)
RenderPass(area::NTuple{4,<:Integer}; samples = 1) = RenderPass(Vk.Rect2D(Vk.Offset2D(area[1:2]...), Vk.Extent2D(area[3:4]...)), samples)


text:
module Utils

using Base: llvmcall

# https://github.com/JuliaCI/BenchmarkTools.jl/pull/92
@inline function clobber()
    llvmcall("""
        call void asm sideeffect "", "~{memory}"()
        ret void
    """, Cvoid, Tuple{})
end

end  # module


text:
using JUGRNModelGenerator
using Test

# -- Model creation tests ------------------------------------------------------ #
function build_grn_model_test()
    return true
end
# ------------------------------------------------------------------------------- #

@testset "model_creation_test_set" begin
    @test build_grn_model_test() == true
end


text:
export getData

import jInv.ForwardShare.getData
function getData(m,pFor::BasicFWIparam,doClear::Bool=false)
    tstart = time_ns();
    # extract pointers
    Mesh  = pFor.Mesh
    omega = pFor.omega
    gamma = pFor.gamma
    Q     = pFor.Sources
    P     = pFor.Receivers

    nrec  = size(P,2)
    nsrc  = size(Q,2)
    nfreq = length(omega)

    # allocate space for data and fields
    D  = zeros(nrec,nsrc,nfreq)
    U  = zeros(ComplexF64,prod(Mesh.n.+1),nsrc,nfreq)
    LU = Array{Any}(undef, nfreq)
    HinvtPs = Array{Any}(undef, nfreq)
    for i=1:length(omega)
        t1 = time_ns();
        H = getHelmholtzOperator(m,gamma,omega[i],Mesh)
        t2 = time_ns();
        LU[i] = lu(H)
        H = nothing
        for k=1:nsrc
            U[:,k,i] = LU[i]\Vector(Q[:,k])
            D[:,k,i] = real(P'*U[:,k,i])
        end
    end
    pFor.Ainv   = LU
    pFor.Fields = U # will need to save to file
    tend = time_ns();
    println("Runtime of getData:");
    println((tend - tstart)/1.0e9);

    return D ,pFor
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.23"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	lambda, N = 1/3, 10^5
	bulbs = [1,10,50]
	xGrid = 0:0.1:10
	C = [:blue :red :green]
	dists = [Gamma(n,1/(n*lambda)) for n in bulbs]
end

# ╔═╡ 72a5caf8-2166-11eb-2bbd-41388c36dafb
function normalizedData(d::Gamma)
    sh = Int64(shape(d))
    data = [sum(-(1/(sh*lambda))*log.(rand(sh))) for _ in 1:N]
end

# ╔═╡ 72a5f870-2166-11eb-1c89-4948a68f9e26
L = [ "Shape = "*string.(shape.(i))*", Scale = "*
    string.(round.(scale.(i),digits=2)) for i in dists ]

# ╔═╡ 72a751e8-2166-11eb-14c1-53fba757ccfc
begin
	stephist( normalizedData.(dists), bins=50,
		normed=:true, c=C, xlims=(0,maximum(xGrid)),ylims=(0,1),
		xlabel="x", ylabel="Density", label="")
	plot!(xGrid, [pdf.(i,xGrid) for i in dists], c=C, label=reshape(L, 1,:))
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.23"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═72a5caf8-2166-11eb-2bbd-41388c36dafb
# ╠═72a5f870-2166-11eb-1c89-4948a68f9e26
# ╠═72a751e8-2166-11eb-14c1-53fba757ccfc
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using Revise
using Layered
using Animations
using Colors

function test()

    duration = 6

    a_y = Animation(
        0, 1,
        polyout(2; n=8, yoyo=true),
        duration, 0)

    a_x = Animation(
        0, 0,
        linear(n=2, yoyo=true),
        duration, 1)

    a_color = Animation(
        0, RGB(0.251, 0.388, 0.847),
        duration * 0.25, RGB(0.22, 0.596, 0.149),
        duration * 0.5, RGB(0.584, 0.345, 0.698),
        duration * 0.75, RGB(0.796, 0.235, 0.2),
        defaulteasing=noease()
    )

    function frame(t)
        c, tl = canvas(2, 2)

        l = layer_in_rect!(tl, c.rect, :w, :norm, margin=30)

        circle!(l, P(a_x(t), a_y(t)), 0.15) + Fill(a_color(t)) + Stroke(nothing)

        c

    end

    filepath = joinpath(dirname(pathof(Animations)), "../misc/example.gif")
    record(frame, filepath, 30, duration)

end; test()


text:
import Base: permutedims

const SparseMatrix = SparseMatrixCSC{ComplexF64, Int}

function gemm_sp_dense_small(alpha::ComplexF64, M::SparseMatrix, B::Matrix{ComplexF64}, result::Matrix{ComplexF64})
    if alpha == ComplexF64(1.)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = M.nzval[i]
                @inbounds for j=1:size(B, 2)
                    result[row, j] += val*B[colindex, j]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = alpha*M.nzval[i]
                @inbounds for j=1:size(B, 2)
                    result[row, j] += val*B[colindex, j]
                end
            end
        end
    end
end

function gemm_sp_dense_big(alpha::ComplexF64, M::SparseMatrix, B::Matrix{ComplexF64}, result::Matrix{ComplexF64})
    if alpha == ComplexF64(1.)
        @inbounds for j=1:size(B, 2)
            @inbounds for colindex = 1:M.n
                m2 = B[colindex, j]
                @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                    row = M.rowval[i]
                    result[row, j] += M.nzval[i]*m2
                end
            end
        end
    else
        @inbounds for j=1:size(B, 2)
            @inbounds for colindex = 1:M.n
                m2 = alpha*B[colindex, j]
                @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                    row = M.rowval[i]
                    result[row, j] += M.nzval[i]*m2
                end
            end
        end
    end
end

function gemm!(alpha::ComplexF64, M::SparseMatrix, B::Matrix{ComplexF64}, beta::ComplexF64, result::Matrix{ComplexF64})
    if beta == ComplexF64(0.)
        fill!(result, beta)
    elseif beta != ComplexF64(1.)
        rmul!(result, beta)
    end
    if nnz(M) > 1000
        gemm_sp_dense_big(alpha, M, B, result)
    else
        gemm_sp_dense_small(alpha, M, B, result)
    end
end

function gemm!(alpha::ComplexF64, B::Matrix{ComplexF64}, M::SparseMatrix, beta::ComplexF64, result::Matrix{ComplexF64})
    if beta == ComplexF64(0.)
        fill!(result, beta)
    elseif beta != ComplexF64(1.)
        rmul!(result, beta)
    end
    dimB = size(result,1)
    if alpha == ComplexF64(1.)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = M.nzval[i]
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[j, colindex] += mi*B[j, mrowvali]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = M.nzval[i]*alpha
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[j, colindex] += mi*B[j, mrowvali]
                end
            end
        end
    end
end

function gemv!(alpha::ComplexF64, M::SparseMatrix, v::Vector{ComplexF64}, beta::ComplexF64, result::Vector{ComplexF64})
    if beta == ComplexF64(0.)
        fill!(result, beta)
    elseif beta != ComplexF64(1.)
        rmul!(result, beta)
    end
    if alpha == ComplexF64(1.)
        @inbounds for colindex = 1:M.n
            vj = v[colindex]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i]] += M.nzval[i]*vj
            end
        end
    else
        @inbounds for colindex = 1:M.n
            vj = alpha*v[colindex]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i]] += M.nzval[i]*vj
            end
        end
    end
end

function gemv!(alpha::ComplexF64, v::Vector{ComplexF64}, M::SparseMatrix, beta::ComplexF64, result::Vector{ComplexF64})
    if beta == ComplexF64(0.)
        fill!(result, beta)
    elseif beta != ComplexF64(1.)
        rmul!(result, beta)
    end
    if alpha == ComplexF64(1.)
        @inbounds for colindex=1:M.n
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[colindex] += M.nzval[i]*v[M.rowval[i]]
            end
        end
    else
        @inbounds for colindex=1:M.n
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[colindex] += M.nzval[i]*alpha*v[M.rowval[i]]
            end
        end
    end
end

function sub2sub(shape1::NTuple{N, Int}, shape2::NTuple{M, Int}, I::CartesianIndex{N}) where {N, M}
    linearindex = LinearIndices(shape1)[I.I...]
    CartesianIndices(shape2)[linearindex]
end

function ptrace(x, shape_nd::Vector{Int}, indices::Vector{Int})
    shape_nd = (shape_nd...,)
    N = div(length(shape_nd), 2)
    shape_2d = (x.m, x.n)
    shape_nd_after = ([i ∈ indices || i-N ∈ indices ? 1 : shape_nd[i] for i=1:2*N]...,)
    shape_2d_after = (prod(shape_nd_after[1:N]), prod(shape_nd_after[N+1:end]))
    I_nd_after_max = CartesianIndex(shape_nd_after...)
    y = spzeros(ComplexF64, shape_2d_after...)
    for I in eachindex(x)
        I_nd = sub2sub(shape_2d, shape_nd, I)
        if I_nd.I[indices] != I_nd.I[indices .+ N]
            continue
        end
        I_after = sub2sub(shape_nd_after, shape_2d_after, min(I_nd, I_nd_after_max))
        y[I_after] += x[I]
    end
    y
end

function permutedims(x, shape, perm)
    shape = (shape...,)
    shape_perm = ([shape[i] for i in perm]...,)
    y = spzeros(ComplexF64, x.m, x.n)
    for I in eachindex(x)
        linear_index = LinearIndices((x.m, x.n))[I.I...]
        cartesian_index = CartesianIndices(shape)[linear_index]
        cartesian_index_perm = [cartesian_index[i] for i=perm]
        linear_index_perm = LinearIndices(shape_perm)[cartesian_index_perm...]
        J = Tuple(CartesianIndices((x.m, x.n))[linear_index_perm])
        y[J...] = x[I.I...]
    end
    y
end


text:
using Test

include("routines/gemm.jl")


text:
@testset "StudentT" begin
    test_interface(
        StudentTLikelihood(3.0, 1.5),
        Distributions.AffineDistribution{Float64,Continuous,<:TDist},
    )
    test_auglik(StudentTLikelihood(3.0, 1.5))
end


text:
# DSMC

using StatsBase: samplepair, percentile
using LinearAlgebra: norm
using Random: MersenneTwister
using Printf: @sprintf
using ProgressMeter: Progress, update!, finish!
using Dates: now

# Storage of dynamic simulation data
mutable struct CloudBuffer
    positions :: Matrix{Float64}
    velocities :: Matrix{Float64}

    accels :: Matrix{Float64}
    potens :: Vector{Float64}
    survival_prob :: Vector{Float64}

    assignments :: Vector{Int64}
    cell_offsets :: Vector{Int64}
    cell_occupancies :: Vector{Int64}
    atom_lookup :: Vector{Int64}

    Nt :: Int64
    F :: Float64

    cellcount :: Int64
    nonempty_count :: Int64
    cellsize :: Vector{Float64}

    function CloudBuffer(conditions)
        # Atom buffers
        pos = copy(conditions.positions)
        vel = copy(conditions.velocities)
        Nt = size(pos, 2)
        acc = zeros(Float64, 3, Nt)
        pot = zeros(Float64, Nt)
        prob = zeros(Float64, Nt)
        assign = zeros(Int64, Nt)
        atoms = zeros(Int64, Nt)

        # Cell buffers (overwritten on first iteration)
        offsets = zeros(Int64, 1)
        occupancies = zeros(Int64, 1)

        # Dynamic cell info (updated on every iteration)
        cellsize = [0.0, 0.0, 0.0]
        nonempty_count = 1
        cellcount = 1

        # Return cloud
        return new(pos, vel, acc, pot, prob, assign, offsets, occupancies, atoms,
                   Nt, conditions.F, cellcount, nonempty_count, cellsize)
    end
end

# Velocty Verlet step + calculate motion-based limit on next timestep
function verlet_step!(cloud, conditions, t, dt)
    Nt = cloud.Nt
    positions = view(cloud.positions, :, 1:Nt)
    velocities = view(cloud.velocities, :, 1:Nt)
    accels = view(cloud.accels, :, 1:Nt)
    # Half timestep
    half_dt = dt / 2.0
    # Approximate position half-way through timestep
    for i in eachindex(positions) # Elementwise iteration
        positions[i] += half_dt * velocities[i]
    end
    # Acceleration halfway through timestep (modifies cloud.accels in-place)
    conditions.acceleration(positions, conditions.species, t + half_dt, accels)
    # Compute final velocities and positions, as well as new timestep
    for atom in 1:Nt
        ax, ay, az = view(accels, :, atom)
        velocities[1, atom] += dt * ax # Update velocity
        velocities[2, atom] += dt * ay
        velocities[3, atom] += dt * az
        vx, vy, vz = view(velocities, :, atom)
        positions[1, atom] += half_dt * vx # Update position
        positions[2, atom] += half_dt * vy
        positions[3, atom] += half_dt * vz
    end
end

# Assign particles to cells
function assign_cells!(cloud, peak_free_path, Nc, Vprev)
    N = cloud.Nt
    positions = view(cloud.positions, :, 1:N)

    # Corners of the bounding box around all atoms
    minx, miny, minz = Inf, Inf, Inf
    maxx, maxy, maxz = -Inf, -Inf, -Inf
    for atom in 1:N
        x, y, z = view(positions, :, atom)
        if x < minx; minx = x; elseif x > maxx; maxx = x; end;
        if y < miny; miny = y; elseif y > maxy; maxy = y; end;
        if z < minz; minz = z; elseif z > maxz; maxz = z; end;
    end

    minpos = [minx, miny, minz]
    maxpos = [maxx, maxy, maxz]

    # Cell dimensions
    V = min(Vprev, prod(maxpos - minpos)) # Volume estimate
    Vc = V * Nc / N
    ds = min(cbrt(Vc), 0.05 * peak_free_path)
    cellsize = [ds, ds, ds]
    cloud.cellsize = cellsize

    # Extend grid slightly beyond the furthest atoms, to
    # avoid edge atoms being placed outside the grid.
    lowerpos = minpos - (0.25 * cellsize)
    upperpos = maxpos + (0.25 * cellsize)
    gridshape = ceil.(Int64, (upperpos - lowerpos) ./ cellsize)
    cellcount = prod(gridshape)

    # Cell occupancies
    occupancy_dict = Dict{Int64, Int64}()

    # Assign atoms to cells and count the number of atoms in each cell.
    assignments = cloud.assignments
    xcount, ycount, _ = gridshape
    lx, ly, lz = lowerpos # Cache lower coordinates
    for atom in 1:N
        x = ceil(Int64, (positions[1, atom] - lx) / ds)
        y = ceil(Int64, (positions[2, atom] - ly) / ds)
        z = ceil(Int64, (positions[3, atom] - lz) / ds)
        cell = x + xcount * ( (y - 1) + ycount * (z - 1) )
        if haskey(occupancy_dict, cell)
            occupancy_dict[cell] += 1
        else
            occupancy_dict[cell] = 1
        end
        assignments[atom] = cell
    end

    cellcount = length(occupancy_dict) # ONLY NONEMPTY CELLS

    # Ensure storage buffers are large enough.
    if cellcount > length(cloud.cell_offsets)
        cloud.cell_offsets = zeros(Int64, cellcount * 2)
        cloud.cell_occupancies = zeros(Int64, cellcount * 2)
    end
    cell_offsets = view(cloud.cell_offsets, 1:cellcount)
    cell_occupancies = view(cloud.cell_occupancies, 1:cellcount)
    fill!(cell_offsets, 0)
    fill!(cell_occupancies, 0)

    # Map cell numbers to nonempty cell index, and
    # compute storage offset for each cell
    index_dict = Dict{Int64, Int64}()
    i = 1
    offset = 1
    for (cell, Nc) in occupancy_dict
        if Nc > 0
            index_dict[cell] = i
            cell_offsets[i] = offset
            cell_occupancies[i] = Nc
            offset += Nc
            i += 1
        end
    end

    # Create lookup index of atoms in each cell.
    # This temporarily modifies the cell offsets.
    atom_lookup = cloud.atom_lookup
    for atom in 1:N
        cell = assignments[atom]
        cell_index = index_dict[cell]
        atom_lookup[ cell_offsets[cell_index] ] = atom
        cell_offsets[cell_index] += 1
    end

    # Restore cell offsets
    for cell in 1:cellcount
        Nc = cell_occupancies[cell]
        cell_offsets[cell] -= Nc
    end

    cloud.nonempty_count = cellcount
    cloud.cellcount = cellcount

    # Compute peak density.
    peak_density = maximum(cell_occupancies) / Vc
    #percentile(cell_occupancies, 97) / Vc

    Vnew = cellcount * Vc # Total volume of non-empty cells

    return peak_density, gridshape, Vnew
end

# Simulate collisions
function collision_step!(cloud, dt, σ, rng=MersenneTwister())
    cellcount = cloud.nonempty_count
    Ncs = cloud.cell_occupancies
    offsets = cloud.cell_offsets
    velocities = cloud.velocities
    atom_lookup = cloud.atom_lookup
    F = cloud.F

    # Check appropriate number of pairs in each cell
    Vc = prod(cloud.cellsize) # Cell volume
    tot_cand, tot_coll = 0, 0
    for cell in 1:cellcount
        Nc = Ncs[cell]
        (Nc < 2) && continue # Skip cells with 1 particle
        offset = offsets[cell]
        # Compute maximum speed
        speeds = zeros(Float64, Nc)
        for i in 1 : Nc
            atom = atom_lookup[i + offset - 1]
            vx, vy, vz = view(velocities, :, atom)
            speeds[i] = sqrt(vx^2 + vy^2 + vz^2)
        end
        max_speed = maximum(speeds)
        # Select appropriate number of pairs
        Mraw = F * (dt * σ / Vc) * Nc * (Nc - 1) * max_speed
        # Note: |v_rel|max = 2 * max_speed, so expression slightly different
        # to theory in report.
        Mcand = ceil(Int64, Mraw)
        Mcoll = 0 # Total collisions in cell
        prob_adjust = Mraw / Mcand # Adjustment for rounding
        Pcoeff = prob_adjust / (2 * max_speed)
        # Check and perform collisions
        for _ in 1:Mcand
            i1, i2 = samplepair(rng, Nc)
            atom1 = atom_lookup[i1 + offset - 1]
            atom2 = atom_lookup[i2 + offset - 1]
            # Get current velocities
            u1x, u1y, u1z = view(velocities, :, atom1)
            u2x, u2y, u2z = view(velocities, :, atom2)
            urel = sqrt((u1x - u2x)^2 + (u1y - u2y)^2 + (u1z - u2z)^2)
            Pij = urel * Pcoeff
            if rand(rng) < Pij # Collision prob.
                # Compute new relative velocity
                ϕ = 2 * π * rand(rng)
                cosθ = 2 * rand(rng) - 1
                sinθ = sqrt(1 - cosθ^2)
                vrx = urel * sinθ * cos(ϕ)
                vry = urel * sinθ * sin(ϕ)
                vrz = urel * cosθ
                # Update stored velocities
                vcmx = 0.5 * (u1x + u2x)
                vcmy = 0.5 * (u1y + u2y)
                vcmz = 0.5 * (u1z + u2z)
                v1x = vcmx + 0.5 * vrx;
                v1y = vcmy + 0.5 * vry;
                v1z = vcmz + 0.5 * vrz;
                v2x = vcmx - 0.5 * vrx;
                v2y = vcmy - 0.5 * vry;
                v2z = vcmz - 0.5 * vrz;
                velocities[1, atom1] = v1x
                velocities[2, atom1] = v1y
                velocities[3, atom1] = v1z
                velocities[1, atom2] = v2x
                velocities[2, atom2] = v2y
                velocities[3, atom2] = v2z
                # Update cell's collision count and maximum speed
                speed1 = sqrt(v1x^2+ v1y^2+v1z^2)
                speed2 = sqrt(v2x^2+ v2y^2+ v2z^2)
                speeds[i1] = speed1
                speeds[i2] = speed2
                max_speed = max(max_speed, speed1, speed2) #maximum(speeds)
                Pcoeff = prob_adjust / (2 * max_speed)
                Mcoll += 1
            end
        end
        tot_cand += Mcand
        tot_coll += Mcoll
    end

    # Calculate collision-based timestep limit
    Γₐ = 2 * tot_coll / (cloud.Nt * dt) # Mean collision rate per atom
    dt = 0.05/Γₐ

    return dt, tot_cand, tot_coll
end

# Atom loss effects: high-energy particles, three-body recombination &
# background collisions.
function atom_loss!(cloud, conditions, t, dt, rng=MersenneTwister())
    τ_bg, K, evap = conditions.τ_bg, conditions.threebodyloss, conditions.evaporate

    N₀, F = cloud.Nt, cloud.F
    dV = prod(cloud.cellsize)
    positions = view(cloud.positions, :, 1:N₀)
    velocities = view(cloud.velocities, :, 1:N₀)

    # Evaporation LOSS probability for each atom (becomes survival later)
    p_survive = view(cloud.survival_prob, 1:N₀)
    evap(positions, velocities, conditions, t, p_survive)

    p_background = 1 - (dt / τ_bg) # Background survival probability (const.)
    for cell in 1:cloud.cellcount
        Nc = cloud.cell_occupancies[cell]
        offset = cloud.cell_offsets[cell]
        # Three-body-loss survival probability (constant within cell)
        p_threebody = 1 - (K * F^2 * Nc * (Nc - 1) / (dV^2) * dt)
        for i in offset : offset + Nc - 1
            # Update survival likelihood to reflect 3-body + background loss
            atom = cloud.atom_lookup[i]
            p_survive[atom] = (1.0 - p_survive[atom]) * p_background * p_threebody
        end
    end

    # Perform atom loss
    N = N₀
    for atom in N₀:-1:1
        if rand(rng) > p_survive[atom]
            # Remove atom by replacing it with atom from the end
            velocities[:,atom] .= view(velocities, :, N)
            positions[:,atom] .= view(positions, :, N)
            N -= 1
        end
    end
    cloud.Nt = N
    return nothing
end

# Repopulate cloud by creating a duplicate of each
# particle, with the opposite velocity.
function duplicate!(cloud)
    N₀ = cloud.Nt
    N₁ = 2 * N₀
    positions = view(cloud.positions, :, 1:N₀)
    velocities = view(cloud.velocities, :, 1:N₀)
    if N₁ > size(cloud.positions, 2) # Buffer requires resizing
        @warn "Cloud duplication triggered above half capacity."
        cloud.positions = hcat(positions, positions)
        cloud.positions = hcat(velocities, velocities)
        # NOTE: switched from reflected velocities to identical velocities.
        accels = view(cloud.accels, :, 1:N₀)
        cloud.accels = hcat(accels, accels)
    else # No resize required
        cloud.positions[:, N₀+1 : N₁] .= positions
        cloud.velocities[:, N₀+1 : N₁] .= velocities
    end

    cloud.Nt = N₁
    cloud.F /= 2

    return nothing
end

# Evolve initial conditions for desired duration
function evolve(conditions, duration; Nc = 1, max_dt = Inf,
    rng=MersenneTwister(), measure = null, meas_dt = 1e-4)
    #= Arguments
    - conditions: a SimulationConditions
    - duration: in virtual time
    - Nc: target average number of atoms per cell.
    - max_dt: upper bound on timestep.
    - rng: allows reproducibility if desired.
    =#

    # Initialise dynamic storage
    cloud = CloudBuffer(conditions)

    σ = conditions.species.σ
    max_dt = time_parametrize(max_dt)

    # Initial values
    peak_density, _, V = assign_cells!(cloud, 1e-5, Nc, Inf)
    dt = min(max_dt(0), 1e-9) # Begin with very short timestep

    # Track progress
    prog_detail = 10000
    progress = Progress(prog_detail, dt = 1, desc = "Simulation progress: ",
                color = :green, showspeed = false, enabled = true, barlen=25)

    accum_cand = 0 # Accumulated since previous measurement
    accum_coll = 0
    accum_n0 = 0
    accum_V = 0
    meas_iter = 0
    meas_t = -Inf # Previous measurement time

    # Iterate simulation
    t = 0 # Virtual time
    iter_count = 0
    start_time = now()
    while t < duration
        # Collisionless motion
        verlet_step!(cloud, conditions, t, dt)

        # Sort atoms by cell
        peak_free_path = 1 / (4 * peak_density * cloud.F * σ);
        peak_density, gridshape, V = assign_cells!(cloud, peak_free_path, Nc, V)

        # Perform collisions
        coll_dt, cand_count, coll_count = collision_step!(cloud, dt, σ, rng)

        atom_loss!(cloud, conditions, t, dt, rng)

        if cloud.Nt == 0
            @warn "Trap empty"
            return cloud
        end

        while cloud.Nt < size(cloud.positions, 2) / 2
            duplicate!(cloud)
        end

        # Increment time and then update timestep
        t += dt
        dt = min(coll_dt, max_dt(t))

        # Record for measurements
        accum_cand += cand_count
        accum_coll += coll_count
        accum_n0 += peak_density
        accum_V += V

        iter_count += 1

        # External measurements on system at specified intervals
        if t >= meas_t + meas_dt
            elapsed_iter = iter_count - meas_iter
            n0 = accum_n0 / elapsed_iter
            V_avg = accum_V / elapsed_iter
            measure(cloud, conditions, accum_cand, accum_coll, t, n0, V_avg)
            accum_cand, accum_coll = 0, 0
            accum_n0, accum_V = 0, 0
            meas_t = t
            meas_iter = iter_count
        end

        #Update progress
        if (mod(iter_count, 100) == 0) || (t >= duration)
            progressvalues = [
                ("Virtual time", @sprintf("%.3g / %.3g s", t, duration)),
                ("Real time", now() - start_time),
                ("Iterations", iter_count),
                ("Coll. dt", @sprintf("%.3g", coll_dt)),
                ("Max. dt", @sprintf("%.3g", max_dt(t))),
                ("dt", @sprintf("%.3g", dt)),
                ("Cell size", @sprintf("%.3g x %.3g x %.3g", cloud.cellsize...)),
                ("Grid shape", @sprintf("%i x %i x %i", gridshape...)),
                ("Non-empty cell count", cloud.nonempty_count),
                ("Peak density", cloud.F * peak_density),
                ("Np", cloud.Nt * cloud.F),
                ("Nt", cloud.Nt)
            ]
            if t >= duration
                finish!(progress, showvalues = progressvalues)
            else
                update!(progress, ceil(Int64, prog_detail * t / duration),
                        showvalues = progressvalues)
            end
        end
    end

    # Return final positions and velocities
    return cloud
end


text:
#prec = 15 #precision

##Single Particle
eigensp(s::State, prec=15) = sort([round(eigvals(Matrix(rhosp(s)))[i], digits = prec) for i in 1:dim(ope(s))], rev=true)
eigensp(s::State_fixed, prec=15) = sort([round(eigvals(Matrix(rhosp(s)))[i], digits = prec) for i in 1:dim(ope(s))], rev=true)
# I can not compute eigenvalues
#directly from sparse matrices

function ssp(sta::State, prec=15)
    eigen = eigensp(sta, prec)
    lene = length(eigen)
    s = 0
    for i in 1:lene
        if eigen[i] != 0 && eigen[i] != 1
            s = s - (eigen[i]*log(2,eigen[i]) + (1 - eigen[i])*log(2,1-eigen[i]))
        end
    end
    return s/lene
end

function ssp(sta::State_fixed, prec=15)
    eigen = eigensp(sta, prec)
    lene = length(eigen)
    s = 0
    for i in 1:lene
        if eigen[i] != 0 && eigen[i] != 1
            s = s - (eigen[i]*log(2,eigen[i]) + (1 - eigen[i])*log(2,1-eigen[i]))
        end
    end
    return s/lene
end
## Quasi Particles
eigenqsp(s::State, prec=15) = sort([round(eigvals(Matrix(rhoqsp(s)))[i], digits = prec) for i in 1:(2*dim(ope(s)))], rev=true)
# I can not compute eigenvalues
#directly from sparse matrices

function sqsp(sta::State, prec=15)
    eigen = eigenqsp(sta, prec)
    lene = length(eigen)
    s = 0
    for i in 1:lene
        if eigen[i] != 0 && eigen[i] != 1
            s = s - (eigen[i]*log(2,eigen[i]) + (1 - eigen[i])*log(2,1-eigen[i]))
        end
    end
    return s/lene
end

#This function outputs 1 if s1 majorizes s2,
#-1 if s2 majorizes s1,
#0 if there is no majorization relation
function majorization_sp(s1, s2)
    e1 = eigensp(s1)
    e2 = eigensp(s2)
    count1 = 0
    count2 = 0
    d = length(e1)
    for j in 1:(d-1) #el último valor a veces hace problemas
        if sum([round(e1[i]-e2[i], digits = 10) for i in 1:j]) < 0
            count1 = count1 + 1
        elseif sum([round(e2[i]-e1[i], digits = 10) for i in 1:j]) < 0
            count2 = count2 + 1
        end
    end
    if count1 == (d-1)
        #println(@Name(s2)," sp majorizes ", @Name(s1))
        return -1
    elseif count2 == (d-1)
        #println(@Name(s1)," sp majorizes ", @Name(s2))
        return 1
    else
        #println("There is no sp majorization relation")
        return 0
    end
end

function majorization_qsp(s1, s2)
    e1 = eigenqsp(s1)
    e2 = eigenqsp(s2)
    count1 = 0
    count2 = 0
    d = length(e1)
    for j in 1:(d-1) #el último valor a veces hace problemas
        if sum([round(e1[i]-e2[i], digits = 10) for i in 1:j]) < 0
            count1 = count1 + 1
        elseif sum([round(e2[i]-e1[i], digits = 10) for i in 1:j]) < 0
            count2 = count2 + 1
        end
    end
    if count1 == (d-1)
        #println(@Name(s2)," qsp majorizes ", @Name(s1))
        return -1
    elseif count2 == (d-1)
        #println(@Name(s1)," qsp majorizes ", @Name(s2))
        return 1
    else
        #println("There is no qsp majorization relation")
        return 0
    end
end

function n_avg(s::State)
    navg = Real(tr(rhosp(s)))
    return navg
end

####################################################
# Following are the matrices of coefficients in
# the m-(N-m) schmidt like decomposition
function coef(s::State, m::Int64, prec=15)
    d = dim(s)
    num = Int(round(n_avg(s),digits = 8))
    bas, _ = basis_m(d, m)
    estat = st(s)
    ty = typ(s)
    cof = cof_mat(d, num, m, bas, estat, ty, prec)
    return round.(cof,digits=prec)
end

function coef(s::State_fixed, m::Int64, prec=15)
    d = dim(s)
    num = nume(s)
    bas_tot = basis(s)
    bas, _ = basis_m(d, m)
    estat = st(s)
    ty = typ(s)
    cof = cof_mat_fixed(d, num, m, bas, bas_tot, estat, ty, prec)
    return round.(cof,digits=prec)
end

#The following are the rhom matrices without diagonalization
# i.e. in the original basis
# For general m with fixed_state
function rhom(s::Union{State,State_fixed}, m::Int64, prec=15)
    cof = coef(s, m, prec)
    rhom = cof*cof'
    return round.(rhom,digits=prec)
end

#The following are the rhom matrices with diagonalization
function rhomd(s::Union{State,State_fixed}, m::Int64, prec=15)
        r = rhom(s,m,prec)
        dim = size(r)[1]
        d = non_zero(eigvals(Matrix(r)))[2]
        row = [i for i in 1:length(d)]
        col = row
        rhomd = sparse(row, col, d, dim, dim)
        return round.(rhomd,digits=prec)
end

#A more efficient 2 body DM in fixed space
function rhom2(s::State_fixed, prec=15)
    m = 2
    o = ope(s)
    num = nume(s)
    d = dim(o)
    rhomnd = spzeros(typ(s),binomial(d, m), binomial(d, m))
    bas,_ = basis_m(d,m)
    estat = st(s)
    for i in 1:binomial(d, m)
        indi = indx(bas[i,:])
        for j in 1:binomial(d,m)
            indj = indx(bas[j,:])
            #rhomd[i,j] = round(estat'*fixed(and[j]*and[i]',num)*estat, digits = 14)
            if indi[2] != indj[1] && indi[1] != indj[2]
                rhomnd[i,j] = round(estat'*ada(o,indj[1],indi[1])*ada(o,indj[2],indi[2])*estat, digits = prec)
            else
                rhomnd[i,j] = round(-estat'*ada(o,indj[1],indi[2])*ada(o,indj[2],indi[1])*estat, digits = prec)
            end
        end
    end
    return rhomnd
end

# Following function is used for constructing coef and rhom
#for example, for a vector [0,1,0,0,1]->[2,5]
function indx(arr)
    l = length(arr)
    ind = spzeros(Int,0)
    for i=1:l
        if arr[i] != 0
            ind = sparse([ind; i])
        end
    end
    return ind
end

# Core functons of the coeficient matrices
function cof_mat(d, num, m, bas, estat, ty, prec=15)
    cmat = spzeros(ty, binomial(d,m),binomial(d,num-m))
    am = sparse([bas[i,:] for i in 1:binomial(d, m)])
    bas_nume_m, _ = basis_m(d, num - m)
    amn = sparse([bas_nume_m[i,:] for i in 1:binomial(d, num-m)])
    for i in 1:binomial(d, m)
        ai = am[i]
        indi = indx(ai)
        for j in 1:binomial(d, num-m)
            aj = amn[j]
            if ai'*aj != 0
                cmat[i,j] = 0
            else
                permut = 0
                for k in 1:length(indi)
                    permut = permut + sum(aj[1:floor(Int,indi[k]-1)])
                end
                indij = indx(ai+aj)
                elem = Int(sum([2^(d-i) for i in indij]) + 1)
                cmat[i,j] = round(estat[elem] * (-1)^permut, digits=prec)
            end
        end
    end
    return cmat
end

function cof_mat_fixed(d, num, m, bas, bas_tot, estat, ty, prec=15)
    cmat = spzeros(ty,binomial(d,m),binomial(d,num-m))
    am = sparse([bas[i,:] for i in 1:binomial(d, m)])
    bas_nume_m, _ = basis_m(d, num - m)
    amn = sparse([bas_nume_m[i,:] for i in 1:binomial(d, num-m)])
    for i in 1:binomial(d, m)
        ai = am[i]
        indi = indx(ai)
        for j in 1:binomial(d, num-m)
            aj = amn[j]
            if ai'*aj != 0
                cmat[i,j] = 0
            else
                permut = 0
                for k in 1:length(indi)
                    permut = permut + sum(aj[1:floor(Int,indi[k]-1)])
                end
                indij = indx(ai+aj)
                elem = 1
                while indx(bas_tot[elem,:]) != indij
                    elem = elem + 1
                end
                cmat[i,j] = round(estat[elem] * (-1)^permut, digits=prec)
            end
        end
    end
    return cmat
end

#Partial trace in a given basis of modes
function trp(state::State, modos::Array{Int64,1})
    d = dim(ope(state))
    bas = basis(ope(state))
    sta = st(state)
    lm = length(modos)
    zvr = zeros(typ(state),2^lm,2^(d-lm))
    full = [i for i in 1:d]
    lista = sort(modos)
    listb = filter(x->x ∉ lista,full)
    for k in 1:2^d
        indicea = parse(Int,join([Int(bas[k,i]) for i in lista]), base=2)+1
        indiceb = parse(Int,join([Int(bas[k,i]) for i in listb]), base=2)+1
        signi = 0
        for i in lista
            listac = filter(x->x <= i-1,listb)
            signl = 0
            for l in listac
                signl = signl+bas[k,l]
            end
            signi = signi + signl*bas[k,i]
        end
        sign = (-1)^(signi)
        zvr[indicea,indiceb] = sta[k]*sign
    end
    rhoa = zvr*zvr'
    #rhob=zvr'*zvr;
    return rhoa
end

#= This was used for finding the diagonalizing
# operators in m
function non_diag_ops(o, d, nume, m, bas, i)
    vec = indx(bas[i, :])
    #le es nume-1 me parece
    mel = 1
    for j in 1:m
        mel = mel*ad(o, Int(vec[j]))
    end
    return mel
end

function diag_ops(o, d, nume, m, bas, u, and, i)
    adi = spzeros(2^d, 2^d)
    for j in 1:binomial(d,m)
        if round(u[j,i], digits=14) != 0
            adi = adi + u[j, i]*and[j]
        end
    end
    return adi
end
=#

#= Work in Progress: non diagonal operators with fixed states
function rhomnd(s::Union{State_fixed,State_sparse_fixed}, m::Int64)
    o = ope(s)
    num = nume(s)
    d = dim(o) #tengo que pasar a matrix lamentablemente
    rhomd = spzeros(binomial(d, m), binomial(d, m))
    bas, _ = basis_m(d, m)
    estat = st(s)
    and = sparse([non_diag_ops(o, d, num, m, bas, j) for j in 1:binomial(d, m)])
    for i in 1:binomial(d, m)
        for j in 1:binomial(d,m)
            rhomd[i,j] = round(estat'*and[j]*and[i]'*estat, digits = 14)
        end
    end
    return rhomd
end

#main function
function rhomnd(s::Union{State_complex_fixed,State_sparse_complex_fixed}, m::Int64)
    o = ope(s)
    num = nume(s)
    d = dim(o) #tengo que pasar a matrix lamentablemente
    rhomd = spzeros(Complex{Float64},binomial(d, m), binomial(d, m))
    bas, _ = basis_m(d, m)
    estat = st(s)
    and = sparse([non_diag_ops(o, d, num, m, bas, j) for j in 1:binomial(d, m)])
    for i in 1:binomial(d, m)
        for j in 1:binomial(d,m)
            rhomd[i,j] = round(estat'*and[j]*and[i]'*estat, digits = 14)
        end
    end
    return rhomd
end
=#

#= Work in Progress: partial trace for fixed particle number
function trp(state::Union{State_fixed,State_sparse_fixed},modos::Array{Int64,1})
    d = dim(ope(state))
    bas = basis(ope(state))
    sta = st(state)
    lm = length(modos)
    zvr = zeros(Complex{Float64},2^lm,2^(d-lm))
    full = [i for i in 1:d]
    lista = sort(modos)
    listb = filter(x->x ∉ lista,full)
    for k in 1:2^d
        indicea=parse(Int,join([Int(bas[k,i]) for i in lista]), base=2)+1
        indiceb=parse(Int,join([Int(bas[k,i]) for i in listb]), base=2)+1
        signi = 0
        for i in lista
            listac = filter(x->x <= i-1,listb)
            signl = 0
            for l in listac
                signl = signl+bas[k,l]
            end
            signi = signi + signl*bas[k,i]
        end
        sign = (-1)^(signi)
        zvr[indicea,indiceb]=sta[k]*sign
    end
    rhoa=zvr*conj(zvr')
    #rhob=zvr'*zvr;
    return rhoa
end

function trp(state::Union{State_complex_fixed,State_sparse_complex_fixed},modos::Array{Int64,1})
    d = dim(ope(state))
    bas = basis(ope(state))
    sta = st(state)
    lm = length(modos)
    zvr = zeros(Complex{Float64},2^lm,2^(d-lm))
    full = [i for i in 1:d]
    lista = sort(modos)
    listb = filter(x->x ∉ lista,full)
    for k in 1:2^d
        indicea=parse(Int,join([Int(bas[k,i]) for i in lista]), base=2)+1
        indiceb=parse(Int,join([Int(bas[k,i]) for i in listb]), base=2)+1
        signi = 0
        for i in lista
            listac = filter(x->x <= i-1,listb)
            signl = 0
            for l in listac
                signl = signl+bas[k,l]
            end
            signi = signi + signl*bas[k,i]
        end
        sign = (-1)^(signi)
        zvr[indicea,indiceb]=sta[k]*sign
    end
    rhoa=zvr*zvr'
    #rhob=zvr'*zvr;
    return rhoa
end
=#

#= This code is faster for many iterations (majorization)
count = 0
@time for i in 1:repeticiones
    #definimos estado premedida
    state_ran = spzeros(Complex{Float64},l)
    for i in 1:l
        if sum(basis(o)[i,:]) == nume
            state_ran[i] = 2*rand(Complex{Float64},1)[1]-1-im
        end
    end
    state_ran = state_ran/sqrt(state_ran'*state_ran)
    state_ran = State_sparse_complex(state_ran, o);
    #e = sort(eigvals(Matrix(rhoqsp(state_ran))), rev=true)

    #ahora medimos
    state_ran_medido = operator*st(state_ran);
    state_ran_medido = state_ran_medido/sqrt(state_ran_medido'*state_ran_medido)
    state_ran_medido = State_sparse_complex(state_ran_medido, o)
    #e_medido = sort(eigvals(Matrix(rhoqsp(state_ran_medido))), rev=true)

    #y chequeamos mayorización
    global count = count + majorization_qsp(state_ran_medido, state_ran)
end

if count == repeticiones
    println("Measured state qsp majorizes unmeasured state")
elseif count == -repeticiones
    println("Unmeasured state qsp majorizes measured state")
else
    println("There is no qsp majorization")
end

=#

#This macro serves for printing as a string the variable's name
#=
macro Name(arg)
   string(arg)
end
=#

#= function rhom(s::State, m::Int64, prec=15)
#     o = ope(s)
#     num = Int(round(n_avg(s),digits = 8))
#     d = dim(o) #tengo que pasar a matrix lamentablemente
#     rhomd = spzeros(typ(s),binomial(d, m), binomial(d, m))
#     bas,_ = basis_m(d, m)
#     estat = st(s)
#     and = sparse([non_diag_ops(o, d, num, m, bas, j) for j in 1:binomial(d, m)])
#     for i in 1:binomial(d, m)
#         for j in 1:binomial(d,m)
#             rhomd[i,j] = round(estat'*and[j]*and[i]'*estat, digits = prec)
#         end
#     end
#     return rhomd
# # end
#
# #main function
# function rhomd(s::State, m::Int64, prec=15)
#     o = ope(s)
#     num = Int(round(n_avg(s),digits = 8))
#     d = dim(o) #tengo que pasar a matrix lamentablemente
#     rhomd = spzeros(typ(s),binomial(d, m), binomial(d, m))
#     bas, _ = basis_m(d, m)
#     estat = st(s)
#     ty = typ(s)
#     u, _ = svd(cof_mat(d, num, m, bas, estat, ty))
#     and = sparse([non_diag_ops(o, d, num, m, bas, j) for j in 1:binomial(d, m)])
#     adi = sparse([diag_ops(o, d, num, m, bas, u, and, j) for j in 1:binomial(d, m)])
#     #adi[1] es evaluar la funcion diagonal en i=1
#     for i in 1:binomial(d, m)
#         rhomd[i,i] = round(estat'*adi[i]*adi[i]'*estat, digits = prec)
#     end
#     return rhomd
# end

método directo (muy lento)
function rhom2(s::State)
    d = dim(ope(s))
    o = ope(s)
    bas,_ = basis_m(d,2)
    l = Int(binomial(d,2))
    rhom = spzeros(l,l);
    stat = st(s)
    for i in 1:l
        indi = indx(bas[i,:])
        for j in 1:l
            indj = indx(bas[j,:])
                rhom[i,j]=(stat'*ad(o,Int(indi[2]))*ad(o,Int(indi[1]))*a(o,Int(indj[1]))*a(o,Int(indj[2]))*stat)
        end
    end
    return rhom
end

rhom(State(1/sqrt(2)*(ad(o,1)*ad(o,2)*ad(o,3)*ad(o,4)+ad(o,3)*ad(o,4)*ad(o,5)*ad(o,6))*vacuum(o),o),2)

# It works, but it is not efficient,
# we need to optimieze ir for state_fixed
function rhom(s::State_fixed, o::Op, m::Int64, prec=15)
    d = dim(o)
    #o = Op(d)
    num = nume(s)
    rhomd = spzeros(typ(s),binomial(d, m), binomial(d, m))
    bas,_ = basis_m(d, m)
    estat = st(s)
    and = sparse([non_diag_ops(o, d, num, m, bas, j) for j in 1:binomial(d, m)])
    for i in 1:binomial(d, m)
        for j in 1:binomial(d,m)
            rhomd[i,j] = round(estat'*fixed(and[j]*and[i]',num)*estat, digits = prec)
        end
    end
    return rhomd
end

#
# # It works, but it is not efficient,
# # we need to optimieze ir for state_fixed
# function rhomd(s::State_fixed, o::Op, m::Int64, prec=15)
#     d = dim(o)
#     #o = Op(d)
#     num = nume(s)
#     rhomd = spzeros(typ(s),binomial(d, m), binomial(d, m))
#     bas_tot, _ = basis_m(d, num)
#     bas, _ = basis_m(d, m)
#     estat = st(s)
#     ty = typ(s)
#     u, _ = svd(cof_mat_fixed(d, num, m, bas, bas_tot, estat, ty))
#     and = sparse([non_diag_ops(o, d, num, m, bas, i) for i in 1:binomial(d, m)])
#     adi = sparse([diag_ops(o, d, num, m, bas, u, and, i) for i in 1:binomial(d, m)])
#     #ad[1] es evaluar la funcion diagonal en i=1
#     for i in 1:binomial(d, m)
#         rhomd[i,i] = round(estat'*fixed(adi[i]*adi[i]',num)*estat, digits = prec)
#     end
#     return rhomd
# end

=#


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

import Base.*, Base.eltype, Base.size
import LinearAlgebra.issymmetric, LinearAlgebra.adjoint

abstract type AbstractLinearOperator end

function *(op::AbstractLinearOperator, v::AbstractArray)
    # @assert size(v) == op.input_shape
    return op.forward(v)
end

# # This function is for computing eigenvalue by eigs.
# # I am still not sure what is the best way to reduce dimension.
# function A_mul_B!(Y::AbstractArray, A::AbstractLinearOperator, B::AbstractArray)
#     for i in 1:size(B,2)
#         Ab = A * B[:,i] # This sometimes return an array of size (N, 1)
#         Ab_squeezed = length(size(Ab)) == 2 ? squeeze(Ab,2) : Ab
#         Y[:,i] .= Ab[:]#_squeezed
#     end
# end

function size(op::AbstractLinearOperator)
    return (op.output_size, op.input_size)
end

function size(op::AbstractLinearOperator, i::Integer)
    if i == 1
        return op.output_size
    elseif i == 2
        return op.input_size
    else
        error("Not implemented.")
    end
end

function eltype(op::AbstractLinearOperator)
    return op.eltype
end

function issymmetric(op::AbstractLinearOperator)
    return op.symmetric
end

mutable struct LinearOperator <: AbstractLinearOperator
    input_shape::Tuple
    output_shape::Tuple
    input_size::Int
    output_size::Int
    forward::Function
    adjoint::Function
    symmetric::Bool
    eltype::DataType
end

function LinearOperator(
        input_shape::Tuple,
        output_shape::Tuple,
        forward::Function,
        adjoint::Function,
        symmetric::Bool = false,
        eltype::DataType = Float32,
    )
    LinearOperator(
        input_shape,
        output_shape,
        prod( s for s in input_shape ),
        prod( s for s in output_shape ),
        forward,
        adjoint,
        symmetric,
        eltype
    )
end

function adjoint(op::LinearOperator)
    return LinearOperator(
        op.output_shape,
        op.input_shape,
        op.output_size,
        op.input_size,
        op.adjoint,
        op.forward,
        op.symmetric,
        op.eltype
    )
end

function adjointnessCheck(A::LinearOperator)
    # Check forward and adjoint operations are correctly implemented.
    x = randn(eltype(A), A.input_shape)
    forward_output = A.forward(x)

    y = randn(eltype(A), A.output_shape)
    adjoint_output = A.adjoint(y)

    return vecdot(forward_output, y) ≈ vecdot(x, adjoint_output)
end

mutable struct StackedLinearOperator <: AbstractLinearOperator
    linops::Vector{LinearOperator}
    input_shape::Tuple
    output_shape::Vector{Tuple}
    input_size::Int
    output_size::Int
    forward!::Function
    adjoint!::Function
    forward::Function
    adjoint::Function
    symmetric::Bool
end

function StackedLinearOperator(linops::AbstractVector{LinearOperator})

    @assert length(unique(linop.input_shape for linop in linops)) == 1
    input_shape = linops[1].input_shape
    output_shape = [linop.output_shape for linop in linops]

    input_size = linops[1].input_size
    output_size = sum( linop.input_size for linop in linops )

    forward! = function(input, output::Vector)
                # Might be better to change it to pmap
                for j in 1:length(output)
                    output[j] .= linops[j].forward(input)
                end
            end

    adjoint! = function(input::Vector, output)
                output .= sum( linops[j].adjoint(input[j]) for j in 1:length(linops) )
            end

    forward_stack = function(input)
                output = [similar(input, dims=sh) for sh in output_shape]
                forward!(input, output)
                return output
            end

    adjoint_stack = function(input)
                output = similar(input, dims=input_shape)
                adjoint!(input, output)
                return output
            end

    return StackedLinearOperator(
        linops,
        input_shape,
        output_shape,
        input_size,
        output_size,
        forward!,
        adjoint!,
        forward_stack,
        adjoint_stack,
        false
    )
end

################################################################################
function IdentityOperator(data_shape::Tuple)
    return LinearOperator(data_shape, data_shape, copy, copy, true)
end

# 3D Second-order differential operator
function SecondOrderDifferentialOperator3D(data_shape::NTuple{3,Int}, dir::NTuple{2,Int})
    I, J, K = data_shape
    diff = zeros(Float32, I, J, K)
    if dir == (1,1)
        diff[end-1,1,1] =  1.0f0
        diff[end,1,1]   = -2.0f0
        diff[1,1,1]     =  1.0f0
    elseif dir == (2,2)
        diff[1,end-1,1] =  1.0f0
        diff[1,end,1]   = -2.0f0
        diff[1,1,1]     =  1.0f0
    elseif dir == (3,3)
        diff[1,1,end-1] =  1.0f0
        diff[1,1,end]   = -2.0f0
        diff[1,1,1]     =  1.0f0
    elseif dir == (1,2) || dir == (2,1)
        # Absorbing the factor √2 to this operator
        diff[end-1,end-1,1] =  sqrt(2.0f0)
        diff[end,1,1]       = -sqrt(2.0f0)
        diff[1,end,1]       = -sqrt(2.0f0)
        diff[1,1,1]         =  sqrt(2.0f0)
    elseif dir == (1,3) || dir == (3,1)
        diff[end-1,1,end-1] =  sqrt(2.0f0)
        diff[end,1,1]       = -sqrt(2.0f0)
        diff[1,1,end]       = -sqrt(2.0f0)
        diff[1,1,1]         =  sqrt(2.0f0)
    elseif dir == (2,3) || dir == (3,2)
        diff[1,end-1,end-1] =  sqrt(2.0f0)
        diff[1,end,1]       = -sqrt(2.0f0)
        diff[1,1,end]       = -sqrt(2.0f0)
        diff[1,1,1]         =  sqrt(2.0f0)
    else
        error("Direction is wrong.")
    end

    diff = to_gpu_or_not_to_gpu(diff)

    dft_mat = rfft3_operator(diff)
    D = dft_mat * diff
    idft_mat = dft_mat'

    D_operator = LinearOperator(data_shape, data_shape,
                    x -> idft_mat * (D .* (dft_mat * x)),
                    x -> idft_mat * (conj.(D) .* (dft_mat * x)))

    return D_operator, D
end


text:
function drawline!(img::Matrix{T}, x0::Int, y0::Int, x1::Int, y1::Int, col::T) where T
    δx = abs(x1 - x0)
    δy = abs(y1 - y0)
    δe = abs(δy / δx)
    er = 0.0

    y = y0
    for x in x0:x1
        img[x, y] = col
        er += δe
        if er > 0.5
            y  += 1
            er -= 1.0
        end
    end

    return img
end

using Images

img = fill(Gray(255.0), 5, 5);
println("\nImage:")
display(img); println()
drawline!(img, 1, 1, 5, 5, Gray(0.0));
println("\nModified image:")
display(img); println()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module SpaceTimeFields

  import Base: +, -, *, show

  export EldredgeRamp, ColoniusRamp, Sinusoid, Gaussian, DGaussian, ConstantProfile, d_dt
  export EmptySpatialField, SpatialGaussian

  include("1dprofiles.jl")
  include("spatialfields.jl")

end


text:
struct DecisionModelSimulationResults <: OperationModelSimulationResults
    variables::FieldResultsByTime
    duals::FieldResultsByTime
    parameters::FieldResultsByTime
    aux_variables::FieldResultsByTime
    expressions::FieldResultsByTime
    forecast_horizon::Int
    container_key_lookup::Dict{String, OptimizationContainerKey}
end

function SimulationProblemResults(
    ::Type{DecisionModel},
    store::SimulationStore,
    model_name::AbstractString,
    problem_params::ModelStoreParams,
    sim_params::SimulationStoreParams,
    path,
    container_key_lookup;
    kwargs...,
)
    name = Symbol(model_name)
    variables = list_decision_model_keys(store, name, STORE_CONTAINER_VARIABLES)
    parameters = list_decision_model_keys(store, name, STORE_CONTAINER_PARAMETERS)
    duals = list_decision_model_keys(store, name, STORE_CONTAINER_DUALS)
    aux_variables = list_decision_model_keys(store, name, STORE_CONTAINER_AUX_VARIABLES)
    expressions = list_decision_model_keys(store, name, STORE_CONTAINER_EXPRESSIONS)

    return SimulationProblemResults{DecisionModelSimulationResults}(
        store,
        model_name,
        problem_params,
        sim_params,
        path,
        DecisionModelSimulationResults(
            _fill_result_value_container(variables),
            _fill_result_value_container(duals),
            _fill_result_value_container(parameters),
            _fill_result_value_container(aux_variables),
            _fill_result_value_container(expressions),
            get_horizon(problem_params),
            container_key_lookup,
        );
        kwargs...,
    )
end

function Base.empty!(res::SimulationProblemResults{DecisionModelSimulationResults})
    foreach(empty!, _get_dicts(res))
    empty!(res.results_timestamps)
    return
end

Base.isempty(res::SimulationProblemResults{DecisionModelSimulationResults}) =
    all(isempty, _get_dicts(res))

# This returns the number of timestamps stored in all containers.
Base.length(res::SimulationProblemResults{DecisionModelSimulationResults}) =
    mapreduce(length, +, _get_dicts(res))

_get_dicts(res::SimulationProblemResults) =
    (y for x in _get_containers(res) for y in values(x))

function get_forecast_horizon(res::SimulationProblemResults{DecisionModelSimulationResults})
    return res.values.forecast_horizon
end

function _get_store_value(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    container_keys::Vector{<:OptimizationContainerKey},
    timestamps,
    ::Nothing,
)
    simulation_store_path = joinpath(get_execution_path(res), "data_store")
    return open_store(HdfSimulationStore, simulation_store_path, "r") do store
        _get_store_value(res, container_keys, timestamps, store)
    end
end

function _get_store_value(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    container_keys::Vector{<:OptimizationContainerKey},
    timestamps,
    store::SimulationStore,
)
    base_power = get_model_base_power(res)
    results =
        Dict{OptimizationContainerKey, SortedDict{Dates.DateTime, DataFrames.DataFrame}}()
    model_name = Symbol(get_model_name(res))
    resolution = get_resolution(res)
    horizon = get_forecast_horizon(res)
    for key in container_keys
        _results = SortedDict{Dates.DateTime, DataFrames.DataFrame}()
        for ts in timestamps
            out = read_result(DataFrames.DataFrame, store, model_name, key, ts)
            if convert_result_to_natural_units(key)
                out .*= base_power
            end
            time_col = range(ts, length=horizon, step=resolution)
            DataFrames.insertcols!(out, 1, :DateTime => time_col)
            _results[ts] = out
        end
        results[key] = _results
    end

    return results
end

function _process_timestamps(
    res::SimulationProblemResults,
    initial_time::Union{Nothing, Dates.DateTime},
    count::Union{Nothing, Int},
)
    if initial_time === nothing
        initial_time = first(get_timestamps(res))
    end

    if initial_time ∉ res.timestamps
        invalid_timestamps = [initial_time]
    else
        if count === nothing
            requested_range = [v for v in res.timestamps if v >= initial_time]
        else
            requested_range =
                collect(range(initial_time, length=count, step=get_interval(res)))
        end
        invalid_timestamps = [v for v in requested_range if v ∉ res.timestamps]
    end
    if !isempty(invalid_timestamps)
        @error "Timestamps $(invalid_timestamps) not stored" get_timestamps(res)
        throw(IS.InvalidValue("Timestamps not stored"))
    end
    return requested_range
end

function _read_variables(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    variable_keys,
    timestamps,
    store,
)
    isempty(variable_keys) && return FieldResultsByTime()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end
    _validate_keys(keys(get_variables(res)), variable_keys)
    same_timestamps = isempty(setdiff(res.results_timestamps, timestamps))
    keys_with_values = [k for (k, v) in get_variables(res) if !isempty(v)]
    same_keys = isempty([n for n in variable_keys if n ∉ keys_with_values])
    if same_timestamps && same_keys
        @debug "reading variables from SimulationsResults"
        vals = filter(p -> (p.first ∈ variable_keys), get_variables(res))
    else
        @debug "reading variables from data store"
        vals = _get_store_value(res, variable_keys, timestamps, store)
    end
    return vals
end

"""
Return the values for the requested variable. It keeps requests when performing multiple retrievals.

# Arguments

  - `args`: Can be a string returned from [`list_variable_names`](@ref) or args that can be
    splatted into a VariableKey.
  - `initial_time::Dates.DateTime` : initial of the requested results
  - `count::Int`: Number of results
  - `store::SimulationStore`: a store that has been opened for reading

# Examples

```julia
read_variable(results, ActivePowerVariable, ThermalStandard)
read_variable(results, "ActivePowerVariable__ThermalStandard")
```
"""
function read_variable(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    args...;
    initial_time::Union{Nothing, Dates.DateTime}=nothing,
    count::Union{Int, Nothing}=nothing,
    store=nothing,
)
    key = _deserialize_key(VariableKey, res, args...)
    timestamps = _process_timestamps(res, initial_time, count)
    return _read_variables(res, [key], timestamps, store)[key]
end

function _read_duals(res::SimulationProblemResults, dual_keys, timestamps, store)
    isempty(dual_keys) && return FieldResultsByTime()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end
    _validate_keys(keys(get_duals(res)), dual_keys)
    same_timestamps = isempty(setdiff(res.results_timestamps, timestamps))
    keys_with_values = [k for (k, v) in get_duals(res) if !isempty(v)]
    same_keys = isempty([n for n in dual_keys if n ∉ keys_with_values])
    if same_timestamps && same_keys
        @debug "reading duals from SimulationsResults"
        vals = filter(p -> (p.first ∈ dual_keys), get_duals(res))
    else
        @debug "reading duals from data store"
        vals = _get_store_value(res, dual_keys, timestamps, store)
    end
    return vals
end

"""
Return the values for the requested dual. It keeps requests when performing multiple retrievals.

# Arguments

  - `args`: Can be a string returned from [`list_dual_names`](@ref) or args that can be
    splatted into a ConstraintKey.
  - `initial_time::Dates.DateTime` : initial of the requested results
  - `count::Int`: Number of results
  - `store::SimulationStore`: a store that has been opened for reading
"""
function read_dual(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    args...;
    initial_time::Union{Nothing, Dates.DateTime}=nothing,
    count::Union{Int, Nothing}=nothing,
    store=nothing,
)
    key = _deserialize_key(ConstraintKey, res, args...)
    timestamps = _process_timestamps(res, initial_time, count)
    return _read_duals(res, [key], timestamps, store)[key]
end

function _read_parameters(res::SimulationProblemResults, parameter_keys, timestamps, store)
    isempty(parameter_keys) && return FieldResultsByTime()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end
    _validate_keys(keys(get_parameters(res)), parameter_keys)
    same_timestamps = isempty(setdiff(res.results_timestamps, timestamps))
    parameters_with_values = [k for (k, v) in get_parameters(res) if !isempty(v)]
    same_parameters = isempty([n for n in parameter_keys if n ∉ parameters_with_values])
    if same_timestamps && same_parameters
        @debug "reading parameters from SimulationsResults"
        vals = filter(p -> (p.first ∈ parameter_keys), get_parameters(res))
    else
        @debug "reading parameters from data store"
        vals = _get_store_value(res, parameter_keys, timestamps, store)
    end
    return vals
end

"""
Return the values for the requested parameter. It keeps requests when performing multiple retrievals.

# Arguments

  - `args`: Can be a string returned from [`list_parameter_names`](@ref) or args that can be
    splatted into a ParameterKey.
  - `initial_time::Dates.DateTime` : initial of the requested results
  - `count::Int`: Number of results
"""
function read_parameter(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    args...;
    time_series_name=nothing,
    initial_time::Union{Nothing, Dates.DateTime}=nothing,
    count::Union{Int, Nothing}=nothing,
    store=nothing,
)
    key = _deserialize_key(ParameterKey, res, args...)
    timestamps = _process_timestamps(res, initial_time, count)
    return _read_parameters(res, [key], timestamps, store)[key]
end

function _read_aux_variables(
    res::SimulationProblemResults,
    aux_variable_keys,
    timestamps,
    store,
)
    isempty(aux_variable_keys) && return FieldResultsByTime()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end
    _validate_keys(keys(get_aux_variables(res)), aux_variable_keys)
    same_timestamps = isempty(setdiff(res.results_timestamps, timestamps))
    keys_with_values = [k for (k, v) in get_aux_variables(res) if !isempty(v)]
    same_keys = isempty([n for n in aux_variable_keys if n ∉ keys_with_values])
    if same_timestamps && same_keys
        @debug "reading aux variables from SimulationsResults"
        vals = filter(p -> (p.first ∈ aux_variable_keys), get_aux_variables(res))
    else
        @debug "reading aux variables from data store"
        vals = _get_store_value(res, aux_variable_keys, timestamps, store)
    end
    return vals
end

"""
Return the values for the requested auxillary variables. It keeps requests when performing multiple retrievals.

# Arguments

  - `args`: Can be a string returned from [`list_aux_variable_names`](@ref) or args that can be
    splatted into a AuxVarKey.
  - `initial_time::Dates.DateTime` : initial of the requested results
  - `count::Int`: Number of results
"""
function read_aux_variable(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    args...;
    time_series_name=nothing,
    initial_time::Union{Nothing, Dates.DateTime}=nothing,
    count::Union{Int, Nothing}=nothing,
    store=nothing,
)
    key = _deserialize_key(AuxVarKey, res, args...)
    timestamps = _process_timestamps(res, initial_time, count)
    return _read_aux_variables(res, [key], timestamps, store)[key]
end

function _read_expressions(
    res::SimulationProblemResults,
    expression_keys,
    timestamps,
    store,
)
    isempty(expression_keys) && return FieldResultsByTime()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end
    _validate_keys(keys(get_expressions(res)), expression_keys)
    same_timestamps = isempty(setdiff(res.results_timestamps, timestamps))
    keys_with_values = [k for (k, v) in get_expressions(res) if !isempty(v)]
    same_keys = isempty([n for n in expression_keys if n ∉ keys_with_values])
    if same_timestamps && same_keys
        @debug "reading expressions from SimulationsResults"
        vals = filter(p -> (p.first ∈ expression_keys), get_expressions(res))
    else
        @debug "reading expressions from data store"
        vals = _get_store_value(res, expression_keys, timestamps, store)
    end
    return vals
end

"""
Return the values for the requested auxillary variables. It keeps requests when performing multiple retrievals.

# Arguments

  - `args`: Can be a string returned from [`list_expression_names`](@ref) or args that can be
    splatted into a ExpressionKey.
  - `initial_time::Dates.DateTime` : initial of the requested results
  - `count::Int`: Number of results
"""
function read_expression(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    args...;
    time_series_name=nothing,
    initial_time::Union{Nothing, Dates.DateTime}=nothing,
    count::Union{Int, Nothing}=nothing,
    store=nothing,
)
    key = _deserialize_key(ExpressionKey, res, args...)
    timestamps = _process_timestamps(res, initial_time, count)
    return _read_expressions(res, [key], timestamps, store)[key]
end

function get_realized_timestamps(
    res::IS.Results;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    timestamps = get_timestamps(res)
    interval = timestamps.step
    resolution = get_resolution(res)
    horizon = get_forecast_horizon(res)
    start_time = isnothing(start_time) ? first(timestamps) : start_time
    end_time =
        isnothing(len) ? last(timestamps) + interval - resolution :
        start_time + (len - 1) * resolution

    requested_range = start_time:resolution:end_time
    available_range =
        first(timestamps):resolution:(last(timestamps) + (horizon - 1) * resolution)
    invalid_timestamps = setdiff(requested_range, available_range)

    if !isempty(invalid_timestamps)
        msg = "Requested time does not match available results"
        @error msg
        throw(IS.InvalidValue(msg))
    end

    return requested_range
end

function read_variables_with_keys(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    meta = RealizedMeta(res, start_time=start_time, len=len)
    timestamps = _process_timestamps(res, meta.start_time, meta.len)
    result_values = _read_variables(res, variables, timestamps, nothing)
    return get_realization(result_values, meta)
end

function read_parameters_with_keys(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    parameters::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    meta = RealizedMeta(res, start_time=start_time, len=len)
    timestamps = _process_timestamps(res, meta.start_time, meta.len)
    result_values = _read_parameters(res, parameters, timestamps, nothing)
    return get_realization(result_values, meta)
end

function read_duals_with_keys(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    duals::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    meta = RealizedMeta(res, start_time=start_time, len=len)
    timestamps = _process_timestamps(res, meta.start_time, meta.len)
    result_values = _read_duals(res, duals, timestamps, nothing)
    return get_realization(result_values, meta)
end

function read_aux_variables_with_keys(
    res::SimulationProblemResults,
    aux_variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    meta = RealizedMeta(res, start_time=start_time, len=len)
    timestamps = _process_timestamps(res, meta.start_time, meta.len)
    result_values = _read_aux_variables(res, aux_variables, timestamps, nothing)
    return get_realization(result_values, meta)
end

function read_expressions_with_keys(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    expressions::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    meta = RealizedMeta(res, start_time=start_time, len=len)
    timestamps = _process_timestamps(res, meta.start_time, meta.len)
    result_values = _read_expressions(res, expressions, timestamps, nothing)
    return get_realization(result_values, meta)
end

"""
Load the simulation results into memory for repeated reads. Running this function twice
overwrites the previously loaded results. This is useful when loading results from remote
locations over network connections.

For each variable/parameter/dual, etc., each element must be the name encoded as a string,
like `"ActivePowerVariable__ThermalStandard"` or a Tuple with its constituent types, like
`(ActivePowerVariable, ThermalStandard)`.

# Arguments

  - `count::Int`: Number of windows to load.
  - `initial_time::Dates.DateTime` : Initial time of first window to load. Defaults to first.
  - `aux_variables::Vector{Union{String, Tuple}}`: Optional list of aux variables to load.
  - `duals::Vector{Union{String, Tuple}}`: Optional list of duals to load.
  - `expressions::Vector{Union{String, Tuple}}`: Optional list of expressions to load.
  - `parameters::Vector{Union{String, Tuple}}`: Optional list of parameters to load.
  - `variables::Vector{Union{String, Tuple}}`: Optional list of variables to load.
"""
function load_results!(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    count::Int;
    initial_time::Union{Dates.DateTime, Nothing}=nothing,
    variables=Vector{Tuple}(),
    duals=Vector{Tuple}(),
    parameters=Vector{Tuple}(),
    aux_variables=Vector{Tuple}(),
    expressions=Vector{Tuple}(),
)
    initial_time = initial_time === nothing ? first(get_timestamps(res)) : initial_time

    res.results_timestamps = _process_timestamps(res, initial_time, count)

    dual_keys = [_deserialize_key(ConstraintKey, res, x...) for x in duals]
    parameter_keys = [_deserialize_key(ParameterKey, res, x...) for x in parameters]
    variable_keys = [_deserialize_key(VariableKey, res, x...) for x in variables]
    aux_variable_keys = [_deserialize_key(AuxVarKey, res, x...) for x in aux_variables]
    expression_keys = [_deserialize_key(ExpressionKey, res, x...) for x in expressions]
    function merge_results(store)
        merge!(
            get_variables(res),
            _read_variables(res, variable_keys, res.results_timestamps, store),
        )
        merge!(get_duals(res), _read_duals(res, dual_keys, res.results_timestamps, store))
        merge!(
            get_parameters(res),
            _read_parameters(res, parameter_keys, res.results_timestamps, store),
        )
        merge!(
            get_aux_variables(res),
            _read_aux_variables(res, aux_variable_keys, res.results_timestamps, store),
        )
        merge!(
            get_expressions(res),
            _read_expressions(res, expression_keys, res.results_timestamps, store),
        )
    end

    if res.store isa InMemorySimulationStore
        merge_results(res.store)
    else
        simulation_store_path = joinpath(res.execution_path, "data_store")
        open_store(HdfSimulationStore, simulation_store_path, "r") do store
            merge_results(store)
        end
    end

    return nothing
end

function _read_optimizer_stats(
    res::SimulationProblemResults{DecisionModelSimulationResults},
    store::SimulationStore,
)
    return read_optimizer_stats(store, Symbol(res.problem))
end


text:
@testset "Convert" begin
    x = variable()
    for (α, β) in ((1 / 2, 1 / 4), (1 / 4, 1 / 2), (1, 2), (2, 1))
        P, Q = Laguerre{α}, Laguerre{β}
        for _ in 1:10
            p = P(rand(1:5, 5))
            q = P(rand(1:5, 5))
            @test convert(P, convert(Q, p))(x) ≈ p(x)
            @test convert(Q, convert(P, q))(x) ≈ q(x)
        end
    end
end


text:
assets(path::String) = joinpath(@__DIR__, "..", "assets", path)

include("basic.jl")
include("texture.jl")


text:
export Schema, FreeSchema, AttrType, Attr,
  SchemaExpr, AttrTypeExpr, AttrExpr, parse_schema

using MLStyle: @match
using ...Meta: strip_lines

# Schema
########

""" The GAT that parameterizes Attributed C-sets
A schema is comprised of a category C, a discrete category D, and a profunctor
Attr : C^op x D → Set. In GAT form, this is given by extending the theory of
categories with two extra types, AttrType for objects of D, and Attr, for elements
of the sets given by the profunctor.
"""
@theory Schema{Ob,Hom,AttrType,Attr} <: Category{Ob,Hom} begin
  AttrType::TYPE
  Attr(dom::Ob,codom::AttrType)::TYPE

  """ Composition is given by the action of the profunctor on C.
  """
  compose(f::Hom(A,B), g::Attr(B,X))::Attr(A,X) ⊣ (A::Ob, B::Ob, X::AttrType)

  (compose(f, compose(g, a)) == compose(compose(f, g), a)
    ⊣ (A::Ob, B::Ob, C::Ob, X::AttrType, f::Hom(A,B), g::Hom(B,C), a::Attr(C, X)))
  compose(id(A), a) == a ⊣ (A::Ob, X::AttrType, a::Attr(A,X))
end

abstract type SchemaExpr{T} <: GATExpr{T} end
abstract type AttrTypeExpr{T} <: SchemaExpr{T} end
abstract type AttrExpr{T} <: SchemaExpr{T} end

@syntax FreeSchema{ObExpr,HomExpr,AttrTypeExpr,AttrExpr} Schema begin
  # should have a normal representation for precompose of a morphism + a generator attribute
  compose(f::Hom, g::Hom) = associate_unit(new(f,g; strict=true), id)
  compose(f::Hom, x::Attr) = associate_unit(new(f,x; strict=true), id)
end

# Type-level representation of Schema
#####################################

function intify(tuple,s)
  findfirst(tuple .== s)
end

function tuplize(xs,f)
  (;[nameof(x) => nameof(f(x)) for x in xs]...)
end

struct SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}
end

const ACSetSchemaDescType = SchemaDescType
const CSetSchemaDescType{obs,homs,doms,codoms} = SchemaDescType{obs,homs,(),(),doms,codoms}

function SchemaDescTypeType(p::Presentation)
  obs,homs,attrtypes,attrs = map(t -> p.generators[t],[:Ob,:Hom,:AttrType,:Attr])
  ob_syms,hom_syms,attrtype_syms,attr_syms = map(xs -> Tuple(nameof.(xs)),
                                                 [obs,homs,attrtypes,attrs])
  homdoms = map(t -> intify(ob_syms,t),tuplize(homs,dom))
  homcodoms = map(t -> intify(ob_syms,t),tuplize(homs,codom))
  attrdoms = map(t -> intify(ob_syms,t),tuplize(attrs,dom))
  attrcodoms = map(t -> intify(attrtype_syms,t),tuplize(attrs,codom))
  SchemaDescType{
    ob_syms,
    hom_syms,
    attrtype_syms,
    attr_syms,
    (;homdoms...,attrdoms...),
    (;homcodoms...,attrcodoms...)
  }
end

ob(::Type{<:SchemaDescType{obs}}) where {obs} = obs
ob_num(::Type{<:SchemaDescType{obs}}, x::Symbol) where {obs} = findfirst(obs .== x)

hom(::Type{<:SchemaDescType{obs,homs}}) where {obs,homs} = homs

function dom_nums(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
  {obs,homs,attrtypes,attrs,doms,codoms}
  Tuple(map(last, filter(((f,_),) -> f ∈ homs, collect(pairs(doms)))))
end

dom(T::Type{<:SchemaDescType}) = map(i -> ob(T)[i], dom_nums(T))
function dom_num(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}},f::Symbol) where
  {obs,homs,attrtypes,attrs,doms,codoms}
  doms[f]
end
dom(T::Type{<:SchemaDescType},f::Symbol) = ob(T)[dom_num(T,f)]

function codom_nums(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
  {obs,homs,attrtypes,attrs,doms,codoms}
  Tuple(map(last, filter(((f,_),) -> f ∈ homs, collect(pairs(codoms)))))
end

codom(T::Type{<:SchemaDescType}) = map(i -> ob(T)[i], codom_nums(T))
function codom_num(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}},f::Symbol) where
  {obs,homs,attrtypes,attrs,doms,codoms}
  codoms[f]
end
codom(T::Type{<:SchemaDescType},f::Symbol) = ob(T)[codom_num(T,f)]

attr(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
  {obs,homs,attrtypes,attrs,doms,codoms} = attrs
attrtype(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
  {obs,homs,attrtypes,attrs,doms,codoms} = attrtypes

function adom_nums(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
  {obs,homs,attrtypes,attrs,doms,codoms}
  Tuple(map(last, filter(((f,_),) -> f ∈ attrs, collect(pairs(doms)))))
end

adom(T::Type{<:SchemaDescType}) = map(i -> ob(T)[i], adom_nums(T))

function acodom_nums(::Type{<:SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
  {obs,homs,attrtypes,attrs,doms,codoms}
  Tuple(map(last, filter(((f,_),) -> f ∈ attrs, collect(pairs(codoms)))))
end

acodom(T::Type{<:SchemaDescType}) = map(i -> attrtype(T)[i], acodom_nums(T))

struct SchemaDesc
  obs::Vector{Symbol}
  homs::Vector{Symbol}
  attrtypes::Vector{Symbol}
  attrs::Vector{Symbol}
  doms::Dict{Symbol,Symbol}
  codoms::Dict{Symbol,Symbol}
  function SchemaDesc(obs::Vector{Symbol},
                      homs::Vector{Symbol},
                      attrtypes::Vector{Symbol},
                      attrs::Vector{Symbol},
                      doms::Dict{Symbol,Symbol},
                      codoms::Dict{Symbol,Symbol})
    new(obs, homs, attrtypes, attrs, doms, codoms)
  end

  function SchemaDesc(::Type{SchemaDescType{obs,homs,attrtypes,attrs,doms,codoms}}) where
      {obs,homs,attrtypes,attrs,doms,codoms}
    homdoms = [f => obs[doms[f]] for f in homs]
    homcodoms = [f => obs[codoms[f]] for f in homs]
    attrdoms = [a => obs[doms[a]] for a in attrs]
    attrcodoms = [a => attrtypes[codoms[a]] for a in attrs]
    new(
      Symbol[obs...],
      Symbol[homs...],
      Symbol[attrtypes...],
      Symbol[attrs...],
      Dict{Symbol,Symbol}(homdoms...,attrdoms...),
      Dict{Symbol,Symbol}(homcodoms...,attrcodoms...)
    )
  end
end

ob_num(S::SchemaDesc, x::Symbol) where {obs} = findfirst(S.obs .== x)

function SchemaDescTypeType(s::SchemaDesc)
  SchemaDescType{
    Tuple(s.obs),
    Tuple(s.homs),
    Tuple(s.attrtypes),
    Tuple(s.attrs),
    NamedTuple(map([s.doms...]) do (f,d)
                 (f,intify(s.obs, d))
               end),
    NamedTuple(map([s.codoms...]) do (f,cd)
                 if f ∈ s.homs
                   (f,intify(s.obs, cd))
                 else
                   (f,intify(s.attrtypes,cd))
                 end
               end)
  }
end

function Presentation(S::Type{T}) where T <: SchemaDescType
  pres = Presentation(FreeSchema)

  obs = map(x -> Ob(FreeSchema, x), collect(ob(S)))
  homs = map(Hom, collect(hom(S)), obs[collect(dom_nums(S))], obs[collect(codom_nums(S))])
  datas = map(x -> AttrType(FreeSchema.AttrType, x), collect(attrtype(S)))
  attrs = map(FreeSchema.Attr, collect(attr(S)), obs[collect(adom_nums(S))], datas[collect(acodom_nums(S))])

  map(gens -> add_generators!(pres, gens), [obs, homs, datas, attrs])
  return pres
end


text:
export iauPn06
"""
Precession-nutation, IAU 2006 model:  a multi-purpose function,
supporting classical (equinox-based) use directly and CIO-based use
indirectly.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double          TT as a 2-part Julian Date (Note 1)
   dpsi,deps    double          nutation (Note 2)

Returned:
   epsa         double          mean obliquity (Note 3)
   rb           double[3][3]    frame bias matrix (Note 4)
   rp           double[3][3]    precession matrix (Note 5)
   rbp          double[3][3]    bias-precession matrix (Note 6)
   rn           double[3][3]    nutation matrix (Note 7)
   rbpn         double[3][3]    GCRS-to-true matrix (Note 8)

Notes:

   1.  The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

            date1          date2

         2450123.7           0.0       (JD method)
         2451545.0       -1421.3       (J2000 method)
         2400000.5       50123.2       (MJD method)
         2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2.  The caller is responsible for providing the nutation components;
   they are in longitude and obliquity, in radians and are with
   respect to the equinox and ecliptic of date.  For high-accuracy
   applications, free core nutation should be included as well as
   any other relevant corrections to the position of the CIP.

   3.  The returned mean obliquity is consistent with the IAU 2006
   precession.

   4.  The matrix rb transforms vectors from GCRS to J2000.0 mean
   equator and equinox by applying frame bias.

   5.  The matrix rp transforms vectors from J2000.0 mean equator and
   equinox to mean equator and equinox of date by applying
   precession.

   6.  The matrix rbp transforms vectors from GCRS to mean equator and
   equinox of date by applying frame bias then precession.  It is
   the product rp x rb.

   7.  The matrix rn transforms vectors from mean equator and equinox
   of date to true equator and equinox of date by applying the
   nutation (luni-solar + planetary).

   8.  The matrix rbpn transforms vectors from GCRS to true equator and
   equinox of date.  It is the product rn x rbp, applying frame
   bias, precession and nutation in that order.

   9.  The X,Y,Z coordinates of the Celestial Intermediate Pole are
   elements (3,1-3) of the GCRS-to-true matrix, i.e. rbpn[2][0-2].

   10. It is permissible to re-use the same array in the returned
   arguments.  The arrays are filled in the stated order.

Called:
   iauPfw06     bias-precession F-W angles, IAU 2006
   iauFw2m      F-W angles to r-matrix
   iauCr        copy r-matrix
   iauTr        transpose r-matrix
   iauRxr       product of two r-matrices

References:

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 November 14

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPn06(double date1, double date2, double dpsi, double deps,
#              double *epsa,
#              double rb[3][3], double rp[3][3], double rbp[3][3],
#              double rn[3][3], double rbpn[3][3])

function iauPn06(date1::Real, date2::Real, dpsi::Real, deps::Real)
   # Allocate return values
   ref_epsa = Ref{Float64}(0.0)
   rb       = zeros(Float64, 3, 3)
   rp       = zeros(Float64, 3, 3)
   rbp      = zeros(Float64, 3, 3)
   rn       = zeros(Float64, 3, 3)
   rbpn     = zeros(Float64, 3, 3)

   ccall((:iauPn06, libsofa_c), Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
         Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
         convert(Float64, date1),
         convert(Float64, date2),
         convert(Float64, dpsi),
         convert(Float64, deps),
         ref_epsa, rb, rp, rbp, rn, rbpn)

   return ref_epsa[], SMatrix{3,3}(rb'), SMatrix{3,3}(rp'),
      SMatrix{3,3}(rbp'), SMatrix{3,3}(rn'), SMatrix{3,3}(rbpn')
end


text:
module LSDescent

export bfgs, bfgs_Stop, bfgs_StopLS
export L_bfgs, L_bfgs_Stop, L_bfgs_StopLS
export Newton_StopLS, Newton_Stop, Newton_Spectral
export CG_generic, formula_FR, formula_PR, formula_HS, formula_HZ
#export CG_FR, CG_PR, CG_Hs, CG_HZ

using LinearAlgebra

using NLPModels
using OneDmin

using LinearOperators
using Logging
using SolverCore   # Pour avoir les utilitaires d'affichage log_header et log_row

using Stopping

include("Newton/NewtonStopLS.jl")
include("Newton/NewtonStop.jl")
include("Newton/NewtonSolver.jl")
include("BFGS/bfgsSolver.jl")
include("BFGS/bfgsStop.jl")
include("BFGS/bfgsStopLS.jl")
include("BFGS/FormuleN2.jl")

include("ConjugateGradient/CG_generic.jl")
include("ConjugateGradient/formulae.jl")

#include("ConjugateGradient/CG_FR.jl")
#include("ConjugateGradient/CG_PR.jl")
#include("ConjugateGradient/CG_HS.jl")
#include("ConjugateGradient/CG_HZ.jl")

end # module

max_stars_repo_name: tmigot/LSDescent

max_issues_repo_name: tmigot/LSDescent

max_forks_repo_name: tmigot/LSDescent


text:
#include("testS2Cell.jl")
include("testS1Angle.jl")
include("testS1ChordAngle.jl")
include("testS2Point.jl")


text:
isdefined(Base, :__precompile__) && __precompile__()

module FlexPlan

# import Compat
import JuMP
import Memento
import PowerModels
import PowerModelsACDC
const _PM = PowerModels
const _PMACDC = PowerModelsACDC
import InfrastructureModels
#import InfrastructureModels: ids, ref, var, con, sol, nw_ids, nws, optimize_model!, @im_fields
const _IM = InfrastructureModels
import MathOptInterface
const _MOI = MathOptInterface

import JuMP: optimizer_with_attributes
export optimizer_with_attributes

using Printf

# Create our module level logger (this will get precompiled)
const _LOGGER = Memento.getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(PowerModels)`
# NOTE: If this line is not included then the precompiled `_PM._LOGGER` won't be registered at runtime.
__init__() = Memento.register(_LOGGER)

include("prob/storage_tnep.jl")
include("prob/flexible_tnep.jl")
include("prob/stochastic_flexible_tnep.jl")
include("prob/simple_stochastic_flexible_tnep.jl")

include("io/parse.jl")
include("io/scale.jl")
include("io/time_series.jl")
include("io/multinetwork.jl")
include("io/plot_geo_data.jl")

include("core/types.jl")
include("core/dimensions.jl")
include("core/variable.jl")
include("core/variableconv.jl")
include("core/variabledcgrid.jl")
include("core/gen.jl")
include("core/flexible_demand.jl")
include("core/storage.jl")
include("core/objective.jl")
include("core/ref_extension.jl")
include("core/constraint_template.jl")
include("core/constraint.jl")
include("core/line_replacement.jl")
include("core/distribution.jl")
include("core/td_coupling.jl")

include("form/bf.jl")
include("form/bfarad.jl")
include("formconv/dcp.jl")

# Submodules

include("json_converter/json_converter.jl")
using .JSONConverter

include("td_decoupling/td_decoupling.jl")
using .TDDecoupling

end


text:
# Quasi-Newton approximations as linear operators.

import Base.push!
import LinearAlgebra.diag

include("lbfgs.jl")
include("lsr1.jl")


text:
function hello_world_and_no_arguments(;)
    println("Hello World!")
    return (;) # Attention! If no values are to be returned, then return an empty named tuple.
end

fns = (;hello_world_and_no_arguments, )