text:
# Remove legend, axes, grid, and force the aspect ratio for shape plots
@recipe function f(::Type{Val{:shapeplot}}, plt::Plot)
    legend --> false
    ticks --> nothing
    grid --> false
    aspect_ratio --> 1
    seriestype := :shape
    ()
end
@shorthands shapeplot

# ----------------------------------------------------------------
# Shapefile.jl plotting

if Plots.is_installed("Shapefile")
    @eval begin
        import Shapefile
        function shapefile_coords(poly::Shapefile.ESRIShape)
            start_indices = poly.parts+1
            end_indices = vcat(poly.parts[2:end], length(poly.points))
            x, y = zeros(0), zeros(0)
            for (si,ei) in zip(start_indices, end_indices)
                push!(x, NaN)
                push!(y, NaN)
                for pt in poly.points[si:ei]
                    push!(x, pt.x)
                    push!(y, pt.y)
                end
            end
            x, y
        end

        function shapefile_coords{T<:Shapefile.ESRIShape}(polys::AbstractArray{T})
            x, y = [], []
            for poly in polys
                xpart, ypart = shapefile_coords(poly)
                push!(x, xpart)
                push!(y, ypart)
            end
            x, y
        end

        @recipe f(poly::Shapefile.ESRIShape) = (seriestype --> :shape; shapefile_coords(poly))
        @recipe f{T<:Shapefile.ESRIShape}(polys::AbstractArray{T}) = (seriestype --> :shape; shapefile_coords(polys))
    end
end

# ----------------------------------------------------------------


text:
using StaticArrays

const PASS = 0
const BET = 1
const KuhnInfoKey = Tuple{Int, Int, SVector{3,Int}} # [player, player_card, action_hist]

struct Hist
    cards::SVector{2,Int}
    action_hist::SVector{3,Int}
end

function Base.length(h::Hist)
    l = 0
    for a in h.action_hist
        a === -1 && break
        l += 1
    end
    return l
end

struct Kuhn <: Game{Hist, KuhnInfoKey}
    cards::Vector{SVector{2,Int}}
    Kuhn() = new([SVector(i,j) for i in 1:3, j in 1:3 if i != j])
end

CounterfactualRegret.initialhist(::Kuhn) = Hist(SA[0,0], SA[-1,-1,-1])

function CounterfactualRegret.isterminal(::Kuhn, h::Hist) # requires some sequence of actions
    L = length(h)
    h = h.action_hist
    if L > 1
        return h[1] == BET || h[2] == PASS || L > 2
    else
        return false
    end
end

function CounterfactualRegret.utility(::Kuhn, i::Int, h::Hist)
    as = h.action_hist
    cards = h.cards
    L = length(as)
    has_higher_card = cards[i] > cards[other_player(i)]
    if as == SA[PASS, PASS, -1]
        return has_higher_card ? 1. : -1.
    elseif as == SA[PASS, BET, PASS]
        return i==2 ? 1. : -1.
    elseif as == SA[PASS, BET, BET]
        return has_higher_card ? 2. : -2.
    elseif as == SA[BET, PASS, -1]
        return i==1 ? 1. : -1.
    elseif as == SA[BET, BET, -1]
        return has_higher_card ? 2. : -2.
    else
        return 0.
    end
end

function CounterfactualRegret.player(::Kuhn, h::Hist)
    if any(iszero, h.cards)
        return 0
    else
        return length(h)%2 + 1
    end
end

CounterfactualRegret.player(::Kuhn, k::KuhnInfoKey) = first(k)

function CounterfactualRegret.chance_actions(game::Kuhn, h::Hist)
    return game.cards
end

function CounterfactualRegret.chance_action(game::Kuhn, h::Hist)
    return rand(game.cards)
end

function CounterfactualRegret.next_hist(::Kuhn, h, a::SVector{2,Int})
    return Hist(a, h.action_hist)
end

function CounterfactualRegret.next_hist(::Kuhn, h::Hist, a::Int)
    L = length(h)
    action_hist = setindex(h.action_hist, a, L+1)
    return Hist(h.cards, action_hist)
end

function CounterfactualRegret.infokey(g::Kuhn, h::Hist)
    p = player(g,h)
    card = p > 0 ? h.cards[p] : 0
    return (p, card, h.action_hist) # [player, player_card, action_hist]
end

CounterfactualRegret.actions(::Kuhn, h::Hist) = PASS:BET

## Extra

function Base.print(io::IO, solver::AbstractCFRSolver{K,G}) where {K,G<:Kuhn}
    println(io)
    for (k,v) in solver.I
        h = k[3]
        Lp1 = findfirst(==(-1), h)
        h_str = if isnothing(Lp1)
            "_"*3
        else
            rpad(join(h[1:Lp1-1]),3,"_")
        end
        σ = copy(v.s)
        σ ./= sum(σ)
        σ = round.(σ, digits=3)
        println(io, "Player: $(k[1]) \t Card: $(k[2]) \t h: $(join(h_str)) \t σ: $σ")
    end
end

ext: jl

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
"""
    FixedPoint.jl
exports `afps` function, for further help type `?afps` in the REPL.
"""
module FixedPoint
"""
    afps(f, x; iters::Int = 5000, vel::Float64 = 0.9, ep::Float64 = 0.01, tol::Float64 = 1e-12, grad_norm=x->maximum(abs,x))

solve equation `f(x) = x` according to:

    `f` : function to find fixed point for

    `x` : initial condition, ideally it should be close to the final solution

    `vel` : amount of Nesterov acceleration in [0,1]

    `ep` : learning rate, typically in ]0,1[

    `tol` : absolute tolerance on |f(x)-x|

    `grad_norm` : function to evaluate the norm for |f(x)-x|

returns a named tuple (x, error, iters) where:

    `x` : is the solution found for f(x)=x

    `error` : is the norm of f(x)-x at the solution point

    `iters` : total number of iterations performed
"""
function afps(
    f::Fun,
    x::Mat;
    iters::Int = 5000,
    vel::T = 0.9,
    ep::T = 0.01,
    tol::T = 1e-12,
    grad_norm = x -> maximum(abs, x),
) where {T <: Number, Mat<:Union{AbstractArray{T},T}, Fun <: Function}
    x_n = identity.(x)
    v_n = zero(x_n)
    β = vel
    ϵ = ep
    runs = 0
    for _ = 1:iters
        trial = x_n + β * v_n
        g = f(trial) - trial
        v_n = β * v_n + ϵ * g
        x_n = x_n + v_n
        runs += 1
        if grad_norm(g) < tol
            break
        end
    end
    (x = x_n, error = grad_norm(f(x_n) - x_n), iters = runs)
end

"""
    afps!(f!, x; iters::Int = 5000, vel::Float64 = 0.9, ep::Float64 = 0.01, tol::Float64 = 1e-12, grad_norm=x->maximum(abs,x))

solve equation `f(x) = x` according to:

    `f!` : inplace version of function to find fixed point for, calling `f!(out,x)` should amount to writing `out = f(x)`

    `x` : initial condition, ideally it should be close to the final solution

    `vel` : amount of Nesterov acceleration in [0,1]

    `ep` : learning rate, typically in ]0,1[

    `tol` : absolute tolerance on |f(x)-x|

    `grad_norm` : function to evaluate the norm for |f(x)-x|

returns a named tuple (x, error, iters) where:

    `x` : is the solution found for f(x)=x

    `error` : is the norm of f(x)-x at the solution point

    `iters` : total number of iterations performed
"""
function afps!(
    f!::Fun,
    x_n::Mat;
    iters::Int = 5000,
    vel::T = 0.9,
    ep::T = 0.01,
    tol::T = 1e-12,
    grad_norm = x -> maximum(abs, x),
) where {T <: Number, Mat<:Union{AbstractArray{T},T}, Fun <: Function}
    v_n = zero(x_n)
    trial = zero(x_n)
    g = zero(x_n)
    β = vel
    ϵ = ep
    runs = 0
    for _ = 1:iters
        #prediction using velocity
        x_n .+= β .* v_n
        #eval gradient
        f!(g,x_n)
        g .-= x_n
        #correction
        x_n .+= ϵ .* g
        #update velocity
        v_n .= β .* v_n .+ ϵ .* g
        runs += 1
        if grad_norm(g) < tol
            break
        end
    end
    (x = x_n, error = grad_norm(g), iters = runs)
end

export afps,afps!
end


text:
module EstimationStatistics

using Plots
using KernelDensity
using Measures
using StatsBase
using Distributions
using Statistics
using DataFrames
using DataStructures

include("quasirandomScatter.jl")
include("jackknife.jl")
include("bootstrap.jl")
include("estimationPlot.jl")
include("addSummaryStat.jl")

struct ConfidenceInterval
    low::Float64
    high::Float64
end

ConfidenceInterval(arr) = ConfidenceInterval(arr[1], arr[2])

struct Bootstrap
    boot::Array
    ci::ConfidenceInterval
end

function customSort!(df::DataFrame, sortops)
    sortv = []
    sortOptions = []
    if(isa(sortops, Array))
        sortv = sortops
    else
        push!(sortv,sortops)
    end
    for i in sortv
        if(isa(i, Tuple))
            if (isa(i[2], Array)) # The second option is a custom order
                orderArray = Array(collect(union(OrderedSet(i[2]),  OrderedSet(unique(df[!, i[1]])))))
                push!(sortOptions, order(i[1], by = x->Dict(x => i for (i,x) in enumerate(orderArray))[x] ))
            else                  # The second option is a reverse direction flag
                push!(sortOptions, order(i[1], rev = i[2]))
            end
        else
          push!(sortOptions, order(i))
        end
    end
    return sort!(df, sortOptions)
end

export Bootstrap, ConfidenceInterval
export quasirandomScatter, addSummaryStat!, estimationPlot
export meanDiff, medianDiff, BCaBoot

end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Filter records in a DataFrame with ability to clear number of rows before and after some occurrence
"""
module FilterRecords

export filter_by_bool!, all_occurrences

using DataFrameSets
using DataFramesMeta

"""
Same as filter_by_bool! over a single dataframe D, but applies it over all records in a DFSet
"""
function filter_by_bool!(f::Function, Ds::DFSet; kwargs...)
  for D in getrecords(Ds)
    filter_by_bool!(f, D; kwargs...)
  end
end

"""
Filter rows of dataframe D according to callback f.  Rows with f returning true are filtered,
along with n_before rows before occurrence, and n_after rows after occurrence.
Safe to use typemax for n_before and n_after.
f is a callback of the form bool=f(row::DataFrameRow).
"""
function filter_by_bool!(f::Function, D::DataFrame;
                        n_before::Int64=0, #safe to use typemax for all previous
                        n_after::Int64=0) #safe to use typemax for all after
  @assert n_before >= 0 && n_after >= 0

  remrows = IntSet()
  rows = all_occurrences(f, D)
  for r in rows
    istart = max(r - n_before, 1)
    n_after = min(n_after, nrow(D)) #workaround for typemax(Int64) case
    iend = min(r + n_after, nrow(D))
    push!(remrows, istart:iend...)
  end

  if !isempty(remrows)
    deleterows!(D, [remrows...])
  end
end

"""
Returns indices of all rows in dataframe 'D' where f returns true.
f is a callback of the form bool=f(row::DataFrameRow).
"""
function all_occurrences(f::Function, D::DataFrame)
  rows = IntSet()
  for (i, r) in enumerate(eachrow(D))
    f(r) && push!(rows, i)
  end
  return rows
end

end #module


text:
module SpectrumSlicingTestSuite

export convergenceRatio,testIntervals

import Polynomials
import GaussQuadrature
using Printf
using LinearAlgebra

include("./chebyshevpoly.jl")
include("./kpm.jl")
include("./phi.jl")
include("./kpmdos.jl")
include("./endPoints.jl")
include("./testIntervals.jl")
include("./convergenceratio.jl")

end # module


text:
"""
A generic, framework agnostic interface for neural networks.
"""
module Network

export AbstractNetwork, OptimiserSpec, CyclicNesterov, Adam

using ..AlphaZero

using Base: @kwdef
using Statistics: mean

"""
    AbstractNetwork

Abstract base type for a neural network.

# Constructor

Any subtype `Network` must implement `Base.copy` along with
the following constructor:

    Network(game_spec, hyperparams)

where the expected type of `hyperparams` is given by
[`HyperParams(Network)`](@ref HyperParams).
"""
abstract type AbstractNetwork end

#####
##### Interface
#####

"""
    HyperParams(::Type{<:AbstractNetwork})

Return the hyperparameter type associated with a given network type.
"""
function HyperParams end

"""
    hyperparams(::AbstractNetwork)

Return the hyperparameters of a network.
"""
function hyperparams end

"""
    game_spec(::AbstractNetwork)

Return the game specification that was passed to the network's constructor.
"""
function game_spec end

"""
    to_gpu(::AbstractNetwork)

Return a copy of the given network that has been transferred to the GPU
if one is available. Otherwise, return the given network untouched.
"""
function to_gpu end

"""
    to_cpu(::AbstractNetwork)

Return a copy of the given network that has been transferred to the CPU
or return the given network untouched if it is already on CPU.
"""
function to_cpu end

"""
    on_gpu(::AbstractNetwork) :: Bool

Test whether or not a network is located on GPU.
"""
function on_gpu end

"""
    set_test_mode!(mode=true)

Put a network in test mode or in training mode.
This is relevant for networks featuring layers such as
batch normalization layers.
"""
function set_test_mode! end

"""
    convert_input(::AbstractNetwork, input)

Convert an array (or number) to the right format so that it can be used
as an input by a given network.
"""
function convert_input end

function convert_input_tuple(nn::AbstractNetwork, input::Tuple)
  return map(input) do arr
    convert_input(nn, arr)
  end
end

"""
    convert_output(::AbstractNetwork, output)

Convert an array (or number) produced by a neural network
to a standard CPU array (or number) type.
"""
function convert_output end

function convert_output_tuple(nn::AbstractNetwork, output::Tuple)
  return map(output) do arr
    convert_output(nn, arr)
  end
end

"""
    forward(::AbstractNetwork, states)

Compute the forward pass of a network on a batch of inputs.

Expect a `Float32` tensor `states` whose batch dimension is the last one.

Return a `(P, V)` triple where:

  - `P` is a matrix of size `(num_actions, batch_size)`. It is allowed
    to put weight on invalid actions (see [`evaluate`](@ref)).
  - `V` is a row vector of size `(1, batch_size)`
"""
function forward end

"""
    regularized_params(::AbstractNetwork)

Return the collection of regularized parameters of a network.
This usually excludes neuron's biases.
"""
function regularized_params end

"""
    params(::AbstractNetwork)

Return the collection of trainable parameters of a network.
"""
function params end

"""
    gc(::AbstractNetwork)

Perform full garbage collection and empty the GPU memory pool.
"""
function gc end

# Optimizers and training

"""
    OptimiserSpec

Abstract type for an optimiser specification.
"""
abstract type OptimiserSpec end

"""
    CyclicNesterov(; lr_base, lr_high, lr_low, momentum_low, momentum_high)

SGD optimiser with a cyclic learning rate and cyclic Nesterov momentum.

  - During an epoch, the learning rate goes from `lr_low`
    to `lr_high` and then back to `lr_low`.
  - The momentum evolves in the opposite way, from high values
    to low values and then back to high values.
"""
@kwdef struct CyclicNesterov <: OptimiserSpec
  lr_base :: Float32
  lr_high :: Float32
  lr_low  :: Float32
  momentum_low :: Float32
  momentum_high :: Float32
end

"""
    Adam(;lr)

Adam optimiser.
"""
@kwdef struct Adam <: OptimiserSpec
  lr :: Float32
end

"""
    train!(callback, ::AbstractNetwork, opt::OptimiserSpec, loss, batches, n)

Update a given network to fit some data.
  - [`opt`](@ref OptimiserSpec) specifies which optimiser to use.
  - `loss` is a function that maps a batch of samples to a tracked real.
  - `data` is an iterator over minibatches.
  - `n` is the number of minibatches. If `length` is defined on `data`,
     we must have `length(data) == n`. However, not all finite
     iterators implement `length` and thus this argument is needed.
  - `callback(i, loss)` is called at each step with the batch number `i`
     and the loss on last batch.
"""
function train! end

#####
##### Derived functions
#####

"""
    num_parameters(::AbstractNetwork)

Return the total number of parameters of a network.
"""
function num_parameters(nn::AbstractNetwork)
  return sum(length(p) for p in params(nn))
end

"""
    num_regularized_parameters(::AbstractNetwork)

Return the total number of regularized parameters of a network.
"""
function num_regularized_parameters(nn::AbstractNetwork)
  return sum(length(p) for p in regularized_params(nn))
end

"""
    mean_weight(::AbstractNetwork)

Return the mean absolute value of the regularized parameters of a network.
"""
function mean_weight(nn::AbstractNetwork)
  sw = sum(sum(abs.(p)) for p in regularized_params(nn))
  sw = convert_output(nn, sw)
  return sw / num_regularized_parameters(nn)
end

"""
    forward_normalized(network::AbstractNetwork, states, actions_mask)

Evaluate a batch of vectorized states. This function is a wrapper
on [`forward`](@ref) that puts a zero weight on invalid actions.

# Arguments

  - `states` is a tensor whose last dimension has size `bach_size`
  - `actions_mask` is a binary matrix of size `(num_actions, batch_size)`

# Returned value

Return a `(P, V, Pinv)` triple where:

  - `P` is a matrix of size `(num_actions, batch_size)`.
  - `V` is a row vector of size `(1, batch_size)`.
  - `Pinv` is a row vector of size `(1, batch_size)`
     that indicates the total probability weight put by the network
     on invalid actions for each sample.

All tensors manipulated by this function have elements of type `Float32`.
"""
function forward_normalized(nn::AbstractNetwork, state, actions_mask)
  p, v = forward(nn, state)
  p = p .* actions_mask
  sp = sum(p, dims=1)
  p = p ./ (sp .+ eps(eltype(p)))
  p_invalid = 1 .- sp
  return (p, v, p_invalid)
end

to_singletons(x) = reshape(x, size(x)..., 1)
from_singletons(x) = reshape(x, size(x)[1:end-1])

"""
    evaluate(::AbstractNetwork, state)

    (nn::AbstractNetwork)(state) = evaluate(nn, state)

Evaluate the neural network as an MCTS oracle on a single state.

Note, however, that evaluating state positions once at a time is slow and so you
may want to use a `BatchedOracle` along with an inference server that uses
[`evaluate_batch`](@ref).
"""
function evaluate(nn::AbstractNetwork, state)
  gspec = game_spec(nn)
  actions_mask = GI.actions_mask(GI.init(gspec, state))
  x = GI.vectorize_state(gspec, state)
  a = Float32.(actions_mask)
  xnet, anet = to_singletons.(convert_input_tuple(nn, (x, a)))
  net_output = forward_normalized(nn, xnet, anet)
  p, v, _ = from_singletons.(convert_output_tuple(nn, net_output))
  return (p[actions_mask], v[1])
end

(nn::AbstractNetwork)(state) = evaluate(nn, state)

"""
    evaluate_batch(::AbstractNetwork, batch)

Evaluate the neural network as an MCTS oracle on a batch of states at once.

Take a list of states as input and return a list of `(P, V)` pairs as defined in the
MCTS oracle interface.
"""
function evaluate_batch(nn::AbstractNetwork, batch)
  gspec = game_spec(nn)
  X = Util.superpose((GI.vectorize_state(gspec, b) for b in batch))
  A = Util.superpose((GI.actions_mask(GI.init(gspec, b)) for b in batch))
  Xnet, Anet = convert_input_tuple(nn, (X, Float32.(A)))
  P, V, _ = convert_output_tuple(nn, forward_normalized(nn, Xnet, Anet))
  return [(P[A[:,i],i], V[1,i]) for i in eachindex(batch)]
end

"""
    copy(::AbstractNetwork; on_gpu, test_mode)

A copy function that also handles CPU/GPU transfers and
test/train mode switches.
"""
function copy(network::AbstractNetwork; on_gpu, test_mode)
  network = Base.copy(network)
  network = on_gpu ? to_gpu(network) : to_cpu(network)
  set_test_mode!(network, test_mode)
  return network
end

end


text:
############################################################################################
#                              WALL CONSTRAINT          								   #
############################################################################################
"""
	WallConstraint{P,T}

Constraint of the form

	   NO constraint violation
-----------------------------------------
              p1 * ///
                 |////
  NO constraint  |////        constraint
  violation      |------> v   violation
                 |////
                 |////
              p2 * ///
-----------------------------------------
	   NO constraint violation

`` (x - p1)^T v \\leq 0 ``
where ``x = [x[x], x[y]]``, ``p1 = [x1,y1]``, ``p2 = [x2,y2]`` ``v = [xv,yv]``.
The ordering of p1 and p2 does not matter.

# Constructor:
```julia
WallConstraint(n, x1::SVector{P}, y1::SVector{P}, x2::SVector{P}, y2::SVector{P}, xv::SVector{P}, yv::SVector{P}, x=1, y=2)
```
"""
struct WallConstraint{P,T} <: TrajectoryOptimization.StateConstraint
	n::Int
	x1::SVector{P,T}
	y1::SVector{P,T}
	x2::SVector{P,T}
	y2::SVector{P,T}
	xv::SVector{P,T}
	yv::SVector{P,T}
	x::Int  # index of x-state
	y::Int  # index of y-state
	function WallConstraint{P,T}(n::Int, x1::AbstractVector, y1::AbstractVector,
		x2::AbstractVector, y2::AbstractVector, xv::AbstractVector, yv::AbstractVector,
		x=1, y=2) where {P,T}
    	@assert length(x1) == length(y1) == length(x2) == length(y2) == length(xv) == length(yv) == P "Lengths of x1, y1, x2, y2, xv, yv must be equal. Got lengths ($(length(x1)), $(length(y1)), $(length(x2)), $(length(y2)), $(length(xv)), $(length(yv)))"
        new{P,T}(n, x1, y1, x2, y2, xv, yv, x, y)
    end
end
function WallConstraint(n::Int, x1::AbstractVector, y1::AbstractVector,
	x2::AbstractVector, y2::AbstractVector, xv::AbstractVector, yv::AbstractVector,
		x=1, y=2)
    T = promote_type(eltype(x1), eltype(y1), eltype(x2), eltype(y2), eltype(xv), eltype(yv))
    P = length(x1)
    WallConstraint{P,T}(n, x1, y1, x2, y2, xv, yv, x, y)
end
TrajectoryOptimization.state_dim(con::WallConstraint) = con.n

function TrajectoryOptimization.evaluate(con::WallConstraint, X::StaticVector)
	x1 = con.x1
	y1 = con.y1
	x2 = con.x2
	y2 = con.y2
	xv = con.xv
	yv = con.yv
	x = X[con.x]
	y = X[con.y]

	left  = (x .- x1) .* (x2 .- x1) + (y .- y1) .* (y2 .- y1) .> 0
	right = (x .- x2) .* (x1 .- x2) + (y .- y2) .* (y1 .- y2) .> 0
	out = (x .- x1) .* xv + (y .- y1) .* yv
	return out .* left .* right
end

function TrajectoryOptimization.jacobian!(∇c, con::WallConstraint{P}, X::SVector) where P
	x1 = con.x1
	y1 = con.y1
	x2 = con.x2
	y2 = con.y2
	xv = con.xv
	yv = con.yv
	x = X[con.x]
	y = X[con.y]

	left  = (x .- x1) .* (x2 .- x1) + (y .- y1) .* (y2 .- y1) .> 0
	right = (x .- x2) .* (x1 .- x2) + (y .- y2) .* (y1 .- y2) .> 0
	for i = 1:P
		∇c[i,con.x] = left[i] * right[i] * xv[i]
		∇c[i,con.y] = left[i] * right[i] * yv[i]
	end
	return false
end

@inline Base.length(::WallConstraint{P}) where P = P
@inline TrajectoryOptimization.sense(::WallConstraint) = Inequality()

function TrajectoryOptimization.change_dimension(con::WallConstraint, n::Int, m::Int, ix=1:n, iu=1:m)
	WallConstraint(n, con.x1, con.y1, con.x2, con.y2, con.xv, con.yv, ix[con.x], ix[con.y])
end

############################################################################################
#                              WALL3D CONSTRAINT          								   #
############################################################################################
"""
	Wall3DConstraint{P,T}

Constraint of the form

	   NO constraint violation
-----------------------------------------
  *  -------- p1 * ///
                 |////
  NO constraint  |////        constraint
  violation      |------> v   violation
                 |////
                 |////
p3 * -------- p2 * ///

                 -
	         -|  NO cons. viol.
      p3 *c.v.|
     -    |   |               -
 - |  cons.|  |           -
   |  viol. | |       -
   |         ||   -
p1 * -------- * p2
   |  c.v. -
   |  -
  -|     NO cons. viol.
-----------------------------------------
	   NO constraint violation
The points p1, p2, p3 form a parallelepiped in 3D space.

`` (x - p1)^T v \\leq 0 ``
where ``x = [x[x], x[y], x[z]]``, ``p1 = [x1,y1,z1]``, ``p2 = [x2,y2,z2]``, ``p3 = [x3,y3,z3]``, ``v = [xv,yv,zv]``.

# Constructor:
```julia
Wall3DConstraint(n, x1::SVector{P}, y1::SVector{P}, z1::SVector{P}, x2::SVector{P}, y2::SVector{P}, z2::SVector{P},  x3::SVector{P}, y3::SVector{P}, z3::SVector{P}, xv::SVector{P}, yv::SVector{P}, x=1, y=2, z=3)
```
"""
struct Wall3DConstraint{P,T} <: TrajectoryOptimization.StateConstraint
	n::Int
	x1::SVector{P,T}
	y1::SVector{P,T}
	z1::SVector{P,T}
	x2::SVector{P,T}
	y2::SVector{P,T}
	z2::SVector{P,T}
	x3::SVector{P,T}
	y3::SVector{P,T}
	z3::SVector{P,T}
	xv::SVector{P,T}
	yv::SVector{P,T}
	zv::SVector{P,T}
	x::Int  # index of x-state
	y::Int  # index of y-state
	z::Int  # index of z-state
	function Wall3DConstraint{P,T}(n::Int,
		x1::AbstractVector, y1::AbstractVector, z1::AbstractVector,
		x2::AbstractVector, y2::AbstractVector, z2::AbstractVector,
		x3::AbstractVector, y3::AbstractVector, z3::AbstractVector,
		xv::AbstractVector, yv::AbstractVector, zv::AbstractVector,
		x=1, y=2, z=3) where {P,T}
    	@assert length(x1) == length(y1) == length(z1) == length(x2) == length(y2) == length(z2) == length(x3) == length(y3) == length(z3) == length(xv) == length(yv) == length(zv) == P "Lengths of x1, y1, x2, y2, xv, yv must be equal. Got lengths ($(length(x1)), $(length(y1)), $(length(x2)), $(length(y2)), $(length(xv)), $(length(yv)))"
        new{P,T}(n, x1, y1, z1, x2, y2, z2, x3, y3, z3, xv, yv, zv, x, y, z)
    end
end
function Wall3DConstraint(n::Int,
		x1::AbstractVector, y1::AbstractVector, z1::AbstractVector,
		x2::AbstractVector, y2::AbstractVector, z2::AbstractVector,
		x3::AbstractVector, y3::AbstractVector, z3::AbstractVector,
		xv::AbstractVector, yv::AbstractVector, zv::AbstractVector,
		x=1, y=2, z=3)
    T = promote_type(
		eltype(x1), eltype(y1), eltype(z1),
		eltype(x2), eltype(y2), eltype(z2),
		eltype(x3), eltype(y3), eltype(z3),
		eltype(xv), eltype(yv), eltype(zv))
    P = length(x1)
    Wall3DConstraint{P,T}(n, x1, y1, z1, x2, y2, z2, x3, y3, z3, xv, yv, zv, x, y, z)
end
TrajectoryOptimization.state_dim(con::Wall3DConstraint) = con.n

function TrajectoryOptimization.evaluate(con::Wall3DConstraint, X::StaticVector)
	x1 = con.x1
	y1 = con.y1
	z1 = con.z1
	x2 = con.x2
	y2 = con.y2
	z2 = con.z2
	x3 = con.x3
	y3 = con.y3
	z3 = con.z3
	xv = con.xv
	yv = con.yv
	zv = con.zv
	x = X[con.x]
	y = X[con.y]
	z = X[con.z]

	left    = (x .- x1) .* (x2 .- x1) + (y .- y1) .* (y2 .- y1) + (z .- z1) .* (z2 .- z1) .> 0
	right   = (x .- x2) .* (x1 .- x2) + (y .- y2) .* (y1 .- y2) + (z .- z2) .* (z1 .- z2) .> 0
	bottom  = (x .- x3) .* (x2 .- x3) + (y .- y3) .* (y2 .- y3) + (z .- z3) .* (z2 .- z3) .> 0
	top     = (x .- x2) .* (x3 .- x2) + (y .- y2) .* (y3 .- y2) + (z .- z2) .* (z3 .- z2) .> 0
	out = (x .- x1) .* xv + (y .- y1) .* yv + (z .- z1) .* zv
	return out .* left .* right .* bottom .* top
end

function TrajectoryOptimization.jacobian!(∇c, con::Wall3DConstraint{P}, X::SVector) where P
	x1 = con.x1
	y1 = con.y1
	z1 = con.z1
	x2 = con.x2
	y2 = con.y2
	z2 = con.z2
	x3 = con.x3
	y3 = con.y3
	z3 = con.z3
	xv = con.xv
	yv = con.yv
	zv = con.zv
	x = X[con.x]
	y = X[con.y]
	z = X[con.z]

	left    = (x .- x1) .* (x2 .- x1) + (y .- y1) .* (y2 .- y1) + (z .- z1) .* (z2 .- z1) .> 0
	right   = (x .- x2) .* (x1 .- x2) + (y .- y2) .* (y1 .- y2) + (z .- z2) .* (z1 .- z2) .> 0
	bottom  = (x .- x3) .* (x2 .- x3) + (y .- y3) .* (y2 .- y3) + (z .- z3) .* (z2 .- z3) .> 0
	top     = (x .- x2) .* (x3 .- x2) + (y .- y2) .* (y3 .- y2) + (z .- z2) .* (z3 .- z2) .> 0
	for i = 1:P
		∇c[i,con.x] = left[i] * right[i] * bottom[i] * top[i] * xv[i]
		∇c[i,con.y] = left[i] * right[i] * bottom[i] * top[i] * yv[i]
		∇c[i,con.z] = left[i] * right[i] * bottom[i] * top[i] * zv[i]
	end
	return false
end

@inline Base.length(::Wall3DConstraint{P}) where P = P
@inline TrajectoryOptimization.sense(::Wall3DConstraint) = Inequality()

function TrajectoryOptimization.change_dimension(con::Wall3DConstraint, n::Int, m::Int, ix=1:n, iu=1:m)
	Wall3DConstraint(n,
		con.x1, con.y1, con.z1,
		con.x2, con.y2, con.z2,
		con.x3, con.y3, con.z3,
		con.x4, con.y4, con.z4,
		con.xv, con.yv, con.zv,
		ix[con.x], ix[con.y], ix[con.z])
end


text:
"""
    calc_grm()
---
Merge imputed data and calculate GRM.
e.g., two.bed + norge.bed, or, dutch.bed, german.bed, norge.bed
"""
function calc_grm(source)
    title("Data imputed within country")
    cd(work_dir)

    fra = joinpath(work_dir, "data/genotypes/step-8.plk")
    ped = joinpath(work_dir, "data/pedigree")
    tmp = joinpath(work_dir, "tmp")
    empty_dir(tmp)

    item("Merge source beds into one")
    open(joinpath(tmp, "beds.lst"), "w") do lst
        for data in source
            write(lst, joinpath(fra, data), '\n')
        end
    end
    _ = read(`plink --cow
                    --merge-list $tmp/beds.lst
                    --make-bed
                    --out $tmp/one`,
             String)
    done()

    cd(tmp)
    item("Prepare pedigree")
    pds = [joinpath(ped, "dutch.ped"),
           joinpath(ped, "german.ped")]
    _ = run(pipeline(`cat $pds`, `$rdBin/pedsort UUUUUUUUUUUUUUUUUUU`, "two.ped"))
    _ = run(pipeline(joinpath(ped, "norge.ped"), `$rdBin/pedsort 0`, "third.ped"))
    peds = [joinpath(tmp, "two.ped"),
            joinpath(tmp, "third.ped")]
    _ = run(pipeline(`$rdBin/merge-ped $peds`, joinpath(tmp, "all.ped")))
    done()

    item("ID country index")
    idic = Dict()
    for line in eachline(joinpath(tmp, "all.ped"))
        id, name = [split(line)[i] for i in [1, 4]]
        push!(idic, name => id)
    end

    icntry = 0
    cdic = Dict()
    countries = ["dutch", "german", "norge"]
    for country in countries
        icntry += 1
        for line in eachline(joinpath(fra, "$country.fam"))
            name = split(line)[2]
            push!(cdic, name=>icntry)
        end
    end

    # To put some MRY Dutch bulls in group 4

    for line in eachline(joinpath(fra, "dutch.fam"))
        name = split(line)[2]
        if name[1:3] == "MRY"
            cdic[name] = 4
        end
    end
    # fourth group ends here

    mv("one.fam", "tmp.fam", force=true)
    cindx = open("country.idx", "w")
    ifam  = open("one.fam", "w")
    for line in eachline("tmp.fam")
        t = split(line)
        name = t[2]
        t[2] = idic[name]
        write(ifam, join(t, ' '), '\n')
        tvec = [0, 0, 0, 0]
        tvec[cdic[name]] = 1
        write(cindx, idic[name], ' ', join(tvec, ' '), '\n')
    end
    close(cindx)
    close(ifam)
    done()

    item("Convert bed to 012 genotypes.txt")
    #plink_012("one", "tmp")
    _ = run(`plink
    		--cow
    		--bfile one
    		--extract $work_dir/data/genotypes/5160.raw
    		--recode A
    		--out tmp`
            )
    _ = run(pipeline("tmp.raw", `$rdBin/raw2gt`, "genotypes.txt"))
    done()

    item("Calculate GRM")
    nlc = length(split(readline("genotypes.txt"))[2])
    inp = ["$nlc",      "genotypes.txt",
           "genotypes", "4 country.idx",
           "vanraden",  "giv 0.00",
           "G ASReml",  "print_giv=asc",
           "print_geno=no genotypes.dat",
           "12", ""]
    write("calc_grm.inp", join(inp, '\n'))

    _ = run(`calc_grm`)
    _ = read(run(pipeline("G.grm", `gawk '{if($1==$2) print $3}'`, "diag.txt")), String)
    done()

    item("Parents offspring consistancy")
    open("int.ped", "w") do io
        for line in eachline("all.ped")
            id, pa, ma = split(line)
            write(io, "$id $pa $ma\n")
        end
    end

    inp = ["$nlc", "genotypes.txt int.ped", "genotypes", "1", "vanraden",
           "giv 0.00", "G ASReml", "print_giv=asc", "print_geno=no genotypes.dat",
           "12", "!checkparents", ""]
    mv("calc_grm.inp", "cal_grm.inp.G", force=true)
    write(joinpath(tmp, "calc_grm.inp"), join(inp, '\n'))
    _ = run(`calc_grm`)
    done()
    cd(work_dir)
end

"""
Calculate grm on MRY related ID in dutch pedigree
"""
function calc_grm_dutch_MRY()
    title("Data imputed within country")
    cd(work_dir)

    fra = joinpath(work_dir, "data/genotypes/step-8.plk")
    ped = joinpath(work_dir, "data/pedigree")
    tmp = joinpath(work_dir, "tmp")
    empty_dir(tmp)

    item("Prepare Dutch MRY releated data")
    _ = run(pipeline(joinpath(ped, "dutch.ped"),
                     `grep MRY`,
                     `gawk '{print $1}'`,
                     joinpath(tmp, "MRY.lst")))
    open(joinpath(tmp, "id.lst"), "w") do io
        for line in eachline(joinpath(tmp, "MRY.lst"))
            write(io, "0 $line\n")
        end
    end
    _ = run(pipeline(joinpath(ped, "dutch.ped"),
                     `$rdBin/pedsort UUUUUUUUUUUUUUUUUUU`,
                     joinpath(tmp, "dutch.ped")))
    _ = run(`plink
		--cow
		--bfile $fra/dutch
		--keep $tmp/id.lst
		--extract $work_dir/data/genotypes/5160.raw
		--make-bed
		--out $tmp/mry`
            )
    cd(tmp)
    ID = Dict()
    for line in eachline("dutch.ped")
        id, name = [split(line)[i] for i in [1, 4]]
        push!(ID, name=>id)
    end
    mv("mry.fam", "tmp.fam", force=true)
    open("mry.fam", "w") do io
        for line in eachline("tmp.fam")
            x = split(line)
            x[2] = ID[x[2]]
            write(io, join(x, ' '), '\n')
        end
    end
    plink_012("mry", "tmp")

    _ = run(pipeline("tmp.raw", `$rdBin/raw2gt`, "genotypes.txt"))
    done()

    item("Calculate GRM")
    nlc = length(split(readline("genotypes.txt"))[2])
    inp = ["$nlc",      "genotypes.txt",
           "genotypes", "1",
           "vanraden",  "giv 0.00",
           "G ASReml",  "print_giv=asc",
           "print_geno=no genotypes.dat",
           "12", ""]
    write("calc_grm.inp", join(inp, '\n'))

    _ = run(`calc_grm`)
    _ = read(run(pipeline("G.grm", `gawk '{if($1==$2) print $3}'`, "diag.txt")), String)
    done()

    item("Parents offspring consistancy")
    open("int.ped", "w") do io
        for line in eachline("dutch.ped")
            id, pa, ma = split(line)
            write(io, "$id $pa $ma\n")
        end
    end

    inp = ["$nlc", "genotypes.txt int.ped", "genotypes", "1", "vanraden",
           "giv 0.00", "G ASReml", "print_giv=asc", "print_geno=no genotypes.dat",
           "12", "!checkparents", ""]
    mv("calc_grm.inp", "cal_grm.inp.G", force=true)
    write(joinpath(tmp, "calc_grm.inp"), join(inp, '\n'))
    _ = run(`calc_grm`)
    done()
    cd(work_dir)
end

"""
    calc_grm_by_country()
---
Calculate **`G`** by country.  Also QQ-plot diagonals.  This must be run after `calc_grm`
to have the data ready.
"""
function calc_grm_by_country()
    title("Examine GRM by country")
    grm = joinpath(work_dir, "data/genotypes/grm")
    tmp = joinpath(work_dir, "tmp")
    bin = joinpath(work_dir, "bin")
    til = joinpath(work_dir, "data/genotypes/grm")

    cd(tmp)
    for country in ["dutch", "german", "norge"]
        item("G of $country")
        _ = run(pipeline(`$bin/merge4grm $country.raw`,
                         "genotypes.txt"))
        open("calc_grm.inp", "w") do io
            nlc = 0
            open(joinpath(tmp, "genotypes.txt"), "r") do gt
                line = readline(gt)
                nlc = length(split(line)) - 1
            end                     # it's 44037 anyway
            write(io, "$nlc\n")
            write(io, "genotypes.txt\n")
            write(io, "genotypes\n")
            write(io, "1\n")
            write(io, "vanraden\n")
            write(io, "giv 0.00\n")
            write(io, "G ASReml\n")
            write(io, "print_giv=asc\n")
            write(io, "print_geno=no genotypes.dat\n")
            write(io, "12\n")       # number of threads
        end
        _ = run(`calc_grm`)
        _ = read(run(pipeline("G.grm", `gawk '{if($1==$2) print $3}'`, "diag.txt")), String)

        isdir("$til/$country") || mkdir("$til/$country")
        for file in ["calc_grm.inp",
                     "calculated_all_freq.dat",
                     "country.idx",
                     "diag.txt",
                     "G_asreml.giv",
                     "genomic_inbr_coef.dat",
                     "genotypes.txt",
                     "G.grm",
                     "high_grm_coefs.log",
                     "ID.dic",
                     "ID_vs_row_number_G.txt",
                     "WARNING_ERROR_calc_grm.log"]
            isfile(file) && mv(file, "$grm/$country/$file", force=true)
        end
    end
    cd(work_dir)
end

"""
    plot_grm_diag()
---
Sort and plot diagonals of `GRM`.
"""
function plot_grm_diag()
    title("QQ-Plot GRM diagonals")
    fra = joinpath(work_dir, "data/genotypes/grm")
    pp = Any[]
    for country in ["all", "dutch", "german", "norge"]
        item("Figure of $country")
        dd = Float64[]
        for line in eachline(joinpath(fra, "$country/diag.txt"))
            push!(dd, parse(Float64, line))
        end
        p = plot(sort(dd), label = country, dpi=300)
        push!(pp, p)
    end
    plot(pp[1], pp[2], pp[3], pp[4], layout=4)
    savefig("notebooks/fig/diag.png")
end

"""
    calc_grm_w_ped()
---
Previously the elements in Dutch and Norwegian block are weird.  This funciton tries to bend the results
toward the **A** matrix, to see if the problem is alieved or not.
"""
function calc_grm_w_ped()
    title("Calculate GRM with pedigree")
    item("Prepare directories")
    countries = ["dutch", "german", "norge"]
    cd(work_dir)
    fra = joinpath(work_dir, "data/genotypes/step-8.plk")
    tmp = joinpath(work_dir, "tmp")
    bin = joinpath(work_dir, "bin")
    grm = joinpath(work_dir, "data/genotypes/grm")
    ped = joinpath(work_dir, "data/pedigree")
    isdir(grm) || mkdir(grm)
    empty_dir(tmp)
    done()

    item("Create unified pedigree")
    run(pipeline(joinpath(ped, "dg.ped"),
                 `$bin/pedsort UUUUUUUUUUUUUUUUUUU`,
                 joinpath(tmp, "dg.ped")))
    run(pipeline(joinpath(ped, "norge.ped"),
                 `$bin/pedsort 0`,
                 joinpath(tmp, "norge.ped")))
    argv = joinpath.(tmp, ["dg.ped", "norge.ped"])
    run(pipeline(`$bin/merge-ped $argv`, joinpath(tmp,"one.ped")))
    done()

    item("Convert country.bed to raw data")
    for country in countries
        plink_012(joinpath(fra, country), joinpath(tmp, country))
    end
    done()

    item("Merge raw data for `calc_grm`")
    argv = joinpath.(tmp, ["dutch.raw", "german.raw", "norge.raw", "country.idx"])
    _ = run(pipeline(joinpath(tmp, "one.ped"), `$bin/merge4grm $argv`, joinpath(tmp, "genotypes.txt")))
    done()

    item("Create a parameter file for `calc_grm`")
    nlc = begin
        open(joinpath(tmp, "genotypes.txt"), "r") do gt
            line = readline(gt)
            length(split(line)) - 1
        end
    end
    open(joinpath(tmp, "three.ped"), "w") do ped
        for line in eachline(joinpath(tmp, "one.ped"))
            id, pa, ma = split(line)[1:3]
            write(ped, "$id $pa $ma\n")
        end
    end

    inp = ["$nlc", "genotypes.txt three.ped", "genotypes", "1", "vanraden",
           "giv 0.00", "G ASReml", "print_giv=asc", "print_geno=no genotypes.dat",
           "12", "!checkparents", ""]
    write(joinpath(tmp, "calc_grm.inp"), join(inp, '\n'))
    done()

    item("Run `calc_grm`")
    cd(tmp)
    _ = run(`calc_grm`)
    done()

    cd(work_dir)
end

"""
This function create an ID (string) dictionary to (integer).
- Unknown ID are coded as 0
- Known ID start from 1
- ID in sire, dam columns and not in ID column are coded first
- Country were also recorded as 1, 2, 3 for dutch, german and norge, respectively.
"""
function generate_ID_dict()
    title("generate a ID dictionary of all ID in 3 countries")
    dic = open(joinpath(work_dir, "data/pedigree/ID.dict"), "w")
    ped = open(joinpath(work_dir, "data/pedigree/union.ped"), "w")
    idx = open(joinpath(work_dir, "data/pedigree/country.idx"), "w")
    ref = Dict()
    iid = 0

    write(dic, "UUUUUUUUUUUUUUUUUUU 0 0\n") # an unknown parent
    write(dic, "0 0 0\n")                   # another unknown parent
    push!(ref, "UUUUUUUUUUUUUUUUUUU" => 0)
    push!(ref, "0" => 0)
    breed = Dict("dutch" => " 1 0 0\n", "german" => " 0 1 0\n", "norge" => " 0 0 1\n")

    for country in ["dutch", "german", "norge"]
        id, pa, ma = begin
            a = String[]
            b = String[]
            c = String[]
            fra = joinpath(work_dir, "data/pedigree/$country.ped")
            for line in eachline(fra)
                x, y, z = (split(line)[i] for i in [1, 2, 3])
                push!(a, x)
                push!(b, y)
                push!(c, z)
            end
            Set(a), delete!(Set(b), "UUUUUUUUUUUUUUUUUUU"), delete!(Set(c), "UUUUUUUUUUUUUUUUUUU")
        end

        println("N id pa ma in $country: ", length(id), ' ', length(pa), ' ', length(ma))
        for i in setdiff(union(pa, ma), id)
            iid += 1
            write(dic, "$i $iid $country\n")
            push!(ref, i => iid)
            write(ped, iid, " 0 0\n")
            write(idx, iid, breed[country])
        end
        for line in eachline(joinpath(work_dir, "data/pedigree/$country.ped"))
            x, y, z = (split(line)[i] for i in [1, 2, 3])
            iid += 1
            push!(ref, x => iid)
            write(dic, "$x $iid $country\n")
            write(ped, iid, ' ', ref[y], ' ', ref[z], '\n')
            write(idx, iid, breed[country])
        end
    end

    close(ped, idx)
    done()
end

function test_ID()
    title("Test ID")
    #til = joinpath
end
#= Discarded codes
    allid = begin
        id = String[]
        for x in eachline(joinpath(work_dir, "data/pedigree/ids"))
            push!(id, x)
        end
        Set(id)
    end
    for country in ["dutch", "german", "norge"]
        cid = begin
            id = String[]
            for line in eachline(joinpath(work_dir, "data/genotypes/step-8.plk/$country.fam"))
                push!(id, split(line)[2])
            end
            Set(id)
        end
        did = setdiff(cid, allid)
        if length(did) > 0
            println(country, " has ", length(did), " ID not in the pedigree")
            write(joinpath(work_dir, "tmp/$country.diff"), join(did, '\n'), '\n')
        end
    end
=#


text:
# Read the following file containing julia code
include("fancymodule.jl")

# Import the module fancy
import .fancy
# You can now call all objects in fancy with fancy.object
fancy.say_hello()
fancy.unexported_function()
fancy.prices
# This doesnt work
say_hello()

# "using" makes the exported objects accessible in the Main module
using .fancy
# Now it works
say_hello()
prices
# But this still doesn't
unexported_function()

# More information:
# https://docs.julialang.org/en/v1/manual/workflow-tips/


text:
# deprecated in GeoJSON v0.5
@deprecate parse read
@deprecate parsefile(path) read(Base.read(path, String))
@deprecate geojson write


text:
module AssociationRuleMining

using Formatting
using Base.Threads
using DataFrames

import Base.==
import Base.unique
import Base.display
import Base.isempty
import Base.find

# export Rule, Sequence, IDList, PrefixNode, SeqRule, apriori, frequent, make_transactions, spade, make_sequences, sequential_rules, dataset

# include("apriori.jl")
include("apriori_views.jl")     # NOTE: This is broken, but maybe a good idea

include("apriori_utils.jl")
include("parallel_tools.jl")
include("spade.jl")
include("spade_utils.jl")
include("sequential_rules.jl")

# end # module


text:
struct Krr
    X::Array{Float64}
    K::Array{Float64}
    U::Array{Float64}
    UtDY::Array{Float64}
    sv::Vector{Float64}
    D::Array{Float64}
    sqrtD::Array{Float64}
    DKt::Array{Float64}
    vtot::Array{Float64}
    lb::Float64
    ymeans::Vector{Float64}
    weights::Vector{Float64}
    kern
    dots
end

"""
    krr(X, Y, weights = ones(size(X, 1)); lb = .01, kern = "krbf", kwargs...)
Kernel ridge regression (KRR) implemented by SVD factorization.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.
* `lb` : A value of the regularization parameter "lambda".
* 'kern' : Type of kernel used to compute the Gram matrices.
    Possible values are "krbf" of "kpol" (see respective functions `krbf` and `kpol`.
* `kwargs` : Named arguments to pass in the kernel function.

KRR is also referred to as least squared SVM regression (LS-SVMR).
The method is close to the particular case of SVM regression
where there is novmarges excluding the observations (epsilon coefficient
set to zero). The difference is that a L2-norm optimization is done,
instead of L1 in SVM.

The kernel Gram matrices are internally centered.

## References

Bennett, K.P., Embrechts, M.J., 2003. An optimization perspective on kernel partial least squares regression,
in: Advances in Learning Theory: Methods, Models and Applications,
NATO Science Series III: Computer & Systems Sciences. IOS Press Amsterdam, pp. 227-250.

Cawley, G.C., Talbot, N.L.C., 2002. Reduced Rank Kernel Ridge Regression.
Neural Processing Letters 16, 293-302. https://doi.org/10.1023/A:1021798002258

Krell, M.M., 2018. Generalizing, Decoding, and Optimizing Support Vector Machine Classification.
arXiv:1801.04929.

Saunders, C., Gammerman, A., Vovk, V., 1998. Ridge Regression Learning Algorithm in Dual Variables,
in: In Proceedings of the 15th International Conference on Machine Learning. Morgan Kaufmann, pp. 515â521.

Suykens, J.A.K., Lukas, L., Vandewalle, J., 2000. Sparse approximation using
least squares support vector machines. 2000 IEEE International Symposium on Circuits and Systems.
Emerging Technologies for the 21st Century. Proceedings (IEEE Cat No.00CH36353). https://doi.org/10.1109/ISCAS.2000.856439

Welling, M., n.d. Kernel ridge regression. Department of Computer Science,
University of Toronto, Toronto, Canada. https://www.ics.uci.edu/~welling/classnotes/papers_class/Kernel-Ridge.pdf

"""
function krr(X, Y, weights = ones(size(X, 1)); lb = .01, kern = "krbf", kwargs...)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    ymeans = colmeans(Y, weights)
    fkern = eval(Meta.parse(kern))
    K = fkern(X, X; kwargs...)
    D = Diagonal(weights)
    DKt = D * K'
    vtot = sum(DKt, dims = 1)
    Kc = K .- vtot' .- vtot .+ sum(D * DKt')
    # Kd = D^(1/2) * Kc * D^(1/2)
    #    = U * Delta^2 * U'
    sqrtD = sqrt.(D)
    Kd = sqrtD * Kc * sqrtD
    res = LinearAlgebra.svd(Kd)
    U = res.V
    sv = sqrt.(res.S)
    # UtDY = U' * D^(1/2) * Y
    UtDY = U' * sqrtD * Y
    Krr(X, K, U, UtDY, sv, D, sqrtD, DKt,
        vtot, lb, ymeans, weights, kern, kwargs)
end

"""
    coef(object::Krr; lb = nothing)
Compute the b-coefficients of a fitted model.
* `object` : The fitted model.
* `lb` : A value of the regularization parameter "lambda".
    If nothing, it is the parameter stored in the fitted model.
"""
function coef(object::Krr; lb = nothing)
    isnothing(lb) ? lb = object.lb : nothing
    eig = object.sv.^2
    z = 1 ./ (eig .+ lb^2)
    A = object.U * (Diagonal(z) * object.UtDY)
    q = length(object.ymeans)
    int = reshape(object.ymeans, 1, q)
    tr = sum(eig .* z)
    (A = A, int = int, df = 1 + tr)
end

"""
    predict(object::Krr, X; lb = nothing)
Compute Y-predictions from a fitted model.
* `object` : The maximal fitted model.
* `X` : X-data for which predictions are computed.
* `lb` : Regularization parameter, or collection of regularization parameters, "lambda" to consider.
    If nothing, it is the parameter stored in the fitted model.
"""
function predict(object::Krr, X; lb = nothing)
    isnothing(lb) ? lb = object.lb : nothing
    fkern = eval(Meta.parse(object.kern))
    K = fkern(X, object.X; object.dots...)
    DKt = object.D * K'
    vtot = sum(DKt, dims = 1)
    Kc = K .- vtot' .- object.vtot .+ sum(object.D * object.DKt')
    le_lb = length(lb)
    pred = list(le_lb, Matrix{Float64})
    @inbounds for i = 1:le_lb
        z = coef(object; lb = lb[i])
        pred[i] = z.int .+ Kc * (object.sqrtD * z.A)
    end
    le_lb == 1 ? pred = pred[1] : nothing
    (pred = pred,)
end


text:
import .FileIO

function loadtodict!(d::Dict, g::Union{File,Group}, prefix::String="")
    for k in keys(g)
        v = g[k]
        if v isa Group
            loadtodict!(d, v, prefix * k * "/")
        else
            d[prefix * k] = read(g, k)
        end
    end
    return d
end

# load with just a filename returns a flat dictionary containing all the variables
function fileio_load(f::FileIO.File{FileIO.format"HDF5"}; kwargs...)
    h5open(FileIO.filename(f), "r"; kwargs...) do file
        loadtodict!(Dict{String,Any}(), file)
    end
end

# when called with explicitly requested variable names, return each one
function fileio_load(f::FileIO.File{FileIO.format"HDF5"}, varname::AbstractString; kwargs...)
    h5open(FileIO.filename(f), "r"; kwargs...) do file
        read(file, varname)
    end
end

function fileio_load(f::FileIO.File{FileIO.format"HDF5"}, varnames::AbstractString...; kwargs...)
    h5open(FileIO.filename(f), "r"; kwargs...) do file
        map(var -> read(file, var), varnames)
    end
end

# save all the key-value pairs in the dict as top-level variables
function fileio_save(f::FileIO.File{FileIO.format"HDF5"}, dict::AbstractDict; kwargs...)
    h5open(FileIO.filename(f), "w"; kwargs...) do file
        for (k,v) in dict
            if !isa(k, AbstractString)
                throw(ArgumentError("keys must be strings (the names of variables), got $k"))
            end
            write(file, String(k), v)
        end
    end
end


text:
using Genie.Router
using Jockey: JockeyPoll
using Jockey: JSession, JSettings
#==== Jockey Poll Statuses:
critical
request
warn
use as symbols, such as on "/"
====#
# Here we need to load server configuration, assign it into a

route("/") do
    serve_static_file("index.html")
end


text:
using ITensors,
      ITensorsGPU,
      LinearAlgebra, # For tr()
      Combinatorics, # For permutations()
      CuArrays,
      Test

      # gpu tests!
@testset "cuITensor, Dense{$SType} storage" for SType ∈ (Float64, ComplexF64)
  mi,mj,mk,ml,ma = 2,3,4,5,6,7
  i = Index(mi,"i")
  j = Index(mj,"j")
  k = Index(mk,"k")
  l = Index(ml,"l")
  a = Index(ma,"a")
  @testset "Constructor" begin
      A = cuITensor(one(SType),i,j,k)
      @test collect(CuArray(A, i, j, k)) == ones(SType, dim(i), dim(j), dim(k))
      A = randomCuITensor(IndexSet(i,j,k))
      @test inds(A) == IndexSet(i, j, k)
      @test ITensorsGPU.store(A) isa ITensorsGPU.CuDense
      Aarr = rand(SType, dim(i)*dim(j)*dim(k))
      @test collect(ITensor(Aarr, i, j, k)) == collect(cuITensor(Aarr, i, j, k))
      @test cuITensor(SType, i, j, k) isa ITensor{3}
      @test store(cuITensor(SType, i, j, k)) isa ITensorsGPU.CuDense{SType}
      @test vec(collect(CuArray(ITensor(Aarr, i,j,k), i, j, k))) == Aarr
  end
  @testset "Test permute(cuITensor,Index...)" begin
    CA     = randomCuITensor(SType,i,k,j)
    permCA = permute(CA,k,j,i)
    permA  = collect(permCA)
    @test k==inds(permA)[1]
    @test j==inds(permA)[2]
    @test i==inds(permA)[3]
    A = collect(CA)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k(kk),i(ii),j(jj)]==permA[i(ii),j(jj),k(kk)]
    end
  end
  #=@testset "Test scalar(cuITensor)" begin
    x = SType(34)
    A = randomCuITensor(a)
    @test x==scalar(A)
  end=#
  @testset "Test CuVector(cuITensor)" begin
      v = CuVector(ones(SType, dim(a)))
      A = cuITensor(v, a)
      @test v==CuVector(A)
  end
  @testset "Test CuMatrix(cuITensor)" begin
      v = CuMatrix(ones(SType, dim(a), dim(l)))
      A = cuITensor(vec(v), a, l)
      @test v==CuMatrix(A, a, l)
      A = cuITensor(vec(v), a, l)
      @test v==CuMatrix(A)
      A = cuITensor(vec(v), a, l)
      @test v==CuArray(A, a, l)
      @test v==CuArray(A)
  end
  @testset "Test norm(cuITensor)" begin
    A = randomCuITensor(SType,i,j,k)
    B = dag(A)*A
    @test norm(A)≈sqrt(real(scalar(B)))
  end
  @testset "Test complex(cuITensor)" begin
    A  = randomCuITensor(SType,i,j,k)
    cA = complex(A)
    @test complex.(CuArray(A)) == CuArray(cA)
  end
  #@testset "Test exp(cuITensor)" begin
  #  A  = randomCuITensor(SType,i,i')
  #  @test CuArray(exp(A,i,i')) ≈ exp(CuArray(A))
  #end
  @testset "Test add cuITensors" begin
    dA = randomCuITensor(SType,i,j,k)
    dB = randomCuITensor(SType,k,i,j)
    A = collect(dA)
    B = collect(dB)
    C = collect(dA+dB)
    @test CuArray(permute(C,i,j,k))==CuArray(permute(A,i,j,k))+CuArray(permute(B,i,j,k))
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test C[i(ii),j(jj),k(kk)]==A[j(jj),i(ii),k(kk)]+B[i(ii),k(kk),j(jj)]
    end
  end

  @testset "Test factorizations of a cuITensor" begin

    A = randomCuITensor(SType,i,j,k,l)

    @testset "Test SVD of a cuITensor" begin
      U,S,V = svd(A,(j,l))
      u = commonind(U,S)
      v = commonind(S,V)
      @test collect(A)≈collect(U*S*dag(V))
      @test collect(U*dag(prime(U,u)))≈δ(SType,u,u') rtol=1e-14
      @test collect(V*dag(prime(V,v)))≈δ(SType,v,v') rtol=1e-14
    end

    #=@testset "Test SVD truncation" begin
        M = randn(4,4)
        (U,s,V) = svd(M)
        ii = Index(4)
        jj = Index(4)
        S = Diagonal(s)
        T = cuITensor(vec(CuArray(U*S*V')),IndexSet(ii,jj))
        (U,S,V) = svd(T,ii;maxm=2)
        @test norm(U*S*V-T)≈sqrt(s[3]^2+s[4]^2)
    end=#

    @testset "Test QR decomposition of a cuITensor" begin
      Q,R = qr(A,(i,l))
      q = commonind(Q,R)
      @test collect(A)≈collect(Q*R)
      @test collect(Q*dag(prime(Q,q)))≈δ(SType,q,q') atol=1e-14
    end

    @testset "Test polar decomposition of a cuITensor" begin
      U,P = polar(A,(k,l))
      @test collect(A)≈collect(U*P)
    end

  end # End ITensor factorization testset
end # End Dense storage test


text:
import BenchmarkTools

function set_quick_params!(bench)
    bench.params.seconds = 0.001
    bench.params.evals = 1
    bench.params.samples = 1
    bench.params.gctrial = false
    bench.params.gcsample = false
    return bench
end

foreach_benchmark(f!, bench::BenchmarkTools.Benchmark) = f!(bench)
function foreach_benchmark(f!, group::BenchmarkTools.BenchmarkGroup)
    for x in values(group)
        foreach_benchmark(f!, x)
    end
end

function quick!(suite)
    foreach_benchmark(set_quick_params!, suite)
    return suite
end

function maybe_quick!(suite)
    if lowercase(get(ENV, "QUICK", "false")) == "true"
        @info "Use quick-run benchmark parameters"
        quick!(suite)
    end
    return suite
end


text:
import Dates
using Compat

@compat abstract type RDate end
# Calendars
@compat abstract type CalendarManager end
@compat abstract type Calendar end

"""
    HolidayRoundingConvention

The convention used in conjunction with holiday calendars to determine what
to do if an adjustments falls on a holiday.
"""
@compat abstract type HolidayRoundingConvention end

"""
    InvalidDayConvention

The convention used for handling month or year increments which fall on a day
which is not valid.
"""
@compat abstract type InvalidDayConvention end

"""
    MonthIncrementConvention

The convention used for handling month or year increments and how they should
be applied. Should handle invalid days explicitly, as these wil be handled
separately by the InvalidDayConvention
"""
@compat abstract type MonthIncrementConvention end

"""
    NullCalendarManager()

The most primitive calendar manager, that will return an error for any request to
get a calendar. The default calendar manager that is available when applying rdates
using the + operator, without an explicit calendar manager.
"""
struct NullCalendarManager <: CalendarManager end

"""
    is_holiday(calendar::Calendar, date::Dates.Date)::Bool

Determine whether the date requested is a holiday for this calendar or not.
"""
is_holiday(x::Calendar, ::Dates.Date) = error("$(typeof(x)) does not support is_holiday")

"""
    holidays(calendar::Calendar, from::Dates.Date, to::Dates.Date)::Vector{Bool}

Get the vector of each day and whether its a holiday or not, inclusive of from and to.
"""
holidays(::C, ::Dates.Date, ::Dates.Date) where {C<:Calendar} = error("holidays not implemented by $C")

"""
    holidaycount(calendar::Calendar, from::Dates.Date, to::Dates.Date)::Int

Get the number of holidays in the calendar between two dates (inclusive)
"""
holidaycount(cal::Calendar, from::Dates.Date, to::Dates.Date) = Base.sum(holidays(cal, from, to))

"""
    bizdaycount(calendar::Calendar, from::Dates.Date, to::Dates.Date)::Int

Get the number of business days in the calendar between two dates (inclusive)
"""
bizdaycount(cal::Calendar, from::Dates.Date, to::Dates.Date) = 1 + (to-from).value - holidaycount(cal, from, to)

"""
    calendar(calendarmgr::CalendarManager, names::Vector)::Calendar

Given a set of calendar names, request the calendar manager to retrieve the associated
calendar that supports the union of them.
"""
calendar(x::CalendarManager, names::Vector)::Calendar = error("$(typeof(x)) does not support calendar")
calendar(x::CalendarManager, name::String)::Calendar = calendar(x, split(name, "|"))

"""
    apply(rdate::RDate, date::Dates.Date, calendarmgr::CalendarManager)::Dates.Date

The application of an rdate to a specific date, given an explicit calendar manager.

### Examples
```jula-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(rd"1b@WEEKEND", Date(2021,7,9), cal_mgr)
2021-07-12
```
"""
apply(rd::RDate, ::Dates.Date, ::CalendarManager)::Dates.Date = error("$(typeof(rd)) does not support date apply")

"""
    multiply(rdate::RDate, count::Integer)::RDate

An internalised multiplication of an rdate which is generated without reapplication of a relative date multiple
times. To apply an rdate multiple times, then use a `Repeat`.

For example "6m" + "6m" != "12m" for all dates, due to the fact that there are different days in each month
and invalid day conventions will kick in.
"""
multiply(rd::R, count::Integer) where {R<:RDate} = error("multiply not implemented by $R")
Base.:*(rdate::RDate, count::Integer) = multiply(rdate, count)
Base.:*(count::Integer, rdate::RDate) = multiply(rdate, count)

"""
    apply(rdate::RDate, date::Dates.Date)::Dates.Date

The application of an rdate to a specific date, without an explicit calendar manager. This
will use the `NullCalendarManager()`.
"""
apply(rd::RDate, date::Dates.Date) = apply(rd, date, NullCalendarManager())
Base.:+(rd::RDate, date::Dates.Date) = apply(rd, date)
Base.:+(date::Dates.Date, rd::RDate) = apply(rd, date)
Base.:-(date::Dates.Date, rd::RDate) = apply(-rd, date)
Base.:-(x::RDate)::RDate = error("$(typeof(x)) does not support negation")

"""
    apply(rounding::HolidayRoundingConvention, date::Dates.Date, calendar::Calendar)::Dates.Date

Apply the holiday rounding to a given date. There is no strict requirement that the resolved date
will not be a holiday for the given date.
"""
apply(x::HolidayRoundingConvention, date::Dates.Date, calendar::Calendar) = error("$(typeof(x)) does not support apply")

"""
    apply(rounding::InvalidDayConvention, day::Integer, month::Integer, year::Integer)::Dates.Date

Given a day, month and year which do not generate a valid date, adjust them in some form to a valid date.
"""
adjust(x::InvalidDayConvention, day, month, year) = error("$(typeof(x)) does not support adjust")

"""
    apply(rounding::MonthIncrementConvention, from::Dates.Date, new_month::Integer, new_year::Integer, calendar_mgr::CalendarManager)::Tuple{Integer, Integer, Integer}

Given the initial day, month and year that we're moving from and a generated new month and new year, determine the new
day, month and year that should be used. Generally used to handle specific features around preservation of month ends
or days of week, if required.

Note that the generated (day, month, year) do not need to be able to produce a valid date. The invalid day convention
should be applied, if required, after this calculation.

May use a calendar manager if required as well
"""
adjust(x::MonthIncrementConvention, from::Dates.Date, new_month, new_year, calendar_mgr::CalendarManager) = error("$(typeof(x)) does not support adjust")


text:
module AlgamesDriving

greet() = print("Hello World!")

using Algames
using BenchmarkTools
using Colors
using CoordinateTransformations
using FileIO
using GeometryBasics
using LinearAlgebra
using MeshCat
using MeshIO
using Parameters
using REPL
using Random
using Rotations
using StaticArrays

# Struct
# vehicle_state
export
    VehicleState,
    standardize,
    specialize
# starting_area
export
    StartingArea,
    randstate
# lane
export
    CircularWall,
    Lane
# roadway
export
    Roadway,
    RoadwayOptions,
    HighwayRoadwayOptions,
    MergingRoadwayOptions,
    FourIntersectionRoadwayOptions,
    build_roadway,
    add_lane!
# player
export
    Player

# Sceanario
# scenario
export
    Scenario,
    get_state,
    get_control_bounds

# Visualization
# visualize_roadway
export
    VisualizationOptions,
    visualize_roadway!
# visualization_methods
export
    clean!,
    ModifiedMeshFileObject,
    rescale_contents,
    select_material,
    transRot
# Camera
export
    set_env!,
    set_camera!,
    set_camera_birdseye!
# visualize_player
export
    PlayerVisualizationOptions,
    WaypointVisualizationOptions,
    LineVisualizationOptions,
    set_player!,
    build_waypoint!,
    set_line_traj!
# visualize_roadway
export
    RoadwayVisualizationOptions,
    set_roadway!
# visualize_scenario
export
    set_scenario!
# visualize_trajectory
export
    set_state!,
    set_traj!,
    set_waypoint!,
    set_waypoint_traj!

# MPC
# Game
export
    get_char
# mpc_struct
export
    MPCOptions,
    MPCStatistics,
    TrajStatistics
# mpc_methods
export
    first_order_rollout,
    simulate_MPC!

# # Plots
# export

# Struct
include("struct/vehicle_state.jl")
include("struct/starting_area.jl")
include("struct/lane.jl")
include("struct/roadway.jl")
include("struct/highway_roadway.jl")
include("struct/merging_roadway.jl")
include("struct/four_intersection_roadway.jl")
include("struct/player.jl")

# Scenario
include("scenario/scenario.jl")

# Visualization
include("visualization/visualization_methods.jl")
include("visualization/camera.jl")
include("visualization/vehicle/visualize_roadway.jl")
include("visualization/vehicle/visualize_player.jl")
include("visualization/vehicle/visualize_scenario.jl")
include("visualization/vehicle/visualize_trajectory.jl")
include("visualization/drone/static_object.jl")
include("visualization/drone/constraint.jl")
include("visualization/drone/animation.jl")
include("visualization/drone/bicycle.jl")
include("visualization/drone/double_integrator.jl")
include("visualization/drone/quadrotor.jl")
include("visualization/drone/unicycle.jl")

# MPC
include("mpc/mpc_struct.jl")
include("mpc/mpc_methods.jl")
include("mpc/game.jl")

# Plots

# Utils
include("utils.jl")

end # module


text:
function loss(a, b)
    return Flux.mse(a, b)
end

@views split_u(uvT, Nz) = uvT[1:Nz, :]
@views split_v(uvT, Nz) = uvT[Nz+1:2Nz, :]
@views split_T(uvT, Nz) = uvT[2Nz+1:3Nz, :]

@views ∂_∂z(profile, D_face) = hcat([D_face * profile[:,i] for i in 1:size(profile, 2)]...)

function calculate_loss_scalings(losses, fractions, train_gradient)
    velocity_scaling = (1 - fractions.T) / fractions.T * losses.T / (losses.u + losses.v)
    profile_loss = velocity_scaling * (losses.u + losses.v) + losses.T

    if train_gradient
        velocity_gradient_scaling = (1 - fractions.∂T∂z) / fractions.∂T∂z * losses.∂T∂z / (losses.∂u∂z + losses.∂v∂z)
        gradient_loss = velocity_gradient_scaling * (losses.∂u∂z + losses.∂v∂z) + losses.∂T∂z
        total_gradient_scaling = (1 - fractions.profile) / fractions.profile * profile_loss / gradient_loss
    else
        velocity_gradient_scaling = 0
        gradient_loss = 0
        total_gradient_scaling = 0
    end

    return (   u = velocity_scaling,
               v = velocity_scaling,
               T = 1,
            ∂u∂z = total_gradient_scaling * velocity_gradient_scaling,
            ∂v∂z = total_gradient_scaling * velocity_gradient_scaling,
            ∂T∂z = total_gradient_scaling )
end

function apply_loss_scalings(losses, scalings)
    return (
        u = scalings.u .* losses.u,
        v = scalings.v .* losses.v,
        T = scalings.T .* losses.T,
        ∂u∂z = scalings.∂u∂z .* losses.∂u∂z,
        ∂v∂z = scalings.∂v∂z .* losses.∂v∂z,
        ∂T∂z = scalings.∂T∂z .* losses.∂T∂z,
    )
end

@views function loss_per_tstep(a, b)
    return [loss(a[:,i], b[:,i]) for i in 1:size(a, 2)]
end

function plot_loss(losses::NamedTuple, OUTPUT_PATH)
    profile_losses = losses.u .+ losses.v .+ losses.T
    gradient_losses = losses.∂u∂z .+ losses.∂v∂z .+ losses.∂T∂z

    x = 1:length(losses.u)
    fig = Figure(resolution=(1500, 750))
    ax_top = Axis(fig[1, 1], yscale=log10)
    ax_bot = Axis(fig[1, 2], yscale=log10)

    ax_top.xlabel = "Epochs"
    ax_bot.xlabel = "Epochs"

    ax_top.ylabel = "Losses"
    ax_bot.ylabel = "Losses"

    colors = distinguishable_colors(8, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

    u_line = CairoMakie.lines!(ax_top, x, losses.u, color=colors[1])
    v_line = CairoMakie.lines!(ax_top, x, losses.v, color=colors[2])
    T_line = CairoMakie.lines!(ax_top, x, losses.T, color=colors[3])

    ∂u∂z_line = CairoMakie.lines!(ax_top, x, losses.∂u∂z, color=colors[4])
    ∂v∂z_line = CairoMakie.lines!(ax_top, x, losses.∂v∂z, color=colors[5])
    ∂T∂z_line = CairoMakie.lines!(ax_top, x, losses.∂T∂z, color=colors[6])

    profile_line = CairoMakie.lines!(ax_bot, x, profile_losses, color=colors[7])
    gradient_line = CairoMakie.lines!(ax_bot, x, gradient_losses, color=colors[8])

    legend_individual = fig[2,1] = Legend(fig, [u_line, v_line, T_line, ∂u∂z_line, ∂v∂z_line, ∂T∂z_line],
                                ["u", "v", "T", "∂u/∂z", "∂v/∂z", "∂T/∂z"], orientation=:horizontal)

    legend_profile = fig[2,2] = Legend(fig, [profile_line, gradient_line],
                                ["profile", "gradient"], orientation=:horizontal)

    colsize!(fig.layout, 1, CairoMakie.Relative(0.5))
    colsize!(fig.layout, 2, CairoMakie.Relative(0.5))

    rowsize!(fig.layout, 2, CairoMakie.Relative(0.05))

    trim!(fig.layout)

    save(OUTPUT_PATH, fig, pt_per_unit=4, px_per_unit=4)
end

function plot_loss(losses, OUTPUT_PATH)
    x = 1:length(losses)
    fig = Figure()
    ax = Axis(fig[1, 1], yscale=log10)

    ax.xlabel = "Epochs"
    ax.ylabel = "Total Loss"

    line = CairoMakie.lines!(ax, x, losses)

    trim!(fig.layout)

    save(OUTPUT_PATH, fig)
end


text:
#   Linear trapezoidal, Dose 100, Dosetime 0, no tau

refdict = Dict(
:Cmax => [
    190.869
    261.177
    105.345
    208.542
    169.334
    154.648
    153.254
    138.327
    167.347
    125.482
],
:Tmax => [
    1
    1
    1.5
    1
    4
    2.5
    2.5
    4
    3
    2
],
:Cdose => [
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
],
:Clast => [
    112.846
    85.241
    67.901
    97.625
    110.778
    69.501
    58.051
    74.437
    93.44
    42.191
],
:AUClast => [
    9585.4218
    10112.176
    5396.5498
    9317.8358
    9561.26
    6966.598
    7029.5735
    7110.6745
    8315.0803
    5620.8945
],
:AUMClast => [
    333582.48
    298701.39
    186032.06
    313955.9
    315181.56
    226977.06
    219797.71
    240526.05
    277613.98
    154893.06
],
:AUCall => [
    9585.4218
    10112.1760
    5396.5498
    9317.8358
    9561.2600
    6966.5980
    7029.5735
    7110.6745
    8315.0803
    5620.8945
],
:Rsq => [
    0.78607696
    0.99276359
    0.81358898
    0.91885869
    0.85335995
    0.95011904
    0.97031231
    0.94796904
    0.94753789
    0.88092269
],
:ARsq => [
    0.71476928
    0.99035145
    0.77630678
    0.83771737
    0.82891994
    0.92517856
    0.96041642
    0.92195356
    0.92130684
    0.86391165
],
:Kel => [
    0.0033847439
    0.014106315
    0.0032914304
    0.0076953442
    0.0068133279
    0.0076922807
    0.012458956
    0.0089300798
    0.0056458649
    0.017189737
],
:HL => [
    204.78571
    49.137367
    210.59148
    90.073577
    101.73401
    90.109450057666
    55.634451
    77.619371
    122.77077
    40.323315
],
:Clast_pred => [
    117.30578
    82.53669
    66.931057
    100.76793
    105.29832
    71.939942
    61.172702
    75.604277
    93.761762
    38.810857
],
:AUCinf => [
    42925.019
    16154.93
    26026.183
    22004.078
    25820.275
    16001.76
    11688.953
    15446.21
    24865.246
    8075.3242
],
:AUCpct => [
    77.669383
    37.405019
    79.26492
    57.65405
    62.969953
    56.463551
    39.861391
    53.964925
    66.559429
    30.394194
],
:MRTlast => [
    34.801023
    29.538786
    34.472406
    33.69408
    32.964438
    32.58076
    31.267574
    33.826053
    33.386807
    27.556657
],
:MRTinf => [
    293.16224
    71.937917
    305.04073
    130.69968
    149.96684
    128.24114
    79.498252
    114.8571
    176.97811
    58.746446
],
:Clinf => [
    0.0023296437
    0.0061900608
    0.0038422846
    0.0045446122
    0.0038729255
    0.0062493127
    0.0085550864
    0.0064740799
    0.0040216775
    0.012383404
],
:Vzinf => [
    0.68827768
    0.43881487
    1.1673601
    0.59056646
    0.56843374
    0.8124135
    0.68666158
    0.72497447
    0.71232266
    0.72039519
],
)

# Linear up Log down, Dose 100, Dosetime 0.25, tau 9

refdict2 = Dict(
:Cmax => [
190.869
261.177
105.345
208.542
169.334
154.648
153.254
138.327
167.347
125.482
],
:Tmax => [
1
1
1.5
1
4
2.5
2.5
4
3
2
],
:Cdose => [
121.239
62.222
49.849
52.421
0
57.882
19.95
22.724
105.438
13.634
],
:Clast => [
112.846
85.241
67.901
97.625
110.778
69.501
58.051
74.437
93.44
42.191
],
:AUClast => [
9566.59680869131
10054.28647805950
5392.45721941379
9297.09633445033
9519.18087436122
6948.98562111745
6988.77263241364
7058.81896352039
8302.36808633358
5486.83888944199
],
:AUMCtau => [
5477.20423544297
8367.57088170951
3455.34643479800
6014.64604481587
6609.78830163090
5064.72384740413
4976.96365993911
2863.00517022791
5386.88322025614
4713.47970846693
],
:AUCall => [
9566.59680869131
10054.28647805950
5392.45721941379
9297.09633445033
9519.18087436122
6948.98562111745
6988.77263241364
7058.81896352039
8302.36808633358
5486.83888944199
],
:Rsq => [
0.786076957
0.992763591
0.81358898
0.918858685
0.853359952
0.95011904
0.970312315
0.94796904
0.947537895
0.88092269

],
:ARsq => [
0.714769276
0.990351454
0.776306776
0.83771737
0.828919944
0.92517856
0.96041642
0.92195356
0.921306842
0.863911645
],
# LZint
:LZint => [
5.00848559255328
5.42889759540296
4.44064607555325
5.16688496904739
5.14735707974283
4.82967584017057
5.01074587961482
4.96847859724365
4.94725938794774
4.89636108788302
],
:Kel => [
0.00338474394000776
0.01410631494324980
0.00329143037249282
0.00769534422298109
0.00681332791154901
0.00769228066663777
0.01245895597676470
0.00893007980967252
0.00564586491870971
0.01718973683041960
],
:HL => [
204.785706938398
49.137367437811
210.591476080649
90.073577019460
101.734011566509
90.109450057666
55.634451382012
77.619371308325
122.770769499451
40.323315440951
],
:Clast_pred => [
117.3057799
82.53668981
66.93105694
100.7679335
105.2983206
71.93994201
61.17270231
75.60427664
93.76176158
38.81085735
],
:AUCinf => [
42906.1941313004
16097.0411126277
26022.0900281352
21983.3384532182
25778.1957695968
15984.1473646863
11648.1518057779
15394.3547690766
24852.5337997128
7941.2685538530
],
:AUCpct => [
77.7034598328254
37.5395365663056
79.2773861992505
57.7084419901233
63.0727419426760
56.5257660444885
40.0010169085628
54.1467046238288
66.5934743183831
30.9072744205350
],
:MRTtauinf => [
299.791671096989
74.654997085457
305.919973652938
143.538421744963
173.022067431888
124.653434795141
92.735873637166
175.461862330056
178.810514188399
69.516339753006
],
:Cltau => [
0.078847213948573
0.054590500813083
0.132511872732088
0.074823364534525
0.076283206573122
0.089747243392665
0.092646906460213
0.130442680913677
0.081991954283052
0.103060243120434
],
:Vztau => [
23.2948829648816
3.8699335037324
40.2596615257358
9.7231991664617
11.1961742577834
11.6671826317919
7.4361693413954
14.6071125559694
14.5224789228203
5.9954520617241
],
:AUCtau => [
1268.27563070553
1831.82052757491
754.64936037981
1336.48093242129
1310.90451610924
1114.24035123260
1079.36685444479
766.62024499617
1219.63186357018
970.30626915116
],
)

################################################################################
# Linear Trapezoidal with Linear Interpolation, Dose 120, Dosetime 0.0, tau 12

refdict3 = Dict(
:Cmax => [190.869
261.177
105.345
208.542
169.334
154.648
153.254
138.327
167.347
125.482],

:Tmax => [1
1
1.5
1
4
2.5
2.5
4
3
2],

:Cdose => [0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0],

:Clast => [112.846
85.241
67.901
97.625
110.778
69.501
58.051
74.437
93.44
42.191],

:AUClast => [9585.4218
10112.176
5396.5498
9317.8358
9561.26
6966.598
7029.5735
7110.6745
8315.0803
5620.8945],

:AUMCtau => [9984.8168
14630.0690
6024.4953
10299.7210
11466.1230
8467.3568
9003.0193
6457.0058
10095.8180
8367.3005],

:AUCall => [9585.4218
10112.1760
5396.5498
9317.8358
9561.2600
6966.5980
7029.5735
7110.6745
8315.0803
5620.8945],

:Rsq => [0.78607696
0.99276359
0.81358898
0.91885869
0.86367664
0.95011904
0.97031231
0.94796904
0.94753789
0.87969895],

:ARsq => [0.71476928
0.99035145
0.77630678
0.83771737
0.84420187
0.92517856
0.96041642
0.92195356
0.92130684
0.86766884],

# LZint
:LZint => [5.0084856
5.4288976
4.4406461
5.166885
5.1496027
4.8296758
5.0107459
4.9684786
4.9472594
4.8651403],

:Kel => [0.0033847439
0.0141063150
0.0032914304
0.0076953442
0.0068579883
0.0076922807
0.0124589560
0.0089300798
0.0056458649
0.0165437520],

:HL => [204.785706938398
49.1373674378108
210.591476080649
90.0735770194602
101.071502239954
90.109450057666
55.6344513820121
77.6193713083247
122.770769499451
41.8978220179993],

:Clast_pred => [117.30578
82.53669
66.931057
100.76793
105.19623
71.939942
61.172702
75.604277
93.761762
39.408841],

:AUCinf => [42925.019
16154.93
26026.183
22004.078
25714.393
16001.76
11688.953
15446.21
24865.246
8171.1624],

:AUCpct => [77.669383
37.405019
79.26492
57.6540502829908
62.817478
56.463551
39.861391
53.964925
66.559429
31.210589],

:MRTtauinf => [302.40303
75.590599
312.72083
148.34069
172.0933
130.19061
91.908297
161.57402
176.30461
70.260736],

#Cltau, CLss
:Cltau => [0.07185191
0.050414459
0.12240579
0.070132959
0.06902661
0.085106504
0.083532913
0.10859036
0.073251565
0.092756742],

:Vztau => [21.228167
3.5738929
37.18924
9.113687
10.06514
11.063884
6.7046479
12.160066
12.974374
5.6067536],

:AUCtau => [1670.1018
2380.2695
980.34575
1711.0358
1738.46
1409.998
1436.5595
1105.0705
1638.1903
1293.7065],
)

################################################################################
#4 Log trapezoidal ATM, Dose 120, Dosetime 0, tau 12

refdict4 = Dict(
:Cmax => [
190.869
261.177
105.345
208.542
169.334
154.648
153.254
138.327
167.347
125.482
],
:Tmax => [
1
1
1.5
1
4
2.5
2.5
4
3
2
],
:Cdose => [
0
0
0
0
0
0
0
0
0
0
],
:Clast => [
112.846
85.241
67.901
97.625
110.778
69.501
58.051
74.437
93.44
42.191
],
:AUClast => [
9572.8582
10054.0367665966
5391.5322
9296.2179
9518.6531
6948.5757
6987.0645
7064.7816
8298.9634
5485.6538
],
:AUMCtau => [
9973.8062
14631.1197073321
6022.9286
10307.954
11473.081
8471.0956
8982.0378
6271.7444
10040.829690586
8361.7894
],
:AUCall => [
9572.8582
10054.0367665966
5391.5322
9296.2179
9518.6531
6948.5757
6987.0645
7064.7816
8298.9634
5485.6538
],
:Rsq => [
0.78607696
0.99276359
0.81358898
0.91885869
0.85335995
0.95011904
0.97031231
0.94796904
0.94753789
0.88092269
],
:ARsq => [
0.71476928
0.99035145
0.77630678
0.83771737
0.82891994
0.92517856
0.96041642
0.92195356
0.92130684
0.86391165
],
# LZint
:LZint => [
5.0084856
5.4288976
4.4406461
5.166885
5.1473571
4.8296758
5.0107459
4.9684786
4.9472594
4.8963611
],
:Kel => [
0.003384744
0.014106315
0.00329143
0.007695344
0.006813328
0.007692281
0.012458956
0.00893008
0.005645865
0.017189737
],
:HL => [
204.78571
49.137367
210.59148
90.073577
101.73401
90.109450057666
55.634451
77.619371
122.77077
40.323315
],
:Clast_pred => [
117.30578
82.53669
66.931057
100.76793
105.29832
71.939942
61.172702
75.604277
93.761762
38.810857
],
:AUCinf => [
42912.456
16096.791
26021.165
21982.4599914207
25777.668
15983.737
11646.444
15400.317
24849.129
7940.0834
],
:AUCpct => [
77.692122
37.540119
79.280204
57.710748
63.074033
56.527216
40.006884
54.12574
66.602599
30.911888
],
:MRTtauinf => [
302.63508
75.323724
313.06798
148.31081
172.5577
130.22554
91.866692
164.91799
176.98523
68.167555
],
:Cltau => [
0.071927102
0.050429351
0.12256044
0.070184147
0.069035447
0.0852177496596485
0.08379761
0.11110872
0.073575577
0.092819834
],
:Vztau => [
21.250382
3.5749486
37.236223
9.1203389
10.132412
11.078346
6.7258934
12.442074
13.031764
5.399724
],
:AUCtau => [
1668.3558
2379.5666
979.10878
1709.7878
1738.2375
1408.1573
1432.0218
1080.0233
1630.976
1292.8271
],
)

################################################################################

urefdict = Dict{Symbol, Float64}(
#:N_Samples => 5,
#:Dose => 100,
:Rsq => 0.90549162,
:ARsq => 0.81098324, #Rsq_adjusted
#:Corr_XY => -0.95157323,
#:No_points_lambda_z => 3,
:Kel => 0.13445441, #Lambda_z
#:Lambda_z_intercept => 0.79280975,
#:Lambda_z_lower => 4,
#:Lambda_z_upper => 15,
:HL => 5.1552579, #HL_Lambda_z
#:Span => 2.1337439,
#:Tlag => 0,
:Tmax => 1.5, #Tmax_Rate
:Maxrate => 4, #Max_Rate
#:Mid_Pt_last => 15,
:Rlast => 0.33333333, #Rate_last
#:Rate_last_pred => 0.2940497,
:AUClast => 17.125, #AURC_last
#:AURC_last_D => 0.17125,
:Vol => 11, #Vol_UR
:AR => 16, #Amount_Recovered
:Prec => 16, #Percent_Recovered
:AUCall => 17.125, #AURC_all
:AUCinf => 19.604155, #AURC_INF_obs
:AUCpct => 12.646069, #AURC_%Extrap_obs
#:AURC_INF_pred => 19.311984,
#:AURC_%Extrap_pred => 11.324493,
)


text:
struct BlockIterator
    rows::Cint
    cols::Cint
    ni::Cint
    nj::Cint
    n::Cint
    xbsize::Cint
    ybsize::Cint
end

function Base.iterate(obj::BlockIterator, iter::Int = 0)
    iter == obj.n && return nothing
    j = floor(Int, iter / obj.ni)
    i = iter % obj.ni
    nrows = if (i + 1) * obj.ybsize < obj.rows
        obj.ybsize
    else
        obj.rows - i * obj.ybsize
    end
    ncols = if (j + 1) * obj.xbsize < obj.cols
        obj.xbsize
    else
        obj.cols - j * obj.xbsize
    end
    (((i, j), (nrows, ncols)), iter + 1)
end

function blocks(fname::AbstractString, dset_name::AbstractString)
    (chunkrows, chunkcols, _) = get_chunk(fname, dset_name)
    rows, cols, _ = size(fname, dset_name)
    ni = ceil(Cint, rows / chunkrows)
    nj = ceil(Cint, cols / chunkcols)
    BlockIterator(rows, cols, ni, nj, ni * nj, chunkcols, chunkrows)
end

Base.length(b::BlockIterator) = b.n

struct WindowIterator
    blockiter::BlockIterator
end

function Base.iterate(obj::WindowIterator, iter::Int = 0)
    handle = obj.blockiter
    next = Base.iterate(handle, iter)
    next == nothing && return nothing
    (((i, j), (nrows, ncols)), iter) = next
    # (((1:ncols) .+ j * handle.xbsize, (1:nrows) .+ i * handle.ybsize), iter)
    (((1:nrows) .+ i * handle.ybsize, (1:ncols) .+ j * handle.xbsize), iter)
end

windows(fname::AbstractString, dset_name::AbstractString) = WindowIterator(blocks(fname, dset_name))

Base.length(w::WindowIterator) = w.blockiter.n
Base.size(w::WindowIterator) = (length(w),)


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using DocStringExtensions;

"""
    NURBStoolbox

NURBS Toolbox module ported from the Matlab NURBS toolbox of [2000Spink](@cite)
(can currently be found on Mathworks FileExchange[^1]). The Matlab package dates
approximately back to the year 2000. Newer versions are written for Octave and
are hosted as Octave Community package[^2] (those might be compatible with Matlab).

This package is licensed under GNU General Puplic License version 2.

[^1]: Penguian (2010). NURBS Toolbox by D.M. Spink [https://www.mathworks.com/matlabcentral/fileexchange/26390-nurbs-toolbox-by-d-m-spink](https://www.mathworks.com/matlabcentral/fileexchange/26390-nurbs-toolbox-by-d-m-spink), MATLAB Central File Exchange. Retrieved April 14, 2021.
[^2]: M. Spink, D. Claxton, C. de Falco, R. Vazquez (2021-03-09) Nurbs. Octave Forge Community packages. [https://octave.sourceforge.io/nurbs/index.html](https://octave.sourceforge.io/nurbs/index.html)

# Exports
$(DocStringExtensions.EXPORTS)

# Imports
$(DocStringExtensions.IMPORTS)
"""
module NURBStoolbox
# deg2rad: is a julia build in function therefore omitted
# rad2deg: is a julia build in function therefore omitted
# igesout: is a not necessary third party functions with uncertain licence
#          and was therefore not ported
# demos  : seemed to be a matlab interface to create demo views -> not ported

# Evaluation functions/methods
export basisfun;
include("basisfun.jl");
export bspeval;
include("bspeval.jl");
export NURBS, NURBS1D, NURBS2D, nrbmak;
include("nrbmak.jl");
export nrbeval;
include("nrbeval.jl");
export nrbplot, nrbplot!;
include("nrbplot.jl");

# Utility and modification functions/methods
export findspan;
include("findspan.jl")
export vecscale;
include("vecscale.jl")
export vectrans;
include("vectrans.jl")
export vecrotx;
include("vecrotx.jl")
export vecroty;
include("vecroty.jl")
export vecrotz;
include("vecrotz.jl")
export nrbtform;
include("nrbtform.jl");
export nrbextrude;
include("nrbextrude.jl");
export nrbrevolve;
include("nrbrevolve.jl");
export vecdot;
include("vecdot.jl");
export veccross;
include("veccross.jl");
export vecmag;
include("vecmag.jl");
export vecmag2;
include("vecmag2.jl");
export vecnorm_toolbox
include("vecnorm_toolbox.jl");
export vecangle;
include("vecangle.jl");
export nrbtransp;
include("nrbtransp.jl");
export bspdegelev;
include("bspdegelev.jl")
export nrbdegelev;
include("nrbdegelev.jl");
export bspkntins;
include("bspkntins.jl");
export nrbkntins;
include("nrbkntins.jl");
export bspderiv;
include("bspderiv.jl");
export nrbderiv;
include("nrbderiv.jl");
export nrbdeval;
include("nrbdeval.jl");
export nrbreverse;
include("nrbreverse.jl");
export nrbruled;
include("nrbruled.jl");
export nrbcoons;
include("nrbcoons.jl");

# "recipes" of NURBS structures
export nrbline;
include("nrbline.jl");
export nrbrect;
include("nrbrect.jl");
export nrb4surf;
include("nrb4surf.jl");
export nrbcirc;
include("nrbcirc.jl");
export nrbcylind;
include("nrbcylind.jl");

# demo/test functions/methods
export nrbtestcrv;
include("nrbtestcrv.jl");
export nrbtestsrf;
include("nrbtestsrf.jl");
export demoline;
include("demoline.jl");
export demorect;
include("demorect.jl");
export democirc;
include("democirc.jl");
export demoellip;
include("demoellip.jl");
export democurve;
include("democurve.jl");
export demohelix;
include("demohelix.jl");
export democylind;
include("democylind.jl");
export demotorus;
include("demotorus.jl");
export demorevolve;
include("demorevolve.jl");
export demodegelev;
include("demodegelev.jl")
export demokntins;
include("demokntins.jl")
export demodercrv;
include("demodercrv.jl");
export demodersrf;
include("demodersrf.jl");
export demoruled;
include("demoruled.jl");
export demo4surf;
include("demo4surf.jl");
export democoons;
include("democoons.jl");
export demogeom;
include("demogeom.jl");

end # NURBStoolbox


text:
# Deprecated functionality

import Base: @deprecate_binding

@deprecate_binding BLAS CUBLAS
@deprecate_binding FFT CUFFT


text:
using Plots, Images; pyplot()

img = load("../data/stars.png")
gImg = red.(img)*0.299 + green.(img)*0.587 + blue.(img)*0.114
rows, cols = size(img)

println("Highest intensity pixel: ", findmax(gImg))

function boxBlur(image,x,y,d)
    if x<=d || y<=d || x>=cols-d || y>=rows-d
        return image[x,y]
    else
        total = 0.0
        for xi = x-d:x+d
            for yi = y-d:y+d
                total += image[xi,yi]
            end
        end
        return total/((2d+1)^2)
    end
end

blurImg = [boxBlur(gImg,x,y,5) for x in 1:cols, y in 1:rows]

yOriginal, xOriginal = argmax(gImg).I
yBoxBlur, xBoxBlur   = argmax(blurImg).I

p1 = heatmap(gImg, c=:Greys, yflip=true)
p1 = scatter!((xOriginal, yOriginal), ms=60, ma=0, msw=4, msc=:red)
p2 = heatmap(blurImg, c=:Greys, yflip=true)
p2 = scatter!((xBoxBlur, yBoxBlur), ms=60, ma=0, msw=4, msc=:red)

plot(p1, p2, size=(800, 400), ratio=:equal, xlims=(0,cols), ylims=(0,rows),
	colorbar_entry=false, border=:none, legend=:none)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

AG.registerdrivers() do
    AG.read("data/point.geojson") do dataset
        layer = AG.getlayer(dataset, 0)
        AG.getfeature(layer, 0) do f1
            AG.getfeature(layer, 2) do f2
                @test sprint(print, f1) == """
                Feature
                  (index 0) geom => POINT
                  (index 0) FID => 2.0
                  (index 1) pointname => point-a
                """
                @test sprint(print, AG.getgeom(f1)) == "Geometry: POINT (100 0)"
                fid1 = AG.getfid(f1)
                @test fid1 == 0

                @test sprint(print, f2) == """
                Feature
                  (index 0) geom => POINT
                  (index 0) FID => 0.0
                  (index 1) pointname => a
                """
                @test sprint(print, AG.getgeom(f2)) == "Geometry: POINT (100 0)"
                fid2 = AG.getfid(f2)
                @test fid2 == 2

                AG.clone(f1) do f3
                    @test AG.equals(AG.getgeom(f1), AG.getgeom(f3)) == true
                end
                AG.setfid!(f1, fid2); AG.setfid!(f2, fid1)
                @test AG.getfid(f1) == 2
                @test AG.getfid(f2) == 0

                @test AG.getgeomfieldindex(f1, "geom") == -1
                @test AG.getgeomfieldindex(f1, "") == 0
                @test AG.getgeomfieldindex(f2, "geom") == -1
                @test AG.getgeomfieldindex(f2, "") == 0
                @test AG.gettype(AG.getgeomfielddefn(f1, 0)) == GDAL.wkbPoint
                @test AG.gettype(AG.getgeomfielddefn(f2, 0)) == GDAL.wkbPoint
            end
        end

        AG.getfeature(layer, 0) do f
            @test AG.toWKT(AG.getgeomfield(f,0)) == "POINT (100 0)"
            AG.setgeomfielddirectly!(f, 0, AG.unsafe_createpoint(0,100))
            @test AG.toWKT(AG.getgeomfield(f,0)) == "POINT (0 100)"
            AG.createpolygon([(0.,100.),(100.,0.)]) do poly
                AG.setgeomfield!(f, 0, poly)
            end
            @test AG.toWKT(AG.getgeomfield(f,0)) == "POLYGON ((0 100,100 0))"

            AG.setstylestring!(f, "@Name")
            @test AG.getstylestring(f) == "@Name"
            AG.setstylestring!(f, "NewName")
            @test AG.getstylestring(f) == "NewName"

            AG.setstyletabledirectly!(f, AG.unsafe_createstyletable())
            AG.createstyletable() do st
                AG.setstyletable!(f, st)
            end

            AG.setnativedata!(f, "nativedata1")
            @test AG.getnativedata(f) == "nativedata1"
            AG.setnativedata!(f, "nativedata2")
            @test AG.getnativedata(f) == "nativedata2"

            AG.setmediatype!(f, "mediatype1")
            @test AG.getmediatype(f) == "mediatype1"
            AG.setmediatype!(f, "mediatype2")
            @test AG.getmediatype(f) == "mediatype2"

            @test AG.validate(f, GDAL.OGR_F_VAL_NULL, false) == true
            @test AG.validate(f, GDAL.OGR_F_VAL_GEOM_TYPE, false) == false
            @test AG.validate(f, GDAL.OGR_F_VAL_WIDTH, false) == true
            @test AG.validate(f, GDAL.OGR_F_VAL_ALLOW_NULL_WHEN_DEFAULT, false) == true
            @test AG.validate(f, GDAL.OGR_F_VAL_ALLOW_DIFFERENT_GEOM_DIM, false) == true

            @test AG.getfield(f, 1) == "point-a"
            AG.setdefault!(AG.getfielddefn(f, 1),"nope")
            @test AG.getfield(f, 1) == "point-a"
            AG.unsetfield!(f, 1)
            @test AG.getfield(f, 1) == nothing
            AG.fillunsetwithdefault!(f, notnull=false)
            @test AG.getfield(f, 1) == "nope"
        end
    end

    @testset "In-Memory Driver" begin
        AG.create("", "MEMORY") do output
            layer = AG.createlayer(output, "dummy", geom=GDAL.wkbPolygon)
            AG.createfielddefn("int64field", GDAL.OFTInteger64) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("doublefield", GDAL.OFTReal) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("intlistfield", GDAL.OFTIntegerList) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("int64listfield", GDAL.OFTInteger64List) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("doublelistfield", GDAL.OFTRealList) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("stringlistfield", GDAL.OFTStringList) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("binaryfield", GDAL.OFTBinary) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfielddefn("datetimefield", GDAL.OFTDateTime) do fielddefn
                AG.createfield!(layer, fielddefn)
            end
            AG.createfeature(layer) do feature
                AG.setfield!(feature, 0, 1)
                AG.setfield!(feature, 1, 1.0)
                AG.setfield!(feature, 2, Cint[1, 2])
                AG.setfield!(feature, 3, Int64[1, 2])
                AG.setfield!(feature, 4, Float64[1.0, 2.0])
                AG.setfield!(feature, 5, ["1", "2.0"])
                AG.setfield!(feature, 6, GDAL.GByte[1,2,3,4])
                AG.setfield!(feature, 7, Dates.DateTime(2016,9,25,21,17,0))

                AG.createfeature(layer) do newfeature
                    AG.setfrom!(newfeature, feature)
                    @test AG.getfield(newfeature, 0) == 1
                    @test AG.getfield(newfeature, 1) ≈ 1.0
                    @test AG.getfield(newfeature, 2) == Int32[1, 2]
                    @test AG.getfield(newfeature, 3) == Int64[1, 2]
                    @test AG.getfield(newfeature, 4) ≈ [1.0, 2.0]
                    @test AG.getfield(newfeature, 5) == ["1", "2.0"]
                    @test AG.getfield(newfeature, 6) == GDAL.GByte[1,2,3,4]
                    @test AG.getfield(newfeature, 7) == Dates.DateTime(2016,9,25,21,17,0)

                    AG.createfeature(layer) do lastfeature
                        AG.setfrom!(lastfeature, feature)
                        AG.setfield!(lastfeature, 0, 45)
                        AG.setfield!(lastfeature, 1, 18.2)
                        AG.setfield!(lastfeature, 5, ["foo", "bar"])
                        @test AG.getfield(lastfeature, 0) == 45
                        @test AG.getfield(lastfeature, 1) ≈ 18.2
                        @test AG.getfield(lastfeature, 2) == Int32[1, 2]
                        @test AG.getfield(lastfeature, 3) == Int64[1, 2]
                        @test AG.getfield(lastfeature, 4) ≈ [1.0, 2.0]
                        @test AG.getfield(lastfeature, 5) == ["foo", "bar"]
                        @test AG.getfield(lastfeature, 6) == GDAL.GByte[1,2,3,4]
                        @test AG.getfield(lastfeature, 7) == Dates.DateTime(2016,9,25,21,17,0)

                        @test AG.getfield(newfeature, 0) == 1
                        @test AG.getfield(newfeature, 1) ≈ 1.0
                        @test AG.getfield(newfeature, 5) == ["1", "2.0"]
                        AG.setfrom!(newfeature, lastfeature, collect(Cint, 0:7))
                        @test AG.getfield(newfeature, 0) == 45
                        @test AG.getfield(newfeature, 1) ≈ 18.2
                        @test AG.getfield(newfeature, 5) == ["foo", "bar"]
                    end
                end
            end
        end
    end
end


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

name = "spglib"
version = v"1.16.5"

# Collection of sources required to build spglib
sources = [
    ArchiveSource("https://github.com/atztogo/spglib/archive/v$(version).tar.gz",
                  "1bbde03b6b78da756c07f458bd90d84f3c253841b9b0632db5b72c5961e87aef"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/spglib-*/
if [[ ${target} == *-mingw32 ]]; then
    sed -i -e 's/LIBRARY/RUNTIME/' CMakeLists.txt
fi
mkdir _build
cd _build/
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      ..
make -j${nproc}
make install VERBOSE=1
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    # For OpenMP we use libomp from `LLVMOpenMP_jll` where we use LLVM as compiler (BSD
    # systems), and libgomp from `CompilerSupportLibraries_jll` everywhere else.
    Dependency("CompilerSupportLibraries_jll"; platforms=filter(!Sys.isbsd, platforms)),
    Dependency("LLVMOpenMP_jll"; platforms=filter(Sys.isbsd, platforms)),
]

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