text:
__precompile__(false)

module Data

using FileIO
using Images
using Juno
using ColorTypes
using FixedPointNumbers
using DataFrames
using CSVFiles
using DSB2018.Errors

using ImageView

global const STAGE1_TEST  = joinpath(@__DIR__, "../etc/stage1_test")
global const STAGE1_TRAIN = joinpath(@__DIR__, "../etc/stage1_train")

export Image, getTrainImageIds
export loadTrainImage
export getTrainImageMasks
export Prediction
export submission, submit!

function loadImage(id::String, stage_path::String)
    image_file = joinpath(stage_path, id, "images", id * ".png")
    return Images.load(image_file)
end

function loadMasks(id::String, stage_path::String)::Dict{String, Matrix{Bool}}
    out = Dict{String, Matrix{Bool}}()
    masks_root_path = joinpath(stage_path, id, "masks")
    for id in readdir(masks_root_path)
        mask_file = joinpath(masks_root_path, id)
        new_mask = Images.load(mask_file)
        out[id] = new_mask .> 0.0
        if isempty(find(out[id]))
            @show mask_file
            throw("Mask is empty")
        end
    end

    return out
end

struct Image{X}
    id::String
    image::Matrix{X}
    mask::Matrix{Int64}
    mask_labels::Vector{Int64}
end

struct Prediction
    id::String
    score::Float64
    mask::Matrix{Int64}
end

function loadTrainImage(id::String)::Image
    global STAGE1_TRAIN
    img = loadImage(id, STAGE1_TRAIN)
    masks = loadMasks(id, STAGE1_TRAIN)
    labeled_mask, labels = label_masks(collect(values(masks)))
    return Image(id, img, labeled_mask, labels)
end

function loadTestImage(id::String)::Image
    global STAGE1_TEST
    img = loadImage(id, STAGE1_TEST)
    return Image(id, img, zeros(Int64, size(img)), Int64[])
end

function getTrainImageIds()
    global STAGE1_TRAIN
    readdir(STAGE1_TRAIN)
end

function getTestImageIds()
    global STAGE1_TEST
    readdir(STAGE1_TEST)
end

function loadTrainImages(n::Int64=9999999999)
    ids = getTrainImageIds()
    n = min(n, length(ids))
    ids = ids[1:n]
    out = Vector{Image}(length(ids))
    Juno.@progress "Loading images" for k in eachindex(ids)
        out[k] = loadTrainImage(ids[k])
    end
    return out
end

function loadTestImages(n::Int64=999999999)
    ids = getTestImageIds()
    n = min(n, length(ids))
    out = Vector{Image}(length(ids))
    Juno.@progress "Loading images" for k in eachindex(ids)
        out[k] = loadTestImage(ids[k])
    end
    return out
end

function getTrainImageMasks(this::Image{X})::Vector{Matrix{Bool}} where {X <: Any}
    out = Vector{Matrix{Bool}}(length(this.mask_labels))
    for k in eachindex(out)
        out[k] = this.mask .== this.mask_labels[k]
    end

    return out
end

function label_masks(this::Vector{Matrix{Bool}})
    szs = size.(this)
    @assert length(unique(szs)) == 1
    sz = first(szs)
    out = zeros(Int64, sz)

    this_float = [x.*1.0 for x in this]
    num_p      = [sum(x) for x in this]

    s = sum(this_float)
    @assert maximum(s) == 1
    @assert minimum(num_p) > 0

    labels = Int64.(collect(1 : length(this)))
    for j in eachindex(this)
        out[this[j]] .= labels[j]
    end

    return out, labels
end

function encoded_string(this::Vector{Pair{Int64, Int64}})
    return join([join([x.first, x.second], " ") for x in this], " ")
end

function submission(this::Vector{Prediction})::DataFrame
    image_ids = Vector{String}(0)
    encodedpixels = Vector{String}(0)
    for prediction in this
        rle_encoded = Errors.encode(prediction.mask)
        for pairs::Vector{Pair{Int64, Int64}} in values(rle_encoded)
            push!(image_ids, prediction.id)
            push!(encodedpixels, encoded_string(pairs))
        end
    end

    return DataFrame(ImageId = image_ids, EncodedPixels = encodedpixels)
end

function submit(this::DataFrame)
    submit_path = joinpath(@__DIR__, "../submissions")
    submit_file = string(now(Dates.UTC), ".csv")
    save(joinpath(submit_path, submit_file), this, quotechar = nothing)
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NextHopParameters <: SwaggerModel
    targetResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetResourceId
    sourceIPAddress::Any # spec type: Union{ Nothing, String } # spec name: sourceIPAddress
    destinationIPAddress::Any # spec type: Union{ Nothing, String } # spec name: destinationIPAddress
    targetNicResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetNicResourceId

    function NextHopParameters(;targetResourceId=nothing, sourceIPAddress=nothing, destinationIPAddress=nothing, targetNicResourceId=nothing)
        o = new()
        validate_property(NextHopParameters, Symbol("targetResourceId"), targetResourceId)
        setfield!(o, Symbol("targetResourceId"), targetResourceId)
        validate_property(NextHopParameters, Symbol("sourceIPAddress"), sourceIPAddress)
        setfield!(o, Symbol("sourceIPAddress"), sourceIPAddress)
        validate_property(NextHopParameters, Symbol("destinationIPAddress"), destinationIPAddress)
        setfield!(o, Symbol("destinationIPAddress"), destinationIPAddress)
        validate_property(NextHopParameters, Symbol("targetNicResourceId"), targetNicResourceId)
        setfield!(o, Symbol("targetNicResourceId"), targetNicResourceId)
        o
    end
end # type NextHopParameters

const _property_map_NextHopParameters = Dict{Symbol,Symbol}(Symbol("targetResourceId")=>Symbol("targetResourceId"), Symbol("sourceIPAddress")=>Symbol("sourceIPAddress"), Symbol("destinationIPAddress")=>Symbol("destinationIPAddress"), Symbol("targetNicResourceId")=>Symbol("targetNicResourceId"))
const _property_types_NextHopParameters = Dict{Symbol,String}(Symbol("targetResourceId")=>"String", Symbol("sourceIPAddress")=>"String", Symbol("destinationIPAddress")=>"String", Symbol("targetNicResourceId")=>"String")
Base.propertynames(::Type{ NextHopParameters }) = collect(keys(_property_map_NextHopParameters))
Swagger.property_type(::Type{ NextHopParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NextHopParameters[name]))}
Swagger.field_name(::Type{ NextHopParameters }, property_name::Symbol) =  _property_map_NextHopParameters[property_name]

function check_required(o::NextHopParameters)
    (getproperty(o, Symbol("targetResourceId")) === nothing) && (return false)
    (getproperty(o, Symbol("sourceIPAddress")) === nothing) && (return false)
    (getproperty(o, Symbol("destinationIPAddress")) === nothing) && (return false)
    true
end

function validate_property(::Type{ NextHopParameters }, name::Symbol, val)
end

ext: jl


text:
using Luxor
@svg begin
    setopacity(0.85)
    steps = 10
    gap   = 12
    for (n, θ) in enumerate(range(π/4, step=(2π-π/2)/steps, length=steps))
        sethue([Luxor.julia_green,
            Luxor.julia_red,
            Luxor.julia_purple,
            Luxor.julia_blue][mod1(n, 4)])
        sector(Point(-150, 0), 50, 100, θ, θ + 2π/steps - deg2rad(gap), :fill)
    end
     fontsize(90)
     fontface("MyanmarMN-Bold")
     sethue("black")
     text("oefplots.jl", Point(-178, -10) ,halign=:left,valign = :middle)
end 700 200 "logo.svg"


text:
@testset "set spins" begin
    # Use a list of integers.
    @test set_spins(3, 2, σᶻ) == [σ⁰, σᶻ, σ⁰]
    @test_throws MethodError set_spins(3, 2, [σᶻ])
    @test set_spins(3, [2], [σᶻ]) == [σ⁰, σᶻ, σ⁰]
    @test set_spins(4, [2], [σᶻ]) == [σ⁰, σᶻ, σ⁰, σ⁰]
    @test set_spins(4, [2, 3], [σᶻ, σˣ]) == [σ⁰, σᶻ, σˣ, σ⁰]

    # Use a CartesianIndex.
    @test set_spins(4, CartesianIndex(2, 3), [σᶻ, σˣ]) == [σ⁰, σᶻ, σˣ, σ⁰]
    @test_throws AssertionError set_spins(4, CartesianIndex(2), [σᶻ, σˣ])
    @test_throws AssertionError set_spins(4, CartesianIndex(2,3), [σᶻ])
end


text:
# Note that we tried and failed to automate the detection of "train" mode looking at the type
# of argument.  The argument is of type Value only during training and only if it is a value
# influenced by model weights.  However people typically apply dropout to the input which is
# not a Value.  So we are going back to no automation, make sure to supply p=0 during testing
# to stop dropout.

# In 1.1.3 trying to automate again, this time using the @diff context: By default drop if
# AutoGrad is recording.

"""
    dropout(x, p; drop, seed)

Given an array `x` and probability `0<=p<=1` return an array `y` in which each element is 0
with probability `p` or `x[i]/(1-p)` with probability `1-p`. Just return `x` if `p==0`, or
`drop=false`. By default `drop=true` in a `@diff` context, `drop=false` otherwise.  Specify a
non-zero `seed::Number` to set the random number seed for reproducible results. See
[(Srivastava et al. 2014)](http://www.jmlr.org/papers/v15/srivastava14a.html) for a reference.

"""
function dropout(x,p; seed=0, drop=training())
    if !drop
        x
    elseif 0 < p < 1
        if seed != 0; Knet.seed!(seed); end
        dropout!(p,x,similar(x))
    elseif p == 0
        x
    elseif p == 1
        zero(x)
    else
        error("Dropout probability not in [0:1]: $p")
    end
end

# Use a more meaningful name than `recording` for Knet.
"`training()` returns `true` only inside a `@diff` context, e.g. during a training iteration of a model."
training() = AutoGrad.recording()

function dropback(dy,y,x,p)
    if 0 < p < 1
        dropback!(p,x,y,dy,similar(x))
    elseif p == 0
        dy
    elseif p == 1
        zero(x)
    else
        error("Dropout probability not in [0:1]: $p")
    end
end

# Turn dropout into an AutoGrad primitive
@primitive dropout(x,p;seed=0,drop=training()),dy,y dropback(value.((dy,y,x,p))...)

# GPU implementation
for S in (32,64)
    T = Symbol("Float$S")
    forw = Symbol("dropout_$S")
    back = Symbol("dropback_$S")
    @eval begin
        function dropout!(p::Number, x::KnetArray{$T}, y::KnetArray{$T})
            rand!(y)
            @knet8($forw,(Cint,$T,Ptr{$T},Ptr{$T}),length(y),$T(p),x,y)
            return y
        end
        function dropback!(p::Number, x::KnetArray{$T}, y::KnetArray{$T}, dy::KnetArray{$T}, dx::KnetArray{$T})
            @knet8($back,(Cint,$T,Ptr{$T},Ptr{$T},Ptr{$T},Ptr{$T}),length(dx),$T(p),x,y,dy,dx)
            return dx
        end
    end
end

# CPU implementation
function dropout!(p,x,y)
    rand!(y)
    p = convert(eltype(y),p)
    q = 1-p
    @inbounds for i=1:length(y)
        if y[i] > p
            y[i] = x[i] / q
        else
            y[i] = 0
        end
    end
    return y
end

function dropback!(p,x,y,dy,dx)
    p = convert(eltype(y),p)
    q = 1-p
    @inbounds for i=1:length(dx)
        if y[i] == 0
            dx[i] = 0
        else
            dx[i] = dy[i] / q
        end
    end
    return dx
end


text:
using GLVisualize, GeometryTypes
using GLAbstraction
if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(linspace(0f0,1f0,360))
end

function sin_torus(radius, thickness, N, t)
    pirange = linspace(0, 2pi, N^2)
    points = zeros(Point3f0, N^2)
    p,q = 6,9
    for (i,x)=enumerate(pirange)
        r = cos(q*(x+t))+thickness
        points[i] = Point3f0(r*cos(p*x),r*sin(p*x), -sin(q*x))/radius
    end
    points
end

t = const_lift(*, 1000f0, timesignal)
ps1 = const_lift(sin_torus, 2, 6, 10, t)
ps2 = const_lift(sin_torus, 5, 4, 10, t)
ps3 = const_lift(sin_torus, 10, 2, 10, t)
ps = map(vcat, ps1, ps2, ps3)
_view(visualize((Circle, ps), scale=Vec2f0(0.04), billboard=true), camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
# Model
include_model("car")

# Horizon
T = 11

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds

# ul <= u <= uu

# Initial and final states
x1 = [0.0; 0.0; 0.0]
xT = [0.0; 1.0; 0.0]
xl, xu = state_bounds(model, T, x1 = x1)#, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [(t < T ? Diagonal(0.0 * ones(model.n))
            : Diagonal(1000.0 * ones(model.n))) for t = 1:T],
        [Diagonal(0.0 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T], [zeros(model.m) for t = 1:T])

# Problem
prob = trajectory_optimization_problem(
           model,
           obj,
           T,
           h = h,
           xl = xl,
           xu = xu)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [0.001 * ones(model.m) for t = 1:T-1]

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve nominal problem
@time z̄, info = solve(prob, copy(z0))
x̄, ū = unpack(z̄, prob)
@show x̄[end]

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x̄, Δt = h, traj = true)

# learn a policy

function policy(model, x, θ)
    n = model.n
    m = model.m

    W1 = reshape(θ[1:(n * n)], n, n)
    b1 = θ[n * n .+ (1:n)]

    W2 = reshape(θ[n * n + n .+ (1:(m * n))], m, n)
    b2 = θ[n * n + n + m * n .+ (1:m)]

    W3 = reshape(θ[n * n + n + m * n + m .+ (1:(m * m))], m, m)

    z1 = tanh.(W1 * x + b1)
    z2 = W3 * tanh.(W2 * z1 + b2)

    # W1 = reshape(θ[1:(m * n)], m, n)
    # b1 = θ[m * n .+ (1:m)]
    #
    # z2 = W1 * x + b1

    return z2
end

function objective(model, x, u, T)
    J = 0.0

    for t = 1:T
        if t < T
            J += 0.01 * (x[:, t] - xT)' * (x[:, t] - xT)
        else
            J += (x[:, T] - xT)' *Diagonal([1.0, 10.0, 1.0]) * (x[:, T] - xT)
        end
    end

    return J
end

function rollout(model, x1, θ, T, h)
    n = model.n
    m = model.m

    x = zeros(eltype(θ), n, T)
    x[:, 1] = x1

    u = zeros(eltype(θ), m, T-1)
    u[:, 1] = policy(model, x[:, 1], θ)

    for t = 1:T-1
        u[:, t] = policy(model, x[:, t], θ)
        x[:, t+1] = fd(model, x[:, t], u[:, t], 0.01 * randn(model.d), h, t)
    end

    J = objective(model, x, u, T)

    return x, u, J
end

function solve(model, x1)
    n = model.n
    m = model.m

    # policy parameters
    p = n * n + n + m * n + m + m * m
    θ = 0.1 * randn(p)

    # rewards
    r_hist = []

    # samples
    N = length(x1)

    α = 1.0e-3
    for i = 1:1e5
        function reward(z)
            J = sum([rollout(model, _x1, z, T, h)[3] for _x1 in x1])
            return J / N
        end

        r = reward(θ)
        push!(r_hist, r)

        if i % 1000 == 0
            println("reward (iter = $i): $(r_hist[end])")
        end
        ∇r = ForwardDiff.gradient(reward, θ)

        # println("gradient $(∇r)")

        θ .-= α * ∇r
    end

    return θ
end

s1 = resample(x1, Diagonal(0.5 * ones(n)), 1.0)
θ = solve(model, [x1, s1...])
s1[4]
_x̄, ū, J = rollout(model, x1, θ, T, h)
x̄ = [_x̄[:, t] for t = 1:T]
visualize!(vis, model, x̄, Δt = h, traj = true)


text:
@testset "sampling.jl" begin
    if GRAPH_T == :coo
        @testset "sample_neighbors" begin
            # replace = false
            dir = :in
            nodes = 2:3
            g = rand_graph(10, 40, bidirected=false, graph_type=GRAPH_T)
            sg = sample_neighbors(g, nodes; dir)
            @test sg.num_nodes == 10
            @test sg.num_edges == sum(degree(g, i; dir) for i in nodes)
            @test size(sg.edata.EID) == (sg.num_edges,)
            @test length(union(sg.edata.EID)) == length(sg.edata.EID)
            adjlist = adjacency_list(g; dir)
            s, t = edge_index(sg)
            @test all(t .∈ Ref(nodes))
            for i in nodes
                @test sort(neighbors(sg, i; dir)) == sort(neighbors(g, i; dir))
            end

            # replace = true
            dir = :out
            nodes = 2:3
            K = 2
            g = rand_graph(10, 40, bidirected=false, graph_type=GRAPH_T)
            sg = sample_neighbors(g, nodes, K; dir, replace=true)
            @test sg.num_nodes == 10
            @test sg.num_edges == sum(K for i in nodes)
            @test size(sg.edata.EID) == (sg.num_edges,)
            adjlist = adjacency_list(g; dir)
            s, t = edge_index(sg)
            @test all(s .∈ Ref(nodes))
            for i in nodes
                @test issubset(neighbors(sg, i; dir), adjlist[i])
            end

            # dropnodes = true
            dir = :in
            nodes = 2:3
            g = rand_graph(10, 40, bidirected=false, graph_type=GRAPH_T)
            g = GNNGraph(g, ndata=(x1=rand(10),), edata=(e1=rand(40),))
            sg = sample_neighbors(g, nodes; dir, dropnodes=true)
            @test sg.num_edges == sum(degree(g, i; dir) for i in nodes)
            @test size(sg.edata.EID) == (sg.num_edges,)
            @test size(sg.ndata.NID) == (sg.num_nodes,)
            @test sg.edata.e1 == g.edata.e1[sg.edata.EID]
            @test sg.ndata.x1 == g.ndata.x1[sg.ndata.NID]
            @test length(union(sg.ndata.NID)) == length(sg.ndata.NID)
        end
    end
end


text:
# Because julia-vim and YouCompleteMe don't get along
_getsigma(x) = x.σ

_dims(c::Flux.Conv{N}) where {N} = N
_dims(c::Flux.CrossCor{N}) where {N} = N

function _conv_impl(c, x::Node)
    N = _dims(c)

    # We flip standard arrays since nGraph really perform cross-correlation
    n = Node(c.weight)
    cn = NNlib.conv(
        x,
        n;
        stride = reverse(c.stride),
        pad = reverse(c.pad),
        dilation = reverse(c.dilation)
    )

    # Broadcast the bias along the first `N` dimensions and the last
    axis_set = [collect(1:N); N+2]
    bb = broadcast(Node(c.bias), size(cn); axes = axis_set)

    node = _getsigma(c).(cn .+ bb)
    return node
end

# Again, getting around a Cassette issue
_dense_impl(d::Flux.Dense, x::Node) = (d.σ).(d.W * x .+ d.b)

# TODO: ngraph makes the dictinction between training and inference. For now, we will
# assume training, but eventually I can add a parameter to SnoopCtx that will determine
# if we're training or inferring and pass that information here.
function _batchnorm_impl(BN::Flux.BatchNorm, x::Node)
    # Create the batchnorm op and then do activation.
    γ = Node(BN.γ)
    β = Node(BN.β)

    n = batchnorm_training(x, γ, β, BN.ϵ)

    # The batchnorm_training op in ngraph returns a tuple
    # (normalized, gamma, beta). We apply the activation function to the normalized output.
    #
    # We also have call `get_output_element` on the other two outputs so that downstream
    # graph rewriting in ngraph works correctly.
    #
    # Also note that we have to `__register` the nodes so they become hidden outputs of the
   # compiled ngraph graph. Otherwide, things break horribly
    a = get_output_element(n, 1)
    __register(get_output_element(n, 2))
    __register(get_output_element(n, 3))

    return BN.λ.(a)
end

# Need to flip the convolution kernels
# NOTE: nGraph's "convolution" is NNlib's crosscorrelation
#
# Need to flip the W and H dimensions of the filters
function flip!(x::AbstractArray{T,N}) where {T,N}
    x .= view(x, size(x, 1):-1:1, size(x, 2):-1:1, ntuple(_->:, N-2)...)
end


text:
# uniform interface for model estimation

export Estimator, MLEstimator, MAPEstimator
export nsamples, estimate, prior_score

abstract Estimator{D<:Distribution}

nsamples{D<:UnivariateDistribution}(e::Estimator{D}, x::Array) = length(x)
nsamples{D<:MultivariateDistribution}(e::Estimator{D}, x::Matrix) = size(x, 2)

type MLEstimator{D<:Distribution} <: Estimator{D} end
MLEstimator{D<:Distribution}(::Type{D}) = MLEstimator{D}()

estimate{D<:Distribution}(e::MLEstimator{D}, x) = fit_mle(D, x)
estimate{D<:Distribution}(e::MLEstimator{D}, x, w) = fit_mle(D, x, w)

prior_score{D<:Distribution}(e::MLEstimator{D}, d::D) = 0.

immutable MAPEstimator{D<:Distribution,Pri} <: Estimator{D}
	pri::Pri
end
MAPEstimator{D<:Distribution,Pri}(::Type{D}, pri::Pri) = MAPEstimator{D,Pri}(pri)

estimate{D<:Distribution,Pri}(e::MAPEstimator{D,Pri}, x) = fit_map(e.pri, D, x)
estimate{D<:Distribution,Pri}(e::MAPEstimator{D,Pri}, x, w) = fit_map(e.pri, D, x, w)


text:
import ThreeJS
import Colors: distinguishable_colors
main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    ts = linspace(0,8π,250);
    x = map(cos,ts)
    y = 0.1ts .* map(sin,ts)
    z = linspace(0,100,250)
    vertices = Tuple{Float64, Float64, Float64}[(x[i]*25,z[i],y[i]*25) for i in 1:size(x,1)]
    ThreeJS.outerdiv() <<
            (ThreeJS.initscene() <<
                [
                    ThreeJS.line(vertices) <<
                    [
                        ThreeJS.linematerial(Dict(:color=>"blue"))
                    ],
                    ThreeJS.pointcloud(collect(zip(y*25,z,x*25)), vertexcolors=distinguishable_colors(size(x,1))) <<
                    [
                        ThreeJS.pointmaterial(Dict(:color=>"white", :texture=>"/assets/disc.png", :transparent=>true, :size=>5.0, :colorkind=>"vertex", :alphatest=>0.5))
                    ],
                    ThreeJS.ambientlight(),
                    ThreeJS.camera(0.0, 0.0, 400.0)
                ]
            )
end


text:
"""
Wrapper around an objective function (compatible with NLSolversBase).
"""
mutable struct EvolutionaryObjective{TC,TF,TX,TP} <: AbstractObjective
    f::TC
    F::TF
    x_f::TX
    f_calls::Int
end

"""
    EvolutionaryObjective(f, x[, F])

Constructor for an objective function object around the function `f` with initial paramter `x`, and objective value `F`.
"""
function EvolutionaryObjective(f, x::AbstractArray,
                               F::Union{Real, AbstractArray{<:Real}} = zero(f(x));
                               eval::Symbol = :serial)
    defval = default_values(x)
    # convert function into the in-place one
    fn = if funargnum(f) == 2 && F isa AbstractArray
        (Fv,xv) -> (Fv .= f(xv))
    else
        f
    end
    EvolutionaryObjective{typeof(fn),typeof(F),typeof(x),Val{eval}}(fn, F, defval, 0)
end

"""
    EvolutionaryObjective(f, x::Expr[, F])

Constructor for an objective object for a Julia evaluatable expression.
"""
function EvolutionaryObjective(f, x::Expr, F::TF = zero(f(x));
                               eval::Symbol = :serial) where {TF<:Real}
    EvolutionaryObjective{typeof(f),TF,typeof(x),Val{eval}}(f, F, :(), 0)
end

f_calls(obj::EvolutionaryObjective) = obj.f_calls

"""
    ismultiobjective(objfun)

Return `true` if the function is multiobjective objective.
"""
ismultiobjective(obj) = obj.F isa AbstractArray

function value(obj::EvolutionaryObjective{TC,TF,TX,TP}, x::TX) where {TC,TF,TX,TP}
    obj.f_calls += 1
    obj.f(x)::TF
end

function value!(obj::EvolutionaryObjective{TC,TF,TX,TP}, x::TX) where {TC,TF,TX,TP}
    obj.F = value(obj, x)
    obj.F
end

function value!!(obj::EvolutionaryObjective{TC,TF,TX,TP}, x::TX) where {TC,TF,TX,TP}
    copyto!(obj.x_f, x)
    value!(obj, x)
end

function value(obj::EvolutionaryObjective{TC,TF,TX,TP}, F, x::TX) where {TC,TF,TX,TP}
    obj.f_calls += 1
    obj.f(F, x)
end
value(obj::EvolutionaryObjective{TC,TF,TX,TP}, x::TX) where {TC,TF<:AbstractArray,TX,TP} = value(obj, copy(obj.F), x)

function value!(obj::EvolutionaryObjective{TC,TF,TX,TP}, F, x::TX) where {TC,TF,TX,TP}
    obj.F = value(obj, F, x)
end

function value!!(obj::EvolutionaryObjective{TC,TF,TX,TP}, F, x::TX) where {TC,TF,TX,TP}
    copyto!(obj.x_f, x)
    value!(obj, F, x)
end

function value!(obj::EvolutionaryObjective{TC,TF,TX,Val{:serial}},
                F::AbstractVector, xs::AbstractVector{TX}) where {TC,TF<:Real,TX}
    broadcast!(x->value(obj,x), F, xs)
    F
end

function value!(obj::EvolutionaryObjective{TC,TF,TX,Val{:thread}},
                F::AbstractVector, xs::AbstractVector{TX}) where {TC,TF<:Real,TX}
    n = length(xs)
    Threads.@threads for i in 1:n
        F[i] = value(obj, xs[i])
    end
    F
end

function value!(obj::EvolutionaryObjective{TC,TF,TX,Val{:serial}},
                F::AbstractMatrix, xs::AbstractVector{TX}) where {TC,TF,TX}
    n = length(xs)
    for i in 1:n
        fv = view(F, :, i)
        value(obj, fv, xs[i])
    end
    F
end

function value!(obj::EvolutionaryObjective{TC,TF,TX,Val{:thread}},
                F::AbstractMatrix, xs::AbstractVector{TX}) where {TC,TF,TX}
    n = length(xs)
    @Threads.threads for i in 1:n
        fv = view(F, :, i)
        value(obj, fv, xs[i])
    end
    F
end


text:
"""
    listholidays(calendar, dt0::Date, dt1::Date) → Vector{Date}

Returns the list of holidays between `dt0` and `dt1`.
"""
function listholidays(hc::HolidayCalendar, dt0::Date, dt1::Date)
    const d0 = min(dt0, dt1)
    const d1 = max(dt0, dt1)
    const dt_range = d0:d1
    isholiday_vec = [ isholiday(hc, i) for i in dt_range ]
    return dt_range[isholiday_vec]
end

listholidays(calendar, dt0::Date, dt1::Date) = listholidays(convert(HolidayCalendar,calendar), dt0, dt1)

"""
    listbdays(calendar, dt0::Date, dt1::Date) → Vector{Date}

Returns the list of business days between `dt0` and `dt1`.
"""
function listbdays(hc::HolidayCalendar, dt0::Date, dt1::Date)
    d = tobday(hc,min(dt0,dt1))
    const d1 = max(dt0, dt1)

    # empty result
    if d > d1
        return Array(Date,0)
    end

    const n = convert(Int, d1-d)+1
    raw_vec::Array{Date, 1} = Array(Date, n)
    raw_vec[1] = d
    i::Int = 2
    d = advancebdays(hc,d,1)
    while d <= d1
        raw_vec[i] = d
        i+=1
        d = advancebdays(hc,d,1)
    end
    return raw_vec[1:(i-1)]
end

listbdays(calendar, dt0::Date, dt1::Date) = listbdays(convert(HolidayCalendar,calendar), dt0, dt1)


text:
abstract type Adjuster{T <: Union{Supply, Demand}} end

abstract type SupplyAdjuster{Supply <: Adjuster} end
abstract type DemandAdjuster{Demand <: Adjuster} end

struct SafetyStockAdjuster2{SupplyAdjuster}
    quantity::Int

    SafetyStockAdjuster2{SupplyAdjuster}(quantity) = new(quantity)
end

struct DemandValueAdjuster{Demand} <: DemandAdjuster{Demand}
    quantity::Int
end

varSupplyAdjusters = [
    SafetyStockAdjuster2{SupplyAdjuster}(20),
    SafetyStockAdjuster2{SupplyAdjuster}(20)]

function applySupplyAdjusters(adj::Array)
    # println("start - Applying SafetyStock Adjuster $(adj)")

    out = []

    for val in adj
        for adjuster in varSupplyAdjusters
            val -= adjuster.quantity
        end
        push!(out, val)
    end

    # println("End - Applying SafetyStock Adjuster $(out)")

    return out
end

function applyDemandAdjusters(adj::DemandValueAdjuster)
    println("Applying SafetyStock Adjuster")
end

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
using NewtonMethods
using LinearAlgebra: Diagonal, I, Symmetric
using Roots: fzero

include("types.jl")
include("base.jl")
include("invert.jl")

println("Running test...")
println("Make synthetic data...")

m0 = [2,1.2,0.1,0.15,0.03,deg2rad(28), deg2rad(39), deg2rad(58), deg2rad(90),deg2rad(58),deg2rad(90)]
ρ = 1.

C0 = thomsen2moduli(m0[1],m0[2],ρ,m0[3],m0[4],m0[5])

np=4
nsv=0
nsh=2

args = (np, nsv, nsh)

dobs = _gt(m0, np, nsv, nsh)

p = [VMeasure(dobs[2k-1],exp(dobs[2k]),.01,0.01) for k in 1:np]
sv = [VMeasure(dobs[2np+2k-1],exp(dobs[2np+2k]),.01,0.01) for k in 1:nsv]
sh = [VMeasure(dobs[2np+2nsv+2k-1],exp(dobs[2np+2nsv+2k]),.01,0.01) for k in 1:nsh]

# a prior moel parameters
m1 = [2.2,1.5,0.0,0.0,0.,deg2rad(25), deg2rad(40), deg2rad(55), deg2rad(90),deg2rad(90)]

C1 = thomsen2moduli(m1[1],m1[2],ρ,m1[3],m1[4],m1[5])

# run inversion thomsen
println("Run inversion to find Thomsen parameters...")

m, anglep, anglesv, anglesh, cm = findthomsen(p,sv,sh, m1[1], m1[2], m1[3], m1[4], m1[5])

println()
println("Model parameters: input | inverted")
show(IOContext(stdout), "text/plain",[m0[1:5]  m])
println()
println()
dcalc = _gt(vcat(m, anglep, anglesv, anglesh), np, nsv, nsh)
println("Observations: obs | calc")
show(IOContext(stdout), "text/plain",[dobs  dcalc])
println()

# run inversion moduli
println("Run inversion to find moduli...")

C, anglep, anglesv, anglesh, cm = findmoduli(p,sv,sh, ρ, C1[1,1], C1[1,3], C1[3,3], C1[5,5], C1[6,6], 100)

m = vcat([C[1,1], C[1,3], C[3,3], C[5,5], C[6,6]], anglep, anglesv, anglesh)

println()
println("Model parameters: input | inverted")
show(IOContext(stdout), "text/plain",C0)
println("\n------")
show(IOContext(stdout), "text/plain",C)
println()
println()
dcalc = _gc(m, ρ, np, nsv, nsh)
println("Observations: obs | calc")
show(IOContext(stdout), "text/plain",[dobs  dcalc])


text:
using QuBase
using Compat
using Base.Test

include("arraytest.jl")
include("multest.jl")
include("tensortest.jl")
include("labelbasistest.jl")
include("operatortest.jl")


text:
# https://docs.julialang.org/en/stable/manual/types/#man-custom-pretty-printing-1
# Base.show(io::IO, z): single line format used in show, print, inside other objects.
# Base.show(io::IO, ::MIME"text/plain", z): multi-line format used by display.
# Base.show(io::IO, ::MIME"text/html", z): multi-line format for html output.
# get(io, :compact, false), show(IOContext(stdout, :compact=>true),z) for compact (array) printing.
# summary(io::IO, x) = print(io, typeof(x))
# string(z): uses print_to_string.

import Base: show

# Summary string
valstr(x)=_valstr(value(x))
_valstr(::Nothing)="nothing"
_valstr(x)=(hasmethod(size,(typeof(x),)) && !isempty(size(x)) ? "$(typeof(x))($(join(size(x),',')))" : string(x))
_valstr(t::Tuple)="($(join(_valstr.(t),',')))"

# One line show used for show, print, string etc.
show(io::IO, x::Bcasted)  = print(IOContext(io,:compact=>true), "B(", valstr(x), ")")
show(io::IO, x::Param)    = print(IOContext(io,:compact=>true), "P(", valstr(x), ")")
show(io::IO, x::Result)   = print(IOContext(io,:compact=>true), "R(", valstr(x), ")")
show(io::IO, x::Tape)     = print(IOContext(io,:compact=>true), "T(", valstr(x), ")")
show(io::IO, x::Sparse)   = print(IOContext(io,:compact=>true), "Sparse(", valstr(x.container), (size.(x.values)...), ")")

# Multi line show used for display:
show(io::IO, ::MIME"text/plain", x::Tape) = show(io, x)
show(io::IO, ::MIME"text/plain", x::Sparse) = show(io, x)

# Hack to take advantage of array display:
import Base: summary, size, getindex
struct ArrayValue{T,N} <: AbstractArray{T,N}; p; end
show(io::IO, m::MIME"text/plain", x::Value{A}) where {A<:AbstractArray} = show(io, m, ArrayValue{eltype(x),ndims(x)}(x))
size(p::ArrayValue) = size(p.p.value)
getindex(p::ArrayValue,i...) = getindex(p.p.value,i...)
summary(io::IO, x::Value{A}) where {A<:AbstractArray} = print(io, Base.dims2string(size(x)), " ", typeof(x))
summary(io::IO, p::ArrayValue) = summary(io, p.p)

function show(io::IO, n::Node)
    og(n::Node)=(n.outgrad === nothing ? 0 : valstr(n.outgrad))
    io = IOContext(io,:compact=>true)
    if isa(n.Value, Param)
        print(io, "N(", og(n), ", ", n.Value, ")")
    else
        r = n.Value
        print(io, "N(", og(n), ", ", r, " = ", r.func, "(", join(valstr.(r.args),", "),
              isempty(r.kwargs) ? "" : "; "*join(["$(x[1])=$(valstr(x[2]))" for x in r.kwargs], ", "), "))")
    end
end

function show(io::IO, ::MIME"text/plain", ts::Vector{Tape}) # to dump _tapes
    if isempty(ts); show(io, ts); return; end
    og(t::Tape,r::Value)=(n=get(t.dict,r,nothing); n===nothing ? '-' : n.outgrad===nothing ? '0' : valstr(n.outgrad))
    argstr(x)=(n=findfirst(a->(a.Value===x),ts[1].list); n===nothing ? valstr(x) : "R$n")
    io = IOContext(io,:compact=>true)
    for (i,n) in enumerate(ts[1].list)
        r = n.Value
        if isa(r,Result)
            print(io, "$i. ", valstr(r), " = ", r.func, "(", join(argstr.(r.args),", "),
                  isempty(r.kwargs) ? "" : "; "*join(["$(x[1])=$(valstr(x[2]))" for x in r.kwargs], ", "), "))")
        else
            print(io, "$i. ", r)
        end
        println(io, " ∇=", join([og(t,r) for t in ts], '|'))
    end
end


text:
#
function q_Sigma_cd_grid(q::Int, c_grid::Array{Float64,1}, d_grid::Array{Float64,1})
    mkΣ_T = Int(1000)#numerical integration grid
    Tmax = mkΣ_T + 2
# Compute cosine transforms and weights
  fvec = collect(π*(1:1:q))'
  tvec = collect((0.5:1:mkΣ_T)/mkΣ_T)
  v = zeros(Tmax, q)
  v[1:mkΣ_T, 1:q] = sqrt(2)*cos(tvec*fvec)
  v_1 = v
#
  x = (1/mkΣ_T)*collect(1:Tmax-1)
  n_c = length(c_grid)
  n_d = length(d_grid)
  qΣcd_grid = NaN*zeros(n_c, n_d, q, q)
    #
  for ic in 1:n_c
    println("c iter $ic of $n_c")
    c = c_grid[ic]
    (c == 0) && (c = 0.01) # corner case
    cx = c*x
    # Summation by parts formula used in d0 < 0.5;
    vqd = NaN*zeros(Tmax, q)
    vqd[1:Tmax-1, :] = (c/mkΣ_T)*v_1[1:Tmax-1, :] - (v_1[2:Tmax, :] - v_1[1:Tmax-1, :])
    vqd[1, :] = vqd[1, :] - v_1[1, :]
    vqd[Tmax, :] = 0
    v_0 = mkΣ_T*vqd
    for id in 1:n_d
      d = d_grid[id]
      vqd = NaN*zeros(Tmax, q)
      d0 = d;
      (d < 0.5) && (d = d+1)
      (d == 0.5) && (d = 0.51)
      if (d0 < 0.5)
        v = v_0
      else
        v = v_1
      end
      #
      ga_0 = c^(1-2*d)*sqrt(π)*gamma(-0.5+d)/gamma(d)
      bselval = besselk(d-0.5, cx)
      ga = ((2^(1.5-d))*sqrt(π)/gamma(d))*((x/c).^(d-0.5)).*bselval
      vqd[1, :] = [ga_0 ga']*v
      for t in 1:Tmax-2
        vqd[t+1, :] = [flipdim(ga[1:t], 1)' ga_0 ga[1:Tmax-t-1]']*v
      end
      vqd[Tmax, :] = [flipdim(ga, 1)' ga_0]*v
      qΣcd = v'*vqd/(mkΣ_T^2)
      qΣcd = qΣcd/(2*π)
      #Normalize to have average value equal to unity
      tmp = mean(diag(qΣcd))
      qΣcd = qΣcd/tmp
      qΣcd_grid[ic, id, :, :] = qΣcd
    end
  end
  return qΣcd_grid
end
#
function interp_q_Sigma_cd_grid(qΣcd_grid, c, d, c_grid, d_grid)
  q = size(qΣcd_grid, 4)
  qΣcd = NaN*zeros(q,q)
  for i in 1:q
    for j in 1:q
      tmp = qΣcd_grid[:, :, i, j] # these are the values of the covariance between transforms i and j for each value of c and d
      spl = Spline2D(c_grid, d_grid, tmp)
      qΣcd[i, j] = evalgrid(spl, [c], [d])[1]
    end
  end
  tmp = mean(diag(qΣcd))
  qΣcd = qΣcd/tmp
  return qΣcd
end
#
function estimate_d(lr::bcdModel)
  q = lr.q
  cvec = lr.c.grid
  dvec = lr.d.grid
  Sigma_d_mat=NaN*zeros(lr.q, lr.q, length(dvec))
  Sigma_inv_d_mat = Sigma_d_mat
  #"Estimate" d via grid through likelihood calibration
  lvec = NaN*zeros(length(dvec))
  X_norm = lr.cost
  for i in 1:length(dvec)
    tic()
    d = dvec[i]
    Sigma = Sigma_Compute2(0, 0, d, lr.q)
    sig_inv = inv(Sigma)
    lvec[i, :] = den_invariant(X_norm, sig_inv)
    Sigma_d_mat[:, :, i] = Sigma
    Sigma_inv_d_mat[:, :, i] = sig_inv
    toc()
  end
  lr_stat = mean(lvec)./lvec
  pct_vec = 1-[2/3 0.9 0.95]'
  lr_d_cv = NaN*zeros(length(dvec), length(pct_vec))
  nsim = 50000
  Ysim = randn(lr.q, nsim)
  lvec_sim = NaN*zeros(size(dvec,1),nsim);
  for i in 1:length(dvec)
    tic()
    sig_ho = Sigma_d_mat[:, :, i]
    chol_ho = chol(Hermitian(sig_ho))
    Xsim = chol_ho*Ysim
    tmp = sqrt(sum(Xsim.^2, 1))
    Xsim_norm = Xsim./repmat(tmp,q,1);
    lvec_sim = NaN*zeros(size(dvec,1),nsim);
    for j = 1:size(dvec,1);
         sig_inv = Sigma_inv_d_mat[:,:,j];
         lvec_sim[j,:] = den_invariant(Xsim_norm,sig_inv);
     end;
     lr_stat_ho = mean(lvec_sim)./lvec_sim[i,:];
     tmp = pctile(lr_stat_ho',pct_vec);
     lr_d_cv[i, :] = tmp';
     toc()
   end
   tmp, ii = findmax(lvec)
   println("MLE (from grid) for d: $(dvec[ii])")
   condmat = repmat(log(lr_stat), 1, size(lr_d_cv, 2)) .> log(lr_d_cv)
   println("Confidence sets for d")
   cset = repmat(pct_vec, 1, 2)
   for i in 1:length(pct_vec)
     tmp = dvec[condmat[:, i].==true]
     #println(" Level = $(pct_vec[i]): $(tmp[1]) $(tmp[end])")
     #cset[i, :] = [tmp[1] tmp[end]]
   end
   #tmp, ii = findmin(abs(dvec-0))
   #lvec = lvec
   lr.d.val = dvec[ii]
   lr.stuf[:lvec] = lvec
   lr.d.lik = lvec
  #return lvec, lr_stat, Sigma_d_mat, Sigma_inv_d_mat, lr_d_cv, lvec_sim, cset
  return lr
end
#
function Sigma_Compute2(b, c, d, q; iconst = false)
  T_general = 150; T_special = 600;
  cov_Id = 0
  if c < 0.0001
    T = T_special;
    cov_Id = frac_cov(T,d)
  else
    if d == 1
      T = T_special;
      rho = 1 - c/T;
      rho = max(0, rho);
      cov_Id, tmp = ar1_cov(T,rho);
    else
      error("not implemented")
      T = T_general;
      rho = 1 - c/T;
      rho = max(0, rho);
      #cov_Id = frac_dif_cov(T, rho, d);
      cov_Id, tmp = ar1_cov(T,rho);
    end
  end
  #Add I(0) components
  scale_I0 = (b*(T^d))^2;
  #println(cov_Id); println(scale_I0); println(eye(T))
  covmat = cov_Id + scale_I0*eye(T);
  psi = psi_compute(T, q);
  if iconst == true
    psi = [ones(T,1)/T psi];
  end
  # Step 3: compute Sigma
  Sigma = psi'*covmat*psi;
  expoente = 1. - 2*d;
  Sigma = (T^expoente)*Sigma;
  #= Verify Sigma is positive definite
  ier = 0;
  tmp, tmp2 = eig(Sigma);
  if minimum(tmp) < 0
    ier = 1;
    Sigma_2 = Sigma*Sigma';
    V, D2 = eig(Sigma_2);
    D_vec = sqrt(diag(D2));
    D_mat = diagm(D_vec);
    Sigma = V*D_mat*V';
  end =#
  return Sigma
end
#


text:
"""
    LeviCivita(T::Type{<:Number} = Sym)

Builds an Array{T,3} of Levi-Civita Symbol `ϵᵢⱼₖ = (i-j) (j-k) (k-i) / 2`

# Examples
```julia
julia> ε = LeviCivita(Sym)
3×3×3 Array{Sym, 3}:
[:, :, 1] =
 0   0  0
 0   0  1
 0  -1  0

[:, :, 2] =
 0  0  -1
 0  0   0
 1  0   0

[:, :, 3] =
  0  1  0
 -1  0  0
  0  0  0
```
"""
LeviCivita(T::Type{<:Number} = Sym) = [T(T((i - j) * (j - k) * (k - i)) / T(2)) for i = 1:3, j = 1:3, k = 1:3]

"""
    𝐞(i::Integer, dim::Int = 3, T::Type{<:Number} = Sym)

Vector of the canonical basis

# Examples
```julia
julia> 𝐞(1)
Tens{1, 3, Sym, Sym, Vec{3, Sym}, CanonicalBasis{3, Sym}}
# data: 3-element Vec{3, Sym}:
 1
 0
 0
# var: (:cont,)
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 1  0  0
 0  1  0
 0  0  1
```
"""
𝐞(::Val{i}, ::Val{dim} = Val(3), ::Val{T} = Val(Sym)) where {i, dim, T<:Number} =
    Tens(Vec{dim}(j -> j == i ? one(T) : zero(T)))

"""
    init_cartesian(coords = symbols("x y z", real = true))

Returns the coordinates, unit vectors and basis of the cartesian basis

# Examples
```julia
julia> coords, vectors, ℬ = init_cartesian() ; x, y, z = coords ; 𝐞₁, 𝐞₂, 𝐞₃ = vectors ;
```
"""
init_cartesian(coords = symbols("x y z", real = true)) = Tuple(coords),
ntuple(i -> 𝐞(Val(i), Val(length(coords)), Val(eltype(coords))), length(coords)),
CanonicalBasis{length(coords),eltype(coords)}()

init_cartesian(::Val{3}) = init_cartesian(symbols("x y z", real = true))
init_cartesian(::Val{2}) = init_cartesian(symbols("x y", real = true))
init_cartesian(dim::Integer) = init_cartesian(Val(dim))

"""
    𝐞ᵖ(i::Integer, θ::T = zero(Sym); canonical = false)

Vector of the polar basis

# Examples
```julia
julia> θ = symbols("θ", real = true) ;

julia> 𝐞ᵖ(1, θ)
Tens{1, 2, Sym, Sym, Vec{2, Sym}, RotatedBasis{2, Sym}}
# data: 2-element Vec{2, Sym}:
 1
 0
# var: (:cont,)
# basis: 2×2 Tensor{2, 2, Sym, 4}:
 cos(θ)  -sin(θ)
 sin(θ)   cos(θ)
```
"""
𝐞ᵖ(::Val{1}, θ::T = zero(Sym); canonical = false) where {T<:Number} =
    canonical ? Tens(Vec{2}([cos(θ), sin(θ)])) :
    Tens(Vec{2}([one(T), zero(T)]), Basis(θ))
𝐞ᵖ(::Val{2}, θ::T = zero(Sym); canonical = false) where {T<:Number} =
    canonical ? Tens(Vec{2}([-sin(θ), cos(θ)])) :
    Tens(Vec{2}([zero(T), one(T)]), Basis(θ))

"""
    init_polar(coords = (symbols("r", positive = true), symbols("θ", real = true)); canonical = false)

Returns the coordinates, base vectors and basis of the polar basis

# Examples
```julia
julia> coords, vectors, ℬᵖ = init_polar() ; r, θ = coords ; 𝐞ʳ, 𝐞ᶿ = vectors ;
```
"""
init_polar(
    coords = (symbols("r", positive = true), symbols("θ", real = true));
    canonical = false,
) = Tuple(coords),
ntuple(i -> 𝐞ᵖ(Val(i), coords[2]; canonical = canonical), 2),
Basis(coords[2])

"""
    𝐞ᶜ(i::Integer, θ::T = zero(Sym); canonical = false)

Vector of the cylindrical basis

# Examples
```julia
julia> θ = symbols("θ", real = true) ;

julia> 𝐞ᶜ(1, θ)
Tens{1, 3, Sym, Sym, Vec{3, Sym}, RotatedBasis{3, Sym}}
# data: 3-element Vec{3, Sym}:
 1
 0
 0
# var: (:cont,)
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 cos(θ)  -sin(θ)  0
 sin(θ)   cos(θ)  0
      0        0  1
```
"""
𝐞ᶜ(::Val{1}, θ::T = zero(Sym); canonical = false) where {T<:Number} =
    canonical ? Tens(Vec{3}([cos(θ), sin(θ), zero(T)])) :
    Tens(Vec{3}([one(T), zero(T), zero(T)]), CylindricalBasis(θ))
𝐞ᶜ(::Val{2}, θ::T = zero(Sym); canonical = false) where {T<:Number} =
    canonical ? Tens(Vec{3}([-sin(θ), cos(θ), zero(T)])) :
    Tens(Vec{3}([zero(T), one(T), zero(T)]), CylindricalBasis(θ))
𝐞ᶜ(::Val{3}, θ::T = zero(Sym); canonical = false) where {T<:Number} =
    canonical ? Tens(Vec{3}([zero(T), zero(T), one(T)])) :
    Tens(Vec{3}([zero(T), zero(T), one(T)]), CylindricalBasis(θ))

"""
    init_cylindrical(coords = (symbols("r", positive = true), symbols("θ", real = true), symbols("z", real = true)); canonical = false)

Returns the coordinates, base vectors and basis of the cylindrical basis

# Examples
```julia
julia> coords, vectors, ℬᶜ = init_cylindrical() ; r, θ, z = coords ; 𝐞ʳ, 𝐞ᶿ, 𝐞ᶻ = vectors ;
```
"""
init_cylindrical(
    coords = (
        symbols("r", positive = true),
        symbols("θ", real = true),
        symbols("z", real = true),
    );
    canonical = false,
) = Tuple(coords),
ntuple(i -> 𝐞ᶜ(Val(i), coords[2]; canonical = canonical), 3),
CylindricalBasis(coords[2])

"""
    𝐞ˢ(i::Integer, θ::T = zero(Sym), ϕ::T = zero(Sym), ψ::T = zero(Sym); canonical = false)

Vector of the basis rotated with the 3 Euler angles `θ, ϕ, ψ` (spherical if `ψ=0`)

# Examples
```julia
julia> θ, ϕ, ψ = symbols("θ, ϕ, ψ", real = true) ;

Tens{1, 3, Sym, Sym, Vec{3, Sym}, RotatedBasis{3, Sym}}
# data: 3-element Vec{3, Sym}:
 1
 0
 0
# var: (:cont,)
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 -sin(ψ)⋅sin(ϕ) + cos(θ)⋅cos(ψ)⋅cos(ϕ)  -sin(ψ)⋅cos(θ)⋅cos(ϕ) - sin(ϕ)⋅cos(ψ)  sin(θ)⋅cos(ϕ)
  sin(ψ)⋅cos(ϕ) + sin(ϕ)⋅cos(θ)⋅cos(ψ)  -sin(ψ)⋅sin(ϕ)⋅cos(θ) + cos(ψ)⋅cos(ϕ)  sin(θ)⋅sin(ϕ)
                        -sin(θ)⋅cos(ψ)                          sin(θ)⋅sin(ψ)         cos(θ)
```
"""
function 𝐞ˢ(
    ::Val{1},
    θ::T1 = 0,
    ϕ::T2 = 0,
    ψ::T3 = 0;
    canonical = false,
) where {T1<:Number,T2<:Number,T3<:Number}
    if canonical
        return Tens(
            Vec{3}([
                -sin(ψ) * sin(ϕ) + cos(θ) * cos(ψ) * cos(ϕ),
                sin(ψ) * cos(ϕ) + sin(ϕ) * cos(θ) * cos(ψ),
                -sin(θ) * cos(ψ),
            ]),
        )
    else
        T = promote_type(T1, T2, T3)
        return Tens(Vec{3}([one(T), zero(T), zero(T)]), Basis(θ, ϕ, ψ))
    end
end
function 𝐞ˢ(
    ::Val{2},
    θ::T1 = 0,
    ϕ::T2 = 0,
    ψ::T3 = 0;
    canonical = false,
) where {T1<:Number,T2<:Number,T3<:Number}
    if canonical
        return Tens(
            Vec{3}([
                -sin(ψ) * cos(θ) * cos(ϕ) - sin(ϕ) * cos(ψ),
                -sin(ψ) * sin(ϕ) * cos(θ) + cos(ψ) * cos(ϕ),
                sin(θ) * sin(ψ),
            ]),
        )
    else
        T = promote_type(T1, T2, T3)
        return Tens(Vec{3}([zero(T), one(T), zero(T)]), Basis(θ, ϕ, ψ))
    end
end
function 𝐞ˢ(
    ::Val{3},
    θ::T1 = 0,
    ϕ::T2 = 0,
    ψ::T3 = 0;
    canonical = false,
) where {T1<:Number,T2<:Number,T3<:Number}
    if canonical
        return Tens(Vec{3}([sin(θ) * cos(ϕ), sin(θ) * sin(ϕ), cos(θ)]))
    else
        T = promote_type(T1, T2, T3)
        return Tens(Vec{3}([zero(T), zero(T), one(T)]), Basis(θ, ϕ, ψ))
    end
end

for eb in (:𝐞, :𝐞ᵖ, :𝐞ᶜ, :𝐞ˢ)
    @eval $eb(i::Integer, args...; kwargs...) = $eb(Val(i), args...; kwargs...)
end

"""
    init_spherical(coords = (symbols("θ", real = true), symbols("ϕ", real = true), symbols("r", positive = true)); canonical = false)

Returns the coordinates, base vectors and basis of the spherical basis.
Take care that the order of the 3 vectors is `𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ` so that
the basis coincides with the canonical one when the angles are null and in consistency
the coordinates are ordered as `θ, ϕ, r`.

# Examples
```julia
julia> coords, vectors, ℬˢ = init_spherical() ; θ, ϕ, r = coords ; 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ  = vectors ;
```
"""
init_spherical(
    coords = (
        symbols("θ", real = true),
        symbols("ϕ", real = true),
        symbols("r", positive = true),
    );
    canonical = false,
) = Tuple(coords),
ntuple(i -> 𝐞ˢ(Val(i), coords[1:2]...; canonical = canonical), 3),
SphericalBasis(coords[1:2]...)

"""
    init_rotated(coords = symbols("θ ϕ ψ", real = true); canonical = false)

Returns the angles, base vectors and basis of the rotated basis.
Note that here the coordinates are angles and do not represent a valid parametrization of `ℝ³`

# Examples
```julia
julia> angles, vectors, ℬʳ = init_rotated() ; θ, ϕ, ψ = angles ; 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = vectors ;
```
"""
init_rotated(angles = symbols("θ ϕ ψ", real = true); canonical = false) = Tuple(angles),
ntuple(i -> 𝐞ˢ(Val(i), angles...; canonical = canonical), 3),
Basis(angles...)

"""
    rot3(θ, ϕ = 0, ψ = 0)

Returns a rotation matrix with respect to the 3 Euler angles `θ, ϕ, ψ`

# Examples
```julia
julia> cθ, cϕ, cψ, sθ, sϕ, sψ = symbols("cθ cϕ cψ sθ sϕ sψ", real = true) ;

julia> d = Dict(cos(θ) => cθ, cos(ϕ) => cϕ, cos(ψ) => cψ, sin(θ) => sθ, sin(ϕ) => sϕ, sin(ψ) => sψ) ;

julia> subs.(rot3(θ, ϕ, ψ),d...)
3×3 StaticArrays.SMatrix{3, 3, Sym, 9} with indices SOneTo(3)×SOneTo(3):
 cθ⋅cψ⋅cϕ - sψ⋅sϕ  -cθ⋅cϕ⋅sψ - cψ⋅sϕ  cϕ⋅sθ
 cθ⋅cψ⋅sϕ + cϕ⋅sψ  -cθ⋅sψ⋅sϕ + cψ⋅cϕ  sθ⋅sϕ
           -cψ⋅sθ              sθ⋅sψ     cθ
```
"""
rot3(θ, ϕ = 0, ψ = 0) = RotZYZ(ϕ, θ, ψ)

"""
    rot2(θ)

Returns a 2D rotation matrix with respect to the angle `θ`

# Examples
```julia
julia> rot2(θ)
2×2 Tensor{2, 2, Sym, 4}:
 cos(θ)  -sin(θ)
 sin(θ)   cos(θ)
```
"""
rot2(θ) = Tensor{2,2}((cos(θ), sin(θ), -sin(θ), cos(θ)))

"""
    rot6(θ, ϕ = 0, ψ = 0)

Returns a rotation matrix with respect to the 3 Euler angles `θ, ϕ, ψ`

# Examples
```julia
julia> cθ, cϕ, cψ, sθ, sϕ, sψ = symbols("cθ cϕ cψ sθ sϕ sψ", real = true) ;

julia> d = Dict(cos(θ) => cθ, cos(ϕ) => cϕ, cos(ψ) => cψ, sin(θ) => sθ, sin(ϕ) => sϕ, sin(ψ) => sψ) ;

julia> R = Tens(subs.(rot3(θ, ϕ, ψ),d...))
Tens.TensCanonical{2, 3, Sym, Tensor{2, 3, Sym, 9}}
# data: 3×3 Tensor{2, 3, Sym, 9}:
 cθ⋅cψ⋅cϕ - sψ⋅sϕ  -cθ⋅cϕ⋅sψ - cψ⋅sϕ  cϕ⋅sθ
 cθ⋅cψ⋅sϕ + cϕ⋅sψ  -cθ⋅sψ⋅sϕ + cψ⋅cϕ  sθ⋅sϕ
           -cψ⋅sθ              sθ⋅sψ     cθ
# var: (:cont, :cont)
# basis: 3×3 Tens.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1

julia> RR = R ⊠ˢ R
Tens.TensCanonical{4, 3, Sym, SymmetricTensor{4, 3, Sym, 36}}
# data: 6×6 Matrix{Sym}:
                          (cθ*cψ*cϕ - sψ*sϕ)^2                            (-cθ*cϕ*sψ - cψ*sϕ)^2           cϕ^2*sθ^2                      √2⋅cϕ⋅sθ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)                     √2⋅cϕ⋅sθ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)                                   √2⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)
                          (cθ*cψ*sϕ + cϕ*sψ)^2                            (-cθ*sψ*sϕ + cψ*cϕ)^2           sθ^2*sϕ^2                      √2⋅sθ⋅sϕ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)                     √2⋅sθ⋅sϕ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)                                   √2⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)
                                     cψ^2*sθ^2                                        sθ^2*sψ^2                cθ^2                                       √2⋅cθ⋅sθ⋅sψ                                    -√2⋅cθ⋅cψ⋅sθ                                                              -sqrt(2)*cψ*sθ^2*sψ
             -√2⋅cψ⋅sθ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)                √2⋅sθ⋅sψ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)    √2⋅cθ⋅sθ⋅sϕ                    cθ*(-cθ*sψ*sϕ + cψ*cϕ) + sθ^2*sψ*sϕ                   cθ*(cθ*cψ*sϕ + cϕ*sψ) - cψ*sθ^2*sϕ                            -cψ⋅sθ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ) + sθ⋅sψ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)
             -√2⋅cψ⋅sθ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)                √2⋅sθ⋅sψ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)    √2⋅cθ⋅cϕ⋅sθ                    cθ*(-cθ*cϕ*sψ - cψ*sϕ) + cϕ*sθ^2*sψ                   cθ*(cθ*cψ*cϕ - sψ*sϕ) - cψ*cϕ*sθ^2                            -cψ⋅sθ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ) + sθ⋅sψ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)
 √2⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)  √2⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)  sqrt(2)*cϕ*sθ^2*sϕ  cϕ⋅sθ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ) + sθ⋅sϕ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)  cϕ⋅sθ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ) + sθ⋅sϕ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)  (cθ*cψ*cϕ - sψ*sϕ)*(-cθ*sψ*sϕ + cψ*cϕ) + (cθ*cψ*sϕ + cϕ*sψ)*(-cθ*cϕ*sψ - cψ*sϕ)
# var: (:cont, :cont, :cont, :cont)
# basis: 3×3 Tens.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1

julia> R6 = invKM(subs.(KM(rot6(θ, ϕ, ψ)),d...))
Tens.TensCanonical{4, 3, Sym, SymmetricTensor{4, 3, Sym, 36}}
# data: 6×6 Matrix{Sym}:
                          (cθ*cψ*cϕ - sψ*sϕ)^2                            (-cθ*cϕ*sψ - cψ*sϕ)^2           cϕ^2*sθ^2                      √2⋅cϕ⋅sθ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)                     √2⋅cϕ⋅sθ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)                                   √2⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)
                          (cθ*cψ*sϕ + cϕ*sψ)^2                            (-cθ*sψ*sϕ + cψ*cϕ)^2           sθ^2*sϕ^2                      √2⋅sθ⋅sϕ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)                     √2⋅sθ⋅sϕ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)                                   √2⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)
                                     cψ^2*sθ^2                                        sθ^2*sψ^2                cθ^2                                       √2⋅cθ⋅sθ⋅sψ                                    -√2⋅cθ⋅cψ⋅sθ                                                              -sqrt(2)*cψ*sθ^2*sψ
             -√2⋅cψ⋅sθ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)                √2⋅sθ⋅sψ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)    √2⋅cθ⋅sθ⋅sϕ                    cθ*(-cθ*sψ*sϕ + cψ*cϕ) + sθ^2*sψ*sϕ                   cθ*(cθ*cψ*sϕ + cϕ*sψ) - cψ*sθ^2*sϕ                            -cψ⋅sθ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ) + sθ⋅sψ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)
             -√2⋅cψ⋅sθ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)                √2⋅sθ⋅sψ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)    √2⋅cθ⋅cϕ⋅sθ                    cθ*(-cθ*cϕ*sψ - cψ*sϕ) + cϕ*sθ^2*sψ                   cθ*(cθ*cψ*cϕ - sψ*sϕ) - cψ*cϕ*sθ^2                            -cψ⋅sθ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ) + sθ⋅sψ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)
 √2⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ)  √2⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ)  sqrt(2)*cde Liv Lehn ϕ*sθ^2*sϕ  cϕ⋅sθ⋅(-cθ⋅sψ⋅sϕ + cψ⋅cϕ) + sθ⋅sϕ⋅(-cθ⋅cϕ⋅sψ - cψ⋅sϕ)  cϕ⋅sθ⋅(cθ⋅cψ⋅sϕ + cϕ⋅sψ) + sθ⋅sϕ⋅(cθ⋅cψ⋅cϕ - sψ⋅sϕ)  (cθ*cψ*cϕ - sψ*sϕ)*(-cθ*sψ*sϕ + cψ*cϕ) + (cθ*cψ*sϕ + cϕ*sψ)*(-cθ*cϕ*sψ - cψ*sϕ)
# var: (:cont, :cont, :cont, :cont)
# basis: 3×3 Tens.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1

julia> R6 == RR
true
```
"""
function rot6(θ, ϕ = 0, ψ = 0)
    R = TensCanonical(rot3(θ, ϕ, ψ))
    return sboxtimes(R, R)
end

export LeviCivita
export 𝐞, 𝐞ᵖ, 𝐞ᶜ, 𝐞ˢ
export init_cartesian, init_polar, init_cylindrical, init_spherical, init_rotated
export rot2, rot3, rot6


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module AddFields

using ExtractMacro
using ..Interface
using ..Common

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphAF, GraphAddFields, GraphAddSubFields

import ..Interface: energy, delta_energy, neighbors,
                    update_cache!, delta_energy_residual

struct GraphAF{ET} <: SimpleGraph{ET}
    N::Int
    fields::Vector{ET}
    function GraphAF(fields::Vector{ET}) where {ET}
        N = length(fields)
        return new{ET}(N, fields)
    end
end

@doc """
    GraphAF(fields::Vector) <: SimpleGraph

An auxiliary `SimpleGraph` used to implement additional fields on variables.

It is only useful when implementing other graph types; see [`GraphAddFields`](@ref).
""" -> GraphAF(fields::Vector)

function energy(X::GraphAF{ET}, C::Config) where {ET}
    @assert X.N == C.N
    @extract X : N fields
    @extract C : s
    E = zero(ET)
    for i = 1:N
        σi = 2s[i] - 1
        E += fields[i] * σi
    end
    return E
end

function delta_energy(X::GraphAF, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : σ = 2s[move]-1
    @extract X : h = fields[move]

    return -2 * σ * h
end

neighbors(X::GraphAF, i::Int) = ()

# Add external fields to any AbstractGraph

mutable struct GraphAddFields{ET,G<:AbstractGraph} <: DoubleGraph{SimpleGraph{ET},ET}
    N::Int
    X0::GraphAF{ET}
    X1::G
    function GraphAddFields(fields::Vector{ET}, g::G) where {ET,G}
        X0 = GraphAF(fields)
        N = X0.N
        getN(g) == N || throw(ArgumentError("incompatible length, fields size=$N graph size=$(getN(g))"))
        return new{ET,G}(N, X0, g)
    end
end

@doc """
    GraphAddFields(fields::Vector, g::AbstractGraph) <: DoubleGraph

A `DoubleGraph` that implements a model in which additional `fields` are added to the
graph `g`. The additional fields can then be taken care of efficiently with [`rrrMC`](@ref).
""" -> GraphAddFields(fields::Vector, g::AbstractGraph)

function update_cache!(X::GraphAddFields, C::Config, move::Int)
    update_cache!(X.X1, C, move)
end

energy(X::GraphAddFields, C::Config) = energy(X.X0, C) + energy(X.X1, C)

delta_energy_residual(X::GraphAddFields, C::Config, move::Int) = delta_energy(X.X1, C, move)

function delta_energy(X::GraphAddFields, C::Config, move::Int)
    return delta_energy(X.X0, C, move) +
           delta_energy_residual(X, C, move)
end

neighbors(X::GraphAddFields, i::Int) = neighbors(X.X1, i)

mutable struct GraphAddSubFields{ET,G<:AbstractGraph} <: DoubleGraph{SimpleGraph{ET},ET}
    N::Int
    X0::GraphAF{ET}
    X1::G
    function GraphAddSubFields(fields::Vector{ET}, g::G) where {ET,G}
        X0 = GraphAF(fields)
        N = X0.N
        getN(g) == N || throw(ArgumentError("incompatible length, fields size=$N graph size=$(getN(g))"))
        return new{ET,G}(N, X0, g)
    end
end

@doc """
    GraphAddSubFields(fields::Vector, g::AbstractGraph) <: DoubleGraph

A `DoubleGraph` that implements a model in which the given `fields` are added and subtracted to the
graph `g`. The fields can then be taken care of efficiently with [`rrrMC`](@ref).
""" -> GraphAddSubFields(fields::Vector, g::AbstractGraph)

function update_cache!(X::GraphAddSubFields, C::Config, move::Int)
    update_cache!(X.X1, C, move)
end

energy(X::GraphAddSubFields, C::Config) = energy(X.X1, C)

delta_energy_residual(X::GraphAddSubFields, C::Config, move::Int) = delta_energy(X.X1, C, move) - delta_energy(X.X0, C, move)

delta_energy(X::GraphAddSubFields, C::Config, move::Int) = delta_energy(X.X1, C, move)

neighbors(X::GraphAddSubFields, i::Int) = neighbors(X.X1, i)

end # module


text:
module testEstimation

using Base.Test
using DataFrames
using TimeData
using Datetime
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running estimator tests\n")

t = Task(() -> AssetMgmt.subsampleProducer(10, 5))
generatedIndices = Array((Int, Int), 5)
counter = 1
for x in t
    generatedIndices[counter] = x
    counter += 1
end

expIndices = [(1,5),
              (1,6),
              (1,7),
              (1,8),
              (1,9)]

@test isequal(generatedIndices, expIndices)

t = Task(() -> AssetMgmt.subsampleProducer(10, 5, rolling=true))
generatedIndices = Array(Any, 5)
counter = 1
for x in t
    generatedIndices[counter] = x
    counter += 1
end

expIndices = [(1,5),
              (2,6),
              (3,7),
              (4,8),
              (5,9)]

@test isequal(generatedIndices, expIndices)

end


text:
using Dispersal, Test, Unitful, DimensionalData
using Dispersal: auxval, cyclic_index, precalc_auxtimeindex, timestep
using DynamicGrids: Extent, SimData
using Unitful: d

suitlayer1 = [0.1 0.2;
              0.3 0.4]

@testset "layer returns the layer value" begin
    @test auxval(suitlayer1, 1, 1, 34325) == 0.1
    @test auxval(suitlayer1, 2, 2, 7685) == 0.4
end

@testset "sequence cycling" begin
    @test cyclic_index(-4, 2) == 2
    @test cyclic_index(-3, 2) == 1
    @test cyclic_index(-2, 2) == 2
    @test cyclic_index(-1, 2) == 1
    @test cyclic_index(0, 2) == 2
    @test cyclic_index(1, 2) == 1
    @test cyclic_index(2, 2) == 2
    @test cyclic_index(3, 2) == 1
    @test cyclic_index(4, 2) == 2
    @test cyclic_index(20, 10) == 10
    @test cyclic_index(21, 10) == 1
    @test cyclic_index(27, 10) == 7
end

a = cat([0.1 0.2; 0.3 0.4], [1.5 1.6; 1.7 1.8]; dims=3)
dimz = X(1:2), Y(1:2), Ti(0d:10d:10d)
suitseq = DimensionalArray(a, dimz)
init = zero(suitseq)

@testset "sequence of layers" begin
    rule = ExponentialGrowthMap(; auxkey=Val(:suit))
    ruleset = Ruleset(rule; timestep=1d)
    data = SimData(Extent(init=init, aux=(suit=suitseq,), tspan=1:1), ruleset)

    @testset "layers sequences are interpolated over timespans" begin
        ind = ((1, 1), (1, 2), (2, 1), (2, 2))
        times = (0d, 20d), (16d, 14d), (19d, 11d), (5d, 45d)#, (15d, 55d)
        for (t1, t2) in times
            interp1 = precalc_auxtimeindex(suitseq, rule, data, t1)
            interp2 = precalc_auxtimeindex(suitseq, rule, data, t2)
            @test interp1 == interp2
        end
    end

    @testset "auxval returns first frame values at 0.5 through the timespan" begin
        t = 5d
        rule1 = precalc_auxtimeindex(suitseq, rule, data, t)
        t = rule1.auxtimeindex
        @test auxval(suitseq, 1, 1, t) == 0.1
        @test auxval(suitseq, 1, 2, t) == 0.2
        @test auxval(suitseq, 2, 1, t) == 0.3
        @test auxval(suitseq, 2, 2, t) == 0.4
    end

    @testset "auxval returns second frame values at 1.5 times through the timespan" begin
        t = 15d
        rule1 = precalc_auxtimeindex(suitseq, rule, data, t)
        t = rule1.auxtimeindex
        @test auxval(suitseq, 1, 1, t) == 1.5
        @test auxval(suitseq, 1, 2, t) == 1.6
        @test auxval(suitseq, 2, 1, t) == 1.7
        @test auxval(suitseq, 2, 2, t) == 1.8
    end

    @testset "auxval returns the second frame one timestep in" begin
        t = 10d
        rule1 = precalc_auxtimeindex(suitseq, rule, data, t)
        t = rule1.auxtimeindex
        @test auxval(suitseq, 1, 1, t) ≈ 1.5
        @test auxval(suitseq, 1, 2, t) ≈ 1.6
        @test auxval(suitseq, 2, 1, t) ≈ 1.7
        @test auxval(suitseq, 2, 2, t) ≈ 1.8
    end

end

@testset "start and stop times" begin
    A = rand(4, 5, 10)
    @test Dispersal.timestep(A) == 1
    @test Dispersal.starttime(A) == 1
    @test Dispersal.stoptime(A) == 10
    A = DimensionalArray(rand(4, 5, 10), (X, Y, Ti(20:10:110)))
    @test Dispersal.timestep(A) == 10
    @test Dispersal.starttime(A) == 20
    @test Dispersal.stoptime(A) == 110
end

@testset "AuxCopy" begin
    init = [0 0]
    l = cat([1 2], [3 4]; dims=3)
    @test AuxCopy(:l, 1) === AuxCopy(auxkey=:l)
    @test AuxCopy{:a,:b}(auxkey=:l, auxtimeindex=2) === AuxCopy{:a,:b}(:l, 2)
    ruleset = Ruleset(AuxCopy(auxkey=:l))
    output = ArrayOutput(init; tspan=1:3, aux=(l=l,))
    sim!(output, ruleset)
    @test output == [[0 0], [3 4], [1 2]]
end


text:
using CSV, DataFrames, ArgParse, Dates

setting = ArgParseSettings()
@add_arg_table setting begin
    "--add"
        help = "add an entry, within double quotes, comma separated as \"name,birthdate,state,relation,email\" with birthdate as yyyy-mm-dd"
    "--latest"
        action = :store_true
        nargs = 0
        help = "print latest (last) entry"
    "--latestfriend"
        action = :store_true
        nargs = 0
        help = "print last friend listed"
    "--latestfamily"
        action = :store_true
        nargs = 0
        help = "print last family member listed"
    "--listbyage"
        action = :store_true
        nargs = 0
        help = "print all ages and entries in birth order"
end

const filename = "example.csv"
const df = CSV.File(filename, dateformat="yyyy-mm-dd") |> DataFrame
const commands = parse_args(setting)
if length(ARGS) == 0
    ArgParse.show_help(setting)
end
const changeflag = [false]

for (k, v) in commands
    if k == "add" && v != nothing
        newrow = Vector{Any}(split(v, r","))
        if length(newrow) == 5 && tryparse(DateTime, newrow[2]) != nothing
            newrow[2] = DateTime(newrow[2])
            push!(df, newrow)
            changeflag[1] = true
            println("Added entry $newrow.")
        end
    elseif k == "latest" && v
        println("The latest entry is $(df[end, :])")
    elseif k == "latestfriend" && v
        println(df[df.Relation .== "friend", :][end, :])
    elseif k == "latestfamily" && v
        println(df[df.Relation .== "family", :][end, :])
    elseif k == "listbyage" && v
        dobcol = df[:Birthdate]
        age = map(x -> round((now() - DateTime(x)).value /(1000*3600*24*365.25), digits=1), dobcol)
        df2 = deepcopy(df)
        df2 = insert!(df, 1, age, :Age)
        println(sort(df2, (:Age)))
    end
end

if changeflag[1]
    CSV.write(filename, df)
    println("Changes written to file $filename.")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
Solitair
Bhupendra Aole
10/20/2020
=#
using Random

#= cards:
    ♠ 1:13
    ♡ 14:26
    ♣ 27:39
    ♢ 40:52
=#
CARD = ["A♠", "2♠", "3♠", "4♠", "5♠", "6♠", "7♠", "8♠", "9♠", "T♠", "J♠", "Q♠", "K♠",
        "A♡", "2♡", "3♡", "4♡", "5♡", "6♡", "7♡", "8♡", "9♡", "T♡", "J♡", "Q♡", "K♡",
        "A♣", "2♣", "3♣", "4♣", "5♣", "6♣", "7♣", "8♣", "9♣", "T♣", "J♣", "Q♣", "K♣",
        "A♢", "2♢", "3♢", "4♢", "5♢", "6♢", "7♢", "8♢", "9♢", "T♢", "J♢", "Q♢", "K♢"]
SUIT = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,]
BACK = ["▮ ¹", "▮ ²", "▮ ³", "▮ ⁴", "▮ ⁵", "▮ ⁶", "▮ ⁷", "▮ ⁸", "▮ ⁹", "▮ ¹⁰",
        "▮ ¹¹", "▮ ¹²", "▮ ¹³", "▮ ¹⁴", "▮ ¹⁵", "▮ ¹⁶", "▮ ¹⁷", "▮ ¹⁸", "▮ ¹⁹", "▮ ²⁰",
        "▮ ²¹", "▮ ²²", "▮ ²³", "▮ ²⁴"]
# card stacks
closed = []
opened = []
completed = [[], [], [], []] # 4 completed stacks
lane = [[], [], [], [], [], [], []] # 7 stacks for closed in lane
working = [[], [], [], [], [], [], []] # 7 working stacks open over lane

function new_game()
    deck = Array(1:52)
    shuffle!(deck)
    for i in 1:7
        for j in 1:i-1
            push!(lane[i], pop!(deck))
        end
        push!(working[i], pop!(deck))
    end
    global closed = deck
end

function display()
    println()
    println("------------------")
    println("X\tO\t\tA\tB\tC\tD")
    # closed stack
    if length(closed)>0
        print(BACK[length(closed)], "\t")
    else
        print("\t")
    end
    # opened stack
    if length(opened)>0
        print(CARD[opened[length(opened)]], "\t\t")
    else
        print("\t\t")
    end
    for c in completed
        if length(c)>0
            print(CARD[c[length(c)]], "\t")
        else
            print("\t")
        end
    end
    println("\n")
    println("1\t2\t3\t4\t5\t6\t7")
    l = deepcopy(lane)
    w = deepcopy(working)

    # first row
    for i = 1:7
        if length(l[i])>0
            print(BACK[length(l[i])], "\t")
        else
            if length(w[i])>0
                wn = popfirst!(w[i])
                print(CARD[wn], "\t")
            else
                print("\t")
            end
        end
    end
    println()

    # reset of the rows
    conti = true
    while conti
        conti = false
        for i = 1:7
            if length(w[i])>0
                wn = popfirst!(w[i])
                print(CARD[wn], "\t")
                conti = true
            else
                print("\t")
            end
        end
        println()
    end
end

function push_all(s, c)
    if isa(c, Array)
        while length(c)>0
            push!(s, popfirst!(c))
        end
    else
        push!(s, c)
    end
end

function place(cards) # card
    # check if it is a single card or an array with a single card
    if isa(cards, Array) && length(cards)==1
        cards = cards[1]
    end

    if isa(cards, Array)
        c = cards[1]
    else # try completed stack if single card is placed
        c = cards
        if c % 13 == 1 # ACE
            for i = 1:4
                if length(completed[i]) == 0
                    push_all(completed[i], c)
                    return true
                end
            end
        else
            for i = 1:4
                if length(completed[i]) > 0 && completed[i][length(completed[i])]+1 == c
                    push_all(completed[i], c)
                    return true
                end
            end
        end
    end

    # check lanes
    for i = 1:7
        if c%13==0 && length(lane[i])==0 && length(working[i])==0 # place king
            push_all(working[i], cards)
            return true
        elseif length(working[i])>0 && working[i][length(working[i])]%13 == (c+1)%13
            # check suit
            t = working[i][length(working[i])]
            if SUIT[t] != SUIT[c]
                push_all(working[i], cards)
                return true
            end
        end
    end

    return false
end

function is_win()
    for i = 1:7
        if length(lane[i]) != 0
            return false
        end
    end
    return true
end

println("Solitair")
new_game()
while true
    display()
    print(">> ")
    reply = lowercase(strip(readline()))

    if length(reply)==0 || reply == "x"
        if length(closed)>0
            push!(opened, pop!(closed))
        else
            global closed = reverse(opened)
            global opened = []
            push!(opened, pop!(closed))
        end
    elseif reply == "end"
        break
    elseif reply == "o"
        if length(opened)>0
            p = opened[length(opened)]
            if place(p)
                pop!(opened)
            end
        end
    elseif length(reply) == 2
        if occursin(reply[1], "1234567")
            if occursin(reply[2], "123456789")
                l = parse(Int32, reply[1])
                w = working[l]
                if length(w)>0
                    n = parse(Int32, reply[2]) # only take these many cards from the top of the stack
                    if place(w[length(w)-(n-1):length(w)])
                        for i = 1:n
                            pop!(w)
                        end
                        if length(lane[l]) > 0 && length(w) == 0
                            push!(w, pop!(lane[l]))
                        end
                    end
                end
            end
        end
    elseif occursin(reply, "1234567")
        l = parse(Int32, reply)
        w = working[l]
        if length(w)>0
            if place(w)
                empty!(w)
                if length(lane[l])>0
                    push!(w, pop!(lane[l]))
                end
            end
        end
    elseif occursin(reply, "abcd")
        w = completed[Int(reply[1])-96]
        if length(w)>0
            if place(w[length(w)])
                pop!(w)
            end
        end
    end

    if is_win()
        display()
        println("You Win!")
        break
    end
end

max_stars_repo_name: aole/Games-in-Julia

max_issues_repo_name: aole/Games-in-Julia

max_forks_repo_name: aole/Games-in-Julia


text:
#####
##### Export dict: SingleFile
#####

export export_dict

function export_dict(D::Dict, outstyle::SingleFile, ft::AbstractFilterTypes=EntireStruct())
  is_leaf(s, prop) = __is_leaf(s, prop, ft)
  if !isempty([triggers(is_leaf(nothing, v)) for v in values(D)])
    open(outstyle.filename,"w") do io
      for (k,v) in D
        if triggers(is_leaf(nothing, v))

          println(io, "$k = $v")
        end
      end
    end
  end
end

#####
##### Export dict: FolderStructure
#####

function export_dict(D::Dict, outstyle::FolderStructure, ft::AbstractFilterTypes=EntireStruct())
  get_path(prop_chain) = joinpath(split(prop_chain, ".")[1:end-1]...)
  get_prop_chain_root(k) = join(split(k, ".")[1:end-1]...)
  get_varname(prop_chain) = split(prop_chain, ".")[end]
  is_leaf(s, prop) = __is_leaf(s, prop, ft)
  mkpath(outstyle.root_folder)
  D_paths = Dict([get_path(k) => [] for k in keys(D)])
  for (k,v) in D
    kp = get_path(k)
    var = get_varname(k)
    push!(D_paths[kp], (var, D[k]))
  end

  for (path,arr) in D_paths
    _fullpath = joinpath(outstyle.root_folder, path)
    folder_needed = [triggers(is_leaf(nothing, val)) for (varname, val) in arr]
    if any(folder_needed)
      mkpath(_fullpath)
      open(joinpath(_fullpath, outstyle.filename), "w") do io
        for (varname, val) in arr
          if triggers(is_leaf(nothing, val))
            println(io, "$varname = $val")
          end
        end
      end
    end
  end
end


text:
# should be your environment and have for nothing to do whith the speed_test
cd("/Users/frank/.julia/dev/TypeDBClient_Speed")
using Pkg
Pkg.activate(".")

using TypeDBClient
using TypeDBClient: CoreSession, CoreClient

using UUIDs
g = TypeDBClient
client = g.CoreClient("127.0.0.1",1729)
Optional{T} = Union{Nothing,T}

# @info "deleting all databases"
# dbs = g.get_all_databases(client)
# for item in dbs
#     g.delete_database(client, item.name)
# end

@info "create database typedb"

if !g.contains_database(client, "typedb")
    g.create_database(client, "typedb")
end
sess = CoreSession(client, "typedb" , g.Proto.Session_Type.DATA, request_timout=Inf)

function trans_func(task)
    return task
end

function coreTransaction(session::g.CoreSession,
    sessionId::g.Bytes,
    type::g.EnumType,
    options::g.TypeDBOptions,
    grpc_controller;
    request_timout::Real=session.request_timeout)

    type = type
    options = options
    input_channel = Channel{g.Proto.Transaction_Client}(10)
    proto_options = g.copy_to_proto(options, g.Proto.Options)

    res_imm = g.Proto.transaction(session.client.core_stub.asyncStub , grpc_controller, input_channel, trans_func)
    req_result, status = fetch(res_imm)

    output_channel = g.grpc_result_or_error(req_result, status, result->result)
    println(status)

    # open_req = g.TransactionRequestBuilder.open_req(session.sessionID, type, proto_options,session.networkLatencyMillis)

    # bidirectionalStream = g.BidirectionalStream(input_channel, output_channel, status)
    # trans_id = uuid4()
    # result = g.CoreTransaction(type, options, bidirectionalStream, trans_id, sessionId, request_timout, session)

    # # The following is only for warming up Transaction. If we didn't do this
    # # it could happen that a transaction reach a timeout.

    # req_result = g.execute(result, open_req, false)
    # kind_of_result = g.Proto.which_oneof(req_result, :res)
    # getproperty(req_result, kind_of_result)

    return input_channel
end

@info "beginning Threads"
client_in = CoreClient("127.0.0.1", 1729)
sess_in = CoreSession(client_in, "typedb", g.Proto.Session_Type.DATA, request_timout=Inf)
grpc_controller = g.gRPCController(request_timeout=sess_in.request_timeout)
Threads.@threads for i in 1:2
    @info "here I am"
    res = coreTransaction(sess_in,
    sess_in.sessionID,
    g.Proto.Transaction_Type.WRITE,
    g.TypeDBOptions(),
    grpc_controller,
    request_timout = sess_in.request_timeout)
    @info "I'm close to end"
end

# close(sess)


text:
mutable struct MBDXMLElement <: MyXMLElement
    el::XMLOrNothing
    precision::MyXMLElement
    precision_prior::MyXMLElement
    is_random::Bool
end

function MBDXMLElement(precision::AbstractArray{Float64, 2},
                       prior_scale::AbstractArray{Float64, 2},
                       is_random::Bool,
                       diagonal_prec::Bool)
    if diagonal_prec
        if !isdiag(prior_scale)
            error("Matrix must be diagonal.")
        else
            precision_param = DiagonalMatrixParameter(diag(prior_scale))
            prior = NothingXMLElement()
        end
    else
        precision_param = MatrixParameter(precision, bn.DIFF_PREC_ID)
        prior = WishartPriorXMLElement(prior_scale, precision_param)
    end
    if !is_random && diagonal_prec
        error("Not implemented") # TODO: gamma prior?
    end
    return MBDXMLElement(nothing, precision_param, prior, is_random)
end

function MBDXMLElement(p::Int; diagonal_prec::Bool = false)
    return MBDXMLElement(Diagonal(ones(p)),
                         Diagonal(ones(p)),
                         true,
                         diagonal_prec)
end

function MBDXMLElement(dm::DiffusionModel)
    P = inv(dm.Σ)
    return MBDXMLElement(P, Diagonal(ones(size(P, 1))), true, false)
end

function make_xml(ml::MBDXMLElement)
    el = new_element(bn.MULTIVARIATE_DIFFUSION_MODEL)
    set_attribute(el, bn.ID, bn.DIFFUSION_ID)
    prec_el = new_element(bn.PRECISION_MATRIX)

    mat_el = make_xml(ml.precision)
    add_child(prec_el, mat_el)

    add_child(el, prec_el)
    ml.el = el

    if ml.is_random
        make_xml(ml.precision_prior)
    end
end

function get_loggables(el::MBDXMLElement)
    if !el.is_random
        return Vector{XMLElement}(undef, 0)
    end
    make_xml(el)
    return [el.el[bn.PRECISION_MATRIX][1][bn.MATRIX_PARAMETER][1]]

end

function get_matrix_parameter(el::MBDXMLElement)
    return get_loggables(el)[1]
end

function get_priors(xml::MBDXMLElement)
    if !xml.is_random
        return Vector{XMLElement}(undef, 0)
    end
    make_xml(xml)
    return [xml.precision_prior.el]
end

function set_precision(mbd::MBDXMLElement, prec::AbstractArray{Float64, 2})
    check_posdef(prec)
    mbd.precision = prec
end

################################################################################
## Priors
################################################################################


text:
using MPI
using StaticArrays
using Logging
using Printf
using LinearAlgebra
using Test
import KernelAbstractions: CPU

using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.ODESolvers
using ClimateMachine.GenericCallbacks
using ClimateMachine.Atmos
using ClimateMachine.VariableTemplates
using ClimateMachine.Thermodynamics

function run_test1(mpicomm, dim, Ne, N, FT, ArrayType)
    warpfun =
        (ξ1, ξ2, ξ3) -> begin
            x1 = ξ1 + (ξ1 - 1 / 2) * cos(2 * π * ξ2 * ξ3) / 4
            x2 = ξ2 + (ξ2 - 1 / 2) * cos(2 * π * ξ2 * ξ3) / 4
            x3 = ξ3 + ξ1 / 4 + sin(2 * π * ξ1) / 16
            return (x1, x2, x3)
        end

    brickrange = (
        range(FT(0); length = Ne + 1, stop = 1),
        range(FT(0); length = Ne + 1, stop = 1),
        range(FT(0); length = 2, stop = 1),
    )
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = (false, false, false),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = warpfun,
    )

    nrealelem = length(topl.realelems)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    vgeo = grid.vgeo
    localvgeo = array_device(vgeo) isa CPU ? vgeo : Array(vgeo)

    S = zeros(Nqk)
    S1 = zeros(Nqk)

    for e in 1:nrealelem
        for k in 1:Nqk
            for j in 1:Nq
                for i in 1:Nq
                    ijk = i + Nq * ((j - 1) + Nq * (k - 1))
                    S[k] +=
                        localvgeo[ijk, grid.x1id, e] *
                        localvgeo[ijk, grid.MHid, e]
                    S1[k] += localvgeo[ijk, grid.MHid, e]
                end
            end
        end
    end

    Stot = zeros(Nqk)
    S1tot = zeros(Nqk)
    Err = 0.0

    for k in 1:Nqk
        Stot[k] = MPI.Reduce(S[k], +, 0, mpicomm)
        S1tot[k] = MPI.Reduce(S1[k], +, 0, mpicomm)
        Err += (0.5 - Stot[k] / S1tot[k])^2
    end
    Err = sqrt(Err / Nqk)
    @test Err < 2e-15
end

function run_test2(mpicomm, dim, Ne, N, FT, ArrayType)
    warpfun =
        (ξ1, ξ2, ξ3) -> begin
            x1 = sin(2 * π * ξ3) / 16 + ξ1
            x2 = ξ2 + (ξ2 - 1 / 2) * cos(2 * π * ξ2 * ξ3) / 4
            x3 = ξ3 + sin(2π * (ξ1)) / 20
            return (x1, x2, x3)
        end

    brickrange = (
        range(FT(0); length = Ne + 1, stop = 1),
        range(FT(0); length = Ne + 1, stop = 1),
        range(FT(0); length = 2, stop = 1),
    )
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = (false, false, false),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = warpfun,
    )

    nrealelem = length(topl.realelems)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    vgeo = grid.vgeo
    localvgeo = array_device(vgeo) isa CPU ? vgeo : Array(vgeo)

    S = zeros(Nqk)
    S1 = zeros(Nqk)
    K = zeros(Nqk)

    for e in 1:nrealelem
        for k in 1:Nqk
            for j in 1:Nq
                for i in 1:Nq
                    ijk = i + Nq * ((j - 1) + Nq * (k - 1))
                    S[k] +=
                        localvgeo[ijk, grid.x1id, e] *
                        localvgeo[ijk, grid.MHid, e]
                    S1[k] += localvgeo[ijk, grid.MHid, e]
                    K[k] = localvgeo[ijk, grid.x3id, e]
                end
            end
        end
    end

    Stot = zeros(Nqk)
    S1tot = zeros(Nqk)
    Err = 0.0

    for k in 1:Nqk
        Stot[k] = MPI.Reduce(S[k], +, 0, mpicomm)
        S1tot[k] = MPI.Reduce(S1[k], +, 0, mpicomm)
        Err += (0.5 + sin(2 * π * K[k]) / 16 - Stot[k] / S1tot[k])^2
    end
    Err = sqrt(Err / Nqk)
    @test 2e-15 > Err
end

function run_test3(mpicomm, dim, Ne, N, FT, ArrayType)
    base_Nhorz = 4
    base_Nvert = 2
    Rinner = 1 // 2
    Router = 1
    expected_result = [
        -4.5894269717905445e-8 -1.1473566985387151e-8
        -2.0621904184281448e-10 -5.155431637149377e-11
        -8.72191208145523e-13 -2.1715962361668062e-13
    ]
    for l in 1:3
        Nhorz = 2^(l - 1) * base_Nhorz
        Nvert = 2^(l - 1) * base_Nvert
        Rrange = grid1d(FT(Rinner), FT(Router); nelem = Nvert)
        topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
        grid = DiscontinuousSpectralElementGrid(
            topl,
            FloatType = FT,
            DeviceArray = ArrayType,
            polynomialorder = N,
            meshwarp = Topologies.cubedshellwarp,
        )

        nrealelem = length(topl.realelems)
        Nq = N + 1
        Nqk = dimensionality(grid) == 2 ? 1 : Nq
        vgeo = grid.vgeo
        localvgeo = array_device(vgeo) isa CPU ? vgeo : Array(vgeo)

        topology = grid.topology
        nvertelem = topology.stacksize
        nhorzelem = div(nrealelem, nvertelem)

        Surfout = 0
        Surfin = 0

        for ev in 1:nvertelem
            for eh in 1:nhorzelem
                e = ev + (eh - 1) * nvertelem
                for i in 1:Nq
                    for j in 1:Nq
                        for k in 1:Nqk
                            ijk = i + Nq * ((j - 1) + Nqk * (k - 1))
                            if (k == Nqk && ev == nvertelem)
                                Surfout += localvgeo[ijk, grid.MHid, e]
                            end
                            if (k == 1 && ev == 1)
                                Surfin += localvgeo[ijk, grid.MHid, e]
                            end
                        end
                    end
                end
            end
        end
        Surfouttot = MPI.Reduce(Surfout, +, 0, MPI.COMM_WORLD)
        Surfintot = MPI.Reduce(Surfin, +, 0, MPI.COMM_WORLD)
        @test (4 * π * Router^2 - Surfouttot) ≈ expected_result[l, 1] rtol =
            1e-3 atol = eps(FT) * 4 * π * Router^2
        @test (4 * π * Rinner^2 - Surfintot) ≈ expected_result[l, 2] rtol = 1e-3 atol =
            eps(FT) * 4 * π * Rinner^2
    end
end

# Test for 2D integral
function run_test4(mpicomm, dim, Ne, N, FT, ArrayType)
    brickrange = ntuple(j -> range(FT(0); length = Ne + 1, stop = 1), 2)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, 2),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    nrealelem = length(topl.realelems)
    Nq = N + 1
    vgeo = grid.vgeo
    localvgeo = array_device(vgeo) isa CPU ? vgeo : Array(vgeo)

    S = zeros(Nq)

    for e in 1:nrealelem
        for i in 1:Nq
            for j in 1:Nq
                ij = i + Nq * (j - 1)
                S[j] +=
                    localvgeo[ij, grid.x1id, e] * localvgeo[ij, grid.MHid, e]
            end
        end
    end

    Err = 0
    for j in 1:Nq
        Err += (S[j] - 0.5)^2
    end

    Err = sqrt(Err / Nq)
    @test Err <= 1e-15
end

function run_test5(mpicomm, dim, Ne, N, FT, ArrayType)
    warpfun = (ξ1, ξ2, ξ3) -> begin
        x1 = cos(π * ξ2) / 16 + abs(ξ1)
        x2 = ξ2
        x3 = ξ3
        return (x1, x2, x3)
    end

    brickrange = ntuple(j -> range(FT(0); length = Ne + 1, stop = 1), 2)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, 2),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = warpfun,
    )

    nrealelem = length(topl.realelems)
    Nq = N + 1
    vgeo = grid.vgeo
    localvgeo = array_device(vgeo) isa CPU ? vgeo : Array(vgeo)

    S = zeros(Nq)
    J = zeros(Nq)

    for e in 1:nrealelem
        for i in 1:Nq
            for j in 1:Nq
                ij = i + Nq * (j - 1)
                S[j] +=
                    localvgeo[ij, grid.x1id, e] * localvgeo[ij, grid.MHid, e]
                J[j] = localvgeo[ij, grid.x2id, e]
            end
        end
    end

    Err = 0
    for j in 1:Nq
        Err += (S[j] - 0.5 - cos(π * J[j]) / 16)^2
    end

    Err = sqrt(Err / Nq)
    @test Err < 1e-15
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    FT = Float64
    dim = 3
    Ne = 1
    polynomialorder = 4

    @info (ArrayType, FT, dim)

    @testset "horizontal_integral" begin
        run_test1(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
        run_test2(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
        run_test3(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
        run_test4(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
        run_test5(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
    end
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# Create summary for Summary_value (default constructor by
# protobuf is broken).
SummaryCollection(;kwargs...) = Summary(value=Vector{Summary_Value}(); kwargs...)
SummaryCollection(summaries::Vector{Summary_Value}; kwargs...) = Summary(value=summaries; kwargs...)
SummaryCollection(summary::Summary_Value; kwargs...) = Summary(value=[summary]; kwargs...)
SummaryCollection(summary::GraphDef; kwargs...) = summary

# TODO Clean up this nothing crazyness
function make_event(logger::TBLogger, summary::Summary; step=TensorBoardLogger.step(logger))
    step = typeof(step) == Nothing ?  TensorBoardLogger.step(logger) : step
    return Event(wall_time=time(), summary=summary, step=step)
end
function make_event(logger::TBLogger, summary::GraphDef; step=TensorBoardLogger.step(logger))
    step = typeof(step) == Nothing ?  TensorBoardLogger.step(logger) : step
    Event(wall_time=time(), graph_def=serialize_proto(summary), step=step)
end

"""
    write_event(out::IOStream, event::Event)

Serializes the Event `event` to the `out` stream according to the TensorBoard
format. The format follows the following rule (in bytes)

#1 0...8       Int64   size of serialized `event`
#2 9...15      UInt32  masked_CRC of #1
#3 16...N        -     serialized `event` as protobuffer
#4 N..N+8      UInt32  masked_CRC of #3
"""
function write_event(out::IOStream, event::Event)
    data = PipeBuffer();
    writeproto(data, event)

    #header
    header     = collect(reinterpret(UInt8, [data.size]))
    crc_header = reinterpret(UInt8, UInt32[masked_crc32c(header)])
    crc_data   = reinterpret(UInt8, UInt32[masked_crc32c(data.data)])

    write(out, header)
    write(out, crc_header)
    write(out, data.data)
    write(out, crc_data)
    flush(out)
end

write_event(logger::TBLogger, event::Event) = write_event(logger.file, event)


text:
#
#   The joint diagonalization problem on the Stiefel manifold St(n,k)
#
using Manopt, Manifolds, ManifoldsBase, LinearAlgebra, Random
import Manifolds: vector_transport_to!
struct IdentityTransport <: AbstractVectorTransportMethod end
vector_transport_to!(::Stiefel, Y, p, X, q, ::IdentityTransport) = (Y .= project(M, q, X))
Random.seed!(42)
n = 12
k = 8
m = 512
A = randn(n, n, m)

for i in 1:m
    A[:, :, i] = diagm(n:-1:1) + 0.1 * (transpose(A[:, :, i]) + A[:, :, i])
end

M = Stiefel(n, k)
F(X::Array{Float64,2}) = -sum([norm(diag(X' * A[:, :, i] * X))^2 for i in 1:m])
function ∇F(X::Array{Float64,2})
    return project(
        M, X, -4 * sum([A[:, :, i] * X * norm(diag(X' * A[:, :, i] * X)) for i in 1:m])
    )
end
x = random_point(M)
@time quasi_Newton(
    M,
    F,
    ∇F,
    x;
    memory_size=32,
    cautious_update=true,
    vector_transport_method=IdentityTransport(),
    stopping_criterion=StopWhenGradientNormLess(norm(M, x, ∇F(x)) * 10^(-6)),
    debug=[:Iteration, " ", :Cost, "\n", 1, :Stop],
)


text:
using LandscapeChange
using NeutralLandscapes
using DynamicGrids
using Rasters
using Plots
using Test
using StaticArrays

@testset "" begin
    proj = EPSG(4326)
    projdims = X(Projected(20:0.01:22; crs=proj)),
               Y(Projected(9:0.01:11; crs=proj))
    random_landscape = Raster(rand(MidpointDisplacement(), projdims); missingval=nothing)
    c = NamedVector((forest=1, cleared=2, settled=3))

    landscape = Rasters.classify(random_landscape,
         0..0.5 => c.forest, 0.5..0.8 => c.cleared, 0.8..1.0 => c.settled;
         others=typemin(Int), missingval=typemin(Int)
    )
    # plot(plot(random_landscape), plot(landscape; c=:viridis))

    counts = broadcast(_ -> zero(c), landscape)

    lu_count_rule = LandCoverCount{:landscape,:counts}(Window(1), c)
    extent = Extent((; landscape, counts); tspan=1:10)
    sd = DynamicGrids.SimData(extent, lu_count_rule)
    sd = step!(sd)
    sd[:counts]

    suitability_params = (
        riverdist=Aux{:distance_to_rivers}(),
        coastdist=Aux{:distance_to_coast}(),
        elevation=Aux{:elevation}(),
        steepness=Aux{:steepness}(),
    )
    suitability_params = (;)

    suitability_functions = (
        forest=lc -> (forest=0.9, cleared=0.1, settled=0.0)[lc],
        cleared=lc -> (forest=0.1, cleared=0.9, settled=0.01)[lc],
        settlement=lc -> (forest=0.0, cleared=0.05, settled=0.95)[lc],
    )

    lu_suitability_rule = LandCoverSuitability{:landscape,:suitability}(;
        functions=suitability_functions, mode=ChangeProbabilities()
    )

    using BenchmarkTools
    f(xs) = map(x -> x^2, xs)
    @btime f(distances(Window(1)))

    suitability = broadcast(_ -> zero(c) .* 0.0, landscape)
    extent = Extent((; landscape, counts, suitability); tspan=1:10)
    sd = DynamicGrids.SimData(extent, lu_suitability_rule)
    sd = step!(sd)
    sd[:suitability]
end

@testset "NamedVector" begin
    using BenchmarkTools
    NamedVector{(:a,:b,:c),3,Int}(NamedTuple{(:a,:b,:c)}((1,2,3)))
    @btime NamedVector{(:a,:b,:c),3}((1.0,2,3))
    @btime NamedVector{(:a,:b,:c)}((1,2,3))
    # This doesn't completely compile away for some reason
    @btime NamedVector{(:a,:b,:c),3}(NamedTuple{(:a,:b,:c)}((1,2,3)))
    @code_warntype NamedVector(NamedTuple{(:a,:b,:c)}((1,2,3)))
    @btime map(x -> 2x, $c)
end


text:
export imageMode, loadImage, noTint, save

## Loading & Displaying

function imageMode(mode)
	state.imageMode = mode
end

function loadImage(filename::AbstractString)
    image = Images.imread(filename)
end

function noTint()
    state.tintStuff = false
end

function save(filename::AbstractString)
	data = Array{GLfloat}(3, state.width*2, state.height*2)
	glReadPixels(0, 0, state.width*2, state.height*2, GL_RGB, GL_FLOAT, data)
	data[1,:,:] = flipud(squeeze(data[1,:,:],1))
	data[2,:,:] = flipud(squeeze(data[2,:,:],1))
	data[3,:,:] = flipud(squeeze(data[3,:,:],1))
	img = Image(data, colorspace = "RGB", spatialorder = ["x", "y"], colordim = 1)
	Images.imwrite(img, filename)
end

# function tint()
#
# end


text:
# Find the Most Frequent Words with Mismatches in a String

const NUCLEOTIDES = ('A', 'C', 'G', 'T')
const NUCLEOTIDE_ORDER = Dict{Char, Int64}('A' => 0, 'C' => 1, 'G' => 2, 'T' => 3)

convert_to_intlist(line) = map(x -> parse(Int, x), split(line))

pattern_to_number(pattern) = foldl((acc, nucleotide) -> 4 * acc + NUCLEOTIDE_ORDER[nucleotide], pattern; init=0)

function number_to_pattern(encoding, k)
    pattern = Char[]
    for _ in 1:k
        encoding, remainder = divrem(encoding, 4)
        push!(pattern, NUCLEOTIDES[remainder+1])
    end
    join(reverse(pattern))
end

function calc_frequencies_for_approximate_occurrences(text, k, d)
    frequencies = zeros(Int32, 4 ^ k)
    for ix in 1:(length(text)-k+1)
        exact_pattern = text[ix:ix+k-1]
        for pattern in generate_neighbourhood(exact_pattern, d)
            code = pattern_to_number(pattern)
            frequencies[code+1] += 1
        end
    end
    frequencies
end

function calc_immediate_neighbours(pattern)
    neighbours = Set{String}()
    for (ix, letter) in enumerate(pattern)
        for nucleotide in NUCLEOTIDES
            if nucleotide != letter
                pattern_array = collect(pattern)
                pattern_array[ix] = nucleotide
                push!(neighbours, join(pattern_array))
            end
        end
    end
    neighbours
end

function generate_neighbourhood(text, d)
    neighbourhood = Set{String}([text])
    for _ in 1:d
        neighbours_of_neighbours = Set{String}()
        for pattern in neighbourhood
            immediate_neighbours = calc_immediate_neighbours(pattern)
            union!(neighbours_of_neighbours, immediate_neighbours)
        end
        union!(neighbourhood, neighbours_of_neighbours)
    end
    neighbourhood
end

function most_frequent_approximate_k_mers(text, k, d)
    frequencies = calc_frequencies_for_approximate_occurrences(text, k, d)
    max_occurrence = maximum(frequencies)
    patterns = String[]
    for (ix, freq) in enumerate(frequencies)
        if freq == max_occurrence
            pattern = number_to_pattern(ix - 1, k)
            push!(patterns, pattern)
        end
    end
    patterns
end

function main()
    text = readline()
    k, d = convert_to_intlist(readline())
    result = most_frequent_approximate_k_mers(text, k, d)
    println(join(result, " "))
end

main()


text:
meshgrid(x, y) = (repeat(x, outer = length(y)), repeat(y, inner = length(x)))

@userplot ScalarHeatMap
@recipe function f(p::ScalarHeatMap)
    field, = p.args
    @series begin
        seriestype := :heatmap
        aspect_ratio := :equal
        transpose(field)
    end
end

@userplot VectorFieldQuiver
@recipe function f(p::VectorFieldQuiver)
    x, y, vectorfield = p.args
    n_arrows = 20
    skip_x = ceil(Int, length(x) / n_arrows)
    skip_y = ceil(Int, length(y) / n_arrows)
    vx_mesh = vectorfield[1][1:skip_x:end, 1:skip_y:end][:]
    vy_mesh = vectorfield[2][1:skip_x:end, 1:skip_y:end][:]
    x_mesh, y_mesh = meshgrid(x[1:skip_x:end], y[1:skip_y:end])
    Lx = x[end] - x[1]
    max_v = maximum(sqrt.(vx_mesh .^ 2 .+ vy_mesh .^ 2))
    scale = 0.2 * Lx / max_v
    headlength = 0.3
    headwidth = 0.2
    @series begin
        seriestype := :quiver
        aspect_ratio := :equal
        xlims := (x[1], x[end])
        ylims := (y[1], y[end])
        arrow := arrow(:closed, :head, headlength, headwidth)
        quiver := (vx_mesh * scale, vy_mesh * scale)
        x_mesh, y_mesh
    end
end

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
using BenchmarkTools
# using BenchmarkTools: @benchmarkable, BenchmarkGroup

using GenomicFeatures

include(joinpath(@__DIR__, "..", "test", "Utilities.jl"))

import ..Utilities: random_intervals

N = 1000
SEED = 1234
SEQNAMES = "chr" .* string.(1:3)

intervals = random_intervals(SEQNAMES, 1000, N, SEED)
intervals_sorted = sort(intervals)

SUITE = BenchmarkGroup()

let suite = SUITE["accessors"] = BenchmarkGroup()
    s0 = suite["$(typeof(intervals))"] = BenchmarkGroup()
    s0["seqname"] = @benchmarkable(seqname.($intervals))
    s0["leftposition"] = @benchmarkable(leftposition.($intervals))
    s0["rightposition"] = @benchmarkable(rightposition.($intervals))
    s0["strand"] = @benchmarkable(strand.($intervals))
    s0["metadata"] = @benchmarkable(metadata.($intervals))
end

let suite = SUITE["sort"] = BenchmarkGroup()
    suite["$(typeof(intervals))"] = @benchmarkable(sort(i), setup=(i = copy($intervals)))
end

let suite = SUITE["insert"] = BenchmarkGroup()
    suite["shorthand"] = @benchmarkable(IntervalCollection($intervals_sorted))
    suite["type"] = @benchmarkable(IntervalCollection{Int}($intervals_sorted))
end

let suite = SUITE["push"] = BenchmarkGroup()
    suite["$(typeof(intervals))"] = @benchmarkable([push!(col, i) for i in $intervals], setup=(col=IntervalCollection{Int}()))
end

let suite = SUITE["eachoverlap"] = BenchmarkGroup()
    intervals_a = intervals_sorted
    intervals_b = sort(random_intervals(SEQNAMES, 1000, N, SEED+1))

    col_a = IntervalCollection(intervals_a)
    col_b = IntervalCollection(intervals_b)

    As = [intervals_a, col_a]
    Bs = [intervals_b, col_b]

    for (A, B) in Iterators.product(As,Bs)
        str = "$(typeof(A)), $(typeof(B))"
        suite[str] = @benchmarkable(collect(eachoverlap($A,$B)))
    end
end


text:
@symbol_func function P_n(cur_reactor::AbstractReactor)
  cur_P = P_F(cur_reactor)

  cur_P *= (4/5)

  cur_P
end


text:
# smoothing polygons

Using Luxor

@png begin
    background("grey15")
   tiles = Tiler(600, 250, 1, 5, margin=10)
    for (pos, n) in tiles
        p = star(pos, tiles.tilewidth/2 - 2, 5, 0.3, 0, vertices=true)
        setdash("dot")
        sethue("white")
        prettypoly(p, close=true, :stroke)
        setdash("solid")
        sethue("antiquewhite")
        polysmooth(p, n * 2, :fill)
    end
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using Knet
using Knet: stack_isempty
include("readdata.jl")

# Trying an rnn model.

@knet function model05(x; hidden=0, o...)
    h = lstm(x; o..., out=hidden)
    if predict
        return wbf(h; o..., out=2, f=:soft)
    end
end

function rnnforw(f, x)
    # x is a (1,T,D,N) one-hot tensor
    # T:sequence length, D:input dimension, N:minibatch size
    global xdbg = x
    (X,T,D,N) = size(x)
    reset!(f)
    ypred = nothing
    for t=1:T
        xt = reshape(x[1,t,:,:], (D, N))
        ypred = sforw(f, xt; predict=(t==T))
    end
    return ypred
end

function train05(; N=2^10, lr=0.001, adam=true, nbatch=128, gclip=0, hidden=256, seed=0, o...)
    seed > 0 && setseed(seed)
    global f = compile(:model05; hidden=hidden, o...)
    setp(f, lr=lr, adam=adam)
    global d = GulayseData(nbatch)
    sloss = zloss = -1; n = nextn = 1
    for (x,ygold) in d
        ypred = rnnforw(f, x)
        sback(f, ygold, softloss)
        while !stack_isempty(f); sback(f); end
        gscale = (gclip > 0 ? min(1.0, gclip/gnorm(f)) : 1.0)
        update!(f; gscale=gscale)
        sl = softloss(ypred,ygold); sloss = (n==1 ? sl : 0.99 * sloss + 0.01 * sl)
        zl = zeroone(ypred,ygold);  zloss = (n==1 ? zl : 0.99 * zloss + 0.01 * zl)
        n==nextn && (println((n,sloss,1-zloss)); nextn*=2)
        (n += 1) > N && break
    end
end

train05()

max_stars_repo_name: JuliaTagBot/BindingPrediction

max_issues_repo_name: JuliaTagBot/BindingPrediction

max_forks_repo_name: JuliaTagBot/BindingPrediction


text:
function FiniteElements(::DomainStyle,model::DiscreteModel,basis::ReferenceFEName,args...;kwargs...)
  @abstractmethod "The factory function FiniteElements has not been defined for the given arguments"
end

function FiniteElements(
  ::ReferenceDomain,
  model::DiscreteModel,
  basis::ReferenceFEName,
  args...;
  kwargs...)

  cell_reffe = ReferenceFE(model,basis,args...;kwargs...)
  cell_map = get_cell_map(Triangulation(model))
  CellFE(cell_map,cell_reffe)
end

function FiniteElements(
  ::PhysicalDomain,
  model::DiscreteModel,
  basis::Lagrangian,
  args...;
  kwargs...)

  # Reference FEs and cell_map
  cell_reffe = ReferenceFE(model,basis,args...;kwargs...)
  cell_map = get_cell_map(Triangulation(model))
  ctype_reffe, cell_ctype = compress_cell_data(cell_reffe)

  # prebasis # TODO we need a better conditioned cell_prebasis
  ctype_prebasis = map(get_prebasis,ctype_reffe)
  cell_prebasis = expand_cell_data(ctype_prebasis,cell_ctype)

  # dof basis
  ctype_ref_dof_basis = map(get_dof_basis,ctype_reffe)
  ctype_ref_nodes = map(i->i.nodes,ctype_ref_dof_basis)
  cell_ref_nodes = expand_cell_data(ctype_ref_nodes,cell_ctype)
  cell_node_ids = lazy_map(evaluate,cell_map,cell_ref_nodes)
  cell_ref_dof_basis = expand_cell_data(ctype_ref_dof_basis,cell_ctype)
  cell_dof_basis = lazy_map(LagrangianDofBasis,cell_ref_dof_basis,cell_node_ids)

  # Shape functions
  cell_dof_values = lazy_map(evaluate,cell_dof_basis,cell_prebasis)
  cell_change = lazy_map(inv,cell_dof_values)
  cell_shapefuns = lazy_map(linear_combination,cell_change,cell_prebasis)

  # Build the CellFE
  ctype_num_dofs = map(num_dofs,ctype_reffe)
  ctype_ldof_comp = map(reffe->get_dof_to_comp(reffe),ctype_reffe)
  cell_conformity = CellConformity(cell_reffe)
  cell_shapefuns_domain = PhysicalDomain()
  cell_dof_basis_domain = cell_shapefuns_domain
  max_order = maximum(map(get_order,ctype_reffe))

  CellFE(
    cell_ctype,
    ctype_num_dofs,
    ctype_ldof_comp,
    cell_conformity,
    cell_shapefuns,
    cell_dof_basis,
    cell_shapefuns_domain,
    cell_dof_basis_domain,
    max_order)
end


text:
function f(x::Real, σ::Real, μ::Real)
    return exp(-0.5((x-μ)/σ)^2) / (σ * sqrt(2π))
end

f(σ::Real, μ::Real) = x -> f(x, σ, μ)

surface(range::AbstractRange; σ=1, μ=0) = f(σ, μ).(range) * f(σ, μ).(range)'

max_stars_repo_name: foldfelis/Replotter

max_issues_repo_name: foldfelis/Replotter

max_forks_repo_name: foldfelis/Replotter


text:
@testset "Condition StraightLineHomotopy" begin
    @polyvar x
    f = [x^2-1]
    g = [x^2+1]

    z=[1.0]
    t=1.0

    H = StraightLineHomotopy(f,g)
    @test κ(H,z,t) ≈ 1/sqrt(2)
    @test κ_norm(H,z,t) ≈ 1.0
    @test μ_norm(H,z,t) ≈ 1.0
end

@testset "Condition GeodesicOnTheSphere" begin
    @polyvar x
    f = [x^2-1]
    g = [x^2+1]

    z=[1.0, 1.0]
    t=1.0

    H = GeodesicOnTheSphere(f,g)
    @test κ(H,z,t) ≈ 1/sqrt(2)
    @test κ_norm(H,z,t) ≈ 1.0
    @test μ_norm(H,z,t) ≈ 1.0
end


text:
import Base: show
import Oceananigans.Grids: short_show
using Oceananigans.Grids: show_domain

location_str(::Face) = "Face"
location_str(::Cell) = "Cell"
show_location(X, Y, Z) = "($(location_str(X())), $(location_str(Y())), $(location_str(Z())))"
show_location(field::AbstractField{X, Y, Z}) where {X, Y, Z} = show_location(X, Y, Z)

short_show(field::Field) = string("Field located at ", show_location(field))

show(io::IO, field::Field) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))\n",
          "└── boundary conditions: $(short_show(field.boundary_conditions))")

short_show(array::OffsetArray{T, D, A}) where {T, D, A} = string("OffsetArray{$T, $D, $A}")


text:
include("init.jl")
include("utils.jl")
include("wrapper.jl")
include("api.jl")


text:
module CRF
export
    Sequence,
    Features,
    @append!,
    empty!,
    label,
    loglikelihood,
    loglikelihood_gradient,
    update

include("base.jl")
include("features.jl")
include("util.jl")
end