text:
module SIMDDualNumbers

using VectorizationBase, SLEEFPirates, ForwardDiff
using VectorizationBase: AbstractSIMD, AbstractMask
import IfElse: ifelse

@generated function Base.abs(x::ForwardDiff.Dual{TAG,S,N}) where {TAG,S<:AbstractSIMD,N}
  quote
    $(Expr(:meta,:inline))
    val = x.value
    p = x.partials
    cmp = val < zero($S)
    absx = $ifelse(cmp, -val, val)
    Base.Cartesian.@nexprs $N n -> p_n = p[n]
    ForwardDiff.Dual{$TAG}(absx, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> $ifelse(cmp, -p_n, p_n)))
  end
end
@inline function Base.max(
  x::ForwardDiff.Dual{TAG,<:AbstractSIMD,N},
  y::ForwardDiff.Dual{TAG,<:AbstractSIMD,N}
) where {TAG,N}
  vx = ForwardDiff.value(x)
  vy = ForwardDiff.value(y)
  xgy = vx > vy
  z = ifelse(xgy, vx, vy)
  p = VectorizationBase.fmap(ifelse, xgy, ForwardDiff.partials(x).values, ForwardDiff.partials(y).values)
  ForwardDiff.Dual{TAG}(z, ForwardDiff.Partials(p))
end

@inline Base.max(x::T, y::Real) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = max(x, T(y))
@inline Base.max(y::Real, x::T) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = max(x, T(y))
@inline Base.max(x::T, y::Int) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = max(x, T(y))
@inline Base.max(y::Int, x::T) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = max(x, T(y))

@inline function Base.min(
  x::ForwardDiff.Dual{TAG,<:AbstractSIMD,N},
  y::ForwardDiff.Dual{TAG,<:AbstractSIMD,N}
) where {TAG,N}
  vx = ForwardDiff.value(x)
  vy = ForwardDiff.value(y)
  xgy = vx < vy
  z = ifelse(xgy, vx, vy)
  p = VectorizationBase.fmap(ifelse, xgy, ForwardDiff.partials(x).values, ForwardDiff.partials(y).values)
  ForwardDiff.Dual{TAG}(z, ForwardDiff.Partials(p))
end
@inline Base.min(x::T, y::Real) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = min(x, T(y))
@inline Base.min(y::Real, x::T) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = min(x, T(y))
@inline Base.min(x::T, y::Int) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = min(x, T(y))
@inline Base.min(y::Int, x::T) where {N,T<:ForwardDiff.Dual{<:Any,<:AbstractSIMD,N}} = min(x, T(y))

@generated function SLEEFPirates.tanh_fast(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
  quote
    $(Expr(:meta,:inline))
    t = tanh_fast(x.value)
    ∂t = $(VectorizationBase.vfnmadd_fast)(t, t, one(S))
    p = x.partials
    ForwardDiff.Dual{T}(t, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> $(Base.FastMath.mul_fast)(∂t, p[n])))
  end
end
@generated function SLEEFPirates.sigmoid_fast(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
  quote
    $(Expr(:meta,:inline))
    s = sigmoid_fast(x.value)
    ∂s =  $(VectorizationBase.vfnmadd_fast)(s,s,s)
    p = x.partials
    ForwardDiff.Dual{T}(s, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> $(Base.FastMath.mul_fast)(∂s, p[n])))
  end
end
@generated function VectorizationBase.relu(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
  quote
    $(Expr(:meta,:inline))
    v = x.value
    z = zero(v)
    cmp = v < z
    r = ifelse(cmp, z, v)
    p = x.partials
    ForwardDiff.Dual{T}(r, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> ifelse(cmp, z, p[n])))
  end
end

@generated function ifelse(m::AbstractMask, x::ForwardDiff.Dual{TAG,V,P}, y::ForwardDiff.Dual{TAG,V,P}) where {TAG,V,P}
  quote
    $(Expr(:meta,:inline))
    z = $ifelse(m, ForwardDiff.value(x), ForwardDiff.value(y))
    px = ForwardDiff.partials(x)
    py = ForwardDiff.partials(y)
    p = Base.Cartesian.@ntuple $P p -> $ifelse(m, px[p], py[p])
    ForwardDiff.Dual{$TAG}(z, ForwardDiff.Partials(p))
  end
end
@generated function ifelse(m::AbstractMask, x::Number, y::ForwardDiff.Dual{TAG,V,P}) where {TAG,V,P}
  quote
    $(Expr(:meta,:inline))
    z = $ifelse(m, x, ForwardDiff.value(y))
    py = ForwardDiff.partials(y)
    p = Base.Cartesian.@ntuple $P p -> $ifelse(m, zero($V), py[p])
    ForwardDiff.Dual{$TAG}(z, ForwardDiff.Partials(p))
  end
end
@generated function ifelse(m::AbstractMask, x::ForwardDiff.Dual{TAG,V,P}, y::Number) where {TAG,V,P}
  quote
    $(Expr(:meta,:inline))
    z = $ifelse(m, ForwardDiff.value(x), y)
    px = ForwardDiff.partials(x)
    p = Base.Cartesian.@ntuple $P p -> $ifelse(m, px[p], zero($V))
    ForwardDiff.Dual{$TAG}(z, ForwardDiff.Partials(p))
  end
end
@inline function SLEEFPirates.softplus(x::ForwardDiff.Dual{TAG}) where {TAG}
  val = ForwardDiff.value(x)
  expx = exp(val)
  vx = log1p(expx)
  px = Base.FastMath.inv_fast(one(val) + Base.FastMath.inv_fast(expx))
  ForwardDiff.Dual{TAG}(vx, Base.FastMath.mul_fast(ForwardDiff.partials(x),px))
end

end


text:
module Tree

export tree

if VERSION >= v"0.7.0-"
    using Dates
end

function colorprint(s, col)
    if VERSION >= v"0.7.0-"
        printstyled(s, color=col)
    else
        print_with_color(col, s)
    end
end

"""
    tree(path;
        maxdepth=6,
        linecolor=:light_black,
        dircolor=:yellow,
        ages=false)

Show the files in `path` in a basic tree, up to a depth of `maxdepth`.

The default value for maxdepth is 6 (which is way too big if you're high enough up in your
directory hierarchy). But the trouble with this option is that it doesn't show directories if
they're in the same directory as the lowest level files... (I think there are some
off-by-at-least one errors to track down.)

Also the scan takes place to deeper levels anyway, it just doesn't show them.

`linecolor` and `dircolor` are the colors for the connecting lines and directories.
The following work on my terminal:

```
:bold, :black, :blue, :cyan, :green, :light_black, :light_blue, :light_cyan, :light_green,
:light_magenta, :light_red, :light_yellow, :magenta, :red, :white, or :yellow
```

"""
function tree(path=".";
    maxdepth=6,
    linecolor=:light_black,
    dircolor=:yellow,
    ages=false)

    CORNER     = string(Char(0x2514)) # '└'
    HORIZONTAL = string(Char(0x2500)) # '─'

    # quit if not a directory
    !isdir(path) && error("Supply a directory")

    rootdepth = length(filter(!isempty, split(abspath(path), "/"))) - 1

    println("files in $(abspath(path)) | showing $maxdepth levels")

    for (root, dirs, files) in walkdir(abspath(path))

        depth = length(filter(!isempty, split(root, "/"))) - rootdepth

        d, f = splitdir(root)

        if depth > maxdepth
            continue
        else
            # print this directory:

            # leading spaces
            print(join(fill(" " ^ 4,  max(0, depth - 1))))

            # angle piece
            colorprint(string("$CORNER", join(fill("$HORIZONTAL" ^ 3, 1))), linecolor)

            # directory name
            colorprint(f * "\n", dircolor)

            # and its contents
            for f in files
                print(join(fill(" " ^ 4,  depth)))
                colorprint(string("$CORNER", join(fill("$HORIZONTAL" ^ 3, 1))), linecolor)
                print(f)
                if ages
                    status = stat(joinpath(path, root, f))
                    fileage = Dates.now() - Dates.unix2datetime(status.mtime)
                    daysold = round(convert(Float64, fileage.value/(1000 * 60 * 60 * 24)), 2)
                    print("\t($(daysold) days old) ")
                end
                println()
            end
        end
    end
end

end # module

ext: jl


text:
module FrequencyTrajectories

using Parameters
using RecipesBase

import Base: iterate, eltype, length
import Base: in, isless, filter!

export FT
export FrequencyBin, Trajectory
export filter!, fixation_probability

include("objects.jl")
include("interfaces.jl")
include("methods.jl")
include("filtering.jl")

const FT = FrequencyTrajectories

end # module


text:
# ---
# title: 1080. Insufficient Nodes in Root to Leaf Paths
# id: problem1080
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search
# link: <https://leetcode.com/problems/insufficient-nodes-in-root-to-leaf-paths/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, consider all _root to leaf paths_ : paths
# from the root to any leaf.  (A leaf is a node with no children.)
#
# A `node` is _insufficient_ if  **every** such root to leaf path intersecting
# this `node` has sum strictly less than `limit`.
#
# Delete all insufficient nodes simultaneously, and return the root of the
# resulting binary tree.
#
#
#
# **Example 1:**
#
#
#
#     **![](https://assets.leetcode.com/uploads/2019/06/05/insufficient-11.png)Input:** root = [1,2,3,4,-99,-99,7,8,9,-99,-99,12,13,-99,14], limit = 1
#     **![](https://assets.leetcode.com/uploads/2019/06/05/insufficient-2.png)Output:** [1,2,3,4,null,null,7,8,9,null,14]
#
#
# **Example 2:**
#
#
#
#     **![](https://assets.leetcode.com/uploads/2019/06/05/insufficient-3.png)Input:** root = [5,4,8,11,null,17,4,7,1,null,null,5,3], limit = 22
#     **![](https://assets.leetcode.com/uploads/2019/06/05/insufficient-4.png)Output:** [5,4,8,11,null,17,4,7,null,null,null,5]
#
#
#
# **Example 3:**
#
#
#
#     **![](https://assets.leetcode.com/uploads/2019/06/11/screen-shot-2019-06-11-at-83301-pm.png)Input:** root = [1,2,-3,-5,null,4,null], limit = -1
#     ![](https://assets.leetcode.com/uploads/2019/06/11/screen-shot-2019-06-11-at-83517-pm.png)Output: [1,null,-3,4]
#
#
#
# **Note:**
#
#   1. The given tree will have between `1` and `5000` nodes.
#   2. `-10^5 <= node.val <= 10^5`
#   3. `-10^9 <= limit <= 10^9`
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    MapManager(params::Params, frame::Frame, extractor::Extractor)

Map Manager is responsible for managing Keyframes in the map
as well as Mappoints.

# Arguments

- `current_frame::Frame`: Current frame that is shared throughout
    all the components in the system.
- `frames_map::Dict{Int64, Frame}`: Map of the Keyframes (its id → Keyframe).
- `params::Params`: Parameters of the system.
- `extractor::Extractor`: Extractor for finding keypoints in the frames.
- `map_points::Dict{Int64, MapPoint}`: Map of all the map_points
    (its id → MapPoints).
- `current_mappoint_id::Int64`: Id of the current Mappoint to be created.
    It is incremented each time a new Mappoint is added to `map_points`.
- `current_keyframe_id::Int64`: Id of the current Keyframe to be created.
    It is incremented each time a new Keyframe is added to `frames_map`.
- `nb_keyframes::Int64`: Total number of keyframes.
- `nb_mappoints::Int64`: Total number of mappoints.
"""
mutable struct MapManager
    current_frame::Frame
    frames_map::Dict{Int64, Frame}

    params::Params
    extractor::Extractor
    map_points::Dict{Int64, MapPoint}

    current_mappoint_id::Int64
    current_keyframe_id::Int64
    nb_keyframes::Int64
    nb_mappoints::Int64

    mappoint_lock::ReentrantLock
    keyframe_lock::ReentrantLock
    map_lock::ReentrantLock
    optimization_lock::ReentrantLock
end

function MapManager(params::Params, frame::Frame, extractor::Extractor)
    mappoint_lock = ReentrantLock()
    keyframe_lock = ReentrantLock()
    map_lock = ReentrantLock()
    optimization_lock = ReentrantLock()

    MapManager(
        frame, Dict{Int64, Frame}(), params, extractor,
        Dict{Int64, MapPoint}(), 0, 0, 0, 0,
        mappoint_lock, keyframe_lock, map_lock, optimization_lock)
end

function get_keyframe(m::MapManager, kfid)
    lock(m.keyframe_lock) do
        get(m.frames_map, kfid, nothing)
    end
end

function has_keyframe(m::MapManager, kfid)
    lock(m.keyframe_lock) do
        kfid in keys(m.frames_map)
    end
end

function get_mappoint(m::MapManager, mpid)
    lock(m.mappoint_lock) do
        get(m.map_points, mpid, nothing)
    end
end

function create_keyframe!(m::MapManager, image)
    @debug "[MM] Creating new keyframe $(m.current_keyframe_id)."
    prepare_frame!(m)
    extract_keypoints!(m, image)
    add_keyframe!(m)
end

function prepare_frame!(m::MapManager)
    m.current_frame.kfid = m.current_keyframe_id
    @debug "[MM] Adding KF $(m.current_frame.kfid) to Map."

    # Filter if there are too many keypoints.
    # if m.current_frame.nb_keypoints > m.params.max_nb_keypoints
    #     # TODO
    # end

    for kp in values(m.current_frame.keypoints)
        mp = get(m.map_points, kp.id, nothing)
        if mp ≡ nothing
            remove_obs_from_current_frame!(m, kp.id)
        else
            add_keyframe_observation!(mp, m.current_keyframe_id)
        end
    end
end

function extract_keypoints!(m::MapManager, image)
    nb_2_detect = m.params.max_nb_keypoints - m.current_frame.nb_occupied_cells
    nb_2_detect ≤ 0 && return
    # Detect keypoints in the provided `image`
    # using current keypoints to set a mask of regions
    # to avoid detecting features in.
    current_points = [kp.pixel for kp in values(m.current_frame.keypoints)]
    keypoints = detect(m.extractor, image, current_points)
    isempty(keypoints) && return

    if m.params.do_local_matching
        descriptors, keypoints = describe(m.extractor, image, keypoints)
    else
        descriptors = fill(BitVector(), length(keypoints))
    end
    add_keypoints_to_frame!(m, m.current_frame, keypoints, descriptors)
end

function add_keypoints_to_frame!(m::MapManager, frame, keypoints, descriptors)
    lock(m.mappoint_lock) do
        for (kp, dp) in zip(keypoints, descriptors)
            # m.current_mappoint_id is incremented in `add_mappoint!`.
            add_keypoint!(frame, Point2f(kp[1], kp[2]), m.current_mappoint_id)
            add_mappoint!(m, dp)
        end
    end
end

@inline function add_mappoint!(m::MapManager, descriptor)
    mp = MapPoint(m.current_mappoint_id, m.current_keyframe_id, descriptor)
    m.map_points[m.current_mappoint_id] = mp
    m.current_mappoint_id += 1
    m.nb_mappoints += 1
end

"""
Remove MapPoint from the map given its id.

Removing a mappoint, also update covisibility scores of the observer Frames.
If MapPoint is observed by the current Frame, remove its keypoint as well.
"""
function remove_mappoint!(m::MapManager, id)
    lock(m.keyframe_lock)
    lock(m.mappoint_lock)

    if !(id in keys(m.map_points))
        unlock(m.mappoint_lock)
        unlock(m.keyframe_lock)
        return
    end

    mp = m.map_points[id]
    for observer_id in mp.observer_keyframes_ids
        observer_kf = get(m.frames_map, observer_id, nothing)
        observer_kf ≡ nothing && continue

        remove_keypoint!(observer_kf, id)
        for co_observer_id in mp.observer_keyframes_ids
            observer_id == co_observer_id && continue
            decrease_covisible_kf!(observer_kf, co_observer_id)
        end
    end

    mp.is_observed && remove_keypoint!(m.current_frame, id)
    mp.is_3d && (m.nb_mappoints -= 1;)
    pop!(m.map_points, id)

    unlock(m.keyframe_lock)
    unlock(m.mappoint_lock)
end

"""
Copy current MapManager's Frame and add it to the KeyFrame map.
Increase current keyframe id & total number of keyframes.
"""
function add_keyframe!(m::MapManager)
    new_keyframe = deepcopy(m.current_frame)

    lock(m.keyframe_lock) do
        m.frames_map[m.current_keyframe_id] = new_keyframe
        m.current_keyframe_id += 1
        m.nb_keyframes += 1
    end
    @debug "[MM] Added."
end

function remove_keyframe!(m::MapManager, kfid)
    lock(m.keyframe_lock)
    lock(m.mappoint_lock)

    kf = get(m.frames_map, kfid, nothing)
    if kf ≡ nothing
        unlock(m.mappoint_lock)
        unlock(m.keyframe_lock)
        return
    end

    for kp in get_keypoints(kf)
        mp = get(m.map_points, kp.id, nothing)
        mp ≢ nothing && remove_kf_observation!(mp, kfid)
    end
    for cov_kfid in keys(kf.covisible_kf)
        cov_kf = get(m.frames_map, cov_kfid, nothing)
        cov_kf ≢ nothing && remove_covisible_kf!(cov_kf, kfid)
    end

    pop!(m.frames_map, kfid)
    m.nb_keyframes -= 1

    unlock(m.mappoint_lock)
    unlock(m.keyframe_lock)
end

"""
Remove a MapPoint observation from current Frame by `id`.
"""
function remove_obs_from_current_frame!(m::MapManager, id::Int64)
    remove_keypoint!(m.current_frame, id)
    # Set MapPoint as not observable.
    mp = get(m.map_points, id, nothing)
    mp ≢ nothing && (mp.is_observed = false;)
end

"""
Remove KeyFrame observation from MapPoint.
"""
function remove_mappoint_obs!(m::MapManager, kpid::Int, kfid::Int)
    lock(m.keyframe_lock)
    lock(m.mappoint_lock)

    # Remove MapPoint from KeyFrame.
    kf = get(m.frames_map, kfid, nothing)
    kf ≢ nothing && remove_keypoint!(kf, kpid)

    mp = get(m.map_points, kpid, nothing)
    if mp ≡ nothing
        unlock(m.mappoint_lock)
        unlock(m.keyframe_lock)
        return
    end

    # Remove KeyFrame observation from MapPoint.
    remove_kf_observation!(mp, kfid)

    if kf ≢ nothing
        for observer_id in mp.observer_keyframes_ids
            observer_kf = get(m.frames_map, observer_id, nothing)
            observer_kf ≡ nothing && continue

            decrease_covisible_kf!(kf, observer_id)
            decrease_covisible_kf!(observer_kf, kfid)
        end
    end

    unlock(m.keyframe_lock)
    unlock(m.mappoint_lock)
end

"""
    update_mappoint!(m::MapManager, mpid, new_position)

Update position of a MapPoint.
"""
function update_mappoint!(m::MapManager, mpid, new_position)
    lock(m.keyframe_lock)
    lock(m.mappoint_lock)

    if !(mpid in keys(m.map_points))
        unlock(m.mappoint_lock)
        unlock(m.keyframe_lock)
        return
    end

    mp = m.map_points[mpid]
    # If MapPoint is 2D, turn it to 3D and update its observing KeyFrames.
    if !mp.is_3d
        for observer_id in mp.observer_keyframes_ids
            if observer_id in keys(m.frames_map)
                turn_keypoint_3d!(m.frames_map[observer_id], mpid)
            else
                remove_kf_observation!(mp, observer_id)
            end
        end
        if mp.is_observed
            # Because we deepcopy frame before putting it to the frames_map,
            # we need to update current frame as well.
            # Which should also update current frame in the FrontEnd.
            turn_keypoint_3d!(m.current_frame, mpid)
        end
    end
    set_position!(mp, new_position)

    unlock(m.mappoint_lock)
    unlock(m.keyframe_lock)
end

"""
    update_frame_covisibility!(map_manager::MapManager, frame::Frame)

Update covisibility graph for the `frame`.
This is done by going through all of the keypoints in the `frame`.
Getting their corresponding mappoints. And joining sets of observers for
those mappoints.
"""
function update_frame_covisibility!(map_manager::MapManager, frame::Frame)
    covisible_keyframes = Dict{Int64, Int64}()
    local_map_ids = Set{Int64}()
    # For each Keypoint in the `frame`, get its corresponding MapPoint.
    # Get the set of KeyFrames that observe this MapPoint.
    # Add them to the covisible map, which contains all KeyFrames
    # that share visibility with the `frame`.
    for kp in get_keypoints(frame)
        if !(kp.id in keys(map_manager.map_points))
            remove_mappoint_obs!(map_manager, kp.id, frame.kfid)
            remove_obs_from_current_frame!(map_manager, kp.id)
            continue
        end
        mp = get_mappoint(map_manager, kp.id)
        mp_observers = get_observers(mp)
        # Get the set of KeyFrames observing this KeyFrame to update covisibility.
        for kfid in mp_observers
            kfid == frame.kfid && continue
            if kfid in keys(covisible_keyframes)
                covisible_keyframes[kfid] += 1
            else
                covisible_keyframes[kfid] = 1
            end
        end
    end
    # Update covisibility for covisible KeyFrames.
    # For each entry in the covisible map, get its corresponding KeyFrame.
    # Update the covisible score for the `frame` in it.
    # Add all 3D Keypoints that are not in the `frame`
    # to the local map for future tracking.
    bad_kfids = Set{Int64}()
    for (kfid, cov_score) in covisible_keyframes
        if !(kfid in keys(map_manager.frames_map))
            push!(bad_kfids, kfid)
            continue
        end
        cov_frame = get_keyframe(map_manager, kfid)
        add_covisibility!(cov_frame, frame.kfid, cov_score)
        for kp in get_3d_keypoints(cov_frame)
            kp.id in keys(frame.keypoints) || push!(local_map_ids, kp.id)
        end
    end
    for bad_kfid in bad_kfids
        pop!(covisible_keyframes, bad_kfid)
    end
    # Update the set of covisible KeyFrames.
    set_covisible_map!(frame, covisible_keyframes)
    # Update local map of unobserved MapPoints.
    if length(local_map_ids) > 0.5 * length(frame.local_map_ids)
        frame.local_map_ids = local_map_ids
    else
        union!(frame.local_map_ids, local_map_ids)
    end
end

"""
    reset!(m::MapManager)

Reset map manager.
"""
function reset!(m::MapManager)
    m.nb_keyframes = 0
    m.nb_mappoints = 0
    m.current_keyframe_id = 0
    m.current_mappoint_id = 0

    m.map_points |> empty!
    m.frames_map |> empty!
end

"""
    merge_mappoints(m::MapManager, prev_id, new_id)

Merge `prev_id` Mappoint into `new_id` Mappoint.
For "previous" observers, update mappoint and keypoint to the new.
"""
function merge_mappoints(m::MapManager, prev_id, new_id)
    lock(m.mappoint_lock)
    lock(m.keyframe_lock)
    try
        prev_mp = get(m.map_points, prev_id, nothing)
        prev_mp ≡ nothing && return
        new_mp = get(m.map_points, new_id, nothing)
        new_mp ≡ nothing && return
        new_mp.is_3d || return

        prev_observers = get_observers(prev_mp)
        new_observers = get_observers(new_mp)

        # For previous mappoint observers, update keypoint for them.
        # If successfull, then add covisibility link between old and new
        # observer keyframes.
        for prev_observer_id in prev_observers
            prev_observer_kf = get(m.frames_map, prev_observer_id, nothing)
            prev_observer_kf ≡ nothing && continue
            update_keypoint!(
                prev_observer_kf, prev_id, new_id, new_mp.is_3d) || continue

            add_keyframe_observation!(new_mp, prev_observer_id)
            for new_observer_id in new_observers
                new_observer_kf = get(m.frames_map, new_observer_id, nothing)
                new_observer_kf ≡ nothing && continue

                add_covisibility!(new_observer_kf, prev_observer_id)
                add_covisibility!(prev_observer_kf, new_observer_id)
            end
        end

        for (kfid, descriptor) in prev_mp.keyframes_descriptors
            add_descriptor!(new_mp, kfid, descriptor)
        end
        if is_observing_kp(m.current_frame, prev_id)
            update_keypoint!(m.current_frame, prev_id, new_id, new_mp.is_3d)
        end

        # Update nb mappoints and erase old mappoint.
        prev_mp.is_3d && (m.nb_mappoints -= 1;)
        pop!(m.map_points, prev_id)
    catch e
        showerror(stdout, e); println()
        display(stacktrace(catch_backtrace())); println()
    finally
        unlock(m.keyframe_lock)
        unlock(m.mappoint_lock)
    end
end

"""
    optical_flow_matching!(map_manager, frame, from_pyramid, to_pyramid, stereo)

Match keypoints in `frame` from `from_pyramid` to `to_pyramid`.
This function is used when matching keypoints temporally from previous frame
to current frame, or when matching keypoints between stereo image.

If there are 3D keypoints in the `frame`, then try to match respectful
keypoints using displacement guess from motion model (when matching temporally)
or calibration pose (in stereo).

# Arguments

- `map_manager::MapManager`: Map manager, used for retrieving parameters info,
    3D mappoints, removing mappoints.
- `frame`: Frame for which to do matching.
- `from_pyramid::LKPyramid`: Pyramid from which to track.
- `to_pyramid::LKPyramid`: Pyramid to which to track.
- `stereo::Bool`: Set to `true` if doing stereo matching. It will retrieve data
    using mutex and update `stereo` keypoints instead of regular keypoints.
    Otherwise set to `false`.
"""
function optical_flow_matching!(
    map_manager::MapManager, frame, from_pyramid, to_pyramid, stereo,
)
    window_size = map_manager.params.window_size
    max_distance = map_manager.params.max_ktl_distance
    pyramid_levels = map_manager.params.pyramid_levels

    pyramid_levels_3d = 1
    ids = Vector{Int64}(undef, frame.nb_keypoints)
    pixels = Vector{Point2f}(undef, frame.nb_keypoints)

    ids3d = Vector{Int64}(undef, frame.nb_3d_kpts)
    pixels3d = Vector{Point2f}(undef, frame.nb_3d_kpts)
    displacements3d = Vector{Point2f}(undef, frame.nb_3d_kpts)

    scale = 1.0 / 2.0^pyramid_levels_3d
    i, i3d = 1, 1
    n_good = 0

    keypoints = stereo ? get_keypoints(frame) : values(frame.keypoints)
    for kp in keypoints
        if !kp.is_3d
            pixels[i] = kp.pixel
            ids[i] = kp.id
            i += 1
            continue
        end

        mp = stereo ? get_mappoint(map_manager, kp.id) : map_manager.map_points[kp.id]
        if mp ≡ nothing
            remove_mappoint_obs!(map_manager, kp.id, frame.kfid)
            continue
        end

        position = get_position(mp)
        projection = stereo ?
            project_world_to_right_image_distort(frame, position) :
            project_world_to_image_distort(frame, position)

        if stereo
            if in_right_image(frame, projection)
                ids3d[i3d] = kp.id
                pixels3d[i3d] = kp.pixel
                displacements3d[i3d] = scale .* (projection .- kp.pixel)
                i3d += 1
            else
                remove_mappoint_obs!(map_manager, kp.id, frame.kfid)
                continue
            end
        else
            if in_image(frame, projection)
                ids3d[i3d] = kp.id
                pixels3d[i3d] = kp.pixel
                displacements3d[i3d] = scale .* (projection .- kp.pixel)
                i3d += 1
            end
        end
    end

    i3d -= 1
    ids3d = @view(ids3d[1:i3d])
    pixels3d = @view(pixels3d[1:i3d])
    displacements3d = @view(displacements3d[1:i3d])

    failed_3d = true
    if !isempty(ids3d)
        new_keypoints, status = fb_tracking!(
            from_pyramid, to_pyramid, pixels3d;
            displacement=displacements3d,
            pyramid_levels=pyramid_levels_3d,
            window_size, max_distance)

        nb_good = 0
        for j in 1:length(status)
            if status[j]
                if stereo
                    succ = maybe_stereo_update!(frame, ids3d[j], new_keypoints[j])
                    succ && (n_good += 1;)
                else
                    update_keypoint!(frame, ids3d[j], new_keypoints[j])
                    nb_good += 1
                end
            else
                # If failed → add to track with 2d keypoints w/o prior.
                pixels[i] = pixels3d[j]
                ids[i] = ids3d[j]
                i += 1
            end
        end
        @debug "[MM] 3D Points tracked $nb_good. Stereo $stereo."
        failed_3d = nb_good < 0.33 * length(ids3d)
    end

    i -= 1
    pixels = @view(pixels[1:i])
    ids = @view(ids[1:i])

    isempty(pixels) && return nothing
    new_keypoints, status = fb_tracking!(
        from_pyramid, to_pyramid, pixels;
        pyramid_levels, window_size, max_distance)

    for j in 1:length(new_keypoints)
        if stereo
            status[j] && maybe_stereo_update!(frame, ids[j], new_keypoints[j]) &&
                (n_good += 1;)
        else
            status[j] ?
                update_keypoint!(frame, ids[j], new_keypoints[j]) :
                remove_obs_from_current_frame!(map_manager, ids[j])
        end
    end
    nothing
end

"""
    maybe_stereo_update!(
        frame::Frame, kpid, new_position::Point2f; epipolar_error::Float64 = 2.0,
    )

Update stereo keypoint if the vertical distance between matched keypoints is
less than `epipolar_error`. In this case, set y-coordinate equal
to the left keypoint. Otherwise do nothing.

# Returns

`true` if successfully updated, otherwise `false`.
"""
function maybe_stereo_update!(
    frame::Frame, kpid, new_position; epipolar_error = 2.0,
)
    kp = get_keypoint(frame, kpid)
    right_pixel = undistort_point(frame.right_camera, new_position)
    abs(kp.undistorted_pixel[1] - right_pixel[1]) > epipolar_error &&
        return false
    # Make y-coordinate the same.
    corrected = Point2f(kp.pixel[1], new_position[2])
    update_stereo_keypoint!(frame, kpid, corrected)
    true
end


text:
using WriteVTK
using Printf

let

    dir = string("../../erickson/output_files/", ARGS[1])
    v_file = open(string(dir, "/vs.dat"))

    @printf "Opened u and v files...\n"

    vs = eachline(v_file)

    @printf "Read v data\n..."

    x = map(x->parse(Float64, x), split.(Iterators.take(vs, 1))[1])
    y = map(x->parse(Float64, x), split.(Iterators.take(vs, 2))[1])

    count = 1
    @printf "Converting to vtk...\n"
    for line in vs
        @printf "\r%d" count
        vtk = vtk_grid(string(dir, "/vtk_filesV/cycles_", count, ".vtr"), x, y)
        v = map(x-> parse(Float64, x), split.(line))
        vtk["v"] = v
        saved_files = vtk_save(vtk)
        count += 1
    end

    close(v_file)

    @printf ("Done.\n")

end


text:
@testset "Line search errors" begin
    hz = LineSearches.HagerZhang(; delta = 0.1, sigma = 0.9, alphamax = Inf,
                      rho = 5.0, epsilon = 1e-6, gamma = 0.66, linesearchmax = 2)
    for optimizer in (ConjugateGradient, GradientDescent, LBFGS, BFGS, Newton, AcceleratedGradientDescent, MomentumGradientDescent)
        debug_printing && println("Testing $(string(optimizer))")
        prob = MultivariateProblems.UnconstrainedProblems.examples["Exponential"]
        @test optimize(MVP.objective(prob), prob.initial_x, optimizer(alphaguess = LineSearches.InitialPrevious(), linesearch = hz)).ls_success == false
    end
end


text:
"ACSR data from spec table; english units."
struct ACSRSpecsEnglish
    "Catchy name"
    label::String
    "Number of conductors per phase"
    bundle::Int64
    "Conductor diameter [in]"
    D::Float64
    "Aluminum weight [lb/1000ft]"
    Al_m::Float64
    "Steel weight [lb/1000ft]"
    St_m::Float64
    "Resistance assuming AC current and 75 C temperature [ohm/1000ft]"
    R::Float64
end

"ACSR data from spec table; SI units."
struct ACSRSpecsMetric
    "Catchy name"
    label::String
    "Number of conductors per phase"
    bundle::Int64
    "Conductor diameter [m]"
    D::Float64
    "Aluminum weight [kg/m]"
    Al_m::Float64
    "Steel weight [kg/m]"
    St_m::Float64
    "Resistance assuming AC current and 75 C temperature [ohm/m]"
    R::Float64
end

"""
In: nothing

Out:

**ACSR table** with the following values for each conductor type:

* `D` [in],             Conductor diameter
* `Al_m` [lb/1000ft],   Aluminum weight
* `St_m` [lb/1000ft],   Steel weight
* `R` [ohm/1000ft],     Resistance assuming AC current and 75 C temperature
* `bundle` [-],         Number of conductors per phase

**`labels`** [string], Label for each conductor type

Credit to Mads Almassalkhi for compiling data.
"""
function acsr_table_english()
    acsr = [
    0.198   24  12  0.806  105 #'Turkey'   ,   '6/1'
    0.250   39  18  0.515  140 #'Swan'     ,   '6/1'
    0.316   62  29  0.332  184 #'Sparrow'  ,   '6/1'
    0.354   78  37  0.268  212 #'Robin'    ,   '6/1'
    0.398   99  47  0.217  242 #'Raven'    ,   '6/1'
    0.447  124  59  0.176  276 #'Quail'    ,   '6/1'
    0.502  156  74  0.144  315 #'Pigeon'   ,   '6/1'
    0.563  197  93  0.119  357 #'Penguin'  ,   '6/1'
    0.609  250  39 0.0787  449 #'Waxwing'  ,  '18/1'
    0.642  251 115 0.0779  475 #'Partridge',  '26/7'
    0.680  283 130 0.0693  492 #'Ostrich'  ,  '26/7'
    0.684  315  49 0.0625  519 #'Merlin'   ,  '18/1'
    0.720  317 146 0.0618  529 #'Linnet'   ,  '26/7'
    0.741  318 209 0.0613  535 #'Oriole'   ,  '30/7'
    0.743  373  58 0.0529  576 #'Chickadee',  '18/1'
    0.772  374 137 0.0526  584 #'Brant'    ,  '24/7'
    0.783  374 172 0.0523  587 #'Ibis'     ,  '26/7'
    0.806  375 247 0.0519  594 #'Lark'     ,  '30/7'
    0.814  447  70 0.0442  646 #'Pelican'  ,  '18/1'
    0.846  449 164 0.0439  655 #'Flicker'  ,  '24/7'
    0.883  450 296 0.0433  666 #'Hen'      ,  '30/7'
    0.879  522  82 0.0379  711 #'Osprey'   ,  '18/1'
    0.914  524 192 0.0376  721 #'Parakeet' ,  '24/7'
    0.953  525 345 0.0372  734 #'Eagle'    ,  '30/7'
    0.953  570 209 0.0346  760 #'Peacock'  ,  '24/7'
    0.994  571 375 0.0342  774 #'Wood Duck',  '30/7'
    0.977  599 219 0.0330  784 #'Rook'     ,  '24/7'
    1.019  600 395 0.0325  798 #'Scoter'   ,  '30/7'
    1.014  628 289 0.0313  812 #'Gannet'   ,  '26/7'
    1.051  674 310 0.0292  849 #'Starling' ,  '26/7'
    1.081  676 435 0.0290  859 #'Redwing'  , '30/19'
    1.040  745  58 0.0268  884 #'Coot'     ,  '36/1'
    1.107  749 344 0.0263  907 #'Drake'    ,  '26/7'
    1.140  751 483 0.0261  918 #'Mallard'  , '30/19'
    1.162  848 310 0.0241  961 #'Canary'   ,  '54/7'
    1.196  899 329 0.0228  996 #'Cardinal' ,  '54/7'
    1.245  973 356 0.0211 1047 #'Curlew'   ,  '54/7'
    1.292 1053 375 0.0197 1093 #'Finch'    , '54/19'
    1.302 1123 219 0.0182 1139 #'Bunting'  ,  '45/7'
    1.345 1198 234 0.0171 1184 #'Bittern'  ,  '45/7'
    1.386 1273 248 0.0162 1229 #'Dipper'   ,  '45/7'
    1.427 1348 263 0.0153 1272 #'Bobolink' ,  '45/7'
    1.504 1498 292 0.0139 1354 #'Lapwing'  ,  '45/7'
    1.602 1685 386 0.0125 1453 #'Chukar'   , '84/19'
    1.735 2051 249 0.0106 1607 #'Kiwi'     ,  '72/7'
    1.762 2040 468 0.0105 1623 #'Bluebird' , '84/19'
        ]
    labels = [
        "Turkey"; "Swan"; "Sparrow"; "Robin"; "Raven";
        "Quail"; "Pigeon"; "Penguin"; "Waxwing"; "Partridge";
        "Ostrich"; "Merlin"; "Linnet"; "Oriole"; "Chickadee";
        "Brant"; "Ibis"; "Lark"; "Pelican"; "Flicker"; "Hen";
        "Osprey"; "Parakeet"; "Eagle"; "Peacock"; "Wood Duck";
        "Rook"; "Scoter"; "Gannet"; "Starling"; "Redwing";
        "Coot"; "Drake"; "Mallard"; "Canary"; "Cardinal";
        "Curlew"; "Finch"; "Bunting"; "Bittern"; "Dipper";
        "Bobolink"; "Lapwing"; "Chukar"; "Kiwi"; "Bluebird"
    ]
    return acsr, labels
end

"""
In: nothing

Out:

**ACSR table** with the following values for each conductor type:

    * `D` [m],         Conductor diameter
    * `Al_m` [kg/m],   Aluminum weight
    * `St_m` [kg/m],   Steel weight
    * `R` [ohm/m],     Resistance assuming AC current and 75 C temperature
    * `bundle` [-],    Number of conductors per phase

**`labels`** [string],   Label for each conductor type

Credit to Mads Almassalkhi for compiling data.
"""
function acsr_table_si()
    acsr, labels = acsr_table_english()
    # Convert diameter: 1 in = 25.4e-3 m
    acsr[:, 1] = acsr[:, 1] * 25.4e-3

    # Convert weight: 1 lb = 0.453592 kg and 1000 ft = 304.8 m
    acsr[:, 2:3] = acsr[:, 2:3] * 0.453592 / 304.8

    # Convert resistance: 1000ft = 304.8 m
    acsr[:, 4] = acsr[:, 4] / 304.8
    # Scale resistances to better match RTS96 system data.
    # These values are tuned so that estimateLineLength.m will
    # give a result very close to the lengths provided for the RTS96 system.
    resistanceFactor = 2.0
    acsr[:, 4] = acsr[:, 4] / resistanceFactor
    return acsr, labels
end

"""
    acsr_specs = acsr_interpolation(I_lim, V_base)

Estimate bundling configuration and ACSR data based on I_lim and V_base.

In:

* `I_lim` [A],          Line current limit
* `V_base` [V],         Line base voltage
* `metric` (Bool),      Whether to use metric (SI) or English units

Out:

* `D` [m],              Conductor diameter
* `Al_m` [kg/m],        Aluminum weight
* `St_m` [kg/m],        Steel weight
* `R` [ohm/m],          Resistance assuming AC current and 75 C temperature
* `bundle` [-],         Number of conductors per phase
* `label` [string],     Label for each conductor type

Credit to Mads Almassalkhi for original MATLAB implementation. See:

*  Bergen and Vittal "Power Systems Analysis" 2nd ed 2000 p.85
*  Weedy and Cory "Electric Power Systems" 4th ed 1998 p.129
"""
function acsr_interpolation(I_lim::Float64, V_base::Float64; metric::Bool=true)
    acsr, labels = metric ? acsr_table_si() : acsr_table_english()

    maxcurrent = maximum(acsr[:, 5])
    largebundle = 1000

    if V_base <= 138e3
        bundle = (I_lim <= 1000) ? 1 :
            (I_lim <= maxcurrent * 2) ? 2 : largebundle
    elseif V_base <= 220e3
        bundle = (I_lim <= 1000) ? 1 :
            (I_lim <= maxcurrent * 2) ? 2 :
                (I_lim <= maxcurrent * 3) ? 3 : largebundle
    elseif V_base <= 345e3
        bundle = (I_lim <= maxcurrent) ? 1 :
            (I_lim <= maxcurrent * 2) ? 2 :
                (I_lim <= maxcurrent * 3) ? 3 : largebundle
    elseif V_base <= 500e3
        bundle = (I_lim <= maxcurrent * 2) ? 2 :
            (I_lim <= maxcurrent * 4) ? 4 :
                (I_lim <= maxcurrent * 6) ? 6 : largebundle
    else
        bundle = (I_lim <= maxcurrent * 3) ? 3 :
            (I_lim <= maxcurrent * 4) ? 4 :
                (I_lim <= maxcurrent * 6) ? 6 : largebundle
    end

    # Search for ACSR match
    row = 1
    found = false
    while !found && row <= size(acsr, 1)
        if I_lim <= acsr[row, 5] * bundle
            found = true
            row -= 1
            row == 0 && (row = 1)
        else
            row += 1
        end
    end

    # Assign output variables
    if !found
        @warn "No ACSR match for line with I_lim = $I_lim A and V_base = $V_base kV"
        D, Al_m, St_m, R = NaN, NaN, NaN, NaN
        label=""
    else
        D, Al_m, St_m, R = acsr[row, :]
        label = labels[row]
    end
    return metric ? ACSRSpecsMetric(label, bundle, D, Al_m, St_m, R) :
                    ACSRSpecsEnglish(label, bundle, D, Al_m, St_m, R)
end

"""
Estimate ACSR specs for all branches, given a `network_data` dict. The input
dict is the output of `PowerModels.parse_file`.
"""
function acsr_interpolation(network_data::Dict{String, Any}; metric::Bool=true)
    !network_data["per_unit"] && @error "Network data is not in per unit!"

    # even with network_data["per_unit"] == true, basekv values are still in kv
    basekv = branch_basekv(network_data)

    baseMVA = network_data["baseMVA"]
    branch_data = Vector{Tuple}()
    for (k, v) in network_data["branch"]
        # note: if network_data["per_unit"] == true, then rate_a has been
        # divided by baseMVA.
        push!(branch_data, (k, v["rate_a"] * baseMVA))
    end

    T = metric ? ACSRSpecsMetric : ACSRSpecsEnglish
    acsr_specs = Dict{String, T}()

    for i in 1:length(branch_data)
        id, rate_a = branch_data[i]
        V_base = basekv[string(id)] * 1e3 # convert kV to V

        # convert long-term VA rating to current rating
        I_lim = rate_a * 1e6 / V_base

        acsr_specs[id] = acsr_interpolation(I_lim, V_base; metric=metric)
    end
    return acsr_specs
end

"""
Return current limits for all branches, given a `network_data`
dict. The input is the output of `PowerModels.parse_file`.
"""
function return_current_limits(network_data::Dict{String, Any})
    !network_data["per_unit"] && @error "Network data is not in per unit!"

    # even with network_data["per_unit"] == true, basekv values are still in kv
    basekv = branch_basekv(network_data)
    baseMVA = network_data["baseMVA"]

    branch_data = Vector{Tuple}()
    for (k, v) in network_data["branch"]
        # note: if network_data["per_unit"] == true, then rate_a has been
        # divided by baseMVA.
        push!(branch_data, (k, v["rate_a"] * baseMVA))
    end

    current_limits = Dict{String, Float64}()

    for i in 1:length(branch_data)
        id, rate_a = branch_data[i]
        V_base = basekv[string(id)] * 1e3 # convert kV to V

        # convert long-term VA rating to current rating
        I_lim = rate_a * 1e6 / V_base

        current_limits[id] = I_lim
    end
    return current_limits
end


text:
using BVHFiles
using Test

@testset "interpolation" begin
    println("Starting test of 'interpolation'")

    g = load("Example.bvh") |>
        interpolate!(1)

    @test frames(g) == 481

    g = load("Example.bvh") |>
    interpolate!(2)

    @test frames(g) == 721

    g = load("Example.bvh") |>
    interpolate!(3)

    @test frames(g) == 961
end


text: import QMTK.Consts


text:
module EMLbatch_dist
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("batch_dist/Diff_Dist.jl")
end


text:
function test_measures()
    m1 = Lebesgue()
    @test !isdiscrete(m1)
    @test iscontinuous(m1)
    @test !isnormalized(m1)
    @test domaintype(m1) == Float64
    @test DomainIntegrals.unsafe_weightfun(m1, 0.4) == 1

    m2 = DomainIntegrals.LebesgueUnit()
    @test !isdiscrete(m2)
    @test iscontinuous(m2)
    @test isnormalized(m2)
    @test domaintype(m2) == Float64
    @test support(m2) == UnitInterval()
    @test weightfun(m2, 0.4) == 1
    @test weightfun(m2, 1.4) == 0

    m3 = LegendreWeight()
    @test !isdiscrete(m3)
    @test iscontinuous(m3)
    @test !isnormalized(m3)
    @test domaintype(m3) == Float64
    @test support(m3) == ChebyshevInterval()
    @test weightfun(m3, 0.4) == 1
    @test weightfun(m3, -0.4) == 1
    @test weightfun(m3, 1.4) == 0
    @test weightfun(m3, -1.4) == 0

    x = 0.5
    m4 = DiracWeight(x)
    @test !isdiscrete(m4)
    @test iscontinuous(m4)
    @test isnormalized(m4)
    @test domaintype(m4) == Float64
    @test support(m4) == Point(x)
    @test point(m4) == x
    @test weightfun(m4, x) == Inf
    @test weightfun(m4, x+1) == 0
    @test weightfun(m4, big(x)) == Inf
    @test weightfun(m4, big(x+1)) == 0

    m5 = GaussianWeight{SVector{2,Float64}}()
    @test !isdiscrete(m5)
    @test iscontinuous(m5)
    @test isnormalized(m5)
    @test domaintype(m5) == SVector{2,Float64}
    @test support(m5) == FullSpace{SVector{2,Float64}}()
    @test weightfun(m5, SVector(0.0,0.0)) ≈ 1/(2pi)
    @test weightfun(m5, SVector(0,0)) ≈ 1/(2pi)
    @test weightfun(m5, SVector(big(0),big(0))) ≈ 1/(2pi)
    @test weightfun(m5, SVector(big(0.0),big(0.0))) ≈ 1/(2pi)

    m6 = LaguerreWeight(0.0)
    @test !isdiscrete(m6)
    @test iscontinuous(m6)
    @test isnormalized(m6)
    @test !isnormalized(LaguerreWeight(0.1))
    @test domaintype(m6) == Float64
    @test support(m6) == HalfLine()
    @test weightfun(m6, 0.4) == exp(-0.4)
    @test weightfun(m6, -0.4) == 0
    @test weightfun(m6, big(0.4)) == exp(-big(0.4))
end


text:
using JuliaInterpreter
using Test, Random, InteractiveUtils, Distributed, Dates

# Much of this file is taken from Julia's test/runtests.jl file.

if !isdefined(Main, :read_and_parse)
    include("utils.jl")
end

const juliadir = dirname(dirname(Sys.BINDIR))
const testdir = joinpath(juliadir, "test")
if isdir(testdir)
    include(joinpath(testdir, "choosetests.jl"))
else
    @error "Julia's test/ directory not found, can't run Julia tests"
end

function test_path(test)
    t = split(test, '/')
    if t[1] in STDLIBS
        if length(t) == 2
            return joinpath(STDLIB_DIR, t[1], "test", t[2])
        else
            return joinpath(STDLIB_DIR, t[1], "test", "runtests")
        end
    else
        return joinpath(testdir, test)
    end
end

nstmts = 10^4  # very quick, aborts a lot
outputfile = "results.md"
i = 1
while i <= length(ARGS)
    global i
    a = ARGS[i]
    if a == "--nstmts"
        global nstmts = parse(Int, ARGS[i+1])
        deleteat!(ARGS, i:i+1)
    elseif a == "--output"
        global outputfile = ARGS[i+1]
        deleteat!(ARGS, i:i+1)
    else
        i += 1
    end
end

tests, _, exit_on_error, seed = choosetests(ARGS)

function spin_up_worker()
    p = addprocs(1)[1]
    remotecall_wait(include, p, "utils.jl")
    remotecall_wait(configure_test, p)
    return p
end

function spin_up_workers(n)
    procs = addprocs(n)
    @sync begin
        @async for p in procs
            remotecall_wait(include, p, "utils.jl")
            remotecall_wait(configure_test, p)
        end
    end
    return procs
end

# Really, we're just going to skip all the tests that run on node1
const node1_tests = String[]
function move_to_node1(t)
    if t in tests
        splice!(tests, findfirst(isequal(t), tests))
        push!(node1_tests, t)
    end
    nothing
end
move_to_node1("precompile")
move_to_node1("SharedArrays")
move_to_node1("stress")
move_to_node1("Distributed")

@testset "Julia tests" begin
    nworkers = min(Sys.CPU_THREADS, length(tests))
    println("Using $nworkers workers")
    results = Dict{String,Any}()
    tests0 = copy(tests)
    all_tasks = Union{Task,Nothing}[]
    try
        @sync begin
            for i = 1:nworkers
                @async begin
                    push!(all_tasks, current_task())
                    while length(tests) > 0
                        nleft = length(tests)
                        test = popfirst!(tests)
                        println(nleft, " remaining, starting ", test, " on task ", i)
                        local resp
                        fullpath = test_path(test) * ".jl"
                        try
                            resp = disable_sigint() do
                                p = spin_up_worker()
                                result = remotecall_fetch(run_test_by_eval, p, test, fullpath, nstmts)
                                rmprocs(p; waitfor=5)
                                result
                            end
                        catch e
                            if isa(e, InterruptException)
                                println("interrupting ", test)
                                break # rethrow(e)
                            end
                            resp = e
                            if isa(e, ProcessExitedException)
                                println("exited on ", test)
                            end
                        end
                        results[test] = resp
                        if resp isa Exception && exit_on_error
                            skipped = length(tests)
                            empty!(tests)
                        end
                    end
                    println("Task ", i, " complete")
                    all_tasks[i] = nothing
                end
            end
        end
    catch err
        isa(err, InterruptException) || rethrow(err)
        # If the test suite was merely interrupted, still print the
        # summary, which can be useful to diagnose what's going on
        foreach(all_tasks) do task
            try
                if isa(task, Task)
                    println("trying to interrupt ", task)
                    schedule(task, InterruptException(); error=true)
                end
            catch
            end
        end
        foreach(wait, all_tasks)
    end

    open(outputfile, "w") do io
        versioninfo(io)
        println(io, "Test run at: ", now())
        println(io)
        println(io, "Maximum number of statements per lowered expression: ", nstmts)
        println(io)
        println(io, "| Test file | Passes | Fails | Errors | Broken | Aborted blocks |")
        println(io, "| --------- | ------:| -----:| ------:| ------:| --------------:|")
        for test in tests0
            result = get(results, test, "")
            if isa(result, Tuple{Test.AbstractTestSet, Vector})
                ts, aborts = result
                passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = Test.get_test_counts(ts)
                naborts = length(aborts)
                println(io, "| ", test, " | ", passes+c_passes, " | ", fails+c_fails, " | ", errors+c_errors, " | ", broken+c_broken, " | ", naborts, " |")
            elseif isa(result, ProcessExitedException)
                println(io, "| ", test, " | ☠ | ☠ | ☠ | ☠ | ☠ |")
            else
                println(test, " => ", result)
                println(io, "| ", test, " | X | X | X | X | X |")
            end
        end
    end
end


text:
# This file was generated, do not modify it. # hide
train, test = partition(eachindex(y), 0.7, shuffle=true)
test[1:3]


text:
"""Returns a human readable string representing seconds in terms of days, hours, minutes or seconds."""
function readable_secs(secs::Real)
    days = Int(floor(secs/3600/24))
    hours = Int(floor((secs/3600) % 24))
    minutes = Int(floor((secs/60) % 60))
    seconds = Int(floor(secs%3600%60))

    if days > 0
        return "$(days)d, $(hours)h"
    elseif hours > 0
        return "$(hours)h, $(minutes)min"
    elseif minutes > 0
        return "$(minutes)min, $(seconds)s"
    else
        return "$(seconds)s"
    end
end

"""Estimates the total time the model integration will take."""
function duration_estimate(i,t,nt,progrtxt)
    time_per_step = (time()-t) / (i-nadvstep)
    time_total = Int(round(time_per_step*nt))
    time_to_go = Int(round(time_per_step*(nt-i)))

    s1 = "Model integration will take approximately "*readable_secs(time_total)*","
    s2 = "and is hopefully done on "*Dates.format(now() + Dates.Second(time_to_go),Dates.RFC1123Format)

    println(s1)     # print inline
    println(s2)
    if output == 1  # print in txt
        write(progrtxt,"\n"*s1*"\n")
        write(progrtxt,s2*"\n")
        flush(progrtxt)
    end
end

"""Returns a boolean whether the prognostic variables contains a NaN."""
function nan_detection(u::AbstractMatrix,v::AbstractMatrix,η::AbstractMatrix,sst::AbstractMatrix)
    #TODO include a check for Posits, are posits <: AbstractFloat?
    #TODO include check for tracer by other means than nan? (semi-Lagrange is unconditionally stable...)

    n_nan = sum(isnan.(u)) + sum(isnan.(v)) + sum(isnan.(η)) + sum(isnan.(sst))
    if n_nan > 0
        return true
    else
        return false
    end
end

"""Initialises the progress txt file."""
function feedback_ini()
    if output
        progrtxt = open(runpath*"progress.txt","w")
        s = "Starting juls run $run_id on "*Dates.format(now(),Dates.RFC1123Format)
        println(s)
        write(progrtxt,s*"\n")
        write(progrtxt,"Juls will integrate $(Ndays)days at a resolution of $(nx)x$(ny) with Δ=$(Δ/1e3)km\n")
        write(progrtxt,"Initial conditions are ")
        if initial_cond == "rest"
            write(progrtxt,"rest.\n")
        else
            write(progrtxt,"last time step of run $init_run_id.\n")
        end
        write(progrtxt,"Boundary conditions are $bc_x with lbc=$lbc.\n")
        write(progrtxt,"Numtype is "*string(Numtype)*".\n")
        write(progrtxt,"\nAll data will be stored in $runpath\n")
    else
        println("Starting juls on "*Dates.format(now(),Dates.RFC1123Format)*" without output.")
        progrtxt = nothing
    end

    return time(),progrtxt
end

"""Feedback function that calls duration estimate, nan_detection and progress."""
function feedback(u,v,η,sst,i,t,nt,nans_detected,progrtxt)
    if i == nadvstep # measure time after tracer advection executed once
        t = time()
    elseif i == min(2*nadvstep,nadvstep+50)
        # after the tracer advection executed twice or at least 50 steps
        duration_estimate(i,t,nt,progrtxt)
    end

    if !nans_detected
        if i % nout == 0    # only check for nans when output is produced
            nans_detected = nan_detection(u,v,η,sst)
            if nans_detected
                println(" NaNs detected at time step $i")
                if output == 1
                    write(progrtxt," NaNs detected at time step $i")
                    flush(progrtxt)
                end
            end
        end
    end

    if i > 100      # show percentage only after duration is estimated
        progress(i,nt,progrtxt)
    end

    return t,nans_detected
end

"""Finalises the progress txt file."""
function feedback_end(progrtxt,t::Real)
    s = " Integration done in "*readable_secs(time()-t)*"."
    println(s)
    if output
        write(progrtxt,"\n"*s[2:end]*"\n")  # close txt file with last output
        flush(progrtxt)
    end
end

"""Converts time step into percent for feedback."""
function progress(i,nt,progrtxt)
    if ((i+1)/nt*100 % 1) < (i/nt*100 % 1)  # update every 1 percent steps.
        percent = Int(round((i+1)/nt*100))
        print("\r\u1b[K")
        print("$percent%")
        if output && (percent % 5 == 0) # write out only every 5 percent step.
            write(progrtxt,"\n$percent%")
            flush(progrtxt)
        end
    end
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# Semiring definitions
abstract type Semiring end

import Base: zero, one, +, *, /
import Base: isapprox, parse, reverse, convert, get
isapprox(a::T,b::T) where { T <: Semiring } = isapprox(a.x,b.x)
reversetype(::Type{S}) where {S <: Semiring} = S
reverse(a::S) where { S <: Semiring } = a
reverseback(a::S) where { S <: Semiring } = a
get(a::S) where { S <: Semiring } = a.x

# properties
"""
`isleft(::Type{W})`

Check if the semiring type `W` satisfies:

``\\forall a,b,c \\in \\mathbb{W} : c \\otimes(a \\oplus b) = c \\otimes a \\oplus c \\otimes b``
"""
isleft(::Type{W}) where {W} = false             # ∀ a,b,c: c*(a+b) = c*a + c*b

"""
`isright(::Type{W})`

Check if the semiring type `W` satisfies:

``\\forall a,b,c \\in \\mathbb{W} : c \\otimes(a \\oplus b) = a \\otimes c \\oplus b \\otimes c``
"""
isright(::Type{W}) where {W}  = false           # ∀ a,b,c: c*(a+b) = a*c + b*c

"""
`isweaklydivisible(::Type{W})`

Check if the semiring type `W` satisfies:

``\\forall a,b \\in \\mathbb{W} \\ \\text{s.t.} \\ a \\oplus b \\neq \\bar{0} \\ \\exists z : x = (x \\oplus y ) \\otimes z``
"""
isweaklydivisible(::Type{W}) where {W}  = false # a+b ≂̸ 0: ∃ z s.t. x = (x+y)*z

"""
`ispath(::Type{W})`

Check if the semiring type `W` satisfies:

``\\forall a,b \\in \\mathbb{W}: a \\oplus b = a \\lor a \\oplus b = b``
"""
ispath(::Type{W}) where {W}  = false            # ∀ a,b: a+b = a or a+b=b

"""
`isidempotent(::Type{W})`

Check if the semiring type `W` satisfies:

``\\forall a \\in \\mathbb{W}: a \\oplus a = a``
"""
isidempotent(::Type{W}) where {W} = false       # ∀ a: a+a = a

"""
`iscommulative(::Type{W})`

Check if the semiring type `W` satisfies:

``\\forall a,b \\in \\mathbb{W}: a \\otimes b = b \\otimes a``
"""
iscommulative(::Type{W}) where {W} = false      # ∀ a,b: a*b = b*a
iscomplete(::Type{W}) where {W} = false

Base.show(io::IO, T::Semiring) = Base.show(io, T.x)


text:
#InspectDR: Draw legends with Cairo layer.
#-------------------------------------------------------------------------------

#==
===============================================================================#

#NOTE: linelength & linegap are absolute values
function legend_renderitem(ctx::CairoContext, wfrm::DWaveform, x::Float64, y::Float64,
	linelength::Float64, linegap::Float64)
	linestart = x+linegap
	lineend = linestart+linelength

	Cairo.set_source(ctx, wfrm.line.color)
	setlinestyle(ctx, wfrm.line.style, Float64(wfrm.line.width))
	drawline(ctx, Point2D(linestart, y), Point2D(lineend, y))
	drawglyph_safe(ctx, wfrm, Point2D(linestart+linelength/2, y))
	x = lineend+linegap #Compute new x
	Cairo.set_source(ctx, COLOR_BLACK)
	render(ctx, wfrm.id, Point2D(x, y), align=ALIGN_VCENTER|ALIGN_LEFT)
end

function legend_render(canvas::PCanvas2D, plot::Plot2D)
	const ctx = canvas.ctx
	const lstyle = plot.layout.legend

	Cairo.save(canvas.ctx)
	setfont(ctx, lstyle.font)
	(w, h) = text_dims(Cairo.text_extents(ctx, "M"))
	ypitch = h*(1+lstyle.vgap)

	#Compute absolute values:
	linegap = lstyle.linegap*w

	x = canvas.bb.xmax - _width(lstyle)
	y = canvas.graphbb.ymin + h/2

	for d in plot.display_data
		if "" == d.id; continue; end
		legend_renderitem(ctx, d, x, y, lstyle.linelength, linegap)
		y += ypitch
	end
	Cairo.restore(canvas.ctx)
end

#Last line


text:
struct PositionVector{T} <: NoAccessArray{T, 1}
    dim::Int
    vec::CuVector{T}
end

PositionVector(
    pos::Vector{Vector{T}}
) where T = PositionVector(
    length(pos),
    CuArray([
        [p[1] for p in pos];
        [p[2] for p in pos];
        [p[3] for p in pos]
    ])
)

@inline Base.size(v::PositionVector{T}) where T = (v.dim,)

@inline Adapt.adapt_storage(
    a::CUDA.Adaptor,
    v::PositionVector{T}
) where T = CuPositionVector(v.dim, Adapt.adapt_storage(a, v.vec))

struct TriangleVector{T} <: NoAccessArray{T, 1}
    dim::Int
    vec::CuVector{T}
end

TriangleVector(
    elements::Vector{Triangle{T}}
) where T = TriangleVector(
    length(elements),
    CuArray(NESSie.unpack(
        [[e.v1; e.v2; e.v3; e.normal; e.distorig; e.area] for e in elements]
    ))
)

@inline Base.size(v::TriangleVector{T}) where T = (v.dim,)

@inline Adapt.adapt_storage(
    a::CUDA.Adaptor,
    v::TriangleVector{T}
) where T = CuTriangleVector(v.dim, Adapt.adapt_storage(a, v.vec))


text:
using Base.Dates.AbstractTime

# Supertypes for atomic q types.
const SUPERTYPE = Dict(
    :_Bool     => Integer,
    :_Unsigned => Unsigned,
    :_Signed   => Signed,
    :_Float    => AbstractFloat,
    :_Text     => Any,
    :_Time     => Dates.TimeType,
    :_Period   => Dates.TimePeriod,
)
K_CLASSES = Type[]
# Create a parametrized type for each type class.
for (class, super) in SUPERTYPE
    @eval begin
        export $class
        struct $class{t,C,T} <: $(super)
            a::Array{C,0}
            function $class{t,C,T}(x::K_) where {t,C,T}
                a = asarray(x)
                t′ = xt(x)
                if t != -t′
                    throw(ArgumentError("type mismatch: t=$t, t′=$t′"))
                end
                new(a)
            end
        end
        ktypecode(::Type{$class{t,C,T}}) where {t,C,T} = -t
        K_new(x::$class{t,C,T}) where {t,C,T} = r1(kpointer(x))
        Base.convert(::Type{$class{t,C,T}}, x::$class{t,C,T}) where {t,C,T} = x
        function Base.convert(::Type{$class{t,C,T}}, x::C) where {t,C,T}
            r = $class{t,C,T}(ka(-t))
            r.a[] = x
            r
        end
        function Base.convert(::Type{$class{t,C,T}}, x) where {t,C,T}
            r = $class{t,C,T}(ka(-t))
            tmp1 = T(x)::T
            tmp2 = _cast(C, tmp1)::C
            r.a[] = tmp2
            r
        end
        push!(K_CLASSES, $class)
        Base.size(::$class) = ()
        Base.size(x::$class, d) =
            convert(Int, d) < 1 ? throw(BoundsError()) : 1
        Base.ndims(x::$class) = 0
        Base.ndims(x::Type{$class}) = 0
        Base.length(x::$class) = 1
        Base.endof(x::$class) = 1
        # payload
        Base.eltype(::Type{$class{t,C,T}}) where {t,C,T} = T
        Base.eltype(x::$class) = eltype(typeof(x))
        Base.getindex(x::$class) = _cast(eltype(x), x.a[])
        Base.pointer(x::$class{t,C,T}) where {t,C,T} = pointer(x.a)
        Base.show(io::IO, x::$class{t,C,T}) where {t,C,T} = print(io,
            "K(", repr(_cast(T, x.a[])), ")")
        Base.serialize(io::AbstractSerializer, x::$class) =
            _serialize(io, x, typeof(x))
    end
end
@eval const K_Atom = Union{$(K_CLASSES...)}

K_CLASS = Dict{Int8,Type}()
# Aliases for concrete scalar types.
for ti in TYPE_INFO
    ktype = Symbol("K_", ti.name)
    class = ti.class
    T = ti.jl_type
    @eval begin
        export $ktype
        global const $ktype =
              $class{$(ti.number),$(ti.c_type),$T}
        K_CLASS[$(ti.number)] = $class
        # Display type aliases by name in REPL.
        Base.show(io::IO, ::Type{$ktype}) = print(io, $(string(ktype)))
        # Conversion to Julia types
        Base.convert(::Type{$T}, x::$ktype) = _cast($T, x.a[])
        Base.promote_rule(x::Type{$T}, y::Type{$ktype}) = x
        # Disambiguate T -> K type conversions
        #Base.convert(::Type{$ktype}, x::$ktype) = x
    end
    if ti.class in [:_Signed, :_Integer]
        @eval Base.dec(x::$ktype, pad::Int=1) =
            string(dec(x.a[], pad), $(ti.letter))
    elseif ti.class === :_Unsigned
        @eval Base.hex(x::$ktype, pad::Int=1, neg::Bool=false) =
            hex(x.a[], pad, neg)
    end
    if ti.number != KC
        @eval ktypecode(::Type{$ti.jl_type}) = $ti.number
    end
end
ktypecode(::Type{Char}) = KC

function K_Atom(x::K_)
    t = -xt(x)
    ti = typeinfo(t)
    K_CLASS[t]{t,ti.c_type,ti.jl_type}(x)
end
# Chars and symbols are special
Base.convert(Char, x::K_char) = Char(x.a[])
Base.print(io::IO, x::K_char) = print(io, Char(x))
Base.show(io::IO, x::K_char) = print(io, "K(", repr(Char(x)), ")")
Base.show(io::IO, ::MIME"text/plain", x::K_char) = show(io, x)
Base.hash(x::K_char, h::UInt) = hash(x.a[], h)
Base.hash(x::K_symbol, h::UInt) = hash(x.a[], h)

Symbol(x::K_symbol) = convert(Symbol, x)
#Base.print(io::IO, x::K_symbol) = print(io, Symbol(x))
Base.show(io::IO, x::K_symbol) = print(io, "K(", repr(Symbol(x)), ")")
Base.show(io::IO, ::MIME"text/plain", x::K_symbol) = show(io, x)

Base.convert(::Type{K_timespan}, x::Integer) = K_timespan(TimeSpan(x))

include("promote_rules.jl")


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("CESMReader.jl")

using .CESMReader

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Ocn domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

function mreplace(x)
    return replace(x, missing=>NaN)
end

Dataset(parsed["domain-file"], "r") do ds
    global area = mreplace(ds["area"][:])
    global lat  = mreplace(ds["yc"][:])

    # Avoid grid points right on the equator.
    global tropic_NH_idx  = (   0.0 .<= lat .<= 20.0 )
    global tropic_SH_idx  = ( -20.0 .<= lat .<=  0.0 )

    global tropic_idx = tropic_NH_idx .| tropic_SH_idx

    global tropic_area = area[tropic_idx]
    global tropic_NH_area = area[tropic_NH_idx]
    global tropic_SH_area = area[tropic_SH_idx]

    global total_tropic_area    = sum(tropic_area)
    global total_tropic_NH_area = sum(tropic_NH_area)
    global total_tropic_SH_area = sum(tropic_SH_area)
end

let

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    PRECC, PRECL = getData(fh, ["PRECC", "PRECL"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    PREC = PRECC + PRECL

    if any(PREC .< 0)
        throw(ErrorException("OHOH"))
    end

    Nt = end_t - beg_t + 1

    global years = Int(Nt/12)
    global PAI         = zeros(Float64, years)
    global NH_mean     = zeros(Float64, years)
    global SH_mean     = zeros(Float64, years)
    global tropic_mean = zeros(Float64, years)

    for y = 1:years

        PREC_m = mean( PREC[ :, :, ((y-1)*12+1):(y*12)] , dims=3 )[:, :, 1]

        NH_m = sum(PREC_m[tropic_NH_idx] .* tropic_NH_area) / total_tropic_NH_area
        SH_m = sum(PREC_m[tropic_SH_idx] .* tropic_SH_area) / total_tropic_SH_area
        m    = sum(PREC_m[tropic_idx]    .* tropic_area   ) / total_tropic_area

        #P = PREC_m[tropic_idx]
        #println( format("Max: {:e}, min: {:e}, mean: {:e}, std: {:e}", maximum(P), minimum(P), mean(P), std(P)) )

        NH_mean[y] = NH_m
        SH_mean[y] = SH_m
        tropic_mean[y] = m

        PAI[y] = ( NH_m - SH_m ) / m

    end

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "year", years)

    for (varname, vardata, vardim, attrib) in [
        ("PAI", PAI, ("year",), Dict()),
        ("NH_mean", NH_mean, ("year",), Dict()),
        ("SH_mean", SH_mean, ("year",), Dict()),
        ("tropic_mean", tropic_mean, ("year",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
#### DomainDecompFuncs
# Provides a set of helper functions that operate on StateVec.

export domain_average!,
       distribute!,
       diagnose_environment!,
       total_covariance!

"""
    domain_average!(sv::StateVec, weight::StateVec, sv_idx, weight_idx, grid::Grid)

Compute domain average of field `sv_idx`, using weights `weight_idx`,
in state vector `sv`, the grid `grid` respectively.

Formulaically, a domain-averaged variable ``⟨ϕ⟩`` is computed from

``⟨ϕ⟩ = Σ_i a_i \\overline{ϕ}_i``

Where variable ``\\overline{ϕ}_i`` represents ``ϕ`` decomposed across multiple
sub-domains, which are weighted by area fractions ``a_i``.

Note that `domain_average!` is the inverse function of `distribute!`.
"""
function domain_average!(sv::StateVec, weight::StateVec, sv_idx, weight_idx, grid::Grid)
  gm, en, ud, sd, al = allcombinations(DomainIdx(sv))
  @inbounds for k in over_elems_real(grid)
    sv[sv_idx, k, gm] = 0
    @inbounds for i in sd
      sv[sv_idx, k, gm] += sv[sv_idx, k, i]*weight[weight_idx, k, i]
    end
  end
end

"""
    distribute!(sv::StateVec, grid::Grid, var_names)

Distributes values in the state vector `sv`, from the grid-mean
domain to the sub-domains.

Formulaically, a domain-decomposed variable ``\\overline{ϕ}_i`` is computed from

``\\overline{ϕ}_i = ⟨ϕ⟩``

Where variable ``⟨ϕ⟩`` is the domain-averaged variable, computed
across multiple sub-domains.

Note that `distribute!` is the inverse function of `domain_average!`.
"""
function distribute!(sv::StateVec, grid::Grid, var_names)
  !(var_names isa Tuple) && (var_names = (var_names,))
  gm, en, ud, sd, al = allcombinations(DomainIdx(sv))
  @inbounds for k in over_elems(grid)
    @inbounds for i in sd
      @inbounds for ϕ in var_names
        sv[ϕ, k, i] = sv[ϕ, k, gm]
      end
    end
  end
end

"""
    diagnose_environment!(sv::StateVec, grid::Grid, names)

Diagnose environment values in the state vector `sv`, from the grid-mean
domain and updraft variables from the decomposition constraint.

Formulaically:

``ϕ_en = (⟨ϕ⟩-Σ_{i=ud} a_i ϕ_i)/a_en``
"""
function diagnose_environment!(sv::StateVec, grid::Grid, weight_idx, var_names)
  !(var_names isa Tuple) && (var_names = (var_names,))
  gm, en, ud, sd, al = allcombinations(DomainIdx(sv))
  @inbounds for k in over_elems(grid)
    sv[weight_idx, k, en] = sv[weight_idx, k, gm] - sum([sv[weight_idx, k, i] for i in ud]...)
    a_en = sv[weight_idx, k, en]
    @inbounds for ϕ in var_names
      sv[ϕ, k, en] = (sv[ϕ, k, gm] - sum([sv[weight_idx, k, i]*sv[ϕ, k, i] for i in ud]...))/a_en
    end
  end
end

"""
    total_covariance!(dst::StateVec, src::StateVec, cv::StateVec, weights::StateVec,
                      dst_idxs, src_idxs, cv_idxs, weight_idx,
                      grid::Grid, decompose_ϕ_ψ::Function)

Computes the total covariance in state vector `dst`, given
 - `src` source state vector
 - `cv` state vector containing co-variances
 - `weights` state vector containing weights
 - `dst_idxs` indexes for destination state vector
 - `cv_idxs` indexes for state vector containing co-variances
 - `weight_idx` indexes for state vector containing weights
 - `grid` the grid
 - `decompose_ϕ_ψ` a function that receives the covariance index and
                   returns the indexes for each variable. For example:
                   `:ϕ_idx, :ψ_idx = decompose_ϕ_ψ(:cv_ϕ_ψ)`

Formulaically, a total covariance between variables ``ϕ`` and ``ψ`` is computed from

``⟨ϕ^*ψ^*⟩ = Σ_i a_i \\overline{ϕ_i'ψ_i'} + Σ_i Σ_j a_i a_j \\overline{ϕ}_i (\\overline{ψ}_i - \\overline{ψ}_j)``

Where variable ``\\overline{ϕ}_i`` represents ``ϕ`` decomposed across multiple
sub-domains, which are weighted by area fractions ``a_i``.
"""
function total_covariance!(tmp::StateVec, sv::StateVec, cv::StateVec,
                           tcv_ϕψ::Symbol, cv_ϕψ::Symbol,
                           a::Symbol, grid::Grid, decompose_ϕ_ψ::Function)
  gm, en, ud, sd, al = allcombinations(DomainIdx(sv))
  @inbounds for k in over_elems(grid)
    _ϕ, _ψ = decompose_ϕ_ψ(tcv_ϕψ)
    tmp[tcv_ϕψ, k] = 0
    @inbounds for i in sd
      tmp[tcv_ϕψ, k] += sv[a, k, i]*cv[cv_ϕψ, k, i]
      @inbounds for j in sd
        tmp[tcv_ϕψ, k] += sv[a, k, i]*sv[a, k, j]*sv[_ϕ, k, i]*(sv[_ψ, k, i] - sv[_ψ, k, j])
      end
    end
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using Statistics,Serialization,StatsPlots, DataFrames, HypothesisTests

filename_prefix_mutation = "../data/mutation/"
filename_prefix_125pinf = "../data/no_mutation_pinf125//"
filename_prefix_no_mutation = "../data/no_mutation/"

# set intial values as in coivd par ini
datat = 7 # interval of data collection
T = 582 #582  # time periods

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 84
# set periods for covid plots
inf_start_day = 1
inf_end_day = 400 # cut of after period 400
inf_start_week = 1
inf_end_week = 58 # 58*7 = 406 and 42*7 = 294

datapoint = fld(T,datat)
include("$(filename_prefix_mutation)beta_l5//covid_par_ini.jl")

### store data in txt file
io = open("stat_tests_results_btw_scenarios.txt", "w")
println(io, "Statisitcal tests: p-values from Mann-Whitney U tests\n")

# POINT A, B, C, D, E in Table 1
paras = ["beta_l50","beta_l5","beta_l5_alpha_l050", "beta_l5_alpha_l050_alpha_o050"]

################# with mutation
### get point A
p=paras[1]
data_prefix = string("$(filename_prefix_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_A_mut = (results -> results[:togdploss]).(worker_results)
totcas_A_mut = (results -> results[:totcas]).(worker_results)

### get point B
p=paras[2]
data_prefix = string("$(filename_prefix_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_B_mut = (results -> results[:togdploss]).(worker_results)
totcas_B_mut = (results -> results[:totcas]).(worker_results)

### get point C
p=paras[3]
data_prefix = string("$(filename_prefix_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_C_mut = (results -> results[:togdploss]).(worker_results)
totcas_C_mut = (results -> results[:totcas]).(worker_results)

### get point D
p=paras[4]
data_prefix = string("$(filename_prefix_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_D_mut = (results -> results[:togdploss]).(worker_results)
totcas_D_mut = (results -> results[:totcas]).(worker_results)

######### higher pinf 125
### get point A
p=paras[1]
data_prefix = string("$(filename_prefix_125pinf)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_A125 = (results -> results[:togdploss]).(worker_results)
totcas_A125 = (results -> results[:totcas]).(worker_results)

### get point B
p=paras[2]
data_prefix = string("$(filename_prefix_125pinf)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_B125 = (results -> results[:togdploss]).(worker_results)
totcas_B125 = (results -> results[:totcas]).(worker_results)

### get point C
p=paras[3]
data_prefix = string("$(filename_prefix_125pinf)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_C125 = (results -> results[:togdploss]).(worker_results)
totcas_C125 = (results -> results[:totcas]).(worker_results)

### get point D
p=paras[4]
data_prefix = string("$(filename_prefix_125pinf)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_D125 = (results -> results[:togdploss]).(worker_results)
totcas_D125 = (results -> results[:totcas]).(worker_results)

### no mutation
### get point A
p=paras[1]
data_prefix = string("$(filename_prefix_no_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_A = (results -> results[:togdploss]).(worker_results)
totcas_A = (results -> results[:totcas]).(worker_results)

### get point B
p=paras[2]
data_prefix = string("$(filename_prefix_no_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_B = (results -> results[:togdploss]).(worker_results)
totcas_B = (results -> results[:totcas]).(worker_results)

### get point C
p=paras[3]
data_prefix = string("$(filename_prefix_no_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_C = (results -> results[:togdploss]).(worker_results)
totcas_C = (results -> results[:totcas]).(worker_results)

### get point D
p=paras[4]
data_prefix = string("$(filename_prefix_no_mutation)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_D = (results -> results[:togdploss]).(worker_results)
totcas_D = (results -> results[:totcas]).(worker_results)

println(io,"\n MUTATION vs NO MUTATION\n")
println(io,"\n TESTS FOR GDP\n")

HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_A_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_A_mut))
println(io,"\n pvalue for GDP A vs A_mut: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_B, gdploss_B_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_B, gdploss_B_mut))
println(io,"\n pvalue for GDP B vs B_mut: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_C, gdploss_C_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_C, gdploss_C_mut))
println(io,"\n pvalue for GDP C vs C_mut: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_D, gdploss_D_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_D, gdploss_D_mut))
println(io,"\n pvalue for GDP D vs D_mut: ", p)

println(io,"\n \n TESTS FOR TOT CAS Table 1\n")

HypothesisTests.MannWhitneyUTest(totcas_A, totcas_A_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_A, totcas_A_mut))
println(io,"\n pvalue for CAS A vs A_mut: ", p)

HypothesisTests.MannWhitneyUTest(totcas_B, totcas_B_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_B, totcas_B_mut))
println(io,"\n pvalue for CAS B vs B_mut: ", p)

HypothesisTests.MannWhitneyUTest(totcas_C, totcas_C_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_C, totcas_C_mut))
println(io,"\n pvalue for CAS C vs C_mut: ", p)

HypothesisTests.MannWhitneyUTest(totcas_D, totcas_D_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_D, totcas_D_mut))
println(io,"\n pvalue for CAS D vs D_mut: ", p)

println(io,"\n \n RATIO GDP\n")

println(io,"\n ratio for GDP A vs A_mut: ", mean(gdploss_A_mut)/mean(gdploss_A))
println(io,"\n ratio for GDP B vs B_mut: ", mean(gdploss_B_mut)/mean(gdploss_B))
println(io,"\n ratio for GDP C vs C_mut: ", mean(gdploss_C_mut)/mean(gdploss_C))
println(io,"\n ratio for GDP D vs D_mut: ", mean(gdploss_D_mut)/mean(gdploss_D))

println(io,"\n \n RATIO CAS\n")

println(io,"\n ratio for CAS A vs A_mut: ", mean(totcas_A_mut)/mean(totcas_A))
println(io,"\n ratio for CAS B vs B_mut: ", mean(totcas_B_mut)/mean(totcas_B))
println(io,"\n ratio for CAS C vs C_mut: ", mean(totcas_C_mut)/mean(totcas_C))
println(io,"\n ratio for CAS D vs D_mut: ", mean(totcas_D_mut)/mean(totcas_D))

println(io,"\n \n DIFF GDP\n")

println(io,"\n ratio for GDP A vs A_mut: ", (mean(gdploss_A_mut) - mean(gdploss_A)))
println(io,"\n ratio for GDP B vs B_mut: ", (mean(gdploss_B_mut) - mean(gdploss_B)))
println(io,"\n ratio for GDP C vs C_mut: ", (mean(gdploss_C_mut) - mean(gdploss_C)))
println(io,"\n ratio for GDP D vs D_mut: ", (mean(gdploss_D_mut) - mean(gdploss_D)))

println(io,"\n \n DIFF CAS\n")

println(io,"\n ratio for CAS A vs A_mut: ", (mean(totcas_A_mut) - mean(totcas_A)))
println(io,"\n ratio for CAS B vs B_mut: ", (mean(totcas_B_mut) - mean(totcas_B)))
println(io,"\n ratio for CAS C vs C_mut: ", (mean(totcas_C_mut) - mean(totcas_C)))
println(io,"\n ratio for CAS D vs D_mut: ", (mean(totcas_D_mut) - mean(totcas_D)))

println(io,"\n MUTATION vs HIGHER PINF 125\n")
println(io,"\n TESTS FOR GDP\n")

HypothesisTests.MannWhitneyUTest(gdploss_A125, gdploss_A_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A125, gdploss_A_mut))
println(io,"\n pvalue for GDP A125 vs A_mut: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_B125, gdploss_B_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_B125, gdploss_B_mut))
println(io,"\n pvalue for GDP B125 vs B_mut: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_C125, gdploss_C_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_C125, gdploss_C_mut))
println(io,"\n pvalue for GDP C125 vs C_mut: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_D125, gdploss_D_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_D125, gdploss_D_mut))
println(io,"\n pvalue for GDP D125 vs D_mut: ", p)

println(io,"\n \n TESTS FOR TOT CAS Table 1\n")

HypothesisTests.MannWhitneyUTest(totcas_A125, totcas_A_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_A125, totcas_A_mut))
println(io,"\n pvalue for CAS A125 vs A_mut: ", p)

HypothesisTests.MannWhitneyUTest(totcas_B125, totcas_B_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_B125, totcas_B_mut))
println(io,"\n pvalue for CAS B125 vs B_mut: ", p)

HypothesisTests.MannWhitneyUTest(totcas_C125, totcas_C_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_C125, totcas_C_mut))
println(io,"\n pvalue for CAS C125 vs C_mut: ", p)

HypothesisTests.MannWhitneyUTest(totcas_D125, totcas_D_mut)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_D125, totcas_D_mut))
println(io,"\n pvalue for CAS D125 vs D_mut: ", p)

println(io,"\n \n RATIO GDP\n")

println(io,"\n ratio for GDP A125 vs A_mut: ", mean(gdploss_A_mut)/mean(gdploss_A125))
println(io,"\n ratio for GDP B125 vs B_mut: ", mean(gdploss_B_mut)/mean(gdploss_B125))
println(io,"\n ratio for GDP C125 vs C_mut: ", mean(gdploss_C_mut)/mean(gdploss_C125))
println(io,"\n ratio for GDP D125 vs D_mut: ", mean(gdploss_D_mut)/mean(gdploss_D125))

println(io,"\n \n RATIO CAS\n")

println(io,"\n ratio for CAS A125 vs A_mut: ", mean(totcas_A_mut)/mean(totcas_A125))
println(io,"\n ratio for CAS B125 vs B_mut: ", mean(totcas_B_mut)/mean(totcas_B125))
println(io,"\n ratio for CAS C125 vs C_mut: ", mean(totcas_C_mut)/mean(totcas_C125))
println(io,"\n ratio for CAS D125 vs D_mut: ", mean(totcas_D_mut)/mean(totcas_D125))

println(io,"\n \n DIFF GDP\n")

println(io,"\n ratio for GDP A125 vs A_mut: ", (mean(gdploss_A_mut) - mean(gdploss_A125)))
println(io,"\n ratio for GDP B125 vs B_mut: ", (mean(gdploss_B_mut) - mean(gdploss_B125)))
println(io,"\n ratio for GDP C125 vs C_mut: ", (mean(gdploss_C_mut) - mean(gdploss_C125)))
println(io,"\n ratio for GDP D125 vs D_mut: ", (mean(gdploss_D_mut) - mean(gdploss_D125)))

println(io,"\n \n DIFF CAS\n")

println(io,"\n ratio for CAS A125 vs A_mut: ", (mean(totcas_A_mut) - mean(totcas_A125)))
println(io,"\n ratio for CAS B125 vs B_mut: ", (mean(totcas_B_mut) - mean(totcas_B125)))
println(io,"\n ratio for CAS C125 vs C_mut: ", (mean(totcas_C_mut) - mean(totcas_C125)))
println(io,"\n ratio for CAS D125 vs D_mut: ", (mean(totcas_D_mut) - mean(totcas_D125)))

### close txt file
close(io)


text:
module GeometricIntegrators

    using Reexport

    include("Utils.jl")
    using .Utils

    include("Config.jl")
    @reexport using .Config
    include("CommonFunctions.jl")
    @reexport using .CommonFunctions
    include("Equations.jl")
    @reexport using .Equations
    include("Solutions.jl")
    @reexport using .Solutions
    include("Interpolation.jl")
    @reexport using .Interpolation
    include("Solvers.jl")
    @reexport using .Solvers
    include("BasisFunctions.jl")
    @reexport using .BasisFunctions
    include("Quadratures.jl")
    @reexport using .Quadratures
    include("Discontinuities.jl")
    @reexport using .Discontinuities
    include("Integrators.jl")
    @reexport using .Integrators
    include("Simulations.jl")
    @reexport using .Simulations
    include("Tableaus.jl")
    @reexport using .Tableaus

    include("TestProblems.jl")
    export TestProblems

    function __init__()
        add_config(:verbosity, 1)
    end

end


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    l = Exo2SplitBinion(Exo2Real(1,2,3,4))
    r = Exo2SplitBinion(1,0,2,0,3,0,4)
    l == r
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo2SplitBinion(a))
end

@test begin
    a = 1
    !isreal(Exo2SplitBinion(a, a, a, a))
end

@test begin
    z = random(Exo2SplitBinion{Int})
    z == +(z)
end

@test begin
    a = CayleyDickson.randomBigFloat()
    real(Exo2SplitBinion(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, Exo2SplitBinion(1.0, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1.0, s: 2.0, a: 3.0, sa: 4.0, b: 5.0, sb: 6.0, ab: 7.0, (sa)b: 8.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2SplitBinion(1.0, 2, 3, 4, 5, 6, 7))
    l = String(take!(io))
    r = "[1: 1.0, s: 2.0, a: 3.0, sa: 4.0, b: 5.0, sb: 6.0, ab: 7.0, (sa)b: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2SplitBinion(1.0, 2, 3, 4, 5, 6))
    l = String(take!(io))
    r = "[1: 1.0, s: 2.0, a: 3.0, sa: 4.0, b: 5.0, sb: 6.0, ab: 0.0, (sa)b: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2SplitBinion(1.0, 2, 3, 4, 5))
    l = String(take!(io))
    r = "[1: 1.0, s: 2.0, a: 3.0, sa: 4.0, b: 5.0, sb: 0.0, ab: 0.0, (sa)b: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2SplitBinion(1.0, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1.0, s: 2.0, a: 3.0, sa: 4.0, b: 0.0, sb: 0.0, ab: 0.0, (sa)b: 0.0]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    d = CayleyDickson.randomBigFloat()
    l = Exo2SplitBinion(Exo1SplitBinion(a, b, c, d))
    r = Exo2SplitBinion(a, b, c, d)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = Exo2SplitBinion(SplitBinion(a, b))
    r = Exo2SplitBinion(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = Exo2SplitBinion(a, b, c)
    r = Exo2SplitBinion(a, b, c, 0)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigFloat})
    y = random(Exo2SplitBinion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo2SplitBinion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Exo1SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo2SplitBinion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Exo1SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l != r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = alternatorL(x, y)
    r = zero(Exo2SplitBinion{BigInt})
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = alternatorR(x, y)
    r = zero(Exo2SplitBinion{BigInt})
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo2SplitBinion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{Rational{BigInt}})
    y = random(Exo2SplitBinion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{Rational{BigInt}})
    y = random(Exo2SplitBinion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo2SplitBinion{Rational{BigInt}})
    y = random(Exo2SplitBinion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo2SplitBinion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Exo1SplitBinion{BigInt})
    x = random(Exo2SplitBinion{BigInt})
    l = commutator(a, x)
    r = zero(Exo2SplitBinion{BigInt})
    l != r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    z = random(Exo2SplitBinion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    quadrance(x) > big(0) || quadrance(x) < big(0) ||quadrance(x) == big(0)
end

@test begin
    x = random(Exo2SplitBinion{BigInt})
    y = random(Exo2SplitBinion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo2SplitBinion(0, 0.0))
end

@test_throws ErrorException begin
    inv(Exo2SplitBinion(1, 1))
end

@test_throws ErrorException begin
    random(Exo2SplitBinion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo2SplitBinion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo2SplitBinion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo2SplitBinion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo2SplitBinion{BigInt})
    z == z + zero(Exo2SplitBinion{BigInt})
end

@test begin
    z = random(Exo2SplitBinion{BigInt})
    z == z * one(Exo2SplitBinion{BigInt})
end

@test begin
    z = random(Exo2SplitBinion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo2SplitBinion{BigInt})
    z == z * one(z)
end


text:
#=

running tribe model from the command line

Ignacio Quintero Mächler

t(-_-t)

August 17 2017
=#

"""
    parse_commandline()

Parses command line arguments to be used in compete
"""
function parse_commandline()

  s = ArgParseSettings(prog        = "Tapestree",
                       description = "Inference of competition on trait evolution along biogeographic history",
                       version     = "0.1",
                       add_version = true)

  @add_arg_table s begin

    # optional arguments
    "--niter", "-i"
      help = "number of mcmc iterations"
      arg_type = Int64
      default = 500000

    "--nburn", "-b"
      help = "number of iterations in burning and tuning phase"
      arg_type = Int64
      default = 500000

    "--nthin", "-t"
      help = "sample every nthin iterations"
      arg_type = Int64
      default = 1000

    "--min_dt", "-d"
      help = "minimum percentage of discrete sampling"
      arg_type = Float64
      default = 0.01

    "--comp_x_prior"
      help = "gaussian prior for trait competition"
      arg_type = Tuple{Float64,Float64}
      default = 0.0, 10.0
      dest_name = "ωxprior"

    "--comp_gain_prior"
      help = "gaussian prior for effect of competition on area gain"
      arg_type = Tuple{Float64,Float64}
      default = 0.0, 10.0
      dest_name = "ω1prior"

    "--comp_loss_prior"
      help = "gaussian prior for effect of competition on area loss"
      arg_type = Tuple{Float64,Float64}
      default = 0.0, 10.0
      dest_name = "ω0prior"

    "--sigma_prior"
      help = "exponential prior for rate of evolution"
      arg_type = Float64
      default = 0.1
      dest_name = "σ²prior"

    "--lambda_prior"
      help = "exponential prior for gain and loss rates"
      arg_type = Float64
      default = 0.1
      dest_name = "λprior"

    "--stbrl"
      help = "stem branch length"
      arg_type = Float64
      default = 1.0

    # flags
    "--fix_biocomp", "-f"
      help = "fix to 0 competition effect on gain and losses"
      action = :store_true
      dest_name = "fix_ω1_ω0"

    # required arguments
    "tree_file"
      help = "full path to phylogenetic tree file"
      required = true
    "data_file"
      help = "full path to data file"
      required = true
    "out_file"
      help = "full path to write file"
      required = true

  end

  parsed_args = parse_args(s, as_symbols = true)

  println("Parsed args:")
  for (arg,val) in parsed_args
      println("  $arg  =>  $val")
  end

  return parsed_args
end

# parse commands
comp_args = parse_commandline()

using Tapestree

# run compete
tribe(comp_args[:tree_file],
      comp_args[:data_file],
      comp_args[:out_file];
      filter((u,v) -> u != :tree_file &&
                      u != :data_file &&
                      u != :out_file, comp_args)...)

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
function interp_face!(fδ, f::T, ll::T1, lr::T1) where {T<:AbstractVector,T1<:AbstractVector}
    fδ[1] = dot(f, ll)
    fδ[2] = dot(f, lr)
end

function interp_face!(fδ, f::T, ll::T1, lr::T1) where {T<:AbstractMatrix,T1<:AbstractVector}
    @views for i in axes(f, 1)
        interp_face!(fδ[i, :], f[i, :], ll, lr)
    end
end


text:
module HTTP2

const Headers = Vector{Tuple{String,String}}
const Maybe{T} = Union{Nothing, T} where {T}

Base.convert(::Type{Headers}, d::Dict) = [(a,b) for (a,b) in d]

bytearr(a::Vector{UInt8}) = a
bytearr(cs::Base.CodeUnits{UInt8,String}) = convert(Vector{UInt8}, cs)
bytearr(s::String) = bytearr(codeunits(s))

include("Frame.jl")
include("Session.jl")

end # module HTTP2


text:
@testset "geometricflux.jl" begin

    @testset "GraphConv" begin
        gnn = GeometricFlux.GraphConv(3 => 3)

        graphs = Matrix.(adjacency_matrix.([random_regular_graph(10, 4) for i = 1:4]))
        nodeFeatures = [rand(Float32, 3, 10) for i = 1:4]
        featuredGraphs = [FeaturedGraph(g; nf=nf) for (g, nf) in zip(graphs, nodeFeatures)]
        batchedFG = SeaPearl.BatchedFeaturedGraph{Float32}(
            cat(graphs...; dims=3);
            nf=cat(nodeFeatures...;dims=3)
        )

        classicOutput = cat(GeometricFlux.FeatureSelector(:node).(gnn.(featuredGraphs))...;dims=3)
        batchedOutput = gnn(batchedFG).nf

        @test all(abs.(classicOutput .- batchedOutput) .< 1e-5)

        grad = gradient(batchedFG) do inp
            gnn(inp).nf[1, 3, 2]
        end
        grad = grad[1]

        @test all(grad.nf[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.nf[:, :, 2] .== 0)
        @test all(grad.graph[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.graph[:, :, 2] .== 0)
    end

    # TODO find why it won't work on CPU
    """
    @testset "GraphNorm" begin
        gnn = SeaPearl.GraphNorm(3)

        graphs = Matrix.(adjacency_matrix.([random_regular_graph(10, 4) for i = 1:4]))
        nodeFeatures = [rand(Float32, 3, 10) for i = 1:4]
        batchedFG = SeaPearl.BatchedFeaturedGraph{Float32}(
            cat(graphs...; dims=3);
            nf=cat(nodeFeatures...;dims=3)
        )
        batchedOutput = gnn(batchedFG).nf

        @test size(batchedOutput) = (3, 10, 3)

        grad = gradient(batchedFG) do inp
            gnn(inp).nf[1, 1, 1]
        end
        grad = grad[1]

        @test all(grad.nf[:, :, 2:4] .== 0)
        @test all(grad.nf[:, 2:10, 1] .== 0)
        @test !all(grad.nf[:, 1, 1] .== 0)
        @test all(grad.graph .== 0)
    end
    """
end


text:
# The master process listens on a well-known port
# Launched workers connect to the master and redirect their STDOUTs to the same
# Workers can join and leave the cluster on demand.

export ElasticManager, elastic_worker

const HDR_COOKIE_LEN = Distributed.HDR_COOKIE_LEN

struct ElasticManager <: ClusterManager
    active::Dict{Int, WorkerConfig}        # active workers
    pending::Channel{TCPSocket}          # to be added workers
    terminated::Set{Int}             # terminated worker ids
    topology::Symbol
    sockname

    function ElasticManager(;addr=IPv4("127.0.0.1"), port=9009, cookie=nothing, topology=:all_to_all)
        Distributed.init_multi()
        cookie !== nothing && cluster_cookie(cookie)

        l_sock = listen(addr, port)

        lman = new(Dict{Int, WorkerConfig}(), Channel{TCPSocket}(typemax(Int)), Set{Int}(), topology, getsockname(l_sock))

        @async begin
            while true
                let s = accept(l_sock)
                    @async process_worker_conn(lman, s)
                end
            end
        end

        @async process_pending_connections(lman)

        lman
    end
end

ElasticManager(port) = ElasticManager(;port=port)
ElasticManager(addr, port) = ElasticManager(;addr=addr, port=port)
ElasticManager(addr, port, cookie) = ElasticManager(;addr=addr, port=port, cookie=cookie)

function process_worker_conn(mgr::ElasticManager, s::TCPSocket)
    # Socket is the worker's STDOUT
    wc = WorkerConfig()
    wc.io = s

    # Validate cookie
    cookie = read(s, HDR_COOKIE_LEN)
    if length(cookie) < HDR_COOKIE_LEN
        error("Cookie read failed. Connection closed by peer.")
    end
    self_cookie = cluster_cookie()
    for i in 1:HDR_COOKIE_LEN
        if UInt8(self_cookie[i]) != cookie[i]
            println(i, " ", self_cookie[i], " ", cookie[i])
            error("Invalid cookie sent by remote worker.")
        end
    end

    put!(mgr.pending, s)
end

function process_pending_connections(mgr::ElasticManager)
    while true
        wait(mgr.pending)
        try
            addprocs(mgr; topology=mgr.topology)
        catch e
            showerror(stderr, e)
            Base.show_backtrace(stderr, Base.catch_backtrace())
        end
    end
end

function launch(mgr::ElasticManager, params::Dict, launched::Array, c::Condition)
    # The workers have already been started.
    while isready(mgr.pending)
        wc=WorkerConfig()
        wc.io = take!(mgr.pending)
        push!(launched, wc)
    end

    notify(c)
end

function manage(mgr::ElasticManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :register
        mgr.active[id] = config
    elseif  op == :deregister
        delete!(mgr.active, id)
        push!(mgr.terminated, id)
    end
end

function Base.show(io::IO, mgr::ElasticManager)
    iob = IOBuffer()

    println(iob, "ElasticManager:")
    print(iob, "  Active workers : [ ")
    for id in sort(collect(keys(mgr.active)))
        print(iob, id, ",")
    end
    seek(iob, position(iob)-1)
    println(iob, "]")

    println(iob, "  Number of workers to be added  : ", Base.n_avail(mgr.pending))

    print(iob, "  Terminated workers : [ ")
    for id in sort(collect(mgr.terminated))
        print(iob, id, ",")
    end
    seek(iob, position(iob)-1)
    println(iob, "]")

    print(io, String(take!(iob)))
end

# Does not return. If executing from a REPL try
# @async connect_to_cluster(.....)
# addr, port that a ElasticManager on the master processes is listening on.
function elastic_worker(cookie, addr="127.0.0.1", port=9009; stdout_to_master=true)
    c = connect(addr, port)
    write(c, rpad(cookie, HDR_COOKIE_LEN)[1:HDR_COOKIE_LEN])
    stdout_to_master && redirect_stdout(c)
    start_worker(c, cookie)
end


text:
using BinaryBuilder

cuda_version = v"10.0.130"  # NOTE: could be less specific

sources_linux = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.0_20191031/cudnn-10.0-linux-x64-v7.6.5.32.tgz",
                   "28355e395f0b2b93ac2c83b61360b35ba6cd0377e44e78be197b6b61b4b492ba")
]
sources_macos = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.0_20191031/cudnn-10.0-osx-x64-v7.6.5.32.tgz",
                  "6fa0b819374da49102e285ecf7fcb8879df4d0b3cc430cc8b781cdeb41009b47")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.0_20191031/cudnn-10.0-windows10-x64-v7.6.5.32.zip",
                  "2767db23ae2cd869ac008235e2adab81430f951a92a62160884c80ab5902b9e8")
]

include("../common.jl")