text:
function check_dataframe_in_setter(value::DataFrame, expected::Int, message::String)
    if size(value)[1] != expected
        throw(DimensionMismatch("'value' and '" * message *"(x)' should have the same number of rows"))
    end

    if !check_dataframe_has_name(value)
        throw(DimensionMismatch("first column of 'value' should exist and be called 'name'"))
    end

    if !check_dataframe_firstcol(value)
        throw(ArgumentError("first column of 'value' should contain strings or nothings"));
    end
end

"""
    setrowdata!(x, value)

Set the row annotations in `x` to `value`.

If `value` is a `DataFrame`, the first column should be called `"name"` and contain the row names of `x`;
this can either be an `AbstractVector{AbstractString}` or a `Vector{Nothing}` (if no row names are available).

If `value` is `nothing`, this is considered to be equivalent to a `DataFrame` with one `"name"` column containing `nothing`s.

The return value is a reference to the modified `x`.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> # using DataFrames

julia> replacement = copy(rowdata(x));

julia> replacement[!,"foobar"] = [ rand() for i in 1:size(x)[1] ];

julia> setrowdata!(x, replacement);

julia> names(rowdata(x))
3-element Vector{String}:
 "name"
 "Type"
 "foobar"
```
"""
function setrowdata!(x::SummarizedExperiment, value::Union{DataFrame,Nothing})
    if isa(value, Nothing)
        x.rowdata = DataFrame(name = Vector{Nothing}(undef, size(x)[1]))
    else
        check_dataframe_in_setter(value, size(x)[1], "rowdata")
        x.rowdata = value
    end
    return x
end

"""
    setcoldata!(x, value)

Set the column annotations in `x` to `value`.

If `value` is a `DataFrame`, the first column should be called `"name"` and contain the column names of `x`;
this can either be a `Vector{String}` or a `Vector{Nothing}` (if no column names are available).

If `value` is `nothing`, this is considered to be equivalent to a `DataFrame` with one `"name"` column containing `nothing`s.

The return value is a reference to the modified `x`.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> replacement = copy(coldata(x));

julia> replacement[!,"foobar"] = [ rand() for i in 1:size(x)[2] ];

julia> setcoldata!(x, replacement);

julia> names(coldata(x))
4-element Vector{String}:
 "name"
 "Treatment"
 "Response"
 "foobar"
```
"""
function setcoldata!(x::SummarizedExperiment, value::Union{DataFrame,Nothing})
    if isa(value, Nothing)
        x.coldata = DataFrame(name = Vector{Nothing}(undef, size(x)[2]))
    else
        check_dataframe_in_setter(value, size(x)[2], "coldata")
        x.coldata = value
    end
    return x
end

"""
    setassay!(x[, i], value)

Set the requested assay in `x` to any array-like `value`.
The first two dimensions of `value` must have extent equal to those of `x`.

`i` may be an integer specifying an index, in which case it must be positive and no greater than `length(assays(x))`;
or a string containing the name, in which case it may be an existing or new name.
If `i` is not supplied, `value` is set as the first assay of `x`.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> first_sum = sum(assay(x));

julia> second_sum = sum(assay(x, 2));

julia> setassay!(x, assay(x, 2)); # Replacing the first assay with the second.

julia> first_sum == sum(assay(x))
false

julia> second_sum == sum(assay(x))
true

julia> setassay!(x, 1, assay(x, 2)); # More explicit forms of the above.

julia> setassay!(x, "foo", assay(x, 2));
```
"""
function setassay!(x::SummarizedExperiment, value::AbstractArray)
    setassay!(x, 1, value)
    return
end

function setassay!(x::SummarizedExperiment, i::Int64, value::AbstractArray)
    dim = size(value)
    xdim = size(x)
    if dim[1] != xdim[1] || dim[2] != xdim[2]
        throw(DimensionMismatch("dimensions of 'value' should be the same as those of 'x'"))
    end

    counter = 0
    for k in keys(x.assays)
        counter += 1
        if counter == i
            x.assays[k] = value
            return
        end
    end

    throw(BoundsError("'i = " * string(i) * "' is out of range of 'assays(x)'"))
end

function setassay!(x::SummarizedExperiment, i::String, value::AbstractArray)
    dim = size(value)
    xdim = size(x)
    if dim[1] != xdim[1] || dim[2] != xdim[2]
        throw(DimensionMismatch("dimensions of 'value' should be the same as those of 'x'"))
    end
    x.assays[i] = value
    return
end

"""
    setassays!(x, value)

Set assays in `x` to `value`, an `OrderedDict` where the keys are assay names and the values are arrays.
All arrays in `value` should have the same extents as `x` for the first two dimensions.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> length(assays(x))
3

julia> refresh = copy(assays(x));

julia> delete!(refresh, "foo");

julia> setassays!(x, refresh)

julia> length(assays(x))
2
```
"""
function setassays!(x::SummarizedExperiment, value::OrderedDict{String,AbstractArray})
    xdim = size(x)
    for (name, val) in value
        dim = size(val)
        if dim[1] != xdim[1] || dim[2] != xdim[2]
            throw(DimensionMismatch("dimensions of 'value[" * repr(name) * "]' should be the same as those of 'x'"))
        end
    end
    x.assays = value
    return
end

"""
    setmetadata!(x, value)

Set metadata in `x` to `value`, a `Dict` where the keys are the metadata field names.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> setmetadata!(x, Dict{String,Any}("foo" => 200));

julia> metadata(x)
Dict{String, Any} with 1 entry:
  "foo" => 200
```
"""
function setmetadata!(x::SummarizedExperiment, value::Dict{String,Any})
    x.metadata = value;
    return
end


text:
using RecipesBase

x = [x ± x/10 for x in 0:0.1:10];
y = sin.(x);

rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), y)
@test rec[1].d == Dict{Symbol, Any}(:yerror => uncertainty.(y))
@test rec[1].args == (value.(y),)

rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), sin, x)
@test rec[1].d == Dict{Symbol, Any}(:xerror => uncertainty.(x),
                                    :yerror => uncertainty.(y))
@test rec[1].args == (value.(x), value.(y))

rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), x, y)
@test rec[1].d == Dict{Symbol, Any}(:xerror => uncertainty.(x),
                                    :yerror => uncertainty.(y))
@test rec[1].args == (value.(x), value.(y))

rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), x, value.(y))
@test rec[1].d == Dict{Symbol, Any}(:xerror => uncertainty.(x))
@test rec[1].args == (value.(x), value.(y))

rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), value.(x), y)
@test rec[1].d == Dict{Symbol, Any}(:yerror => uncertainty.(y))
@test rec[1].args == (value.(x), value.(y))

ext: jl


text:
"""
The `NomnomlJS` package provides an interface to the `nomnoml` JavaScript
library for rendering UML diagrams from `.noml` strings and via the
`read(filename, Diagram)` overload.

The public interface for this package is the `Diagram` type, which constructs
a `nomnoml` diagram. This object can then be written to several different
file formats via `write(filename, diagram)`.

The following list of support syntax and features is quoted from the upstream
repository, [nomnoml](https://github.com/skanaar/nomnoml):

### Association types

```plaintext
-    association
->   association
<->  association
-->  dependency
<--> dependency
-:>  generalization
<:-  generalization
--:> implementation
<:-- implementation
+-   composition
+->  composition
o-   aggregation
o->  aggregation
--   note
-/-  hidden
_>   weightless edge
__   weightless dashed edge
```

### Classifier types

```plaintext
[name]
[<abstract> name]
[<instance> name]
[<reference> name]
[<note> name]
[<package> name]
[<frame> name]
[<database> name]
[<start> name]
[<end> name]
[<state> name]
[<choice> name]
[<sync> name]
[<input> name]
[<sender> name]
[<receiver> name]
[<transceiver> name]
[<actor> name]
[<usecase> name]
[<label> name]
[<hidden> name]
[<table> name| a | 5 || b | 7]
```

### Directives

```plaintext
#import: my-common-styles.nomnoml
#arrowSize: 1
#bendSize: 0.3
#direction: down | right
#gutter: 5
#edgeMargin: 0
#gravity: 1
#edges: hard | rounded
#background: transparent
#fill: #eee8d5; #fdf6e3
#fillArrows: false
#font: Calibri
#fontSize: 12
#leading: 1.25
#lineWidth: 3
#padding: 8
#spacing: 40
#stroke: #33322E
#title: filename
#zoom: 1
#acyclicer: greedy
#ranker: network-simplex | tight-tree | longest-path
```

### Custom classifier styles

A directive that starts with "." define a classifier style. The style is
written as a space separated list of modifiers and key/value pairs.

```plaintext
#.box: fill=#8f8 dashed
#.blob: visual=ellipse title=bold
[<box> GreenBox]
[<blob> HideousBlob]
```

Modifiers

```plaintext
dashed
empty
```

Key/value pairs

```plaintext
fill=(any css color)

stroke=(any css color)

align=center
align=left

direction=right
direction=down

visual=actor
visual=class
visual=database
visual=ellipse
visual=end
visual=frame
visual=hidden
visual=input
visual=none
visual=note
visual=package
visual=receiver
visual=rhomb
visual=roundrect
visual=sender
visual=start
visual=table
visual=transceiver
```

Style title and text body with a comma separated list of text modifiers

```plaintext
title=left,italic,bold
body=center,italic,bold
```

Text modifiers

```plaintext
bold
center
italic
left
underline
```
"""
module NomnomlJS

using Artifacts, Librsvg_jll, NodeJS

export Diagram

"""
    Diagram(str::String)

Construct a *nomnoml* diagram from the given source code `str`. To create a
`Diagram` from a file use `read(filename, Diagram)` instead.
"""
struct Diagram
    src::String
end

Base.read(file::AbstractString, ::Type{Diagram}) = Diagram(read(file, String))
Base.write(filename::String, d::Diagram) = open(io -> show(io, mimetype(filename)(), d), filename, "w")

Base.show(io::IO, d::Diagram) = print(io, "$Diagram(...)")

### Keep next section in sync.

const SVG = MIME"image/svg+xml"
const PNG = MIME"image/png"
const PDF = MIME"application/pdf"
const EPS = MIME"application/postscript"

const SUPPORTED_MIMES = Union{SVG,PNG,PDF,EPS}

extension(::SVG) = "svg"
extension(::PNG) = "png"
extension(::PDF) = "pdf"
extension(::EPS) = "eps"
extension(other) = error("unknown extension '$other'.")

function mimetype(filename::AbstractString)
    _, ext = splitext(filename)
    return mimetype(Val{Symbol(ext[2:end])}())
end
mimetype(::Val{:svg}) = SVG
mimetype(::Val{:png}) = PNG
mimetype(::Val{:pdf}) = PDF
mimetype(::Val{:eps}) = EPS
mimetype(::Val{s}) where s = error("unknown format '$s'.")

###

Base.show(io::IO, mime::SUPPORTED_MIMES, d::Diagram) = write(io, converter(d.src, mime))

nomnoml_bin() = joinpath(artifact"nomnoml", "index.js")
nomnoml() = `$(nodejs_cmd()) $(nomnoml_bin())`

function exec(cmd::Cmd, input::IOBuffer)
    output, errors = IOBuffer(), IOBuffer()
    yes = success(pipeline(cmd; stdout = output, stdin = input, stderr = errors))
    return (yes = yes, io = yes ? output : errors)
end
exec(cmd::Cmd, input="") = exec(cmd, IOBuffer(input))

converter(src, ::SVG) = take!(check(exec(nomnoml(), src)).io)
function converter(src, m::SUPPORTED_MIMES)
    Librsvg_jll.rsvg_convert() do bin
        take!(check(exec(`$bin -f $(extension(m))`, converter(src, SVG()))).io)
    end
end
check(result) = result.yes ? result : error(String(take!(result.io)))

end # module


text:
"""
    with_gil(f, [c=true])

Compute `f()` with the GIL enabled.

This may need a `try-finally` block to ensure the GIL is released again. If you know that `f` cannot throw, pass `c=false` to avoid this overhead.
"""
function with_gil(f, c::Bool = true)
    if !CONFIG.isembedded
        f()
    elseif c
        g = C.PyGILState_Ensure()
        try
            f()
        finally
            C.PyGILState_Release(g)
        end
    else
        g = C.PyGILState_Ensure()
        r = f()
        C.PyGILState_Release(g)
        r
    end
end


text:
include("aes-code.jl")

global const BLOCK_BYTES = (Nb * WORDLENGTH)

# Carries out AES in ECB mode on the given blocks using the given key and
# returns the result.
# ECB mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, ECB mode is carried out in decryption
# direction.
# Both the blocks and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned result is a string consisting of only hexadecimal
# characters.
function AESECB(blocks::String, key::String, encrypt::Bool)
	bytes2hex(AESECB(hex2bytes(blocks), hex2bytes(key), encrypt))
end

# Carries out AES in ECB mode on the given blocks using the given key.
# CBC mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, ECB mode is carried out in decryption
# direction.
function AESECB(blocks::Array{UInt8, 1}, key::Array{UInt8, 1}, encrypt::Bool)
	noBlocks = paddedCheck(blocks, key)
	o = Array{UInt8}(undef, length(blocks))

	for i=1:noBlocks
		indices = blockIndices(blocks, i)
		o[indices] = encrypt ? AESEncrypt(blocks[indices], key) : AESDecrypt(blocks[indices], key)
	end

	return o
end

# Carries out AES in CBC mode on the given blocks using the given key and
# returns the result.
# CBC mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, CBC mode is carried out in decryption
# direction.
# Both the blocks and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned result is a string consisting of only hexadecimal
# characters.
function AESCBC(blocks::String, key::String,
	iv::String, encrypt::Bool)
	bytes2hex(AESCBC(hex2bytes(blocks), hex2bytes(key), hex2bytes(iv),
	encrypt))
end

# Carries out AES in CBC mode on the given blocks using the given key.
# CBC mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, CBC mode is carried out in decryption
# direction.
function AESCBC(blocks::Array{UInt8, 1}, key::Array{UInt8, 1},
	iv::Array{UInt8, 1}, encrypt::Bool)
	noBlocks = paddedCheck(blocks, key)
	if length(iv) != BLOCK_BYTES
		error("IV does not have 16 bytes!")
	end
	o = Array{UInt8}(undef, length(blocks))
	prev = iv

	for i=1:noBlocks
		indices = blockIndices(blocks, i)
		if encrypt
			curr = AESEncrypt(xor.(prev , blocks[indices]), key)
			o[indices] = curr
			prev = curr
		else
			# decrypt
			curr = xor.(AESDecrypt(blocks[indices], key) , prev)
			o[indices] = curr
			prev = blocks[indices]
		end
	end

	return o
end

# Carries out AES in CFB mode on the given blocks using the given key and
# returns the result.
# CFB mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, CFB mode is carried out in decryption
# direction.
# Both the blocks and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned result is a string consisting of only hexadecimal
# characters.
function AESCFB(blocks::String, key::String,
	iv::String, encrypt::Bool)
	bytes2hex(AESCFB(hex2bytes(blocks), hex2bytes(key), hex2bytes(iv),
	encrypt))
end

# Carries out AES in CFB mode on the given blocks using the given key.
# CFB mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, CFB mode is carried out in decryption
# direction.
function AESCFB(blocks::Array{UInt8, 1}, key::Array{UInt8, 1},
	iv::Array{UInt8, 1}, encrypt::Bool)
	noBlocks = keyStreamCheck(blocks, key, iv)
	o = Array{UInt8}(undef, length(blocks))
	curr = iv

	for i=1:noBlocks
		indices = blockIndices(blocks, i)
		o[indices] = xor.(AESEncrypt(curr, key)[1:length(indices)] , blocks[indices])
		if encrypt
			curr = o[indices]
		else
			# decrypt
			curr = blocks[indices]
		end
	end

	return o
end

# Carries out AES in OFB mode on the given blocks using the given key and
# returns the result.
# Both the blocks and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned result is a string consisting of only hexadecimal
# characters.
function AESOFB(blocks::String, key::String, iv::String)
	bytes2hex(AESOFB(hex2bytes(blocks), hex2bytes(key), hex2bytes(iv)))
end

# Carries out AES in OFB mode on the given blocks using the given key.
function AESOFB(blocks::Array{UInt8, 1}, key::Array{UInt8, 1},
	iv::Array{UInt8, 1})
	noBlocks = keyStreamCheck(blocks, key, iv)
	o = Array{UInt8}(undef, length(blocks))
	prev = iv

	for i=1:noBlocks
		indices = blockIndices(blocks, i)
		eb = AESEncrypt(prev, key)
		o[indices] = xor.(eb[1:length(indices)] , blocks[indices])
		prev = eb
	end

	return o
end

# Carries out AES in CTR mode on the given blocks using the given key and
# returns the result.
# Both the blocks and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned result is a string consisting of only hexadecimal
# characters. Treats the low eight bytes of the iv array as the little endian
# counter.
function AESCTR(blocks::String, key::String, iv::String)
	bytes2hex(AESCTR(hex2bytes(blocks), hex2bytes(key), hex2bytes(iv)))
end

# Carries out AES in CTR mode on the given blocks using the given key.
# CTR mode is carried out in encryption direction if the parameter
# encrypt is true. Otherwise, CTR mode is carried out in decryption
# direction. Treats the least-significant eight bytes of the iv array
# as the big endian counter.
function AESCTR(blocks::Array{UInt8, 1}, key::Array{UInt8, 1},
	iv::Array{UInt8, 1})
	noBlocks = keyStreamCheck(blocks, key, iv)
	curr = copy(iv)
	o = Array{UInt8}(undef, length(blocks))

	for i=1:noBlocks
		indices = blockIndices(blocks, i)
		eb = AESEncrypt(curr, key)
		o[indices] = xor.(eb[1:length(indices)] , blocks[indices])
		for bi=reverse((length(curr) - 7):length(curr))
			tmp = curr[bi]
			curr[bi] = UInt8(mod(Int(tmp) + 1, 256))
			if curr[bi] > tmp
				# no byte overflow
				break
			end
		end
	end

	return o
end

# Checks whether the parameters are OK for ECB/CBC mode.
# Returns the number of blocks including the last block whose length must
# be as long as the standard block length (16 bytes).
function paddedCheck(blocks::Array{UInt8, 1}, key::Array{UInt8, 1})
	noBlocks = div(length(blocks), BLOCK_BYTES)
	if (noBlocks < 1) || ((noBlocks * BLOCK_BYTES) != length(blocks))
		error("No blocks or length of blocks is not a multplile of ",
		"16!")
	end
	# Check if key is OK
	AESParameters(key)
	return noBlocks
end

# Checks whether the parameters are OK for CFB/OFB/CTR mode.
# Returns the number of blocks including the last block whose length may
# be less than standard block length (16 bytes).
function keyStreamCheck(blocks::Array{UInt8, 1}, key::Array{UInt8, 1},
	iv::Array{UInt8, 1})
	if length(blocks) < 1
		error("blocks must contain at least one byte!")
	end
	noBlocks = div(length(blocks) + BLOCK_BYTES - 1, BLOCK_BYTES)
	# Check if key is OK
	AESParameters(key)
	if length(iv) != BLOCK_BYTES
		error("iv does not have 16 bytes!")
	end
	return noBlocks
end

# Returns the indices of the bytes of a block given a one-based
# block number. The last block may be incomplete.
function blockIndices(blocks::Array{UInt8, 1}, blockNumber::Int)
	@assert(blockNumber >= 1)
	((blockNumber - 1) * BLOCK_BYTES + 1):(min(blockNumber * BLOCK_BYTES, length(blocks)))
end


text:
## plot range tracking error

Gadfly.set_default_plot_size(10cm, 5cm)
pl = Gadfly.plot(dferr,
  x=:x, y=:y, color=:Legend, Geom.line,
  style(line_width=0.5mm, point_size=0.2mm),
  Guide.xlabel("time [s]"),
  Guide.ylabel("range err [m]"),
)

pl |> SVG(joinpath(dirname(@__FILE__),"exports","rangeTrackingError.svg"),10cm,5cm)

##  Plot Ranges

Gadfly.set_default_plot_size(10cm, 5cm)
pl = Gadfly.plot(df,
  x=:x, y=:y, color=:Legend, Geom.line,
  style(line_width=0.5mm, point_size=0.2mm),
  Guide.xlabel("time [s]"),
  Guide.ylabel("range [m]"),
)

pl |> SVG(joinpath(dirname(@__FILE__),"exports","trackingRange.svg"),10cm,5cm)

##  Plot angles

Gadfly.set_default_plot_size(10cm, 5cm)
pl = Gadfly.plot(dfa,
  x=:x, y=:y, color=:Legend, Geom.line,
  style(line_width=0.5mm, point_size=0.2mm),
  Guide.xlabel("time [s]"),
  Guide.ylabel("angle [rad]"),
)

pl |> SVG(joinpath(dirname(@__FILE__),"exports","trackingAngle.svg"),10cm,5cm)

## Plot Cartesian figure

Gadfly.set_default_plot_size(10cm, 5cm)
pl = Gadfly.plot(dfcm,
  x=:x, y=:y, color=:Legend, Geom.path,
  style(line_width=0.5mm, point_size=0.2mm),
  Guide.xlabel("horizontal [m]"),
  Guide.ylabel("vertical [m]"),
)

pl.coord=Gadfly.Coord.Cartesian(xmin=-20,xmax=10, ymin=-60,ymax=20)

# Gadfly.set_default_plot_size(10cm, 8cm)
pl |> SVG(joinpath(dirname(@__FILE__),"exports","cartesian.svg"),10cm,6cm)

##  HEATMAP

n = 100
x = LinRange(-5.0, 60.0, n)
y = LinRange(-pi, pi, n)
z = zeros(n,n)

xy = (x,y) -> sum( [TT[i](collect([x y]'))[1] for i in [1:3;5:10;]] )

for i in 1:n, j in 1:n
  z[i,j] = xy(x[i], y[j])
end

z[1,1] -= 0.15

# scene = AbstractPlotting.surface(x, y, fill(0f0, N, N), color = z, shading = false)
scene = AbstractPlotting.contour(x, y, -z, levels = 0, linewidth = 0, fillrange = true )

swit = 301
lines!(scene, randense[1:(swit-1)], angdense[1:(swit-1)], color = :red, linewidth=2.5, linestyle=:dash)
lines!(scene, randense[(swit+1):end], angdense[(swit+1):end], color = :red, linewidth=2.5, linestyle=:dash)

axis = scene[Axis] # get the axis object from the scene
# axis[:grid][:linecolor] = ((:red, 0.5), (:blue, 0.5))
# axis[:names][:textcolor] = ((:red, 1.0), (:blue, 1.0))
axis[:names][:axisnames] = ("range [m]", "angle [rad]")
scene

##

save(joinpath(dirname(@__FILE__),"exports","plot.png"), scene)

## Look at cylindrical product data

sym = :t0
for sym in [Symbol("t$i") for i in 0:(nposes-1)]

pll, plc = plotLocalProductCylinder(fg, sym, show=false, scale=0.1);

# Gadfly.set_default_plot_size(8cm, 6cm)
pll |> SVG(joinpath(dirname(@__FILE__),"exports","localProductLinear$sym.svg"),8cm, 6cm)
plc |> SVG(joinpath(dirname(@__FILE__),"exports","localProductCirc$sym.svg"),8cm, 6cm)

end


text:
"""
Function to pool the univariate estimators
Arguments:
unibetas::Array{Unibeta, 1} -> Array which contains the univariate estimators which have to be pooled
"""
function pool_unibetas(unibetas::Array{Unibeta, 1})
	pooledbetas = unibetas[1].n .* unibetas[1].unibeta
	sumN = unibetas[1].n
	@simd for i = 2 : length(unibetas)
		@inbounds pooledbetas = pooledbetas .+ unibetas[i].n .* unibetas[i].unibeta
		@inbounds sumN = sumN .+ unibetas[i].n
	end
	pooledbetas = pooledbetas ./ sumN
	return Unibeta(pooledbetas, sumN, unibetas[1].labels)
end
"""
Function to pool the covariance matrices
Arguments:
covarmat::Array{Covarmat, 1} -> Array which contains the covariances which have to be pooled
"""
function pool_covarmats(unibeta::Unibeta, covarmat::Array{Covarmat, 1})
	pooledcovarmat = zeros(size(covarmat[1].covarmat))
	sumN = unibeta.n[1:size(covarmat[1].covarmat)[1]]
	@simd for i = 1 : size(covarmat[1].covarmat)[2]
		@inbounds pooledcovarmat[:,i] = covarmat[1].covarmat[:,i] .* unibeta.n[1:size(covarmat[1].covarmat)[1]]
	end
	@simd for i = 2 : length(covarmat)
		@simd for j = 1 : size(covarmat[1].covarmat)[2]
			@inbounds pooledcovarmat[:,j] = pooledcovarmat[:,j] .+ covarmat[i].covarmat[:,j] .* unibeta.n[1:size(covarmat[1].covarmat)[1]]
		end
		@inbounds sumN = sumN .+ unibeta.n[1:size(covarmat[1].covarmat)[1]]
	end
	@simd for i = 1 : size(covarmat[1].covarmat)[2]
		@inbounds pooledcovarmat[:,i] = pooledcovarmat[:,i] ./ sumN
	end
	return pooledcovarmat
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
using MiniGNN
using Test


text:
# Adapted from https://github.com/terasakisatoshi/sysimage_creator/
import Pkg

Pkg.add(["PackageCompiler", "IJulia"])

using PackageCompiler

sysimage_path = joinpath(@__DIR__, "sysimage.so")

@info "SysImage path: " sysimage_path

# Do not include too many packages because memory on CI machines is limited
PackageCompiler.create_sysimage(
    ["Agents"];
    project=".",
    sysimage_path=sysimage_path,
    cpu_target=PackageCompiler.default_app_cpu_target()
)

nthreads = Sys.CPU_THREADS

using IJulia
kernelpath = IJulia.installkernel("Julia-sys-$(nthreads)-threads",
    "--project=@.", "--sysimage=$(sysimage_path)";
    specname="julia", env=Dict("JULIA_NUM_THREADS" => "$(nthreads)",)
)

@info "IJulia kernal path: " kernelpath

Pkg.rm("PackageCompiler")
Pkg.gc()

max_stars_repo_name: sosiristseng/juliabook-abm

max_issues_repo_name: sosiristseng/juliabook-abm

max_forks_repo_name: sosiristseng/juliabook-abm


text:
module OptimTestProblems

export MultivariateProblems, UnivariateProblems

include("optim_tests/multivariate/MultivariateProblems.jl")
include("optim_tests/univariate/bounded.jl")

# Deprecation stuff
UnconstrainedProblems = OptimTestProblems.MultivariateProblems.UnconstrainedProblems
export UnconstrainedProblems

end # module


text:
module Layers

export Linear, Sequential, forward, backward, parms, ∇parms

include("./tensor_dt.jl")
using ..TensorDT: Tensor

include("./abstract_layers.jl")
import ..AbstractLayers: AbstractLayer, AL, forward, backward, parms, ∇parms, idims, odims

include("./initializations.jl")
using ..Initializations: init_xavier, INIT_FNs

## ======================================================================
## Linear (or Dense) Layer
## ======================================================================
struct Linear <: AbstractLayer
  idim::Integer   ## input dim
  odim::Integer   ## output dim
  ## parameters
  w::Tensor
  b::Tensor
  #
  _type::Symbol
  store::Dict # {Any. Any}

  ## Assume type is DT == Float64
  function Linear(idim::Integer, odim::Integer;
      DT=Float64, init_fn=init_xavier, seed=42, kwargs...)
    ##
    @assert 1 ≤ idim
    @assert 1 ≤ odim
    @assert init_fn ∈ INIT_FNs
    ##
    w = init_fn((odim, idim); seed, DT, kwargs...)
    b = zeros(DT, (odim, 1))
    new(idim, odim, w, b, :dense, Dict())
  end
end

# gen_key(self::Linear, prefix::Symbol) = Symbol(prefix, self._id[1])

function forward(self::Linear, input::Tensor)::Tensor
  self.store[:input] = input  ## Storing the input for backward pass
  self.w * input + self.b
end

function backward(self::Linear, ∇p::Tensor)::Tensor
  ## as each bᵢ is added to output oᵢ, ∇b is the same as output ∇
  self.store[:∇b] = ∇p
  self.store[:∇w] = ∇p * self.store[:input]'
  r = sum.(self.w' * ∇p)
  @assert(size(self.w) == size(self.store[:∇w]),
          "w and ∇w should have same shape: $(size(self.w)) == $(size(self.store[:∇w]))")
  @assert(size(self.b) == size(self.store[:∇b]),
          "b and ∇b should have same shape: $(size(self.b)) == $(size(self.store[:∇b]))")
  r
end

parms(self::Linear) = [self.w, self.b]
∇parms(self::Linear) = [self.store[:∇w], self.store[:∇b]]

## Extra
idims(self::Linear) = self.idim
odims(self::Linear) = self.odim

## ======================================================================
## Neural Network as a sequence of Layers
## ======================================================================

const CHECK_DIMS_EXCLUSION = [:activation, :regularization] #:embedding]

struct Sequential <: AbstractLayer
  layers::Vector{AbstractLayer}
  _type::Symbol

  function Sequential(layers::Vector{AbstractLayer})
    @assert length(layers) > 0
    check_dim_layer(layers)
    #
    new(layers, :sequential)
  end
end

function forward(self::Sequential, input::Tensor)::Tensor
  for (ix, l) ∈ enumerate(self.layers)
    input = forward(l, input)
  end
  input
end

function backward(self::Sequential, ∇p::Tensor)::Tensor
  for (ix, l) ∈ enumerate(reverse(self.layers))
    ∇p = backward(l, ∇p)
  end
  ∇p
end

parms(self::Sequential) = [p for l ∈ self.layers for p ∈ parms(l)]
∇parms(self::Sequential) = [∇p for l ∈ self.layers for ∇p ∈ ∇parms(l)]

idims(self::Sequential) = nothing
odims(self::Sequential) = nothing

##
## Internal helpers
##

function check_dim_layer(layers::Vector{AbstractLayer})
  """
  output of prev. layer == input of curr layer
  """
  @assert length(layers) > 0
  pl = layers[1]._type ∈ CHECK_DIMS_EXCLUSION ? layers[1] : nothing

  for cl ∈ layers[2:end]
    cl._type ∈ CHECK_DIMS_EXCLUSION && continue
    isnothing(pl) && (pl = cl; continue)

    odims(pl) ≠ idims(cl) &&
      throw(ArgumentError("incompatible shape between pl.odim; $(odims(pl)) vs cl.idim: $(idims(cl))"))
    pl = cl
  end
end

end ## Module

max_stars_repo_path: Data Science From Scratch/src/layers.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/src/layers.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/src/layers.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
using AMDB
using GZip

# data_dir = "/Data/AMDB/RawTables/"
data_dir = expanduser("~/research/AMDB/data/")

"""
Parse lines from `io`, and keep the `id` if the spell status is in
`status_types`. For the AMP database. Returns a `Set` of ids.
"""
function filtered_ids(io::IO, status_types, ids = Set{Int}())
    counter = 0
    readlines(io) do line
        (id_field, empnum, status_field, other_fields) = split(line, ';'; limit=4)
        id = parse(Int, id_field)
        status = tryparse(AMP.spell_dict, status_field)
        if isnull(status)
            warn("Could not parse line:\n$line\n")
        else
            if get(status) ∈ status_types
                push!(ids, id)
            end
        end
        counter += 1
        if mod(counter, 1000000) == 0
            print(".")
        end
    end
    println()
    ids
end

# keep these for people giving birth
keep_status = Set([AMP.maternity_active,
                   AMP.maternity_inactive,
                   AMP.parental_leave_inactive,
                   AMP.parental_leave_active])

maternity_ids =
    process_compressed_database(joinpath(data_dir, "hv_epi_uni_roh.csv.gz")) do io
        filtered_ids(io, keep_status)
    end


text:
using ACE1
using IPFitting

al = IPFitting.Data.read_xyz(@__DIR__() * "/TiAl_tutorial_DB.xyz", energy_key="energy", force_key="force", virial_key="virial")[1:10:end];

r0 = 0.5*(rnn(:Ti) + rnn(:Al))

ACE_B = ACE1.Utils.rpi_basis(species = [:Ti, :Al],
                              N = 3,
                              r0 = r0,
                              rin = 0.6 * r0,
                              rcut = 5.5,
                              maxdeg = 6)

Bpair = pair_basis(species = [:Ti, :Al],
      r0 = r0,
      maxdeg = 6,
      rcut = 7.0,
      pcut = 1,
      pin = 0)

B = JuLIP.MLIPs.IPSuperBasis([Bpair, ACE_B]);

dB = LsqDB("", B, al)

Vref = OneBody(:Ti => -1586.0195, :Al => -105.5954)

weights = Dict(
        "FLD_TiAl" => Dict("E" => 30.0, "F" => 1.0 , "V" => 1.0 ),
        "TiAl_T5000" => Dict("E" => 5.0, "F" => 1.0 , "V" => 1.0 ))

solver_type = :ard
laplace_precon = true

if solver_type == :lsqr
	solver = Dict(
        	"solver" => :lsqr,
        	"lsqr_damp" => 1e-2,
        	"lsqr_atol" => 1e-6)
elseif solver_type == :rrqr
	solver = Dict(
        	"solver" => :rrqr,
        	"rrqr_tol" => 1e-5)
elseif solver_type == :brr
	solver = Dict(
        	"solver" => :brr,
		"brr_tol" => 1e-3)
elseif solver_type == :ard
	solver= Dict(
         	"solver" => :ard,
         	"ard_tol" => 1e-3,
         	"ard_threshold_lambda" => 10000)
end

if laplace_precon
	using LinearAlgebra
	rlap_scal = 3.0
	P = Diagonal(vcat(ACE1.scaling.(dB.basis.BB, rlap_scal)...))
	solver["P"] = P
end

IP, lsqinfo = lsqfit(dB, solver=solver, weights=weights, Vref=Vref, error_table = true);

rmse_table(lsqinfo["errors"])

save_dict("./TiAl_tutorial_pot.json", Dict("IP" => write_dict(IP), "info" => lsqinfo))

ACE1.ExportMulti.export_ACE("./TiAl_tutorial_pot.yace", IP)


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("../module/CustomStencil.jl")
include("misc.jl")
include("cpu_stencils.jl")

## Stencil Definition
coefs = Float64.([1;0.5;0.25])
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 2, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def, m_step=2)

## Input Data size Definition
radius = 2
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 4
## GPU Implementation
data = CreateData(dx,dy,dz)

dout = ApplyStencil(st_inst, data, t_steps)
gpu_out = Array(dout)
gpu_out = gpu_out[3:38, 3:38, :]
## CPU Implementation
cpu_out =  cpu_star_stencil(data, radius, coefs, t_steps)

## Check Result
mean_er = sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out)))
@assert(sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out))) < 0.001)

using Plots

shift = 0
Int(round(dz/2))+shift
heatmap(cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(gpu_out[:,:,Int(round(dz/2))+shift])

d2 = ApplyMultiGPU(4, st_inst, t_steps, data)
gpu_out  = d2
heatmap(gpu_out[:,:,Int(round(dz/2))+shift])
heatmap(gpu_out[:,:,Int(round(dz/2))+shift]-cpu_out[:,:,Int(round(dz/2))+shift])
## Visual Check
using Plots
shift = 0
Int(round(dz/2))+shift
heatmap(cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(gpu_out[:,:,Int(round(dz/2))+shift])

for i = 1:size(original_data,3)

    mean_er = sum(abs.(cpu_out[:,:,i] .- gpu_out[:,:,i]))./(reduce(*, size(gpu_out[:,:,1])))
    if mean_er > 0.01
        println(i, " ", mean_er)
    end
end

st_def = @def_stencil_expression c[0]*D[x,y,z] + c[1]*D[x+1,y,z]
c = Float64.([1; 0.5])
st_inst = CreateStencilDefinition(st_def,c)
st = NewStencilInstance(st_inst, m_step=3)
st_sym, r = CombineStencils(st_inst)
new_sym = st_sym
for i = 2:3
    global new_sym
    new_sym = CombineTimeSteps(st_sym, 1, new_sym, i)
end
new_sym


text:
function LineAxis(parent::Scene; kwargs...)

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

    decorations = Dict{Symbol, Any}()

    @extract attrs (endpoints, limits, flipped, ticksize, tickwidth,
        tickcolor, tickalign, ticks, ticklabelalign, ticklabelrotation, ticksvisible,
        ticklabelspace, ticklabelpad, labelpadding,
        ticklabelsize, ticklabelsvisible, spinewidth, spinecolor, label, labelsize, labelcolor,
        labelfont, ticklabelfont,
        labelvisible, spinevisible, trimspine, flip_vertical_label)

    pos_extents_horizontal = lift(endpoints) do endpoints
        if endpoints[1][2] == endpoints[2][2]
            horizontal = true
            extents = (endpoints[1][1], endpoints[2][1])
            position = endpoints[1][2]
            return (position, extents, horizontal)
        elseif endpoints[1][1] == endpoints[2][1]
            horizontal = false
            extents = (endpoints[1][2], endpoints[2][2])
            position = endpoints[1][1]
            return (position, extents, horizontal)
        else
            error("Axis endpoints $(endpoints[1]) and $(endpoints[2]) are neither on a horizontal nor vertical line")
        end
    end

    ticksnode = Node(Point2f0[])
    ticklines = linesegments!(
        parent, ticksnode, linewidth = tickwidth, color = tickcolor,
        show_axis = false, visible = ticksvisible
    )[end]
    decorations[:ticklines] = ticklines

    ticklabelannosnode = Node{Vector{Tuple{String, Point2f0}}}([("temp", Point2f0(0, 0))])
    ticklabels = annotations!(
        parent,
        ticklabelannosnode,
        align = ticklabelalign,
        rotation = ticklabelrotation,
        textsize = ticklabelsize,
        font = ticklabelfont,
        show_axis = false,
        visible = ticklabelsvisible)[end]

    ticklabel_ideal_space = lift(ticklabelannosnode, ticklabelalign, ticklabelrotation, ticklabelfont, ticklabelsvisible, typ=Float32) do args...
        maxwidth = if pos_extents_horizontal[][3]
                # height
                ticklabelsvisible[] ? height(FRect2D(boundingbox(ticklabels))) : 0f0
            else
                # width
                ticklabelsvisible[] ? width(FRect2D(boundingbox(ticklabels))) : 0f0
        end
    end

    attrs[:actual_ticklabelspace] = 0f0
    actual_ticklabelspace = attrs[:actual_ticklabelspace]

    onany(ticklabel_ideal_space, ticklabelspace) do idealspace, space
        s = if space == AbstractPlotting.automatic
            idealspace
        else
            space
        end
        if s != actual_ticklabelspace[]
            actual_ticklabelspace[] = s
        end
    end

    decorations[:ticklabels] = ticklabels

    tickspace = lift(ticksvisible, ticksize, tickalign) do ticksvisible,
            ticksize, tickalign

        ticksvisible ? max(0f0, ticksize * (1f0 - tickalign)) : 0f0
    end

    labelgap = lift(spinewidth, tickspace, ticklabelsvisible, actual_ticklabelspace,
        ticklabelpad, labelpadding) do spinewidth, tickspace, ticklabelsvisible,
            actual_ticklabelspace, ticklabelpad, labelpadding

        spinewidth + tickspace +
            (ticklabelsvisible ? actual_ticklabelspace + ticklabelpad : 0f0) +
            labelpadding
    end

    labelpos = lift(pos_extents_horizontal, flipped, labelgap) do (position, extents, horizontal), flipped, labelgap

        # fullgap = tickspace[] + labelgap

        middle = extents[1] + 0.5f0 * (extents[2] - extents[1])

        x_or_y = flipped ? position + labelgap : position - labelgap

        if horizontal
            Point2(middle, x_or_y)
        else
            Point2(x_or_y, middle)
        end
    end

    labelalign = lift(pos_extents_horizontal, flipped, flip_vertical_label) do (position, extents, horizontal), flipped, flip_vertical_label
        if horizontal
            (:center, flipped ? :bottom : :top)
        else
            (:center, if flipped
                    flip_vertical_label ? :bottom : :top
                else
                    flip_vertical_label ? :top : :bottom
                end
            )
        end
    end

    labelrotation = lift(pos_extents_horizontal, flip_vertical_label) do (position, extents, horizontal), flip_vertical_label
        if horizontal
            0f0
        else
            if flip_vertical_label
                Float32(-0.5pi)
            else
                Float32(0.5pi)
            end
        end
    end

    labeltext = text!(
        parent, label, textsize = labelsize, color = labelcolor,
        position = labelpos, show_axis = false, visible = labelvisible,
        align = labelalign, rotation = labelrotation, font = labelfont,
    )[end]

    decorations[:labeltext] = labeltext

    tickvalues = Node(Float32[])

    tickvalues = lift(pos_extents_horizontal, limits, ticks) do (position, extents, horizontal),
            limits, ticks
        px_width = extents[2] - extents[1]
        compute_tick_values(ticks, limits..., px_width)
    end

    tickpositions = Node(Point2f0[])
    tickstrings = Node(String[])

    onany(tickvalues) do tickvalues
        # limoy = limits[].origin[2]
        # limh = limits[].widths[2]
        # px_ox = pixelarea(scene)[].origin[1]
        # px_oy = pixelarea(scene)[].origin[2]
        # px_w = pixelarea(scene)[].widths[1]
        # px_h = pixelarea(scene)[].widths[2]

        position, extents, horizontal = pos_extents_horizontal[]

        px_o = extents[1]
        px_width = extents[2] - extents[1]

        lim_o = limits[][1]
        lim_w = limits[][2] - limits[][1]

        tick_fractions = (tickvalues .- lim_o) ./ lim_w
        tick_scenecoords = px_o .+ px_width .* tick_fractions

        tickpos = if horizontal
            [Point(x, position) for x in tick_scenecoords]
        else
            [Point(position, y) for y in tick_scenecoords]
        end

        # now trigger updates
        tickpositions[] = tickpos

        tickstrings[] = get_tick_labels(ticks[], tickvalues)
    end

    onany(tickstrings, labelgap, flipped) do tickstrings, labelgap, flipped
        # tickspace is always updated before labelgap
        # tickpositions are always updated before tickstrings
        # so we don't need to lift those

        position, extents, horizontal = pos_extents_horizontal[]

        nticks = length(tickvalues[])

        ticklabelgap = spinewidth[] + tickspace[] + ticklabelpad[]

        shift = if horizontal
            Point2f0(0f0, flipped ? ticklabelgap : -ticklabelgap)
        else
            Point2f0(flipped ? ticklabelgap : -ticklabelgap, 0f0)
        end

        ticklabelpositions = tickpositions[] .+ Ref(shift)
        ticklabelannosnode[] = collect(zip(tickstrings, ticklabelpositions))
    end

    onany(tickpositions, tickalign, ticksize, spinewidth) do tickpositions,
            tickalign, ticksize, spinewidth

        position, extents, horizontal = pos_extents_horizontal[]

        if horizontal
            tickstarts = [tp + (flipped[] ? -1f0 : 1f0) * Point2f0(0f0, tickalign * ticksize - 0.5f0 * spinewidth) for tp in tickpositions]
            tickends = [t + (flipped[] ? -1f0 : 1f0) * Point2f0(0f0, -ticksize) for t in tickstarts]
            ticksnode[] = interleave_vectors(tickstarts, tickends)
        else
            tickstarts = [tp + (flipped[] ? -1f0 : 1f0) * Point2f0(tickalign * ticksize - 0.5f0 * spinewidth, 0f0) for tp in tickpositions]
            tickends = [t + (flipped[] ? -1f0 : 1f0) * Point2f0(-ticksize, 0f0) for t in tickstarts]
            ticksnode[] = interleave_vectors(tickstarts, tickends)
        end
    end

    linepoints = lift(pos_extents_horizontal, flipped, spinewidth, trimspine, tickpositions, tickwidth) do (position, extents, horizontal),
            flipped, sw, trimspine, tickpositions, tickwidth

        if !trimspine
            if horizontal
                y = position + (flipped ? 1f0 : -1f0) * 0.5f0 * sw
                p1 = Point2f0(extents[1] - sw, y)
                p2 = Point2(extents[2] + sw, y)
                [p1, p2]
            else
                x = position + (flipped ? 1f0 : -1f0) * 0.5f0 * sw
                p1 = Point2f0(x, extents[1] - sw)
                p2 = Point2f0(x, extents[2] + sw)
                [p1, p2]
            end
        else
            [tickpositions[1], tickpositions[end]] .+ [
                (horizontal ? Point2f0(-0.5f0 * tickwidth, 0) : Point2f0(0, -0.5f0 * tickwidth)),
                (horizontal ? Point2f0(0.5f0 * tickwidth, 0) : Point2f0(0, 0.5f0 * tickwidth)),
            ]
        end
    end

    decorations[:axisline] = lines!(parent, linepoints, linewidth = spinewidth, visible = spinevisible,
        color = spinecolor, raw = true)[end]

    protrusion = lift(ticksvisible, label, labelvisible, labelpadding, labelsize, tickalign, spinewidth,
            tickspace, ticklabelsvisible, actual_ticklabelspace, ticklabelpad, labelfont, ticklabelfont) do ticksvisible,
            label, labelvisible, labelpadding, labelsize, tickalign, spinewidth, tickspace, ticklabelsvisible,
            actual_ticklabelspace, ticklabelpad, labelfont, ticklabelfont

        position, extents, horizontal = pos_extents_horizontal[]

        real_labelsize = if iswhitespace(label)
            0f0
        else
            horizontal ? boundingbox(labeltext).widths[2] : boundingbox(labeltext).widths[1]
        end

        labelspace = (labelvisible && !iswhitespace(label)) ? real_labelsize + labelpadding : 0f0
        spinespace = spinewidth
        # tickspace = ticksvisible ? max(0f0, xticksize * (1f0 - xtickalign)) : 0f0
        ticklabelgap = ticklabelsvisible ? actual_ticklabelspace + ticklabelpad : 0f0

        together = spinespace + tickspace + ticklabelgap + labelspace
    end

    # extents[] = extents[] # trigger

    LineAxis(parent, protrusion, attrs, decorations, tickpositions, tickvalues, tickstrings)
end

function tight_ticklabel_spacing!(la::LineAxis)

    horizontal = if la.attributes.endpoints[][1][2] == la.attributes.endpoints[][2][2]
        true
    elseif la.attributes.endpoints[][1][1] == la.attributes.endpoints[][2][1]
        false
    else
        error("endpoints not on a horizontal or vertical line")
    end

    tls = la.decorations[:ticklabels]
    maxwidth = if horizontal
            # height
            tls.visible[] ? height(FRect2D(boundingbox(tls))) : 0f0
        else
            # width
            tls.visible[] ? width(FRect2D(boundingbox(tls))) : 0f0
    end
    la.attributes.ticklabelspace = maxwidth
end

function iswhitespace(str)
    match(r"^\s+$", str) !== nothing
end

compute_tick_values(ct::CustomTicks, args...) = ct.f_tickvalues(args...)
get_tick_labels(ct::CustomTicks, values) = ct.f_ticklabels(values)


text:
@testset "distances.jl" begin
    p1 = Point(30.0, 100.0)
    p2 = Point(50.0, 210.0)
    p3 = Point(40.0, 180.0)

    ls = Arc(p1,p2)
    lss = Arcs([p1, p2, p3])

    @test (angular_distance(p1, p2)) ≈ 8773e3/6371008.7714 atol = 1e3
    @test (angular_length(ls)) ≈ 8773e3/6371008.7714 atol = 1e3
    @test (angular_distance(p1, p1)) == 0.0
    @test (angular_length(lss)) ≈ (8773e3+2585e3)/6371008.7714 atol = 1e3

    # Example http://edwilliams.org/avform.htm Cross track error
    pLAX = Point(33.95,-118.4)
    pJFK = Point(40.63333, -73.78333)
    pD = Point(34.5, -116.5)

    dist_LAX2D = angular_distance(pLAX, pD)
    @test dist_LAX2D ≈ rad2deg(0.02905) atol = 0.1
    bearing_LAX2D = azimuth(pLAX, pD)
    @test bearing_LAX2D ≈ 70.17 atol = 0.01
    bearing_LAX2JFK = azimuth(pLAX, pJFK)

    @test angular_distance(dist_LAX2D, bearing_LAX2D, bearing_LAX2JFK) ≈
    rad2deg(0.00216747) atol = 0.0001
    @test angular_distance(pD, pLAX, bearing_LAX2JFK) ≈
    rad2deg(0.00216747) atol = 0.0001
    @test angular_distance(pD, Line(pLAX, bearing_LAX2JFK)) ≈
    rad2deg(0.00216747) atol = 0.0001

    @test along_line_angular_distance(rad2deg(0.02905), rad2deg(0.00216747)) ≈
    rad2deg(0.0289691) atol = 0.00001
    @test along_line_angular_distance(pD, Line(pLAX, bearing_LAX2JFK)) ≈
    rad2deg(0.0289691) atol = 0.00001

    px1 = Point(-10.0, -20.0)
    px2 = Point(10.0, -20.0)
    px3a = Point(0.0, -10.0)
    @test angular_distance(px3a, px1, px2) ≈ 10.0 atol = 0.1
    px3b = Point(15.0, -20.0)
    @test angular_distance(px3b, px1, px2) ≈ 5.0 atol = 0.1
    px3c = Point(-17.0, -20.0)
    @test angular_distance(px3c, px1, px2) ≈ 7.0 atol = 0.1
    @test angular_distance(px3c, Arc(px1, px2)) ≈ 7.0 atol = 0.1

    px4 = Point(0.0, 0.0)
    px5 = Point(10.0, 0.0)
    arc = Arc(px4, px5)
    @test (angular_distance(Point(-5.0, 0.0), arc)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(-5.0, 5.0), arc)) ≈ 785.8*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(0.0, 5.0), arc)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(5.0, 5.0), arc)) ≈ 553.9*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(10.0, 5.0), arc)) ≈ 547.5*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(15.0, 0.0), arc)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(15.0, -5.0), arc)) ≈ 776.9*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(10.0, -5.0), arc)) ≈ 547.5*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(5.0, -5.0), arc)) ≈ 553.9*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(0.0, -5.0), arc)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(5.0, 0.0), arc)) ≈ 0.0*1000/6371008.7714 atol = 100.0

    @test (angular_distance(Point(-5.0, 0.0), px4)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(-5.0, 5.0), px4)) ≈ 785.8*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(0.0, 5.0), px4)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(5.0, 5.0), arc)) ≈ 553.9*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(10.0, 5.0), px5)) ≈ 547.5*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(15.0, 0.0), px5)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(15.0, -5.0), px5)) ≈ 776.9*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(10.0, -5.0), px5)) ≈ 547.5*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(5.0, -5.0), arc)) ≈ 553.9*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(0.0, -5.0), px4)) ≈ 556.0*1000/6371008.7714 atol = 100.0
    @test (angular_distance(Point(5.0, 0.0), arc)) ≈ 0.0*1000/6371008.7714 atol = 100.0

    @test (angular_length(Polygon(p1, [p1, p2, p3, p1]))) ≈ (8773e3+
    2585e3)/6371008.7714+(angular_distance(p1, p3)) atol = 1e3

end


text:
const β_min, β_step = 0.01, 0.01

function x_upper(β)
    q = -Inf
    x = 0.01
    while cdf(GeneralizedKolmogorov(β), x; precomputed_interpolation = false) < 0.9999
        x *= 2
    end
    return x
end

function make_cdfs()
    β_grid = β_min:β_step:1.0
    δ = 0.001

    cdf_inter_results = []
    @showprogress for β in β_grid
        upper_value = x_upper(β)
        x_grid = 0:δ:upper_value
        cdf_vals = vcat([cdf(GeneralizedKolmogorov(β), x; precomputed_interpolation = false) for x in x_grid], 1.0)
        x_grid = union(x_grid, 2*upper_value)
        push!(cdf_inter_results, LinearInterpolation(x_grid, cdf_vals, extrapolation_bc=Flat()))
    end
    return cdf_inter_results
end


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("Utils.jl")
using LinearAlgebra
using RecursiveArrayTools
using Distributions
using Distances
using Random
using Profile
using JLD
using MAT
using BlockArrays
Random.seed!(123)

NTrain = 2000
NTest = 5
N = NTest + NTrain
d =  2
x = randn(d, N)

l = 1.0
σ = Diagonal(2.0 * ones(N))
cov = r -> exp(-r/l)

ρ = 4.0
λ = 1.3

# reordering points and forming skeletons
@time colptr, rowval, P, revP, distances = sortSparseRev(x, ρ, 1)
x = x[:, P]
skeletons = construct_skeletons(colptr, rowval, distances, λ)

# Creating
Θ = cov.(pairwise(Euclidean(), x, dims=2))
Θ = Θ + σ

L = assembleL(skeletons, x, cov)
U = sparse(triu(L * L') + min.(inv(σ), 1e30))
icholU_high_level!(U)

ΘApprox = inv(Matrix(L' * L)) * U' * U * σ

@show KL(Θ, ΘApprox)
@show KL(ΘApprox, Θ)


text:
using Calculus
import Base:+,-,*,/,vec

function vec(ex::Expr)
	nex = :([])
	nc = length(ex.args)
	nr = length(ex.args[1].args)
	@assert ex.head==:vcat

	for r in ex.args
		@assert r.head == :row && length(r.args) == nr
	end

    for ri = 1:nr
    	for ci = 1:nc
			push!(nex.args,ex.args[ci].args[ri])
		end
	end
	return nex
end

+{T<:Union{Expr,Symbol,Number}}(x::T) = x

function +{T<:Union{Expr,Symbol,Number},S<:Union{Expr,Symbol,Number}}(x1::T,x2::S)
    :(+($x1,$x2))
end

function +{T<:Union{Expr,Symbol},S<:Float64}(x1::T,x2::Vector{S})
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]+x1))
    end
    return out
end

function +{T<:Union{Expr,Symbol},S<:Float64}(x2::Vector{S},x1::T)
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]+x1))
    end
    return out
end

-{T<:Union{Expr,Symbol,Number}}(x::T) = -1*x

function -{T<:Union{Expr,Symbol,Number},S<:Union{Expr,Symbol,Number}}(x1::T,x2::S)
    :(-($x1,$x2))
end

function -{T<:Union{Expr,Symbol},S<:Float64}(x1::T,x2::Vector{S})
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]-x1))
    end
    return out
end

function -{T<:Union{Expr,Symbol},S<:Float64}(x2::Vector{S},x1::T)
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]-x1))
    end
    return out
end

function *{T<:Union{Expr,Symbol,Number},S<:Union{Expr,Symbol,Number}}(x1::T,x2::S)
    :(*($x1,$x2))
end

function *{T<:Union{Expr,Symbol},S<:Float64}(x1::T,x2::Vector{S})
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]*x1))
    end
    return out
end

function *{T<:Union{Expr,Symbol},S<:Float64}(x2::Vector{S},x1::T)
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]*x1))
    end
    return out
end

function /{T<:Union{Expr,Symbol,Number},S<:Union{Expr,Symbol,Number}}(x1::T,x2::S)
    :(/($x1,$x2))
end

function /{T<:Union{Expr,Symbol},S<:Float64}(x1::T,x2::Vector{S})
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x1/x2[i]))
    end
    return out
end

function /{T<:Union{Expr,Symbol},S<:Float64}(x2::Vector{S},x1::T)
    out = :([])
    for i = 1:length(x2)
        push!(out.args,simplify(x2[i]/x1))
    end
    return out
end


text:
using Revise,ApproxOperator

elements, nodes = importmsh("./msh/cantilever.msh")
np = length(nodes[:x])

set𝓖!(elements["Ω"],:QuadGI2)
set𝓖!(elements["Γᵗ"],:SegGI2)
set𝓖!(elements["Γᵍ"],:SegGI2)

P = 1000.0
Ē = 3e6
ν̄ = 0.499999999
E = Ē/(1.0-ν̄^2)
ν = ν̄/(1.0-ν̄)
L = 48.0
D = 12.0
I = D^3/12
EI = E*I
t₁(x,y,z) = 0.0
t₂(x,y,z) = P/2/I*(D^2/4-y^2)
g₁(x,y,z) = -P*y/6/EI*((6*L-3x)*x + (2+ν)*(y^2-D^2/4))
g₂(x,y,z) = P/6/EI*(3*ν*y^2*(L-x) + (4+5*ν)*D^2*x/4 + (3*L-x)*x^2)
prescribe!(elements["Γᵗ"],:t₁,t₁)
prescribe!(elements["Γᵗ"],:t₂,t₂)
prescribe!(elements["Γᵍ"],:g₁,g₁)
prescribe!(elements["Γᵍ"],:g₂,g₂)
prescribe!(elements["Γᵍ"],:n₁₁,(x,y,z)->1.0)
prescribe!(elements["Γᵍ"],:n₂₂,(x,y,z)->1.0)

op_Ω = Operator(:∫∫εᵢⱼσᵢⱼdxdy,:E=>E,:ν=>ν)
op_Ωᵛ = Operator(:∫∫εᵛᵢⱼσᵛᵢⱼdxdy,:E=>Ē,:ν=>ν̄)
op_Ωᵈ = Operator(:∫∫εᵈᵢⱼσᵈᵢⱼdxdy,:E=>Ē,:ν=>ν̄)
op_Γᵗ = Operator(:∫vᵢtᵢds)
op_Γᵍ = Operator(:∫vᵢgᵢds,:α=>1e7*Ē)

k = zeros(2*np,2*np)
f = zeros(2*np)

# op_Ω(elements["Ω"],k)
op_Ωᵈ(elements["Ω"],k)
set𝓖!(elements["Ω"],:QuadGI1)
op_Ωᵛ(elements["Ω"],k)
op_Γᵗ(elements["Γᵗ"],f)
op_Γᵍ(elements["Γᵍ"],k,f)

d = k\f

nodes[:d₁] = d[1:2:2*np-1]
nodes[:d₂] = d[2:2:2*np]

set𝓖!(elements["Ω"],:TriGI3)
prescribe!(elements["Ω"],:u,g₁)
prescribe!(elements["Ω"],:v,g₂)
prescribe!(elements["Ω"],:∂u∂x,(x,y,z)->-P/EI*(L-x)*y)
prescribe!(elements["Ω"],:∂u∂y,(x,y,z)->-P/6/EI*((6*L-3*x)*x + (2+ν)*(3*y^2-D^2/4)))
prescribe!(elements["Ω"],:∂v∂x,(x,y,z)->P/6/EI*((6*L-3*x)*x - 3*ν*y^2 + (4+5*ν)*D^2/4))
prescribe!(elements["Ω"],:∂v∂y,(x,y,z)->P/EI*(L-x)*y*ν)
op = Operator(:Hₑ_PlaneStress,:E=>E,:ν=>ν)
h1,l2 = op(elements["Ω"])

l2 = log10(l2)
h1 = log10(h1)

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
# faststreamingmodetest.jl

using PicoScope
using Base.Test

print("testing ps3000_open_unit...")
handle = ps3000_open_unit()
@test handle != 0
# make sure scope closes even if we have errors
close_unit_atexit() = ps3000_close_unit(handle)
atexit(close_unit_atexit)
println("done")

print("testing ps3000_set_channel...")
enabled = true
dc_coupling = true
ps3000_set_channel(handle, PS3000Channel.A, enabled, dc_coupling, PS3000Range.V1)
println("done")

print("testing ps3000_set_trigger...")
threshold = 0
delay = -50
autotriggerms = 100
ps3000_set_trigger(handle, PS3000Source.CHANNEL_A, threshold,
                    PS3000Direction.RISING, delay, autotriggerms)
println("done")

println("testing fast streaming mode...")
fast_streaming_buffer = zeros(Int16,(8,10000))
function get_overview_buffers(overviewbuffers::Ptr{Ptr{Int16}},
                              overflow::Int16,
                              triggeredat::UInt32,
                              triggered::Int16,
                              auto_stop::Int16,
                              nvalues::UInt32)
  fast_streaming_buffer = unsafe_load(overviewbuffers)
  return nothing
end
const get_overview_buffers_c =
  cfunction(get_overview_buffers, Void,
  (Ptr{Ptr{Int16}},Int16,UInt32,Int16,Int16,UInt32))
sample_interval = 10
time_units = PS3000TimeUnits.MS  # miliseconds
max_samples = 1000
auto_stop = true
no_of_samples_to_aggergate = 1
overview_buffer_size = 15000
print("  call ps3000_runStreaming_ns...")
ps3000_run_streaming_ns(handle, sample_interval, time_units, max_samples,
                        auto_stop, no_of_samples_to_aggergate, overview_buffer_size)
println("return")
print("  call ps3000_get_streaming_last_values...")
number_of_data_values_per_channel =
  ps3000_get_streaming_last_values(handle, get_overview_buffers_c)
println("return")
print("  call ps3000_overview_buffer_status...")
  @test ~ps3000_overview_buffer_status(handle)
println("  return")

print("  call ps3000_stop...")
ps3000_stop(handle)
println("return")
println("done fast streaming mode")

print("closing scope...")
ps3000_close_unit(handle)
println("done")


text:
# Autogenerated wrapper script for tblis_jll for x86_64-w64-mingw32-cxx11
export tblis, tci

using Hwloc_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("tblis")
JLLWrappers.@declare_library_product(tblis, "libtblis-0.dll")
JLLWrappers.@declare_library_product(tci, "libtci-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Hwloc_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        tblis,
        "bin\\libtblis-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        tci,
        "bin\\libtci-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
fetchvcat(J::judiVector) = J
fetchvcat(x) = vcat(x...)

function JUDI.time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, dm, srcnum::UnitRange{Int64}, op::Char, mode::Int64, options)

    # Broadcast common parameters
    _model = @bcast model
    _dm = isnothing(dm) ? dm : @bcast dm
    opts = make_opt([model, dm, op, srcGeometry])
    iter = make_parts(srcnum)
    # Run on azure
    results = @batchexec pmap(sx -> time_modeling_azure(_model, subsample(srcGeometry, sx), subsample(srcData, sx),
                                                        subsample(recGeometry, sx), subsample(recData, sx), _dm,
                                                        op, mode, subsample(options, sx)), iter) opts
    # Gather results
    if op=='F' || (op=='J' && mode==1)
        temp = fetch(results)
        argout1 = fetchvcat(temp)
    elseif op=='J' && mode==-1
        argout1 = fetchreduce(results; op=+, remote=true)
    else
        error("operation no defined")
    end
    return argout1
end

# FWI
function JUDI.fwi_objective(model::Model, source::judiVector, dObs::judiVector; options=JUDI.Options())
    # fwi_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.
    # Broadcast common parameters
    _model = @bcast model
    opts = make_opt([model, dObs, source])
    iter = make_parts(1:dObs.nsrc)
    results = @batchexec pmap(j -> fwi_objective_azure(_model, source[j], dObs[j], subsample(options, j)), iter) opts

    # Collect and reduce gradients
    obj, gradient = fetchreduce(results; op=+, remote=true)
    # first value corresponds to function value, the rest to the gradient
    return obj, gradient
end

# lsrtm

function JUDI.lsrtm_objective(model::Model, source::judiVector, dObs::judiVector, dm; options=JUDI.Options(), nlind=false)
    # lsrtm_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.

    # Broadcast common parameters
    _model = @bcast model
    _dm = isnothing(dm) ? dm : @bcast dm

    opts = make_opt([model, dObs, source, dm])
    iter = make_parts(1:dObs.nsrc)
    results = @batchexec pmap(j -> lsrtm_objective_azure(_model, source[j], dObs[j], _dm, subsample(options, j); nlind=nlind), iter) opts
    # Collect and reduce gradients
    obj, gradient = fetchreduce(results; op=+, remote=false)
    delete_job(results)

    # first value corresponds to function value, the rest to the gradient
    return obj, gradient
end

function JUDI.lsrtm_objective(model::Array{Model,1}, source::Array{judiVector{T,Array{T,2}},1}, dObs::Array{judiVector{T,Array{T,2}},1}, dm::Union{Array{Array{T,1},1}, Array{PhysicalParameter{T},1}}; options=Options(), nlind=false) where T
# lsrtm_objective function for multiple sources and multiple vintages. The function distributes the sources and the input data amongst the available workers.

    println("serial lsrtm_objective on each vintage now")
    ctrl = Array{Any, 1}(undef, length(dObs))
    obj = Array{T, 1}(undef, length(dObs))
    gradient = Array{Array{T,2}, 1}(undef, length(dObs))

    for i = 1:length(dObs)
        _model = model[i]
        _source = source[i]
        _dObs = dObs[i]
        _dm = dm[i]
        _model1 = @bcast _model
        _dm1 = @bcast _dm
        opts = make_opt([_model, dObs, source, _dm])
        ctrl[i] = @batchexec pmap(j -> lsrtm_objective_azure(_model1, _source[j], _dObs[j], _dm1, subsample(options, j); nlind=nlind), 1:dObs[i].nsrc) opts
    end
    for i = 1:length(dObs)
        results = fetch(ctrl[i])
        delete_job(ctrl[i])
        obj[i] = sum([results[j][1] for j =  1:dObs[i].nsrc])
        gradient[i] = sum([results[j][2] for j = 1:dObs[i].nsrc])
    end

    obj1 = sum(obj)

    # first value corresponds to function value, the rest to the gradient
    return obj1, gradient
end

# Make options for unique batch ids
make_opt(a::Array) = AzureClusterlessHPC.Options(;task_name="task_$(objectid(a))")

# Batchexec function
slicec(i::Integer, j::Integer) = i==j ? i : (i:j)
chunk(arr, n) = [arr[slicec(i, min(i + n - 1, length(arr)))] for i in 1:n:length(arr)]
make_parts(iter) = _njpi > 1 ? chunk(iter, _njpi) : iter


text:
function elements2df(els::Vector{Element}; rtod=180.0/pi)
  DataFrame(
    length  = [els[i].length for i in 1:length(els)],
    id      = [els[i].id for i in 1:length(els)],
    od      = [els[i].od for i in 1:length(els)],
    radius  = [els[i].radius for i in 1:length(els)],
    inclination    = [els[i].inclination*rtod for i in 1:length(els)],
    heading = [els[i].heading*rtod for i in 1:length(els)],
    fc      = [els[i].fc for i in 1:length(els)],
    mass    = [els[i].mass for i in 1:length(els)],
    ea      = [els[i].ea for i in 1:length(els)],
    er      = [els[i].er for i in 1:length(els)],
    gj      = [els[i].gj for i in 1:length(els)]
  )
end

function nodes2df(nodes::Vector{Node}; rtod=180.0/pi)
  yvals = [nodes[i].y for i in 1:length(nodes)]
  for i in 1:length(yvals)
    if abs(yvals[i]) < 1.0e-13
      yvals[i] = 0.0
    end
  end
  DataFrame(
    x  = [nodes[i].x for i in 1:length(nodes)],
    y  = yvals,
    z  = [nodes[i].z for i in 1:length(nodes)],
    tvd  = [nodes[i].tvd for i in 1:length(nodes)],
    xvert  = [nodes[i].xvert for i in 1:length(nodes)],
    dispew  = [nodes[i].dispew for i in 1:length(nodes)],
    dispns  = [nodes[i].dispns for i in 1:length(nodes)],
    inclination  = [nodes[i].inclination*rtod for i in 1:length(nodes)],
    dip  = [nodes[i].dip*rtod for i in 1:length(nodes)],
    heading  = [nodes[i].heading*rtod for i in 1:length(nodes)],
    holeradius  = [nodes[i].holeradius for i in 1:length(nodes)],
    stringradius  = [nodes[i].stringradius for i in 1:length(nodes)],
    clearance  = [nodes[i].clearance for i in 1:length(nodes)],
    fc  = [nodes[i].fc for i in 1:length(nodes)],
    weight  = [nodes[i].weight for i in 1:length(nodes)],
    mass  = [nodes[i].mass for i in 1:length(nodes)],
    massinertia  = [nodes[i].massinertia for i in 1:length(nodes)],
    axialinertia  = [nodes[i].axialinertia for i in 1:length(nodes)],
    radialinertia  = [nodes[i].radialinertia for i in 1:length(nodes)],
    deltainclination  = [nodes[i].deltainclination*rtod for i in 1:length(nodes)],
    deltaheading  = [nodes[i].deltaheading*rtod for i in 1:length(nodes)],
    externalxforce  = [nodes[i].externalxforce for i in 1:length(nodes)],
    externalyforce  = [nodes[i].externalyforce for i in 1:length(nodes)],
    fixedxcoord  = [nodes[i].fixedxcoord for i in 1:length(nodes)],
    fixedycoord  = [nodes[i].fixedycoord for i in 1:length(nodes)],
    fixedxtheta  = [nodes[i].fixedxtheta for i in 1:length(nodes)],
    fixedytheta  = [nodes[i].fixedytheta for i in 1:length(nodes)],
  )
end

function segments2df(data::Dict{Symbol, Any})
  syms = [:segment_type, :material, :length, :id, :od, :fc, :angle, :fx, :fy]
  df = DataFrame(Any[DataType[], Symbol[], Float64[], Float64[],
    Float64[], Float64[], Float64[], Float64[], Float64[]], syms)
  for seg in data[:segments]
    if typeof(seg) in [Bit, Stabilizer]
      push!(df, [typeof(seg) seg.material 0.0 seg.id seg.od seg.fc 0.0 0.0 0.0])
    elseif typeof(seg) in [Collar, Pipe]
      push!(df, [typeof(seg) seg.material seg.length seg.id seg.od seg.fc 0.0 0.0 0.0])
    elseif typeof(seg) in [Rig]
      push!(df, [typeof(seg) seg.material 0.0 seg.id seg.od 0.0 0.0 0.0 0.0])
    elseif typeof(seg) in [Bentsub]
      push!(df, [typeof(seg) :none 0.0 0.0 0.0 0.0 seg.angle 0.0 0.0])
    elseif typeof(seg) in [ExtForce]
      push!(df, [typeof(seg) :none 0.0 0.0 0.0 0.0 0.0 seg.fx seg.fy])
    else
      println("Unknown segment_type found.")
    end
  end
  df
end

function sol2df(fname::AbstractString)
  tmp_df = DataFrames.readtable(fname)

	syms = [:depth, :dtb, :radius, :incl, :hdng, :x, :y, :delx, :dely, :fx, :fy, :fz]
	df = DataFrame(Any[Float64[], Float64[], Float64[], Float64[], Float64[], Float64[],
    Float64[], Float64[], Float64[], Float64[], Float64[], Float64[]], syms
	)

  (rows, colums) = size(tmp_df)
  maxdepth = maximum(tmp_df[:depth])

  for i in 1:rows
  	push!(df, [
      tmp_df[:depth][i], tmp_df[:distanceToBit][i],
      tmp_df[:holeRadius][i], tmp_df[:inclination][i], tmp_df[:heading][i],
      round(tmp_df[:x][i], digits=4),
      round(tmp_df[:y][i], digits=4),
      round(tmp_df[:delX][i], digits=4),
      round(tmp_df[:delY][i], digits=5),
      round(tmp_df[:fX][i], digits=5),
      round(tmp_df[:fY][i], digits=5),
      round(tmp_df[:fZ][i], digits=5)
      ]
    )
  end
	df
end

function dyn2df(params::BHAtp.Parameters)
	syms = [:depth, :xampl, :xphase, :yampl, :yphase, :fxamp, :fxphase, :fyampl, :fyphase]
	df = DataFrame(Any[Float64[], Float64[], Float64[], Float64[], Float64[], Float64[],
    Float64[], Float64[], Float64[]], syms
	)
  for i in 1:params.noofnodes
    xampl = round(sqrt(params.del[1,2,i]^2 + params.del[1,3,i]^2), digits=4)
    xphase = xampl > 1.0e-20 ? atan2(params.del[1,2,i], params.del[1,3,i])/params.dtor : 0.0
    yampl = round(sqrt(params.del[2,2,i]^2 + params.del[2,3,i]^2), digits=4)
    yphase = yampl > 1.0e-20 ? atan2(-params.del[2,2,i], -params.del[2,3,i])/params.dtor : 0.0
    fxampl = round(sqrt(params.bf[1,2,i]^2 + params.bf[1,3,i]^2), digits=4)
    fxphase = fxampl > 1.0e-20 ? atan2(params.bf[1,2,i], params.bf[1,3,i])/params.dtor : 0.0
    fyampl = round(sqrt(params.bf[2,2,i]^2 + params.bf[2,3,i]^2), digits=4)
    fyphase = fyampl > 1.0e-20 ? atan2(-params.bf[2,2,i], -params.bf[2,3,i])/params.dtor : 0.0
  	push!(df, [
      params.nodes[i].z/12.0,
      xampl,  xphase,
      yampl,  yphase,
      fxampl,  fxphase,
      fyampl,  fyphase
      ]
    )
  end
	df
end

function tp2df(tptable)
  syms = [:wob, :incl, :tp, :fbit, :dtft, :ft]
	df = DataFrame(Any[Float64[], Float64[], Float64[], Float64[], Float64[], Float64[]],
    syms)
  for i in 1:length(tptable)
    wob = round(tptable[i].wob, digits=4)/1000.0
    inclination = round(round(tptable[i].inclination, digits=4), digits=4)
    dalpha = -round(tptable[i].dalpha, digits=4)
    fbit = round(round(tptable[i].fbit, digits=4), digits=4)
    firsttouch = round(round(tptable[i].firsttouch, digits=4), digits=4)
    f_firsttouch = round(round(tptable[i].f_firsttouch, digits=4), digits=4)
  	push!(df, [wob, inclination, dalpha, fbit,  firsttouch, f_firsttouch])
  end
	df
end

function sol2df(params::BHAtp.Parameters)
	syms = [:depth, :dtb, :radius, :od, :incl, :hdng, :x, :y, :dx, :dy, :dz,
    :fx, :fy, :fz, :friction, :fn, :fs]
	df = DataFrame(Any[Float64[], Float64[], Float64[], Float64[], Float64[],
    Float64[], Float64[], Float64[], Float64[], Float64[], Float64[], Float64[],
    Float64[], Float64[], Float64[], Float64[], Float64[]], syms
	)
  nnodes = params.noofnodes
  maxdepth = maximum(params.nodes[nnodes].z)
  for i in 1:nnodes
  	push!(df, [
      (maxdepth - params.nodes[i].z)/12.0, params.nodes[i].z/12.0,
      params.nodes[i].holeradius, params.nodes[i].stringradius,
      (params.nodes[i].inclination/params.dtor),
      (params.nodes[i].heading/params.dtor),
      params.nodes[i].x, params.nodes[i].y,
      round((params.del[1,1,i] - params.nodes[i].x), digits=4),
      round((params.del[2,1,i] - params.nodes[i].y), digits=4),
      round(params.del[3,1,i], digits=4),
      round(params.f[1, 1, i], digits=5),
      round(params.f[2, 1, i], digits=5),
      round(params.f[3, 1, i], digits=5),
      round(params.friction[i], digits=5),
      round(params.fn[i], digits=5),
      round(params.fs[i], digits=5)
      ]
    )
  end
	df
end


text:
"""
# Polling
"""

"""
Global queue of file descriptors to poll.

`register_for_events()` puts new FDs in `fd_new`.

`poll_task()` calls `poll_wait(::PollQueue)`.

`poll_wait(::PollQueue)` copies new FDs from `fd_new` to `poll_vector`
then passes `poll_vector` (containing C `struct pollfd`s) to
[`poll(2)`](https://man7.org/linux/man-pages/man2/poll.2.html) to wait
for IO activity.

`register_for_events()` uses `wakeup_pipe` to wake `poll(2)` when there is
a new FD waiting in `fd_new`.
"""
struct PollQueue
    fd_new::Channel{FD}
    fd_del::Channel{FD}
    fd_vector::Vector{FD}
    poll_vector::Vector{C.pollfd}
    wakeup_pipe::Vector{Cint}
end

const poll_queue = PollQueue(Channel{FD}(Inf),
                             Channel{FD}(Inf),
                             [],
                             [],
                             [])

@db function poll_queue_init()

    # Pipe for asking `poll(2)` to return before timeout.
    copy!(poll_queue.wakeup_pipe, wakeup_pipe())
    fd = poll_queue.wakeup_pipe[1]
    push!(poll_queue.poll_vector, C.pollfd(fd, C.POLLIN, 0))
    push!(poll_queue.fd_vector, FD{In}(fd))

    # Global Task to run `poll(2)`.
    Threads.@spawn poll_task(poll_queue)              ;@db "@spawn poll_task()"

    if Sys.islinux()
        epoll_queue_init()
    end
end

function wakeup_pipe()
    wakeup_pipe = fill(Cint(-1), 2)
    @cerr C.pipe(wakeup_pipe)
    for fd in wakeup_pipe
        fcntl_setfl(fd, C.O_NONBLOCK)
    end
    return wakeup_pipe
end

"""
Write a byte to the pipe to wake up `poll(2)` before its timeout.
"""
function wakeup_poll(q::PollQueue);                       @db 3 "wakeup_poll()"
    @cerr C.write(q.wakeup_pipe[2], Ref(0), 1)
    nothing
end

"""
Wait for an event to occur on `fd`.
"""
@db 2 function wait_for_event(queue, fd::FD)
    assert_havelock(fd)

    timer = register_timer(fd.deadline) do
        @dblock fd notify(fd, :timeout)
    end
    fd.nwaiting += 1
    event = try
        register_for_events(queue, fd)
        @db 2 "$(fd.nwaiting) waiting for $fd..."
        wait(fd.ready) # Wait for: `poll_task()`
    finally
        close(timer)
        @assert fd.nwaiting > 0
        fd.nwaiting -= 1
        if fd.nwaiting == 0
            unregister_for_events(queue, fd)
        end
    end
    @db 2 return event "Ready: $fd"
end

@db 4 function register_for_events(q::PollQueue, fd)
    put!(q.fd_new, fd)
    wakeup_poll(q)                                           ;@db 5 poll_queue
end

@db 4 function unregister_for_events(q::PollQueue, fd)
    put!(q.fd_del, fd)
    wakeup_poll(q)
end

"""
Run `poll_wait()` in a loop.
"""
@db 1 function poll_task(q)

    if Threads.threadid() == 1
        @warn "UnixIO.poll_task() is running on thread No. 1!\n" *
              "Other Tasks on thread 1 may be blocked for up to 100ms while " *
              "poll_task() is waiting for IO.\n" *
              "Consider increasing JULIA_NUM_THREADS (`julia --threads N`)."

        timeout_ms = 100
    else
        timeout_ms = 60_000
    end

    while true
        try
            poll_wait(q, timeout_ms) do events, fd
                if events & (C.POLLHUP | C.POLLNVAL) != 0
                    fd.gothup = true
                    @db 1 "$(db_c(events,"POLL")) -> $fd 💥"
                end
                if fd.nwaiting <= 0
                    @db 1 "$(db_c(events,r"POLL[A-Z]")) -> $fd None Waiting!"
                else
                    @db 2 "$(db_c(events,r"POLL[A-Z]")) -> notify($fd)"
                    @dblock fd notify(fd, events); # Wake: `wait_for_event()`
                end
            end
        catch err
            notify_error(q, err)
            exception=(err, catch_backtrace())
            @error "Error in poll_task()" exception
        end
        yield()
    end
    @assert false
end

function notify_error(q::PollQueue, err)
    for fd in q.fd_vector
        @lock fd.ready Base.notify_error(fd.ready, err)
    end
end

gc_safe_poll(fds, nfds, timeout_ms) = @gc_safe C.poll(fds, nfds, timeout_ms)

"""
- Update `poll_vector` based on `fd_del` and `fd_new`.
- Wait for events.
- Call `f(events, fd)` for each event.
"""
@db 4 function poll_wait(f::Function, q::PollQueue, timeout_ms::Int)

    fdv = q.fd_vector
    pollv = q.poll_vector
    timeout_ms = next_timer_deadline_ms(timeout_ms)
    @db 4 " [ poll($timeout_ms ms)..."

    # Remove old entries from fd_vector and poll_vector.
    while !isempty(q.fd_del)
        fd = take!(q.fd_del)
        i = findfirst(isequal(fd), fdv)
        if i != nothing                                        ;@db 4 "del $fd"
            deleteat!(pollv, i)
            deleteat!(fdv, i)
        end
    end

    # Add new entries to fd_vector and poll_vector.
    while !isempty(q.fd_new)
        fd = take!(q.fd_new)                                   ;@db 4 "add $fd"
        pushfirst!(fdv, fd)
        pushfirst!(pollv, C.pollfd(fd, poll_event_type(fd), 0))
    end

    # Wait for events
    n = @cerr(allow=C.EINTR,
              gc_safe_poll(pollv, length(pollv), timeout_ms))          ;@db 5 n
                                                                   ;@db 4 pollv
    # Check poll vector for events.
    indexes_to_delete = Int[]
    for (i, e) in enumerate(pollv)
        if e.revents != 0
            if e.fd == q.wakeup_pipe[1]
                C.read(q.wakeup_pipe[1], Ref(0), 1)        ;@db 4 "got wakeup!"
            else
                push!(indexes_to_delete, i)
                @assert RawFD(e.fd) == fdv[i].fd
                f(e.revents, fdv[i])
            end
        end
    end
    deleteat!(pollv, indexes_to_delete)                        ;@db 6 del pollv
    deleteat!(fdv, indexes_to_delete)

    # Check for timeouts.
    notify_timers()
    nothing
end

# Linux epoll(7)             https://man7.org/linux/man-pages/man7/epoll.7.html

"""
See `struct epoll_event` in "sys/epoll.h"
"""
struct epoll_event
    events::UInt32
    data::UInt64
end

"""
See PollQueue above.
"""
mutable struct EPollQueue
    fd::RawFD
    dict::Dict{RawFD,FD}
    cvector::Vector{epoll_event}
    lock::Threads.SpinLock
end

const epoll_queue = EPollQueue(RawFD(-1),
                               Dict{RawFD,FD}(),
                               Vector{epoll_event}(undef, 10),
                               Threads.SpinLock())

@db function epoll_queue_init()

    @assert Sys.ARCH != :x86_64 """
        FIXME: UnixIO does not support epoll on x86_64 because
        `struct epoll_event` has `__attribute__((packed))` on x86.
        See https://git.io/JCGMK and https://git.io/JCGDz
        This is not difficult to work around but is not handled yet.
    """
    @assert C.EPOLLIN == C.POLLIN
    @assert C.EPOLLOUT == C.POLLOUT
    @assert C.EPOLLHUP == C.POLLHUP

    # Global FD interface to Linux epoll(7).
    epoll_queue.fd = RawFD(@cerr C.epoll_create1(C.EPOLL_CLOEXEC))

    # Global Task to run `epoll_wait(7)`.
    Threads.@spawn poll_task(epoll_queue)
end

"""
Add, modify or delete epoll target FDs.
See [epoll_ctl(7)(https://man7.org/linux/man-pages/man7/epoll_ctl.7.html)
"""
@db 4 function epoll_ctl(fd, op, events=0, data=fd; allow=())
                                                  @db 4 fd db_c(op,"EPOLL_CTL")
    e = [epoll_event(events, fd)]
    GC.@preserve e @cerr(allow=allow,
                         C.epoll_ctl(epoll_queue.fd, op, fd, pointer(e)))
end

epoll_ctl(fd::FD, op, events=0; kw...) =
    epoll_ctl(convert(Cint, fd), op, events; kw...)

"""
Register `fd` to wake up `epoll_wait(7)` on `event`:
"""
@db 4 function register_for_events(q::EPollQueue, fd::FD)
    @dblock q.lock begin
        if !haskey(q.dict, fd.fd)
            q.dict[fd.fd] = fd
            epoll_ctl(fd, C.EPOLL_CTL_ADD, poll_event_type(fd))
        else
            @db 2 "Already registered!"
        end
    end
end

@db 4 function unregister_for_events(q::EPollQueue, fd::FD)
    @dblock q.lock begin
        if haskey(q.dict, fd.fd)
            delete!(q.dict, fd.fd)
            epoll_ctl(fd, C.EPOLL_CTL_DEL)#; allow=C.ENOENT)
        end
    end
end

poll_event_type(::FD{In}) = C.POLLIN
poll_event_type(::FD{Out}) = C.POLLOUT

"""
Call `epoll_wait(7)` to wait for events.
Call `f(events, fd)` for each event.
"""
@db 4 function poll_wait(f::Function, q::EPollQueue, timeout_ms::Int)
    v = q.cvector                                                 ;@db 6 q.dict
    n = @cerr(allow=C.EINTR,
              gc_safe_epoll_wait(q.fd, v, length(v), timeout_ms))    ;@db 6 n v
    if n >= 0
        for i in 1:n
            fd = get(q.dict, RawFD(v[i].data), nothing)             ;@db 6 i fd
            if fd == nothing
                @error "Ignoring epoll_wait() notification v[$i] = $(v[i]) " *
                       "for unknown FD (already deleted?)." v[i] q.dict
            else
                unregister_for_events(q, fd)
                f(v[i].events, fd)
            end
        end
    end
    nothing
end

function notify_error(q::EPollQueue, err)
    for fd in values(q.dict)
        @lock fd Base.notify_error(fd.ready, err)
    end
end

gc_safe_epoll_wait(epfd, events, maxevents, timeout_ms) =
    @gc_safe C.epoll_wait(epfd, events, maxevents, timeout_ms)

# Pretty Printing.

Base.show(io::IO, fd::C.pollfd) =
    print(io, "(", fd.fd, ", ", fd.events, ", ", fd.revents, ")")

Base.show(io::IO, q::PollQueue) =
    dbprint(io, "(", q.fd_vector, ", ", q.poll_vector, ")")

# End of file: poll.h


text:
module SelfOrganizingMaps

# package code goes here

end # module

max_stars_repo_name: nathandesdouits/SelfOrganizingMaps

max_issues_repo_name: nathandesdouits/SelfOrganizingMaps

max_forks_repo_name: nathandesdouits/SelfOrganizingMaps


text:
using HybridArrays, ArrayInterface, Test, StaticArrays

@testset "ArrayInterface compatibility" begin
    M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2; 4 5])
    @test ArrayInterface.ismutable(M)
    @test ArrayInterface.can_setindex(M)
    @test ArrayInterface.parent_type(M) === Matrix{Int}
    @test ArrayInterface.restructure(M, [2, 4, 6, 8]) == HybridMatrix{2, StaticArrays.Dynamic()}([2 6; 4 8])
    @test isa(ArrayInterface.restructure(M, [2, 4, 6, 8]), HybridMatrix{2, StaticArrays.Dynamic()})
end


text:
################################
# split-Cholesky Factorization #
################################

#
# A split-Cholesky factorization of an SPD banded matrix A is:
#
#   S = [ U    ]
#       [ M  L ],
#
# where U and M are upper-triangular and L is lower-triangular so that SᵀS = A.
#
# if n = size(A, 1), n == size(A, 2) and kd = bandwidth(A), then if m = (n+b)÷2,
#
# size(U) == (m, m),  size(L) == (n-m, n-m),  and size(M) == (n-m, m).
#

struct SplitCholesky{T,M<:AbstractMatrix{T}} <: Factorization{T}
    factors::M
    uplo::Char

    function SplitCholesky{T,M}(factors, uplo) where {T,M<:AbstractMatrix}
        require_one_based_indexing(factors)
        new{T,M}(factors, uplo)
    end
end

SplitCholesky(factors::AbstractMatrix{T}, uplo::Char) where T = SplitCholesky{T,typeof(factors)}(factors, uplo)

size(S::SplitCholesky) = size(S.factors)
size(S::SplitCholesky, i::Integer) = size(S.factors, i)

splitcholesky!(A::Symmetric{T,<:BandedMatrix{T}}) where T = splitcholesky!(A, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
splitcholesky!(A::Symmetric{T,<:BandedMatrix{T}}, ::Type{Tr}) where {T,Tr} = splitcholesky!(MemoryLayout(A), A, Tr)

function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor},
                       A::AbstractMatrix{T}, ::Type{UpperTriangular}) where T<:BlasFloat
    _, info = pbstf!('U', size(A, 1), bandwidth(A,2), symbandeddata(A))
    SplitCholesky(A, A.uplo)
end

function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor},
                       A::AbstractMatrix{T}, ::Type{LowerTriangular}) where T<:BlasFloat
    _, info = pbstf!('L', size(A, 1), bandwidth(A,1), symbandeddata(A))
    SplitCholesky(A, A.uplo)
end

adjoint(S::SplitCholesky) = Adjoint(S)

function lmul!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
    require_one_based_indexing(B)
    n, nrhs = size(B, 1), size(B, 2)
    if size(S, 1) != n
        throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))
    end
    A = S.factors
    b = bandwidth(A, 1)
    m = (n+b)÷2
    @inbounds for l = 1:nrhs
        for j = n:-1:m+1
            t = zero(T)
            @simd for k = j-b:j
                t += A[j,k]*B[k,l]
            end
            B[j,l] = t
        end
        for j = 1:m
            t = zero(T)
            @simd for k = j:min(j+b,m)
                t += A[j,k]*B[k,l]
            end
            B[j,l] = t
        end
    end
    B
end

function ldiv!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
    require_one_based_indexing(B)
    n, nrhs = size(B, 1), size(B, 2)
    if size(S, 1) != n
        throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))
    end
    A = S.factors
    b = bandwidth(A, 1)
    m = (n+b)÷2
    @inbounds for l = 1:nrhs
        for j = m:-1:1
            t = zero(T)
            @simd for k = j+1:min(j+b,m)
                t += A[j,k]*B[k,l]
            end
            B[j,l] = (B[j,l]-t)/A[j,j]
        end
        for j = m+1:n
            t = zero(T)
            @simd for k = j-b:j-1
                t += A[j,k]*B[k,l]
            end
            B[j,l] = (B[j,l]-t)/A[j,j]
        end
    end
    B
end

function lmul!(S::Adjoint{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
    require_one_based_indexing(B)
    n, nrhs = size(B, 1), size(B, 2)
    if size(S, 1) != n
        throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))
    end
    A = S.parent.factors
    b = bandwidth(A, 1)
    m = (n+b)÷2
    @inbounds for l = 1:nrhs
        for j = m:-1:1
            t = zero(T)
            @simd for k = max(1,j-b):j
                t += A[k,j]*B[k,l]
            end
            B[j,l] = t
        end
        for j = m-b+1:m
            t = zero(T)
            @simd for k = m+1:j+b
                t += A[k,j]*B[k,l]
            end
            B[j,l] += t
        end
        for j = m+1:n
            t = zero(T)
            @simd for k = j:min(j+b,n)
                t += A[k,j]*B[k,l]
            end
            B[j,l] = t
        end
    end
    B
end

function ldiv!(S::Adjoint{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
    require_one_based_indexing(B)
    n, nrhs = size(B, 1), size(B, 2)
    if size(S, 1) != n
        throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))
    end
    A = S.parent.factors
    b = bandwidth(A, 1)
    m = (n+b)÷2
    @inbounds for l = 1:nrhs
        for j = n:-1:m+1
            t = zero(T)
            @simd for k = j+1:min(j+b,n)
                t += A[k,j]*B[k,l]
            end
            B[j,l] = (B[j,l]-t)/A[j,j]
        end
        for j = m:-1:m-b+1
            t = zero(T)
            @simd for k = m+1:j+b
                t += A[k,j]*B[k,l]
            end
            B[j,l] -= t
        end
        for j = 1:m
            t = zero(T)
            @simd for k = max(1,j-b):j-1
                t += A[k,j]*B[k,l]
            end
            B[j,l] = (B[j,l]-t)/A[j,j]
        end
    end
    B
end


text:
@testset "cfr" begin

function test_nash_kuhn_poker(game, policy)
    game_values = expected_returns(new_initial_state(game), policy, -1)

    # 1/18 is the Nash value. See https://en.wikipedia.org/wiki/Kuhn_poker
    nash_value = 1.0 / 18.0
    eps = 1e-3
    @test length(game_values) == 2
    @test isapprox(game_values[1], -nash_value, atol=eps)
    @test isapprox(game_values[2], nash_value, atol=eps)
end

test_exploitability_kuhn_poker(game, policy) = @test exploitability(game, policy) <= 0.05

@testset "CFRSolver" begin
    game = load_game("kuhn_poker")
    solver = CFRSolver(game)
    for _ in 1:300
        evaluate_and_update_policy(solver)
    end
    avg_policy = average_policy(solver)
    test_nash_kuhn_poker(game, avg_policy)
    test_exploitability_kuhn_poker(game, avg_policy)
end

@testset "CFRPlusSolver" begin
    game = load_game("kuhn_poker")
    solver = CFRPlusSolver(game)
    for _ in 1:200
        evaluate_and_update_policy(solver)
    end
    avg_policy = average_policy(solver)
    test_nash_kuhn_poker(game, avg_policy)
    test_exploitability_kuhn_poker(game, avg_policy)
end

end


text:
using HorizonSideRobots

function make_cross(r::Robot)
    for side in (Sud,West,Ost,Nord)
        putmarkers!(r,side)
        move_markers(r,inverse(side))
    end
end
function putmarkers!(r::Robot,side::HorizonSide)
    while  isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

function move_markers(r::Robot,side::HorizonSide)
    while ismarker(r)==true
        move!(r,side)
    end
end

function inverse(side::HorizonSide)
    new_side=HorizonSide(mod(Int(side)+2,4))
    return new_side
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


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

name = "ZPares"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/cometscome/zpares_mirror.git", "1e84be7cd0f8368da09ac0c262f815583ce04b7d")
    FileSource("https://zpares.cs.tsukuba.ac.jp/?download=242","3c34257d249451b0b984abc985e296ebb73ae5331025f1b8ea08d50301c7cf9a",filename="zpares_0.9.6a.tar.gz")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
tar -xvf zpares_0.9.6a.tar.gz
cd zpares_mirror/
install_license ./LICENSE
cd ../zpares_0.9.6a
cp ../zpares_mirror/wrapper/zpares_wrapper.f90 ./

if [[ $target == *"aarch64-apple-darwin"* ]]; then
    Rankmismatch="-fallow-argument-mismatch"
fi

cp Makefile.inc/make.inc.gfortran.seq ./make.inc
make BLAS="-L${libdir} -lopenblas" USE_MPI="0" FFLAG="-O3 ${Rankmismatch} -shared -fPIC -L${libdir} -lopenblas" LAPACK="-L./"
gfortran -O3 -shared -fPIC zpares_wrapper.f90 -I./include -L./lib -lzpares -L${libdir} -lopenblas -o "${libdir}/libzpares.${dlext}"
cp include/zpares.mod "${includedir}"
cp zpares_wrapper.mod "${includedir}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)
platforms = expand_gfortran_versions(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))

]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using  DataStructures

t = range(0.0, stop=10.0, length=100)

result = OrderedDict{String,Any}()

result["time"] = t
result["phi"]  = sin.(t)

plot(result, "phi", heading="Sine(time)")

println("\n... result info:")
printResultInfo(result)


text:
function generate_kernel_functions(dist::Distributions.Distribution)
    kernel(x) = pdf.(dist, x)
    kernel_int(x) = cdf.(dist, x)
    return kernel, kernel_int
end


text:
"""
    nipals(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing, AbstractString,AbstractArray}="", dim::Number=3, numepoch::Number=3, lower::Number=0, upper::Number=1.0f+38, offsetStoch::Number=1f-20, logdir::Union{Nothing,AbstractString}=nothing)
Online CCA solved by NIPALS.
Input Arguments
---------
- `input` : Julia Binary file generated by `OnlineCCA.csv2bin` function. Only two files can be specified (e.g. X.zst and Y.zst)
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `colmeanlist` : The mean of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `colvarlist` : The variance of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `dim` : The number of dimension of CCA.
- `numepoch` : The number of epoch.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `evalfreq` : Evaluation Frequency of Reconstruction Error
- `offsetStoch` : Off set value for avoding overflow when calculating stochastic gradient
---------
- `T` : = (t1,t2,...,tK) (No. columns of X × dim), the score vectors of X
- `P` : = (p1,p2,...,pK) (No. rows of X × dim), the loading of X
- `U` : = (u1,u2,...,uK) (No. columns of Y × dim), the score vectors of Y
- `Q` : = (q1,q2,...,qK) (No. rows of Y × dim), the loading of Y
- `β` : The Regression of Coefficient (dim × dim), Y is predicted by the new X by Y_new = X_new P β Q'
"""
function nipals(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing,AbstractString,AbstractArray}="", dim::Number=3, numepoch::Number=3, lower::Number=0, upper::Number=1.0f+38, offsetStoch::Number=1f-20, logdir::Union{Nothing,AbstractString}=nothing)
    # Initial Setting
    cca = NIPALS()
    pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch = init(input, pseudocount, dim, colmeanlist, colvarlist, cca, lower, upper, offsetStoch, scale)
    # Perform CCA
    out = nipals(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, numepoch, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch)
    if outdir isa String
        output(outdir, out, cca)
    end
    return out
end

function nipals(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, numepoch, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch)
    @assert length(input) == 2
    l = 2
    N, M = nm(input)
    for i=1:l
        @assert 0 < dim ≤ min(N[i], M[i]) - 1
    end
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    X = input[1]
    Y = input[2]
    x = zeros(UInt32, M[1])
    normx = zeros(UInt32, M[1])
    defx = zeros(UInt32, M[1])
    y = zeros(UInt32, M[2])
    normy = zeros(UInt32, M[2])
    defy = zeros(UInt32, M[2])
    T = zeros(Float32, N[1], dim)
    P = zeros(Float32, M[1], dim)
    U = zeros(Float32, N[2], dim)
    Q = zeros(Float32, M[2], dim)
    tmp = tempdir()
    old_t = zeros(Float32, N[1])
    # NIPALS Iteration
    for k=1:dim
        println("#######################################")
        println(string(k)*"-th u, t, p, and q are calculated...")
        println("#######################################")
        outsvd = oocmcca(input, outdir, scale, pseudocount, colmeanlist, colvarlist, 1, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, offsetStoch, false)
        u = outsvd[3][1][:,1]
        t = outsvd[3][2][:,1]
        p = zeros(Float32, M[1])
        q = zeros(Float32, M[2])
        iter = 1
        relchange = Any[]
        push!(relchange, 10^20)
        while (relchange[iter] > lower) && (iter <= numepoch)
            println(string(iter)*"-th iteration : ||t_old - t_new|| = "*string(relchange[iter]))
            old_t .= t
            ####################
            ##### p = X'u ######
            ####################
            open(input[1]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[1]
                    # Row vector of data matrix
                    read!(stream, x)
                    if k == 1
                        normx = normalizex(x, n, stream, scale[1], pseudocount[1], colmeanlist[1], colmeanvec[1], colvarlist[1], colvarvec[1])
                    else
                        normx = x
                    end
                    # p = X'u
                    p .+= normx[:,1] .* u[n]
                end
                close(stream)
            end
            p = p ./ norm(p)

            ####################
            ###### t = Xp ######
            ####################
            open(input[1]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[1]
                    # Row vector of data matrix
                    read!(stream, x)
                    if k == 1
                        normx = normalizex(x, n, stream, scale[1], pseudocount[1], colmeanlist[1], colmeanvec[1], colvarlist[1], colvarvec[1])
                    else
                        normx = x
                    end
                    # t = Xp
                    t[n] = normx[:,1]'p
                end
                close(stream)
            end

            ####################
            ###### q = Y't #####
            ####################
            open(input[2]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[2]
                    # Row vector of data matrix
                    read!(stream, y)
                    if k == 1
                        normy = normalizex(y, n, stream, scale[2], pseudocount[2], colmeanlist[2], colmeanvec[2], colvarlist[2], colvarvec[2])
                    else
                        normy = y
                    end
                    # q = Y't
                    q .+= normy[:,1] .* t[n]
                end
                close(stream)
            end
            q = q ./ norm(q)

            ####################
            ###### u = Yq ######
            ####################
            open(input[2]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[2]
                    # Row vector of data matrix
                    read!(stream, y)
                    if k == 1
                        normy = normalizex(y, n, stream, scale[2], pseudocount[2], colmeanlist[2], colmeanvec[2], colvarlist[2], colvarvec[2])
                    else
                        normy = y
                    end
                    # u = Yq
                    u[n] = normy[:,1]'q
                end
                close(stream)
            end
            # Relative Change
            iter = iter + 1
            push!(relchange, abs(norm(old_t - t)))
        end
        # Update
        T[:,k] = t
        P[:,k] = p
        U[:,k] = u
        Q[:,k] = q
        ####################
        #### Deflation #####
        ####################
        if k < dim
            defXfile = tmp*"/defX"*string(k)*".zst"
            defYfile = tmp*"/defY"*string(k)*".zst"
            open(defXfile, "w") do deffile
                defstream = ZstdCompressorStream(deffile)
                open(input[1]) do file
                    stream = ZstdDecompressorStream(file)
                    read!(stream, tmpN)
                    read!(stream, tmpM)
                    write(defstream, tmpN)
                    write(defstream, tmpM)
                    # Each step n
                    for n=1:N[1]
                        # Row vector of data matrix
                        read!(stream, x)
                        normx = normalizex(x, n, stream, scale[1], pseudocount[1], colmeanlist[1], colmeanvec[1], colvarlist[1], colvarvec[1])
                        # X = X - tp'
                        write(defstream, normx[:,1]' - t[n]*p')
                    end
                    close(stream)
                end
                close(defstream)
            end

            open(defYfile, "w") do deffile
                defstream = ZstdCompressorStream(deffile)
                open(input[2]) do file
                    stream = ZstdDecompressorStream(file)
                    read!(stream, tmpN)
                    read!(stream, tmpM)
                    write(defstream, tmpN)
                    write(defstream, tmpM)
                    # Each step n
                    for n=1:N[2]
                        # Row vector of data matrix
                        read!(stream, y)
                        normy = normalizex(y, n, stream, scale[2], pseudocount[2], colmeanlist[2], colmeanvec[2], colvarlist[2], colvarvec[2])
                        # Y = Y - uq'
                        write(defstream, normy[:,1]' - u[n]*q')
                    end
                    close(stream)
                end
                close(defstream)
            end
            input[1] = defXfile
            input[2] = defYfile
        end
    end
    # Regression Coefficient
    W,Σ,V = svd(T'T)
    β = V*Diagonal(1/Σ)*W'T'U
    return (T, P, U, Q, β)
end