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

name = "libfacedetection"
# I guess it's v3 https://github.com/ShiqiYu/libfacedetection/blob/master/ChangeLog#L3
# but doesn't really have any versions
version = v"3"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/ShiqiYu/libfacedetection.git", "36783e346c596968f6648bdfdc76efc8c6f082a9"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/libfacedetection
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
	-DBUILD_SHARED_LIBS=ON \
	-DENABLE_AVX2=OFF \
	-DDEMO=OFF \
	..
make -j${nproc}
make install
install_license ../LICENSE
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = []

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module AlgebraicDecisionDiagrams

include("DecisionDiagrams.jl")
include("Operations.jl")
include("MultilinearExpressions.jl")

end # module

ext: jl


text:
# # Machine Learning in Julia

# An introduction to the
# [MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
# toolbox.

# ### Set-up

# Inspect Julia version:

VERSION

# The following instantiates a package environment.

# The package environment has been created using **Julia 1.6** and may not
# instantiate properly for other Julia versions.

using Pkg
Pkg.activate("env")
Pkg.instantiate()

# ## General resources

# - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
# - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
# - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
# - [Data Science Tutorials in Julia](https://juliaai.github.io/DataScienceTutorials.jl/)

# ## Part 1 - Data Representation

# > **Goals:**
# > 1. Learn how MLJ specifies it's data requirements using "scientific" types
# > 2. Understand the options for representing tabular data
# > 3. Learn how to inspect and fix the representation of data to meet MLJ requirements

# ### Scientific types

# To help you focus on the intended *purpose* or *interpretation* of
# data, MLJ models specify data requirements using *scientific types*,
# instead of machine types. An example of a scientific type is
# `OrderedFactor`. The other basic "scalar" scientific types are
# illustrated below:

# ![](scitypes.png)

# A scientific type is an ordinary Julia type (so it can be used for
# method dispatch, for example) but it usually has no instances. The
# `scitype` function is used to articulate MLJ's convention about how
# different machine types will be interpreted by MLJ models:

using ScientificTypes
scitype(3.141)

#-

time = [2.3, 4.5, 4.2, 1.8, 7.1]
scitype(time)

# To fix data which MLJ is interpreting incorrectly, we use the
# `coerce` method:

height = [185, 153, 163, 114, 180]
scitype(height)

#-

height = coerce(height, Continuous)

# Here's an example of data we would want interpreted as
# `OrderedFactor` but isn't:

exam_mark = ["rotten", "great", "bla",  missing, "great"]
scitype(exam_mark)

#-

exam_mark = coerce(exam_mark, OrderedFactor)

#-

levels(exam_mark)

# Use `levels!` to put the classes in the right order:

levels!(exam_mark, ["rotten", "bla", "great"])
exam_mark[1] < exam_mark[2]

# When sub-sampling, no levels are lost:

levels(exam_mark[1:2])

# **Note on binary data.** There is no separate scientific type for
# binary data. Binary data is `OrderedFactor{2}` or
# `Multiclass{2}`. If a binary measure like `truepositive` is a
# applied to `OrderedFactor{2}` then the "positive" class is assumed
# to appear *second* in the ordering. If such a measure is applied to
# `Multiclass{2}` data, a warning is issued. A single `OrderedFactor`
# can be coerced to a single `Continuous` variable, for models that
# require this, while a `Multiclass` variable can only be one-hot
# encoded.

# ### Two-dimensional data

# Whenever it makes sense, MLJ Models generally expect two-dimensional
# data to be *tabular*. All the tabular formats implementing the
# [Tables.jl API](https://juliadata.github.io/Tables.jl/stable/) (see
# this
# [list](https://github.com/JuliaData/Tables.jl/blob/master/INTEGRATIONS.md))
# have a scientific type of `Table` and can be used with such models.

# Probably the simplest example of a table is the julia native *column
# table*, which is just a named tuple of equal-length vectors:

column_table = (h=height, e=exam_mark, t=time)

#-

scitype(column_table)

#-

# Notice the `Table{K}` type parameter `K` encodes the scientific
# types of the columns. (This is useful when comparing table scitypes
# with `<:`). To inspect the individual column scitypes, we use the
# `schema` method instead:

schema(column_table)

# Here are five other examples of tables:

dict_table = Dict(:h => height, :e => exam_mark, :t => time)
schema(dict_table)

# (To control column order here, instead use `LittleDict` from
# OrderedCollections.jl.)

row_table = [(a=1, b=3.4),
             (a=2, b=4.5),
             (a=3, b=5.6)]
schema(row_table)

#-

import DataFrames
df = DataFrames.DataFrame(column_table)

#-

schema(df) == schema(column_table)

#-

using UrlDownload, CSV
csv_file = urldownload("https://raw.githubusercontent.com/ablaom/"*
                   "MachineLearningInJulia2020/"*
                   "for-MLJ-version-0.16/data/horse.csv");
schema(csv_file)

# Most MLJ models do not accept matrix in lieu of a table, but you can
# wrap a matrix as a table:

using Tables
matrix_table = Tables.table(rand(2,3))
schema(matrix_table)

# The matrix is *not* copied, only wrapped. Some models may perform
# better if one wraps the adjoint of the transpose - see
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Observations-correspond-to-rows,-not-columns).

# **Manipulating tabular data.** In this workshop we assume
# familiarity with some kind of tabular data container (although it is
# possible, in principle, to carry out the exercises without this.)
# For a quick start introduction to `DataFrames`, see [this
# tutorial](https://juliaai.github.io/DataScienceTutorials.jl/data/dataframe/).

# ### Fixing scientific types in tabular data

# To show how we can correct the scientific types of data in tables,
# we introduce a cleaned up version of the UCI Horse Colic Data Set
# (the cleaning work-flow is described
# [here](https://juliaai.github.io/DataScienceTutorials.jl/end-to-end/horse/#dealing_with_missing_values)).
# We already downloaded this data set immediately above.q

horse = DataFrames.DataFrame(csv_file); # convert to data frame
first(horse, 4)

#-

# From [the UCI
# docs](http://archive.ics.uci.edu/ml/datasets/Horse+Colic) we can
# surmise how each variable ought to be interpreted (a step in our
# work-flow that cannot reliably be left to the computer):

# variable                    | scientific type (interpretation)
# ----------------------------|-----------------------------------
# `:surgery`                  | Multiclass
# `:age`                      | Multiclass
# `:rectal_temperature`       | Continuous
# `:pulse`                    | Continuous
# `:respiratory_rate`         | Continuous
# `:temperature_extremities`  | OrderedFactor
# `:mucous_membranes`         | Multiclass
# `:capillary_refill_time`    | Multiclass
# `:pain`                     | OrderedFactor
# `:peristalsis`              | OrderedFactor
# `:abdominal_distension`     | OrderedFactor
# `:packed_cell_volume`       | Continuous
# `:total_protein`            | Continuous
# `:outcome`                  | Multiclass
# `:surgical_lesion`          | OrderedFactor
# `:cp_data`                  | Multiclass

# Let's see how MLJ will actually interpret the data, as it is
# currently encoded:

schema(horse)

# As a first correction step, we can get MLJ to "guess" the
# appropriate fix, using the `autotype` method:

autotype(horse)

#-

# Okay, this is not perfect, but a step in the right direction, which
# we implement like this:

coerce!(horse, autotype(horse));
schema(horse)

# All remaining `Count` data should be `Continuous`:

coerce!(horse, Count => Continuous);
schema(horse)

# We'll correct the remaining truant entries manually:

coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
schema(horse)

# ### Resources for Part 1
#
# - From the MLJ manual:
#    - [A preview of data type specification in
#   MLJ](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#A-preview-of-data-type-specification-in-MLJ-1)
#    - [Data containers and scientific types](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Data-containers-and-scientific-types-1)
#    - [Working with Categorical Data](https://alan-turing-institute.github.io/MLJ.jl/dev/working_with_categorical_data/)
# - [Summary](https://juliaai.github.io/ScientificTypes.jl/dev/#Summary-of-the-default-convention) of the MLJ convention for representing scientific types
# - [ScientificTypes.jl](https://juliaai.github.io/ScientificTypes.jl/dev/)
# - From Data Science Tutorials:
#     - [Data interpretation: Scientific Types](https://juliaai.github.io/DataScienceTutorials.jl/data/scitype/)
#     - [Horse colic data](https://juliaai.github.io/DataScienceTutorials.jl/end-to-end/horse/)
# - [UCI Horse Colic Data Set](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)

# ### Exercises for Part 1

# #### Exercise 1

# Try to guess how each code snippet below will evaluate:

scitype(42)

#-

questions = ["who", "why", "what", "when"]
scitype(questions)

#-

elscitype(questions)

#-

t = (3.141, 42, "how")
scitype(t)

#-

A = rand(2, 3)

# -

scitype(A)

#-

elscitype(A)

#-

using SparseArrays
Asparse = sparse(A)

#-

scitype(Asparse)

#-

C = coerce(A, Multiclass)

#-

scitype(C)

#-

elscitype(C)

#-

v = [1, 2, missing, 4]
scitype(v)

#-

elscitype(v)

#-

scitype(v[1:2])

# Can you guess at the general behavior of
# `scitype` with respect to tuples, abstract arrays and missing
# values? The answers are
# [here](https://github.com/juliaai/ScientificTypesBase.jl#2-the-scitype-and-scitype-methods)
# (ignore "Property 1").

# #### Exercise 2

# Coerce the following vector to make MLJ recognize it as a vector of
# ordered factors (with an appropriate ordering):

quality = ["good", "poor", "poor", "excellent", missing, "good", "excellent"]

#-

# #### Exercise 3 (fixing scitypes in a table)

# Fix the scitypes for the [House Prices in King
# County](https://mlr3gallery.mlr-org.com/posts/2020-01-30-house-prices-in-king-county/)
# dataset:

house_csv = urldownload("https://raw.githubusercontent.com/ablaom/"*
                        "MachineLearningInJulia2020/for-MLJ-version-0.16/"*
                        "data/house.csv");
house = DataFrames.DataFrame(house_csv)
first(house, 4)

# (Two features in the original data set have been deemed uninformative
# and dropped, namely `:id` and `:date`. The original feature
# `:yr_renovated` has been replaced by the `Bool` feature `is_renovated`.)

#


text:
# Set up Problem
model = Dynamics.car
costfun = Dynamics.car_costfun
xf = [0,1,0]
N = 51
n,m = model.n, model.m
bnd = BoundConstraint(n,m, x_min=[-0.5, -0.01, -Inf], x_max=[0.5, 1.01, Inf], u_min=[0.1,-2], u_max=1.5)
bnd1 = BoundConstraint(n,m, u_min=bnd.u_min)
bnd_x = BoundConstraint(n,m, x_min=[-0.5, -0.01, -Inf], x_max=[0.5, 1.01, Inf])
goal = goal_constraint(xf)
obs = (([0.2, 0.6], 0.25),
       ([-0.5, 0.5], 0.4))
obs1 = TO.planar_obstacle_constraint(n,m, obs[1]..., :obstacle1)
obs2 = TO.planar_obstacle_constraint(n,m, obs[2]..., :obstacle2)
con = Constraints(N)
con[1] += bnd1
for k = 2:N-1
    con[k] += bnd  + obs1 + obs2
end
con[N] += goal
prob = Problem(rk4(model), Objective(costfun, N), constraints=con, tf=3)

# Solve with ALTRO
initial_controls!(prob, ones(m,N-1))
ilqr = iLQRSolverOptions()
al = AugmentedLagrangianSolverOptions(opts_uncon=ilqr)
al.constraint_tolerance = 1e-2
al.constraint_tolerance_intermediate = 1e-1
al.kickout_max_penalty = false
al.verbose = true
opts = ALTROSolverOptions{Float64}(opts_al=al,
    projected_newton=false,
    projected_newton_tolerance=1e-2)
solver0 = solve!(prob, opts)
@test 1e-7 < max_violation(prob) < 1e-4
max_violation(prob)
solver0.stats[:time]

opts.projected_newton = true
opts.projected_newton_tolerance = 1e-2
opts.opts_pn.feasibility_tolerance = 1e-10
opts.opts_pn.active_set_tolerance = 1e-4
solver = SequentialNewtonSolver(prob, opts.opts_pn)
begin
    copyto!(solver.V.X, prob.X)
    copyto!(solver.V.U, prob.U)
    TO.newton_step!(prob, solver)
end
solver0 = ProjectedNewtonSolver(prob, opts.opts_pn)
TO.newton_step!(prob, solver)
TO.newton_step!(prob, solver0)

TO.update!(res, solver0)
Y,y = TO.active_constraints(res, solver0)
Hinv = inv(Diagonal(solver0.H))
S = Y*Hinv*Y'
cholesky(Symmetric(S))

res,solver = solve(prob, opts)
@test max_violation(res) < 1e-10z
solver.stats[:time]
@test solver.stats[:time] < solver0.stats[:time]

opts_pn = ProjectedNewtonSolverOptions{Float64}()
solver = SequentialNewtonSolver(res, opts_pn)
TO.newton_step!(res, solver)

opts = ALTROSolverOptions{Float64}(projected_newton=false, projected_newton_tolerance=1e-1, opts_al=opts_al)
res, = solve(prob, opts)
max_violation(res)

solver = SequentialNewtonSolver(res, opts.opts_pn)
solve(res, solver)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
include("guesses.jl")
include("answers.jl")

# list of possible answers
ANS_LIST = make_answers() # place to hold list of words so we only load it once

# list of allowable guesses
# GUESS_LIST = Set(ANS_LIST) ∪ Set(make_guesses()) |> collect |> sort

GUESS_LIST = vcat(ANS_LIST, make_guesses())

# reverse lookup dictionaries. Given a word, find where in the appropriate list
# to find that word.

ANS_DICT = Dict{String,Int}()
for j = 1:length(ANS_LIST)
    w = ANS_LIST[j]
    ANS_DICT[w] = j
end

NA = length(ANS_DICT)

@info "Answer dictionary has $NA entries"

GUESS_DICT = Dict{String,Int}()
for j = 1:length(GUESS_LIST)
    w = GUESS_LIST[j]
    GUESS_DICT[w] = j
end

NG = length(GUESS_DICT)
@info "Guess dictionary has $NG entries"

include("pre_compute_tools.jl")

@info "Precomputing scoring functions"

SCORE_MATRIX = all_pairs_compute()


text:
"""
This module uses classical laminate theory to calculate stresses and strains in
a thin plate laminate.  Analysis is based on chapter 3 of "Design and Analysis
of Composite Structures" by Kassapoglou.
"""
module Composites

import LinearAlgebra

# Some modules structs
include("structs.jl")
# ABD matrix calculations
include("abd.jl")
# Stress and strain calculations
include("stressstrain.jl")
# Material Failure Theories
include("failuretheories.jl")
# Buckling related functions
include("buckling.jl")

include("deprecated.jl")

end #module


text: Dict{String,Any}("equivalent_two" => Dict{String,Any}("value" => "The quick brown fox jumps over the lazy dog.","type" => "string"),"multiline_empty_four" => Dict{String,Any}("value" => "","type" => "string"),"multiline_empty_one" => Dict{String,Any}("value" => "","type" => "string"),"equivalent_three" => Dict{String,Any}("value" => "The quick brown fox jumps over the lazy dog.","type" => "string"),"equivalent_one" => Dict{String,Any}("value" => "The quick brown fox jumps over the lazy dog.","type" => "string"),"multiline_empty_two" => Dict{String,Any}("value" => "","type" => "string"),"multiline_empty_three" => Dict{String,Any}("value" => "","type" => "string"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
#=
day1:
- Julia version: 1.5.2
- Author: Paul.Mealor
- Date: 2020-12-01
=#

import AoC

function countIncreases(depths)
    increases = 0
    prev=typemax(eltype(depths))
    for depth in depths
        if depth > prev
            increases += 1
        end
        prev = depth
    end
    return increases
end
function countIncreases2(depths)
    return countIncreases((sum(depths[i:i+2]) for i in 1:length(depths)-2))
end

show(AoC.exampleInts(1,1) |> x -> @time countIncreases(x))
show(AoC.ints(1) |> x -> @time countIncreases(x))
show(AoC.exampleInts(1,1) |> x -> @time countIncreases2(x))
show(AoC.ints(1) |> x -> @time countIncreases2(x))


text:
module ForecastAccuracy

using Statistics

export me, mde, mae, gmae, mdae, mse, mdse, rmse, mpe, mdpe, mape, mdape, smape,
    smdape, maape, mase, mdase, rmsse, pl

err(real::Vector{T}, forecast::Vector{T}) where T = real .- forecast
abs_err(real::Vector{T}, forecast::Vector{T}) where T = abs.(err(real, forecast))
sq_err(real::Vector{T}, forecast::Vector{T}) where T = err(real, forecast) .^ 2
percent_err(real::Vector{T}, forecast::Vector{T}) where T = err(real, forecast) ./ real
abs_percent_err(real::Vector{T}, forecast::Vector{T}) where T = abs.(percent_err(real, forecast))

function scaled_err(real::Vector{T}, forecast::Vector{T}, steps_ahead::Int) where T
    denom = zero(T)
    n = length(real)
    for i in 1 + steps_ahead : n
        @inbounds denom += abs(real[i] - real[i - 1])
    end
    denom = denom/(n - steps_ahead)
    return err(real, forecast)/denom
end

function scaled_sq_err(real::Vector{T}, forecast::Vector{T}, steps_ahead::Int) where T
    denom = zero(T)
    n = length(real)
    for i in 1 + steps_ahead : n
        @inbounds denom += (real[i] - real[i - 1])^2
    end
    denom = denom/(n - steps_ahead)
    return mean(sq_err(real, forecast))/denom
end

function geomean(v::Vector{T}) where T
    s = zero(T)
    n = length(v)
    for i in 1 : n
        @inbounds s += log(v[i])
    end
    return exp(s / n)
end

"""
    me(real::Vector{T}, forecast::Vector{T}) where T

Mean Error.
"""
function me(real::Vector{T}, forecast::Vector{T}) where T
    return mean(err(real, forecast))
end

"""
    mde(real::Vector{T}, forecast::Vector{T}) where T

Median Error.
"""
function mde(real::Vector{T}, forecast::Vector{T}) where T
    return median(err(real, forecast))
end

"""
    mae(real::Vector{T}, forecast::Vector{T}) where T

Mean Absolute Error.
"""
function mae(real::Vector{T}, forecast::Vector{T}) where T
    return mean(abs_err(real, forecast))
end

"""
    gmae(real::Vector{T}, forecast::Vector{T}) where T

Geometric Mean Absolute Error.
"""
function gmae(real::Vector{T}, forecast::Vector{T}) where T
    return geomean(abs_err(real, forecast))
end

"""
    mdae(real::Vector{T}, forecast::Vector{T}) where T

Median Absolute Error.
"""
function mdae(real::Vector{T}, forecast::Vector{T}) where T
    return median(abs_err(real, forecast))
end

"""
    mse(real::Vector{T}, forecast::Vector{T}) where T

Mean Square Error.
"""
function mse(real::Vector{T}, forecast::Vector{T}) where T
    return mean(sq_err(real, forecast))
end

"""
    mdse(real::Vector{T}, forecast::Vector{T}) where T

Median Square Error.
"""
function mdse(real::Vector{T}, forecast::Vector{T}) where T
    return median(sq_err(real, forecast))
end

"""
    rmse(real::Vector{T}, forecast::Vector{T}) where T

Root Mean Squared Error.
"""
function rmse(real::Vector{T}, forecast::Vector{T}) where T
    return mean(sqrt.(sq_err(real, forecast)))
end

"""
    mpe(real::Vector{T}, forecast::Vector{T}) where T

Mean Percentual Error.
"""
function mpe(real::Vector{T}, forecast::Vector{T}) where T
    return mean(percent_err(real, forecast))
end

"""
    mdpe(real::Vector{T}, forecast::Vector{T}) where T

Median Percentual Error.
"""
function mdpe(real::Vector{T}, forecast::Vector{T}) where T
    return median(percent_err(real, forecast))
end

"""
    mape(real::Vector{T}, forecast::Vector{T}) where T

Mean Absolute Percentual Error.
"""
function mape(real::Vector{T}, forecast::Vector{T}) where T
    return mean(abs_percent_err(real, forecast))
end

"""
    mdape(real::Vector{T}, forecast::Vector{T}) where T

Median Absolute Percentual Error.
"""
function mdape(real::Vector{T}, forecast::Vector{T}) where T
    return median(abs_percent_err(real, forecast))
end

"""
    smape(real::Vector{T}, forecast::Vector{T}) where T

Symmetric Mean Absolute Percentual Error.
"""
function smape(real::Vector{T}, forecast::Vector{T}) where T
    return mean((2 .* abs_err(real, forecast) ./ (real .+ forecast)))
end

"""
    smdape(real::Vector{T}, forecast::Vector{T}) where T

Symmetric Median Absolute Percentual Error.
"""
function smdape(real::Vector{T}, forecast::Vector{T}) where T
    return median((2 .* abs_err(real, forecast) ./ (real .+ forecast)))
end

"""
    maape(real::Vector{T}, forecast::Vector{T}) where T

Mean Arctangent Absolute Percentage Error. As developed in https://www.sciencedirect.com/science/article/pii/S0169207016000121
"""
function maape(real::Vector{T}, forecast::Vector{T}) where T
    return mean(atan.(abs_percent_err(real, forecast)))
end

"""
    mase(real::Vector{T}, forecast::Vector{T}; steps_ahead::Int = 1) where T

Mean Absolute Scaled Error. As developed in https://robjhyndman.com/papers/mase.pdf
"""
function mase(real::Vector{T}, forecast::Vector{T}; steps_ahead::Int = 1) where T
    return mean(abs.(scaled_err(real, forecast, steps_ahead)))
end

"""
    mdase(real::Vector{T}, forecast::Vector{T}; steps_ahead::Int = 1) where T

Median Absolute Scaled Error. As developed in https://robjhyndman.com/papers/mase.pdf
"""
function mdase(real::Vector{T}, forecast::Vector{T}; steps_ahead::Int = 1) where T
    return median(abs.(scaled_err(real, forecast, steps_ahead)))
end

"""
    rmsse(real::Vector{T}, forecast::Vector{T}; steps_ahead::Int = 1) where T

Root Mean Squared Scaled Error.
"""
function rmsse(real::Vector{T}, forecast::Vector{T}; steps_ahead::Int = 1) where T
    return sqrt(scaled_sq_err(real, forecast, steps_ahead))
end

"""
    pl(real::Vector{T}, quantile_forecast::Vector{T}, tau::T) where T

Pinball loss function defined in terms of the quantile `tau`
"""
function pl(real::Vector{T}, quantile_forecast::Vector{T}; tau::T = 0.5) where T
    @assert ((tau <= 1) && (tau >= 0)) ""
    pinball_loss = zero(T)
    for i in 1:length(real)
        if real[i] >= quantile_forecast[i]
            pinball_loss += tau * (real[i] - quantile_forecast[i])
        else
            pinball_loss += (1 - tau) * (quantile_forecast[i] - real[i])
        end
    end
    return pinball_loss
end

const METRICS = [
    me,
    mde,
    mae,
    gmae,
    mdae,
    mse,
    mdse,
    rmse,
    mpe,
    mdpe,
    mape,
    mdape,
    smape,
    smdape,
    maape,
    mase,
    mdase,
    rmsse,
    pl
]

end # module


text:
using Distributions

lRange = 1; uRange = 9;         unidd    = DiscreteUniform(lRange,uRange)
p = 0.8;                         bernd    = Bernoulli(p)
n = 20; p = 0.8;                 bind     = Binomial(n,p)
ps = [0.2,0.5,0.2,0.1];          catd     = Categorical(ps)
n=20; ps = [0.2,0.5,0.3];        multind  = Multinomial(n, ps)
p = 0.2;                         geomd    = Geometric(p)
sAv = 10; fAv = 20; nTrials = 3; hypgeomd = Hypergeometric(sAv, fAv, nTrials)
# hypergeometrical multivariate form not available
rate = 10;                       poissd   = Poisson(rate)
nSucc =10; p = 0.2;              pascd    = NegativeBinomial(nSucc,p) # Return the number of failures before n successes instead of total trials to n successes

lRange = 1; uRange = 9;          unicd    = Uniform(lRange,uRange)
rate = 2.5;                      expd     = Exponential(rate)
μ = 2; σsq = 5 ;                 normd    = Normal(μ,sqrt(σsq))  # parametrised by sd!
n = 5; rate = 2.5;               erld     = Erlang(n,rate)
μ = 2.5; σ = 2.5;                cuachyd  = Cauchy(μ, σ)
df = 3;                          chisqd   = Chisq(df)
df = 4.5;                        td       = TDist(df)
df1 = 6 ; df2 = 7;               fd       = FDist(df1, df2)
shapeα = 2.5; shapeβ=4.5;        betad    = Beta(shapeα,shapeβ)
shapeα = 2.5; rateβ = 4.5;       gammad   = Gamma(shapeα,1/rateβ)

uvds = [unidd, bernd, bind, catd, geomd, hypgeomd, poissd, pascd,
        unicd, expd, normd, erld, cuachyd, chisqd,  td, fd, betad, gammad]
mvds = [multind]

x = 3
for d in uvds
    println("--------")
    println(d)
    println("Mean: $(mean(d))")
    println("Variance: $(var(d))")
    d != chisqd && println("Median: $(median(d))")
    println("Sample: $(rand(d))")
    println("PDF: $(pdf(d,x))")
    println("CDF: $(cdf(d,x))")
    println("20% quantile: $(quantile(d,0.2))")
end

xs = [3,12,5]
for d in mvds
    println("--------")
    println(d)
    println("Mean: $(mean(d))")
    println("Variance: $(var(d))")
    #println("Median: $(median(d))")
    println("Sample: $(rand(d))")
    println("PDF: $(pdf(d,xs))")
    #println("CDF: $(cdf(d,xs))")
    #println("20% quantile: $(quantile(d,0.2))")
end


text:
amean(A) = sum(A)/length(A)

gmean(A) = prod(A)^(1/length(A))

hmean(A) = length(A)/sum(1./A)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
dao_palette = [
    colorant"#d77255",
    colorant"#009afa",
    colorant"#707070",
    colorant"#21ab74",
    colorant"#ba3030",
    colorant"#9467bd"
]

_themes[:dao] = PlotTheme(
    background = :white,
    framestyle = :box,
    grid=true,
    gridalpha=0.4,
    linewidth=1.4,
    markerstrokewidth=0,
    fontfamily = "Computer Modern",
    colorgradient = :magma,
    guidefontsize=12,
    titlefontsize=12,
    tickfontsize=8,
    palette = dao_palette,
    minorgrid=true,
    minorticks = 5,
    gridlinewidth = 0.7,
    minorgridalpha=0.06,
    legend=:outertopright,
    # marker=:auto
)


text:
"""
    Diagnostics

Accumulate mean fields and covariance statistics on the computational grid.

"""
module Diagnostics

export DiagnosticsGroup,
    setup_atmos_default_diagnostics,
    setup_atmos_core_diagnostics,
    setup_dump_state_and_aux_diagnostics,
    VecGrad,
    compute_vec_grad,
    SphericalCoord,
    compute_vec_grad_spherical,
    Vorticity,
    compute_vorticity

using Dates
using FileIO
using JLD2
using MPI
using OrderedCollections
using Printf
using StaticArrays

using ClimateMachine
using ..DGmethods
using ..DGmethods:
    number_state_conservative,
    vars_state_conservative,
    number_state_auxiliary,
    vars_state_auxiliary,
    vars_state_gradient_flux,
    number_state_gradient_flux
using ..Mesh.Interpolation
using ..MPIStateArrays
using ..VariableTemplates
using ..Writers

Base.@kwdef mutable struct Diagnostic_Settings
    mpicomm::MPI.Comm = MPI.COMM_WORLD
    dg::Union{Nothing, DGModel} = nothing
    Q::Union{Nothing, MPIStateArray} = nothing
    starttime::Union{Nothing, String} = nothing
    output_dir::Union{Nothing, String} = nothing
end
const Settings = Diagnostic_Settings()

"""
    init(mpicomm, dg, Q, starttime, output_dir)

Initialize the diagnostics collection module -- save the parameters into `Settings`.
"""
function init(
    mpicomm::MPI.Comm,
    dg::DGModel,
    Q::MPIStateArray,
    starttime::String,
    output_dir::String,
)
    Settings.mpicomm = mpicomm
    Settings.dg = dg
    Settings.Q = Q
    Settings.starttime = starttime
    Settings.output_dir = output_dir
end

"""
    DiagnosticsGroup

Holds a set of diagnostics that share a collection interval, a filename prefix,
and an interpolation.

TODO: to be completed; will be restructured.
"""
mutable struct DiagnosticsGroup
    name::String
    init::Function
    fini::Function
    collect::Function
    interval::String
    out_prefix::String
    writer::AbstractWriter
    interpol::Union{Nothing, InterpolationTopology}
    project::Bool
    num::Int

    DiagnosticsGroup(
        name,
        init,
        fini,
        collect,
        interval,
        out_prefix,
        writer,
        interpol,
        project,
    ) = new(
        name,
        init,
        fini,
        collect,
        interval,
        out_prefix,
        writer,
        interpol,
        project,
        0,
    )
end

function (dgngrp::DiagnosticsGroup)(currtime; init = false, fini = false)
    if init
        dgngrp.init(dgngrp, currtime)
    else
        dgngrp.collect(dgngrp, currtime)
    end
    if fini
        dgngrp.fini(dgngrp, currtime)
    end
    dgngrp.num += 1
end

"""
    setup_atmos_default_diagnostics(
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
        project = true,
    )

Create and return a `DiagnosticsGroup` containing the "AtmosDefault"
diagnostics. All the diagnostics in the group will run at the specified
`interval`, be interpolated to the specified boundaries and resolution,
and will be written to files prefixed by `out_prefix` using `writer`.
"""
function setup_atmos_default_diagnostics(
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
    project = true,
)
    return DiagnosticsGroup(
        "AtmosDefault",
        Diagnostics.atmos_default_init,
        Diagnostics.atmos_default_fini,
        Diagnostics.atmos_default_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        project,
    )
end

"""
    setup_atmos_core_diagnostics(
            interval::Int,
            out_prefix::String;
            writer::AbstractWriter,
            interpol = nothing,
            project  = true)

Create and return a `DiagnosticsGroup` containing the "AtmosCore"
diagnostics. All the diagnostics in the group will run at the
specified `interval`, be interpolated to the specified boundaries
and resolution, and will be written to files prefixed by `out_prefix`
using `writer`.
"""
function setup_atmos_core_diagnostics(
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
    project = true,
)
    return DiagnosticsGroup(
        "AtmosCore",
        Diagnostics.atmos_core_init,
        Diagnostics.atmos_core_fini,
        Diagnostics.atmos_core_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        project,
    )
end

"""
    setup_dump_state_and_aux_diagnostics(
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
        project = true,
    )

Create and return a `DiagnosticsGroup` containing a diagnostic that
simply dumps the state and aux variables at the specified `interval`
after being interpolated, into NetCDF files prefixed by `out_prefix`.
"""
function setup_dump_state_and_aux_diagnostics(
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
    project = true,
)
    return DiagnosticsGroup(
        "DumpStateAndAux",
        Diagnostics.dump_state_and_aux_init,
        Diagnostics.dump_state_and_aux_fini,
        Diagnostics.dump_state_and_aux_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        project,
    )
end

"""
    visitQ()

Helper macro to iterate over the DG grid. Generates the needed loops
and indices: `eh`, `ev`, `e`, `k,`, `j`, `i`, `ijk`.
"""
macro visitQ(nhorzelem, nvertelem, Nqk, Nq, expr)
    return esc(quote
        for eh in 1:nhorzelem
            for ev in 1:nvertelem
                e = ev + (eh - 1) * nvertelem
                for k in 1:Nqk
                    for j in 1:Nq
                        for i in 1:Nq
                            ijk = i + Nq * ((j - 1) + Nq * (k - 1))
                            $expr
                        end
                    end
                end
            end
        end
    end)
end

# Helpers to extract data from the various state arrays
function extract_state_conservative(dg, state_conservative, ijk, e)
    bl = dg.balance_law
    FT = eltype(state_conservative)
    num_state_conservative = number_state_conservative(bl, FT)
    local_state_conservative = MArray{Tuple{num_state_conservative}, FT}(undef)
    for s in 1:num_state_conservative
        local_state_conservative[s] = state_conservative[ijk, s, e]
    end
    return Vars{vars_state_conservative(bl, FT)}(local_state_conservative)
end
function extract_state_auxiliary(dg, state_auxiliary, ijk, e)
    bl = dg.balance_law
    FT = eltype(state_auxiliary)
    num_state_auxiliary = number_state_auxiliary(bl, FT)
    local_state_auxiliary = MArray{Tuple{num_state_auxiliary}, FT}(undef)
    for s in 1:num_state_auxiliary
        local_state_auxiliary[s] = state_auxiliary[ijk, s, e]
    end
    return Vars{vars_state_auxiliary(bl, FT)}(local_state_auxiliary)
end
function extract_state_gradient_flux(dg, state_gradient_flux, ijk, e)
    bl = dg.balance_law
    FT = eltype(state_gradient_flux)
    num_state_gradient_flux = number_state_gradient_flux(bl, FT)
    local_state_gradient_flux =
        MArray{Tuple{num_state_gradient_flux}, FT}(undef)
    for s in 1:num_state_gradient_flux
        local_state_gradient_flux[s] = state_gradient_flux[ijk, s, e]
    end
    return Vars{vars_state_gradient_flux(bl, FT)}(local_state_gradient_flux)
end

include("atmos_common.jl")
include("thermo.jl")
include("atmos_default.jl")
include("atmos_core.jl")
include("dump_state_and_aux.jl")
include("diagnostic_fields.jl")

end # module Diagnostics


text:
using BINOD
using Test

@testset "BINOD.jl" begin
    @test sprint(print, "BINOD") == "BINOD"
end


text:
using Gadfly

Gadfly.push_theme(Theme(major_label_font="Helvetica", major_label_font_size=16pt,
                        minor_label_font="Helvetica", minor_label_font_size=14pt,
                        key_title_font="Helvetica", key_title_font_size=16pt,
                        key_label_font="Helvetica", key_label_font_size=14pt,
                        line_width=0.8mm, point_size=0.5mm,
                        lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2),
                        default_color=colorant"#000000"))

colors = [colorant"#e41a1c", colorant"#377eb8", colorant"#4daf4a",
          colorant"#984ea3", colorant"#ff7f00", colorant"#ffff33"]

function plot_spikes(n::Network, filename::String)
    spikes = n.history["spikes"]
    s = map(x->ind2sub(spikes, x), find(spikes))
    ntype = ["exc" for i in 1:length(n.exc)]
    ntype[n.inputs] = "input"
    ntype[n.outputs] = "output"
    ntype[n.inh] = "inh"
    y = map(x->x[1], s)
    color = map(x->ntype[x], y)
    plt = plot(x=map(x->x[2], s), y=y, color=color, Geom.point,
              Guide.xlabel(nothing), Guide.ylabel(nothing),
              Scale.color_discrete_manual(colors...))

    draw(PDF(filename, 8inch, 6inch), plt)
end


text:
using IcanProj
using NetCDF
using DataFrames
using ProgressMeter
using CSV
using Statistics

function results_table(path)

    df_all = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

    for variable in ["gsurf", "hatmo", "latmo", "melt", "rnet", "rof", "snowdepth", "swe"], #, "tsoil", "tsurf"],
        spaceres in ["1km", "50km"],
        cfg in 1:32

        @show variable, spaceres, cfg

        file = joinpath(path, "results_$(cfg)", "$(variable)_$(spaceres).nc")

        idname = Symbol(ncgetatt(file, "id", "id"))

        dataname = Symbol("$(variable)_cfg$(cfg)_$(spaceres)_mean")

        ids = convert(Array{Int64}, ncread(file, "id"))

        data = ncread(file, variable)

        df_nc = DataFrame()

        df_nc[idname] = ids

        df_nc[dataname] = mean(data, 2)[:]

        df_all = join(df_all, df_nc, on = idname)

        ncclose()

    end

    return df_all

end

path = "/data04/jmg/fsm_simulations/netcdf/fsmres"

df_all = results_table(path)

file = joinpath(dirname(pathof(IcanProj)), "..", "data", "table_results.txt")

df_all |> CSV.write(file)


text:
export bsplsurf, bsplsurfu, dbsurf

"""
	bsplsurf(B[], ordx, ordy, npts, mpts, p1, p2)

Calculate a Cartesian product B-spline surface using open uniform knot vectors.

---

# Arguments
- `B::Array{Float64}`: array containing the control net vertices
- `ordx::Int64`: order in the x direction
- `ordy::Int64`: order in the y direction
- `npts::Int64`: the number of control net vertices in x direction
- `mpts::Int64`: the number of control net vertices in y direction
- `p1::Int64`: number of parametric lines in the x direction
- `p2::Int64`: number of parametric lines in the y direction

---

# Examples
```jldoctest
julia> b = [
        0.0 0.0 0.0 1.0 1.0 1.0 2.0 2.0 2.0;
        0.0 2.0 4.0 0.0 2.0 4.0 0.0 2.0 4.0;
        2.0 0.0 2.0 0.0 1.0 0.0 2.0 0.0 2.0
       ]
3×9 Array{Float64,2}:
 0.0  0.0  0.0  1.0  1.0  1.0  2.0  2.0  2.0
 0.0  2.0  4.0  0.0  2.0  4.0  0.0  2.0  4.0
 2.0  0.0  2.0  0.0  1.0  0.0  2.0  0.0  2.0

julia> bsplsurf(b, 2, 2, 3, 3, 5, 5)[1]
3×25 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.5  0.5   0.5  0.5   0.5  1.0  …  1.0  1.5  1.5   1.5  1.5   1.5  2.0  2.0  2.0  2.0  2.0
 0.0  1.0  2.0  3.0  4.0  0.0  1.0   2.0  3.0   4.0  0.0     4.0  0.0  1.0   2.0  3.0   4.0  0.0  1.0  2.0  3.0  4.0
 2.0  1.0  0.0  1.0  2.0  1.0  0.75  0.5  0.75  1.0  0.0     0.0  1.0  0.75  0.5  0.75  1.0  2.0  1.0  0.0  1.0  2.0

julia> bsplsurf(b, 2, 2, 3, 3, 5, 5)[2]
44-element Array{Array{Int64,1},1}:
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 5]
 [5, 6]
 [6, 7]
 [7, 8]
 [8, 9]
 [9, 10]
 [10, 11]
 ⋮
 [12, 17]
 [13, 18]
 [14, 19]
 [15, 20]
 [16, 21]
 [17, 22]
 [18, 23]
 [19, 24]
 [20, 25]

julia> bsplsurf(b, 2, 2, 3, 3, 5, 5)[3]
16-element Array{Array{Int64,1},1}:
 [1, 2, 6, 7]
 [2, 3, 7, 8]
 [3, 4, 8, 9]
 [4, 5, 9, 10]
 [6, 7, 11, 12]
 [7, 8, 12, 13]
 [8, 9, 13, 14]
 [9, 10, 14, 15]
 [11, 12, 16, 17]
 [12, 13, 17, 18]
 [13, 14, 18, 19]
 [14, 15, 19, 20]
 [16, 17, 21, 22]
 [17, 18, 22, 23]
 [18, 19, 23, 24]
 [19, 20, 24, 25]
```

---

_By Paolo Macciacchera, Elia Onofri_
"""

function bsplsurf(B::Array{Float64,2}, ordx::Int64, ordy::Int64, npts::Int64, mpts::Int64, p1::Int64, p2::Int64)::Tuple{Array{Float64, 2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}

    nplusc = npts + ordx
    mplusc = mpts + ordy
    x = zeros(nplusc)
    y = zeros(mplusc)
    nbasis = zeros(1, npts)
    mbasis = zeros(1, mpts)
    q = zeros(3, p1*p2)

    #generate the open uniform knot vectors
    x = knot(npts, ordx)
    y = knot(mpts, ordy)

    icount = 0
    #calculate the points on the B-spline surface
    stepu = x[nplusc] / (p1 - 1)
    stepw = y[mplusc] / (p2 - 1)
    for u = 0 : stepu : x[nplusc]
        nbasis = basis(npts, ordx, u, x)
        for w = 0 : stepw : y[mplusc]
            mbasis = basis(mpts, ordy, w, y)
            icount = icount + 1
            for i = 1 : npts
                for j = 1 : mpts
                    j1 = mpts * (i - 1) + j
                    for s = 1 : 3
                        q[s, icount] = q[s, icount] + B[s, j1] * nbasis[i] * mbasis[j]
                    end
                end
            end
        end
    end

    EV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for i = 1 : p1 * p2 - 1
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+1
        push!(EV,C)
    end

    for i = 1 : p2 * (p1 - 1)
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+p2
        push!(EV,C)
    end

    FV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for j = 0 : p1 - 2
        for i = 1 : p2 - 1
            C = Array{Int64}(4)
            C[1] = j * p2 + i
            C[2] = j * p2 + i + 1
            C[3] = (j + 1) * p2 + i
            C[4] = (j + 1) * p2 + i + 1
            push!(FV,C)
        end
    end
    return(q, EV, FV)
end

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

"""
	bsplsurfu(B[], ordx, ordy, npts, mpts, p1, p2)

Calculate a Cartesian product B-spline surface using periodic uniform knot vectors.

---

# Arguments
- `B::Array{Float64}`: array containing the control net vertices
- `ordx::Int64`: order in the x direction
- `ordy::Int64`: order in the y direction
- `npts::Int64`: the number of control net vertices in x direction
- `mpts::Int64`: the number of control net vertices in y direction
- `p1::Int64`: number of parametric lines in the x direction
- `p2::Int64`: number of parametric lines in the y direction

---

# Examples
```jldoctest
julia> b = [
        0.0 0.0 0.0 1.0 1.0 1.0 2.0 2.0 2.0;
        0.0 2.0 4.0 0.0 2.0 4.0 0.0 2.0 4.0;
        2.0 0.0 2.0 0.0 1.0 0.0 2.0 0.0 2.0
       ]
3×9 Array{Float64,2}:
 0.0  0.0  0.0  1.0  1.0  1.0  2.0  2.0  2.0
 0.0  2.0  4.0  0.0  2.0  4.0  0.0  2.0  4.0
 2.0  0.0  2.0  0.0  1.0  0.0  2.0  0.0  2.0

julia> bsplsurfu(b, 2, 2, 3, 3, 5, 5)[1]
3×25 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.5  0.5   …  1.5   1.5  2.0  2.0  2.0  2.0  2.0
 0.0  1.0  2.0  3.0  4.0  0.0  1.0      3.0   4.0  0.0  1.0  2.0  3.0  4.0
 2.0  1.0  0.0  1.0  2.0  1.0  0.75     0.75  1.0  2.0  1.0  0.0  1.0  2.0

julia> bsplsurfu(b, 2, 2, 3, 3, 5, 5)[2]
44-element Array{Array{Int64,1},1}:
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 5]
 [5, 6]
 [6, 7]
 [7, 8]
 [8, 9]
 [9, 10]
 [10, 11]
 ⋮
 [12, 17]
 [13, 18]
 [14, 19]
 [15, 20]
 [16, 21]
 [17, 22]
 [18, 23]
 [19, 24]
 [20, 25]

julia> bsplsurfu(b, 2, 2, 3, 3, 5, 5)[3]
16-element Array{Array{Int64,1},1}:
 [1, 2, 6, 7]
 [2, 3, 7, 8]
 [3, 4, 8, 9]
 [4, 5, 9, 10]
 [6, 7, 11, 12]
 [7, 8, 12, 13]
 [8, 9, 13, 14]
 [9, 10, 14, 15]
 [11, 12, 16, 17]
 [12, 13, 17, 18]
 [13, 14, 18, 19]
 [14, 15, 19, 20]
 [16, 17, 21, 22]
 [17, 18, 22, 23]
 [18, 19, 23, 24]
 [19, 20, 24, 25]
```

---

_By Paolo Macciacchera, Elia Onofri_
"""

function bsplsurfu(B::Array{Float64,2}, ordx::Int64, ordy::Int64, npts::Int64, mpts::Int64, p1::Int64, p2::Int64)::Tuple{Array{Float64, 2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}

    nplusc = npts + ordx
    mplusc = mpts + ordy
    x = zeros(nplusc)
    y = zeros(mplusc)
    nbasis = zeros(1, npts)
    mbasis = zeros(1, mpts)
    q = zeros(3, p1 * p2)

    #generate the open uniform knot vectors
    x = knotu(npts, ordx)
    y = knotu(mpts, ordy)

    icount = 0
    #calculate the points on the B-spline surface
    stepu = (npts - ordx + 1) / (p1 - 1)
    stepw = (mpts - ordy + 1) / (p2 - 1)
    for u = ordx-1 : stepu : npts
        nbasis = basis(npts, ordx, u, x)
        for w = ordy-1 : stepw : mpts
            mbasis = basis(mpts, ordy, w, y)
            icount = icount + 1
            for i = 1 : npts
                for j = 1 : mpts
                    for s = 1 : 3
                        j1 = mpts * (i - 1) + j
                        q[s, icount] = q[s, icount] + B[s, j1] * nbasis[i] * mbasis[j]
                    end
                end
            end
        end
    end

    EV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for i = 1 : p1 * p2 - 1
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+1
        push!(EV,C)
    end

    for i = 1 : p2 * (p1 - 1)
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+p2
        push!(EV,C)
    end

    FV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for j = 0 : p1 - 2
        for i = 1 : p2 - 1
            C = Array{Int64}(4)
            C[1] = j * p2 + i
            C[2] = j * p2 + i + 1
            C[3] = (j + 1) * p2 + i
            C[4] = (j + 1) * p2 + i + 1
            push!(FV,C)
        end
    end
    return(q, EV, FV)
end

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

"""
	dbsurf(B[], ordx, ordy, npts, mpts, p1, p2)

Calculate a Cartesian product B-spline surface and its derivatives using open uniform knot vectors.

---

# Arguments
- `B::Array{Float64}`: array containing the control net vertices
- `ordx::Int64`: order in the x direction
- `ordy::Int64`: order in the y direction
- `npts::Int64`: the number of control net vertices in x direction
- `mpts::Int64`: the number of control net vertices in y direction
- `p1::Int64`: number of parametric lines in the x direction
- `p2::Int64`: number of parametric lines in the y direction

---

# Examples
```jldoctest
julia> b = [
        0.0 0.0 0.0 1.0 1.0 1.0 2.0 2.0 2.0;
        0.0 2.0 4.0 0.0 2.0 4.0 0.0 2.0 4.0;
        2.0 0.0 2.0 0.0 1.0 0.0 2.0 0.0 2.0
       ]
3×9 Array{Float64,2}:
 0.0  0.0  0.0  1.0  1.0  1.0  2.0  2.0  2.0
 0.0  2.0  4.0  0.0  2.0  4.0  0.0  2.0  4.0
 2.0  0.0  2.0  0.0  1.0  0.0  2.0  0.0  2.0

julia> bsplsurf(arr, 4, 4, 5, 5, 100, 100)[1] == dbsurf(arr, 4, 4, 5, 5, 100, 100)[1]
true
julia> bsplsurf(arr, 4, 4, 5, 5, 100, 100)[2] == dbsurf(arr, 4, 4, 5, 5, 100, 100)[7]
true
julia> bsplsurf(arr, 4, 4, 5, 5, 100, 100)[3] == dbsurf(arr, 4, 4, 5, 5, 100, 100)[8]
true
```

---

_By Paolo Macciacchera, Elia Onofri_
"""

function dbsurf(B::Array{Float64,2}, ordx::Int64, ordy::Int64, npts::Int64, mpts::Int64, p1::Int64, p2::Int64)::Tuple{Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}
    nplusc = npts + ordx
    mplusc = mpts + ordy
    x = zeros(nplusc)
    y = zeros(mplusc)
    nbasis = zeros(1, npts)
    mbasis = zeros(1, mpts)
    d1nbasis = zeros(1, npts)
    d1mbasis = zeros(1, mpts)
    d2nbasis = zeros(1, npts)
    d2mbasis = zeros(1, mpts)
    q = zeros(3, p1 * p2)
    qu = zeros(3, p1 * p2)
    qw = zeros(3, p1 * p2)
    quu = zeros(3, p1 * p2)
    quw = zeros(3, p1 * p2)
    qww = zeros(3, p1 * p2)

    x = knot(npts, ordx)
    y = knot(mpts, ordy)

    icount = 0
    #calculate the points on the B-spline surface
    stepu = x[nplusc] / (p1 - 1)
    stepw = y[mplusc] / (p2 - 1)
    for u = 0 : stepu : x[nplusc]
        nbasis, d1nbasis, d2nbasis = dbasis(npts, ordx, u, x)
        for w = 0 : stepw : y[mplusc]
            mbasis, d1mbasis, d2mbasis = dbasis(mpts, ordy, w, y)
            icount = icount + 1
            for i = 1 : npts
                for j = 1 : mpts
                    for s = 1 : 3
                        j1 = mpts * (i - 1) + j
                        q[s, icount] = q[s, icount] + B[s, j1] * nbasis[i] * mbasis[j]
                        qu[s, icount] = qu[s, icount] + B[s, j1] * d1nbasis[i] * mbasis[j]
                        qw[s, icount] = qw[s, icount] + B[s, j1] * nbasis[i] * d1mbasis[j]
                        quu[s, icount] = quu[s, icount] + B[s, j1] * d1nbasis[i] * d1mbasis[j]
                        quw[s, icount] = quw[s, icount] + B[s, j1] * d2nbasis[i] * mbasis[j]
                        qww[s, icount] = qww[s, icount] + B[s, j1] * nbasis[i] * d2mbasis[j]
                    end
                end
            end
        end
    end

    EV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for i = 1 : p1 * p2 - 1
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+1
        push!(EV,C)
    end

    for i = 1 : p2 * (p1 - 1)
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+p2
        push!(EV,C)
    end

    FV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for j = 0 : p1 - 2
        for i = 1 : p2 - 1
            C = Array{Int64}(4)
            C[1] = j * p2 + i
            C[2] = j * p2 + i + 1
            C[3] = (j + 1) * p2 + i
            C[4] = (j + 1) * p2 + i + 1
            push!(FV,C)
        end
    end

    return(q, qu, qw, quu, quw, qww, EV, FV)

end


text:
"""
    waveguide_domains = planar_waveguides(;width, index, direction, x0, y0, waveguide_number)

returns vector of Domains
"""
function planar_waveguide_domains(width::Number, index::Number, direction::Symbol, x0::Number, y0::Number, waveguide_number::Int, background_index::Number)
    if direction ∈ [:ew, :lr, :h, :horizontal, :x]
        bulk_domain_type = :bulk_planar_waveguide_x
        which_asymptote1 = :left
        which_asymptote2 = :right
        a = 1e5
        domain = Rectangle(a=a, b=width, x0=x0-a/2, y0=y0-width/2, Dict(:n₁=>index))
    elseif direction ∈ [:ns, :ud, :tb, :v, :vertical, :y]
        bulk_domain_type = :bulk_planar_waveguide_y
        which_asymptote1 = :top
        which_asymptote2 = :bottom
        b=1e5
        domain = Rectangle(a=width, b=b, x0=x0-width/2, y0=y0-b/2, Dict(:n₁=>index))
    else
        throw(ArgumentError("invalid waveguide direction $(direction)"))
    end

    bgd1 = Domain(; :domain_type => :planar_waveguide_background, :n₁ => background_index, :which_waveguide => waveguide_number, :which_asymptote => which_asymptote1)
    bgd2 = Domain(; :domain_type => :planar_waveguide_background, :n₁ => background_index, :which_waveguide => waveguide_number+1, :which_asymptote => which_asymptote2)
    wvg1 = Domain(domain; :domain_type => :planar_waveguide, :which_waveguide => waveguide_number, :which_asymptote => which_asymptote1)
    wvg2 = Domain(domain; :domain_type => :planar_waveguide, :which_waveguide => waveguide_number+1, :which_asymptote => which_asymptote2)
    wvg_bulk = Domain(domain; :domain_type => bulk_domain_type)

    return vcat(bgd1, bgd2, wvg1, wvg2, wvg_bulk)
end

"""
    new_sys = add_planar_waveguide(sys; width, index, direction, position)

    new_sim = add_planar_waveguide(sim; width, index, direction, position)
"""
function add_planar_waveguide(sys::System; width::Number, index::Number, direction::Symbol, x0::Number, y0::Number, background_index::Number=1)
    waveguide_number = length(sys.waveguides)+1
    waveguide_domains = planar_waveguide_domains(width, index, direction, x0, y0, waveguide_number, background_index)
    return System(vcat(waveguide_domains,sys.domains))
end
function add_planar_waveguide(sim::Simulation; width::Number, index::Number, direction::Symbol, x0::Number, y0::Number, background_index::Number=sim.sys.domains[end].domain_params[:n₁])
    sys = add_planar_waveguide(sim.sys; width=width, index=index, direction=direction, x0=x0, y0=y0, background_index=background_index)
    return Simulation(sim; :sys => sys)
end

"""
    new_sim = add_planar_waveguides(sim, dim; width, index, position)
"""
function add_planar_waveguides(sim::Simulation, dim::Int; width::Number, index::Number, direction::Symbol, x0::Number, y0::Number, background_index::Number=sim.sys.domains[end].params[:n₁])
    if dim == 1
        return add_horizontal_planar_waveguides(sim; width=width, index=index, y0=position, background_index=background_index)
    elseif dim == 2
        return add_vertical_planar_waveguides(sim; width=width, index=index, x0=position, background_index=background_index)
    else
        throw(ArgumentError("invalid dimension $dim"))
    end
end
function add_planar_waveguides(sim::Simulation, dir::Symbol; width::Number, index::Number, direction::Symbol, x0::Number, y0::Number, background_index::Number=sim.sys.domains[end].params[:n₁])
    @assert dir ∈ [:v, :V, :vertical, :Vertical, :vert, :Vert, :h, :H, :horizontal, :Horizontal, :horiz, :Horiz] "unrecognized direction $direction"
    if dir ∈ [:v, :V, :vertical, :Vertical, :vert, :Vert]
        return add_planar_waveguides(sim, 2; width=width, index=index, direction=direction, x0=x0, y0=y0, background_index=background_index)
    else
        return add_planar_waveguides(sim, 1; width=width, index=index, direction=direction, x0=x0, y0=y0, background_index=background_index)
    end
end

"""
    new_sim = add_horizontal_planar_waveguides(sim; width, index, y0)

    new_sys = add_horizontal_planar_waveguides(sys; width, index, y0)
"""
function add_horizontal_planar_waveguides(sys::System; width, index, y0, background_index=1)
    sys = add_planar_waveguide(sys; width=width, index=index, direction=:h, x0=0, y0=y0, background_index=background_index)
    return sys
end
function add_horizontal_planar_waveguides(sim::Simulation; width, index, y0, background_index=sim.sys.domains[end][:n₁])
    sys = add_horizontal_planar_waveguides(sim.sys, width=width, index=index, y0=y0, background_index=background_index)
    return Simulation(sim; :sys => sys)
end

"""
    new_sim = add_vertical_planar_waveguides(sim; width, index, x0)

    new_sys = add_vertical_planar_waveguides(sys; width, index, x0)
"""
function add_vertical_planar_waveguides(sys::System; width::Number, index::Number, x0::Number, background_index::Number=1)
    sys = add_planar_waveguide(sys; width=width, index=index, direction=:v, x0=x0, y0=0, background_index=background_index)
    return sys
end
function add_vertical_planar_waveguides(sim::Simulation; width::Number, index::Number, x0::Number, background_index::Number=sim.sys.domains[end][:n₁])
    sys = add_vertical_planar_waveguides(sim.sys; width=width, index=index, x0=x0, background_index=background_index)
    return Simulation(sim; :sys => sys)
end

"""
    new_sim = add_planar_waveguide(sim, side::Symbol; index=1)

    new_sys = add_planar_waveguide(sys, side::Symbol; index=1)
"""
function add_planar_waveguide(sim::Simulation, side::Symbol; index::Number=1)
    throw(ErrorException("this is intended for metallic waveguides, haven't done yet"))
    if side == :top
        direction = :n
        x0 = sim.bnd.∂Ω[1,1]
        y0 = sim.bnd.∂Ω[2,2]
        width = sim.bnd.∂Ω[2,1]-sim.bnd.∂Ω[1,1]
    elseif side == :bottom
        direction = :s
        x0 = sim.bnd.∂Ω[1,1]
        y0 = sim.bnd.∂Ω[1,2]
        width = sim.bnd.∂Ω[2,1]-sim.bnd.∂Ω[1,1]
    elseif side == :left
        direction = :w
        x0 = sim.bnd.∂Ω[1,1]
        y0 = sim.bnd.∂Ω[1,2]
        width = sim.bnd.∂Ω[2,2]-sim.bnd.∂Ω[1,2]
    elseif side == :right
        direction = :e
        x0 = sim.bnd.∂Ω[2,1]
        y0 = sim.bnd.∂Ω[1,2]
        width = sim.bnd.∂Ω[2,2]-sim.bnd.∂Ω[1,2]
    else
        throw(ArgumentError("invalide side $side, must be one of :top, :bottom, :left, :right"))
    end
    return add_planar_waveguide(simsys; width=width, index=index, direction=direction, x0=x0, y0=y0)
end


text:
module JungleHelperClimbableOneWayPlatform

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/ClimbableOneWayPlatform" ClimbableOneWayPlatform(x::Integer, y::Integer, height::Integer=Maple.defaultBlockHeight, left::Bool=true, texture::String="wood", animationDelay::Number=0.0,
	surfaceIndex::Int=-1, staminaBehavior::String="None", sameDirectionJumpBoost::Bool=false, momentumJumpDecayTime::Number=0.0, momentumJumpDecayCurvature::Number=1.0)

textures = ["wood", "dream", "temple", "templeB", "cliffside", "reflection", "core", "moon"]
staminaBehaviors = ["None", "Conserve", "Regain"]

const placements = Ahorn.PlacementDict(
    "Climbable One-Way Platform (Left) (Jungle Helper)" => Ahorn.EntityPlacement(
        ClimbableOneWayPlatform,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => true
        )
    ),
    "Climbable One-Way Platform (Right) (Jungle Helper)" => Ahorn.EntityPlacement(
        ClimbableOneWayPlatform,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => false
        )
    )
)

quads = Tuple{Integer, Integer, Integer, Integer}[
    (0, 0, 8, 7) (8, 0, 8, 7) (16, 0, 8, 7);
    (0, 8, 8, 5) (8, 8, 8, 5) (16, 8, 8, 5)
]

Ahorn.editingOptions(entity::ClimbableOneWayPlatform) = Dict{String, Any}(
    "texture" => textures,
    "surfaceIndex" => Maple.tileset_sound_ids,
    "staminaBehavior" => staminaBehaviors
)

Ahorn.minimumSize(entity::ClimbableOneWayPlatform) = 0, 8
Ahorn.resizable(entity::ClimbableOneWayPlatform) = false, true

function Ahorn.selection(entity::ClimbableOneWayPlatform)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, 8, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ClimbableOneWayPlatform, room::Maple.Room)
    texture = get(entity.data, "texture", "wood")
    texture = texture == "default" ? "wood" : texture

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 8))
    left = get(entity.data, "left", true)

    startX = div(x, 8) + 1
    startY = div(y, 8) + 1
    stopY = startY + div(height, 8) - 1
    animated = Number(get(entity.data, "animationDelay", 0)) > 0

    Ahorn.Cairo.save(ctx)

    Ahorn.rotate(ctx, pi / 2)

    if left
        Ahorn.scale(ctx, 1, -1)
    end

    len = stopY - startY
    for i in 0:len
        if animated
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)00", 8 * i, left ? 0 : -8)
        else
            connected = false
            qx = 2
            if i == 0
                connected = get(room.fgTiles.data, (startY - 1, startX), false) != '0'
                qx = 1

            elseif i == len
                connected = get(room.fgTiles.data, (stopY + 1, startX), false) != '0'
                qx = 3
            end

            quad = quads[2 - connected, qx]
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)", 8 * i, left ? 0 : -8, quad...)
        end
    end

    Ahorn.Cairo.restore(ctx)
end

end


text:
# ----------- #
# Calibration #
# ----------- #

struct FlatCalibration <: AbstractDict{Symbol,Float64}
    d::OrderedDict{Symbol,Float64}
end

FlatCalibration(pairs::Pair{Symbol,Float64}...) =
    FlatCalibration(OrderedDict(pairs))

FlatCalibration() = FlatCalibration(OrderedDict{Symbol,Float64}())

struct GroupedCalibration <: AbstractDict{Symbol,Vector{Float64}}
    d::Dict{Symbol,Vector{Float64}}
end

GroupedCalibration() = GroupedCalibration(Dict{Symbol,Vector{Float64}}())

GroupedCalibration(pairs::Pair{Symbol,Vector{Float64}}...) =
    GroupedCalibration(Dict(pairs))

# use let block so `Calib` isn't a member of the Dolo module
let
    Calib = Union{FlatCalibration,GroupedCalibration}

    Base.getindex(c::Calib, n::Symbol) = c.d[n]
    Base.getindex(c::Calib, nms::Symbol...) = [c[n] for n in nms]

    # define the rest of the AbstractDict interface by forwarding to c.d

    @eval Base.iterate(c::$Calib, args...) = iterate(c.d, args...)

    # 1-arg functions
    for f in [:length, :keys, :values, :keytype, :valtype, :eltype, :isempty]
        @eval Base.$(f)(c::$(Calib)) = $(f)(c.d)
    end

    # 2-arg functions
    for f in [:haskey, :delete!, :pop!, :sizehint!]
        @eval Base.$(f)(c::$(Calib), arg) = $(f)(c.d, arg)
    end

    # 3-arg functions
    for f in [:get, :get!, :getkey, :pop!]
        @eval Base.$(f)(c::$(Calib), arg1, arg2) = $(f)(c.d, arg1, arg2)
    end

    # methods that didn't match signatures above
    Base.get!(f::Function, c::Calib, key) = get!(f, c.d, key)
    Base.get(f::Function, c::Calib, key) = get(f, c.d, key)
    Base.merge!(c::T, others::T...) where {T<:Calib} =
        T(merge!(c.d, [o.d for o in others]...))
    Base.copy(c::T) where {T<:Calib} = T(copy(c.d))
    Base.deepcopy(c::T) where {T<:Calib} = T(deepcopy(c.d))
end

# setting a single value
Base.setindex!(fc::FlatCalibration, v::Real, k::Symbol) =
    fc.d[k] = convert(Float64, v)

# setting multiple values
function Base.setindex!(fc::FlatCalibration, vs::Union{Tuple,AbstractVector},
                        ks::Symbol...)
    if length(vs) != length(ks)
        throw(DimensionMismatch("length of keys and values must be the same"))
    end

    for (v, k) in zip(vs, ks)
        fc[k] = v
    end
    fc
end

# setting a single group
function Base.setindex!(gc::GroupedCalibration, v::AbstractVector, k::Symbol)
    # use get with default of v in case we are adding a new key
    if length(v) != length(get(gc, k, v))
        msg = string("length $(length(v)) of new values for key $k does",
                     "not match length of existing value ($(length(gc[k])))")
        throw(DimensionMismatch(msg))
    end
    gc.d[k] = convert(Vector{Float64}, v)
end

# setting multiple groups
function Base.setindex!(gc::GroupedCalibration, vs::Tuple, ks::Symbol...)
    if length(vs) != length(ks)
        throw(DimensionMismatch("length of keys and values must be the same"))
    end

    for (v, k) in zip(vs, ks)
        gc[k] = v
    end
    gc
end

# Contains a FlatCalibration and GroupedCalibration, plus information on symbols
struct ModelCalibration
    flat::FlatCalibration
    grouped::GroupedCalibration
    symbol_table::Dict{Symbol,Tuple{Symbol,Int}}
    symbol_groups::OrderedDict{Symbol,Vector{Symbol}}
end

function ModelCalibration(calib::OrderedDict{Symbol,Real}, symbols::OrderedDict{Symbol,Array{Symbol,1}})
    flat = FlatCalibration(calib)
    grouped = GroupedCalibration()
    for (k, nms) in symbols
        grouped[k] = Float64[flat[nm] for nm in nms]
    end
    # grouped[:definitions] = Float64[flat[nm] for nm in keys(sm.definitions)]
    symbol_table = Dict{Symbol,Tuple{Symbol,Int}}()
    for (grp, vals) in symbols
        for (i, v) in enumerate(vals)
            symbol_table[v] = (grp, i)
        end
    end

    # for (i, v) in enumerate(keys(sm.definitions))
    #     symbol_table[v] = (:definitions, i)
    # end

    # # make sure we documented where in grouped every symbol is
    # @assert sort(collect(keys(symbol_table))) == sort(collect(keys(flat)))

    ModelCalibration(flat, grouped, symbol_table, deepcopy(symbols))
end

for f in (:copy, :deepcopy)
    @eval Base.$(f)(mc::ModelCalibration) =
        ModelCalibration($(f)(mc.flat), $(f)(mc.grouped),
                         $(f)(mc.symbol_table), $(f)(mc.symbol_groups))
end

# getindex and setindex! should be forwarded on to the grouped field
Base.getindex(mc::ModelCalibration, n::Symbol) = mc.grouped[n]
Base.getindex(mc::ModelCalibration, nms::Symbol...) = [mc[n] for n in nms]

function _setindex_group!(mc::ModelCalibration, v::AbstractVector, grp::Symbol)
    ks = mc.symbol_groups[grp]
    if length(v) != length(ks)
        msg = string("Calibration has $(length(ks)) symbols in $grp, ",
        "but passed $(length(v)) values")
        throw(DimensionMismatch(msg))
    end

    # update grouped
    mc.grouped[grp] = v

    # update flat
    for (vi, ki) in zip(v, ks)
        mc.flat[ki] = vi
    end
    mc
end

function _setindex_flat!(mc::ModelCalibration, v::Real, k::Symbol)
    # update flat
    mc.flat[k] = v

    # update grouped
    grp, ix = mc.symbol_table[k]
    mc.grouped[grp][ix] = v
    mc
end

function Base.setindex!(mc::ModelCalibration, v, k::Symbol)
    if haskey(mc.symbol_table, k)       # setting a variable
        _setindex_flat!(mc, v, k)
    elseif haskey(mc.symbol_groups, k)  # setting a group
        _setindex_group!(mc, v, k)
    else
        throw(KeyError("ModelCalibration has no variable or group named $k"))
    end
end

function Base.setindex!(mc::ModelCalibration, vs::Tuple, ks::Symbol...)
    if length(vs) != length(ks)
        throw(DimensionMismatch("length of keys and values must be the same"))
    end

    for (k, v) in zip(ks, vs)
        mc[k] = v
    end
    mc
end

# tries to replace a symbol if the key is in the calibration, otherwise just
# keeps the symbol in place
_replace_me(mc::ModelCalibration, s::Symbol) = get(mc.flat, s, s)
_replace_me(mc, o) = o

# eval_with will work on
function eval_with(mc::ModelCalibration, ex::Expr)
    # put in let block to allow us to define intermediates in expr and not
    # have them become globals in `Dolo`
    new_ex = MacroTools.prewalk(s->_replace_me(mc, s), ex)
    Core.eval(Dolo, :(
    let
        $new_ex
    end))
end

eval_with(mc::ModelCalibration, s::AbstractString) = eval_with(mc, _to_expr(s))
eval_with(mc::ModelCalibration, s::Symbol) = _replace_me(mc, s)
eval_with(mc::ModelCalibration, x::Number) = x
eval_with(mc::ModelCalibration, x::AbstractArray) = map(y->eval_with(mc, y), x)

function eval_with(mc::ModelCalibration, d::AbstractDict)
    out = Dict{Symbol,Any}()
    for (k, v) in d
        sk = Symbol(k)
        if sk == :tag
            out[sk] = v
        else
            out[sk] = eval_with(mc, v)
        end
    end
    out
end


text:
import MPI
include("02-broadcast-impl.jl")

function main()
    MPI.Init()

    do_broadcast()

    MPI.Finalize()
end

main()


text:
using Documenter, Literate, GXBeam

# Pre-install matplotlib
import Plots; Plots.pyplot()

# Generate examples
include("generate.jl")

# Build documentation
makedocs(;
    modules = [GXBeam],
    pages = [
        "Home" => "index.md",
        "Getting Started" => joinpath("examples", "guide.md"),
        "Examples" => [
            joinpath("examples", "cantilever.md"),
            joinpath("examples", "overdetermined.md"),
            joinpath("examples", "tipforce.md"),
            joinpath("examples", "tipmoment.md"),
            joinpath("examples", "curved.md"),
            joinpath("examples", "rotating.md"),
            joinpath("examples", "wind-turbine-blade.md"),
            joinpath("examples", "static-joined-wing.md"),
            joinpath("examples", "dynamic-joined-wing.md"),
            joinpath("examples", "vertical-axis-wind-turbine.md"),
        ],
        "Using GXBeam with DifferentialEquations.jl" => joinpath("examples", "diffeq.md"),
        "API Reference" => joinpath("reference", "reference.md"),
    ],
    sitename = "GXBeam.jl",
    authors = "Taylor McDonnell <taylormcd@byu.edu>",
)

deploydocs(
    repo = "github.com/byuflowlab/GXBeam.jl.git",
)


text:
function plotmodel(modelnames, project)
    data_to_plot = []
    key_labels =[]

	rpoB_eq_table = project["rpoB_equiv_table"]
	reads_table = project["read_counts_table"]
	hmm_database_info = project["hmm_database_info"]
	sample_names = map(string,names(project["rpoB_equiv_table"])[2:end])

    #Getting data out of the big table and putting it into a small data frame using the melt command to pivot that dataframe
    for (i, modelname) in enumerate(modelnames)
        single_model = melt(rpoB_eq_table[rpoB_eq_table[:PFAM_Model] .== modelname, :], :PFAM_Model)
        single_model_reads = melt(reads_table[reads_table[:PFAM_Model] .== modelname, :], :PFAM_Model)

        single_reads = []
        for readnum in single_model_reads[:value]
            push!(single_reads,string(readnum))
        end
        single_model[:reads] = single_reads
        single_model[:sample] = sample_names
		single_model[:PFAM_full] = hmm_database_info[hmm_database_info[:NAME] .== modelname, :DESC][1]
    relative_plot_position = 1:length(sample_names)
    relative_plot_position = relative_plot_position + 0.5 - i*(1/length(modelnames))
    single_model[:read_plot_pos] = relative_plot_position
        push!(data_to_plot, single_model)
    	push!(key_labels, hmm_database_info[hmm_database_info[:NAME] .== modelname, :DESC][1])
	end

    all_data_to_plot = data_to_plot[1]
    if length(data_to_plot) > 1
        for i in 2:length(data_to_plot)
            all_data_to_plot = vcat(all_data_to_plot, data_to_plot[i])
        end
    end

    #plotting
    modelplot = plot(all_data_to_plot, x=:value,y=:sample, colour=:PFAM_full, Geom.bar(position=:dodge,orientation=:horizontal),
    Guide.ylabel(""),
    Guide.xlabel("rpoB equivalents"),
    	Theme(bar_highlight=color(colorant"black"),
    	key_position=:bottom,
    	default_color=color(colorant"black"),
    	panel_stroke=color(colorant"black"),
    	grid_color=color(colorant"gray"),
    	major_label_font="Helvetica",
    	major_label_color=color(colorant"black"),
    	key_title_color=color(colorant"white"),
    	minor_label_font="Helvetica",
    	key_label_font="Helvetica",
    	minor_label_color=color(colorant"black"),
		point_label_font_size=6pt),
	Guide.annotation(compose(context(),
    text(all_data_to_plot[:value]+0.01,all_data_to_plot[:read_plot_pos],all_data_to_plot[:reads],[hleft])))
    )

	return modelplot
end

function perform_pca(project)

	rpoB_eq_table = project["rpoB_equiv_table"]
	sample_names = map(string,names(project["rpoB_equiv_table"])[2:end])

	rpoB_eq_matrix = convert(Array, rpoB_eq_table)
	rpoB_eq_matrix_converted = convert(Array{Float64,2},rpoB_eq_matrix[:,2:end])
	PCA_model = fit(PCA, rpoB_eq_matrix_converted)

	PC1 = transform(PCA_model, rpoB_eq_matrix_converted)[1,:]
	PC2 = transform(PCA_model, rpoB_eq_matrix_converted)[2,:]

	pca_plot = plot(x=PC1,y=PC2, Guide.annotation(compose(context(),text(PC1,PC2,map(string,sample_names),[hleft]))))

	return pca_plot
end

function plotcorrelation(modelnames,parameter,project)

	rpoB_eq_table = project["rpoB_equiv_table"]
	metadata_table = project["project_table"]
	hmm_database_info = project["hmm_database_info"]

	data_to_plot = []
    key_labels =[]

    #Getting data out of the big table and putting it into a small data frame using the melt command to pivot that dataframe
    for modelname in modelnames
        single_model = melt(rpoB_eq_table[rpoB_eq_table[:PFAM_Model] .== modelname, :], :PFAM_Model)
        single_model[:param] = reverse(metadata_table[parameter])
        single_model[:PFAM_full] = hmm_database_info[hmm_database_info[:NAME] .== modelname, :DESC][1]
        push!(data_to_plot, single_model)
	end

    #Putting all the data frames for each PFAM model together
    all_data_to_plot = data_to_plot[1]
    if length(data_to_plot) > 1
        for i in 2:length(data_to_plot)
            all_data_to_plot = vcat(all_data_to_plot, data_to_plot[i])
        end
    end

    #Plotting all of that data together
    plot(all_data_to_plot, x=:param,y=:value, color=:PFAM_full,
    Theme(default_color=color(colorant"black"),
    panel_stroke=color(colorant"black"),
    grid_color=color(colorant"gray"),
    major_label_font="Helvetica",
    major_label_color=color(colorant"black"),
    key_title_color=color(colorant"white"),
    minor_label_font="Helvetica",
    key_label_font="Helvetica",
    minor_label_color=color(colorant"black"),
    key_position=:bottom),
    Guide.xlabel(string(parameter)),
    Guide.ylabel("rpoB equivalents")
    #Guide.annotation(compose(context(), text(rpoB_eq,metadata_subset,sample_names)))
    )
end

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
module BPlot
    export plotb

    using PyPlot
    import Main.BMath

    function plotb(points::Vector; precision=10, Title="Math Plot", xlbl="X", ylbl="Y")
        xvals = Float64[]
        yvals = Float64[]
        for point in points
            if BMath.dim_count(point) != 1
                error("Cannot Plot Due to way to many dimensions! Needs to be 1")
            end
            c = BMath.cartesian(point)
            push!(xvals, round(c[1], digits=precision))
            push!(yvals, round(c[2], digits=precision))
        end

        pygui(true)
        scatter(xvals, yvals, 1)
        title(Title)
        xlabel(xlbl)
        ylabel(ylbl)
    end

    function plotb(inputs::Vector, points::Vector; precision=10, Title="Math Plot", xlbl="Variable", ylbl="X_Comp", zlbl="Y_Comp")
        xvals = Float64[]
        yvals = Float64[]
        zvals = Float64[]
        for i in 1:length(points)
            point = points[i]
            if BMath.dim_count(point) != 1
                error("Cannot Plot Due to way to many dimensions! Needs to be 1")
            end
            c = BMath.cartesian(point)
            push!(xvals, round(inputs[i], digits=precision))
            push!(yvals, round(c[1], digits=precision))
            push!(zvals, round(c[2], digits=precision))
        end

        pygui(true)
        using3D()
        fig = plt.figure()
        ax = fig.add_subplot(projection="3d")
        ax.set_title(Title)
        ax.set_xlabel(xlbl)
        ax.set_ylabel(ylbl)
        ax.set_zlabel(zlbl)
        ax.scatter(xvals, yvals, zvals, s=1)
    end
end


text:
module _EspEn2D
export EspEn2D
    """
        Esp2D,  = EspEn2D(Mat)

    Returns the bidimensional Espinosa entropy estimate (`Esp2D`)
    estimated for the data matrix (`Mat`) using the default parameters:
    time delay = 1, tolerance threshold = 20, percentage similarity = 0.7
    logarithm = natural, matrix template size = [floor(H/10) floor(W/10)]
    (where H and W represent the height (rows) and width (columns) of
    the data matrix `Mat`)
    ** The minimum number of rows and columns of `Mat` must be > 10.

        Esp2D = EspEn2D(Mat::AbstractArray{T,2} where T<:Real; m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
                                        tau::Int=1, r::Real=20, ps::Float=.7, Logx::Real=exp(1), Lock::Bool=true)

    Returns the bidimensional Espinosa entropy (`Esp2D`) estimates for the data
    matrix (`Mat`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Template submatrix dimensions, an integer scalar (i.e. the same
              height and width) or a two-element vector of integers [height, width] with a minimum value > 1.
              (default: [floor(H/10) floor(W/10)]) \n
    `tau`   - Time Delay, a positive integer       (default: 1) \n
    `r`     - Tolerance threshold, a positive scalar  (default: 20) \n
    `ps`    - Percentage similarity, a value in range [0 1]  (default: 0.7) \n
    `Logx`  - Logarithm base, a positive scalar    (default: natural)  \n
    `Lock`  - By default, EspEn2D only permits matrices with a maximum
              size of 128 x 128 to prevent memory errors when storing data
              on RAM. e.g. For Mat = [200 x 200], m = 3, and tau = 1,
              EspEn2D creates a vector of 753049836 elements.
              To enable matrices greater than [128 x 128] elements,
              set `Lock` to false.  (default: true)  \n
              `WARNING: unlocking the permitted matrix size may cause your Julia
              IDE to crash.`

    #  See also `SampEn2D`, `FuzzEn2D`, `DispEn2D`, `DistEn2D`, `PermEn2D`

    # References:
        [1] Ricardo Espinosa, et al.,
            "Two-dimensional EspEn: A New Approach to Analyze Image Texture
            by Irregularity."
            Entropy,
            23:1261 (2021)

    """
    function EspEn2D(Mat::AbstractArray{T,2} where T<:Real;
        m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
        tau::Int=1, r::Real=20, ps::Real=0.7, Logx::Real=exp(1), Lock::Bool=true)

    NL, NW = size(Mat)
    ((NL > 128 || NW > 128) && Lock) ?
        error("To prevent memory errors, matrix width & length must have <= 128 elements.
            To estimate EspEn2D  for the current matrix ($NL,$NW) change Lock to 'false'.
            Caution: unlocking the safe matrix size may cause the Julia IDE to crash.") :
            nothing

    length(m)==1 ?  (mL = m; mW = m) : (mL = m[1];  mW = m[2])
    (NL > 10 && NW > 10) ? nothing :
        error("Number of rows and columns in Mat must be > 10")
    (minimum(m)>1) ? nothing :
        error("m:     must be an integer > 1, or a 2 element tuple of integer values > 1")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (r >= 0) ?  nothing : error("r:     must be a positive value")
    ((ps >= 0) && (ps <= 1)) ?  nothing : error("ps:   must be a value in range [0 1]")
    (Logx>0) ? nothing :  error("Logx:     must be a positive number > 0")

    NL = NL - (mL-1)*tau
    NW = NW - (mW-1)*tau
    X = zeros(mL,mW,NL*NW)
    p = 0
    for k = 1:NL
        for n = 1:NW
            p += 1
            X[:,:,p] = Mat[k:tau:(mL-1)*tau+k,n:tau:(mW-1)*tau+n]
        end
    end

    p = size(X,3)
    p != NL*NW ? @warn("Potential error with submatrix division.") : nothing
    Ny = p*(p-1)/2
    Ny > 300000000 ? @warn("Number of pairwise distance calculations is $Ny") : nothing

    Cij = -ones(p-1,p-1)
    for k = 1:p-1
        Temp = abs.(X[:,:,k+1:p] .- X[:,:,k]) .<= r
        Cij[1:(end-k+1),k] = sum(Temp,dims=(1,2))
    end

    Dm = sum((Cij[:]/(mL*mW)).>=ps)/(p*(p-1)/2)
    Esp2D = -log(Logx, Dm)

    return Esp2D
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using Plots

f(x) = sin(x)

function trapezoidal_rule(a, b, n)
    h = (b - a) / n
    res = 0.0
    for i in 1:n
        if (i==0 || i==n)
            res += f(a+i*h)
        else
            res += 2 * f(a+i*h)
        end
    end

    res*h/2
end

println("Analytical=", 2)
println("Numerical=", trapezoidal_rule(0, pi, 100))


text:
using ArarForecast
using Test

@testset "ArarForecast.jl" begin
    # Write your tests here.
end


text:
using Test, PSDUtils, LegendTextIO, MJDSigGen, DSP, Distributions
using Random: seed!, MersenneTwister

@testset "PSDUtils.jl" begin
	@testset "todetcoords" begin
		x, y, z = rand(3) .* 50 .- 25
		zlen = rand() * 100

		@test todetcoords(x, y, z, zlen) == ((z + 200) * 10, x * 10, -10 * y + zlen / 2)

		@test todetcoords((x, y, z), zlen) == todetcoords(x, y, z, zlen)

		@test todetcoords(0, 5, -200, 100) == (0, 0, 0)

		event = read(DarioHitsFile("test.root.hits"))
		xtal_length = SigGenSetup("GWD6022_01ns.config").xtal_length

		pos = todetcoords.(event.pos, xtal_length)
		todetcoords!(event, xtal_length)

		@test event.pos == pos
	end

	@testset "get_signal!" begin
		setup = SigGenSetup("GWD6022_01ns.config")
		event = rand(collect(DarioHitsFile("test.root.hits")))

		signal = get_signal!(setup, event.pos, event.E)
		signal2 = get_signal!(setup, event)

		@test signal == signal2

		E = sum(event.E)

		@test all(-0.01 * E .< signal .< 1.01 * E)

		@test all(-0.01 * E .< diff(signal))

		@test signal[end] ≈ E

		@test abs(signal[1]) < 0.05 * E
	end

	@testset "charge_cloud_size" begin
		@test charge_cloud_size(1500) == -0.03499440089585633 + 0.0003359462486002238 * 1500

		@test charge_cloud_size(0) == 0.01
	end

	@testset "Signal Processing" begin
		signal = cumsum(rand(100))

		@test getA(signal) == maximum(diff(signal))

		let idx1 = rand(1:50), idx2 = rand(51:100), Δ = idx2 - idx1
			@test get_rise_time(signal, signal[idx1], signal[idx2]) == Δ

			@test get_rise_time(signal, signal[idx1], signal[idx2], 5.3) == 5.3Δ
		end

		@test get_rise_time(signal, signal[9], 101) == (100 + 1) - 9

		@test all(moving_average(signal, 5)[1:3] .≈ sum(signal[1:5]) / 5)
		@test moving_average(signal, 5)[79] ≈ sum(signal[79-2:79+2]) / 5
		@test all(moving_average(signal, 5)[99:100] .≈ sum(signal[96:100]) / 5)

		@test moving_average(signal, 9, 2) == moving_average(moving_average(signal, 9), 9)
	end

	@testset "apply_coll_effects" begin
		w = PSDUtils.gausswindow(0.53, 7, 0.13)

		@test w[1] == w[end] == exp(-0.5 * 7^2)

		@test w[cld(length(w), 2)] == maximum(w)

		@test length(w) == round(Int, 2 * 7 * 0.53 / 0.13)

		signal = cumsum(rand(1000))

		pad = zeros(Float64, 1099)

		PSDUtils.pad_signal!(pad, signal, 100)

		@test pad[1:1000] == signal

		@test all(pad[1001:1099] .== signal[end])

		padcoll   = apply_coll_effects(signal, 0.59, 0.1,  true; nσ = 5)
		nopadcoll = apply_coll_effects(signal, 0.59, 0.1, false; nσ = 5)

		wcoll = PSDUtils.gausswindow(0.59 / (2 * √(2 * log(2))), 5, 0.1)

		@test conv(signal, wcoll) == nopadcoll

		wlen = length(wcoll)
		padsignal = PSDUtils.pad_signal!(zeros(Float64, wlen + 1000 - 1), signal, wlen)

		@test conv(padsignal, wcoll)[1:1000 + wlen - 1] == padcoll
	end

	@testset "noise" begin
		let mt1 = MersenneTwister(123), mt2 = MersenneTwister(123)
			get_noisy_energy(mt1, 1500, 3) == rand(mt2, Normal(1500, 3))
		end

		let mt = MersenneTwister(538)
			seed!(538)
			get_noisy_energy(1500, 3) == rand(mt, Normal(1500, 3))
		end

		noise = randn(1000)

		let signal = cumsum(rand(100)), signal2 = copy(signal)
			addnoise!(signal, noise, 57)
			@test signal == signal2 .+ noise[57:57+99]
		end

		let signal = cumsum(rand(100)), signal2 = copy(signal)
			addnoise!(signal, noise, 970)
			@test signal == signal2 .+ [noise[970:end]; noise[1:69]]
		end

		let signal = cumsum(rand(100)), signal2 = copy(signal)
			addnoise!(MersenneTwister(11235), signal, noise)
			addnoise!(signal2, noise, rand(MersenneTwister(11235), 1:1000))
			@test signal == signal2
		end

		let signal = cumsum(rand(100)), signal2 = copy(signal)
			seed!(23579)
			addnoise!(signal, noise)
			addnoise!(signal2, noise, rand(MersenneTwister(23579), 1:1000))
			@test signal == signal2
		end

		@test_throws BoundsError addnoise!(rand(100), rand(1000), 9970)
	end
end


text:
using Documenter, ProjectiveVectors
using LinearAlgebra

makedocs(
    format = Documenter.HTML(),
    sitename = "ProjectiveVectors",
    pages = [
        "Index" => "index.md"
        ],
    doctest=false,
    modules=[ProjectiveVectors],
    checkdocs=:exports
)

deploydocs(
    repo   = "github.com/JuliaHomotopyContinuation/ProjectiveVectors.jl.git"
)


text:
module BSDiff

export bsdiff, bspatch, bsindex

using ArgTools
using SuffixArrays
using TranscodingStreams, CodecBzip2
using TranscodingStreams: Codec
using BufferedStreams

# abstract Patch format type
# specific formats defined below
abstract type Patch end

# control over compression details

compressor() = Bzip2Compressor(blocksize100k=9)
decompressor() = Bzip2Decompressor()

# specific format implementations

include("classic.jl")
include("endsley.jl")

# format names, patch types, auto detection

const DEFAULT_FORMAT = :classic
const FORMATS = Dict{Symbol,Type{<:Patch}}()
const MAGICS = Vector{Pair{Symbol,String}}()

function register_format!(format::Symbol, type::Type{<:Patch})
    magic = format_magic(type)
    FORMATS[format] = type
    push!(MAGICS, format => magic)
    sort!(MAGICS, by = ncodeunits∘last)
    return
end

register_format!(:classic, ClassicPatch)
register_format!(:endsley, EndsleyPatch)

function patch_type(format::Symbol)
    type = get(FORMATS, format, nothing)
    type !== nothing && return type
    throw(ArgumentError("unknown patch format: $format"))
end

function detect_format(patch_io::IO)
    data = UInt8[]
    for (format, magic) in MAGICS
        n = ncodeunits(magic)
        m = n - length(data)
        m > 0 && append!(data, read(patch_io, m))
        view(data, 1:n) == codeunits(magic) && return format
    end
    return :unknown
end
detect_format(path::AbstractString) = open(detect_format, path)

const INDEX_HEADER = "SUFFIX ARRAY\0"

## high-level API (similar to the C tool) ##

"""
    bsdiff(old, new, [ patch ]; format = [ :classic | :endsley ]) -> patch

Compute a binary patch that will transform the content of `old` into the content
of `new`. All arguments can be strings or IO handles. If no `patch` argument is
provided, the patch data is written to a temporary file whose path is returned.

The `old` argument can also be a 2-tuple of strings and/or IO handles, in which
case the first is used as the old data and the second is used as a precomputed
index of the old data, as computed by [`bsindex`](@ref). Since indexing the old
data is the slowest part of generating a diff, precomputing this and reusing it
can significantly speed up generting diffs from the same old file to multiple
different new files.

The `format` keyword argument allows selecting a patch format to generate. The
value must be one of the symbols `:classic` or `:endsley` indicating a bsdiff
patch format. The classic patch format is generated by default, but the Endsley
format can be selected with `bsdiff(old, new, patch, format = :endsley)`.
"""
function bsdiff(
    old::Union{ArgRead, NTuple{2, ArgRead}},
    new::ArgRead,
    patch::Union{ArgWrite, Nothing} = nothing;
    format::Symbol = DEFAULT_FORMAT,
)
    type = patch_type(format)
    old_data, index = data_and_index(old)
    new_data = arg_read(read, new)
    arg_write(patch) do patch_io
        write(patch_io, format_magic(type))
        patch_obj = write_start(type, patch_io, old_data, new_data)
        generate_patch(patch_obj, old_data, new_data, index)
        write_finish(patch_obj)
    end
end

"""
    bspatch(old, [ new, ] patch; format = [ :classic | :endsley ]) -> new

Apply a binary patch given by the `patch` argument to the content of `old` to
produce the content of `new`. All arguments can be strings or IO handles. If no
`new` argument is provided, the new data is written to a temporary file whose
path is returned.

Note that the optional argument is the middle argument, which is a bit unusual
but makes the argument order when passing all three paths consistent with the
`bspatch` command and with the `bsdiff` function.

By default `bspatch` auto-detects the patch format, so the `format` keyword
argument is usually unnecessary. If you wish to restrict the format of patch
that will be accepted, however, you can use this keyword argument: `bspatch`
will raise an error unless the patch file has indicated format.
"""
function bspatch(
    old::ArgRead,
    new::Union{ArgWrite, Nothing},
    patch::ArgRead;
    format::Symbol = :auto,
)
    format == :auto || format in keys(FORMATS) ||
        error("unknown patch format: $format")
    old_data = arg_read(read, old)
    arg_read(patch) do patch_io
        detected = detect_format(patch_io)
        detected == :unknown && error("unrecognized/corrupt patch file")
        format == :auto || format == detected ||
            error("patch has $detected format, expected $format format")
        type = patch_type(detected)
        patch_obj = read_start(type, patch_io)
        arg_write(new) do new_io
            new_io = BufferedOutputStream(new_io)
            apply_patch(patch_obj, old_data, new_io)
            finalize(patch_obj)
            flush(new_io)
        end
    end
end

function bspatch(
    old::ArgRead,
    patch::ArgRead;
    format::Symbol = :auto,
)
    bspatch(old, nothing, patch; format = format)
end

"""
    bsindex(old, [ index ]) -> index

Save index data (a sorted suffix array) for the content of `old` into `index`.
All arguments can be strings or IO handles. If no `index` argument is provided,
the index data is saved to a temporary file whose path is returned.

The index can be passed to `bsdiff` to speed up the diff computation by passing
`(old, index)` as the first argument instead of just `old`. Since indexing the
old data is the slowest part of generating a diff, precomputing this and reusing
it can significantly speed up generting diffs from the same old file to multiple
different new files.
"""
function bsindex(
    old::ArgRead,
    index::Union{ArgWrite, Nothing} = nothing,
)
    old_data = arg_read(read, old)
    arg_write(index) do index_io
        write(index_io, INDEX_HEADER)
        index_data = generate_index(old_data)
        write(index_io, UInt8(sizeof(eltype(index_data))))
        write(index_io, index_data)
    end
end

## loading data and index ##

const IndexType{T<:Integer} = Vector{T}

function data_and_index(data_path::ArgRead)
    data = arg_read(read, data_path)
    data, generate_index(data)
end

function data_and_index(
    (data_path, index_path)::NTuple{2, ArgRead},
)
    data = arg_read(read, data_path)
    index = arg_read(index_path) do index_io
        hdr = String(read(index_io, ncodeunits(INDEX_HEADER)))
        hdr == INDEX_HEADER || error("corrupt bsdiff index file")
        unit = Int(read(index_io, UInt8))
        T = unit == 1 ? UInt8 :
            unit == 2 ? UInt16 :
            unit == 4 ? UInt32 :
            unit == 8 ? UInt64 :
            error("invalid unit size for bsdiff index file: $unit")
        read!(index_io, Vector{T}(undef, length(data)))
    end
    return data, index
end

## generic patch generation and application logic ##

generate_index(data::AbstractVector{<:UInt8}) = suffixsort(data, 0)

# transform used to serialize integers to avoid lots of
# high bytes being emitted for small negative values
int_io(x::Signed) = ifelse(x == abs(x), x, typemin(x) - x)
write_int(io::IO, x::Signed) = write(io, int_io(Int64(x)))
read_int(io::IO) = Int(int_io(read(io, Int64)))

"""
Return lexicographic order and length of common prefix.
"""
function strcmplen(p::Ptr{UInt8}, m::Int, q::Ptr{UInt8}, n::Int)
    i = 0
    while i < min(m, n)
        a = unsafe_load(p + i)
        b = unsafe_load(q + i)
        a ≠ b && return (a - b) % Int8, i
        i += 1
    end
    return (m - n) % Int8, i
end

"""
Search for the longest prefix of new[t:end] in old.
Uses the suffix array of old to search efficiently.
"""
function prefix_search(
    index::IndexType, # suffix array
    old::AbstractVector{UInt8}, # old data to search in
    new::AbstractVector{UInt8}, # new data to search for
    t::Int, # search for longest match of new[t:end]
)
    old_n = length(old)
    new_n = length(new) - t + 1
    old_p = pointer(old)
    new_p = pointer(new, t)
    # invariant: longest match is in index[lo:hi]
    lo, hi = 1, old_n
    c = lo_c = hi_c = 0
    while hi - lo ≥ 2
        m = (lo + hi) >>> 1
        s = index[m]
        x, l = strcmplen(new_p+c, new_n+c, old_p+s+c, old_n-s-c)
        if 0 < x
            lo, lo_c = m, c+l
        else
            hi, hi_c = m, c+l
        end
        c = min(lo_c, hi_c)
    end
    lo_c > hi_c ? (index[lo]+1, lo_c) : (index[hi]+1, hi_c)
end

"""
Computes and emits the diff of the byte vectors `new` versus `old`.
The `index` array is a zero-based suffix array of `old`.
"""
function generate_patch(
    patch::Patch,
    old::AbstractVector{UInt8},
    new::AbstractVector{UInt8},
    index::IndexType = generate_index(old),
)
    oldsize, newsize = length(old), length(new)
    scan = len = pos = lastscan = lastpos = lastoffset = 0

    while scan < newsize
        oldscore = 0
        scsc = scan += len
        while scan < newsize
            pos, len = prefix_search(index, old, new, scan+1)
            pos -= 1 # zero-based
            while scsc < scan + len
                oldscore += scsc + lastoffset < oldsize &&
                    old[scsc + lastoffset + 1] == new[scsc + 1]
                scsc += 1
            end
            if len == oldscore && len ≠ 0 || len > oldscore + 8
                break
            end
            oldscore -= scan + lastoffset < oldsize &&
                old[scan + lastoffset + 1] == new[scan + 1]
            scan += 1
        end
        if len ≠ oldscore || scan == newsize
            i = s = Sf = lenf = 0
            while lastscan + i < scan && lastpos + i < oldsize
                s += old[lastpos + i + 1] == new[lastscan + i + 1]
                i += 1
                if 2s - i > 2Sf - lenf
                    Sf = s
                    lenf = i
                end
            end
            lenb = 0
            if scan < newsize
                s = Sb = 0
                i = 1
                while scan ≥ lastscan + i && pos ≥ i
                    s += old[pos - i + 1] == new[scan - i + 1]
                    if 2s - i > 2Sb - lenb
                        Sb = s
                        lenb = i
                    end
                    i += 1
                end
            end
            if lastscan + lenf > scan - lenb
                overlap = (lastscan + lenf) - (scan - lenb)
                i = s = Ss = lens = 0
                while i < overlap
                    s += new[lastscan + lenf - overlap + i + 1] ==
                         old[lastpos + lenf - overlap + i + 1]
                    s -= new[scan - lenb + i + 1] ==
                         old[pos - lenb + i + 1]
                    if s > Ss
                        Ss = s
                        lens = i + 1;
                    end
                    i += 1
                end
                lenf += lens - overlap
                lenb -= lens
            end

            diff_size = lenf
            copy_size = (scan - lenb) - (lastscan + lenf)
            skip_size = (pos - lenb) - (lastpos + lenf)

            # skip if both blocks are empty
            diff_size == copy_size == 0 && continue

            encode_control(patch, diff_size, copy_size, skip_size)
            encode_diff(patch, diff_size, new, lastscan, old, lastpos)
            encode_data(patch, copy_size, new, lastscan + diff_size)

            lastscan = scan - lenb
            lastpos = pos - lenb
            lastoffset = pos - scan
        end
    end
end

"""
Apply a patch stream to the `old` data buffer, emitting a `new` data stream.
"""
function apply_patch(
    patch::Patch,
    old::AbstractVector{UInt8},
    new::IO,
    new_size::Int = hasfield(typeof(patch), :new_size) ? patch.new_size : typemax(Int),
)
    old_pos = new_pos = 0
    old_size = length(old)
    while true
        ctrl = decode_control(patch)
        ctrl == nothing && break
        diff_size, copy_size, skip_size = ctrl

        # sanity checks
        0 ≤ diff_size && 0 ≤ copy_size &&                # block sizes are non-negative
        new_pos + diff_size + copy_size ≤ new_size &&    # don't write > new_size bytes
        0 ≤ old_pos && old_pos + diff_size ≤ old_size || # bounds check for old data
            error("corrupt bsdiff patch")

        decode_diff(patch, diff_size, new, old, old_pos)
        decode_data(patch, copy_size, new)

        new_pos += diff_size + copy_size
        old_pos += diff_size + skip_size
    end
    return new_pos
end

end # module


text:
function solve_modelKARC(nlp_stop, X :: PDataK, α:: T) where T
    # target value should be close to satisfy αλ=||d||
    target = [ ( abs( α*X.shifts[i] - X.norm_dirs[i]^X.τ) )   for i = 1 : X.nshifts ];

    # pick the closest shift to the target within positive definite H+λI
    X.indmin = max(1,X.indmin)
    start = X.indmin
    bidon,indmin = findmin(target[X.positives[start:end]])
    X.indmin = start + indmin - 1
    p_imin = X.positives[X.indmin]
    X.d = X.xShift[:,p_imin]
    X.λ = X.shifts[p_imin]

    return X.d, NaN * rand(length(X.d)), X.λ
end