text:
#!/usr/bin/env julia

Fract_KelvinVoigt =  RheoModelClass(
        # Model name
        name="fractKV",
        # Model parameters,
        p = (:cₐ, :a, :cᵦ, :β),
        # Relaxation modulus
        G = quote
              cₐ*t^(-a)/gamma(1 - a) + cᵦ*t^(-β)/gamma(1 - β)
            end,
        # Creep modulus
        J = quote
              (t^(a)/cₐ)*mittleff(a - β, 1 + a, -cᵦ*t^(a - β)/cₐ)
            end,
        # Storage modulus
        Gp = quote
                cₐ*ω^a*cos(a*π/2) + cᵦ*ω^β*cos(β*π/2)
             end,
        # Loss modulus
        Gpp = quote
                cₐ*ω^a*sin(a*π/2) + cᵦ*ω^β*sin(β*π/2)
              end,
        # Constraints
        constraint = quote
                 all([  (a<1) & (a>0)
                         (β<1) & (β>0)
                          -a+β < 0])
                end,
        # Network
        info= "
                ________ ╱╲ ________
               |         ╲╱  cₐ, a  |
           ____|                    |____
               |                    |
               |________ ╱╲ ________|
                         ╲╱  cᵦ, β
                "
        )

FractS_KelvinVoigt =  RheoModelClass(
        # Model name
        name="fractSpringKV",
        # Model parameters,
        p = (:cₐ, :a, :k),
        # Relaxation modulus
        G = quote
              cₐ*t^(-a)/gamma(1 - a) + k
            end,
        # Creep modulus
        J = quote
              (t^(a)/cₐ)*mittleff(a, 1 + a, -k*t^a/cₐ)
            end,
        # Storage modulus
        Gp = quote
                #if a!=1.0        #  This is to be removed as the condition contradicts the constraints
                    cₐ*ω^a*cos(a*π/2) + k
                #else
                #    k
                #end
             end,
        # Loss modulus
        Gpp = quote
                cₐ*ω^a*sin(a*π/2)
              end,
        # Constraints
        constraint = quote
                 (a<1) & (a>0)
                end,
        # Network
        info= "
                ________ ╱╲ ________
               |         ╲╱  cₐ, a  |
           ____|                    |____
               |                    |
               |____╱╲  ╱╲  ╱╲  ____|
                      ╲╱  ╲╱  ╲╱
                                k
                "
        )

FractD_KelvinVoigt =  RheoModelClass(
        # Model name
        name="fractDashpotKV",
        # Model parameters,
        p = (:η, :cᵦ, :β),
        # Relaxation modulus
        G = quote
            cᵦ*t^(-β)/gamma(1 - β)
            #t!=0.0 ? cᵦ*t^(-β)/gamma(1 - β) : Inf
        end,
        # Creep modulus
        J = quote
              (t/η)*mittleff(1 - β, 1 + 1, -cᵦ*t^(1.0 - β)/η)
            end,
        # Storage modulus
        Gp = quote
                # cosine floating point error work-around
                if β==1.0
                    0.0
                else
                    cᵦ*ω^β*cos(β*π/2)
                end
             end,
        # Loss modulus
        Gpp = quote
                η*ω + cᵦ*ω^β*sin(β*π/2)
              end,
        # Constraints
        constraint = quote
                 (β<1) & (β>0)
                end,
        # Network
        info= "
                        ___
                _________| |________
               |        _|_| η      |
           ____|                    |____
               |                    |
               |________ ╱╲ ________|
                         ╲╱  cᵦ, β
                "
                )

KelvinVoigt =  RheoModelClass(
        # Model name
        name="KelvinVoigt",
        # Model parameters,
        p = (:η, :k),
        # Relaxation modulus  -    Should be singular
        G = quote
            k
        end,
        # Creep modulus
        J = quote
              (1 - exp(-k*t/η))/k
            end,
        # Storage modulus
        Gp = quote
                k
             end,
        # Loss modulus
        Gpp = quote
                η*ω
              end,
        # Network
        info= "
                        ___
                _________| |________
               |        _|_| η      |
           ____|                    |____
               |                    |
               |____╱╲  ╱╲  ╱╲  ____|
                      ╲╱  ╲╱  ╲╱
                                k
                "
                )


text:
using GAP_pkg_ferret
using Test

import GAP_pkg_ferret
println("TODO")

ext: jl


text:
using PastaQ
using ITensors
using ITensorsGPU
using Random

Random.seed!(1234)

using CUDA: allowscalar
allowscalar(false)

# Set to identity to run on CPU
gpu = cu

N = 40     # Number of qubits
depth = 18 # Depth of the quantum circuit

# Starting state on GPU
# `|ψ₀⟩ = |0,0,…,0⟩`
println("Make the starting state |ψ₀⟩ = |0,0,…,0⟩...")
ψ₀ = gpu(qubits(N))
@show maxlinkdim(ψ₀)
println()

# Generate random quantum circuit built out of
# layers of single-qubit random rotations + `CX`
# gates, alternating between even and of odd layers.
println("Random circuit of depth $depth on $N qubits:")
gates = randomcircuit(N, depth)
display(gates)
println()

# Build the circuit as tensors and move them to GPU
gate_tensors = gpu.(complex.(buildcircuit(ψ₀, gates)))

# Obtain the approximate circuit evolution
println("Approximating random circuit evolution Û|0,0,…,0⟩...")
ψ = @time runcircuit(ψ₀, gate_tensors)
@show maxlinkdim(ψ)
println()


text:
module PauliParticle3d

    include("pauli_particle_3d/symmetric_field.jl")
    include("pauli_particle_3d/theta_pinch.jl")
    include("pauli_particle_3d/solovev_iter.jl")
    include("pauli_particle_3d/solovev_iter_xpoint.jl")
    include("pauli_particle_3d/tokamak_iter_cylindrical.jl")
    include("pauli_particle_3d/tokamak_small_cartesian.jl")
    include("pauli_particle_3d/tokamak_small_cylindrical.jl")
    include("pauli_particle_3d/tokamak_small_toroidal.jl")

end


text:
f = open("jawiki-country.txt", "r")

for line in readlines(f)
  if ismatch(r"(\[+)*(ファイル|File)",line)
    m = match(r"(ファイル|File)\:(.*\.\w+)\|", line)
    println(m.captures[2])
  end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "cohomCalg"
version = v"0.32"

sources = [
    GitSource("https://github.com/BenjaminJurke/cohomCalg",
              "c663c8e37cceab3cc0b2bcc57d35cb895930ab1f"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd cohomCalg*

atomic_patch -p1 ../patches/Makefile.patch
make -j${nproc}

# install
mkdir -p ${bindir}
cp bin/cohomcalg  ${bindir}
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("cohomcalg", :cohomcalg),
]

# Dependencies that must be installed before this package can be built
dependencies = [
#    Dependency("GMP_jll", v"6.2.0"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ---
# title: 283. Move Zeroes
# id: problem283
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/move-zeroes/description/>
# hidden: true
# ---
#
# Given an array `nums`, write a function to move all `0`'s to the end of it
# while maintaining the relative order of the non-zero elements.
#
# **Example:**
#
#
#
#     Input: [0,1,0,3,12]
#     Output: [1,3,12,0,0]
#
# **Note** :
#
#   1. You must do this **in-place** without making a copy of the array.
#   2. Minimize the total number of operations.
#
#
## @lc code=start
using LeetCode

## use double pointers
function move_zeros!(nums::Vector{Int})::Nothing
    p2 = 1
    for p1 in 1:length(nums)
        if nums[p1] != 0
            nums[p1], nums[p2] = nums[p2], nums[p1]
            p2 += 1
        end
    end
    return nothing
end

## @lc code=end


text:
using Documenter, Measurements

makedocs(
    modules = [Measurements],
    format = :html,
    sitename = "Measurements",
    pages    = Any[
        "Introduction" => "index.md",
        "Installation" => "installation.md",
        "Usage"        => "usage.md",
        "Examples"     => "examples.md",
        "Performance"  => "performance.md",
        "Development"  => "todo.md",
        "Appendix"     => "appendix.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaPhysics/Measurements.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
immutable LogNormal <: ContinuousUnivariateDistribution
    meanlog::Float64
    sdlog::Float64
    function LogNormal(ml::Real, sdl::Real)
    	sdl > zero(sdl) || error("sdlog must be positive")
    	new(float64(ml), float64(sdl))
    end
end

LogNormal(ml::Real) = LogNormal(ml, 1.0)
LogNormal() = LogNormal(0.0, 1.0)

## Support
@continuous_distr_support LogNormal 0.0 Inf

## Properties
mean(d::LogNormal) = exp(d.meanlog + d.sdlog^2 / 2)

median(d::LogNormal) = exp(d.meanlog)

mode(d::LogNormal) = exp(d.meanlog - d.sdlog^2)

function var(d::LogNormal)
    sigsq = d.sdlog^2
    (exp(sigsq) - 1) * exp(2d.meanlog + sigsq)
end

function skewness(d::LogNormal)
    (exp(d.sdlog^2) + 2.0) * sqrt(exp(d.sdlog^2) - 1.0)
end

function kurtosis(d::LogNormal)
   exp(4.0 * d.sdlog^2) + 2.0 * exp(3.0 * d.sdlog^2) +
        3.0 * exp(2.0 * d.sdlog^2) - 6.0
end

entropy(d::LogNormal) = 0.5 + 0.5 * log(2.0 * pi * d.sdlog^2) + d.meanlog

## Functions
pdf(d::LogNormal, x::Real) = pdf(Normal(d.meanlog,d.sdlog),log(x))/x
logpdf(d::LogNormal, x::Real) = (lx = log(x); logpdf(Normal(d.meanlog,d.sdlog),lx)-lx)

cdf(d::LogNormal, q::Real) = q <= zero(q) ? 0.0 : cdf(Normal(d.meanlog,d.sdlog),log(q))
ccdf(d::LogNormal, q::Real) = q <= zero(q) ? 1.0 : ccdf(Normal(d.meanlog,d.sdlog),log(q))
logcdf(d::LogNormal, q::Real) = q <= zero(q) ? -Inf : logcdf(Normal(d.meanlog,d.sdlog),log(q))
logccdf(d::LogNormal, q::Real) = q <= zero(q) ? 0.0 : logccdf(Normal(d.meanlog,d.sdlog),log(q))

quantile(d::LogNormal, p::Real) = exp(quantile(Normal(d.meanlog,d.sdlog),p))
cquantile(d::LogNormal, p::Real) = exp(cquantile(Normal(d.meanlog,d.sdlog),p))
invlogcdf(d::LogNormal, p::Real) = exp(invlogcdf(Normal(d.meanlog,d.sdlog),p))
invlogccdf(d::LogNormal, p::Real) = exp(invlogccdf(Normal(d.meanlog,d.sdlog),p))

function gradloglik(d::LogNormal, x::Float64)
  insupport(LogNormal, x) ? - ((log(x) - d.meanlog) / (d.sdlog^2) + 1.0) / x : 0.0
end

# mgf(d::LogNormal)
# cf(d::LogNormal)

## Sampling
rand(d::LogNormal) = exp(rand(Normal(d.meanlog,d.sdlog)))

## Fitting
function fit_mle{T <: Real}(::Type{LogNormal}, x::Array{T})
    lx = log(x)
    LogNormal(mean(lx), std(lx))
end


text:
import Base: ==, +, -, *, /, Broadcast
using Base.Cartesian

"""
    DenseOperator(b1[, b2, data])

Dense array implementation of Operator.

The matrix consisting of complex floats is stored in the `data` field.
"""
mutable struct DenseOperator{BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} <: DataOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    data::Matrix{ComplexF64}
    function DenseOperator{BL,BR,T}(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}}
        if !(length(b1) == size(data, 1) && length(b2) == size(data, 2))
            throw(DimensionMismatch())
        end
        new(b1, b2, data)
    end
end

DenseOperator{BL,BR}(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{BL,BR,T}(b1, b2, data)
DenseOperator(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{BL,BR,T}(b1, b2, data)
DenseOperator(b1::Basis, b2::Basis, data) = DenseOperator(b1, b2, convert(Matrix{ComplexF64}, data))
DenseOperator(b::Basis, data) = DenseOperator(b, b, data)
DenseOperator(b1::Basis, b2::Basis) = DenseOperator(b1, b2, zeros(ComplexF64, length(b1), length(b2)))
DenseOperator{B1,B2}(b1::B1, b2::B2) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(b1, b2, zeros(ComplexF64, length(b1), length(b2)))
DenseOperator(b::Basis) = DenseOperator(b, b)
DenseOperator(op::AbstractOperator) = dense(op)

Base.copy(x::T) where T<:DataOperator = T(x.basis_l, x.basis_r, copy(x.data))

"""
    dense(op::AbstractOperator)

Convert an arbitrary Operator into a [`DenseOperator`](@ref).
"""
dense(x::DenseOperator) = copy(x)

==(x::DenseOperator, y::DenseOperator) = false
==(x::T, y::T) where T<:DenseOperator = (x.data == y.data)

# Arithmetic operations
+(a::T, b::T) where T<:DenseOperator = T(a.basis_l, a.basis_r, a.data+b.data)
+(a::DenseOperator, b::DenseOperator) = throw(IncompatibleBases())

-(a::T) where T<:DenseOperator = T(a.basis_l, a.basis_r, -a.data)
-(a::T, b::T) where T<:DenseOperator = T(a.basis_l, a.basis_r, a.data-b.data)
-(a::DenseOperator, b::DenseOperator) = throw(IncompatibleBases())

*(a::DenseOperator{BL,BR}, b::Ket{BR}) where {BL<:Basis,BR<:Basis} = Ket{BL}(a.basis_l, a.data*b.data)
*(a::DenseOperator, b::Ket) = throw(IncompatibleBases())
*(a::Bra{BL}, b::DenseOperator{BL,BR}) where {BL<:Basis,BR<:Basis} = Bra{BR}(b.basis_r, transpose(b.data)*a.data)
*(a::Bra, b::DenseOperator) = throw(IncompatibleBases())
*(a::DenseOperator{B1,B2,T}, b::DenseOperator{B2,B3,T}) where {B1<:Basis,B2<:Basis,B3<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{B1,B3,T}(a.basis_l, b.basis_r, a.data*b.data)
*(a::DenseOperator, b::DenseOperator) = throw(IncompatibleBases())
*(a::DenseOperator, b::Number) = DenseOperator(a.basis_l, a.basis_r, complex(b)*a.data)
*(a::Number, b::DenseOperator) = DenseOperator(b.basis_l, b.basis_r, complex(a)*b.data)
function *(op1::AbstractOperator{B1,B2}, op2::DenseOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}
    result = DenseOperator{B1,B3}(op1.basis_l, op2.basis_r)
    gemm!(Complex(1.), op1, op2, Complex(0.), result)
    return result
end
function *(op1::DenseOperator{B1,B2}, op2::AbstractOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}
    result = DenseOperator{B1,B3}(op1.basis_l, op2.basis_r)
    gemm!(Complex(1.), op1, op2, Complex(0.), result)
    return result
end
function *(op::AbstractOperator{BL,BR}, psi::Ket{BR}) where {BL<:Basis,BR<:Basis}
    result = Ket{BL}(op.basis_l)
    gemv!(Complex(1.), op, psi, Complex(0.), result)
    return result
end
function *(psi::Bra{BL}, op::AbstractOperator{BL,BR}) where {BL<:Basis,BR<:Basis}
    result = Bra{BR}(op.basis_r)
    gemv!(Complex(1.), psi, op, Complex(0.), result)
    return result
end

/(a::DenseOperator, b::Number) = DenseOperator(a.basis_l, a.basis_r, a.data/complex(b))

dagger(x::DenseOperator) = DenseOperator(x.basis_r, x.basis_l, x.data')

ishermitian(A::DenseOperator) = false
ishermitian(A::DenseOperator{B,B}) where B<:Basis = ishermitian(A.data)

tensor(a::DenseOperator, b::DenseOperator) = DenseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))

conj(a::DenseOperator) = DenseOperator(a.basis_l, a.basis_r, conj(a.data))
conj!(a::DenseOperator) = conj!(a.data)

transpose(op::DenseOperator{BL,BR,T}) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{BR,BL,T}(op.basis_r, op.basis_l, T(transpose(op.data)))

"""
    tensor(x::Ket, y::Bra)

Outer product ``|x⟩⟨y|`` of the given states.
"""
tensor(a::Ket, b::Bra) = DenseOperator(a.basis, b.basis, reshape(kron(b.data, a.data), prod(a.basis.shape), prod(b.basis.shape)))

tr(op::DenseOperator{B,B}) where B<:Basis = tr(op.data)

function ptrace(a::DenseOperator, indices::Vector{Int})
    check_ptrace_arguments(a, indices)
    rank = length(a.basis_l.shape)
    result = _ptrace(Val{rank}, a.data, a.basis_l.shape, a.basis_r.shape, indices)
    return DenseOperator(ptrace(a.basis_l, indices), ptrace(a.basis_r, indices), result)
end

function ptrace(psi::Ket, indices::Vector{Int})
    check_ptrace_arguments(psi, indices)
    b = basis(psi)
    b_ = ptrace(b, indices)
    rank = length(b.shape)
    result = _ptrace_ket(Val{rank}, psi.data, b.shape, indices)
    return DenseOperator(b_, b_, result)
end
function ptrace(psi::Bra, indices::Vector{Int})
    check_ptrace_arguments(psi, indices)
    b = basis(psi)
    b_ = ptrace(b, indices)
    rank = length(b.shape)
    result = _ptrace_bra(Val{rank}, psi.data, b.shape, indices)
    return DenseOperator(b_, b_, result)
end

normalize!(op::DenseOperator) = (rmul!(op.data, 1.0/tr(op)); nothing)

function expect(op::DenseOperator{B,B}, state::Ket{B}) where B<:Basis
    state.data' * op.data * state.data
end

function expect(op::DenseOperator{B1,B2}, state::AbstractOperator{B2,B2}) where {B1<:Basis,B2<:Basis}
    result = ComplexF64(0.)
    @inbounds for i=1:size(op.data, 1), j=1:size(op.data,2)
        result += op.data[i,j]*state.data[j,i]
    end
    result
end

function exp(op::T) where {B<:Basis,T<:DenseOperator{B,B}}
    return T(op.basis_l, op.basis_r, exp(op.data))
end

function permutesystems(a::DenseOperator{B1,B2}, perm::Vector{Int}) where {B1<:CompositeBasis,B2<:CompositeBasis}
    @assert length(a.basis_l.bases) == length(a.basis_r.bases) == length(perm)
    @assert isperm(perm)
    data = reshape(a.data, [a.basis_l.shape; a.basis_r.shape]...)
    data = permutedims(data, [perm; perm .+ length(perm)])
    data = reshape(data, length(a.basis_l), length(a.basis_r))
    DenseOperator(permutesystems(a.basis_l, perm), permutesystems(a.basis_r, perm), data)
end

identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:DenseOperator} = DenseOperator(b1, b2, Matrix{ComplexF64}(I, length(b1), length(b2)))

"""
    projector(a::Ket, b::Bra)

Projection operator ``|a⟩⟨b|``.
"""
projector(a::Ket, b::Bra) = tensor(a, b)
"""
    projector(a::Ket)

Projection operator ``|a⟩⟨a|``.
"""
projector(a::Ket) = tensor(a, dagger(a))
"""
    projector(a::Bra)

Projection operator ``|a⟩⟨a|``.
"""
projector(a::Bra) = tensor(dagger(a), a)

"""
    dm(a::StateVector)

Create density matrix ``|a⟩⟨a|``. Same as `projector(a)`.
"""
dm(x::Ket) = tensor(x, dagger(x))
dm(x::Bra) = tensor(dagger(x), x)

# Partial trace implementation for dense operators.
function _strides(shape::Vector{Int})
    N = length(shape)
    S = zeros(Int, N)
    S[1] = 1
    for m=2:N
        S[m] = S[m-1]*shape[m-1]
    end
    return S
end

# Dense operator version
@generated function _ptrace(::Type{Val{RANK}}, a::Matrix{ComplexF64},
                            shape_l::Vector{Int}, shape_r::Vector{Int},
                            indices::Vector{Int}) where RANK
    return quote
        a_strides_l = _strides(shape_l)
        result_shape_l = copy(shape_l)
        result_shape_l[indices] .= 1
        result_strides_l = _strides(result_shape_l)
        a_strides_r = _strides(shape_r)
        result_shape_r = copy(shape_r)
        result_shape_r[indices] .= 1
        result_strides_r = _strides(result_shape_r)
        N_result_l = prod(result_shape_l)
        N_result_r = prod(result_shape_r)
        result = zeros(ComplexF64, N_result_l, N_result_r)
        @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1))
        @nloops $RANK ir (d->1:shape_r[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides_r[d]; if !(d in indices) Jr_d+=result_strides_r[d] end)) begin
            @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1))
            @nloops $RANK il (k->1:shape_l[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides_l[k]; continue end)) (k->(Il_k+=a_strides_l[k]; if !(k in indices) Jl_k+=result_strides_l[k] end)) begin
                result[Jl_0, Jr_0] += a[Il_0, Ir_0]
            end
        end
        return result
    end
end

@generated function _ptrace_ket(::Type{Val{RANK}}, a::Vector{ComplexF64},
                            shape::Vector{Int}, indices::Vector{Int}) where RANK
    return quote
        a_strides = _strides(shape)
        result_shape = copy(shape)
        result_shape[indices] .= 1
        result_strides = _strides(result_shape)
        N_result = prod(result_shape)
        result = zeros(ComplexF64, N_result, N_result)
        @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1))
        @nloops $RANK ir (d->1:shape[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides[d]; if !(d in indices) Jr_d+=result_strides[d] end)) begin
            @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1))
            @nloops $RANK il (k->1:shape[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides[k]; continue end)) (k->(Il_k+=a_strides[k]; if !(k in indices) Jl_k+=result_strides[k] end)) begin
                result[Jl_0, Jr_0] += a[Il_0]*conj(a[Ir_0])
            end
        end
        return result
    end
end

@generated function _ptrace_bra(::Type{Val{RANK}}, a::Vector{ComplexF64},
                            shape::Vector{Int}, indices::Vector{Int}) where RANK
    return quote
        a_strides = _strides(shape)
        result_shape = copy(shape)
        result_shape[indices] .= 1
        result_strides = _strides(result_shape)
        N_result = prod(result_shape)
        result = zeros(ComplexF64, N_result, N_result)
        @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1))
        @nloops $RANK ir (d->1:shape[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides[d]; if !(d in indices) Jr_d+=result_strides[d] end)) begin
            @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1))
            @nloops $RANK il (k->1:shape[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides[k]; continue end)) (k->(Il_k+=a_strides[k]; if !(k in indices) Jl_k+=result_strides[k] end)) begin
                result[Jl_0, Jr_0] += conj(a[Il_0])*a[Ir_0]
            end
        end
        return result
    end
end

# Fast in-place multiplication
gemm!(alpha, a::Matrix{ComplexF64}, b::Matrix{ComplexF64}, beta, result::Matrix{ComplexF64}) = BLAS.gemm!('N', 'N', convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result)
gemv!(alpha, a::Matrix{ComplexF64}, b::Vector{ComplexF64}, beta, result::Vector{ComplexF64}) = BLAS.gemv!('N', convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result)
gemv!(alpha, a::Vector{ComplexF64}, b::Matrix{ComplexF64}, beta, result::Vector{ComplexF64}) = BLAS.gemv!('T', convert(ComplexF64, alpha), b, a, convert(ComplexF64, beta), result)

gemm!(alpha, a::DenseOperator{B1,B2}, b::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = gemm!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data)
gemv!(alpha, a::DenseOperator{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis} = gemv!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data)
gemv!(alpha, a::Bra{B1}, b::DenseOperator{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis} = gemv!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data)

# Multiplication for Operators in terms of their gemv! implementation
function gemm!(alpha, M::AbstractOperator{B1,B2}, b::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}
    for i=1:size(b.data, 2)
        bket = Ket(b.basis_l, b.data[:,i])
        resultket = Ket(M.basis_l, result.data[:,i])
        gemv!(alpha, M, bket, beta, resultket)
        result.data[:,i] = resultket.data
    end
end

function gemm!(alpha, b::DenseOperator{B1,B2}, M::AbstractOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}
    for i=1:size(b.data, 1)
        bbra = Bra(b.basis_r, vec(b.data[i,:]))
        resultbra = Bra(M.basis_r, vec(result.data[i,:]))
        gemv!(alpha, bbra, M, beta, resultbra)
        result.data[i,:] = resultbra.data
    end
end

# Broadcasting
Base.size(A::DataOperator) = size(A.data)
@inline Base.axes(A::DataOperator) = axes(A.data)
Base.broadcastable(A::DataOperator) = A

# Custom broadcasting styles
abstract type DataOperatorStyle{BL<:Basis,BR<:Basis} <: Broadcast.BroadcastStyle end
struct DenseOperatorStyle{BL<:Basis,BR<:Basis} <: DataOperatorStyle{BL,BR} end

# Style precedence rules
Broadcast.BroadcastStyle(::Type{<:DenseOperator{BL,BR}}) where {BL<:Basis,BR<:Basis} = DenseOperatorStyle{BL,BR}()
Broadcast.BroadcastStyle(::DenseOperatorStyle{B1,B2}, ::DenseOperatorStyle{B3,B4}) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(IncompatibleBases())

# Out-of-place broadcasting
@inline function Base.copy(bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:DenseOperatorStyle{BL,BR},Axes,F,Args<:Tuple}
    bcf = Broadcast.flatten(bc)
    bl,br = find_basis(bcf.args)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    # TODO: remove convert
    return DenseOperator{BL,BR}(bl, br, convert(Matrix{ComplexF64}, copy(bc_)))
end
find_basis(a::DataOperator, rest) = (a.basis_l, a.basis_r)

const BasicMathFunc = Union{typeof(+),typeof(-),typeof(*)}
function Broadcasted_restrict_f(f::BasicMathFunc, args::Tuple{Vararg{<:DataOperator}}, axes)
    args_ = Tuple(a.data for a=args)
    return Broadcast.Broadcasted(f, args_, axes)
end
function Broadcasted_restrict_f(f, args::Tuple{Vararg{<:DataOperator}}, axes)
    throw(error("Cannot broadcast function `$f` on type `$(eltype(args))`"))
end

# In-place broadcasting
@inline function Base.copyto!(dest::DataOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:DataOperatorStyle{BL,BR},Axes,F,Args}
    axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && isa(bc.args, Tuple{<:DataOperator{BL,BR}}) # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    # Get the underlying data fields of operators and broadcast them as arrays
    bcf = Broadcast.flatten(bc)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    copyto!(dest.data, bc_)
    return dest
end
@inline Base.copyto!(A::DataOperator{BL,BR},B::DataOperator{BL,BR}) where {BL<:Basis,BR<:Basis} = (copyto!(A.data,B.data); A)
@inline Base.copyto!(dest::DataOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:DataOperatorStyle,Axes,F,Args} =
    throw(IncompatibleBases())


text:
# https://github.com/JuliaLang/Cairo.jl/blob/master/samples/sample_arc.jl

using Cairo

c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

## original example, following here
xc = 128.0;
yc = 128.0;
radius = 100.0;
angle1 = 45.0  * (pi/180.0);  # angles are specified
angle2 = 180.0 * (pi/180.0);  # in radians

set_line_width(cr, 10.0);
arc(cr, xc, yc, radius, angle1, angle2);
stroke(cr);

# draw helping lines
set_source_rgba(cr, 1, 0.2, 0.2, 0.6);
set_line_width(cr, 6.0);

arc(cr, xc, yc, 10.0, 0, 2*pi);
fill(cr);

arc(cr, xc, yc, radius, angle1, angle1);
line_to(cr, xc, yc);
arc(cr, xc, yc, radius, angle2, angle2);
line_to(cr, xc, yc);
stroke(cr);

## mark picture with current date
move_to(cr,0.0,12.0);
set_source_rgb (cr, 0,0,0);
show_text(cr,strftime(time()));
write_to_png(c,"sample_arc.png");


text:
module GraphTheoryMips

using JuMP
using Cbc
#using Graphs
using LightGraphs

include("graphs/structural.jl")
include("graphs/flows.jl")
include("util.jl")

export

max_independent_set,
min_dominating_set,
min_vertex_cover,
distances,
bipartite_set,
bin_decisions

end


text:
module MVCalculusTest
using BenchmarkTools, LinearAlgebra, Random, Test, WiSER

@testset "kron_axpy!" begin
Random.seed!(123)
m, n = 100, 50
p, q = 30, 20
A = randn(m, n)
X = randn(p, q)

Y1 = zeros(m * p, n * q)
kron_axpy!(A, X, Y1)
@test Y1 == kron(A, X)
bm = @benchmark kron_axpy!($A, $X, $Y1) setup=(fill!($Y1, 0))
display(bm); println()
@test allocs(bm) == 0

Y2 = zeros(n * p, m * q)
kron_axpy!(transpose(A), X, Y2)
@test Y2 == kron(transpose(A), X)
bm = @benchmark kron_axpy!($(transpose(A)), $X, $Y2) setup=(fill!($Y2, 0))
display(bm); println()
@test allocs(bm) == 0

Y3 = zeros(m * q, n * p)
kron_axpy!(A, transpose(X), Y3)
@test Y3 == kron(A, transpose(X))
bm = @benchmark kron_axpy!($A, $(transpose(X)), $Y3) setup=(fill!($Y3, 0))
display(bm); println()
@test allocs(bm) == 0

Y4 = zeros(n * q, m * p)
kron_axpy!(transpose(A), transpose(X), Y4)
@test Y4 == kron(transpose(A), transpose(X))
bm = @benchmark kron_axpy!($(transpose(A)), $(transpose(X)), $Y4) setup=(fill!($Y4, 0))
display(bm); println()
@test allocs(bm) == 0
end

@testset "kr_axpy!" begin
Random.seed!(123)
m, n, p = 100, 50, 30
A = randn(m, n)
X = randn(p, n)

Y1 = zeros(m * p, n)
kr_axpy!(A, X, Y1)
@test Y1 == reshape([A[i1, j] * X[i2, j] for i2 in 1:p, i1 in 1:m, j in 1:n], m * p, n)
bm = @benchmark kr_axpy!($A, $X, $Y1) setup=(fill!($Y1, 0))
display(bm); println()
@test allocs(bm) == 0
end

@testset "CopyMatrix" begin
Random.seed!(123)
n = 50
C = CopyMatrix(n)

@test size(C) == (n * n, (n * (n + 1)) >> 1)
# display(C); println()
# display(Matrix{Float64}(C)); println()
M = randn(n, n)
@test vec(LowerTriangular(M)) == C * vech(LowerTriangular(M))
@test vec(M)' * C == vech(M)'
@test C' * vec(M) == vech(M)

result1 = similar(vec(M)' * C)
bm = @benchmark mul!($result1, $(transpose(vec(M))), $C)
display(bm); println()
@test allocs(bm) == 0

result2 = similar(C' * vec(M))
bm = @benchmark mul!($result2, $(transpose(C)), $(vec(M)))
display(bm); println()
@test_skip allocs(bm) == 0
end

@testset "Ct_At_kron_A_KC" begin
Random.seed!(123)
q = 5
q◺ = WiSER.◺(q)
A = randn(q, q); A = A'A * LowerTriangular(randn(q, q))
H = Ct_At_kron_A_KC(A)
Cq  = CopyMatrix(q)
Kqq = commutation(q, q)
@test size(H) == (q◺, q◺)
@test issymmetric(H)
@test H == Cq' * kron(A', A) * Kqq * Cq
bm = @benchmark Ct_At_kron_A_KC!($H, $A) setup=(fill!($H, 0))
display(bm); println()
@test allocs(bm) == 0
end

@testset "Ct_A_kron_B_C" begin
Random.seed!(123)
q = 5
q◺ = WiSER.◺(q)
A = randn(q, q); A = A'A
B = randn(q, q); B = B'B
H = Ct_A_kron_B_C(A, B)
Cq  = CopyMatrix(q)
@test size(H) == (q◺, q◺)
@test issymmetric(H)
@test all(eigvals(Symmetric(H)) .≥ 0)
@test H == Cq' * kron(A, B) * Cq
bm = @benchmark Ct_A_kron_B_C!($H, $A, $B) setup=(fill!($H, 0))
display(bm); println()
@test allocs(bm) == 0
end

@testset "Ct_A_kr_B" begin
Random.seed!(123)
q, n = 5, 10
q◺ = ◺(q)
A = randn(q, n)
B = randn(q, n)
H = Ct_A_kr_B(A, B)
Cq  = CopyMatrix(q)
@test size(H) == (q◺, n)
@test H == Cq' * kr(A, B)
bm = @benchmark Ct_A_kr_B!($H, $A, $B) setup=(fill!($H, 0))
display(bm); println()
@test allocs(bm) == 0
end
end


text:
using Base.Broadcast: broadcastable
using Swizzles
using Swizzles.ShallowArrays
using Swizzles.WrapperArrays

struct SimpleShallowArray{T, N, P} <: ShallowArray{T, N, P}
    arg::P
end

function SimpleShallowArray(arg)
    arg = broadcastable(arg)
    SimpleShallowArray{eltype(arg), ndims(typeof(arg)), typeof(arg)}(arg)
end
SimpleShallowArray(arg::Tuple) = SimpleShallowArray{eltype(arg), 1, typeof(arg)}(arg)
Base.parent(arr::SimpleShallowArray) = arr.arg
WrapperArrays.adopt(arr::SimpleShallowArray, arg) = SimpleShallowArray(arg)
Base.size(arr::SimpleShallowArray{<:Any, <:Any, <:Tuple}) = length(arr.arg)
Base.axes(arr::SimpleShallowArray{<:Any, <:Any, <:Tuple}) = (Base.OneTo(length(arr.arg)),)

foo(x) = x

@testset "ShallowArrays" begin
    for arg in ((1, 2, 3.0), (1, 2, 3), (), [1, 2, 3.0], [1, 2, 3], [], [1 2; 3 4], transpose([1, 2]))

        x = foo.(arg)
        y = foo.(SimpleShallowArray(arg))
        @test x == y
        @test typeof(x) == typeof(y)

        x = arg .+ arg
        y = SimpleShallowArray(arg) .+ arg
        @test x == y
        @test typeof(x) == typeof(y)

        x = arg .+ arg
        y = SimpleShallowArray(arg) .+ SimpleShallowArray(arg)
        @test x == y
        @test typeof(x) == typeof(y)
    end
end


text:
using Caesar, Caesar.ZmqCaesar
using Test

@testset "ZMQ Callback Compatibility Tests" begin

    # Bring in the callbacks symbols
    include(joinpath(dirname(@__FILE__), "..", "..", "src", "zmq", "services", "ZmqServer.jl"))

    @testset "ZMQ Callback Compatibility Tests: System Callbacks" begin
        for callback in systemverbs
            @test typeof(getfield(Caesar.ZmqCaesar, callback)) <: Function
        end
    end

    @testset "ZMQ Callback Compatibility Tests: Config Callbacks" begin
        for callback in configverbs
            @test typeof(getfield(Caesar.ZmqCaesar, callback)) <: Function
        end
    end

    @testset "ZMQ Callback Compatibility Tests: Session Callbacks" begin
        for callback in sessionverbs
            @test typeof(getfield(Caesar.ZmqCaesar, callback)) <: Function
        end
    end
end


text:
# io.jl
### Files for saving and reading mh results

abstract type Fields end

struct Result_Fields <: Fields
    name::String
    label::String
    cond::String
    gene::String
    model::String
    nalleles::String
end

struct Summary_Fields <: Fields
    name::String
    label::String
    cond::String
    model::String
end

const raterow_dict = Dict([("ml", 1),("mean",2),("median",3),("last",4)])
const statrow_dict = Dict([("mean",1),("SD", 2),("median",3),("MAD",4)])

"""
  write_dataframes(resultfolder::String,datafolder::String,measure=:AIC)
  write_dataframes(resultfolder::String,measure)

  collates run results into a csv file

Arguments
- `resultfolder`: name of folder with result files
- `datafolder`: name of folder where data is stored
- `measure`: measure used to assess winner
- `assemble`: if true then assemble results into summary files
"""

function write_dataframes(resultfolder::String,datafolder::String;measure::Symbol=:AIC,assemble::Bool=true)
    write_dataframes_only(resultfolder,datafolder,assemble=assemble)
    write_winners(resultfolder,measure)
end

function write_dataframes_only(resultfolder::String,datafolder::String;assemble::Bool=true)
    dfs = make_dataframes(resultfolder,datafolder,assemble)
    for df in dfs
        for dff in dfs
            for dfff in dff
                csvfile = joinpath(resultfolder,dfff[1])
                CSV.write(csvfile,dfff[2])
            end
        end
    end
    nothing
end

"""
write_winners(resultfolder,measure)

Write best performing model for measure

"""
function write_winners(resultfolder,measure)
    df = best_measure(resultfolder,measure)
    for i in eachindex(df)
        csvfile = joinpath(resultfolder,df[i][1])
        CSV.write(csvfile,df[i][2])
    end
    nothing
end

"""
write_augmented(summaryfile::String,resultfolder,datafolder;fishdata=false)

Augment summary file with G=2 burst size, model predicted moments, and fit measures

"""
write_augmented(summaryfile::String,resultfolder,datafolder;fishdata=false) = CSV.write(summaryfile,augment_dataframe(read_dataframe(summaryfile),resultfolder,datafolder))

"""
read_dataframe(csvfile::String)
"""
read_dataframe(csvfile::String) = DataFrame(CSV.File(csvfile))

"""
get_suffix(file::String)

"""
get_suffix(file::String) = chop(file,tail=4), last(file,3)

# does not account for csv files with less than 4 fields
function fields(file::String)
    file,suffix = get_suffix(file)
    v = split(file,"_")
    if suffix == "csv"
        if length(v) == 4
            s = Summary_Fields(v[1],v[2],v[3],v[4])
        else
            println(file)
            throw("Incorrect file name format")
        end
    else
        if length(v) == 6
            s = Result_Fields(v[1],v[2],v[3],v[4],v[5],v[6])
        else
            println(file)
            throw("Incorrect file name format")
        end
    end
    return s
end

function isratefile(folder::String)
    files=readdir(folder)
    any(occursin.(".csv",files) .& occursin.("rates",files))
end

isfish(string::String) = occursin("FISH",string)

function get_genes(file::String)
    r,header = readdlm(file,',',header=true)
    return r[:,1]
end

get_genes(root,cond,datafolder) = get_genes(cond,joinpath(root,datafolder))

function get_genes(cond,datafolder)
    genes = Vector{String}(undef,0)
    files = readdir(datafolder)
    for file in files
        if occursin(cond,file)
            push!(genes,split(file,"_")[1])
        end
    end
    return genes
end
"""
    get_genes(folder,type,label,cond,model)

"""
# not working
function get_genes(folder,type,label,cond,model)
    genes = Array{String,1}(undef,0)
    files = get_files(folder,type,label,cond,model)
    for file in files
        push!(genes,get_gene(file))
    end
    return genes
end

get_files(folder::String,resultname,label,cond,model) = get_files(get_resultfiles(folder),resultname,label,cond,model)

function get_files(files::Vector,resultname,label,cond,model)
    parts = fields.(files)
    files[(getfield.(parts,:name) .== resultname) .& (getfield.(parts,:label) .== label) .& (getfield.(parts,:cond) .== cond) .& (getfield.(parts,:model) .== model)]
end

get_gene(file::String) = fields(file).gene
get_model(file::String) = fields(file).model
get_label(file::String) = fields(file).label
get_cond(file::String) = fields(file).cond
get_nalleles(file::String) = fields(file).nalleles

get_fields(parts::Vector{T}, field::Symbol) where T <: Fields = unique(getfield.(parts,field))

get_models(parts::Vector{T}) where T <: Fields = get_fields(parts,:model)

get_genes(parts::Vector{T}) where T <: Fields = get_fields(parts,:gene)

get_conds(parts::Vector{T}) where T <: Fields = get_fields(parts,:cond)

get_labels(parts::Vector{T}) where T <: Fields = get_fields(parts,:label)

get_names(parts::Vector{T}) where T <: Fields = get_fields(parts,:name)

get_nalleles(parts::Vector{T}) where T <: Fields = get_fields(parts,:nalleles)

get_resultfiles(folder::String) = get_resultfiles(readdir(folder))
get_resultfiles(files::Vector) = files[occursin.(".txt",files) .& occursin.("_",files)]

get_summaryfiles(folder::String) = get_summaryfiles(readdir(folder))
get_summaryfiles(files::Vector) = files[occursin.(".csv",files) .& occursin.("_",files)]
get_summaryfiles(files::Vector,name) = files[occursin.(".csv",files) .& occursin.(name,files)]

get_ratesummaryfiles(files::Vector) = get_summaryfiles(files,"rates")
get_ratesummaryfiles(folder::String) = get_ratesummaryfiles(get_summaryfiles(folder))

get_measuresummaryfiles(files::Vector) = get_summaryfiles(files,"measures")
get_measuresummaryfiles(folder::String) = get_measuresummaryfiles(get_summaryfiles(folder))

get_burstsummaryfiles(files::Vector) = get_summaryfiles(files,"burst")
get_burstsummaryfiles(folder::String) = get_burstsummaryfiles(get_summaryfiles(folder))

"""
write_moments(outfile,genelist,cond,datafolder,fish,root)

"""
function write_moments(outfile,genelist,cond,datafolder,fish,root)
    f = open(outfile,"w")
    writedlm(f,["Gene" "Expression Mean" "Expression Variance"],',')
    for gene in genelist
        h = get_histogram_rna(gene,cond,datafolder,fish,root)
        writedlm(f,[gene mean_histogram(h) var_histogram(h)],',')
    end
    close(f)
end

"""
    write_burst_stats(outfile,infile::String,G::String,cell,folder,cond,root)

"""
function write_burst_stats(outfile,infile::String,G::String,cell,folder,cond,root)
    folder = joinpath(root,folder)
    condarray = split(cond,"-")
    g = parse(Int,G)
    lr = 2*g
    lc = 2*g-1
    freq = Array{Float64,1}(undef,2*length(condarray))
    burst = similar(freq)
    f = open(joinpath(folder,outfile),"w")
    contents,head = readdlm(joinpath(folder,infile),',',header=true)
    label = Array{String,1}(undef,0)
    for c in condarray
        label = vcat(label, "Freq " * c, "sd","Burst Size " * c, "sd")
    end
    writedlm(f,["gene" reshape(label,1,length(label))],',')
    for r in eachrow(contents)
        gene = String(r[1])
        rates = r[2:end]
        rdecay = decay(root,cell,gene)
        cov = read_covparam(joinpath(folder,getfile("param-stats",gene,G,folder,cond)[1]))
        # mu = readmean(joinpath(folder,getfile("param-stats",gene,G,folder,cond)[1]))
        if size(cov,2) < 2
            println(gene)
        end
        for i in eachindex(condarray)
            j = i-1
            freq[2*i-1], freq[2*i] = frequency(rates[1+lr*(i-1)],sqrt(cov[1+lc*j,1+lc*j]),rdecay)
            burst[2*i-1], burst[2*i] = burstsize(rates[3+lr*j],rates[2+lr*j],cov[3+lc*j,3+lc*j],cov[2+lc*j,2+lc*j],cov[2+lc*j,3+lc*j])
        end
        writedlm(f,[gene freq[1] freq[2] burst[1] burst[2] freq[3] freq[4] burst[3] burst[4]],',')
        flush(f)
    end
    close(f)
end
"""
    write_histograms(resultfolder,ratefile,cell,datacond,G::Int,datafolder::String,fish,root,outfolder = "histograms")

"""
function write_histograms(resultfolder,ratefile,cell,datacond,G::Int,datafolder::String,fish,root,outfolder = "histograms")
    ratefile = joinpath(resultfolder,ratefile)
    rates,head = readdlm(ratefile,',',header=true)
    outfolder = joinpath(resultfolder,outfolder)
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    cond = string.(split(datacond,"-"))
    for r in eachrow(rates)
        h = histograms(r,cell,cond,G,datafolder,fish,root)
        for i in eachindex(cond)
            f = open(joinpath(outfolder,string(r[1]) * cond[i] * ".txt"),"w")
            writedlm(f,h[i])
            close(f)
        end
    end
end

"""
    assemble_all(folder)

"""
function assemble_all(folder::String)
    files = get_resultfiles(folder)
    parts = fields.(files)
    labels = get_labels(parts)
    conds = get_conds(parts)
    models = get_models(parts)
    assemble_all(folder,files,labels,conds,models)
end

function assemble_all(folder::String,files::Vector,labels::Vector,conds::Vector,models::Vector)
    for l in labels, c in conds, g in models
        assemble_all(folder,files,l,c,g,isfish(l))
    end
end

function assemble_all(folder::String,files::Vector,label::String,cond::String,model::String,fish::Bool)
    assemble_rates(folder,files,label,cond,model,fish)
    assemble_measures(folder,files,label,cond,model)
    assemble_stat(folder,files,label,cond,model)
    if model == "2"
        assemble_burst_model2(folder,files,label,cond)
    end
end

function assemble_files(folder::String,files::Vector,outfile::String,header,readfunction)
    f = open(outfile,"w")
    writedlm(f,header,',')
    for file in files
        gene = get_gene(file)
        r = readfunction(joinpath(folder, file))
        writedlm(f,[gene r],',')
    end
    close(f)
end

function assemble_rates(folder::String,files::Vector,label::String,cond::String,model::String,fish)
    outfile = joinpath(folder,"rates_" * label * "_" * cond * "_" * model * ".csv")
    assemble_files(folder,get_files(files,"rates",label,cond,model),outfile,ratelabels(model,length(split(cond,"-"))),readmean)
end

function assemble_measures(folder::String,files,label::String,cond::String,model::String)
    outfile = joinpath(folder,"measures_" * label * "_" * cond * "_" * model *  ".csv")
    header = ["Gene" "Nalleles" "Deviance" "LogMaxLikelihood" "WAIC" "SD" "AIC" "Acceptance" "Temperature" "Rhat"]
    # assemble_files(folder,get_files(files,"measures",label,cond,model),outfile,header,readmeasures)
    files = get_files(files,"measures",label,cond,model)
    f = open(outfile,"w")
    writedlm(f,header,',')
    for file in files
        gene = get_gene(file)
        nalleles = get_nalleles(file)
        r = readmeasures(joinpath(folder, file))
        writedlm(f,[gene nalleles r],',')
    end
    close(f)
end

function assemble_mad(folder::String,files,label::String,cond::String,model::String)
    outfile = joinpath(folder,"mad_" * label * "_" * cond * "_" * model *  ".csv")
    assemble_files(folder,get_files(files,"param-stats",label,cond,model),outfile,ratelabels(model,length(split(cond,"-")))[:,1:end-1],readmad)
end

function assemble_sd(folder::String,files,label::String,cond::String,model::String)
    outfile = joinpath(folder,"sd_" * label * "_" * cond * "_" * model *  ".csv")
    assemble_files(folder,get_files(files,"param-stats",label,cond,model),outfile,ratelabels(model,length(split(cond,"-")),"SD")[:,1:end-1],readsd)
end

function assemble_stat(folder::String,files,label::String,cond::String,model::String)
    outfile = joinpath(folder,"stats_" * label * "_" * cond * "_" * model *  ".csv")
    assemble_files(folder,get_files(files,"param-stats",label,cond,model),outfile,statlabels(model,length(split(cond,"-"))),readstats)
end

function assemble_burst_model2(folder::String,files,label::String,cond::String,model::String = "2")
    outfile = joinpath(folder,"burst_" * label * "_" * cond * "_" * model *  ".csv")
    assemble_files(folder,get_files(files,"param-stats",label,cond,model),outfile,["Gene" "BurstSize" "BurstSD" "BurstVar" "E" "O" "VEE" "VOO" "VEO"],read_burst_model2)
end

function Gratelabels(model,nsets)
    G = parse(Int,model)
    n = G-1
    Grates = Array{String,2}(undef,1,2*n)
    for i = 0:n-1
        Grates[1,2*i+1] = "Rate$i$(i+1)"
        Grates[1,2*i+2] = "Rate$(i+1)$i"
    end
    return Grates
end

function ratelabels(model,nsets)
    Grates = Gratelabels(model,nsets)
    rates = [Grates "Eject" "Decay"]
    for i = 2:nsets
        rates = [rates rates]
    end
    return ["Gene" rates]
end

function statlabels(model,nsets)
    label = ["SD","Median","MAD"]
    Grates = Gratelabels(model,nsets)
    Grates = [Grates "Eject"]
    for i = 2:nsets
        Grates = [Grates Grates]
    end
    rates = Matrix{String}(undef,1,0)
    for i in 1:3
        rates = [rates Grates .* label[i]]
    end
    return ["Gene" rates]
end

# function ratelabels(model,nsets,sd::Bool=false)
#     G = parse(Int,model)
#     n = G-1
#     Grates = Array{String,2}(undef,1,2*n)
#     for i = 0:n-1
#         Grates[1,2*i+1] = "Rate$i$(i+1)"
#         Grates[1,2*i+2] = "Rate$(i+1)$i"
#     end
#     rates = [Grates "Eject" "Decay"]
#     for i = 2:nsets
#         rates = [rates Grates "Eject" "Decay"]
#     end
#     if typeof(model) <: AbstractGMlossmodel
#         rates = [rates "Yield"]
#         lenr = 2*G*nsets + 1
#     else
#         lenr = 2*G*nsets
#     end
#     if sd
#         rates = reshape([rates; fill("SD",1,lenr)],1,2*(lenr))
#     end
#     return ["Gene" rates]
# end

function get_all_rates(file::String,header::Bool)
    r = readdlm(file,',',header=header)
    if header
        r = r[1]
    end
    return r
end

filename(data,model::AbstractGRMmodel) = filename(data.name,data.gene,model.G,model.R,model.nalleles)
filename(data,model::AbstractGMmodel) = filename(data.name,data.gene,model.G,model.nalleles)
filename(label::String,gene::String,G::Int,R::Int,nalleles::Int) = filename(label,gene,"$G"*"$R","$(nalleles)")
filename(label::String,gene,G::Int,nalleles::Int) = filename(label,gene,"$G","$(nalleles)")
filename(label::String,gene::String,model::String,nalleles::String) = "_" * label  * "_" * gene *  "_" * model * "_" * nalleles * txtstr

"""
write_results(file::String,x)
"""
function writeall(path::String,fit,stats,measures,data,temp,model::StochasticGRmodel)
    if ~isdir(path)
        mkpath(path)
    end
    name = filename(data,model)
    write_rates(joinpath(path,"rates" * name ),fit,stats,model)
    write_measures(joinpath(path,"measures" * name),fit,measures,deviance(fit,data,model),temp)
    write_param_stats(joinpath(path,"param-stats" * name),stats)

end

"""
write_rates(file::String,fit)

Write rate parameters, rows in order are
maximum likelihood
mean
median
last accepted
"""
function write_rates(file::String,fit::Fit,stats,model)
    f = open(file,"w")
    writedlm(f,[get_rates(fit.parml,model)],',')
    writedlm(f,[get_rates(stats.meanparam,model)],',')
    writedlm(f,[get_rates(stats.medparam,model)],',')
    writedlm(f,[get_rates(fit.param[:,end],model)],',')
    close(f)
end
"""
write_measures(file,fit,waic,dev)
"""
function write_measures(file::String,fit::Fit,measures::Measures,dev,temp)
    f = open(file,"w")
    writedlm(f,[fit.llml mean(fit.ll) std(fit.ll) quantile(fit.ll,[.025;.5;.975])' measures.waic[1] measures.waic[2] aic(fit)],',')
    writedlm(f,dev,',')
    writedlm(f,[fit.accept fit.total],',')
    writedlm(f,temp,',')
    writedlm(f,measures.rhat',',')
    writedlm(f,maximum(measures.rhat),',')
    close(f)
end
"""
write_param_stats(stats,waic,data,model)

"""
function write_param_stats(file,stats::Stats)
    f = open(file,"w")
    writedlm(f,stats.meanparam',',')
    writedlm(f,stats.stdparam',',')
    writedlm(f,stats.medparam',',')
    writedlm(f,stats.madparam',',')
    writedlm(f,stats.qparam,',')
    writedlm(f,stats.corparam,',')
    writedlm(f,stats.covparam,',')
    writedlm(f,stats.covlogparam,',')
    close(f)
end

"""
readrates(file::String,row::Int)

row
1       maximum likelihood
2       mean
3       median
4       last value of previous run
"""

readrates(file::String) = readrates(file,3)
readrates(file::String,row::Int) = readrow(file,row)

function readrow(file::String,row)
    if isfile(file) && ~isempty(read(file))
        contents = readdlm(file,',')
        m = contents[row,:]
        return m[.~isempty.(m)]
    else
        println(file)
    end
end

function readrow_flip(file,row)
    m = readrow(file,row)
    reshape(m,1,length(m))
end

function readmeasures(file::String)
    d = readdeviance(file)
    w = readwaic(file)
    a = readaccept(file)
    t = readtemp(file)
    r = readrhat(file)
    [d[1] w[1] w[7] w[8] w[9] a t[1] r[1]]
end

readdeviance(file::String) = readrow(file,2)

readwaic(file::String) = readrow(file,1)

function readaccept(file::String)
    a = readrow(file,3)
    a[1]/a[2]
end

readtemp(file::String) = readrow(file,4)

readrhat(file::String) = readrow(file,6)

function readmean(ratefile::String)
    m = readrow(ratefile,2)
    reshape(m,1,length(m))
end

function readsd(ratefile::String)
    m = readrow(ratefile,2)
    reshape(m,1,length(m))
end

function readstats(statfile::String)
    sd = readrow_flip(statfile,2)
    median = readrow_flip(statfile,3)
    mad = readrow_flip(statfile,4)
    [sd median mad]
end

function readmedian(statfile::String)
    m = readrow(statfile,3)
    reshape(m,1,length(m))
end

function readmad(statfile::String)
    m = readrow(file,4)
    reshape(m,1,length(m))
end

function read_corparam(file::String)
    c = readdlm(file,',')
    n = length(c[1,:])
    # c[5+n:4+2*n,1:n]
    c[8:7+n,1:n]
end

function read_covparam(file::String)
    c = readdlm(file,',')
    read_covparam(c)
end

function read_covparam(c::Matrix)
    n = length(c[1,:])
    # c[5+2*n:4+3*n,1:n]
    c[8+n:7+2*n,1:n]
end

function read_covlogparam(file::String)
    c = readdlm(file,',')
    n = length(c[1,:])
    c[8+2*n:7+3*n,1:n]
end

read_crosscov(statfile::String) = read_crosscov(read_covparam(statfile))

function read_crosscov(C::Matrix)
    c = Float64[]
    N = size(C,1)
    for i in 1:N
        for j in i+1:N
            push!(c,C[i,j])
        end
    end
    c
end

function read_burst_model2(file::String)
    c = readdlm(file,',')
    b = c[1,end]/c[1,end-1]
    cov = read_covparam(c)
    v = var_ratio(c[1,end],c[1,end-1],cov[end,end],cov[end-1,end-1],cov[end-1,end])
    [b sqrt(abs(v)) v c[1,end] c[1,end-1] cov[end,end] cov[end-1,end-1] cov[end-1,end]]
end

function write_residency_G(fileout::String,filein::String,G,header)
    rates = get_all_rates(filein,header)
    n = G-1
    f = open(fileout,"w")
    writedlm(f,["gene" collect(0:n)'],',')
    for r in eachrow(rates)
        writedlm(f,[r[1] residenceprob_G(r[2:2*n+1],n)],',')
    end
    close(f)
end

# Functions for saving and loading data and models

"""
write_log(file,datafile,data,model)
write all information necessary for rerunning
"""
function save_data(file::String,data::TransientRNAData)
    f = open(file,"w")
    writedlm(f, [typeof(data)])
    writedlm(f,[data.name])
    writedlm(f,[data.gene])
    writedlm(f,[data.nRNA])
    writedlm(f,[data.time])
    writedlm(f,[data.histRNA])
    close(f)
end

function load_data(file::String,model::AbstractGMmodel)

end

function save_model(file::String,model::AbstractGMmodel)
    f = open(file,"w")
    write(f, model.G)
    writedlm(f,model.nalleles)
    writedlm(f,model.ratepriors)
    writedlm(f,model.proposal)
    writedlm(f,model.fittedparam)
    writedlm(f,model.method)
    close(f)

end

function load_model(file::String,model::AbstractGRMmodel)

end

"""
    makestring(v::Vector)

    turn a vector of strings into a single string
"""
function makestring(v)
    s = ""
    for i in v
        s *= i
    end
    return s
end


text:
type Bus
   bindex::Int
   btype::Int
   Pd::Float64
   Qd::Float64
   Gs::Float64
   Bs::Float64
   area::Int
   Vm::Float64
   Va::Float64
   baseKV::Float64
   bzone::Int
   Vmax::Float64
   Vmin::Float64
   children::Vector{Int}
   ancestor::Vector{Int}
   inline::Vector{Int}
   outline::Vector{Int}
   function Bus(bindex, btype, Pd, Qd, Gs, Bs, area, Vm, Va, baseKV, bzone, Vmax, Vmin)
      b = new(bindex, btype)
      b.Pd = Pd
      b.Qd = Qd
      b.Gs = Gs
      b.Bs = Bs
      b.area = area
      b.Vm = Vm
      b.Va = Va
      b.baseKV = baseKV
      b.bzone = bzone
      b.Vmax = Vmax
      b.Vmin = Vmin
      b.children = Int[]
      b.ancestor = Int[]
      b.inline = Int[]
      b.outline = Int[]
      return b
   end
end

type Line
   lindex::Int
   fbus::Int # the "from" node
   tbus::Int # the "to" node
   r::Float64 # the resistance value
   x::Float64 # the reactance value
   b::Float64 # the susceptance value
   u::Float64 # the capacity of the line
   children::Vector{Int}
   ancestor::Vector{Int}
   function Line(lindex, fbus, tbus, r, x, b, u)
      line = new(lindex, fbus, tbus)
      line.r = r
      line.x = x
      line.b = b
      line.u = u
      line.children = Int[]
      line.ancestor = Int[]
      return line
   end
end

type Generator
   gindex::Int
   gtype::String
   location::Int
   Pg::Float64
   Qg::Float64
   Qmax::Float64
   Qmin::Float64
   Vg::Float64
   mBase::Float64
   status::Int
   Pmax::Float64
   Pmin::Float64
   cost::Vector{Float64}
   SUcost::Float64
   SDcost::Float64
   RU::Float64
   RD::Float64
   UPtime::Float64
   DNtime::Float64
   function Generator(gindex, gtype, location, Pg, Qg, Qmax, Qmin, Vg, mBase, status, Pmax, Pmin, cost, SUcost, SDcost, RU, RD, UPtime, DNtime)
      g = new(gindex, gtype, location)
      g.Pg = Pg
      g.Qg = Qg
      g.Qmax = Qmax
      g.Qmin = Qmin
      g.Vg = Vg
      g.mBase = mBase
      g.status = status
      g.Pmax = Pmax
      g.Pmin = Pmin
      g.cost = cost
      g.SUcost = SUcost
      g.SDcost = SDcost
      g.RU = RU
      g.RD = RD
      g.UPtime = UPtime
      g.DNtime = DNtime
      return g
   end
end


text:
using DifferentialEquations, Flux, DiffEqFlux
using Plots, NPZ

true_data = npzread("data.npy")[1:20:4000,:]
data_len = size(true_data)[1]
x_data = true_data[:,1]
state_data = true_data[:,2:3]

u0 = [true_data[1,2],true_data[1,3]]
tspan = (0.0,true_data[data_len,1])
tsteps = 0.0:true_data[data_len,1]-true_data[data_len-1,1]:true_data[data_len,1]

my_nn = Chain(
  Dense(2, 32, relu),
  Dense(32, 2)) |> f32

p, re = Flux.destructure(my_nn)

function right_hand_side(du,u,p,t)
  m = re(p)
  nn_output = m(u)
  du[1] = nn_output[1]
  du[2] = nn_output[2]
end

prob = ODEProblem(right_hand_side,u0,tspan,p)
sol = solve(prob, Tsit5())

# Plot the solution
plot(sol)
savefig("Initial_NN_ODE.png")

function predict_ode(p)
  tmp_prob = remake(prob, p = p)
  pred = Array(solve(tmp_prob, Tsit5(), saveat = tsteps))
  return pred
end

function loss_ode(p)
    pred = predict_ode(p)
    loss = sum(abs2, transpose(state_data)-pred)
    return loss, pred
end

callback = function (p, l, pred)
  display(l)
  plt = plot(x_data,hcat(transpose(pred),state_data),
    line=(4, [:solid :solid :dash :dash]), ylim = (-6, 6),
    label=["y1-pred" "y2-pred" "y1-true" "y2-true"])
  display(plt)
  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  return false
end

result_ode = DiffEqFlux.sciml_train(loss_ode, p,
                                    cb = callback,
                                    maxiters = 100)

prob_final = ODEProblem(right_hand_side, u0, tspan, result_ode)
sol_final = solve(prob_final,Tsit5())
plot(sol_final)
savefig("NN_ode_optimized.png") # Both state components close to 1

max_stars_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_issues_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_forks_repo_name: Romit-Maulik/Tutorials-Demos-Practice


text:
import Base.==
function ==(i1::Index, i2::Index)
    i1.meta == i2.meta
end

function buff(size=3000)
    io = IOBuffer()
    write(io, [UInt8(1) for _ in 1:size])
    seek(io, 0)
    io
end

function rw_leafnode()
    io = buff()
    ni = Index1024.NodeInfo(Index1024.tag(Index1024.leaf, 0x01), Index1024.Leaf(1, 0))
    @assert write(io, ni) == 24
    seek(io, 0)
    ni == read(io, Index1024.NodeInfo)
end

function rw_empty_index()
    io = buff()
    idx = Index(io)
    @assert write(io, idx) > 0
    seek(io,0)
    idx == read(io, Index)
end

function rw_index1p()
    io = buff()
    idx = Index(io)
    push!(idx.meta, "file1")
    push!(idx.meta, "file2")
    push!(idx.meta, "file3")
    @assert write(io, idx) > 0
    seek(io,0)
    idx == read(io, Index)
end

function rw_index()
    io = buff()
    idx = Index(io)
    for i in 1:256
        push!(idx.meta, "file$i")
    end
    @assert write(io, idx) > 0
    seek(io,0)
    idx == read(io, Index)
end

function nb(n=0)
    io = buff()
    if n > 0
        write(io, [1 for _ in 0:n])
    end
    Index1024.nextblock(io)
    position(io)
end

kval(k) = UInt64(k+mod(k,7))
kdata(k) = UInt64(10k)
kaux(k) = UInt64(1000*(1+mod(k,3)))
ktup(k) = (data=kdata(k), aux=kaux(k))

egtree_entries(n=16) = Dict([kval(k)=>(data=kdata(k), aux=kaux(k)) for k in 1:n])

function egtree(n=16; padding=0)
    io = buff()
    print(io, repeat("p", padding))
    tree_pos = position(io)
    entries = egtree_entries(n)
    build_index_file(io, entries)
    open_index(io)
end

tsearch_exists(n; padding=0) = search(egtree(n; padding), kval(n))
tsearch_fails(n) = search(egtree(n), 0)

tget_exists(n; padding=0) = get(egtree(n; padding), kval(n), 0)
tget_fails(n) = get(egtree(n), 0x0, 0)

function tmeta(n=16)
    io = buff()
    entries = egtree_entries(n)
    meta = String["one", "two", "three"]
    build_index_file(io, entries; meta)
    seekstart(io)
    idx = open_index(io)
    search(idx, kval(n)) == ktup(n)
    idx.meta == meta
end


text:
using QuantumOpticsBase
using QuantumOpticsBase: check_samebases, check_multiplicable
import ..timeevolution: recast!, QO_CHECKS, pure_inference

import DiffEqCallbacks, StochasticDiffEq, OrdinaryDiffEq

const DiffArray{T} = Union{AbstractArray{T,1}, AbstractArray{T, 2}}

"""
    integrate_stoch(tspan::Vector, df::Function, dg::Vector{Function}, x0::Vector{ComplexF64},
            state::T, dstate::T, fout::Function; kwargs...)

Integrate using StochasticDiffEq
"""
function integrate_stoch(tspan::Vector, df::Function, dg::Function, x0::Vector,
            state::T, dstate::T, fout::Function, n::Int;
            save_everystep = false, callback=nothing, saveat=tspan,
            alg::StochasticDiffEq.StochasticDiffEqAlgorithm=StochasticDiffEq.EM(),
            noise_rate_prototype = nothing,
            noise_prototype_classical = nothing,
            noise=nothing,
            ncb=nothing,
            kwargs...) where T

    function df_(dx::T, x::T, p, t) where T
        recast!(x, state)
        recast!(dx, dstate)
        df(t, state, dstate)
        recast!(dstate, dx)
    end

    function dg_(dx, x, p, t) where T
        recast!(x, state)
        dg(dx, t, state, dstate, n)
    end

    function fout_(x::Vector, t, integrator)
        recast!(x, state)
        fout(t, state)
    end

    nc = isa(noise_prototype_classical, Nothing) ? 0 : size(noise_prototype_classical)[2]
    if isa(noise, Nothing) && n > 0
        if n + nc == 1
            noise_ = StochasticDiffEq.RealWienerProcess(0.0, 0.0)
        else
            noise_ = StochasticDiffEq.RealWienerProcess!(0.0, zeros(n + nc))
        end
    else
        noise_ = noise
    end
    if isa(noise_rate_prototype, Nothing)
        if n > 1 || nc > 1 || (n > 0 && nc > 0)
            noise_rate_prototype = zeros(eltype(x0), length(x0), n + nc)
        end
    end

    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})

    out = DiffEqCallbacks.SavedValues(eltype(tspan),out_type)

    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=saveat,
                                         save_everystep=save_everystep,
                                         save_start = false)

    full_cb = OrdinaryDiffEq.CallbackSet(callback, ncb, scb)

    prob = StochasticDiffEq.SDEProblem{true}(df_, dg_, x0,(tspan[1],tspan[end]),
                    noise=noise_,
                    noise_rate_prototype=noise_rate_prototype)

    sol = StochasticDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-3,
                abstol = 1.0e-3,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=full_cb, kwargs...)

    out.t,out.saveval
end

"""
    integrate_stoch

Define fout if it was omitted.
"""
function integrate_stoch(tspan::Vector, df::Function, dg::Function, x0::Vector,
    state::T, dstate::T, ::Nothing, n::Int; kwargs...) where T
    function fout(t, state::T)
        copy(state)
    end
    integrate_stoch(tspan, df, dg, x0, state, dstate, fout, n; kwargs...)
end


text:
#######
##FEM Stochastic Heat Animation Test
#######

#Generates an animation for a solution of the heat equation
#Uses Plots.jl, requires matplotlib >=1.5
using FiniteElementDiffEq#, ImageMagick
prob = prob_femheat_stochasticbirthdeath

sol = solve(prob,FEMDiffEqHeatSemiImplicitCrankNicholson(),save_everystep=true,solver=:LU)

println("Generating Animation")
#animate(sol::FEMSolution;zlims=(0,3),cbar=false)

# Returns true if nothing error'd
true


text:
# Generating new terms in a sequence

function generate!(S::Sequence,n::Int)
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n s[k] = a[k] end
end

function generate!(S::Cumsum,n::Int)
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n s[k] = a[k]+s[k-1] end
end

function generate!(S::Cumprod,n::Int)
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n s[k] = a[k]*s[k-1] end
end

function generate!(S::Difference{T,D,1},n::Int) where {T,D}
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n s[k] = a[k+1]-a[k] end
end

function generate!(S::Difference{T,D,2},n::Int) where {T,D}
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n s[k] = a[k+2]-2a[k+1]+a[k] end
end

function generate!(S::Shift{T,D,O},n::Int) where {T,D,O}
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n s[k] = a[k+O] end
end

# Aitken's Δ² process

function generate!(S::Aitken,n::Int)
    s,a = sequence(S),generator(S)
    i = length(s)
    i ≤ n && resize!(s,n)
    for k = 1+i:n
        s[k] = a[k+2]-(a[k+2]-a[k+1])^2/(a[k+2]-2a[k+1]+a[k])
    end
end

# Wynn's ϵ- and ρ-algorithms

for trans in (:ϵ,:ρ)
    rec = Meta.parse(string(trans)*"recurrence!")
    @eval begin
        function generate!(S::$trans,n::Int)
            data,s = sequence(S),generator(S)
            aux = auxiliary(S)
            i = length(data)
            i ≤ n && resize!(data,n);resize!(aux,n)
            for k = 1+i:n
                aux[k] = s[k]
                $rec(k,aux,data)
            end
        end
    end
end

function ϵrecurrence!(k::Int,aux::Vector{T},data::Vector{T}) where T
    aux2 = zero(T)
    for j = k:-1:2
        aux1 = aux2
        aux2 = aux[j-1]
        diff = aux[j]-aux2
        aux[j-1] = aux1 + one(T)/diff
    end
    data[k] = isodd(k) ? aux[1] : aux[2]
end

function ρrecurrence!(k::Int,aux::Vector{T},data::Vector{T}) where T
    aux2 = zero(T)
    for j = k:-1:2
        aux1 = aux2
        aux2 = aux[j-1]
        diff = aux[j]-aux2
        aux[j-1] = aux1 + j/diff
    end
    data[k] = isodd(k) ? aux[1] : aux[2]
end

# Drummond, Levin, and Weniger transformations

for trans in (:Drummond,:Levin,:Weniger)
    rec = Meta.parse(string(trans)*"recurrence!")
    @eval begin
        function generate!(S::$trans{T},n::Int) where T
            ℓ,s = sequence(S),generator(S)
            N,D = numerator(S),denominator(S)
            i = length(ℓ)
            i ≤ n && resize!(ℓ,n);resize!(N,n);resize!(D,n)
            for k = 1+i:n
                ω = remainder(S,k)
                N[k],D[k] = s[k]/ω,one(T)/ω
                $rec(k,ω,N,D,ℓ)
            end
        end
    end
end

function Drummondrecurrence!(k::Int,ω::T,N::Vector{T},D::Vector{T},ℓ::Vector{T}) where T
    for j = 1:k-1
        N[k-j] = N[k-j+1]-N[k-j]
        D[k-j] = D[k-j+1]-D[k-j]
    end
    ℓ[k] = N[1]/D[1]
end

function Levinrecurrence!(k::Int,ω::T,N::Vector{T},D::Vector{T},ℓ::Vector{T}) where T
    kp1 = k+one(T)
    kdkp1 = k/kp1
    kdkp1jm2 = 1/kdkp1
    for j = 1:k-1
        cst = -(kp1-j)*kdkp1jm2/kp1
        N[k-j] = muladd(cst,N[k-j],N[k-j+1])
        D[k-j] = muladd(cst,D[k-j],D[k-j+1])
        kdkp1jm2 *= kdkp1
    end
    ℓ[k] = N[1]/D[1]
end

function Wenigerrecurrence!(k::Int,ω::T,N::Vector{T},D::Vector{T},ℓ::Vector{T}) where T
    km1 = k-one(T)
    for j = 1:k-1
        kjm1 = km1+j
        kjm2 = kjm1-1
        cst = -k*(km1/kjm1)/kjm2
        N[k-j] = muladd(cst,N[k-j],N[k-j+1])
        D[k-j] = muladd(cst,D[k-j],D[k-j+1])
    end
    ℓ[k] = N[1]/D[1]
end

max_stars_repo_name: MikaelSlevinsky/SequenceTransformations

max_issues_repo_name: MikaelSlevinsky/SequenceTransformations

max_forks_repo_name: MikaelSlevinsky/SequenceTransformations


text:
module Services

using ..AWSCorePrototype

const migration_hub = AWSCorePrototype.JSONService("mgh", "2017-05-31", "1.1", "AWSMigrationHub")
const accessanalyzer = AWSCorePrototype.RestXMLService("access-analyzer", "2019-11-01")
const acm = AWSCorePrototype.JSONService("acm", "2015-12-08", "1.1", "CertificateManager")
const acm_pca = AWSCorePrototype.JSONService("acm-pca", "2017-08-22", "1.1", "ACMPrivateCA")
const alexa_for_business = AWSCorePrototype.JSONService("a4b", "2017-11-09", "1.1", "AlexaForBusiness")
const amplify = AWSCorePrototype.RestXMLService("amplify", "2017-07-25")
const api_gateway = AWSCorePrototype.RestXMLService("apigateway", "2015-07-09")
const apigatewaymanagementapi = AWSCorePrototype.RestXMLService("execute-api", "2018-11-29")
const apigatewayv2 = AWSCorePrototype.RestXMLService("apigateway", "2018-11-29")
const appconfig = AWSCorePrototype.RestXMLService("appconfig", "2019-10-09")
const application_auto_scaling = AWSCorePrototype.JSONService("application-autoscaling", "2016-02-06", "1.1", "AnyScaleFrontendService")
const application_insights = AWSCorePrototype.JSONService("applicationinsights", "2018-11-25", "1.1", "EC2WindowsBarleyService")
const app_mesh = AWSCorePrototype.RestXMLService("appmesh", "2018-10-01")
const app_mesh = AWSCorePrototype.RestXMLService("appmesh", "2019-01-25")
const appstream = AWSCorePrototype.JSONService("appstream2", "2016-12-01", "1.1", "PhotonAdminProxyService")
const appsync = AWSCorePrototype.RestXMLService("appsync", "2017-07-25")
const athena = AWSCorePrototype.JSONService("athena", "2017-05-18", "1.1", "AmazonAthena")
const auto_scaling = AWSCorePrototype.RestXMLService("autoscaling", "2011-01-01")
const auto_scaling_plans = AWSCorePrototype.JSONService("autoscaling-plans", "2018-01-06", "1.1", "AnyScaleScalingPlannerFrontendService")
const backup = AWSCorePrototype.RestXMLService("backup", "2018-11-15")
const batch = AWSCorePrototype.RestXMLService("batch", "2016-08-10")
const budgets = AWSCorePrototype.JSONService("budgets", "2016-10-20", "1.1", "AWSBudgetServiceGateway")
const cost_explorer = AWSCorePrototype.JSONService("ce", "2017-10-25", "1.1", "AWSInsightsIndexService")
const chime = AWSCorePrototype.RestXMLService("chime", "2018-05-01")
const cloud9 = AWSCorePrototype.JSONService("cloud9", "2017-09-23", "1.1", "AWSCloud9WorkspaceManagementService")
const clouddirectory = AWSCorePrototype.RestXMLService("clouddirectory", "2016-05-10")
const clouddirectory = AWSCorePrototype.RestXMLService("clouddirectory", "2017-01-11")
const cloudformation = AWSCorePrototype.RestXMLService("cloudformation", "2010-05-15")
const cloudfront = AWSCorePrototype.RestXMLService("cloudfront", "2016-11-25")
const cloudfront = AWSCorePrototype.RestXMLService("cloudfront", "2017-03-25")
const cloudfront = AWSCorePrototype.RestXMLService("cloudfront", "2017-10-30")
const cloudfront = AWSCorePrototype.RestXMLService("cloudfront", "2018-06-18")
const cloudfront = AWSCorePrototype.RestXMLService("cloudfront", "2018-11-05")
const cloudfront = AWSCorePrototype.RestXMLService("cloudfront", "2019-03-26")
const cloudhsm = AWSCorePrototype.JSONService("cloudhsm", "2014-05-30", "1.1", "CloudHsmFrontendService")
const cloudhsm_v2 = AWSCorePrototype.JSONService("cloudhsmv2", "2017-04-28", "1.1", "BaldrApiService")
const cloudsearch = AWSCorePrototype.RestXMLService("cloudsearch", "2011-02-01")
const cloudsearch = AWSCorePrototype.RestXMLService("cloudsearch", "2013-01-01")
const cloudsearch_domain = AWSCorePrototype.RestXMLService("cloudsearchdomain", "2013-01-01")
const cloudtrail = AWSCorePrototype.JSONService("cloudtrail", "2013-11-01", "1.1", "com.amazonaws.cloudtrail.v20131101.CloudTrail_20131101")
const codebuild = AWSCorePrototype.JSONService("codebuild", "2016-10-06", "1.1", "CodeBuild_20161006")
const codecommit = AWSCorePrototype.JSONService("codecommit", "2015-04-13", "1.1", "CodeCommit_20150413")
const codedeploy = AWSCorePrototype.JSONService("codedeploy", "2014-10-06", "1.1", "CodeDeploy_20141006")
const codeguru_reviewer = AWSCorePrototype.RestXMLService("codeguru-reviewer", "2019-09-19")
const codeguruprofiler = AWSCorePrototype.RestXMLService("codeguru-profiler", "2019-07-18")
const codepipeline = AWSCorePrototype.JSONService("codepipeline", "2015-07-09", "1.1", "CodePipeline_20150709")
const codestar = AWSCorePrototype.JSONService("codestar", "2017-04-19", "1.1", "CodeStar_20170419")
const codestar_connections = AWSCorePrototype.JSONService("codestar-connections", "2019-12-01", "1.0", "com.amazonaws.codestar.connections.CodeStar_connections_20191201")
const codestar_notifications = AWSCorePrototype.RestXMLService("codestar-notifications", "2019-10-15")
const cognito_identity = AWSCorePrototype.JSONService("cognito-identity", "2014-06-30", "1.1", "AWSCognitoIdentityService")
const cognito_identity_provider = AWSCorePrototype.JSONService("cognito-idp", "2016-04-18", "1.1", "AWSCognitoIdentityProviderService")
const cognito_sync = AWSCorePrototype.RestXMLService("cognito-sync", "2014-06-30")
const comprehend = AWSCorePrototype.JSONService("comprehend", "2017-11-27", "1.1", "Comprehend_20171127")
const comprehendmedical = AWSCorePrototype.JSONService("comprehendmedical", "2018-10-30", "1.1", "ComprehendMedical_20181030")
const compute_optimizer = AWSCorePrototype.JSONService("compute-optimizer", "2019-11-01", "1.0", "ComputeOptimizerService")
const config_service = AWSCorePrototype.JSONService("config", "2014-11-12", "1.1", "StarlingDoveService")
const connect = AWSCorePrototype.RestXMLService("connect", "2017-08-08")
const connectparticipant = AWSCorePrototype.RestXMLService("participant.connect", "2018-09-07")
const cost_and_usage_report_service = AWSCorePrototype.JSONService("cur", "2017-01-06", "1.1", "AWSOrigamiServiceGatewayService")
const dataexchange = AWSCorePrototype.RestXMLService("dataexchange", "2017-07-25")
const data_pipeline = AWSCorePrototype.JSONService("datapipeline", "2012-10-29", "1.1", "DataPipeline")
const datasync = AWSCorePrototype.JSONService("datasync", "2018-11-09", "1.1", "FmrsService")
const dax = AWSCorePrototype.JSONService("dax", "2017-04-19", "1.1", "AmazonDAXV3")
const detective = AWSCorePrototype.RestXMLService("api.detective", "2018-10-26")
const device_farm = AWSCorePrototype.JSONService("devicefarm", "2015-06-23", "1.1", "DeviceFarm_20150623")
const direct_connect = AWSCorePrototype.JSONService("directconnect", "2012-10-25", "1.1", "OvertureService")
const application_discovery_service = AWSCorePrototype.JSONService("discovery", "2015-11-01", "1.1", "AWSPoseidonService_V2015_11_01")
const dlm = AWSCorePrototype.RestXMLService("dlm", "2018-01-12")
const database_migration_service = AWSCorePrototype.JSONService("dms", "2016-01-01", "1.1", "AmazonDMSv20160101")
const docdb = AWSCorePrototype.RestXMLService("rds", "2014-10-31")
const directory_service = AWSCorePrototype.JSONService("ds", "2015-04-16", "1.1", "DirectoryService_20150416")
const dynamodb = AWSCorePrototype.JSONService("dynamodb", "2011-12-05", "1.0", "DynamoDB_20111205")
const dynamodb = AWSCorePrototype.JSONService("dynamodb", "2012-08-10", "1.0", "DynamoDB_20120810")
const ebs = AWSCorePrototype.RestXMLService("ebs", "2019-11-02")
const ec2 = AWSCorePrototype.RestXMLService("ec2", "2016-11-15")
const ec2_instance_connect = AWSCorePrototype.JSONService("ec2-instance-connect", "2018-04-02", "1.1", "AWSEC2InstanceConnectService")
const ecr = AWSCorePrototype.JSONService("api.ecr", "2015-09-21", "1.1", "AmazonEC2ContainerRegistry_V20150921")
const ecs = AWSCorePrototype.JSONService("ecs", "2014-11-13", "1.1", "AmazonEC2ContainerServiceV20141113")
const eks = AWSCorePrototype.RestXMLService("eks", "2017-11-01")
const elastic_inference = AWSCorePrototype.RestXMLService("api.elastic-inference", "2017-07-25")
const elasticache = AWSCorePrototype.RestXMLService("elasticache", "2015-02-02")
const elastic_beanstalk = AWSCorePrototype.RestXMLService("elasticbeanstalk", "2010-12-01")
const efs = AWSCorePrototype.RestXMLService("elasticfilesystem", "2015-02-01")
const elastic_load_balancing = AWSCorePrototype.RestXMLService("elasticloadbalancing", "2012-06-01")
const elastic_load_balancing_v2 = AWSCorePrototype.RestXMLService("elasticloadbalancing", "2015-12-01")
const emr = AWSCorePrototype.JSONService("elasticmapreduce", "2009-03-31", "1.1", "ElasticMapReduce")
const elastic_transcoder = AWSCorePrototype.RestXMLService("elastictranscoder", "2012-09-25")
const ses = AWSCorePrototype.RestXMLService("email", "2010-12-01")
const marketplace_entitlement_service = AWSCorePrototype.JSONService("entitlement.marketplace", "2017-01-11", "1.1", "AWSMPEntitlementService")
const elasticsearch_service = AWSCorePrototype.RestXMLService("es", "2015-01-01")
const eventbridge = AWSCorePrototype.JSONService("events", "2015-10-07", "1.1", "AWSEvents")
const cloudwatch_events = AWSCorePrototype.JSONService("events", "2015-10-07", "1.1", "AWSEvents")
const firehose = AWSCorePrototype.JSONService("firehose", "2015-08-04", "1.1", "Firehose_20150804")
const fms = AWSCorePrototype.JSONService("fms", "2018-01-01", "1.1", "AWSFMS_20180101")
const forecast = AWSCorePrototype.JSONService("forecast", "2018-06-26", "1.1", "AmazonForecast")
const forecastquery = AWSCorePrototype.JSONService("forecastquery", "2018-06-26", "1.1", "AmazonForecastRuntime")
const frauddetector = AWSCorePrototype.JSONService("frauddetector", "2019-11-15", "1.1", "AWSHawksNestServiceFacade")
const fsx = AWSCorePrototype.JSONService("fsx", "2018-03-01", "1.1", "AWSSimbaAPIService_v20180301")
const gamelift = AWSCorePrototype.JSONService("gamelift", "2015-10-01", "1.1", "GameLift")
const glacier = AWSCorePrototype.RestXMLService("glacier", "2012-06-01")
const global_accelerator = AWSCorePrototype.JSONService("globalaccelerator", "2018-08-08", "1.1", "GlobalAccelerator_V20180706")
const glue = AWSCorePrototype.JSONService("glue", "2017-03-31", "1.1", "AWSGlue")
const greengrass = AWSCorePrototype.RestXMLService("greengrass", "2017-06-07")
const groundstation = AWSCorePrototype.RestXMLService("groundstation", "2019-05-23")
const guardduty = AWSCorePrototype.RestXMLService("guardduty", "2017-11-28")
const health = AWSCorePrototype.JSONService("health", "2016-08-04", "1.1", "AWSHealth_20160804")
const iam = AWSCorePrototype.RestXMLService("iam", "2010-05-08")
const imagebuilder = AWSCorePrototype.RestXMLService("imagebuilder", "2019-12-02")
const importexport = AWSCorePrototype.RestXMLService("importexport", "2010-06-01")
const inspector = AWSCorePrototype.JSONService("inspector", "2016-02-16", "1.1", "InspectorService")
const iot = AWSCorePrototype.RestXMLService("iot", "2015-05-28")
const iot_data_plane = AWSCorePrototype.RestXMLService("data.iot", "2015-05-28")
const iot_jobs_data_plane = AWSCorePrototype.RestXMLService("data.jobs.iot", "2017-09-29")
const iot_1click_devices_service = AWSCorePrototype.RestXMLService("devices.iot1click", "2018-05-14")
const iot_1click_projects = AWSCorePrototype.RestXMLService("projects.iot1click", "2018-05-14")
const iotanalytics = AWSCorePrototype.RestXMLService("iotanalytics", "2017-11-27")
const iot_events = AWSCorePrototype.RestXMLService("iotevents", "2018-07-27")
const iot_events_data = AWSCorePrototype.RestXMLService("data.iotevents", "2018-10-23")
const iotsecuretunneling = AWSCorePrototype.JSONService("api.tunneling.iot", "2018-10-05", "1.1", "IoTSecuredTunneling")
const iotthingsgraph = AWSCorePrototype.JSONService("iotthingsgraph", "2018-09-06", "1.1", "IotThingsGraphFrontEndService")
const kafka = AWSCorePrototype.RestXMLService("kafka", "2018-11-14")
const kendra = AWSCorePrototype.JSONService("kendra", "2019-02-03", "1.1", "AWSKendraFrontendService")
const kinesis = AWSCorePrototype.JSONService("kinesis", "2013-12-02", "1.1", "Kinesis_20131202")
const kinesis_video_archived_media = AWSCorePrototype.RestXMLService("kinesisvideo", "2017-09-30")
const kinesis_video_media = AWSCorePrototype.RestXMLService("kinesisvideo", "2017-09-30")
const kinesis_video_signaling = AWSCorePrototype.RestXMLService("kinesisvideo", "2019-12-04")
const kinesis_analytics = AWSCorePrototype.JSONService("kinesisanalytics", "2015-08-14", "1.1", "KinesisAnalytics_20150814")
const kinesis_analytics_v2 = AWSCorePrototype.JSONService("kinesisanalytics", "2018-05-23", "1.1", "KinesisAnalytics_20180523")
const kinesis_video = AWSCorePrototype.RestXMLService("kinesisvideo", "2017-09-30")
const kms = AWSCorePrototype.JSONService("kms", "2014-11-01", "1.1", "TrentService")
const lakeformation = AWSCorePrototype.JSONService("lakeformation", "2017-03-31", "1.1", "AWSLakeFormation")
const lambda = AWSCorePrototype.RestXMLService("lambda", "2014-11-11")
const lambda = AWSCorePrototype.RestXMLService("lambda", "2015-03-31")
const lex_model_building_service = AWSCorePrototype.RestXMLService("models.lex", "2017-04-19")
const license_manager = AWSCorePrototype.JSONService("license-manager", "2018-08-01", "1.1", "AWSLicenseManager")
const lightsail = AWSCorePrototype.JSONService("lightsail", "2016-11-28", "1.1", "Lightsail_20161128")
const cloudwatch_logs = AWSCorePrototype.JSONService("logs", "2014-03-28", "1.1", "Logs_20140328")
const machine_learning = AWSCorePrototype.JSONService("machinelearning", "2014-12-12", "1.1", "AmazonML_20141212")
const macie = AWSCorePrototype.JSONService("macie", "2017-12-19", "1.1", "MacieService")
const managedblockchain = AWSCorePrototype.RestXMLService("managedblockchain", "2018-09-24")
const marketplace_catalog = AWSCorePrototype.RestXMLService("catalog.marketplace", "2018-09-17")
const marketplace_commerce_analytics = AWSCorePrototype.JSONService("marketplacecommerceanalytics", "2015-07-01", "1.1", "MarketplaceCommerceAnalytics20150701")
const mediaconnect = AWSCorePrototype.RestXMLService("mediaconnect", "2018-11-14")
const mediaconvert = AWSCorePrototype.RestXMLService("mediaconvert", "2017-08-29")
const medialive = AWSCorePrototype.RestXMLService("medialive", "2017-10-14")
const mediapackage = AWSCorePrototype.RestXMLService("mediapackage", "2017-10-12")
const mediapackage_vod = AWSCorePrototype.RestXMLService("mediapackage-vod", "2018-11-07")
const mediastore = AWSCorePrototype.JSONService("mediastore", "2017-09-01", "1.1", "MediaStore_20170901")
const mediastore_data = AWSCorePrototype.RestXMLService("data.mediastore", "2017-09-01")
const mediatailor = AWSCorePrototype.RestXMLService("api.mediatailor", "2018-04-23")
const marketplace_metering = AWSCorePrototype.JSONService("metering.marketplace", "2016-01-14", "1.1", "AWSMPMeteringService")
const migrationhub_config = AWSCorePrototype.JSONService("migrationhub-config", "2019-06-30", "1.1", "AWSMigrationHubMultiAccountService")
const mobile = AWSCorePrototype.RestXMLService("mobile", "2017-07-01")
const mobile_analytics = AWSCorePrototype.RestXMLService("mobileanalytics", "2014-06-05")
const cloudwatch = AWSCorePrototype.RestXMLService("monitoring", "2010-08-01")
const mq = AWSCorePrototype.RestXMLService("mq", "2017-11-27")
const mturk = AWSCorePrototype.JSONService("mturk-requester", "2017-01-17", "1.1", "MTurkRequesterServiceV20170117")
const neptune = AWSCorePrototype.RestXMLService("rds", "2014-10-31")
const networkmanager = AWSCorePrototype.RestXMLService("networkmanager", "2019-07-05")
const opsworks = AWSCorePrototype.JSONService("opsworks", "2013-02-18", "1.1", "OpsWorks_20130218")
const opsworkscm = AWSCorePrototype.JSONService("opsworks-cm", "2016-11-01", "1.1", "OpsWorksCM_V2016_11_01")
const organizations = AWSCorePrototype.JSONService("organizations", "2016-11-28", "1.1", "AWSOrganizationsV20161128")
const outposts = AWSCorePrototype.RestXMLService("outposts", "2019-12-03")
const personalize = AWSCorePrototype.JSONService("personalize", "2018-05-22", "1.1", "AmazonPersonalize")
const personalize_events = AWSCorePrototype.RestXMLService("personalize-events", "2018-03-22")
const personalize_runtime = AWSCorePrototype.RestXMLService("personalize-runtime", "2018-05-22")
const pi = AWSCorePrototype.JSONService("pi", "2018-02-27", "1.1", "PerformanceInsightsv20180227")
const pinpoint = AWSCorePrototype.RestXMLService("pinpoint", "2016-12-01")
const pinpoint_email = AWSCorePrototype.RestXMLService("email", "2018-07-26")
const polly = AWSCorePrototype.RestXMLService("polly", "2016-06-10")
const pricing = AWSCorePrototype.JSONService("api.pricing", "2017-10-15", "1.1", "AWSPriceListService")
const qldb = AWSCorePrototype.RestXMLService("qldb", "2019-01-02")
const qldb_session = AWSCorePrototype.JSONService("session.qldb", "2019-07-11", "1.0", "QLDBSession")
const quicksight = AWSCorePrototype.RestXMLService("quicksight", "2018-04-01")
const ram = AWSCorePrototype.RestXMLService("ram", "2018-01-04")
const rds = AWSCorePrototype.RestXMLService("rds", "2013-01-10")
const rds = AWSCorePrototype.RestXMLService("rds", "2013-02-12")
const rds = AWSCorePrototype.RestXMLService("rds", "2013-09-09")
const rds = AWSCorePrototype.RestXMLService("rds", "2014-09-01")
const rds = AWSCorePrototype.RestXMLService("rds", "2014-10-31")
const rds_data = AWSCorePrototype.RestXMLService("rds-data", "2018-08-01")
const redshift = AWSCorePrototype.RestXMLService("redshift", "2012-12-01")
const rekognition = AWSCorePrototype.JSONService("rekognition", "2016-06-27", "1.1", "RekognitionService")
const resource_groups = AWSCorePrototype.RestXMLService("resource-groups", "2017-11-27")
const resource_groups_tagging_api = AWSCorePrototype.JSONService("tagging", "2017-01-26", "1.1", "ResourceGroupsTaggingAPI_20170126")
const robomaker = AWSCorePrototype.RestXMLService("robomaker", "2018-06-29")
const route_53 = AWSCorePrototype.RestXMLService("route53", "2013-04-01")
const route_53_domains = AWSCorePrototype.JSONService("route53domains", "2014-05-15", "1.1", "Route53Domains_v20140515")
const route53resolver = AWSCorePrototype.JSONService("route53resolver", "2018-04-01", "1.1", "Route53Resolver")
const lex_runtime_service = AWSCorePrototype.RestXMLService("runtime.lex", "2016-11-28")
const sagemaker_runtime = AWSCorePrototype.RestXMLService("runtime.sagemaker", "2017-05-13")
const s3 = AWSCorePrototype.RestXMLService("s3", "2006-03-01")
const s3_control = AWSCorePrototype.RestXMLService("s3-control", "2018-08-20")
const sagemaker = AWSCorePrototype.JSONService("api.sagemaker", "2017-07-24", "1.1", "SageMaker")
const sagemaker_a2i_runtime = AWSCorePrototype.RestXMLService("a2i-runtime.sagemaker", "2019-11-07")
const savingsplans = AWSCorePrototype.RestXMLService("savingsplans", "2019-06-28")
const schemas = AWSCorePrototype.RestXMLService("schemas", "2019-12-02")
const simpledb = AWSCorePrototype.RestXMLService("sdb", "2009-04-15")
const secrets_manager = AWSCorePrototype.JSONService("secretsmanager", "2017-10-17", "1.1", "secretsmanager")
const securityhub = AWSCorePrototype.RestXMLService("securityhub", "2018-10-26")
const serverlessapplicationrepository = AWSCorePrototype.RestXMLService("serverlessrepo", "2017-09-08")
const service_quotas = AWSCorePrototype.JSONService("servicequotas", "2019-06-24", "1.1", "ServiceQuotasV20190624")
const service_catalog = AWSCorePrototype.JSONService("servicecatalog", "2015-12-10", "1.1", "AWS242ServiceCatalogService")
const servicediscovery = AWSCorePrototype.JSONService("servicediscovery", "2017-03-14", "1.1", "Route53AutoNaming_v20170314")
const sesv2 = AWSCorePrototype.RestXMLService("email", "2019-09-27")
const shield = AWSCorePrototype.JSONService("shield", "2016-06-02", "1.1", "AWSShield_20160616")
const signer = AWSCorePrototype.RestXMLService("signer", "2017-08-25")
const sms = AWSCorePrototype.JSONService("sms", "2016-10-24", "1.1", "AWSServerMigrationService_V2016_10_24")
const pinpoint_sms_voice = AWSCorePrototype.RestXMLService("sms-voice.pinpoint", "2018-09-05")
const snowball = AWSCorePrototype.JSONService("snowball", "2016-06-30", "1.1", "AWSIESnowballJobManagementService")
const sns = AWSCorePrototype.RestXMLService("sns", "2010-03-31")
const sqs = AWSCorePrototype.RestXMLService("sqs", "2012-11-05")
const ssm = AWSCorePrototype.JSONService("ssm", "2014-11-06", "1.1", "AmazonSSM")
const sso = AWSCorePrototype.RestXMLService("portal.sso", "2019-06-10")
const sso_oidc = AWSCorePrototype.RestXMLService("oidc", "2019-06-10")
const sfn = AWSCorePrototype.JSONService("states", "2016-11-23", "1.0", "AWSStepFunctions")
const storage_gateway = AWSCorePrototype.JSONService("storagegateway", "2013-06-30", "1.1", "StorageGateway_20130630")
const dynamodb_streams = AWSCorePrototype.JSONService("streams.dynamodb", "2012-08-10", "1.0", "DynamoDBStreams_20120810")
const sts = AWSCorePrototype.RestXMLService("sts", "2011-06-15")
const support = AWSCorePrototype.JSONService("support", "2013-04-15", "1.1", "AWSSupport_20130415")
const swf = AWSCorePrototype.JSONService("swf", "2012-01-25", "1.0", "SimpleWorkflowService")
const textract = AWSCorePrototype.JSONService("textract", "2018-06-27", "1.1", "Textract")
const transcribe = AWSCorePrototype.JSONService("transcribe", "2017-10-26", "1.1", "Transcribe")
const transfer = AWSCorePrototype.JSONService("transfer", "2018-11-05", "1.1", "TransferService")
const translate = AWSCorePrototype.JSONService("translate", "2017-07-01", "1.1", "AWSShineFrontendService_20170701")
const waf = AWSCorePrototype.JSONService("waf", "2015-08-24", "1.1", "AWSWAF_20150824")
const waf_regional = AWSCorePrototype.JSONService("waf-regional", "2016-11-28", "1.1", "AWSWAF_Regional_20161128")
const wafv2 = AWSCorePrototype.JSONService("wafv2", "2019-07-29", "1.1", "AWSWAF_20190729")
const workdocs = AWSCorePrototype.RestXMLService("workdocs", "2016-05-01")
const worklink = AWSCorePrototype.RestXMLService("worklink", "2018-09-25")
const workmail = AWSCorePrototype.JSONService("workmail", "2017-10-01", "1.1", "WorkMailService")
const workmailmessageflow = AWSCorePrototype.RestXMLService("workmailmessageflow", "2019-05-01")
const workspaces = AWSCorePrototype.JSONService("workspaces", "2015-04-08", "1.1", "WorkspacesService")
const xray = AWSCorePrototype.RestXMLService("xray", "2016-04-12")

end

max_stars_repo_name: nicoleepp/JuliaCloud-AWS-State

max_issues_repo_name: nicoleepp/JuliaCloud-AWS-State

max_forks_repo_name: nicoleepp/JuliaCloud-AWS-State


text:
# changing color saturation by channel
using Images, ImageView
img = load("sample-images/cats-3061372_640.jpg");
img_ch_view = channelview(img);
img_ch_view = permuteddimsview(img_ch_view, (2, 3, 1));

x_coords = 320:640

img_ch_view[:, x_coords, 1] = min.(img_ch_view[:, x_coords, 1] .* 1.1, 1);
img_ch_view[:, x_coords, 2] = min.(img_ch_view[:, x_coords, 2] .* 1.2, 1);
img_ch_view[:, x_coords, 3] = min.(img_ch_view[:, x_coords, 3] .* 1.4, 1);
imshow(img)


text:
export vertices, num_vertices, vertex_dims

"""
    vertices( scenario :: BlackBox;
        rep = "normalized" :: String
    ) :: Vector{Vector{Int64}}

Generates the Local Polytope vertices for a BlackBox scenario. Valid represenations
are:

*`rep == "normalized"` or `rep == "generalized"`.
"""
function vertices(
    scenario :: BlackBox;
    rep = "normalized" :: String
) :: Vector{Vector{Int64}}
    if !(rep in ["normalized","generalized"])
        throw(DomainError(rep, "Argument `rep` must be either 'normalized' or 'generalized'."))
    end

    strategies = deterministic_strategies(scenario)

    vertices = (rep == "normalized") ? map(
        s -> s[1:(end-1),:][:], strategies
    ) : map(
        s -> s[:], strategies
    )

    vertices
end

"""
    vertices( scenario :: LocalSignaling;
        rep = "normalized" :: String
        rank_d_only = false :: Bool
    ) ::  Vector{Vector{Int64}}

Generates the deterministic strategies for the local polytope of `LocalSignaling`
scenarios.  The `rank_d_only` keyword arg specifies  whether to  exclude vertices
which use fewer dits  of communication and thus have a  matrix rank less than d.

!!! warning
    The vertices computed in this method are vectorized directly from a strategy
    matrix by column-majorization. These vertices are distinct from those produced
    by older `LocalPolytope.vertices()` methods which are row-majorized.
"""
function vertices(scenario :: LocalSignaling;
    rep = "normalized" :: String,
    rank_d_only = false :: Bool
) ::  Vector{Vector{Int64}}

    if !(rep in ("normalized", "generalized"))
        throw(DomainError(rep, "Argument `rep` must be either 'normalized' or 'generalized'."))
    end

    num_rows = (rep == "normalized") ? scenario.Y - 1 : scenario.Y
    lower_dits_bound = rank_d_only ? scenario.d : 1

    vertices = Vector{Vector{Int64}}(undef, num_vertices(scenario, rank_d_only = rank_d_only))
    vertex_id = 1
    for dits in lower_dits_bound:scenario.d
        d_perms = permutations(1:dits)
        X_partitions = stirling2_partitions(scenario.X, dits)
        Y_combinations = combinations(1:scenario.Y, dits)

        for Y in Y_combinations
            for d in d_perms
                for X in X_partitions
                    # construct matrix m from permutation ids. This is faster
                    # than performing matrix multiplication
                    m = zeros(Int64, num_rows, scenario.X)
                    for i in 1:dits
                        if Y[i] > num_rows
                            continue
                        end

                        row_id = Y[i]
                        partition_ids = X[d[i]]

                        m[row_id, partition_ids] .= 1
                    end
                    vertices[vertex_id] = m[:]

                    vertex_id +=  1
                end
            end
        end
    end

    vertices
end

"""
    vertices( scenario :: BipartiteNonSignaling,
        rep="non-signaling" :: String
    ) :: Vector{Vector{Int64}}

Enumerates the LocalPolytope vertices for the [`BipartiteNonSignaling`](@ref) scenario.
Valid representations for the vertices include:
* `"non-signaling"`, `"normalized"`, `"generalized"`

A `DomainError` is thrown if a valid representation is not specified.
"""
function vertices(scenario :: BipartiteNonSignaling, rep="non-signaling" :: String) :: Vector{Vector{Int64}}
    alice_strategies = deterministic_strategies(scenario.A, scenario.X)
    bob_strategies = deterministic_strategies(scenario.B, scenario.Y)

    if rep == "non-signaling"
        α_strategies = map(s -> s[1:end-1,:], alice_strategies)
        β_strategies = map(s -> s[1:end-1,:], bob_strategies)

        dim_α = scenario.X*(scenario.A - 1)
        dim_β = scenario.Y*(scenario.B - 1)
        dim_v = dim_α + dim_β + dim_α*dim_β

        strategies = map( (α, β) for α in α_strategies for β in β_strategies) do (α,β)
            v = zeros(Int64, dim_v)

            v[1:dim_α+dim_β] = vcat(α[:],β[:])
            v[dim_α+dim_β+1:end] = kron(α,β)[:]

            v
        end

        return strategies

    elseif rep in ("normalized", "generalized")

        strategies = (rep == "normalized") ? map(
            (α, β) for α in alice_strategies for β in bob_strategies) do (α,β)
                kron(α,β)[1:end-1,:][:]
            end : map(
            (α, β) for α in alice_strategies for β in bob_strategies) do (α,β)
                kron(α,β)[:]
            end

        return strategies
    else
        throw(DomainError(rep, "`rep in (\"non-signaling\",\"normalized\",\"generalized\")` must be satisfied" ))
    end
end

"""
    num_vertices( scenario :: BlackBox ) :: Int64

For ``n`` outputs and ``m`` inputs the number of vertices ``|\\mathcal{V}|`` are counted:

```math
|\\mathbf{V}| = n^m
```
"""
function num_vertices(scenario :: BlackBox) :: Int64
    scenario.num_out^scenario.num_in
end

"""
    num_vertices( scenario :: LocalSignaling;
        rank_d_only = false :: Bool
    ) :: Int64

If `rank_d_only = true`, then only strategies using  `d`-dits are counted. For
``X`` inputs and ``Y`` outputs the number of vertices ``|\\mathcal{V}|`` are counted:

```math
|\\mathbf{V}| = \\sum_{c=1}^d \\left\\{X \\atop c \\right\\}\\binom{Y}{c}c!
```
"""
function num_vertices(scenario :: LocalSignaling; rank_d_only = false :: Bool) :: Int64
    lower_dits_bound = rank_d_only ? scenario.d : 1
    sum(map(i -> stirling2(scenario.X, i)*binomial(scenario.Y, i)*factorial(i), lower_dits_bound:scenario.d))
end

"""
    num_vertices( scenario :: BipartiteNonSignaling ) :: Int64

For two non-signaling black-boxes with ``X`` and ``Y`` inputs and ``A`` and ``B``
outputs respectively, the number of vertices ``|\\mathcal{V}|`` are counted:

```math
|\\mathbf{V}| = A^X B^Y
```
"""
function num_vertices(scenario :: BipartiteNonSignaling) :: Int64
    black_box_A = BlackBox(scenario.A, scenario.X)
    black_box_B = BlackBox(scenario.B, scenario.Y)

    num_vertices(black_box_A)*num_vertices(black_box_B)
end

"""
For the given [`Scenario`](@ref), returns the length of the vertex in the representation
specified by `rep`. A `DomainError` is thrown if the `rep` is invalid.

    vertex_dims(scenario :: Union{BlackBox,LocalSignaling}, rep :: String) :: Int64

Valid values of `rep` are `"normalized"` and `"generalized"`.
"""
function vertex_dims(scenario :: Union{BlackBox,LocalSignaling}, rep :: String) :: Int64
    if !(rep in ["normalized", "generalized"])
        throw(DomainError(rep, "Argument `rep` must be either \"normalized\" or \"generalized\"."))
    end

    strat_dims = strategy_dims(scenario)

    (rep == "normalized") ? strat_dims[2] * (strat_dims[1] - 1) : strat_dims[2] * strat_dims[1]
end

"""
    vertex_dims( scenario:: BipartiteNonSignaling, rep :: String ) :: Int64

Valid values for `rep` include:
* "non-signaling"
* "normalized"
* "generalized"
"""
function vertex_dims(scenario :: BipartiteNonSignaling, rep :: String) :: Int64
    if !(rep in ["non-signaling", "normalized", "generalized"])
        throw(DomainError(rep, "Argument `rep` must be either \"non-signaling\", \"normalized\", or \"generalized\"."))
    end

    dim_α = scenario.X*(scenario.A - 1)
    dim_β = scenario.Y*(scenario.B - 1)

    (rep == "non-signaling") ? dim_α + dim_β + dim_α*dim_β : vertex_dims(
        BlackBox(strategy_dims(scenario)...), rep)
end


text:
#'
#' # Testing exercise generation
#' #'
#' This would be the text for the exercise... The main document that we maintain is the complete notebook and in each code cell, the parts that should be removed in the skeleton are separated with lines containing only #%
#' $Ax = b$
#'
#+

println("Hello!")

#TODO
#% # Solution
#% using JuMP
#% m = Model()

#+

using Cbc
function solve_model(m)
    #TODO
#%     optimize!(m)
end


text:
"""
    struct Trajectory

Representation of a single crane trajectory as part of as solution to the GCSPET. The crane movements are modeled as a sequence of crane-positions through time.
Positions relevant to the execution of a job are marked with the id of that particular job. Each job-id thus occurs twice: once for the starting point of the job and once more for finishing.
"""
struct Trajectory
    id::Int
    X::Vector{Int}
    T::Vector{Int}
    J::Vector{Int}
end

id(t::Trajectory) = t.id

jobset(t::Trajectory) = unique(filter( x -> x > 0, t.J))

function starttimes(t::Trajectory)
    d = Dict{Int,Int}()
    for j in jobset(t)
        idx = findfirst(j, t.J)
        d[j] = t.T[idx]
    end
    return d
end

function endtimes(t::Trajectory)
    for j in jobset(t)
        idx = findlast(j, t.J)
        d[j] = t.T[idx]
    end
    return d
end

function durations(t::Trajectory)
    d = Dict{Int,Int}
    S = starttimes(t)
    E = endtimes(t)
    for j in jobset(t)
        d[j] = E[j] - S[j]
    end
    return d
end

function Base.show(io::IO, ::MIME"text/plain", t::Trajectory)
    println(io, "Crane $(t.id)")
    println(io, "t\tx\tjob")
    for (x,y,z) in Iterators.zip(t.T, t.X, t.J)
        if z > 0
            println(io, "$x\t$y\t$z")
        else
            println(io, "$x\t$y\t")
        end
    end
end

##

"""
    struct Solution

Representation of a GCSPET-solution, consisting of all the crane trajectories and job executions.
"""
struct Solution
    name::String
    jobs::Vector{Job}
    cranes::Vector{Crane}
    trajectories::Vector{Trajectory}
end

name(s::Solution) = s.name
jobs(s::Solution) = s.jobs
cranes(s::Solution) = s.cranes
trajectories(s::Solution) = s.trajectories

njobs(s::Solution) = length(jobs(s))
ncranes(s::Solution) = length(cranes(s))
load(s::Solution) = mapreduce(istruck, +, jobs(s))/njobs(s)

crane_starting_pos(s::Solution) = map(starting_position, cranes(s))

objval(s::Solution) = makespan(s) + truckwaitingtime(s)

makespan(s::Solution) = maximum(t -> last(t.T) , trajectories(s))

function truckwaitingtime(s::Solution)
    twt = 0
    for j in Iterators.filter(istruck, jobs(s))
        tc = t_completion(s, id(j))
        twt += tc - t_arrival(j)
    end
    return twt
end

function t_completion(s::Solution, id)
    for t in trajectories(s)
        idx = findlast(x -> x == id ,t.J)
        !isnothing(idx) && return t.T[idx]
    end
    return error("The provided job-id is not in the schedule!")
end

 function Base.show(io::IO, ::MIME"text/plain", s::Solution)
    njobs = length(jobs(s))
    ncranes = length(cranes(s))
    cmax = makespan(s)
    twt = truckwaitingtime(s)
    println(io, "GCSPET-solution: $(name(s))")
    println(io, "njobs: $njobs, ncranes: $ncranes")
    println(io, "CMAX: $cmax, TWT: $twt")
    println(io, "--------------")
    for t in trajectories(s)
        println(io, t)
    end
 end