text:
using CorpusLoaders
using MultiResolutionIterators
using TextAnalysis
using TextAnalysis: CRF, crf_loss
using WordTokenizers
using Embeddings
using Flux
using Flux: onehot, param, onehotbatch, Conv, LSTM, flip, Momentum, reset!, onecold, batchseq
using Tracker
using BSON: @save, @load
using LinearAlgebra
using Base.Iterators: partition
using CorpusLoaders: NERTaggedWord

println("Dependencies Loaded")

device = :gpu

if device == :gpu
    using CuArrays
end

CHAR_EMBED_DIMS = 25
WORD_EMBED_DIMS = 100

train_set = load(CoNLL(), "train") # training set
# test_set = load(CoNLL(), "test") # test set
dev_set = load(CoNLL(), "dev") # dev set

dataset = flatten_levels(train_set, lvls(CoNLL, :document)) |> full_consolidate
dev_dataset = flatten_levels(dev_set, lvls(CoNLL, :document)) |> full_consolidate

get_pos(w::NERTaggedWord) = w.pos

X_train = [CorpusLoaders.word.(sent) for sent in dataset]
Y_train = [get_pos.(sent) for sent in dataset]
@assert length.(X_train) == length.(Y_train)

X_dev = [CorpusLoaders.word.(sent) for sent in dev_dataset]
Y_dev = [get_pos.(sent) for sent in dev_dataset]
@assert length.(X_dev) == length.(Y_dev)

words_vocab = unique(vcat(X_train...))
alphabets = unique(vcat(collect.(words_vocab)...))
labels = unique(vcat(Y_train...))
DENSE_OUT_SIZE = num_labels = length(labels)

#### Word Embeddings
# using these indices to generate the flux one hot input word embeddings vectors.

embtable = load_embeddings(GloVe, 2)
get_word_index = Dict(word => ii for (ii, word) in enumerate(embtable.vocab))
get_word_from_index = Dict(value => key for (key, value) in get_word_index)
W_word_Embed = (embtable.embeddings)

# One Vec for unknown chars
UNK_Word = "<UNK>"
UNK_Word_Idx = length(get_word_index) + 1
embedding_vocab_length = length(get_word_index) + 1
@assert UNK_Word ∉ collect(keys(get_word_index))
get_word_index[UNK_Word] = UNK_Word_Idx
get_word_from_index[UNK_Word_Idx] = UNK_Word
W_word_Embed = hcat(W_word_Embed, rand(WORD_EMBED_DIMS))

# One Vec for PAD chars
PAD_Word = "<pad>"
PAD_Word_Idx = length(get_word_index) + 1
embedding_vocab_length = length(get_word_index) + 1
@assert PAD_Word ∉ collect(keys(get_word_index))
get_word_index[PAD_Word] = PAD_Word_Idx
get_word_from_index[PAD_Word_Idx] = PAD_Word
W_word_Embed = hcat(W_word_Embed, rand(WORD_EMBED_DIMS))

W_word_Embed = Float32.(W_word_Embed)
@assert size(W_word_Embed, 2) == embedding_vocab_length

##### Char Embeddings
UNK_char = '¿'
PAD_char = '█'
@assert UNK_char ∉ alphabets
@assert PAD_char ∉ alphabets
push!(alphabets, UNK_char)
push!(alphabets, PAD_char)
sort!(unique!(alphabets))

# using these indices to generate the flux one hot input Char embeddings vectors.

W_Char_Embed = rand(CHAR_EMBED_DIMS, length(alphabets))
W_Char_Embed = Float32.(W_Char_Embed)
W_Char_Embed = param(W_Char_Embed)
get_char_index = Dict(char => ii for (ii, char) in enumerate(alphabets))
get_char_from_index = Dict(value => key for (key, value) in get_char_index)

println("Embeddings Made")

############# Creating onehoot vectors (useful for embeddings lookup), convenient
BATCH_SIZE = 10
X_train = sort(X_train, by=length, alg=MergeSort)
Y_train = sort(Y_train, by=length, alg=MergeSort) # For stable sorting

onehotword(word) = onehot(get(get_word_index, lowercase(word), get_word_index[UNK_Word]), 1:embedding_vocab_length)
onehotbatchword(words) = onehotbatch([get(get_word_index, lowercase(word), get_word_index[UNK_Word]) for word in words], 1:embedding_vocab_length)
onehotchars(word) = onehotbatch([get(get_char_index, c, get_char_index[UNK_char]) for c in word], 1:length(alphabets))
onehotlabel(label) = onehot(label, labels)
onehotinput(word) = (onehot(get(get_word_index, lowercase(word), get_word_index[UNK_Word]), 1:embedding_vocab_length),
                onehotbatch([get(get_char_index, c, get_char_index[UNK_char]) for c in word], 1:length(alphabets)))

# align(k) = [[j[i]  for j in k] for i in 1:length(k[1])]
# batch_char_input(ohed_chars, p) = [align(padded(b, p)) for b in collect(partition(ohed_chars, BATCH_SIZE))]
# batches(xs, p) = [batchseq(b, p) for b in partition(xs, BATCH_SIZE)]

oh_seq(arr, f) = [f(element) for element in arr]

padded(b, p, m) = [length(elem) < m ? [elem..., repeat([p], m - length(elem))...] : elem for elem in b]
padded(b, p) = padded(b, p, maximum(length.(b)))

align_chars_new(b, p, m) = [[[p, p, elem[i]..., repeat([p], m[i] - length(elem[i]))..., p, p] for i in 1:length(elem)] for elem in b]
align_chars_new(b, p) = align_chars_new(b, p, [maximum([length(j[i]) for j in b])   for i in 1:length(b[1])])

batch_oh_chars(o) = [[char_w_seq[i] for char_w_seq in o] for i in 1:length(o[1])]
gpuit(x) = cu.(x)
batch_new_chars(xs, p) = [gpuit.(batch_oh_chars(oh_seq.(align_chars_new(padded(b, p), '█'), onehotchars))) for b in partition(xs, BATCH_SIZE)]
# X_chars_train = batch_char_input([oh_seq(sentence, onehotchars) for sentence in X_train], onehotchars("█")) # A Bunch Sequences of Array of Chars, done to prevent repeated computations.
X_newchars_train = batch_new_chars([collect.(sentence) for sentence in X_train], collect("█"))

padded_word_batches(b, p, m) = [[length(elem) < i ? p : elem[i] for elem in b] for i in 1:m]
padded_word_batches(b, p) = padded_word_batches(b, p, maximum(length.(b)))
X_words_train = [cu.(oh_seq(padded_word_batches(b, PAD_Word), onehotbatchword)) for b in partition(X_train, BATCH_SIZE)] # A Bunch of Sequences of words, i.e. sentences

X_input_train = [[(k,l) for (k,l) in zip(i,j)] for (i,j) in zip(X_words_train, X_newchars_train)]
Y_oh_train = [cu.(oh_seq.(b, onehotlabel)) for b in partition(Y_train, BATCH_SIZE)]
deleteat!(Y_oh_train, length(Y_oh_train))
deleteat!(X_input_train, length(X_input_train))

X_dev = sort(X_dev, by=length, alg=MergeSort)
Y_dev = sort(Y_dev, by=length, alg=MergeSort) # For stable sorting

X_newchars_dev = batch_new_chars([collect.(sentence) for sentence in X_dev], collect("█"))
X_words_dev = [cu.(oh_seq(padded_word_batches(b, PAD_Word), onehotbatchword)) for b in partition(X_dev, BATCH_SIZE)] # A Bunch of Sequences of words, i.e. sentences
X_input_dev = [[(k,l) for (k,l) in zip(i,j)] for (i,j) in zip(X_words_dev, X_newchars_dev)]
Y_oh_dev = [oh_seq.(b, onehotlabel) for b in partition(Y_dev, BATCH_SIZE)]
# deleteat!(Y_oh_dev, length(Y_oh_dev))
# deleteat!(X_input_dev, length(X_input_dev))

# if device == :gpu
#     X_input_train = [cu.(oh_seq(sentence, onehotinput)) for sentence in X_train]
# else
#     X_input_train = [oh_seq(sentence, onehotinput) for sentence in X_train]
# end
# Y_oh_train = [oh_seq(tags_sequence, onehotlabel) for tags_sequence in Y_train]
#
# if device == :gpu
#     Y_oh_train = cu.(Y_oh_train)
# end
#
# X_input_dev = [oh_seq(sentence, onehotinput) for sentence in X_dev]
# Y_oh_dev = [oh_seq(tags_sequence, onehotlabel) for tags_sequence in Y_dev]
#
# @assert length.(X_train) == length.(X_input_train) == length.(Y_oh_train)

#################### MODEL ######################
# 1. Character Embeddings with CNNs
# 2. Obtain Word Embeddings, concatenate Char & Word Embeddings
# 3. Bi-LSTM
# 4. Softmax or CRF on each unit in sequence

CNN_OUTPUT_SIZE = 30
CONV_WINDOW_LENGTH = 3
LSTM_STATE_SIZE = 200
DROPOUT_RATE_CNN = DROPOUT_INPUT_EMBEDDING = DROPOUT_OUT_LAYER = Float32(0.5)

# 1. Dropout before conv, Max poll layer, PADDING on both sides, hyperparams = window size, output vector size.
W_Char_Embed = W_Char_Embed |> gpu
conv1 = Conv((CHAR_EMBED_DIMS, CONV_WINDOW_LENGTH), 1=>CNN_OUTPUT_SIZE, pad=(0,2)) |> gpu
dropout1 = Dropout(DROPOUT_RATE_CNN) |> gpu

slice_embed(x) = W_Char_Embed * x
char_features = Chain(x -> slice_embed.(x),
                      x -> cat(x..., dims=4),
                      dropout1,
                      conv1,
                      x -> maximum(x, dims=2),
                      x -> reshape(x, CNN_OUTPUT_SIZE, BATCH_SIZE))
# 2. Input embeddings

W_word_Embed = param(W_word_Embed) |> gpu # For trainable
get_word_embedding(w) = W_word_Embed * w |> gpu # works coz - onehot

# Dropout before LSTM
dropout_embed = Dropout(DROPOUT_INPUT_EMBEDDING) |> gpu
input_embeddings((w, cs)) = dropout_embed(vcat(get_word_embedding(w), char_features(cs))) |> gpu

# 3. Bi-LSTM

if device == :gpu
    forward_lstm = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE) |> gpu
    backward = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE) |> gpu
else
    forward_lstm = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE)
    backward = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE)
end
backward_lstm(x) = reverse(backward.(reverse(x)))

bilstm_layer(x) = vcat.(forward_lstm.(x), backward_lstm(x))

# 4. Softmax / CRF Layer
# Dropout after LSTM

if device == :gpu
    dropout2 = Dropout(DROPOUT_OUT_LAYER) |> gpu
    d_out = Dense(LSTM_STATE_SIZE * 2, DENSE_OUT_SIZE + 2) |> gpu
else
    dropout2 = Dropout(DROPOUT_OUT_LAYER)
    d_out = Dense(LSTM_STATE_SIZE * 2, DENSE_OUT_SIZE + 2)
end

m = Chain(x -> input_embeddings.(x),
                bilstm_layer,
                x -> dropout2.(x),
                x -> d_out.(x))

if device == :gpu
    m = gpu(m)
end
# dropout.(bilstm_layer(input_embeddings.(w_cs))) |> gpu

using TextAnalysis: crf_loss, CRF
Flux.@treelike TextAnalysis.CRF

c = TextAnalysis.CRF(num_labels)

init_α = Float32.(fill(-10000, (c.n + 2, 1)))
init_α[c.n + 1] = 0
if device == :gpu
    init_α = cu(init_α)
    c = c |> gpu
end

geti(x, i) = x[:, i]
loss(x, y) = sum([crf_loss(c, [geti(x[j], i) for j in 1:length(y[i])], y[i], init_α) for i in 1:BATCH_SIZE])

η = 0.01 # learning rate
β = 0.05 # rate decay
ρ = 0.9 # momentum

opt = Flux.Optimiser(ExpDecay(β), Momentum(η, ρ))
data = zip(X_input_train, Y_oh_train)

println("Model Built")

function try_outs()
    Flux.testmode!(dropout1)
    Flux.testmode!(dropout2)
    Flux.testmode!(dropout_embed)

    confusion_matrix = zeros(Int, (num_labels, num_labels))
    for (x, y) in zip(X_input_dev, Y_oh_dev)
        reset!(forward_lstm)
        reset!(backward)
        x_seq = m(x)
        sum([crf_loss(c, [geti(x_seq[j], i) for j in 1:length(y[i])], y[i], init_α) for i in 1:BATCH_SIZE])

        ohs = [TextAnalysis.viterbi_decode(cpu(c), cpu.([geti(x_seq[j], i) for j in 1:length(y[i])]), cpu(init_α)) for i in 1:BATCH_SIZE]

        for i in 1:BATCH_SIZE
            for d in zip(ohs[i], y[i])
                confusion_matrix[d[1].ix, d[2].ix] += 1
            end
        end
    end

    s1 = sum(confusion_matrix, dims = 2)
    s2 = sum(confusion_matrix, dims = 1)
    dg = diag(confusion_matrix)

    a = sum([s1[i] == 0 ? 0 : dg[i]/s1[i] for i in eachindex(dg)]) / num_labels
    b = sum(dg ./ s2') / num_labels

    println("Precision and recall are:", a, " ", b)
    println("F1 is:", (2 * a * b) / (a + b))

    Flux.testmode!(dropout1, false)
    Flux.testmode!(dropout2, false)
    Flux.testmode!(dropout_embed, false)
end

function save_weights(epoch, conv1, W_word_Embed, W_Char_Embed, forward_lstm, backward, d_out, c)
    conv_cpu = conv1 |> cpu
    W_word_cpu = W_word_Embed |> cpu
    W_char_cpu = W_Char_Embed |> cpu
    forward_lstm_cpu = forward_lstm |> cpu
    backward_lstm_cpu = backward |> cpu
    crf_cpu = c |> cpu
    d_cpu = d_out |> cpu

    @save "./weights/conv_cpu$(epoch).bson" conv_cpu
    @save "./weights/W_word_cpu$(epoch).bson" W_word_cpu
    @save "./weights/W_char_cpu$(epoch).bson" W_char_cpu
    @save "./weights/forward_lstm$(epoch).bson" forward_lstm_cpu
    @save "./weights/backward_lstm$(epoch).bson" backward_lstm_cpu
    @save "./weights/d_cpu$(epoch).bson" d_cpu
    @save "./weights/crf$(epoch).bson" crf_cpu
end

UPPER_BOUND = Float32(5.0)

# Gradient Clipping
grad_clipping(g) = min(g, UPPER_BOUND)

# Backward
function back!(p::Flux.Params, l, opt)
    l = Tracker.hook(x -> grad_clipping(x), l)
    grads = Tracker.gradient(() -> l, p)
    Tracker.update!(opt, p, grads)
    return
end

function train(EPOCHS)
    reset!(forward_lstm)
    reset!(backward)

    ps = params(params(conv1)..., params(W_word_Embed)..., params(W_Char_Embed)..., params(forward_lstm, backward)..., params(d_out)..., params(c)...)
    for epoch in 1:EPOCHS
        if epoch % 2 == 1
            save_weights(epoch, conv1, W_word_Embed, W_Char_Embed, forward_lstm, backward, d_out, c)
        end
        reset!(forward_lstm)
        reset!(backward)
        try_outs()
        println("----------------------- Starting with epoch : $(epoch) ----------------------")

        for d in data
            reset!(forward_lstm)
            reset!(backward)

            l = loss(m(d[1]), d[2])
            back!(ps, l, opt)
        end
        reset!(forward_lstm)
        reset!(backward)
    end
end

train(100)
reset!(forward_lstm)
reset!(backward)
# save_weights(51, conv_cpu, W_word_Embed, W_Char_Embed, forward_lstm, backward, d_out, c)
try_outs()


text:
using ImageCore, Colors, FixedPointNumbers, ColorVectorSpace
using Base.Test

@testset "Image traits" begin
    for (B,S) in ((rand(UInt16(1):UInt16(20), 3, 5),"Gray"),
                  (rand(Gray{Float32}, 3, 5),"Gray"),
                  (rand(RGB{Float16}, 3, 5),"RGB"),
                  (bitrand(3, 5),"Binary"),
                  (rand(UInt32, 3, 5),"RGB24"))
        @test pixelspacing(B) == (1,1)
        @test spacedirections(B) == ((1,0),(0,1))
        @test sdims(B) == 2
        @test coords_spatial(B) == (1,2)
        @test nimages(B) == 1
        @test size_spatial(B) == (3,5)
        @test indices_spatial(B) == (Base.OneTo(3), Base.OneTo(5))
        assert_timedim_last(B)
        @test width(B) == 3
        @test height(B) == 5
    end
end

@testset "Utilities" begin
    using ImageCore.permutation
    @test permutation((1,2,3), (1,2,3)) == [1,2,3]
    @test permutation((3,1,2), (1,2,3)) == [3,1,2]
    @test permutation(["b", "c", "a"], ["a", "b", "c"]) == [2,3,1]
    @test_throws ArgumentError permutation(["b", "c", "a"], ["a", "bb", "c"])
    @test permutation(["b", "c", "a"], ["a", "bb", "b", "c"]) == [3,4,1,2]
    @test permutation(["a", "bb", "b", "c"], ["b", "c", "a"]) == [3,0,1,2]
    @test_throws ArgumentError permutation(["a", "b"], [:cat, :dog])
    @test_throws ArgumentError permutation(["a", "b"], ["a", "a", "b"])
end

nothing

ext: jl


text:
if basename(pwd()) == "aoc"
    cd("2021/25")
end

function loadfile(filename::AbstractString)
    permutedims(reshape(
        [c for c in readchomp(filename) if !isspace(c)],
        length(readline(filename)),
        :
    ))
end

step(src) = first(step!(similar(src), src))

function step!(dest::AbstractMatrix{Char}, src::AbstractMatrix{Char})
    copyto!(dest, src)
    h, w = size(dest)
    for y in 1:h
        for x in 1:w-1
            if src[y, x] == '>' && src[y, x+1] == '.'
                dest[y, x+1] = '>'
                dest[y, x] = '.'
            end
        end
        if src[y, w] == '>' && src[y, 1] == '.'
            dest[y, 1] = '>'
            dest[y, w] = '.'
        end
    end
    for x in 1:w
        for y in 1:h-1
            if src[y, x] == 'v' && dest[y+1, x] == '.'
                dest[y+1, x] = 'v'
                dest[y, x] = '.'
            end
        end
        if src[h, x] == 'v' && src[1, x] != 'v' && dest[1, x] == '.'
            dest[1, x] = 'v'
            dest[h, x] = '.'
        end
    end
    dest, src
end

function prettyprint(grid, msg = "")
    println(msg)
    println.(join(row) for row in eachrow(grid))
    println()
    nothing
end

function part1(filename::AbstractString; verbose::Bool = false)
    previous = loadfile(filename)
    current = step(previous)
    total = 1
    if verbose
        prettyprint(previous, "Initial state:")
        prettyprint(current, "After 1 step:")
    end
    while previous != current
        current, previous = step!(previous, current)
        total += 1
        if verbose
            prettyprint(current, "After $(total) steps:")
        end
    end
    total
end

@assert part1("example.txt") == 58
println("Part 1: ", part1("input.txt"))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using XGBoost

###
# advanced: cutomsized loss function
#

dtrain = DMatrix("../data/agaricus.txt.train")
dtest = DMatrix("../data/agaricus.txt.test")

# note: for customized objective function, we leave objective as default
# note: what we are getting is margin value in prediction
# you must know what you are doing

param = ["max_depth"=>2, "eta"=>1, "silent"=>1]
watchlist  = [(dtest,"eval"), (dtrain,"train")]
num_round = 2

function logregobj(preds::Vector{Float32}, dtrain::DMatrix)
    labels = get_info(dtrain, "label")
    preds = 1.0 ./ (1.0 + exp(-preds))
    grad = preds - labels
    hess = preds .* (1.0-preds)
    return (grad, hess)
end

# user defined evaluation function, return a pair metric_name, result
# NOTE: when you do customized loss function, the default prediction value is margin
# this may make buildin evalution metric not function properly
# for example, we are doing logistic loss, the prediction is score before logistic transformation
# the buildin evaluation error assumes input is after logistic transformation
# Take this in mind when you use the customization, and maybe you need write customized evaluation function
function evalerror(preds::Vector{Float32}, dtrain::DMatrix)
    labels = get_info(dtrain, "label")
    # return a pair metric_name, result
    # since preds are margin(before logistic transformation, cutoff at 0)
    return ("self-error", sum((preds .> 0.0) .!= labels) / float(size(preds)[1]))
end

# training with customized objective, we can also do step by step training
# simply look at xgboost_lib.jl's implementation of train
bst = xgboost(dtrain, num_round, param=param, watchlist=watchlist,
              obj=logregobj, feval=evalerror)


text:
function sumbinom(K::Int)
    res = 0
    for k = 0:K
        res += binomial(K, k)^2
    end
    return res, binomial(2K, K)
end

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
# CUDA pointer types

export CuPtr, CU_NULL, PtrOrCuPtr, CuArrayPtr, CUA_NULL, CuRef, RefOrCuRef

#
# CUDA device pointer
#

# FIXME: should be called CuDevicePtr...

"""
    CuPtr{T}

A memory address that refers to data of type `T` that is accessible from the GPU. A `CuPtr`
is ABI compatible with regular `Ptr` objects, e.g. it can be used to `ccall` a function that
expects a `Ptr` to GPU memory, but it prevents erroneous conversions between the two.
"""
CuPtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type CuPtr{T} 64 end
else
    primitive type CuPtr{T} 32 end
end

# constructor
CuPtr{T}(x::Union{Int,UInt,CuPtr}) where {T} = Base.bitcast(CuPtr{T}, x)

const CU_NULL = CuPtr{Cvoid}(0)

## getters

Base.eltype(::Type{<:CuPtr{T}}) where {T} = T

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::CuPtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{CuPtr{T}}, x::Union{Int,UInt}) where {T} = CuPtr{T}(x)
Int(x::CuPtr)  = Base.bitcast(Int, x)
UInt(x::CuPtr) = Base.bitcast(UInt, x)

# between regular and CUDA pointers
Base.convert(::Type{<:Ptr}, p::CuPtr) =
    throw(ArgumentError("cannot convert a GPU pointer to a CPU pointer"))

# between CUDA pointers
Base.convert(::Type{CuPtr{T}}, p::CuPtr) where {T} = Base.bitcast(CuPtr{T}, p)

# defer conversions to unsafe_convert
Base.cconvert(::Type{<:CuPtr}, x) = x

# fallback for unsafe_convert
Base.unsafe_convert(::Type{P}, x::CuPtr) where {P<:CuPtr} = convert(P, x)

# from arrays
Base.unsafe_convert(::Type{CuPtr{S}}, a::AbstractArray{T}) where {S,T} =
    convert(CuPtr{S}, Base.unsafe_convert(CuPtr{T}, a))
Base.unsafe_convert(::Type{CuPtr{T}}, a::AbstractArray{T}) where {T} =
    error("conversion to pointer not defined for $(typeof(a))")

## limited pointer arithmetic & comparison

Base.isequal(x::CuPtr, y::CuPtr) = (x === y)
Base.isless(x::CuPtr{T}, y::CuPtr{T}) where {T} = x < y

Base.:(==)(x::CuPtr, y::CuPtr) = UInt(x) == UInt(y)
Base.:(<)(x::CuPtr,  y::CuPtr) = UInt(x) < UInt(y)
Base.:(-)(x::CuPtr,  y::CuPtr) = UInt(x) - UInt(y)

Base.:(+)(x::CuPtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::CuPtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::CuPtr) = y + x

#
# Host or device pointer
#

"""
    PtrOrCuPtr{T}

A special pointer type, ABI-compatible with both `Ptr` and `CuPtr`, for use in `ccall`
expressions to convert values to either a GPU or a CPU type (in that order). This is
required for CUDA APIs which accept pointers that either point to host or device memory.
"""
PtrOrCuPtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type PtrOrCuPtr{T} 64 end
else
    primitive type PtrOrCuPtr{T} 32 end
end

function Base.cconvert(::Type{PtrOrCuPtr{T}}, val) where {T}
    # `cconvert` is always implemented for both `Ptr` and `CuPtr`, so pick the first result
    # that has done an actual conversion

    gpu_val = Base.cconvert(CuPtr{T}, val)
    if gpu_val !== val
        return gpu_val
    end

    cpu_val = Base.cconvert(Ptr{T}, val)
    if cpu_val !== val
        return cpu_val
    end

    return val
end

function Base.unsafe_convert(::Type{PtrOrCuPtr{T}}, val) where {T}
    ptr = if Core.Compiler.return_type(Base.unsafe_convert,
                                       Tuple{Type{Ptr{T}}, typeof(val)}) !== Union{}
        Base.unsafe_convert(Ptr{T}, val)
    elseif Core.Compiler.return_type(Base.unsafe_convert,
                                     Tuple{Type{CuPtr{T}}, typeof(val)}) !== Union{}
        Base.unsafe_convert(CuPtr{T}, val)
    else
        throw(ArgumentError("cannot convert to either a CPU or GPU pointer"))
    end

    return Base.bitcast(PtrOrCuPtr{T}, ptr)
end

#
# CUDA array pointer
#

if sizeof(Ptr{Cvoid}) == 8
    primitive type CuArrayPtr{T} 64 end
else
    primitive type CuArrayPtr{T} 32 end
end

# constructor
CuArrayPtr{T}(x::Union{Int,UInt,CuArrayPtr}) where {T} = Base.bitcast(CuArrayPtr{T}, x)

## getters

Base.eltype(::Type{<:CuArrayPtr{T}}) where {T} = T

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::CuArrayPtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{CuArrayPtr{T}}, x::Union{Int,UInt}) where {T} = CuArrayPtr{T}(x)
Int(x::CuArrayPtr)  = Base.bitcast(Int, x)
UInt(x::CuArrayPtr) = Base.bitcast(UInt, x)

# between regular and CUDA pointers
Base.convert(::Type{<:Ptr}, p::CuArrayPtr) =
    throw(ArgumentError("cannot convert a GPU array pointer to a CPU pointer"))

# between CUDA array pointers
Base.convert(::Type{CuArrayPtr{T}}, p::CuArrayPtr) where {T} = Base.bitcast(CuArrayPtr{T}, p)

# defer conversions to unsafe_convert
Base.cconvert(::Type{<:CuArrayPtr}, x) = x

# fallback for unsafe_convert
Base.unsafe_convert(::Type{P}, x::CuArrayPtr) where {P<:CuArrayPtr} = convert(P, x)

#
# CUDA reference objects
#

if sizeof(Ptr{Cvoid}) == 8
    primitive type CuRef{T} 64 end
else
    primitive type CuRef{T} 32 end
end

# general methods for CuRef{T} type
Base.eltype(x::Type{<:CuRef{T}}) where {T} = @isdefined(T) ? T : Any

Base.convert(::Type{CuRef{T}}, x::CuRef{T}) where {T} = x

# conversion or the actual ccall
Base.unsafe_convert(::Type{CuRef{T}}, x::CuRef{T}) where {T} = Base.bitcast(CuRef{T}, Base.unsafe_convert(CuPtr{T}, x))
Base.unsafe_convert(::Type{CuRef{T}}, x) where {T} = Base.bitcast(CuRef{T}, Base.unsafe_convert(CuPtr{T}, x))

# CuRef from literal pointer
Base.convert(::Type{CuRef{T}}, x::CuPtr{T}) where {T} = x

# indirect constructors using CuRef
CuRef(x::Any) = CuRefArray(CuArray([x]))
CuRef{T}(x) where {T} = CuRefArray{T}(CuArray(T[x]))
CuRef{T}() where {T} = CuRefArray(CuArray{T}(undef, 1))
Base.convert(::Type{CuRef{T}}, x) where {T} = CuRef{T}(x)

## CuRef object backed by a CUDA array at index i

struct CuRefArray{T,A<:AbstractArray{T}} <: Ref{T}
    x::A
    i::Int
    CuRefArray{T,A}(x,i) where {T,A<:AbstractArray{T}} = new(x,i)
end
CuRefArray{T}(x::AbstractArray{T}, i::Int=1) where {T} = CuRefArray{T,typeof(x)}(x, i)
CuRefArray(x::AbstractArray{T}, i::Int=1) where {T} = CuRefArray{T}(x, i)
Base.convert(::Type{CuRef{T}}, x::AbstractArray{T}) where {T} = CuRefArray(x, 1)

function Base.unsafe_convert(P::Type{CuPtr{T}}, b::CuRefArray{T}) where T
    return pointer(b.x, b.i)
end
function Base.unsafe_convert(P::Type{CuPtr{Any}}, b::CuRefArray{Any})
    return convert(P, pointer(b.x, b.i))
end
Base.unsafe_convert(::Type{CuPtr{Cvoid}}, b::CuRefArray{T}) where {T} =
    convert(CuPtr{Cvoid}, Base.unsafe_convert(CuPtr{T}, b))

## Union with all CuRef 'subtypes'

const CuRefs{T} = Union{CuPtr{T}, CuRefArray{T}}

## RefOrCuRef

if sizeof(Ptr{Cvoid}) == 8
    primitive type RefOrCuRef{T} 64 end
else
    primitive type RefOrCuRef{T} 32 end
end

Base.convert(::Type{RefOrCuRef{T}}, x::Union{RefOrCuRef{T}, Ref{T}, CuRef{T}, CuRefs{T}}) where {T} = x

# prefer conversion to CPU ref: this is generally cheaper
Base.convert(::Type{RefOrCuRef{T}}, x) where {T} = Ref{T}(x)
Base.unsafe_convert(::Type{RefOrCuRef{T}}, x::Ref{T}) where {T} =
    Base.bitcast(RefOrCuRef{T}, Base.unsafe_convert(Ptr{T}, x))
Base.unsafe_convert(::Type{RefOrCuRef{T}}, x) where {T} =
    Base.bitcast(RefOrCuRef{T}, Base.unsafe_convert(Ptr{T}, x))

# support conversion from GPU ref
Base.unsafe_convert(::Type{RefOrCuRef{T}}, x::CuRefs{T}) where {T} =
    Base.bitcast(RefOrCuRef{T}, Base.unsafe_convert(CuPtr{T}, x))

# support conversion from arrays
Base.convert(::Type{RefOrCuRef{T}}, x::Array{T}) where {T} = convert(Ref{T}, x)
Base.convert(::Type{RefOrCuRef{T}}, x::AbstractArray{T}) where {T} = convert(CuRef{T}, x)
Base.unsafe_convert(P::Type{RefOrCuRef{T}}, b::CuRefArray{T}) where T =
    Base.bitcast(RefOrCuRef{T}, Base.unsafe_convert(CuRef{T}, b))


text:
export
    alg3  # Algorithm for learning neural networks

"""
    Algorithm for learning neural networks (we only implement one layer yet.)

arguments:

  features,     column vectors
  labels
  m,            # of layers
  k,            sample k points (10)
  s,            # of hidden nodes of last layer, not relevant if m = 1
  T,            repeat T rounds
  B,            [-B, B]^k cube (10)
  weight,       default to []

"""
function alg3(features, labels; m=1, k=10, s=1, T=1, B=10, weight=[])
  d, n = size(features) # feature dim

  size(labels) == (1, n) || error("alg2: column number of features should equal to number of labels")
  n >= k || error("alg3: sample size (=$n) should > k (=$k)")
  m == 1 || error("alg3: higher layer network not implemented yet")
  T == 1 || error("alg3: support T = 1 only")

  # sample k points w.r.t. weight
  ind = weight == [] ? sample(1:n, k, replace=false) :
    sample(1:n, WeightVec(weight), k, replace=false)

  x = features[:, ind]

  # draw u
  u = B * rand(k)'

  # find v
  v = opt(x, u, B)
end

max_stars_repo_name: innerlee/Halfspace

max_issues_repo_name: innerlee/Halfspace

max_forks_repo_name: innerlee/Halfspace


text:
######################################################################
# descender.jl
# provides descent engine
######################################################################
export Descender

export InverseStep
export InverseRootStep
export ArmijoLineSearch
export ConstantStep

struct InverseStep <: Descender
    alpha::Float64
end

function get_step_size(ir::InverseStep, em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF, zd::VVF, vd::VVF, ad::MF, bd::MF, K::Vector{MF}, i::Int)
    return ir.alpha / i
end

struct InverseRootStep <: Descender
    alpha::Float64
end

#get_step_size(ir::InverseRootStep, i::Int) = ir.alpha / sqrt(i)
function get_step_size(ir::InverseRootStep, em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF, zd::VVF, vd::VVF, ad::MF, bd::MF, K::Vector{MF}, i::Int)
    return ir.alpha / sqrt(i)
end

struct ConstantStep <: Descender
    alpha::Float64
end

function get_step_size(cs::ConstantStep, em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF, zd::VVF, vd::VVF, ad::MF, bd::MF, K::Vector{MF}, i::Int)
    return cs.alpha
end

include("armijo_line_search.jl")


text:
using Catalyst
using DiffEqJump
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

sir_model = @reaction_network begin
  β*c/(s+i+r), s + i --> 2i
  γ, i --> r
end β c γ;

tmax = 40.0
tspan = (0.0,tmax);

δt = 0.1
t = 0:δt:tmax;

u0 = [990,10,0]; # S,I,R

p = [0.05,10.0,0.25];

Random.seed!(1234);

prob_discrete = DiscreteProblem(sir_model,u0,tspan,p);

prob_jump = JumpProblem(sir_model,prob_discrete,Direct());

sol_jump = solve(prob_jump,SSAStepper());

out_jump = sol_jump(t);

df_jump = DataFrame(out_jump')
df_jump[!,:t] = out_jump.t;

@df df_jump plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_jump,SSAStepper())


text:
using ImplicitArrays
using Test

include("../src/compat.jl")

@testset "ImplicitArrays.jl" begin
    @testset "BlockMatrix" begin include("block.jl") end
    @testset "InteractionMatrix" begin include("interaction.jl") end
    @testset "FixedValueArray" begin include("fixedvalue.jl") end
    @testset "RowProjectionMatrix" begin include("rowprojection.jl") end
end


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module DataUtils
include("../shared/load.jl")
using GMMData
using BAData
using HandData
using LSTMData
using Printf

export create_empty_output_for, save_output_to_file
export objective_file_name, jacobian_file_name, times_file_name, save_time_to_file, save_value_to_file, save_vector_to_file, save_matrix_to_file
export save_gmm_output_to_file

"""
    create_empty_output_for(input)

For given input object creates a corresponding empty output object.
"""
function create_empty_output_for end
create_empty_output_for(::GMMInput) = empty_gmm_output()
create_empty_output_for(::BAInput) = empty_ba_output()
create_empty_output_for(::HandInput) = empty_hand_output()
create_empty_output_for(::LSTMInput) = empty_lstm_output()

"""
    save_output_to_file(output, output_prefix, input_name, module_name)

Saves provided output to a correctly named file.
"""
function save_output_to_file end
save_output_to_file(output::GMMOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = save_gmm_output_to_file(output, output_prefix, input_name, module_name)
save_output_to_file(output::BAOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = save_ba_output_to_file(output, output_prefix, input_name, module_name)
save_output_to_file(output::HandOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = save_hand_output_to_file(output, output_prefix, input_name, module_name)
save_output_to_file(output::LSTMOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = save_lstm_output_to_file(output, output_prefix, input_name, module_name)

format(x::Float64) = @sprintf "%.12g" x
format(x::Int) = @sprintf "%d" x

objective_file_name(output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = "$(output_prefix)$(input_name)_F_$(module_name).txt"
jacobian_file_name(output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = "$(output_prefix)$(input_name)_J_$(module_name).txt"
times_file_name(output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString) = "$(output_prefix)$(input_name)_times_$(module_name).txt"

function save_time_to_file(filepath::AbstractString, objective_time::Float64, derivative_time::Float64)
    open(filepath, "w") do io
        println(io, format(objective_time))
        print(io, format(derivative_time))
    end
end

function save_value_to_file(filepath::AbstractString, x::T) where T
    open(filepath, "w") do io
        print(io, format(x))
    end
end

function save_vector_to_file(filepath::AbstractString, v::Vector{T}) where T
    open(filepath, "w") do io
        for x ∈ v
            println(io, format(x))
        end
    end
end

function save_matrix_to_file(filepath::AbstractString, m::Matrix{T}) where T
    open(filepath, "w") do io
        for i ∈ 1:size(m, 1)
            join(io, map(j -> format(m[i, j]), 1:size(m, 2)), "\t")
            println(io)
        end
    end
end

function save_errors_to_file(filepath::AbstractString, reproj_err::Matrix{Float64}, w_err::Vector{Float64})
    open(filepath, "w") do io
        println(io, "Reprojection error:")
        for i ∈ reproj_err
            println(io, format(i))
        end
        println(io, "Zach weight error:")
        for i ∈ w_err
            println(io, format(i))
        end
    end
end

function save_sparse_j_to_file(filepath::AbstractString, jacobian::BASparseMatrix)
    open(filepath, "w") do io
        println(io, format(jacobian.nrows), " ", format(jacobian.ncols))
        println(io, format(size(jacobian.rows, 1)))
        join(io, map(format, jacobian.rows), " ")
        println(io)
        println(io, format(size(jacobian.cols, 1)))
        join(io, map(format, jacobian.cols), " ")
        println(io)
        join(io, map(format, jacobian.vals), " ")
    end
end

function save_gmm_output_to_file(output::GMMOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString)
    save_value_to_file(objective_file_name(output_prefix, input_name, module_name), output.objective)
    save_vector_to_file(jacobian_file_name(output_prefix, input_name, module_name), output.gradient)
end

function save_ba_output_to_file(output::BAOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString)
    save_errors_to_file(objective_file_name(output_prefix, input_name, module_name), output.reproj_err, output.w_err)
    save_sparse_j_to_file(jacobian_file_name(output_prefix, input_name, module_name), output.jacobian)
end

function save_hand_output_to_file(output::HandOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString)
    save_vector_to_file(objective_file_name(output_prefix, input_name, module_name), output.objective)
    save_matrix_to_file(jacobian_file_name(output_prefix, input_name, module_name), output.jacobian)
end

function save_lstm_output_to_file(output::LSTMOutput, output_prefix::AbstractString, input_name::AbstractString, module_name::AbstractString)
    save_value_to_file(objective_file_name(output_prefix, input_name, module_name), output.objective)
    save_vector_to_file(jacobian_file_name(output_prefix, input_name, module_name), output.gradient)
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
# Objects in the game and their supporting functions (update!, collide!, ...)

"""
    WorldPos(5.0,-200.0)
x,y float coordinates in the game world (not necessarily the same as pixel
coordinates on the screen).
"""
struct WorldPos  # 0,0 == middle
    x::Float64
    y::Float64
end
"""
    Vector2D(-2.5,1.0)
x,y vector representing direction in the game world. Could represent a velocity,
a distance, etc. Subtracting two `WorldPos`itions results in a `Vector2D`.
"""
struct Vector2D
    x::Float64
    y::Float64
end
import Base.*, Base./, Base.-, Base.+
+(a::Vector2D, b::Vector2D) = Vector2D(a.x+b.x, a.y+b.y)
-(a::Vector2D, b::Vector2D) = Vector2D(a.x-b.x, a.y-b.y)
*(a::Vector2D, x::Number) = Vector2D(a.x*x, a.y*x)
*(x::Number, a::Vector2D) = a*x
/(a::Vector2D, x::Number) = Vector2D(a.x/x, a.y/x)
+(a::WorldPos, b::Vector2D) = WorldPos(a.x+b.x, a.y+b.y)
-(a::WorldPos, b::Vector2D) = WorldPos(a.x-b.x, a.y-b.y)
+(a::Vector2D, b::WorldPos) = WorldPos(a.x+b.x, a.y+b.y)
-(a::Vector2D, b::WorldPos) = WorldPos(a.x-b.x, a.y-b.y)
-(a::WorldPos, b::WorldPos) = Vector2D(a.x-b.x, a.y-b.y)
-(x::WorldPos) = WorldPos(-x.x, -x.y)
-(x::Vector2D) = Vector2D(-x.x, -x.y)

mutable struct Ball
    pos::WorldPos
    vel::Vector2D
end
mutable struct Paddle
    pos::WorldPos
    vel::Vector2D
    length
end

pingSound = nothing
scoreSound = nothing

collide!(a::Ball, b::Paddle) = collide!(b,a)
function collide!(p::Paddle, b::Ball)
    xIncr = 100
    xSign = 1
    if b.pos.x - p.pos.x > p.length/4; # In right quarter
        if b.vel.x < 0 ; xSign=-1; end
    elseif b.pos.x - p.pos.x < -p.length/4; # In left
        xIncr *=-1; if b.vel.x > 0 ; xSign=-1; end
    else xIncr = 0;
    end
    b.vel = Vector2D(b.vel.x * xSign + xIncr, -b.vel.y)

    audioEnabled && SDL2.Mix_PlayChannel( Int32(-1), pingSound, Int32(0) )
end

"""
    LineSegment(WorldPos(0,0), WorldPos(1,1))
Line segment used for collision detection algorithm.
"""
struct LineSegment
     a::WorldPos
     b::WorldPos
end
""" Will they collide on the next update? """
willCollide(b::Ball, p::Paddle, dt) = willCollide(p,b,dt)
function willCollide(p::Paddle, b::Ball, dt)
     # If the ball is in the right X-axis vicinity of the paddle.
    if (abs(p.pos.x - b.pos.x) <= (p.length/2.0+ballWidth/2.0+abs(b.vel.x)*ballWidth))
        l = LineSegment(b.pos, b.pos+b.vel*dt) # If next update will bring collision.
        return isColliding(p, l, ballWidth)  # will it cross the paddle in the Y-axis
    else
        return false
    end
end
function isColliding(p::Paddle, l::LineSegment, width)
    v = l.b - l.a
    if v.y != 0
        c0 = (p.pos.y - l.a.y) / v.y
        if !(0.0 <= c0 <= 1.0) return false end
        lc0 = (c0*v + l.a)
    else
        lc0 = l.a
        if (l.a.y != p.pos.y) return false end
    end
    return abs(lc0.x - p.pos.x) <= (p.length/2.0 + width/2.0)
end

""" Perform game updates for `b` given `dt` seconds since last update. """
function update!(b::Ball, dt)
    global scoreA, scoreB
    b.pos = b.pos + (b.vel * dt)

    if b.pos.y > winHeight[]/2.0
        scoreB += 1
        audioEnabled && SDL2.Mix_PlayChannel( Int32(-1), scoreSound, Int32(0) )
        b.pos = WorldPos(0,0)
        b.vel = Vector2D(rand(-ballSpeed:ballSpeed), rand([ballSpeed,-ballSpeed]))
    elseif b.pos.y < -winHeight[]/2.0
        scoreA += 1
        audioEnabled && SDL2.Mix_PlayChannel( Int32(-1), scoreSound, Int32(0) )
        b.pos = WorldPos(0,0)
        b.vel = Vector2D(rand(-ballSpeed:ballSpeed), rand([ballSpeed,-ballSpeed]))
    end
    if b.pos.x > winWidth[]/2.0
        b.vel = Vector2D(-abs(b.vel.x), b.vel.y)
    elseif b.pos.x < -winWidth[]/2.0
        b.vel = Vector2D(abs(b.vel.x), b.vel.y)
    end
end
""" Game updates for `p` with `keys` pressed and `dt` seconds since last update. """
function update!(p::Paddle, keys, dt)
    # Apply velocity
    p.pos = p.pos + (p.vel * dt)

    # TODO: another way i could calculate this, which might give more control is
    # to set the position manually as a function of how long they've been
    # holding the button. That way, if I want, i could implement normal
    # acceleration as I've done here, but OR if i want, i could also do like t^3
    # instead of t^2 to make it more dramatic or something.

    # Calculate acceleration from user input.
    accel = 0
    decelerating = false  # keep track of whether speeding up or slowing down.
    cur_vel_sign = 0
    # If both keys are down or neither key is down, decelerate.
    if (!xor(keys.leftDown, keys.rightDown))
        if abs(p.vel.x) > 0
            cur_vel_sign = p.vel.x / abs(p.vel.x)
            accel = paddleDeccel * -1 * cur_vel_sign
            decelerating = true
        end
    elseif (keys.leftDown)
        accel = -paddleAccel
        # Boost if switching directions
        if p.vel.x > 0;
            accel += -paddleDeccel
        end
    elseif (keys.rightDown)
        accel = paddleAccel
        # Boost if switching directions
        if p.vel.x < 0;
            accel += paddleDeccel
        end
    end

    # Apply accel
    p.vel = p.vel + Vector2D(accel * dt, 0)
    if decelerating
        # If decelerating pushed it "past" 0, bring it to a stop.
        if (cur_vel_sign < 0 && p.vel.x > 0
           || cur_vel_sign > 0 && p.vel.x < 0)
            p.vel = Vector2D(0, p.vel.y)
        end
    else
        # If accelerating pushed it past max, bring it to a stop.
        if p.vel.x > paddleSpeed
            p.vel = Vector2D(paddleSpeed, p.vel.y)
        elseif p.vel.x < -paddleSpeed
            p.vel = Vector2D(-paddleSpeed, p.vel.y)
        end
    end

    # Check position bounds.
    if p.pos.x > winWidth[]/2.0
        p.pos = WorldPos(winWidth[]/2.0, p.pos.y)
    elseif p.pos.x < -winWidth[]/2.0
        p.pos = WorldPos(-winWidth[]/2.0, p.pos.y)
    end
end

max_stars_repo_name: NHDaly/pongClone

max_issues_repo_name: NHDaly/pongClone

max_forks_repo_name: NHDaly/pongClone


text:
"""
`module  JAC.Semiempirical`
    ... a submodel of JAC that contains all methods to set-up and process (simple)  semiempirical estimations of
        atomic properties that cannot be calculated so easily in a rigorious manner.
"""
module Semiempirical

    using ..Basics, ..Defaults, ..ManyElectron, ..PeriodicTable    ## ..ImpactIonization

    """
    `@enum   AtomicCrossSection`  ... defines a enumeration for (supported) cross section estimation.
    """
    @enum   AtomicCrossSection    NoCrossSection    ImpactIonizationCS

    """
    `Semiempirical.AtomicCrossSection(sa::String)`  ... constructor for a given String.
    """
    function AtomicCrossSection(sa::String)
        if       sa == "No cross section"      wa = NoCrossSection
        elseif   sa == "impact-ionization"     wa = ImpactIonizationCS
        else
            error("stop a")
        end

        return( wa )
    end

    # `Base.show(io::IO, cs::AtomicCrossSection)`  ... prepares a proper printout of the variable cs::AtomicCrossSection.
    ## function Base.show(io::IO, cs::AtomicCrossSection)
    ## print(io, string(cs) )
    ## end

    # `Base.string(cs::AtomicCrossSection)`  ... provides a proper printout of the variable cs::AtomicCrossSection.
    function Base.string(cs::AtomicCrossSection)
        if       cs == NoCrossSection       return("No cross section")
        elseif   cs == ImpactIonizationCS   return("impact-ionization cross section")
        else     error("stop a")
        end
    end

    #=========================================
    """
    `struct  Estimation`
        ... defines a type for defining  (a framework for) semiempirical estimates of different properties that cannot be calculated
            rigoriously in JAC.

        + calcCrossSections    ::Bool                        ... True, if atomic cross sections of some kind are to be calculated
        + crossSection         ::AtomicCrossSection          ... Selected cross sections type.
        + config               ::Configuration               ... A non-relativistic configurations for which computations are to be carried out.
        + crossSectionSettings ::Union{ImpactIonization.Settings} ... Provides the settings for the selected cross section estimates.
    """
    struct  Estimation
        calcCrossSections      ::Bool
        crossSection           ::AtomicCrossSection
        config                 ::Configuration
        crossSectionSettings   ::Union{ImpactIonization.Settings}
    end

    """
    `Semiempirical.Estimation()`  ... constructor for an 'empty' instance::Semiempirical.Estimation.
    """
    function Estimation()
        Estimation(false, NoCrossSection, Configuration(), ImpactIonization.Settings() )
    end

    # `Base.show(io::IO, estimation::Semiempirical.Estimation)`  ... prepares a proper printout of the variable estimation::Semiempirical.Estimation.
    function Base.show(io::IO, estimation::Semiempirical.Estimation)
        println(io, "calcCrossSections:        $(estimation.calcCrossSections)  ")
        println(io, "crossSection:             $(estimation.crossSection)  ")
        println(io, "config:                   $(estimation.config)  ")
        println(io, "crossSectionSettings:     $(estimation.crossSectionSettings)  ")
    end

    """
    `Semiempirical.Estimation("interactive")`
        ... constructor to generate a new instance of Semiempirical.Estimation interactively by replying to some detailed dialog.
            **Not yet fully implemented !**
    """
    function Semiempirical.Estimation(sa::String)
        sa != "interactive"   &&    error("Unsupported keystring = $sa.")
    end
    =====================================================#

    """
    `Semiempirical.estimate("ionization potential: inner-shell", shell::Shell, Z::Int64)`
        ... to estimate the ionization potential = mean binding energy of an electron in the given shell; an energy value::Float64
            is returned. These ionization potentials are taken from the tabulation ....
    """
    function  estimate(sa::String, shell::Shell, Z::Int64)
        if      sa == "ionization potential: inner-shell"
            shellIndex = Dict( Subshell("1s_1/2") => 1, Subshell("2s_1/2") => 2, Subshell("2p_1/2") => 3, Subshell("2p_3/2") => 4,
                            Subshell("3s_1/2") => 5, Subshell("3p_1/2") => 6, Subshell("3p_3/2") => 7,
                            Subshell("3d_3/2") => 8, Subshell("3d_3/2") => 9  )
            idx = shellIndex[shell]
            return( PeriodicTable.store("binding energies: Williams (2000)", Z)[idx] )
        else
            error("Unsupported keystring = $sa")
        end
    end

    """
    `Semiempirical.estimate("binding energy", Z::Float64, sh::Subshell)`
        ... to provide the binding energy of a subshell electron, taken from a semi-empirical tabulations by Williams et al.,
            https://userweb.jlab.org/~gwyn/ebindene.html. A energy::Float64 in  Hartree is returned.
    """
    function estimate(sa::String, Z::Int64, sh::Subshell)
        if     sa == "binding energy"
            wa = PeriodicTable.store_Williams2000(Z)
            if      sh == Subshell("1s_1/2")    wb = wa[1]
            elseif  sh == Subshell("2s_1/2")    wb = wa[2]
            elseif  sh == Subshell("2p_1/2")    wb = wa[3]
            elseif  sh == Subshell("2p_3/2")    wb = wa[4]
            elseif  sh == Subshell("3s_1/2")    wb = wa[5]
            elseif  sh == Subshell("3p_1/2")    wb = wa[6]
            elseif  sh == Subshell("3p_3/2")    wb = wa[7]
            elseif  sh == Subshell("3d_3/2")    wb = wa[8]
            elseif  sh == Subshell("3d_5/2")    wb = wa[9]
            elseif  sh == Subshell("4s_1/2")    wb = wa[10]
            elseif  sh == Subshell("4p_1/2")    wb = wa[11]
            elseif  sh == Subshell("4p_3/2")    wb = wa[12]
            else    error("No binding energy available for Z = $Z and subshell $sh ")
            end
        else   error("Unsupported keystring")
        end
        #
        if     wb == -1.   error("No binding energy available for Z = $Z and subshell $sh ")
        else   wb = Defaults.convertUnits("energy: from eV to atomic", wb)
        end
        #
        return( wb )
    end

    """
    `Semiempirical.estimate("binding energy", Z::Float64, conf::Configuration)`
        ... to provide an approximate binding energy of a given electron configuration. This estimate adds the binding
            energies of all subshell, taken frogm a semi-empirical tabulations by Williams et al.,
            https://userweb.jlab.org/~gwyn/ebindene.html. If no binding energy is provided by this table, it simply scales the
            binding energies (e/ 12.) from the next lower shell of the same symmetry
            No relaxation effects are included if several hole states
            occur with regard to the neutral atom. An energy::Float64 in  Hartree is returned.
    """
    function estimate(sa::String, Z::Int64, conf::Configuration)
        if     sa == "binding energy"
            wa = PeriodicTable.bindingEnergies_Williams2000(Z);    wb = 0.
            for (sh,v) in  conf.shells
                if      sh == Shell("1s")    if wa[1]  == -1.   error("stop aa")          else    wb = wb + v * wa[1]    end
                elseif  sh == Shell("2s")    if wa[2]  == -1.   error("stop ab")          else    wb = wb + v * wa[2]    end
                elseif  sh == Shell("2p")    if wa[3]  == -1.   wb = wb + v * wa[2]/12.   else    wb = wb + v * wa[3]    end
                elseif  sh == Shell("3s")    if wa[5]  == -1.   wb = wb + v * wa[2]/12.   else    wb = wb + v * wa[5]    end
                                             if wa[2]  == -1.   error("stop ad")   end
                elseif  sh == Shell("3p")    if wa[6]  == -1.   wb = wb + v * wa[3]/12.   else    wb = wb + v * wa[6]    end
                                             if wa[3]  == -1.   error("stop af")   end
                elseif  sh == Shell("3d")    if wa[8]  == -1.                             else    wb = wb + v * wa[8]    end
                elseif  sh == Shell("4s")    if wa[10] == -1.                             else    wb = wb + v * wa[10]   end
                elseif  sh == Shell("4p")    if wa[11] == -1.                             else    wb = wb + v * wa[11]   end
                else    ## error("No binding energy available for Z = $Z and subshell $sh ")
                end
            end
        else   error("Unsupported keystring")
        end
        #
        wb = Defaults.convertUnits("energy: from eV to atomic", wb)
        #
        return( wb )
    end

end # module


text:
######################################################################
# plotting.jl
# Handles all the calls to PyPlot to generate plots
######################################################################

"""
`plot_b(m::SearchDomain, X::VehicleSet)`

Plots the belief, jammer, and vehicles.
"""
function plot_b(m::SearchDomain, X::VehicleSet)
	#hold(true)
	#plot_b(m.b, X, m.theta)
	plot_theta(m)
	plot_vehicles(X)
	imshow(m.b', interpolation="none", cmap="Greys", origin="lower")
	labels()
end

# This is needed for plot_sim
function plot_b(m::SearchDomain, b::Belief, X::VehicleSet)
	plot_theta(m)
	hold(true)
	plot_vehicles(X)
	imshow(b', interpolation="none", cmap="Greys", origin="lower")
	labels()
end

"""
`plot_eid(m::SearchDomain, X::VehicleSet)`

Plots the eid using the Fisher information matrix.
"""
function plot_eid(m::SearchDomain, X::VehicleSet)
	eid = EID(m,m.b)
	imshow(eid', interpolation="none", cmap="Greys", origin="lower")
	plot_contour(m, eid)
	plot_theta(m)
	plot_vehicles(X)
	labels()
end

"""
`plot_mi(m::SearchDomain, X::VehicleSet)`

Plots the mutual information.
"""
function plot_mi(m::SearchDomain, X::VehicleSet)
	mut_info = mutual_information(m)
	plot_vehicles(X)
	imshow(mut_info', interpolation="none", cmap="Greys", origin="lower")
end

"""
`plot_sim(m::SearchDomain, s::Simulation)`

Steps through a simulation.
"""
function plot_sim(m::SearchDomain, s::Simulation)
	for t = 0:s.T
		hold(false)
		pause(1)
		plot_b(m, s.belief_list[t+1], s.state_list[t+1])
	end
end

######################################################################
# Helper functions
######################################################################
# Plots locations of the vehicles
function plot_vehicles(X::VehicleSet)
	mark_size = 12
	for xi in X
		plot(xi.x, xi.y, "b*", markersize=mark_size)
	end
end

# Plots jammer location
function plot_theta(m::SearchDomain)
	mark_size = 12
	plot(m.theta[1], m.theta[2], "r^", markersize=mark_size)
end

# Plots contours of some distribution `d` (a matrix).
function plot_contour(m::SearchDomain, d)
	X,Y = meshgrid(0:m.num_cells-1, 0:m.num_cells-1)
	contour(X, Y, d')
end

# Sets the appropriate plot labels
function labels()
	xlabel("x")
	ylabel("y")
end

######################################################################
# Old, deprecated, or unused
######################################################################
function plot_eid2(m, b, x, theta)
	mark_size = 12
	eid = EID(m,theta[1],theta[2])
	plot(x[1], x[2], "b*", markersize=mark_size)
	plot(theta[1], theta[2], "r^", markersize=mark_size)
	imshow(eid', interpolation="none", cmap="Greys", origin="lower")
	xlabel("x")
	ylabel("y")
end

function plot_eid3(m, b, x, theta)
	mark_size = 12
	eid = EID(m,theta[1],theta[2])
	#plot(x[1], m.num_cells - 1 - x[2], "b*", markersize=mark_size)
	#plot(theta[1], m.num_cells - 1 - theta[2], "r^", markersize=mark_size)
	#eid_plot = (eid')[end:-1:1, :]

	#imshow(eid_plot, interpolation="none", cmap="Greys")
	plot_contour(m, eid)
	xlabel("x")
	ylabel("y")
end

"""
`meshgrid(x, y)`

Returns X and Y.
"""
function meshgrid(x, y)
	nx = length(x)
	ny = length(y)

	X = zeros(ny, nx)
	Y = zeros(ny, nx)
	for xi = 1:nx
		for xj = 1:ny
			X[xj, xi] = x[xi]
			Y[xj, xi] = y[xj]
		end
	end
	return X, Y
end


text:
export TR_generic_ls
"""
A generic one dimensionnal Trust Region algorithm adapted for line search.
In this generic version we assume a second order polynomial is used to
approximate our h function. A third degree approximation is used in TR_Cub_ls
We use q(d) = h(θ) + h'(θ)⋆d + 0.5 * h''(θ) ⋆ d². Also h''(θ) is sometimes
approximated by a secant (Sec) or an improved secant (SecA). That approximation
is decided with hyper parameter "direction"

Some specific hyper parameter to this line search:
symmetrical :: The method is designed to used an asymmetrical  trust region
to avoid having a negative step size. If symmetrical true then the trust region
becomes symmetrical

Some parameters are used for TR and ARC methods
eps1 and eps2 :: treshold to determine if q is a bad or a good++ approximation
of h
red and aug :: Controls the speed at which we reduce or augment the size of
the trust region depending on the previously mentionned treshold.

All the parameters mentionned above can be found in the LS_Function_Meta
"""
function TR_generic_ls(h :: LineModel,  stop_ls :: LS_Stopping;
                       f_meta = LS_Function_Meta(),
                       φ_dφ        :: Function = (x, y) -> phi_dphi(x, y),
                       verboseLS   :: Bool = false, symmetrical :: Bool = false,
                       kwargs...)

  state = stop_ls.current_state

  tmp = obj(h, 1.0); tmpg = grad(h, 1.0)
  # # Check if 1.0 is an admissible step size.
  OK = update_and_start!(stop_ls, ht = tmp, gt = tmpg, tmps = time())

  # If 1.0 isn't an admissible step size we start our line search algorithm
  iter = 0

  # If 1.0 is a good step size we don't want to compute  unnecessary function
  # evaluations
  H = NaN; φt = NaN; dφt = NaN                   # only for the verbose

  #if OK == false
    φt, dφt = φ_dφ(h, state) #Tangi: Ici on recalcule la fct objectif ?

    tprec = NaN; φtestTR = NaN; dφtestTR = NaN; φtprec = NaN; dφtprec = NaN;

    # H will denote the approximation to φ'' hereafter
    if f_meta.dir == :Nwt
      H = hess(h, state.x)
    elseif f_meta.dir == :Sec || f_meta.dir == :SecA
      H = 1.0
    end
  #end # if OK == false

  # Quadratic approximation of h
  # q(d) = φt + dφt * d + 0.5 * H * d^2

  verboseLS &&
    @printf("  iter   t       φt        dφt        ddφt        Δn         Δp")
  verboseLS &&
    @printf("   Successful        dN        ratio\n")
  verboseLS &&
    @printf("%4d %9.2e  %9.2e %9.2e %9.2e %9.2e %9.2e\n", iter, state.x, φt, dφt, H,
                                                          f_meta.Δn, f_meta.Δp);

  while !OK #admissible: Armijo & Wolfe, tired: iterations
      # We select our descent direction. We compute the classical Newton
      # direction and then we check we the bounds of our trust region
      dN = -dφt / H
      d = TR_ls_step_computation(H, dφt, dN, f_meta.Δn, f_meta.Δp)

      # We see where the direction d gets us
      candidate_state = copy(state)
      update!(candidate_state, x = state.x + d)
      φtestTR, dφtestTR = φ_dφ(h, candidate_state; kwargs...)

      # depending on the approximation of the second derivate we need different
      # information
      if f_meta.dir == :Sec || f_meta.dir == :SecA
        tprec = state.x
        dφtprec = dφt
      end
      if f_meta.dir == :SecA φtprec = φt end

      # We compute the predicted and actual reduction as well as the ratio
      # between the two.
      (pred, ared, ratio) = pred_ared_computation(dφt, φt, H, d, φtestTR,
                                                  dφtestTR)
      # mettre le if qui suit dans une fonction pour aléger le code.
      if ratio < f_meta.eps1  # Unsuccessful approximation => reduce interval
        f_meta.Δp = f_meta.red * f_meta.Δp
        f_meta.Δn = f_meta.red * f_meta.Δn
        iter += 1
        verboseLS && @printf("%4d %9.2e %9.2e  %9.2e %9.2e %9.2e %9.2e %9e %9.2e %9e\n",
                              iter, state.x, φt, dφt, H, f_meta.Δn, f_meta.Δp, 0, dN, ratio);
      else             # Successful approximation => we move towards an
                       # admissible step size
        stop_ls.current_state = copy(candidate_state)
        state = stop_ls.current_state
        φt, dφt = φtestTR, dφtestTR
        H = update_H(f_meta.dir, h, state.x, tprec, φt, dφt, φtprec, dφtprec)

        if symmetrical        # We adjust our interval if we want it symmetrical
                              # or not
          if ratio > f_meta.eps2 # Very Successful iteration => augment size of
                                 # trust region
            f_meta.Δp = f_meta.aug * f_meta.Δp
            f_meta.Δn = f_meta.aug * f_meta.Δn
          end
        else
          if ratio > f_meta.eps2  # Very Successful iteration=> augment size of
                                  # trust region
            f_meta.Δp = f_meta.aug * f_meta.Δp
            f_meta.Δn = min(-state.x, f_meta.Δn)
          else
            f_meta.Δn = min(-state.x, f_meta.Δn)
          end
        end # à mettre dans une fonction

          iter += 1 # à mettre dans ls_at_t a un certain point

          verboseLS &&
            @printf("%4d %9.2e %9.2e  %9.2e %9.2e %9.2e %9.2e %9e %9.2e  %9e\n",
                    iter, state.x, φt, dφt, H, f_meta.Δn, f_meta.Δp, 1, dN, ratio);
      end;
      OK = stop!(stop_ls)
  end;

  return (state, stop_ls.meta.optimal)

end # function


text:
# Chapter 4.3

# CDF of a Gaussian mixture

## Julia code to generate the PDF and CDF
using Distributions
using Plots

gmm = MixtureModel(Normal[Normal(0, 1), Normal(5, 1),], [0.3, 0.7])

p1 = plot(x -> 0, -5, 1, fillrange= x->pdf(gmm, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false) # style
plot!(p1, x -> pdf(gmm, x), -5, 10,
    color=1, linewidth=6, label="PDF") # style

p2 = plot(x -> cdf(gmm, x), -5, 10,
    linewidth=6, label="CDF")

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

# CDF of a uniform random variable

# Julia code to generate the PDF and CDF of a uniform random variable
using Distributions
using Plots

u = Uniform(-3, 4)

# pdf
p1 = plot(x -> 0, -5:0.01:1, fillrange=x->pdf(u, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false)
plot!(p1, x -> pdf(u, x), -5, 10,
    color=1, linewidth=6, ylims=(0, 0.4), label="PDF")

# cdf
p2 = plot(x -> cdf(u, x), -5:0.01:10,
    linewidth=6, label="CDF")
vline!(p2, [-3, 4],
    linestyle=:dash, color=:green, label=false)

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

# CDF of an exponential random variable

# Julia code to generate the PDF and CDF of an exponential random variable
using Distributions
using Plots

u = Exponential(2)

# pdf
p1 = plot(x -> 0, -5:0.01:1, fillrange=x->pdf(u, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false)
plot!(p1, x -> pdf(u, x), -5, 10,
    color=1, linewidth=6, ylims=(0, 0.6), label="PDF")

# cdf
p2 = plot(x -> cdf(u, x), -5:0.01:10,
    linewidth=6, label="CDF")

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

# Chapter 4.5

# Generate a uniform random variable

# Julia code to generate 1000 uniform random numbers
using Distributions
using Plots

u = Uniform(0, 1) # same as Uniform()
X = rand(u, 1000)
histogram(X)

# Mean, variance, median, mode of a uniform random variable

# Julia code to compute empirical mean, var, median, mode
using Distributions
u = Uniform(0, 1)
X = rand(u, 1000)

M = mean(X)
V = var(X)
Med = median(X)
Mod = mode(X)

# Probability of a uniform random variable

# Julia code to compute the probability P(0.2 < X < 0.3)
using Distributions

u = Uniform(0, 1)
F = cdf(u, 0.3) - cdf(u, 0.2)

# PDF of an exponential random variable

# Julia code to generate PDF and CDF of an exponential random variable
using Distributions
using Plots

exp1 = Exponential(1/2)
exp2 = Exponential(1/5)

# Plotting the pdfs
p1 = plot(x -> pdf(exp1, x), 0, 1,
    linewidth=4, linestyle=:dashdot, color=RGB(0,0.2,0.8), label="λ = 2")
plot!(p1, x -> pdf(exp2, x), 0, 1,
    linewidth=4, color=RGB(0.8,0.2,0), label="λ = 5")

# Plotting the cdfs
p2 = plot(x -> cdf(exp1, x), 0, 1,
    linewidth=4, linestyle=:dashdot, color=RGB(0,0.2,0.8), label="λ = 2", legend=:topleft)
plot!(p2, x -> cdf(exp2, x), 0, 1,
    linewidth=4, color=RGB(0.8,0.2,0), label="λ = 5")

# Figure with both plots
plot(p1, p2, layout=(1, 2),
    size=(600, 400))

# Chapter 4.6

# PDF and CDF of a Gaussian random variable

# Julia code to generate standard Gaussian PDF and CDF
using Distributions
using Plots

normdist = Normal()

p1 = plot(x -> 0, -5:0.01:-1, fillrange=x->pdf(normdist, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false)
plot!(p1, x -> pdf(normdist, x), -5, 5,
    color=1, linewidth=6, ylims=(0, 0.6), label="PDF")

p2 = plot(x -> cdf(normdist, x), -5:0.01:5,
    linewidth=6, label="CDF")

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

# Skewness and kurtosis of a random variable

# Julia code to plot a Gamma distribution
using Distributions
using Plots

θ = 1

p1 = plot(size=(600, 300))
plot!(p1, x -> pdf(Gamma(2, θ), x), 0, 30, lw=4, c=RGB(0,0,0), label="k = 2")
plot!(p1, x -> pdf(Gamma(5, θ), x), 0, 30, lw=4, c=RGB(0.2,0.2,0.2), label="k = 5")
plot!(p1, x -> pdf(Gamma(10, θ), x), 0, 30, lw=4, c=RGB(0.4,0.4,0.4), label="k = 10")
plot!(p1, x -> pdf(Gamma(15, θ), x), 0, 30, lw=4, c=RGB(0.6,0.6,0.6), label="k = 15")
plot!(p1, x -> pdf(Gamma(20, θ), x), 0, 30, lw=4, c=RGB(0.8,0.8,0.8), label="k = 20")

# Julia code to compute skewness and kurtosis
X = rand(Gamma(3, 5), 10_000)
s = skewness(X)
k = kurtosis(X)

# Chapter 4.8

# Generating Gaussians from uniform

# Julia code to generate Gaussian from uniform
using Distributions
using Plots

mu = 3
sigma = 2
U = rand(10_000)
gU = sigma * quantile(Normal(), U) .+ mu

p1 = histogram(U, label="U")
p2 = histogram(gU, label="gU")
plot(p1, p2, layout=(1, 2), size=(600, 300), legend=:topleft)


text:
function viewVort2D(tev::Matrix{Float64}, lev::Matrix{Float64}, bv::Matrix{Float64})
    scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none")
    sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none")
    sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none")
    #plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0)
end

#For quick views while debugging
function viewVort2D(tev::Vector{TwoDVort}, lev::Vector{TwoDVort}, bv::Vector{TwoDVort})
    tev = hcat(map(q->q.s, tev), map(q->q.x, tev), map(q->q.z, tev))
    lev = hcat(map(q->q.s, lev), map(q->q.x, lev), map(q->q.z, lev))
    bv = hcat(map(q->q.s, bv), map(q->q.x, bv), map(q->q.z, bv))

    scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none")
    sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none")
    sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none")
    #plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0)
end

function viewVortConnect2D(tev::Matrix{Float64}, lev::Matrix{Float64}, bv::Matrix{Float64})
    scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none")
    sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none")
    sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none")

    plot(tev[:,2], tev[:,3], color = "red", linewidth=1.0)
    plot(lev[:,2], lev[:,3], color = "blue", linewidth=1.0)
    plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0)
end

function makeVortPlots2D()

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
    #Determine axis limits
    dirmax = maximum(dirresults)
    cd("$(dirmax)")

    multsurfflag = 0
    if isfile("boundv-1") == true  #only 1 surface
        multsurfflag = 1
    end

    if multsurfflag == 0 #single surface

        tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
        lev =   try
            DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
        catch
            zeros(0,3)
        end
        bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true)

        xmin = minimum([tev[:,2];lev[:,2];bv[:,2];])
        zmin = minimum([tev[:,3];lev[:,3];bv[:,3];])
        xmax = maximum([tev[:,2];lev[:,2];])
        zmax = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];])

        cd("..")

        if "vortPlots" in dirvec
            rm("vortPlots", recursive=true)
        end
        mkdir("vortPlots")
        for i=1:length(dirresults)
            if dirresults[i] != 0
                dirstr="$(dirresults[i])"
                cd(dirstr)
                tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
                lev =   try
                    DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
                catch
                    zeros(0,3)
                end
                bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true)

                viewVort2D(tev, lev, bv)
                axis([xmin-1, xmax+1, zmin-1, zmax+1])
                savefig("../vortPlots/$(dirresults[i]).png")
                close()
                cd("..")
            end
        end

    else #multiple surfaces

        tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
        lev =   try
            DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
        catch
            zeros(0,3)
        end

        nsurf = 0
        for i = 1:1000
            if isfile("boundv-$i") == true
                nsurf += 1

            end
        end
        bv, _ = DelimitedFiles.readdlm("boundv-1", '\t', Float64, header=true)
        for i = 2:nsurf
            bvt, _ = DelimitedFiles.readdlm("boundv-$i", '\t', Float64, header=true)
            bv = [bv;bvt;]
        end

        xmin = minimum([tev[:,2];lev[:,2];bv[:,2];])
        zmin = minimum([tev[:,3];lev[:,3];bv[:,3];])
        xmax = maximum([tev[:,2];lev[:,2];])
        zmax = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];])

        cd("..")

        if "vortPlots" in dirvec
            rm("vortPlots", recursive=true)
        end
        mkdir("vortPlots")
        for i=1:length(dirresults)
            if dirresults[i] != 0
                dirstr="$(dirresults[i])"
                cd(dirstr)
                tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
                lev =   try
                    DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
                catch
                    zeros(0,3)
                end
                bv, _ = DelimitedFiles.readdlm("boundv-1", '\t', Float64, header=true)
                for i = 2:nsurf
                    bvt, _ = DelimitedFiles.readdlm("boundv-$i", '\t', Float64, header=true)
                    bv = [bv;bvt;]
                end
                viewVort2D(tev, lev, bv)
                axis([xmin-1, xmax+1, zmin-1, zmax+1])
                savefig("../vortPlots/$(dirresults[i]).png")
                close()
                cd("..")
            end
        end
    end
end

function makeForcePlots2D()
    dirvec = readdir()
    if "forcePlots" in dirvec
        rm("forcePlots", recursive=true)
    end
    mkdir("forcePlots")

    mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true)

    if length(mat[1,:]) == 8 #only 1 surface
        t = mat[:,1]
        alpha = mat[:,2]*180/pi
        h = mat[:,3]
        cl = mat[:,6]
        cd = mat[:,7]
        cm = mat[:,8]

        plot(t, alpha)
        len = length(t)
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(alpha) - 0.1*abs(minimum(alpha))
        zmax = maximum(alpha) + 0.1*abs(maximum(alpha))
        axis([xmin, xmax, zmin, zmax])
        xlabel(L"$t^*$")
        ylabel(L"$\alpha$ (deg)")
        savefig("forcePlots/pitch.png")
        close()

        plot(t, h)
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(h) - 0.1*abs(minimum(h))
        zmax = maximum(h) + 0.1*abs(maximum(h))
        axis([xmin, xmax, zmin, zmax])
        xlabel(L"$t^*$")
        ylabel(L"$h/c$")
        savefig("forcePlots/plunge.png")
        close()

        plot(t, cl)
        range = Int(round(0.05*len)):Int(round(0.95*len))
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
        zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
        axis([xmin, xmax, zmin, zmax])
        xlabel(L"$t^*$")
        ylabel(L"$C_l$")
        savefig("forcePlots/cl.png")
        close()

        plot(t, cd)
        range = Int(round(0.05*len)):Int(round(0.95*len))
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(cd[range]) - 0.1*abs(minimum(cd[range]))
        zmax = maximum(cd[range]) + 0.1*abs(maximum(cd[range]))
        axis([xmin, xmax, zmin, zmax])
        xlabel(L"$t^*$")
        ylabel(L"$C_d$")
        savefig("forcePlots/cd.png")
        close()

        plot(t, cm)
        range = Int(round(0.05*len)):Int(round(0.95*len))
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(cm[range]) - 0.1*abs(minimum(cm[range]))
        zmax = maximum(cm[range]) + 0.1*abs(maximum(cm[range]))
        axis([xmin, xmax, zmin, zmax])
        xlabel(L"$t^*$")
        ylabel(L"$C_m$")
        savefig("forcePlots/cm.png")
        close()

    else #multiple surfaces

        nsurf = (length(mat[1,:]) - 1)/7
        t = mat[:,1]

        for i = 1:nsurf
            alpha = mat[:,Int((i-1)*7+2)]*180/pi
            h = mat[:,Int((i-1)*7+3)]
            cl = mat[:,Int((i-1)*7+6)]
            cd = mat[:,Int((i-1)*7+7)]
            cm = mat[:,Int((i-1)*7+8)]

            plot(t, alpha)
            len = length(t)
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha) - 0.1*abs(minimum(alpha))
            zmax = maximum(alpha) + 0.1*abs(maximum(alpha))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$\alpha$ (deg)")
            savefig("forcePlots/pitch-$i.png")
            close()

            plot(t, h)
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h) - 0.1*abs(minimum(h))
            zmax = maximum(h) + 0.1*abs(maximum(h))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$h/c$")
            savefig("forcePlots/plunge-$i.png")
            close()

            plot(t, cl)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
            zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$C_l$")
            savefig("forcePlots/cl-$i.png")
            close()

            plot(t, cd)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cd[range]) - 0.1*abs(minimum(cd[range]))
            zmax = maximum(cd[range]) + 0.1*abs(maximum(cd[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$C_d$")
            savefig("forcePlots/cd-$i.png")
            close()

            plot(t, cm)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cm[range]) - 0.1*abs(minimum(cm[range]))
            zmax = maximum(cm[range]) + 0.1*abs(maximum(cm[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$C_m$")
            savefig("forcePlots/cm-$i.png")
            close()
        end

    end

end

function checkConverge(k::Float64)

    mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true)

    T = pi/k

    end_cycle = mat[end,1]

    #Find number of cycles
    ncyc = 0
    for i = 1:1000
        t_l = real(i)*T
        if t_l > end_cycle
            break
        end
        ncyc = ncyc + 1
    end

    #Lift convergence
    figure
    for i = 1:ncyc
        start_t = real(i-1)*T
        end_t = real(i)*T
        start_ind = argmin(abs.(mat[:,1] .- start_t))
        end_ind = argmin(abs.(mat[:,1] .- end_t))
        plot((mat[start_ind:end_ind,1] .- start_t)/T, mat[start_ind:end_ind,6])
    end
    xmin = 0.
    xmax = 1.
    zmin = minimum(mat[:,6])
    zmax = maximum(mat[:,6])
    axis([xmin, xmax, zmin, zmax])

    xlabel(L"$t^*$")
    ylabel(L"$C_l$")
    savefig("forcePlots/cl-convergence.png")
    close()

    #Drag convergence
    figure
    for i = 1:ncyc
        start_t = real(i-1)*T
        end_t = real(i)*T
        start_ind = argmin(abs.(mat[:,1] .- start_t))
        end_ind = argmin(abs.(mat[:,1] .- end_t))
        plot((mat[start_ind:end_ind,1] .- start_t)/T, mat[start_ind:end_ind,7])
    end
    xmin = 0.
    xmax = 1.
    zmin = minimum(mat[:,7])
    zmax = maximum(mat[:,7])
    axis([xmin, xmax, zmin, zmax])

    xlabel(L"$t^*$")
    ylabel(L"$C_d$")
    savefig("forcePlots/cd-convergence.png")
    close()

    #Pitching moment convergence
    figure
    for i = 1:ncyc
        start_t = real(i-1)*T
        end_t = real(i)*T
        start_ind = argmin(abs.(mat[:,1] .- start_t))
        end_ind = argmin(abs.(mat[:,1] .- end_t))
        plot((mat[start_ind:end_ind,1] .- start_t)/T, mat[start_ind:end_ind,8])
    end
    xmin = 0.
    xmax = 1.
    zmin = minimum(mat[:,8])
    zmax = maximum(mat[:,8])
    axis([xmin, xmax, zmin, zmax])

    xlabel(L"$t^*$")
    ylabel(L"$C_m$")
    savefig("forcePlots/cm-convergence.png")
    close()
end

function makeKinemClVortPlots2D()

    mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true)

    t = mat[:,1]
    len = length(t)
    cl = mat[:,6]
    alpha = mat[:,2]*180/pi
    h = mat[:,3]
    u = mat[:,4]

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)

    #Determine axis limits
    dirmax = maximum(dirresults)

    cd("$(dirmax)")

    multsurfflag = 0
    if isfile("boundv-1") == true  #only 1 surface
        error("this plot function is only written for single surface")
    end

    tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
    lev =   try
        DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
    catch
        zeros(0,3)
    end
    bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true)

    xminv = minimum([tev[:,2];lev[:,2];bv[:,2];])
    zminv = minimum([tev[:,3];lev[:,3];bv[:,3];])
    xmaxv = maximum([tev[:,2];lev[:,2];])
    zmaxv = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];])

    cd("..")

    if "infoPlots" in dirvec
        rm("infoPlots", recursive=true)
    end
    mkdir("infoPlots")

    for i=1:length(dirresults)
        if dirresults[i] != 0
            dirstr="$(dirresults[i])"
            cd(dirstr)

            t_cur = dirresults[i]

            figure(figsize=(8,8))

            subplot2grid((6, 2), (0, 0))
            plot(t, alpha)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha[range]) - 0.1*abs(minimum(alpha[range]))
            zmax = maximum(alpha[range]) + 0.1*abs(maximum(alpha[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel(L"$t^*$")
            xticks([])
            ylabel(L"$\alpha (deg)$")

            subplot2grid((6, 2), (1, 0))
            plot(t, h)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h[range]) - 0.1*abs(minimum(h[range]))
            zmax = maximum(h[range]) + 0.1*abs(maximum(h[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel(L"$t^*$")
            xticks([])
            ylabel(L"$h$")

            subplot2grid((6, 2), (2, 0))
            plot(t, u)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(u[range]) - 0.1*abs(minimum(u[range]))
            zmax = maximum(u[range]) + 0.1*abs(maximum(u[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$u$")

            subplot2grid((6, 2), (0, 1), rowspan=3)
            plot(t, cl)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
            zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$C_l$")

            subplot2grid((6, 2), (3, 0), rowspan=3, colspan=2)
            tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])

            tight_layout()

            savefig("../infoPlots/$(dirresults[i]).png")
            close()
            cd("..")
        end
    end
end

function makeKinemVelVortPlots2D()

    mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true)

    t = mat[:,1]
    len = length(t)
    cl = mat[:,6]
    alpha = mat[:,2]*180/pi
    h = mat[:,3]
    u = mat[:,4]

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)

    #Determine axis limits
    dirmax = maximum(dirresults)

    cd("$(dirmax)")

    multsurfflag = 0
    if isfile("boundv-1") == true  #only 1 surface
        error("this plot function is only written for single surface")
    end

    tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
    lev =   try
        DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
    catch
        zeros(0,3)
    end
    bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true)

    xminv = minimum([tev[:,2];lev[:,2];bv[:,2];])
    zminv = minimum([tev[:,3];lev[:,3];bv[:,3];])
    xmaxv = maximum([tev[:,2];lev[:,2];])
    zmaxv = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];])

    cd("..")

    if "infoPlots" in dirvec
        rm("infoPlots", recursive=true)
    end
    mkdir("infoPlots")

    for i=1:length(dirresults)
        if dirresults[i] != 0
            dirstr="$(dirresults[i])"
            cd(dirstr)

            t_cur = dirresults[i]

            figure(figsize=(8,8))

            subplot2grid((6, 2), (0, 0))
            plot(t, alpha)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha[range]) - 0.1*abs(minimum(alpha[range]))
            zmax = maximum(alpha[range]) + 0.1*abs(maximum(alpha[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel(L"$t^*$")
            xticks([])
            ylabel(L"$\alpha (deg)$")

            subplot2grid((6, 2), (1, 0))
            plot(t, h)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h[range]) - 0.1*abs(minimum(h[range]))
            zmax = maximum(h[range]) + 0.1*abs(maximum(h[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel(L"$t^*$")
            xticks([])
            ylabel(L"$h$")

            subplot2grid((6, 2), (2, 0))
            plot(t, u)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(u[range]) - 0.1*abs(minimum(u[range]))
            zmax = maximum(u[range]) + 0.1*abs(maximum(u[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$u$")

            subplot2grid((6, 2), (0, 1), rowspan=3)
            plot(t, cl)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
            zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$C_l$")

            subplot2grid((6, 2), (3, 0), rowspan=3, colspan=2)
            tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])

            tight_layout()

            savefig("../infoPlots/$(dirresults[i]).png")
            close()
            cd("..")
        end
    end
end

max_stars_repo_name: hjabird/UNSflow

max_issues_repo_name: hjabird/UNSflow

max_forks_repo_name: hjabird/UNSflow


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# runtime glass cats
const MIL_GLASSES = Dict{Int,Glass}()
const MODEL_GLASSES = Vector{Glass}(undef, 0)

"""
    glassfromMIL(glasscode::Union{Float64,Int}) -> Glass

Generates a glass object for the given glass code based on U.S. military standard MIL-G-174, see [the MIL specification](http://www.newportglass.com/GeneCd.htm) for further details.

The glass code is a six-digit number specifying the glass according to its refractive index `Nd` at d-light (587.5618nm), and its Abbe number `Vd` also taken at d-light.
The resulting glass code is the value of `Nd - 1` rounded to three digits, followed by `Vd` rounded to three digits, with all decimal points ignored.
For example, N_BK7 has `Nd = 1.5168` and `Vd = 64.17`, giving a six-digit glass code of `517642`.

For `Nd > 1.999` the format `1.123642` can be used representing `Nd = 2.123` and `Vd = 64.2`.

**Accuracy is poor given the low precision of the input parameters**, the mean error to measured data may be significant.
Behavior may differ from other optical simulation tools when using MIL glasses.
The approximate dispersion calculation used these glasses is generally only valid for visible wavelengths, in this case a limit of 360nm to 750nm is imposed.

# Examples
```julia-repl
julia> index(glassfromMIL(517642), 0.5875618)
1.5170003960064509

julia> index(glassfromMIL(1.134642), 0.5875618)
2.1340008686098946
```
"""
function glassfromMIL(glasscode::Int)
    Nd = floor(Int, glasscode / 1000) / 1000 + 1
    Vd = (glasscode - floor(Int, glasscode / 1000) * 1000) / 10
    g = _modelglass(Nd, Vd, 0.0, GlassID(MIL, glasscode))
    MIL_GLASSES[glasscode] = g
    return g
end

function glassfromMIL(glasscode::Float64)
    @assert glasscode > 1.0
    glasscodeid = round(Int, glasscode * 1000000)
    Nd = floor(Int, glasscode * 1000) / 1000 + 1
    Vd = round((glasscode * 1000 - floor(Int, glasscode * 1000)) * 100, digits = 1)
    g = _modelglass(Nd, Vd, 0.0, GlassID(MIL, glasscodeid))
    MIL_GLASSES[glasscodeid] = g
    return g
end

"""
    modelglass(Nd::Float64, Vd::Float64, ΔPgF::Float64) -> Glass

Generates a glass object for the given refractive index at d-light (587.5618nm), `Nd`, the Abbe number also at d-light, `Vd`, and partial dispersion, `ΔPgF`.
The mean error to measured data for these models is typically small - usually < 0.0001.
Behavior may differ from other optical simulation tools when using model glasses.

The approximate dispersion calculation used for these glasses is generally only valid for visible wavelengths, in this case a limit of 360nm to 750nm is imposed.

# Examples
```julia-repl
julia> index(modelglass(1.5168, 64.17, 0.0), 0.5875618)
1.5168003970108495

julia> index(modelglass(1.2344, 61.57, 0.003), 0.678)
1.2329425902693352
```
"""
function modelglass(Nd::Float64, Vd::Float64, ΔPgF::Float64)
    g = _modelglass(Nd, Vd, ΔPgF, GlassID(MODEL, length(MODEL_GLASSES) + 1))
    push!(MODEL_GLASSES, g)
    return g
end

function _modelglass(Nd::Float64, Vd::Float64, ΔPgF::Float64, ID::GlassID)
    # from Schott "TIE-29: Refractive Index and Dispersion"
    a = ΔPgF + 0.6438 - 0.001682 * Vd
    # Using fitting results from https://www.gnu.org/software/goptical/manual/Material_Abbe_class_reference.html
    C1 = a * -6.11873891971188577088 + 1.17752614766485175224
    C2 = a * 18.27315722388047447566 + -8.93204522498095698779
    C3 = a * -14.55275321129051135927 + 7.91015964461522003148
    C4 = a * 3.48385106908642905310 + -1.80321117937358499361
    return Glass(ID, -1, C1, C2, C3, C4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.36, 0.75, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, TEMP_REF, ΔPgF, -1.0, -1.0, 0.0, -1.0, 0, -1.0, nothing, Nd, -1.0, -1.0, 0, Vd, 1, 0.0, 0.0)
end


text:
using CellularAutomata

const rule0 = DCA(0)
@test rule0.ruleset == [0, 0, 0, 0, 0, 0, 0, 0] #http://atlas.wolfram.com/01/01/0/

const rule1 = DCA(1)
@test rule1.ruleset == [1, 0, 0, 0, 0, 0, 0, 0] #http://atlas.wolfram.com/01/01/1/

const rule2 = DCA(2)
@test rule2.ruleset == [0, 1, 0, 0, 0, 0, 0, 0] #http://atlas.wolfram.com/01/01/2/

const rule3 = DCA(3)
@test rule3.ruleset == [1, 1, 0, 0, 0, 0, 0, 0] #http://atlas.wolfram.com/01/01/2/

const rule30 = DCA(30)
@test rule30.ruleset == [0, 1, 1, 1, 1, 0, 0, 0] #http://atlas.wolfram.com/01/01/30/


text:
struct SNumber
    digits::Vector{Float32}
    point::Int16
    base::Float32

    function SNumber(digits::Vector{Float32}, point::Int16; base=10.0, num_digits=10)
        new(digits, point, base)
    end

    function SNumber(v::Number; base::Float64=10.0, num_digits=10)
        if v == 0
            return new([0], 0, base)
        end
        digits = zeros(Float32, num_digits)
        point::Int16 = 0
        start_digit = Int32(floor(log(base, v))) + 1
        point = start_digit - 1

        counter = 1
        for i in start_digit:-1:(start_digit - num_digits + 1)
            b = base ^ (i - point)
            d = v / b
            res = isapprox(d, ceil(d)) ? ceil(d) : floor(d)
            if res > 0
                v -= b * res
                digits[counter] = res
            end
            counter += 1
        end
        SNumber(digits, point, base=base, num_digits=num_digits)
    end

    Base.:+(n::SNumber, i::Number) = n + SNumber(i)
    Base.getindex(n::SNumber, i::Integer) = i > length(n.digits) || i <= 0 ? 0 : n.digits[i]
    digits(n::SNumber) = digits
    Base.length(n::SNumber) = length(n.digits)
    Base.print(io::IO, n::SNumber) = Base.string(n)

    function Base.string(n::SNumber)
        str = ""
        for i in 1:n.point
            str *= i != n.point ? string(n.digits[i]) * ", " : string(n.digits[i])
        end
        str *= ":"
        for i in (n.point + 1):length(n)
                str *= i != length(n) ? string(n.digits[i]) * ", " : string(n.digits[i])
        end
        str * "    Base:" * string(n.base)
    end

    function Base.:+(n::SNumber, n2::SNumber)
        b = max(n.base, n2.base)
        SNumber(basesum(n, b) + basesum(n2, b), base=n.base, num_digits=length(n))
    end

end

function basesum(n::SNumber; base_num::Float64=base)
    sum::Float64 = 0
    for d in 1:length(n)
        sum += n[d] * base_num ^ (d - n.point)
    end
    sum
end


text:
import Base.position

export SPK, position, velocity, state, segments, print_segments

const SECONDS_PER_DAY = 86400
const SIZE_FLOAT64 = sizeof(Float64)

type OutOfRangeError <: Exception
    date::Float64
    startdate::Float64
    finaldate::Float64
end

Base.showerror(io::IO, err::OutOfRangeError) = print(io, "The requested date $(err.date) is outside the intervall ($(err.startdate), $(err.finaldate)).")

type Segment
    name::ASCIIString
    firstsec::Float64
    lastsec::Float64
    firstdate::Float64
    lastdate::Float64
    target::Int32
    center::Int32
    frame::Int32
    spktype::Int32
    firstaddr::Int32
    lastaddr::Int32
    firstword::Int32
    lastword::Int32
    initialsecond::Float64
    intlen::Float64
    rsize::Int32
    n_records::Int32
    cache::Matrix{Float64}
    cached_record::Int
end

jd(sec) = 2451545 + sec/SECONDS_PER_DAY
seconds(jd) = (jd - 2451545)*SECONDS_PER_DAY

function Segment(daf, name, record)
    firstsec, lastsec = reinterpret(Float64, record[1:16], daf.little)
    target, center, frame, spktype, firstaddr, lastaddr = reinterpret(Int32, record[17:end], daf.little)
    if spktype != 2
        error("Type $spktype SPK file detected. Only Type 2 SPK files are supported.")
    end
    init, intlen, rsize, n_records = reinterpret(Float64, daf.array[lastaddr*SIZE_FLOAT64-4*SIZE_FLOAT64+1:lastaddr*SIZE_FLOAT64], daf.little)
    n_records = round(Int32, n_records)
    Segment(
        name,
        firstsec,
        lastsec,
        jd(firstsec),
        jd(lastsec),
        target,
        center,
        frame,
        spktype,
        firstaddr,
        lastaddr,
        firstaddr*SIZE_FLOAT64 - SIZE_FLOAT64 + 1,
        lastaddr*SIZE_FLOAT64 - SIZE_FLOAT64*4,
        init,
        intlen,
        round(Int32, rsize),
        n_records,
        Matrix{Float64}(),
        -1,
    )
end

type SPK
    daf::DAF
    segments::Dict{Int, Dict{Int, Segment}}
end

function SPK(filename)
    daf = DAF(filename)
    segments = Dict{Int, Dict{Int, Segment}}()
    for (name, summary) in getsummaries(daf)
        seg = Segment(daf, name, summary)
        if haskey(segments, seg.center)
            merge!(segments[seg.center], Dict(seg.target=>seg))
        else
            merge!(segments, Dict(seg.center=>Dict(seg.target=>seg)))
        end
    end
    SPK(daf, segments)
end

segments(spk::SPK) = spk.segments

function list_segments(spk::SPK)
    s = ASCIIString[]
    for (k,v) in spk.segments
        for l in keys(v)
            push!(s, "$(name_from_naifid(k)) ($k) => $(name_from_naifid(l)) ($l)")
        end
    end
    return sort!(s, lt=segstrlt)
end

function print_segments(spk::SPK)
    s = list_segments(spk)
    println(join(s, "\n"))
end

function segstrlt(a::ASCIIString, b::ASCIIString)
   rex = r"\([0-9]*\)$"
   ma = match(rex, a)
   mb = match(rex, b)
   ia = parse(Int, a[ma.offset+1:end-1])
   ib = parse(Int, b[mb.offset+1:end-1])
   return ia < ib
end

function checkdate(seg::Segment, tdb::Float64)
    if !(seg.firstdate <= tdb <= seg.lastdate)
        throw(OutOfRangeError(tdb, seg.firstdate, seg.lastdate))
    end
end

function getcoefficients(spk::SPK, seg::Segment, tdb::Float64, tdb2::Float64=0.0)
    checkdate(seg, tdb+tdb2)
    components = 3
    order = (seg.rsize - 2) ÷ 3
    secs = (seconds(tdb) - seg.initialsecond) + tdb2*SECONDS_PER_DAY
    recordnum, frac = divrem(secs, seg.intlen)
    recordnum = round(Int, recordnum)
    if recordnum == seg.n_records
        recordnum -= 1
        frac = seg.intlen
    end
    if recordnum == seg.cached_record
        c = seg.cache
    else
        # Drop the MID and RADIUS values
        first = seg.firstword + SIZE_FLOAT64*seg.rsize*recordnum + SIZE_FLOAT64*2
        last = seg.firstword + SIZE_FLOAT64*seg.rsize*(recordnum+1)
        cv = reinterpret(Float64, spk.daf.array[first:last], spk.daf.little)
        c = reshape(cv, (order, components))'
        seg.cache = c
        seg.cached_record = recordnum
    end
    x = Array(Float64, order)
    tc = 2.0 * frac/seg.intlen - 1.0
    x[1] = 1.0
    x[2] = tc
    twotc = tc + tc
    @inbounds for i = 3:order
        x[i] = twotc*x[i-1] - x[i-2]
    end
    c, x, seg.intlen, twotc, order
end

function position(c::Matrix, x::Vector, order::Int)
    r = zeros(3)
    @inbounds @simd for i = 1:3
        for j = 1:order
            r[i] += c[i, j] * x[j]
        end
    end
    return r
end

function position(spk::SPK, seg::Segment, tdb::Float64, tdb2::Float64=0.0)
    c, x, dt, twotc, order = getcoefficients(spk, seg, tdb, tdb2)
    position(c, x, order)
end

function velocity(c::Matrix, x::Vector, dt::Float64, twotc::Float64, order::Int)
    v = zeros(Float64, 3)
    t = zeros(Float64, order)
    t[2] = 1.0
    if order > 2
        t[3] = twotc + twotc
        for i = 4:order
            t[i] = twotc*t[i-1] - t[i-2] + x[i-1] + x[i-1]
        end
    end
    t *= 2.0
    t /= dt
    @inbounds @simd for i = 1:3
        for j = 1:order
            v[i] += c[i, j] * t[j]
        end
    end
    return v
end

function velocity(spk::SPK, seg::Segment, tdb::Float64, tdb2::Float64=0.0)
    c, x, dt, twotc, order = getcoefficients(spk, seg, tdb, tdb2)
    velocity(c, x, dt, twotc, order)
end

function state(spk::SPK, seg::Segment, tdb::Float64, tdb2::Float64=0.0)
    c, x, dt, twotc, order = getcoefficients(spk, seg, tdb, tdb2)
    r = position(c, x, order)
    v = velocity(c, x, dt, twotc, order)
    [r;v]
end

for (f, n) in zip((:state, :velocity, :position), (6, 3, 3))
    @eval begin
        function ($f)(spk::SPK, center::Int, target::Int, tdb::Float64, tdb2::Float64=0.0)
            seg = spk.segments[center][target]
            ($f)(spk, seg, tdb, tdb2)
        end

        function ($f)(spk::SPK, target::Int, tdb::Float64, tdb2::Float64=0.0)
            seg = spk.segments[0][target]
            ($f)(spk, seg, tdb, tdb2)
        end

        function ($f)(spk::SPK, target::AbstractString, tdb::Float64, tdb2::Float64=0.0)
            ($f)(spk, naifid(target), tdb, tdb2)
        end

        function ($f)(spk::SPK, center::AbstractString, target::AbstractString, tdb::Float64, tdb2::Float64=0.0)
            ($f)(spk, naifid(center), naifid(target), tdb, tdb2)
        end

        function ($f)(spk::SPK, center::Int, target::AbstractString, tdb::Float64, tdb2::Float64=0.0)
            ($f)(spk, center, naifid(target), tdb, tdb2)
        end

        function ($f)(spk::SPK, center::AbstractString, target::Int, tdb::Float64, tdb2::Float64=0.0)
            ($f)(spk, naifid(center), target, tdb, tdb2)
        end

        function ($f)(spk::SPK, center, target, tdb::AbstractArray, tdb2::AbstractArray=zeros(length(tdb)))
            m = length(tdb)
            if m != length(tdb2)
                error("'tdb' and 'tdb2' must have the same length.")
            end
            out = zeros(m, $n)
            for (i, t, t2) in zip(1:m, tdb, tdb2)
                out[i,:] = ($f)(spk, center, target, t, t2)
            end
            return out
        end

        function ($f)(spk::SPK, target, tdb::AbstractArray, tdb2::AbstractArray=zeros(length(tdb)))
            ($f)(spk, 0, target, tdb, tdb2)
        end
    end
end


text:
#=
Written by Rasmus Kyng.

Given an adjacency matrix mat
where A[i,j] != 0 indicates a directed edge from i to j,
compute a topological sorting of the vertices, s.t.

Please check that the edge direction agrees with what you want before usage!
=#

function toposort{Tv,Ti}(mat::SparseMatrixCSC{Tv,Ti})
    # each COL i gives the outgoing edges of vertex i
    n = mat.n
    colptr::Array{Ti,1} = mat.colptr
    rowval::Array{Ti,1} = mat.rowval # the row index corresponding to a given colptr entry
    popSorted::Array{Ti,1} = zeros(Ti,n);
    seen::Array{Bool,1} = falses(n);
    popCount::Ti = 0
    s = DataStructures.Stack(Tuple{Ti,Ti})
    for i = 1:n
        if seen[i]
            continue
        end
        seen[i] = true
        DataStructures.push!(s,(i,colptr[i]))
        while !isempty(s)
            vertNbrTuple = DataStructures.pop!(s)
            v = vertNbrTuple[1]
            edgeNum = vertNbrTuple[2]
            if edgeNum < colptr[v+1]
                DataStructures.push!(s,(v,edgeNum+1))
                nbr = rowval[edgeNum]
                if !seen[nbr]
                    seen[nbr] = true
                    DataStructures.push!(s,(nbr,colptr[nbr]))
                end
            else
                popCount += 1
                popSorted[popCount] = v
            end
        end
    end
    return popSorted
end

"""The signed edge-vertex adjacency matrix"""
function dirEdgeVertexMat(A::SparseMatrixCSC)
    (ai,aj) = findnz(A)
    m = length(ai)
    n = size(A)[1]
    return sparse(collect(1:m),ai,1.0,m,n) - sparse(collect(1:m),aj,1.0,m,n)
end


text:
using AbstractTrees

mutable struct BinaryNode{T}
    data::T
    parent::Union{Nothing,BinaryNode{T}}
    left::Union{Nothing,BinaryNode{T}}
    right::Union{Nothing,BinaryNode{T}}

    function BinaryNode{T}(data, parent=nothing, l=nothing, r=nothing) where T
        new{T}(data, parent, l, r)
    end
end
BinaryNode(data) = BinaryNode{typeof(data)}(data)

function leftchild!(parent::BinaryNode, data)
    isnothing(parent.left) || error("left child is already assigned")
    node = typeof(parent)(data, parent)
    parent.left = node
end
function rightchild!(parent::BinaryNode, data)
    isnothing(parent.right) || error("right child is already assigned")
    node = typeof(parent)(data, parent)
    parent.right = node
end

## Things we need to define
function AbstractTrees.children(node::BinaryNode)
    if isnothing(node.left) && isnothing(node.right)
        ()
    elseif isnothing(node.left) && !isnothing(node.right)
        (node.right,)
    elseif !isnothing(node.left) && isnothing(node.right)
        (node.left,)
    else
        (node.left, node.right)
    end
end

AbstractTrees.nodevalue(n::BinaryNode) = n.data

AbstractTrees.ParentLinks(::Type{<:BinaryNode}) = StoredParents()

AbstractTrees.parent(n::BinaryNode) = n.parent

AbstractTrees.NodeType(::Type{<:BinaryNode{T}}) where {T} = HasNodeType()
AbstractTrees.nodetype(::Type{<:BinaryNode{T}}) where {T} = BinaryNode{T}

function binarynode_example()
    n₀ = BinaryNode(0)
    l₁ = leftchild!(n₀, 1)
    r₁ = rightchild!(n₀, 2)
    r₂ = rightchild!(l₁, 3)

    n₀
end