text:
module Laplace

## Package imports
#==========================================================================================#

using StaticArrays
using LinearAlgebra

import ..MathTools: Point2D, Point3D, magnitude, angle

import ..AeroMDAO: velocity

abstract type AbstractLaplace end

## Legacy (to be removed?)
#==========================================================================================#

# Performs velocity and potential computations for an object on a grid
grid_data(object :: AbstractLaplace, xs) = velocity(object, xs), potential(object, xs)
velocity(object :: AbstractLaplace, xs) = map(x -> velocity(object, x...), xs)
potential(object :: AbstractLaplace, xs) = map(x -> potential(object, x...), xs)

# Performs velocity and potential calculations on a grid
function grid_data(objects :: Vector{<: AbstractLaplace}, xs)
    vels = foldl((v1, v2) -> [ u .+ v for (u, v) ∈ zip(v1, v2) ], velocity(object, xs) for object ∈ objects)
    pots = foldl((v1, v2) -> v1 + v2, potential(object, xs) for object ∈ objects)

    vels, pots
end

## 2D singularities
#============================================#

struct Singularity2D{T <: Real} <: AbstractLaplace
    strength :: T
    r        :: Point2D{T}
end

# Getters
strength(s :: Singularity2D) = s.strength
x(s :: Singularity2D)        = s.r.x
y(s :: Singularity2D)        = s.r.y

source_velocity(src :: Singularity2D, x, y) = SVector(strength(src) / (2π) * (x - y(src)) / ((x - y(src))^2 + (y - x(src))^2), str / (2π) * (y - x(src)) / ((x - y(src))^2 + (y - x(src))^2))
source_potential(src :: Singularity2D, x, y) = strength(src) / (4π) * log((x - y(src))^2 + (y - x(src))^2)
source_stream(src :: Singularity2D, x, y) = strength(src) / (2π) * atan(y - x(src), x - y(src))

doublet_velocity(dub :: Singularity2D, x, y) = SVector(strength(dub) / (2π) * ((x - y(dub))^2 - (y - x(dub))^2) / ((x - y(dub))^2 + (y - x(dub))^2)^2, - strength(dub) / (2π) * 2 * (x - y(dub)) * (y - x(dub)) / ((x - y(dub))^2 + (y - x(dub))^2)^2)
doublet_potential(dub :: Singularity2D, x, y) = -strength(dub) / (2π) * (y - x(dub)) / ((x - y(dub))^2 + (y - x(dub))^2)
doublet_stream(dub :: Singularity2D, x, y) = -strength(dub) / (2π) * (y - x(dub)) / ((x - y(dub))^2 + (y - x(dub))^2)

vortex_velocity(vor :: Singularity2D, x, y) = SVector(-strength(vor) / (2π) * (y - x(vor)) / ((x - y(vor))^2 + (y - x(vor))^2), str / (2π) * (x - y(vor)) / ((x - y(vor))^2 + (y - x(vor))^2))
vortex_potential(vor :: Singularity2D, x, y) = strength(vor) / (2π) * atan(y - x(vor), x - y(vor))
vortex_stream(vor :: Singularity2D, x, y) = -strength(vor) / (4π) * log((x - y(vor))^2 + (y - x(vor))^2)

struct Uniform2D{T <: Real} <: AbstractLaplace
    magnitude :: T
    angle     :: T
    Uniform2D{T}(mag, ang) where T <: Real = new(mag, deg2rad(ang))
end

magnitude(uni :: Uniform2D) = uni.magnitude
angle(uni :: Uniform2D)     = uni.angle

Uniform2D(mag :: T, ang :: T) where T <: Real = Uniform2D{T}(mag, ang)
Uniform2D(mag, ang) = Uniform2D(promote(mag, ang)...)
Uniform2D(velocity) = Uniform2D((sqrt ∘ sum)(velocity.^2), atan(velocity[2], velocity[1]))
Uniform2D(; angle)  = Uniform2D(1., angle)

velocity(uni :: Uniform2D) = let (sa, ca) = sincos(uni.angle); uni.magnitude * SVector(ca, sa) end
potential(uni :: Uniform2D, x, y) = uni.magnitude * (x * cos(uni.angle) + y * sin(uni.angle))
stream(uni :: Uniform2D, x, y)    = uni.magnitude * (y * cos(uni.angle) - x * sin(uni.angle))

## 3D singularities
#============================================#

# struct Singularity3D{T <: Real} <: AbstractLaplace
#       str :: T
#       r   :: Point3D{T}
# end

# source_velocity(src :: Source2D, x, y, z)
# source_potential(src :: Source2D, x, y, z)
# source_stream(src :: Source2D, x, y, z)

struct DoubletLine3D{T <: Real} <: AbstractLaplace
    strength :: T
    r1       :: SVector{3,T}
    r2       :: SVector{3,T}
    eta      :: SVector{3,T}
end

function doublet_influence(r, φ, η)
    r_φ = dot(r, φ)
    r_η = dot(r, η)
    den = (norm(r)^2 - dot(r, φ)^2 ) * r

    ((r_φ * η + r_η * φ) * den - (den * r / norm(r)^2 + 2 * (r - r_φ * η) * r) * r_φ * r_η) / den^2
end

function velocity(src :: DoubletLine3D)
    l = normalize(src.r2 - src.r1)

    f(x) = src.strength / 4π * (doublet_influence(x - src.r2, l, src.eta) - doublet_influence(x - src.r1, l, src.eta))
end

## Freestream
#============================================#
"""
    freestream_to_cartesian(r, θ, φ)

Convert freestream flow (spherical polar) coordinates to Cartesian coordinates.
"""
freestream_to_cartesian(r, θ, φ) = r * SVector(cos(θ) * cos(φ), -sin(φ), sin(θ) * cos(φ))

"""
    cartesian_to_freestream(U)

Convert Cartesian coordinates to freestream (spherical polar) flow coordinates.
"""
cartesian_to_freestream(U) = SVector(norm(U), -atand(U[3], U[1]), -atand(U[2], √(U[1]^2 + U[3]^2)))

# 2D versions
cartesian_to_freestream(u, w) = magnitude(u, w), angle(u, w)
freestream_to_cartesian(V, α) = V * cos(α), V * sin(α)

end


text:
using BinaryBuilder

# Collection of sources required to build Ogg
sources = [
    "JuliaInterop"
]

# Bash recipe for building across all platforms
function getscript(version)
    shortversion = version[1:3]
    return """
    Julia_ROOT=/usr/local

    # Download julia
    cd /usr/local
    curl -L "https://github.com/JuliaPackaging/JuliaBuilder/releases/download/$version/julia-$version-\$target.tar.gz" | tar -zx --strip-components=1

    # Build libcxxwrap
    cd \$WORKSPACE/srcdir/libcxxwrap-julia*
    mkdir build && cd build
    cmake -DJulia_ROOT=\$Julia_ROOT -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=/opt/\$target/\$target.toolchain -DCMAKE_CXX_FLAGS="-march=x86-64" -DCMAKE_INSTALL_PREFIX=\${prefix} ..
    VERBOSE=ON cmake --build . --config Release --target install
    if [[ "\$target" == "x86_64-linux-gnu" ]]; then
        ctest -V
    fi
    """
end

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
    #Linux(:i686)
]

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libcxxwrap", :libcxxwrap),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
version_number = get(ENV, "TRAVIS_TAG", "")
if version_number == ""
    version_number = "v0.99"
end
build_tarballs(ARGS, "libcxxwrap-julia-0.7", VersionNumber(version_number), sources, getscript("0.7.0"), platforms, products, dependencies)
build_tarballs(ARGS, "libcxxwrap-julia-1.0", VersionNumber(version_number), sources, getscript("1.0.0"), platforms, products, dependencies)

ext: jl

max_stars_repo_name: RelationalAI-oss/libcxxwrap-julia

max_issues_repo_name: RelationalAI-oss/libcxxwrap-julia

max_forks_repo_name: RelationalAI-oss/libcxxwrap-julia


text:
##### Packages #####

## Import package manager
import Pkg

## Activate project
Pkg.activate(".")

## Add necessary packages to project
Pkg.add("Plots")
Pkg.add("DifferentialEquations")
Pkg.add("SparseArrays")
Pkg.add("DataFrames")
Pkg.add("CSV")

## check package status
Pkg.status()

## load packages
using Plots, DifferentialEquations, SparseArrays, DataFrames, CSV


text:
module BlockCoordinateDescent

export multinomialBCD

include("bcd.jl")

end # module


text:
function addnode!{V,W}(g::Graph{V,W}, u::V)
    if !haskey(g, u)
        g[u] = valtype(g)()
    end
    nothing
end

function deletenode!{V,W}(g::Graph{V,W}, u::V)
    for v in keys(g[u])
        delete!(g[v], u)
    end
    delete!(g, u)
end

function addedge!{V,W}(g::Graph{V,W}, u::V, v::V, ew::W)
    addnode!(g, u)
    addnode!(g, v)
    g[u][v] = ew
    g[v][u] = ew
end
addedge!{V,W}(g::Graph{V,W}, u::V, v::V) = addedge!(g, u, v, zero(W))

function deleteedge!{V,W}(g::Graph{V,W}, u::V, v::V)
    delete!(g[u], v)
    delete!(g[v], u)
end

function ne{V,W}(g::Graph{V,W})
	m = 0
	for u in keys(g)
		m += length(g[u])
	end
	return div(m, 2)
end


text:
# This file was generated, do not modify it. # hide
using MLJ, StatsBase, Random, PyPlot, CategoricalArrays, PrettyPrinting, DataFrames
X, y = @load_crabs
X = DataFrame(X)
@show size(X)
@show y[1:3]
first(X, 3) |> pretty

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using OptimalTransport

using Distances
using PythonOT: PythonOT

using LinearAlgebra
using Random
using Test

const POT = PythonOT

Random.seed!(100)

@testset "sinkhorn_barycenter.jl" begin
    # set up support
    support = range(-1; stop=1, length=500)
    N = 10
    μ = hcat([normalize!(exp.(-(support .+ rand()) .^ 2 ./ 0.1^2), 1) for _ in 1:N]...)

    # create cost matrix
    C = pairwise(SqEuclidean(), support'; dims=2)

    # regularisation parameter
    ε = 0.05

    # weights
    w = ones(N) / N

    @testset "example" begin
        α = sinkhorn_barycenter(μ, C, ε, w, SinkhornGibbs())

        # compare with POT
        # need to use a larger tolerance here because of a quirk with the POT solver
        α_pot = POT.barycenter(μ, C, ε; weights=w, stopThr=1e-9)
        @test α ≈ α_pot rtol = 1e-6
    end

    # different element type
    @testset "Float32" begin
        μ32 = map(Float32, μ)
        ε32 = map(Float32, ε)
        C32 = map(Float32, C)
        w32 = map(Float32, w)
        α = sinkhorn_barycenter(μ32, C32, ε32, w32, SinkhornGibbs())

        α_pot = POT.barycenter(μ32, C32, ε32; weights=w32, stopThr=1e-9)
        @test α ≈ α_pot rtol = 1e-6
    end
end


text:
function stats(id::String)
  res = HTTP.get(ENDPOINT * "products/$id/stats")
  rawData = JSON.parse(String(res.body))
  DataFrame(rawData)
end

function stats(ccy1::String, ccy2::String)
  stats(ccy1 * "-" * ccy2)
end


text:
# test FourierOperators
function testFT(N=16)
  # random image
  x = zeros(ComplexF64,N,N)
  for i=1:N,j=1:N
    x[i,j] = rand()
  end

  # FourierMatrix
  idx = CartesianIndices((N,N))[collect(1:N^2)]
  F = [ exp(-2*pi*im*((idx[j][1]-1)*(idx[k][1]-1)+(idx[j][2]-1)*(idx[k][2]-1))/N) for j=1:N^2, k=1:N^2 ]
  F_adj = F'

  # Operators
  tr = CartesianTrajectory(N,N)
  F_nfft = NFFTOp((N,N),tr,symmetrize=false)
  F_exp = ExplicitOp((N,N),tr,zeros(ComplexF64,N,N),symmetrize=false)

  # test agains FourierOperators
  y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N)) )
  y_adj = vec( ifftshift(reshape(F_adj*vec(fftshift(x)),N,N)) )

  y_nfft = F_nfft*vec(x)
  y_adj_nfft = adjoint(F_nfft) * vec(x)

  y_exp = F_exp*vec(x)
  y_adj_exp = adjoint(F_exp) * vec(x)

  @test (norm(y-y_nfft)/norm(y)) < 1e-2
  @test (norm(y_adj-y_adj_nfft)/norm(y_adj)) < 1e-2
  @test (norm(y-y_exp)/norm(y)) < 1e-2
  @test (norm(y_adj-y_adj_exp)/norm(y_adj)) < 1e-2
end

function testFT3d(N=12)
  # random image
  x = zeros(ComplexF64,N,N,N)
  for i=1:N,j=1:N,k=1:N
    x[i,j,k] = rand()
  end

  # FourierMatrix
  idx = CartesianIndices((N,N,N))[collect(1:N^3)]
  F = [ exp(-2*pi*im*((idx[j][1]-1)*(idx[k][1]-1)+(idx[j][2]-1)*(idx[k][2]-1)+(idx[j][3]-1)*(idx[k][3]-1))/N) for j=1:N^3, k=1:N^3 ]
  F_adj = F'

  # Operators
  tr = CartesianTrajectory3D(N,N,numSlices=N)
  F_nfft = NFFTOp((N,N,N),tr,symmetrize=false)
  F_exp = ExplicitOp((N,N,N),tr,zeros(ComplexF64,N,N,N),symmetrize=false)

  # test agains FourierOperators
  y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N,N)) )
  y_adj = vec( ifftshift(reshape(F_adj*vec(fftshift(x)),N,N,N)) )

  y_nfft = F_nfft*vec(x)
  y_adj_nfft = adjoint(F_nfft) * vec(x)

  y_exp = F_exp*vec(x)
  y_adj_exp = adjoint(F_exp) * vec(x)

  relError = (norm(y-y_exp)/norm(y))
  println("Relative error ExplicitOp: ", relError)
  @test  relError< 1e-2
  relError = (norm(y_adj-y_adj_exp)/norm(y_adj))
  println("Relative error adj. ExplicitOp: ", relError)
  @test  relError< 1e-2
  relError = (norm(y-y_nfft)/norm(y))
  println("Relative error NFFT: ", relError)
  @test relError < 1e-2
  relError = (norm(y_adj-y_adj_nfft)/norm(y_adj))
  println("Relative error adj. NFFT: ", relError)
  @test  relError< 1e-2
end

# test FieldmapNFFTOp
function testFieldmapFT(N=16)
  # random image
  x = zeros(ComplexF64,N,N)
  for i=1:N,j=1:N
    x[i,j] = rand()
  end

  tr = CartesianTrajectory(N,N;TE=0.0,AQ=0.01)
  times = readoutTimes(tr)
  nodes = kspaceNodes(tr)
  cmap = im*quadraticFieldmap(N,N)[:,:,1]

  # FourierMatrix
  idx = CartesianIndices((N,N))[collect(1:N^2)]
  F = [exp(-2*1im*pi*(nodes[1,k]*(idx[l][1]-size(x,1)/2-1)+nodes[2,k]*(idx[l][2]-size(x,2)/2-1))-cmap[idx[l][1],idx[l][2]]*times[k]) for k=1:size(nodes,2), l=1:length(x)]
  F_adj = F'

  # Operators
  F_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false)

  # test agains FourierOperators
  y = F*vec(x)
  y_adj = F_adj*vec(x)

  y_nfft = F_nfft*vec(x)
  y_adj_nfft = adjoint(F_nfft) * vec(x)

  @test (norm(y-y_nfft)/norm(y)) < 1e-2
  @test (norm(y_adj-y_adj_nfft)/norm(y_adj)) < 1e-2
end

function testOperators()
  @testset "Linear Operator" begin
    testFT()
    testFT3d()
    testFieldmapFT()
  end
end

testOperators()


text:
include("src/TransitRouter.jl")

using .TransitRouter
using .TransitRouter.OSRM
using ArgParse

function main()
    parser = ArgParseSettings()
    @add_arg_table! parser begin
        "--osrm-network", "-n"
            help = "path to a .osrm file to use for transfer finding"
        "--max-transfer-distance", "-t"
            help = "maximum transfer distance (meters)"
            arg_type = Float64
            default = 1000.0
        "output"
            help = "Output file to save network (ends in .trjl)"
            required = true
        "gtfs"
            help = "GTFS input files"
            nargs = '+'
            required = true
    end

    parsed_args = parse_args(parser)

    gtfs::Vector{String} = convert(Vector{String}, parsed_args["gtfs"])::Vector{String}
    output::String = parsed_args["output"]::String

    if haskey(parsed_args, "osrm-network")
        @info "Starting OSRM to route through the street network"
        # TODO don't hardwire mld
        osrm = start_osrm(parsed_args["osrm-network"]::String, "mld")
        network = build_network(gtfs, osrm, parsed_args["max-transfer-distance"])
        stop_osrm!(osrm)
        save_network(network, output)
    else
        network = build_network(gtfs, max_transfer_distance_meters=parsed_args["max-transfer-distance"])
        save_network(network, output)
    end
end

main()


text:
using Test, GeneticAlgorithm

@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 2) == [1, 2, 30, 40]
@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 1) == [1, 20, 30, 40]
@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 0) == [10, 20, 30, 40]
@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 4) == [1, 2, 3, 4]

include("runtests_example01.jl")
include("runtests_example02.jl")

msg = "!! ALL DONE !!"
printstyled("-" ^ length(msg), color=:green)
printstyled(msg, color=:green)
printstyled("-" ^ length(msg), color=:green)


text:
module Optimizers

import PyCall: PyObject, pycall

import ..Keras
import ..Keras: PyDoc

const keras_optimizers = [
    "SGD",
    "RMSprop",
    "Adagrad",
    "Adadelta",
    "Adam",
    "Adamax",
    "Nadam",
]

for o in keras_optimizers
    opt_name = Symbol(o)

    @eval begin
        struct $opt_name
            obj::PyObject

            @doc PyDoc(Keras._optimizers, Symbol($o)) function $opt_name(args...; kwargs...)
                new(Keras._optimizers[Symbol($o)](args...; kwargs...))
            end
        end

        # convert(::Type{$(opt_name)}, obj::PyObject) = $opt_name(obj)
        PyObject(opt::$(opt_name)) = opt.obj
    end
end

end


text:
# ---
# title: 1001. Grid Illumination
# id: problem1001
# author: Indigo
# date: 2021-01-30
# difficulty: Hard
# categories: Hash Table
# link: <https://leetcode.com/problems/grid-illumination/description/>
# hidden: true
# ---
#
# You are given a `grid` of size `N x N`, and each cell of this grid has a lamp
# that is initially **turned off**.
#
# You are also given an array of lamp positions `lamps`, where `lamps[i] =
# [rowi, coli]` indicates that the lamp at `grid[rowi][coli]` is **turned on**.
# When a lamp is turned on, it **illuminates its cell** and **all  other cells**
# in the same **row, column, or diagonal**.
#
# Finally, you are given a query array `queries`, where `queries[i] = [rowi,
# coli]`. For the `ith` query, determine whether `grid[rowi][coli]` is
# illuminated or not. After answering the `ith` query,  **turn off** the lamp at
# `grid[rowi][coli]` and its **8 adjacent lamps** if they exist. A lamp is
# adjacent if its cell shares either a side or corner with `grid[rowi][coli]`.
#
# Return _an array of integers_`ans` _,_ _  where _`ans[i]` _should be_`1` _  if
# the lamp in the _`ith` _query was illuminated, or  _`0` _  if the lamp was
# not._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/19/illu_1.jpg)
#
#
#
#     Input: N = 5, lamps = [[0,0],[4,4]], queries = [[1,1],[1,0]]
#     Output: [1,0]
#     Explanation: We have the initial grid with all lamps turned off. In the above picture we see the grid after turning on the lamp at grid[0][0] then turning on the lamp at grid[4][4].
#     The 0th query asks if the lamp at grid[1][1] is illuminated or not (the blue square). It is illuminated, so set ans[0] = 1. Then, we turn off all lamps in the red square.
#     ![](https://assets.leetcode.com/uploads/2020/08/19/illu_step1.jpg)
#     The 1st query asks if the lamp at grid[1][0] is illuminated or not (the blue square). It is not illuminated, so set ans[1] = 1. Then, we turn off all lamps in the red rectangle.
#     ![](https://assets.leetcode.com/uploads/2020/08/19/illu_step2.jpg)
#
#
# **Example 2:**
#
#
#
#     Input: N = 5, lamps = [[0,0],[4,4]], queries = [[1,1],[1,1]]
#     Output: [1,1]
#
#
# **Example 3:**
#
#
#
#     Input: N = 5, lamps = [[0,0],[0,4]], queries = [[0,4],[0,1],[1,4]]
#     Output: [1,1,0]
#
#
#
#
# **Constraints:**
#
#   * `1 <= N <= 109`
#   * `0 <= lamps.length <= 20000`
#   * `lamps[i].length == 2`
#   * `0 <= lamps[i][j] < N`
#   * `0 <= queries.length <= 20000`
#   * `queries[i].length == 2`
#   * `0 <= queries[i][j] < N`
#
#
## @lc code=start
using LeetCode

function grid_illumination(N, lamps::Vector{Tuple{Int,Int}}, query::Vector{Tuple{Int,Int}})
    function update!(lamp, q)
        col_illum[lamp[2]] += q
        row_illum[lamp[1]] += q
        main_diag_illum[lamp[2]-lamp[1]] += q
        d_diag_illum[lamp[2]+lamp[1]] += q
    end
    dirs = [
        -1 -1 -1 0 0 0 1 1 1
        1 0 -1 1 0 -1 1 0 -1
    ]
    col_illum = DefaultDict{Int,Int}(0)
    row_illum = DefaultDict{Int,Int}(0)
    main_diag_illum = DefaultDict{Int,Int}(0)
    d_diag_illum = DefaultDict{Int,Int}(0)
    illum_pos = Set(lamps)
    res = fill(0, length(query))
    for lamp in lamps
        update!(lamp, 1)
    end
    for (idx, q) in enumerate(query)
        if col_illum[q[2]] > 0 ||
           row_illum[q[1]] > 0 ||
           main_diag_illum[q[2]-q[1]] > 0 ||
           main_diag_illum[q[2]+q[1]] > 0
            res[idx] = 1
        end
        for dir in 1:9
            new_pos = (q[1] + dirs[1, dir], q[2] + dirs[2, dir])
            if new_pos in illum_pos
                pop!(illum_pos)
                update!(new_pos, -1)
            end
        end
    end
    res
end
## @lc code=end


text:
using SymbolicMDPs
using Test

@testset "SymbolicMDPs.jl" begin
    # Write your tests here.
end


text:
using Iris
using Test

@testset "Domains" begin
    @testset "1D" begin
        shape = Interval(-1,1)
        bnd = Boundary(shape)
        @testset "Nondispersive" begin
            ε = @test_nowarn NondispersiveDomain(shape,1+.1im,:test,:only)
            @test_nowarn NondispersiveDomain(shape,1)
            @test_nowarn NondispersiveDomain(shape,1im)
            @test_nowarn NondispersiveDomain(shape,1,.1)
            @test_nowarn NondispersiveDomain(shape,1,0)
            @test_nowarn NondispersiveDomain(shape,0,1)
            @test_nowarn NondispersiveDomain(bnd)
        end
        @testset "Dispersive" begin
            disp = TwoLevelSystem(3,.1,1)
            @test_nowarn DispersiveDomain(bnd,disp;type=:test,name=:only)
            @test_nowarn DispersiveDomain(bnd,disp,.5;type=:test,name=:only)
            ε = @test_nowarn DispersiveDomain(shape,disp;type=:test,name=:only)
        end
        @testset "Lattice" begin
            @testset "Cartesian" begin
                lat = Lattices.Cartesian(.1)
                dom = @test_nowarn LatticeDomain(bnd,lat)
                @test_nowarn LatticeDomain(lat,bnd)
                @test_nowarn LatticeDomain(bnd,lat,2)
                @test_nowarn LatticeDomain(lat,bnd,2)
                @test_nowarn LatticeDomain(bnd,2,lat)
                @test_nowarn LatticeDomain(lat,2,bnd)
                @test_nowarn LatticeDomain(2,bnd,lat)
                @test_nowarn LatticeDomain(2,lat,bnd)
            end
            # @testset "Polar" begin
            #     lat = Lattices.Cartesian(.1)
            #     dom = @test_nowarn LatticeDomain(bnd,lat)
            #     @test_nowarn LatticeDomain(lat,bnd)
            #     @test_nowarn LatticeDomain(bnd,lat,2)
            #     @test_nowarn LatticeDomain(lat,bnd,2)
            #     @test_nowarn LatticeDomain(bnd,2,lat)
            #     @test_nowarn LatticeDomain(lat,2,bnd)
            #     @test_nowarn LatticeDomain(2,bnd,lat)
            #     @test_nowarn LatticeDomain(2,lat,bnd)
            # end
            # @testset "Spherical" begin
            #     lat = Lattices.Cartesian(.1)
            #     dom = @test_nowarn LatticeDomain(bnd,lat)
            #     @test_nowarn LatticeDomain(lat,bnd)
            #     @test_nowarn LatticeDomain(bnd,lat,2)
            #     @test_nowarn LatticeDomain(lat,bnd,2)
            #     @test_nowarn LatticeDomain(bnd,2,lat)
            #     @test_nowarn LatticeDomain(lat,2,bnd)
            #     @test_nowarn LatticeDomain(2,bnd,lat)
            #     @test_nowarn LatticeDomain(2,lat,bnd)
            # end
        end
    end
    @testset "2D" begin
        shape = Square(1,0,0)
        bnd = Boundary(shape)
        @testset "Nondispersive" begin
            ε = @test_nowarn NondispersiveDomain(shape,1+.1im,:test,:only)
            @test_nowarn NondispersiveDomain(shape,1)
            @test_nowarn NondispersiveDomain(shape,1im)
            @test_nowarn NondispersiveDomain(shape,1,.1)
            @test_nowarn NondispersiveDomain(shape,1,0)
            @test_nowarn NondispersiveDomain(shape,0,1)
            @test_nowarn NondispersiveDomain(bnd)
        end
        @testset "Dispersive" begin
            disp = TwoLevelSystem(3,.1,1)
            @test_nowarn DispersiveDomain(bnd,disp;type=:test,name=:only)
            @test_nowarn DispersiveDomain(bnd,disp,.5;type=:test,name=:only)
            ε = @test_nowarn DispersiveDomain(shape,disp;type=:test,name=:only)
        end
        @testset "Lattice" begin
            @testset "Cartesian" begin
                lat = Lattices.Cartesian(.1,.1)
                dom = @test_nowarn LatticeDomain(bnd,lat)
                @test_nowarn LatticeDomain(lat,bnd)
                @test_nowarn LatticeDomain(bnd,lat,2)
                @test_nowarn LatticeDomain(lat,bnd,2)
                @test_nowarn LatticeDomain(bnd,2,lat)
                @test_nowarn LatticeDomain(lat,2,bnd)
                @test_nowarn LatticeDomain(2,bnd,lat)
                @test_nowarn LatticeDomain(2,lat,bnd)
                lat = Lattices.Cartesian(.05,.05,angle=.05)
                @test_nowarn typeof(LatticeDomain(bnd,lat))<:LatticeDomain{2,Domains.Unsymmetric}
            end
            # @testset "Polar" begin
            #     lat = Lattices.Cartesian(.1)
            #     dom = @test_nowarn LatticeDomain(bnd,lat)
            #     @test_nowarn LatticeDomain(lat,bnd)
            #     @test_nowarn LatticeDomain(bnd,lat,2)
            #     @test_nowarn LatticeDomain(lat,bnd,2)
            #     @test_nowarn LatticeDomain(bnd,2,lat)
            #     @test_nowarn LatticeDomain(lat,2,bnd)
            #     @test_nowarn LatticeDomain(2,bnd,lat)
            #     @test_nowarn LatticeDomain(2,lat,bnd)
            # end
            # @testset "Spherical" begin
            #     lat = Lattices.Cartesian(.1)
            #     dom = @test_nowarn LatticeDomain(bnd,lat)
            #     @test_nowarn LatticeDomain(lat,bnd)
            #     @test_nowarn LatticeDomain(bnd,lat,2)
            #     @test_nowarn LatticeDomain(lat,bnd,2)
            #     @test_nowarn LatticeDomain(bnd,2,lat)
            #     @test_nowarn LatticeDomain(lat,2,bnd)
            #     @test_nowarn LatticeDomain(2,bnd,lat)
            #     @test_nowarn LatticeDomain(2,lat,bnd)
            # end
        end
    end
end


text:
#    Unit test.jl, part of PosDefManifold Package for julia language
#
#    MIT License
#    Copyright (c) 2019-21, Marco Congedo, CNRS, Grenobe, France:
#    https://sites.google.com/site/marcocongedo/home
#
#    DESCRIPTION
#    This Unit tests all functions in PosDefManifold.
#    Some functions are fully tested, the others are just executed.
#    Unce you ran it, for each method of each function,
#    a ⭐ sign is printed if the test is succesful, while
#    a ⛔ sign is printed if the test is not succesful.
#    a ☆ sign is printed if the function has been executed correctly.
#    If there are fails, the concerned functions will be listed as Warnings
#    and returned by the testall() function as an array of strings

#    RUN TESTS
#    just hit CRTL+A (seect all) and Shift+Enter (run),
#    then invoke 'testall()' in the REPL

#    USE IT FOR TESTING YOUR OWN PACKAGES
#    - Copy this unit in the same directory where the Main
#      module of your library is.
#    - Replace 'PosDefManifold' by the name of your module(s) in the
#      'using' command here below.
#    - Overwrite the tests() function with your own tests.

push!(LOAD_PATH, pwd())
using LinearAlgebra, Statistics, PosDefManifold

function tests();
    metrics=[Metric(i) for i in 1:10]
    n=10
    m=20
    P=randP(n)      # an SPD Hermitian matrix
    Q=randP(n)      # an SPD Hermitian matrix
    PC=randP(ComplexF64, n) # a complex Hermitian matrix
    QC=randP(ComplexF64, n) # a complex Hermitian matrix
    T=randn(m, n)   # a real tall Matrix
    TC=randn(ComplexF64, m, n)   # a complex tall Matrix
    T2=randn(m, n)   # a real tall Matrix
    TC2=randn(ComplexF64, m, n)   # a complex tall Matrix
    W=randn(n, m)   # a wide Matrix
    X=randn(n, n)   # a square matrix
    XC=randn(ComplexF64, n, n) # a complex square matrix
    λ=[0.0, 0.1, 0.2, 0.3] # a real vector
    λc=[0.0+0.0im, 0.1+0.1im, 0.2+0.2im, 0.3+0.3im] # a complex vector
    Λ=Diagonal(λ)   # a real Diagonal matrix
    ΛC=Diagonal(λc) # a complex Diagonal matrix

    P_=ℍ([  3. 2. 1.;
            2. 4. 2.;
            1. 2. 5.])
    PC_=ℍ([3.0+0.0im 2.0-0.2im 1.0-0.1im;
            2.0+0.2im 4.0+0.0im 2.0-0.2im;
            1.0+0.1im 2.0+0.2im 5.0+0.0im])

    Q_=ℍ(P_+[ 2. 1. 0.;
            1. 3. 1.;
            0. 1. 4.])
    QC_=ℍ(PC_+[2.0+0.0im 1.0-0.2im 0.0-0.1im;
            1.0+0.2im 3.0+0.0im 1.0-0.2im;
            0.0+0.1im 1.0+0.2im 4.0+0.0im])

    D_=Diagonal(P_)
    DC_=Diagonal(PC_)
    E_=Diagonal(Q_)
    EC_=Diagonal(QC_)

    PC_small=randP(ComplexF64, 3)

    A3x2=[4. 3.; 2. 5.; 1. 2.]

    𝐏=randP(8, 4)
    𝐏C=randP(ComplexF64, 8, 4)
    𝐐=randP(8, 4)
    𝐐C=randP(ComplexF64, 8, 4)
    𝐃=randΛ(8, 4)
    𝐄=randΛ(8, 4)
    weights=[0.1, 0.2, 0.3, 0.4]

    # functions in LinearAlgebrainP.jl
    print("- Unit 'linearAlgebra.jl'")

    ## 1. Utilities
    name="typeofMatrix"; newTest(name)
    typeofMatrix(𝐏); RUN()

    name="dim"; newTest(name)
    dim(𝐏); RUN()

    ## 2. Matrix Normalizations and approximations

    name="det1"; newTest(name)
    det(det1(P))  ≈ 1 ?  OK() : OH(name*" real case")
    det(det1(PC)) ≈ 1 ?  OK() : OH(name*" complex case")

    name="tr1"; newTest(name)
    tr(tr1(P))  ≈ 1 ?    OK() : OH(name*" real case")
    tr(tr1(PC)) ≈ 1 ?    OK() : OH(name*" complex case")

    name="nearestPosDef"; newTest(name)
    nearestPosDef(X); RUN()
    nearestPosDef(XC); RUN()

    name="normalizeCol!"; newTest(name)
    j=rand(1:n)
    normalizeCol!(T, j)
    norm(T[:, j])       ≈ 1 ? OK() : OH(name*" Method 1 real case")
    normalizeCol!(T, j, 2)
    norm(T[:, j])       ≈ 0.5 ? OK() : OH(name*" Method 2 real case")
    normalizeCol!(T, 2:3)
    norm(T[:, 2])+norm(T[:, 3]) ≈ 2 ? OK() : OH(name*" Method 3 real case")
    normalizeCol!(T, 2:3, 2)
    norm(T[:, 2])+norm(T[:, 3]) ≈ 1 ? OK() : OH(name*" Method 4 real case")
    normalizeCol!(TC, j)
    norm(TC[:, j])       ≈ 1 ? OK() : OH(name*" Method 1 complex case")
    normalizeCol!(TC, j, 2)
    norm(TC[:, j])       ≈ 0.5 ? OK() : OH(name*" Method 2 complex case")
    normalizeCol!(TC, 2:3)
    norm(TC[:, 2])+norm(TC[:, 3]) ≈ 2 ? OK() : OH(name*" Method 3 complex case")
    normalizeCol!(TC, 2:3, 2)
    norm(TC[:, 2])+norm(TC[:, 3]) ≈ 1 ? OK() : OH(name*" Method 4 complex case")

    ## 3. Boolean functions of matrices

    name="ispos"; newTest(name)
    ispos(λ, 🔔=false) == false ? OK() : OH(name*" Method 1 real case")
    ispos(Λ, 🔔=false) == false ? OK() : OH(name*" Method 2 real case")

    ## 4. Scalar Functions of Matrices

    name="colProd"; newTest(name)
    j1=1; j2=rand(2:n);
    s=𝚺(T[:, j1].*T[:, j2])
    colProd(T, j1, j2) ≈ s ? OK() : OH(name*" Method 1 real case")
    s=𝚺(conj(TC[:, j1]).*TC[:, j2])
    colProd(TC, j1, j2) ≈ s ? OK() : OH(name*" Method 1 complex case")
    s=𝚺(T[:, j1].*T2[:, j2])
    colProd(T, T2, j1, j2) ≈ s ? OK() : OH(name*" Method 2 real case")
    s=𝚺(conj(TC[:, j1]).*TC2[:, j2])
    colProd(TC, TC2, j1, j2) ≈ s ? OK() : OH(name*" Method 2 complex case")

    name="colNorm"; newTest(name)
    colNorm(X, 2)      ≈ norm(X[:, 2]) ? OK()  : OH(name*" Method 1 real case")
    colNorm(XC, 2)     ≈ norm(XC[:, 2]) ? OK() : OH(name*" Method 1 complex case")

    name="sumOfSqr"; newTest(name)
    sumOfSqr(Matrix(P_))  ≈ 68    ? OK() : OH(name*" Method 1 real case")
    sumOfSqr(P_)          ≈ 68    ? OK() : OH(name*" Method 2 real case")
    sumOfSqr(LowerTriangular(P_)) ≈ 59 ? OK() : OH(name*" Method 3 real case")
    sumOfSqr(Diagonal(P_)) ≈ 50 ? OK() : OH(name*" Method 4 real case")
    sumOfSqr(P_, 2)       ≈ 24    ? OK() : OH(name*" Method 5 real case")
    sumOfSqr(P_, 1:2)     ≈ 38    ? OK() : OH(name*" Method 6 real case")
    sumOfSqr(Matrix(PC_)) ≈ 68.18 ? OK() : OH(name*" Method 1 complex case")
    sumOfSqr(PC_)         ≈ 68.18 ? OK() : OH(name*" Method 2 complex case")
    sumOfSqr(LowerTriangular(PC_)) ≈ 59.089999999999996 ? OK() : OH(name*" Method 3 complex case")
    sumOfSqr(Diagonal(PC_)) ≈ 50 ? OK() : OH(name*" Method 4 complex case")
    sumOfSqr(PC_, 2)      ≈ 24.08 ? OK() : OH(name*" Method 5 complex case")
    sumOfSqr(PC_, 1:3)    ≈ 68.18 ? OK() : OH(name*" Method 6 complex case")

    name="sumOfSqrDiag"; newTest(name)
    sumOfSqrDiag(P_)    ≈ 50 ? OK() : OH(name*" Method 1 real case")
    A=randn(ComplexF64, 3, 3)
    s=sum(abs2(A[i, i]) for i in 1:size(A, 1))
    sumOfSqrDiag(A)   ≈ s  ? OK() : OH(name*" Method 1 complex case")
    A=real(A)
    s=sum(A[i, i]^2 for i in 1:size(A, 1))
    sumOfSqrDiag(A)     ≈ s  ? OK() : OH(name*" Method 2")

    name="sumOfSqrTril"; newTest(name)
    sumOfSqrTril(A3x2, -1)≈ 9 ? OK() : OH(name*" Method 1 real case")
    A=randn(ComplexF64, 3, 2)
    s=sumOfSqr(A)-abs2(A[1, 2])
    sumOfSqrTril(A, 0)≈ s ? OK() : OH(name*" Method 1 complex case")

    name="tr"; newTest(name)
    tr(P, Q) ≈ tr(P*Q) ? OK() : OH(name*" Method 1 real case")
    tr(PC, QC) ≈ tr(PC*QC) ? OK() : OH(name*" Method 1 complex case")
    tr(P, X) ≈ tr(P*X) ? OK() : OH(name*" Method 2 real case")
    tr(PC, XC) ≈ tr(PC*XC) ? OK() : OH(name*" Method 2 complex case")
    tr(D_, P_) ≈ tr(D_*P_) ? OK() : OH(name*" Method 3 real case")
    tr(DC_, PC_) ≈ tr(DC_*PC_) ? OK() : OH(name*" Method 3 complex case")
    tr(P_, D_) ≈ tr(D_*P_) ? OK() : OH(name*" Method 4 real case")
    tr(PC_, DC_) ≈ tr(DC_*PC_) ? OK() : OH(name*" Method 4 complex case")

    name="quadraticForm"; newTest(name)
    v=randn(n)
    vC=randn(ComplexF64, n)
    quadraticForm(v, P) ≈ v'*P*v ? OK() : OH(name*" Method 1 real case")
    quadraticForm(v, LowerTriangular(Matrix(P))) ≈ v'*LowerTriangular(Matrix(P))*v ? OK() : OH(name*" Method 2 real case")
    quadraticForm(v, Matrix(P), true) ≈ v'*P*v ? OK() : OH(name*" Method 3a real case")
    quadraticForm(v, Matrix(P), false) ≈ v'*P*v ? OK() : OH(name*" Method 3b real case")
    quadraticForm(vC, PC) ≈ vC'*PC*vC ? OK() : OH(name*" Method 4a complex case")
    quadraticForm(vC, Matrix(PC)) ≈ vC'*PC*vC ? OK() : OH(name*" Method 4b complex case")
    quadraticForm(vC, LowerTriangular(Matrix(PC))) ≈ vC'*LowerTriangular(Matrix(PC))*vC ? OK() : OH(name*" Method 4c complex case")

    name="fidelity"; newTest(name);
    # Test compilation only
    f=fidelity(P, Q); RUN()
    f=fidelity(PC, QC); RUN()

    ## 5. Diagonal functions of matrices

    name="fDiag"; newTest(name)
    D=fDiag(x->x^2, P_)
    D≈Diagonal(diagm(0 => [9.,16.,25.])) ? OK() : OH(name)

    name="DiagOfProd"; newTest(name)
    DiagOfProd(P, Q)≈Diagonal(P*Q) ? OK() : OH(name*" Real Input")
    DiagOfProd(PC, QC)≈Diagonal(PC*QC) ? OK() : OH(name*" Complex Input")

    ## 6. Unitary functions of matrices

    name="mgs"; newTest(name)
    U=mgs(X)
    U'*U≈I ? OK() : OH(name*" Real Input")
    U=mgs(XC)
    U'*U≈I ? OK() : OH(name*" Complex Input")

    ## 7. Matrix function of matrices

    name="fVec"; newTest(name)
    Pset=randP(4, 64)
    mean(Pset)≈fVec(mean, Pset) ? OK() : OH(name*" simple mean, real")
    mean(sqrt, Pset)≈fVec(mean, sqrt, Pset) ? OK() : OH(name*" mean, sqrt, real")
    fVecw=(randn(64)).^2
    fVecw=fVecw./sum(fVecw) # generate normalized random weights
    mean(Pset[i]*fVecw[i] for i=1:64)≈fVec(mean, Pset; w=fVecw) ? OK() : OH(name*" simple weighted mean, real")
    mean(sqrt(Pset[i])*fVecw[i] for i=1:64)≈fVec(mean, sqrt, Pset; w=fVecw) ? OK() : OH(name*" weighted mean, sqrt, real")

    name="congruence (I)"; newTest(name)
    M=randn(n, n)
    congruence(M, P, ℍ)≈M*P*M' ? OK() : OH(name)

    name="congruence (II)"; newTest(name)
    Pset=randP(n, n*4)
    Qset=cong(M, Pset, ℍVector)
    Qset≈[M*Pset[i]*M' for i=1:n*4] ? OK() : OH(name)

    name="congruence (III)"; newTest(name)
    Pset1=randP(4, 10); # generate 10 positive definite 4x4 matrices
    Pset2=randP(4, 8);
    Pset=ℍVector₂([Pset1, Pset2]);
    M=randn(4, 4)
    Qset=cong(M, Pset, MatrixVector₂)
    Qset[1][1]≈M*Pset[1][1]*M' &&
    Qset[1][5]≈M*Pset[1][5]*M' &&
    Qset[2][1]≈M*Pset[2][1]*M' &&
    Qset[2][4]≈M*Pset[2][4]*M' ? OK() : OH(name)

    name="congruence (IV)"; newTest(name)
    Pset1=randP(4, 2); # generate 2 positive definite 4x4 matrices
    Pset2=randP(4, 2);
    Pset=ℍVector₂([Pset1, Pset2]);
    U=𝕄Vector([randU(4), randU(4)])
    Qset=cong(U, Pset, MatrixVector₂)
    Qset[1][1]≈U[1]*Pset[1][1]*U[1]' &&
    Qset[1][2]≈U[1]*Pset[1][2]*U[2]' &&
    Qset[2][1]≈U[2]*Pset[2][1]*U[1]' &&
    Qset[2][2]≈U[2]*Pset[2][2]*U[2]' ? OK() : OH(name)

    ## 8. Spectral decompositions of positive matrices

    name="evd"; newTest(name)
    (Λ, U) = evd(P)
    U*Λ*U'≈P ? OK() : OH(name*" Real Input")
    (Λ, U) = evd(PC)
    U*Λ*U'≈PC ? OK() : OH(name*" Complex Input")

    name="frf"; newTest(name)
    F = frf(P)
    F*F'≈P ? OK() : OH(name*" Real Input")
    F = frf(PC)
    F*F'≈PC ? OK() : OH(name*" Complex Input")

    name="invfrf"; newTest(name)
    F = invfrf(P)
    F*P*F'≈I ? OK() : OH(name*" Real Input")
    F = invfrf(PC)
    F*PC*F'≈I ? OK() : OH(name*" Complex Input")

    name="spectralFunctions"; newTest(name);
    spectralFunctions(P, x->x+1); RUN()
    spectralFunctions(PC, abs2); RUN()
    spectralFunctions(D_, x->x+1); RUN()

    name="pow"; newTest(name)
    P½, P½ⁱ=pow(P_, 0.5, -0.5)
    P½*P½ⁱ≈I && P½*P½≈P_ ?  OK() : OH(name*" Real Input method 1")
    P½, P½ⁱ=pow(PC_, 0.5, -0.5)
    P½*P½ⁱ≈I && P½*P½≈PC_ ? OK() : OH(name*" Complex Input method 1")
    D½, D½ⁱ=pow(D_, 0.5, -0.5)
    D½*D½ⁱ≈I && D½*D½≈D_ ?  OK() : OH(name*" Real Input method 2")

    name="invsqrt"; newTest(name)
    P½ⁱ=invsqrt(P_)
    P½ⁱ*P_*P½ⁱ'≈I ? OK() : OH(name*" Real Input method 1")
    P½ⁱ=invsqrt(PC_)
    P½ⁱ*PC_*P½ⁱ'≈I ? OK() : OH(name*" Complex Input method 1")
    D½ⁱ=invsqrt(D_)
    D½ⁱ*D_*D½ⁱ'≈I ? OK() : OH(name*" Real Input method 2")

    name="sqr"; newTest(name)
    P²=sqr(P_)
    P² ≈ P_*P_' ? OK() : OH(name*" Real Input method 1")
    P²=sqr(PC_)
    P² ≈ PC_*PC_' ? OK() : OH(name*" Complex Input method 1")
    D²=sqr(D_)
    D² ≈ D_*D_' ? OK() : OH(name*" Real Input method 2")

    name="powerIterations"; newTest(name)
    Λ, U, iterations, convergence=powIter(P_, size(P_, 2), evalues=true)
    sort(diag(Λ))≈eigvals(P_) && U*Λ*U'≈P_ ? OK() : OH(name*" Real Input Method 1")
    Λ, U, iterations, convergence=powIter(PC_, size(PC_, 2), evalues=true)
    sort(diag(Λ))≈eigvals(PC_) && U*Λ*U'≈PC_ ? OK() : OH(name*" Complex Input Method 1")
    Λ, U, iterations, convergence=powIter(ℍ(P_), size(P_, 2), evalues=true)
    sort(diag(Λ))≈eigvals(P_) && U*Λ*U'≈P_ ? OK() : OH(name*" Real Input Method 2")
    Λ, U, iterations, convergence=powIter(ℍ(PC_), size(PC_, 2), evalues=true)
    sort(diag(Λ))≈eigvals(PC_) && U*Λ*U'≈PC_ ? OK() : OH(name*" Complex Input Method 2")
    Λ, U, iterations, convergence=powIter(LowerTriangular(P_), size(P_, 2), evalues=true)
    sort(diag(Λ))≈eigvals(P_) && U*Λ*U'≈P_ ? OK() : OH(name*" Real Input Method 3")

    name="choL"; newTest(name)
    L=choL(P_)
    L*L'≈P_ ? OK() : OH(name*" Real Input method 1")
    L=choL(PC_)
    L*L'≈PC_ ? OK() : OH(name*" Complex Input method 2")
    L=choL(D_)
    L*L'≈D_ ? OK() : OH(name*" Real Input method 2")

    name="choInv"; newTest(name)
    L, Linv=choInv(P_)
    invP_=inv(P_)
    L*L'≈P_ ? OK() : OH(name*" Real Input test L")
    Linv*Linv'≈invP_ ? OK() : OH(name*" Real Input test Linv")
    L, D, Linv=choInv(P_; kind=:LDLt)
    L*D*L'≈P_ ? OK() : OH(name*" Real Input, kind :LDLt test L")
    Linv*inv(D)*Linv'≈invP_ ? OK() : OH(name*" Real Input, kind :LDLt test Linv")

    L, Linv=choInv(PC_)
    invPC_=inv(PC_)
    L*L'≈PC_ ? OK() : OH(name*" Complex Input test L")
    Linv*Linv'≈invPC_ ? OK() : OH(name*" Complex Input test Linv")
    L, D, Linv=choInv(PC_; kind=:LDLt)
    L*D*L'≈PC_ ? OK() : OH(name*" Complex Input, kind :LDLt test L")
    Linv*inv(D)*Linv'≈invPC_ ? OK() : OH(name*" Complex Input, kind :LDLt test Linv")

    # 9. Decompositions involving triangular matrices

    # functions in SignalProcessinginP.jl
    println(" ")
    print("\n- Unit 'signalProcessing.jl'")

    name="randλ"; newTest(name);
    randλ(10); RUN()

    name="randΛ"; newTest(name);
    randΛ(10); RUN()
    randΛ(10, 2); RUN()

    name="randU"; newTest(name);
    U=randU(10)
    U'*U≈I ? OK() : OH(name*" Real Input")
    U=randU(ComplexF64, 10)
    U'*U≈I ? OK() : OH(name*" Complex Input")

    name="randP"; newTest(name);
    randP(3); RUN()
    randP(ComplexF64, 3); RUN()

    name="regularize!"; newTest(name)
    signalVar=tr(P)
    regularize!(P, SNR=10)
    signalPlusNoiseVar=tr(P)
    signalVar/(signalPlusNoiseVar-signalVar) ≈ 10 ? OK() : OH(name*" Real Input Method 1")

    signalVar=tr(PC)
    regularize!(PC, SNR=10)
    signalPlusNoiseVar=tr(PC)
    signalVar/(signalPlusNoiseVar-signalVar) ≈ 10 ? OK() : OH(name*" Complex Input Method 1")

    𝐏2=randP(5, 20)
    signalVar=𝚺(tr(P) for P in 𝐏2)
    regularize!(𝐏2, SNR=10)
    signalPlusNoiseVar=𝚺(tr(P) for P in 𝐏2)
    signalVar/(signalPlusNoiseVar-signalVar) ≈ 10 ? OK() : OH(name*" Real Input Method 2")

    𝐏C2=randP(ComplexF64, 5, 20)
    signalVar=𝚺(tr(P) for P in 𝐏C2)
    regularize!(𝐏C2, SNR=10)
    signalPlusNoiseVar=𝚺(tr(P) for P in 𝐏C2)
    signalVar/(signalPlusNoiseVar-signalVar) ≈ 10 ? OK() : OH(name*" Real Input Method 2")

    name="gram"; newTest(name);
    gram(T); RUN()
    gram(W); RUN()

    name="trade"; newTest(name);
    trade(P); RUN()
    trade(PC); RUN()

    # functions in RiemannianGeometryinP.jl
    println(" ")
    print("\n- Unit 'riemannianGeometry.jl'")

    name="geodesic"; newTest(name);
    (geodesic(m, P, Q, 0.5) for m in metrics if m≠9); RUN()
    (geodesic(m, PC, QC, 0.5) for m in metrics if m≠9); RUN()
    (geodesic(m, 𝐃[1], 𝐃[2], 0.5) for m in metrics if m≠9); RUN()

    name="distanceSqr (I)"; newTest(name);
    for m in metrics distanceSqr(m, P) end; RUN()
    for m in metrics distanceSqr(m, P, Q) end; RUN()
    for m in metrics distanceSqr(m, PC) end; RUN()
    for m in metrics distanceSqr(m, PC, QC) end; RUN()

    name="distanceSqr (II)"; newTest(name);
    for m in metrics
         distanceSqr(m, D_)≈distanceSqr(m, ℍ(𝕄(D_))) ? OK() : OH(name*", metric "*string(m))
    end

    name="distanceSqr (III)"; newTest(name);
    for m in metrics
         distanceSqr(m, D_, E_)≈distanceSqr(m, ℍ(𝕄(D_)), ℍ(𝕄(E_))) ? OK() : OH(name*", metric "*string(m))
    end

    name="distance (I)"; newTest(name);
    # since it calls distanceSqr just check the call works
    for m in metrics distance(m, P) end; RUN()
    for m in metrics distance(m, P, Q) end; RUN()
    for m in metrics distance(m, PC) end; RUN()
    for m in metrics distance(m, PC, QC); end; RUN()

    name="distance (II)"; newTest(name);
    for m in metrics distance(m, D_) end; RUN()
    for m in metrics distance(m, D_, E_) end; RUN()

    𝐏=randP(8, 4)
    𝐏C=randP(ComplexF64, 8, 4)
    k=length(𝐏)
    kC=length(𝐏C)

    name="distanceSqrMat (I)"; newTest(name);
    for m in metrics
            L1=distanceSqrMat(Float64, m, 𝐏, ⏩=false)
            manualL1=𝕃{Float64}(diagm(0 => zeros(k)))
            for j=1:k-1, i=j+1:k manualL1[i, j]=distanceSqr(m, 𝐏[i], 𝐏[j]) end
            manualL1≈L1 ? OK() : OH(name*" Real Input, metric "*string(m))
    end

    name="distanceSqrMat (I ⏩ )"; newTest(name);
    for m in metrics
            L2=distanceSqrMat(Float64, m, 𝐏)
            manualL2=𝕃{Float64}(diagm(0 => zeros(k)))
            for j=1:k-1, i=j+1:k manualL2[i, j]=distanceSqr(m, 𝐏[i], 𝐏[j]) end
            manualL2≈L2 ? OK() : OH(name*" Real Input, metric "*string(m))
    end

    name="distanceSqrMat (II)"; newTest(name);
    for m in metrics
            L3=distanceSqrMat(Float64, m, 𝐏C, ⏩=false)
            manualL3=𝕃{Float64}(diagm(0 => zeros(kC)))
            for j=1:kC-1, i=j+1:kC manualL3[i, j]=distanceSqr(m, 𝐏C[i], 𝐏C[j]) end
            manualL3≈L3 ? OK() : OH(name*" Complex Input, metric "*string(m))
    end

    name="distanceSqrMat (II ⏩ )"; newTest(name);
    for m in metrics
            L4=distanceSqrMat(Float64, m, 𝐏C)
            manualL4=𝕃{Float64}(diagm(0 => zeros(kC)))
            for j=1:kC-1, i=j+1:kC manualL4[i, j]=distanceSqr(m, 𝐏C[i], 𝐏C[j]) end
            manualL4≈L4 ? OK() : OH(name*" Complex Input, metric "*string(m))
    end

    name="distanceMat (I)"; newTest(name);
    for m in metrics
            L5=distanceMat(Float64, m, 𝐏, ⏩=false)
            manualL5=𝕃{Float64}(diagm(0 => zeros(k)))
            for j=1:k-1, i=j+1:k manualL5[i, j]=distance(m, 𝐏[i], 𝐏[j]) end
            manualL5≈L5 ? OK() : OH(name*" Real Input, metric "*string(m))
    end

    name="distanceMat (I ⏩ )"; newTest(name);
    for m in metrics
            L6=distanceMat(Float64, m, 𝐏)
            manualL6=𝕃{Float64}(diagm(0 => zeros(k)))
            for j=1:k-1, i=j+1:k manualL6[i, j]=distance(m, 𝐏[i], 𝐏[j]) end
            manualL6≈L6 ? OK() : OH(name*" Real Input, metric "*string(m))
    end

    name="distanceMat (II)"; newTest(name);
    for m in metrics
            L7=distanceMat(Float64, m, 𝐏C, ⏩=false)
            manualL7=𝕃{Float64}(diagm(0 => zeros(kC)))
            for j=1:kC-1, i=j+1:kC manualL7[i, j]=distance(m, 𝐏C[i], 𝐏C[j]) end
            manualL7≈L7 ? OK() : OH(name*" Complex Input, metric "*string(m))
    end

    name="distanceMat (II ⏩ )"; newTest(name);
    for m in metrics
            L8=distanceMat(Float64, m, 𝐏C)
            manualL8=𝕃{Float64}(diagm(0 => zeros(kC)))
            for j=1:kC-1, i=j+1:kC manualL8[i, j]=distance(m, 𝐏C[i], 𝐏C[j]) end
            manualL8≈L8 ? OK() : OH(name*" Complex Input, metric "*string(m))
    end

    name="laplacian"; newTest(name);
    Dsqr=distanceSqrMat(logEuclidean, 𝐏)
    lap=laplacian(Dsqr); RUN()

    name="laplacianEigenMaps"; newTest(name);
    laplacianEM(lap, 2); RUN()

    name="spectralEmbedding"; newTest(name);
    spectralEmbedding(logEuclidean, 𝐏, 2); RUN()

    name="mean (I)"; newTest(name);
    for m=1:length(metrics)
            if m ∉ (7, 9) mean(metrics[m], P, Q) end end; RUN()
    for m=1:length(metrics)
            if m ∉ (6, 7, 9, 10) mean(metrics[m], 𝐏) end end; RUN()
    for m=1:length(metrics)
            if m ∉ (7, 9) mean(metrics[m], PC, QC) end end; RUN()
    for m=1:length(metrics)
            if m ∉ (6, 7, 9, 10) mean(metrics[m], 𝐏C) end end; RUN()
    for m=1:length(metrics)
            if m ∉ (6, 7, 9, 10) mean(metrics[m], 𝐃) end end; RUN()

    name="mean (II)"; newTest(name);
    k=length(𝐃)
    for m=1:length(metrics)
        if m ∉ (6, 7, 9, 10)
            D1=mean(metrics[m], 𝐃, ⏩=false)
            𝐃H=Vector{Hermitian}(undef, k)
            for i=1:k 𝐃H[i]=Hermitian(Matrix(𝐃[i])) end
            D2=mean(metrics[m], 𝐃H, ⏩=false)
            norm(𝕄(D1)-𝕄(D2))/k<0.0001 ? OK() : OH(name*" Real Diagonal Input, metric "*string(m))
        end
    end

    name="mean (⏩ )"; newTest(name);
    for m=1:length(metrics)
            if m ∉ (6, 7, 9, 10) mean(metrics[m], 𝐏) end end; RUN()
    for m=1:length(metrics)
            if m ∉ (7, 9) mean(metrics[m], 𝐃) end end; RUN()

    name="means"; newTest(name);
    means(logEuclidean, ℍVector₂([𝐏, 𝐐]); ⏩=false); RUN()
    means(logEuclidean, ℍVector₂([𝐏C, 𝐐C]); ⏩=false); RUN()
    means(logEuclidean, 𝔻Vector₂([𝐃, 𝐄]); ⏩=false); RUN()

    name="means (⏩ )"; newTest(name);
    means(logEuclidean, ℍVector₂([𝐏, 𝐐])); RUN()
    means(logEuclidean, ℍVector₂([𝐏C, 𝐐C])); RUN()
    means(logEuclidean, 𝔻Vector₂([𝐃, 𝐄])); RUN()

    name="generalizedMean"; newTest(name);
    𝐏2=ℍVector([P_, Q_])
    w=[0.2, 0.8]
    p=0.5
    ℍ( (P_^p+Q_^p)/2) ^(1/p) ≈ generalizedMean(𝐏2, p) ? OK() : OH(name*" Real Input 1")
    ℍ( (ℍ(0.2*P_^p)+ℍ(0.8*Q_^p))  )^(1/p) ≈ generalizedMean(𝐏2, p; w=w, ✓w=false) ? OK() : OH(name*" Real Input 2")
    w=w.*2.0
    ℍ( (ℍ(0.2*P_^p)+ℍ(0.8*Q_^p))  )^(1/p) ≈ generalizedMean(𝐏2, p; w=w) ? OK() : OH(name*" Real Input 3")
    ℍ( (ℍ(0.2*P_^p)+ℍ(0.8*Q_^p))  )^(1/p) ≉ generalizedMean(𝐏2, p; w=w, ✓w=false) ? OK() : OH(name*" Real Input 4")
    ℍ( (ℍ(0.4*P_^p)+ℍ(1.6*Q_^p))  )^(1/p) ≈ generalizedMean(𝐏2, p; w=w, ✓w=false) ? OK() : OH(name*" Real Input 5")
    𝐏2=ℍVector([PC_, QC_])
    w=[0.2, 0.8]
    ℍ( (PC_^p+QC_^p)/2) ^(1/p) ≈ generalizedMean(𝐏2, p) ? OK() : OH(name*" Complex Input 1")
    ℍ( (ℍ(0.2*PC_^p)+ℍ(0.8*QC_^p))  )^(1/p) ≈ generalizedMean(𝐏2, p; w=w, ✓w=false) ? OK() : OH(name*" Complex Input 2")
    w=w.*2.0
    ℍ( (ℍ(0.2*PC_^p)+ℍ(0.8*QC_^p))  )^(1/p) ≈ generalizedMean(𝐏2, p; w=w) ? OK() : OH(name*" Complex Input 3")
    ℍ( (ℍ(0.2*PC_^p)+ℍ(0.8*QC_^p))  )^(1/p) ≉ generalizedMean(𝐏2, p; w=w, ✓w=false) ? OK() : OH(name*" Complex Input 4")
    ℍ( (ℍ(0.4*PC_^p)+ℍ(1.6*QC_^p))  )^(1/p) ≈ generalizedMean(𝐏2, p; w=w, ✓w=false) ? OK() : OH(name*" Complex Input 5")
    ((𝐃[1]^p+𝐃[2]^p)/2)^(1/p) ≈ generalizedMean(𝔻Vector([𝐃[1], 𝐃[2]]), p) ? OK() : OH(name*" Real Diagonal Input")

    name="generalizedMean"; newTest(name);
    generalizedMean(𝐏, 0.5; ⏩=false); RUN()
    generalizedMean(𝐏, 0.5; w=weights, ✓w=false, ⏩=false); RUN()
    generalizedMean(𝐏C, 0.5; ⏩=false); RUN()
    generalizedMean(𝐏C, 0.5; w=weights, ✓w=false, ⏩=false); RUN()
    generalizedMean(𝐃, 0.5; ⏩=false); RUN()
    generalizedMean(𝐃, 0.5; w=weights, ✓w=false, ⏩=false); RUN()

    name="geometricMean(⏩ )"; newTest(name);
    geometricMean(𝐏); RUN()
    geometricMean(𝐏, w=weights, ✓w=false);  RUN()
    geometricMean(𝐏C); RUN()
    geometricMean(𝐏C, w=weights, ✓w=false); RUN()
    geometricMean(𝐃); RUN()
    geometricMean(𝐃, w=weights, ✓w=false); RUN()

    name="geometricMean"; newTest(name);
    geometricMean(𝐏; ⏩=false); RUN()
    geometricMean(𝐏; w=weights, ✓w=false, ⏩=false); RUN()
    geometricMean(𝐏C; ⏩=false); RUN()
    geometricMean(𝐏C; w=weights, ✓w=false, ⏩=false); RUN()
    geometricMean(𝐃; ⏩=false); RUN()
    geometricMean(𝐃, w=weights, ✓w=false, ⏩=false); RUN()

    name="logdet0Mean"; newTest(name);
    w=[0.5, 0.5]
    P½, P½ⁱ=pow(P_, 0.5, -0.5)
    GM=P½*(P½ⁱ*Q_*P½ⁱ)^0.5*P½  # Fisher mean for k=2
    ldG, iter, conv = logdet0Mean(ℍVector([P_, Q_]); ⏩=false) # logdet0 mean for k=2
    GM ≈ ldG ? OK() : OH(name*" Real Input 1")
    ldG, iter, conv = logdet0Mean(ℍVector([P_, Q_]); w=w, ⏩=false) # weighted logdet0 mean for k=2
    GM ≈ ldG ? OK() : OH(name*" Real Input 2")
    P½, P½ⁱ=pow(PC_, 0.5, -0.5)
    GM=P½*(P½ⁱ*QC_*P½ⁱ)^0.5*P½  # Fisher mean for k=2
    ldG, iter, conv = logdet0Mean(ℍVector([PC_, QC_]); ⏩=false) # logdet0 mean for k=2
    GM ≈ ldG ? OK() : OH(name*" Complex Input 1")
    ldG, iter, conv = logdet0Mean(ℍVector([PC_, QC_]); w=w, ⏩=false) # weighted logdet0 mean for k=2
    GM ≈ ldG ? OK() : OH(name*" Complex Input 2")
    GM=(inv(𝐃[1])*𝐃[2])^0.5*𝐃[1]  # Fisher mean for k=2
    ldG, iter, conv = logdet0Mean(𝔻Vector([𝐃[1], 𝐃[2]]); ⏩=false) # logdet0 mean for k=2
    GM ≈ ldG ? OK() : OH(name*" Real Diagonal Input")

    name="logdet0Mean(⏩ )"; newTest(name);
    logdet0Mean(𝐏); RUN()
    logdet0Mean(𝐏; w=weights, ✓w=false); RUN()
    logdet0Mean(𝐏C); RUN()
    logdet0Mean(𝐏C; w=weights, ✓w=false); RUN()
    logdet0Mean(𝐃); RUN()
    logdet0Mean(𝐃, w=weights, ✓w=false); RUN()

    name="wasMean(⏩ )"; newTest(name);
    wasMean(𝐏); RUN()
    wasMean(𝐏; w=weights); RUN()
    wasMean(𝐏C); RUN()
    wasMean(𝐏C; w=weights); RUN()
    wasMean(𝐃); RUN()
    wasMean(𝐃; w=weights); RUN()

    name="wasMean"; newTest(name);
    wasMean(𝐏; ⏩=false); RUN()
    wasMean(𝐏; w=weights, ✓w=false, ⏩=false); RUN()
    wasMean(𝐏C; ⏩=false); RUN()
    wasMean(𝐏C; w=weights, ✓w=false, ⏩=false); RUN()
    wasMean(𝐃; ⏩=false); RUN()
    wasMean(𝐃, w=weights, ✓w=false, ⏩=false); RUN()

    name="powerMean(⏩ )"; newTest(name);
    powerMean(𝐏, 0.5); RUN()
    powerMean(𝐏, 0.5; w=weights); RUN()
    powerMean(𝐏C, 0.5); RUN()
    powerMean(𝐏C, 0.5; w=weights); RUN()
    powerMean(𝐃, 0.5); RUN()
    powerMean(𝐃, 0.5; w=weights); RUN()

    name="powerMean(⏩ )"; newTest(name);
    powerMean(𝐏, 0.5; ⏩=false); RUN()
    powerMean(𝐏, 0.5; w=weights, ✓w=false, ⏩=false); RUN()
    powerMean(𝐏C, 0.5; ⏩=false); RUN()
    powerMean(𝐏C, 0.5; w=weights, ✓w=false, ⏩=false); RUN()
    powerMean(𝐃, 0.5; ⏩=false); RUN()
    powerMean(𝐃, 0.5, w=weights, ✓w=false, ⏩=false); RUN()

    name="logMap"; newTest(name);
    logMap(Fisher, P, Q); RUN()
    logMap(Fisher, PC, QC); RUN()

    name="expMap"; newTest(name);
    expMap(Fisher, P, Q); RUN()
    expMap(Fisher, PC, QC); RUN()

    G=mean(Fisher, P, Q)
    S=logMap(Fisher, P, G)
    H=expMap(Fisher, S, G)
    if P≉H
        println("")
        @error("either logMap or expMap or both do not give the expected output in the real case")
    end

    G=mean(Fisher, PC, QC)
    S=logMap(Fisher, PC, G)
    H=expMap(Fisher, S, G)
    if PC≉H
        println("")
        @error("either logMap or expMap or both do not give the expected output in the complex case")
    end

    name="vecP"; newTest(name);
    v=vecP(P); RUN()
    vC=vecP(PC); RUN()

    name="matP"; newTest(name);
    Pnew=matP(v); RUN()
    PCnew=matP(vC); RUN()

    if P≉Pnew
        println("")
        @warn("either vecP or matP or both do not give the expected output in the real case")
    end

    if PC≉PCnew
        println("")
        @warn("either vecP or matP or both do not give the expected output in the complex case")
    end

    name="procrustes"; newTest(name);
    U=randU(size(P, 1))
    rotatedP=Hermitian(U'*P*U)
    procrustesU=procrustes(P, rotatedP)
    newP=Hermitian(procrustesU'*rotatedP*procrustesU)
    newP ≈ P ? OK() : OH(name*" Real Input")

    U=randU(ComplexF64, size(PC, 1))
    rotatedP=Hermitian(U'*PC*U)
    procrustesU=procrustes(PC, rotatedP)
    newP=Hermitian(procrustesU'*rotatedP*procrustesU)
    newP ≈ PC ? OK() : OH(name*" Complex Input")

    # functions in classification.jl
    println(" ")
    print("\n- Unit 'statistics.jl'")

    name="softmax"; newTest(name);
    g=[1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0]
    h= [0.0236405430216, 0.0642616585105, 0.17468129856, 0.47483299974438,
        0.0236405430216, 0.0642616585105, 0.17468129856]
    hh=softmax(g)
    hh ≈ h ? OK() : OH(name)

    name="mean (scalar version)"; newTest(name);
    g=[1.0, 2.0, 3.0]
    for m=1:length(metrics)
            if m ∉ (7, 9) mean(metrics[m], g) end end; RUN()

    name="std (scalar version)"; newTest(name);
    g=[1.0, 2.0, 3.0]
    for m=1:length(metrics)
            if m in (Fisher, Euclidean) std(metrics[m], g) end end; RUN()

end # function tests

function newTest(name::String)
    sleep(0.025)
    println(" ")
    print(rpad(name*":", 24))
end

OK()=print("⭐ ")

RUN()=print("☆ ")

function OH(name::String)
    print("⛔ ")
    push!(failing_tests, name)
end

SKIP()=print(" skypped")

failing_tests=[]

function testall()
    println("\n⭐","  PosDefManifold testing utility ", "⭐\n")
    println("Starting tests...\n")
    tests()
    # print out the tests that have failed (if any)
    println("\n")
    length(failing_tests)==0 ? @info("All tests were succesful!") :
        for s in failing_tests @warn("Test of $s failed") end
    for i=1:length(failing_tests) pop!(failing_tests) end
    return nothing
end # function testall

#clipboard("testall()")
#@info("\nhit CTRL+V+ENTER on the REPL for running the tests of PosDefManifold.")
#fails=testAll()


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module CarmichaelNumbers

using Nemo, NumberTheory, Counts

export ModuleCarmichaelNumbers
export isCarmichael, I002997, F002997, L002997
export isweakCarmichael, I225498, F225498, L225498

"""

* isCarmichael, I002997, F002997, L002997
* isweakCarmichael, I225498, F225498, L225498
"""
const ModuleCarmichaelNumbers = ""

"""

Is ``n`` a Carmichael (or Šimerka) number?
"""
function isCarmichael(n)
    (n == 1 || isEven(n) || isPrime(n)) && return false
    for f ∈ Factors(n)
        (f[2] > 1 || (n - 1) % (f[1] - 1) != 0) && return false
    end
    return true
end

"""

Iterate over the first ``n`` Carmichael (or Šimerka) numbers.
"""
I002997(n) = takeFirst(isCarmichael, n)

"""

Iterate over the Carmichael (or Šimerka) numbers which do not exceed ``n``.
"""
F002997(n) = filter(isCarmichael, 1:n)

"""

Return the first ``n`` Carmichael (or Šimerka) numbers in an array.
"""
L002997(n) = collect(I002997(n))

"""

Is ``n`` a weak Carmichael number?
"""
function isweakCarmichael(n)
    (n == 1 || isEven(n) || isPrime(n)) && return false
    for f ∈ Factors(n)
        (n - 1) % (f[1] - 1) != 0 && return false
    end
    return true
end

"""

Iterate over the first ``n`` weak Carmichael numbers.
"""
I225498(n) = takeFirst(isweakCarmichael, n)

"""

Iterate over the weak Carmichael numbers which do not exceed ``n``.
"""
F225498(n) = filter(isweakCarmichael, 1:n)

"""

Return the first ``n`` weak Carmichael numbers in an array.
"""
L225498(n) = collect(I225498(n))

#START-TEST-########################################################

using Test, SeqTests

function test()
    @testset "Carmichael" begin
        @test !isCarmichael(560)
        @test isCarmichael(561)
        @test !isCarmichael(563)
        @test isweakCarmichael(561)
        @test !isweakCarmichael(563)
        @test isweakCarmichael(625)

        if data_installed()
            L = [L002997, L225498]
            SeqTest(L, 'L')
        end

    end
end

function demo()
    for n ∈ 1:30000
        isCarmichael(n) && println(n)
    end

    println()
    for v ∈ I225498(10)
        println(v)
    end
end

"""

@time L = L002997(30)
0.967192 seconds (4.01 M allocations: 248.863 MiB, 28.97% gc time)
"""
function perf()
    @time L = L002997(30)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
using LinearAlgebra

function cheb(N)
    if N==0
        D=0; x=1
        return D,x
    else
        x= cos.(pi * collect(0:N) ./ N)
        c=[2; ones(N-1, 1); 2] .* (-1).^collect(0:N);
        X = repeat(x, 1, N+1);
        dX=X .- X';
        D = (c .* (1 ./ c)')  ./ (dX .+ Diagonal( ones(N+1) ) );
        D = D .- Diagonal( vec( sum(D', dims=1) ) );
        return D,x
    end
end


text:
"""
    crawl_macro(Val(:macroname)[, l::LineNumberNode], exprs...)

Crawl a `:macrocall` expression.
This allows certain macro names to be handled specially by JSExpr (such as
`@var` and `@new`).
"""
function crawl_macrocall end

# Allow override-able parsing of macros
function crawl(h::Val{:macrocall}, m::Symbol, l::LineNumberNode, b...)
    return crawl_macrocall(Val(m), l, b...)
end

function crawl_macrocall(m::Val{M}, ::LineNumberNode, b...) where {M}
    # Remove LineNumberNode info to allow for more convenient crawl_macrocall
    # definitions that don't need it.
    return crawl_macrocall(m, b...)
end

function crawl_macrocall(m::Val{Symbol("@js_str")}, l::LineNumberNode, body...)
    return :(JSTerminal(@js_str $(body...)))
end


text:
"""
```
simulate_switching(m, scen::SwitchingScenario; verbose = :low)
```

Simulate switching in and out of a default scenario for the
SwitchingScenario scen. See `SwitchingScenario` for more
info. Returns a dictionary of results and the proportion
of times a switch actually occured.
"""
function simulate_switching(m::AbstractDSGEModel, scen::SwitchingScenario;
                            verbose::Symbol = :low)

    info_print(verbose, :low, "Simulating switching for " * string(scen.key) * "...")
    println(verbose, :low, "Start time: " * string(now()))
    println(verbose, :low, "Outputs will be saved in " * rawpath(m, "scenarios"))

    tic = time_ns()
    # Revert model alt policy to historical rule
    m <= Setting(:alternative_policy, AltPolicy(:historical, solve, eqcond), false, "apol",
                 "Alternative policy")

    # Get output file names
    original_output_files = get_scenario_output_files(m, scen.original,
                                                      [:forecastobs, :forecastpseudo])
    default_output_files = get_scenario_output_files(m, scen.default,
                                                     [:forecastobs, :forecastpseudo])

    results = Dict{Symbol, Array{Float64}}()
    switching_results = [0.0, 0.0]

    for (i, output_var) in enumerate([:forecastobs, :forecastpseudo])
        # Read in original and default draws
        original_draws = JLD2.jldopen(original_output_files[output_var], "r") do file
            read(file, "arr")
        end
        default_draws  = JLD2.jldopen(default_output_files[output_var], "r") do file
            read(file, "arr")
        end

        n_draws, n_vars, n_periods = size(original_draws)
        n_default_draws = size(default_draws, 1)

        @assert n_vars == size(default_draws, 2)
        @assert length(scen.probs_enter) == n_periods <= size(default_draws, 3)

        # Simulate switching
        results[output_var] = zeros(n_draws, n_vars, n_periods)
        proportion_switched = 0.0
        for j = 1:n_draws
            k = rand(1:n_default_draws)
            results[output_var][j, :, :], switched = switch(original_draws[j, :, :], default_draws[k, :, :],
                                                            scen.probs_enter, scen.probs_exit)
            proportion_switched += Float64(switched)
        end
        switching_results[i] = proportion_switched / n_draws
    end

    results[:proportion_switched] = [mean(switching_results)]

    # Write outputs
    output_files = get_scenario_output_files(m, scen, [:forecastobs, :forecastpseudo])
    write_scenario_forecasts(m, output_files, results, verbose = verbose)

    # Print
    switching_time = (time_ns() - tic)/1e9
    switching_time_min = switching_time/60
    println(verbose, :low, "\nTime elapsed: " * string(switching_time_min) * " minutes")
    println(verbose, :low, "Switching complete: " * string(now()))

    return results, switching_results
end

#for switching with multiple models
function simulate_switching(m_original::AbstractDSGEModel, m_default::AbstractDSGEModel, scen::SwitchingScenario;
                            verbose::Symbol = :low)

    info(verbose, :low, "Simulating switching for " * string(scen.key) * "...")
    println(verbose, :low, "Start time: " * string(now()))
    println(verbose, :low, "Outputs will be saved in " * rawpath(m_original, "scenarios"))
    tic()

    # Revert model alt policy to historical rule
    m_original <= Setting(:alternative_policy, AltPolicy(:historical, solve, eqcond), false, "apol",
                 "Alternative policy")
    m_default <= Setting(:alternative_policy, AltPolicy(:historical, solve, eqcond), false, "apol",
                 "Alternative policy")

    # Get output file names
    original_output_files = get_scenario_output_files(m_original, scen.original,
                                                      [:forecastobs, :forecastpseudo])
    default_output_files = get_scenario_output_files(m_default, scen.default,
                                                     [:forecastobs, :forecastpseudo])

    results = Dict{Symbol, Array{Float64}}()
    switching_results = [0.0, 0.0]

    for (i, output_var) in enumerate([:forecastobs, :forecastpseudo])
        if output_var==:forecastobs
            include_obs = intersect(keys(m_default.observables), keys(m_original.observables))
            inds_default = Vector(0)
            inds_original = Vector(0)
            for key in include_obs
                inds_default = vcat(inds_default, m_default.observables[key])
                inds_original = vcat(inds_original, m_original.observables[key])
            end
        elseif output_var==:forecastpseudo
            include_obs = intersect(keys(m_default.pseudo_observables), keys(m_original.pseudo_observables))
            inds_default = Vector(0)
            inds_original = Vector(0)
            for key in include_obs
                inds_default = vcat(inds_default, m_default.pseudo_observables[key])
                inds_original = vcat(inds_original, m_original.pseudo_observables[key])
            end
        end

        # Read in original and default draws
        original_draws = load(original_output_files[output_var], "arr")[:, inds_original, :]
        default_draws  = load(default_output_files[output_var], "arr")[:, inds_default, :]

        n_draws, n_vars, n_periods = size(original_draws)
        n_default_draws = size(default_draws, 1)

        @assert n_vars == size(default_draws, 2)
        @assert length(scen.probs_enter) == n_periods <= size(default_draws, 3)

        # Simulate switching
        results[output_var] = zeros(n_draws, n_vars, n_periods)
        proportion_switched = 0.0
        for j = 1:n_draws
            k = rand(1:n_default_draws)
            results[output_var][j, :, :], switched = switch(original_draws[j, :, :], default_draws[k, :, :],
                                                            scen.probs_enter, scen.probs_exit)
            proportion_switched += Float64(switched)
        end
        switching_results[i] = proportion_switched / n_draws
    end

    results[:proportion_switched] = [mean(switching_results)]

    # Write outputs
    output_files = get_scenario_output_files(m_original, scen, [:forecastobs, :forecastpseudo])
    write_scenario_forecasts(m_original, output_files, results, verbose = verbose)

    # Print
    switching_time = toq()
    switching_time_min = switching_time/60
    println(verbose, :low, "\nTime elapsed: " * string(switching_time_min) * " minutes")
    println(verbose, :low, "Switching complete: " * string(now()))

    return results, switching_results
end

"""
```
switch(original, default, probs_enter, probs_exit)
```

Simulate entry and exit from the original vector into the default one
according to the specified probabilities.
"""
function switch(original::Matrix{Float64}, default::Matrix{Float64},
                probs_enter::Vector{Float64}, probs_exit::Vector{Float64})

    n_vars, n_periods = size(original)
    output = zeros(n_vars, n_periods)

    # Find period of switch from default to original
    period_in = choose_switch_period(probs_enter, 1)

    # Record whether original is actually switched into
    switched = period_in < n_periods

    # Find period of switch from original to default
    # Start in period_in + 1 s.t. at least 1 period is spent in the original scenario
    period_out = choose_switch_period(probs_exit, period_in + 1)

    output[:, 1:(period_in-1)]          = default[:, 1:(period_in-1)]
    output[:, period_in:(period_out-1)] = original[:, 1:(period_out-period_in)]
    output[:, period_out:n_periods]     = default[:, period_out:n_periods]

    return output, switched
end

function choose_switch_period(probs::Vector{Float64}, first_period::Int)
    n_periods = length(probs)

    for t = first_period:n_periods
        if rand() < probs[t]
            return t
        end
    end

    return n_periods + 1
end


text:
"""
    Search(dataSource::DataSource, intervals::Vector{<:Interval}, query::SearchQuerySpec; <keyword arguments>)

Returns a Druid.SearchResult which is compatible with the Tables.jl interface.

A method with all arguments as keyword arguments is also provided.

# Arguments
- granularity::Granularity = nothing
- filter::Filter = nothing
- virtualColumns::Vector{<:VirtualColumn} = nothing
- sort::String = nothing
- limit::Integer = nothing
- context::Dict = nothing
"""
mutable struct Search <: Query
    queryType::String
    dataSource::DataSource
    intervals::Vector{<:Interval}
    query::SearchQuerySpec
    searchDimensions
    granularity
    filter
    virtualColumns
    sort
    limit
    context
    function Search(
        dataSource, intervals, query;
        searchDimensions=nothing, granularity=nothing, filter=nothing,
        virtualColumns=nothing, sort=nothing, limit=nothing, context=nothing
    )
        nothing_or_type(searchDimensions, Vector{<:DimensionSpec})
        nothing_or_type(granularity, Granularity)
        nothing_or_type(filter, Filter)
        nothing_or_type(virtualColumns, Vector{<:VirtualColumn})
        sort === nothing || (sort = lowercase(sort)) ∈ ["lexicographic", "alphanumeric", "numeric", "strlen"] || error("Invalid sort value")
        limit === nothing || (isa(limit, Integer) && limit >= 0) || error("limit must be a non-negative integer")
        nothing_or_type(context, Dict)
        new("search", dataSource, intervals, query, searchDimensions, granularity, filter, virtualColumns, Dict("type" => sort), limit, context)
    end
end
Search(
    ; dataSource, intervals, query,
    searchDimensions=nothing, granularity=nothing, filter=nothing, virtualColumns=nothing, sort=nothing, limit=nothing, context=nothing
) = Search(
    dataSource, intervals, query;
    searchDimensions, granularity, filter, virtualColumns, sort, limit, context
)

function execute(client::Client, query::Search; pretty=false)
    res = convert(Vector{Dict}, JSON.parse(execute_native_query(client, query; pretty)))
    names = [:timestamp, Symbol.(keys(res[1]["result"][1]))...]
    num_rows = sum(length(subres["result"]) for subres ∈ res)
    SearchResult(names, num_rows, res)
end

struct SearchResult <: QueryResult
    names::Vector{Symbol}
    num_rows::Int
    inner_array::Vector{Dict}
end

names(sr::SearchResult) = getfield(sr, :names)

Base.summary(io::IO, sr::SearchResult) = print(io, string(length(sr)) * "-element " * string(typeof(sr)))
Base.getindex(sr::SearchResult, i::Int) = (i <= length(sr) || throw(BoundsError(sr, i))) && SearchRow(i, sr)

Tables.rowaccess(::SearchResult) = true
Tables.rows(sr::SearchResult) = sr

Base.eltype(::SearchResult) = SearchRow
Base.length(sr::SearchResult) = getfield(sr, :num_rows)
Base.iterate(sr::SearchResult, state=1) = state > length(sr) ? nothing : (SearchRow(state, sr), state + 1)

function Tables.schema(sr::SearchResult)
    if length(sr) != 0
        types = [typeof(cl) for cl in sr[1]]
    else
        types = DataType[]
    end
    Tables.Schema(names(sr), types)
end

struct SearchRow <: Tables.AbstractRow
    row::Int
    source::SearchResult
end

function Tables.getcolumn(sr::SearchRow, name::Symbol)
    elems = 0
    i = getfield(sr, :row)
    for subres ∈ getfield(getfield(sr, :source), :inner_array)
        subelems = length(subres["result"])
        if i <= elems + subelems
            row = Dict(subres["result"][i - elems])
            row["timestamp"] = subres["timestamp"]
            return row
        end
        elems += subelems
    end
end

Tables.getcolumn(sr::SearchRow, i::Int) = Tables.getcolumn(sr, names(getfield(sr, :source))[i])
Tables.columnnames(sr::SearchRow) = names(getfield(sr, :source))


text:
include("../roblib.jl")

function find_marker!(robot::Robot)
    robot = interface_protected_robot(robot)
    line = interface_line(robot.move!)

    function move_ifnomarker!(side)::Bool
        if robot.ismarker()
            return false
        end
        move_unlimited_line!(side)
        return true
    end

    function move_unlimited_line!(side)
        coordinate = 0 # - текущая координата Робота относительно начального положения
        num_steps = 1 # - текущее число шагов при поиске края перегородки
        orthogonal = right(side) # - текущее направление поиска края перегородки
        back_side = orthogonal # - первоначальное направление поиска края перегородки
        while robot.isborder(side) # поиск края полубесконечной прямолинейной перегородки
            line.movements!(orthogonal,num_steps)
            if orthogonal == back_side
                coordinate += num_steps
            else
                coordinate -= num_steps
            end
            num_steps+=1
            orthogonal=inverse(orthogonal)
        end
            #УТВ: Робот - за краем перегородки
        robot.move!(side)
        if coordinate > 0
            back_side = inverse(back_side)
        else
            coordinate = -coordinate
        end
        line.movements!(back_side,coordinate)
            #УТВ: Робот - в соседней, по отношению к исходной, позиции с другой стороны от перегородки
    end

    trajectories = interface_trajectories(robot)
    trajectories.spiral!(move_ifnomarker!)
end


text:
"Use schedule_MATTER_watcher to launch in a task"
function MATTER_watcher(masschannels, exitchannel, setup_channel_func, maxchannels)
	println("starting MATTER_watcher")
	ljhname, writingbool = "",false
	last_noise_filename = ""
	last_started_analyzing_fname = "" # not neccesary, but I could see using the info in the future
	have_planned_to_end = true
	i=0
	while true
		watch_file(LJHUtil.sentinel_file_path,10) #blocks task until file changes, or timeout given in seconds
		oldljhname, oldwritingbool = ljhname, writingbool
		ljhname, writingbool = LJHUtil.matter_writing_status()
		if (ljhname != oldljhname)
			if contains(ljhname, "noise") || contains(ljhname,".noi")
				last_noise_filename = ljhname
				println("New noise file $last_noise_filename")
			else
				println("new ljh file observed: $ljhname")
				if last_noise_filename == ""
					println("must set last_noise_filename before processing data")
					continue
				end
				last_started_analyzing_fname = ljhname
				println("Starting analysis of $ljhname with noise from $last_noise_filename")
				sleep(2) # wait for LJH files to exist, a bit hacky
				t0 = time()
				channums = LJHUtil.allchannels(ljhname)
				channums = channums[1:min(maxchannels,length(channums))]
				ljh_filenames = [LJHUtil.fnames(ljhname,channum) for channum in channums]
				noise_filenames = [LJHUtil.fnames(last_noise_filename,channum) for channum in channums]
				for (i, channum) in enumerate(channums)
					masschannels[channums[i]] = setup_channel_func(ljh_filenames[i], noise_filenames[i])
				end
				tf = time()
				println("setup $(length(channums)) channels in $(tf-t0) seconds")
				tschedule = @elapsed map(schedule, values(masschannels))
				have_planned_to_end = false
				println("launched $(length(channums)) channels in $tschedule seconds")
			end
		end
		if !writingbool && !have_planned_to_end
			# stop tasks for previous file once they've finished their work
			println("allowing analysis tasks for $oldljhname to end")
			map(plantoend, values(masschannels))
			have_planned_to_end=true
		end
		if isready(exitchannel)
			println("exiting MATTER_watcher, due to ready exitchannel")
			return
		end
	end
end

"`masschannels, watcher_exitchannel, watcher_task = schedule_MATTER_watcher(setup_channel_func, maxchannels=1000)`
`setup_channel_func(ljh_filename, noise_filename)` should be a function that returns a `MassChannel`, designed to analyze the given files.
`maxchannels` is an `Int` that limits the maximum number of channels to analyze
Launches MATTER_watcher in a task. It will monitor sentinel_file_path for changes. It stores the name
of files that contain the text `noise` or `.noi` as noise files, then starts analysis on any other files, using
the last seen noise file as the noise. When writing stops, or when a new file starts being written it
plans to end the analysis steps. This means that when all work has finished, the analysis tasks will end.
`masschannels` is a Dictionary channel_number->`MassChannel` for the currently analyzing files
`watcher_exitchannel` will cause the MATTER_watcher task to end upon `put!(watcher_exitchannel,1)` (or any other `Int`)
`watcher_task` is the task executing MATTER_watcher"
function schedule_MATTER_watcher(setup_channel_func, maxchannels=1000)
	masschannels = Dict()
	exitchannel = Channel{Int}(1)
	task = @schedule MATTER_watcher(masschannels, exitchannel, setup_channel_func, maxchannels)
	masschannels, exitchannel, task
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[1.23603, 1.34652, 1.31271, 1.00791, 1.48861, 1.21097, 1.95192, 1.9999, 1.25166, 1.98667],
     y=[1.55575, 1.43711, 1.42472, 1.77322, 1.28119, 1.20947, 1.25138, 1.02037, 1.2877, 1.85951],
     Scale.x_log10(minvalue=1.0, maxvalue=10),
     Scale.y_log10(minvalue=1.0, maxvalue=10))


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/IainNZ/Coverage.jl
#######################################################################
module Coverage

    import JuliaParser.Parser
    using Compat

    # process_cov
    # Given a .cov file, return the counts for each line, where the
    # lines that can't be counted are denoted with a -1
    export process_cov, amend_coverage_from_src!, coverage_file, coverage_folder, analyze_malloc
    function process_cov(filename)
        if !isfile(filename)
            srcname, ext = splitext(filename)
            lines = open(srcname) do fp
                readlines(fp)
            end
            coverage = Array(Union(Nothing,Int), length(lines))
            return fill!(coverage, nothing)
        end
        fp = open(filename, "r")
        lines = readlines(fp)
        num_lines = length(lines)
        coverage = Array(Union(Nothing,Int), num_lines)
        for i = 1:num_lines
            cov_segment = lines[i][1:9]
            coverage[i] = cov_segment[9] == '-' ? nothing : int(cov_segment)
        end
        close(fp)
        return coverage
    end
    function amend_coverage_from_src!(coverage, srcname)
        # To make sure things stay in sync, parse the file position corresonding to each new line
        linepos = Int[]
        open(srcname) do io
            while !eof(io)
                push!(linepos, position(io))
                readline(io)
            end
            push!(linepos, position(io))
        end
        open(srcname) do io
            while !eof(io)
                pos = position(io)
                linestart = minimum(searchsorted(linepos, pos))
                ast = Parser.parse(io)
                isa(ast, Expr) || continue
                flines = function_body_lines(ast)
                if !isempty(flines)
                    flines += linestart-1
                    for l in flines
                        if coverage[l] == nothing
                            coverage[l] = 0
                        end
                    end
                end
            end
        end
        coverage
    end
    function coverage_file(filename)
        results = Coveralls.process_file(filename)
        coverage = results["coverage"]
        tot = sum(x->x!=nothing, coverage)
        covered = sum(x->x!=nothing && x>0, coverage)
        covered, tot
    end
    function coverage_folder(folder="src")
        results = Coveralls.process_folder(folder)
        tot = covered = 0
        for item in results
            coverage = item["coverage"]
            tot += sum(x->x!=nothing, coverage)
            covered += sum(x->x!=nothing && x>0, coverage)
        end
        covered, tot
    end

    function_body_lines(ast) = function_body_lines!(Int[], ast, false)
    function_body_lines!(flines, arg, infunction) = flines
    function function_body_lines!(flines, node::LineNumberNode, infunction)
        line = node.line
        if infunction
            push!(flines, line)
        end
        flines
    end
    function function_body_lines!(flines, ast::Expr, infunction)
        if ast.head == :line
            line = ast.args[1]
            if infunction
                push!(flines, line)
            end
            return flines
        end
        infunction |= Base.Cartesian.isfuncexpr(ast)
        for arg in ast.args
            flines = function_body_lines!(flines, arg, infunction)
        end
        flines
    end

    export Coveralls
    module Coveralls
        using Requests
        using Coverage
        using JSON
        using Compat

        # coveralls_process_file
        # Given a .jl file, return the Coveralls.io dictionary for this
        # file by reading in the file and its matching .cov. Don't convert
        # to JSON yet, just return dictionary.
        # https://coveralls.io/docs/api
        # {
        #   "name" : "$filename"
        #   "source": "...\n....\n...."
        #   "coverage": [null, 1, null]
        # }
        export process_file
        function process_file(filename)
            return @compat Dict("name" => filename,
                    "source" => readall(filename),
                    "coverage" => amend_coverage_from_src!(process_cov(filename*".cov"), filename))
        end

        # coveralls_process_src
        # Recursively walk through a Julia package's src/ folder
        # and collect coverage statistics
        export process_folder
        function process_folder(folder="src")
            source_files=Any[]
            filelist = readdir(folder)
            for file in filelist
                fullfile = joinpath(folder,file)
                println(fullfile)
                if isfile(fullfile)
                    try
                        new_sf = process_file(fullfile)
                        push!(source_files, new_sf)
                    catch e
#                         if !isa(e,SystemError)
#                             rethrow(e)
#                         end
                        # Skip
                        println("Skipped $fullfile")
                    end
                else isdir(fullfile)
                    append!(source_files, process_folder(fullfile))
                end
            end
            return source_files
        end

        # submit
        # Submit coverage to Coveralls.io
        # https://coveralls.io/docs/api
        # {
        #   "service_job_id": "1234567890",
        #   "service_name": "travis-ci",
        #   "source_files": [
        #     {
        #       "name": "example.rb",
        #       "source": "def four\n  4\nend",
        #       "coverage": [null, 1, null]
        #     },
        #     {
        #       "name": "lib/two.rb",
        #       "source": "def seven\n  eight\n  nine\nend",
        #       "coverage": [null, 1, 0, null]
        #     }
        #   ]
        # }
        export submit, submit_token
        function submit(source_files)
            data = @compat Dict("service_job_id" => ENV["TRAVIS_JOB_ID"],
                    "service_name" => "travis-ci",
                    "source_files" => source_files)
            r = Requests.post(URI("https://coveralls.io/api/v1/jobs"), files =
                [FileParam(JSON.json(data),"application/json","json_file","coverage.json")])
            dump(r.data)
        end

        function submit_token(source_files)
            data = @compat Dict("repo_token" => ENV["REPO_TOKEN"],
                    "source_files" => source_files)
            r = post(URI("https://coveralls.io/api/v1/jobs"), files =
                [FileParam(JSON.json(data),"application/json","json_file","coverage.json")])
            dump(r.data)
        end
    end  # module Coveralls

    ## Analyzing memory allocation
    immutable MallocInfo
        bytes::Int
        filename::UTF8String
        linenumber::Int
    end

    sortbybytes(a::MallocInfo, b::MallocInfo) = a.bytes < b.bytes

    function analyze_malloc_files(files)
        bc = MallocInfo[]
        for filename in files
            open(filename) do file
                for (i,ln) in enumerate(eachline(file))
                    tln = strip(ln)
                    if !isempty(tln) && isdigit(tln[1])
                        s = split(tln)
                        b = parseint(s[1])
                        push!(bc, MallocInfo(b, filename, i))
                    end
                end
            end
        end
        sort(bc, lt=sortbybytes)
    end

    function find_malloc_files(dirs)
        files = ByteString[]
        for dir in dirs
            filelist = readdir(dir)
            for file in filelist
                file = joinpath(dir, file)
                if isdir(file)
                    append!(files, find_malloc_files(file))
                elseif endswith(file, "jl.mem")
                    push!(files, file)
                end
            end
        end
        files
    end
    find_malloc_files(file::ByteString) = find_malloc_files([file])

    analyze_malloc(dirs) = analyze_malloc_files(find_malloc_files(dirs))
    analyze_malloc(dir::ByteString) = analyze_malloc([dir])

    # Support Unix command line usage like `julia Coverage.jl $(find ~/.julia/v0.3 -name "*.jl.mem")`
    if !isinteractive()
        bc = analyze_malloc_files(ARGS)
        println(bc)
    end
end


text:
export Vertex, DVertex, IVertex, vertex, dvertex

import Base: copy, hash, ==, <, <<

abstract Vertex{T}

Base.eltype{T}(::Vertex{T}) = T

Base.show{T<:Vertex}(io::IO, ::Type{T}) =
  print(io, T.name.name, (T.parameters[1] == Any ? [] : ["{", T.parameters[1], "}"])...)

include("set.jl")
include("dlgraph.jl")
include("ilgraph.jl")
include("conversions.jl")

thread!(to::Vertex, from) = thread!(to, convert(typeof(to), from))

thread!(v::Vertex, xs...) = foldl(thread!, v, xs)

(::Type{T}){T<:Vertex}(x, args...) = thread!(T(x), args...)

Base.getindex(v::Vertex, i::Integer) = inputs(v)[i]
Base.getindex(v::Vertex, is::Integer...) = foldl(getindex, v, is)

function collectv(v::Vertex, vs = OASet{typeof(v)}())
  v ∈ vs && return collect(vs)
  push!(vs, v)
  foreach(v′ -> collectv(v′, vs), inputs(v))
  foreach(v′ -> collectv(v′, vs), outputs(v))
  return collect(vs)
end

function topo_up(v::Vertex, vs, seen)
  v ∈ seen && return vs
  push!(seen, v)
  foreach(v′ -> topo_up(v′, vs, seen), inputs(v))
  push!(vs, v)
end

function topo(v::Vertex)
  seen, vs = OSet{typeof(v)}(), typeof(v)[]
  for v in sort!(collectv(v), by = x -> x ≡ v)
    topo_up(v, vs, seen)
  end
  return vs
end

function isreaching(from::Vertex, to::Vertex, seen = OSet())
  to ∈ seen && return false
  push!(seen, to)
  any(v -> v ≡ from || isreaching(from, v, seen), inputs(to))
end

Base.isless(a::Vertex, b::Vertex) = isreaching(a, b)

<<(a::Vertex, b::Vertex) = a < b && !(a > b)

↺(v::Vertex) = v < v
↺(a::Vertex, b::Vertex) = a < b && b < a


text:
import Plots
import PyPlot
import StatsPlots

"""
    plotbench(bench_result, plot_filename; kwargs...)
"""
function plotbench(bench_result::BenchmarkResult,
                   plot_filename::AbstractString;
                   plotting_backend = Plots.pyplot)
    df = benchmark_result_df(bench_result)
    df[!, :Size] = Float64.(df[!, :Size])
    df[!, :TFLOPS] = Float64.(df[!, :TFLOPS])
    if !(plotting_backend isa Nothing)
        plotting_backend()
    end
    p = StatsPlots.@df df StatsPlots.plot(
        :Size,
        :TFLOPS;
        group = :Library,
        legend = :outertopright,
        xscale = :log2,
        xlabel = "Size",
        ylabel = "TFLOPS",
    );
    Plots.savefig(p, plot_filename)
    return nothing
end


text:
abstract type Model end
abstract type TerraceOnly <:Model end
abstract type TerraceKnickpoint <:Model end
abstract type TerraceDiffusion <:Model end
abstract type TerraceKnickpointDiffusion <:Model end

struct ModelOpts{T} end

function solve(solver::Type{ModelOpts{TerraceOnly}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer,
    ScratchTerrace_FEM, Scratch_FEM, to)

    @timeit to "newton terrace" terracenewton2!(
        TerracePhysics,
        (Δt * yr),
        Terrace.z,
        h_sea,
        id_shore.terrace)

    return Terrace, River,to
end

function solve(solver::Type{ModelOpts{TerraceKnickpoint}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer,
    ScratchTerrace_FEM, Scratch_FEM, to)

    # RIVER PROFILE
    @timeit to "newton river" terracenewton2!(
        TerracePhysics,
        Δt * yr,
        River.z,
        h_sea,
        id_shore.river,
        )

    # TERRACE PROFILE
    @timeit to "newton terrace" terracenewton2!(
        TerracePhysics,
        (Δt * yr),
        Terrace.z,
        h_sea,
        id_shore.terrace
        )

    # KNICKPOINT MIGRATION ADVECTION -> FLUX LIMITER TVD
    @timeit to "update r" updater!(RiverPhysics, RiverArrays, River)
    @timeit to "TVD" TVD!(River, buffer, RiverPhysics, RiverArrays, Δt, id_shore.river.i1)

    return Terrace, River,to
end

function solve(solver::Type{ModelOpts{TerraceKnickpointDiffusion}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer,
    ScratchTerrace_FEM, Scratch_FEM, to)

    # -- RIVER PROFILE
    @timeit to "newton river" terracenewton2!(
        TerracePhysics,
        Δt * yr,
        River.z,
        h_sea,
        id_shore.river,
    )

    # -- TERRACE PROFILE
    @timeit to "newton terrace" terracenewton2!(
        TerracePhysics,
        (Δt * yr),
        Terrace.z,
        h_sea,
        id_shore.terrace
    )

    # KNICKPOINT MIGRATION ADVECTION -> FLUX LIMITER TVD
    @timeit to "update r" updater!(RiverPhysics, RiverArrays, River)
    @timeit to "TVD" TVD!(River, buffer, RiverPhysics, RiverArrays, Δt, id_shore.river.i1)

    # Diffusion over Δt (FEM)
    @timeit to "FEM River" River.z = femsolver(Δt, Δx, nn, River, h_sea,
                                               Scratch_FEM,
                                               HillSlopePhysics,
                                               SubmarinePhysics)

    # TERRACE
    # Diffusion over Δt
    @timeit to "FEM Terrace" Terrace.z = femsolver(
            Δt,
            Terrace,
            h_sea,
            ScratchTerrace_FEM,
            Scratch_FEM,
            HillSlopePhysics
    )

    return Terrace, River,to
end

function solve(solver::Type{ModelOpts{TerraceDiffusion}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer,
    ScratchTerrace_FEM, Scratch_FEM, to)

    # -- RIVER PROFILE
    @timeit to "newton river" terracenewton2!(
        TerracePhysics,
        Δt * yr,
        River.z,
        h_sea,
        id_shore.river,
    )

    # TERRACE
    @timeit to "FEM Terrace" Terrace.z = femsolver(
            Δt,
            Terrace,
            h_sea,
            ScratchTerrace_FEM,
            Scratch_FEM,
            HillSlopePhysics
    )

    return Terrace, River,to
end

function old_solver!(Terrace, TerracePhysics, h_sea, id_shore, dt)
    P, h_wb, βz = TerracePhysics.P0, TerracePhysics.h_wb, TerracePhysics.βz
    for i in id_shore.terrace.i1: id_shore.terrace.i2
        @inbounds Terrace.z[i] += - βz .* P .* exp(-4*(h_sea-Terrace.z[i])/h_wb) * dt
    end
end

#-------------------------------------------------------------------------
function terracenewton!(TerracePhysics::TerraceParameters,dt::Float64,
    sol::Array{Float64},h_sea::Float64,infected_nodes::NamedTuple,sol_n,sol_old)
    #=
        Add derivation of newton-raphson
    =#
    P, h_wb, βz = TerracePhysics.P0, TerracePhysics.h_wb, TerracePhysics.βz
    cte1 = βz * P * dt
    cte2 = 4.0/h_wb
    cte3 = cte1 * cte2
    copyto!(sol_old, sol)
    residual = 1.0
    i1, i2 = infected_nodes.i1, infected_nodes.i2
    while residual > 1e-8
        copyto!(sol_n,sol)
        @tturbo for ii in i1:i2
            sol[ii] -= (sol[ii] + cte1 * exp(-cte2*(h_sea-sol[ii])) - sol_old[ii]) * (1/ (1 + cte3 * exp(-cte2*(h_sea-sol[ii]))))
        end
        residual = norm2(sol, sol_n, infected_nodes)
    end
end

function _newton(sol,h_sea, cte1, cte2)
    h0, h, hprevious = sol, sol, sol
    residual = 1.0
    its = 0
    while its < 21 && residual > 1e-8
        cte3 = cte1 * exp(-cte2*(h_sea-h))
        up = h - hprevious + cte3
        down = 1/(1 - cte3)
        h -= up * down
        residual = abs(h-h0)/abs(h0)
        h0 = h
        its += 1
    end
    h
end

function terracenewton2!(TerracePhysics::TerraceParameters,dt::Float64,
    sol::Array{Float64},h_sea::Float64,infected_nodes::NamedTuple)
    #=
        Add derivation of newton-raphson
    =#
    P, h_wb, βz = TerracePhysics.P0, TerracePhysics.h_wb, TerracePhysics.βz
    cte1 = βz * P * dt
    cte2 = 4.0/h_wb
    i1, i2 = infected_nodes.i1, infected_nodes.i2
    for ii in i1:i2
        @inbounds sol[ii] = _newton(sol[ii] ,h_sea, cte1, cte2)
    end

end

#-------------------------------------------------------------------------
function terracenewton!(βz::Float64,P::Float64,h_wb::Float64,dt::Float64,
    sol::Array{Float64},h_sea::Float64,i1::Int64,sol_n,sol_old)
    #=
        Add derivation of newton-raphson
    =#
    cte1 = βz * P * dt
    cte2 = 4.0/h_wb
    cte3 = cte1 * cte2
    copyto!(sol_old, sol)
    residual    = 1.0
    while residual > 1e-3
        copyto!(sol_n,sol)
        @turbo for ii = i1:length(sol)
            cte4    = exp(-cte2*(h_sea-sol[ii]))
            sol[ii]-= (sol[ii] + cte1 * cte4 - sol_old[ii]) * (1/ (1.0 + cte3 * cte4))
        end
        residual = norm2(sol, sol_n)/norm(sol_n)
    end
end

# CLIFF RETREAT ================================================================
function cliffretreat!(Terrace,h_sea,P0,h_wb,βx,Δt)
    shelf_depth         = 10; # max depth of continental shelf
    bool1               = @. (abs(h_sea-Terrace.z)<shelf_depth)
    bool2               = @. h_sea>=Terrace.z
    id_shelf            = [bool1[kk] == true && bool2[kk] == true ? true : false for kk = 1:length(bool1)]
    z_sea_terrace       = Terrace.z[id_shelf];
    h2sea_terrace       = @. abs(h_sea - z_sea_terrace);

    # -- Cliff retreat
    a                   = Terrace.x[id_shelf]
    b                   = @. P0*exp(-4*h2sea_terrace/h_wb)
    integral            = integrate(a,b,TrapezoidalEvenFast())
    # integral            = integrate1D(a,b)
    dx_cliff            = (P0 - integral) * Δt * βx
    # -- Modify geometry
    x_terrace               = deepcopy(Terrace.x)
    id_cliff                = findlast(x->x>h_sea,Terrace.z)  # find clif position
    x_terrace[1:id_cliff] .-= abs(dx_cliff)
    # iz_correct              = id_cliff + 1
    idx                     = id_cliff-25:id_cliff+25
    Terrace.z[idx]         .= interpolate((x_terrace[idx],), Terrace.z[idx], Gridded(Linear()))(Terrace.x[idx])
    Terrace.z               = interpolate((x_terrace,), Terrace.z, Gridded(Linear()))(Terrace.x)
    Terrace.z[id_cliff]     = Terrace.z[id_cliff+1]

end
# ==============================================================================

@inline integrate1D(x,f) = sum(f*(x[2]-x[1]))

@inline function fixterracex!(x,id)
    x0 = x[id]
    for i ∈ 1:id-1
        @inbounds if x[i] >= x0
            x[i] -= x0 - 0.01*i
        end
    end
end