text:
mutable struct Tree
    state
    action
    value
    count
    depth
    terminal
    children
    parent
end

EPSILON = 0.01

function transition(state, action)

    next_state = nothing
    reward = nothing
    done = nothing

    # next_state
    if state < 0.5 - EPSILON && action == 1
        next_state = 0.6 * rand()  # [0, 0.6]
        reward = 1.
    elseif state < 0.5 - EPSILON && action == 2
        next_state = 0.75 * rand()  # [0, 0.75]
        reward = 0.
    elseif state > 0.5 + EPSILON && action == 2
        next_state = 0.4 + 0.5 * rand()  # [0.4, 1.0]
        reward = 1.
    elseif state > 0.5 + EPSILON && action == 1
        next_state = 0.25 + 0.75 * rand()  # [0.25, 1.0]
        reward = 0.
    end

    # done
    if next_state > 0.5 - EPSILON && next_state < 0.5 + EPSILON
        done = true
    else
        done = false
    end

    # return
    return next_state, reward, done
end

ACTIONS = [1, 2]

default_policy(state) = rand(ACTIONS)

function score(node, alpha = 1.0)
    avg_value = node.value / node.count
    rel_freq = 2 * alpha * sqrt(2 * log(node.parent.count) / node.count)
    return avg_value + rel_freq
end

function uct_search(init_state, max_depth = 20, max_steps = 100)
    root = Tree(init_state, nothing, 0.0, 0, 0, false, Array{Tree,1}(), nothing)
    steps = 0
    node = root
    while node.depth < max_depth && steps < max_steps
        value = evaluate(node, default_policy)
        backup(node, value)
        node = tree_policy(root)
        steps += 1
        println("step = $steps, depth = $(node.depth)")
    end
    return argmax([n.value for n in root.children])
end

function tree_policy(node)
    """Find the next leaf to evaluate.

        If the current node is expandable, return a node for an untried action. Else, choose the "best" child node and repeat until you arrive at a terminal node.
    """
    while ~node.terminal
        child = expand(node)
        if child == nothing
            node = best_child(node)
        else
            return child
        end
    end
    return node
end

function expand(node)
    available_actions = setdiff(ACTIONS, Set([n.action for n in node.children]))
    length(available_actions) == 0 && return nothing  # node *not* expandable
    action = rand(available_actions)
    state, _, done = transition(node.state, action)
    child = Tree(state, action, 0.0, 0, node.depth + 1, done, Array{Tree,1}(), node)
    push!(node.children, child)
    return child
end

best_child(node) = node.children[argmax(score.(node.children))]

function evaluate(node, policy)
    """Evaluate a node by running `policy` until end of episode."""
    value = 0.0
    state = node.state
    done = node.terminal
    while ~done
        action = policy(state)
        state, reward, done = transition(state, action)
        value += reward
    end
    return value
end

function backup(node, value)
    """Update value of each node by adding `value` and incrementing count."""
    while node != nothing
        node.count += 1
        node.value += value
        node = node.parent
    end
end

init_state = rand()
opt_action = uct_search(init_state)


text:
# ------------------------------------------------------------------------------------------
# # Julia 很快
#
# 基准数值测试程序（benchmarks）通常被用来比较编程语言的性能。
# 这些 benchmarks 可能引起很多的讨论，首先是作为基准测试的到底是什么，再有就是究竟是什么导致不同语言间的差异。
# 这些简单的问题有时比你乍想之下复杂得多。
#
# 这个 notebook 目的在于让你实现一个简单的基准测试。
# 通过这个 notebook 可以看到作者的配有四核 Macbook Pro 表现如何，或者你可以自己跑跑看。
#
# （材料来自于 MIT 的 Steven Johnson 的一场演讲：[18S096/Boxes-and-registers.ipynb at master ·
# mitmath/18S096](https://github.com/mitmath/18S096/blob/master/lectures/lecture1/Boxes-and-
# registers.ipynb)）
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 目录
#
# - [`sum` 函数的定义](#sum:-一个简单易懂的函数)
# - [不同语言中的各种基准测试](#不同语言中的各种基准测试)
#     1. [C (手写)](#1.-C-语言)
#     2. [C (手写 + `-ffast-math`)](#2.-使用--ffast-math-的-C-语言)
#     3. [python (内置)](#3.-Python-的内置-sum)
#     4. [python (numpy)](#4.-Python:-numpy)
#     5. [python (手写)](#5.-Python（手写）)
#     6. [Julia (内置)](#6.-Julia（内置）)
#     7. [Julia (手写)](#7.-Julia（手写）)
#     8. [Julia (手写 + SIMD)](#8.-Julia（手写-+--@simd）)
# - [性能测试小结](#性能测试小结)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## `sum`: 一个简单易懂的函数
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 考虑**求和**函数 `sum(a)`，它是用来计算
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i
# $$
# 其中 $n$ 是 `a` 的长度。
# ------------------------------------------------------------------------------------------

a = rand(10^7)  # 在[0,1)上均匀分布的随机数组成的1维向量

sum(a)

# ------------------------------------------------------------------------------------------
# 预期结果应该在 $0.5 \times 10^7$ 左右，因为所有元素的平均值应该是 0.5 左右
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 不同语言中的各种基准测试
# ------------------------------------------------------------------------------------------

@time sum(a)

@time sum(a)

@time sum(a)

# ------------------------------------------------------------------------------------------
# `@time` 宏会产生受干扰的结果，所以不是做基准测试的最佳选择！
#
# 好在 Julia 有个 `BenchmarkTools.jl` 包来简单并准确地进行基准测试：
# ------------------------------------------------------------------------------------------

using Pkg
Pkg.add("BenchmarkTools")

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# ###  1. C 语言
#
# C 语言对人类不友好但是对机器友好，经常作为一个黄金标准。
# 在 C 语言的两倍以内通常来说就比较理想了。
# 然而，即使仅考虑 C 语言，不同的 C 语言程序员也可能采用或不采用不同的优化方式。
#
# 原作者不会 C 语言，所以他看不懂下面这个代码块中的代码，
# 但是他知道 Julia 中可以编译、运行 C 语言代码。注意 `"""` 所包围的是个多行字符串。
# ------------------------------------------------------------------------------------------

using Libdl
C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()  # 创建一个临时文件

# 通过将 C_code 导入 gcc 来编译成共享链接库
# (只有装了 gcc 才行):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# 定义一个调用这个 C 语言函数的 Julia 函数：
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a)  # 输入\approx后按<TAB>键得到≈符号

c_sum(a) - sum(a)

≈  # `isapprox` 函数的别名

?isapprox

# ------------------------------------------------------------------------------------------
# 现在我们可以在 Julia 中直接对 C 语言代码进行基准测试了：
# ------------------------------------------------------------------------------------------

c_bench = @benchmark c_sum($a)

println("C: 最快 $(minimum(c_bench.times) / 1e6) 毫秒")

d = Dict()  # 字典，也就是关联数组（键值对）
d["C"] = minimum(c_bench.times) / 1e6  # 毫秒
d

using Plots
gr()

using Statistics  # 导入统计包以计算标准差
t = c_bench.times / 1e6  # 毫秒
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ------------------------------------------------------------------------------------------
# ### 2. 使用 `-ffast-math` 的 C 语言
#
# 如果我们允许 C 重新安排浮点操作，它就会使用 SIMD 指令(single instruction, multiple data)进行矢量化。
# ------------------------------------------------------------------------------------------

const Clib_fastmath = tempname()   # 创建临时文件

# 同上但是增加了 -ffast-math 标志
open(`gcc -fPIC -O3 -msse3 -xc -shared -ffast-math -o $(Clib_fastmath * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# 定义一个调用这个 C 语言函数的 Julia 函数：
c_sum_fastmath(X::Array{Float64}) = ccall(("c_sum", Clib_fastmath), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_fastmath_bench = @benchmark $c_sum_fastmath($a)

d["C -ffast-math"] = minimum(c_fastmath_bench.times) / 1e6  # 毫秒

# ------------------------------------------------------------------------------------------
# ### 3. Python 的内置 `sum`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `PyCall` 包给 Julia 提供了 Python 接口：
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("PyCall")
using PyCall

# 获得 Python 的内置 "sum" 函数：
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($a)

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 4. Python: `numpy`
#
# **条件允许的情况下利用硬件 "SIMD"**
#
# `numpy` 是一个 Python 的优化过的 C 语言库。
# 在 Julia 版 Python 中的安装方法为：
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("Conda")
using Conda

Conda.add("numpy")

numpy_sum = pyimport("numpy")["sum"]

py_numpy_bench = @benchmark $numpy_sum($a)

numpy_sum(a)

numpy_sum(a) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 5. Python（手写）
# ------------------------------------------------------------------------------------------

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

py_hand = @benchmark $sum_py($a)

sum_py(a)

sum_py(a) ≈ sum(a)

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 6. Julia（内置）
#
# **是直接用 Julia 写的，不是 C！**
# ------------------------------------------------------------------------------------------

@which sum(a)

j_bench = @benchmark sum($a)

d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 7. Julia（手写）
# ------------------------------------------------------------------------------------------

function mysum(A)
    s = 0.0  # s = zero(eltype(a))
    for a in A
        s += a
    end
    s
end

j_bench_hand = @benchmark mysum($a)

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 8. Julia（手写 +  `@simd`）
# ------------------------------------------------------------------------------------------

function mysum_simd(A)
    s = 0.0  # s = zero(eltype(A))
    @simd for a in A
        s += a
    end
    s
end

j_bench_hand_simd = @benchmark mysum_simd($a)

mysum_simd(a)

d["Julia hand-written simd"] = minimum(j_bench_hand_simd.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ## 性能测试小结
# ------------------------------------------------------------------------------------------

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value; digits=1), 6, "."))
end

ext: jl

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
# Robust Statistics

#############################
#
#   Trimming outliers
#
#############################

# Trimmed set
"Return the upper and lower bound elements used by `trim` and `winsor`"
function uplo(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    n = length(x)
    n > 0 || throw(ArgumentError("x can not be empty."))

    if count == 0
        0 <= prop < 0.5 || throw(ArgumentError("prop must satisfy 0 ≤ prop < 0.5."))
        count = floor(Int, n * prop)
    else
        prop == 0 || throw(ArgumentError("prop and count can not both be > 0."))
        0 <= count < n/2 || throw(ArgumentError("count must satisfy 0 ≤ count < length(x)/2."))
    end

    # indices for lowest count values
    x2 = copy(x)
    lo = partialsort!(x2, 1:count+1)[end]
    # indices for largest count values
    up = partialsort!(x2, n-count:n)[1]

    up, lo
end

"""
    trim(x::AbstractVector; prop=0.0, count=0)

Return an iterator of all elements of `x` that omits either `count` or proportion
`prop` of the highest and lowest elements.

The number of trimmed elements could be smaller than specified if several
elements equal the lower or upper bound.

To compute the trimmed mean of `x` use `mean(trim(x))`;
to compute the variance use `trimvar(x)` (see [`trimvar`](@ref)).

# Example
```julia
julia> collect(trim([5,2,4,3,1], prop=0.2))
3-element Array{Int64,1}:
 2
 4
 3
```
"""
function trim(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    up, lo = uplo(x; prop=prop, count=count)

    (xi for xi in x if lo <= xi <= up)
end

"""
    trim!(x::AbstractVector; prop=0.0, count=0)

A variant of [`trim`](@ref) that modifies `x` in place.
"""
function trim!(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    up, lo = uplo(x; prop=prop, count=count)
    ix = (i for (i,xi) in enumerate(x) if lo > xi || xi > up)
    deleteat!(x, ix)
    return x
end

"""
    winsor(x::AbstractVector; prop=0.0, count=0)

Return an iterator of all elements of `x` that replaces either `count` or
proportion `prop` of the highest elements with the previous-highest element
and an equal number of the lowest elements with the next-lowest element.

The number of replaced elements could be smaller than specified if several
elements equal the lower or upper bound.

To compute the Winsorized mean of `x` use `mean(winsor(x))`.

# Example
```julia
julia> collect(winsor([5,2,3,4,1], prop=0.2))
5-element Array{Int64,1}:
 4
 2
 3
 4
 2
```
"""
function winsor(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    up, lo = uplo(x; prop=prop, count=count)

    (clamp(xi, lo, up) for xi in x)
end

"""
    winsor!(x::AbstractVector; prop=0.0, count=0)

A variant of [`winsor`](@ref) that modifies vector `x` in place.
"""
function winsor!(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    copyto!(x, winsor(x; prop=prop, count=count))
    return x
end

#############################
#
#   Other
#
#############################

# Variance of a trimmed set.
"""
    trimvar(x; prop=0.0, count=0)

Compute the variance of the trimmed mean of `x`. This function uses
the Winsorized variance, as described in Wilcox (2010).
"""
function trimvar(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    n = length(x)
    n > 0 || throw(ArgumentError("x can not be empty."))

    if count == 0
        0 <= prop < 0.5 || throw(ArgumentError("prop must satisfy 0 ≤ prop < 0.5."))
        count = floor(Int, n * prop)
    else
        0 <= count < n/2 || throw(ArgumentError("count must satisfy 0 ≤ count < length(x)/2."))
        prop = count/n
    end

    return var(winsor(x, count=count)) / (n * (1 - 2prop)^2)
end


text:
import ..Filter
import ..Wirreducible  #remove all polarization subdiagrams
import ..Girreducible  #remove all self-energy inseration
import ..NoHatree
import ..NoFock
import ..NoBubble  # true to remove all bubble subdiagram
import ..Proper  #ver4, ver3, and polarization diagrams may require to be irreducible along the transfer momentum/frequency

import ..DiagramType
import ..GreenDiag
import ..SigmaDiag
import ..PolarDiag
import ..Ver3Diag
import ..Ver4Diag
import ..GnDiag
import ..GcDiag

import ..Composite
import ..ChargeCharge
import ..SpinSpin
import ..UpUp
import ..UpDown
import ..Response

import ..Instant
import ..Dynamic
import ..D_Instant
import ..D_Dynamic
import ..AnalyticProperty

import ..symbol
import ..short

import ..Interaction

import ..GenericPara

import ..innerTauNum

# unique id
# uid() = abs(rand(Int)) % 10000

# let z = 0
#     global function uid()
#         z += 1
#         return z
#     end
# end
_counter = 0
function uid()
    global _counter += 1
    return _counter
end

function uidreset()
    global _counter = 0
end

function getK(loopNum, loopIdx)
    k = zeros(loopNum)
    k[loopIdx] = 1.0
    return k
end


text:
function getKelleyPrimalCuts(thePortfolio::SparsePortfolioData, useCopyOfVariables::Bool, minReturnConstraint::Bool, theStabilizationPoint::Array{Float64}, kelleyPrimalEpochs::Int64)

    rootModel = Model(solver=CplexSolver(CPX_PARAM_SCRIND=0))
    @variable(rootModel, sRoot[1:thePortfolio.n]>=0.0)
    @constraint(rootModel, sRoot.<=ones(thePortfolio.n))
    @variable(rootModel, tRoot>=-1e12)

    # Objective
    @objective(rootModel, Min, tRoot)

    # Constraints
    @constraint(rootModel, sum(sRoot)<=thePortfolio.k)
    @constraint(rootModel, sum(sRoot)>=1)

    # Max investment constraint
    if useCopyOfVariables
      @constraint(rootModel, sum(sRoot[i]*(thePortfolio.A[1,i].>thePortfolio.l[1]) for i=1:thePortfolio.n)>=1.0) #ensures feasibility of problem
      @variable(rootModel,   xRoot[1:thePortfolio.n]>=0.0)
      @constraint(rootModel, thePortfolio.A*xRoot.>=thePortfolio.l)
      @constraint(rootModel, thePortfolio.A*xRoot.<=thePortfolio.u)
      @constraint(rootModel, minInvestment[i=1:thePortfolio.n], xRoot[i]>=sRoot[i]*thePortfolio.min_investment[i])
      @constraint(rootModel, bigMMaster[i=1:thePortfolio.n], xRoot[i]<=thePortfolio.u[i+1].*sRoot[i])
      @constraint(rootModel, sumsToOne, sum(xRoot[i] for i=1:thePortfolio.n)==1.0)
    end

    if minReturnConstraint #This line assumes that minimum return constraint comes first whenever it is included
      @constraint(rootModel, sum(sRoot[i]*(thePortfolio.A[1,i].>thePortfolio.l[1]) for i=1:thePortfolio.n)>=1.0) #ensures feasibility of problem
    end

    UB = Inf; LB = -Inf
    rootStabilizationTrick = :inOut
    stabilizationPoint=theStabilizationPoint
    rootCutsSense=1
    ε= 1e-10
    λ = (rootStabilizationTrick == :inOut) ? .1 : 1.
    δ = (rootStabilizationTrick == :inOut || rootStabilizationTrick == :twoEps) ? 2*ε : 0.
    rootCutsLim=kelleyPrimalEpochs
    rootCutCount = 0
    oaRootCutCount=0
    consecutiveNonImprov_1 = 0
    consecutiveNonImprov_2 = 0
    theCutPool=Cut[]
    for epoch in 1:kelleyPrimalEpochs
      #@show epoch, LB, UB, oaRootCutCount
      solve(rootModel)
      zstar = min.(max.(getvalue(sRoot), 0.), 1.) #cope with numerical errors
      stabilizationPoint += zstar; stabilizationPoint /= 2

      if LB >= getobjectivevalue(rootModel) - eps()
        if consecutiveNonImprov_1 == 5
          consecutiveNonImprov_2 += 1
        else
          consecutiveNonImprov_1 += 1
        end
      else
        if consecutiveNonImprov_1 < 5
          consecutiveNonImprov_1 = 0
        end
        consecutiveNonImprov_2 = 0
      end
      LB = max(LB, getobjectivevalue(rootModel))

      if consecutiveNonImprov_1 == 5
        λ = 1
      elseif consecutiveNonImprov_2 == 5
        δ = 0.
      end

      z0 = λ*zstar + (1-λ)*stabilizationPoint .+ δ
      z0=min.(ones(thePortfolio.n), max.(z0, zeros(thePortfolio.n)))

      Cut  = portfolios_objective2(thePortfolio, z0)
      addOACut=true
      if Cut.status==:Optimal
        UB=min(UB, Cut.p)
      else
        z0 = zstar
        z0=min.(ones(thePortfolio.n), max.(z0, zeros(thePortfolio.n)))
        @suppress Cut  = portfolios_objective2(thePortfolio, z0)
        addOACut=false
      end

      if Cut.status==:Optimal
        @constraint(rootModel, tRoot >= Cut.p + dot(Cut.∇s,sRoot-z0))
        if ( (rootCutsSense > 0)&(epoch <= rootCutsLim) )||( (rootCutsSense < 0)&(epoch >= kelleyPrimalEpochs-rootCutsLim+1) ) & addOACut
            Cut.p+=dot(Cut.∇s,-z0)
            push!(theCutPool, Cut)
        end
      else
        @constraint(rootModel, sum(z0[i]*(1.0-sRoot[i])+sRoot[i]*(1.0-z0[i]) for i=1:thePortfolio.n)>=1.0)
        stabilizationPoint=theStabilizationPoint
      end

      if abs(UB-LB)/abs(UB) <= 1e-4 || consecutiveNonImprov_2 >= 10
        #println("Root node gap: ", abs(UB-LB))
      end
    end
    #println("Kelley Primal: Generated $kelleyPrimalEpochs cuts")

  return theCutPool
end


text:
#!/usr/local/bin/julia
#
# This script demonstrates how to submit multiple OrdinalGWAS runs from multiple sets of
# Plink files on UCLA Hoffman2 cluster. It assumes that a demo data is available by
# running `julia cluster_preparedata.jl` at current folder.
#

using OrdinalGWAS, Serialization

# Step 1: fit null model and save result to file `fittednullmodel.jls`
nm = ordinalgwas(@formula(trait ~ sex), "covariate.txt", nothing)
open("fittednullmodel.jls", "w") do io
    Serialization.serialize(io, nm)
end

# Step 2: GWAS for each chromosome
for chr in 1:23
    println("submit job for chromosome=$chr")
    jcode = "using OrdinalGWAS, Serialization;
    nm = open(deserialize, \"fittednullmodel.jls\");
    bedfile = \"hapmap3.chr.\" * string($chr) * \".bed.gz\";
    bimfile = \"hapmap3.chr.\" * string($chr) * \".bim.gz\";
    pvalfile = \"hapmap3.chr.\" * string($chr) * \".pval.txt\";
    ordinalgwas(nm, bedfile, bimfile, 324; pvalfile=pvalfile);"
    # prepare sh file for qsub
    open("tmp.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l h_rt=0:30:00,h_data=2G") # request runtime and memory
        println(io, "#\$ -pe shared 2") # request # shared-memory nodes
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.0.1") # available Julia version
        println(io)
        println(io, "# run julia code")
        println(io, "julia -e '$jcode' > output.\$JOB_ID 2>&1")
    end
    # submit job
    run(`qsub tmp.sh`)
end

max_stars_repo_name: Hua-Zhou/PolrGWAS

max_issues_repo_name: Hua-Zhou/PolrGWAS


text:
abstract type PointTriangulationPartition{D, T} <: AbstractTriangulationPartition{D, T} end
abstract type PointTriangulationPartitionFull{D, T} <: AbstractTriangulationPartitionFull{D, T} end

abstract type MutablePointTriangulationPartition{D, T} <: AbstractMutableTriangulationPartition{D, T} end
abstract type MutablePointTriangulationPartitionFull{D, T} <: AbstractMutableTriangulationPartitionFull{D, T} end

# Minimal
struct PointTriangulation{D, T} <: PointTriangulationPartition{D, T}
	data::Vector{Vector{T}}
	simplexindices::DelaunayTriangulation
end

mutable struct MutablePointTriangulation{D, T} <: MutablePointTriangulationPartition{D, T}
	data::Vector{Vector{T}}
	simplexindices::DelaunayTriangulation
end

# Full info
struct PointTriangulationFull{D, T} <: PointTriangulationPartitionFull{D, T}
	data::Vector{Vector{T}}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

# Full info
struct MutablePointTriangulationFull{D, T} <: MutablePointTriangulationPartitionFull{D, T}
	data::Vector{Vector{T}}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

PT = Union{PointTriangulationPartition, PointTriangulationPartitionFull,
    MutablePointTriangulationPartition, MutablePointTriangulationPartitionFull}
getpoints(pt::PT) = [pt.data[i] for i = 1:length(pt.data)]

getsimplices(pt::Union{PointTriangulationPartitionFull, MutablePointTriangulationPartitionFull}) =
    pt.simplices
getsimplices(pt::Union{PointTriangulationPartition, MutablePointTriangulationPartition}) =
    [Simplex(pt.data[pt.simplexindices[i]]) for i = 1:length(pt.simplexindices)]

####################
# Pretty printing
####################
function summarise(tp::Union{PointTriangulation, MutablePointTriangulation})
    _type = typeof(tp)
    n_simplices = nsimplices(tp)
    n_pts = length(tp.data)
    D = dimension(tp)
    summary = "$_type: $n_simplices simplices constructed from $n_pts points"
end

Base.show(io::IO, tp::Union{PointTriangulation, MutablePointTriangulation}) =
    println(io, summarise(tp))

export
PointTriangulation, MutablePointTriangulation,
PointTriangulationFull, MutablePointTriangulationFull


text:
using Formatting
using ArgParse
using DataStructures, JSON

s = ArgParseSettings()
@add_arg_table s begin

    "--casename"
        help = "Casename"
        arg_type = String
        required = true

    "--year-rng"
        help = "How many years?"
        arg_type = Int64
        nargs = 2
        required = true

end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

function runOneCmd(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

function pleaseRun(cmd)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i])
        end
    else
        runOneCmd(cmd)
    end
end

casename = parsed["casename"]
beg_year = parsed["year-rng"][1]
end_year = parsed["year-rng"][2]

println("Begin year: $(beg_year)")
println("End   year: $(end_year)")

layers = 33
in_dir = "/glade/u/home/tienyiao/scratch-tienyiao/archive/$(casename)/ocn/hist"

coord = "z_t,z_w,z_w_top,z_w_bot,TAREA"
varnames = "HMXL,TEMP,SALT,TAUX,TAUY,SHF,SHF_QSW,SFWF"
year_rng      = format( "{:04d}-{:04d}",  beg_year, end_year )
year_rng_eval = format( "{:04d}..{:04d}",  beg_year, end_year )

ref_file="$(in_dir)/$(casename).pop.h.$(format("{:04d}", beg_year))-01.nc"

out_dir="./output_$(year_rng)_layers$(layers)"
coord_file="$(out_dir)/coord.nc"
time_file="$(out_dir)/time.nc"

if isdir(out_dir)
    throw(ErrorException("ERROR: directory $(out_dir) already exists."))
end

#pleaseRun(`rm -rf $(out_dir)`)

pleaseRun(`mkdir -p $(out_dir)`)
pleaseRun(`julia make_monthly_time.jl --output $(time_file) --years $( end_year - beg_year + 1 )`)

println("Output directory: $(out_dir)")

pleaseRun(`ncks -O -F -v $coord -d z_t,1,$(layers) -d z_w_top,1,$(layers) -d z_w_bot,1,$(layers) -d z_w,1,$(layers) $ref_file $coord_file`)
pleaseRun(`ncap2 -O -s 'z_t=-z_t/100.0;z_w_top=-z_w_top/100.0;z_w_bot=-z_w_bot/100.0;z_w=-z_w/100.0' $coord_file $coord_file`)

output_files = Dict()

for varname in split(varnames, ",")
    output_file = "$(out_dir)/$(varname).nc"
    output_files[varname] = output_file

    println("Concatnating var: $varname to $(output_file)")

    if varname != "HMXL"
        pleaseRun(`bash -c "ncrcat -O -F -v $(varname) -d z_t,1,$(layers) -d z_w_top,1,$(layers) -d z_w_bot,1,$(layers) $(in_dir)/$(casename).pop.h.{$(year_rng_eval)}-{01..12}.nc $(output_file)"`)
    else

        # HMXL was output in daily. So we have to average it to monthly

        mkpath("tmp_HMXL")
        for y=beg_year:end_year, m=1:12
            timestr = format("{:04d}-{:02d}", y, m)
            pleaseRun(`ncra -O -F -v HMXL $(in_dir)/$(casename).pop.h.nday1.$(timestr)-01.nc tmp_HMXL/HMXL_$(timestr).nc`)
        end
        pleaseRun(`bash -c "ncrcat -O -F -v HMXL tmp_HMXL/HMXL_{$(year_rng_eval)}-{01..12}.nc $(output_file)"`)

    end

end

println("Converting forcing for IOM...")

for varname in ["VSFLX", "HMXL", "SWFLX", "NSWFLX"]
    output_files[varname] = "$(out_dir)/$(varname).nc"
end

pleaseRun(`ncap2 -O -v -s 'VSFLX=SFWF;'      $(output_files["SFWF"]) $(output_files["VSFLX"])`)
pleaseRun(`ncap2 -O -v -s 'HMXL=HMXL/100.0;' $(output_files["HMXL"]) $(output_files["HMXL"])`)
pleaseRun(`ncap2 -O -v -s 'TAUX=TAUX/10.0;'  $(output_files["TAUX"]) $(output_files["TAUX"])`)
pleaseRun(`ncap2 -O -v -s 'TAUY=TAUY/10.0;'  $(output_files["TAUY"]) $(output_files["TAUY"])`)
pleaseRun(`ncap2 -O -v -s 'SWFLX=-SHF_QSW;'  $(output_files["SHF_QSW"]) $(output_files["SWFLX"])`)

pleaseRun(`ncks  -O    -v SHF                      $(output_files["SHF"])     $(output_files["NSWFLX"])`)
pleaseRun(`ncks     -A -v SHF_QSW                  $(output_files["SHF_QSW"]) $(output_files["NSWFLX"])`)
pleaseRun(`ncap2 -O -v -s 'NSWFLX=-(SHF-SHF_QSW);' $(output_files["NSWFLX"])  $(output_files["NSWFLX"])`)

for (_, output_file) in output_files
    println(coord_file, "; ", output_file)
    pleaseRun(`ncks -A -v $coord         $coord_file $output_file`)
    pleaseRun(`ncks -A -v time,time_bound $time_file $output_file`)
end

println("Doing monthly mean:")

out_dir_monthly = "./$(out_dir)/monthly"
tmp_dir = "$(out_dir_monthly)/tmp"
mkpath(tmp_dir)

a_year_time_file = "$(tmp_dir)/time.nc"

pleaseRun(`julia make_monthly_time.jl --output $(a_year_time_file) --years 1`)

output_files_monthly = Dict()
for varname in keys(output_files)

    println("Doing variable $(varname)")

    output_files_monthly[varname] = "$(out_dir_monthly)/$(varname)_monthly.nc"

    for m = 1:12
        pleaseRun(`ncra -O -F -d time,$(m),,12 $(output_files[varname]) $tmp_dir/$(varname)_$(format("{:02d}", m)).nc`)
    end

    pleaseRun(`bash -c "ncrcat -O -F $tmp_dir/$(varname)_{01..12}.nc $(output_files_monthly[varname])"`)
    pleaseRun(`ncks -A -v time,time_bound $a_year_time_file $(output_files_monthly[varname])`)

end

rm(tmp_dir, recursive=true, force=true)

println("All done.")

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function bisect_right(data::Array{Int64}, item::Int64)::Int64
    low = 0
    high = length(data)
    while low < high
        middle = i32::from(((low + high) / 2))
        if item < data[middle]
            high = middle
        else

            low = (middle + 1)
        end
    end
    return low
end

function bin_it(limits::Array{Int64}, data::Array{Int64})::Array{Int64}
    bins = [0]
    for x in limits
        bins.push(0)
    end
    for d in data
        bins[bisect_right(limits, d)] += 1
    end
    return bins
end


text:
GridLayout(; kwargs...) = GridLayout(1, 1; kwargs...)

GridLayout(scene::Scene, args...; kwargs...) = GridLayout(args...; bbox = lift(x -> BBox(x), pixelarea(scene)), parentscene = scene, kwargs...)

function GridLayout(nrows::Int, ncols::Int;
        rowsizes = nothing,
        colsizes = nothing,
        addedrowgaps = nothing,
        addedcolgaps = nothing,
        alignmode = Inside(),
        equalprotrusiongaps = (false, false),
        bbox = nothing,
        parentscene = nothing,
        kwargs...)

    if isnothing(rowsizes)
        rowsizes = [Auto() for _ in 1:nrows]
    # duplicate a single row size into a vector for every row
    elseif rowsizes isa ContentSize
        rowsizes = [rowsizes for _ in 1:nrows]
    elseif !(typeof(rowsizes) <: Vector{<: ContentSize})
        error("Row sizes must be one size or a vector of sizes, not $(typeof(rowsizes))")
    end

    if isnothing(colsizes)
        colsizes = [Auto() for _ in 1:ncols]
    # duplicate a single col size into a vector for every col
    elseif colsizes isa ContentSize
        colsizes = [colsizes for _ in 1:ncols]
    elseif !(typeof(colsizes) <: Vector{<: ContentSize})
        error("Column sizes must be one size or a vector of sizes, not $(typeof(colsizes))")
    end

    if isnothing(addedrowgaps)
        addedrowgaps = [Fixed(20) for _ in 1:nrows-1]
    elseif addedrowgaps isa GapSize
        addedrowgaps = [addedrowgaps for _ in 1:nrows-1]
    elseif !(typeof(addedrowgaps) <: Vector{<: GapSize})
        error("Row gaps must be one size or a vector of sizes, not $(typeof(addedrowgaps))")
    end

    if isnothing(addedcolgaps)
        addedcolgaps = [Fixed(20) for _ in 1:ncols-1]
    elseif addedcolgaps isa GapSize
        addedcolgaps = [addedcolgaps for _ in 1:ncols-1]
    elseif !(typeof(addedcolgaps) <: Vector{<: GapSize})
        error("Column gaps must be one size or a vector of sizes, not $(typeof(addedcolgaps))")
    end

    needs_update = Node(true)

    content = []

    attrs = merge!(Attributes(kwargs), default_attributes(GridLayout))

    sizeattrs = sizenode!(attrs.width, attrs.height)

    alignment = lift(tuple, attrs.halign, attrs.valign)

    autosizenode = Node{NTuple{2, Optional{Float32}}}((nothing, nothing))

    computedsize = computedsizenode!(sizeattrs, autosizenode)

    suggestedbbox = create_suggested_bboxnode(bbox)

    finalbbox = alignedbboxnode!(suggestedbbox, computedsize, alignment,
        sizeattrs, autosizenode)

    protrusions = Node(RectSides(0f0, 0f0, 0f0, 0f0))

    layoutnodes = LayoutNodes{GridLayout, GridLayout}(suggestedbbox, protrusions, computedsize, autosizenode, finalbbox, nothing)

    gl = GridLayout(
        content, nrows, ncols, rowsizes, colsizes, addedrowgaps,
        addedcolgaps, alignmode, equalprotrusiongaps, needs_update, layoutnodes, attrs, parentscene)

    on(finalbbox) do bb
        align_to_bbox!(gl, bb)
    end

    on(needs_update) do u
        # TODO: is this correct? or should the bbox change somehow when a member
        # size changes

        w = determinedirsize(gl, Col())
        h = determinedirsize(gl, Row())

        new_autosize = (w, h)
        new_protrusions = RectSides{Float32}(
            protrusion(gl, Left()),
            protrusion(gl, Right()),
            protrusion(gl, Bottom()),
            protrusion(gl, Top()),
        )

        if autosizenode[] == new_autosize &&
                gl.layoutnodes.protrusions[] == new_protrusions

            gl.layoutnodes.suggestedbbox[] = gl.layoutnodes.suggestedbbox[]
        else
            # otherwise these values will not already be up to date when adding the
            # gridlayout into the next one

            # TODO: this is a double update?
            gl.layoutnodes.protrusions[] = new_protrusions
            autosizenode[] = new_autosize

            if isnothing(gl.layoutnodes.gridcontent)
                gl.layoutnodes.suggestedbbox[] = gl.layoutnodes.suggestedbbox[]
            end
        end

        nothing
    end

    gl
end

computedsizenode(gridlayout::GridLayout) = gridlayout.layoutnodes.computedsize
protrusionnode(gridlayout::GridLayout) = gridlayout.layoutnodes.protrusions

function validategridlayout(gl::GridLayout)
    if gl.nrows < 1
        error("Number of rows can't be smaller than 1")
    end
    if gl.ncols < 1
        error("Number of columns can't be smaller than 1")
    end

    if length(gl.rowsizes) != gl.nrows
        error("There are $nrows rows but $(length(rowsizes)) row sizes.")
    end
    if length(gl.colsizes) != gl.ncols
        error("There are $ncols columns but $(length(colsizes)) column sizes.")
    end
    if length(gl.addedrowgaps) != gl.nrows - 1
        error("There are $nrows rows but $(length(addedrowgaps)) row gaps.")
    end
    if length(gl.addedcolgaps) != gl.ncols - 1
        error("There are $ncols columns but $(length(addedcolgaps)) column gaps.")
    end
end

function with_updates_suspended(f::Function, gl::GridLayout)
    gl.block_updates = true
    f()
    gl.block_updates = false
    gl.needs_update[] = true
end

function connect_layoutnodes!(gc::GridContent)

    disconnect_layoutnodes!(gc::GridContent)

    gc.protrusions_handle = on(protrusionnode(gc.content)) do p
        gc.needs_update[] = true
    end
    gc.computedsize_handle = on(computedsizenode(gc.content)) do c
        gc.needs_update[] = true
    end
end

function disconnect_layoutnodes!(gc::GridContent)
    if !isnothing(gc.protrusions_handle)
        Observables.off(protrusionnode(gc.content), gc.protrusions_handle)
        gc.protrusions_handle = nothing
    end
    if !isnothing(gc.computedsize_handle)
        Observables.off(computedsizenode(gc.content), gc.computedsize_handle)
        gc.computedsize_handle = nothing
    end
end

function add_to_gridlayout!(g::GridLayout, gc::GridContent)

    # to be safe
    remove_from_gridlayout!(gc)

    push!(g.content, gc)

    # let the gridcontent know that it's inside a gridlayout
    gc.parent = g

    on(gc.needs_update) do update
        g.needs_update[] = true
    end

    # trigger relayout
    g.needs_update[] = true
end

function remove_from_gridlayout!(gc::GridContent)
    if isnothing(gc.parent)
        return
    end

    i = findfirst(x -> x === gc, gc.parent.content)
    if isnothing(i)
        error("""GridContent had a parent but was not in the content array.
        This must be a bug.""")
    end
    deleteat!(gc.parent.content, i)

    gc.parent = nothing

    # remove all listeners from needs_update because they could be pointing
    # to previous parents if we're re-nesting layout objects
    empty!(gc.needs_update.listeners)
end

function convert_contentsizes(n, sizes)::Vector{ContentSize}
    if isnothing(sizes)
        [Auto() for _ in 1:n]
    elseif sizes isa ContentSize
        [sizes for _ in 1:n]
    elseif sizes isa Vector{<:ContentSize}
        length(sizes) == n ? sizes : error("$(length(sizes)) sizes instead of $n")
    else
        error("Illegal sizes value $sizes")
    end
end

function convert_gapsizes(n, gaps)::Vector{GapSize}
    if isnothing(gaps)
        [Fixed(20) for _ in 1:n]
    elseif gaps isa GapSize
        [gaps for _ in 1:n]
    elseif gaps isa Vector{<:GapSize}
        length(gaps) == n ? gaps : error("$(length(gaps)) gaps instead of $n")
    else
        error("Illegal gaps value $gaps")
    end
end

function appendrows!(gl::GridLayout, n::Int; rowsizes=nothing, addedrowgaps=nothing)

    rowsizes = convert_contentsizes(n, rowsizes)
    addedrowgaps = convert_gapsizes(n, addedrowgaps)

    with_updates_suspended(gl) do
        gl.nrows += n
        append!(gl.rowsizes, rowsizes)
        append!(gl.addedrowgaps, addedrowgaps)
    end
end

function appendcols!(gl::GridLayout, n::Int; colsizes=nothing, addedcolgaps=nothing)

    colsizes = convert_contentsizes(n, colsizes)
    addedcolgaps = convert_gapsizes(n, addedcolgaps)

    with_updates_suspended(gl) do
        gl.ncols += n
        append!(gl.colsizes, colsizes)
        append!(gl.addedcolgaps, addedcolgaps)
    end
end

function prependrows!(gl::GridLayout, n::Int; rowsizes=nothing, addedrowgaps=nothing)

    rowsizes = convert_contentsizes(n, rowsizes)
    addedrowgaps = convert_gapsizes(n, addedrowgaps)

    newcontent = map(gl.content) do gc
        span = gc.span
        newspan = Span(span.rows .+ n, span.cols)
        GridContent(gc.content, newspan, gc.side)
    end
    # because of too specific typing with the new parametric GridContent
    gl.content = convert(Vector{GridContent}, newcontent)

    with_updates_suspended(gl) do
        gl.nrows += n
        prepend!(gl.rowsizes, rowsizes)
        prepend!(gl.addedrowgaps, addedrowgaps)
    end
end

function prependcols!(gl::GridLayout, n::Int; colsizes=nothing, addedcolgaps=nothing)

    colsizes = convert_contentsizes(n, colsizes)
    addedcolgaps = convert_gapsizes(n, addedcolgaps)

    newcontent = map(gl.content) do gc
        span = gc.span
        newspan = Span(span.rows, span.cols .+ n)
        GridContent(gc.content, newspan, gc.side)
    end
    # because of too specific typing with the new parametric GridContent
    gl.content = convert(Vector{GridContent}, newcontent)

    with_updates_suspended(gl) do
        gl.ncols += n
        prepend!(gl.colsizes, colsizes)
        prepend!(gl.addedcolgaps, addedcolgaps)
    end
end

function deleterow!(gl::GridLayout, irow::Int)
    if !(1 <= irow <= gl.nrows)
        error("Row $irow does not exist.")
    end

    if gl.nrows == 1
        error("Can't delete the last row")
    end

    # new_content = GridContent[]
    to_remove = GridContent[]
    for c in gl.content
        rows = c.span.rows
        newrows = if irow in rows
            # range is one shorter now
            rows.start : rows.stop - 1
        elseif irow > rows.stop
            # content before deleted row stays the same
            rows
        else
            # content completely after is moved forward 1 step
            rows .- 1
        end
        if isempty(newrows)
            # the row span was just one row and now zero, remove the element
            push!(to_remove, c)
        else
            c.span = Span(newrows, c.span.cols)
        end
    end

    for c in to_remove
        remove_from_gridlayout!(c)
    end
    # gl.content = new_content
    deleteat!(gl.rowsizes, irow)
    deleteat!(gl.addedrowgaps, irow == 1 ? 1 : irow - 1)
    gl.nrows -= 1
    gl.needs_update[] = true
end

function deletecol!(gl::GridLayout, icol::Int)
    if !(1 <= icol <= gl.ncols)
        error("Col $icol does not exist.")
    end

    if gl.ncols == 1
        error("Can't delete the last col")
    end

    to_remove = GridContent[]
    for c in gl.content
        cols = c.span.cols
        newcols = if icol in cols
            # range is one shorter now
            cols.start : cols.stop - 1
        elseif icol > cols.stop
            # content before deleted col stays the same
            cols
        else
            # content completely after is moved forward 1 step
            cols .- 1
        end
        if isempty(newcols)
            # the col span was just one col and now zero, remove the element
            push!(to_remove, c)
        else
            c.span = Span(c.span.rows, newcols)
        end
    end

    for c in to_remove
        remove_from_gridlayout!(c)
    end
    # gl.content = new_content
    deleteat!(gl.colsizes, icol)
    deleteat!(gl.addedcolgaps, icol == 1 ? 1 : icol - 1)
    gl.ncols -= 1
    gl.needs_update[] = true
end

function Base.isempty(gl::GridLayout, dir::GridDir, i::Int)
    !any(gl.content) do c
        span = dir isa Row ? c.span.rows : c.span.cols
        i in span
    end
end

function trim!(gl::GridLayout)
    irow = 1
    while irow <= gl.nrows && gl.nrows > 1
        if isempty(gl, Row(), irow)
            deleterow!(gl, irow)
        else
            irow += 1
        end
    end

    icol = 1
    while icol <= gl.ncols && gl.ncols > 1
        if isempty(gl, Col(), icol)
            deletecol!(gl, icol)
        else
            icol += 1
        end
    end
end

function gridnest!(gl::GridLayout, rows::Indexables, cols::Indexables)

    newrows, newcols = adjust_rows_cols!(gl, rows, cols)

    subgl = GridLayout(
        length(newrows), length(newcols);
        parent = nothing,
        colsizes = gl.colsizes[newcols],
        rowsizes = gl.rowsizes[newrows],
        addedrowgaps = gl.addedrowgaps[newrows.start:(newrows.stop-1)],
        addedcolgaps = gl.addedcolgaps[newcols.start:(newcols.stop-1)],
    )

    # remove the content from the parent that is completely inside the replacement grid
    subgl.block_updates = true
    i = 1
    while i <= length(gl.content)
        gc = gl.content[i]

        if (gc.span.rows.start >= newrows.start && gc.span.rows.stop <= newrows.stop &&
            gc.span.cols.start >= newcols.start && gc.span.cols.stop <= newcols.stop)

            subgl[gc.span.rows .- (newrows.start - 1), gc.span.cols .- (newcols.start - 1), gc.side] = gc.content
            continue
            # don't advance i because there's one piece of content less in the queue
            # and the next item is in the same position as the old removed one
        end
        i += 1
    end
    subgl.block_updates = false

    gl[newrows, newcols] = subgl

    subgl
end

function find_in_grid(obj, container::GridLayout)
    for i in 1:length(container.content)
        if container.content[i].content === obj
            return i
        end
    end
    nothing
end

function topmost_grid(gl::GridLayout)
    candidate = gl
    while true
        if candidate.parent isa GridLayout
            candidate = candidate.parent
        else
            return candidate
        end
    end
end

function find_in_grid_and_subgrids(obj, container::GridLayout)
    for i in 1:length(container.content)
        candidate = container.content[i].content

        if candidate === obj
            return container, i
        elseif candidate isa GridLayout
            return find_in_grid_and_subgrids(obj, candidate)
        end
    end
    nothing, nothing
end

function find_in_grid_tree(obj, container::GridLayout)
    topmost = topmost_grid(container)
    find_in_grid_and_subgrids(obj, topmost)
end

function detachfromgridlayout!(obj, gl::GridLayout)
    i = find_in_grid(obj, gl)
    if !isnothing(i)
        deleteat!(gl.content, i)
    end
end

function Base.show(io::IO, ::MIME"text/plain", gl::GridLayout)

    function spaceindent(str, n, downconnection)
        joinstr = if downconnection
            "\n" * (" " ^ 1) * "┃" * (" " ^ (n-2))
        else
            "\n" * (" " ^ n)
        end
        join(split(str, "\n"), joinstr)
    end

    println(io, "GridLayout[$(gl.nrows), $(gl.ncols)] with $(length(gl.content)) children")

    for (i, c) in enumerate(gl.content)
        rows = c.span.rows
        cols = c.span.cols
        content = c.content

        connector = i == length(gl.content) ? " ┗━ " : " ┣━ "

        if content isa GridLayout
            downconnection = i < length(gl.content)
            str = spaceindent(repr(MIME"text/plain"(), content), 2, downconnection)
            println(io, connector * "[$rows | $cols] $str")
        else
            println(io, connector * "[$rows | $cols] $(typeof(content))")
        end
    end
end

function Base.show(io::IO, gl::GridLayout)
    print(io, "GridLayout[$(gl.nrows), $(gl.ncols)] ($(length(gl.content)) children)")
end

function colsize!(gl::GridLayout, i::Int, s::ContentSize)
    if !(1 <= i <= gl.ncols)
        error("Can't set size of invalid column $i.")
    end
    gl.colsizes[i] = s
    gl.needs_update[] = true
end

function rowsize!(gl::GridLayout, i::Int, s::ContentSize)
    if !(1 <= i <= gl.nrows)
        error("Can't set size of invalid row $i.")
    end
    gl.rowsizes[i] = s
    gl.needs_update[] = true
end

function colgap!(gl::GridLayout, i::Int, s::GapSize)
    if !(1 <= i <= (gl.ncols - 1))
        error("Can't set size of invalid column gap $i.")
    end
    gl.addedcolgaps[i] = s
    gl.needs_update[] = true
end

function colgap!(gl::GridLayout, s::GapSize)
    gl.addedcolgaps .= s
    gl.needs_update[] = true
end

function colgap!(gl::GridLayout, r::Real)
    gl.addedcolgaps .= Ref(Fixed(r))
    gl.needs_update[] = true
end

function rowgap!(gl::GridLayout, i::Int, s::GapSize)
    if !(1 <= i <= (gl.nrows - 1))
        error("Can't set size of invalid row gap $i.")
    end
    gl.addedrowgaps[i] = s
    gl.needs_update[] = true
end

function rowgap!(gl::GridLayout, s::GapSize)
    gl.addedrowgaps .= s
    gl.needs_update[] = true
end

function rowgap!(gl::GridLayout, r::Real)
    gl.addedrowgaps .= Ref(Fixed(r))
    gl.needs_update[] = true
end

"""
This function solves a grid layout such that the "important lines" fit exactly
into a given bounding box. This means that the protrusions of all objects inside
the grid are not taken into account. This is needed if the grid is itself placed
inside another grid.
"""
function align_to_bbox!(gl::GridLayout, suggestedbbox::BBox)

    # compute the actual bbox for the content given that there might be outside
    # padding that needs to be removed
    bbox = if gl.alignmode isa Outside
        pad = gl.alignmode.padding
        BBox(
            left(suggestedbbox) + pad.left,
            right(suggestedbbox) - pad.right,
            bottom(suggestedbbox) + pad.bottom,
            top(suggestedbbox) - pad.top)
    elseif gl.alignmode isa Inside
        suggestedbbox
    elseif gl.alignmode isa Mixed
        pad = gl.alignmode.padding
        BBox(
            left(suggestedbbox) + ifnothing(pad.left, 0f0),
            right(suggestedbbox) - ifnothing(pad.right, 0f0),
            bottom(suggestedbbox) + ifnothing(pad.bottom, 0f0),
            top(suggestedbbox) - ifnothing(pad.top, 0f0))
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode))")
    end

    # first determine how big the protrusions on each side of all columns and rows are
    maxgrid = RowCols(gl.ncols, gl.nrows)
    # go through all the layout objects placed in the grid
    for c in gl.content
        idx_rect = side_indices(c)
        mapsides(idx_rect, maxgrid) do side, idx, grid
            grid[idx] = max(grid[idx], protrusion(c, side))
        end
    end

    # for the outside alignmode
    topprot = maxgrid.tops[1]
    bottomprot = maxgrid.bottoms[end]
    leftprot = maxgrid.lefts[1]
    rightprot = maxgrid.rights[end]

    # compute what size the gaps between rows and columns need to be
    colgaps = maxgrid.lefts[2:end] .+ maxgrid.rights[1:end-1]
    rowgaps = maxgrid.tops[2:end] .+ maxgrid.bottoms[1:end-1]

    # determine the biggest gap
    # using the biggest gap size for all gaps will make the layout more even, but one
    # could make this aspect customizable, because it might waste space
    if gl.equalprotrusiongaps[2]
        colgaps = ones(gl.ncols - 1) .* (gl.ncols <= 1 ? 0.0 : maximum(colgaps))
    end
    if gl.equalprotrusiongaps[1]
        rowgaps = ones(gl.nrows - 1) .* (gl.nrows <= 1 ? 0.0 : maximum(rowgaps))
    end

    # determine the vertical and horizontal space needed just for the gaps
    # again, the gaps are what the protrusions stick into, so they are not actually "empty"
    # depending on what sticks out of the plots
    sumcolgaps = (gl.ncols <= 1) ? 0.0 : sum(colgaps)
    sumrowgaps = (gl.nrows <= 1) ? 0.0 : sum(rowgaps)

    # compute what space remains for the inner parts of the plots
    remaininghorizontalspace = if gl.alignmode isa Inside
        width(bbox) - sumcolgaps
    elseif gl.alignmode isa Outside
        width(bbox) - sumcolgaps - leftprot - rightprot
    elseif gl.alignmode isa Mixed
        rightal = getside(gl.alignmode, Right())
        leftal = getside(gl.alignmode, Left())
        width(bbox) - sumcolgaps -
            (isnothing(leftal) ? 0 : leftprot) -
            (isnothing(rightal) ? 0 : rightprot)
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode))")
    end

    remainingverticalspace = if gl.alignmode isa Inside
        height(bbox) - sumrowgaps
    elseif gl.alignmode isa Outside
        height(bbox) - sumrowgaps - topprot - bottomprot
    elseif gl.alignmode isa Mixed
        topal = getside(gl.alignmode, Top())
        bottomal = getside(gl.alignmode, Bottom())
        height(bbox) - sumrowgaps -
            (isnothing(bottomal) ? 0 : bottomprot) -
            (isnothing(topal) ? 0 : topprot)
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode))")
    end

    # compute how much gap to add, in case e.g. labels are too close together
    # this is given as a fraction of the space used for the inner parts of the plots
    # so far, but maybe this should just be an absolute pixel value so it doesn't change
    # when resizing the window
    addedcolgaps = map(gl.addedcolgaps) do cg
        if cg isa Fixed
            return cg.x
        elseif cg isa Relative
            return cg.x * remaininghorizontalspace
        else
            return 0.0 # for float type inference
        end
    end
    addedrowgaps = map(gl.addedrowgaps) do rg
        if rg isa Fixed
            return rg.x
        elseif rg isa Relative
            return rg.x * remainingverticalspace
        else
            return 0.0 # for float type inference
        end
    end

    # compute the actual space available for the rows and columns (plots without protrusions)
    spaceforcolumns = remaininghorizontalspace - ((gl.ncols <= 1) ? 0.0 : sum(addedcolgaps))
    spaceforrows = remainingverticalspace - ((gl.nrows <= 1) ? 0.0 : sum(addedrowgaps))

    colwidths, rowheights = compute_col_row_sizes(spaceforcolumns, spaceforrows, gl)

    # don't allow smaller widths than 1 px even if it breaks the layout (better than weird glitches)
    colwidths = max.(colwidths, ones(length(colwidths)))
    rowheights = max.(rowheights, ones(length(rowheights)))

    # this is the vertical / horizontal space between the inner lines of all plots
    finalcolgaps = colgaps .+ addedcolgaps
    finalrowgaps = rowgaps .+ addedrowgaps

    # compute the resulting width and height of the gridlayout and compute
    # adjustments for the grid's alignment (this will only matter if the grid is
    # bigger or smaller than the bounding box it occupies)

    gridwidth = sum(colwidths) + sum(finalcolgaps) +
        (gl.alignmode isa Outside ? (leftprot + rightprot) : 0.0)
    gridheight = sum(rowheights) + sum(finalrowgaps) +
        (gl.alignmode isa Outside ? (topprot + bottomprot) : 0.0)

    # compute the x values for all left and right column boundaries
    xleftcols = if gl.alignmode isa Inside
        left(bbox) .+ cumsum([0; colwidths[1:end-1]]) .+
            cumsum([0; finalcolgaps])
    elseif gl.alignmode isa Outside
        left(bbox) .+ cumsum([0; colwidths[1:end-1]]) .+
            cumsum([0; finalcolgaps]) .+ leftprot
    elseif gl.alignmode isa Mixed
        leftal = getside(gl.alignmode, Left())
        left(bbox) .+ cumsum([0; colwidths[1:end-1]]) .+
            cumsum([0; finalcolgaps]) .+ (isnothing(leftal) ? 0 : leftprot)
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode))")
    end
    xrightcols = xleftcols .+ colwidths

    # compute the y values for all top and bottom row boundaries
    ytoprows = if gl.alignmode isa Inside
        top(bbox) .- cumsum([0; rowheights[1:end-1]]) .-
            cumsum([0; finalrowgaps])
    elseif gl.alignmode isa Outside
        top(bbox) .- cumsum([0; rowheights[1:end-1]]) .-
            cumsum([0; finalrowgaps]) .- topprot
    elseif gl.alignmode isa Mixed
        topal = getside(gl.alignmode, Top())
        top(bbox) .- cumsum([0; rowheights[1:end-1]]) .-
            cumsum([0; finalrowgaps]) .- (isnothing(topal) ? 0 : topprot)
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode))")
    end
    ybottomrows = ytoprows .- rowheights

    # now we can solve the content thats inside the grid because we know where each
    # column and row is placed, how wide it is, etc.
    # note that what we did at the top was determine the protrusions of all grid content,
    # but we know the protrusions before we know how much space each plot actually has
    # because the protrusions should be static (like tick labels etc don't change size with the plot)

    gridboxes = RowCols(
        xleftcols, xrightcols,
        ytoprows, ybottomrows
    )

    for c in gl.content
        idx_rect = side_indices(c)
        bbox_cell = mapsides(idx_rect, gridboxes) do side, idx, gridside
            gridside[idx]
        end

        solving_bbox = bbox_for_solving_from_side(maxgrid, bbox_cell, idx_rect, c.side)

        c.content.layoutnodes.suggestedbbox[] = solving_bbox
    end

    nothing
end

dirlength(gl::GridLayout, c::Col) = gl.ncols
dirlength(gl::GridLayout, r::Row) = gl.nrows

function dirgaps(gl::GridLayout, dir::GridDir)
    starts = zeros(Float32, dirlength(gl, dir))
    stops = zeros(Float32, dirlength(gl, dir))
    for c in gl.content
        span = getspan(c, dir)
        start = span.start
        stop = span.stop
        starts[start] = max(starts[start], protrusion(c, startside(dir)))
        stops[stop] = max(stops[stop], protrusion(c, stopside(dir)))
    end
    starts, stops
end

dirsizes(gl::GridLayout, c::Col) = gl.colsizes
dirsizes(gl::GridLayout, r::Row) = gl.rowsizes

"""
Determine the size of a grid layout along one of its dimensions.
`Row` measures from bottom to top and `Col` from left to right.
The size is dependent on the alignmode of the grid, `Outside` includes
protrusions and paddings.
"""
function determinedirsize(gl::GridLayout, gdir::GridDir)
    sum_dirsizes = 0

    sizes = dirsizes(gl, gdir)

    for idir in 1:dirlength(gl, gdir)
        # width can only be determined for fixed and auto
        sz = sizes[idir]
        dsize = determinedirsize(idir, gl, gdir)

        if isnothing(dsize)
            # early exit if a colsize can not be determined
            return nothing
        end
        sum_dirsizes += dsize
    end

    dirgapsstart, dirgapsstop = dirgaps(gl, Col())

    forceequalprotrusiongaps = gl.equalprotrusiongaps[gdir isa Row ? 1 : 2]

    dirgapsizes = if forceequalprotrusiongaps
        innergaps = dirgapsstart[2:end] .+ dirgapsstop[1:end-1]
        m = maximum(innergaps)
        innergaps .= m
    else
        innergaps = dirgapsstart[2:end] .+ dirgapsstop[1:end-1]
    end

    inner_gapsizes = dirlength(gl, gdir) > 1 ? sum(dirgapsizes) : 0

    addeddirgapsizes = gdir isa Row ? gl.addedrowgaps : gl.addedcolgaps

    addeddirgaps = dirlength(gl, gdir) == 1 ? 0 : sum(addeddirgapsizes) do c
        if c isa Fixed
            c.x
        elseif c isa Relative
            error("Auto grid size not implemented with relative gaps")
        end
    end

    inner_size_combined = sum_dirsizes + inner_gapsizes + addeddirgaps
    return if gl.alignmode isa Inside
        inner_size_combined
    elseif gl.alignmode isa Outside
        paddings = if gdir isa Row
            gl.alignmode.padding.top + gl.alignmode.padding.bottom
        else
            gl.alignmode.padding.left + gl.alignmode.padding.right
        end
        inner_size_combined + dirgapsstart[1] + dirgapsstop[end] + paddings
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode))")
    end
end

"""
Determine the size of one row or column of a grid layout.
"""
function determinedirsize(idir, gl, dir::GridDir)

    sz = dirsizes(gl, dir)[idir]

    if sz isa Fixed
        # fixed dir size can simply be returned
        return sz.x
    elseif sz isa Relative
        # relative dir size can't be inferred
        return nothing
    elseif sz isa Auto
        # auto dir size can either be determined or not, depending on the
        # trydetermine flag
        !sz.trydetermine && return nothing

        dirsize = nothing
        for c in gl.content
            # content has to be single span to be determinable in size
            singlespanned = getspan(c, dir).start == getspan(c, dir).stop == idir

            # content has to be placed with Inner side, otherwise it's protrusion
            # content
            is_inner = c.side isa Inner

            if singlespanned && is_inner
                s = determinedirsize(c.content, dir, c.side)
                if !isnothing(s)
                    dirsize = isnothing(dirsize) ? s : max(dirsize, s)
                end
            end
        end
        return dirsize
    end
    nothing
end

function compute_col_row_sizes(spaceforcolumns, spaceforrows, gl)
    # the space for columns and for rows is divided depending on the sizes
    # stored in the grid layout

    # rows / cols with Auto size must have single span content that
    # can report its own size, otherwise they can't determine their own size
    # alternatively, if there is only one Auto row / col it can get the rest
    # that remains after the other specified sizes are subtracted

    # rows / cols with Fixed size get that size

    # rows / cols with Relative size get that fraction of the available space

    colwidths = zeros(gl.ncols)
    rowheights = zeros(gl.nrows)

    # store sizes in arrays with indices so we can keep track of them
    icsizes = collect(enumerate(gl.colsizes))
    irsizes = collect(enumerate(gl.rowsizes))

    filtersize(T::Type) = isize -> isize[2] isa T

    fcsizes = filter(filtersize(Fixed), icsizes)
    frsizes = filter(filtersize(Fixed), irsizes)

    relcsizes = filter(filtersize(Relative), icsizes)
    relrsizes = filter(filtersize(Relative), irsizes)

    acsizes = filter(filtersize(Auto), icsizes)
    arsizes = filter(filtersize(Auto), irsizes)

    aspcsizes = filter(filtersize(Aspect), icsizes)
    asprsizes = filter(filtersize(Aspect), irsizes)

    determined_acsizes = map(acsizes) do (i, c)
        (i, determinedirsize(i, gl, Col()))
    end
    det_acsizes = filter(tup -> !isnothing(tup[2]), determined_acsizes)
    nondets_c = filter(tup -> isnothing(tup[2]), determined_acsizes)
    nondet_acsizes = map(nondets_c) do (i, noth)
        i, gl.colsizes[i]
    end

    determined_arsizes = map(arsizes) do (i, c)
        (i, determinedirsize(i, gl, Row()))
    end
    det_arsizes = filter(tup -> !isnothing(tup[2]), determined_arsizes)
    nondets_r = filter(tup -> isnothing(tup[2]), determined_arsizes)
    nondet_arsizes = map(nondets_r) do (i, noth)
        i, gl.rowsizes[i]
    end

    # assign fixed sizes first
    map(fcsizes) do (i, f)
        colwidths[i] = f.x
    end
    map(frsizes) do (i, f)
        rowheights[i] = f.x
    end

    # next relative sizes
    map(relcsizes) do (i, r)
        colwidths[i] = r.x * spaceforcolumns
    end
    map(relrsizes) do (i, r)
        rowheights[i] = r.x * spaceforrows
    end

    # next determinable auto sizes
    map(det_acsizes) do (i, x)
        colwidths[i] = x
    end
    map(det_arsizes) do (i, x)
        rowheights[i] = x
    end

    # next aspect sizes
    map(aspcsizes) do (i, asp)
        index = asp.index
        ratio = asp.ratio
        rowsize = gl.rowsizes[index]
        rowheight = if rowsize isa Union{Fixed, Relative, Auto}
            if rowsize isa Auto
                if !isempty(nondet_arsizes) && any(x -> x[1] == index, nondet_arsizes)
                    error("Can't use aspect ratio with an undeterminable Auto size")
                end
            end
            rowheights[index]
        else
            error("Aspect size can only work in conjunction with Fixed, Relative, or determinable Auto, not $(typeof(gl.rowsizes[index]))")
        end
        colwidths[i] = rowheight * ratio
    end

    map(asprsizes) do (i, asp)
        index = asp.index
        ratio = asp.ratio
        colsize = gl.colsizes[index]
        colwidth = if colsize isa Union{Fixed, Relative, Auto}
            if colsize isa Auto
                if !isempty(nondet_acsizes) && any(x -> x[1] == index, nondet_acsizes)
                    error("Can't use aspect ratio with an undeterminable Auto size")
                end
            end
            colwidths[index]
        else
            error("Aspect size can only work in conjunction with Fixed, Relative, or determinable Auto, not $(typeof(gl.rowsizes[index]))")
        end
        rowheights[i] = colwidth * ratio
    end

    # next undeterminable auto sizes
    if !isempty(nondet_acsizes)
        remaining = spaceforcolumns - sum(colwidths)
        nondet_acsizes
        sumratios = sum(nondet_acsizes) do (i, c)
            c.ratio
        end
        map(nondet_acsizes) do (i, c)
            colwidths[i] = remaining * c.ratio / sumratios
        end
    end
    if !isempty(nondet_arsizes)
        remaining = spaceforrows - sum(rowheights)
        sumratios = sum(nondet_arsizes) do (i, r)
            r.ratio
        end
        map(nondet_arsizes) do (i, r)
            rowheights[i] = remaining * r.ratio / sumratios
        end
    end

    colwidths, rowheights
end

function Base.setindex!(g::GridLayout, content, rows::Indexables, cols::Indexables, side::Side = Inner())
    add_content!(g, content, rows, cols, side)
    content
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray{T, 2}) where T
    rowrange = 1:size(content_array, 1)
    colrange = 1:size(content_array, 2)
    g[rowrange, colrange] = content_array
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray{T, 1}) where T
    error("""
        You can only assign a one-dimensional content AbstractArray if you also specify the direction in the layout.
        Valid options are :h for horizontal and :v for vertical.
        Example:
            layout[:h] = contentvector
    """)
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray{T, 1}, h_or_v::Symbol) where T
    if h_or_v == :h
        g[1, 1:length(content_array)] = content_array
    elseif h_or_v == :v
        g[1:length(content_array), 1] = content_array
    else
        error("""
            Invalid direction specifier $h_or_v.
            Valid options are :h for horizontal and :v for vertical.
        """)
    end
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray, rows::Indexables, cols::Indexables)

    rows, cols = to_ranges(g, rows, cols)

    if rows.start < 1
        error("Can't prepend rows using array syntax so far, start row $(rows.start) is smaller than 1.")
    end
    if cols.start < 1
        error("Can't prepend columns using array syntax so far, start column $(cols.start) is smaller than 1.")
    end

    nrows = length(rows)
    ncols = length(cols)
    ncells = nrows * ncols

    if ndims(content_array) == 2
        if size(content_array) != (nrows, ncols)
            error("Content array size is size $(size(content_array)) for $nrows rows and $ncols cols")
        end
        # put the array content into the grid layout in order
        for (i, r) in enumerate(rows), (j, c) in enumerate(cols)
            g[r, c] = content_array[i, j]
        end
    elseif ndims(content_array) == 1
        if length(content_array) != nrows * ncols
            error("Content array size is length $(length(content_array)) for $nrows * $ncols cells")
        end
        # put the content in the layout along columns first, because that is more
        # intuitive
        for (i, (c, r)) in enumerate(Iterators.product(cols, rows))
            g[r, c] = content_array[i]
        end
    else
        error("Can't assign a content array with $(ndims(content_array)) dimensions, only 1 or 2.")
    end
    content_array
end

function GridContent(content::T, span::Span, side::Side) where T
    needs_update = Node(false)
    # connect the correct nodes
    protrusions_handle = on(protrusionnode(content)) do p
        needs_update[] = true
    end
    computedsize_handle = on(computedsizenode(content)) do c
        needs_update[] = true
    end
    GridContent{GridLayout, T}(nothing, content, span, side, needs_update,
        protrusions_handle, computedsize_handle)
end

function add_content!(g::GridLayout, content, rows, cols, side::Side)
    rows, cols = adjust_rows_cols!(g, rows, cols)

    gc = if !isnothing(content.layoutnodes.gridcontent)
        # take the existing gridcontent, remove it from its gridlayout if it has one,
        # and modify it with the new span and side
        gridc = content.layoutnodes.gridcontent
        remove_from_gridlayout!(gridc)
        gridc.span = Span(rows, cols)
        gridc.side = side
        gridc
    else
        # make a new one if none existed
        GridContent(content, Span(rows, cols), side)
    end

    content.layoutnodes.gridcontent = gc

    connect_layoutnodes!(gc)

    add_to_gridlayout!(g, gc)
end

function Base.lastindex(g::GridLayout, d)
    if d == 1
        g.nrows
    elseif d == 2
        g.ncols
    else
        error("A grid only has two dimensions, you're indexing dimension $d.")
    end
end

function is_range_within(inner::UnitRange, outer::UnitRange)
    inner.start >= outer.start && inner.stop <= outer.stop
end

function Base.getindex(g::GridLayout, rows::Indexables, cols::Indexables)
    GridPosition(g, rows, cols)
end

function Base.setindex!(gp::GridPosition, element)
    gp.layout[gp.rows, gp.cols] = element
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Capstone_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Capstone")
JLLWrappers.@generate_main_file("Capstone", UUID("0bfc08ce-5f35-58d0-8256-ec6a7315737c"))
end  # module Capstone_jll


text:
include("../src/UnitConversions.jl")
include("../src/Fluid.jl")
include("../src/Tubing.jl")
include("../src/Trajectory.jl")

using Shot
using Test

f = Fluid(50.0,0.65,650.0)
f2 = Fluid(35.0, 0.77, 650.0)

@testset "Shot" begin
    # Write your tests here.
    @test oilGravity(f) == 0.7796143250688705
    @test oilGravity(62.4) == 1.0
    @test fahrenhetToRankine(60) == 520
    @test rankineToFahrenheit(520) == 60
    @test oilDensity(f,180.0)== 39.961609859419426
    @test adiabaticOilBulkModulus(f2, 80.0,1000.0 )==231170.2088159817
    @test pressureGradient(62.4) == 0.43333333333333335
    @test adiabaticWaterBulkModulus(150.0, 1000.0) == 351400.00
    @test gasDensity(f, 1500.0, 150.0) == 4.3139351901364344
    @test psiFactor(0.1,0.01,0.3,1) == 1.1236363636363635
    @test psiFactor(0.1,0.01,0.3,2) == 1.069090909090909
    @test psiFactor(0.1,0.01,0.3,3) == 1.1690909090909092
    @test getTVDSS(10.0,[1.0, 5.0, 8.0, 15.0],[0., 1., 2., 3.])==2.2857142857142856
    @test getTemperature(10.0,[1.0, 5.0, 8.0, 15.0],[0., 1., 2., 3.])==2.2857142857142856
end

max_stars_repo_name: leytzher/Shot

max_issues_repo_name: leytzher/Shot

max_forks_repo_name: leytzher/Shot


text:
abstract type Dynamics end
abstract type Population end


text:
# PostgreSQL

@info "LibPQ.jl is available - including functionality to read / write PostgreSQL databases"

using .LibPQ
using CSV

"""
    read_table(conn:: LibPQ.Connection, tablename:: AbstractString)

Reads the content of a whole table and returns a Tables.jl compatible object.
"""
function read_table(conn:: LibPQ.Connection, tablename:: AbstractString; kwargs...)
    _checktablename(tablename)
   return execute(conn, "select * from $tablename") # SQL parameters cannot be used for table names
end

read_sql(conn:: LibPQ.Connection, sql:: AbstractString) = execute(conn, sql)

"""
    write_table(conn:: LibPQ.Connection, tablename:: AbstractString, table)

Writes data into an existing PostgreSQL table.
The table columns must have the same names as in the input table and the types must be compliant. It is OK to have more types in the PostgreSQL table
than in the input table if these columns are nullable.

Note that this method does not create a non-existing table (in contrast to the corresponding SQLite method). This is a design decision because PostgreSQL databases are usually more persistant
than (often "throw-away") SQLite databases.

This method is using `COPY FROM STDIN` on CSV data, which is much faster than uploading using SQL statements.
"""
function write_table!(conn:: LibPQ.Connection, tablename:: AbstractString, table; kwargs...)
    # Uploading data to P
    _checktable(table)
    _checktablename(tablename)
    iter = CSV.RowWriter(table)
    column_names = first(iter)
    copyin = LibPQ.CopyIn("COPY $tablename ($column_names) FROM STDIN (FORMAT CSV, HEADER);", iter)
    execute(conn, copyin)
    nothing
end


text:
# Concise JSON representations

jsonfmt(x::TextNode) = Dict(:txt => x.text)

function jsonfmt{ns, tag}(x::Elem{ns, tag})
    dict = Dict{Any, Any}('t' => tag)
    if ns !== :xhtml
        dict['n'] = ns
    end
    if hasproperties(x)
        dict['p'] = properties(x)
    end
    if haschildren(x)
        dict['c'] = [jsonfmt(c) for c in x.children]
    end
    dict
end
jsonfmt(::Void) = nothing

# values from vtree/vpatch.js
const VPATCH_NONE = 0
const VPATCH_VTEXT = 1
const VPATCH_VNODE = 2
# const VPATCH_WIDGET = 3
const VPATCH_PROPS = 4
const VPATCH_ORDER = 5
const VPATCH_INSERT = 6
const VPATCH_REMOVE = 7
# const VPATCH_THUNK = 8

jsonfmt{T <: Elem}(p::Overwrite{T}) = Dict(VPATCH_VNODE =>jsonfmt(p.b))
jsonfmt(p::Overwrite{TextNode}) = Dict(VPATCH_VTEXT => p.b.text)
jsonfmt(p::DictDiff)        = Dict(VPATCH_PROPS => p.updates)
jsonfmt(p::Reorder)         = Dict(VPATCH_ORDER => p.moves)
jsonfmt(p::Insert)          = Dict(VPATCH_INSERT => jsonfmt(p.b))
jsonfmt(p::Delete)          = Dict(VPATCH_REMOVE => nothing)

jsonfmt(ps::AbstractArray{Patch}) =
    length(ps) == 1 ?
        jsonfmt(ps[1]) :
        [jsonfmt(p) for p in ps]

function jsonfmt(ps::Dict)
    d=Dict()
    for (k,p) in ps
        d[k]=jsonfmt(p)
    end
    d
end


text:
using Test

using FlightMechanicsSimulator

@testset "ISA1976" begin
    # Check constructor
    atmosphere = ISA1976(0.0)

    # Check return type is the expected one.
    @test isa(atmosphere, ISA1976)
    # Check inheritance.
    @test isa(atmosphere, FlightMechanicsSimulator.Atmosphere)

    # Check getters
    @test isapprox(get_temperature(atmosphere), 288.15)
    @test isapprox(get_pressure(atmosphere), 101325)
    @test isapprox(get_density(atmosphere), 1.225)
    @test isapprox(get_sound_velocity(atmosphere), 340.293988026089)
    @test isapprox(get_height(atmosphere), 0.0)

    # Check against other constructor
    atmosphere2 = ISA1976(
        atmosphere.temperature,
        atmosphere.pressure,
        atmosphere.density,
        atmosphere.sound_velocity,
        atmosphere.height,
    )
    @test atmosphere2 === atmosphere
end

@testset "F16StevensAtmosphere" begin
    # Check constructor
    atmosphere = F16StevensAtmosphere(0.0)

    # Check return type is the expected one.
    @test isa(atmosphere, F16StevensAtmosphere)
    # Check inheritance.
    @test isa(atmosphere, FlightMechanicsSimulator.Atmosphere)

    # Check getters
    @test isapprox(get_temperature(atmosphere), 288, atol=1)
    @test isapprox(get_pressure(atmosphere), 101325, atol=100)
    @test isapprox(get_density(atmosphere), 1.225, atol=0.001)
    @test isapprox(get_sound_velocity(atmosphere), 340, atol=0.5)
    @test isapprox(get_height(atmosphere), 0.0)

    # Check against other constructor
    atmosphere2 = F16StevensAtmosphere(
        atmosphere.temperature,
        atmosphere.pressure,
        atmosphere.density,
        atmosphere.sound_velocity,
        atmosphere.height,
    )
    @test atmosphere2 === atmosphere
end


text:
using Test
using AxisArrays
using SampledSignals
using Statistics
using Unitful
using DSP
# using PlotAxes
using CorticalSpectralTemporalResponses

x = #= SampleBuf( =#sin.(2π .* 1000 .* range(0,stop=1,length=8000)) #,8000)
X = filt(audiospect,x,fs=8000)

err = 0.05
x_hat = filt(inv(audiospect,target_error=err,max_iterations=1000),X)
as_x_hat = filt(audiospect,x_hat,fs=8000)

@testset "CorticalSpectralTemporalResponses" begin

@testset "Spectrogram" begin
  @test_throws ErrorException filt(audiospect,collect(1:10),fs=4000)
  @test_logs (:warn, r"Unknown framerate") filt(audiospect,collect(1:10))
  @test_logs filt(audiospect,collect(1:10),fs=8000)

  @test mean(X[:,0.9kHz ..1.1kHz]) > mean(X[:,1.9kHz .. 2.1kHz])
  @test sum((as_x_hat .* (mean(X) / mean(as_x_hat)) .- X).^2) ./
    sum(X.^2) <= err
  X_small = filt(Audiospect(freq_step=2),x,fs=8000)
  @test nfrequencies(X_small) == 64
  @test mean(X_small[:,0.9kHz ..1.1kHz]) > mean(X_small[:,1.9kHz .. 2.1kHz])
  msg = r"It's recommended that you have a frame length of at least 8.0 samples"
  @test_logs((:warn,msg), Audiospect(Δt = 0.5ms))

  data = AxisArray(rand(2000),Axis{:time}(range(0,step=1/8000,length=2000)))
  @test filt(audiospect,data,fs=8000) isa CorticalSpectralTemporalResponses.AuditorySpectrogram
  # data = AxisArray(rand(2000),Axis{:time}(range(0,step=1/3000,length=2000)))
  # @test_throws ErrorException("Expected samplerate of 8000 Hz.") filt(audiospect,data)

  @test eltype(filt(audiospect,collect(1:10),fs=8000)) == float(Int)

  @test length(frequencies(X)) > 0
  @test minimum(frequencies(X)) > 0Hz
  @test maximum(frequencies(X)) > 3kHz
  @test nfrequencies(X) == 128
  @test nfrequencies(X[:,1:20]) == 20
  @test Δt(X) ≈ times(X)[2] - times(X)[1]
  @test Δf(X) ≈ frequencies(X)[2]/frequencies(X)[1]
  @test 3 < length(freq_ticks(X)[1]) < 7
  @test floor.(freq_ticks(X)[1]./10) == freq_ticks(X)[1]./10
end

# TODO: add tests for plotting cortical data as well
scalef = scalefilter()
S_cr = filt(scalef,X)
S_X̂ = filt(inv(scalef),S_cr)

ratef = ratefilter()
R_cr = filt(ratef,X)
R_X̂ = filt(inv(ratef),R_cr)

cort = cortical()
cr = filt(cort,X)

@testset "Single dimension cortical model" begin
  @test mean(abs,S_cr[:,:,0.9kHz ..1.1kHz]) >
    mean(abs,S_cr[:,:,1.9kHz .. 2.1kHz])
  @test mean(abs,R_cr[:,:,0.9kHz ..1.1kHz]) >
    mean(abs,R_cr[:,:,1.9kHz .. 2.1kHz])
  @test quantile(vec((S_X̂ .- X).^2 ./ mean(abs2,X)),0.75) < 5e-3
  @test quantile(vec((R_X̂ .- X).^2 ./ mean(abs2,X)),0.75) < 5e-3
  @test isempty(setdiff(scales(S_cr),default_scales))
  @test isempty(setdiff(rates(R_cr),default_rates))
  @test freq_ticks(S_cr) == freq_ticks(X)
  @test freq_ticks(R_cr) == freq_ticks(X)
end

@testset "Multi dimensional cortical model" begin
  X̂ = filt(inv(cort),cr)

  # TODO: fix tests
  @test mean(abs,cr[:,:,:,0.9kHz ..1.1kHz]) >
    mean(abs,cr[:,:,:,1.9kHz .. 2.1kHz])
  @test quantile(vec((X̂ .- X).^2 ./ mean(abs2,X)),0.75) < 5e-3
  @test isempty(setdiff(scales(cr),default_scales))
  @test isempty(setdiff(rates(cr),default_rates))
  @test freq_ticks(cr) == freq_ticks(X)
end

@testset "Position and keyword arugment support" begin
  sc = [1cycoct, 2cycoct]
  rt = [1Hz, 2Hz]

  @test scalefilter(sc).data ==
    scalefilter(scales=sc).data
  @test ratefilter(rt).data == ratefilter(rates=rt).data
  cr1 = cortical(sc, rt)
  cr2 = cortical(scales=sc, rates=rt)
  @test cr1.scales.data == cr2.scales.data
  @test cr1.rates.data == cr2.rates.data

  msg = "Cannot specify both a position and keyword value for scales."
  @test_throws ErrorException(msg) scalefilter(sc.+1cycoct,scales=sc)
  msg = "Cannot specify both a position and keyword value for rates."
  @test_throws ErrorException(msg) ratefilter(rt.+1Hz,rates=rt)
  msg = "Cannot specify both a position and keyword value for scales."
  @test_throws ErrorException(msg) cortical(sc.+1cycoct,scales=sc)
  msg = "Cannot specify both a position and keyword value for rates."
  @test_throws ErrorException(msg) cortical(default_scales,rt.+1Hz,rates=rt)
end

@testset "Handles implicit and explicit unit." begin
  @test scalefilter([1cycoct]).data == scalefilter([1]).data
  @test ratefilter([1Hz]).data == ratefilter([1]).data
end

# @testset "Data is plottable" begin
#   # TODO: fix tests (PlotAxes needs imporvements first)
#   @test size(PlotAxes.asplotable(X)[1],1) == length(X)
#   @test size(PlotAxes.asplotable(S_cr,quantize=(1000,1000,20,20))[1],1) == length(S_cr)
#   @test size(PlotAxes.asplotable(R_cr,quantize=(1000,1000,20,20))[1],1) == length(R_cr)
#   @test size(PlotAxes.asplotable(cr,quantize=(1000,1000,20,20))[1],1) == length(cr)

#   @test size(PlotAxes.asplotable(X,quantize=(5,5,5,5))[1],1) < length(X)
#   @test size(PlotAxes.asplotable(S_cr,quantize=(5,5,5,5))[1],1) < length(S_cr)
#   @test size(PlotAxes.asplotable(R_cr,quantize=(5,5,5,5))[1],1) < length(R_cr)
#   @test size(PlotAxes.asplotable(cr,quantize=(5,5,5,5))[1],1) < length(cr)

#   @test logrange(666) == log(666)
#   @test PlotAxes.fn_prefix(logrange) == "log"
#   @test logrange.(1:10) == log.(1:10)
# end

@testset "Data is printed in console" begin
  iobuf = IOBuffer()
  display(TextDisplay(iobuf), X)
  @test length(String(take!(iobuf))) > 0

  iobuf = IOBuffer()
  display(TextDisplay(iobuf), S_cr)
  @test length(String(take!(iobuf))) > 0

  iobuf = IOBuffer()
  display(TextDisplay(iobuf), R_cr)
  @test length(String(take!(iobuf))) > 0

  iobuf = IOBuffer()
  display(TextDisplay(iobuf), cr)
  @test length(String(take!(iobuf))) > 0
end

@testset "Repeated axis cortical model" begin

  ratef = ratefilter(axis=:rateslow)
  crR = filt(ratef,cr)
  @test ndims(crR) == 5

  ratef = ratefilter()
  @test_throws(ArgumentError("axis name :rate is used more than once"),
    filt(ratef,cr))
end

end