text:
using MimiFUND
using Mimi
using Test
using DataFrames
using CSVFiles

@testset "fund" begin

#------------------------------------------------------------------------------
#   1. Run the whole model
#------------------------------------------------------------------------------

@testset "fund-model" begin

#default model exported by fund module
default_nsteps = 1050
m = MimiFUND.get_model()
run(m)
@test Mimi.time_labels(m) == collect(1950:1:1950+default_nsteps)

#default model created by MimiFUND.get_model()
m1 = MimiFUND.get_model()
run(m1)
@test Mimi.time_labels(m1) == collect(1950:1:1950+default_nsteps)

#use optional args for MimiFUND.get_model()
new_nsteps = 10
@test_throws ErrorException m2 = MimiFUND.get_model(nsteps = new_nsteps) #should error because parameter lenghts won't match time dim

end #fund-model testset

#------------------------------------------------------------------------------
#   2. Run tests to make sure integration version (Mimi v0.5.0)
#   values match Mimi 0.4.0 values
#------------------------------------------------------------------------------
@testset "test-integration" begin

Mimi.reset_compdefs()

m = MimiFUND.get_model()
run(m)

missingvalue = -999.999
err_number = 1.0e-9
err_array = 0.0

for c in map(name, Mimi.compdefs(m)), v in Mimi.variable_names(m, c)

    #load data for comparison
    filename = joinpath(@__DIR__, "../contrib/validation_data_v040/$c-$v.csv")
    results = m[c, v]

    df = load(filename) |> DataFrame
    if typeof(results) <: Number
        validation_results = df[1,1]
        @test results ≈ validation_results atol = err_number #slight imprecision with these values due to rounding

    else
        validation_results = convert(Matrix, df)

        #replace missings with missingvalue so they can be compared
        results[ismissing.(results)] .= missingvalue
        validation_results[ismissing.(validation_results)] .= missingvalue

        #match dimensions
        if size(validation_results,1) == 1
            validation_results = validation_results'
        end

        @test results ≈ validation_results atol = err_array

    end
end

end #fund-integration testset

#------------------------------------------------------------------------------
# 3. Test marginal damages functions (test that each function does not error)
#------------------------------------------------------------------------------

@testset "test-marginaldamages" begin

# Test functions from file "new_marginaldamages.jl"

# Test the compute_scc function with various keyword arguments
scc1 = MimiFUND.compute_scc(year = 2020)
@test scc1 isa Float64   # test that it's not missing or a NaN
scc2 = MimiFUND.compute_scc(year = 2020, last_year=2300)
@test scc2 < scc1  # test that shorter horizon made it smaller
scc3 = MimiFUND.compute_scc(year = 2020, last_year=2300, equity_weights=true)
@test scc3 > scc2  # test that equity weights made itbigger
scc4 = MimiFUND.compute_scc(year = 2020, last_year=2300, equity_weights=true, eta=.8, prtp=0.01)
@test scc4 > scc3   # test that lower eta and prtp make scc higher
scc5 = MimiFUND.compute_scc(year = 2020, gas=:CH4)
@test scc5 > scc1   # test social cost of methane is higher

# Test with a modified model
m = MimiFUND.get_model()
update_param!(m, :climatesensitivity, 5)
scc6 = MimiFUND.compute_scc(m, year=2020, last_year=2300)
@test scc6 > scc2 # test that it's higher than the default because of a higher climate sensitivity

# Test get_marginal_model
mm = MimiFUND.get_marginal_model(year=2020, gas=:CH4)
run(mm)

# Test compute_scc_mm
result = MimiFUND.compute_scc_mm(year=2050)
@test result.scc isa Float64
@test result.mm isa Mimi.MarginalModel

# Test old exported versions of the functions
scc = MimiFUND.get_social_cost()
md = MimiFUND.getmarginaldamages()

end #marginaldamages testset

#------------------------------------------------------------------------------
# 4. Run basic test of Marginal Damages and MCS functionality
#------------------------------------------------------------------------------

@testset "test-mcs" begin

# mcs
MimiFUND.run_fund_mcs(10)        # Run 10 trials of basic FUND MCS
MimiFUND.run_fund_scc_mcs(10)    # Run 10 trials of FUND MCS SCC calculations

end #test-mcs testset

end #fund testset

nothing


text:
module CUDAGpuUtils

using CUDA, StaticArrays

export syncThreadsAnd,gridDimX,atomicallySetValueTrreeDim,atomicallyAddOneInt,clearLocArrdefineIndicies,computeBlocksFromOccupancy,reduce_warp,getKernelContants,assignWorkToCooperativeBlocks,getMaxBlocksPerMultiproc,reduce_warp_max,reduce_warp_min,reduce_warp_min,reduce_warp_or,reduce_warp_and,blockIdxZ,blockIdxY,blockIdxX,blockDimZ, blockDimY, blockDimX, threadIdxX, threadIdxY, threadIdxZ
export @unroll, @ifX, @ifY, @ifXY, @widL, @wid, @lan,getThreadsAndBlocksNumbForKernel

"""
convinience macro that will execute only if it has given thread Id X
"""
macro ifX(x, ex)
    return esc(:(if threadIdxX()==$x
        $ex
    end))

end

"""
convinience macro that will execute only if it has given thread Id Y
"""
macro ifY(y, ex)
    return esc(:(if threadIdxY()==$y
        $ex
    end))

end

"""
convinience macro that will execute only if it has given thread Id Y
"""
macro ifXY(x,y, ex)
        return esc(:(if threadIdxY()==$y && threadIdxX()==$x
            $ex
        end))

end
"""
convinience macro that will execute only if it has given wid and lane
"""
macro widL(innerWid,innerLane, ex)
    return esc(:(if wid==$innerWid && lane==$innerLane
        $ex
    end))

end

"""
convinience macro that will execute only if it has given wid
"""
macro wid(innerWid, ex)
    return esc(:(if wid==$innerWid
        $ex
    end))

end
"""
convinience macro that will execute only if it has given lane
"""
macro lan(innerWid,innerLane, ex)
    return esc(:(if lane==$innerLane
        $ex
    end))

end

"""
Type{maskNumb}  - type of the numbers hold in mask
G - 3 dimensional array holding ground truth segmentation
T - 3 dimensional array holding segmentation that we want to compare to ground truth
isVariedSlice - if true it will mark that  slices have varying thickness  - hence we will need later corrextion...

"""
function defineBlocks(::Type{maskNumb}
                    ,G::Array{maskNumb, 3}
                    ,T::Array{maskNumb, 3}
                    ,isVariedSlice::Bool
                    ) where{maskNumb}

end#defineBlocks

"""
wrapper to get x coordinates of thread
"""
function threadIdxX()::UInt32
    threadIdx().x
end
"""
wrapper to get y coordinates of thread
"""
function threadIdxY()::UInt32
    threadIdx().y
end
"""
wrapper to get x coordinates of thread
"""
function threadIdxZ()::UInt32
    threadIdx().z
end

"""
wrapper to get x dimension of thread block
"""
function blockDimX()::UInt32
    blockDim().x
end
"""
wrapper to get y dimension of thread block
"""
function blockDimY()::UInt32
    blockDim().y
end
"""
wrapper to get x dimension of thread block
"""
function blockDimZ()::UInt32
    blockDim().z
end

"""
wrapper to get x coordinates of thread block
"""
function blockIdxX()::UInt32
    blockIdx().x
end
"""
wrapper to get y coordinates of thread block
"""
function blockIdxY()::UInt32
    blockIdx().y
end
"""
wrapper to get x coordinates of thread block
"""
function blockIdxZ()::UInt32
    blockIdx().z
end

"""
wrapper to get x length of grid - how many thread blocks in x dimension
"""
function gridDimX()::UInt32
    gridDim().x
end

"""
wrapper for and synchronization in whole block
"""
function syncThreadsAnd(pred)
return sync_threads_and(pred)

end

"""
defining basic indicies for 3 dimensional case
"""
function defineIndicies()
    i= (blockIdx().x-1) * blockDim().x + threadIdxX()
    j = (blockIdx().y-1) * blockDim().y + threadIdxY()
    z = (blockIdx().z-1) * blockDim().z + threadIdx().z
    return (i,j,z)

end#defineIndicies

"""
calculates for getBlockTpFpFn optimal number of blocks and thread blocks
    also it poins out to maximum number of blocks that we can squeeze on device ..
args - tupple with arguments for kernel
int32Shemm per warp - we are assuming we get some shared memory and some number of it per warp
    """
function computeBlocksFromOccupancy(args, int32Shemm)
    wanted_threads =1000000
    function compute_threads(max_threads)
        if wanted_threads > max_threads
            true ? prevwarp(device(), max_threads) : max_threads
        else
            wanted_threads
        end
    end
    compute_shmem(threads) = Int64((threads/32)*int32Shemm*sizeof(Int32) )

       kernel = @cuda launch=false getBlockTpFpFn(args...)
       kernel_config = launch_configuration(kernel.fun; shmem=compute_shmem∘compute_threads)
       blocks =  kernel_config.blocks
       threads =  kernel_config.threads
       maxBlocks = attribute(device(), CUDA.DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT)

return blocks,threads,maxBlocks
end

"""
copied from https://github.com/JuliaGPU/CUDA.jl/blob/afe81794038dddbda49639c8c26469496543d831/perf/volumerhs.jl
"""
function loopinfo(name, expr, nodes...)
    if expr.head != :for
        error("Syntax error: pragma $name needs a for loop")
    end
    push!(expr.args[2].args, Expr(:loopinfo, nodes...))
    return expr
end

"""
copied from https://github.com/JuliaGPU/CUDA.jl/blob/afe81794038dddbda49639c8c26469496543d831/perf/volumerhs.jl
"""
macro unroll(expr)
    expr = loopinfo("@unroll", expr, (Symbol("llvm.loop.unroll.full"),))
    return esc(expr)
end

"""
Reduce a value across a warp and sum
"""
@inline function reduce_warp( vall::T, lanesNumb) where T

    offset = UInt16(1)
    while(offset <lanesNumb)
        vall+=shfl_down_sync(FULL_MASK, vall, offset)
        offset<<= 1
    end

    return vall
end

"""
Reduce a value across a warp and return max
"""
@inline function reduce_warp_max( vall, lanesNumb::UInt8)
    offset = UInt16(1)
    while(offset <lanesNumb)
        vall=max(vall, shfl_down_sync(FULL_MASK, vall, offset))
        offset<<= 1
    end
    return vall
end

"""
Reduce a value across a warp and return min
"""
@inline function reduce_warp_min( vall::T, lanesNumb::UInt8) where T
    offset = UInt16(1)
    while(offset <lanesNumb)
        #if(vall==0)  vall=T(10000) end
        vall=min(vall, shfl_down_sync(FULL_MASK, vall, offset))
        offset<<= 1
    end
    return vall
end

"""
Reduce a value across a warp and return or (so if any value is true it will return true)
"""
@inline function reduce_warp_or( vall::Bool, lanesNumb::UInt8)
    offset = UInt16(1)
    while(offset <lanesNumb)
        #if(vall==0)  vall=T(10000) end
        vall=vall | shfl_down_sync(FULL_MASK, vall, offset)
        offset<<= 1
    end
    return vall
end

"""
Reduce a value across a warp and return and (so if all values are true it will return true)
"""
@inline function reduce_warp_and( vall::Bool, lanesNumb::UInt8)
    offset = UInt16(1)
    while(offset <lanesNumb)
        #if(vall==0)  vall=T(10000) end
        vall=vall & shfl_down_sync(FULL_MASK, vall, offset)
        offset<<= 1
    end
    return vall
end

"""
generally  we want to get one block per slice as dimensions of slices in the medical images are friendly - so 256x256 ; 512x512 and 1024x1024 - it should all be possible
In order to avoid edge cases we will keep number of threads to some even multiply of 32 for ecxample 512
    arguments
        threadnum - number of threads per block
        slicesNumber - number of slices of our data
        pixelNumberPerSlice - number of pixels in a slice
    output
        blockNum - number of blocks we need  - generally will return number of slices
        loopNumb - number of  iteration of each single line it need to do so a single block will cover whole slice
        indexCorr - as one lane will get access to multiple data elements we need to take correction for it
"""
function getKernelContants(threadnum::Int,pixelNumberPerSlice::Int  )

indexCorr =  Int64(ceil(pixelNumberPerSlice/threadnum))
loopNumb= Int64(indexCorr-1)

return ( loopNumb, indexCorr )
end#getKernelContants

"""
Calculates  maximum number of
args - tuple with arguments needed to execute a kernel (it will not be executed it is only needed for occupancy API)
kernelFun - the function describing kernel

"""
function getMaxBlocksPerMultiproc(args,kernelFun )
 kernel = @cuda launch=false kernelFun(args...)
return CUDA.active_blocks(kernel.fun, 512)
end#getMaxBlocksPerMultiproc

"""
now plan is to create a matrix on the CPU that will be pushed onto the GPU
matrix will tell each block which slices it should manage in order to cover all slices
number of blocks will be calclulated in a way to fit them all in a single grid and enable full utilization of cooperative groups
matrix will be n x m  where
    n is number of blocks we can fit in the grid
    m is maximum number of slices managed per block - of course in most cases some blocks will not be needed in final round - in this case this will evaluate to 0 ...
  code adapted from
    https://stackoverflow.com/questions/63929929/processing-shared-work-queue-using-cuda-atomic-operations-and-grid-synchronizati/63930239#63930239
    https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#grid-synchronization-cg
    https://github.com/JuliaGPU/CUDA.jl/blob/afe81794038dddbda49639c8c26469496543d831/test/execution.jl

    generally the number of threads to avoid corner cases  needs to be 256 or 512 or 1024

    slicesNumb - number of slices to manage, generally it is assumed that each block controls one slice
    numberOfBlocksPerMultprocessor- how many block can be run in a single SM
    """
function assignWorkToCooperativeBlocks(slicesNumb, numberOfBlocksPerMultprocessor=1 )
   numberOfBlocks = attribute(device(), CUDA.DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT)*numberOfBlocksPerMultprocessor
   maxSlicesPerBlock = Int64(ceil(slicesNumb/numberOfBlocks))
   sliceAssignMatrix = ones(Int8,numberOfBlocks, maxSlicesPerBlock ).-2
   # first filling first 4 columns
   for i in 0:numberOfBlocks-1
    for j in 1:maxSlicesPerBlock-1
        z = i*(maxSlicesPerBlock-1) +j
        if(z<=slicesNumb)
          sliceAssignMatrix[i+1,j]=z
        end#if
    end#for
   end#for
   index =0
   #filling last column with what is left
   for i in (maximum(sliceAssignMatrix)+1):slicesNumb
    index+=1
    sliceAssignMatrix[index,maxSlicesPerBlock]= i
   end#for

    zz = CuArray(sliceAssignMatrix)
   return (maxSlicesPerBlock, zz,numberOfBlocks)
end #assignWorkToCooperativeBlocks

####### profiling
#   nv-nsight-cu-cli --mode=launch julia

"""
function to be invoked on CPU
use the occupancy API to calculate optimal number of threads and blocks
it assumes that we always want to have x dimension of thread block to be 32
the Y dimension will vary depending on kernel and GPU
get_shmem - function accepting  integer that is representing number of threads and returning number of bytes that the kernel will utilize with given number of threads
kernelFun - reference to kernel function that is to be analyzed by occupancy API
args - arguments tuple  for kernel function
reurn 2 tuple with optimal threads and blocks number (threads,blocks)
"""
function getThreadsAndBlocksNumbForKernel(get_shmemm,kernelFun,args)
    # calculate the amount of dynamic shared memory for a 2D block size
    #get_shmem(threads) = (sizeof(UInt32)*3*4)

    function get_threads(threads)
        threads_x = 32
        threads_y = cld(threads,threads_x )
        return (threads_x, threads_y)
    end

    kernel = @cuda launch=false kernelFun(args...)

    config = launch_configuration(kernel.fun, shmem=threads->get_shmemm(get_threads(threads)))
    threads = get_threads(config.threads)
    blocks = UInt32(config.blocks)
    return (threads,blocks)
end

end #CUDAGpuUtils

ext: jl


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Origins
export Point, RectUniform, RectGrid, Hexapolar

using ....OpticSim
using ...Emitters
using ...Geometry
using LinearAlgebra
using Distributions
using Random

abstract type AbstractOriginDistribution{T<:Real} end

Base.iterate(a::AbstractOriginDistribution, state = 1) = state > length(a) ? nothing : (generate(a, state - 1), state + 1)
Base.getindex(a::AbstractOriginDistribution, index) = generate(a, index)
Base.firstindex(a::AbstractOriginDistribution) = 0
Base.lastindex(a::AbstractOriginDistribution) = length(a) - 1
Base.copy(a::AbstractOriginDistribution) = a # most don't have any heap allocated stuff so don't really need copying

"""
    Point{T} <: AbstractOriginDistribution{T}

Encapsulates a single point origin.

```julia
Point(position::Vec3{T}) where {T<:Real}
Point(x::T, y::T, z::T) where {T<:Real}
Point(::Type{T} = Float64) where {T<:Real}
```
"""
struct Point{T} <: AbstractOriginDistribution{T}
    origin::Vec3{T}

    function Point(position::Vec3{T}) where {T<:Real}
        return new{T}(position)
    end

    function Point(x::T, y::T, z::T) where {T<:Real}
        return new{T}(Vec3{T}(x, y, z))
    end

    function Point(::Type{T} = Float64) where {T<:Real}
        return new{T}(zero(Vec3))
    end
end

Base.length(::Point) = 1
Emitters.visual_size(::Point) = 1
Emitters.generate(o::Point, ::Int64) = o.origin

"""
    RectUniform{T} <: AbstractOriginDistribution{T}

Encapsulates a uniformly sampled rectangle with user defined number of samples.

```julia
RectUniform(width::T, height::T, samples_count::Int64) where {T<:Real}
```
"""
struct RectUniform{T} <: AbstractOriginDistribution{T}
    width::T
    height::T
    samples_count::Int64
    rng::Random.AbstractRNG

    function RectUniform(width::T, height::T, samples_count::Int64; rng=Random.GLOBAL_RNG) where {T<:Real}
        return new{T}(width, height, samples_count, rng)
    end
end

Base.length(o::RectUniform) = o.samples_count
Emitters.visual_size(o::RectUniform) = max(o.width, o.height)

# generate origin on the agrid
function Emitters.generate(o::RectUniform{T}, n::Int64) where {T<:Real}
    n = mod(n, length(o))
    u = rand(o.rng, Distributions.Uniform(-one(T), one(T)))
    v = rand(o.rng, Distributions.Uniform(-one(T), one(T)))
    return zero(Vec3{T}) + ((o.width / 2) * u * unitX3(T)) + ((o.height/2) * v * unitY3(T))
end

"""
    RectGrid{T} <: AbstractOriginDistribution{T}

Encapsulates a rectangle sampled in a grid fashion.

```julia
RectGrid(width::T, height::T, usamples::Int64, vsamples::Int64) where {T<:Real}
```
"""
struct RectGrid{T} <: AbstractOriginDistribution{T}
    width::T
    height::T
    usamples::Int64
    vsamples::Int64
    ustep::T
    vstep::T

    function RectGrid(width::T, height::T, usamples::Int64, vsamples::Int64) where {T<:Real}
        return new{T}(width, height, usamples, vsamples, width / (usamples - 1), height / (vsamples - 1))
    end
end

Base.length(o::RectGrid) = o.usamples * o.vsamples
Emitters.visual_size(o::RectGrid) = max(o.width, o.height)

# generate origin on the grid
function Emitters.generate(o::RectGrid{T}, n::Int64) where {T<:Real}
    n = mod(n, length(o))
    v = o.vsamples == 1 ? zero(T) : 2 * Int64(floor(n / o.usamples)) / (o.vsamples - 1) - 1.0
    u = o.usamples == 1 ? zero(T) : 2 * mod(n, o.usamples) / (o.usamples - 1) - 1.0
    return zeros(Vec3{T}) + ((o.width / 2) * u * unitX3(T)) + ((o.height/2) * v * unitY3(T))
end

"""
    Hexapolar{T} <: AbstractOriginDistribution{T}

Encapsulates an ellipse (or a circle where halfsizeu=halfsizev) sampled in an hexapolar fashion (rings).

```julia
Hexapolar(nrings::Int64, halfsizeu::T, halfsizev::T) where {T<:Real}
```
"""
struct Hexapolar{T} <: AbstractOriginDistribution{T}
    halfsizeu::T
    halfsizev::T
    nrings::Int64

    function Hexapolar(nrings::Int64, halfsizeu::T, halfsizev::T) where {T<:Real}
        return new{T}(halfsizeu, halfsizev, nrings)
    end
end

Base.length(o::Hexapolar) = 1 + round(Int64, (o.nrings * (o.nrings + 1) / 2) * 6)
Emitters.visual_size(o::Hexapolar) = max(o.halfsizeu*2, o.halfsizev*2)

function Emitters.generate(o::Hexapolar{T}, n::Int64) where {T<:Real}
    n = mod(n, length(o))
    if n == 0
        return zeros(Vec3{T})
    else
        t = 1
        ringi = 1
        for i in 1:(o.nrings)
            t += 6 * i
            if n < t
                ringi = i
                break
            end
        end
        ρ = ringi / o.nrings
        pind = n - (t - 6 * ringi)
        ϕ = (pind / (6 * ringi)) * 2π
        u = cos(ϕ) * o.halfsizeu
        v = sin(ϕ) * o.halfsizev
        return zeros(Vec3{T}) + ρ * (u * unitX3(T) + v * unitY3(T))
    end
end

end # module Origins


text:
export SymComplex

struct SymComplex <: Number
    ex::Expression
end

SymComplex(name::Symbol) = SymComplex(Variable(name))
SymComplex(x::SymComplex) = x
SymComplex(x::Number) = SymComplex(Numeric(x))
SymComplex(x::Complex) = SymComplex(Numeric(real(x)) + Numeric(imag(x)) * Im())
SymComplex(x::Irrational{sym}) where sym = SymComplex(Constant(sym))
SymComplex(x::SymReal) = SymComplex(term(x))
SymComplex(x::SymInteger) = SymComplex(term(x))

DOMAIN_TYPES[:Complex] = SymComplex

term(x::SymComplex) = x.ex
Base.show(io::IO, t::SymComplex) = print(io, t.ex)

for (M, f, arity) in DiffRules.diffrules()
    op = :($M.$f)

    if arity == 1
        @eval $op(x::SymComplex) = track($op, x)
    elseif arity == 2
        @eval $op(x::SymComplex, y::SymComplex) = track($op, x, y)
        @eval $op(x::SymReal, y::Complex{Bool}) = track($op, promote(x, y)...)
        @eval $op(x::Complex{Bool}, y::SymReal) = track($op, promote(x, y)...)
    end
end

Base.promote_rule(::Type{<:SymComplex}, ::Type{<:Number}) = SymComplex
Base.promote_rule(::Type{<:SymReal}, ::Type{<:Complex}) = SymComplex

Base.convert(::Type{<:SymComplex}, x::SymComplex) = x
Base.convert(::Type{<:SymComplex}, x::Number) = SymComplex(x)
Base.convert(::Type{<:SymComplex}, x::Complex) = SymComplex(x)
Base.convert(::Type{<:SymComplex}, x::Irrational) = SymComplex(x)


text:
function AlongTrackAiming(arg0::OneAxisEllipsoid, arg1::Orbit, arg2::jboolean)
    return AlongTrackAiming((OneAxisEllipsoid, Orbit, jboolean), arg0, arg1, arg2)
end

function along_tile_direction(obj::AlongTrackAiming, arg0::Vector3D, arg1::GeodeticPoint)
    return jcall(obj, "alongTileDirection", Vector3D, (Vector3D, GeodeticPoint), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_singular_points(obj::AlongTrackAiming)
    return jcall(obj, "getSingularPoints", List, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using Pkg
pkg"add https://github.com/kalmarek/SymbolicWedderburn.jl#bl/sos"

import MutableArithmetics
const MA = MutableArithmetics
using MultivariatePolynomials
const MP = MultivariatePolynomials
using MultivariateBases
const MB = MultivariateBases

using GroupsCore
using SymbolicWedderburn
using PermutationGroups
using Cyclotomics
using SumOfSquares

function SymbolicWedderburn.ExtensionHomomorphism(basis::MB.MonomialBasis, action)
    monos = collect(basis.monomials)
    mono_to_index = Dict(monos[i] => i for i in eachindex(monos))
    return SymbolicWedderburn.ExtensionHomomorphism(monos, mono_to_index, action)
end

function MP.polynomialtype(::Type{<:MB.AbstractPolynomialVectorBasis{PT}}, T::Type) where PT
    C = MP.coefficienttype(PT)
    U = MA.promote_operation(*, C, T)
    V = MA.promote_operation(+, U, U)
    return MP.polynomialtype(PT, V)
end

struct SymmetricIdeal{C, GT, AT} <: Certificate.AbstractIdealCertificate
    certificate::C
    group::GT
    action::AT
end
SumOfSquares.matrix_cone_type(::Type{<:SymmetricIdeal{C}}) where {C} = SumOfSquares.matrix_cone_type(C)
Certificate.get(::Type{<:SymmetricIdeal}, ::SumOfSquares.Certificate.GramBasisType) = Vector{Vector{MB.FixedPolynomialBasis}}
Certificate.zero_basis_type(::Type{<:SymmetricIdeal}) = MB.MonomialBasis
Certificate.zero_basis(::SymmetricIdeal) = MB.MonomialBasis
function Certificate.get(certificate::SymmetricIdeal, attr::Certificate.ReducedPolynomial, poly, domain)
    return Certificate.get(certificate.certificate, attr, poly, domain)
end
_type(::Type{MOI.PositiveSemidefiniteConeTriangle}) = Float64
_type(::Type{SumOfSquares.COI.HermitianPositiveSemidefiniteConeTriangle}) = Complex{Float64}
function Certificate.get(cert::SymmetricIdeal, attr::Certificate.GramBasis, poly)
    basis = Certificate.get(cert.certificate, attr, poly)
    T = _type(SumOfSquares.matrix_cone_type(typeof(cert)))
    Rs, ms = SymbolicWedderburn.symmetry_adapted_basis(T, cert.group, basis, cert.action)
    return map(zip(Rs, ms)) do (R, m)
        F = convert(Matrix{T}, R)
        N = size(R, 1)
        d = div(N, m)
        if d > 1
            if m > 1
                ps = R * basis.monomials
                S = map(gens(cert.group)) do g
                    Si = Matrix{T}(undef, N, N)
                    for i in eachindex(ps)
                        p = ps[i]
                        q = cert.action(p, g)
                        coefs = coefficients(q, basis.monomials)
                        col = row_echelon_linsolve(R, coefs)
                        Si[:, i] = col
                    end
                    return Si
                end
                U = ordered_block_diag(S, d)
            else
                U = Matrix(1.0I, N, N)
            end
            map(1:d) do i
                FixedPolynomialBasis((transpose(U[:, i:d:(i+d*(m-1))]) * F) * basis.monomials)
            end
        else
            [FixedPolynomialBasis(F * basis.monomials)]
        end
    end
end

include(joinpath(dirname(dirname(pathof(SumOfSquares))), "examples", "block_diag.jl"))


text:
using Random, Distributions
using ProgressBars
using Plots

ENV["JULIA_DEBUG"] = "all"

# Partial HARQ-IR
# Hybrid Automatic Repeat Request (Incremental Redundancy) for the delay-tolerant
# user.

const tmax = 3 # number of frames wrong
const expectedrates = [
    0.5, # delay-critical user u
    1.8, # delay-tolerant user v
]
const ζ = 2 # path loss exponent
const a = [
    0.7, # power for user u
    0.3,  # power for user v
]
const d = [
    2, # distance of user u
    1,  # distance of user v
]
const total_frames = 2000000

channel_distribution = Normal(0, √2 / 2)
outdated_state_distribution = Normal(0, √2 / 2)
plot_points = ([], [])

Threads.@threads for SNR ∈ tqdm(0:2:40)
    # for SNR in tqdm(0:1:20)
    for τ ∈ [0.862]
        successful_frames = zeros((2, 1))

        progress_bar = 1:total_frames
        # progress_bar = ProgressBar(1:total_frames)
        # set_description(progress_bar, "($SNR dB, τ = $τ)")
        for frame ∈ progress_bar
            h = rand(channel_distribution, (2, 2)) * [1, 1im]
            # Step 1.5: Apply time-correlated small fading coefficient

            h = (
                begin
                    t = repeat([tmax], 1, tmax)
                    ωⱼ = rand(channel_distribution, (tmax, 2)) * [1, 1im]
                    h_channel = @. τ^t * h[1] + √(1 - τ^(t)) * ωⱼ
                    # g_channel = @. h_channel / √(1 + d[1] ^ ζ)
                end,
                begin
                    t = 1:tmax
                    ωⱼ =
                        rand(outdated_state_distribution, (tmax, 2)) * [1, 1im]
                    h_channel = @. τ^t * h[2] + √(1 - τ^(2t)) * ωⱼ
                    # g_channel = @. h_channel / √(1 + d[2] ^ ζ)
                end,
            )

            # for t ∈ 1:tmax
            #     if abs2(h[1][t]) > abs2(h[2][t])
            #         h[2][t], h[1][t] = h[1][t], h[2][t]
            #     end
            # end

            g = (@. h[1] / √(1 + d[1]^ζ), @. h[2] / √(1 + d[2]^ζ))

            # Step 2: Calculate AWGN
            σ = 10^(-SNR / 20)
            n = rand(Normal(0, σ), 2)

            # Step 3: Calculate reception SNR
            # Note: We are using the actual noise value n, instead of just the SNR-based σ
            SNRuu = @. a[1] * abs2(g[1]) / (a[2] * abs2(g[1]) + σ^2)
            SNRvv = @. a[2] * abs2(g[2]) / σ^2
            SNRvu = @. a[1] * abs2(g[2]) / (a[2] * abs2(g[2]) + σ^2)

            success_frame = false
            # Feedback from user to base carries ACK signal if the user decodes
            # it successfully by combining the previously received signals
            for t ∈ tmax:tmax # Only assume final transmission (3 ACKs), not intermediate ones
                success_transmission_vv =
                    (1 / t * sum(@. log2(1 + SNRvv[1:t]))) >= expectedrates[2]
                # success_transmission_vu = all(rate -> rate >= expectedrates[1], @. log2(1 + SNRvu[1:t]))
                success_transmission_vu = true # Assume that u is always successfully decoded

                if success_transmission_vv && success_transmission_vu
                    success_frame = true
                    break
                end
            end

            successful_frames +=
                [log2(1 + SNRuu[tmax]) >= expectedrates[1], success_frame]
        end

        outage_rates = @. 1 - successful_frames / total_frames

        for rate ∈ outage_rates
            if rate > 0
                push!(plot_points[1], SNR)
                push!(plot_points[2], rate)
            end
        end
    end
end

## THEORETICAL PLOTS

using SymPy

# Calculate theoretical values
λ = 1 ./ (1 .+ d .^ ζ)
r = 2 .^ expectedrates .- 1
φ = max(r[1] / (a[1] - a[2] * r[1]), r[2] / a[2])

SNR = LinRange(0, 45, 100)
σ = @. 10^(-SNR / 20)
ρ = @. 1 / σ^2

L = 10 # Gaver-Stehfest procedure parameter

FY(y, τ, σ) = begin
    @info "Called with $y $σ"
    FYₖ = Array{Float64}(undef, L)
    for k ∈ 1:L
        l = floor(Int, (k + 1) / 2)
        l = l:min(k, L ÷ 2) # typo in paper
        ŵₖ =
            (-1)^(L / 2 + k) * sum(@. (
                l^(L ÷ 2 + 1) / factorial(L ÷ 2) *
                binomial(L ÷ 2, l) *
                binomial(2l, l) *
                binomial(l, k - l)
            ))
        t = 1:tmax
        FYₖ[k] =
            ŵₖ / k *
            (
                prod(@. 1 + a[2] * λ[2] * (1 - τ^(2t)) * k * log(2) / y / σ^2) *
                (
                    1 + (sum(@. a[2] * λ[2] * τ^(2t) * k * log(2) / (
                        y * σ^2 + a[2] * λ[2] * (1 - τ^(2t)) * k * log(2)
                    )))
                )
            )^(-1)
    end
    sum(FYₖ)
end

plotlyjs()
plot()
τ = 0.8621
# for τ ∈ [0.8621]
    PoutTheoreticalNOMACC = @. FY(2^(tmax * expectedrates[2]) - 1, τ, σ)
    # PoutTheoreticalNOMAIR = @. FY(σ, τ, σ)
    PoutTheoreticalNOMAIR = @. FY((tmax * 2^(expectedrates[2]) - 1) / 1, τ, σ)
    # for σ₀ ∈ σ
    #     push!(PoutTheoreticalNOMA, )
    # end
    t = 1:tmax

    meijerg = convert(Float64, sympy.functions.special.hyper.meijerg(
        reshape([1 repeat([2], 1, tmax)], tmax + 1),
        [],
        [],
        reshape([repeat([1], 1, tmax) 0], tmax + 1),
        2^(tmax * expectedrates[2]),
    ).evalf())

    PoutUpperNOMAIR = vec(meijerg .* (1 + sum(@. τ^(2t) / (1 - τ^(2t))))^(-1) ./ prod((@. a[2] * ρ' * λ[2] * (1 - τ^(2t))), dims = 1))

    # OMA - Outdated CSI
    PoutTheoreticalOMA = @. [
        1 - exp(
            -2 * (2^(2 * expectedrates[1]) - 1) / ((2 - τ^(2tmax)) * ρ * λ[1]),
        ),
        1 - 2 * exp(-(2^(2 * expectedrates[2]) - 1) / ρ / λ[2]) + exp(
            -2 * (2^(2 * expectedrates[2]) - 1) / ((2 - τ^(2tmax)) * ρ * λ[2]),
        ),
    ]
    PoutTheoreticalNOMAInit = @. [
        1 - exp(
            -2 * r[1] * σ^2 / (2 - τ^(2tmax)) / (a[1] - a[2] * r[1]) / λ[1],
        ),
        1 - 2 * exp(-φ * σ^2 / λ[2]) +
        exp(-2 * φ * σ^2 / (2 - τ^(2tmax)) / λ[2]),
    ]

    plot!(
        SNR,
        PoutTheoreticalNOMACC,
        yscale = :log10,
        label = "NOMA-HARQ-CC, User v, τ = $τ",
        lw = 2,
    )
    plot!(
        SNR,
        PoutTheoreticalNOMAIR,
        yscale = :log10,
        label = "NOMA-HARQ-IR, User v, τ = $τ, Lower Bound",
        lw = 2,
    )
    plot!(
        SNR,
        PoutUpperNOMAIR,
        yscale = :log10,
        label = "NOMA-HARQ-IR, User v, τ = $τ, Upper Bound",
        lw = 2,
    )
    plot!(
        SNR,
        PoutTheoreticalNOMAInit,
        yscale = :log10,
        label = "Correct",
        lw = 4,
    )
    plot!(
        SNR,
        PoutTheoreticalOMA,
        yscale = :log10,
        label = ["OMA, User u, τ = $τ" "OMA, User v, τ = $τ"],
        lw = 2,
    )
# end

plot!(
    plot_points[1],
    plot_points[2],
    seriestype = :scatter,
    yscale = :log10,
    label = "NOMA Simulation results",
)

xlabel!("SNR (dB)")
ylabel!("Outage rate")

display(plot!(legend = :best, size = (1000, 800)))


text:
@recipe function f(b::Backtest, name::String)
    xguide := "lead times"
    @series begin
        seriestype := :path
        label := "MAE " * name
        marker := :circle
        b.mae
    end
    @series begin
        seriestype := :path
        label := "Mean CRPS " * name
        marker := :circle
        b.mean_crps
    end
end


text:
"""
$(TYPEDEF)

Houbolt's integration scheme with given step-size.

### Fields

- `Δt` -- step-size

### References

See [[HOU50]](@ref).
"""
struct Houbolt{N} <: AbstractSolver
    Δt::N
end

Houbolt(; Δt::N) where N = Houbolt(Δt)

function _init(alg::Houbolt, M, C, K)
    Δt = alg.Δt

    # compute integration constants
    a₀ = 2/Δt^2
    a₁ = 11/(6*Δt)
    a₂ = 5/Δt^2
    a₃ = 3/Δt
    a₄ = -2*a₀
    a₅ = -a₃/2
    a₆ = a₀/2
    a₇ = a₃/9

    K̂ = K + a₀ * M + a₁ * C

    return a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇, K̂
end

function _solve(alg::Houbolt{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int; kwargs...) where {N, VT, ST, XT<:Tuple{VT, VT}}

    sys = system(ivp)
    (U₀, U₀′) = initial_state(ivp)

    IMAX = NSTEPS + 1
    M, C, K, R = _unwrap(sys, IMAX)

    U₀′′ = M \ (R[1] - C * U₀′ - K * U₀)
    a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇, K̂ = _init(alg, M, C, K)
    K̂⁻¹ = factorize(K̂)

    U = Vector{VT}(undef, IMAX)
    U[1] = U₀

    # special starting procedure to calculate U(Δt) and U(2Δt)
    V = _solve(CentralDifference(Δt=alg.Δt), ivp, 2) |> displacements
    U[2] = V[2]
    U[3] = V[3]

    @inbounds for i in 3:NSTEPS
        mᵢ = M * (a₂ * U[i] + a₄ * U[i-1] + a₆ * U[i-2])
        cᵢ = C * (a₃ * U[i] + a₅ * U[i-1] + a₇ * U[i-2])
        R̂ᵢ₊₁ = R[i+1] + mᵢ + cᵢ
        U[i+1] = K̂⁻¹ \ R̂ᵢ₊₁
    end

    return _build_solution(alg, U, NSTEPS)
end

function _build_solution(alg::Houbolt{N}, U, NSTEPS) where {N}
    t = range(zero(N), step=alg.Δt, length=(NSTEPS+1))
    return Solution(alg, U, nothing, nothing, t)
end


text:
using BayesianNonparametrics

(X, Y) = bloobs(randomize = false)

μ0 = vec(mean(X, 1))
κ0 = 5.0
ν0 = 9.0
Σ0 = cov(X)

H = WishartGaussian(μ0, κ0, ν0, Σ0)

model = DPM(H)
initialisation = KMeansInitialisation(k = 10)

modelBuffer = init(X, model, initialisation)
models = train(modelBuffer, DPMHyperparam(), Gibbs(maxiter = 500))

densities = Float64[m.energy for m in models]
activeComponents = Int[sum(m.weights .> 0) for m in models]
assignments = [m.assignments for m in models]

A = reduce(hcat, assignments)
(N, D) = size(X)

PSM = zeros(N, N)
M = size(A, 2)
for i in 1:N
  for j in 1:i
    PSM[i, j] = sum(A[i,:] .== A[j,:]) / M
    PSM[j, i] = sum(A[i,:] .== A[j,:]) / M
  end
end

mink = minimum([length(m.weights) for m in models])
maxk = maximum([length(m.weights) for m in models])

(peassignments, _) = pointestimate(PSM, method = :average, mink = mink, maxk = maxk)


text:
module ConvexFit

using LinearAlgebra: mul!, dot
using Printf
using Random: AbstractRNG, randperm!

import Base: ==, show, push!, getindex

export convexfit,

       ModelSelection,
       GridSearch,
       grid,
       OptimSearch,
       optim,
       ModelSelectionResult,
       GridSearchResult,
       kfoldcv,
       loocv

if VERSION < v"1.3"
    import LinearAlgebra: mul!
    using LinearAlgebra: lmul!
    function mul!(C::Vector, A::AbstractMatrix, B::Vector, α::Bool, β::Real)
        lmul!(β, C)
        C .= C + A*B
    end
    using Random: GLOBAL_RNG
    default_rng() = GLOBAL_RNG
else
    using Random: default_rng
end

include("utils.jl")
include("result.jl")
include("fit.jl")
include("select.jl")

end # module


text:
# Ref FluxML/Flux.jl#1209 1D input
@testset "BPTT-1D" begin
  seq = [rand(Float32, 2) for i = 1:3]
  for r ∈ [RNN,]
    rnn = r(2, 3)
    Flux.reset!(rnn)
    grads_seq = gradient(Flux.params(rnn)) do
        sum(rnn.(seq)[3])
    end
    Flux.reset!(rnn);
    bptt = gradient(Wh -> sum(tanh.(rnn.cell.Wi * seq[3] + Wh *
                                  tanh.(rnn.cell.Wi * seq[2] + Wh *
                                        tanh.(rnn.cell.Wi * seq[1] +
                                            Wh * rnn.cell.state0
                                        + rnn.cell.b)
                                  + rnn.cell.b)
                            + rnn.cell.b)),
                    rnn.cell.Wh)
    @test grads_seq[rnn.cell.Wh] ≈ bptt[1]
end
end

# Ref FluxML/Flux.jl#1209 2D input
@testset "BPTT-2D" begin
  seq = [rand(Float32, (2, 1)) for i = 1:3]
  for r ∈ [RNN,]
    rnn = r(2, 3)
    Flux.reset!(rnn)
    grads_seq = gradient(Flux.params(rnn)) do
        sum(rnn.(seq)[3])
    end
    Flux.reset!(rnn);
    bptt = gradient(Wh -> sum(tanh.(rnn.cell.Wi * seq[3] + Wh *
                                  tanh.(rnn.cell.Wi * seq[2] + Wh *
                                        tanh.(rnn.cell.Wi * seq[1] +
                                            Wh * rnn.cell.state0
                                        + rnn.cell.b)
                                  + rnn.cell.b)
                            + rnn.cell.b)),
                    rnn.cell.Wh)
    @test grads_seq[rnn.cell.Wh] ≈ bptt[1]
end
end

@testset "RNN-shapes" begin
    @testset for R in [RNN, GRU, LSTM]
        m1 = R(3, 5)
        m2 = R(3, 5)
        x1 = rand(Float32, 3)
        x2 = rand(Float32,3,1)
        Flux.reset!(m1)
        Flux.reset!(m2)
        @test size(m1(x1)) == (5,)
        @test size(m1(x1)) == (5,) # repeat in case of effect from change in state shape
        @test size(m2(x2)) == (5,1)
        @test size(m2(x2)) == (5,1)
    end
end

@testset "RNN-input-state-eltypes" begin
  @testset for R in [RNN, GRU, LSTM]
      m = R(3, 5)
      x = rand(Float64, 3, 1)
      Flux.reset!(m)
      @test_throws MethodError m(x)
  end
end


text:
"""
```
function impulse_responses(m::AbstractDSGEModel, paras::Union{Vector{S}, Matrix{S}},
                           input_type::Symbol, method::Symbol,
                           lags::Int, observables::Vector{Symbol},
                           shocks::Vector{Symbol},
                           n_obs_shock::Int; parallel::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           flip_shocks::Bool = false,
                           use_intercept::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false,
                           minimize::Bool = true,
                           forecast_string::String = "",
                           verbose::Symbol = :high) where {S<:Real}
```
computes the impulse responses of a VAR(p) approximation to a DSGE.

### Inputs
* `m::Union{AbstractDSGEModel,AbstractDSGEVARModel}`: DSGE/DSGEVAR model object
* `paras::Matrix{S}` or `paras::Vector{S}`: parameters to calibrate the model
* `input_type::Symbol`: `:mode` specifies a modal impulse response, and
    `:full` specifies a full-distribution forecast if `paras` is not given.
    This argument is also used to construct the file names of computed `MeansBands`.
* `method::Symbol`: type of impulse response to compute. The options are
    `:cholesky`, `:maximum_business_cycle_variance` or `:maxBC`,
    and `:cholesky_long_run` or `:choleskyLR`. See `?cholesky_shock`,
    `?maxBC_shock`, and `?choleskyLR_shock`.
* `lags::Int`: number of lags in the VAR(p) approximation, i.e. p = lags
* `observables::Vector{Symbol}`: observables to be used in the VAR. These can be
    any of the observables or pseudo-observables in `m`.
* `shocks::Vector{Symbol}`: (structural) exogenous shocks to be used in the DSGE-VAR.
    These shocks must be in `m`.
* `n_obs_shock::Int`: the index of the observable corresponding to the orthogonalized shock causing the impulse response.

### Keywords
* `parallel::Bool`: use parallel workers or not
* `frequency_band::Tuple{S,S}`: See `?maxBC_shock`.
* `flip_shocks::Bool`: impulse response shocks are negative by default. Set to `true` for
    a positive signed shock.
* `density_bands::Vector{Float64}`: bands for full-distribution IRF computations
* `create_meansbands::Bool`: set to `true` to save output as a `MeansBands` object.
* `minimize::Bool`: choose shortest interval if true, otherwise just chop off lowest and
    highest (percent/2)
* `forecast_string::String`: string tag for identifying this impulse response
* `verbose::Symbol`: quantity of output desired

"""
function impulse_responses(m::AbstractDSGEModel, paras::Union{Vector{S}, Matrix{S}},
                           input_type::Symbol, method::Symbol,
                           lags::Int, observables::Vector{Symbol},
                           shocks::Vector{Symbol},
                           n_obs_shock::Int; parallel::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           flip_shocks::Bool = false,
                           use_intercept::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false,
                           minimize::Bool = true,
                           forecast_string::String = "",
                           verbose::Symbol = :high) where {S<:Real}

    # Compute VAR coefficients implied by DSGE
    if isa(m, AbstractDSGEModel)
        _dsgevar = DSGEVAR(m, shocks, "ss0")
        DSGE.update!(_dsgevar, lags = lags, observables = observables)
    else
        _dsgevar = m

    end

    return impulse_responses(_dsgevar, paras, input_type, method, n_obs_shock;
                             parallel = parallel, frequency_band = frequency_band,
                             flip_shocks = flip_shocks, use_intercept = use_intercept,
                             density_bands = density_bands, save_as_DSGE = isa(m, AbstractDSGEModel),
                             create_meansbands = create_meansbands, minimize = minimize,
                             forecast_string = forecast_string, verbose = verbose)
end

function impulse_responses(m::AbstractDSGEVARModel, paras::Union{Vector{S}, Matrix{S}},
                           input_type::Symbol, method::Symbol,
                           lags::Int, observables::Vector{Symbol},
                           shocks::Vector{Symbol},
                           n_obs_shock::Int; parallel::Bool = false,
                           use_intercept::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           flip_shocks::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false,
                           minimize::Bool = true,
                           forecast_string::String = "",
                           save_as_DSGE::Bool = false,
                           verbose::Symbol = :high) where {S<:Real}
    DSGE.update!(m, lags = lags, observables = observables, shocks = shocks)
    return impulse_responses(m, paras, input_type, method, n_obs_shock;
                             parallel = parallel,
                             use_intercept = use_intercept,
                             frequency_band = frequency_band,
                             flip_shocks = flip_shocks,
                             density_bands = density_bands,
                             create_meansbands = create_meansbands,
                             minimize = minimize,
                             forecast_string = forecast_string,
                             save_as_DSGE = save_as_DSGE, verbose = verbose)
end

function impulse_responses(m::AbstractDSGEVARModel, paras::Union{Vector{S}, Matrix{S}},
                           input_type::Symbol, method::Symbol,
                           n_obs_shock::Int; parallel::Bool = false,
                           use_intercept::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           flip_shocks::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false, test_meansbands::Bool = false,
                           minimize::Bool = true, save_as_DSGE::Bool = false,
                           forecast_string::String = "",
                           verbose::Symbol = :high) where {S<:Real}
    if n_obs_shock <= 0
        error("To use method $method, user must specify the index of" *
              " the target observable with keyword `n_obs_shock`.")
    end
    if isa(paras, Vector{S})
        paras = reshape(paras, 1, length(paras))
    end

    # Set up computation method
    mapfcn = parallel ? pmap : map
    h = impulse_response_horizons(m)

    # Compute IRFs
    paras = mapslices(x -> [vec(x)], paras, dims = 2)
    function get_irf!(para)
        DSGE.update!(m, para)
        return impulse_responses(m, method, n_obs_shock;
                                 horizon = h, use_intercept = use_intercept,
                                 flip_shocks = flip_shocks,
                                 frequency_band = frequency_band)
    end

    irf_output = mapfcn(para -> get_irf!(para), paras)

    if create_meansbands
        # Set up metadata and output from IRFs computation
        metadata = Dict{Symbol,Any}()
        metadata[:para] = input_type
        metadata[:cond_type] = :none
        metadata[:product] = :dsgevarirf
        metadata[:class] = :obs # We default everything to an observable
        metadata[:date_inds] = OrderedDict()

        # Set up for loop over variable names
        means = DataFrame()
        bands = Dict{Symbol,DataFrame}()
        observables = collect(keys(get_observables(m)))
        metadata[:indices] =
            OrderedDict{Symbol,Int}(name => name_i
                                    for (name_i, name) in enumerate(observables))

        # Means and Bands for each variable in a class
        for (name_i,name) in enumerate(observables)
            # irf_output is Vector{nperiod x nobs} -> for each observable,
            # we want to select its specific IRF, i.e. map(x -> x[:,obs_index]).
            # This creates a nperiod x ndraws matrix, which we want to transpose
            # to get a ndraws x nperiod matrix
            single_var = Matrix(reduce(hcat, map(x -> x[:,name_i], irf_output))')
            means[!,name] = vec(mean(single_var, dims = 1))
            bands[name]   = find_density_bands(single_var, density_bands;
                                               minimize = minimize)
        end
        mb = MeansBands(metadata, means, bands)

        # Save MeansBands
        if !test_meansbands
            tail = if method == :cholesky
                :cholesky
            elseif method == :maxBC || method == :maximum_business_cycle_variance
                :maxBC
            else
                :choleskyLR
            end

            var_names = save_as_DSGE ?
                Symbol("_" * join(string.(DSGE.detexify(observables)), "_") * "_") : Symbol("_")
            fp = get_meansbands_output_file(save_as_DSGE ? get_dsge(m) : m, input_type, :none,
                                            Symbol(:dsgevarirf, :obs,
                                                   var_names, tail),
                                            forecast_string = forecast_string)
            dirpath = dirname(fp)
            isdir(dirpath) || mkpath(dirpath)
            JLD2.jldopen(fp, true, true, true, IOStream) do file
                write(file, "mb", mb)
            end
            println(verbose, :high, "  " * "wrote " * basename(fp))
        end
        return mb
    else
        # Reshape irf_output to nobs x nperiod x ndraw
        # return cat(map(x -> x', irf_output)..., dims = 3)
        return cat(irf_output..., dims = 3)
    end
end


text:
using Test, Random, QPDAS
using LinearAlgebra, SparseArrays

@testset "CholeskySpecial" begin
    include("testCholeskySpecial.jl")
end

@testset "testCholeskySpecialShifted" begin
    include("testCholeskySpecialShifted.jl")
end

@testset "testCholeskySpecialShifted" begin
    include("testCholeskySpecialShiftedProjection.jl")
end

@testset "testBoxQP" begin
    include("testBoxQP.jl")
end

@testset "testQPNoScaling" begin
    include("testQPNoScaling.jl")
end

@testset "testQP" begin
    include("testQP.jl")
end

@testset "testQPsemidefinite" begin
    include("testQPsemidefinite.jl")
end

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
using Base.Test

# test_meshwarp()
img = reshape(float(collect(1:121).%2), 11, 11) # 11x11 checkerboard
triangles = [1 2 3]

tform = [1 0 0;
        0 1 0;
        0 0 1]
src = [0.0 0.0;
        0.0 10.0;
        10.0 0.0]
dst = warp_pts(tform, src)
offset = [0,0]
img_warped, warped_offset = meshwarp(img, src, dst, triangles, offset)
# @test_approx_eq_eps img_warped img 1e-5
@test warped_offset == [0,0]

tform = [1 0 0;
        0 1 0;
        0 0 1]
src = [0.0 0.0;
        0.0 10.0;
        10.0 0.0]
dst = warp_pts(tform, src)
offset = [5,10]
img_warped, warped_offset = meshwarp(img, src, dst, triangles, offset)
@test_approx_eq_eps img_warped zeros(11,11) 1e-10
@test warped_offset == [0,0]

tform = [1 0 0;
        0 1 0;
        0 0 1]
src = [0.0 0.0;
        0.0 20.0;
        20.0 0.0]
dst = warp_pts(tform, src)
offset = [5,10]
img_warped, warped_offset = meshwarp(img, src, dst, triangles, offset)
@test img_warped[6,11] == 1.0
@test warped_offset == [0,0]

src = [2.0 2.0;
        10.0 2.0;
        6.0 10.0]
dst = [4.0 2.0;
        8.0 2.0;
        6.0 10.0]
offset = [0, 0]
img_warped, warped_offset = meshwarp(img, src, dst, triangles, offset)
@test warped_offset == [4,2]

# test_poly2mask()
# vertices=[1.5 2.5; 2.5 1.5; 3.5 4.5]
# img, off = poly2mask(vertices[:,1],vertices[:,2])
# """ off should be [1,1] and img should be:
# false  false  false  false  false
# false   true   true  false  false
# false  false   true   true  false
# false  false  false  false  false
# """
# vertices=[1.5 2.4; 2.5 1.5; 3.5 4.5]
# img, off = poly2mask(vertices[:,1],vertices[:,2])
# """ off should be [1,1] and img should be:
# false  false  false  false  false
# false   true  false  false  false
# false  false   true  false  false
# false  false  false  false  false,
# """
# vertices=[1 2; 2.5 1.5; 3.5 4.5]
# img, off = poly2mask(vertices[:,1],vertices[:,2])
# """ corner case: off should be [1,1] and img should be:
# false   true  false  false  false
# false   true   true  false  false
# false  false   true   true  false
# false  false  false  false  false
# """
# p=plot(x=vertices[:,1],y=vertices[:,2])
# display(p)
# img,off


text:
function prepare_coordinates(frame_num::Int64)::Matrix{Coordinate{Float32}}
    coordinates = Matrix{Coordinate{Float32}}(undef, 3, frame_num)
    for atom_idx in 1:3, frame_idx in 1:frame_num
        coordinates[atom_idx, frame_idx] =
            Coordinate([atom_idx * 10.0f0 + frame_idx * 1.0f0 + 0.1f0,
                        atom_idx * 10.0f0 + frame_idx * 1.0f0 + 0.2f0,
                        atom_idx * 10.0f0 + frame_idx * 1.0f0 + 0.3f0])
    end
    coordinates
end

@testset "read_dcd function" begin
    trj = read_dcd("data/test_position.dcd")
    @test isapprox(Array(trj.coordinates[1,1]),    [15.308, 14.180, -2.955], atol = 1e-3)
    @test isapprox(Array(trj.coordinates[5,1001]), [-2.577, 88.384, -7.513], atol = 1e-3)

    trj_frame = read_dcd("data/test_position.dcd", frame_indices = 1:10:1001)
    @test isapprox(Array(trj_frame.coordinates[1,  1]), Array(trj.coordinates[1,    1]), atol = 1e-3)
    @test isapprox(Array(trj_frame.coordinates[5,101]), Array(trj.coordinates[5, 1001]), atol = 1e-3)

    trj_step = read_dcd("data/test_position.dcd", step = 10)
    @test isapprox(Array(trj_step.coordinates[1,  1]), Array(trj.coordinates[1,    1]), atol = 1e-3)
    @test isapprox(Array(trj_step.coordinates[5,101]), Array(trj.coordinates[5, 1001]), atol = 1e-3)

end

@testset "write_dcd function" begin
    trj = read_dcd("data/test_position.dcd")
    write_dcd("data/write_test.dcd", trj)
    writed_trj = read_dcd("data/write_test.dcd")
    @test trj.nframe == writed_trj.nframe
    @test trj.natom  == writed_trj.natom
    @test isapprox(trj.coordinates[3, 10].x, writed_trj.coordinates[3, 10].x, atol = 1e-3)
end

@testset "read_pdb function" begin
    trj = read_pdb("data/1aki.pdb")
    @test isapprox(Array(trj.coordinates[1,1]), [35.365, 22.342, -11.980], atol = 1e-3)
    @test isapprox(Array(trj.coordinates[1079, 1]), [43.755, 23.843, 8.038], atol = 1e-3)
    @test trj.attributes[123].atomname == "N"
    @test trj.attributes[200].resid == 26
    @test trj.attributes[231].resname == "VAL"

    trj = read_pdb("data/1aki.pdb", model = :AA)
    @test isapprox(Array(trj.coordinates[1,1]), [35.365, 22.342, -11.980], atol = 1e-3)
    @test isapprox(Array(trj.coordinates[1079, 1]), [43.755, 23.843, 8.038], atol = 1e-3)
    @test trj.attributes[123].atomname == "N"
    @test trj.attributes[200].resid == 26
    @test trj.attributes[231].resname == "VAL"
    @test isapprox(trj.attributes[123].mass, 14.0069, atol = 1e-3)

    trj = read_pdb("data/1aki.pdb", model = :CA)
    @test isapprox(Array(trj.coordinates[1,1]), [35.365, 22.342, -11.980], atol = 1e-3)
    @test isapprox(Array(trj.coordinates[1079, 1]), [43.755, 23.843, 8.038], atol = 1e-3)
    @test trj.attributes[123].atomname == "N"
    @test trj.attributes[200].resid == 26
    @test trj.attributes[231].resname == "VAL"
    @test isapprox(trj.attributes[231].mass, 99.06841, atol = 1e-3)

    # test for multi frame pdb case
    trj = read_pdb("data/test_position.pdb")
    @test trj.nframe == 1001
    @test trj.natom  == 5
    @test isapprox(Array(trj.coordinates[1,   1]), [15.308, 14.180, -2.955], atol = 1e-3)
    @test isapprox(Array(trj.coordinates[5,1001]), [-2.578, 88.385, -7.513] ,atol = 1e-3)
end

@testset "write_pdb function" begin
    trj = read_pdb("data/1aki.pdb")
    write_pdb("data/write_test.pdb", trj)
    writed_trj = read_pdb("data/write_test.pdb")
    @test isapprox(trj.coordinates[1, 1].x, writed_trj.coordinates[1, 1].x, atol = 1e-3)
    @test isapprox(trj.coordinates[1079, 1].x, writed_trj.coordinates[1079, 1].x, atol = 1e-3)
    @test trj.attributes[123].atomname == writed_trj.attributes[123].atomname
    @test trj.attributes[200].resid    == writed_trj.attributes[200].resid
    @test trj.attributes[231].resname  == writed_trj.attributes[231].resname

    write_pdb("data/write_test.pdb", trj, tempfactor = 0.00)
    @test isapprox(trj.coordinates[1, 1].x, writed_trj.coordinates[1, 1].x, atol = 1e-3)
    @test isapprox(trj.coordinates[1079, 1].x, writed_trj.coordinates[1079, 1].x, atol = 1e-3)
    @test trj.attributes[123].atomname == writed_trj.attributes[123].atomname
    @test trj.attributes[200].resid    == writed_trj.attributes[200].resid
    @test trj.attributes[231].resname  == writed_trj.attributes[231].resname

    temp_arr = []
    open("data/write_test.pdb", "r") do fp
        for line in eachline(fp)
            if occursin(r"^ATOM", line)
                push!(temp_arr, line[63:66])
            end
        end
    end
    @test temp_arr[5]   == "0.00"
    @test temp_arr[756] == "0.00"
end

@testset "read_xyz function" begin
    trj = read_xyz("data/test_position.xyz")
    @test isapprox(Array(trj.coordinates[1, 1]),    [15.308, 14.180, -2.955], atol = 1e-3)
    @test isapprox(Array(trj.coordinates[5, 1001]), [-2.577, 88.384, -7.513], atol = 1e-3)
    @test trj.attributes[3].atomname == "CA"
end

@testset "write_xyz function" begin
    trj = read_xyz("data/test_position.xyz")
    write_xyz("data/write_test.xyz", trj)
    writed_trj = read_xyz("data/write_test.xyz")
    @test trj.nframe == writed_trj.nframe
    @test trj.natom  == writed_trj.natom
    @test writed_trj.attributes[3].atomname == trj.attributes[3].atomname
    @test isapprox(trj.coordinates[3, 10].x, writed_trj.coordinates[3, 10].x, atol = 1e-3)

    trj = read_dcd("data/test_position.dcd")
    write_xyz("data/write_test.xyz", trj)
    writed_trj = read_xyz("data/write_test.xyz")
    @test writed_trj.attributes[3].atomname == "UNK"
end

@testset "get_frame" begin
    coordinates = Matrix{Coordinate{Float32}}(undef, 2, 2)
    coordinates[1,1] = Coordinate([1.1f0, 1.2f0, 1.3f0])
    coordinates[1,2] = Coordinate([2.1f0, 2.2f0, 2.3f0])
    coordinates[2,1] = Coordinate([3.1f0, 3.2f0, 3.2f0])
    coordinates[2,2] = Coordinate([4.1f0, 4.2f0, 4.1f0])
    attributes  = [Attribute(atomname = "hoge"), Attribute(atomname = "huga")]
    trj = Trajectory(coordinates, attributes)
    frame1 = get_frame(1, trj)
    frame2 = get_frame(2, trj)
    @test isapprox(Array(frame1.coordinates[1])[1], 1.1f0, atol = 1e-7)
    @test frame2.attributes[2].atomname   == "huga"
end

@testset "get_atom" begin
    coordinates = Matrix{Coordinate{Float32}}(undef, 2, 2)
    coordinates[1,1] = Coordinate([1.1f0, 1.2f0, 1.3f0])
    coordinates[1,2] = Coordinate([2.1f0, 2.2f0, 2.3f0])
    coordinates[2,1] = Coordinate([3.1f0, 3.2f0, 3.2f0])
    coordinates[2,2] = Coordinate([4.1f0, 4.2f0, 4.1f0])
    attributes  = [Attribute(atomname = "hoge"), Attribute(atomname = "huga")]
    trj = Trajectory(coordinates, attributes)
    atom1_series = get_atom(1, trj)
    atom2_series = get_atom(2, trj)
    @test isapprox(atom1_series[1].coordinate.x, 1.1f0, atol=1e-7)
    @test atom2_series[2].attribute.atomname == "huga"

    frame = Frame(coordinates[:,1], attributes)
    atom1 = get_atom(1, frame)
    atom2 = get_atom(2, frame)
    @test isapprox(atom1.coordinate.x, 1.1f0, atol = 1e-7)
    @test atom2.attribute.atomname == "huga"
end

@testset "clip_trajectory(query::Integer)" begin
    coordinates = prepare_coordinates(3)
    attributes = [Attribute(atomname = "hoge"), Attribute(atomname = "huga"), Attribute(atomname = "piyo")]
    trj = Trajectory(coordinates, attributes)

    cliped_frame_int = clip_trajectory(2, trj)
    @test isapprox(cliped_frame_int.coordinates[2, 1].x, trj.coordinates[2, 2].x, atol = 1e-3)
    @test cliped_frame_int.attributes   == trj.attributes

    cliped_atom_int = clip_trajectory(2, trj, query_key = :atom)
    @test isapprox(cliped_atom_int.coordinates[1, 2].x,  trj.coordinates[2, 2].x, atol = 1e-3)
    @test cliped_atom_int.attributes   == [trj.attributes[2]]

end

@testset "clip_trajectory(query::Union{Array, StepRange})" begin
    coordinates = prepare_coordinates(3)
    attributes = [Attribute(atomname = "hoge"), Attribute(atomname = "huga"), Attribute(atomname = "piyo")]
    trj = Trajectory(coordinates, attributes)

    cliped_frame_array = clip_trajectory([2, 3], trj)
    @test isapprox(cliped_frame_array.coordinates[2, 2].x, trj.coordinates[2, 3].x, atol = 1e-3)
    @test cliped_frame_array.attributes  == trj.attributes

    cliped_frame_range = clip_trajectory(1:2:3, trj)
    @test isapprox(cliped_frame_range.coordinates[2, 2].y, trj.coordinates[2, 3].y, atol = 1e-3)
    @test cliped_frame_range.attributes  == trj.attributes

    cliped_atom_array = clip_trajectory([2, 3], trj, query_key = :atom)
    @test isapprox(cliped_atom_array.coordinates[2, 2].z,  trj.coordinates[3, 2].z, atol = 1e-3)
    @test cliped_atom_array.attributes  == trj.attributes[[2, 3]]

    cliped_atom_range = clip_trajectory(1:2:3, trj, query_key = :atom)
    @test isapprox(cliped_atom_range.coordinates[2, 2].x,  trj.coordinates[3, 2].x, atol = 1e-3)
    @test cliped_atom_range.attributes  == trj.attributes[1:2:3]
end

@testset "geometric_center_of_mass" begin
    coordinates = [Coordinate(0.0, 1.0, 2.0),
                   Coordinate(1.0, 0.0, 3.0),
                   Coordinate(2.0, 2.0, 4.0)]
    @test isapprox(Array(geometric_center_of_mass(coordinates)), [1.0, 1.0, 3.0], atol = 1e-3)
    @test isapprox(Array(geometric_center_of_mass(coordinates, atom_indices = 1:2:3)), [1.0, 1.5, 3.0], atol = 1e-3)
end

@testset "center_of_mass for Coordinate vector" begin
    coord_vec = [Coordinate(11.1f0, 11.2f0, 11.3f0),
                 Coordinate(21.1f0, 21.2f0, 21.3f0),
                 Coordinate(31.1f0, 31.2f0, 31.3f0)]
    mass_vec  = [2.0f0, 3.0f0, 4.0f0]

    @test isapprox(Array(center_of_mass(coord_vec, mass_vec)), [23.3222f0, 23.4222f0, 23.5222f0],
                   atol=1e-3)
end

@testset "center_of_mass for trajectory" begin
    coordinates = prepare_coordinates(3)
    trj_wo_attr = Trajectory(coordinates)

    @test_throws ArgumentError center_of_mass(trj_wo_attr)

    attributes = [Attribute(mass = 2.0f0), Attribute(mass = 3.0f0), Attribute(mass = 4.0f0)]
    trj = Trajectory(coordinates, attributes)

    @test isapprox(Array(center_of_mass(trj)[1]), [23.3222f0, 23.4222f0, 23.5222f0], atol = 1e-3)
    @test isapprox(Array(center_of_mass(trj, atom_indices = 1:2:3)[2]), [25.4333f0, 25.5333f0, 25.6333f0],
                   atol = 1e-3)
    @test isapprox(Array(center_of_mass(trj, atom_indices = 1:2:3, frame_indices = 2:3)[1]),
                   [25.4333f0, 25.5333f0, 25.6333f0], atol = 1e-3)

    @test isapprox(Array(center_of_mass(trj, geometric = true)[1]), [21.1f0, 21.2f0, 21.3f0], atol = 1e-3)
    @test isapprox(Array(center_of_mass(trj, geometric = true, atom_indices = 1:2:3)[2]),
                   [22.1f0, 22.2f0, 22.3f0], atol = 1e-3)
    @test isapprox(Array(center_of_mass(trj, geometric = true, atom_indices = 1:2:3, frame_indices = 2:3)[1]),
                   [22.1f0, 22.2f0, 22.3f0], atol = 1e-3)

end

@testset "pair_length_matrix" begin
    coordinates = prepare_coordinates(3)
    length_matrix = pair_length_matrix(coordinates[:, 1], coordinates[:, 3])

    @test size(length_matrix) == (3, 3)
    @test isapprox(38.1051f0, length_matrix[1, 3], atol = 1e-3)

    trj = Trajectory(coordinates)
    length_matrices = pair_length_matrix(trj, frame_indices = 1:2,
                                         first_atom_indices = 1:2:3,
                                         second_atom_indices = :)

    @test size(length_matrices[1]) == (2, 3)
    @test isapprox(34.6410f0, length_matrices[2][2, 1], atol = 1e-3)
end

@testset "pair_length_matrix_pbc" begin
    coordinates = prepare_coordinates(3)
    box_size    = Coordinate(40.0, 40.0, 40.0)
    length_matrix = pair_length_matrix_pbc(coordinates[:, 1], coordinates[:, 3], box_size)

    @test size(length_matrix) == (3, 3)
    @test isapprox(31.1769, length_matrix[1, 3], atol = 1e-3)

    trj = Trajectory(coordinates)
    length_matrices = pair_length_matrix_pbc(trj, box_size,
                                             frame_indices = 1:2:3,
                                             first_atom_indices = 1:2:3,
                                             second_atom_indices = :)

    @test size(length_matrices[1]) == (2, 3)
    @test isapprox(distance_pbc(coordinates[3, 3], coordinates[1, 3], box_size),
                   length_matrices[2][2, 1], atol = 1e-3)
end

@testset "pair_length_matrix_parallel" begin
    coordinates = prepare_coordinates(3)
    trj = Trajectory(coordinates)

    @test isapprox(pair_length_matrix(trj), pair_length_matrix_parallel(trj), atol = 1e-3)
end

@testset "contact_bool_matrix" begin
    coordinates = prepare_coordinates(2)
    bool_matrix = contact_bool_matrix(20.0f0, coordinates[:, 1], coordinates[:, 2])

    @test size(bool_matrix) == (3, 3)
    @test bool_matrix == [true true false; true true true; false true true]

    trj = Trajectory(coordinates)
    contact_bool_mat_arr = contact_bool_matrix(18.0f0, trj, first_atom_indices = 1:2, second_atom_indices = 2:3)

    @test contact_bool_mat_arr[1][1, 1] == true
    @test contact_bool_mat_arr[2][1, 2] == false
end

@testset "contact_bool_matrix_pbc" begin
    first_coord  = Coordinate(1.0, 2.0, 2.0)
    second_coord = Coordinate(2.5, 2.0, 2.0)
    third_coord  = Coordinate(5.0, 2.0, 2.0)
    coordinates  = [first_coord  first_coord ;
                    second_coord second_coord;
                    third_coord  third_coord ]
    trj          = Trajectory(coordinates)
    box_size     = Coordinate(5.5, 4.0, 4.0)

    trj = Trajectory(coordinates)
    contact_bool_mat_arr = contact_bool_matrix_pbc(2.0f0, trj, box_size)

    @test contact_bool_mat_arr[1][1, 2] == true
    @test contact_bool_mat_arr[2][3, 1] == true
    @test contact_bool_mat_arr[2][2, 3] == false
end

@testset "contact_bool_matrix_parallel" begin
    coordinates = prepare_coordinates(3)
    trj = Trajectory(coordinates)

    @test isapprox(contact_bool_matrix(18.0f0, trj), contact_bool_matrix_parallel(18.0f0, trj), atol = 1e-3)
end

@testset "contact_probability_matrix" begin
    coordinates = prepare_coordinates(3)
    trj = Trajectory(coordinates)
    contact_prob_mat = contact_probability_matrix(18.0f0, trj,
                                                  frame_indices = 1:2:3,
                                                  first_atom_indices = 1:2,
                                                  second_atom_indices = 2:3)

    @test size(contact_prob_mat) == (2, 2)
    @test isapprox(contact_prob_mat[1, 1], 1.0f0, atol = 1e-3)
    @test isapprox(contact_prob_mat[1, 2], 0.0f0, atol = 1e-3)
end

@testset "contact_probability_matrix_parallel" begin
    coordinates = prepare_coordinates(3)
    trj = Trajectory(coordinates)

    @test isapprox(contact_probability_matrix(18.0f0, trj), contact_probability_matrix_parallel(18.0f0, trj), atol = 1e-3)
end

@testset "radius_of_gyration" begin
    coordinates = prepare_coordinates(3)
    attributes = [Attribute(mass = 2.0f0), Attribute(mass = 3.0f0), Attribute(mass = 4.0f0)]
    trj = Trajectory(coordinates, attributes)

    rg_arr     = Array{Float32, 1}(undef, 3)
    geo_rg_arr = Array{Float32, 1}(undef, 3)
    for frame_idx in 1:size(coordinates, 2)
        # com calculation
        com     = Coordinate([0.0, 0.0, 0.0])
        geo_com = Coordinate([0.0, 0.0, 0.0])
        for atom_idx in 1:size(coordinates, 1)
            com     += coordinates[atom_idx, frame_idx] * attributes[atom_idx].mass
            geo_com += coordinates[atom_idx, frame_idx]
        end
        com     /= sum([attr.mass for attr in attributes])
        geo_com /= size(coordinates, 1)

        # rg calculation
        rg     = 0.0
        geo_rg = 0.0
        for atom_idx in 1:size(coordinates, 1)
            rg     += norm(coordinates[atom_idx, frame_idx] - com)^2 * attributes[atom_idx].mass
            geo_rg += norm(coordinates[atom_idx, frame_idx] - geo_com)^2
        end
        rg     = sqrt(rg     / sum([attr.mass for attr in attributes]))
        geo_rg = sqrt(geo_rg / size(coordinates, 1))
        rg_arr[frame_idx]     = rg
        geo_rg_arr[frame_idx] = geo_rg
    end

    @test isapprox(radius_of_gyration(trj, frame_indices = 1:2:3),                   rg_arr[1:2:3],     atol = 1e-3)
    @test isapprox(radius_of_gyration(trj, frame_indices = 1:2:3, geometric = true), geo_rg_arr[1:2:3], atol = 1e-3)
end

@testset "atom_mass" begin
    @test isapprox(atom_mass("HB3"), 1.00798, atol = 1e-5)
    @test isapprox(atom_mass("CA"),  12.0106, atol = 1e-5)
    @test isapprox(atom_mass("NB"),  14.0069, atol = 1e-5)
    @test isapprox(atom_mass("NA"),  22.9898, atol = 1e-5)
    @test isapprox(atom_mass("CL"),  35.4520, atol = 1e-5)
end

@testset "residue_mass" begin
    @test isapprox(residue_mass("ALA"), 71.03711, atol = 1e-5)
    @test isapprox(residue_mass("HOH"), 18.01528, atol = 1e-5)
end

@testset "distance_pbc" begin
    first_coord  = Coordinate(1.0, 2.0, 2.0)
    second_coord = Coordinate(2.5, 2.0, 2.0)
    third_coord  = Coordinate(3.5, 2.0, 2.0)
    box_size     = Coordinate(4.0, 4.0, 4.0)
    first_second =
        distance_pbc(first_coord, second_coord, box_size)
    first_third  =
        distance_pbc(first_coord, third_coord,  box_size)
    @test isapprox(first_second, 1.5, atol = 1e-3)
    @test isapprox(first_third,  1.5, atol = 1e-3)
end

@testset "fix_pbc" begin
    dcd = read_dcd("data/test_position.dcd")
    groupid_vec = [1, 1, 2, 2, 2]
    box_size  = Coordinate(10.0f0, 10.0f0, 10.0f0)
    @test_nowarn fix_pbc(dcd.coordinates[:, 1], groupid_vec, box_size)
    @test_nowarn fix_pbc(dcd.coordinates[:, 1], groupid_vec, box_size, z = false)

    pdb = read_pdb("data/test_position.pdb")
    @test_nowarn fix_pbc(pdb, groupid_vec, box_size)
    @test_nowarn fix_pbc(pdb, groupid_vec, box_size, z = false)
end

@testset "fix_pbc!" begin
    dcd = read_dcd("data/test_position.dcd")
    pdb = read_pdb("data/test_position.pdb")
    dcd.attributes = pdb.attributes
    box_size = Coordinate(10.0, 10.0, 10.0)
    @test_nowarn fix_pbc!(dcd, box_size)
    @test_nowarn fix_pbc!(dcd, box_size, z = false)
end

@testset "move_pbc_center" begin
    first_coord  = Coordinate(1.0, 2.0, 2.0)
    second_coord = Coordinate(2.5, 2.0, 2.0)
    third_coord  = Coordinate(3.5, 2.0, 2.0)
    coordinates  = [first_coord, second_coord, third_coord]
    new_center   = Coordinate(1.0, 2.0, 2.0)
    box_size     = Coordinate(4.0, 4.0, 4.0)
    new_coords = move_pbc_center(coordinates, new_center, box_size)
    @test isapprox(Array(first_coord), Array(new_coords[1]), atol = 1e-3)
    @test isapprox([-0.5, 2.0, 2.0],   Array(new_coords[3]), atol = 1e-3)
end


text:
# Test skipchars for IOStreams
mktemp() do path, file
    function append_to_file(str)
        mark(file)
        println(file, str)
        flush(file)
        reset(file)
    end

    # test it doesn't error on eof
    @test skipchars(file, isspace) == file

    # test if it correctly skips
    append_to_file("    ")
    @test eof(skipchars(file, isspace))

    # test it correctly detects comment lines
    append_to_file("#    \n   ")
    @test eof(skipchars(file, isspace, linecomment='#'))

    # test it stops at the appropriate time
    append_to_file("   not a space")
    @test skipchars(file, isspace) == file
    @test !eof(file) && read(file, Char) == 'n'

    # test it correctly ignores the contents of comment lines
    append_to_file("  #not a space \n   not a space")
    @test skipchars(file, isspace, linecomment='#') == file
    @test !eof(file) && read(file, Char) == 'n'
end


text:
using jPCA
using Base.Test

# main
X = rand(100,100)
Y = fit(JPCA, X, 10)

@test size(Y.M,1) == size(Y.M,2)

@test_throws ErrorException fit(JPCA, rand(5,5), 6)
@test_throws ErrorException fit(JPCA, rand(5,10))

# skew sym
import jPCA.reshape_skew
# create a random skew sym matrix
Z = Z = rand(100,100)
Z .= 0.5(Z - Z')
@test Z == reshape_skew(reshape_skew(Z))

@test_throws ErrorException reshape_skew(Z[:,1:99])
Z[1,1] += 0.1
@test_throws ErrorException reshape_skew(Z)

max_stars_repo_name: alexmorley/jPCA

max_issues_repo_name: alexmorley/jPCA

max_forks_repo_name: alexmorley/jPCA


text:
using BernoulliFactory
using Random
import Random.GLOBAL_RNG
import Statistics: mean, var

μ() = rand()*3*π
m = 100000

vs = Vector{Float64}(undef, m)
flips = Vector{Int64}(undef, m)
calls = Vector{Int64}(undef, m)
for i in 1:m
  vs[i], flips[i], calls[i] = BernoulliFactory.signedEstimate(μ, sin, 1.0, 0.1, 1)
end
println("actual value = ", 2/(3*π))
println("estimated value = ", mean(vs))
println("variance of estimates = ", var(vs))
println("should be approx. standard normal: ", (mean(vs)-2/(3*π))/sqrt(var(vs)/m))


text:
"""
    blosc2_meta_exists(schunk::Ptr{blosc2_schunk}, name)

Find whether the schunk has a metalayer or not.

# Arguments

- `schunk` – The super-chunk from which the metalayer will be checked.

- `name` – The name of the metalayer to be checked.

# Returns

If successful, return the index of the metalayer. Else, return a negative value.
"""
function blosc2_meta_exists(schunk::Ptr{blosc2_schunk}, name)
    @ccall lib.blosc2_meta_exists(schunk::Ptr{blosc2_schunk}, name::Cstring)::Cint
end

"""
    blosc2_meta_add(schunk::Ptr{blosc2_schunk}, name, content, content_len)

Add content into a new metalayer.

# Arguments

- `schunk` – The super-chunk to which the metalayer should be added.

- `name` – The name of the metalayer.

- `content` – The content of the metalayer.

- `content_len` – The length of the content.

# Returns

If successful, the index of the new metalayer. Else, return a negative value.
"""
function blosc2_meta_add(schunk::Ptr{blosc2_schunk}, name, content, content_len)
    @ccall lib.blosc2_meta_add(schunk::Ptr{blosc2_schunk}, name::Cstring, content::Ptr{UInt8}, content_len::UInt32)::Cint
end

"""
    blosc2_meta_update(schunk::Ptr{blosc2_schunk}, name, content, content_len)

Update the content of an existing metalayer.

# Arguments

- `schunk` – The frame containing the metalayer.

- `name` – The name of the metalayer to be updated.

- `content` – The new content of the metalayer.

- `content_len` – The length of the content.

# Returns

If successful, the index of the metalayer. Else, return a negative value.

# Note

Contrarily to `blosc2_meta_add` the updates to metalayers are automatically serialized into a possible attached frame.
"""
function blosc2_meta_update(schunk::Ptr{blosc2_schunk}, name, content, content_len)
    @ccall lib.blosc2_meta_update(schunk::Ptr{blosc2_schunk}, name::Cstring, content::Ptr{UInt8}, content_len::UInt32)::Cint
end

"""
    blosc2_meta_get(schunk::Ptr{blosc2_schunk}, name)

Get the content out of a  metalayer.

# Arguments

- `schunk` – The frame containing the metalayer.

- `name` – The name of the metalayer to be updated.

- `content` – The new content of the metalayer.

- `content_len` – The length of the content.

# Returns

If successful, the content of the metalayer. Else, return `nothing`.
"""
function blosc2_meta_get(schunk::Ptr{blosc2_schunk}, name)
   content_ptr = Ref{Ptr{UInt8}}()
   content_len = Ref{UInt32}()
   r = @ccall lib.blosc2_meta_get(schunk::Ptr{blosc2_schunk}, name::Cstring, content_ptr::Ptr{Ptr{UInt8}}, content_len::Ref{UInt32})::Cint
   r < 0 && return nothing
   return unsafe_wrap(Array, content_ptr[], content_len[], own = true)
end
"""
    blosc2_vmeta_exists(schunk::Ptr{blosc2_schunk}, name)

Find whether the schunk has a variable-length metalayer or not.

# Arguments

- `schunk` – The super-chunk from which the metalayer will be checked.

- `name` – The name of the metalayer to be checked.

# Returns

If successful, return the index of the metalayer. Else, return a negative value.
"""
function blosc2_vlmeta_exists(schunk::Ptr{blosc2_schunk}, name)
    @ccall lib.blosc2_vlmeta_exists(schunk::Ptr{blosc2_schunk}, name::Cstring)::Cint
end

"""
    blosc2_vmeta_add(schunk::Ptr{blosc2_schunk}, name, content, content_len, cparams::blosc2_cparams)

Add content into a new variable-length metalayer.

# Arguments

- `schunk` – The super-chunk to which the metalayer should be added.

- `name` – The name of the metalayer.

- `content` – The content of the metalayer.

- `content_len` – The length of the content.

- `cparams` – The parameters for compressing the variable-length metalayer content.

# Returns

If successful, the index of the new metalayer. Else, return a negative value.
"""
function blosc2_vlmeta_add(schunk::Ptr{blosc2_schunk}, name, content, content_len, cparams::blosc2_cparams = blosc2_cparams())
    cparams_ref = Ref(cparams)
    @ccall lib.blosc2_vlmeta_add(
        schunk::Ptr{blosc2_schunk},
        name::Cstring,
        content::Ptr{UInt8},
        content_len::UInt32,
        cparams_ref::Ptr{blosc2_cparams}
    )::Cint
end

"""
    blosc2_vlmeta_update(schunk::Ptr{blosc2_schunk}, name, content, content_len)

Update the content of an existing variable-length metalayer.

# Arguments

- `schunk` – The frame containing the metalayer.

- `name` – The name of the metalayer to be updated.

- `content` – The new content of the metalayer.

- `content_len` – The length of the content.

- `cparams` – The parameters for compressing the variable-length metalayer content.

# Returns

If successful, the index of the metalayer. Else, return a negative value.

# Note

Contrarily to `blosc2_meta_add` the updates to metalayers are automatically serialized into a possible attached frame.
"""
function blosc2_vlmeta_update(schunk::Ptr{blosc2_schunk}, name, content, content_len, cparams::blosc2_cparams = blosc2_cparams())
    cparams_ref = Ref(cparams)
    @ccall lib.blosc2_vlmeta_update(
        schunk::Ptr{blosc2_schunk},
        name::Cstring,
        content::Ptr{UInt8},
        content_len::UInt32,
        cparams_ref::Ptr{blosc2_cparams}
    )::Cint
end

"""
    blosc2_vlmeta_get(schunk::Ptr{blosc2_schunk}, name)

Get the content out of a variable-length metalayer.

# Arguments

- `schunk` – The frame containing the metalayer.

- `name` – The name of the metalayer to be updated.

- `content` – The new content of the metalayer.

- `content_len` – The length of the content.

# Returns

If successful, the content of the metalayer. Else, return `nothing`.
"""
function blosc2_vlmeta_get(schunk::Ptr{blosc2_schunk}, name)
   content_ptr = Ref{Ptr{UInt8}}()
   content_len = Ref{UInt32}()
   r = @ccall lib.blosc2_vlmeta_get(schunk::Ptr{blosc2_schunk}, name::Cstring, content_ptr::Ptr{Ptr{UInt8}}, content_len::Ref{UInt32})::Cint
   r < 0 && return nothing
   return unsafe_wrap(Array, content_ptr[], content_len[], own = true)
end


text:
using Random
using PlotUtils
using Makie
using AbstractPlotting
import Showoff
using Printf
using AbstractPlotting: Rect2D
import AbstractPlotting: IRect2D

const BBox = Rect2D{Float64}

left(rect::Rect2D) = minimum(rect)[1]
right(rect::Rect2D) = maximum(rect)[1]

bottom(rect::Rect2D) = minimum(rect)[2]
top(rect::Rect2D) = maximum(rect)[2]

abstract type Side end

struct Left <: Side end
struct Right <: Side end
struct Top <: Side end
struct Bottom <: Side end

Base.getindex(bbox::Rect2D, ::Left) = left(bbox)
Base.getindex(bbox::Rect2D, ::Right) = right(bbox)
Base.getindex(bbox::Rect2D, ::Bottom) = bottom(bbox)
Base.getindex(bbox::Rect2D, ::Top) = top(bbox)

mutable struct LayoutedAxis
    parent::Scene
    scene::Scene
    xlabel::Node{String}
    ylabel::Node{String}
    limits::Node{FRect2D} # these should steer the camera, not used yet
end

width(rect::Rect2D) = right(rect) - left(rect)
height(rect::Rect2D) = top(rect) - bottom(rect)

function BBox(left::Number, right::Number, top::Number, bottom::Number)
    mini = (left, bottom)
    maxi = (right, top)
    return BBox(mini, maxi .- mini)
end

function IRect2D(bbox::Rect2D)
    return IRect2D(
        round.(Int, minimum(bbox)),
        round.(Int, widths(bbox))
    )
end

struct RowCols{T <: Union{Number, Vector{Float64}}}
    lefts::T
    rights::T
    tops::T
    bottoms::T
end

function RowCols(ncols::Int, nrows::Int)
    return RowCols(
        zeros(ncols),
        zeros(ncols),
        zeros(nrows),
        zeros(nrows)
    )
end

Base.getindex(rowcols::RowCols, ::Left) = rowcols.lefts
Base.getindex(rowcols::RowCols, ::Right) = rowcols.rights
Base.getindex(rowcols::RowCols, ::Top) = rowcols.tops
Base.getindex(rowcols::RowCols, ::Bottom) = rowcols.bottoms

"""
    eachside(f)
Calls f over all sides (Left, Right, Top, Bottom), and creates a BBox from the result of f(side)
"""
function eachside(f)
    return BBox(map(f, (Left(), Right(), Top(), Bottom()))...)
end

"""
mapsides(
       f, first::Union{Rect2D, RowCols}, rest::Union{Rect2D, RowCols}...
   )::BBox
Maps f over all sides of the rectangle like arguments.
e.g.
```
mapsides(BBox(left, right, top, bottom)) do side::Side, side_val::Number
    return ...
end::BBox
```
"""
function mapsides(
        f, first::Union{Rect2D, RowCols}, rest::Union{Rect2D, RowCols}...
    )
    return eachside() do side
        f(side, getindex.((first, rest...), (side,))...)
    end
end

abstract type Alignable end

"""
Used to specify space that is occupied in a grid. Like 1:1|1:1 for the first square,
or 2:3|1:4 for a rect over the 2nd and 3rd row and the first four columns.
"""
struct Span
    rows::UnitRange{Int64}
    cols::UnitRange{Int64}
end

"""
An object that can be aligned that also specifies how much space it occupies in
a grid via its span.
"""
struct SpannedAlignable{T <: Alignable}
    al::T
    sp::Span
end
function side_indices(c::SpannedAlignable)
    return RowCols(
        c.sp.cols.start,
        c.sp.cols.stop,
        c.sp.rows.start,
        c.sp.rows.stop,
    )
end

"""
These functions tell whether an object in a grid touches the left, top, etc. border
of the grid. This means that it is relevant for the grid's own protrusion on that side.
"""
ismostin(sp::SpannedAlignable, grid, ::Left) = sp.sp.cols.start == 1
ismostin(sp::SpannedAlignable, grid, ::Right) = sp.sp.cols.stop == grid.ncols
ismostin(sp::SpannedAlignable, grid, ::Bottom) = sp.sp.rows.stop == grid.nrows
ismostin(sp::SpannedAlignable, grid, ::Top) = sp.sp.cols.start == 1

isleftmostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Left())
isrightmostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Right())
isbottommostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Bottom())
istopmostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Top())

struct SolvedAxisLayout <: Alignable
    inner::BBox
    outer::BBox
    axis::LayoutedAxis
end

struct AxisLayout <: Alignable
    decorations::BBox
    axis::LayoutedAxis
end

struct SolvedGridLayout <: Alignable
    bbox::BBox
    content::Vector{SpannedAlignable}
    nrows::Int
    ncols::Int
    grid::RowCols{Vector{Float64}}
end

struct GridLayout <: Alignable
    content::Vector{SpannedAlignable}
    nrows::Int
    ncols::Int
    colratios::Vector{Float64}
    rowratios::Vector{Float64}
    colgapfraction::Float64
    rowgapfraction::Float64
end

"""
All the protrusion functions calculate how much stuff "sticks out" of a layoutable object.
This is so that collisions are avoided, while what is actually aligned is the
"important" edges of the layout objects.
"""
leftprotrusion(x) = protrusion(x, Left())
rightprotrusion(x) = protrusion(x, Right())
bottomprotrusion(x) = protrusion(x, Bottom())
topprotrusion(x) = protrusion(x, Top())

protrusion(u::AxisLayout, side::Side) = u.decorations[side]
protrusion(sp::SpannedAlignable, side::Side) = protrusion(sp.al, side)

function protrusion(gl::GridLayout, side::Side)
    return mapreduce(max, gl.content, init = 0.0) do elem
        # we use only objects that stick out on this side
        # And from those we use the maximum protrusion
        ismostin(elem, gl, side) ? protrusion(elem, side) : 0.0
    end
end

protrusion(s::SolvedAxisLayout, ::Left) = left(s.inner) - left(s.outer)
function protrusion(s::SolvedAxisLayout, side::Side)
    return s.outer[side] - s.inner[side]
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 solve(gl::GridLayout, bbox::BBox)

    # 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.al, side))
        end
    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
    maxcolgap = maximum(colgaps)
    maxrowgap = maximum(rowgaps)

    # 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 = maxcolgap * (gl.ncols - 1)
    sumrowgaps = maxrowgap * (gl.nrows - 1)

    # compute what space remains for the inner parts of the plots
    remaininghorizontalspace = width(bbox) - sumcolgaps
    remainingverticalspace = height(bbox) - sumrowgaps

    # 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
    addedcolgap = gl.colgapfraction * remaininghorizontalspace
    addedrowgap = gl.rowgapfraction * remainingverticalspace

    # compute the actual space available for the rows and columns (plots without protrusions)
    spaceforcolumns = remaininghorizontalspace - addedcolgap * (gl.ncols - 1)
    spaceforrows = remainingverticalspace - addedrowgap * (gl.nrows - 1)

    # compute the column widths and row heights using the specified row and column ratios
    colwidths = gl.colratios ./ sum(gl.colratios) .* spaceforcolumns
    rowheights = gl.rowratios ./ sum(gl.rowratios) .* spaceforrows

    # this is the vertical / horizontal space between the inner lines of all plots
    colgap = maxcolgap + addedcolgap
    rowgap = maxrowgap + addedrowgap

    # compute the x values for all left and right column boundaries
    xleftcols = [left(bbox) + sum(colwidths[1:i-1]) + (i - 1) * colgap for i in 1:gl.ncols]
    xrightcols = xleftcols .+ colwidths

    # compute the y values for all top and bottom row boundaries
    ytoprows = [top(bbox) - sum(rowheights[1:i-1]) - (i - 1) * rowgap for i in 1:gl.nrows]
    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
    )
    solvedcontent = map(gl.content) do c
        idx_rect = side_indices(c)
        bbox_cell = mapsides(idx_rect, gridboxes) do side, idx, gridside
            gridside[idx]
        end
        solved = solve(c.al, bbox_cell)
        return SpannedAlignable(solved, c.sp)
    end
    # return a solved grid layout in which all objects are also solved layout objects
    return SolvedGridLayout(
        bbox, solvedcontent,
        gl.nrows, gl.ncols,
        gridboxes
    )
end

"""
This function solves a grid layout so that it fits exactly inside a bounding box.
Exactly means that the protrusions of all other objects inside this grid layout
also have to fit into the bounding box. This is needed if the grid is the outermost
object in the layout, the bounding box would then be the scene boundary.
"""
function outersolve(gl::GridLayout, bbox::BBox)
    maxgrid = RowCols(gl.ncols, gl.nrows)
    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.al, side))
        end
    end

    topprot = maxgrid.tops[1]
    bottomprot = maxgrid.bottoms[end]
    leftprot = maxgrid.lefts[1]
    rightprot = maxgrid.rights[end]

    colgaps = maxgrid.lefts[2:end] .+ maxgrid.rights[1:end-1]
    rowgaps = maxgrid.tops[2:end] .+ maxgrid.bottoms[1:end-1]

    maxcolgap = gl.ncols <= 1 ? 0 : maximum(colgaps)
    maxrowgap = gl.nrows <= 1 ? 0 : maximum(rowgaps)

    sumcolgaps = maxcolgap * (gl.ncols - 1)
    sumrowgaps = maxrowgap * (gl.nrows - 1)

    remaininghorizontalspace = width(bbox) - sumcolgaps - leftprot - rightprot
    remainingverticalspace = height(bbox) - sumrowgaps - topprot - bottomprot

    addedcolgap = gl.colgapfraction * remaininghorizontalspace
    addedrowgap = gl.rowgapfraction * remainingverticalspace

    spaceforcolumns = remaininghorizontalspace - addedcolgap * (gl.ncols - 1)
    spaceforrows = remainingverticalspace - addedrowgap * (gl.nrows - 1)

    colwidths = gl.colratios ./ sum(gl.colratios) .* spaceforcolumns
    rowheights = gl.rowratios ./ sum(gl.rowratios) .* spaceforrows

    colgap = maxcolgap + addedcolgap
    rowgap = maxrowgap + addedrowgap

    xleftcols = map(1:gl.ncols) do i
        left(bbox) + leftprot + sum(colwidths[1:i-1]) + (i - 1) * colgap
    end
    xrightcols = xleftcols .+ colwidths

    ytoprows = map(1:gl.nrows) do i
        top(bbox) - topprot - sum(rowheights[1:i-1]) - (i - 1) * rowgap
    end

    ybottomrows = ytoprows .- rowheights
    gridboxes = RowCols(
        xleftcols, xrightcols,
        ytoprows, ybottomrows
    )
    solvedcontent = map(gl.content) do c
        idx_rect = side_indices(c)
        bbox_cell = mapsides(idx_rect, gridboxes) do side, idx, gridside
            gridside[idx]
        end
        solved = solve(c.al, bbox_cell)
        return SpannedAlignable(solved, c.sp)
    end
    # solvedcontent = solve.(gl.content)
    SolvedGridLayout(
        bbox, solvedcontent, gl.nrows, gl.ncols,
        gridboxes
    )
end

function solve(ua::AxisLayout, innerbbox)
    bbox = mapsides(innerbbox, ua.decorations) do side, iside, decside
        op = side isa Union{Left, Top} ? (-) : (+)
        return op(iside, decside)
    end
    SolvedAxisLayout(innerbbox, bbox, ua.axis)
end

const Indexables = Union{UnitRange, Int, Colon}

"""
This function allows indexing syntax to add a layout object to a grid.
You can do:

grid[1, 1] = obj
grid[1, :] = obj
grid[1:3, 2:5] = obj

and all combinations of the above
"""
function Base.setindex!(g, a::Alignable, rows::Indexables, cols::Indexables)
    if rows isa Int
        rows = rows:rows
    elseif rows isa Colon
        rows = 1:g.nrows
    end
    if cols isa Int
        cols = cols:cols
    elseif cols isa Colon
        cols = 1:g.ncols
    end

    if !((1 <= rows.start <= g.nrows) || (1 <= rows.stop <= g.nrows))
        error("invalid row span $rows for grid with $(g.nrows) rows")
    end
    if !((1 <= cols.start <= g.ncols) || (1 <= cols.stop <= g.ncols))
        error("invalid col span $cols for grid with $(g.ncols) columns")
    end

    push!(g.content, SpannedAlignable(a, Span(rows, cols)))
end

function axislines!(scene, rect)
    points = lift(rect) do r
        p1 = Point2(r.origin[1], r.origin[2] + r.widths[2])
        p2 = Point2(r.origin[1], r.origin[2])
        p3 = Point2(r.origin[1] + r.widths[1], r.origin[2])
        [p1, p2, p3]
    end
    lines!(scene, points, linewidth = 2, show_axis = false)
end

"""
A cheaper function that tries to come up with usable tick locations for a given value range
"""
function locateticks(xmin, xmax)
    # whats the distance?
    d = xmax - xmin
    # which order of magnitude is the distance in?
    ex = log10(d)
    # round the exponent to the closest one
    exrounded = round(ex)
    # this factor is used so we can find integer steps that then map back to
    # nice numbers in the given value range
    factor = 1 / 10 ^ (exrounded - 1)

    # minimum needs to be at least the lower value times the scaling factor
    xminf = ceil(xmin * factor)
    # maximum needs to be at most the higher value times the scaling factor
    xmaxf = floor(xmax * factor)

    # xminf and xmaxf are now both integers that are in an order of magnitude around ten steps apart
    df = xmaxf - xminf

    # step sizes we like
    steps = [5, 4, 2, 1]

    # from the highest to the lowest step size, choose the first that fits at least
    # two times between the end values (gives three ticks including the end value)
    for s in steps
        n, remainder = divrem(df, s)
        if n >= 2
            rang = 1:n
            # calculate the ticks by dividing with the factor from above
            ticks = [xminf; [xminf + x * s for x in rang]] ./ factor
            return ticks
        end
    end
end

function LayoutedAxis(parent::Scene)
    scene = Scene(parent, Node(IRect(0, 0, 100, 100)))
    limits = Node(FRect(0, 0, 100, 100))
    xlabel = Node("x label")
    ylabel = Node("y label")

    cam = cam2d!(scene)

    ticksnode = Node(Point2f0[])
    ticks = linesegments!(
        parent, ticksnode, linewidth = 2, show_axis = false
    )[end]

    # the algorithm from above seems to not give more than 7 ticks with the step sizes I chose
    nmaxticks = 7

    xticklabelnodes = [Node("0") for i in 1:nmaxticks]
    xticklabelposnodes = [Node(Point(0.0, 0.0)) for i in 1:nmaxticks]
    xticklabels = map(1:nmaxticks) do i
        text!(
            parent,
            xticklabelnodes[i],
            position = xticklabelposnodes[i],
            align = (:center, :top),
            textsize = 20,
            show_axis = false
        )[end]
    end

    yticklabelnodes = [Node("0") for i in 1:nmaxticks]
    yticklabelposnodes = [Node(Point(0.0, 0.0)) for i in 1:nmaxticks]
    yticklabels = map(1:nmaxticks) do i
        text!(
            parent,
            yticklabelnodes[i],
            position = yticklabelposnodes[i],
            align = (:center, :bottom),
            rotation = pi/2,
            textsize = 20,
            show_axis = false
        )[end]
    end

    on(cam.area) do a

        pxa = scene.px_area[]
        px_aspect = pxa.widths[1] / pxa.widths[2]

        # @printf("cam %.1f, %.1f, %.1f, %.1f\n", a.origin..., a.widths...)
        # @printf("pix %.1f, %.1f, %.1f, %.1f\n", pxa.origin..., pxa.widths...)

        width = px_aspect > 1 ? a.widths[1] / px_aspect : a.widths[1]
        xrange = (a.origin[1], a.origin[1] + width)

        if width == 0 || !isfinite(xrange[1]) || !isfinite(xrange[2])
            return
        end

        xtickvals = locateticks(xrange...)
        # xtickvals, vminbest, vmaxbest = optimize_ticks(xrange...)

        # this code here tries to transform between values given in pixels of the
        # scene and the camera area, but this is incorrect right now and everything
        # should just be determined by the now unused limits that are saved in the
        # LayoutedAxis object
        xfractions = (xtickvals .- xrange[1]) ./ width
        xrange_scene = (pxa.origin[1], pxa.origin[1] + pxa.widths[1])
        width_scene = xrange_scene[2] - xrange_scene[1]
        xticks_scene = xrange_scene[1] .+ width_scene .* xfractions
        ticksize = 10 # px
        y = pxa.origin[2]
        xtickstarts = [Point(x, y) for x in xticks_scene]
        xtickends = [t + Point(0.0, -ticksize) for t in xtickstarts]

        height = px_aspect < 1 ? a.widths[2] * px_aspect : a.widths[2]
        yrange = (a.origin[2], a.origin[2] + height)

        ytickvals = locateticks(yrange...)
        # ytickvals, vminbest, vmaxbest = optimize_ticks(yrange...)
        yfractions = (ytickvals .- yrange[1]) ./ height
        yrange_scene = (pxa.origin[2], pxa.origin[2] + pxa.widths[2])
        height_scene = yrange_scene[2] - yrange_scene[1]
        yticks_scene = yrange_scene[1] .+ height_scene .* yfractions
        ticksize = 10 # px
        x = pxa.origin[1]
        ytickstarts = [Point(x, y) for y in yticks_scene]
        ytickends = [t + Point(-ticksize, 0.0) for t in ytickstarts]

        # set and position tick labels
        xtickstrings = Showoff.showoff(xtickvals, :plain)
        nxticks = length(xtickvals)
        for i in 1:nmaxticks
            if i <= nxticks
                xticklabelnodes[i][] = xtickstrings[i]
                xticklabelposnodes[i][] = xtickends[i] + Point(0.0, -10.0)
                xticklabels[i].visible = true
            else
                xticklabels[i].visible = false
            end
        end

        ytickstrings = Showoff.showoff(ytickvals, :plain)
        nyticks = length(ytickvals)
        for i in 1:nmaxticks
            if i <= nyticks
                yticklabelnodes[i][] = ytickstrings[i]
                yticklabelposnodes[i][] = ytickends[i] + Point(-10.0, 0.0)
                yticklabels[i].visible = true
            else
                yticklabels[i].visible = false
            end
        end

        # set tick mark positions
        ticksnode[] = collect(Iterators.flatten(zip(
            [xtickstarts; ytickstarts],
            [xtickends; ytickends]
        )))
    end

    labelgap = 50

    xlabelpos = lift(scene.px_area) do a
        Point2(a.origin[1] + a.widths[1] / 2, a.origin[2] - labelgap)
    end

    ylabelpos = lift(scene.px_area) do a
        Point2(a.origin[1] - labelgap, a.origin[2] + a.widths[2] / 2)
    end

    tx = text!(
        parent, xlabel, textsize = 20, position = xlabelpos, show_axis = false
    )[end]
    tx.align = [0.5, 1]
    ty = text!(
        parent, ylabel, textsize = 20,
        position = ylabelpos, rotation = pi/2, show_axis = false
    )[end]
    ty.align = [0.5, 0]

    axislines!(parent, scene.px_area)

    LayoutedAxis(parent, scene, xlabel, ylabel, limits)
end

function applylayout(sg::SolvedGridLayout)
    for c in sg.content
        applylayout(c.al)
    end
end

function applylayout(sa::SolvedAxisLayout)
    sa.axis.scene.px_area[] = IRect2D(sa.inner)
end

function shrinkbymargin(rect, margin)
    IRect((rect.origin .+ margin), (rect.widths .- 2 .* margin))
end

begin
    scene = Scene(resolution=(600, 600));
    screen = display(scene)
    campixel!(scene);

    la1 = LayoutedAxis(scene)
    la2 = LayoutedAxis(scene)
    la3 = LayoutedAxis(scene)
    la4 = LayoutedAxis(scene)
    la5 = LayoutedAxis(scene)

    # lines!(la1.scene, rand(200, 2) .* 100, color=:black, show_axis=false)
    img = rand(100, 100)
    image!(la1.scene, img, show_axis=false)
    lines!(la2.scene, rand(200, 2) .* 100, color=:blue, show_axis=false)
    scatter!(la3.scene, rand(200, 2) .* 100, markersize=3, color=:red, show_axis=false)
    lines!(la4.scene, rand(200, 2) .* 100, color=:orange, show_axis=false)
    lines!(la5.scene, rand(200, 2) .* 100, color=:pink, show_axis=false)
    update!(scene)

    gl = GridLayout([], 2, 2, [1, 1], [1, 1], 0.01, 0.01)

    gl[2, 1:2] = AxisLayout(BBox(75, 0, 0, 75), la1)
    gl[1, 2] = AxisLayout(BBox(75, 0, 0, 75), la2)

    gl2 = GridLayout([], 2, 2, [0.8, 0.2], [0.2, 0.8], 0.01, 0.01)
    gl2[2, 1] = AxisLayout(BBox(75, 0, 0, 75), la3)
    gl2[1, 1] = AxisLayout(BBox(75, 0, 0, 75), la4)
    gl2[2, 2] = AxisLayout(BBox(75, 0, 0, 75), la5)

    gl[1, 1] = gl2

    sg = outersolve(gl, BBox(shrinkbymargin(pixelarea(scene)[], 30)))
    applylayout(sg)
    # when the scene is resized, apply the outersolve'd outermost grid layout
    # this recursively updates all layout objects that are contained in the grid
    on(scene.events.window_area) do area
        sg = outersolve(gl, BBox(shrinkbymargin(pixelarea(scene)[], 30)))
        applylayout(sg)
    end
    # @profiler wait(screen)
end

max_stars_repo_name: SimonDanisch/MakieLayout

max_issues_repo_name: SimonDanisch/MakieLayout

max_forks_repo_name: SimonDanisch/MakieLayout