text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using BinaryBuilderBase

include(joinpath(@__DIR__, "..", "..", "platforms", "microarchitectures.jl"))

name = "finufft"
version = v"2.0.4"
julia_compat = "1.6"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/flatironinstitute/finufft/archive/v$(version).zip", "2434f694b4fbdbeb65c77f65d784a1712852130b9c61e15999555a2e2cf1a9fa")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/finufft*/

CFLAGS="-fopenmp -fPIC -O3 -funroll-loops -Iinclude"
if [[ "${target}" != *-freebsd* ]] && [[ "${target}" != *-apple-* ]]; then
    CFLAGS="${CFLAGS} -fcx-limited-range"
fi

# Overwrite LIBSFFT such that we do not require fftw3_threads or fftw3_omp for OMP support. Since the libraries in FFTW_jll already provide for threading, we do not loose anything.
# Make use of the -DFFTW_PLAN_SAFE flag to allow for multiple threads using finufft at the same time.
make lib \
    CC=${CC} \
    CXX=${CXX} \
    CFLAGS="${CFLAGS}" \
    CXXFLAGS="${CFLAGS} -std=c++14 -DFFTW_PLAN_SAFE" \
    LIBSFFT="-lfftw3 -lfftw3f -lm" \
    DYNLIB="lib/libfinufft.${dlext}"
install -Dvm 0755 "lib/libfinufft.${dlext}" "${libdir}/libfinufft.${dlext}"
"""

# Expand for microarchitectures on x86_64 (library doesn't have CPU dispatching)
# Tests on Linux/x86_64 yielded a slow binary with avx512 for some reason, so disable that
platforms = expand_cxxstring_abis(expand_microarchitectures(supported_platforms(), ["x86_64", "avx", "avx2"]); skip=!Sys.iswindows)

augment_platform_block = """
    $(MicroArchitectures.augment)

    function augment_platform!(platform::Platform)
        # We augment only x86_64
        @static if Sys.ARCH === :x86_64
            augment_microarchitecture!(platform)
        else
            platform
        end
    end
    """

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a")),
    # For OpenMP we use libomp from `LLVMOpenMP_jll` where we use LLVM as compiler (BSD
    # systems), and libgomp from `CompilerSupportLibraries_jll` everywhere else.
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
]

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

max_stars_repo_name: onetonfoot/Yggdrasil

max_issues_repo_name: onetonfoot/Yggdrasil

max_forks_repo_name: onetonfoot/Yggdrasil


text:
using TaskGraphs

## set up the environment
vtx_grid = initialize_dense_vtx_grid(4,4) # 4 x 4 grid world
#  1   2   3   4
#  5   6   7   8
#  9  10  11  12
# 13  14  15  16
env = construct_factory_env_from_vtx_grid(vtx_grid)

## Define the initial conditions of the robots
robot_ics = [
    ROBOT_AT(1,2), # robot 1 starts at vertex 2
    ROBOT_AT(2,9), # robot 2 starts at vertex 9
]

## Define the manufacturing project
spec = ProjectSpec()
# set initial conditions of "raw materials"
set_initial_condition!(spec,OBJECT_AT(1,4)) # object 1 starts at vertex 4
set_initial_condition!(spec,OBJECT_AT(2,16))  # object 2 starts at vertex 16
# define the operations that need to take place
op1 = Operation(Δt=2) # operation 1 has a duration of 2 time steps
# inputs
set_precondition!(op1,OBJECT_AT(1,8)) # object 1 must be at vertex 8 before op1 can begin
set_precondition!(op1,OBJECT_AT(2,12)) # object 2 must be at vertex 12 before op1 can begin
# outputs
set_postcondition!(op1,OBJECT_AT(3,7)) # object 3 appears at vertex 7 when op1 is completed
add_operation!(spec,op1)
# add a terminal operation with no outputs
op2 = Operation(Δt=0)
set_precondition!(op2,OBJECT_AT(3,13))
add_operation!(spec,op2)

## define solver
solver = NBSSolver()
# finalize problem construction (the solver is passed as an argument here
# because it determines what cost model is used for the problem)
prob = pctapf_problem(solver,spec,env,robot_ics)
# solve the problem
solution, cost = solve!(solver,prob)
# check if the problem was solved to optimality
@show feasible_status(solver)
@show optimal_status(solver)

ext: jl


text:
# ---
# title: 304. Range Sum Query 2D - Immutable
# id: problem304
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/range-sum-query-2d-immutable/description/>
# hidden: true
# ---
#
# Given a 2D matrix _matrix_ , find the sum of the elements inside the rectangle
# defined by its upper left corner ( _row_ 1, _col_ 1) and lower right corner (
# _row_ 2, _col_ 2).
#
# ![Range Sum Query 2D](https://assets.leetcode.com/static/images/courses/range_sum_query_2d.png)
# The above rectangle (with the red border) is defined by (row1, col1) = **(2,
# 1)** and (row2, col2) = **(4, 3)** , which contains sum = **8**.
#
# **Example:**
#
#
#
#     Given matrix = [
#       [3, 0, 1, 4, 2],
#       [5, 6, 3, 2, 1],
#       [1, 2, 0, 1, 5],
#       [4, 1, 0, 1, 7],
#       [1, 0, 3, 0, 5]
#     ]
#
#     sumRegion(2, 1, 4, 3) -> 8
#     sumRegion(1, 1, 2, 2) -> 11
#     sumRegion(1, 2, 2, 4) -> 12
#
#
# **Note:**
#
#   1. You may assume that the matrix does not change.
#   2. There are many calls to _sumRegion_ function.
#   3. You may assume that _row_ 1 ≤ _row_ 2 and _col_ 1 ≤ _col_ 2.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#!/usr/bin/env julia

using Test

include("day03.jl")
using .Wires: parse_wire, parse_direction, gen_segment, problem1, problem2

@testset "parse_direction" begin
  for direction = ["U0", "L0", "R0", "D0"]
    @test (0,0) == parse_direction(direction)
  end

  @test (1,0) == parse_direction("R1")
  @test (0,1) == parse_direction("U1")
  @test (-1,0) == parse_direction("L1")
  @test (0,-1) == parse_direction("D1")

  @test (45,0) == parse_direction("R45")
  @test (0,10) == parse_direction("U10")
  @test (-3,0) == parse_direction("L3")
  @test (0,-6) == parse_direction("D6")
end

@testset "gen_segment" begin
  origin = (0, 0)
  @test Dict([]) == gen_segment(origin, origin, 1)

  for direction = [(1,0), (0,1), (-1,0), (0,-1)]
    @test Dict(direction => 1) == gen_segment(origin, direction, 1)
  end

  for y = [0, -1, 1, 5]
    @test Dict((1,y)=>1, (2,y)=>2, (3,y)=>3) == gen_segment((0,y), (3,0), 1)
    @test Dict((-1,y)=>1, (-2,y)=>2, (-3,y)=>3) == gen_segment((0,y), (-3,0), 1)
  end

  for x = [0, -1, 1, 5]
    @test Dict((x,1)=>1, (x,2)=>2, (x,3)=>3) == gen_segment((x,0), (0,3), 1)
    @test Dict((x,-1)=>1, (x,-2)=>2, (x,-3)=>3) == gen_segment((x,0), (0,-3), 1)
  end

  @test Dict((1,1)=>1) == gen_segment((0,1), (1,0), 1)
  @test Dict((1,1)=>1) == gen_segment((1,0), (0,1), 1)
end

@testset "parse_wire" begin
  wire = parse_wire("U7")
  @test Dict((0,1)=>1, (0,2)=>2, (0,3)=>3, (0,4)=>4, (0,5)=>5, (0,6)=>6, (0,7)=>7) == wire

  wire = parse_wire("U1,R1")
  @test Dict((0,1)=>1, (1,1)=>2) == wire

  wire = parse_wire("R1,U1")
  @test Dict((1,0)=>1, (1,1)=>2) == wire
end

@testset "problem1" begin
  @test 6 == problem1("R8,U5,L5,D3", "U7,R6,D4,L4")
  @test 159 == problem1("R75,D30,R83,U83,L12,D49,R71,U7,L72",
                        "U62,R66,U55,R34,D71,R55,D58,R83")
  @test 135 == problem1("R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51",
                        "U98,R91,D20,R16,D67,R40,U7,R15,U6,R7")
end

@testset "problem2" begin
  @test 30 == problem2("R8,U5,L5,D3", "U7,R6,D4,L4")
  @test 610 == problem2("R75,D30,R83,U83,L12,D49,R71,U7,L72",
                        "U62,R66,U55,R34,D71,R55,D58,R83")
  @test 410 == problem2("R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51",
                        "U98,R91,D20,R16,D67,R40,U7,R15,U6,R7")
end

max_stars_repo_name: bfontaine/advent-of-code

max_issues_repo_name: bfontaine/advent-of-code

max_forks_repo_name: bfontaine/advent-of-code


text:
include("entities/AbstractDataStore.jl")
include("entities/AbstractBigDataEntries.jl")
include("entities/InMemoryDataStore.jl")
include("entities/FileDataStore.jl")

include("services/AbstractDataStore.jl")
include("services/AbstractBigDataEntries.jl")
include("services/InMemoryDataStore.jl")
include("services/FileDataStore.jl")

export AbstractDataStore

export AbstractBigDataEntry, GeneralBigDataEntry, MongodbBigDataEntry, FileBigDataEntry
export InMemoryDataStore, FileDataStore

export getBigData, addBigData!, updateBigData!, deleteBigData!, listStoreEntries
export copyStore


text:
using SnoopCompile

SnoopCompile.@snoop "glp_compiles.csv" begin
    using GLPlot;GLPlot.init()
    using Colors, GeometryTypes
    glplot(rand(Float32, 32,32))
    glplot(rand(Float32, 32,32), :surface)
    glplot(rand(Point3f0,32))
    glplot(rand(Point3f0,32), :lines)
    glplot(rand(Point2f0,32), :lines)
    glplot(rand(Point2f0,32))
    glplot(RGBA{Float32}[RGBA{Float32}(rand(), rand(), rand(), rand()) for i=1:512, j=1:512])
    while isopen(GLPlot.viewing_screen)
        yield()
    end
end

using GLPlot
data = SnoopCompile.read("glp_compiles.csv")
blacklist = ["MIME"]
pc = SnoopCompile.format_userimg(data[end:-1:1,2], blacklist=blacklist)
SnoopCompile.write(Pkg.dir("GLPlot", "src", "glp_userimg.jl"), pc)


text:
import RRStrace
using Cxx
sess = icxx"""rr::ReplaySession::create("");"""
RRStrace.strace(sess)


text:
"""
    GranulatedAggregation

Functor object for using partial aggregation in an L-shaped algorithm. Create by supplying a [`GranulatedAggregate`](@ref) object through `aggregate` in `LShaped.Optimizer` or by setting the [`Aggregator`](@ref) attribute.

...
# Parameters
- `size::Int`: Number of cuts in each aggregate
...
"""
struct GranulatedAggregation{T <: AbstractFloat, Agg <: AbstractAggregation} <: AbstractAggregation
    size::Int
    num_subproblems::Int
    start_id::Int
    collection::CutCollection{T}
    inner::Agg

    function GranulatedAggregation(size::Integer, start_id::Integer, num_subproblems::Integer, aggregator::AbstractAggregation, ::Type{T}) where T <: AbstractFloat
        Agg = typeof(aggregator)
        return new{T,Agg}(min(size, num_subproblems), num_subproblems, start_id, CutCollection(T, start_id), aggregator)
    end
end

function aggregate_cut!(lshaped::AbstractLShaped, aggregation::GranulatedAggregation, cut::HyperPlane)
    added = passthrough!(lshaped, aggregation, cut)
    add_to_collection!(aggregation.collection, cut, lshaped.x)
    if considered(aggregation.collection) == aggregation.size
        if collection_size(aggregation.collection) == aggregation.size
            granulated_cut = aggregate(aggregation.collection)
            added |= aggregate_cut!(lshaped, aggregation.inner, granulated_cut)
        end
        renew!(aggregation.collection, aggregation.collection.id % num_thetas(lshaped) + 1)
    end
    return added
end

function passthrough!(lshaped::AbstractLShaped, aggregation::GranulatedAggregation, cut::HyperPlane{H}) where H <: HyperPlaneType
    return add_cut!(lshaped, HyperPlane(cut.δQ, cut.q, aggregation.collection.id, H))
end

function passthrough!(lshaped::AbstractLShaped, aggregation::GranulatedAggregation, cut::HyperPlane{OptimalityCut})
    return false
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::GranulatedAggregation, metadata::MetaData, t::Integer, cut::HyperPlane, x::AbstractArray)
    passthrough!(cutqueue, aggregation, cut, t, x)
    add_to_collection!(aggregation.collection, cut, x)
    if considered(aggregation.collection) == aggregation.size
        if collection_size(aggregation.collection) == aggregation.size
            granulated_cut = aggregate(aggregation.collection)
            aggregate_cut!(cutqueue, aggregation.inner, metadata, t, granulated_cut, x)
        end
        new_id = aggregation.start_id + (aggregation.collection.id - aggregation.start_id + 1) % ceil(Int, aggregation.num_subproblems/aggregation.size)
        renew!(aggregation.collection, new_id)
    end
    return nothing
end

function passthrough!(cutqueue::CutQueue, aggregation::GranulatedAggregation, cut::HyperPlane{H}, t::Integer, x::AbstractVector) where H <: HyperPlaneType
    put!(cutqueue, (t, HyperPlane(cut.δQ, cut.q, aggregation.collection.id, H)))
    return nothing
end

function passthrough!(::CutQueue, ::GranulatedAggregation, ::HyperPlane{OptimalityCut}, ::Integer, ::AbstractVector)
    return nothing
end

function num_thetas(num_subproblems::Integer, aggregation::GranulatedAggregation)
    return ceil(Int, num_subproblems / aggregation.size)
end

function num_thetas(num_subproblems::Integer, aggregation::GranulatedAggregation, sp::StochasticPrograms.ScenarioProblems)
    jobsize = ceil(Int, num_subproblems / nworkers())
    n = ceil(Int, jobsize/aggregation.size)
    remainder = num_subproblems - (nworkers() - 1) * jobsize
    aggregationem = ceil(Int, remainder/aggregation.size)
    return n * (nworkers() - 1) + aggregationem
end

function num_thetas(::Integer, aggregation::GranulatedAggregation, sp::DistributedScenarioProblems)
    return sum([ceil(Int, nscen/aggregation.size) for nscen in sp.scenario_distribution])
end

function flush!(lshaped::AbstractLShaped, aggregation::GranulatedAggregation)
    added = false
    if collection_size(aggregation.collection) > 0 && collection_size(aggregation.collection) == considered(aggregation.collection)
        granulated_cut = aggregate(aggregation.collection)
        added |= aggregate_cut!(lshaped, aggregation.inner, granulated_cut)
    end
    renew!(aggregation.collection, aggregation.start_id)
    added |= flush!(lshaped, aggregation.inner)
    return added
end

function flush!(cutqueue::CutQueue, aggregation::GranulatedAggregation, metadata::MetaData, t::Integer, x::AbstractArray)
    if collection_size(aggregation.collection) > 0 && collection_size(aggregation.collection) == considered(aggregation.collection)
        granulated_cut = aggregate(aggregation.collection)
        aggregate_cut!(cutqueue, aggregation.inner, metadata, t, granulated_cut, x)
    end
    renew!(aggregation.collection, aggregation.start_id)
    flush!(cutqueue, aggregation.inner, metadata, t, x)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    GranulatedAggregate

Factory object for [`GranulatedAggregation`](@ref). Pass to `aggregate` in `LShaped.Optimizer` or by setting the [`Aggregator`](@ref) attribute.  See ?GranulatedAggregation for parameter descriptions.

"""
mutable struct GranulatedAggregate <: AbstractAggregator
    size::Int
    start_id::Int
    inner::AbstractAggregator

    function GranulatedAggregate(size::Integer, aggregator::AbstractAggregator)
        return new(size, 1, aggregator)
    end

    function GranulatedAggregate(size::Integer, start_id::Int, aggregator::AbstractAggregator)
        return new(size, start_id, aggregator)
    end
end

function (aggregator::GranulatedAggregate)(num_subproblems::Integer, ::Type{T}) where T <: AbstractFloat
    aggregator.size == 1 && return NoAggregation()
    aggregator.size == num_subproblems && return FullAggregation(1, num_subproblems, T)
    inner = aggregator.inner(num_subproblems, T)
    return GranulatedAggregation(aggregator.size, aggregator.start_id, num_subproblems, inner, T)
end

function remote_aggregator(aggregation::GranulatedAggregation, sp::ScenarioProblems, w::Integer)
    (nscen, extra) = divrem(num_subproblems(sp), nworkers())
    prev = map(2:(w - 1)) do p
        jobsize = nscen + (extra + 2 - p > 0)
        ceil(Int, jobsize / aggregation.size)
    end
    start_id = isempty(prev) ? 1 : sum(prev) + 1
    return GranulatedAggregate(aggregation.size, start_id, remote_aggregator(aggregation.inner, sp, w))
end

function remote_aggregator(aggregation::GranulatedAggregation, sp::DistributedScenarioProblems, w::Integer)
    prev = map(2:(w - 1)) do p
        ceil(Int, sp.scenario_distribution[p-1] / aggregation.size)
    end
    start_id = isempty(prev) ? 1 : sum(prev) + 1
    nscen = sp.scenario_distribution[w-1]
    return GranulatedAggregate(aggregation.size, start_id, remote_aggregator(aggregation.inner, sp, w))
end

function str(aggregator::GranulatedAggregate)
    return "Granulated cut aggregation of size $(aggregator.size) with inner aggregator: $(str(aggregator.inner))"
end


text:
import Theta: ellipsoid_pointwise, ellipsoid_uniform

@testset "Ellipsoid" begin
    T = [1 0; 0 1];
    shift = [0; 0];
    @test Set(ellipsoid_pointwise(T, 2, shift)) == Set([[0,-1], [-1,0], [1,0], [0,1], [0,0]])
    @test issubset(ellipsoid_pointwise(T, 2, shift), ellipsoid_uniform(T, 2))
end


text:
using JpegTurbo
using Documenter

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true"
)

makedocs(;
    modules=[JpegTurbo],
    sitename="JpegTurbo.jl",
    format=format,
    pages=[
        "Home" => "index.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/johnnychen94/JpegTurbo.jl",
    devbranch="master",
)


text:
struct S a; b end
struct S c; d end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function tmul_threads!(C, A, B)
    @tullio C[m,n] = A[m,k] * B[k,n]
end
function tmul_no_threads!(C, A, B)
    @tullio C[m,n] = A[m,k] * B[k,n] threads=false
end
function lvmul_threads!(C, A, B)
    @avxt for n ∈ indices((C,B), 2), m ∈ indices((C,A), 1)
        Cmn = zero(eltype(C))
        for k ∈ indices((A,B), (2,1))
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
end
function lvmul_no_threads!(C, A, B)
    @avx for n ∈ indices((C,B), 2), m ∈ indices((C,A), 1)
        Cmn = zero(eltype(C))
        for k ∈ indices((A,B), (2,1))
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
end

function generic_matmul!(C, A, B)
    istransposed(C) === 'N' || (generic_matmul!(untransposed(C), _transpose(B), _transpose(A)); return C)
    transA = istransposed(A)
    transB = istransposed(B)
    pA     = untransposed(A);
    pB     = untransposed(B)
    LinearAlgebra.generic_matmatmul!(C, transA, transB, pA, pB)
end

function getfuncs(libs::Vector{Symbol}, threaded::Bool)::Vector{Function}
    map(libs) do i
        if i === :MKL
            gemmmkl!
        elseif i === :MKL_DIRECT || i === :MKL_direct
            gemmmkl_direct!
        elseif i === :OpenBLAS
            gemmopenblas!
        elseif i === :BLIS || i === :blis
            gemmblis!
        elseif i === :Octavian
            threaded ? matmul! : matmul_serial!
        elseif i === :Tullio
            threaded ? tmul_threads! : tmul_no_threads!
        elseif i === :Gaius
            threaded ? Gaius.mul! : Gaius.mul_serial!
        elseif i === :LoopVectorization
            threaded ? lvmul_threads! : lvmul_no_threads!
        elseif i === :generic || i === :Generic || i === :GENERIC
            generic_matmul!
        else
            throw("Library $i not reognized.")
        end
    end
end


text:
__precompile__()

module HemirealNumbers

import Base: +, -, *, /, \, ^, abs, abs2, conj, convert, isfinite, promote_op, promote_rule, real, show, zero
import Base.LinAlg: A_mul_Bt, A_mul_Bc, At_mul_B, Ac_mul_B, At_mul_Bt, Ac_mul_Bc

export PureHemi, Hemireal, μ, ν, mu, nu

immutable PureHemi{T<:Real} <: Number
    m::T
    n::T
end
PureHemi(m::Real, n::Real) = PureHemi(promote(m,n)...)

immutable Hemireal{T<:Real} <: Number
    r::T
    h::PureHemi{T}
end

Hemireal(r::Real, m::Real, n::Real) = Hemireal(r, PureHemi(m, n))
function Hemireal{R,H}(r::R, h::PureHemi{H})
    T = promote_type(R,H)
    Hemireal{T}(r, PureHemi{T}(h.m, h.n))
end
Hemireal{R}(h::PureHemi{R}) = Hemireal(zero(R), h)
Hemireal{R}(r::R) = Hemireal(r, zero(PureHemi{R}))

const μ = PureHemi(true,false)
const ν = PureHemi(false,true)

## PureHemi implementation
convert{R}(::Type{PureHemi{R}}, x::PureHemi) = PureHemi{R}(x.m, x.n)
convert{R}(::Type{PureHemi{R}}, x::Real) = x == 0 ? zero(PureHemi{R}) : throw(DomainError())  # error("Non-zero reals cannot be converted to pure-hemi numbers")

(-)(x::PureHemi) = PureHemi(-x.m, -x.n)

(+)(x::PureHemi, y::PureHemi) = PureHemi(x.m+y.m, x.n+y.n)
(-)(x::PureHemi, y::PureHemi) = PureHemi(x.m-y.m, x.n-y.n)
(*)(x::PureHemi, y::PureHemi) = x.m*y.n + x.n*y.m
(*)(c::Bool, x::PureHemi)     = PureHemi(c*x.m, c*x.n)
(*)(c::Real, x::PureHemi)     = PureHemi(c*x.m, c*x.n)
(*)(x::PureHemi, c::Real)     = c*x
(/)(x::PureHemi, c::Real)     = PureHemi(x.m/c, x.n/c)
(\)(c::Real, x::PureHemi)     = x/c
function (/)(x::PureHemi, y::PureHemi)
    @assert x.m*y.n == x.n*y.m
    y.n == 0 ? x.m/y.m : x.n/y.n
end
(\)(x::PureHemi, y::PureHemi) = y/x
(^)(x::PureHemi, p::Integer)  = (x*x)^(p/2)  # note x*x is much faster than x^2
(^)(x::PureHemi, p::Rational) = (x*x)^(p/2)
(^)(x::PureHemi, p::Real)     = (x*x)^(p/2)

# Symmetric division (there are other solutions y to y*x = c)
(/)(c::Real, x::PureHemi) = PureHemi(c/(2*x.n), c/(2*x.m))
(\)(x::PureHemi, c::Real) = c/x

real{R}(::Type{PureHemi{R}}) = R
real{R}(x::PureHemi{R}) = zero(R)
conj(x::PureHemi) = x
mu(x::PureHemi) = x.m
nu(x::PureHemi) = x.n

zero{R}(::Type{PureHemi{R}}) = PureHemi{R}(0, 0)
zero{R}(::PureHemi{R}) = zero(PureHemi{R})

isfinite(x::PureHemi) = isfinite(x.m) && isfinite(x.n)
# I'm not entirely certain whether the next two should even be defined.
# (note that abs2(x) != x*x)
abs2(x::PureHemi) = x.m*x.m + x.n*x.n
abs(x::PureHemi) = sqrt(abs2(x))

## Hemireal implementation
convert{R}(::Type{Hemireal{R}}, x::Hemireal) = Hemireal{R}(x.r, x.h)
convert{R}(::Type{Hemireal{R}}, x::PureHemi) = Hemireal{R}(0, x)
convert{R}(::Type{Hemireal{R}}, x::Real) = Hemireal{R}(x, zero(PureHemi{R}))
convert{R}(::Type{Hemireal{R}}, r::Real, m::Real, n::Real) = Hemireal{R}(r, PureHemi{R}(m, n))

(-)(x::Hemireal) = Hemireal(-x.r, -x.h)

(+)(x::Hemireal, y::Hemireal) = Hemireal(x.r+y.r, x.h+y.h)
(-)(x::Hemireal, y::Hemireal) = Hemireal(x.r-y.r, x.h-y.h)
(*)(x::Hemireal, y::Hemireal) = Hemireal(x.r*y.r + x.h*y.h, x.r*y.h + x.h*y.r)

(+)(x::Hemireal, y::PureHemi) = x+Hemireal(y)
(+)(x::PureHemi, y::Hemireal) = y+x
(-)(x::Hemireal, y::PureHemi) = x-Hemireal(y)
(-)(x::PureHemi, y::Hemireal) = Hemireal(x)-y
(*)(x::Hemireal, y::PureHemi) = x*Hemireal(y)
(*)(x::PureHemi, y::Hemireal) = y*x

(+)(x::Hemireal, y::Real) = x+Hemireal(y)
(+)(x::Real, y::Hemireal) = y+x
(-)(x::Hemireal, y::Real) = x-Hemireal(y)
(-)(x::Real, y::Hemireal) = Hemireal(x)-y
(*)(c::Bool, x::Hemireal) = Hemireal(c*x.r, c*x.h)
(*)(c::Real, x::Hemireal) = Hemireal(c*x.r, c*x.h)
(*)(x::Hemireal, c::Real) = c*x
(/)(x::Hemireal, c::Real) = Hemireal(x.r/c, x.h/c)
(\)(c::Real, x::Hemireal) = x/c

real{R}(::Type{Hemireal{R}}) = R
real{R}(x::Hemireal{R}) = x.r
conj(x::Hemireal) = x
mu(x::Hemireal) = mu(x.h)
nu(x::Hemireal) = nu(x.h)

zero{R}(::Type{Hemireal{R}}) = Hemireal{R}(0, PureHemi{R}(0, 0))
zero{R}(::Hemireal{R}) = zero(Hemireal{R})

isfinite(x::Hemireal) = isfinite(x.r) && isfinite(x.h)
# ?
abs2(x::Hemireal) = x.r*x.r + abs2(x.h)
abs(x::Hemireal) = sqrt(abs2(x))

for f in (:mu, :nu)
    @eval begin
        function $f(X::AbstractArray)
            out = Array(real(eltype(X)), size(X))
            for I in eachindex(X)
                out[I] = $f(X[I])
            end
            out
        end
    end
end

(*){H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVector{H2}) = A_mul_B!(Array(promote_type(real(H1),real(H2)), size(A,1)), A, B)
(*){H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedMatrix{H2}) = A_mul_B!(Array(promote_type(real(H1),real(H2)), size(A,1), size(B,2)), A, B)
A_mul_Bt{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVecOrMat{H2}) = A_mul_Bt!(Array(promote_type(real(H1),real(H2)), size(A,1), size(B,1)), A, B)
A_mul_Bc{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVecOrMat{H2}) = A_mul_Bt!(Array(promote_type(real(H1),real(H2)), size(A,1), size(B,1)), A, B)
At_mul_B{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVector{H2}) = At_mul_B!(Array(promote_type(real(H1),real(H2)), size(A,2)), A, B)
At_mul_B{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedMatrix{H2}) = At_mul_B!(Array(promote_type(real(H1),real(H2)), size(A,2), size(B,2)), A, B)
Ac_mul_B{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVector{H2}) = At_mul_B!(Array(promote_type(real(H1),real(H2)), size(A,2)), A, B)
Ac_mul_B{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedMatrix{H2}) = At_mul_B!(Array(promote_type(real(H1),real(H2)), size(A,2), size(B,2)), A, B)
At_mul_Bt{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVecOrMat{H2}) = At_mul_Bt!(Array(promote_type(real(H1),real(H2)), size(A,2), size(B,1)), A, B)
Ac_mul_Bc{H1<:PureHemi,H2<:PureHemi}(A::StridedVecOrMat{H1}, B::StridedVecOrMat{H2}) = At_mul_Bt!(Array(promote_type(real(H1),real(H2)), size(A,2), size(B,1)), A, B)

promote_op{H1,H2}(::Base.MulFun, ::Type{PureHemi{H1}}, ::Type{PureHemi{H2}}) = promote_type(H1,H2)
promote_op{H1,H2}(::Base.DotMulFun, ::Type{PureHemi{H1}}, ::Type{PureHemi{H2}}) = promote_type(H1,H2)
promote_op{R<:Real,H}(::Base.MulFun, ::Type{R}, ::Type{PureHemi{H}}) = PureHemi{promote_type(R,H)}
promote_op{R<:Real,H}(::Base.DotMulFun, ::Type{R}, ::Type{PureHemi{H}}) = PureHemi{promote_type(R,H)}
promote_op{R<:Real,H}(::Base.MulFun, ::Type{PureHemi{H}}, ::Type{R}) = promote_op(Base.MulFun(), R, PureHemi{H})
promote_op{R<:Real,H}(::Base.DotMulFun, ::Type{PureHemi{H}}, ::Type{R}) = promote_op(Base.MulFun(), R, PureHemi{H})

promote_rule{H}(::Type{Bool}, ::Type{PureHemi{H}}) = Hemireal{promote_type(Bool,H)}
promote_rule{R,H}(::Type{Irrational{R}}, ::Type{PureHemi{H}}) = Hemireal{promote_type(Irrational{R},H)}
promote_rule{R<:Real,H}(::Type{R}, ::Type{PureHemi{H}}) = Hemireal{promote_type(R,H)}
promote_rule{H1,H2}(::Type{PureHemi{H1}}, ::Type{PureHemi{H2}}) = PureHemi{promote_type(H1,H2)}
promote_rule{H}(::Type{Bool}, ::Type{Hemireal{H}}) = Hemireal{promote_type(Bool,H)}
promote_rule{R,H}(::Type{Irrational{R}}, ::Type{Hemireal{H}}) = Hemireal{promote_type(Irrational{R},H)}
promote_rule{R<:Real,H}(::Type{R}, ::Type{Hemireal{H}}) = Hemireal{promote_type(R,H)}
promote_rule{H1,H2}(::Type{PureHemi{H1}}, ::Type{Hemireal{H2}}) = Hemireal{promote_type(H1,H2)}
promote_rule{H1,H2}(::Type{Hemireal{H1}}, ::Type{Hemireal{H2}}) = Hemireal{promote_type(H1,H2)}

show(io::IO, x::PureHemi) = x.n >= 0 ? print(io, x.m, "μ + ", x.n, 'ν') : print(io, x.m, "μ - ", abs(x.n), 'ν')
show(io::IO, x::Hemireal) = x.h.m >= 0 ? print(io, x.r, " + ", x.h) : print(io, x.r, " - ", PureHemi(-x.h.m, x.h.n))

end # module


text:
using Test: include
using HighDimPDE, Test

@testset "HighDimPDE" begin
    include("reflect.jl")
    include("MLP.jl")
    include("DeepSplitting.jl")
    include("MCSample.jl")
end


text:
using GLVisualize

import GLVisualize: get_docs, all_docs

# prints the documenation of parameters for visualizing a float matrix
get_docs(rand(Float32, 32, 32))
# prints the documenation of parameters for visualizing a float matrix in surface style
get_docs(rand(Float32, 32, 32), :surface)

# prints all visualization methods available with some documenation if available
io = IOBuffer()
all_docs(io)
String(take!(io)) # TODO test it returns the correct string!


text:
abstract type Objective end

mutable struct TrackingObjective{Q,U,C,B} <: Objective
    q::Vector{Q}
    u::Vector{U}
    γ::Vector{C}
    b::Vector{B}
end

function TrackingObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.nq})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.nu})) for t = 1:H],
    γ = [Diagonal(zeros(SizedVector{model.nc})) for t = 1:H],
    b = [Diagonal(zeros(SizedVector{model.nc * friction_dim(env)})) for t = 1:H])
    return TrackingObjective(q, u, γ, b)
end

mutable struct TrackingVelocityObjective{Q,V,U,C,B,VT} <: Objective
    q::Vector{Q}
    v::Vector{V}
    u::Vector{U}
    γ::Vector{C}
    b::Vector{B}
    v_target::Vector{VT}
end

function TrackingVelocityObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.nq})) for t = 1:H],
    v = [Diagonal(zeros(SizedVector{model.nq})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.nu})) for t = 1:H],
    γ = [Diagonal(zeros(SizedVector{model.nc})) for t = 1:H],
    b = [Diagonal(zeros(SizedVector{model.nc * friction_dim(env)})) for t = 1:H],
    v_target = [zeros(SizedVector{model.nq}) for t = 1:H])
    return TrackingVelocityObjective(q, v, u, γ, b, v_target)
end


text:
################################
### Aerodynamic conductance ####
################################

#' Aerodynamic Conductance
#'
#' Bulk aerodynamic conductance, including options for the boundary layer conductance
#'              formulation and stability correction functions.
#'
#' - data              DataFrame or matrix containing all required variables
#' - Tair              Air temperature (deg C)
#' - pressure          Atmospheric pressure (kPa)
#' - wind              Wind speed (m s-1)
#' - ustar             Friction velocity (m s-1)
#' - H                 Sensible heat flux (W m-2)
#' - zr                Instrument (reference) height (m)
#' - zh                Canopy height (m)
#' - d                 Zero-plane displacement height (m)
#' - z0m               Roughness length for momentum (m), optional; if not provided, it is estimated from `roughness_parameters`
#'                          (method=Val(:wind_profile)). Only used if `wind_profile = true` and/or `Gb_model` = `Val(:Su_2001)` or
#'                          `Val(:Choudhury_1988)`.
#' - Dl                Characteristic leaf dimension (m) (if `Gb_model` = `Val(:Su_2001)`)
#'                          or leaf width (if `Gb_model` = `Val(:Choudhury_1988)`); ignored otherwise.
#' - N                 Number of leaf sides participating in heat exchange (1 or 2); only used if `Gb_model = Val(:Su_2001)`.
#'                          Defaults to 2.
#' - fc                Fractional vegetation cover (-); only used if `Gb_model = Val(:Su_2001)`. See Details.
#' - LAI               One-sided leaf area index (m2 m-2); only used if `Gb_model` = `Val(:Choudhury_1988)` or `Val(:Su_2001)`.
#' - Cd                Foliage drag coefficient (-); only used if `Gb_model = Val(:Su_2001)`.
#' - hs                Roughness length of bare soil (m); only used if `Gb_model = Val(:Su_2001)`.
#' - wind_profile      Should Ga for momentum be calculated based on the logarithmic wind profile equation?
#'                          Defaults to `false`.
#' - stab_correction   Should stability correction be applied? Defaults to `true`. Ignored if `wind_profile = false`.
#' - stab_formulation  Stability correction function. Either `Val(:Dyer_1970)` (default) or
#'                          `Val(:Businger_1971)`. Ignored if `wind_profile = false` or if `stab_correction = false`.
#' - Gb_model          Boundary layer resistance formulation. One of `Val(:Thom_1972),Val(:Choudhury_1988),Val(:Su_2001),Val(:constant_kB1)`.
#' - kB_h              kB-1 value for heat transfer; only used if `Gb_model = Val(:constant_kB1)`
#' - Sc                Optional: Schmidt number of additional quantities to be calculated
#' - Sc_name           Optional: Name of the additonal quantities, has to be of same length than
#'                          `Sc_name`
#' - constants         k - von Karman constant
#'                          cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                          Kelvin - conversion degree Celsius to Kelvin
#'                          g - gravitational acceleration (m s-2)
#'                          pressure0 - reference atmospheric pressure at sea level (Pa)
#'                          Tair0 - reference air temperature (K)
#'                          Sc_CO2 - Schmidt number for CO2
#'                          Pr - Prandtl number (if `Sc` is provided)
#'
#' # Details

#'
#' Aerodynamic conductance for heat (Ga_h) is calculated as:
#'
#'     ``Ga_h = 1 / (Ra_m + Rb_h)``
#'
#'  where Ra_m is the aerodynamic resistance for momentum and Rb the (quasi-laminar) canopy boundary
#'  layer resistance ('excess resistance').
#'
#'  The aerodynamic resistance for momentum Ra_m is given by:
#'
#'     ``Ra_m = u/ustar^2``
#'
#'  Note that this formulation accounts for changes in atmospheric stability, and does not require an
#'  additional stability correction function.
#'
#'  An alternative method to calculate Ra_m is provided
#'  (calculated if `wind_profile = true`):
#'
#'     ``Ra_m = (ln((zr - d)/z0m) - psi_h) / (k ustar)``
#'
#'  If the roughness parameters z0m and d are unknown, they can be estimated using
#'  [`roughness_parameters`](@ref). The argument `stab_formulation` determines the stability correction function used
#'  to account for the effect of atmospheric stability on Ra_m (Ra_m is lower for unstable
#'  and higher for stable stratification). Stratification is based on a stability parameter zeta (z-d/L),
#'  where z = reference height, d the zero-plane displacement height, and L the Monin-Obukhov length,
#'  calculated with [`Monin_Obukhov_length`](@ref)
#'  The stability correction function is chosen by the argument `stab_formulation`. Options are
#'  `Val(:Dyer_1970)` and `Val(:Businger_1971)`.
#'
#'  The model used to determine the canopy boundary layer resistance for heat (Rb_h) is specified by
#'  the argument `Gb_model`. The following options are implemented:
#'  `Val(:Thom_1972)` is an empirical formulation based on the friction velocity (ustar) (Thom 1972):
#'
#'    ``Rb_h = 6.2ustar^-0.667``
#'
#'  The model by Choudhury & Monteith 1988 (`Gb_model = Val(:Choudhury_1988)`),
#'  calculates Rb_h based on leaf width, LAI and ustar (Note that function argument `Dl`
#'  represents leaf width (w) and not characteristic leaf dimension (Dl)
#'  if `Gb_model` = `Val(:Choudhury_1988)`):
#'
#'     ``Gb_h = LAI((0.02/\\alpha)*sqrt(u(zh)/w)*(1-exp(-\\alpha/2)))``
#'
#'  where ``\\alpha`` is a canopy attenuation coefficient modeled in dependence on LAI,
#'  u(zh) is wind speed at canopy height (calculated from [`wind_profile`](@ref)),
#'  and w is leaf width (m). See [`Gb_Choudhury`](@ref) for further details.
#'
#'  The option `Gb_model = Val(:Su_2001)` calculates Rb_h based on the physically-based Rb model by Su et al. 2001,
#'  a simplification of the model developed by Massman 1999:
#'
#'     ``kB-1 = (k Cd fc^2) / (4Ct ustar/u(zh)) + kBs-1(1 - fc)^2``
#'
#'  where Cd is a foliage drag coefficient (defaults to 0.2), fc is fractional
#'  vegetation cover, Bs-1 is the inverse Stanton number for bare soil surface,
#'  and Ct is a heat transfer coefficient. See [`Gb_Su`](@ref) for
#'  details on the model.
#'
#'
#'  The models calculate the parameter kB-1, which is related to Rb_h:
#'
#'     ``kB-1 = Rb_h * (k * ustar)``
#'
#'  Rb (and Gb) for water vapor and heat are assumed to be equal in this package.
#'  Gb for other quantities x is calculated as (Hicks et al. 1987):
#'
#'     ``Gb_x = Gb / (Sc_x / Pr)^0.67``
#'
#'  where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).
#'
#' # Value
 a dataframe with the following columns:
#'         - Ga_m: Aerodynamic conductance for momentum transfer (m s-1)
#'         - Ra_m: Aerodynamic resistance for momentum transfer (s m-1)
#'         - Ga_h: Aerodynamic conductance for heat transfer (m s-1)
#'         - Ra_h: Aerodynamic resistance for heat transfer (s m-1)
#'         - Gb_h: Canopy boundary layer conductance for heat transfer (m s-1)
#'         - Rb_h: Canopy boundary layer resistance for heat transfer (s m-1)
#'         - kB_h: kB-1 parameter for heat transfer
#'         - zeta: Stability parameter 'zeta' (missing if `wind_profile = false`)
#'         - psi_h: Integrated stability correction function (missing if `wind_profile = false`)
#'         - Ra_CO2: Aerodynamic resistance for CO2 transfer (s m-1)
#'         - Ga_CO2: Aerodynamic conductance for CO2 transfer (m s-1)
#'         - Gb_CO2: Canopy boundary layer conductance for CO2 transfer (m s-1)
#'         \item{Ga_Sc_name}{Aerodynamic conductance for `Sc_name` (m s-1). Only added if `Sc_name` and
#'                           the respective `Sc` are provided}
#'         \item{Gb_Sc_name}{Boundary layer conductance for `Sc_name` (m s-1). Only added if `Sc_name` and
#'                           the respective `Sc` are provided}
#'
#' # Note
#' The roughness length for water and heat (z0h) is not returned by this function, but
#'       it can be calculated from the following relationship (e.g. Verma 1989):
#'
#'       ``kB-1 = ln(z0m/z0h)``
#'
#'       it follows:
#'
#'       ``z0h = z0m / exp(kB-1)``
#'
#'       `kB-1` is an output of this function.
#'
#'       Input variables such as LAI, Dl, or zh can be either constants, or
#'       vary with time (i.e. vectors of the same length as `data`).
#'
#'       Note that boundary layer conductance to water vapor transfer (Gb_w) is often
#'       assumed to equal Gb_h. This assumption is also made in this R package, for
#'       example in the function [`surface_conductance`](@ref).
#'
#'       If the roughness length for momentum (`z0m`) is not provided as input, it is estimated
#'       from the function `roughness_parameters` within `wind_profile` if `wind_profile = true`
#'       and/or `Gb_model` = `Val(:Su_2001)` or `Val(:Choudhury_1988)` The `roughness_parameters`
#'       function estimates a single `z0m` value for the entire time period! If a varying `z0m` value
#'       (e.g. across seasons or years) is required, `z0m` should be provided as input argument.
#'
#'
#' #References
#' Verma, S., 1989: Aerodynamic resistances to transfers of heat, mass and momentum.
#'             In: Estimation of areal evapotranspiration, IAHS Pub, 177, 13-20.
#'
#'             Verhoef, A., De Bruin, H., Van Den Hurk, B., 1997: Some practical notes on the parameter kB-1
#'             for sparse vegetation. Journal of Applied Meteorology, 36, 560-572.
#'
#'             Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
#'             A preliminary multiple resistance routine for deriving dry deposition velocities
#'             from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
#'
#'             Monteith, J_L., Unsworth, M_H., 2008: Principles of environmental physics.
#'             Third Edition. Elsevier Academic Press, Burlington, USA.
#'
#' #See also
#' [`Gb_Thom`](@ref), [`Gb_Choudhury`](@ref), [`Gb_Su`](@ref) for calculations of Rb / Gb only
#'
#' ```@example; output = false
#' ```
#' df = DataFrame(Tair=25,pressure=100,wind=c(3,4,5),ustar=c(0.5,0.6,0.65),H=c(200,230,250))
#'
#' # simple calculation of Ga
#' aerodynamic_conductance(df,Gb_model=Val(:Thom_1972))
#'
#' # calculation of Ga using a model derived from the logarithmic wind profile
#' aerodynamic_conductance(df,Gb_model=Val(:Thom_1972),zr=40,zh=25,d=17.5,z0m=2,wind_profile=true)
#'
#' # simple calculation of Ga, but a physically based canopy boundary layer model
#' aerodynamic_conductance(df,Gb_model=Val(:Su_2001),zr=40,zh=25,d=17.5,Dl=0.05,N=2,fc=0.8)
#'
"""
"""
function aerodynamic_conductance(data,Tair="Tair",pressure="pressure",wind="wind",ustar="ustar",H="H",
                                    zr,zh,d,z0m=nothing,Dl,N=2,fc=nothing,LAI,Cd=0.2,hs=0.01,wind_profile=false,
                                    stab_correction=true,stab_formulation=c(Val(:Dyer_1970),Val(:Businger_1971)),
                                    Gb_model=c(Val(:Thom_1972),Val(:Choudhury_1988),Val(:Su_2001),Val(:constant_kB1)),
                                    kB_h=nothing,Sc=nothing,Sc_name=nothing,constants=BigleafConstants())

  Gb_model         = match_arg(Gb_model)
  stab_formulation = match_arg(stab_formulation)

  check_input(data,list(Tair,pressure,wind,ustar,H))

  ## calculate canopy boundary layer conductance (Gb)
  if (Gb_model in c(Val(:Thom_1972),Val(:Choudhury_1988),Val(:Su_2001)))

    if (Gb_model == Val(:Thom_1972))

      Gb_mod = Gb_Thom(ustar=ustar,Sc=Sc,Sc_name=Sc_name,constants=constants)

elseif (Gb_model == Val(:Choudhury_1988))

      Gb_mod = Gb_Choudhury(data,Tair=Tair,pressure=pressure,wind=wind,ustar=ustar,
                             H=H,leafwidth=Dl,LAI=LAI,zh=zh,zr=zr,d=d,z0m=z0m,
                             stab_formulation=stab_formulation,Sc=Sc,Sc_name=Sc_name,
                             constants=constants)

elseif (Gb_model == Val(:Su_2001))

      Gb_mod = Gb_Su(data=data,Tair=Tair,pressure=pressure,ustar=ustar,wind=wind,
                      H=H,zh=zh,zr=zr,d=d,z0m=z0m,Dl=Dl,N=N,fc=fc,LAI=LAI,Cd=Cd,hs=hs,
                      stab_formulation=stab_formulation,Sc=Sc,Sc_name=Sc_name,
                      constants=constants)

end

    kB_h = Gb_mod[,"kB_h"]
    Rb_h = Gb_mod[,"Rb_h"]
    Gb_h = Gb_mod[,"Gb_h"]
    Gb_x = DataFrame(Gb_mod[,grep(colnames(Gb_mod),pattern="Gb_")[-1]])
    colnames(Gb_x) = grep(colnames(Gb_mod),pattern="Gb_",value=true)[-1]

elseif (Gb_model == Val(:constant_kB1))

    if(isnothing(kB_h))
      stop("value of kB-1 has to be specified if Gb_model is set to 'constant_kB-1'!")
else
      Rb_h = kB_h/(constants.k * ustar)
      Gb_h = 1/Rb_h

      if (!isnothing(Sc) | !isnothing(Sc_name))
        if (length(Sc) != length(Sc_name))
          stop("arguments 'Sc' and 'Sc_name' must have the same length")
end
        if (!is_numeric(Sc))
          stop("argument 'Sc' must be numeric")
end
end

      Sc   = c(constants.Sc_CO2,Sc)
      Gb_x = DataFrame(lapply(Sc,function(x) Gb_h / (x/constants.Pr)^0.67))
      colnames(Gb_x) = paste0("Gb_",c("CO2",Sc_name))

end

end

  ## calculate aerodynamic conductance for momentum (Ga_m)
  if (wind_profile)

    if (isnothing(z0m) & Gb_model in c(Val(:constant_kB1),Val(:Thom_1972)))
      stop("z0m must be provided if wind_profile=true!")
elseif (isnothing(z0m) & Gb_model in c(Val(:Choudhury_1988),Val(:Su_2001)))
      # z0m estimated as in Choudhury_1988 or Su_2001
      z0m = roughness_parameters(method=Val(:wind_profile),zh=zh,zr=zr,d=d,data=data,
                                  Tair=Tair,pressure=pressure,wind=wind,ustar=ustar,H=H,
                                  stab_roughness=true,stab_formulation=stab_formulation,
                                  constants=constants)[,"z0m"]
end

    if (stab_correction)

      zeta  =  stability_parameter(data=data,Tair=Tair,pressure=pressure,ustar=ustar,
                                    H=H,zr=zr,d=d,constants=constants)

      if (stab_formulation in c(Val(:Dyer_1970),Val(:Businger_1971)))

        psi_h = stability_correction(zeta,formulation=stab_formulation)[,"psi_h"]

else
        stop("'stab_formulation' has to be one of 'Dyer_1970' or 'Businger_1971'.
             Choose 'stab_correction = false' if no stability correction should be applied.")
end

      Ra_m  = pmax((log((zr - d)/z0m) - psi_h),0) / (constants.k*ustar)

else

        Ra_m  = pmax((log((zr - d)/z0m)),0) / (constants.k*ustar)
        zeta = psi_h = rep(NA_integer_,length=length(Ra_m))

end

else

    if ((!missing(zr) | !missing(d) | !missing(z0m)) & Gb_model in c(Val(:constant_kB1),Val(:Thom_1972)))
      @warn"Provided roughness length parameters (zr,d,z0m) are not used if 'wind_profile = false' (the default). Ra_m is calculated as wind / ustar^2")
end

    Ra_m = wind / ustar^2
    zeta = psi_h = rep(NA_integer_,length=length(Ra_m))

end

  Ga_m   = 1/Ra_m
  Ra_h   = Ra_m + Rb_h
  Ga_h   = 1/Ra_h
  Ga_x   = 1/(Ra_m + 1/Gb_x)
  Ra_CO2 = 1/Ga_x[,1]
  colnames(Ga_x) = paste0("Ga_",c("CO2",Sc_name))

  Gab_x = cbind(Ga_x,Gb_x)
  Gab_x = Gab_x[rep(c(1,ncol(Gab_x)-(ncol(Gab_x)/2-1)),ncol(Gab_x)/2) + sort(rep(0:(ncol(Gab_x)/2-1),2))] # reorder columns

  return(DataFrame(Ga_m,Ra_m,Ga_h,Ra_h,Gb_h,Rb_h,kB_h,zeta,psi_h,Ra_CO2,Gab_x))

end


text:
using Pkg
Pkg.activate(".")

using DataFrames
using CSV
using TravelingSalesmanHeuristics

distance_data = CSV.read("distance_data.csv",header=[:i, :j, :dist])
N = length(unique(distance_data.i))
const distance_mx = zeros(Float64,(N,N))
for r in 1:nrow(distance_data)
   distance_mx[distance_data.i[r],distance_data.j[r]] = distance_data.dist[r]
end

sol = TravelingSalesmanHeuristics.solve_tsp(distance_mx,quality_factor =100)
sol
sbws_la = CSV.read("Subway_LV.csv")

using PyPlot
ids = sol[1]
PyPlot.cla()
PyPlot.ioff()
p = PyPlot.plot(sbws_la.long[ids], sbws_la.latt[ids])
PyPlot.display_figs()
PyPlot.show()


text:
### Apply the inverse of an HSS matrix to another HSS matrix
#
# as seen in
# Chandrasekaran, S., Gu, M., & Pals, T. (2006). A Fast $ULV$ Decomposition Solver for Hierarchically Semiseparable Representations.
# SIAM Journal on Matrix Analysis and Applications, 28(3), 603–622. https://doi.org/10.1137/S0895479803436652
#
# Massei, S., Robol, L., & Kressner, D. (n.d.). hm-toolbox: MATLAB SOFTWARE FOR HODLR AND HSS MATRICES.
# Retrieved from http://scg.ece.ucsb.edu/software.html
#
# Written by Boris Bonev, Feb. 2021

# function for acess that imitate the behavior in LinearAlgebra.jl
ldiv!(hssA::HssMatrix, hssB::HssMatrix) = _ldiv!(copy(hssA), hssB)
# lazy implementation of rdiv!
function rdiv!(hssA::HssMatrix, hssB::HssMatrix)
  hssA = adjoint(_ldiv!(adjoint(hssB), adjoint(hssA)))
end

## ULV divide algorithm to apply the inverse to another HSS matrix
# temporary name for function that actually just computes the ULV factorization
# the cluster structure in hssA and hssB should be compatible w/e that means...
function _ldiv!(hssA::HssMatrix, hssB::HssMatrix)
  if isleaf(hssA)
    D, U, V = _hssleaf(hssB)
    D = full(hssA) \ D
    return HssMatrix(D, U, V, isroot(hssB))
  else
    # bottom-up stage of the ULV solution algorithm
    hssL, QU, QL, QV, mk, nk, ktree  = _ulvfactor_leaves!(hssA, 0)
    hssB = _utransforms!(hssB, QU)
    hssQB = _extract_crows(hssB, nk)
    hssY0 = _ltransforms!(hssB, QL)

    # early exit if all remaining blocks are 0.
    if size(hssQB,1) == 0; return hssB = _vtransforms!(hssB, QV); end
    # TODO: this still fails if only one hss block gets fully eliminated - fix that!

    hssQB = hssQB - hssL * hssY0 # multiply triangularized part with solved part and substract from the
    hssQB = recompress!(hssQB)
    hssQB = prune_leaves!(hssQB)

    # reduce to the remainder block (and regain sqare HSS matrix for recursive division)
    hssL = _extract_ccols(hssL, nk) # extract uncompressed rows to form Matrix with one less level
    hssL = prune_leaves!(hssL)

    # recursively call mldivide
    hssY1 = _ldiv!(hssL, hssQB)

    # do the unpacking
    hssB = _unpackadd_rows!(hssY0, hssY1, ktree)
    hssB = _vtransforms!(hssB, QV)
    hssB = recompress!(hssB)
  end
end

## The core routines which make up the division follow here
# core routine to reduce the rows and triangularize the diagonal block via QR/LQ decompositions
function _ulvreduce!(D::Matrix{T}, U::Matrix{T}, V::Matrix{T}) where T
  T <: Complex ? adj = 'C' : adj = 'T'
  m, n = size(D)
  k = size(U, 2)
  nk = min(m-k,n)
  ind = 1:m-k
  cind = m-k+1:m
  # can't be compressed, exit early
  if k >= m
    #@warn "Encountered a full-rank block with k=$(k). Clustering might not yield best performance!"
    L1 = Matrix{T}(undef, 0, 0)
    qlf = (Matrix{T}(undef, m, 0), Vector{T}(undef, 0))
    lqf = (Matrix{T}(undef, 0, n), Vector{T}(undef, 0))
  else
    # form QL decomposition of the row generators and apply it
    qlf = geqlf!(U);
    U = tril(U[end-k+1:end,:]) # k x k block
    ormql!('L', adj, qlf..., D) # transform the diagonal block
    # Form the LQ decomposition of the first m-k rows of D
    lqf = gelqf!(D[1:end-k,:])
    L1 = tril(lqf[1])
    L1 = L1[:,1:nk]
    L2 = ormlq!('R', adj, lqf..., D[end-k+1:end,:]) # update the bottom block of the diagonal block
    V = ormlq!('L', 'N', lqf..., V) # compute the updated off-diagonal generators on the right
    D = L2
  end
  return D, U, V, L1, qlf, lqf, m-k, nk, k
end

# recursive function that operates only on the leaves
function _ulvfactor_leaves!(hssA::HssMatrix{T}, co::Int) where T
  if isleaf(hssA)
    cols = collect(co .+ (1:size(hssA,2)))
    m, n = size(hssA.D)
    D, U, V, L1, qlf, lqf, mk, nk, k = _ulvreduce!(hssA.D, hssA.U, hssA.V)
    # remember which columns got reduced
    QU = BinaryNode(qlf)
    QL = BinaryNode(L1)
    QV = BinaryNode(lqf) # check whether indices are correct
    return HssMatrix(D, U, V, hssA.rootnode), QU, QL, QV, BinaryNode(mk), BinaryNode(nk), BinaryNode(k)
  else
    m1, n1 = hssA.sz1; m2, n2 = hssA.sz2
    hssA.A11, QU1, QL1, QV1, mk1, nk1, k1 = _ulvfactor_leaves!(hssA.A11, co)
    hssA.A22, QU2, QL2, QV2, mk2, nk2, k2 = _ulvfactor_leaves!(hssA.A22, co+n1)
    # update sizes
    hssA.sz1 = size(hssA.A11); hssA.sz2 = size(hssA.A22)
    QU = BinaryNode(QU1, QU2)
    QL = BinaryNode(QL1, QL2)
    QV = BinaryNode(QV1, QV2)
    mk = BinaryNode(mk1, mk2)
    nk = BinaryNode(nk1, nk2)
    ktree = BinaryNode(k1, k2)
    return hssA, QU, QL, QV, mk, nk, ktree # not sure mk, nk are needed
  end
end

## functions to apply U, L and V transforms at the leaf level to hssB
function _utransforms_kernel!(hssA::HssMatrix, Q::BinaryNode)
  eltype(hssA) <: Complex ? adj = 'C' : adj = 'T'
  if size(Q.data[1], 1) == 0 return hssA end
  hssA.D = ormql!('L', adj, Q.data..., hssA.D)
  hssA.U = ormql!('L', adj, Q.data..., hssA.U)
  return hssA
end
function _ltransforms_kernel!(hssA::HssMatrix, Q::BinaryNode)
  nk = size(Q.data,1)
  hssA.D[1:nk, :] = trsm('L', 'L', 'N', 'N', 1., Q.data, hssA.D[1:nk,:])
  hssA.U[1:nk, :] = trsm('L', 'L', 'N', 'N', 1., Q.data, hssA.U[1:nk,:])
  hssA.D[nk+1:end, :] .= eltype(hssA)(0)
  hssA.U[nk+1:end, :] .= eltype(hssA)(0)
  return hssA
end
function _vtransforms_kernel!(hssA::HssMatrix, Q::BinaryNode)
  eltype(hssA) <: Complex ? adj = 'C' : adj = 'T'
  hssA.D = ormlq!('L', adj, Q.data..., hssA.D)
  hssA.U = ormlq!('L', adj, Q.data..., hssA.U)
  return hssA
end
# function _rapplyv!(hssA::HssMatrix, Q::BinaryNode)
#   eltype(hssA) <: Complex ? adj = 'C' : adj = 'T'
#   hssA.D = ormlq!('R', adj, Q.data..., hssA.D)
#   hssA.V = ormlq!('R', adj, Q.data..., hssA.U)
#   return hssA
# end

# generate branch level recursive function for each transform
for (f,g) in zip((:_utransforms!, :_ltransforms!, :_vtransforms!), (:_utransforms_kernel!, :_ltransforms_kernel!, :_vtransforms_kernel!))
  @eval begin
    function $f(hssA::HssMatrix, Qtree::BinaryNode)
      if isleaf(hssA)
        hssA = $g(hssA, Qtree)
      else
        hssA.A11 = $f(hssA.A11, Qtree.left)
        hssA.A22 = $f(hssA.A22, Qtree.right)
        hssA.sz1 = size(hssA.A11); hssA.sz2 = size(hssA.A22)
        return hssA
      end
    end
  end
end

## extractor routines
_extract_nrows_kernel(hssA::HssMatrix, ntree::BinaryNode{Int}) = HssMatrix(hssA.D[1:ntree.data,:], hssA.U[1:ntree.data,:], hssA.V, hssA.rootnode)
_extract_crows_kernel(hssA::HssMatrix, ntree::BinaryNode{Int}) = HssMatrix(hssA.D[ntree.data+1:end,:], hssA.U[ntree.data+1:end,:], hssA.V, hssA.rootnode)
_extract_ncols_kernel(hssA::HssMatrix, ntree::BinaryNode{Int}) = HssMatrix(hssA.D[:,1:ntree.data], hssA.U, hssA.V[1:ntree.data,:], hssA.rootnode)
_extract_ccols_kernel(hssA::HssMatrix, ntree::BinaryNode{Int}) = HssMatrix(hssA.D[:,ntree.data+1:end], hssA.U, hssA.V[ntree.data+1:end,:], hssA.rootnode)

# generate branch routines via metaprogramming
for (f,g) in zip((:_extract_nrows, :_extract_crows, :_extract_ncols, :_extract_ccols), (:_extract_nrows_kernel, :_extract_crows_kernel, :_extract_ncols_kernel, :_extract_ccols_kernel))
  @eval begin
    function $f(hssA::HssMatrix, ntree::BinaryNode{Int})
      if isleaf(hssA)
        hssA = $g(hssA, ntree)
      else
        A11 = $f(hssA.A11, ntree.left)
        A22 = $f(hssA.A22, ntree.right)
        return HssMatrix(A11, A22, hssA.B12, hssA.B21, hssA.R1, hssA.W1, hssA.R2, hssA.W2, hssA.rootnode)
      end
    end
  end
end

## unpacking routine: unpacks hssB and adds it to hssA in place to contribute solution from lower levels
function _unpackadd_rows!(hssA::HssMatrix{T}, hssB::HssMatrix{T}, ktree::BinaryNode{Int}) where T
  if isbranch(hssA)
    if isbranch(hssB)
      hssA.B12 = blkdiagm(hssA.B12, hssB.B12)
      hssA.B21 = blkdiagm(hssA.B21, hssB.B21)

      hssA.R1 = blkdiagm(hssA.R1, hssB.R1)
      hssA.W1 = blkdiagm(hssA.W1, hssB.W1)
      hssA.R2 = blkdiagm(hssA.R2, hssB.R2)
      hssA.W2 = blkdiagm(hssA.W2, hssB.W2)

      hssA.A11 = _unpackadd_rows!(hssA.A11, hssB.A11, ktree.left)
      hssA.A22 = _unpackadd_rows!(hssA.A22, hssB.A22, ktree.right)
      hssA.sz1 == size(hssA.A11) || error("dimensions got screwed up!")
      hssA.sz2 == size(hssA.A22) || error("dimensions got screwed up!")
    else
      isbranch(ktree) || throw(ArgumentError("didn't expect ktree to be a leaf node"))
      k1 = ktree.left.data; k2 = ktree.right.data
      k1 + k2 == size(hssB,1) || throw(DimensionMismatch("first dimension of D does not match the expected k1+k2 rows"))
      m1, n1 = hssA.sz1
      r, w = gensize(hssB)

      # perform the unpacking on the leaves of hssA
      hssA.B12 = blkdiagm(hssA.B12, zeros(T, r, w), Matrix{T}(I, k1, k1))
      hssA.B21 = blkdiagm(hssA.B21, zeros(T, r, w), Matrix{T}(I, k2, k2))

      hssA.R1 = [blkdiagm(hssA.R1, Matrix{T}(I,r,r)); zeros(T, k1, size(hssA.R1,2)+r)]
      hssA.W1 = [blkdiagm(hssA.W1, Matrix{T}(I,w,w)); zeros(T, k2, size(hssA.W1,2)+w)]
      hssA.R2 = [blkdiagm(hssA.R2, Matrix{T}(I,r,r)); zeros(T, k2, size(hssA.R2,2)+r)]
      hssA.W2 = [blkdiagm(hssA.W2, Matrix{T}(I,w,w)); zeros(T, k1, size(hssA.W2,2)+w)]

      # perform the extend-add operation on the leaves
      hssA.A11 = _unpackadd_rows!(hssA.A11, hssB.D[1:k1,1:n1], hssB.U[1:k1,:], hssB.V[1:n1,:], hssB.D[k1+1:end,1:n1]) # probably should add @view
      hssA.A22 = _unpackadd_rows!(hssA.A22, hssB.D[k1+1:end,n1+1:end], hssB.U[k1+1:end,:], hssB.V[n1+1:end,:], hssB.D[1:k1,n1+1:end])

      hssA.sz1 == size(hssA.A11) || error("dimensions got screwed up!")
      hssA.sz2 == size(hssA.A22) || error("dimensions got screwed up!")
    end
    return hssA
  else
    if isbranch(hssB)
      error("Something went wrong. Expected hssB to also be a leaf node")
    else
      k, rk = size(hssB.U);
      k <= size(hssA,1) || throw(DimensionMismatch("expected hssB to have less rows than hssA."))
      hssA.U = [hssA.U zeros(size(hssA.U,1), rk)]
      hssA.U[end-k+1:end, end-rk+1:end] .= hssB.U
      hssA.V = [hssA.V hssB.V]
      hssA.D[end-k+1:end, :] .= hssA.D[end-k+1:end, :] .+ hssB.D
      return hssA
    end
  end
end

function _unpackadd_rows!(hssA::HssMatrix, D::Matrix, U::Matrix, V::Matrix, Brow::Matrix)
  if isleaf(hssA)
    k = size(D,1);
    size(D,1) == size(U,1) || throw(DimensionMismatch("first dimension of U does not match first dimension of D. Expected $(size(D,1)), got $(size(U,1))"))
    size(D,2) == size(V,1)|| throw(DimensionMismatch("first dimension of V does not match second dimension of D. Expected $(size(D,2)), got $(size(V,1))"))
    m, n = size(hssA.D)

    # perform the extend-add operation on the leaves
    hssA.D[end-k+1:end, :] .= hssA.D[end-k+1:end, :] .+ D
    hssA.U = [hssA.U zeros(size(hssA.U,1), size(U,2)+k)]
    hssA.U[end-k+1:end, end-size(U,2)-k+1:1:end] = [U Matrix(I,k,k)]
    hssA.V = [hssA.V V Brow']
    return hssA
  else
    error("Expected hssA to be a leaf node")
  end
end


text:
using Documenter, DataIndexer

makedocs(;
    modules=[DataIndexer],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/lizz/DataIndexer.jl/blob/{commit}{path}#L{line}",
    sitename="DataIndexer.jl",
    authors="lizz",
    assets=String[],
)


text:
#

using Distributions
using Base.Test

dg = Gamma(1,1)

d_s = EdgeworthSum(dg,10)
d_m = EdgeworthMean(dg,10)
d_z = EdgeworthZ(dg,10)

dg_s = Gamma(10,1)
dg_m = Gamma(10,0.1)
dg_za = Gamma(10,1/std(dg_s))

for i = 0.01:0.01:0.99

    q = quantile(dg_s,i)
    @test_approx_eq_eps quantile(d_s,i) q 0.02
    @test_approx_eq_eps cquantile(d_s,1-i) q 0.02
    @test_approx_eq_eps cdf(d_s,q) i 0.002
    @test_approx_eq_eps ccdf(d_s,q) 1-i 0.002
    @test_approx_eq_eps pdf(d_s,q) pdf(dg_s,q) 0.005

    q = quantile(dg_m,i)
    @test_approx_eq_eps quantile(d_m,i) q 0.01
    @test_approx_eq_eps cquantile(d_m,1-i) q 0.01
    @test_approx_eq_eps cdf(d_m,q) i 0.002
    @test_approx_eq_eps ccdf(d_m,q) 1-i 0.002
    @test_approx_eq_eps pdf(d_m,q) pdf(dg_m,q) 0.05

    q = quantile(dg_za,i) - mean(dg_za)
    @test_approx_eq_eps quantile(d_z,i) q 0.01
    @test_approx_eq_eps cquantile(d_z,1-i) q 0.01
    @test_approx_eq_eps cdf(d_z,q) i 0.002
    @test_approx_eq_eps ccdf(d_z,q) 1-i 0.002
    @test_approx_eq_eps pdf(d_z,q) pdf(dg_za,q+mean(dg_za)) 0.02

end


text:
# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

# spaces:
#    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 # hallway
#    12, 13 # Room for As
#    14, 15 # Room for Bs
#    16, 17 # Room for Cs
#    18, 19 # Room for Ds

const MULTIPLIER = Dict(
    :A1 => 1,
    :A2 => 1,
    :B1 => 10,
    :B2 => 10,
    :C1 => 100,
    :C2 => 100,
    :D1 => 1000,
    :D2 => 1000
)

const ENTRY_POINTS = Dict(
    :A1 => 3, :A2 => 3,
    :B1 => 5, :B2 => 5,
    :C1 => 7, :C2 => 7,
    :D1 => 9, :D2 => 9
)

const SIBLING = Dict(
    k1 => k2 for k1 in keys(ENTRY_POINTS) for k2 in keys(ENTRY_POINTS)
    if k1 ≠ k2 && ENTRY_POINTS[k1] == ENTRY_POINTS[k2]
)

# initial configuration for the test

#############
#...........#
###B#C#B#D###
  #A#D#C#A#
  #########

initial_positions_test = Dict(
    :A1 => 13,
    :A2 => 19,
    :B1 => 12,
    :B2 => 16,
    :C1 => 14,
    :C2 => 17,
    :D1 => 15,
    :D2 => 18
)

# An easy one for testing 4 x 1 + 4 x 10 + 2 x 1 = 46 points

#############
#...........#
###B#A#C#D###
  #A#B#C#D#
  #########

initial_positions_test_almostgameoverAB = Dict(
    :A1 => 14,
    :A2 => 13,
    :B1 => 12,
    :B2 => 15,
    :C1 => 16,
    :C2 => 17,
    :D1 => 18,
    :D2 => 19
)

initial_positions_test_almostgameoverCD = Dict(
    :A1 => 12,
    :A2 => 13,
    :B1 => 14,
    :B2 => 15,
    :C1 => 18,
    :C2 => 17,
    :D1 => 16,
    :D2 => 19
)

initial_positions_test_almostgameoverAD = Dict(
    :A1 => 18,
    :A2 => 13,
    :B1 => 14,
    :B2 => 15,
    :C1 => 16,
    :C2 => 17,
    :D1 => 12,
    :D2 => 19
)

#############
#...........#
###D#B#C#A###
  #D#B#C#A#
  #########

initial_positions_test_almostgameoverAADD = Dict(
    :A1 => 18,
    :A2 => 19,
    :B1 => 14,
    :B2 => 15,
    :C1 => 16,
    :C2 => 17,
    :D1 => 12,
    :D2 => 13
)

# initial configuration for the challenge

#############
#...........#
###B#D#C#A###
  #C#D#B#A#
  #########

initial_positions_challenge = Dict(
    :A1 => 18,
    :A2 => 19,
    :B1 => 12,
    :B2 => 17,
    :C1 => 13,
    :C2 => 16,
    :D1 => 14,
    :D2 => 15
)

# Part 1

gameover(positions) = all(
    Set(get.(identity, Ref(positions), keys)) == Set(pos)
    for (keys, pos) in zip(
        ((:A1, :A2), (:B1, :B2), (:C1, :C2), (:D1, :D2)),
        ((i, i+1) for i in 12:2:18)
    )
)

function allowedmoves(positions, k)
    pk = positions[k]
    mvs = []
    e = ENTRY_POINTS[k]

    e9 = e + 9
    e10 = e + 10
    if pk ≤ 11
        if isdisjoint(values(positions), union(e:pk-1, pk+1:e, e9:e10))
            push!(mvs, (e10, length(e:pk-1) + length(pk+1:e) + 2))
        elseif isdisjoint(values(positions), union(e:pk-1, pk+1:e, e9)) &&
            positions[SIBLING[k]] == e10
            push!(mvs, (e9, length(e:pk-1) + length(pk+1:e) + 1))
        end
    elseif pk != e10 && pk:positions[SIBLING[k]] != e9:e10
        if pk in 12:2:18 && isdisjoint(values(positions), union(pk-9:e, e:pk-9, e9:e10))
            push!(mvs, (e10, length(pk-9:e) + length(e:pk-9) + 2))
        elseif pk in 12:2:18 && isdisjoint(values(positions), union(pk-9:e, e:pk-9, e9)) &&
            positions[SIBLING[k]] == e10
            push!(mvs, (e9, length(pk-9:e) + length(e:pk-9) + 1))
        elseif pk in 13:2:19 && isdisjoint(values(positions), union(pk-1, pk-10:e, e:pk-10, e9:e10))
            push!(mvs, (e10, length(pk-10:e) + length(e:pk-10) + 3))
        elseif pk in 13:2:19 && isdisjoint(values(positions), union(pk-1, pk-10:e, e:pk-10, e9)) &&
            positions[SIBLING[k]] == e10
            push!(mvs, (e9, length(pk-10:e) + length(e:pk-10) + 2))
        else
            for p in (1, 2, 4, 6, 8, 10, 11)
                if pk in 12:2:18 && isdisjoint(values(positions), union(p:pk-9, pk-9:p))
                    push!(mvs, (p, length(p:pk-9) + length(pk-9:p)))
                elseif isdisjoint(values(positions), union(pk-1, p:pk-10, pk-10:p))
                    push!(mvs, (p, length(p:pk-10) + length(pk-10:p) + 1))
                end
            end
        end
    end
    return mvs
end

function move(positions, points = 0, minpoints = typemax(Int),
    moves = Vector{Tuple{Symbol, Pair{Int, Int}, Int}}(), minmoves = Vector{Tuple{Symbol, Pair{Int, Int}, Int}}())
    if points < minpoints && gameover(positions)
        minpoints = points
        minmoves = moves
    else
        for k in keys(positions)
            for (pos, steps) in allowedmoves(positions, k)
                step_pts = MULTIPLIER[k] * steps
                newpoints = points + step_pts
                if newpoints < minpoints
                    newpositions = copy(positions)
                    newpositions[k] = pos
                    newmoves = [moves; (k, positions[k] => pos, step_pts)]
                    minmoves, minpoints = move(newpositions, newpoints, minpoints, newmoves, minmoves)
                end
            end
        end
    end
    return minmoves, minpoints
end

function part1(initial_positions)
    move(initial_positions)
end

testit("Test Part 1 AB", part1(initial_positions_test_almostgameoverAB)[2], ==(46))
testit("Test Part 1 CD", part1(initial_positions_test_almostgameoverCD)[2], ==(4600))
testit("Test Part 1 DC", part1(initial_positions_test_almostgameoverAD)[2], ==(8010))
testit("Test Part 1 AADD", part1(initial_positions_test_almostgameoverAADD)[2], ==(18022))
testit("Test Part 1", part1(initial_positions_test)[2], ==(12521))
testit("Challenge Part 1", part1(initial_positions_challenge)[2], ==(15322))


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

module TestCholesky

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted,
    PosDefException, RankDeficientException, chkfullrank

function unary_ops_tests(a, ca, tol; n=size(a, 1))
    @test inv(ca)*a ≈ Matrix(I, n, n)
    @test a*inv(ca) ≈ Matrix(I, n, n)
    @test abs((det(ca) - det(a))/det(ca)) <= tol # Ad hoc, but statistically verified, revisit
    @test logdet(ca) ≈ logdet(a)
    @test logdet(ca) ≈ log(det(ca))  # logdet is less likely to overflow
    @test isposdef(ca)
    @test_throws ErrorException ca.Z
    @test size(ca) == size(a)
    @test Array(copy(ca)) ≈ a
end

function factor_recreation_tests(a_U, a_L)
    c_U = cholesky(a_U)
    c_L = cholesky(a_L)
    cl  = c_L.U
    ls = c_L.L
    @test Array(c_U) ≈ Array(c_L) ≈ a_U
    @test ls*ls' ≈ a_U
    @test triu(c_U.factors) ≈ c_U.U
    @test tril(c_L.factors) ≈ c_L.L
    @test istriu(cl)
    @test cl'cl ≈ a_U
    @test cl'cl ≈ a_L
end

@testset "core functionality" begin
    n = 10

    # Split n into 2 parts for tests needing two matrices
    n1 = div(n, 2)
    n2 = 2*n1

    Random.seed!(12344)

    areal = randn(n,n)/2
    aimg  = randn(n,n)/2
    a2real = randn(n,n)/2
    a2img  = randn(n,n)/2
    breal = randn(n,2)/2
    bimg  = randn(n,2)/2

    for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)

        ε = εa = eps(abs(float(one(eltya))))

        # Test of symmetric pos. def. strided matrix
        apd  = a'*a
        @inferred cholesky(apd)
        capd  = factorize(apd)
        r     = capd.U
        κ     = cond(apd, 1) #condition number

        unary_ops_tests(apd, capd, ε*κ*n)
        if eltya != Int
            @test Factorization{eltya}(capd) === capd
            if eltya <: Real
                @test Array(Factorization{complex(eltya)}(capd)) ≈ Array(factorize(complex(apd)))
                @test eltype(Factorization{complex(eltya)}(capd)) == complex(eltya)
            end
        end
        @testset "throw for non-square input" begin
            A = rand(eltya, 2, 3)
            @test_throws DimensionMismatch cholesky(A)
            @test_throws DimensionMismatch cholesky!(A)
        end

        #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1

        #these tests were failing on 64-bit linux when inside the inner loop
        #for eltya = ComplexF32 and eltyb = Int. The E[i,j] had NaN32 elements
        #but only with Random.seed!(1234321) set before the loops.
        E = abs.(apd - r'*r)
        for i=1:n, j=1:n
            @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j]))
        end
        E = abs.(apd - Matrix(capd))
        for i=1:n, j=1:n
            @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j]))
        end
        @test LinearAlgebra.issuccess(capd)
        @inferred(logdet(capd))

        apos = apd[1,1]
        @test all(x -> x ≈ √apos, cholesky(apos).factors)

        # Test cholesky with Symmetric/Hermitian upper/lower
        apds  = Symmetric(apd)
        apdsL = Symmetric(apd, :L)
        apdh  = Hermitian(apd)
        apdhL = Hermitian(apd, :L)
        if eltya <: Real
            capds = cholesky(apds)
            unary_ops_tests(apds, capds, ε*κ*n)
            if eltya <: BlasReal
                capds = cholesky!(copy(apds))
                unary_ops_tests(apds, capds, ε*κ*n)
            end
            ulstring = sprint((t, s) -> show(t, "text/plain", s), capds.UL)
            @test sprint((t, s) -> show(t, "text/plain", s), capds) == "$(typeof(capds))\nU factor:\n$ulstring"
        else
            capdh = cholesky(apdh)
            unary_ops_tests(apdh, capdh, ε*κ*n)
            capdh = cholesky!(copy(apdh))
            unary_ops_tests(apdh, capdh, ε*κ*n)
            capdh = cholesky!(copy(apd))
            unary_ops_tests(apd, capdh, ε*κ*n)
            ulstring = sprint((t, s) -> show(t, "text/plain", s), capdh.UL)
            @test sprint((t, s) -> show(t, "text/plain", s), capdh) == "$(typeof(capdh))\nU factor:\n$ulstring"
        end

        # test cholesky of 2x2 Strang matrix
        S = Matrix{eltya}(SymTridiagonal([2, 2], [-1]))
        @test Matrix(cholesky(S).U) ≈ [2 -1; 0 sqrt(eltya(3))] / sqrt(eltya(2))

        # test extraction of factor and re-creating original matrix
        if eltya <: Real
            factor_recreation_tests(apds, apdsL)
        else
            factor_recreation_tests(apdh, apdhL)
        end

        #pivoted upper Cholesky
        if eltya != BigFloat
            cpapd = cholesky(apdh, Val(true))
            unary_ops_tests(apdh, cpapd, ε*κ*n)
            @test rank(cpapd) == n
            @test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing

            @test cpapd.P*cpapd.L*cpapd.U*cpapd.P' ≈ apd
        end

        for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
            b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
            εb = eps(abs(float(one(eltyb))))
            ε = max(εa,εb)

            for b in (b, view(b, 1:n, 1)) # Array and SubArray

                # Test error bound on linear solver: LAWNS 14, Theorem 2.1
                # This is a surprisingly loose bound
                x = capd\b
                @test norm(x-apd\b,1)/norm(x,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                @test norm(apd*x-b,1)/norm(b,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ

                @test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit

                if eltya != BigFloat && eltyb != BigFloat
                    lapd = cholesky(apdhL)
                    @test norm(apd * (lapd\b) - b)/norm(b) <= ε*κ*n
                    @test norm(apd * (lapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n
                end

                if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted Cholesky decomposition in julia

                    cpapd = cholesky(apdh, Val(true))
                    @test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
                    @test norm(apd * (cpapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n

                    lpapd = cholesky(apdhL, Val(true))
                    @test norm(apd * (lpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
                    @test norm(apd * (lpapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n
                end
            end
        end

        for eltyb in (Float64, ComplexF64)
            Breal = convert(Matrix{BigFloat}, randn(n,n)/2)
            Bimg  = convert(Matrix{BigFloat}, randn(n,n)/2)
            B = (eltya <: Complex || eltyb <: Complex) ? complex.(Breal, Bimg) : Breal
            εb = eps(abs(float(one(eltyb))))
            ε = max(εa,εb)

            for B in (B, view(B, 1:n, 1:n)) # Array and SubArray

                # Test error bound on linear solver: LAWNS 14, Theorem 2.1
                # This is a surprisingly loose bound
                BB = copy(B)
                ldiv!(capd, BB)
                @test norm(apd \ B - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                @test norm(apd * BB - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                if eltya != BigFloat
                    cpapd = cholesky(apdh, Val(true))
                    BB = copy(B)
                    ldiv!(cpapd, BB)
                    @test norm(apd \ B - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                    @test norm(apd * BB - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                end
            end
        end

        @testset "solve with generic Cholesky" begin
            Breal = convert(Matrix{BigFloat}, randn(n,n)/2)
            Bimg  = convert(Matrix{BigFloat}, randn(n,n)/2)
            B = eltya <: Complex ? complex.(Breal, Bimg) : Breal
            εb = eps(abs(float(one(eltype(B)))))
            ε = max(εa,εb)

            for B in (B, view(B, 1:n, 1:n)) # Array and SubArray

                # Test error bound on linear solver: LAWNS 14, Theorem 2.1
                # This is a surprisingly loose bound
                cpapd = cholesky(eltya <: Complex ? apdh : apds)
                BB = copy(B)
                rdiv!(BB, cpapd)
                @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                cpapd = cholesky(eltya <: Complex ? apdhL : apdsL)
                BB = copy(B)
                rdiv!(BB, cpapd)
                @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                if eltya != BigFloat
                    cpapd = cholesky(eltya <: Complex ? apdh : apds, Val(true))
                    BB = copy(B)
                    rdiv!(BB, cpapd)
                    @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                    @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                    cpapd = cholesky(eltya <: Complex ? apdhL : apdsL, Val(true))
                    BB = copy(B)
                    rdiv!(BB, cpapd)
                    @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                    @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
                end
            end
        end
        if eltya <: BlasFloat
            @testset "generic cholesky!" begin
                if eltya <: Complex
                    A = complex.(randn(5,5), randn(5,5))
                else
                    A = randn(5,5)
                end
                A = convert(Matrix{eltya}, A'A)
                @test Matrix(cholesky(A).L) ≈ Matrix(invoke(LinearAlgebra._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular)[1])
                @test Matrix(cholesky(A).U) ≈ Matrix(invoke(LinearAlgebra._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular)[1])
            end
        end
    end
end

@testset "behavior for non-positive definite matrices" for T in (Float64, ComplexF64)
    A = T[1 2; 2 1]
    B = T[1 2; 0 1]
    # check = (true|false)
    for M in (A, Hermitian(A), B)
        @test_throws PosDefException cholesky(M)
        @test_throws PosDefException cholesky!(copy(M))
        @test_throws PosDefException cholesky(M; check = true)
        @test_throws PosDefException cholesky!(copy(M); check = true)
        @test !LinearAlgebra.issuccess(cholesky(M; check = false))
        @test !LinearAlgebra.issuccess(cholesky!(copy(M); check = false))
    end
    for M in (A, Hermitian(A), B)
        @test_throws RankDeficientException cholesky(M, Val(true))
        @test_throws RankDeficientException cholesky!(copy(M), Val(true))
        @test_throws RankDeficientException cholesky(M, Val(true); check = true)
        @test_throws RankDeficientException cholesky!(copy(M), Val(true); check = true)
        @test !LinearAlgebra.issuccess(cholesky(M, Val(true); check = false))
        @test !LinearAlgebra.issuccess(cholesky!(copy(M), Val(true); check = false))
        C = cholesky(M, Val(true); check = false)
        @test_throws RankDeficientException chkfullrank(C)
        C = cholesky!(copy(M), Val(true); check = false)
        @test_throws RankDeficientException chkfullrank(C)
    end
    @test !isposdef(A)
    str = sprint((io, x) -> show(io, "text/plain", x), cholesky(A; check = false))
end

@testset "Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices" begin
    X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3]
    L = Matrix(LinearAlgebra._chol!(X*X', LowerTriangular)[1])
    U = Matrix(LinearAlgebra._chol!(X*X', UpperTriangular)[1])
    XX = Matrix(X*X')

    @test sum(sum(norm, L*L' - XX)) < eps()
    @test sum(sum(norm, U'*U - XX)) < eps()
end

struct WrappedVector{T} <: AbstractVector{T}
    data::Vector{T}
end
Base.copy(v::WrappedVector) = WrappedVector(copy(v.data))
Base.size(v::WrappedVector) = size(v.data)
Base.getindex(v::WrappedVector, i::Integer) = getindex(v.data, i)
Base.setindex!(v::WrappedVector, val, i::Integer) = setindex!(v.data, val, i)

@testset "cholesky up- and downdates" begin
    A = complex.(randn(10,5), randn(10, 5))
    v = complex.(randn(5), randn(5))
    w = WrappedVector(v)
    for uplo in (:U, :L)
        AcA = A'*A
        BcB = AcA + v*v'
        BcB = (BcB + BcB')/2
        F = cholesky(Hermitian(AcA, uplo))
        G = cholesky(Hermitian(BcB, uplo))
        @test getproperty(lowrankupdate(F, v), uplo) ≈ getproperty(G, uplo)
        @test getproperty(lowrankupdate(F, w), uplo) ≈ getproperty(G, uplo)
        @test_throws DimensionMismatch lowrankupdate(F, Vector{eltype(v)}(undef,length(v)+1))
        @test getproperty(lowrankdowndate(G, v), uplo) ≈ getproperty(F, uplo)
        @test getproperty(lowrankdowndate(G, w), uplo) ≈ getproperty(F, uplo)
        @test_throws DimensionMismatch lowrankdowndate(G, Vector{eltype(v)}(undef,length(v)+1))
    end
end

@testset "issue #13243, unexpected nans in complex cholesky" begin
    apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.3824869f0im 2.619998f0 + 1.8532984f0im -1.8306153f0 - 1.2336911f0im 0.32275113f0 + 0.015575029f0im 2.1968813f0 + 1.0640624f0im 0.27894387f0 + 0.97911835f0im 3.0476584f0 + 0.18548489f0im 0.3842994f0 + 0.7050991f0im
        -0.79540455f0 - 0.7066077f0im 8.313246f0 + 0.0f0im -1.8076122f0 - 0.8882447f0im 0.47806996f0 + 0.48494184f0im 0.5096429f0 - 0.5395974f0im -0.7285097f0 - 0.10360408f0im -1.1760061f0 - 2.7146957f0im -0.4271084f0 + 0.042899966f0im -1.7228563f0 + 2.8335886f0im 1.8942566f0 + 0.6389735f0im
        0.98274714f0 - 1.3824869f0im -1.8076122f0 + 0.8882447f0im 9.367975f0 + 0.0f0im -0.1838578f0 + 0.6468568f0im -1.8338387f0 + 0.7064959f0im 0.041852742f0 - 0.6556877f0im 2.5673025f0 + 1.9732997f0im -1.1148382f0 - 0.15693812f0im 2.4704504f0 - 1.0389464f0im 1.0858271f0 - 1.298006f0im
        2.619998f0 - 1.8532984f0im 0.47806996f0 - 0.48494184f0im -0.1838578f0 - 0.6468568f0im 3.1117508f0 + 0.0f0im -1.956626f0 + 0.22825956f0im 0.07081801f0 - 0.31801307f0im 0.3698375f0 - 0.5400855f0im 0.80686307f0 + 1.5315914f0im 1.5649154f0 - 1.6229297f0im -0.112077385f0 + 1.2014246f0im
        -1.8306153f0 + 1.2336911f0im 0.5096429f0 + 0.5395974f0im -1.8338387f0 - 0.7064959f0im -1.956626f0 - 0.22825956f0im 3.6439795f0 + 0.0f0im -0.2594722f0 + 0.48786148f0im -0.47636223f0 - 0.27821827f0im -0.61608654f0 - 2.01858f0im -2.7767487f0 + 1.7693765f0im 0.048102796f0 - 0.9741874f0im
        0.32275113f0 - 0.015575029f0im -0.7285097f0 + 0.10360408f0im 0.041852742f0 + 0.6556877f0im 0.07081801f0 + 0.31801307f0im -0.2594722f0 - 0.48786148f0im 3.624376f0 + 0.0f0im -1.6697118f0 + 0.4017511f0im -1.4397877f0 - 0.7550918f0im -0.31456697f0 - 1.0403451f0im -0.31978557f0 + 0.13701046f0im
        2.1968813f0 - 1.0640624f0im -1.1760061f0 + 2.7146957f0im 2.5673025f0 - 1.9732997f0im 0.3698375f0 + 0.5400855f0im -0.47636223f0 + 0.27821827f0im -1.6697118f0 - 0.4017511f0im 6.8273163f0 + 0.0f0im -0.10051322f0 + 0.24303961f0im 1.4415971f0 + 0.29750675f0im 1.221786f0 - 0.85654986f0im
        0.27894387f0 - 0.97911835f0im -0.4271084f0 - 0.042899966f0im -1.1148382f0 + 0.15693812f0im 0.80686307f0 - 1.5315914f0im -0.61608654f0 + 2.01858f0im -1.4397877f0 + 0.7550918f0im -0.10051322f0 - 0.24303961f0im 3.4057708f0 + 0.0f0im -0.5856801f0 - 1.0203559f0im 0.7103452f0 + 0.8422135f0im
        3.0476584f0 - 0.18548489f0im -1.7228563f0 - 2.8335886f0im 2.4704504f0 + 1.0389464f0im 1.5649154f0 + 1.6229297f0im -2.7767487f0 - 1.7693765f0im -0.31456697f0 + 1.0403451f0im 1.4415971f0 - 0.29750675f0im -0.5856801f0 + 1.0203559f0im 7.005772f0 + 0.0f0im -0.9617417f0 - 1.2486815f0im
        0.3842994f0 - 0.7050991f0im 1.8942566f0 - 0.6389735f0im 1.0858271f0 + 1.298006f0im -0.112077385f0 - 1.2014246f0im 0.048102796f0 + 0.9741874f0im -0.31978557f0 - 0.13701046f0im 1.221786f0 + 0.85654986f0im 0.7103452f0 - 0.8422135f0im -0.9617417f0 + 1.2486815f0im 3.4629636f0 + 0.0f0im]
    b = [-0.905011814118756 + 0.2847570854574069im -0.7122162951294634 - 0.630289556702497im
        -0.7620356655676837 + 0.15533508334193666im 0.39947219167701153 - 0.4576746001199889im
        -0.21782716937787788 - 0.9222220085490986im -0.727775859267237 + 0.50638268521728im
        -1.0509472322215125 + 0.5022165705328413im -0.7264975746431271 + 0.31670415674097235im
        -0.6650468984506477 - 0.5000967284800251im -0.023682508769195098 + 0.18093440285319276im
        -0.20604111555491242 + 0.10570814584017311im 0.562377322638969 - 0.2578030745663871im
        -0.3451346708401685 + 1.076948486041297im 0.9870834574024372 - 0.2825689605519449im
        0.25336108035924787 + 0.975317836492159im 0.0628393808469436 - 0.1253397353973715im
        0.11192755545114 - 0.1603741874112385im 0.8439562576196216 + 1.0850814110398734im
        -1.0568488936791578 - 0.06025820467086475im 0.12696236014017806 - 0.09853584666755086im]
    cholesky(Hermitian(apd, :L), Val(true)) \ b
    r = factorize(apd).U
    E = abs.(apd - r'*r)
    ε = eps(abs(float(one(ComplexF32))))
    n = 10
    for i=1:n, j=1:n
        @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j]))
    end
end

@testset "fail for non-BLAS element types" begin
    @test_throws ArgumentError cholesky!(Hermitian(rand(Float16, 5,5)), Val(true))
end

@testset "cholesky Diagonal" begin
    # real
    d = abs.(randn(3)) .+ 0.1
    D = Diagonal(d)
    CD = cholesky(D)
    CM = cholesky(Matrix(D))
    @test CD isa Cholesky{Float64}
    @test CD.U ≈ Diagonal(.√d) ≈ CM.U
    @test D ≈ CD.L * CD.U
    @test CD.info == 0

    # real, failing
    @test_throws PosDefException cholesky(Diagonal([1.0, -2.0]))
    Dnpd = cholesky(Diagonal([1.0, -2.0]); check = false)
    @test Dnpd.info == 2

    # complex
    D = complex(D)
    CD = cholesky(D)
    CM = cholesky(Matrix(D))
    @test CD isa Cholesky{ComplexF64}
    @test CD.U ≈ Diagonal(.√d) ≈ CM.U
    @test D ≈ CD.L * CD.U
    @test CD.info == 0

    # complex, failing
    D[2, 2] = 0.0 + 0im
    @test_throws PosDefException cholesky(D)
    Dnpd = cholesky(D; check = false)
    @test Dnpd.info == 2

    # InexactError for Int
    @test_throws InexactError cholesky!(Diagonal([2, 1]))
end

@testset "constructor with non-BlasInt arguments" begin

    x = rand(5,5)
    chol = cholesky(x'x)

    factors, uplo, info = chol.factors, chol.uplo, chol.info

    @test Cholesky(factors, uplo, Int32(info)) == chol
    @test Cholesky(factors, uplo, Int64(info)) == chol

    cholp = cholesky(x'x, Val(true))

    factors, uplo, piv, rank, tol, info =
        cholp.factors, cholp.uplo, cholp.piv, cholp.rank, cholp.tol, cholp.info

    @test CholeskyPivoted(factors, uplo, Vector{Int32}(piv), rank, tol, info) == cholp
    @test CholeskyPivoted(factors, uplo, Vector{Int64}(piv), rank, tol, info) == cholp

    @test CholeskyPivoted(factors, uplo, piv, Int32(rank), tol, info) == cholp
    @test CholeskyPivoted(factors, uplo, piv, Int64(rank), tol, info) == cholp

    @test CholeskyPivoted(factors, uplo, piv, rank, tol, Int32(info)) == cholp
    @test CholeskyPivoted(factors, uplo, piv, rank, tol, Int64(info)) == cholp

end

@testset "issue #33704, casting low-rank CholeskyPivoted to Matrix" begin
    A = randn(1,8)
    B = A'A
    C = cholesky(B, Val(true), check=false)
    @test B ≈ Matrix(C)
end

@testset "CholeskyPivoted and Factorization" begin
    A = randn(8,8)
    B = A'A
    C = cholesky(B, Val(true), check=false)
    @test CholeskyPivoted{eltype(C)}(C) === C
    @test Factorization{eltype(C)}(C) === C
    @test Array(CholeskyPivoted{complex(eltype(C))}(C)) ≈ Array(cholesky(complex(B), Val(true), check=false))
    @test Array(Factorization{complex(eltype(C))}(C)) ≈ Array(cholesky(complex(B), Val(true), check=false))
    @test eltype(Factorization{complex(eltype(C))}(C)) == complex(eltype(C))
end

@testset "REPL printing of CholeskyPivoted" begin
    A = randn(8,8)
    B = A'A
    C = cholesky(B, Val(true), check=false)
    cholstring = sprint((t, s) -> show(t, "text/plain", s), C)
    rankstring = "$(C.uplo) factor with rank $(rank(C)):"
    factorstring = sprint((t, s) -> show(t, "text/plain", s), C.uplo == 'U' ? C.U : C.L)
    permstring   = sprint((t, s) -> show(t, "text/plain", s), C.p)
    @test cholstring == "$(summary(C))\n$rankstring\n$factorstring\npermutation:\n$permstring"
end

@testset "destructuring for Cholesky[Pivoted]" begin
    for val in (true, false)
        A = rand(8, 8)
        B = A'A
        C = cholesky(B, Val(val), check=false)
        l, u = C
        @test l == C.L
        @test u == C.U
    end
end

@testset "issue #37356, diagonal elements of hermitian generic matrix" begin
    B = Hermitian(hcat([one(BigFloat) + im]))
    @test Matrix(cholesky(B)) ≈ B
    C = Hermitian(hcat([one(BigFloat) + im]), :L)
    @test Matrix(cholesky(C)) ≈ C
end

@testset "constructing a Cholesky factor from a triangular matrix" begin
    A = [1.0 2.0; 3.0 4.0]
    let
        U = UpperTriangular(A)
        C = Cholesky(U)
        @test C isa Cholesky{Float64}
        @test C.U == U
        @test C.L == U'
    end
    let
        L = LowerTriangular(A)
        C = Cholesky(L)
        @test C isa Cholesky{Float64}
        @test C.L == L
        @test C.U == L'
    end
end

@testset "adjoint of Cholesky" begin
    A = randn(5, 5)
    A = A'A
    F = cholesky(A)
    b = ones(size(A, 1))
    @test F\b == F'\b
end

@testset "Float16" begin
    A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
    B = cholesky(A)
    B32 = cholesky(Float32.(A))
    @test B isa Cholesky{Float16, Matrix{Float16}}
    @test B.U isa UpperTriangular{Float16, Matrix{Float16}}
    @test B.L isa LowerTriangular{Float16, Matrix{Float16}}
    @test B.UL isa UpperTriangular{Float16, Matrix{Float16}}
    @test B.U ≈ B32.U
    @test B.L ≈ B32.L
    @test B.UL ≈ B32.UL
end

@testset "det and logdet" begin
    A = [4083 3825 5876 2048 4470 5490;
         3825 3575 5520 1920 4200 5140;
         5876 5520 8427 2940 6410 7903;
         2048 1920 2940 1008 2240 2740;
         4470 4200 6410 2240 4875 6015;
         5490 5140 7903 2740 6015 7370]
    B = cholesky(A, Val(true), check=false)
    @test det(B)  ==  0.0
    @test det(B)  ≈  det(A) atol=eps()
    @test logdet(B)  ==  -Inf
 end

end # module TestCholesky

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
# TODO: use take and drop instead of firstn given in Julia base.

function firstn(a,n)
    b=Any[]
    for i in a
        push!(b,i)
        length(b)==n && break
    end
    b
end

function firstn(a::SketchEngine,n)
    D=length(a.dict)
    d=cell(D)
    for (w,i) in a.dict
        d[i]=w
    end
    b=Any[]
    for s in a
        push!(b,map(i->(i>D?"<unk>":d[i]), s))
        length(b)==n && break
    end
    b
end

function onerows(x::SparseMatrixCSC,dict::Dict)
    D=length(dict)+1
    d=cell(D)
    for (w,i) in dict
        d[i]=w
    end
    d[end]="<unk>"
    map(1:size(x,2)) do i
        x.nzval[i] == 0 ? "0" : d[x.rowval[i]]
    end
end

function onerows(x::SparseMatrixCSC)
    map(1:size(x,2)) do i
        x.nzval[i] == 0 ? 0 : x.rowval[i]
    end
end

function firstn(a::TagData,n)
    b=Any[]
    for i in a
        if i==nothing
            push!(b,i)
        else
            (x,y,m) = i
            !isa(a.x.sgen,Array) && isdefined(a.x.sgen,:dict) ?
            push!(b, (convert(Vector{Int},m), onerows(x,a.x.sgen.dict), onerows(y,a.y.sgen.dict))) :
            push!(b, map(z->convert(Vector{Int},z), (m, onerows(x), onerows(y))))
        end
        length(b)==n && break
    end
    b
end