text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of row name column.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Show row names" begin
    expected = """
\\begin{table}
  \\begin{tabular}{rrrrr}
    \\hline\\hline
     & \\textbf{C1} & \\textbf{C2} & \\textbf{C3} & \\textbf{C4} \\\\
     & \\texttt{Int} & \\texttt{Bool} & \\texttt{Float} & \\texttt{Hex} \\\\\\hline
    Row 1 & 1 & false & 1.0 & 1 \\\\
    Row 2 & 2 & true & 2.0 & 2 \\\\
    Row 3 & 3 & false & 3.0 & 3 \\\\
    Row 4 & 4 & true & 4.0 & 4 \\\\
    Row 5 & 5 & false & 5.0 & 5 \\\\
    Row 6 & 6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data;
                          backend = Val(:latex),
                          header = (["C1",  "C2",   "C3",    "C4"],
                                    ["Int", "Bool", "Float", "Hex"]),
                          row_names = ["Row $i" for i = 1:6])

    @test result == expected

    expected = """
\\begin{table}
  \\begin{tabular}{rrrrr}
    \\hline\\hline
    \\textbf{Row names} & \\textbf{C1} & \\textbf{C2} & \\textbf{C3} & \\textbf{C4} \\\\
     & \\texttt{Int} & \\texttt{Bool} & \\texttt{Float} & \\texttt{Hex} \\\\\\hline
    Row 1 & 1 & false & 1.0 & 1 \\\\
    Row 2 & 2 & true & 2.0 & 2 \\\\
    Row 3 & 3 & false & 3.0 & 3 \\\\
    Row 4 & 4 & true & 4.0 & 4 \\\\
    Row 5 & 5 & false & 5.0 & 5 \\\\
    Row 6 & 6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data;
                          backend = Val(:latex),
                          header = (["C1",  "C2",   "C3",    "C4"],
                                    ["Int", "Bool", "Float", "Hex"]),
                          row_names = ["Row $i" for i = 1:6],
                          row_name_column_title = "Row names")

    @test result == expected

    expected = """
\\begin{table}
  \\begin{tabular}{rcrrrr}
    \\hline\\hline
    \\textbf{Row} & \\textbf{Row names} & \\textbf{C1} & \\textbf{C2} & \\textbf{C3} & \\textbf{C4} \\\\
     &  & \\texttt{Int} & \\texttt{Bool} & \\texttt{Float} & \\texttt{Hex} \\\\\\hline
    1 & Row 1 & 1 & false & 1.0 & 1 \\\\
    2 & Row 2 & 2 & true & 2.0 & 2 \\\\
    3 & Row 3 & 3 & false & 3.0 & 3 \\\\
    4 & Row 4 & 4 & true & 4.0 & 4 \\\\
    5 & Row 5 & 5 & false & 5.0 & 5 \\\\
    6 & Row 6 & 6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data;
                          backend = Val(:latex),
                          header = (["C1",  "C2",   "C3",    "C4"],
                                    ["Int", "Bool", "Float", "Hex"]),
                          row_names = ["Row $i" for i = 1:6],
                          row_name_column_title = "Row names",
                          row_name_alignment = :c,
                          show_row_number = true)

    @test result == expected
end


text:
import Base.show
import Base.print

const spceStr = " " #string(spce)

abstract GraphicsArea
abstract Canvas <: GraphicsArea

function print(io::IO, c::GraphicsArea)
  for row in 1:nrows(c)
    printRow(io, c, row)
    print(io, "\n")
  end
end

function show(io::IO, c::GraphicsArea)
  b = borderDashed
  borderLength = ncols(c)
  drawBorderTop(io, "", borderLength, :solid)
  print(io, "\n")
  for row in 1:nrows(c)
    print(io, b[:l])
    printRow(io, c, row)
    print(io, b[:r], "\n")
  end
  drawBorderBottom(io, "", borderLength, :solid)
  print(io, "\n")
end

nrows(c::Canvas) = size(c.grid,2)
ncols(c::Canvas) = size(c.grid,1)

function printRow(io::IO, c::Canvas, row::Int)
  nunrows = nrows(c)
  0 < row <= nunrows || throw(ArgumentError("Argument row out of bounds: $row"))
  y = row
  for x in 1:ncols(c)
    printColor(c.colors[x,y], io, c.grid[x,y])
  end
end

ext: jl


text:
# this file depends on jacobi.jl from JuliaOptics/OpticsPolynomials that
# should be found in the same directory

export legendre, legendre_series, legendre_sum

"""
    legendre(n, x)

Compute the Legendre polynomial of order n at point x.

The Legendre polynomials are a special case of the Jacobi polynomials with
α, β = 0.

See also: [`legendre_series`](@ref), [`legendre_sum`](@ref)
"""
function legendre(n, x)
    return jacobi(n, 0, 0, x)
end

"""
    legendre_series(ns, α, β, x)

Compute a series of Legendre polynomials of orders n.
Returns an array with shape (size(x)..., length(ns)).
That is, the _final_ dimension contains the modes and the first dimension(s)
are spatial.

See also: [`legendre`](@ref), [`legendre_sum`](@ref)
"""
function legendre_series(ns, x)
    return jacobi_series(ns, 0., 0., x)
end
"""
    legendre_sum(ns, weights, x)

Compute a sum of Legendre polynomial of order n weighted by weights.

See also: [`legendre`](@ref), [`legendre_series`](@ref)
"""
function legendre_sum(ns, weights, x)
    return jacobi_sum(ns, weights, 0., 0., x)
end


text:
#=
Idmrg environments are only to be used internally.
They have to be updated manually, without any kind of checks
=#

struct IDMRGEnv{H,V<:AbstractTensorMap}
    opp :: H
    lw :: PeriodicArray{V,2}
    rw :: PeriodicArray{V,2}
end

function IDMRGEnv(state::Union{MPSMultiline,InfiniteMPS},env)
    state === env.dependency || recalculate!(env,state);
    IDMRGEnv(env.opp,deepcopy(env.lw),deepcopy(env.rw));
end

leftenv(envs::IDMRGEnv,pos::Int) = envs.lw[pos,:];
leftenv(envs::IDMRGEnv,row::Int,col::Int) = envs.lw[row,col];
leftenv(envs::IDMRGEnv,pos::Int,state::InfiniteMPS) = envs.lw[pos,:];
function setleftenv!(envs::IDMRGEnv,pos,lw)
    envs.lw[pos,:] = lw[:]
end
function setleftenv!(envs::IDMRGEnv,row,col,val)
    envs.lw[row,col] = val;
end

rightenv(envs::IDMRGEnv,row::Int,col::Int) = envs.rw[row,col];
rightenv(envs::IDMRGEnv,pos::Int) = envs.rw[pos,:];
rightenv(envs::IDMRGEnv,pos::Int,state::InfiniteMPS) = envs.rw[pos,:];
function setrightenv!(envs::IDMRGEnv,pos,rw)
    envs.rw[pos,:] = rw[:]
end
function setrightenv!(envs::IDMRGEnv,row,col,val)
    envs.rw[row,col] = val;
end


text:
function solve()
    n = parse(Int, readline())
    sum([x for x in 1:n if x%3≠0 && x%5≠0])
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# This file was generated, do not modify it. # hide
KNNC = @load KNNClassifier
MNC = @load MultinomialClassifier pkg=MLJLinearModels;

KnnPipe = Standardizer |> KNNC
MnPipe = Standardizer |> MNC

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
mutable struct AbstractParametersVSPARK{IT,DT,TT,D,S,R,P,FT,PT,UT,GT,ϕT} <: Parameters{DT,TT}
    f_f::FT
    f_p::PT
    f_u::UT
    f_g::GT
    f_ϕ::ϕT

    Δt::TT

    t_q::CoefficientsARK{TT}
    t_p::CoefficientsARK{TT}
    t_q̃::CoefficientsPRK{TT}
    t_p̃::CoefficientsPRK{TT}
    t_λ::CoefficientsMRK{TT}
    t_ω::Matrix{TT}
    t_δ::Matrix{TT}
    d_v::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function AbstractParametersVSPARK{IT,DT,TT,D,S,R,P,FT,PT,UT,GT,ϕT}(f_f, f_p, f_u, f_g, f_ϕ, Δt, t_q, t_p, t_q̃, t_p̃, t_λ, t_ω, t_δ, d_v) where {IT,DT,TT,D,S,R,P,FT,PT,UT,GT,ϕT}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new(f_f, f_p, f_u, f_g, f_ϕ, Δt,
            t_q, t_p, t_q̃, t_p̃, t_λ, t_ω, t_δ, d_v,
            zero(TT), q, p, λ)
    end
end

function update_params!(params::AbstractParametersVSPARK, cache::IntegratorCacheSPARK)
    # set time for nonlinear solver and copy previous solution
    params.t  = cache.t
    params.q .= cache.q
    params.p .= cache.p
    params.λ .= cache.λ
end

equation(int::AbstractIntegratorVSPARK) = int.equation
timestep(int::AbstractIntegratorVSPARK) = int.params.Δt
tableau(int::AbstractIntegratorVSPARK) = int.tableau
nstages(int::AbstractIntegratorVSPARK) = int.tableau.s
pstages(int::AbstractIntegratorVSPARK) = int.tableau.r

function create_integrator_cache(int::AbstractIntegratorVSPARK{DT,TT}) where {DT,TT}
    IntegratorCacheSPARK{DT, TT, ndims(int), nstages(int), pstages(int)}()
end

function initialize!(int::AbstractIntegratorVSPARK, cache::IntegratorCacheSPARK)
    cache.t̅ = cache.t - timestep(int)

    equation(int).v(cache.t, cache.q, cache.p, cache.v)
    equation(int).f(cache.t, cache.q, cache.p, cache.f)

    initialize!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f,
                            cache.t̅, cache.q̅, cache.p̅, cache.v̅, cache.f̅)
end

function initial_guess!(int::AbstractIntegratorVSPARK, cache::IntegratorCacheSPARK)
    for i in 1:nstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in 1:ndims(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+3] = cache.ṽ[k]
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = 0
        end
    end

    if int.params.t_λ.c[1] == 0
        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(k-1)+3] = cache.λ[k]
        end
    end

    if isdefined(tableau(int), :d)
        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*ndims(int)*pstages(int)+k] = 0
        end
    end
end


text:
@doc doc"""
Representation of a specific conditional density decomposition. A
`Vine` basically stores an `Array` of `Tree` instances, where each
`Tree` represents all univariate conditional information sets of a
given variable. For reasons of memory efficiency, individual trees are
not stored as `Array{Tree, 1}`, but each tree is stored as a single
column in a `Array{Int, 2}` matrix as what I call parent notation.
That is, entry [ii, jj] denotes the parent of variable ii in the
conditional information set tree of variable jj.

For future extension a `Vine` instance also comprises field `names`.

- trees::Array{Int, 2}

- names::Array{T, 1}

Matrix `trees` must be square.
"""->
type Vine
    trees::Array{Int, 2}

    function Vine(tr::Array{Int, 2})
        ## test symmetry
        if !(size(tr,1) == size(tr,2))
            error("tree matrix columns and rows must be of same size")
        end
        return new(tr)
    end
end

##################
## constructors ##
##################

function Vine(trs::Array{AbstractCTree, 1})
    nTrees = length(trs)
    parNot = Array(Int, nTrees, nTrees)
    for ii=1:nTrees
        currTree = convert(CTreeParRef, trs[ii])
        parNot[:, ii] = currTree.tree
    end
    return Vine(parNot)
end

function Vine(trs::Array{CTreePaths, 1})
    nTrees = length(trs)
    parNot = Array(Int, nTrees, nTrees)
    for ii=1:nTrees
        currTree = convert(CTreeParRef, trs[ii], nTrees)
        parNot[:, ii] = currTree.tree
    end
    return Vine(parNot)
end

function Vine(trs::Array{CTreeParRef, 1})
    nTrees = length(trs)
    parNot = Array(Int, nTrees, nTrees)
    for ii=1:nTrees
        currTree = convert(CTreeParRef, trs[ii])
        parNot[:, ii] = currTree.tree
    end
    return Vine(parNot)
end

## dg function
##------------

function dg(vn::Vine)
    return vn.trees
end

#############
## display ##
#############

import Base.Multimedia.display
@doc doc"""
Transform columns of parent notation matrix to tree instances and
display individual conditioning trees.
"""->
function display(vn::Vine)
    ## transform to Tree array

    nVars = size(dg(vn), 1)

    vineTrees = Array(CTreePaths, nVars)
    for ii=1:nVars
        trPar = CTreeParRef(vn.trees[:, ii])
        vineTrees[ii] = convert(CTreePaths, trPar)
    end
    println("VINE of $nVars variables:")
    println("")
    display(dg(vn))
    println("In TREE notation:")
    display(vineTrees)
end

##############
## equality ##
##############

function ==(vn1::Vine, vn2::Vine)
    return dg(vn1) == dg(vn2)
end

################
## conversion ##
################

import Base.convert
function convert(::Type{CTreePaths}, vn::Vine)
    nVars = size(vn.trees, 1)
    trees = Array(CTreePaths, nVars)
    for ii=1:nVars
        trees[ii] = convert(CTreePaths, vn.trees[:, ii])
    end
    return trees
end

function convert(::Type{CTreeParRef}, vn::Vine)
    nVars = size(vn.trees, 1)
    trees = Array(CTreeParRef, nVars)
    for ii=1:nVars
        trees[ii] = convert(CTreeParRef, vn.trees[:, ii])
    end
    return trees
end

function convert(::Type{Vine}, trArr::Array{CTreePaths, 1})
    nVars = length(trArr)
    vnTrees = Array(Int, nVars, nVars)
    for ii=1:nVars
        trPar = convert(CTreeParRef, trArr[ii])
        vnTrees[:, ii] = trPar.tree
    end
    return Vine(vnTrees)
end

function convert(::Type{Vine}, trArr::Array{CTreeParRef, 1})
    nVars = length(trArr)
    vnTrees = Array(Int, nVars, nVars)
    for ii=1:nVars
        trPar = convert(CTreeParRef, trArr[ii])
        vnTrees[:, ii] = trPar.tree
    end
    return Vine(vnTrees)
end

##############
## analysis ##
##############

@doc doc"""
Array{Array{Int, 1}, 1} of all possible simulation sequences given
that the variables of `condSet` are already known.

Applies condSetChk to each conditioning set and each tree. If
condSetChk returns true for a given tree, then the tree's root is
added to the conditioning set.

"""->
function getSimSequences(vn::Vine, condSet::Array{Int, 1})
    ## find possible simulation variable sequences

    ## transform vine to array of c-trees in parent reference notation
    nVars = size(dg(vn), 1)
    trees = CTreePaths[CTreeParRef(dg(vn)[:, ii]) for ii=1:nVars]

    nVars = size(vn.trees, 1)

    ## preallocation
    condSetList = Array(Array{Int, 1}, 1)
    condSetList[1] = condSet

    ## partialSequences = Array(Array{Int, 1}, 0)
    simSequences = Array(Array{Int, 1}, 0)

    ## algorithm:
    ## - take last element of list
    ## - get missing variables
    ## - run through missing variables
    ## - add matching variables to condList
    ## - finished?

    while !isempty(condSetList)
        currentCondSet = pop!(condSetList)
        println("Current condition set: $currentCondSet")

        ## get remaining variables to be simulated
        toSim = setdiff([1:nVars], currentCondSet)
        nRemaining = length(toSim)

        for ii=1:nRemaining
            missingVar = toSim[nRemaining+1-ii]
            if condSetChk(currentCondSet, trees[missingVar])
                ## add to list
                newCondSet = [currentCondSet, missingVar]
                if length(newCondSet) == nVars
                    push!(simSequences, [newCondSet])
                    println("New simulation sequence found!")
                else
                    println("New condition set: $newCondSet")
                    push!(condSetList, [newCondSet])
                    ## display(condSetList)
                end
            end
        end
    end
    return simSequences
end

@doc doc"""
`Array{Array{Int, 1}, 1}` of all possible simulation sequences for a
given vine.

WARNING: no unit tests.
"""->
function getSimSequences(vn::Vine)
    ## get all variable sequences for simulation
    nVars = size(vn.trees, 1)
    simSequences = Array(Array{Int, 1}, 0)
    for ii=1:nVars
        newSequ = getSimSequences(vn, [ii])
        simSequences = [simSequences, newSequ]
    end
    return simSequences
end

@doc doc"""
Evaluating the result of `getSimSequences`: for each variable find its
average position in the sequences of simulation. This could be an
indicator of how central a given variable is, and where it should be
plotted in a visualization of vines. Also, variables that can be
observed earlier should probably get an earlier average position in
the simulation sequences.

WARNING: no unit tests.
"""->
function avgSimSequPosition(simSequ::IntArrays)
    ## whats the average position of a variable during simulation?
    nVars = length(simSequ[1])
    nSequ = length(simSequ)
    avgPos = Array(Float64, nVars)
    println("")
    for ii=1:nVars
        val = 0
        for jj=1:nSequ
            val = val + find(ii .== simSequ[jj])[1]
        end
        avgVal = val/nSequ
        avgPos[ii] = avgVal
        println("$ii: $avgVal")
    end
    println("")
    return avgPos
end

@doc doc"""
`Array{Int, 2}` square matrix with level of layers where each pairs of
variables are connected. A number of 1 denotes pairs that are
connected unconditionally, while a number of 3 denotes pairs that are
connected with a conditioning set of size 2.

WARNING: no unit tests.
"""->
function linkLayers(vn::Vine)
    ## at which layer are links between variables modelled?
    nVars = size(vn.trees, 1)
    links = zeros(Int, nVars, nVars)
    for ii=1:(nVars-1)
        root = ii
        for jj=(ii+1):nVars
            # follow jj in column ii to root
            notAtRoot = true
            steps = 0
            currNode = jj
            while notAtRoot
                steps = steps + 1
                currNode = vn.trees[currNode, ii]
                if currNode == root
                    notAtRoot = false
                end
            end
            links[jj, ii] = steps
            links[ii, jj] = steps
        end
    end
    return links
end

#######################
## conditioning sets ##
#######################

@doc doc"""
Get the conditioning set for variables `var1` and `var2`. More
specific, if var1, var2|k, l, m, then return conditioning set {k,l,m}
such that the order matches the requirement to calculate var1|k,l,m.
"""->
function getCondSet(vn::Vine, var1::Int, var2::Int)
    parNot = vn.trees[:, var1]
    trPar = CTreeParRef(parNot)
    pathToVar1 = getPathToRoot(trPar, var2)
    return pathToVar1
end

@doc doc"""
Get the conditioning sets for all pairs of variables. Conditioning
sets are given in the order of linear triangular indexing. Each
conditioning set is sorted such that it represents the conditioning
set for the variable with smaller index.
"""->
function getAllCopCondSets(vn::Vine)
    nVars = size(vn.trees, 1)
    nSets = sum((nVars-1):-1:1)

    condSets = Array(Array{Int, 1}, nSets)
    ind = 1
    for ii=1:(nVars-1)
        root = ii
        for jj=(ii+1):nVars
            currCondSet = getPathToRoot(CTreeParRef(vn.trees[:, ii]), jj)
            condSets[ind] = currCondSet
            ind += 1
        end
    end
    return condSets
end

@doc doc"""
Get both conditioning sets for all pairs of variables. That is, the
first output returns the conditioning sets in correct order for the
smaller variable, while the second output has them ordered for the
larger variable.
"""->
function getAllVarCondSets(vn::Vine)
    nVars = size(vn.trees, 1)
    nSets = sum((nVars-1):-1:1)

    condSets1 = Array(Array{Int, 1}, nSets)
    condSets2 = Array(Array{Int, 1}, nSets)
    ind = 1
    for ii=1:(nVars-1)
        root = ii
        for jj=(ii+1):nVars
            path = getPathToRoot(CTreeParRef(vn.trees[:, ii]), jj)
            condSets1[ind] = path
            path = getPathToRoot(CTreeParRef(vn.trees[:, jj]), ii)
            condSets2[ind] = path
            ind += 1
        end
    end
    return (condSets1, condSets2)
end


text:
"""
Generate a population of size `n` distributed spatially according to
a bivariate distribution, `d`
"""
function create_pop_db_bivariate(n::Int64, d::ContinuousMultivariateDistribution)
  coordinates = rand(d, n)'
  DataFrame(ind_id = 1:n, x = coordinates[:,1], y = coordinates[:,2])
end

"""
Generate a population of size `n` distributed spatially according to
a univariate distributions, `d1` and `d2`
"""
function create_pop_db_univariates(n::Int64,
                                   d1::ContinuousUnivariateDistribution,
                                   d2::ContinuousUnivariateDistribution)
  DataFrame(ind_id = 1:n, x = rand(d1, n), y = rand(d2, n))
end

"""
Construct a matrix containing the euclidean distance between every
pair of individuals contained in `pop_db`
"""
function create_dist_mtx(pop_db::DataFrame)
  distance_mat = fill(0.0, (pop_db[end, 1], pop_db[end, 1]))
  for i = 1:(size(distance_mat)[1]-1)
    for j = (i+1):(size(distance_mat)[1])
      distance_mat[i,j] = euclidean([pop_db[i, 2],pop_db[i, 3]], [pop_db[j, 2],pop_db[j, 3]])
    end
  end
  distance_mat' + distance_mat
end


text:
using WindMixing

DATA_NAMES = [
  "NDE_18sim_windcooling_windheating_18simBFGST0.8grad_divide1f5_gradient_smallNN_leakyrelu_rate_2e-4_T0.8"
]

EXTRACTED_DATA_DIR = joinpath(pwd(), "extracted_training_output")
OUTPUT_DIR = joinpath("C:\\Users\\xinle\\Documents\\OceanParameterizations.jl")

test_files = [
  # "wind_-5e-4_diurnal_5e-8"
  # "wind_-5e-4_diurnal_3e-8"
  # "wind_-5e-4_diurnal_1e-8"

  # "wind_-3.5e-4_diurnal_5e-8"
  # "wind_-3.5e-4_diurnal_3e-8"
  # "wind_-3.5e-4_diurnal_1e-8"

  # "wind_-2e-4_diurnal_5e-8"
  # "wind_-2e-4_diurnal_3e-8"
  # "wind_-2e-4_diurnal_1e-8"

  # "wind_-2e-4_diurnal_2e-8"
  # "wind_-2e-4_diurnal_3.5e-8"
  # "wind_-3.5e-4_diurnal_2e-8"
  # "wind_-3.5e-4_diurnal_3.5e-8"
  # "wind_-5e-4_diurnal_2e-8"
  # "wind_-5e-4_diurnal_3.5e-8"

  "wind_-5e-4_cooling_3e-8_new",
  # "wind_-5e-4_cooling_1e-8_new",
  # "wind_-2e-4_cooling_3e-8_new",
  # "wind_-2e-4_cooling_1e-8_new",
  # "wind_-5e-4_heating_-3e-8_new",
  # "wind_-2e-4_heating_-1e-8_new",
  # "wind_-2e-4_heating_-3e-8_new",
  # "wind_-5e-4_heating_-1e-8_new",

  # "wind_-3.5e-4_cooling_2e-8_new",
  # "wind_-3.5e-4_heating_-2e-8_new",

  # "wind_-5e-4_cooling_2e-8_new",
  # "wind_-3.5e-4_cooling_3e-8_new",
  # "wind_-3.5e-4_cooling_1e-8_new",
  # "wind_-2e-4_cooling_2e-8_new",
  # "wind_-3.5e-4_heating_-3e-8_new",
  # "wind_-3.5e-4_heating_-1e-8_new",
  # "wind_-2e-4_heating_-2e-8_new",
  # "wind_-5e-4_heating_-2e-8_new",
  ]

for DATA_NAME in DATA_NAMES
    animate_training_results(test_files, DATA_NAME,
                         EXTRACTED_DATA_DIR=EXTRACTED_DATA_DIR, OUTPUT_DIR=OUTPUT_DIR, convective_adjustment=true)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using BenchmarkTools
using Base.Sort

# %%
function f!(a, b)
    for x in b
        push!(a, x)
        sort!(a)
    end
    a
end

function g!(a, b)
    append!(a, b)
    sort!(a)
end

function h!(a, b)
    for x in b
        i = Base.Sort.searchsortedfirst(a, x)
        insert!(a, i, x)
    end
    a
end

# %%
a, b = sort(rand(1000)), rand(10)

# %%
f!(copy(a), b) == g!(copy(a), b) == h!(copy(a), b)

# %%
@btime f!(A, $b) setup=(A = copy(a));

# %%
@btime g!(A, $b) setup=(A = copy(a));

# %%
@btime h!(A, $b) setup=(A = copy(a));

# %%
@btime f!(A, $b) setup=(A = copy(a); sizehint!(A, length(a)+length(b)));

# %%
@btime g!(A, $b) setup=(A = copy(a); sizehint!(A, length(a)+length(b)));

# %%
@btime h!(A, $b) setup=(A = copy(a); sizehint!(A, length(a)+length(b)));

# %%
A = copy(a)
@btime begin
    for i in 1:10
        push!(A, b[i])
        sort!(A)
    end
    A
end

# %%
@btime begin
    for i in 1:10
        push!(B, b[i])
        sort!(B)
    end
    B
end setup=(B = copy(a); sizehint!(B, length(a)+length(b)));

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@sumProductRule(:node_type     => Nonlinear{Sampling},
                :outbound_type => Message{SampleList},
                :inbound_types => (Nothing, Message),
                :name          => SPNonlinearSOutNM)

@sumProductRule(:node_type     => Nonlinear{Sampling},
                :outbound_type => Message{Function},
                :inbound_types => (Message, Nothing),
                :name          => SPNonlinearSIn1MN)

mutable struct SPNonlinearSOutNGX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSOutNGX}) = Message{SampleList}
function isApplicable(::Type{SPNonlinearSOutNGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false

    for input_type in input_types[2:end]
        matches(input_type, Message{Gaussian}) || return false
    end

    return true
end

mutable struct SPNonlinearSInGX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSInGX}) = Message{GaussianWeightedMeanPrecision}
function isApplicable(::Type{SPNonlinearSInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] != Nothing) || return false

    nothing_inputs = 0
    factorfunction_input = false
    gaussian_inputs = 0
    if matches(input_types[1], Message{FactorFunction})
        factorfunction_input = true
    end
    for input_type in input_types[2:end]
        if input_type == Nothing
            nothing_inputs += 1
        elseif matches(input_type, Message{Gaussian})
            gaussian_inputs += 1
        end
    end

    return (nothing_inputs == 1) && (gaussian_inputs == total_inputs-2) && factorfunction_input
end

mutable struct SPNonlinearSOutNFactorX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSOutNFactorX}) = Message{SampleList}
function isApplicable(::Type{SPNonlinearSOutNFactorX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false

    factorNode_input = false
    gaussian_inputs = 0
    for input_type in input_types[2:end]
        if matches(input_type, Message{SampleList})
            return false
        elseif matches(input_type, Message{FactorNode})
            factorNode_input += 1
            if matches(input_type, Message{Gaussian})
                gaussian_inputs += 1
            end
        else
            return false
        end
    end

    return (gaussian_inputs < total_inputs-1) && (factorNode_input == total_inputs-1)
end

mutable struct SPNonlinearSInFactorX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSInFactorX}) = Message{Function}
function isApplicable(::Type{SPNonlinearSInFactorX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] != Nothing) || return false

    nothing_inputs = 0
    factorfunction_input = false
    factorNode_input = 0
    gaussian_inputs = 0
    if matches(input_types[1], Message{FactorFunction})
        factorfunction_input = true
    end
    for input_type in input_types[2:end]
        if input_type == Nothing
            nothing_inputs += 1
        elseif matches(input_type, Message{PointMass})
            return false
        elseif matches(input_type, Message{SampleList})
            return false
        elseif matches(input_type, Message{FactorNode})
            factorNode_input += 1
            if matches(input_type, Message{Gaussian})
                gaussian_inputs += 1
            end
        else
            return false
        end
    end

    return (nothing_inputs == 1) && (gaussian_inputs < total_inputs-2) && (factorNode_input == total_inputs-2) && factorfunction_input
end

mutable struct MNonlinearSInGX <: MarginalRule{Nonlinear{Sampling}} end
function isApplicable(::Type{MNonlinearSInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false # Indicates marginalization over outbound variable

    for input_type in input_types[2:end]
        matches(input_type, Message{Gaussian}) || return false
    end

    return true
end


text:
using Documenter
using AxisIndices
using LinearAlgebra
using Metadata
using Statistics

makedocs(;
    modules=[AxisIndices],
    format=Documenter.HTML(),
    pages=[
        "index.md",
        "references.md",
    ],
    repo="https://github.com/Tokazama/AxisIndices.jl/blob/{commit}{path}#L{line}",
    sitename="AxisIndices.jl",
    authors="Zachary P. Christensen",
)

deploydocs(
    repo = "github.com/Tokazama/AxisIndices.jl.git",
)


text:
using Categores, Test

function setup(::Type{K}=Float64) where {K}
    A,B,C,D = Category.([:A,:B,:C,:D],Symbol)

    W,X,Y,Z = Ob.([A,B,C,D],[:W,:X,:Y,:Z])

    f = Hom(:f,W,X)
    g = Hom(:g,X,Y)
    h = Hom(:h,Y,Z)

    tuv = Tensor(U⊗V,[1 2 3;4 5 6])
    tvv = Tensor(V⊗V,[1 2 3;4 5 6;7 8 9])
    tvw = Tensor(V⊗W,[1 2 3 4;5 6 7 8;9 10 11 12])
    auvw = Array{Float64,3}(undef,2,3,4)
    auvw[:,:,1] = [1 2 3;4 5 6]
    auvw[:,:,2] = [7 8 9;10 11 12]
    auvw[:,:,3] = [13 14 15;16 17 18]
    auvw[:,:,4] = [19 20 21;22 23 24]
    tuvw = Tensor(U⊗V⊗W,auvw)

    uα = u⊗α
    uvw = u⊗v⊗w
    αβγ = α⊗β⊗γ

    U,V,W,u,v,w,α,β,γ,tuv,tvv,tvw,tuvw,uα,uvw,αβγ
end

U,V,W,u,v,w,α,β,γ,tuv,tvv,tvw,tuvw,uα,uvw,αβγ = setup(Float64)

@testset "Dimension errors" begin
    @test_throws DimensionMismatch Vector(V,[1,2])
    @test_throws DimensionMismatch Vector(V,[1,2,3,4])
end

@testset "Domain errors" begin
    @test_throws DomainError α(v)
    @test_throws DomainError (α⊗β)(-,u)
    @test_throws DomainError (α⊗β)(v,-)
    @test_throws DomainError α⊗β+β⊗α
    @test_throws DomainError α⊗β-β⊗α
end

@testset "Inclusion" begin
    @test (v ∈ U) === false
    @test (v ∈ V) === true
    @test (v ∈ W) === false
    @test (v ∈ U^*) === false
    @test (v ∈ V^*) === false
    @test (v ∈ W^*) === false
    @test (w ∈ U) === false
    @test (w ∈ V) === false
    @test (w ∈ W) === true
    @test (w ∈ U^*) === false
    @test (w ∈ V^*) === false
    @test (w ∈ W^*) === false
    @test (α ∈ U) === false
    @test (α ∈ V) === false
    @test (α ∈ W) === false
    @test (α ∈ U^*) === true
    @test (α ∈ V^*) === false
    @test (α ∈ W^*) === false
    @test (β ∈ U) === false
    @test (β ∈ V) === false
    @test (β ∈ W) === false
    @test (β ∈ U^*) === false
    @test (β ∈ V^*) === true
    @test (β ∈ W^*) === false
    @test (α⊗β ∈ (U⊗V)^*) === true
end

@testset "Product spaces" begin
    @test V×V === ProductSpace(V,V)
    @test V×W === ProductSpace(V,W)
    @test V×(V^*) === ProductSpace(V,dual(V))
    @test V×(W^*) === ProductSpace(V,dual(W))
    @test V⊗V === TensorSpace(V,V)
    @test V⊗W === TensorSpace(V,W)
    @test V⊗(V^*) === TensorSpace(V,dual(V))
    @test V⊗(W^*) === TensorSpace(V,dual(W))
end

@testset "Addition, subtraction and negation" begin
    @test u+u ≈ 2*u
    @test u-u ≈ 0*u
    @test -u ≈ -1*u
end

@testset "Tensor products" begin
    @test (domain(uα) === U⊗U) === false
    @test (domain(uα) === (U^*)⊗U) === true
    @test (domain(uα) === U⊗U^*) === false
    @test (domain(uα) === (U⊗U)^*) === false
    @test (domain(uvw) === U⊗V⊗W) === false
    @test (domain(uvw) === (U^*)⊗V⊗W) === false
    @test (domain(uvw) === U⊗(V^*)⊗W) === false
    @test (domain(uvw) === U⊗V⊗W^*) === false
    @test (domain(uvw) === ((U⊗V)^*)⊗W) === false
    @test (domain(uvw) === U⊗(V⊗W)^*) === false
    @test (domain(uvw) === (U⊗V⊗W)^*) === true
    @test (domain(αβγ) === U⊗V⊗W) === true
end

@testset "Evaluation" begin
    @test α(u) ≈ 5.0
    @test u(α) ≈ 5.0
    @test (α⊗β)(u,v) ≈ α(u)*β(v)
    @test (α⊗β)(u⊗v) ≈ α(u)*β(v)
    @test (α⊗β)(-,v)(u) ≈ α(u)*β(v)
    @test (α⊗β)(u,-)(v) ≈ α(u)*β(v)
    @test (u⊗v)(α,β) ≈ u(α)*v(β)
    @test (u⊗v)(α⊗β) ≈ u(α)*v(β)
    @test (u⊗v)(-,β)(α) ≈ u(α)*v(β)
    @test (u⊗v)(α,-)(β) ≈ u(α)*v(β)
    @test (α⊗β⊗γ)(u,v,w) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(u⊗v⊗w) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(-,v,w)(u) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(u,-,w)(v) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(u,v,-)(w) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(-,-,w)(u,v) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(-,v,-)(u,w) ≈ α(u)*v(β)*γ(w)
    @test (α⊗β⊗γ)(u,-,-)(v,w) ≈ α(u)*v(β)*γ(w)
    @test tuv(-,v)(u) ≈ tuv(u,v)
    @test tuv(u,-)(v) ≈ tuv(u,v)
    @test tuvw(-,v,w)(u) ≈ tuvw(u,v,w)
    @test tuvw(u,-,w)(v) ≈ tuvw(u,v,w)
    @test tuvw(u,v,-)(w) ≈ tuvw(u,v,w)
    @test tuvw(-,-,w)(u,v) ≈ tuvw(u,v,w)
    @test tuvw(-,v,-)(u,w) ≈ tuvw(u,v,w)
    @test tuvw(u,-,-)(v,w) ≈ tuvw(u,v,w)
    @test (α⊗β⊗γ⊗tuvw)(u,v,w,u,v,w) ≈ α(u)*v(β)*γ(w)*tuvw(u,v,w)
end


text:
using .ArchGDAL

const AG = ArchGDAL

const GDAL_LON_INDEX = ForwardIndex()
const GDAL_LAT_INDEX = ReverseIndex()
const GDAL_BAND_INDEX = ForwardIndex()
const GDAL_LON_ARRAY = ForwardArray()
const GDAL_LAT_ARRAY = ReverseArray()
const GDAL_BAND_ARRAY = ForwardArray()
const GDAL_RELATION = ForwardRelation()
const GDAL_LON_LOCUS = Start()
const GDAL_LAT_LOCUS = Start()

export GDALarray, GDALstack, GDALarrayMetadata, GDALdimMetadata

# Metadata ########################################################################

"""
    GDALdimMetadata(val::Union{Dict,NamedTuple})
    GDALdimMetadata(pairs::Pair...) => GDALdimMetadata{Dict}
    GDALdimMetadata(; kw...) => GDALdimMetadata{NamedTuple}

`Metadata` wrapper for `GDALarray` dimensions.
"""
struct GDALdimMetadata{T} <: AbstractDimMetadata{T}
    val::T
end

"""
    GDALarrayMetadata(val::Union{Dict,NamedTuple})
    GDALarrayMetadata(pairs::Pair...) => GDALarrayMetadata{Dict}
    GDALarrayMetadata(; kw...) => GDALarrayMetadata{NamedTuple}

`Metadata` wrapper for `GDALarray`.
"""
struct GDALarrayMetadata{T} <: AbstractArrayMetadata{T}
    val::T
end

# Array ########################################################################

"""
    GDALarray(filename;
              crs=nothing,
              mappedcrs=nothing,
              name="",
              dims=nothing,
              refdims=(),
              metadata=GDALdimMetadata(),
              missingval=nothing)

Load a file lazily using gdal. `GDALarray` will be converted to [`GeoArray`](@ref)
after indexing or other manipulations. `GeoArray(GDALarray(filename))` will do this
immediately.

`GDALarray`s are always 3 dimensional, and have [`Lat`](@ref), [`Lon`](@ref) and
[`Band`](@ref) dimensions.

## Arguments

- `filename`: `String` pointing to a grd file. Extension is optional.

## Keyword arguments

- `usercrs`: CRS format like `EPSG(4326)` used in `Selectors` like `Between` and `At`, and
  for plotting. Can be any CRS `GeoFormat` from GeoFormatTypes.jl, like `WellKnownText`.
- `name`: `Symbol` name for the array.
- `dims`: `Tuple` of `Dimension`s for the array. Detected automatically, but can be passed in.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from.
- `missingval`: Value reprsenting missing values. Detected automatically when possible, but
  can be passed it.
- `metadata`: [`Metadata`](@ref) object for the array. Detected automatically as
  [`GDALarrayMetadata`](@ref), but can be passed in.

# Example

```julia
A = GDALarray("folder/file.tif"; usercrs=EPSG(4326))
# Select Australia using lat/lon coords, whatever the crs is underneath.
A[Lat(Between(-10, -43), Lon(Between(113, 153)))
```
"""
struct GDALarray{T,N,F,D<:Tuple,R<:Tuple,Na<:Symbol,Me,Mi,S
                } <: DiskGeoArray{T,N,D,LazyArray{T,N}}
    filename::F
    dims::D
    refdims::R
    name::Na
    metadata::Me
    missingval::Mi
    size::S
end
GDALarray(filename::AbstractString; kwargs...) = begin
    isfile(filename) || error("file not found: $filename")
    gdalread(filename) do raster
        GDALarray(raster, filename; kwargs...)
    end
end
GDALarray(raster::AG.RasterDataset, filename, key=nothing;
          crs=nothing,
          mappedcrs=nothing,
          dims=dims(raster, crs, mappedcrs),
          refdims=(),
          name=Symbol(""),
          metadata=metadata(raster),
          missingval=missingval(raster)) = begin
    sze = size(raster)
    T = eltype(raster)
    N = length(sze)
    name = Symbol(name)
    GDALarray{T,N,typeof.((filename,dims,refdims,name,metadata,missingval,sze))...
             }(filename, dims, refdims, name, metadata, missingval, sze)
end

# AbstractGeoArray methods

"""
    Base.write(filename::AbstractString, ::Type{GDALarray}, A::AbstractGeoArray;
               driver="GTiff", compress="DEFLATE", tiled=true)

Write a [`GDALarray`](@ref) to file, `.tiff` by default, but other GDAL drivers also work.

GDAL flags `driver`, `compress` and `tiled` can be passed in as keyword arguments.

Returns `filename`.
"""
Base.write(filename::AbstractString, ::Type{<:GDALarray}, A::AbstractGeoArray{T,2}; kwargs...) where T = begin
    all(hasdim(A, (Lon, Lat))) || error("Array must have Lat and Lon dims")

    correctedA = permutedims(A, (Lon(), Lat())) |>
        a -> reorder(a, (Lon(GDAL_LON_INDEX), Lat(GDAL_LAT_INDEX))) |>
        a -> reorder(a, GDAL_RELATION)
    checkarrayorder(correctedA, (GDAL_LON_ARRAY, GDAL_LAT_ARRAY))

    nbands = 1
    indices = 1
    gdalwrite(filename, correctedA, nbands, indices; kwargs...)
end
Base.write(filename::AbstractString, ::Type{<:GDALarray}, A::AbstractGeoArray{T,3}, kwargs...) where T = begin
    all(hasdim(A, (Lon, Lat))) || error("Array must have Lat and Lon dims")
    hasdim(A, Band()) || error("Must have a `Band` dimension to write a 3-dimensional array")

    correctedA = permutedims(A, (Lon(), Lat(), Band())) |>
        a -> reorder(a, (Lon(GDAL_LON_INDEX), Lat(GDAL_LAT_INDEX), Band(GDAL_BAND_INDEX))) |>
        a -> reorder(a, GDAL_RELATION)
    checkarrayorder(correctedA, (GDAL_LON_ARRAY, GDAL_LAT_ARRAY, GDAL_BAND_ARRAY))

    nbands = size(correctedA, Band())
    indices = Cint[1:nbands...]
    gdalwrite(filename, correctedA, nbands, indices; kwargs...)
end

# AbstractGeoStack methods

"""
    GDALstack(filenames; keys, kwargs...)
    GDALstack(filenames...; keys, kwargs...)
    GDALstack(filenames::NamedTuple;
              window=(),
              metadata=NoMetadata(),
              childkwargs=(),
              refdims=())

Convenience method to create a DiskStack  of [`GDALarray`](@ref) from `filenames`.

Load a stack of files lazily from disk.

## Arguments

- `filenames`: A NamedTuple of stack keys and `String` filenames, or a `Tuple`,
  `Vector` or splatted arguments of `String` filenames.

## Keyword arguments

- `keys`: Used as stack keys when a `Tuple`, `Vector` or splat of filenames are passed in.
- `window`: A `Tuple` of `Dimension`/`Selector`/indices that will be applied to the
  contained arrays when they are accessed.
- `metadata`: a [`GDALstackMetadata`](@ref) object.
- `childkwargs`: A `NamedTuple` of keyword arguments to pass to the `childtype` constructor.
- `refdims`: `Tuple` of  position `Dimension` the array was sliced from.

## Example

Create a `GDALstack` from four files, that sets the child arrays `usercrs` value
when they are loaded.

```julia
files = (:temp="temp.tif", :pressure="pressure.tif", :relhum="relhum.tif")
stack = GDALstack(files; childkwargs=(usercrs=EPSG(4326),))
stack[:relhum][Lat(Contains(-37), Lon(Contains(144))
```
"""
GDALstack(args...; kwargs...) =
    DiskStack(args...; childtype=GDALarray, kwargs...)

withsource(f, ::Type{<:GDALarray}, filename::AbstractString, key...) =
    gdalread(f, filename)

# DimensionalData methods for ArchGDAL types ###############################

dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing) = begin
    gt = try
        AG.getgeotransform(raster) catch GDAL_EMPTY_TRANSFORM end
    lonsize, latsize = size(raster)

    nbands = AG.nraster(raster)
    band = Band(1:nbands, mode=Categorical(Ordered()))
    crs = crs isa Nothing ? GeoData.crs(raster) : crs

    lonlat_metadata = GDALdimMetadata()

    # Output Sampled index dims when the transformation is lat/lon alligned,
    # otherwise use Transformed index, with an affine map.
    if isalligned(gt)
        lonstep = gt[GDAL_WE_RES]
        lonmin = gt[GDAL_TOPLEFT_X]
        lonmax = gt[GDAL_TOPLEFT_X] + lonstep * (lonsize - 1)
        lonindex = LinRange(lonmin, lonmax, lonsize)

        latstep = gt[GDAL_NS_RES] # A negative number
        latmax = gt[GDAL_TOPLEFT_Y] + latstep
        latmin = gt[GDAL_TOPLEFT_Y] + latstep * latsize
        latindex = LinRange(latmax, latmin, latsize)

        # Spatial data defaults to area/inteval
        lonsampling, latsampling = if gdalmetadata(raster.ds, "AREA_OR_POINT") == "Point"
            Points(), Points()
        else
            # GeoTiff uses the "pixelCorner" convention
            Intervals(GDAL_LON_LOCUS), Intervals(GDAL_LAT_LOCUS)
        end

        lonmode = Projected(
            order=Ordered(GDAL_LON_INDEX, GDAL_LON_ARRAY, GDAL_RELATION),
            span=Regular(step(lonindex)),
            sampling=lonsampling,
            crs=crs,
            mappedcrs=mappedcrs,
        )
        latmode = Projected(
            order=Ordered(GDAL_LAT_INDEX, GDAL_LAT_ARRAY, GDAL_RELATION),
            sampling=latsampling,
            # Use the range step as is will be different to latstep due to float error
            span=Regular(step(latindex)),
            crs=crs,
            mappedcrs=mappedcrs,
        )

        lon = Lon(lonindex; mode=lonmode, metadata=lonlat_metadata)
        lat = Lat(latindex; mode=latmode, metadata=lonlat_metadata)

        DimensionalData._formatdims(map(Base.OneTo, (lonsize, latsize, nbands)), (lon, lat, band))
    else
        error("Rotated/transformed dimensions are not handled yet. Open a github issue for GeoData.jl if you need this.")
        # affinemap = geotransform2affine(geotransform)
        # x = X(affinemap; mode=TransformedIndex(dims=Lon()))
        # y = Y(affinemap; mode=TransformedIndex(dims=Lat()))

        # formatdims((lonsize, latsize, nbands), (x, y, band))
    end
end

missingval(raster::AG.RasterDataset, args...) = begin
    # We can only handle data where all bands have
    # the same missingval
    band = AG.getband(raster.ds, 1)
    missingval = AG.getnodatavalue(band)
    T = AG.pixeltype(band)
    try
        missingval = convert(T, missingval)
    catch
        @warn "No data value from GDAL $(missingval) is not convertible to data type $T. `missingval` is probably incorrect."
    end
    missingval
end

metadata(raster::AG.RasterDataset, args...) = begin
    band = AG.getband(raster.ds, 1)
    # color = AG.getname(AG.getcolorinterp(band))
    scale = AG.getscale(band)
    offset = AG.getoffset(band)
    # norvw = AG.noverview(band)
    path = first(AG.filelist(raster))
    units = AG.getunittype(band)
    upair = units == "" ? () : (:units=>units,)
    GDALarrayMetadata(Dict(:filepath=>path, :scale=>scale, :offset=>offset, upair...))
end

# metadata(raster::RasterDataset, key) = begin
#     regex = Regex("$key=(.*)")
#     i = findfirst(f -> occursin(regex, f), meta)
#     if i isa Nothing
#         nothing
#     else
#         match(regex, meta[i])[1]
#     end
# end

crs(raster::AG.RasterDataset, args...) =
    WellKnownText(GeoFormatTypes.CRS(), string(AG.getproj(raster.ds)))

# Utils ########################################################################

gdalmetadata(dataset::AG.Dataset, key) = begin
    meta = AG.metadata(dataset)
    regex = Regex("$key=(.*)")
    i = findfirst(f -> occursin(regex, f), meta)
    if i isa Nothing
        nothing
    else
        match(regex, meta[i])[1]
    end
end

gdalread(f, filename::AbstractString) =
    AG.readraster(filename) do raster
        f(raster)
    end

gdalwrite(filename, A, nbands, indices; driver="GTiff", compress="DEFLATE", tiled=true) = begin
    tiledstring = tiled isa Bool ? (tiled ? "YES" : "NO") : tiled
    options = driver == "GTiff" ? ["COMPRESS=$compress", "TILED=$tiledstring"] : String[]

    AG.create(filename;
        driver=AG.getdriver(driver),
        width=size(A, Lon()),
        height=size(A, Lat()),
        nbands=nbands,
        dtype=eltype(A),
        options=options,
    ) do dataset
        gdalsetproperties!(dataset, A)
        AG.write!(dataset, data(A), indices)
    end
    return filename
end

function gdalsetproperties!(dataset, A)
    # Convert the dimensions to `Projected` if they are `Converted`
    # This allows saving NetCDF to Tiff
    # Set the index loci to the start of the cell for the lat and lon dimensions.
    # NetCDF or other formats use the center of the interval, so they need conversion.
    lon = shiftindexloci(GDAL_LON_LOCUS, dims(A, Lon))
    lat = shiftindexloci(GDAL_LAT_LOCUS, dims(A, Lat))
    lon = convertmode(Projected, lon)
    lat = convertmode(Projected, lat)
    @assert indexorder(lat) == GDAL_LAT_INDEX
    @assert indexorder(lon) == GDAL_LON_INDEX
    # Get the geotransform from the updated lat/lon dims
    geotransform = dims2geotransform(lat, lon)
    # Convert projection to a string of well known text
    proj = convert(String, convert(WellKnownText, crs(lon)))

    # Write projection, geotransform and data to GDAL
    AG.setproj!(dataset, proj)
    AG.setgeotransform!(dataset, geotransform)

    # Set the nodata value. GDAL can't handle missing
    # We could choose a default, but we would need to do this for all
    # possible types.
    if missingval(A) !== missing
        bands = if hasdim(A, Band)
            index(A, Band)
        else
            1
        end
        for i in bands
            AG.setnodatavalue!(AG.getband(dataset, i), missingval(A))
        end
    end

    dataset
end

# Create a GeoArray from a memory-backed dataset
function GeoArray(dataset::AG.Dataset;
                  crs=nothing,
                  mappedcrs=nothing,
                  dims=dims(AG.RasterDataset(dataset), crs, mappedcrs),
                  refdims=(),
                  name=Symbol(""),
                  metadata=metadata(AG.RasterDataset(dataset)),
                  missingval=missingval(AG.RasterDataset(dataset)))
    GeoArray(AG.read(dataset), dims, refdims, name, metadata, missingval)
end

# Create a memory-backed GDAL dataset from any AbstractGeoArray
function unsafe_ArchGDALdataset(A::AbstractGeoArray)
    width = size(A, Lon)
    height = size(A, Lat)
    nbands = size(A, Band)

    dataset = AG.unsafe_create("tmp",
                               driver=AG.getdriver("MEM"),
                               width=width,
                               height=height,
                               nbands=nbands,
                               dtype=eltype(A))
    # write bands to dataset
    AG.write!(dataset,
              data(permutedims(A, (Lon, Lat, Band))),
              Cint[1:nbands...])
    gdalsetproperties!(dataset, A)
    dataset
end

function AG.Dataset(f::Function, A::AbstractGeoArray)
    dataset = unsafe_ArchGDALdataset(A)
    try
        f(dataset)
    finally
        AG.destroy(dataset)
    end
end

#= Geotranforms ########################################################################

See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html

"In the particular, but common, case of a “north up” image without any rotation or
shearing, the georeferencing transform takes the following form" :
adfGeoTransform[0] /* top left x */
adfGeoTransform[1] /* w-e pixel resolution */
adfGeoTransform[2] /* 0 */
adfGeoTransform[3] /* top left y */
adfGeoTransform[4] /* 0 */
adfGeoTransform[5] /* n-s pixel resolution (negative value) */
=#

const GDAL_EMPTY_TRANSFORM = [0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
const GDAL_TOPLEFT_X = 1
const GDAL_WE_RES = 2
const GDAL_ROT1 = 3
const GDAL_TOPLEFT_Y = 4
const GDAL_ROT2 = 5
const GDAL_NS_RES = 6

isalligned(geotransform) =
    geotransform[GDAL_ROT1] == 0 && geotransform[GDAL_ROT2] == 0

geotransform2affine(gt) =
    AffineMap([gt[GDAL_WE_RES] gt[GDAL_ROT1]; gt[GDAL_ROT2] gt[GDAL_NS_RES]],
              [gt[GDAL_TOPLEFT_X], gt[GDAL_TOPLEFT_Y]])

function dims2geotransform(lat::Lat, lon::Lon)
    gt = zeros(6)
    gt[GDAL_TOPLEFT_X] = first(lon)
    gt[GDAL_WE_RES] = step(lon)
    gt[GDAL_ROT1] = 0.0
    gt[GDAL_TOPLEFT_Y] = first(lat) - step(lat)
    gt[GDAL_ROT2] = 0.0
    gt[GDAL_NS_RES] = step(lat)
    return gt
end


text:
# This file contains functions for the Directional Multiplier DEA model
"""
    DirectionalMultiplierDEAModel
An data structure representing a directional distance function multiplier DEA model.
"""
struct DirectionalMultiplierDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    Gx::Symbol
    Gy::Symbol
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    v::Matrix
    u::Matrix
    omega::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaddfm(X, Y; Gx, Gy)
Compute data envelopment analysis directional distance function multiplier model for inputs
`X` and outputs `Y`, using directions `Gx` and `Gy`.

# Direction specification:

The directions `Gx` and `Gy` can be one of the following symbols.
- `:Zeros`: use zeros.
- `:Ones`: use ones.
- `:Observed`: use observed values.
- `:Mean`: use column means.

Alternatively, a vector or matrix with the desired directions can be supplied.

# Optional Arguments
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.
"""
function deaddfm(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    Gx::Union{Symbol, Matrix, Vector}, Gy::Union{Symbol, Matrix, Vector},
    rts::Symbol = :CRS,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::DirectionalMultiplierDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    # Build or get user directions
    if typeof(Gx) == Symbol
        Gxsym = Gx

        if Gx == :Zeros
            Gx = zeros(size(X))
        elseif Gx == :Ones
            Gx = ones(size(X))
        elseif Gx == :Observed
            Gx = X
        elseif Gx == :Mean
            Gx = repeat(mean(X, dims = 1), size(X, 1))
        else
            throw(ArgumentError("Invalid `Gx`"));
        end

    else
        Gxsym = :Custom
    end

    if typeof(Gy) == Symbol
        Gysym = Gy

        if Gy == :Zeros
            Gy = zeros(size(Y))
        elseif Gy == :Ones
            Gy = ones(size(Y))
        elseif Gy == :Observed
            Gy = Y
        elseif Gy == :Mean
            Gy = repeat(mean(Y, dims = 1), size(Y, 1))
        else
            throw(ArgumentError("Invalid `Gy`"));
        end

    else
        Gysym = :Custom
    end

    if (size(Gx, 1) != size(X, 1)) | (size(Gx, 2) != size(X, 2))
        throw(DimensionMismatch("size of Gx and X ($(size(Gx)), $(size(X))) are not equal"));
    end
    if (size(Gy, 1) != size(Y, 1)) | (size(Gy, 2) != size(Y, 2))
        throw(DimensionMismatch("size of Gy and Y ($(size(Gy)), $(size(Y))) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    vi = zeros(n, m)
    ui = zeros(n, s)
    omegai = zeros(n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Directions to use
        Gx0 = Gx[i,:]
        Gy0 = Gy[i,:]

        # Solve if any direction is different from zero
        if any(Gx0 .!= 0) | any(Gy0 .!= 0)
            # Create the optimization model
            deamodel = newdeamodel(optimizer)

            @variable(deamodel, v[1:m] >= 0)
            @variable(deamodel, u[1:s] >= 0)
            @variable(deamodel, omega)

            @objective(deamodel, Min, - sum(u[r] * y0[r] for r in 1:s) + sum(v[i] * x0[i] for i in 1:m) + omega)

            @constraint(deamodel, [j in 1:nref], sum(u[r] * Yref[j,r] for r in 1:s) - sum(v[i] * Xref[j,i] for i in 1:m) - omega <= 0)
            @constraint(deamodel, sum(u[r] * Gy0[r] for r in 1:s) + sum(v[i] * Gx0[i] for i in 1:m) == 1)

            # Add return to scale constraints
            if rts == :CRS
                @constraint(deamodel, omega == 0)
            elseif rts == :VRS
                # No contraint to add for variable returns to scale
            else
                throw(ArgumentError("`rts` must be :CRS or :VRS"));
            end

            # Optimize and return results
            JuMP.optimize!(deamodel)

            effi[i]  = JuMP.objective_value(deamodel)
            vi[i, :] = JuMP.value.(v)
            ui[i, :] = JuMP.value.(u)
            omegai[i] = JuMP.value(omega)

            # Check termination status
            if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
                @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
            end
        else
            effi[i]  = 0.0
            vi[i,:] .= 0.0
            ui[i,:] .= 0.0
            omegai[i] = 0.0
        end

    end

    # Get first-stage X and Y targets
    Xtarget = X .- effi .* Gx
    Ytarget = Y .+ effi .* Gy

    if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
    if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end

    return DirectionalMultiplierDEAModel(n, m, s, Gxsym, Gysym, rts, names, effi, vi, ui, omegai, Xtarget, Ytarget)

end

function Base.show(io::IO, x::DirectionalMultiplierDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    v = multipliers(x, :X)
    u = multipliers(x, :Y)

    if !compact
        print(io, "Directional DF DEA Model (Multiplier form)\n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Gx = ", string(x.Gx), "; Gy = ", string(x.Gy))
        print(io, "\n")

        show(io, CoefTable(hcat(eff, v, u), ["efficiency"; ["v$i" for i in 1:m ]; ["u$i" for i in 1:s ]], dmunames))
    end
end


text:
# @everywhere ENV["JULIA_PKGDIR"] = "/work/fredrikb/.julia/"
@everywhere using Robotlib
@everywhere using Plots
@everywhere include("pf_state_estimator.jl")
@everywhere import PFstateEstimator
@everywhere include("simfuns.jl")
@everywhere import PFsimulator
@everywhere pyplot()
@everywhere default(size=(1000,700))
import DSP
@everywhere using Robotlib: Rt2T, ad, adi, T2R, T2t

MCruns          = 100
@assert MCruns >= 2
# @everywhere srand(1)

# Particle filter parameters
Npart       = round(Int,logspace(3,2,MCruns))#round(Int,logspace(3,2.3,MCruns))
h           = 0.01
sigmaW      = [0.3,  0.02π/180, 1e-6] # t, R, f
sigmaV      = [0.1, 0.1, 0.1π/180] # x, y, theta
dFK         = 1.2 # plateau width
kd          = 1e-3 # extension to plateau width due to force (essentialy a super simple compliance model)
drFK        = 1π/180
sigmaVFK    = [0.15, 1π/180, dFK, kd, drFK] # xyz, theta, plateau width, norm(f)
NpartInit   = 100
force_filter= 0.1 # in [0,1], lower value = more filtering

# Simulation parameters
T           = 500
σt          = 0.4
σR          = 0.02π/180
σq          = 0*0.0001π/180 # TODO: Verkar vara något vajsing på denna, även detta sjukt lilla värdet är högst significant, kanske ska man inte köra J *= 1000 utan /1000 på orienteringselementen
σm          = 0.1sigmaV
σf          = 1
force_ref   = [1,100,500,1,1,1] #
Kreact      = zeros(6,6)
Kcomp       = diagm([1e-3,1e-3,1e-3, 0.1π/180000, 0.1π/180000, 0.1π/180000]) # Compliance matrix
Kcomp[4,2]  = 0.5*π/180000 # Force in negative y dir creates deflection around x axis
R_TF_M      = [
-1 0 0;
0  0 1;
0  1 0]
# Kcomp      += 1e-6randn(size(Kcomp))
T_TF_M      = Rt2T(R_TF_M, [0,20.,0]) # Tool flange to meta
T_TF_M      = cat(3,T_TF_M, Rt2T(R_TF_M, [0.,40.,0]))

# Construct nominal trajectories
triwave     = cumsum(sign(sin(0.05*(1:T))))
triwave   ./= maximum(triwave)
triwave     = DSP.filtfilt(ones(15),[15.], triwave)
sinwave     = sin(0.005*(1:T))
traj_nom    = PFsimulator.traj2frames3D([5triwave linspace(0,200,T) zeros(T)]')
traj_nomyz  = PFsimulator.traj2frames3D([zeros(T) linspace(0,200,T) 100sinwave]', "x")

# Setup parameter objects
simparams   = PFsimulator.SimParams(h, T, σt, σR, σq, σm, σf, force_ref, Kreact, Kcomp, T_TF_M)

# Setup main loop
seeds               = rand(1:1000,MCruns)
simresult           = Array(Any,MCruns)
simresultNoMeas     = Array(Any,MCruns)
simresult1Sensor    = Array(Any,MCruns)

simresultyz         = Array(Any,MCruns)
simresultyzNoMeas   = Array(Any,MCruns)
simresultyz1Sensor  = Array(Any,MCruns)
savetrace           = MCruns <= 20
factor              = 1

tic()
# Main loop, perform all simulations asynchronously
@sync for mc = 1:MCruns

    # Generate simulation data
    simparamsi           = PFsimulator.perturb_simparams(simparams, factor)
    simdata              = PFsimulator.generate_simdata(traj_nom, simparamsi, optimize = false)
    simdataNoMeas        = deepcopy(simdata)
    simdataNoMeas.meas  *= 0
    simparams1Sensor     = deepcopy(simparamsi)
    simparams1Sensor.T_TF_M = T_TF_M[:,:,1]
    pfparams             = PFstateEstimator.PFparams(sigmaW, sigmaV, sigmaVFK, Npart[mc], NpartInit, h, force_filter)
    PFsimulator.perturb_pfparams!(pfparams, factor)

    # Generate simulation data
    simdatayz              = PFsimulator.generate_simdata(traj_nomyz, simparamsi, optimize = false)
    simdatayzNoMeas        = deepcopy(simdatayz)
    simdatayzNoMeas.meas  *= 0

    # Run simulations
    simresult[mc]        = @spawn PFsimulator.run_tracking("xy 2 sens", simdata,       simparamsi,
    pfparams,   seeds[mc], savetrace = savetrace)
    simresultNoMeas[mc]  = @spawn PFsimulator.run_tracking("xy 0 sens",     simdataNoMeas, simparamsi,
    pfparams,   seeds[mc], savetrace = savetrace)
    simresult1Sensor[mc] = @spawn PFsimulator.run_tracking("xy 1 sens", simdata, simparams1Sensor,
    pfparams,   seeds[mc], savetrace = savetrace)

    simresultyz[mc]        = @spawn PFsimulator.run_tracking("yz 2 sens", simdatayz,       simparamsi,
    pfparams,   seeds[mc], savetrace = savetrace)
    simresultyzNoMeas[mc]  = @spawn PFsimulator.run_tracking("yz 0 sens",     simdatayzNoMeas, simparamsi,
    pfparams,   seeds[mc], savetrace = savetrace)
    simresultyz1Sensor[mc] = @spawn PFsimulator.run_tracking("yz 1 sens", simdatayz, simparams1Sensor,
    pfparams,   seeds[mc], savetrace = savetrace)

    println("Spawned MCrun ", mc)
end

# Fetch result of simulations
for mc = 1:MCruns
    simresult[mc]        = fetch(simresult[mc])
    simresultNoMeas[mc]  = fetch(simresultNoMeas[mc])
    simresult1Sensor[mc] = fetch(simresult1Sensor[mc])

    simresultyz[mc]        = fetch(simresultyz[mc])
    simresultyzNoMeas[mc]  = fetch(simresultyzNoMeas[mc])
    simresultyz1Sensor[mc] = fetch(simresultyz1Sensor[mc])
end
time = toq()
println("Done! It took $time seconds")

# Plot results

PFsimulator.plot_errors(simresultNoMeas, simresult1Sensor, simresult, simresultyzNoMeas, simresultyz1Sensor, simresultyz)
# PFsimulator.savetikz("/local/home/fredrikb/papers/PFseamTracking/errorbox.tex", PyPlot.gcf(), ["x tick label style={rotate=45,anchor=east}"])

# PFsimulator.plot(simresult[1], :tracking, :resample)
# PFsimulator.plot(simresult1Sensor[1], :tracking, :resample, :angle)
# PFsimulator.plot(simresultNoMeas[1], :tracking, :resample, :angle)

# plot_traj_sub(simresult[1].simdata.traj_nom-simresult[1].simdata.traj_act, title="Noise in traj_act")
# plot_traj_sub(simresult[1].simdata.traj_nom-simresult[1].simdata.traj_meas, title="Noise in traj_meas")

# Linear modeling
using GLM, DataFrames
import ExperimentalAnalysis
df      = PFsimulator.get_dataframe(simresultNoMeas, simresult1Sensor, simresult, simresultyzNoMeas, simresultyz1Sensor, simresultyz);
# ExperimentalAnalysis.scattermatrix(df[(df[:sens] .>= 1) & (df[:traj_type] .== "yz"),:],ex +ey+ez+er +sens ~ Npart + σW1 + σW2 + σW3 + σV1+ σV2+ σV3)
# ExperimentalAnalysis.scattermatrix(df[df[:sens] .>= 1,:],ex +ey+ez+er ~ σVFK1+ σVFK2+ σVFK3+ σVFK4+ σVFK5 )
# ExperimentalAnalysis.scattermatrix(df[df[:sens] .>= 1,:],ex +ey+ez+er ~ force_r1 + force_r2 + force_r3 + force_r4 + force_r5 + force_r6 )
# ExperimentalAnalysis.scattermatrix(df[df[:sens] .>= 1,:],ex +ey+ez+er+sens ~  σm1 + σm2 + σm3 + σf + sens)

modelx  = lm(ex~Npart + sens + et0 + σW1 + σW2 + σW3+ σm1+ σm2+ σm3 + σV1+ σV2+ σV3+ σVFK1+ σVFK2+ σVFK3+ σVFK4+ σVFK5 + force_r1 + force_r2 + force_r3 + force_r4 + force_r5 + force_r6 + σf, df)
modely  = lm(ey~Npart + sens + et0 + σW1 + σW2 + σW3+ σm1+ σm2+ σm3 + σV1+ σV2+ σV3+ σVFK1+ σVFK2+ σVFK3+ σVFK4+ σVFK5 + force_r1 + force_r2 + force_r3 + force_r4 + force_r5 + force_r6 + σf, df)
modelz  = lm(ez~Npart + sens + et0 + σW1 + σW2 + σW3+ σm1+ σm2+ σm3 + σV1+ σV2+ σV3+ σVFK1+ σVFK2+ σVFK3+ σVFK4+ σVFK5 + force_r1 + force_r2 + force_r3 + force_r4 + force_r5 + force_r6 + σf, df)
modelr  = lm(er~Npart + sens + et0 + σW1 + σW2 + σW3+ σm1+ σm2+ σm3 + σV1+ σV2+ σV3+ σVFK1+ σVFK2+ σVFK3+ σVFK4+ σVFK5 + force_r1 + force_r2 + force_r3 + force_r4 + force_r5 + force_r6 + σf, df)

# modelx  = lm(ex~Npart + sens + et0 + traj_type, df)
# modely  = lm(ey~Npart + sens + et0 + traj_type, df)
# modelz  = lm(ez~Npart + sens + et0 + traj_type, df)
# modelr  = lm(er~Npart + sens + et0 + traj_type, df)

modeld = []
for sym in [:ex, :ey, :ez, :er]
    f = $sym + (sens .== 1) - (sens .== 2) ~ Npart + et0 + traj_type + σW1*σV1 + σW2 + σW3 + σm1+ σm2+ σm3 + σV2+ σV3+ σVFK1+ σVFK2+ σVFK3+ σVFK4+ σVFK5 + σf
    push!(modeld, PFsimulator.dglm(f,df, Normal()))
    println("=== Linear modeling of differences for $sym ===")
    display(modeld[end])
end

fig = ExperimentalAnalysis.modelheatmap(["x", "y", "z", "R"], modelx, modely, modelz, modelr)
# PFsimulator.savetikz("/local/home/fredrikb/papers/PFseamTracking/heatmap.tex", PyPlot.gcf(), ["x tick label style={rotate=45,anchor=east}"])

ExperimentalAnalysis.modelheatmap(["ex", "dy", "dz", "dR"], modeld...)
# PFsimulator.savetikz("/local/home/fredrikb/papers/PFseamTracking/heatmapdiff.tex", PyPlot.gcf(), ["x tick label style={rotate=45,anchor=east}"])

scatter(df[:Npart][(df[:sens2].==true) & (df[:traj_type].=="xy")],df[:ey][(df[:sens2].==true) & (df[:traj_type].=="xy")])

# file = open("/work/fredrikb/savefile_varyN","w")
# serialize(file,(simresultNoMeas, simresult1Sensor, simresult, simresultyzNoMeas, simresultyz1Sensor, simresultyz))
# close(file)

# file = open("/work/fredrikb/savefile_varyN","r")
# (simresultNoMeas, simresult1Sensor, simresult, simresultyzNoMeas, simresultyz1Sensor, simresultyz) = deserialize(file);
# close(file)
#

savetrace && PFsimulator.plot_trace(simresult[1], xIndices = 1:3, yIndices = 1:3)
pyplot()


text:
using LinearAlgebra
@testset "topo order" begin
    g = Compgraph()
    add_lincomb!(g, :X1, 0.2, :I, 0.2, :A)
    add_lincomb!(g, :X2, 0.2, :I, 0.22, :A)
    add_lincomb!(g, :X3, 0.3, :I, 0.1, :A)

    add_mult!(g, :X4, :I, :X3)

    add_mult!(g, :Z1, :X1, :X2)
    add_mult!(g, :Z2, :X2, :X3)
    add_mult!(g, :Z3, :Z1, :Z2)
    add_mult!(g, :Z4, :Z3, :Z3)

    helper = Dict()
    helper[:X1] = -10000 + 1
    helper[:X2] = -10000 + 2
    helper[:X3] = -10000 + 3
    helper[:X4] = -3

    add_output!(g, :Z4)
    order = get_topo_order(g, priohelp = helper, will_not_dealloc = [])[1]
    # Check that 1) helper prioritizes x1,x2,x3 first
    # 2) that x4 comes next since it "frees" :I (which is not in will_not_dealloc)
    @test order[1:4] == [:X1; :X2; :X3; :X4]
    # 3) that z1 comes next since it "frees" if will_not_dealloc=:I
    #    (and before z2 since it has better priority)
    helper[:Z1] = -2
    helper[:Z2] = -1
    order = get_topo_order(g, priohelp = helper, will_not_dealloc = [:I])[1]
    @test order[1:4] == [:X1; :X2; :X3; :Z1]
end


text:
"""
A datastructure which in addition to a list of faces stores connectivity information for each vertex.
"""
struct CachedDS
    faces
    connectivity
end

Faces(t::CachedDS) = t.faces
Edges(t::CachedDS) = filter(x->x[1]<x[2],decompose(Face{2,Int},t.faces))

"""
    CacheDS(t)

Constructs cached face based datastructure `CachedDS` from arbitrary topology `t` which provides `EdgeRing` iterator.
"""
function CachedDS(faces)
    vmax = maximum(maximum(faces))
    connectivity = Array{Int,1}[]
    for i in 1:vmax
        v1 = Array{Int}(undef,0)
        v2 = Array{Int}(undef,0)
        for (v1i,v2i) in EdgeRing(i,faces)
            push!(v1,v1i)
            push!(v2,v2i)
        end

        vj = v2[1]
        coni = [vj]
        for j in 2:length(v1)
            vj, = v2[v1.==vj]
            push!(coni,vj)
        end
        push!(connectivity,coni)
    end
    return CachedDS(faces,connectivity)
end

VertexRing(vi::Int64,t::CachedDS) = t.connectivity[vi]
EdgeRing(vi::Int64,t::CachedDS) = PairIterator(t.connectivity[vi])
FaceRing(vi::Int64,t::CachedDS) = error("Face ring is not suitable for this datastructure.")


text:
# Cycle through regularized evolution many times,
# printing the fittest equation every 10% through
function SRCycle(dataset::Dataset{T}, baseline::T,
        pop::Population,
        ncycles::Int,
        curmaxsize::Int,
        frequencyComplexity::AbstractVector{T};
        verbosity::Int=0,
        options::Options
        )::Population where {T<:Real}

    top = convert(T, 1)
    allT = LinRange(top, convert(T, 0), ncycles)
    for temperature in 1:size(allT)[1]
        if options.annealing
            pop = regEvolCycle(dataset, baseline, pop, allT[temperature], curmaxsize, frequencyComplexity, options)
        else
            pop = regEvolCycle(dataset, baseline, pop, top, curmaxsize, frequencyComplexity, options)
        end

        if verbosity > 0 && (temperature % verbosity == 0)
            bestPops = bestSubPop(pop)
            bestCurScoreIdx = argmin([bestPops.members[member].score for member=1:bestPops.n])
            bestCurScore = bestPops.members[bestCurScoreIdx].score
            debug(verbosity, bestCurScore, " is the score for ", stringTree(bestPops.members[bestCurScoreIdx].tree, options, varMap=dataset.varMap))
        end
    end

    return pop
end


text:
immutable PointwiseLayer{P<:AbstractPointwise} <: Mocha.Layer
    ρ::P
    bottoms::Vector{Symbol}
    name::AbstractString
    tops::Vector{Symbol}
end

Mocha.can_do_bp(::PointwiseLayer) = false
Mocha.has_neuron(::PointwiseLayer) = false
Mocha.has_param(::PointwiseLayer) = false

function PointwiseLayer( ;
        ρ::AbstractPointwise = Identity(),
        bottoms::Vector{Symbol} = Symbol[],
        name::AbstractString = "pointwise",
        tops::Vector{Symbol} = Symbol[])
    @assert length(bottoms) == 1
    @assert length(tops) == 1
    PointwiseLayer(ρ, bottoms, name, tops)
end


text:
using LinearAlgebra
import LinearAlgebra.Eigen
using SparseArrays
#this lib will be used to estimate the executon time

struct SVDCM
	U::Matrix{Float64}
	S::SparseMatrixCSC{Float64, Int64}
	S_inv::SparseMatrixCSC{Float64, Int64}
	V::Matrix{Float64}

	function SVDCM(A::Matrix{Float64})
		U, S, V,S_inv = svdCm(A)
		return new( U, S, S_inv,V)
	end
	function svdCm(A::Matrix{Float64})

		r,c = size(A)
		shorter=min(r,c)
		longer=max(r,c)

		Asq = zeros(shorter,shorter)
		mul!(Asq, A', A)
		a, V = eigen!(Asq)

		S = sparse(zeros(longer,shorter))
		Sin = sparse(zeros(shorter,longer))

		S[diagind(S)] .= (sqrt.(a))
		Sin[diagind(Sin)] .= 1 ./ S[diagind(S)]

		U1 = zeros(r,c)
		U = zeros(longer,longer)

		mul!(U1, A, V)
		mul!(U,U1,Sin)
		return U, S, V, Sin
	end
end


text:
export  PlaneStrain

@doc raw"""
    PlaneStrain

Creates a plane strain element

- `H`: Linear elasticity matrix, $3\times3$
- `E`: Young's modulus
- `ν`: Poisson's ratio
- `ρ`: density
- `σ0`: stress at the **last** time step
- `σ0_`: (for internal use), stress to be updated in `commitHistory`
- `ε0`: strain at the **last** time step
- `ε0_`: (for internal use), strain to be updated in `commitHistory`

# Example
```julia
prop = Dict("name"=> "PlaneStrain", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
mat = PlaneStrain(prop)
```
"""
mutable struct PlaneStrain
    H::Array{Float64}
    E::Float64
    ν::Float64
    ρ::Float64
    σ0::Array{Float64}
    σ0_::Array{Float64}
    ε0::Array{Float64}
    ε0_::Array{Float64}
end

"""
    PlaneStrain(prop::Dict{String, Any})

`prop` should contain at least the following three fields: `E`, `nu`, `rho`
"""
function PlaneStrain(prop::Dict{String, Any})
    E = prop["E"]; ν = prop["nu"]; ρ = prop["rho"]
    H = zeros(3,3)
    H[1,1] = E*(1. -ν)/((1+ν)*(1. -2. *ν));
    H[1,2] = H[1,1]*ν/(1-ν);
    H[2,1] = H[1,2];
    H[2,2] = H[1,1];
    H[3,3] = H[1,1]*0.5*(1. -2. *ν)/(1. -ν);
    σ0 = zeros(3); σ0_ = zeros(3); ε0 = zeros(3); ε0_ = zeros(3)
    PlaneStrain(H, E, ν, ρ, σ0, σ0_,ε0,ε0_)
end

function getStress(self::PlaneStrain, strain::Array{Float64}, Dstrain::Array{Float64}, Δt::Float64 = 0.0)
    sigma = self.H * strain
    self.σ0_ = copy(sigma)
    self.ε0_ = copy(strain)
    return sigma, self.H
end

function getTangent(self::PlaneStrain)
    self.H
end

function commitHistory(self::PlaneStrain)
    self.σ0 = self.σ0_
    self.ε0 = self.ε0_
end


text:
"""
    FiniteDiffBackend(method::Val{Symbol} = Val{:central})

Differentiation backend based on FiniteDiff package.
"""
struct FiniteDiffBackend{TM<:Val} <: AbstractDiffBackend
    method::TM
end

FiniteDiffBackend() = FiniteDiffBackend(Val(:central))

push!(_diff_backends, FiniteDiffBackend())

function _derivative(f, p, backend::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_derivative(f, p, Method)
end

function _gradient(f, p, backend::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_gradient(f, p, Method)
end

function _gradient!(f, X, p, backend::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_gradient!(X, f, p, Method)
end


text:
include("TotalVariation.jl")
using TotalVariation
using PyPlot

ground_truth = vcat(ones(1000),
-10sin.(linspace(0,pi,400)), 1ones(950))
noise        = 10randn(length(ground_truth))
combined     = ground_truth .+ noise
gstv_out     = TotalVariation.gstv(combined, 40, 15.0)
tv_out       = TotalVariation.tv(combined, 200.0)

figure(1)
PyPlot.clf();
subplot(2,2,1)
title("original signal")
plot(ground_truth)

subplot(2,2,2)
title("original + noise")
plot(combined)

subplot(2,2,3)
title("tv correction")
plot(tv_out)
plot(ground_truth, color="g")

subplot(2,2,4)
title("gstv correction")
plot(gstv_out)
plot(ground_truth, color="g")

rms(s) = sqrt(mean((ground_truth.-s).^2))

println("Initial Error        = ", rms(combined))
println("Output  Error [tv]   = ", rms(tv_out))
println("Output  Error [gstv] = ", rms(gstv_out))


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using .Gherkin: DataTable, DataTableRow

@testset "Executor Datatables  " begin
    @testset "Scenario step has a data table; Context contains the data tables" begin
        table = DataTable(
            DataTableRow[
                ["header1", "header2"],
                ["value 11", "value 12"],
                ["value 21", "value 22"],
            ]
        )
        given = Given("Some precondition", datatable=table)
        function check_datatable_step_definition(context::StepDefinitionContext, _args)
            expectedtable = DataTable(
                DataTableRow[
                    ["header1", "header2"],
                    ["value 11", "value 12"],
                    ["value 21", "value 22"],
                ]
            )

            if context.datatable == expectedtable
                Behavior.SuccessfulStepExecution()
            else
                Behavior.StepFailed("")
            end
        end
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => check_datatable_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
    end

    @testset "Scenario step has no data table; Context table is empty" begin
        given = Given("Some precondition")
        function check_datatable_step_definition(context::StepDefinitionContext, _args)
            if context.datatable == []
                Behavior.SuccessfulStepExecution()
            else
                Behavior.StepFailed("")
            end
        end
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => check_datatable_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
    end

    @testset "First step has table but not second; Second context table is empty" begin
        given1 = Given("Some precondition")
        given2 = Given("Some other precondition")
        function check_datatable_step_definition(context::StepDefinitionContext, _args)
            if context.datatable == []
                Behavior.SuccessfulStepExecution()
            else
                Behavior.StepFailed("")
            end
        end
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(
            given1 => successful_step_definition,
            given2 => check_datatable_step_definition
        ))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given1, given2])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[2], Behavior.SuccessfulStepExecution)
    end

    @testset "Scenario Outline step has a data table; Context contains the data tables" begin
        table = DataTable(
            DataTableRow[
                ["header1", "header2"],
            ]
        )
        given = Given("Some precondition <foo>", datatable=table)
        when = When("Some precondition", datatable=table)
        then = Then("Some precondition", datatable=table)
        function check_datatable_step_definition(context::StepDefinitionContext, _args)
            expectedtable = DataTable(
                DataTableRow[
                    ["header1", "header2"],
                ]
            )

            if context.datatable == expectedtable
                Behavior.SuccessfulStepExecution()
            else
                Behavior.StepFailed("")
            end
        end
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(
            Given("Some precondition 1") => check_datatable_step_definition,
            Given("Some precondition 2") => check_datatable_step_definition,
            Given("Some precondition 3") => check_datatable_step_definition,
            when => check_datatable_step_definition,
            then => check_datatable_step_definition,
        ))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = ScenarioOutline(
            "Description",
            String[],                        # Tags
            ScenarioStep[given, when, then], # Steps
            ["foo"],                         # Placeholders
            [ "1" "2"]                       # Examples
        )

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult[1].steps[1], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult[1].steps[2], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult[1].steps[3], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult[2].steps[1], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult[2].steps[2], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult[2].steps[3], Behavior.SuccessfulStepExecution)
    end

end


text:
"""
    _SubSolver <: AbstractSolver

An internal solver type called by MetaSolver when multithreading is enabled.

# Arguments:
- `id::Int`: subsolver id for debugging
- `model::Model`: a ref to the model of the main solver
- `state::_State`: a `deepcopy` of the main solver that evolves independently
- `options::Options`: a ref to the options of the main solver
"""
mutable struct _SubSolver <: AbstractSolver
    meta_local_id::Tuple{Int, Int}
    model::_Model
    options::Options
    pool::Pool
    state::State
    strategies::MetaStrategy
end

function solver(mlid, model, options, pool, ::RemoteChannel, ::RemoteChannel, ::RemoteChannel, strats, ::Val{:sub})
    sub_options = deepcopy(options)
    set_option!(options, "print_level", :silent)
    return _SubSolver(mlid, model, sub_options, pool, state(), strats)
end

_init!(s::_SubSolver) = _init!(s, :local)

stop_while_loop(::_SubSolver, stop, ::Int, ::Float64) = !(stop[])


text:
#--- JalerkinType

"""Abstract type JalerkinType: root of the type hierarchy"""
abstract type JalerkinType end

#--- Floating point precision

"""Floating point precision to be used by Jalerkin"""
fp_precision() = Float64

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
using Test
using Distributions: Categorical, kldivergence
@test kldivergence(Categorical([0.0, 0.1, 0.9]), Categorical([0.1, 0.1, 0.8])) ≥ 0
@test kldivergence(Categorical([0.0, 0.1, 0.9]), Categorical([0.1, 0.1, 0.8])) ≈
    kldivergence([0.0, 0.1, 0.9], [0.1, 0.1, 0.8])


text:
module Seglib

export
    randindex,
    randerror,
    contingencymatrix,
    seg2boundary,
    seg2boundaryindex,
    seg2boundarystrength,
    spdist

include("superpixel.jl")
include("metrics.jl")

matlab_homepath = begin
    try
        dirname(dirname(realpath(chomp(readall(`which matlab`)))))
    catch
        nothing
    end
end

#if !is(matlab_homepath, nothing)
#    include("matlab.jl")
#
#    export
#        @mquiet,
#        slic,
#        thinboundary,
#        correspondpixels,
#        gPb,
#        ucm_regions,
#        ucm_strength
#end

end # module


text:
module SimplexifyTests

using Gridap

grid = CartesianGrid(domain=(0.0,1.0,-1.0,2.0),partition=(3,4))
tgrid = simplexify(grid)
test_grid(tgrid,20,24)
#writevtk(tgrid,"tgrid")

grid = CartesianGrid(partition=(3,4,3))
tgrid = simplexify(grid)
test_grid(tgrid,80,216)
#writevtk(tgrid,"tgrid")

model = CartesianDiscreteModel(partition=(3,4,2))
tmodel = simplexify(model)

model = CartesianDiscreteModel(partition=(3,4))
tmodel = simplexify(model)
#writevtk(tmodel,"tmodel")

end # module


text:
mutable struct Iterator
    val::Any
    st::Any
end
Iterator(x) = Iterator(x, nothing)
Base.length(x::Iterator) = length(x.val)

const PyJuliaIteratorValue_Type__ref = Ref(PyNULL)
PyJuliaIteratorValue_Type() = begin
    ptr = PyJuliaIteratorValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.IteratorValue",
                base = base,
                iter = pyjliter_iter,
                iternext = pyjliter_iternext,
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaIteratorValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaIteratorValue_New(x::Iterator) = PyJuliaValue_New(PyJuliaIteratorValue_Type(), x)

pyjliter_iter(xo::PyPtr) =
    PyJuliaIteratorValue_New(Iterator((PyJuliaValue_GetValue(xo)::Iterator).val))

pyjliter_iternext(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::Iterator
        val = x.val
        st = x.st
        if st === nothing
            z = iterate(val)
        else
            z = iterate(val, something(st))
        end
        if z === nothing
            PyNULL
        else
            r, newst = z
            x.st = Some(newst)
            PyObject_From(r)
        end
    catch err
        if err isa MethodError && err.f === iterate
            PyErr_SetStringFromJuliaError(PyExc_TypeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end


text:
### LINEAR SOLVE FUNCTIONS ####
# Direct solvers
function linear_solve!(u,K,RHS,::BackslashSolver)
    u .= K \ RHS
    return nothing
end
linear_solve(K,RHS,::BackslashSolver) = K \ RHS

function linear_solve!(u,K,RHS,::MUMPS)
    u .= MUMPSjInv.solveMUMPS(K, RHS, 1, 1)
    return nothing
end

# indirect solvers :
function linear_solve!(δu,K,res,linear_solver::ConjugateGradients)
    if linear_solver.package == :KrylovMethods
        #!isdefined(KrylovMethods) && @eval(import KrylovMethods)
        δu, flag, relres, iter, resvec = KrylovMethods.cg(K, res, linear_solver.kwargs...)
        @assert flag == 0
        return nothing
    elseif linear_solver.package == :IterativeSolvers
        #!isdefined(IterativeSolvers) && @eval(import IterativeSolvers)
        #Preconditioners.UpdatePreconditioner!(p, K)
        IterativeSolvers.cg!(δu, K, res; linear_solver.kwargs...)
        return nothing
    else
        @error "linear solver package $(linear_solver.package) is not implemented."
    end
end

function linear_solve!(δu,K,res,linear_solver::MINRES)
    if linear_solver.package == :IterativeSolvers
        IterativeSolvers.minres!(δu, K, res; linear_solver.kwargs...)
        return nothing
    else
        @error "linear solver package $(linear_solver.package) is not interfaced."
    end
end


text:
@testset "Bounding volumes" begin

@testset "Create bounding box" begin

    p1 = Point3(0, 0, 0)
    p2 = Point3(1, 2, 3)
    p3 = Point3(-3, -2, -1)

    @testset "Degenerate bounding box" begin
        BB = GRAYTR.BoundingBox(p1, p1)
        @test BB.pMin == p1
        @test BB.pMax == p1

        BB = GRAYTR.BoundingBox(p2, p2)
        @test BB.pMin == p2
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p3, p3)
        @test BB.pMin == p3
        @test BB.pMax == p3
    end

    @testset "Bounding box from two points" begin
        BB = GRAYTR.BoundingBox(p1, p2)
        @test BB.pMin == p1
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p2, p1)
        @test BB.pMin == p1
        @test BB.pMax == p2
    end

    @testset "Bounding box from three points" begin
        BB = GRAYTR.BoundingBox(p1, p2, p3)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p1, p3, p2)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p2, p1, p3)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p2, p3, p1)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p3, p1, p2)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p3, p2, p1)
        @test BB.pMin == p3
        @test BB.pMax == p2
    end

    @testset "Bounding box from array of points" begin
        BB = GRAYTR.BoundingBox([p1, p2, p3])
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox([p1, p3, p2])
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox([p2, p1, p3])
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox([p2, p3, p1])
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox([p3, p1, p2])
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox([p3, p2, p1])
        @test BB.pMin == p3
        @test BB.pMax == p2
    end

    @testset "Bounding box from two bounding boxes" begin
        BB1 = GRAYTR.BoundingBox(p1, p2)
        BB2 = GRAYTR.BoundingBox(p3, p2)

        BB = GRAYTR.BoundingBox(BB1, BB2)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(BB2, BB1)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox([BB1, BB2])
        @test BB.pMin == p3
        @test BB.pMax == p2
    end

    @testset "Bounding box from a bounding box and a point" begin
        BB1 = GRAYTR.BoundingBox(p1, p2)

        BB = GRAYTR.BoundingBox(BB1, p3)
        @test BB.pMin == p3
        @test BB.pMax == p2

        BB = GRAYTR.BoundingBox(p3, BB1)
        @test BB.pMin == p3
        @test BB.pMax == p2
    end
end

@testset "Area" begin
    p1 = Point3(0, 0, 0)
    p2 = Point3(1, 1, 1)
    p3 = Point3(-3, -2, -1)

    BB = GRAYTR.BoundingBox(p1, p2)
    @test GRAYTR.area(BB) ≈ 6.0

    BB = GRAYTR.BoundingBox(p1, p1)
    @test GRAYTR.area(BB) ≈ 0.0

    BB = GRAYTR.BoundingBox(p1, p3)
    @test GRAYTR.area(BB) ≈ 2*2.0 + 2*3.0 + 2*6.0
end

@testset "Max extent" begin
    p1 = Point3(0, 0, 0)
    p2 = Point3(1, 2, 3)
    p3 = Point3(-3, -2, -1)
    p4 = Point3(1, 3, 2)

    BB = GRAYTR.BoundingBox(p1, p2)
    @test GRAYTR.max_extent(BB) == 3

    BB = GRAYTR.BoundingBox(p1, p3)
    @test GRAYTR.max_extent(BB) == 1

    BB = GRAYTR.BoundingBox(p2, p3)
    @test GRAYTR.max_extent(BB) == 3

    BB = GRAYTR.BoundingBox(p1, p1)
    @test GRAYTR.max_extent(BB) == 3

    BB = GRAYTR.BoundingBox(p1, p4)
    @test GRAYTR.max_extent(BB) == 2
end

@testset "True/false ray intersection test" begin
    p1 = Point3(0, 0, 0)
    p2 = Point3(2, 2, 1)
    BB = GRAYTR.BoundingBox(p1, p2)

    @testset "Ray from inside box" begin
        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(1, 0, 0))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, 1, 0))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, 0,  1))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(-1, 0, 0))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, -1, 0))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, 0, -1))
        @test GRAYTR.intersectP(R, BB)
    end

    @testset "Ray travelling in -Z direction" begin
        R = GRAYTR.Ray(Point3(0.0, 0.0, 2.0), Vector3(0, 0, -1))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(1.0, 1.0, 2.0), Vector3(0, 0, -1))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(2.0, 2.0, 2.0), Vector3(0, 0, -1))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(-1.0, 1.0, 2.0), Vector3(0, 0, -1))
        @test !GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(3.0, 1.0, 2.0), Vector3(0, 0, -1))
        @test !GRAYTR.intersectP(R, BB)
    end

    @testset "Ray travelling +X direction" begin
        R = GRAYTR.Ray(Point3(-1.0, 1.0, 0.5), Vector3(1, 0, 0))
        @test GRAYTR.intersectP(R, BB)

        R = GRAYTR.Ray(Point3(-1.0, 3.0, 0.5), Vector3(1, 0, 0))
        @test !GRAYTR.intersectP(R, BB)
    end

end

@testset "Full ray intersection test" begin
    p1 = Point3(0, 0, 0)
    p2 = Point3(2, 2, 1)
    BB = GRAYTR.BoundingBox(p1, p2)

    @testset "Ray from inside box" begin
        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(1, 0, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 0.0
        @test t1 ≈ 1.0

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, 1, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 0.0
        @test t1 ≈ 1.0

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, 0,  1))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 0.0
        @test t1 ≈ 0.5

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(-1, 0, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 0.0
        @test t1 ≈ 1.0

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, -1, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 0.0
        @test t1 ≈ 1.0

        R = GRAYTR.Ray(Point3(1.0, 1.0, 0.5), Vector3(0, 0, -1))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 0.0
        @test t1 ≈ 0.5
    end

    @testset "Ray travelling in -Z direction" begin
        R = GRAYTR.Ray(Point3(0.0, 0.0, 2.0), Vector3(0, 0, -1))
        isect, t0, t1 =  GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 1.0
        @test t1 ≈ 2.0

        R = GRAYTR.Ray(Point3(1.0, 1.0, 2.0), Vector3(0, 0, -1))
        isect, t0, t1 =  GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 1.0
        @test t1 ≈ 2.0

        R = GRAYTR.Ray(Point3(-1.0, 1.0, 2.0), Vector3(0, 0, -1))
        isect, t0, t1 =  GRAYTR.intersect(R, BB)
        @test !isect
        @test isnan(t0)
        @test isnan(t1)

        R = GRAYTR.Ray(Point3(3.0, 1.0, 2.0), Vector3(0, 0, -1))
        isect, t0, t1 =  GRAYTR.intersect(R, BB)
        @test !isect
        @test isnan(t0)
        @test isnan(t1)
    end

    @testset "Ray travelling +X direction" begin
        R = GRAYTR.Ray(Point3(-1.0, 1.0, 0.5), Vector3(1, 0, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 1.0
        @test t1 ≈ 3.0

        R = GRAYTR.Ray(Point3(-1.0, 3.0, 0.5), Vector3(1, 0, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test !isect
        @test isnan(t0)
        @test isnan(t1)
    end

    @testset "Ray travelling -X direction" begin
        R = GRAYTR.Ray(Point3(3.0, 1.0, 0.5), Vector3(-1, 0, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test isect
        @test t0 ≈ 1.0
        @test t1 ≈ 3.0

        R = GRAYTR.Ray(Point3(3.0, 3.0, 0.5), Vector3(-1, 0, 0))
        isect, t0, t1 = GRAYTR.intersect(R, BB)
        @test !isect
        @test isnan(t0)
        @test isnan(t1)
    end

end

@testset "Bounding Sphere" begin
    @testset "Creating" begin
        BSph = GRAYTR.BoundingSphere(Point3(0), 1.0)
        @test BSph.center ≈ Point3(0)
        @test BSph.radius ≈ 1.0

        p1 = Point3(0, 0, 0)
        p2 = Point3(1, 2, 4)

        d = sqrt(1 + 4 + 16)

        BB = GRAYTR.BoundingBox(p1, p2)

        BSph = GRAYTR.BoundingSphere(BB)
        @test BSph.center ≈ Point3(0.5, 1.0, 2.0)
        @test isapprox(BSph.radius, d/2 ; atol=1e-14)
    end

    @testset "Ray intersection" begin
        BSph = GRAYTR.BoundingSphere(Point3(0), 1.0)

        R = GRAYTR.Ray(Point3(0.0, 0.0, 3.0), Vector3(0.0, 0.0, -1.0))
        @test GRAYTR.intersectP(R, BSph)

        R = GRAYTR.Ray(Point3(0.0, 0.0, -3.0), Vector3(0.0, 0.0, 1.0))
        @test GRAYTR.intersectP(R, BSph)

        R = GRAYTR.Ray(Point3(0.0, 0.0, -3.0), Vector3(0.0, 0.0, -1.0))
        @test !GRAYTR.intersectP(R, BSph)

        R = GRAYTR.Ray(Point3(3.0, 0.0, 3.0), Vector3(0.0, 0.0, -1.0))
        @test !GRAYTR.intersectP(R, BSph)

    end

end

end # main testset

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR