text:
## Testing the New Parameter API

import Mimi: model_param, is_shared

#
# Section 1. update_param!, add_shared_param! and connect_param!
#

@defcomp A begin

    p1 = Parameter{Symbol}()
    p2 = Parameter(default = 2)
    p3 = Parameter()
    p4 = Parameter(unit = "dollars")
    p5 = Parameter(unit = "\$")
    p6 = Parameter(index = [time])
    p7 = Parameter(index = [regions, time])

    function run_timestep(p,v,d,t)
    end
end

function _get_model()
    m = Model()
    set_dimension!(m, :time, 1:5);
    set_dimension!(m, :regions, [:R1, :R2, :R3])
    add_comp!(m, A)
    return m
end

# General Functionality
m = _get_model()

@test_throws MethodError update_param!(m, :A, :p1, 3) # can't convert
update_param!(m, :A, :p1, :hello)
@test model_param(m, :A, :p1).value == :hello
add_shared_param!(m, :p1_fail, 3)
@test_throws ErrorException connect_param!(m, :A, :p1, :p1_fail) # we throw specific error here
add_shared_param!(m, :p1, :goodbye)
connect_param!(m, :A, :p1, :p1)
@test model_param(m, :A, :p1).value == :goodbye
@test_throws ErrorException update_param!(m, :A, :p1, :foo) # can't call this method on a shared parameter
update_param!(m, :p1, :foo)
@test model_param(m, :A, :p1).value == :foo
disconnect_param!(m, :A, :p1)
update_param!(m, :A, :p1, :foo) # now we can update :p1 with this syntax since it was disconnected
update_param!(m, :p1, :bar) # this is the shared parameter named :p1
@test model_param(m, :A, :p1).value == :foo
@test model_param(m, :p1).value == :bar

m = _get_model()

add_shared_param!(m, :shared_param, 100)
connect_param!(m, :A, :p2, :shared_param)
connect_param!(m, :A, :p3, :shared_param)
@test model_param(m, :A, :p2).value == model_param(m, :A, :p3).value == 100

# still error because they are connected to a shared parameter
@test_throws ErrorException update_param!(m, :A, :p2, 1)
@test_throws ErrorException update_param!(m, :A, :p3, 2)

disconnect_param!(m, :A, :p2)
update_param!(m, :A, :p2, 1)
model_param(m, :A, :p2).value == 1
model_param(m, :A, :p3).value == model_param(m, :shared_param).value == 100

# Defaults
m = _get_model()

@test model_param(m, :A, :p2).value == 2
@test !(is_shared(model_param(m, :A, :p2)))
update_param!(m, :A, :p2, 100)
@test !(is_shared(model_param(m, :A, :p2)))

# Dimensions
m = _get_model()

@test_throws ErrorException add_shared_param!(m, :x, [1:10]) # need dimensions to be specified
@test_throws ErrorException add_shared_param!(m, :x, [1:10], dims = [:time]) # wrong dimensions
add_shared_param!(m, :x, 1:5, dims = [:time])

@test_throws ErrorException add_shared_param!(m, :y, fill(1, 3, 5)) # need dimensions to be specified
@test_throws ErrorException add_shared_param!(m, :y, fill(1, 3, 5), dims = [:time, :regions]) # need dimensions to be specified
add_shared_param!(m, :y, fill(1, 5, 3), dims = [:time, :regions])

@test_throws ErrorException connect_param!(m, :A, :p7, :y) # wrong dimensions, flipped around

# Units and Datatypes
m = _get_model()

add_shared_param!(m, :myparam, 100)
connect_param!(m, :A, :p3, :myparam)
@test_throws ErrorException connect_param!(m, :A, :p4, :myparam) # units error
connect_param!(m, :A, :p4, :myparam; ignoreunits = true)
@test model_param(m, :A, :p3).value == model_param(m, :A, :p4).value == 100
@test_throws ErrorException update_param!(m, :myparam, :boo) # cannot convert
update_param!(m, :myparam, 200)
@test model_param(m, :A, :p3).value == model_param(m, :A, :p4).value == 200
@test_throws ErrorException connect_param!(m, :A, :p3, :myparam) # units error

#
# Section 2. add_shared_param! defaults
#

@defcomp A begin
    pA1 = Parameter{Symbol}()   # type will by Symbol
    pA2 = Parameter()           # type will be Number
    function run_timestep(p,v,d,t)
    end
end

@defcomp B begin
    pB1 = Parameter{Number}()   # type will be Number
    pB2 = Parameter{Int64}()    # type will be Int64
    function run_timestep(p,v,d,t)
    end
end

function _get_model()
    m = Model()
    set_dimension!(m, :time, 1:5);
    add_comp!(m, A)
    add_comp!(m, B)
    return m
end

# typical behavior
m = _get_model()
add_shared_param!(m, :myparam, 5)
@test model_param(m, :myparam) isa Mimi.ScalarModelParameter{Float64} # by default same as model, which defaults to number_type(m) == Float64

exp = :(connect_param!(m, :A, :pA1, :myparam)) # pA1 should have a specified parameter type of Symbol and !(Float64 <: Symbol)
myerr1 = try eval(exp) catch err err end
@test occursin("Mismatched datatype of parameter connection", sprint(showerror, myerr1))

connect_param!(m, :A, :pA2, :myparam) # pA2 should have a parameter type of Number by default and Float64 <: Number
connect_param!(m, :B, :pB1, :myparam) # pB1 should have a specified parameter type of Number and Float64 <: Number

exp = :(connect_param!(m, :B, :pB2, :myparam)) # pB2 should have a specified parameter type of Int64 and !(Float64 <: Int64)
myerr2 = try eval(exp) catch err err end
@test occursin("Mismatched datatype of parameter connection", sprint(showerror, myerr2))

# try data_type keyword argument
m = _get_model() # number_type(m) == Float64

exp = :(add_shared_param!(m, :myparam, :foo; data_type = Int64)) # !(:foo isa Int64)
myerr3 = try eval(exp) catch err err end
@test occursin("Mismatched datatypes:", sprint(showerror, myerr3))

add_shared_param!(m, :myparam, 5; data_type = Int64)
@test model_param(m, :myparam) isa Mimi.ScalarModelParameter{Int64} # 5 is convertible to Int64

connect_param!(m, :B, :pB2, :myparam) # pB2 should have a specified parameter type of Int64 and Int64 <: Int64
connect_param!(m, :A, :pA2, :myparam) # we allow pB2 and pA2 types to conflict as long as they both passed compatibilty with the model parameter

#
# Section 2. update_leftover_params! and set_leftover_params!
#

@defcomp A begin

    p1 = Parameter{Symbol}()
    p2 = Parameter(default = 100)
    p3 = Parameter()

    function run_timestep(p,v,d,t)
    end
end

@defcomp B begin

    p1 = Parameter{Symbol}()
    p2 = Parameter()
    p3 = Parameter()
    p4 = Parameter()
    function run_timestep(p,v,d,t)
    end
end

function _get_model()
    m = Model()
    set_dimension!(m, :time, 1:5);
    add_comp!(m, A)
    add_comp!(m, B)
    return m
end

#
# set_leftover_params!
#

m = _get_model()

# wrong type (p1 must be a Symbol)
m = _get_model()
parameters = Dict("p1" => 1, "p2" => 2, "p3" => 3, "p4" => 4)
fail_expr1 = :(set_leftover_params!(m, parameters))
err1 = try eval(fail_expr1) catch err err end
@test occursin("Cannot `convert`", sprint(showerror, err1))

# missing entry (missing p4)
m = _get_model()
parameters = Dict("p1" => :foo, "p2" => 2, "p3" => 3)
fail_expr2 = :(set_leftover_params!(m, parameters))
err2 = try eval(fail_expr2) catch err err end
@test occursin("not found in provided dictionary", sprint(showerror, err2))

# successful calls
m = _get_model()
parameters = Dict(:p1 => :foo, "p2" => 2, :p3 => 3, "p4" => 4) # keys can be Symbols or Strings
set_leftover_params!(m, parameters)
run(m)
@test m[:A, :p1] == m[:B, :p1] == :foo
@test model_param(m, :p1).is_shared

@test m[:A, :p2] == 100 # remained default value
@test !model_param(m, :A, :p2).is_shared # remains its default so is not shared

@test m[:B, :p2] == 2 # took on shared value
@test model_param(m, :p2).is_shared

@test m[:A, :p3] == m[:B, :p3] == 3
@test model_param(m, :p3).is_shared

@test m[:B, :p4] == 4
@test model_param(m, :p4).is_shared

#
# update_leftover_params!
#

# wrong type (p1 must be a Symbol)
m = _get_model()
parameters = Dict(  (:A, :p1) => 1, (:B, :p1) => 10,
                    (:B, :p2) => 20,
                    (:A, :p3) => 3, (:B, :p3) => 30,
                    (:A, :p4) => 4, (:B, :p4) => 40
                )
fail_expr3 = :(update_leftover_params!(m, parameters))
err3 = try eval(fail_expr3) catch err err end
@test occursin("Cannot `convert`", sprint(showerror, err3))

# missing entry (missing B's p4)
m = _get_model()
parameters = Dict(  (:A, :p1) => :foo, (:B, :p1) => :bar,
                    (:B, :p2) => 20,
                    (:A, :p3) => 3, (:B, :p3) => 30,
                    (:A, :p4) => 4
                )
fail_expr4 = :(update_leftover_params!(m, parameters))
err4 = try eval(fail_expr4) catch err err end
@test occursin("not found in provided dictionary", sprint(showerror, err4))

# successful calls
m = _get_model()
parameters = Dict(  (:A, :p1) => :foo, (:B, "p1") => :bar,
                    (:B, :p2) => 20,
                    (:A, :p3) => 3, (:B, :p3) => 30,
                    (:A, "p4") => 4, (:B, :p4) => 40
                )
update_leftover_params!(m, parameters)
run(m)
@test m[:A, :p1] == :foo && m[:B, :p1] == :bar
@test !model_param(m, :A, :p1).is_shared && !model_param(m, :B, :p1).is_shared
@test isnothing(model_param(m, :p1, missing_ok = true)) # no shared model parameter created

@test m[:A, :p2] == 100 # remained default value
@test !model_param(m, :A, :p2).is_shared # remains its default so is not shared

@test m[:B, :p2] == 20 # took on shared value
@test !model_param(m, :B, :p2).is_shared

@test isnothing(model_param(m, :p2, missing_ok = true)) # no shared model parameter created

@test m[:A, :p3] == 3 && m[:B, :p3] == 30
@test !model_param(m, :A, :p3).is_shared && !model_param(m, :B, :p3).is_shared
@test isnothing(model_param(m, :p3, missing_ok = true)) # no shared model parameter created

@test m[:B, :p4] == 40
@test !model_param(m, :B, :p4).is_shared
@test isnothing(model_param(m, :p4, missing_ok = true)) # no shared model parameter created

#
# Section 3. update_params!
#

@defcomp A begin

    p1 = Parameter(default = 0)
    p2 = Parameter(default = 0)
    p3 = Parameter()
    p4 = Parameter()
    p5 = Parameter()
    p6 = Parameter()

    function run_timestep(p,v,d,t)
    end
end

function _get_model()
    m = Model()
    set_dimension!(m, :time, 1:5);
    add_comp!(m, A)

    add_shared_param!(m, :shared_param, 0)
    connect_param!(m, :A, :p3, :shared_param)
    connect_param!(m, :A, :p4, :shared_param)

    return m
end

# update the shared parameters and unshared parameters separately
m = _get_model()

shared_dict = Dict(:shared_param => 1)
update_params!(m, shared_dict)

unshared_dict = Dict((:A, :p5) => 2, (:A, :p6) => 3)
update_params!(m, unshared_dict)

run(m)
@test m[:A, :p3] == m[:A, :p4] == 1
@test m[:A, :p5] == 2
@test m[:A, :p6] == 3

# update both at the same time
m = _get_model()

dict = Dict(:shared_param => 1, (:A, :p5) => 2, (:A, :p6) => 3)
update_params!(m, dict)

run(m)
@test m[:A, :p3] == m[:A, :p4] == 1
@test m[:A, :p5] == 2
@test m[:A, :p6] == 3

# test failures

m = _get_model()

shared_dict = Dict(:shared_param => :foo)
@test_throws ErrorException update_params!(m, shared_dict) # units failure
shared_dict = Dict(:p3 => 3)
@test_throws ErrorException update_params!(m, shared_dict) # can't find parameter

unshared_dict = Dict((:A, :p5) => :foo, (:A, :p6) => 3)
@test_throws MethodError update_params!(m, unshared_dict) # units failure
unshared_dict = Dict((:B, :p5) => 5)
@test_throws ErrorException update_params!(m, unshared_dict) # can't find component
unshared_dict = Dict((:B, :missing) => 5)
@test_throws ErrorException update_params!(m, unshared_dict) # can't find parameter

nothing


text:
#=
using Cuas
using Test
using ExternalBallistics

@testset "Cuas.jl" begin
    # Write your tests here.
    drone = createTargetRect(15.0,25.0,800.0)
    ϵ = RectErrorB(25.0,100.0,3.0,-10.0)
    @test pHit(drone,ϵ)== 0.0881314035058921

end
=#

ext: jl


text:
using NumericalAnalysis: SEq1
using Test

@testset "SEq1 test" begin
    @test 3.13<=SEq1.Bisection(sin, π/2, 3π/2)<=3.15
    @test 1.3<=SEq1.fixed_point(x -> x - (x^3+4x^2-10)/(3x^2+8x), 1.5)<=1.37
    @test 0.73<=SEq1.Newton(x->cos(x)-x, 0.74)<=0.74
    @test 0.73<=SEq1.Secant(x -> cos(x) - x, 0.5, π/4)<=0.74
    @test 0.73<=SEq1.FalsePos(x->cos(x)-x, 0.5, π/4)<=0.74
end


text:
# from test/criteria.jl:
losses2 = [9.5, 9.3, 10,            # 10.8     0     0     1
           9.3, 9.1, 8.9, 8,        # 11.2     0     0     2
           8.3, 8.4, 9,             # 6.02     12.5  2.08  3
           9.9, 9.5, 10,            # 21.2     25.0  1.18  4
           10.6, 10.4, 11]          # 9.61     37.5  3.90  5

stopper = EarlyStopper(PQ(alpha=3.8, k=2), NotANumber())

@test !done!(stopper, losses2[1], training=true)
@test !done!(stopper, losses2[2], training=true)
@test !done!(stopper, losses2[3])

@test !done!(stopper, losses2[4], training=true)
@test !done!(stopper, losses2[5], training=true)
@test !done!(stopper, losses2[6], training=true)
@test !done!(stopper, losses2[7])

@test !done!(stopper, losses2[8], training=true)
@test !done!(stopper, losses2[9], training=true)
@test !done!(stopper, losses2[10])

@test !done!(stopper, losses2[11], training=true)
@test !done!(stopper, losses2[12], training=true)
@test !done!(stopper, losses2[13])

@test !done!(stopper, losses2[14], training=true)
@test !done!(stopper, losses2[15], training=true)
@test done!(stopper, losses2[16])

message(stopper) == "Early stop triggered by "*
    "PQ(3.8, 2, 2.220446049250313e-16) stopping criterio
n. "

# verbose case:
stopper = EarlyStopper(NotANumber(), verbosity=1)

@test_logs (:info, r"training loss: 1.0") done!(stopper, 1.0, training=true)
@test_logs (:info, r"loss: 2.0")  done!(stopper, 2.0)


text:
type Word
	id::Int64
	value::AbstractString
	freq::Int64
	doc_freq::Int64
	used_in_docs::Vector{Any}
	function Word(id, value)
		new(id, value, 0, 0, Vector{Any}())
	end
end

typealias Document Dict{Word,Int64}
typealias DocumentList Vector{Document}
typealias DateList Vector{DateTime}

type Lexicon
	words::Vector{Word}
	string_dict::Dict{AbstractString, Int64}
	function Lexicon()
		new(Vector{Word}(), Dict{String, Word}())
	end
	function Lexicon{T<:AbstractString}(strings::Vector{T})
		l = Lexicon()
		for (i, s) in enumerate(strings)
			w = Word(i, s)
			push!(l.words,w)
			l.string_dict[s] = i
		end
		l
	end
end

type Corpus{T<:Union{DocumentList, DynamicData.DynamicDataSet{Document, DateTime}}}
	documents::T
	lexicon::Lexicon

	function Corpus()
		new(T(), CorpusUtils.Lexicon())
	end

	function Corpus(docs::T, lex::Lexicon)
		new(docs, lex)
	end
end

typealias DynamicCorpus Corpus{DynamicData.DynamicDataSet{Document, DateTime}}


text:
export Circularised, lb, ub

"""
    Circularised{K<:Kernel} <: Kernel

A kernel `k` which only accepts data on the domain `lb(k)` to `ub(k)`. This condition is not
checked at run-time though, as it would imply an unacceptably large overhead. Note that this
is a hack that doesn't work in all situations as a covariance matrix generated by it is not
necessarily positive definite. You have been warned...
"""
struct Circularised{Tk<:Kernel, Tx<:AbstractFloat} <: Kernel
    k::Tk
    lb::Tx
    ub::Tx
    d::Tx
    twod::Tx
    Circularised(k::Tk, lb::Tx, ub::Tx) where {Tk<:Kernel, Tx<:AbstractFloat} =
        new{Tk, Tx}(k, lb, ub, (ub - lb) / 2, ub - lb)
end

lb(k::Circularised) = k.lb
ub(k::Circularised) = k.ub
(k::Circularised)(x::Real, y::Real) = k.k(mod(x - y + k.d, k.twod) - k.d, 0.0) # Definitely doing slightly more work than necessary here.
==(a::Circularised, b::Circularised) = (lb(a) == lb(b)) && (ub(a) == ub(b)) && (a.k == b.k)


text:
using BinaryBuilder

name = "Kaleido"
version = v"0.2.1"

url_prefix = "https://github.com/plotly/Kaleido/releases/download/v$(version)/kaleido"
sources = [
    ArchiveSource("$(url_prefix)_linux_x64.zip", "3eaa3efd41a00db05dd71add7f54ad4e9d4259552e3f7470ac4c1a8594cea214"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)_linux_arm64.zip", "a6040fa4d95692b7047d3880391f539ea7af9353b0c7b0558ea7b8e280013703"; unpack_target = "aarch64-linux-gnu"),
    ArchiveSource("$(url_prefix)_mac_x64.zip", "57e3e1a1d98f1c25f565a1d37cbc2baf0509c13e37e6b92f7d3cb89c53b28f27"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("$(url_prefix)_mac_arm64.zip", "e5cc7bceb096db135384928d1f11b115c5697dcd1f66e7064ddb91562a61d55f"; unpack_target = "aarch64-apple-darwin20"),
    ArchiveSource("$(url_prefix)_win_x64.zip", "22857d6696fe348fea166b2950263e8c2fa41e465a1a354f2e09a6f90dfe6df3"; unpack_target = "x86_64-w64-mingw32"),
    ArchiveSource("$(url_prefix)_win_x86.zip", "771d0372e003393a313adb38a6a3fed3eaf9a9ad46b1453f92cb5489efcee376"; unpack_target = "i686-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
cp -r ${target}/* ${prefix}/.
chmod +x ${bindir}/*
if [[ "${target}" == *-mingw* ]]; then
    chmod +x ${prefix}/kaleido.cmd
fi
if [[ "${target}" != *-apple-darwin* ]]; then
    chmod +x ${bindir}/swiftshader/*.${dlext}
fi
LIC_DIR="${prefix}/share/licenses/${SRC_NAME}"
mkdir -p "${LIC_DIR}"
mv "${prefix}/LICENSE.txt" "${LIC_DIR}/LICENSE.txt"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"),
    Platform("aarch64", "linux"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

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

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

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


text:
#
println("Cfg) Test of the LandeZeeman module with ASF from an internally generated multiplet.")
#
setDefaults("print summary: open", "zzz-LandeZeeman.sum")
wa = Atomic.Computation("xx",  Nuclear.Model(26., "Fermi", 58., 3.75, AngularJ64(5//2), 1.0, 2.0);
                        properties=[JAC.LandeJ],
                        configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                        zeemanSettings=LandeZeeman.Settings(true, true, true, true, 0., true, false, Int64[] ) )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
function get_date(obj::TurnSpan)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_turn_duration(obj::TurnSpan)
    return jcall(obj, "getTurnDuration", jdouble, ())
end

function get_turn_end_date(obj::TurnSpan)
    return jcall(obj, "getTurnEndDate", AbsoluteDate, ())
end

function get_turn_start_date(obj::TurnSpan)
    return jcall(obj, "getTurnStartDate", AbsoluteDate, ())
end

function in_turn_time_range(obj::TurnSpan, arg0::AbsoluteDate)
    return jcall(obj, "inTurnTimeRange", jboolean, (AbsoluteDate,), arg0)
end

function update_end(obj::TurnSpan, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "updateEnd", void, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function update_start(obj::TurnSpan, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "updateStart", void, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
abstract type NoAlphaModel <: AlphaModel end

struct NoAlphaParam <: EoSParam
end

@newmodelsimple NoAlpha NoAlphaModel NoAlphaParam

export NoAlpha
function NoAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    model = NoAlpha(NoAlphaParam())
    return model
end

function α_function(model::CubicModel,V,T,z,alpha_model::NoAlphaModel)
   return FillArrays.Ones{Float64}(length(z))
end

is_splittable(::NoAlpha) = false


text:
using LinearAlgebra

println("Hi")

#
# Do work
#

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Test of the function `_view`.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Views" begin
    # Create the pager structure.
    a   = [(i, j) for i = 1:99, j = 1:9]
    str = sprint(show, MIME"text/plain"(), a)
    pagerd = _create_pagerd(str)

    # Initial view
    # ==========================================================================

    expected = """
99×9 Matrix{Tuple{Int64, In
 (1, 1)   (1, 2)   (1, 3)
 (2, 1)   (2, 2)   (2, 3)
 (3, 1)   (3, 2)   (3, 3)
 (4, 1)   (4, 2)   (4, 3)
 (5, 1)   (5, 2)   (5, 3)
 (6, 1)   (6, 2)   (6, 3)
 (7, 1)   (7, 2)   (7, 3)
 (8, 1)   (8, 2)   (8, 3)
 (9, 1)   (9, 2)   (9, 3)
"""

    buf = IOBuffer()
    pagerd.display_size = (11, 9*3)
    TerminalPager._view!(pagerd)
    result = String(take!(pagerd.buf.io))

    @test expected == result
    @test pagerd.lines_cropped == 90
    @test pagerd.columns_cropped == 52

    # Moving the view
    # ==========================================================================

    expected = """
2, 1)   (2, 2)   (2, 3)   (
3, 1)   (3, 2)   (3, 3)   (
4, 1)   (4, 2)   (4, 3)   (
5, 1)   (5, 2)   (5, 3)   (
6, 1)   (6, 2)   (6, 3)   (
7, 1)   (7, 2)   (7, 3)   (
8, 1)   (8, 2)   (8, 3)   (
9, 1)   (9, 2)   (9, 3)   (
10, 1)  (10, 2)  (10, 3)  (
11, 1)  (11, 2)  (11, 3)  (
"""

    buf = IOBuffer()
    pagerd.display_size = (11, 9*3)
    pagerd.start_row = 3
    pagerd.start_col = 3
    TerminalPager._view!(pagerd)
    result = String(take!(pagerd.buf.io))

    @test expected == result
    @test pagerd.lines_cropped == 88
    @test pagerd.columns_cropped == 51
end


text:
@testset "traces" begin
    @testset "Trajectory" begin
        t = Trajectory(; state = Vector{Int}(), reward = Vector{Bool}())
        @test (:state, :reward) == keys(t)
        @test haskey(t, :state)
        @test haskey(t, :reward)
        push!(t; state = 3, reward = true)
        push!(t; state = 4, reward = false)
        @test t[:state] == [3, 4]
        @test t[:reward] == [true, false]
        pop!(t)
        @test t[:state] == [3]
        @test t[:reward] == [true]
        empty!(t)
        @test t[:state] == Int[]
        @test t[:reward] == Bool[]
    end

    @testset "CircularArraySARTTrajectory" begin
        t = CircularArraySARTTrajectory(;
            capacity = 3,
            state = Vector{Int} => (4,),
            action = Int => (),
            reward = Float32 => (),
            terminal = Bool => (),
        )

        @test length(t) == 0
        push!(t; state = ones(Int, 4), action = 1)
        @test length(t) == 0
        push!(t; reward = 1.0f0, terminal = false, state = 2 * ones(Int, 4), action = 2)
        @test length(t) == 1

        @test t[:state] == hcat(ones(Int, 4), 2 * ones(Int, 4))

        push!(t; reward = 2.0f0, terminal = false, state = 3 * ones(Int, 4), action = 3)
        @test length(t) == 2

        push!(t; reward = 3.0f0, terminal = false, state = 4 * ones(Int, 4), action = 4)
        @test length(t) == 3
        @test t[:state] == [j for i in 1:4, j in 1:4]
        @test t[:reward] == [1, 2, 3]

        # test circle works as expected
        push!(t; reward = 4.0f0, terminal = true, state = 5 * ones(Int, 4), action = 5)
        @test length(t) == 3
        @test t[:state] == [j for i in 1:4, j in 2:5]
        @test t[:reward] == [2, 3, 4]
    end

    @testset "ReservoirTrajectory" begin
        # test length
        t = ReservoirTrajectory(3; a = Array{Float64,2}, b = Bool)
        push!(t; a = rand(2, 3), b = rand(Bool))
        @test length(t) == 1
        push!(t; a = rand(2, 3), b = rand(Bool))
        @test length(t) == 2
        push!(t; a = rand(2, 3), b = rand(Bool))
        @test length(t) == 3

        for _ in 1:100
            push!(t; a = rand(2, 3), b = rand(Bool))
        end

        @test length(t) == 3

        # test distribution

        Random.seed!(110)
        k, n, N = 3, 10, 10000
        stats = Dict(i => 0 for i in 1:n)
        for _ in 1:N
            t = ReservoirTrajectory(k; a = Array{Int,2}, b = Int)
            for i in 1:n
                push!(t; a = i .* ones(Int, 2, 3), b = i)
            end

            for i in 1:length(t)
                stats[t[:b][i]] += 1
            end
        end

        for v in values(stats)
            @test isapprox(v / N, k / n; atol = 0.03)
        end
    end
end


text:
"""
Shuffles the residues in each column, keeping fixed the gap positions
"""
function shuffle_residues_columnwise!(aln::Matrix{Residue})
    nseq, nres = size(aln)
    for i in 1:nres
        @inbounds for j in 1:nseq
            a = aln[j,i]
            if a != GAP
                k = rand(1:nseq)
                b = aln[k,i]
                while b == GAP
                    k = rand(1:nseq)
                    b = aln[k,i]
                end
                aln[k,i] = a
                aln[j,i] = b
            end
        end
    end
    aln
end

# 0.00084 seconds faster than an implemetation similar to shuffle_residues_columnwise (PF00085)
"""
Shuffles the residues in each sequence, keeping fixed the gap positions
"""
function shuffle_residues_sequencewise!(aln::Matrix{Residue})
    taln = transpose(aln)
    shuffle_residues_columnwise!(taln)
    transpose!(aln, taln)
end

"""
Shuffles the residues in each sequence
"""
function shuffle_sequencewise!(aln::Matrix{Residue})
    nseq, nres = size(aln)
    for i in 1:nseq
        @inbounds for j in 1:nres
            k = rand(1:nres)
            aln[i,k], aln[i,j] = aln[i,j], aln[i,k]
        end
    end
    aln
end

"""
Shuffles the residues in each column
"""
function shuffle_columnwise!(aln::Matrix{Residue})
    nseq, nres = size(aln)
    for i in 1:nres
        @inbounds for j in 1:nseq
            k = rand(1:nseq)
            aln[k,i], aln[j,i] = aln[j,i], aln[k,i]
        end
    end
    aln
end

for fun in [ :shuffle_columnwise!, :shuffle_sequencewise!, :shuffle_residues_sequencewise!, :shuffle_residues_columnwise! ]
    @eval $(fun)(aln::AbstractMultipleSequenceAlignment) = $(fun)(aln.msa)
end


text:
# ------------------------- #
# Testing normal form games #
# ------------------------- #

# NOTE: We include `@inferred` at least once for each function name. We do
#       that multiple times for the same function if we have particular reason
#       to believe there might be a type stability with that function.

using MathOptInterface
const MOI = MathOptInterface
using Clp
using CDDLib

@testset "Testing normal_form_game.jl" begin

    # Player #

    @testset "Player with 1 opponent" begin
        coordination_game_matrix = [4 0; 3 2]
        player = Player(coordination_game_matrix)

        @test @inferred(best_response(player, 2)) == 2
        @test @inferred(best_response(player, [1/2, 1/2])) == 2
        @test @inferred(best_response(player, [1/2, 1/2], BROptions())) == 2
        @test sort(@inferred(best_responses(player, [2/3, 1/3]))) ==
              sort([1, 2])
        @test best_response(
            player, [2/3, 1/3], tie_breaking=:random
            ) in [1, 2]
        @test best_response(
            player, [2/3, 1/3], BROptions(tol=1e-5, tie_breaking=:random)
            ) in [1, 2]
        @test_throws ArgumentError best_response(
            player, [2/3, 1/3], tie_breaking=:invalid_symbol
            )
        @test @inferred(is_best_response(player, 1, 1))
        @test @inferred(is_best_response(player, [1/2, 1/2], [2/3, 1/3]))

        # Perturbed best response
        @test best_response(player, [2/3, 1/3], [0., 0.1]) == 2
        @test best_response(player, [2, 1], [0., 0.1]) == 2

        # Dominated actions
        @test !(@inferred(is_dominated(player, 1)))
        @test @inferred(dominated_actions(player)) == Int[]
    end

    @testset "Player with 2 opponents" begin
        payoffs_2opponents = Array{Int}(undef, 2, 2, 2)
        payoffs_2opponents[:, 1, 1] = [3, 1]
        payoffs_2opponents[:, 1, 2] = [6, 0]
        payoffs_2opponents[:, 2, 1] = [4, 5]
        payoffs_2opponents[:, 2, 2] = [2, 7]
        player = @inferred Player(payoffs_2opponents)

        @test @inferred(payoff_vector(player, (1, 2))) == [6, 0]
        @test !(@inferred(is_best_response(player, 1, (2, 1))))
        @test @inferred(best_response(player, (2, 2))) == 2
        @test sort(@inferred(best_responses(player, ([3/7, 4/7], [1/2, 1/2])))) ==
              sort([1, 2])

        @test_throws MethodError best_response(player, (1, [1/2, 1/2]))

        @test !(@inferred(is_dominated(player, 1)))
        @test @inferred(dominated_actions(player)) == Int[]
    end

    @testset "convert for Player" begin
        T1 = Int
        T2 = Float64
        player = Player(T1[1 2; 3 4])
        N = ndims(player.payoff_array)

        for T in [T1, T2]
            player_new = Player{N,T}(player)
            @test eltype(player_new.payoff_array) == T
            @test player_new.payoff_array == player.payoff_array

            # Constructor always makes a copy
            @test player_new.payoff_array !== player.payoff_array
        end

        for T in [T1, T2]
            player_new = Player(T, player)
            @test eltype(player_new.payoff_array) == T
            @test player_new.payoff_array == player.payoff_array

            # Constructor always makes a copy
            @test player_new.payoff_array !== player.payoff_array
        end

        player_new = @inferred Player(player)
        @test eltype(player_new.payoff_array) == T1
        @test player_new.payoff_array == player.payoff_array

        # Constructor always makes a copy
        @test player_new.payoff_array !== player.payoff_array

        for T in [T1, T2]
            player_new = convert(Player{N,T}, player)
            @test eltype(player_new.payoff_array) == T
            @test player_new.payoff_array == player.payoff_array

            @test (player_new.payoff_array === player.payoff_array) ==
                  (T == T1)
        end
    end

    @testset "repr(Player)" begin
        A = [1 2; 3 4]
        player = Player(A)
        r = repr("text/plain", A)
        @test repr(player) ==
            replace(r, string(typeof(A)) =>
                       split(string(typeof(player)), ".")[end])
    end

    @testset "Tests on delete_action for Player" begin
        shapley_game = [0 1 0; 0 0 1; 1 0 0]
        player = Player(shapley_game)

        player_new_1 = @inferred delete_action(player, 1, 1)
        player_new_2 = @inferred delete_action(player, [1, 2], 1)
        @test player_new_1.payoff_array == Player([0 0 1; 1 0 0]).payoff_array
        @test player_new_2.payoff_array == Player([1 0 0]).payoff_array
    end

    # NormalFormGame #

    @testset "symmetric NormalFormGame with 2 players" begin
        coordination_game_matrix = [4 0; 3 2]
        g = @inferred(NormalFormGame(coordination_game_matrix))

        # NOTE: getindex(g, 1, 2) is equivalent to `g[1, 2]`. I use the former
        #       so we can test the type stability of the get index function
        @test @inferred(getindex(g, 1, 2)) == [0, 3]
        @test @inferred(getindex(g, CartesianIndex(1, 2))) == [0, 3]
        @test @inferred is_nash(g, (1, 1))
        @test @inferred is_nash(g, ([2/3, 1/3], [2/3, 1/3]))
    end

    @testset "asymmetric NormalFormGame with 2 players" begin
        matching_pennies_bimatrix = Array{Float64}(undef, 2, 2, 2)
        matching_pennies_bimatrix[:, 1, 1] = [1, -1]
        matching_pennies_bimatrix[:, 1, 2] = [-1, 1]
        matching_pennies_bimatrix[:, 2, 1] = [-1, 1]
        matching_pennies_bimatrix[:, 2, 2] = [1, -1]
        g = @inferred(NormalFormGame(matching_pennies_bimatrix))

        @test g[2, 1] == [-1, 1]
        @test !(is_nash(g, (1, 1)))
        @test is_nash(g, ([1/2, 1/2], [1/2, 1/2]))
    end

    @testset "asymmetric NormalFormGame with 3 players" begin
        payoffs_2opponents = Array{Int}(undef, 2, 2, 2)
        payoffs_2opponents[:, 1, 1] = [3, 1]
        payoffs_2opponents[:, 1, 2] = [6, 0]
        payoffs_2opponents[:, 2, 1] = [4, 5]
        payoffs_2opponents[:, 2, 2] = [2, 7]
        player = @inferred Player(payoffs_2opponents)
        g = @inferred NormalFormGame(tuple([player for i in 1:3]...))

        @test @inferred(getindex(g, 1, 1, 2)) == [6, 4, 1]
        @test @inferred(getindex(g, CartesianIndex(1, 1, 2))) == [6, 4, 1]
        @test @inferred is_nash(g, (1, 1, 1))
        @test @inferred !(is_nash(g, (1, 1, 2)))

        p = (1 + sqrt(65)) / 16
        @test is_nash(g, ([1-p, p], [1-p, p], [1-p, p]))
    end

    @testset "NormalFormGame input action sizes" begin
        g = @inferred NormalFormGame((2, 3, 4))

        @test @inferred(num_players(g)) == 3
        @test g.players[1].payoff_array == zeros((2, 3, 4))
        @test g.players[2].payoff_array == zeros((3, 4, 2))
        @test g.players[3].payoff_array == zeros((4, 2, 3))
    end

    @testset "NormalFormGame setindex" begin
        g = @inferred NormalFormGame((2, 2))
        g[1, 1] = [0, 10]
        g[1, 2] = [0, 10]
        g[2, 1] = [3, 5]
        g[CartesianIndex(2, 2)] = [-2, 0]

        @test g.players[1].payoff_array == [0 0; 3 -2]
        @test g.players[2].payoff_array == [10 5; 10 0]
    end

    @testset "NormalFormGame constant payoffs" begin
        g = NormalFormGame((2, 2))

        @test @inferred is_nash(g, (1, 1))
        @test is_nash(g, (1, 2))
        @test is_nash(g, (2, 1))
        @test is_nash(g, (2, 2))
    end

    @testset "convert for NormalFormGame" begin
        T1 = Int
        T2 = Float64
        players = (Player(T1[1 2 3; 4 5 6]), Player(T1[7 8; 9 10; 11 12]))
        g = NormalFormGame(players)
        N = num_players(g)

        for T in [T1, T2]
            g_new = @inferred NormalFormGame{N,T}(g)
            for (player_new, player) in zip(g_new.players, players)
                @test eltype(player_new.payoff_array) == T
                @test player_new.payoff_array == player.payoff_array

                # Constructor always makes a copy
                @test player_new.payoff_array !== player.payoff_array
            end
        end

        for T in [T1, T2]
            g_new = @inferred NormalFormGame(T, g)
            for (player_new, player) in zip(g_new.players, players)
                @test eltype(player_new.payoff_array) == T
                @test player_new.payoff_array == player.payoff_array

                # Constructor always makes a copy
                @test player_new.payoff_array !== player.payoff_array
            end
        end

        g_new = @inferred NormalFormGame(g)
        for (player_new, player) in zip(g_new.players, players)
            @test eltype(player_new.payoff_array) == T1
            @test player_new.payoff_array == player.payoff_array

            # Constructor always makes a copy
            @test player_new.payoff_array !== player.payoff_array
        end

        for T in [T1, T2]
            g_new = @inferred convert(NormalFormGame{N,T}, g)
            for (player_new, player) in zip(g_new.players, players)
                @test eltype(player_new.payoff_array) == T
                @test player_new.payoff_array == player.payoff_array

                @test (player_new.payoff_array === player.payoff_array) ===
                      (T == T1)
            end
        end
    end

    @testset "Tests on delete_action for NormalFormGame" begin
        shapley_game = Array{Int}(undef, 3, 3, 2)
        shapley_game[:, 1, 1] = [0, 0, 1]
        shapley_game[:, 2, 1] = [1, 0, 0]
        shapley_game[:, 3, 1] = [0, 1, 0]
        shapley_game[:, 1, 2] = [0, 1, 0]
        shapley_game[:, 2, 2] = [0, 0, 1]
        shapley_game[:, 3, 2] = [1, 0, 0]
        g = NormalFormGame(shapley_game)

        deleted_game_1 = Array{Int}(undef, 2, 3, 2)
        deleted_game_1[:, 1, 1] = [0, 1]
        deleted_game_1[:, 2, 1] = [0, 0]
        deleted_game_1[:, 3, 1] = [1, 0]
        deleted_game_1[:, 1, 2] = [1, 0]
        deleted_game_1[:, 2, 2] = [0, 1]
        deleted_game_1[:, 3, 2] = [0, 0]

        deleted_game_2 = Array{Int}(undef, 1, 3, 2)
        deleted_game_2[1, 1, 1] = 1
        deleted_game_2[1, 2, 1] = 0
        deleted_game_2[1, 3, 1] = 0
        deleted_game_2[1, 1, 2] = 0
        deleted_game_2[1, 2, 2] = 1
        deleted_game_2[1, 3, 2] = 0

        g_new_1 = @inferred delete_action(g, 1, 1)
        g_new_2 = @inferred delete_action(g, [1, 2], 1)
        @test g_new_1.players[1].payoff_array ==
            NormalFormGame(deleted_game_1).players[1].payoff_array
        @test g_new_1.players[2].payoff_array ==
            NormalFormGame(deleted_game_1).players[2].payoff_array
        @test g_new_2.players[1].payoff_array ==
            NormalFormGame(deleted_game_2).players[1].payoff_array
        @test g_new_2.players[2].payoff_array ==
            NormalFormGame(deleted_game_2).players[2].payoff_array
    end

    # Trivial cases with one player #

    @testset "Player with 0 opponents" begin
        payoffs = [0, 1]
        player = Player(payoffs)

        @test @inferred(payoff_vector(player, nothing)) == [0, 1]
        @test @inferred is_best_response(player, 2, nothing)
        @test @inferred(best_response(player, nothing)) == 2
        @test @inferred is_dominated(player, 1)
        @test !is_dominated(player, 2)

        payoffs = [0, 1, -1]
        player = Player(payoffs)
        dom_actions = [1, 3]
        @test @inferred(dominated_actions(player)) == dom_actions

    end

    @testset "NormalFormGame with 1 player" begin
        payoffs = [0, 1, 1]
        g = @inferred NormalFormGame(Player(payoffs))
        @test num_players(g) == 1
        @test g.players[1].payoff_array == [0, 1, 1]
        @test g[1] == 0
        @test is_nash(g, 2)
        @test !(is_nash(g, 1))
        @test is_nash(g, [0, 1/2, 1/2])

        g = NormalFormGame((2,))
        @test num_players(g) == 1
        @test g.players[1].payoff_array == zeros(2)
        g[1] = 10
        @test g.players[1].payoff_array == [10, 0]
    end

    # Invalid inputs #

    @testset "Player 0 actions" begin
        @test_throws ArgumentError p = Player(Float64[])
    end

    @testset "NormalFormGame invalid players shape inconsistent" begin
        p1 = Player(zeros((2, 3)))
        p2 = Player(zeros((2, 3)))
        @test_throws ArgumentError g = NormalFormGame((p1, p2))
    end

    @testset "NormalFormGame invalid players number inconsistent" begin
        p1 = Player(zeros((2, 2, 2)))
        p2 = Player(zeros((2, 2, 2)))
        @test_throws MethodError g = NormalFormGame((p1, p2))
    end

    @testset "NormalFormGame invalid nonsquare matrix" begin
        @test_throws ArgumentError g = NormalFormGame(zeros((2, 3)))

        A = [0, 1, 1]
        A = reshape(A, (size(A)..., 1))
        @test_throws ArgumentError g = NormalFormGame(A)
    end

    @testset "NormalFormGame invalid payoff profiles" begin
        @test_throws ArgumentError g = NormalFormGame(zeros((2, 2, 1)))
    end

    @testset "NormalFormGame empty tuple" begin
        @test_throws ArgumentError g = NormalFormGame(tuple())
    end

    @testset "NormalFormGame 0 actions" begin
        @test_throws ArgumentError g = NormalFormGame((0, 2))
    end

    @testset "payoff_vector empty tuple" begin
        p1 = Player(zeros((2, 2, 2)))
        @test_throws ArgumentError payoff_vector(p1, tuple())
    end

    @testset "is_dominated linprog error" begin
        player = Player([1. 1.; 0. -1.; -1. 0.])

        function clp_optimizer_silent_maxiter1()
            optimizer = Clp.Optimizer()
            MOI.set(optimizer, MOI.Silent(), true)
            MOI.set(optimizer, MOI.RawParameter("MaximumIterations"), 1)
            return optimizer
        end
        lp_solver = clp_optimizer_silent_maxiter1
        @test_throws ErrorException is_dominated(player, 1,
                                                 lp_solver=lp_solver)
    end

    # Utility functions #

    @testset "pure2mixed" begin
        num_actions = 3
        pure_action = 1
        mixed_action = [1., 0., 0.]
        @test @inferred(pure2mixed(num_actions, pure_action)) == mixed_action
    end

    # Pareto efficiency & Pareto dominance #

    @testset "Tests on Pareto efficiency and dominance" begin
        coordination_game_matrix = [4 0;
                                    3 2]

        equal_po_p1_bimatrix = Array{Float64}(undef, 2, 2, 2)
        equal_po_p1_bimatrix[1, 1, :] = [1, -1]
        equal_po_p1_bimatrix[1, 2, :] = [1, 1]
        equal_po_p1_bimatrix[2, 1, :] = [1, 1]
        equal_po_p1_bimatrix[2, 2, :] = [1, -1]

        three_p_equal_po_array = Array{Int}(undef, 2, 2, 2)
        three_p_equal_po_array[:, :, 1] = [2 0; 0 2]
        three_p_equal_po_array[:, :, 2] = [2 0; 0 2]

        p1 = p2 = p3 = Player(three_p_equal_po_array)

        games_dict = [NormalFormGame(coordination_game_matrix),
                      NormalFormGame(equal_po_p1_bimatrix),
                      NormalFormGame((p1, p2, p3))]

        act_prof_dict = [[(1, 1), (1, 2), (2, 1), (2, 2)],
                         [(1, 1), (1, 2), (2, 1), (2, 2)],
                         [(1, 1, 1), (2, 1, 1), (1, 2, 1), (2, 2, 1),
                          (1, 1, 2), (2, 1, 2), (1, 2, 2), (2, 2, 2)]]

        @testset "Testing is_pareto_efficient" begin
            output_dict = [[true, false, false, false],
                           [false, true, true, false],
                           [true, false, false, false, false, false, false,
                            true]]

            for i = 1:length(games_dict)
                for j in 1:length(act_prof_dict[i])
                    @test @inferred is_pareto_efficient(games_dict[i],
                                                       act_prof_dict[i][j]) ==
                                    output_dict[i][j]
                end
            end
        end

        @testset "Testing Pareto dominance" begin
            output_dict = [[true, false, false, false],
                           [false, false, false, false],
                           [false, false, false, false, false, false, false,
                            false]]

            for i = 1:length(games_dict)
                for j in 1:length(act_prof_dict[i])
                    @test @inferred is_pareto_dominant(games_dict[i],
                                                       act_prof_dict[i][j]) ==
                                    output_dict[i][j]
                end
            end
        end

        @testset "Test is_dominated" begin
            coordination_game_matrix = [4 0; 3 2]
            player = Player(coordination_game_matrix)
            for action = 1:num_actions(player)
                @test !is_dominated(player, action)
            end

            payoffs_2opponents = Array{Int}(undef, 2, 2, 2)
            payoffs_2opponents[:, 1, 1] = [3, 1]
            payoffs_2opponents[:, 1, 2] = [6, 0]
            payoffs_2opponents[:, 2, 1] = [4, 5]
            payoffs_2opponents[:, 2, 2] = [2, 7]
            player = Player(payoffs_2opponents)

            for i = 1:num_actions(player)
                @test !is_dominated(player, i)
            end

        end

        @testset "Test player corner cases" begin
            n, m = 3, 4
            player = Player(zeros((n, m)))
            for action = 1:n
                @test is_best_response(player, action, ones(m) * 1/m)
                @test !is_dominated(player, action)
            end

            e = 1e-8
            player = Player([-e -e;
                             1 -1;
                             -1 1])
            action = 1
            @test is_best_response(player, action, [1/2, 1/2], tol=e)
            @test !is_best_response(player, action, [1/2, 1/2], tol=e/2)
            @test !is_dominated(player, action, tol=e+1e-16)
            @test dominated_actions(player, tol=e+1e-16) == Int[]
            @test is_dominated(player, action, tol=e/2)
            @test dominated_actions(player, tol=e/2) == [action]
        end

        @testset "Test rational input game" begin
            T = Rational{BigInt}
            lp_solver = CDDLib.Optimizer{T}

            # Corner cases
            e = 1//(2^25)
            player = Player([-e -e;
                             1//1 -1//1;
                             -1//1 1//1])

            action = 1
            @test !is_dominated(T, player, action, tol=e, lp_solver=lp_solver)
            @test is_dominated(T, player, action, tol=e//2, lp_solver=lp_solver)

            player.payoff_array[1, 1:2] .= 0;

            @test !is_dominated(T, player, action, tol=0, lp_solver=lp_solver)

            # Simple game
            game_matrix = [2//3 1//3;
                           1//3 2//3]
            player = Player(game_matrix)

            @test dominated_actions(T, player, lp_solver=lp_solver) == Int[]
        end
    end

end


text:
typealias BondCoeff Tuple{Int, Int, Int, Int, Float64}
typealias BondCoeffs Vector{BondCoeff}

function bond_coeffs(S2::Integer)
  S = 0.5S2
  coeff = BondCoeff[]

  ## diagonal
  for si in 0:S2
    mi = si-S
    for sj in 0:S2
      mj = sj-S
      JZ = mi*mj
      if JZ == 0.0
        continue
      end
      push!(coeff, (si, si, sj, sj, JZ))
    end
  end

  ## offdiagonal
  cp = zeros(S2+1)
  cm = zeros(S2+1)
  for s in 0:S2
    m = s - S
    cp[s+1] = sqrt((S-m)*(S+m+1))
    cm[s+1] = sqrt((S+m)*(S-m+1))
  end
  for si in 0:S2
    mi = si-S

    for sj in 0:S2
      mj = sj-S

      ## S^+ S^-
      si2 = si+1
      sj2 = sj-1
      if si2 <= S2 && sj2 >= 0
        c = 0.5cp[si+1] * cm[sj+1]
        push!(coeff, (si2, si, sj2, sj, c))
      end

      ## S^- S^+
      si2 = si-1
      sj2 = sj+1
      if si2 >= 0 && sj2 <= S2
        c = 0.5cm[si+1] * cp[sj+1]
        push!(coeff, (si2, si, sj2, sj, c))
      end
    end
  end

  return coeff
end

@inline function bond(io, i::Integer,j::Integer,z::Real,xy::Real, coeffs::BondCoeffs)
  bond(io, i,j,z,complex(xy,zero(xy)), coeffs)
end

@inline function bond(io, i::Integer, j::Integer, z::Real, xy::Complex, coeffs::BondCoeffs)
  for coeff in coeffs
    s1,s2,s3,s4,c = coeff
    Js = [c*xy', complex(z*c, 0.0), c*xy]
    J = Js[s1-s2+2]
    #            i  s1 i  s2 j  s3 j  s4 J.re J.im
    @printf(io, "%d %d %d %d %d %d %d %d %llg %llg \n", i, s1, i, s2, j, s3, j, s4, J.re, J.im)
  end
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
export dtm,dtm_debug,dtm_single_threaded

# note: these functions (in cross_validation.jl have major redundancies with runmodel_multi.jl
# we may want to improve this

function dtm_single_threaded(dtmtable::DTMTable, sett::ModelSettings, cvo::CVOptions;file::String=joinpath(mktempdir(), defaultModelNameWtihCSVext))
    @info "DTM: (bk) You may want to consider the multi threaded verison of this function: simply call dtm(...)"
# if folds<0 then we consider n disjoint training sets

# currently we have three possible CVs
# k disjoint samples
#    KfoldDisjoint(length(weight),k)
# random sampling, k samples of size samplesize
#    RandomSub(length(weight),samplesize,k)
# k samples where the k validation data sets are a partition
#    Kfold(length(weight),k)

    # during CV  we will overwrite trnidx and validx!
    # store orig indices
    trnidx_orig = copy(dtmtable.trnidx)
    trnsize_orig = length(trnidx_orig)
    fullsize = length(dtmtable.weight)
    validx_orig = copy(dtmtable.validx)
    minw_orig = copy(sett.minWeight)
    seed_orig = copy(sett.seed)
    minw_prop = minw_orig / trnsize_orig

# set rnd state to make this function reporducible (irrespective of trn/val idx)
# Random.srand should not depend on sett.seed as we do not 'store' the original seed in the resulting Excel file.
	intDatahash = floor(Int, .25 * hash(2231, hash(dtmtable.features, hash(dtmtable.numerator, hash(dtmtable.denominator, hash(dtmtable.weight))))))
    Random.seed!(intDatahash)

# 1. sample Data
    size_which_is_sampled = 0
    # local cvsampler
    if cvo.use_all_data
        # this is the default case
        # trn and val datat is considered for the "sampled trn data"
        size_which_is_sampled = fullsize
    else
        size_which_is_sampled = trnsize_orig
    end
    if cvo.training_proportion > 0
        @info "DTM: training_proportion > 0 provided. Performing random sampling without replacement ($(cvo.folds) samples)"
        sample_size = max(1, round(cvo.training_proportion * size_which_is_sampled))
        @warn "DTM: RandomSub cvsampler is currently not working! Algorithm may crash."
        # need to fix line 14 of boosting where we have 'obs=obstrn+obsval' this is not meaningful if val and trn do not 'fill out' size(features,1)
        cvsampler = RandomSub(size_which_is_sampled, sample_size, cvo.folds)
    else
        if cvo.folds < 0
            @info "DTM: Performing 'DISJOINT training sets' k-fold cross validation (k=$(cvo.folds))"
            cvsampler = KfoldDisjoint(size_which_is_sampled, abs(cvo.folds))
        else
            # we have folds>=0
            @info "DTM: Performing k-fold cross validation (disjoint validation sets, k=$(cvo.folds))"
            cvsampler = Kfold(size_which_is_sampled, cvo.folds)
        end
    end

    # initialize variables
    path_and_fn_wo_extension = "some"
    path_and_fn_wo_extension, ext = splitext(file)
    header = Vector{String}(undef, 1)
    header_settings = Vector{String}(undef, 1)
    allstats = Array{Float64,2}
    allsettings = Array{Any,2}
    allmodels = Vector{Any}
    i = 1
    n_folds = abs(cvo.folds)

    # 2. run models
    for this_sample in cvsampler
        sett.seed = floor(Int, 0.31 * hash(intDatahash, hash(99812, hash(i))))
        if cvo.use_all_data
            # size_which_is_sampled=fullsize
            this_trnidx = this_sample
            dtmtable.validx = setdiff(1:size_which_is_sampled, this_sample)
        else
            # size_which_is_sampled=trnsize_orig
            # here the training data is a subset of the trn orig data
            this_trnidx = trnidx_orig[this_sample]
            # here we leave the val data as it is
            # notably trn and val "do not fill out" all data, there is third part of the data which remains unused here
        end
        dtmtable.trnidx = this_trnidx
        if sett.minWeight > 0 # if sett.minWeight is specified as a negative number then it will be interpreted as a percentage anyway (no need to update it)
            sett.minWeight = minw_prop * length(dtmtable.trnidx) # note minWeight is adjusted for each run!
        end

        path_and_fn_wo_extension_mod = string(path_and_fn_wo_extension, "_", i)
        fnmod = string(path_and_fn_wo_extension_mod, ext)
        somestrings, model = run_model_actual(dtmtable, sett, fnmod)
        desc, numbrs, desc_settingsvec, settingsvec, selectedPerformanceMeasure = get_stats(model, perfMeasure=sett.performanceMeasure)

        # add index as first column
        pushfirst!(numbrs, i)
        pushfirst!(settingsvec, string(i))
        pushfirst!(desc, "Number")
        pushfirst!(desc_settingsvec, "Number")
        if i == 1
            header = deepcopy(desc)
            header_settings = deepcopy(desc_settingsvec)
            allstats = Array{Float64,2}(length(numbrs), n_folds)
            allsettings = Array{Any,2}(length(settingsvec), n_folds)
            allmodels = Vector{Any}(n_folds)
        else
            if !all(header[2:end] .== desc[2:end])
                @warn("Model run $(i) returned an unexpected results vector:");@show desc, header
            end
            if !all(header_settings[2:end] .== desc_settingsvec[2:end])
                @warn("Model run $(i) returned an unexpected results vector:");@show desc_settingsvec, header_settings
            end
        end
        allstats[:,i] .= deepcopy(numbrs)
        allsettings[:,i] .= deepcopy(settingsvec)
        allmodels[i] = model
        i += 1
    end # this_sample in cvsampler

    # 3. aggregate some statistics
    statsdf = DataFrame(transpose(allstats), :auto)
    settsdf = DataFrame(permutedims(allsettings, [2,1]), :auto)
    DataFrames.names!(settsdf, Symbol.(header_settings))

    fld, namestr = splitdir(path_and_fn_wo_extension)
    filen = string(path_and_fn_wo_extension, "_multistats.xlsx")
    if isfile(filen)
        @info "Deleting $(filen)."
        rm(filen)
    end

    # calc some 'stats of stats'
    col_rng_stats = 2:size(statsdf, 2)
    means = map(x->Statistics.mean(statsdf[x]), col_rng_stats)
    medians = map(x->Statistics.median(statsdf[x]), col_rng_stats)
    mins = map(x->minimum(statsdf[x]), col_rng_stats)
    maxs = map(x->maximum(statsdf[x]), col_rng_stats)
    stds = map(x->StatsBase.std(statsdf[x]), col_rng_stats)
    sums = map(x->sum(statsdf[x]), col_rng_stats)

    stats_of_stats = hcat(means, medians, mins, maxs, stds, sums)
    titles = ["Mean","Median","Min","Max","Std Dev","Sum"]
    stats_of_stats = hcat(titles, transpose(stats_of_stats))
    stats_of_stats = vcat(repeat([""], 1, size(stats_of_stats, 2)), stats_of_stats)

    allstats_with_stats = vcat(transpose(allstats), stats_of_stats)
    # NOTE! statsdf is RE defined here!
    statsdf = DataFrame(allstats_with_stats, :auto)
    DataFrames.names!(statsdf, Symbol.(header))

    # define Exceldata
    sh1 = ExcelSheet("settings", settsdf)
    sh2 = ExcelSheet("stats", statsdf)
    xld = ExcelData(ExcelSheet[sh1,sh2], Array{Chart,1}(undef, 0))
    # write data
    @info "DTM: Writing multistats to file:\r\n$(filen)"
    try
		@time writeStatistics(xld, filen, true, false)
	catch e
        @show e
		@warn("DTM: Failed to create Excel Statistics file. You may want to check the PyCall installation and whether the required Python packages are installed. \r\n $(filen)")
	end
    # restore indices
    dtmtable.trnidx = copy(trnidx_orig)
    dtmtable.validx = copy(validx_orig)
    sett.minWeight = deepcopy(minw_orig)
    sett.seed = seed_orig

    return statsdf, settsdf, allmodels
end

function dtm(dtmtable::DTMTable, sett::ModelSettings, cvo::CVOptions;file::String=joinpath(mktempdir(), defaultModelNameWtihCSVext))
    # warn("This requires testing especially when Distributed.nprocs()>2")
    # if folds<0 then we consider n disjoint training sets

    # currently we have three possible CVs
    # k disjoint samples
    #    KfoldDisjoint(length(weight),k)
    # random sampling, k samples of size samplesize
    #    RandomSub(length(weight),samplesize,k)
    # k samples where the k validation data sets are a partition
    #    Kfold(length(weight),k)

    # during CV  we will overwrite trnidx and validx!
    # store orig indices
    trnidx_orig = copy(dtmtable.trnidx)
    trnsize_orig = length(trnidx_orig)
    fullsize = length(dtmtable.weight)
    validx_orig = copy(dtmtable.validx)
    minw_orig = copy(sett.minWeight)
    seed_orig = copy(sett.seed)
    minw_prop = minw_orig / trnsize_orig

    # set rnd state to make this function reporducible (irrespective of trn/val idx)
    # Random.srand should not depend on sett.seed as we do not 'store' the original seed in the resulting Excel file.
    intDatahash = floor(Int, .25 * hash(2231, hash(dtmtable.features, hash(dtmtable.numerator, hash(dtmtable.denominator, hash(dtmtable.weight))))))
    Random.seed!(intDatahash)
    # 1. sample Data
    size_which_is_sampled = 0
        # local cvsampler
    if cvo.use_all_data
            # this is the default case
            # trn and val datat is considered for the "sampled trn data"
        size_which_is_sampled = fullsize
    else
        size_which_is_sampled = trnsize_orig
    end
    if cvo.training_proportion > 0
        @info "DTM: training_proportion > 0 provided. Performing random sampling without replacement ($(cvo.folds) samples)"
        sample_size = max(1, round(cvo.training_proportion * size_which_is_sampled))
        @warn "DTM: RandomSub cvsampler is currently not working! Algorithm may crash."
            # need to fix line 14 of boosting where we have 'obs=obstrn+obsval' this is not meaningful if val and trn do not 'fill out' size(features,1)
        cvsampler = RandomSub(size_which_is_sampled, sample_size, cvo.folds)
    else
        if cvo.folds < 0
            @info "DTM: Performing 'DISJOINT training sets' k-fold cross validation (k=$(cvo.folds))"
            cvsampler = KfoldDisjoint(size_which_is_sampled, abs(cvo.folds))
        else
                # we have folds>=0
            @info "DTM: Performing k-fold cross validation (disjoint validation sets, k=$(cvo.folds))"
            cvsampler = Kfold(size_which_is_sampled, cvo.folds)
        end
    end

        # initialize variables
    path_and_fn_wo_extension = "some"
    path_and_fn_wo_extension, ext = splitext(file)
    header = Vector{String}(undef, 1)
    header_settings = Vector{String}(undef, 1)
    i = 1
    n_folds = abs(cvo.folds)
        # run an initial dummy model to get the 'structure' (ie header and size) of the output statistics
    dummySett = deepcopy(sett)
    dummySett.print_details = false;dummySett.writeTree = false;dummySett.writeSasCode = false;dummySett.writeIterationMatrix = false;dummySett.writeResult = false;dummySett.writeStatistics = false;dummySett.boolCreateZipFile = false;dummySett.writeCsharpCode = false;dummySett.writeVbaCode = false;dummySett.saveJLDFile = false;dummySett.saveResultAsJLDFile = false;dummySett.showProgressBar_time = false;dummySett.prroduceEstAndLeafMatrices = false;dummySett.write_dot_graph = false;
    path_and_fn_wo_extension_mod = string(path_and_fn_wo_extension, "_0")
    fnmod = string(path_and_fn_wo_extension_mod, ext)
    somestrings, model = run_model_actual(dtmtable, dummySett, fnmod)
    desc, numbrs, desc_settingsvec, settingsvec, selectedPerformanceMeasure = get_stats(model, perfMeasure=sett.performanceMeasure)
            # add first column
    pushfirst!(numbrs, i)
    pushfirst!(settingsvec, string(i))
    pushfirst!(desc, "Number")
    pushfirst!(desc_settingsvec, "Number")
    header = deepcopy(desc)
    header_settings = deepcopy(desc_settingsvec)
    allstats = Array{Float64,2}(undef, length(numbrs), n_folds)
    allsettings = Array{Any,2}(undef, length(settingsvec), n_folds)
    allmodels = Vector{Any}(undef, n_folds)
        # define default stats if the model fails
    defaulted_stats = deepcopy(numbrs)
    defaulted_settings = deepcopy(settingsvec)
    fill!(defaulted_stats, 0.0)

        # create a dictionary which contains all data
    di = Dict("trnidx_orig" => trnidx_orig, "ext" => ext, "minw_prop" => minw_prop, "size_which_is_sampled" => size_which_is_sampled, "defaulted_settings" => defaulted_settings, "header_settings" => header_settings, "header" => header, "cvsampler" => cvsampler, "intDatahash" => intDatahash, "sett" => deepcopy(sett), "cvo" => cvo, "dtmtable" => dtmtable, "path_and_fn_wo_extension" => path_and_fn_wo_extension, "defaulted_stats" => defaulted_stats)
        # send data to all workers: without this command we will have a lot of overhead to send the data to the processes; but the data is constant and only needs to be sent ONCE!)
    sendto_module(DecisionTrees, Distributed.workers(), local_data_dict=deepcopy(di))
    local_data_dict = di

        # run all models in parallel
        # @warn("this is currently terribly slow as the local_data_dict might be transferred to each worker for each iteration -> improve this!.... ? global const variable....?")
    pmapresult = Distributed.pmap(iLoop->run_cvsample_on_a_process(iLoop, local_data_dict), 1:length(cvsampler))

        # 2. run models
    for pmap_entry in pmapresult# this_sample in cvsampler
        numbrs = pmap_entry[1]
        settingsvec = pmap_entry[2]
        i = pmap_entry[3]
        mdl = pmap_entry[4]
            # append stats
        allstats[:,i] .= deepcopy(numbrs)
        allsettings[:,i] .= deepcopy(settingsvec)
        allmodels[i] = mdl
    end

        # 3. aggregate some statistics
    statsdf = DataFrame(transpose(allstats), :auto)
    settsdf = DataFrame(permutedims(allsettings, [2,1]), :auto)
    DataFrames.rename!(settsdf, Symbol.(header_settings))

    fld, namestr = splitdir(path_and_fn_wo_extension)
    filen = string(path_and_fn_wo_extension, "_multistats.xlsx")
    if isfile(filen)
        @info "Deleting $(filen)."
        rm(filen)
    end

        # calc some 'stats of stats'
    col_rng_stats = 2:size(statsdf, 2)
    means = map(x->Statistics.mean(statsdf[!,x]), col_rng_stats)
    medians = map(x->Statistics.median(statsdf[!,x]), col_rng_stats)
    mins = map(x->minimum(statsdf[!,x]), col_rng_stats)
    maxs = map(x->maximum(statsdf[!,x]), col_rng_stats)
    stds = map(x->StatsBase.std(statsdf[!,x]), col_rng_stats)
    sums = map(x->sum(statsdf[!,x]), col_rng_stats)

    stats_of_stats = hcat(means, medians, mins, maxs, stds, sums)
    titles = ["Mean","Median","Min","Max","Std Dev","Sum"]
    stats_of_stats = hcat(titles, transpose(stats_of_stats))
    stats_of_stats = vcat(repeat([""], 1, size(stats_of_stats, 2)), stats_of_stats)

    allstats_with_stats = vcat(transpose(allstats), stats_of_stats)
        # NOTE! statsdf is RE defined here!
    statsdf = DataFrame(allstats_with_stats, :auto)
    DataFrames.rename!(statsdf, Symbol.(header))

        # define Exceldata
    sh1 = ExcelSheet("settings", settsdf)
    sh2 = ExcelSheet("stats", statsdf)
    xld = ExcelData(ExcelSheet[sh1,sh2], Array{Chart,1}(undef, 0))
        # write data
    try
        @time writeStatistics(xld, filen, true, false)
    catch e
        @show e
        @warn("DTM: Failed to create Excel Statistics file. You may want to check the PyCall installation and whether the required Python packages are installed. \r\n $(filen)")
    end
        # restore indices
    dtmtable.trnidx = copy(trnidx_orig)
    dtmtable.validx = copy(validx_orig)
    sett.minWeight = deepcopy(minw_orig)
    sett.seed = seed_orig

    return statsdf, settsdf, allmodels
end

function run_cvsample_on_a_process(i::Int, local_data_dict::Dict)
    di = local_data_dict
    defaulted_settings = di["defaulted_settings"]
    defaulted_stats = di["defaulted_stats"]

    try
        ext = di["ext"]
        minw_prop = di["minw_prop"]
        size_which_is_sampled = di["size_which_is_sampled"]
        header_settings = di["header_settings"]
        header = di["header"]
        cvsampler = di["cvsampler"]
        intDatahash = di["intDatahash"]
        sett = di["sett"]
        cvo = di["cvo"]
        dtmtable = di["dtmtable"]
        trnidx_orig = di["trnidx_orig"]
        path_and_fn_wo_extension = di["path_and_fn_wo_extension"]

        # find ith sample
        jj = 1
        local this_sample
        for a_sample in cvsampler
            jj += 1
            if jj >= i
					this_sample = a_sample # this line is new; on 0.7alpha Julia was complainig that I use the value of the iterator variable which may not exist outside the loop in later verions (although I clearly defined this_sample outside the loop....)
                break
            end
        end
        # define trnidx
        if cvo.use_all_data
                # size_which_is_sampled=fullsize
            this_trnidx = this_sample
            dtmtable.validx = setdiff(1:size_which_is_sampled, this_sample)
        else
                # here the training data is a subset of the trn orig data
            this_trnidx = deepcopy(trnidx_orig[this_sample])
                # here we leave the val data as it is;notably trn and val "do not fill out" all data, there is third part of the data which remains unused here
        end
        dtmtable.trnidx = this_trnidx

        # update minWeight
        if sett.minWeight > 0 # if sett.minWeight is specified as a negative number then it will be interpreted as a percentage anyway (no need to update it)
            sett.minWeight = minw_prop * length(dtmtable.trnidx) # note minWeight is adjusted for each run!
        end
        # set seed
        sett.seed =  floor(Int, 0.31 * hash(intDatahash, hash(99812, hash(i))))
        # define path
        path_and_fn_wo_extension_mod = string(path_and_fn_wo_extension, "_", i)
        fnmod = string(path_and_fn_wo_extension_mod, ext)
        # run model
        somestrings, model = run_model_actual(dtmtable, sett, fnmod)
        desc, numbrs, desc_settingsvec, settingsvec, selectedPerformanceMeasure = get_stats(model, perfMeasure=sett.performanceMeasure)
        # add index as first column
        pushfirst!(numbrs, i)
        pushfirst!(settingsvec, string(i))
        pushfirst!(desc, "Number")
        pushfirst!(desc_settingsvec, "Number")
        # check if header and settings header are as expected
        if !all(header[2:end] .== desc[2:end])
            @warn("Model run $(i) returned an unexpected results vector:")
            @show desc
            @show header
        end
        if !all(header_settings[2:end] .== desc_settingsvec[2:end])
            @warn("Model run $(i) returned an unexpected results vector:")
            @show desc_settingsvec
            @show header_settings
        end
        return numbrs, settingsvec, i, model
    catch eri
        @warn("DTM: CV model $(i) failed.")
        println(eri)
        return defaulted_stats, defaulted_settings, i, deepcopy(EmtpyDTModel())
    end

    return nothing
end	 # run_cvsample_on_a_process


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.BalanceLaws: vars_state, Prognostic, Auxiliary
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates: flattenednames
using ClimateMachine.Ocean.SplitExplicit01
using ClimateMachine.GenericCallbacks
using ClimateMachine.VTK
using ClimateMachine.Checkpoint

using Test
using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.Ocean.SplitExplicit01:
    ocean_init_aux!,
    ocean_init_state!,
    ocean_boundary_state!,
    CoastlineFreeSlip,
    CoastlineNoSlip,
    OceanFloorFreeSlip,
    OceanFloorNoSlip,
    OceanSurfaceNoStressNoForcing,
    OceanSurfaceStressNoForcing,
    OceanSurfaceNoStressForcing,
    OceanSurfaceStressForcing
import ClimateMachine.DGMethods:
    update_auxiliary_state!, update_auxiliary_state_gradient!, VerticalDirection
# using GPUifyLoops

const ArrayType = ClimateMachine.array_type()

struct SimpleBox{T, BC} <: AbstractOceanProblem
    Lˣ::T
    Lʸ::T
    H::T
    τₒ::T
    λʳ::T
    θᴱ::T
    boundary_conditions::BC
end

function ocean_init_state!(p::SimpleBox, Q, A, localgeo, t)
    coords = localgeo.coord
    @inbounds y = coords[2]
    @inbounds z = coords[3]
    @inbounds H = p.H

    Q.u = @SVector [-0, -0]
    Q.η = -0
    Q.θ = (5 + 4 * cos(y * π / p.Lʸ)) * (1 + z / H)

    return nothing
end

function ocean_init_aux!(m::OceanModel, p::SimpleBox, A, geom)
    FT = eltype(A)
    @inbounds A.y = geom.coord[2]

    # not sure if this is needed but getting weird intialization stuff
    A.w = -0
    A.pkin = -0
    A.wz0 = -0
    A.u_d = @SVector [-0, -0]
    A.ΔGu = @SVector [-0, -0]

    return nothing
end

# A is Filled afer the state
function ocean_init_aux!(m::BarotropicModel, P::SimpleBox, A, geom)
    @inbounds A.y = geom.coord[2]

    A.Gᵁ = @SVector [-0, -0]
    A.U_c = @SVector [-0, -0]
    A.η_c = -0
    A.U_s = @SVector [-0, -0]
    A.η_s = -0
    A.Δu = @SVector [-0, -0]
    A.η_diag = -0
    A.Δη = -0

    return nothing
end

function main(; restart = 0)
    mpicomm = MPI.COMM_WORLD

    ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
    loglevel =
        ll == "DEBUG" ? Logging.Debug :
        ll == "WARN" ? Logging.Warn :
        ll == "ERROR" ? Logging.Error : Logging.Info
    logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
    global_logger(ConsoleLogger(logger_stream, loglevel))

    if restart == 0 && MPI.Comm_rank(mpicomm) == 0 && isdir(vtkpath)
        @info @sprintf("""Remove old dir: %s and make new one""", vtkpath)
        rm(vtkpath, recursive = true)
    end

    brickrange_2D = (xrange, yrange)
    topl_2D =
        BrickTopology(mpicomm, brickrange_2D, periodicity = (false, false))
    grid_2D = DiscontinuousSpectralElementGrid(
        topl_2D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    brickrange_3D = (xrange, yrange, zrange)
    topl_3D = StackedBrickTopology(
        mpicomm,
        brickrange_3D;
        periodicity = (false, false, false),
        boundary = ((1, 1), (1, 1), (2, 3)),
    )
    grid_3D = DiscontinuousSpectralElementGrid(
        topl_3D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    BC = (
        ClimateMachine.Ocean.SplitExplicit01.CoastlineNoSlip(),
        ClimateMachine.Ocean.SplitExplicit01.OceanFloorNoSlip(),
        ClimateMachine.Ocean.SplitExplicit01.OceanSurfaceStressForcing(),
    )
    prob = SimpleBox{FT, typeof(BC)}(Lˣ, Lʸ, H, τₒ, λʳ, θᴱ, BC)
    gravity::FT = grav(param_set)

    #- set model time-step:
    dt_fast = 240
    dt_slow = 5400
    # dt_fast = 300
    # dt_slow = 300
    if t_chkp > 0
        n_chkp = ceil(Int64, t_chkp / dt_slow)
        dt_slow = t_chkp / n_chkp
    else
        n_chkp = ceil(Int64, runTime / dt_slow)
        dt_slow = runTime / n_chkp
        n_chkp = 0
    end
    n_outp =
        t_outp > 0 ? floor(Int64, t_outp / dt_slow) :
        ceil(Int64, runTime / dt_slow)
    ivdc_dt = numImplSteps > 0 ? dt_slow / FT(numImplSteps) : dt_slow

    model = OceanModel{FT}(
        prob,
        grav = gravity,
        cʰ = cʰ,
        add_fast_substeps = add_fast_substeps,
        numImplSteps = numImplSteps,
        ivdc_dt = ivdc_dt,
        κᶜ = FT(0.1),
    )
    # model = OceanModel{FT}(prob, cʰ = cʰ, fₒ = FT(0), β = FT(0) )
    # model = OceanModel{FT}(prob, cʰ = cʰ, νʰ = FT(1e3), νᶻ = FT(1e-3) )
    # model = OceanModel{FT}(prob, cʰ = cʰ, νʰ = FT(0), fₒ = FT(0), β = FT(0) )

    barotropicmodel = BarotropicModel(model)

    minΔx = min_node_distance(grid_3D, HorizontalDirection())
    minΔz = min_node_distance(grid_3D, VerticalDirection())
    #- 2 horiz directions
    gravity_max_dT = 1 / (2 * sqrt(gravity * H) / minΔx)
    # dt_fast = minimum([gravity_max_dT])

    #- 2 horiz directions + harmonic visc or diffusion: 2^2 factor in CFL:
    viscous_max_dT = 1 / (2 * model.νʰ / minΔx^2 + model.νᶻ / minΔz^2) / 4
    diffusive_max_dT = 1 / (2 * model.κʰ / minΔx^2 + model.κᶻ / minΔz^2) / 4
    # dt_slow = minimum([diffusive_max_dT, viscous_max_dT])

    @info @sprintf(
        """Update
           Gravity Max-dT = %.1f
           Timestep       = %.1f""",
        gravity_max_dT,
        dt_fast
    )

    @info @sprintf(
        """Update
       Viscous   Max-dT = %.1f
       Diffusive Max-dT = %.1f
       Timestep      = %.1f""",
        viscous_max_dT,
        diffusive_max_dT,
        dt_slow
    )

    if restart > 0
        direction = EveryDirection()
        Q_3D, _, t0 =
            read_checkpoint(vtkpath, "baroclinic", ArrayType, mpicomm, restart)
        Q_2D, _, _ =
            read_checkpoint(vtkpath, "barotropic", ArrayType, mpicomm, restart)

        dg = OceanDGModel(
            model,
            grid_3D,
            RusanovNumericalFlux(),
            CentralNumericalFluxSecondOrder(),
            CentralNumericalFluxGradient(),
        )
        barotropic_dg = DGModel(
            barotropicmodel,
            grid_2D,
            RusanovNumericalFlux(),
            CentralNumericalFluxSecondOrder(),
            CentralNumericalFluxGradient(),
        )

        Q_3D = restart_ode_state(dg, Q_3D; init_on_cpu = true)
        Q_2D = restart_ode_state(barotropic_dg, Q_2D; init_on_cpu = true)

    else
        t0 = 0
        dg = OceanDGModel(
            model,
            grid_3D,
            # CentralNumericalFluxFirstOrder(),
            RusanovNumericalFlux(),
            CentralNumericalFluxSecondOrder(),
            CentralNumericalFluxGradient(),
        )
        barotropic_dg = DGModel(
            barotropicmodel,
            grid_2D,
            # CentralNumericalFluxFirstOrder(),
            RusanovNumericalFlux(),
            CentralNumericalFluxSecondOrder(),
            CentralNumericalFluxGradient(),
        )

        Q_3D = init_ode_state(dg, FT(0); init_on_cpu = true)
        Q_2D = init_ode_state(barotropic_dg, FT(0); init_on_cpu = true)

    end
    timeend = runTime + t0

    lsrk_ocean = LSRK54CarpenterKennedy(dg, Q_3D, dt = dt_slow, t0 = t0)
    lsrk_barotropic =
        LSRK54CarpenterKennedy(barotropic_dg, Q_2D, dt = dt_fast, t0 = t0)

    odesolver = SplitExplicitLSRK2nSolver(lsrk_ocean, lsrk_barotropic)

    #-- Set up State Check call back for config state arrays, called every ntFrq_SC time steps
    cbcs_dg = ClimateMachine.StateCheck.sccreate(
        [
            (Q_3D, "oce Q_3D"),
            (dg.state_auxiliary, "oce aux"),
            # (dg.diffstate,"oce diff",),
            # (lsrk_ocean.dQ,"oce_dQ",),
            # (dg.modeldata.tendency_dg.state_auxiliary,"tend Int aux",),
            # (dg.modeldata.conti3d_Q,"conti3d_Q",),
            (Q_2D, "baro Q_2D"),
            (barotropic_dg.state_auxiliary, "baro aux"),
        ],
        ntFrq_SC;
        prec = 12,
    )
    # (barotropic_dg.diffstate,"baro diff",),
    # (lsrk_barotropic.dQ,"baro_dQ",)
    #--

    cb_ntFrq = [n_outp, n_chkp]
    outp_nb = round(Int64, restart * n_chkp / n_outp)
    step = [outp_nb, outp_nb, restart + 1]
    cbvector = make_callbacks(
        vtkpath,
        step,
        cb_ntFrq,
        timeend,
        mpicomm,
        odesolver,
        dg,
        model,
        Q_3D,
        barotropic_dg,
        barotropicmodel,
        Q_2D,
    )

    eng0 = norm(Q_3D)
    @info @sprintf """Starting
    norm(Q₀) = %.16e
    ArrayType = %s""" eng0 ArrayType

    # slow fast state tuple
    Qvec = (slow = Q_3D, fast = Q_2D)
    # solve!(Qvec, odesolver; timeend = timeend, callbacks = cbvector)
    cbv = (cbvector..., cbcs_dg)
    solve!(Qvec, odesolver; timeend = timeend, callbacks = cbv)

    ## Enable the code block below to print table for use in reference value code
    ## reference value code sits in a file named $(@__FILE__)_refvals.jl. It is hand
    ## edited using code generated by block below when reference values are updated.
    regenRefVals = false
    if regenRefVals
        ## Print state statistics in format for use as reference values
        println(
            "# SC ========== Test number ",
            1,
            " reference values and precision match template. =======",
        )
        println("# SC ========== $(@__FILE__) test reference values ======================================")
        ClimateMachine.StateCheck.scprintref(cbcs_dg)
        println("# SC ====================================================================================")
    end

    ## Check results against reference if present
    checkRefVals = true
    if checkRefVals
        include("../refvals/simple_box_ivd_refvals.jl")
        refDat = (refVals[1], refPrecs[1])
        checkPass = ClimateMachine.StateCheck.scdocheck(cbcs_dg, refDat)
        checkPass ? checkRep = "Pass" : checkRep = "Fail"
        @info @sprintf("""Compare vs RefVals: %s""", checkRep)
        @test checkPass
    end

    return nothing
end

function make_callbacks(
    vtkpath,
    step,
    ntFrq,
    timeend,
    mpicomm,
    odesolver,
    dg_slow,
    model_slow,
    Q_slow,
    dg_fast,
    model_fast,
    Q_fast,
)
    n_outp = ntFrq[1]
    n_chkp = ntFrq[2]
    mkpath(vtkpath)
    mkpath(vtkpath * "/slow")
    mkpath(vtkpath * "/fast")

    function do_output(span, step, model, dg, Q)
        outprefix = @sprintf(
            "%s/%s/mpirank%04d_step%04d",
            vtkpath,
            span,
            MPI.Comm_rank(mpicomm),
            step
        )
        @info "doing VTK output" outprefix
        statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
        auxnames = flattenednames(vars_state(model, Auxiliary(), eltype(Q)))
        writevtk(outprefix, Q, dg, statenames, dg.state_auxiliary, auxnames)

        mycomm = Q.mpicomm
        ## Generate the pvtu file for these vtk files
        if MPI.Comm_rank(mpicomm) == 0 && MPI.Comm_size(mpicomm) > 1
            ## name of the pvtu file
            pvtuprefix = @sprintf("%s/%s/step%04d", vtkpath, span, step)
            ## name of each of the ranks vtk files
            prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
                @sprintf("mpirank%04d_step%04d", i - 1, step)
            end
            writepvtu(
                pvtuprefix,
                prefixes,
                (statenames..., auxnames...),
                eltype(Q),
            )
            @info "Done writing VTK: $pvtuprefix"
        end

    end

    do_output("slow", step[1], model_slow, dg_slow, Q_slow)
    step[1] += 1
    cbvtk_slow =
        GenericCallbacks.EveryXSimulationSteps(n_outp) do (init = false)
            do_output("slow", step[1], model_slow, dg_slow, Q_slow)
            step[1] += 1
            nothing
        end

    do_output("fast", step[2], model_fast, dg_fast, Q_fast)
    step[2] += 1
    cbvtk_fast =
        GenericCallbacks.EveryXSimulationSteps(n_outp) do (init = false)
            do_output("fast", step[2], model_fast, dg_fast, Q_fast)
            step[2] += 1
            nothing
        end

    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q_slow)
            @info @sprintf(
                """Update
                simtime = %8.2f / %8.2f
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(odesolver),
                timeend,
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end

    if n_chkp > 0
        # Note: write zeros instead of Aux vars (not needed to restart); would be
        # better just to write state vars (once write_checkpoint() can handle it)
        cb_checkpoint = GenericCallbacks.EveryXSimulationSteps(n_chkp) do
            write_checkpoint(
                Q_slow,
                zero(Q_slow),
                odesolver,
                vtkpath,
                "baroclinic",
                mpicomm,
                step[3],
            )

            write_checkpoint(
                Q_fast,
                zero(Q_fast),
                odesolver,
                vtkpath,
                "barotropic",
                mpicomm,
                step[3],
            )

            step[3] += 1
            nothing
        end
        return (cbvtk_slow, cbvtk_fast, cbinfo, cb_checkpoint)
        # return (cbinfo, cb_checkpoint)
    else
        return (cbvtk_slow, cbvtk_fast, cbinfo)
        # return (cbinfo)
    end

end

#################
# RUN THE TESTS #
#################
FT = Float64
vtkpath = "vtk_split"

const runTime = 5 * 24 * 3600 # s
const t_outp = 24 * 3600 # s
const t_chkp = runTime  # s
#const runTime = 6 * 3600 # s
#const t_outp = 6 * 3600 # s
#const t_chkp = 0
const ntFrq_SC = 1 # frequency (in time-step) for State-Check output

const N = 4
const Nˣ = 20
const Nʸ = 20
const Nᶻ = 20
const Lˣ = 4e6  # m
const Lʸ = 4e6  # m
const H = 1000  # m

xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
zrange = range(FT(-H); length = Nᶻ + 1, stop = 0)

#const cʰ = sqrt(gravity * H)
const cʰ = 1  # typical of ocean internal-wave speed
const cᶻ = 0

#- inverse ratio of additional fast time steps (for weighted average)
#  --> do 1/add more time-steps and average from: 1 - 1/add up to: 1 + 1/add
# e.g., = 1 --> 100% more ; = 2 --> 50% more ; = 3 --> 33% more ...
add_fast_substeps = 2

#- number of Implicit vertical-diffusion sub-time-steps within one model full time-step
# default = 0 : disable implicit vertical diffusion
numImplSteps = 5

const τₒ = 2e-1  # (Pa = N/m^2)
const λʳ = 20 // 86400 # m/s
#- since we are using old BC (with factor of 2), take only half:
#const τₒ = 1e-1
#const λʳ = 10 // 86400
const θᴱ = 10    # deg.C

@testset "$(@__FILE__)" begin
    main(restart = 0)
end


text:
using StatsBase
using Base.Test


text:
using KDTrees
using StatsBase
using Plotly
function run_bench_knn_points(dim, knn, exps, rounds)
    println("Running KNN benchmark for 10^(", exps, ") points in ", dim, " dimensions with k = ", knn, ".")
    n_points = [10^i for i in exps]

    times = zeros(length(n_points), rounds)
    gc_disable()
    timer = 0.0
    for (i, n_point) in enumerate(n_points)
        for (j , round) in enumerate(1:rounds)
            n_iters = 25
            print("Round ", j, " out of ", rounds, " for ", dim, "x", int(n_point), "...\r")
            data = rand(dim, int(n_point))
            tree = KDTree(data, 10)
            while true
                gc_disable()
                timer = time_ns()
                for k in 1:n_iters
                	  p = rand(dim)
                    knn(tree, p, knn)
                end
                timer = (time_ns() - float(timer)) / 10^9 # To seconds
                gc_enable()
                if timer < 1.0
                    n_iters *= 3
                    continue
                end
                break # Ends this round
            end
        times[i, j] = timer / n_iters
        end
        print("\n")
    end
    println("\nDone!")
    return mean_and_std(1./times, 2)
end

data = Dict[]
for knn in [1, 5, 10]
  dim = 3
  exps = 3:0.5:6
  rounds = 5

  speeds, stds = run_bench_knn_points(dim, knn, exps, rounds)
  speeds = vec(speeds)
  stds = vec(stds)

  stderr = stds / sqrt(length(stds))
  ymins = speeds - 1.96*stderr
  ymaxs = speeds + 1.96*stderr
  sizes = [10^i for i in exps]
  n_points = [10^i for i in exps]
  trace = [
    "x" => sizes,
    "y" => speeds,
    "type" => "scatter",
    "mode" => "lines+markers",
    "name" => string("k = ", knn),
    "error_y" => [
        "type" => "data",
        "array" => stderr*1.96*2,
        "visible" => true
    ]
  ]
  push!(data, trace)
end

####################################################################
# Plotting
####################################################################

layout = [
  "title" => "KNN search speed (dim = 3)",
  "xaxis" => [
      "title" => "Number of points",
      "type" => "log"
  ],
  "yaxis" => [
      "title" => "KNN searches / second [1/s]",
      "type" => "log"
  ],
  "type" => "log",
  "autorange" => true
]

Plotly.signin("kcarlsson89", "lolololoololololo")

response = Plotly.plot(data, ["layout" => layout,
                              "filename" => "bench_linux_knn_3",
                              "fileopt" => "overwrite"])
plot_url = response["url"]

#=
p = plot(x = sizes, y = speeds,  Geom.point, Geom.line, Geom.errorbar, Scale.x_log10, Scale.y_log10,
         ymin = ymins,
         ymax = ymaxs,
         Guide.xticks(ticks=collect(exps)),
         Guide.yticks(ticks=collect(log10(speeds))),
         Guide.xlabel("Number of points"),
         Guide.ylabel("KNN searches / second [1/s]"),
         Guide.title("KNN search speed (dim = 3, k = 5)"))
g = render(p)
img = PNG(string("knn_plot_", knn, ".png"), 10inch, 8inch)
draw(img, g)
=#


text:
include("constraint_struct.jl")
include("settings_struct.jl")


text:
module PkgServerCacheWarmer

include("defaults.jl")
include("types.jl")

include("commands.jl")
include("public.jl")
include("temporary_environment.jl")
include("warmup_cache.jl")

end # module


text:
using NNlib: glu
using Zygote

@testset "glu" begin
    x = [1 2 3; 4 5 6; 7 8 9; 10 11 12]
    @test ceil.(glu(x, 1)) == [1 2 3; 4 5 6]
    @test_throws AssertionError glu(x, 2)
end

@testset "AutoDiff" begin
    local rng = StableRNG(17)
    gradtest(glu, rand(rng, 4, 3))
end


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

function withlocales(f, newlocales)
    # save current locales
    locales = Dict{Int,String}()
    for cat in 0:9999
        cstr = ccall(:setlocale, Cstring, (Cint, Cstring), cat, C_NULL)
        if cstr != C_NULL
            locales[cat] = unsafe_string(cstr)
        end
    end
    timestrs = String[]
    try
        # change to each of given locales
        for lc in newlocales
            set = true
            for (cat, _) in locales
                set &= ccall(:setlocale, Cstring, (Cint, Cstring), cat, lc) != C_NULL
            end
            set && f()
        end
    finally
        # recover locales
        for (cat, lc) in locales
            cstr = ccall(:setlocale, Cstring, (Cint, Cstring), cat, lc)
        end
    end
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
"""
Naive approach to calculating forces (scales very poorly, O(n^2)).
    Takes as input an array of floats X with dimensions (n,d) for n particles
    in d dimensions

    Returns a vector of forces F (n,d) when collisions are detected.
"""
function Fcol(X,radius=1.0)
    F = zeros(size(X,1),2)
    r2 = 4.0*radius*radius # distance threshold
    for i in 1:size(X,1)
        for j in i+1:size(X,1) # use symmetry
            if j != i
                rij = X[j,1:2] .- X[i,1:2] # relative dist from i to j
                d2 = rij[1]*rij[1]+rij[2]*rij[2] # can avoid one sqrt
                if d2 < r2
                    d = sqrt(d2)
                    F[i,:] += (d-2*radius)*rij/d
                    F[j,:] -= (d-2*radius)*rij/d # use symmetry
                end
            end
        end
    end
    return F
end

"""
    Takes as input the positions of N particles X ~ (N,2), and returns
    a new positions in the next time frame.
"""
function Step(X;radius=1.0,v0=1.0,dt=0.001,Dr=0.1,Dt=0.0,L=10,k=300.)
    Y = copy(X)
    F = Fcol(X,radius) # get forces
    for i in 1:size(X,1)
        d = [cos(X[i,3]),sin(X[i,3])] # i's directions vector
        Y[i,3] += sqrt(2.0*Dr/dt)*randn()*dt # increment of rotational noise
        # self propulsion + forces + spatial diffusion
        Y[i,1:2] += dt*(v0*d + k*F[i,:] + sqrt(2.0*Dt/dt)*randn(2))
    end
    # periodic box
    Y[:,1:2] = mod.(Y[:,1:2],L)
    return Y
end

max_stars_repo_name: Jerboa-app/ActiveBrownianParticles

max_issues_repo_name: Jerboa-app/ActiveBrownianParticles

max_forks_repo_name: Jerboa-app/ActiveBrownianParticles


text:
using TheNatureOfComputation
using Documenter

makedocs(; modules=[TheNatureOfComputation], sitename="TheNatureOfComputation", authors="Beacon Biosignals, Inc.",
         pages=["API Documentation" => "index.md"])

deploydocs(; repo="github.com/beacon-biosignals/TheNatureOfComputation.jl.git", push_preview=true)