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

@doc raw"""DaemonSetStatus represents the current status of a daemon set.

    IoK8sApiAppsV1beta2DaemonSetStatus(;
        collisionCount=nothing,
        conditions=nothing,
        currentNumberScheduled=nothing,
        desiredNumberScheduled=nothing,
        numberAvailable=nothing,
        numberMisscheduled=nothing,
        numberReady=nothing,
        numberUnavailable=nothing,
        observedGeneration=nothing,
        updatedNumberScheduled=nothing,
    )

    - collisionCount::Int32 : Count of hash collisions for the DaemonSet. The DaemonSet controller uses this field as a collision avoidance mechanism when it needs to create the name for the newest ControllerRevision.
    - conditions::Vector{IoK8sApiAppsV1beta2DaemonSetCondition} : Represents the latest available observations of a DaemonSet's current state.
    - currentNumberScheduled::Int32 : The number of nodes that are running at least 1 daemon pod and are supposed to run the daemon pod. More info: https://kubernetes.io/docs/concepts/workloads/controllers/daemonset/
    - desiredNumberScheduled::Int32 : The total number of nodes that should be running the daemon pod (including nodes correctly running the daemon pod). More info: https://kubernetes.io/docs/concepts/workloads/controllers/daemonset/
    - numberAvailable::Int32 : The number of nodes that should be running the daemon pod and have one or more of the daemon pod running and available (ready for at least spec.minReadySeconds)
    - numberMisscheduled::Int32 : The number of nodes that are running the daemon pod, but are not supposed to run the daemon pod. More info: https://kubernetes.io/docs/concepts/workloads/controllers/daemonset/
    - numberReady::Int32 : The number of nodes that should be running the daemon pod and have one or more of the daemon pod running and ready.
    - numberUnavailable::Int32 : The number of nodes that should be running the daemon pod and have none of the daemon pod running and available (ready for at least spec.minReadySeconds)
    - observedGeneration::Int64 : The most recent generation observed by the daemon set controller.
    - updatedNumberScheduled::Int32 : The total number of nodes that are running updated daemon pod
"""
mutable struct IoK8sApiAppsV1beta2DaemonSetStatus <: SwaggerModel
    collisionCount::Any # spec type: Union{ Nothing, Int32 } # spec name: collisionCount
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiAppsV1beta2DaemonSetCondition} } # spec name: conditions
    currentNumberScheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: currentNumberScheduled
    desiredNumberScheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: desiredNumberScheduled
    numberAvailable::Any # spec type: Union{ Nothing, Int32 } # spec name: numberAvailable
    numberMisscheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: numberMisscheduled
    numberReady::Any # spec type: Union{ Nothing, Int32 } # spec name: numberReady
    numberUnavailable::Any # spec type: Union{ Nothing, Int32 } # spec name: numberUnavailable
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration
    updatedNumberScheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: updatedNumberScheduled

    function IoK8sApiAppsV1beta2DaemonSetStatus(;collisionCount=nothing, conditions=nothing, currentNumberScheduled=nothing, desiredNumberScheduled=nothing, numberAvailable=nothing, numberMisscheduled=nothing, numberReady=nothing, numberUnavailable=nothing, observedGeneration=nothing, updatedNumberScheduled=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("collisionCount"), collisionCount)
        setfield!(o, Symbol("collisionCount"), collisionCount)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("currentNumberScheduled"), currentNumberScheduled)
        setfield!(o, Symbol("currentNumberScheduled"), currentNumberScheduled)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("desiredNumberScheduled"), desiredNumberScheduled)
        setfield!(o, Symbol("desiredNumberScheduled"), desiredNumberScheduled)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("numberAvailable"), numberAvailable)
        setfield!(o, Symbol("numberAvailable"), numberAvailable)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("numberMisscheduled"), numberMisscheduled)
        setfield!(o, Symbol("numberMisscheduled"), numberMisscheduled)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("numberReady"), numberReady)
        setfield!(o, Symbol("numberReady"), numberReady)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("numberUnavailable"), numberUnavailable)
        setfield!(o, Symbol("numberUnavailable"), numberUnavailable)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        validate_property(IoK8sApiAppsV1beta2DaemonSetStatus, Symbol("updatedNumberScheduled"), updatedNumberScheduled)
        setfield!(o, Symbol("updatedNumberScheduled"), updatedNumberScheduled)
        o
    end
end # type IoK8sApiAppsV1beta2DaemonSetStatus

const _property_map_IoK8sApiAppsV1beta2DaemonSetStatus = Dict{Symbol,Symbol}(Symbol("collisionCount")=>Symbol("collisionCount"), Symbol("conditions")=>Symbol("conditions"), Symbol("currentNumberScheduled")=>Symbol("currentNumberScheduled"), Symbol("desiredNumberScheduled")=>Symbol("desiredNumberScheduled"), Symbol("numberAvailable")=>Symbol("numberAvailable"), Symbol("numberMisscheduled")=>Symbol("numberMisscheduled"), Symbol("numberReady")=>Symbol("numberReady"), Symbol("numberUnavailable")=>Symbol("numberUnavailable"), Symbol("observedGeneration")=>Symbol("observedGeneration"), Symbol("updatedNumberScheduled")=>Symbol("updatedNumberScheduled"))
const _property_types_IoK8sApiAppsV1beta2DaemonSetStatus = Dict{Symbol,String}(Symbol("collisionCount")=>"Int32", Symbol("conditions")=>"Vector{IoK8sApiAppsV1beta2DaemonSetCondition}", Symbol("currentNumberScheduled")=>"Int32", Symbol("desiredNumberScheduled")=>"Int32", Symbol("numberAvailable")=>"Int32", Symbol("numberMisscheduled")=>"Int32", Symbol("numberReady")=>"Int32", Symbol("numberUnavailable")=>"Int32", Symbol("observedGeneration")=>"Int64", Symbol("updatedNumberScheduled")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2DaemonSetStatus }) = collect(keys(_property_map_IoK8sApiAppsV1beta2DaemonSetStatus))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2DaemonSetStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta2DaemonSetStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2DaemonSetStatus }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2DaemonSetStatus[property_name]

function check_required(o::IoK8sApiAppsV1beta2DaemonSetStatus)
    (getproperty(o, Symbol("currentNumberScheduled")) === nothing) && (return false)
    (getproperty(o, Symbol("desiredNumberScheduled")) === nothing) && (return false)
    (getproperty(o, Symbol("numberMisscheduled")) === nothing) && (return false)
    (getproperty(o, Symbol("numberReady")) === nothing) && (return false)
    true
end

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


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Broadcast styles
import Base.Broadcast
using Base.Broadcast: DefaultArrayStyle, Broadcasted, tail

struct StructuredMatrixStyle{T} <: Broadcast.AbstractArrayStyle{2} end
StructuredMatrixStyle{T}(::Val{2}) where {T} = StructuredMatrixStyle{T}()
StructuredMatrixStyle{T}(::Val{N}) where {T,N} = Broadcast.DefaultArrayStyle{N}()

const StructuredMatrix = Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UnitLowerTriangular,UpperTriangular,UnitUpperTriangular}
Broadcast.BroadcastStyle(::Type{T}) where {T<:StructuredMatrix} = StructuredMatrixStyle{T}()

# Promotion of broadcasts between structured matrices. This is slightly unusual
# as we define them symmetrically. This allows us to have a fallback to DefaultArrayStyle{2}().
# Diagonal can cavort with all the other structured matrix types.
# Bidiagonal doesn't know if it's upper or lower, so it becomes Tridiagonal
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Diagonal}, ::StructuredMatrixStyle{<:Diagonal}) =
    StructuredMatrixStyle{Diagonal}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Diagonal}, ::StructuredMatrixStyle{<:Union{Bidiagonal,SymTridiagonal,Tridiagonal}}) =
    StructuredMatrixStyle{Tridiagonal}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Diagonal}, ::StructuredMatrixStyle{<:Union{LowerTriangular,UnitLowerTriangular}}) =
    StructuredMatrixStyle{LowerTriangular}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Diagonal}, ::StructuredMatrixStyle{<:Union{UpperTriangular,UnitUpperTriangular}}) =
    StructuredMatrixStyle{UpperTriangular}()

Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Bidiagonal}, ::StructuredMatrixStyle{<:Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal}}) =
    StructuredMatrixStyle{Tridiagonal}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:SymTridiagonal}, ::StructuredMatrixStyle{<:Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal}}) =
    StructuredMatrixStyle{Tridiagonal}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Tridiagonal}, ::StructuredMatrixStyle{<:Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal}}) =
    StructuredMatrixStyle{Tridiagonal}()

Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:LowerTriangular}, ::StructuredMatrixStyle{<:Union{Diagonal,LowerTriangular,UnitLowerTriangular}}) =
    StructuredMatrixStyle{LowerTriangular}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:UpperTriangular}, ::StructuredMatrixStyle{<:Union{Diagonal,UpperTriangular,UnitUpperTriangular}}) =
    StructuredMatrixStyle{UpperTriangular}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:UnitLowerTriangular}, ::StructuredMatrixStyle{<:Union{Diagonal,LowerTriangular,UnitLowerTriangular}}) =
    StructuredMatrixStyle{LowerTriangular}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:UnitUpperTriangular}, ::StructuredMatrixStyle{<:Union{Diagonal,UpperTriangular,UnitUpperTriangular}}) =
    StructuredMatrixStyle{UpperTriangular}()

Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Union{LowerTriangular,UnitLowerTriangular}}, ::StructuredMatrixStyle{<:Union{UpperTriangular,UnitUpperTriangular}}) =
    StructuredMatrixStyle{Matrix}()
Broadcast.BroadcastStyle(::StructuredMatrixStyle{<:Union{UpperTriangular,UnitUpperTriangular}}, ::StructuredMatrixStyle{<:Union{LowerTriangular,UnitLowerTriangular}}) =
    StructuredMatrixStyle{Matrix}()

# Make sure that `StructuredMatrixStyle{<:Matrix}` doesn't ever end up falling
# through and give back `DefaultArrayStyle{2}`
Broadcast.BroadcastStyle(T::StructuredMatrixStyle{<:Matrix}, ::StructuredMatrixStyle) = T
Broadcast.BroadcastStyle(::StructuredMatrixStyle, T::StructuredMatrixStyle{<:Matrix}) = T
Broadcast.BroadcastStyle(T::StructuredMatrixStyle{<:Matrix}, ::StructuredMatrixStyle{<:Matrix}) = T

# All other combinations fall back to the default style
Broadcast.BroadcastStyle(::StructuredMatrixStyle, ::StructuredMatrixStyle) = DefaultArrayStyle{2}()

# And a definition akin to similar using the structured type:
structured_broadcast_alloc(bc, ::Type{<:Diagonal}, ::Type{ElType}, n) where {ElType} =
    Diagonal(Array{ElType}(undef, n))
# Bidiagonal is tricky as we need to know if it's upper or lower. The promotion
# system will return Tridiagonal when there's more than one Bidiagonal, but when
# there's only one, we need to make figure out upper or lower
merge_uplos(::Nothing, ::Nothing) = nothing
merge_uplos(a, ::Nothing) = a
merge_uplos(::Nothing, b) = b
merge_uplos(a, b) = a == b ? a : 'T'

find_uplo(a::Bidiagonal) = a.uplo
find_uplo(a) = nothing
find_uplo(bc::Broadcasted) = mapreduce(find_uplo, merge_uplos, bc.args, init=nothing)

function structured_broadcast_alloc(bc, ::Type{<:Bidiagonal}, ::Type{ElType}, n) where {ElType}
    uplo = find_uplo(bc)
    if uplo == 'T'
        return Tridiagonal(Array{ElType}(undef, n-1), Array{ElType}(undef, n), Array{ElType}(undef, n-1))
    end
    return Bidiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n-1), uplo)
end
structured_broadcast_alloc(bc, ::Type{<:SymTridiagonal}, ::Type{ElType}, n) where {ElType} =
    SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n-1))
structured_broadcast_alloc(bc, ::Type{<:Tridiagonal}, ::Type{ElType}, n) where {ElType} =
    Tridiagonal(Array{ElType}(undef, n-1),Array{ElType}(undef, n),Array{ElType}(undef, n-1))
structured_broadcast_alloc(bc, ::Type{<:LowerTriangular}, ::Type{ElType}, n) where {ElType} =
    LowerTriangular(Array{ElType}(undef, n, n))
structured_broadcast_alloc(bc, ::Type{<:UpperTriangular}, ::Type{ElType}, n) where {ElType} =
    UpperTriangular(Array{ElType}(undef, n, n))
structured_broadcast_alloc(bc, ::Type{<:UnitLowerTriangular}, ::Type{ElType}, n) where {ElType} =
    UnitLowerTriangular(Array{ElType}(undef, n, n))
structured_broadcast_alloc(bc, ::Type{<:UnitUpperTriangular}, ::Type{ElType}, n) where {ElType} =
    UnitUpperTriangular(Array{ElType}(undef, n, n))
structured_broadcast_alloc(bc, ::Type{<:Matrix}, ::Type{ElType}, n) where {ElType} =
    Matrix(Array{ElType}(undef, n, n))

# A _very_ limited list of structure-preserving functions known at compile-time. This list is
# derived from the formerly-implemented `broadcast` methods in 0.6. Note that this must
# preserve both zeros and ones (for Unit***erTriangular) and symmetry (for SymTridiagonal)
const TypeFuncs = Union{typeof(round),typeof(trunc),typeof(floor),typeof(ceil)}
isstructurepreserving(bc::Broadcasted) = isstructurepreserving(bc.f, bc.args...)
isstructurepreserving(::Union{typeof(abs),typeof(big)}, ::StructuredMatrix) = true
isstructurepreserving(::TypeFuncs, ::StructuredMatrix) = true
isstructurepreserving(::TypeFuncs, ::Ref{<:Type}, ::StructuredMatrix) = true
function isstructurepreserving(::typeof(Base.literal_pow), ::Ref{typeof(^)}, ::StructuredMatrix, ::Ref{Val{N}}) where N
    return N isa Integer && N > 0
end
isstructurepreserving(f, args...) = false

_iszero(n::Number) = iszero(n)
_iszero(x) = x == 0
fzeropreserving(bc) = (v = fzero(bc); !ismissing(v) && _iszero(v))
# Like sparse matrices, we assume that the zero-preservation property of a broadcasted
# expression is stable.  We can test the zero-preservability by applying the function
# in cases where all other arguments are known scalars against a zero from the structured
# matrix. If any non-structured matrix argument is not a known scalar, we give up.
fzero(x::Number) = x
fzero(::Type{T}) where T = T
fzero(r::Ref) = r[]
fzero(t::Tuple{Any}) = t[1]
fzero(S::StructuredMatrix) = zero(eltype(S))
fzero(x) = missing
function fzero(bc::Broadcast.Broadcasted)
    args = map(fzero, bc.args)
    return any(ismissing, args) ? missing : bc.f(args...)
end

function Base.similar(bc::Broadcasted{StructuredMatrixStyle{T}}, ::Type{ElType}) where {T,ElType}
    inds = axes(bc)
    if isstructurepreserving(bc) || (fzeropreserving(bc) && !(T <: Union{SymTridiagonal,UnitLowerTriangular,UnitUpperTriangular}))
        return structured_broadcast_alloc(bc, T, ElType, length(inds[1]))
    end
    return similar(convert(Broadcasted{DefaultArrayStyle{ndims(bc)}}, bc), ElType)
end

function copyto!(dest::Diagonal, bc::Broadcasted{<:StructuredMatrixStyle})
    !isstructurepreserving(bc) && !fzeropreserving(bc) && return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    axs = axes(dest)
    axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
    for i in axs[1]
        dest.diag[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i))
    end
    return dest
end

function copyto!(dest::Bidiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
    !isstructurepreserving(bc) && !fzeropreserving(bc) && return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    axs = axes(dest)
    axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
    for i in axs[1]
        dest.dv[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i))
    end
    if dest.uplo == 'U'
        for i = 1:size(dest, 1)-1
            dest.ev[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i+1))
        end
    else
        for i = 1:size(dest, 1)-1
            dest.ev[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i+1, i))
        end
    end
    return dest
end

function copyto!(dest::SymTridiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
    !isstructurepreserving(bc) && !fzeropreserving(bc) && return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    axs = axes(dest)
    axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
    for i in axs[1]
        dest.dv[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i))
    end
    for i = 1:size(dest, 1)-1
        v = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i+1))
        v == Broadcast._broadcast_getindex(bc, CartesianIndex(i+1, i)) || throw(ArgumentError("broadcasted assignment breaks symmetry between locations ($i, $(i+1)) and ($(i+1), $i)"))
        dest.ev[i] = v
    end
    return dest
end

function copyto!(dest::Tridiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
    !isstructurepreserving(bc) && !fzeropreserving(bc) && return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    axs = axes(dest)
    axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
    for i in axs[1]
        dest.d[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i))
    end
    for i = 1:size(dest, 1)-1
        dest.du[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i+1))
        dest.dl[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i+1, i))
    end
    return dest
end

function copyto!(dest::LowerTriangular, bc::Broadcasted{<:StructuredMatrixStyle})
    !isstructurepreserving(bc) && !fzeropreserving(bc) && return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    axs = axes(dest)
    axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
    for j in axs[2]
        for i in j:axs[1][end]
            dest.data[i,j] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, j))
        end
    end
    return dest
end

function copyto!(dest::UpperTriangular, bc::Broadcasted{<:StructuredMatrixStyle})
    !isstructurepreserving(bc) && !fzeropreserving(bc) && return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    axs = axes(dest)
    axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
    for j in axs[2]
        for i in 1:j
            dest.data[i,j] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, j))
        end
    end
    return dest
end

# We can also implement `map` and its promotion in terms of broadcast with a stricter dimension check
function map(f, A::StructuredMatrix, Bs::StructuredMatrix...)
    sz = size(A)
    all(map(B->size(B)==sz, Bs)) || throw(DimensionMismatch("dimensions must match"))
    return f.(A, Bs...)
end

ext: jl

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
#= Top-llevel file to run nearly neutral simulation from the command line.
Examples: (run from directory evotech/nearly-neutral/src):
> julia run.jl examples/ia_example1
> julia run.jl examples/is_example1
=#
# The following are top-level commands run when the file is executed by julia.
using Random
using Distributed
include("dfe_functions.jl")   # These need to be defined to read in the infinite sites configuration files
#push!(LOAD_PATH,"../src/")
#include("InfAlleles.jl")
#using Main.InfAlleles

if length(ARGS) == 0
  simname = "../experiments/examples/ia_example1"
else
  simname = ARGS[1]
  if length(ARGS) >= 2   # second command-line argument is random number seed
    seed = parse(Int,ARGS[2])
    println("random number seed: ",seed)
    Random.seed!(seed)    # Replaces function srand()
  end
end
println("reading configuration file: ","$(simname).jl")
include("$(simname).jl")
println("opening output CSV file: ","$(simname).csv")
stream = open("$(simname).csv","w")

if !isdefined(Main,:mu_list_flag)
  mu_list_flag=false
end
if nn_simtype==1
  using Main.InfAlleles
  println("infinite alleles model")
  include("run_infalleles.jl")
  ia_run_trials( mu_list_flag )
elseif nn_simtype == 2
  println("infinite sites model")
  include("run_infsites.jl")
  is_run_trials( mu_list_flag )
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using .APIUtils

using CEnum: @cenum

using Printf

import Libdl

using LazyArtifacts

# low-level wrappers
const CUdeviceptr = CuPtr{Cvoid}
const CUarray = CuArrayPtr{Cvoid}
const GLuint = Cuint    # FIXME: get these from somewhere
const GLenum = Cuint
include("libcuda_common.jl")
include("error.jl")
include("libcuda.jl")
include("libcuda_deprecated.jl")

const _libcuda = Ref{String}()
function libcuda()
    if !isassigned(_libcuda)
        # find the local system library
        system_driver = if Sys.iswindows()
            Libdl.find_library("nvcuda")
        else
            Libdl.find_library(["libcuda.so.1", "libcuda.so"])
        end
        if name == ""
            if Sys.iswindows()
                error("""Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
                         If you're sure it's installed, look for `libcuda.so` in your system and make sure it's discoverable by the linker.
                         Typically, that involves an entry in '/etc/ld.so.conf', or setting LD_LIBRARY_PATH.""")
            else
                error("""Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
                         If you're sure it's installed, look for `nvcuda.dll` in your system and make sure it's discoverable by the linker.
                         Typically, that involves adding an entry to PATH.""")
            end
        end

        # NOTE: we can't re-use the function wrappers from libcuda.so,
        #       as those use libcuda() whose value is cached by ccall.

        # query the system driver version
        function get_version(driver)
            library_handle = Libdl.dlopen(driver)
            try
                function_handle = Libdl.dlsym(library_handle, "cuDriverGetVersion")
                version_ref = Ref{Cint}()
                @check ccall(function_handle, CUresult, (Ptr{Cint},), version_ref)
                major, ver = divrem(version_ref[], 1000)
                minor, patch = divrem(ver, 10)
                return VersionNumber(major, minor, patch)
            finally
                Libdl.dlclose(library_handle)
            end
        end
        _system_version[] = get_version(system_driver)
        _libcuda[] = system_driver
        # XXX: apparently cuDriverGetVersion can be used before cuInit,
        #      despite the docs stating "any function [...] will return
        #      CUDA_ERROR_NOT_INITIALIZED"; is this a recent change?

        # if we're using an older driver; consider using forward compatibility
        function do_init(driver)
            library_handle = Libdl.dlopen(driver)
            try
                function_handle = Libdl.dlsym(library_handle, "cuInit")
                @check ccall(function_handle, CUresult, (UInt32,), 0)
            finally
                Libdl.dlclose(library_handle)
            end
            return
        end
        if _system_version[] < v"11.4"
            platform = Base.BinaryPlatforms.HostPlatform()
            platform.tags["cuda"] = "11.4"

            artifact = try
                @artifact_str("CUDA_compat", platform)
            catch ex
                @debug "Could not download forward compatibility package" exception=(ex,catch_backtrace())
                nothing
            end

            if artifact !== nothing
                # TODO: do we need to dlopen the JIT compiler library for it to
                #       be discoverable? doesn't that clash with a system one
                #       if compat cuInit fails? or should we load it _after_
                #       the compat driver initialization succeeds?
                #compat_compiler = joinpath(artifact, "lib", "libnvidia-ptxjitcompiler.so")
                #Libdl.dlopen(compat_compiler)

                compat_driver = joinpath(artifact, "lib", "libcuda.so")
                try
                    # just calling cuDriverGetVersion doesn't trigger errors,
                    # so perform driver intialization here
                    do_init(compat_driver)
                    _libcuda[] = compat_driver
                catch ex
                    @debug "Could not use forward compatibility package" exception=(ex,catch_backtrace())
                end
            end
        end

        cuInit(0)
    end

    # memoized because otherwise each ccall would perform discovery again
    _libcuda[]
end

# high-level wrappers
include("types.jl")
include("version.jl")
include("devices.jl")
include("context.jl")
include("stream.jl")
include("pool.jl")
include("memory.jl")
include("module.jl")
include("events.jl")
include("execution.jl")
include("profile.jl")
include("occupancy.jl")
include("graph.jl")

# TODO: figure out if these wrappers may use the runtime-esque state (stream(), context()).
#       it's inconsitent now: @finalize_in_ctx doesn't, memory.jl does use stream(), etc.


text:
function KalmanFilter(T::Matrix, Z::Matrix, Q::Matrix, H::Matrix, n::Integer, α0, y)

    # Initialize state variables
    αhat = repeat([zeros(size(α0))], n);
    αhat[1] = T*α0;

    println("I'm here 1")

    yhat = repeat([zeros(size(y[1]))], n);
    v = copy(yhat);

    P0 = [2.0]; # Initial guess given in the problem

    P = repeat([zeros(size(Q))], n);
    P[1] = T * P0 * T' + Q;

    println("I'm here 2, P[1] = $(P[1])")

    F = repeat([zeros(size(H))], n)
    G = repeat([zeros(size(Z * Q))], n)
    K = repeat([zeros(size(T * G[1]' * ones(size(H))))], n)

    for i in 1:n-1
        yhat[i] = Z * αhat[i];
        v[i] = y[i] - yhat[i];

        F[i] = Z * P[i] * Z' + H;
        G[i] = Z * P[i];
        K[i] = T * G[i]' * inv(F[i]);

        αhat[i + 1] = T * αhat[i] + K[i] * v[i];
        P[i + 1] = T * (P[i] - G[i]' * inv(F[i]) * G[i]) * T' + Q;
    end

    yhat[n] = Z * αhat[n];
    v[n] = y[n] - yhat[n];

    F[n] = Z * P[n] * Z' + H;
    G[n] = Z * P[n];
    K[n] = T * G[n]' * inv(F[n]);

    return αhat, v, yhat, P, F, G, K
end


text:
############################################################################################
"""
$(TYPEDEF)
Contains maximum memory for latent and observed data. Relevant for number of
propagation steps from initial particle distribution and for first time particles
are weighted with data point.

# Fields
$(TYPEDFIELDS)
"""
struct ParticleFilterMemory
    "Latent variable memory. Number of times that initial particle distribution is used."
    latent::Int64
    "Observed data memory. First particle weighting taken at this point."
    data::Int64
    "Number of times when sampled from initial distribution. Per default equal to latent."
    initial::Int64
    function ParticleFilterMemory(latent::Int64, data::Int64, initial::Int64)
        ArgCheck.@argcheck latent >= 0 "Latent memory cannot be negative."
        ArgCheck.@argcheck data >= 0 "Data memory cannot be negative."
        ArgCheck.@argcheck initial >= 0 "initial number of sampling steps cannot be negative."
        return new(latent, data, initial)
    end
end

function max(memory::ParticleFilterMemory)
    return max(memory.latent, memory.data)
end

############################################################################################
"""
$(SIGNATURES)
Propagate forward a single trajectory given ParticleKernel.

# Examples
```julia
```

"""
function propagate(
    _rng, kernel::ParticleKernel, memory::ParticleFilterMemory, reference::AbstractArray{P}
) where {P}
    trajectory = Vector{P}(undef, size(reference, 1))
    Ninitial = memory.initial #max(1, memory.latent)
    ## Initialize particle
    for iter in 1:Ninitial
        trajectory[iter] = initial(_rng, kernel)
    end
    ## Propoagate particle forward
    for iter in (Ninitial + 1):size(trajectory, 1)
        trajectory[iter] = transition(_rng, kernel, trajectory, iter)
    end
    return trajectory
end

############################################################################################
"""
$(SIGNATURES)
Check if function is exectuable.

# Examples
```julia
```

"""
function _try(f::Function, input...)
    try
        f(input...)
    catch err
        return err
    end
end

"""
$(SIGNATURES)
Return first iteration at which function is executable.

# Examples
```julia
```

"""
function _checkmemory(f::Function; maxiter=1000)
    for iter in Base.OneTo(maxiter)
        val = _try(f, iter)
        if !isa(val, Exception)
            return iter
        end
    end
    return "Memory for Objective exceeds 1000 - check your dynamics(objective) implementation."
end

"""
$(SIGNATURES)
Guess memory of ParticleKernel.

# Examples
```julia
```

"""
function _guessmemory(
    _rng::Random.AbstractRNG, kernel::ParticleKernel, reference::AbstractArray{T}
) where {T}
    ## Initialize state and evidence transition distribution as a function of the current iteration
    guess_state(iter) = transition(_rng, kernel, reference, iter)
    #!NOTE: We dont need to use more general ℓevidence here as evidence distribution independent of Particle type
    guess_evidence(iter) = kernel.evidence(reference, iter)
    ## Compute first iteration at which kernel works
    Mparticle = _checkmemory(guess_state)
    Mevidence = _checkmemory(guess_evidence)
    ## Check if dependency structure could be determined
    ArgCheck.@argcheck isa(Mparticle, Integer) "Could not determine particle dependency - check kernel.transition in dynamics, or consider writing an issue in BaytesFilters."
    ArgCheck.@argcheck isa(Mevidence, Integer) "Could not determine data dependency - check kernel.observation in dynamics, or consider writing an issue in BaytesFilters."
    ## Return ParticleFilterMemory - subtract 1 from memory as iterations start with 1, which would mean memoryless transition.
    #=
    println(
    "Latent memory set to: ", Mparticle - 1,
    ". Data memory set to: ", Mevidence - 1,
    ". Initial distribution sampling steps: ", max(1, Mparticle - 1),
    ". If either of those is not as intended, you can manually define the memory in the 'ParticleFilterDefault' container as
    ParticleFilterDefault(; memory = ParticleFilterMemory(latent, data, initial))")
    =#
    return ParticleFilterMemory(Mparticle - 1, Mevidence - 1, max(1, Mparticle - 1))
end

################################################################################
#export
export ParticleFilterMemory


text:
# PPO loss
function ppo_loss(π, 𝒫, 𝒟; info = Dict())
    new_probs = logpdf(π, 𝒟[:s], 𝒟[:a])
    r = exp.(new_probs .- 𝒟[:logprob])

    A = 𝒟[:advantage]
    p_loss = -mean(min.(r .* A, clamp.(r, (1f0 - 𝒫[:ϵ]), (1f0 + 𝒫[:ϵ])) .* A))
    e_loss = -mean(entropy(π, 𝒟[:s]))

    # Log useful information
    ignore() do
        info[:entropy] = -e_loss
        info[:kl] = mean(𝒟[:logprob] .- new_probs)
        info[:clip_fraction] = sum((r .> 1 + 𝒫[:ϵ]) .| (r .< 1 - 𝒫[:ϵ])) / length(r)
    end
    𝒫[:λp]*p_loss + 𝒫[:λe]*e_loss
end

function PPO(;π::ActorCritic,
     ϵ::Float32 = 0.2f0,
     λp::Float32 = 1f0,
     λe::Float32 = 0.1f0,
     target_kl = 0.012f0,
     a_opt::NamedTuple=(;),
     c_opt::NamedTuple=(;),
     log::NamedTuple=(;),
     required_columns=[],
     kwargs...)

     OnPolicySolver(;agent=PolicyParams(π),
                    𝒫=(ϵ=ϵ, λp=λp, λe=λe),
                    log = LoggerParams(;dir = "log/ppo", log...),
                    a_opt = TrainingParams(;loss = ppo_loss, early_stopping = (infos) -> (infos[end][:kl] > target_kl), name = "actor_", a_opt...),
                    c_opt = TrainingParams(;loss = (π, 𝒫, D; kwargs...) -> Flux.mse(value(π, D[:s]), D[:return]), name = "critic_", c_opt...),
                    post_sample_callback = (𝒟; kwargs...) -> (𝒟[:advantage] .= whiten(𝒟[:advantage])),
                    required_columns = unique([required_columns..., :return, :logprob, :advantage]),
                    kwargs...)
end

# PPO loss with a penalty
function lagrange_ppo_loss(π, 𝒫, 𝒟; info = Dict())
    new_probs = logpdf(π, 𝒟[:s], 𝒟[:a])
    r = exp.(new_probs .- 𝒟[:logprob])

    A = 𝒟[:advantage]
    p_loss = -mean(min.(r .* A, clamp.(r, (1f0 - 𝒫[:ϵ]), (1f0 + 𝒫[:ϵ])) .* A))
    e_loss = -mean(entropy(π, 𝒟[:s]))

    #update the cost penalty
    penalty = ignore() do
        # 𝒫[:penalty_param][1] = clamp(𝒫[:penalty_param][1], -7, 10)
        # Flux.softplus(𝒫[:penalty_param][1])

        # Average cost
        Jc = sum(𝒟[:cost]) / sum(𝒟[:episode_end])
        # Jc = maximum(𝒟[:cost])

        # Compute the error
        Δ = Jc - 𝒫[:target_cost]

        # Update integral term
        𝒫[:I][1] = clamp(𝒫[:I][1] + 𝒫[:Ki]*Δ, 0, 𝒫[:Ki_max])

        # Smooth out the values
        α = 𝒫[:ema_α]
        𝒫[:smooth_Δ][1] = α * 𝒫[:smooth_Δ][1] + (1 - α)*Δ
        𝒫[:smooth_Jc][1] = α * 𝒫[:smooth_Jc][1] + (1 - α)*Jc

        # Compute the derivative term
        ∂ = max(0, 𝒫[:smooth_Jc][1] - 𝒫[:Jc_prev][1])

        # Update the previous cost
        𝒫[:Jc_prev][1] = 𝒫[:smooth_Jc][1]

        # PID update
        penalty = clamp(𝒫[:Kp] * 𝒫[:smooth_Δ][1] + 𝒫[:I][1] + 𝒫[:Kd]*∂, 0, 𝒫[:penalty_max])

        info["penalty"] = penalty
        info["cur_cost"] = Jc
        info["prop_term"] = 𝒫[:Kp] * 𝒫[:smooth_Δ][1]
        info["deriv_term"] = ∂
        info["integral term"] = 𝒫[:I][1]

        penalty
    end

    # cost_loss = 𝒫[:penalty_scale] * penalty * mean(r .* 𝒟[:cost_advantage])
    cost_loss = penalty * mean(max.(r .* 𝒟[:cost_advantage], clamp.(r, (1f0 - 𝒫[:ϵ]), (1f0 + 𝒫[:ϵ])) .* 𝒟[:cost_advantage]))

    # Log useful information
    ignore() do
        info[:entropy] = -e_loss
        info[:kl] = mean(𝒟[:logprob] .- new_probs)
        info[:clip_fraction] = sum((r .> 1 + 𝒫[:ϵ]) .| (r .< 1 - 𝒫[:ϵ])) / length(r)
        info["p_loss"] = 𝒫[:λp]*p_loss
        info["cost_loss"] = cost_loss
    end
    (𝒫[:λp]*p_loss + 𝒫[:λe]*e_loss + cost_loss) / (1 + penalty)
end

function LagrangePPO(;π::ActorCritic,
     Vc::ContinuousNetwork, # value network for estimating cost
     ϵ::Float32 = 0.2f0,
     λp::Float32 = 1f0,
     λe::Float32 = 0.1f0,
     λ_gae = 0.95f0,
     target_kl = 0.012f0,
     target_cost = 0.025f0,
     penalty_scale = 1f0,
     penalty_max = Inf32,
     Ki_max = 10f0,
     Ki = 1f-3,
     Kp = 1,
     Kd = 0,
     ema_α = 0.95,
     a_opt::NamedTuple=(;),
     c_opt::NamedTuple=(;),
     cost_opt::NamedTuple=(;),
     log::NamedTuple=(;),
     required_columns=[],
     kwargs...)

     𝒫=(ϵ=ϵ, λp=λp, λe=λe,
        target_cost=target_cost,
        penalty_scale=penalty_scale,
        penalty_max=penalty_max,
        Ki_max=Ki_max,
        I = [0f0],
        Jc_prev = [0f0],
        Ki=Ki,
        Kp=Kp,
        Kd=Kd,
        ema_α=ema_α,
        smooth_Δ = [0f0],
        smooth_Jc = [0f0]
        )

     OnPolicySolver(;agent=PolicyParams(π),
                    𝒫=𝒫,
                    Vc=Vc,
                    log = LoggerParams(;dir = "log/lagrange_ppo", log...),
                    a_opt = TrainingParams(;loss = lagrange_ppo_loss, early_stopping = (infos) -> (infos[end][:kl] > target_kl), name = "actor_", a_opt...),
                    c_opt = TrainingParams(;loss = (π, 𝒫, D; kwargs...) -> Flux.mse(value(π, D[:s]), D[:return]), name = "critic_", c_opt...),
                    cost_opt = TrainingParams(;loss = (π, 𝒫, D; kwargs...) -> Flux.mse(value(π, D[:s]), D[:cost_return]), name = "cost_critic_", cost_opt...),
                    required_columns = unique([required_columns..., :return, :advantage, :logprob, :cost_advantage, :cost, :cost_return]),
                    post_batch_callback = (𝒟; kwargs...) -> (𝒟[:advantage] .= whiten(𝒟[:advantage])),
                    kwargs...)
end


text:
import SDtoolbox: ct, cvsolve, CANTERA_MECH_FILETYPE
using PyCall
using StatsBase
using Interpolations

P₁ = 2e5
T₁ = 1500

ϕ = 1

X₁ = "H2:$(42*ϕ), O2:21,N2:79"
mech = "gri30."*CANTERA_MECH_FILETYPE

gas = ct.Solution(mech)
gas.TPX = T₁,P₁,X₁

out = cvsolve(gas,t_end=1e-3)

# plot!(out["time"],out["T"])

# compare with original sdtoolbox
pyimport_conda("scipy","scipy")
pyimport_conda("matplotlib","matplotlib")
pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
cvsolve_py = pyimport("sdtoolbox.cv").cvsolve
gas.TPX = T₁,P₁,X₁
out_py = cvsolve_py(gas,t_end=out["time"][end])
# plot!(out_py["time"],out_py["T"])

# calculate cross correlation to quantify error
itp = LinearInterpolation(out["time"][1:(end-1)], out["T"][1:(end-1)])
sol_itp = itp(out_py["time"])

err = maximum(crosscor(sol_itp,out_py["T"]))


text:
export TransferFunction, sampleTF, setTF

mutable struct TransferFunction
  freq::Vector{Float64}
  data::Matrix{ComplexF64}
  inductionFactor::Vector{Float64}

  function TransferFunction(freq_, datain::Array{T}, inductionFactor=ones(size(datain,2))) where {T<:Complex}
    freq = freq_[1]:(freq_[2]-freq_[1]):freq_[end]
    data=reshape(deepcopy(datain),size(datain,1), size(datain,2))
    return new(freq_, data, inductionFactor)
  end
end

function TransferFunction(freq_, ampdata, phasedata, args...)
  data = ampdata.*exp.(im.*phasedata)
  return TransferFunction(freq_, data, args...)
end

function TransferFunction(filename::String; kargs...)
    filenamebase, ext = splitext(filename)
    if  ext == ".h5"
      tf = load_tf(filename)
    else #if ext == "s1p" || ext == "s2p"
      tf = load_tf_fromVNA(filename; kargs...)
    end
    return tf
end

function getindex(tmf::TransferFunction, x::UnitRange, chan::Integer=1)
  a = tmf.data[x,chan]
  return a
end

#function getindex(tmf::TransferFunction, x::Real, chan::Integer=1)
#  a = tmf.interp[chan](x)
#  return a
#end

function getindex(tmf::TransferFunction, X::Union{Vector,AbstractRange}, chan::Integer=1)
  I = extrapolate(interpolate((tmf.freq,), tmf.data[:,chan], Gridded(Linear())), Interpolations.Flat())

  return [I(x) for x in X]
end

function getindex(tmf::TransferFunction, X::Union{Vector,AbstractRange}, chan::AbstractRange)
  out = zeros(ComplexF64, length(X), length(chan))
  for d=1:length(chan)
    out[:,d] = tmf[X,d]
  end
  return out
end

function load_tf(filename::String)
  tf = h5read(filename,"/transferFunction")
  freq = h5read(filename,"/frequencies")
  inductionFactor = h5read(filename,"/inductionFactor")
  return TransferFunction(freq,tf,inductionFactor)
end

function combine(tf1, tf2)
  freq = tf1.freq
  data = cat(tf1.data,tf2.data, dims=2)
  inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1)
  return TransferFunction(freq, data, inductionFactor)
end

function save(filename::String, tf::TransferFunction)
  h5write(filename, "/transferFunction", tf.data)
  h5write(filename, "/frequencies", tf.freq)
  h5write(filename, "/inductionFactor", tf.inductionFactor)
  return nothing
end

function load_tf_fromVNA(filename::String;
    frequencyWeighting=false,
    R = 50.0, #Ω
    N = 10, #5# Turns
    A = 7.4894*10.0^-4) #m^2 #1.3e-3^2*pi;)

  file = open(filename)
  lines = readlines(file)
  apdata = Float64[]
  aϕdata = Float64[]
  freq = Float64[]
  for i=4:length(lines)
      tmp = split(strip(lines[i])," ")
      tmp=tmp[tmp.!=""]
      f = parse(Float64,strip(tmp[1]))
      if occursin(lines[3],"# kHz S MA R 50")
          ap = parse(Float64,strip(tmp[2]))
          aphi = parse(Float64,strip(tmp[3]))
          f=f*1000
       elseif occursin(lines[3],"# Hz S DB R 50")
            ap = 10^(parse(Float64,strip(tmp[2]))/20)
            aphi = parse(Float64,strip(tmp[3]))*pi/180
            f=f
       elseif occursin(lines[3],"# kHz S DB R 50")
            ap = 10^(parse(Float64,strip(tmp[2]))/20)
            aphi = parse(Float64,strip(tmp[3]))*pi/180
            f=f*1000
        elseif occursin(lines[3],"# MHz S DB R 50")
            ap = 10^(parse(Float64,strip(tmp[2]))/20)
            aphi = parse(Float64,strip(tmp[3]))*pi/180
            f=f*1000000
      elseif occursin(lines[3],"# kHz S RI R 50")
          tf_complex=parse(Float64,strip(tmp[2]))+im*parse(Float64,strip(tmp[3]))
          ap=abs.(tf_complex);
          aphi=angle(tf_complex)
          f=f*1000
      else
	      error("Wrong data Format! Please export in kHz domain S21 parameter with either Magnitude/Phase, DB/Phase or Real/Imaginary!")
      end
      push!(apdata, ap)
      push!(aϕdata, aphi)
      push!(freq, f)
  end
  close(file)
  if frequencyWeighting
  	apdata ./= (freq.*2*pi) # As TF is defined as u_ADC = u_coil *TF the derivative from magnetic moment is applied as the division of the TF by w
  end
  return TransferFunction(freq, apdata, aϕdata)
end

function sampleTF(tmf::TransferFunction, f::MPIFile)
  freq = rxFrequencies(f)
  numChan = rxNumChannels(f)
  numFreq = length(freq)
  return tmf[freq,1:numChan]
end

function _writeBrukerParamFile(b::BrukerFile, A, filename)
  mv(joinpath(filepath(b),filename), joinpath(filepath(b),filename*"_bac"), force=true )

  open(joinpath(filepath(b),filename),"w") do fd
    for u in A
      println(fd,u)
    end
  end
  return
end

function setTF(b::BrukerFile, filenameTF::AbstractString)
  A = readlines(joinpath(filepath(b),"acqp"))
  l = findfirst( s->occursin("ACQ_comment", s), A)

  if l == nothing
   ll = findfirst( s->occursin("ACQ_experiment_mode", s), A)
   insert!(A, ll, "##\$ACQ_comment=( 2048 )")
   insert!(A, ll+1, "<$(filenameTF)>")
  else
    A[l+1] = "<$(filenameTF)>"
  end

  _writeBrukerParamFile(b, A, "acqp")

  # For some reason ACQ_comment is also in the method file -> also change that
  B = readlines(joinpath(filepath(b),"method"))
  g = findfirst( s->occursin("ACQ_comment", s), B)
  if g != nothing
    B[g+1] = "<$(filenameTF)>"
    _writeBrukerParamFile(b, B, "method")
  end

  if isfile(joinpath(filepath(b),"mdf"))
    MPIFile(filepath(b)) do f
      setTF(f, filenameTF)
    end
  end
  return
end

function setTF(f::MDFFile, filenameTF::AbstractString)
  tmf = TransferFunction(filenameTF)
  tf = sampleTF(tmf, f)

  # We need to close the HDF5 file handle before we can write to it
  close(f.file)

  h5open(filepath(f), "r+") do file
	if haskey(file, "/acquisition/receiver/transferFunctionFileName")
      delete_object(file, "/acquisition/receiver/transferFunctionFileName")
    end
    write(file, "/acquisition/receiver/transferFunctionFileName", filenameTF)
    if haskey(file, "/acquisition/receiver/transferFunction")
      delete_object(file, "/acquisition/receiver/transferFunction")
    end
    write(file, "/acquisition/receiver/transferFunction", tf)
    if haskey(file, "/acquisition/receiver/inductionFactor")
      delete_object(file, "/acquisition/receiver/inductionFactor")
    end
    write(file, "/acquisition/receiver/inductionFactor", tmf.inductionFactor)
  end
  return
end


text:
#
# Provides Restricted Two-body Problem equations
#

"""
Returns the conic section, as specified by eccentricity `e`.
"""
function conic(e::T) where T<:Real

    !isnan(e) || throw(ArgumentError("Provided eccentricity is not a number (NaN)."))

    if e ≈ zero(T)
        return Circular
    elseif e ≈ one(T)
        return Parabolic
    elseif zero(T) < e && e < one(T)
        return Elliptical
    elseif e > one(T)
        return Hyperbolic
    else
        throw(ArgumentError("Provided eccentricity, $e, is not valid."))
    end

end

"""
Returns a Keplarian representation of a Cartesian orbital state.
Algorithm taught in ENAE601.
"""
function keplerian(rᵢ, vᵢ, μ)

    safe_acos(num) = isapprox(num, one(num)) ?
                            acos(one(num)) :
                                isapprox(num, -one(num)) ?
                                    acos(-one(num)) :
                                        acos(num)

    î = SVector{3}(1, 0, 0)
    ĵ = SVector{3}(0, 1, 0)
    k̂ = SVector{3}(0, 0, 1)

    h̅ = specific_angular_momentum_vector(rᵢ, vᵢ)
    a = semimajor_axis(norm(rᵢ), norm(vᵢ), μ)
    n̅ = k̂ × specific_angular_momentum_vector(rᵢ, vᵢ)
    e̅ = eccentricity_vector(rᵢ, vᵢ, μ)
    e = norm(e̅)

    i = safe_acos((h̅ ⋅ k̂) / norm(h̅))

    Ω = ustrip(n̅ ⋅ ĵ) > 0 ?
            safe_acos((n̅ ⋅ î) / norm(n̅)) :
            2π - safe_acos((n̅ ⋅ î) / norm(n̅))

    ω = ustrip(e̅ ⋅ k̂) > 0 ?
            safe_acos((n̅ ⋅ e̅) / (norm(n̅) * e)) :
            2π - safe_acos((n̅ ⋅ e̅) / (norm(n̅) * e))

    ν = ustrip(rᵢ ⋅  vᵢ) > 0 ?
            safe_acos((e̅ ⋅ rᵢ) / (e * norm(rᵢ))) :
            2π - safe_acos((e̅ ⋅ rᵢ) / (e * norm(rᵢ)))

    return upreferred(e), upreferred(a), upreferred(i),
           upreferred(Ω), upreferred(ω),
           upreferred(ν)

end

"""
Returns a Cartesian representation of a Keplerian two-body orbital state
in an inertial frame, centered at the center of mass of the central body.
Algorithm taught in ENAE601.
"""
function cartesian(e, a, i, Ω, ω, ν, μ)
    rᵢ, vᵢ = spatial(i, Ω, ω, perifocal(a, e, ν, μ)...)
    return  upreferred.(rᵢ),
            upreffered.(vᵢ)
end

"""
Returns a spatial representation of the provied Perifocal state.
"""
function spatial(i, Ω, ω, rₚ, vₚ)

    # Set up Perifocal ⟶ Cartesian conversion
    R_3Ω =  SMatrix{3,3}(
            [cos(Ω)          -sin(Ω)            0.;
             sin(Ω)           cos(Ω)            0.;
             0.               0.                1.])
    R_1i = SMatrix{3,3}(
            [1.               0.                0.;
             0.               cos(i)           -sin(i);
             0.               sin(i)            cos(i)])
    R_3ω = SMatrix{3,3}(
            [cos(ω)          -sin(ω)            0.
             sin(ω)           cos(ω)            0.
             0.               0.                1.])

    ᴵTₚ = R_3Ω * R_1i * R_3ω

    return ᴵTₚ * rₚ, ᴵTₚ * vₚ

end

"""
Returns position and velocity vectors in the Perifocal frame.
"""
function perifocal(a, e, ν, μ)

        p = semi_parameter(a, e)
        r = distance(p, e, ν)

        P̂=SVector{3, Float64}(1, 0, 0)
        Q̂=SVector{3, Float64}(0, 1, 0)
        # Ŵ=SVector{3, Float64}(0, 0, 1)

        rₚ = (r * cos(ν) .* P̂ .+ r * sin(ν) .* Q̂)
        vₚ = √(μ/p) * ((-sin(ν) * P̂) .+ ((e + cos(ν)) .* Q̂))

        return rₚ, vₚ

end

function perifocal(i, Ω, ω, rᵢ, vᵢ)

    # Set up Cartesian ⟶ Perifocal conversion
    R_3Ω =  SMatrix{3,3}(
            [cos(Ω)          -sin(Ω)            0.;
             sin(Ω)           cos(Ω)            0.;
             0.               0.                1.])
    R_1i = SMatrix{3,3}(
            [1.               0.                0.;
             0.               cos(i)           -sin(i);
             0.               sin(i)            cos(i)])
    R_3ω = SMatrix{3,3}(
            [cos(ω)          -sin(ω)            0.
             sin(ω)           cos(ω)            0.
             0.               0.                1.])

    ᵖTᵢ = inv(R_3Ω * R_1i * R_3ω)

    return ᵖTᵢ*rᵢ, ᵖTᵢ*vᵢ

end

"""
Returns semimajor axis parameter, a.
"""
semimajor_axis(r, v, μ) = inv( (2 / r) - (v^2 / μ) )

"""
Returns specific angular momentum vector, h̅.
"""
specific_angular_momentum_vector(rᵢ, vᵢ) = rᵢ × vᵢ

"""
Returns scalar specific angular momentum vector, h.
"""
specific_angular_momentum(rᵢ, vᵢ) = norm(specific_angular_momentum_vector(rᵢ, vᵢ))

"""
Returns specific orbital energy, ϵ.
"""
specific_energy(a, μ) = ( -μ / (2 * a) )
specific_energy(r, v, μ) = (v^2 / 2) - (μ / r)

"""
Returns C3 value.
"""
C3(r, v, μ) = v^2 - 2μ/r

"""
Returns v∞.
"""
v_infinity(r, v, μ) = √C3(r, v, μ)

"""
Returns potential energy for an orbit about a `RestrictedTwoBodySystem`.
"""
specific_potential_energy(r, μ) = (μ/r)
specific_potential_energy(r, μ, R, J₂, ϕ) = (μ/r) * (1 - J₂ * (R/r)^2 * ((3/2) * (sin(ϕ))^2 - (1/2)))

"""
Returns orbital eccentricity vector e̅.
"""
function eccentricity_vector(rᵢ, vᵢ, μ)

    return map(x-> abs(x) < eps(typeof(x)) ? zero(x) : x, (1 / μ) * ((vᵢ × specific_angular_momentum_vector(rᵢ, vᵢ)) - μ * rᵢ / norm(rᵢ)))

end

"""
Returns orbital eccentricity, e.
"""
eccentricity(rᵢ, vᵢ, μ) = norm(eccentricity_vector(rᵢ, vᵢ, μ)) |> upreferred

"""
Returns semilatus parameter, p.
"""
semi_parameter(a, e) = a * (1 - e^2)

"""
Returns distance, r.
"""
distance(p, e, ν) = upreferred(p / (1 + e * cos(ν)))

"""
Returns instantaneous velocity, v, for any orbital representation.
"""
speed(r, a, μ) =  upreferred(√( (2 * μ / r) - (μ / a)))

"""
Returns the instantaneous velocity, `v`, for any orbital representation.
"""
speed(p, e, ν, a, μ) = speed(distance(p,e,ν), a, μ)

"""
Returns periapsis distance, rₚ.
"""
periapsis_radius(a, e) = a * (1 - e)

"""
Returns apoapsis distance, rₐ.
"""
apoapsis_radius(a, e) = a * (1 + e)

"""
Returns the orbital period.
"""
period(a, μ) = 2π * √(upreferred(a^3 / μ))

"""
Returns true anomoly, ν.
"""
function true_anomoly(r, h, e, μ)
    val = (h^2 - μ * r) / (μ * r * e)
    acos(isapprox(val, one(val)) ? one(val) : val)
end

"""
Returns mean motion, n.
"""
mean_motion(a, μ) = √(μ / a^3)

"""
Returns time since periapsis, t.
"""
time_since_periapsis(n, e, E) = (E - e * sin(E)) / (n)

"""
Sphere of influence.
"""
SOI(a, m, M) = a * (m / M)^(2/5)

"""
Sphere of activity.
"""
SOA(a, m, M) = a * (m / 3M)^(1/3)

"""
Computes a Hohmann transfer, and returns the departure and
arrival velocity vectors.
"""
function hohmann(r₁, r₂, μ)

    vₐ = √((2μ/r₁) - (2μ/(r₁+r₂)))
    vₚ = √((2μ/r₂) - (2μ/(r₁+r₂)))

    return vₐ, vₚ

end


text:
zerosum(A::AbstractArray; dims = 1) = A .- mean(A; dims)
zerosum!(A::AbstractArray; dims = 1) = A .= zerosum(A; dims)

"""
    zerosum(rbm)

Returns an equivalent `rbm` in zerosum gauge. Only affects Potts layers. If the `rbm`
doesn't have `Potts` layers, does nothing.
"""
function zerosum(rbm::RBM)
    vdims = ntuple(identity, ndims(rbm.visible))
    hdims = ntuple(d -> d + ndims(rbm.visible), ndims(rbm.hidden))
    if rbm.visible isa Potts && rbm.hidden isa Potts
        ωv = mean(rbm.w; dims = 1)
        ωh = mean(rbm.w; dims = 1 + ndims(rbm.visible))
        ω = mean(rbm.w; dims = (1, 1 + ndims(rbm.visible)))
        visible = Potts(rbm.visible.θ .- mean(rbm.visible.θ; dims = 1) .+ reshape(sum(ωh .- ω; dims = hdims), size(rbm.visible)))
        hidden = Potts(rbm.hidden.θ .- mean(rbm.hidden.θ; dims = 1) .+ reshape(sum(ωv .- ω; dims = vdims), size(rbm.hidden)))
        return oftype(rbm, RBM(visible, hidden, rbm.w .- ωv .- ωh .+ ω))
    elseif rbm.visible isa Potts
        ωv = mean(rbm.w; dims = 1)
        visible = Potts(rbm.visible.θ .- mean(rbm.visible.θ; dims = 1))
        hidden = shift_fields(rbm.hidden, reshape(sum(ωv; dims = vdims), size(rbm.hidden)))
        return oftype(rbm, RBM(visible, hidden, rbm.w .- ωv))
    elseif rbm.hidden isa Potts
        ωh = mean(rbm.w, dims = ndims(rbm.visible) + 1)
        visible = shift_fields(rbm.visible, reshape(sum(ωh; dims = hdims), size(rbm.visible)))
        hidden = Potts(rbm.hidden.θ .- mean(rbm.hidden.θ; dims = 1))
        return oftype(rbm, RBM(visible, hidden, rbm.w .- ωh))
    else
        # if the RBM doesn't have Potts layers, do nothing
        return rbm
    end
end

shift_fields(l::Binary, a::AbstractArray) = Binary(l.θ .+ a)
shift_fields(l::Spin, a::AbstractArray) = Spin(l.θ .+ a)
shift_fields(l::Potts, a::AbstractArray) = Potts(l.θ .+ a)
shift_fields(l::Gaussian, a::AbstractArray) = Gaussian(l.θ .+ a, l.γ)
shift_fields(l::ReLU, a::AbstractArray) = ReLU(l.θ .+ a, l.γ)
shift_fields(l::dReLU, a::AbstractArray) = dReLU(l.θp .+ a, l.θn .+ a, l.γp, l.γn)
shift_fields(l::pReLU, a::AbstractArray) = pReLU(l.θ .+ a, l.γ, l.Δ, l.η)
shift_fields(l::xReLU, a::AbstractArray) = xReLU(l.θ .+ a, l.γ, l.Δ, l.ξ)

function shift_fields!(l::Union{Binary,Spin,Potts,Gaussian,ReLU,pReLU,xReLU}, a::AbstractArray)
    l.θ .+= a
    return l
end

function shift_fields!(l::dReLU, a::AbstractArray)
    l.θp .+= a
    l.θn .+= a
    return l
end

"""
    zerosum!(rbm)

In-place zero-sum gauge on `rbm`.
"""
function zerosum!(rbm::RBM)
    if rbm.visible isa Potts
        zerosum!(rbm.visible.θ; dims = 1)
        ωv = mean(rbm.w; dims = 1)
        rbm.w .-= ωv
        dims = ntuple(identity, ndims(rbm.visible))
        shift_fields!(rbm.hidden, reshape(sum(ωv; dims), size(rbm.hidden)))
    end
    if rbm.hidden isa Potts
        zerosum!(rbm.hidden.θ; dims = 1)
        ωh = mean(rbm.w; dims = 1 + ndims(rbm.visible))
        rbm.w .-= ωh
        dims = ntuple(d -> d + ndims(rbm.visible), ndims(rbm.hidden))
        shift_fields!(rbm.visible, reshape(sum(ωh; dims), size(rbm.visible)))
    end
    return rbm
end

"""
    zerosum!(∂, rbm)

Projects the gradient so that it doesn't modify the zerosum gauge.
"""
function zerosum!(∂::NamedTuple, rbm::RBM)
    if rbm.visible isa Potts
        zerosum!(∂.visible.θ; dims = 1)
        zerosum!(∂.w; dims = 1)
    end
    if rbm.hidden isa Potts
        zerosum!(∂.hidden.θ; dims = 1)
        zerosum!(∂.w; dims = ndims(rbm.visible) + 1)
    end
    return ∂
end

function zerosum_weights(weights::AbstractArray, rbm::RBM)
    @assert size(weights) == size(rbm.w)
    if rbm.visible isa Potts && rbm.hidden isa Potts
        ωv = mean(weights; dims = 1)
        ωh = mean(weights; dims = 1 + ndims(rbm.visible))
        ω = mean(weights; dims = (1, 1 + ndims(rbm.visible)))
        return oftype(weights, weights .- ωv .- ωh .+ ω)
    elseif rbm.visible isa Potts
        ωv = mean(weights; dims = 1)
        return oftype(weights, weights .- ωv)
    elseif rbm.hidden isa Potts
        ωh = mean(weights, dims = ndims(rbm.visible) + 1)
        return oftype(weights, rbm.w .- ωh)
    else
        # if the RBM doesn't have Potts layers, do nothing
        return weights
    end
end


text:
module ServiceMeshLearning

export create_experiment, run_experiment, create_env, create_wrapped_env, create_agent

using Dates
using Distributions
using Flux
using ReinforcementLearning
using Logging, TensorBoardLogger
using Random, StableRNGs
using IntervalSets
using Statistics
using BSON

using ServiceMesh

include("agents/guided_policy.jl") # Hooks need this and agents need hooks...
include("hooks.jl")
include("agents/agents.jl")
include("env.jl")

function create_experiment(;
    alg::Symbol = :SAC,
    tag::String = "test",
    basepath::String = joinpath(homedir(), "servicemesh_results"),
    timesteps::Int = 1_000_000,
    n_env = 1,
    log_every = 60,
    logging = true,
    running_reward_steps = 3600*24*2, # 2 days?
    verbose=true,
    env=:default,
    kwargs...
)
    wrapped_env = default_wrappers(create_env(; env, kwargs...); kwargs...)
    if alg in [:SAC, :Guided, :LoadedAgent]
        wrapped_env = rl_wrappers(default_wrappers(wrapped_env; kwargs...); kwargs...)
    end

    agent = create_agent(Val(alg); env=wrapped_env, n_env=n_env, kwargs...)

    hooks = AbstractHook[RunningReward(running_reward_steps; envmap = env -> env[!])]
    if logging
        curr_time = Dates.format(now(), "yymmdd_HHMMSS")
        save_dir = mkpath(joinpath(basepath, "$(env)", tag, "$(alg)", "$curr_time"))
        verbose && @show save_dir
        push!(hooks, LogEveryNStep(save_dir; n=log_every))
        push!(hooks, DoOnExit(
            function(agent, env)
                # Log all hyperparam data
                open("$(save_dir)/params.txt", "w") do io
                    params = (; alg, tag, timesteps, n_env, log_every, logging, running_reward_steps, verbose, env, kwargs...)
                    join(io, map(k -> "$(k) = $(params[k])", keys(params)), "\n")
                end
                # Log policy if agent has policy field
                if hasproperty(agent, :policy)
                    BSON.@save "$(save_dir)/policy.bson" policy=agent.policy
                end
            end
        ))
    end
    hooks = ComposedHook(hooks...)

    stop_condition = StopAfterStep(timesteps; is_show_progress=verbose)

    Experiment(agent, wrapped_env, stop_condition, hooks, "## ServiceMesh $(tag)")
end

function run_experiment(; verbose=true, kwargs...)
    e = create_experiment(; verbose=verbose, kwargs...)
    run(e; describe=verbose)
    return mean(e.hook.hooks[1].rewards)
end

end # module

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
__precompile__()

module Nabla
    using ChainRules
    using ChainRulesCore
    using ChainRulesOverloadGeneration
    using ExprTools: ExprTools
    using ForwardDiff: ForwardDiff
    using LinearAlgebra
    using Random
    using SpecialFunctions
    using Statistics

    # Some aliases used repeatedly throughout the package.
    export ∇Scalar, ∇Array, SymOrExpr, ∇ArrayOrScalar
    const ∇Scalar = Number
    const ∇Array = AbstractArray{<:∇Scalar}
    const ∇AbstractVector = AbstractVector{<:∇Scalar}
    const ∇AbstractMatrix = AbstractMatrix{<:∇Scalar}
    const ∇ArrayOrScalar = Union{AbstractArray{<:∇Scalar}, ∇Scalar}
    const SymOrExpr = Union{Symbol, Expr}

    # ones/zeros(::AbstractArray) is deprecated in 0.7 and removed in 1.0, but it's a
    # pretty useful method, so we'll define our own for internal use
    for f in (:ones, :zeros)
        like = Symbol(f, "like")
        @eval begin
            $(like)(a::AbstractArray) = $(f)(eltype(a), size(a))
            $(like)(n::Integer) = $(f)(n)
        end
    end

    # Meta-programming utilities specific to Nabla.
    include("code_transformation/util.jl")
    include("code_transformation/differentiable.jl")

    # Functionality for constructing computational graphs.
    include("core.jl")

    # Functionality for defining new sensitivities.
    include("sensitivity.jl")

    # Finite differencing functionality - only used in tests. Would be good to move this
    # into a separate module at some point.
    include("finite_differencing.jl")

    # Sensitivities via ChainRules
    include("sensitivities/chainrules.jl")

    # Overloads we added to solve ambiguitities (probably caused by ChainRules)
    include("sensitivities/ambiguity_resolution.jl")

    # Sensitivities for the basics.
    include("sensitivities/indexing.jl")
    include("sensitivities/scalar.jl")

    # Sensitivities for functionals.
    include("sensitivities/functional/functional.jl")
    include("sensitivities/functional/reduce.jl")
    include("sensitivities/functional/reducedim.jl")

    # Linear algebra optimisations.
    include("sensitivities/linalg/generic.jl")
    include("sensitivities/linalg/symmetric.jl")
    include("sensitivities/linalg/blas.jl")
    include("sensitivities/linalg/diagonal.jl")
    include("sensitivities/linalg/factorization/cholesky.jl")
    include("sensitivities/linalg/factorization/svd.jl")

    # Checkpointing
    include("checkpointing.jl")

    # Link up to ChainRulesCore so rules are generated when new rrules are declared.
    on_new_rule(generate_overload, rrule)

end # module Nabla


text:
export GitRevWalker, hide!, reset!, walk, sortby!

type GitRevWalker
    repo::GitRepo
    ptr::Ptr{Void}

    function GitRevWalker(repo::GitRepo, ptr::Ptr{Void})
        @assert ptr != C_NULL
        w = new(repo, ptr)
        finalizer(w, free!)
        return w
    end
end

GitRevWalker(r::GitRepo) = begin
    wptr = Ptr{Void}[0]
    @check ccall((:git_revwalk_new, libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}), wptr, r)
    return GitRevWalker(r, wptr[1])
end

free!(w::GitRevWalker) = begin
    if w.ptr != C_NULL
        ccall((:git_revwalk_free, libgit2), Void, (Ptr{Void},), w.ptr)
        w.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, w::GitRevWalker) = w.ptr

Base.start(w::GitRevWalker) = begin
    id_ptr = [Oid()]
    err = ccall((:git_revwalk_next, libgit2), Cint, (Ptr{Oid}, Ptr{Void}), id_ptr, w)
    if err == GitErrorConst.ITEROVER
        return (nothing, true)
    elseif err != GitErrorConst.GIT_OK
        throw(LibGitError(err))
    end
    return (lookup_commit(w.repo, id_ptr[1]), false)
end

Base.done(w::GitRevWalker, state) = state[2]::Bool

Base.next(w::GitRevWalker, state) = begin
    id_ptr = [Oid()]
    err = ccall((:git_revwalk_next, libgit2), Cint,
                (Ptr{Oid}, Ptr{Void}), id_ptr, w)
    if err == GitErrorConst.ITEROVER
        return (state[1], (nothing, true))
    elseif err != GitErrorConst.GIT_OK
        throw(LibGitError(err))
    end
    return (state[1], (lookup_commit(w.repo, id_ptr[1]), false))
end

Base.push!(w::GitRevWalker, cid::Oid) = begin
    @check ccall((:git_revwalk_push, libgit2), Cint,
                 (Ptr{Void}, Ptr{Oid}), w, &cid)
    return w
end

function sortby!(w::GitRevWalker, sort_mode::Symbol; rev::Bool=false)
    s = sort_mode === :none ? GitConst.SORT_NONE :
        sort_mode === :topo ? GitConst.SORT_TOPOLOGICAL :
        sort_mode === :date ? GitConst.SORT_TIME :
        throw(ArgumentError("unknown git sort flag :$s"))
    rev && (s |= GitConst.SORT_REVERSE)
    sortby!(w, s)
    return
end

function sortby!(w::GitRevWalker, sort_mode::Cint)
    ccall((:git_revwalk_sorting, libgit2), Void,
           (Ptr{Void}, Cint), w, sort_mode)
    return w
end

function hide!(w::GitRevWalker, id::Oid)
    @check ccall((:git_revwalk_hide, libgit2), Cint,
                 (Ptr{Void}, Ptr{Oid}), w, &id)
    return w
end

function reset!(w::GitRevWalker)
    ccall((:git_revwalk_reset, libgit2), Void, (Ptr{Void},), w)
    return w
end

function walk(r::GitRepo, from::Oid, sort_mode::Symbol=:date, rev::Bool=false)
    walker = GitRevWalker(r)
    sortby!(walker, sort_mode, rev=rev)
    push!(walker, from)
    return (@task for c in walker
        produce(c)
    end)
end

# walk(repo, oid) do commit
      # do something with commit
# end
function walk(f::Function, r::GitRepo, from::Oid, sort_mode::Symbol=:date, rev::Bool=false)
    walker = GitRevWalker(r)
    sortby!(walker, sort_mode, rev=rev)
    push!(walker, from)
    for c in walker
        f(c)
    end
    return
end


text:
M(P) = [Int(sign(P[i]-P[j])) for i=1:length(P), j=1:length(P)]
preference(d::Dict{String,Int}) = [d[k] for k in sort([keys(d)...])]

"""
    ⦸(d₁, d₂[, dₛ...])

Return "light" preference matrix of given rankings.
A binary (infix) operation on two or more rankings.

# Examples
```jldoctest
julia> ranking"a₂,a₄,a₁,a₃" ⦸
       ranking"a₁,a₃∼a₄,a₂" ⦸
       ranking"a₂∼a₃,a₄,a₁" ⦸
       ranking"a₃,a₂,a₁∼a₄"
4×4 Array{Int64,2}:
  0   2  0   1
 -2   0  1  -2
  0  -1  0  -1
 -1   2  1   0

julia> ⦸(ranking"123",ranking"231",ranking"312")
3×3 Array{Int64,2}:
  0  -1   1
  1   0  -1
 -1   1   0
```
"""
function  ⦸(d₁::Dict{String,Int}, d₂::Dict{String,Int}, dₛ::Dict{String,Int}...)
    @assert length(d₁) == length(d₂)
    @assert sort([keys(d₁)...]) == sort([keys(d₂)...])
    R =  -M(preference(d₁))-M(preference(d₂))
    for d ∈ dₛ
        @assert size(R)[1] == size(R)[2] == length(d)
        @assert sort([keys(d₁)...]) == sort([keys(d)...])
        R = ⦸(R, d)
    end
    return R
end

⦸(R::Matrix{Int}, d::Dict{String,Int}) = R - M(preference(d))


text:
using Test
using FileIO, Random, Pkg
using MakieGallery
using AbstractPlotting, GLMakie
using Statistics

database = MakieGallery.load_test_database()
# filter!(database) do x
#     x.title == "Axis theming"
# end
# function has_no_recording(x)
#     isempty(readdir(joinpath(test_record_path, string(x.unique_name), "media")))
# end
# filter!(has_no_recording, database)

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")

isdir(tested_diff_path) && rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)

isdir(test_record_path) && rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)

examples = MakieGallery.record_examples(test_record_path)

@test length(examples) == length(database)

printstyled("Running ", color = :green, bold = true)
println("visual regression tests")

MakieGallery.run_comparison(test_record_path, tested_diff_path)
#
# cd(test_record_path) do
#   MakieGallery.generate_preview(test_record_path, "index.html")
# end
# MakieGallery.generate_thumbnails(test_record_path)
# MakieGallery.gallery_from_recordings(test_record_path, joinpath(test_record_path, "index.html"))


text:
function perform_wavelet_transf(f, Jmin, dir, filter = "9-7",separable = 0, ti = 0)
  """""
  perform_wavelet_transf - peform fast lifting transform

  y = perform_wavelet_transf(x, Jmin, dir, filter = "9-7",separable = 0, ti = 0);

  Implement 1D and 2D symmetric wavelets with symmetric boundary treatements, using
  a lifting implementation.

  filter gives the coefficients of the lifting filter.
  You can use h='linear' or h='7-9' to select automatically biorthogonal
  transform with 2 and 4 vanishing moments.

  You can set ti=1 to compute a translation invariant wavelet transform.

  You can set separable=1 to compute a separable 2D wavelet
  transform.

  Copyright (c) 2008 Gabriel Peyre
  """

  # Copy f.
  x = copy(f)

  # Convert Jmin to Int.
  Jmin = Int(Jmin)

  # Number of dimensions of x.
  d = ndims(x)

  if filter in ["linear","5-3"]
    h = [1/2, 1/4, sqrt(2)]

  elseif filter in ["9-7","7-9"]
    h = [1.586134342, -.05298011854, -.8829110762, .4435068522, 1.149604398]

  else
    error("Unknown filter")   # Raising an error if the input filter is unknown.

  end

  if (d == 2) & (separable == 1)
    ti = 0
    if ti == 1
      warn("Separable does not works for translation invariant transform")

      # perform a separable wavelet transform
    end

    n = size(x)[1]

    if dir == 1
      for i in 1:n
        x[:,i] = perform_wavelet_transf(x[:,i], Jmin, dir, filter, separable, ti)
      end
      for i in 1:n
        x[i,:] = transpose(perform_wavelet_transf(np.transpose(x[i,:]), Jmin, dir, filter, separable, ti))
      end

    else
      for i in 1:n
        x[i,:] = transpose(perform_wavelet_transf(np.transpose(x[i,:]), Jmin, dir, filter, separable, ti))
      end
      for i in 1:n
        x[:,i] = perform_wavelet_transf(x[:,i], Jmin, dir, filter, separable, ti)
      end
    end

    # number of lifting steps
    if ndims(x) == 1
      n = length(x)
    else
      n = size(x)[2]
    end
    m = Int((length(h) - 1)/2)
    Jmax = Int(log2(n) - 1)
    jlist = Jmax:-1:Jmin

    if dir == -1
      jlist = Jmin:Jmax
    end

    if ti == 0
      # subsampled
      for j in jlist
        if d == 1
          x[1:2^(j+1), :] = lifting_step(x[1:2^(j+1)], h, dir) # A vérifier les indices.
        else
          x[1:2^(j+1), 1:2^(j+1)] = lifting_step(x[1:2^(j+1), 1:2^(j+1)], h, dir) # A vérifier les indices.
          x[1:2^(j+1), 1:2^(j+1)] = transpose(lifting_step(transpose(x[1:2^(j+1), 1:2^(j+1)]), h, dir)) # A vérifier les indices.
        end
      end

    else
      # TI
      nJ = Jmax - Jmin + 1
      if (dir == 1) & (d == 1)
        x = repeat(x, outer = [nJ + 1, 1, 1])   # A vérifier.
      elseif (dir == 1) & (d == 2)
        x = repeat(x, outer = [3*nJ + 1,1,1])
      end
      for j in jlist
        dist = 2^(Jmax - j)

        if d == 1
          if dir == 1
            x[1:(j - Jmin + 2), :, :] = lifting_step_ti(x[0, :, :], h, dir, dist)
          else
            x[1, :, :] = lifting_step_ti(x[1:(j - Jmin + 2), :, :], h, dir, dist)
          end
        else
          dj = 3*(j - Jmin)

          if dir == 1
#              x[1 : dj + 1, :, :] = lifting_step_ti(x[1, :, :], h, dir, dist) # A vérifier avec précautions ligne 104
              x[[1, dj + 1], :, :] = lifting_step_ti(x[1, :, :], h, dir, dist)

#              x[1:dj + 2, :, :] = lifting_step_ti(transpose(x[1, :, :]), h, dir, dist) # A vérifier avec précautions ligne 106
              x[[1, dj + 2],:,:] = lifting_step_ti(transpose(x[1, :, :]), h, dir, dist)
              x[1, :, :] = transpose(x[1, :, :])
              x[dj + 3, :, :] = transpose(x[dj + 3, :, :])

#              x[2 + dj : 3 + dj, :, :] = lifting_step_ti(transpose(x[dj, :, :]), h, dir, dist)
              x[[2 + dj, 3 + dj], :, :] = lifting_step_ti(transpose(x[dj + 2, :, :]), h, dir, dist)
              x[dj + 2, :, :] = transpose(x[dj + 2, :, :])
              x[dj + 4, :, :] = transpose(x[dj + 4, :, :])
          else

              x[dj + 2, :, :] = transpose(x[dj + 2, :, :])
              x[dj + 4, :, :] = transpose(x[dj + 4, :, :])

#              x[dj + 2, :, :] = transpose(lifting_step_ti(x[[2 + dj : 3 + dj], :, :], h, dir, dist))
              x[dj + 2, :, :] = transpose(lifting_step_ti(x[[2 + dj, 3 + dj], :, :], h, dir, dist))

              x[1, :, :] = transpose(x[1, :, :])
              x[dj + 3, :, :] = transpose(x[dj + 3, :, :])
#              x[1, :, :] = transpose(lifting_step_ti(x[[1 : dj + 2], :, :], h, dir, dist))
              x[1, :, :] = transpose(lifting_step_ti(x[[1, dj+2],:,:], h, dir, dist))

#              x[1, :, :] = lifting_step_ti(x[[1, dj+1], :, :], h, dir, dist)
              x[1, :, :] = lifting_step_ti(x[[1 : dj+1], :, :], h, dir, dist)
            end
          end
        end
      end
      if dir == -1
        x = x[1, :, :]
      end
    end
    return x
end

function lifting_step(x0, h, dir)

    #copy x
    x = copy(x0)

    # number of lifting steps
    m = Int((length(h) - 1)/2)

    if dir == 1
        # split
#        d = x[1::2,]
        d = x[collect(i for i in 1:length(x) if i%2 == 0)]
#        x = x[0::2,]
        x = x[collect(i for i in 1:length(x) if i%2 == 1)]

        for i in 1 : m
            d = d - h[2*i] * (x + [x[1 : end, :]'; x[end, :]']) # I think we have to transpose.
            x = x + h[2*i + 1] * (d + [d[1, :]'; d[1 : end - 1,]']) #Be very cautiaus about transposing.
        end
        x = [x*h[end]; d/h[end]]

    else
        # retrieve detail coefs
        fin = length(x)
        d = x[fin/2 + 1 : end, :]*h[end]
        x = x[1:fin/2, :]/h[end]
        for i in m:-1:1
            x = x - h[2*i - 1] * (d + [d[1, :]; d[1:end - 1, :]])
            d = d + h[2*i - 2] * (x + [x[2:end, :]; x[end, :]])
        end
        # merge
        x1 = [x; x]
#        x1[::2,] = x
        x1[collect(i for i in 1:length(x[:, 1]) if i%2 == 1)] = x
#        x1[1::2,] = d
        x1[collect(i for i in 1:length(x[:, 1]) if i%2 == 0)] = d
        x = x1
    end

    return x
end

function lifting_step_ti(x0, h, dir, dist)

    #copy x
    x = copy(x0)

    # number of lifting steps
    m = Int((len(h) - 1)/2)
    n = size(x[0])[0]

    s1 = collect(i for i in 1:n) + dist
    s2 = collect(i for i in 1:n) - dist

    # boundary conditions
    s1[collect(i for i in 1:length(s1) if s1[i] > n)] = 2*n - s1[collect(i for i in 1:length(s1) if s1[i] > n)]
    s1[collect(i for i in 1:length(s1) if s1[i] < 1)] = 2 - s1[collect(i for i in 1:length(s1) if s1[i] < 1)]

    s2[collect(i for i in 1:length(s2) if s2[i] > n)] = 2*n - s2[collect(i for i in 1:length(s2) if s2[i] > n)]
    s2[s2 < 1] = 2   - s2[s2 < 1]
    s2[collect(i for i in 1:length(s2) if s2[i] < 1)] = 2 - s2[collect(i for i in 1:length(s2) if s2[i] < 1)]

    if dir == 1
        # split
        d = x
        for i in 1:m
            if ndims(x) == 2
                x = repeat(x, outer = [1, 1, 1])
            end
            d = d - h[2*i] * (x[:, s1, :] + x[:, s2, :])
            x = x + h[2*i + 1] * (d[:, s1, :] + d[:, s2, :])
        end

        #merge
#        x = np.concatenate((x*h[end],d/h[end]))
        x = [x*h[end]; d/h[end]]

    else
        # retrieve detail coefs

        d = x[1, :, :]*h[end]
        x = x[0, :, :]/h[end]

        for i in m:-1:1
            x = x - h[2*i - 1] * (d[s1, :] + d[s2, :])
            d = d + h[2*i - 2] * (x[s1, :] + x[s2, :])
        end

        # merge
        x = (x + d)/2
    end

    return x
end


text:
using RRTMGP.Gases
using RRTMGP.GasConcentrations
using RRTMGP.GasOptics
using RRTMGP.ReferenceStates
using NCDatasets

function read_char_vec(ds, var_name)
  tmp = ds[var_name][:]
  return String[strip(join(tmp[:,i])) for i = 1:size(tmp, 2)]
end

read_gases(ds, var_name) = convert.(AbstractGas, read_char_vec(ds, var_name))

"""
    load_and_init(ds, ::Type{FT}, gases_prescribed::Vector{AbstractGas}) where {FT<:AbstractFloat}

Initialize the gas optics class with data. The calls look slightly different depending
  on whether the radiation sources are internal to the atmosphere (longwave) or external (shortwave)
"""
function load_and_init(ds, ::Type{FT}, gases_prescribed::Vector{AbstractGas}) where {FT<:AbstractFloat}

  # Reading the properties from the NetCDF file

  I = Int
  kmajor                          = Array{FT}(ds["kmajor"][:])
  key_species                     = Array{I}(ds["key_species"][:])
  band2gpt                        = Array{I}(ds["bnd_limits_gpt"][:])
  band_lims_wavenum               = Array{FT}(ds["bnd_limits_wavenumber"][:])
  gases_in_database               = read_gases(ds, "gas_names")
  gas_minor                       = read_gases(ds, "gas_minor")
  identifier_minor                = read_gases(ds, "identifier_minor")

  rayl_lower = haskey(ds,"rayl_lower") ? Array{FT}(ds["rayl_lower"][:]) : nothing
  rayl_upper = haskey(ds,"rayl_upper") ? Array{FT}(ds["rayl_upper"][:]) : nothing
  @assert haskey(ds,"rayl_lower") == haskey(ds,"rayl_upper")

  lower = GasOpticsVars{FT,I}(Array{I}(ds["minor_limits_gpt_lower"][:]),
                              Array{Bool}(ds["minor_scales_with_density_lower"][:]),
                              Array{Bool}(ds["scale_by_complement_lower"][:]),
                              Array{I}(ds["kminor_start_lower"][:]),
                              Array{FT}(ds["kminor_lower"][:]),
                              read_gases(ds, "scaling_gas_lower"),
                              read_gases(ds, "minor_gases_lower"))
  upper = GasOpticsVars{FT,I}(Array{I}(ds["minor_limits_gpt_upper"][:]),
                              Array{Bool}(ds["minor_scales_with_density_upper"][:]),
                              Array{Bool}(ds["scale_by_complement_upper"][:]),
                              Array{I}(ds["kminor_start_upper"][:]),
                              Array{FT}(ds["kminor_upper"][:]),
                              read_gases(ds, "scaling_gas_upper"),
                              read_gases(ds, "minor_gases_upper"))

  ref = ReferenceState(Array{FT}(ds["press_ref"][:]),
                       Array{FT}(ds["temp_ref"][:]),
                       FT(ds["press_ref_trop"][:]),
                       Array{FT}(ds["vmr_ref"][:]),
                       gases_prescribed,
                       gases_in_database)

  optical_props = OpticalPropsBase("GasOptics optical props", band_lims_wavenum, band2gpt)

  args = (gases_prescribed,
          gases_in_database,
          key_species,
          optical_props,
          kmajor,
          gas_minor,
          identifier_minor,
          lower,
          upper)

  if haskey(ds,"totplnk")
    totplnk     = Array{FT}(ds["totplnk"][:])
    planck_frac = Array{FT}(ds["plank_fraction"][:])
    kdist = get_k_dist_lw(totplnk, planck_frac, rayl_lower, rayl_upper, ref, args...)
  else
    solar_src = ds["solar_source"][:]
    kdist = get_k_dist_sw(solar_src, rayl_lower, rayl_upper, ref, args...)
  end
  return kdist

end


text:
import BangBang
import OnlineStats
using Documenter
using Transducers
# Not importing `JSON` and `Literate` here so that this file can be
# included from test without adding them as dependencies.

EXAMPLE_PAGES = [
    "Tutorial: Missing values" => "tutorials/tutorial_missings.md",
    "Tutorial: Parallelism" => "tutorials/tutorial_parallel.md",
    "Parallel word count" => "tutorials/words.md",
    "Empty result handling" => "howto/empty_result_handling.md",
    "Writing transducers" => "howto/transducers.md",
    "Writing reducibles" => "howto/reducibles.md",
]

function transducers_rm_examples()
    rm(joinpath(@__DIR__, "src/tutorials"); force = true, recursive = true)
    rm(joinpath(@__DIR__, "src/howto"); force = true, recursive = true)
end

function transducers_literate(;
        inputbase = joinpath(@__DIR__, "..", "examples"),
        outputbase = joinpath(@__DIR__, "src"),
        examples = EXAMPLE_PAGES,
        kwargs...)
    transducers_rm_examples()
    for (_, outpath) in examples
        name, = splitext(basename(outpath))
        inputfile = joinpath(inputbase, "$name.jl")
        outputdir = joinpath(outputbase, dirname(outpath))
        Literate.markdown(
            inputfile, outputdir;
            documenter = true,
            kwargs...)
    end
end

function transducers_makedocs(;
        examples = EXAMPLE_PAGES,
        kwargs...)
    if isempty(examples)
        # Make some dummy examples
        examples = String[]
        for outdir in ["tutorials", "howto"]
            dummypage = joinpath(outdir, "dummy.md")
            path = joinpath(@__DIR__, "src", dummypage)
            mkpath(dirname(path))
            write(
                path,
                """
            # Dummy page
            ## Empty result handling
            """
            )
            push!(examples, dummypage)
        end
    end
    tutorials = filter(((_, path),) -> startswith(path, "tutorials/"), examples)
    howto = filter(((_, path),) -> startswith(path, "howto/"), examples)
    @assert issetequal(union(tutorials, howto), examples)
    makedocs(;
        modules = [Transducers],
        pages = [
            "Home" => "index.md",
            "Reference" => [
                "Manual" => "reference/manual.md",
                "Interface" => "reference/interface.md",
            ],
            "Tutorials" => tutorials,
            "How-to guides" => howto,
            "Explanation" => [
                "Parallelism" => "parallelism.md",  # TODO: merge this to index.md
                "Comparison to iterators" => "explanation/comparison_to_iterators.md",
                "Glossary" => "explanation/glossary.md",
                hide("Internals" => "explanation/internals.md"),
            ],
        ],
        repo = "https://github.com/JuliaFolds/Transducers.jl/blob/{commit}{path}#L{line}",
        sitename = "Transducers.jl",
        authors = "Takafumi Arakaki",
        root = @__DIR__,
        strict = true,
        kwargs...)
end

function transducers_redirection_mapping()
    mapping = [
        # old page => new page
        "manual" => "reference/manual",
        "interface" => "reference/interface",
    ]

    old_examples = [
        "tutorials/tutorial_missings",
        "tutorials/tutorial_parallel",
        "tutorials/words",
        "howto/empty_result_handling",
        "howto/transducers",
        "howto/reducibles",
    ]

    for e in old_examples
        push!(mapping, joinpath("examples", basename(e)) => e)
    end

    return mapping
end

function transducers_make_redirections(;
    mapping = transducers_redirection_mapping(),
    build = joinpath((@__DIR__), "build"),
)
    for (old, new) in mapping
        oldpath = joinpath(build, old)
        newpath = joinpath(build, new)
        relurl = relpath(newpath, oldpath)
        html = """<meta http-equiv="refresh" content="0; url=$relurl"/>"""
        mkpath(oldpath)
        write(joinpath(oldpath, "index.html"), html)
    end
end

function should_push_preview(event_path = get(ENV, "GITHUB_EVENT_PATH", nothing))
    event_path === nothing && return false
    event = JSON.parsefile(event_path)
    pull_request = get(event, "pull_request", nothing)
    pull_request === nothing && return false
    labels = [x["name"] for x in pull_request["labels"]]
    # https://developer.github.com/v3/activity/events/types/#pullrequestevent
    yes = "push_preview" in labels
    if yes
        @info "Trying to push preview as label `push_preview` is specified." labels
    else
        @info "Not pushing preview as label `push_preview` is not specified." labels
    end
    return yes
end


text:
#######################
## Timematr get mean ##
#######################

import Base.mean
function mean(tm::AbstractTimematr, dim::Int)
    ## output: DataFrame, since date dimension is lost
    if dim == 2
        error("For rowwise mean use rowmeans function")
    end
    meanVals = DataFrame()
    for (nam, col) in eachcol(tm.vals)
        meanVals[nam] = mean(col)
    end
    return meanVals
end

function mean(tn::AbstractTimenum, dim::Int)
    ## output: DataFrame, since date dimension is lost
    ## NAs are ignored!
    if dim == 2
        error("For rowwise mean use rowmeans function")
    end
    nObs, nAss = size(tn)
    meanVals = DataFrame()
    nVals = zeros(Int, nObs)
    for (nam, col) in eachcol(tn.vals)
        if isa(col, DataArray)
            meanVals[nam] = mean(col[!col.na])
        else
            meanVals[nam] = mean(col)
        end
    end
    return meanVals
end

## rowmeans
##---------

function rowmeans(tm::AbstractTimematr)
    ## output: Timematr
    nObs, nAss = size(tm)
    meanVals = zeros(Float64, nObs)
    for (nam, col) in eachcol(tm.vals)
        meanVals += col
    end
    return Timematr(DataFrame(mean_values = meanVals./nAss),
                    idx(tm))
end

function rowmeans(tn::Timenum)
    ## output: Timematr, NAs are ignored
    nObs, nAss = size(tn)
    meanVals = DataArray(zeros(Float64, nObs))
    nVals = zeros(Int, nObs)
    for (nam, col) in eachcol(tn.vals)
        if isa(col, DataArray)
            nVals[!col.na] += 1
            meanVals[!col.na] += col.data[!col.na]
        else
            nVals += 1
            meanVals += col
        end
    end
    return Timematr(DataFrame(mean_values = meanVals./nVals),
                    idx(tn))
end

#######################
## Timematr get prod ##
#######################

import Base.prod
function prod(tm::AbstractTimematr, dim::Int = 1)
    ## output: DataFrame, since date dimension is lost
    if dim == 2
        error("For rowwise prod use rowprods function")
    end
    prodVals = DataFrame()
    for (nam, col) in eachcol(tm.vals)
        prodVals[nam] = prod(col)
    end
    return prodVals
end

function rowprods(tm::AbstractTimematr)
    ## output: Timematr
    prodVals = prod(core(tm), 2)
    prods = Timematr(prodVals, idx(tm))
end
## function rowprods2(tm::AbstractTimematr)
##     ## output: Timematr
##     nObs, nAss = size(tm)
##     prodVals = ones(Float64, nObs)
##     for (nam, col) in eachcol(tm.vals)
##         prodVals = prodVals .* col
##     end
##     return Timematr(composeDataFrame(prodVals, [:prod_values]),
##                     idx(tm))
## end

######################################
## Timematr get row and column sums ##
######################################

import Base.sum
function sum(tm::AbstractTimematr, dim::Int = 1)
    ## output: DataFrame, since date dimension is lost
    if dim == 2
        error("For rowwise sum use rowsums function")
    end
    sumVals = sum(core(tm), dim)
    sums = composeDataFrame(sumVals, names(tm))
end

function rowsums(tm::AbstractTimematr)
    ## output: Timematr
    sumVals = sum(core(tm), 2)
    sums = Timematr(sumVals, idx(tm))
end

#########################
## Timematr covariance ##
#########################

import Base.cov
function cov(tm::AbstractTimematr)
    ## output: DataFrame
    covDf = composeDataFrame(cov(core(tm)), names(tm.vals))
    return covDf
end

##########################
## Timematr correlation ##
##########################

import Base.cor
function cor(tm::AbstractTimematr)
    ## output: DataFrame
    corDf = composeDataFrame(cor(core(tm)), names(tm.vals))
    return corDf
end

##########################
## Timematr std ##
##########################

import Base.std
function std(tm::AbstractTimematr)
    ## output: DataFrame
    stdDf = composeDataFrame(std(core(tm), 1), names(tm.vals))
    return stdDf
end

function std(tm::AbstractTimematr, dim::Integer)
    ## output: DataFrame
    if dim == 2
        error("For rowwise standard deviations use rowstds")
    end
    stdDf = composeDataFrame(std(core(tm), 1), names(tm.vals))
    return stdDf
end

function rowstds(tm::AbstractTimematr)
    rowStd = Timematr(std(core(tm), 2), idx(tm))
    return rowStd
end

#########################
## minimum and maximum ##
#########################

import Base.minimum
function minimum(tm::AbstractTimematr)
    return minimum(core(tm))
end
function minimum(tm::AbstractTimematr, dim::Integer)
    if dim == 2
        error("For rowwise minimum use rowmin function")
    end
    return composeDataFrame(minimum(core(tm), 1), names(tm))
end

############
## cumsum ##
############

import Base.cumsum
function cumsum(tm::AbstractTimematr, dim::Integer)
    cumulated = cumsum(core(tm), dim)
    return Timematr(composeDataFrame(cumulated, names(tm)), idx(tm))
end

import Base.cumprod
function cumprod(tm::AbstractTimematr, dim::Integer)
    cumulated = cumprod(core(tm), dim)
    return Timematr(composeDataFrame(cumulated, names(tm)), idx(tm))
end

####################
## geometric mean ##
####################

function geomMean(x::AbstractTimematr; percent = true)
    nObs = size(x, 1)
    if percent
        vals = (prod(1 + core(x)./100, 1).^(1/nObs) - 1)*100
    else
        vals = prod(1 + core(x), 1).^(1/nObs) - 1
    end
    return composeDataFrame(vals, names(x))
end

function geomMean(x; percent = true)
    nObs = size(x, 1)
    if percent
        vals = (prod(1 + x./100, 1).^(1/nObs) - 1)*100
    else
        vals = prod(1 + x./100, 1).^(1/nObs) - 1
    end
    return vals
end

#####################
## moving averages ##
#####################

function movAvg(tm::AbstractTimematr, nPeriods::Integer)
    (nObs, nAss) = size(tm)
    nMAs = nObs - nPeriods + 1
    vals = core(tm)

    ## get start and ending indices
    startInd = ceil(nPeriods/2)
    finInd = startInd + nMAs - 1

    ## preallocation
    movAvgs = ones(nMAs, nAss)

    movAvgs[1, :] = sum(vals[1:nPeriods, :], 1)

    for ii=2:nMAs
        movAvgs[ii, :] = movAvgs[ii-1, :] - vals[(ii-1), :] +
        vals[(ii-1)+nPeriods, :]
    end

    ## divide
    scaledVals = movAvgs / nPeriods
    dats = idx(tm)[startInd:finInd]

    movAvgTm = Timematr(scaledVals, names(tm), dats)
    return movAvgTm
end