text:
# ------------------------------------------
#   solver

using IPFitting, ACE1, LinearAlgebra, JuLIP

export solver_params

"""
`solver_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct one of the solvers.
All parameters are passed as keyword argument and the kind of
parameters required depend on "type".

### LSQR Parameters
* `type = "lsqr"`
* `lsqr_damp = 5e-3`
* `lsqr_atol = 1e-6`
* `lsqr_colim = 1e8`
* `lsqr_maxiter = 1e5`
* `lsqr_verbose = false`

### RRQR Parameters
* `type = "rrqr"`
* `rrqr_tol = 1e-5`

### BBR
* `type = "bbr"`
* `brr_n_iter = 300`
* `brr_tol = 1e-3`

### ARD
* `type = "ard"`
* `ard_n_iter = 300`
* `ard_tol = 1e-3`
* `ard_threshold_lambda = 1e4`

"""
function solver_params(; type = nothing, kwargs...)
    # TODO error message
    solver_type = _solver_to_params(type)
    @assert solver_type in keys(_solvers_params)
    return _solvers_params[solver_type](; kwargs...)
end

"""
`lsqr_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct a lsqr solver.
All parameters are passed as keyword argument.

### Parameters
* `lsqr_damp = 5e-3`
* `lsqr_atol = 1e-6`
* `lsqr_colim = 1e8`
* `lsqr_maxiter = 1e5`
* `lsqr_verbose = false`
"""
lsqr_params(; lsqr_damp = 5e-3, lsqr_atol = 1e-6, lsqr_conlim = 1e8, lsqr_maxiter = Integer(1e5), lsqr_verbose = false) =
    Dict("type" => "lsqr", "lsqr_damp" => lsqr_damp, "lsqr_atol" => lsqr_atol,
         "lsqr_conlim" => lsqr_conlim, "lsqr_maxiter" => lsqr_maxiter, "lsqr_verbose" => lsqr_verbose)

"""
`rrqr_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct a rrqr solver.
All parameters are passed as keyword argument.

### Parameters
* `rrqr_tol = 1e-5`
"""
rrqr_params(; rrqr_tol = 1e-5) =
    Dict("type" => "rrqr", "rrqr_tol" => rrqr_tol)

"""
`brr_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct a scikit-learn
Bayesian ridge regression solver. All parameters are passed as
keyword argument.

### Parameters
* `brr_n_iter = 300`
* `brr_tol = 1e-3`
"""
brr_params(; brr_n_iter = 300, brr_tol = 1e-3) =
    Dict("type" => "brr", "brr_n_iter" => brr_n_iter, "brr_tol" => brr_tol)

"""
`ard_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct a scikit-learn
Automatic Relevance Detemrination solver. All parameters are
passed as keyword argument.

### Parameters
* `ard_n_iter = 300`
* `ard_tol = 1e-3`
* `ard_threshold_lambda = 1e4`
"""
ard_params(; ard_n_iter = 300, ard_tol = 1e-3, ard_threshold_lambda = 1e4) =
    Dict("type" => "ard", "ard_n_iter" => ard_n_iter, "ard_tol" => ard_tol, "ard_threshold_lambda" => ard_threshold_lambda)

_solver_to_params(solver_type::Union{Symbol, AbstractString}) =
    string(solver_type)

_solvers_params = Dict("lsqr" => lsqr_params, "rrqr" => rrqr_params,
                        "brr" => brr_params, "ard" => ard_params)

_params_to_solver(solver_type::AbstractString) =
    Symbol.(solver_type)

function generate_solver(params::Dict)
    params = copy(params)
    params["solver"] = _params_to_solver(pop!(params, "type"))
    return params
end


text:
import CategoricalArrays: recode!

"""
    recode!(dat, col, recode_lookup, typ)
This function performs the re-coding specified by the `recode_lookup`
dictionary and modifies the `dat` argument in place. Note that we must
specify the type of the column using the `typ` argument. This should be
the type of recoded values.
"""

function recode!(dat::DataTable, col::Symbol, recode_lookup::Dict, typ::DataType)
    tmp = dat[col]
    n = size(dat, 1)
    dat[col] = NullableArray{typ, 1}(n)
    for i = 1:n
        if tmp[i] ∈ keys(recode_lookup)
            dat[i, col] = recode_lookup[tmp[i]]
        else
            dat[i, col] = tmp[i]
        end
    end
end

"""
    clean_string(x)
Given a single string `x`, this function returns a copy of the string
where all non-alphabetic characters have been replaced with a blank space.
The exception is the hyphen, which is not removed or altered.
"""
function clean_string(x::T) where {T <: AbstractString}
    xlwr = lowercase(x)
    xcln = replace(xlwr, r"[^-a-z]", " ")     # keep only letters and hyphens
    return xcln
end

function cut_stopwords(x_arr, stopwords)
    res = [z for z in x_arr if z ∉ stopwords]
    return res
end

function clean_all_strings(x::NullableArray{T, 1}, stopwords::Array{String, 1}) where {T <: AbstractString}
    n = length(x)
    res = Array{String,1}(n)
    for i = 1:n
        xi_cln = clean_string(x[i].value)
        res[i] = join(cut_stopwords(split(xi_cln), stopwords), " ")
    end
    res
end

function clean_all_strings(x::Array{T, 1}, stopwords::Array{String, 1}) where {T <: AbstractString}
    n = length(x)
    res = Array{String,1}(n)
    for i = 1:n
        xi_cln = clean_string(x[i])
        res[i] = join(cut_stopwords(split(xi_cln), stopwords), " ")
    end
    res
end

"""
    document_frequency(x)
Give a vector of strings, `x`, this function return a dictionary where the
keys are the unique words in the entire vector `x` and the values are integers
representing the number of documents that contained the given word.
"""
function document_frequency(x::NullableArray{String, 1})
    n = length(x)
    res = Dict{String, Int}()
    for i = 1:n
        unq_words = unique(split(x[i].value))
        for word in unq_words
            res[word] = get(res, word, 0) + 1
        end
    end
    res
end

function document_frequency(x::Array{String, 1})
    n = length(x)
    res = Dict{String, Int}()
    for i = 1:n
        unq_words = unique(split(x[i]))
        for word in unq_words
            res[word] = get(res, word, 0) + 1
        end
    end
    res
end

"""
    tf_idf_datatable(x, stopwords)
Given a vector (of strings, typically), `x`, this function returns a
DataTable where columns contain the TF/IDF values for the words appearing
in the original vector `x`. In particular, the rows in the output DataTable
correspond to the indices of the input vector `x`. And the columns correspond
to a sorted list of unique terms in the the vocabulary of the input vector (after
the stop words have been removed).
"""
function tf_idf_datatable(x::AbstractArray, stopwords::Array{String,1})
     x_cln = clean_all_strings(x, stopwords)
     doc_freq = document_frequency(x_cln)
     unq_words = collect(keys(doc_freq))
     sort!(unq_words)

     n = length(x)
     p = length(unq_words)

     word_idx = Dict{String, Int}()
     for j = 1:p
         word_idx[unq_words[j]] = j
     end

     mat = zeros(n, p)
     for i = 1:n
         word_cnts = countmap(split(x_cln[i]))
         for (word, cnt) in word_cnts
             j = word_idx[word]
             mat[i, j] = cnt * log(n / doc_freq[word])
         end
     end
     res = convert(DataTable, mat)
     col_names = map(Symbol, unq_words)
     names!(res, col_names)
     res
 end

ext: jl


text:
"""
Calculates the Hamming distance (numbers of different bits) between two numbers.
"""
hamming_distance(x::N1, y::N2) where {N1<:Integer, N2<:Integer} = count_ones(xor(x,y))

"""
A node in the B-K tree.
"""
mutable struct Node{T}
    item::Union{Nothing, T}
    children::Dict{Float64, Node{T}}
end

# Constructors
Node(item::T) where T = Node{T}(item, Dict{Float64, Node{T}}())

Node{T}() where T = Node{T}(nothing, Dict{Float64, Node{T}}())

"""
The B-K tree structure.
"""
mutable struct BKTree{T}
    f::Function     # distance function: df(x,y) should return the distance
    root::Node{T}    # node
end

# Constructors
# Remarks:
#   If no function `f` is provided, the `hamming_distance` will
#   be used however, the hamming distance works only for ::Integer data
BKTree{T}(f::Function) where {T} = BKTree(f, Node{T}())

BKTree(root::Node{T}) where {T} = BKTree(hamming_distance, root)

BKTree(item::T) where {T} = BKTree(hamming_distance, Node(item))

BKTree{T}() where T = BKTree(Node{T}())

BKTree(f::Function, items::Vector{T}) where {T} = begin
    bkt = BKTree{T}(f)
    for item in items
        add!(bkt, item)
    end
    return bkt
end

BKTree(items::Vector{T}) where {T} = begin
    bkt = BKTree{T}()
    for item in items
        add!(bkt, item)
    end
    return bkt
end

# Show methods
show(io::IO, bkt::BKTree{T}) where {T} =
    print(io, "B-K Tree{$T}, function=$(bkt.f), ",
          "$(length(bkt.root.children)) branches")

show(io::IO, node::Node{T}) where T = begin
    if node.item != nothing
        print(io, "B-K Node, item=$(node.item), $(length(node.children)) branches")
    else
        print(io, "B-K Node, no items.")
    end
end

# Useful functions
"""
Determines whether a Node in the BKTree is empty or not.
"""
is_empty_node(node::Node) = node.item == nothing ? true : false

"""
Adds a node (i.e. `Node{T}(item::T)`) to the tree.
"""
function add!(tree::BKTree{T}, item::T) where T
    if is_empty_node(tree.root)
        tree.root = Node(item)
    else
        node = tree.root
        while true
            parent, children = node.item, node.children
            distance = tree.f(item, parent)
            node = get(children, distance, nothing)
            if node == nothing
                children[distance] = Node(item)
                break
            end
        end
    end
    return tree
end

"""
Find items in the tree whose distance is less or equal to
`n` and returns the top `k` items, ordered ascending according
to the distance.
"""
function find(tree::BKTree{T}, item::T, n::Int; k::Int=1) where {T}
    found = Vector{Tuple{Float64,T}}()
    candidates = Deque{Node{T}}()
    push!(candidates, tree.root)
    while !isempty(candidates)
        candidate = popfirst!(candidates)
        distance = tree.f(candidate.item, item)
        if distance <= n
            push!(found, (distance, candidate.item))
        end
        if !isempty(candidate.children)
            lower = distance - n
            upper = distance + n
            for (d, c) in candidate.children
                if lower <= d <= upper
                    push!(candidates, c)
                end
            end
        end
    end
    # Sort by distance
    k = min(k, length(found))
    sfound = sort(found, by=x->x[1], rev=false)[1:k]
    return sfound
end


text:
using JuMP
using CPLEX
include("parser.jl")
include("exercice1.jl")
include("exercice2.jl")
include("exercice3.jl")
include("exercice4.jl")
include("exercice5.jl")

timee = Array{Float64,2}(zeros(10,5))
for i in 10000:10000:50000
  size = i
  j = i÷10000
  time_1 = @timed Benders("instance1.txt", size) # exo 1
  timee[j,1] = time_1.time

  time_2 = @timed Normal("instance1.txt", size) # exo 2
  timee[j,2] = time_2.time

  time_3 = @timed AutoBenders("instance1.txt", size) # exo 3
  timee[j,3] = time_3.time

  time_4 = @timed FastBenders("instance1.txt", size) # exo 4
  timee[j,4] = time_4.time

  time_5 = @timed FastBendersGurobi("instance1.txt", size) # exo 5
  timee[j,5] = time_5.time
end

for j in 1:5
  println("size of instance = ", j*10000)
  println(timee[j,1], " seconds for exercise 1")
  println(timee[j,2], " seconds for exercise 2")
  println(timee[j,3], " seconds for exercise 3")
  println(timee[j,4], " seconds for exercise 4")
  println(timee[j,5], " seconds for exercise 5")
end

max_stars_repo_name: boumpteryx/MLA

max_issues_repo_name: boumpteryx/MLA

max_forks_repo_name: boumpteryx/MLA


text:
using BenchmarkTools, SymbolicUtils
using SymbolicUtils: is_literal_number

using Random

SUITE = BenchmarkGroup()

@syms a b c d; Random.seed!(123);

let r = @rule(~x => ~x), rs = RuleSet([r]),
    acr = @rule(~x::is_literal_number + ~y => ~y)

    overhead = SUITE["overhead"]  = BenchmarkGroup()
    overhead["rule"]  = BenchmarkGroup()

    overhead["rule"]["noop:Int"]  = @benchmarkable $r(1)
    overhead["rule"]["noop:Sym"]  = @benchmarkable $r($a)
    overhead["rule"]["noop:Term"] = @benchmarkable $r($(a+2))

    overhead["acrule"]  = BenchmarkGroup()
    overhead["acrule"]["noop:Int"]  = @benchmarkable $acr(1)
    overhead["acrule"]["noop:Sym"]  = @benchmarkable $acr($a)
    overhead["acrule"]["a+2"] = @benchmarkable $acr($(a+2))
    overhead["acrule"]["a+b"] = @benchmarkable $acr($(a+b))
    overhead["acrule"]["a+2+b"] = @benchmarkable $acr($(a+2+b))

    overhead["ruleset"]  = BenchmarkGroup()
    overhead["ruleset"]["noop:Int"]  = @benchmarkable $rs(1)
    overhead["ruleset"]["noop:Sym"]  = @benchmarkable $rs($a)
    overhead["ruleset"]["noop:Term"] = @benchmarkable $rs($(a+2))

    overhead["simplify"]  = BenchmarkGroup()
    overhead["simplify"]["noop:Int"]  = @benchmarkable simplify(1)
    overhead["simplify"]["noop:Sym"]  = @benchmarkable simplify($a)
    overhead["simplify"]["noop:Term"] = @benchmarkable simplify($(a+2))

    function random_term(len; atoms, funs, fallback_atom=1)
        xs = rand(atoms, len)
        while length(xs) > 1
            xs = map(Iterators.partition(xs, 2)) do xy
                x = xy[1]; y = get(xy, 2, fallback_atom)
                rand(funs)(x, y)
            end
        end
        xs[]
    end
    ex1 = random_term(1000, atoms=[a, b, c, d, a^(-1), b^(-1), 1, 2.0], funs=[+, *])
    ex2 = random_term(1000, atoms=[a, b, c, d, a^(-1), b^(-1), 1, 2.0], funs=[/, *])

    overhead["simplify"]["randterm (+, *):serial"] = @benchmarkable simplify($ex1, threaded=false)
    overhead["simplify"]["randterm (/, *):serial"] = @benchmarkable simplify($ex2, threaded=false)
    overhead["simplify"]["randterm (+, *):thread"] = @benchmarkable simplify($ex1, threaded=true)
    overhead["simplify"]["randterm (/, *):thread"] = @benchmarkable simplify($ex2, threaded=true)

    overhead["substitute"] = BenchmarkGroup()

    overhead["substitute"]["a"] = @benchmarkable substitute(subs_expr, $(Dict(a=>1))) setup=begin
        subs_expr = (sin(a+b) + cos(b+c)) * (sin(b+c) + cos(c+a)) * (sin(c+a) + cos(a+b))
    end

    overhead["substitute"]["a,b"] = @benchmarkable substitute(subs_expr, $(Dict(a=>1, b=>2))) setup=begin
        subs_expr = (sin(a+b) + cos(b+c)) * (sin(b+c) + cos(c+a)) * (sin(c+a) + cos(a+b))
    end

    overhead["substitute"]["a,b,c"] = @benchmarkable substitute(subs_expr, $(Dict(a=>1, b=>2, c=>3))) setup=begin
        subs_expr = (sin(a+b) + cos(b+c)) * (sin(b+c) + cos(c+a)) * (sin(c+a) + cos(a+b))
    end

end

let
    pform = SUITE["polyform"]  = BenchmarkGroup()
    @syms a b c d e f g h i
    ex = (f + ((((g*(c^2)*(e^2)) / d - e*h*(c^2)) / b + (-c*e*f*g) / d + c*e*i) /
              (i + ((c*e*g) / d - c*h) / b + (-f*g) / d) - c*e) / b +
         ((g*(f^2)) / d + ((-c*e*f*g) / d + c*f*h) / b - f*i) /
         (i + ((c*e*g) / d - c*h) / b + (-f*g) / d)) / d

    o = (d + (e*((c*(g + (-d*g) / d)) / (i + (-c*(h + (-e*g) / d)) / b + (-f*g) / d))) / b +
         (-f*(g + (-d*g) / d)) / (i + (-c*(h + (-e*g) / d)) / b + (-f*g) / d)) / d
    pform["simplify_fractions"] = @benchmarkable simplify_fractions($ex)
    pform["iszero"] = @benchmarkable SymbolicUtils.fraction_iszero($ex)
    pform["isone"] = @benchmarkable SymbolicUtils.fraction_isone($o)
    pform["easy_iszero"] = @benchmarkable SymbolicUtils.fraction_iszero($((b*(h + (-e*g) / d)) / b + (e*g) / d - h))
end


text:
#### fWn-MLE ####

# This files contains a version of MLE for a single scale fWn.
# For legacy only, all functions can be realized via fWn bank MLE.

"""
    fWn_log_likelihood_H(X, H, ψ, G)

Log-likelihood of a fWn model with the optimal volatility.
"""
function fWn_log_likelihood_H(X::AbstractVecOrMat{<:Real}, H::Real, ψ::AbstractVector{<:Real}, G::AbstractVector{<:Integer})
    @assert length(G) == size(X,1)
    proc = FractionalWaveletNoise(H, ψ)
    return log_likelihood_H(covmat(proc, G), X)

    # Σ::AbstractMatrix = if length(G)>0
    #     @assert length(G) == size(X,1)
    #     covmat(proc, G)
    # else
    #     covmat(proc, size(X,1))
    # end
    # return log_likelihood_H(Σ, X)
end

"""
Maximum likelihood estimation of Hurst exponent and volatility for fractional Wavelet noise.

# Args
- X: sample vector or matrix. For matrix each column is a sample.
- ψ: wavelet filter used for computing `X`.
- G: integer time grid of `X`, by default the regular grid `0:size(X,1)-1` is used.
- method: `:optim` for optimization based or `:table` for lookup table based procedure

# Returns
- H, σ, L, obj: estimation of Hurst exponent, of volatility, log-likelihood, object of optimizer

# Notes
- The MLE is known for its sensitivity to mis-specification of model, as well as to missing value (NaN) and outliers.
- The starting point of the grid `G` has no importance since fWn is stationary.
"""
function fWn_MLE_estim(X::AbstractVecOrMat{<:Real}, ψ::AbstractVector{<:Real}, G::AbstractVector{<:Integer}; method::Symbol=:optim)
    if length(X) == 0 || any(isnan.(X))  # for empty input or input containing nans
        return (hurst=NaN, σ=NaN, loglikelihood=NaN, optimizer=nothing)
    else
        # Force the zero-mean condition: this might be dangerous and gives lower estimations
        # X = (ndims(X) == 1) ? X : X .- mean(X, dims=2)

        @assert length(G) == size(X,1)  "Mismatched dimensions."
        @assert minimum(abs.(diff(sort(G)))) > 0  "All elements of the grid must be distinct."
        # if length(G)>0
        #     @assert length(G) == size(X,1)
        #     @assert minimum(abs.(diff(sort(G)))) > 0  # all elements are distinct
        # end

        func = h -> -fWn_log_likelihood_H(X, h, ψ, G)

        opm = nothing
        hurst = nothing
        ε::Real=1e-2   # search hurst in the range [ε, 1-ε]

        if method == :optim
            # Gradient-free constrained optimization
            opm = Optim.optimize(func, ε, 1-ε, Optim.Brent())
            # # Gradient-based optimization
            # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
            # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
            hurst = Optim.minimizer(opm)[1]
        elseif method == :table
            Hs = collect(ε:ε:1-ε)
            hurst = Hs[argmin([func(h) for h in Hs])]
        else
            throw("Unknown method: ", method)
        end

        proc = FractionalWaveletNoise(hurst, ψ)

        # Σ = covmat(proc, length(G)>0 ? G : size(X,1))
        Σ = covmat(proc, G)

        # Estimation of volatility
        σ = sqrt(xiAx(Σ, X) / length(X))
        # σ *= mean(diff(G))^(hurst)  # <- why this?

        # Log-likelihood
        L = log_likelihood_H(Σ, X)
        # # or equivalently
        # L = log_likelihood(σ^2*Σ, X)

        L /= length(X)  # normalization
        return (hurst=hurst, σ=σ, loglikelihood=L, optimizer=opm)
    end
end

"""
Accelerated fWn-MLE by dividing a large vector of samples into smaller ones.

The MLE method can be computationally expensive on data of large dimensions due to the inversion of covariance matrix. This function accelerates the MLE method by applying rolling vectorization on the large vector `X` with the parameter `(s,u,l)`. In this way the original vector is divided into smaller vectors which are then treated by MLE as i.i.d. samples.

# Args
- X: sample path of a fWn.
- ψ: wavelet filter used for computing `X`.
- s: sub window size
- u: downsampling factor
- l: length of decorrelation

# Notes
The rolling vectorization returns
- empty, if `s>size(X)[end]`
- `X` in matrix form (i.e. same as `reshape(X, :, 1)`), if `s==size(X)[end]`
"""
function fWn_MLE_estim(X::AbstractVector{<:Real}, ψ::AbstractVector{<:Real}, s::Integer, u::Integer, l::Integer; kwargs...)
    t, V = rolling_vectorize(X, s, u, l; mode=:causal)

    fWn_MLE_estim(V, ψ, u*(0:s-1); kwargs...)
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
# Optional flat/non-resource MVC folder structure
# push!(LOAD_PATH,  "controllers", "views", "views/layouts", "models")


text:
# Torsion QuadraticForm
#
# Example:
# A = matrix(ZZ, [[2,0,0,-1],[0,2,0,-1],[0,0,2,-1],[-1,-1,-1,2]])
# L = Zlattice(gram = A)
# T = Hecke.discriminant_group(T)

# We representation torsion quadratic modules as quotients of Z-lattices
# by a full rank sublattice.
#
# We store them as a Z-lattice M together with a projection p : M -> A
# onto an abelian group A. The bilinear structure of A is induced via p,
# that is <a, b> = <p^-1(a), p^-1(a)> with values in Q/nZ, where n
# is the modulus and depends on the kernel of p.
#
# Elements of A are basically just elements of the underlying abelian group.
# To move between M and A, we use the lift function lift : M -> A
# and coercion A(m).
#
# N.B. Since there are no elements of Z-latties, we think of elements of M as
# elements of the ambient vector space. Thus if v::Vector is such an element
# then the coordinates with respec to the basis of M are given by
# v * inv(basis_matrix(M)).
mutable struct TorQuadMod
  ab_grp::GrpAbFinGen             # underlying abelian group
  cover::ZLat                     # ZLat -> ab_grp, x -> x * proj
  proj::fmpz_mat                  # is a projection and respects the forms
  gens_lift::Vector{Vector{fmpz}}
  gens_lift_ambient::Vector{Vector{fmpq}}
  gens_lift_mat::fmpz_mat          # integer matrix
  gens_lift_mat_ambient::fmpq_mat
  d::fmpz
  rels::fmpz_mat
  modulus::fmpq
  modulus_qf::fmpq
  value_module::QmodnZ
  value_module_qf::QmodnZ
  gram_matrix_bilinear::fmpq_mat
  gram_matrix_quadratic::fmpq_mat

  TorQuadMod() = new()
end

################################################################################
#
#  Construction
#
################################################################################

# compute the torsion quadratic module M/N
function torsion_quadratic_module(M::ZLat, N::ZLat)
  @req ambient_space(M) === ambient_space(N) "Lattices must have same ambient space"
  _rels = basis_matrix(N) * inv(basis_matrix(M))
  @req isone(denominator(_rels)) "Second lattice must be a submodule of first lattice"
  rels = change_base_ring(FlintZZ, _rels)
  A = abelian_group(rels)
  S, mS = snf(A)
  gens_lift = [collect(mS(s).coeff) for s in gens(S)]

  num = basis_matrix(M) * gram_matrix(ambient_space(M)) * basis_matrix(N)'
  modulus = reduce(gcd, [a for a in num], init = zero(fmpq))
  norm = reduce(gcd, diagonal(gram_matrix(N)), init = zero(fmpq))
  modulus_qf = gcd(norm, 2 * modulus)

  T = TorQuadMod()
  T.cover = M
  T.ab_grp = S
  T.proj = inv(mS).map
  T.gens_lift = gens_lift
  T.gens_lift_mat = matrix(ZZ, length(gens_lift), ngens(A), reduce(vcat, gens_lift))
  T.gens_lift_mat_ambient = change_base_ring(FlintQQ, T.gens_lift_mat) * basis_matrix(M)
  T.modulus = modulus
  T.modulus_qf = modulus_qf
  T.value_module = QmodnZ(modulus)
  T.value_module_qf = QmodnZ(modulus_qf)
  return T
end

# compute M^#/M
function discriminant_group(L::ZLat)
  # I need to check that M is integral
  return torsion_quadratic_module(dual(L), L)
end

@doc Markdown.doc"""
    order(T::TorQuadMod) -> fmpz

Return the order of `T`
"""
function order(T::TorQuadMod)
  return order(abelian_group(T))
end

@doc Markdown.doc"""
    exponent(T::TorQuadMod) -> fmpz

Returns the exponent of `T`
"""
function exponent(T::TorQuadMod)
  return exponent(abelian_group(T))
end

################################################################################
#
#  Basic field access
#
################################################################################

abelian_group(T::TorQuadMod) = T.ab_grp

cover(T::TorQuadMod) = T.cover

value_module(T::TorQuadMod) = T.value_module

value_module_quadratic_form(T::TorQuadMod) = T.value_module_qf

################################################################################
#
#  Gram matrices
#
################################################################################

function gram_matrix_bilinear(T::TorQuadMod)
  if isdefined(T, :gram_matrix_bilinear)
    return T.gram_matrix_bilinear
  end
  g = gens(T)
  G = zero_matrix(FlintQQ, length(g), length(g))
  for i in 1:length(g)
    for j in 1:i
      G[i, j] = G[j, i] = lift(g[i] * g[j])
    end
  end
  T.gram_matrix_bilinear = G
  return G
end

function gram_matrix_quadratic(T::TorQuadMod)
  if isdefined(T, :gram_matrix_quadratic)
    return T.gram_matrix_quadratic
  end
  g = gens(T)
  r = length(g)
  G = zero_matrix(FlintQQ, r, r)
  for i in 1:r
    for j in 1:(i - 1)
      G[i, j] = G[j, i] = lift(g[i] * g[j])
    end
    G[i, i] = lift(quadratic_product(g[i]))
  end
  T.gram_matrix_quadratic = G
  return G
end

################################################################################
#
#  I/O
#
################################################################################

function Base.show(io::IO, T::TorQuadMod)
  print(io, "Finite quadratic module over Integer Ring with invariants ")
  println(io, elementary_divisors(abelian_group(T)))
  print(io, "Gram matrix of the quadratic form with values in ")
  println(io, value_module_quadratic_form(T))
  print(io, gram_matrix_quadratic(T))
end

################################################################################
#
#  Elements
#
################################################################################

mutable struct TorQuadModElem
  a::GrpAbFinGenElem
  parent::TorQuadMod

  TorQuadModElem(T::TorQuadMod, a::GrpAbFinGenElem) = new(a, T)
end

# TODO: Check the parents ...
(T::TorQuadMod)(a::GrpAbFinGenElem) = TorQuadModElem(T, a)

function (T::TorQuadMod)(v::Vector{fmpq})
  vv = change_base_ring(FlintZZ, matrix(FlintQQ, 1, length(v), v) * inv(basis_matrix(cover(T))))
  return T(abelian_group(T)(vv * T.proj))
end

# TODO: Cache this
gens(T::TorQuadMod) = [T(g) for g in gens(abelian_group(T))]

parent(a::TorQuadModElem) = a.parent

################################################################################
#
#  Inner product
#
################################################################################

function Base.:(*)(a::TorQuadModElem, b::TorQuadModElem)
  T = parent(a)
  z = inner_product(ambient_space(cover(T)), lift(a), lift(b))
  return value_module(T)(z)
end

################################################################################
#
#  Quadratic product
#
################################################################################

function quadratic_product(a::TorQuadModElem)
  T = parent(a)
  al = lift(a)
  z = inner_product(ambient_space(cover(T)), al, al)
  return value_module_quadratic_form(T)(z)
end

################################################################################
#
#  Lift
#
################################################################################

# Lift an element to the ambient space of cover(parent(a))
function lift(a::TorQuadModElem)
  T = parent(a)
  z = change_base_ring(FlintQQ, a.a.coeff) * T.gens_lift_mat_ambient
  return fmpq[z[1, i] for i in 1:ncols(z)]
end

# this is broken
#function TorQuadMod(q::fmpq_mat)
#  @req issquare(q) "Matrix must be a square matrix"
#  @req issymmetric(q) "Matrix must be symmetric"
#
#  d = denominator(q)
#  Q = change_base_ring(FlintZZ, d * q)
#  S, U, V = snf_with_transform(Q)
#  D = change_base_ring(FlintQQ, U) * q * change_base_ring(FlintQQ, V)
#  @show D
#  L = Zlattice(gram = d^2 * q)
#  denoms = [denominator(D[i, i]) for i in 1:ncols(D)]
#  rels = diagonal_matrix(denoms) * U
#  _A = abelian_group(rels)
#  S, T = snf(_A)
#  @show T
#  @show T(gens(S)[1])
#  value_module = QmodnZ()
#  return TorQuadMod(S, L, d, rels, value_module)
#end

#        if modulus is None or check:
#           # The inner product of two elements `b(v1+W,v2+W)`
#           # is defined `mod (V,W)`
#           num = V.basis_matrix() * V.inner_product_matrix() * W.basis_matrix().T
#           max_modulus = gcd(num.list())
#
#       if modulus is None:
#           modulus = max_modulus
#       elif check and max_modulus / modulus not in V.base_ring():
#           raise ValueError("the modulus must divide (V, W)")
#
#       if modulus_qf is None or check:
#           # The quadratic_product of an element `q(v+W)` is defined
#           # `\mod 2(V,W) + ZZ\{ (w,w) | w in w\}`
#           norm = gcd(W.gram_matrix().diagonal())
#           max_modulus_qf = gcd(norm, 2 * modulus)
#
#       if modulus_qf is None:
#           modulus_qf = max_modulus_qf
#       elif check and max_modulus_qf / modulus_qf not in V.base_ring():
#           raise ValueError("the modulus_qf must divide (V, W)")
#       return super(TorsionQuadraticModule, cls).__classcall__(cls, V, W, gens, modulus, modulus_qf)


text:
using .TurbulenceClosures: ▶xy_cfa, ▶xy_fca

#####
##### Functions for non-rotating models
#####

@inline x_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT)
@inline y_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT)
@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT)

#####
##### The 'FPlane' approximation. This is equivalent to a model with a constant
##### rotation rate around its vertical axis.
#####

"""
    FPlane{FT} <: AbstractRotation

A parameter object for constant rotation around a vertical axis.
"""
struct FPlane{FT} <: AbstractRotation
    f :: FT
end

"""
    FPlane([FT=Float64;] f=nothing, rotation_rate=nothing, latitude=nothing)

Returns a parameter object for constant rotation at the angular frequency
`f/2`, and therefore with background vorticity `f`, around a vertical axis.
If `f` is not specified, it is calculated from `rotation_rate` and
`latitude` according to the relation `f = 2*rotation_rate*sind(latitude).

Also called `FPlane`, after the "f-plane" approximation for the local effect of
Earth's rotation in a planar coordinate system tangent to the Earth's surface.
"""
function FPlane(FT=Float64; f=nothing, rotation_rate=nothing, latitude=nothing)

    if f == nothing && rotation_rate != nothing && latitude != nothing
        return FPlane{FT}(2rotation_rate*sind(latitude))
    elseif f != nothing && rotation_rate == nothing && latitude == nothing
        return FPlane{FT}(f)
    else
        throw(ArgumentError("Either both keywords rotation_rate and
                             latitude must be specified, *or* only f
                             must be specified."))
    end
end

@inline x_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = - coriolis.f * ▶xy_fca(i, j, k, grid, U.v)
@inline y_f_cross_U(i, j, k, grid, coriolis::FPlane, U) =   coriolis.f * ▶xy_cfa(i, j, k, grid, U.u)
@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::FPlane, U) where FT = zero(FT)

#####
##### The Beta Plane
#####

"""
    BetaPlane{T} <: AbstractRotation

A parameter object for meridionally increasing Coriolis
parameter (`f = f₀ + βy`).
"""
struct BetaPlane{T} <: AbstractRotation
    f₀ :: T
     β :: T
end

"""
    BetaPlane([T=Float64;] f₀=nothing, β=nothing,
                           rotation_rate=nothing, latitude=nothing, radius=nothing)

A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`).

The user may specify both `f₀` and `β`, or the three parameters
`rotation_rate`, `latitude`, and `radius` that specify the rotation rate and radius
of a planet, and the central latitude at which the `β`-plane approximation is to be made.
"""
function BetaPlane(T=Float64; f₀=nothing, β=nothing,
                              rotation_rate=nothing, latitude=nothing, radius=nothing)

    f_and_β = f₀ != nothing && β != nothing
    planet_parameters = rotation_rate != nothing && latitude != nothing && radius != nothing

    (f_and_β && all(Tuple(p === nothing for p in (rotation_rate, latitude, radius)))) ||
    (planet_parameters && all(Tuple(p === nothing for p in (f₀, β)))) ||
        throw(ArgumentError("Either both keywords f₀ and β must be specified,
                            *or* all of rotation_rate, latitude, and radius."))

    if planet_parameters
        f₀ = 2rotation_rate * sind(latitude)
         β = 2rotation_rate * cosd(latitude) / radius
     end

    return BetaPlane{T}(f₀, β)
end

@inline x_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) =
    @inbounds - (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, U.v)
@inline y_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) =
    @inbounds   (coriolis.f₀ + coriolis.β * grid.yF[j]) * ▶xy_cfa(i, j, k, grid, U.u)
@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::BetaPlane, U) where FT = zero(FT)


text:
function repair()
	sol = Bool[1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1]
	cost,matrix = loadSPP("instances/pb_200rnd0900.dat")
	println("avant repair : z = ", transpose(sol) * cost)

	constr_violated = Int64[]
	cst = matrix*sol
	for i = 1:length(cst)
		if cst[i] > 1
			push!(constr_violated, i)
		end
	end
	# println(constr_violated)

	var_1 = Int64[]
	for i = 1:length(sol)
		if sol[i] == 1
			push!(var_1, i)
		end
	end
	# println(var_1)
	var_1_conflit = Int64[]
	for i = 1:length(sol)
		if sol[i] == 1
			for cst in constr_violated
				if matrix[cst,i] == 1
					push!(var_1_conflit, i)
				end
			end
		end
	end
	# println(var_1_conflit)
	while length(constr_violated) != 0
		# println("--------------------")
		inu = inutilite(var_1_conflit, constr_violated, cost, matrix)
		# println("inutility = ", inu)
		e = getBestInutil(var_1_conflit, inu)
		# println("variable a mettre a zero :", e)
		sol[e] = 0

		# ind = findfirst((x -> x == e),var_1_conflit)
		# deleteat!(var_1_conflit,ind)

		var_1_conflit = Int64[]
		for i = 1:length(sol)
			if sol[i] == 1
				for cst in constr_violated
					if matrix[cst,i] == 1
						push!(var_1_conflit, i)
					end
				end
			end
		end

		# println("var_1_conflit = ", var_1_conflit)
		constr_violated = updateCstV(constr_violated, sol, matrix)
		# println("constr_violated", constr_violated)
	end
	println("apres repair : z = ", transpose(sol) * cost)
	return sol
end

function getBestInutil(var_1::Array{Int64,1}, inutility::Array{Float64,1})
	max_x = 0
	max_inu = 0.0
	for j in 1:length(var_1)
		if inutility[j] > max_inu
			max_inu = inutility[j]
			max_x = var_1[j]
		end
	end
	return max_x
end

function inutilite(var_1, constr, cost, matrix)
	inu = zeros(Float64, length(var_1))
	i = 1
	for var in var_1
		sum = 0
		for cst in constr
			sum = sum + matrix[cst,var]
		end
		inu[i] = sum / cost[i]
		i = i + 1
	end
	return inu
end

function updateCstV(constr_violated, sol, matrix)
	# res = copy(constr_violated)
	# for constr in constr_violated
	# 	if transpose(matrix[:,constr])*sol <= 1
	# 		ind = findfirst((x -> x == constr),res)
	# 		deleteat!(res,ind)
	# 	end
	# end
	# return res
	constr_violated = Int64[]
	cst = matrix*sol
	for i = 1:length(cst)
		if cst[i] > 1
			push!(constr_violated, i)
		end
	end
	return constr_violated
end


text:
__precompile__()

module IRKGaussLegendre

using Reexport
@reexport using DiffEqBase

using LinearAlgebra
using Parameters
using OrdinaryDiffEq
using RecursiveArrayTools

const CompiledFloats = Union{Float32,Float64}

include("IRKCoefficients.jl")
include("IRKGL16AuxFunctions.jl")
include("IRKGL16Solver.jl")

export IRKGL16,IRKAlgorithm
export tcoeffs, CompiledFloats

end # module


text:
struct IntDisjointSet{T<:Integer}
    parents::Vector{T}
    normalized::Ref{Bool}
end

IntDisjointSet{T}() where {T<:Integer} = IntDisjointSet{T}(Vector{T}[], Ref(true))
Base.length(x::IntDisjointSet) = length(x.parents)

function Base.push!(x::IntDisjointSet{T}) where {T}
    push!(x.parents, convert(T, -1))
    convert(T, length(x))
end

function find_root(x::IntDisjointSet{T}, i::T) where {T}
    while x.parents[i] >= 0
        i = x.parents[i]
    end
    return convert(T, i)
end

function in_same_set(x::IntDisjointSet{T}, a::T, b::T) where {T}
    find_root(x, a) == find_root(x, b)
end

function Base.union!(x::IntDisjointSet{T}, i::T, j::T) where {T}
    pi = find_root(x, i)
    pj = find_root(x, j)
    if pi != pj
        x.normalized[] = false
        isize = -x.parents[pi]
        jsize = -x.parents[pj]
        if isize > jsize # swap to make size of i less than j
            pi, pj = pj, pi
            isize, jsize = jsize, isize
        end
        x.parents[pj] -= isize # increase new size of pj
        x.parents[pi] = pj # set parent of pi to pj
    end
    return convert(T, pj)
end

function normalize!(x::IntDisjointSet{T}) where {T}
    for i in convert(T, length(x))
        pi = find_root(x, i)
        if pi != i
            x.parents[i] = convert(T, pi)
        end
    end
    x.normalized[] = true
end

# If normalized we don't even need a loop here.
function _find_root_normal(x::IntDisjointSet{T}, i::T) where {T}
    pi = x.parents[i]
    if pi < 0 # Is `i` a root?
        return i
    else
        return pi
    end
    # return pi
end

function _in_same_set_normal(x::IntDisjointSet, a::Int64, b::Int64)
    _find_root_normal(x, a) == _find_root_normal(x, b)
end

function find_root_if_normal(x::IntDisjointSet, i::Int64)
    if x.normalized[]
        _find_root_normal(x, i)
    else
        find_root(x, i)
    end
end


text:
include("gc.jl")
include("test.jl")
include("compile.jl")
include("examples.jl")


text:
doc"""
    GeneralizedExtremeValue(μ, σ, ξ)

The *Generalized extreme value distribution* with shape parameter `ξ`, scale `σ` and location `μ` has probability density function

$f(x; \xi, \sigma, \mu) = \begin{cases}
        \frac{1}{\sigma} \left[ 1+\left(\frac{x-\mu}{\sigma}\right)\xi\right]^{-1/\xi-1} \exp\left\{-\left[ 1+ \left(\frac{x-\mu}{\sigma}\right)\xi\right]^{-1/\xi} \right\} & \text{for } \xi \neq 0 \\
        \frac{1}{\sigma} \exp\left\{-\frac{x-\mu}{\sigma}\right\} \exp\left\{-\exp\left[-\frac{x-\mu}{\sigma}\right]\right\} & \text{for } \xi = 0
    \end{cases}$

for

$x \in \begin{cases}
        \left[ \mu - \frac{\sigma}{\xi}, + \infty \right) & \text{for } \xi > 0 \\
        \left( - \infty, + \infty \right) & \text{for } \xi = 0 \\
        \left( - \infty, \mu - \frac{\sigma}{\xi} \right] & \text{for } \xi < 0
    \end{cases}$

```julia
GeneralizedExtremeValue(k, s, m)      # Generalized Pareto distribution with shape k, scale s and location m.

params(d)       # Get the parameters, i.e. (k, s, m)
shape(d)        # Get the shape parameter, i.e. k (sometimes called c)
scale(d)        # Get the scale parameter, i.e. s
location(d)     # Get the location parameter, i.e. m
```

External links

* [Generalized extreme value distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution)

"""

immutable GeneralizedExtremeValue{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    ξ::T

    function GeneralizedExtremeValue(μ::T, σ::T, ξ::T)
        σ > zero(σ) || error("Scale must be positive")
        new(μ, σ, ξ)
    end
end

GeneralizedExtremeValue{T<:Real}(μ::T, σ::T, ξ::T) = GeneralizedExtremeValue{T}(μ, σ, ξ)
GeneralizedExtremeValue(μ::Real, σ::Real, ξ::Real) = GeneralizedExtremeValue(promote(μ, σ, ξ)...)
function GeneralizedExtremeValue(μ::Integer, σ::Integer, ξ::Integer)
    GeneralizedExtremeValue(Float64(μ), Float64(σ), Float64(ξ))
end

#### Conversions
function convert{T<:Real}(::Type{GeneralizedExtremeValue{T}}, μ::Real, σ::Real, ξ::Real)
    GeneralizedExtremeValue(T(μ), T(σ), T(ξ))
end
function convert{T <: Real, S <: Real}(::Type{GeneralizedExtremeValue{T}}, d::GeneralizedExtremeValue{S})
    GeneralizedExtremeValue(T(d.μ), T(d.σ), T(d.ξ))
end

minimum{T<:Real}(d::GeneralizedExtremeValue{T}) =
        d.ξ > 0 ? d.μ - d.σ / d.ξ : -T(Inf)
maximum{T<:Real}(d::GeneralizedExtremeValue{T}) =
        d.ξ < 0 ? d.μ - d.σ / d.ξ : T(Inf)

#### Parameters

shape(d::GeneralizedExtremeValue) = d.ξ
scale(d::GeneralizedExtremeValue) = d.σ
location(d::GeneralizedExtremeValue) = d.μ
params(d::GeneralizedExtremeValue) = (d.μ, d.σ, d.ξ)
@inline partype{T<:Real}(d::GeneralizedExtremeValue{T}) = T

#### Statistics

testfd(d::GeneralizedExtremeValue) = d.ξ^3
g(d::GeneralizedExtremeValue, k::Real) = gamma(1 - k * d.ξ) # This should not be exported.

function median(d::GeneralizedExtremeValue)
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps() # ξ == 0
        return μ - σ * log(log(2))
    else
        return μ + σ * (log(2) ^ (- ξ) - 1) / ξ
    end
end

function mean{T<:Real}(d::GeneralizedExtremeValue{T})
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        return μ + σ * γ
    elseif ξ < 1
        return μ + σ * (gamma(1 - ξ) - 1) / ξ
    else
        return T(Inf)
    end
end

function mode(d::GeneralizedExtremeValue)
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        return μ
    else
        return μ + σ * ((1 + ξ) ^ (-ξ) - 1) / ξ
    end
end

function var{T<:Real}(d::GeneralizedExtremeValue{T})
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        return σ ^2 * π^2 / 6
    elseif ξ < 1/2
        return σ^2 * (g(d, 2) - g(d, 1) ^ 2) / ξ^2
    else
        return T(Inf)
    end
end

function skewness{T<:Real}(d::GeneralizedExtremeValue{T})
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        return 12sqrt(6) * zeta(3) / pi ^ 3 * one(T)
    elseif ξ < 1/3
        g1 = g(d, 1)
        g2 = g(d, 2)
        g3 = g(d, 3)
        return sign(ξ) * (g3 - 3g1 * g2 + 2g1^3) / (g2 - g1^2) ^ (3/2)
    else
        return T(Inf)
    end
end

function kurtosis{T<:Real}(d::GeneralizedExtremeValue{T})
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        return T(12)/5
    elseif ξ < 1 / 4
        g1 = g(d, 1)
        g2 = g(d, 2)
        g3 = g(d, 3)
        g4 = g(d, 4)
        return (g4 - 4g1 * g3 + 6g2 * g1^2 - 3 * g1^4) / (g2 - g1^2)^2 - 3
    else
        return T(Inf)
    end
end

function entropy(d::GeneralizedExtremeValue)
    (μ, σ, ξ) = params(d)
    return log(σ) + γ * ξ + (1 + γ)
end

function quantile(d::GeneralizedExtremeValue, p::Real)
    (μ, σ, ξ) = params(d)

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        return μ + σ * (-log(-log(p)))
    else
        return μ + σ * ((-log(p))^(-ξ) - 1) / ξ
    end
end

#### Support

insupport(d::GeneralizedExtremeValue, x::Real) = minimum(d) <= x <= maximum(d)

#### Evaluation

function logpdf{T<:Real}(d::GeneralizedExtremeValue{T}, x::Real)
    if x == -Inf || x == Inf || ! insupport(d, x)
      return -T(Inf)
    else
        (μ, σ, ξ) = params(d)

        z = (x - μ) / σ # Normalise x.
        if abs(ξ) < eps(one(ξ)) # ξ == 0
            t = z
            return -log(σ) - t - exp(-t)
        else
            if z * ξ == -1 # Otherwise, would compute zero to the power something.
                return -T(Inf)
            else
                t = (1 + z * ξ) ^ (-1/ξ)
                return - log(σ) + (ξ + 1) * log(t) - t
            end
        end
    end
end

function pdf{T<:Real}(d::GeneralizedExtremeValue{T}, x::Real)
    if x == -Inf || x == Inf || ! insupport(d, x)
        return zero(T)
    else
        (μ, σ, ξ) = params(d)

        z = (x - μ) / σ # Normalise x.
        if abs(ξ) < eps(one(ξ)) # ξ == 0
            t = exp(-z)
            return (t * exp(-t)) / σ
        else
            if z * ξ == -1 # In this case: zero to the power something.
                return zero(T)
            else
                t = (1 + z*ξ)^(- 1 / ξ)
                return (t^(ξ + 1) * exp(- t)) / σ
            end
        end
    end
end

function logcdf{T<:Real}(d::GeneralizedExtremeValue{T}, x::Real)
    if insupport(d, x)
        (μ, σ, ξ) = params(d)

        z = (x - μ) / σ # Normalise x.
        if abs(ξ) < eps(one(ξ)) # ξ == 0
            return -exp(- z)
        else
            return - (1 + z * ξ) ^ ( -1/ξ)
        end
    elseif x <= minimum(d)
        return -T(Inf)
    else
        return zero(T)
    end
end

function cdf{T<:Real}(d::GeneralizedExtremeValue{T}, x::Real)
    if insupport(d, x)
        (μ, σ, ξ) = params(d)

        z = (x - μ) / σ # Normalise x.
        if abs(ξ) < eps(one(ξ)) # ξ == 0
            t = exp(-z)
        else
            t = (1 + z * ξ) ^ (-1/ξ)
        end
        return exp(- t)
    elseif x <= minimum(d)
        return zero(T)
    else
        return one(T)
    end
end

logccdf(d::GeneralizedExtremeValue, x::Real) = log1p(- cdf(d, x))
ccdf(d::GeneralizedExtremeValue, x::Real) = - expm1(logcdf(d, x))

#### Sampling

function rand(d::GeneralizedExtremeValue)
    (μ, σ, ξ) = params(d)

    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand()

    if abs(ξ) < eps(one(ξ)) # ξ == 0
        rd = - log(- log(u))
    else
        rd = expm1(- ξ * log(- log(u))) / ξ
    end

    return μ + σ * rd
end


text:
struct DPModel{U<:ObjectiveFunction}
    variables               ::Array{String}
    objective               ::U
    constraints             ::Dict{Union{DataType,UnionAll}, Array{Any}}
    domains                 ::Dict{String, Union{IntervalTreeNode,Nothing}}
end

"""
    deepCopy(domains::Union{Dict{String, IntervalTreeNode}, Dict{String, Union{Nothing, IntervalTreeNode}}})

Make a deep copy of a domain object

# Arguments
- `domains::Union{Dict{String, IntervalTreeNode}, Dict{String, Union{Nothing, IntervalTreeNode}}}`: the domains to copy
"""
function deepCopy(domains::Union{Dict{String, IntervalTreeNode}, Dict{String, Union{Nothing, IntervalTreeNode}}})
    newDomains = Dict{String, Union{IntervalTreeNode,Nothing}}()
    for (key, value) in domains
        if !isnothing(value)
            newDomains[key] = deepCopyIntervalTreeNodes(value, false)
        end
    end
    return newDomains
end


text:
using NLOptControl
using VehicleModels
using Parameters
using DataFrames
using MAVs

try
  close(x.s1)
  close(x.s2)
end
mdl=0;n=0;r=0;params=0;c=0;x=0;d=0;pa=0;s=0;

# THINGS TO THINK ABOUT TODO
# 6) intialize obstacle field in intial solve
# 7) get the vehicle states from MATLAB so that I can plot things as well

c=defineCase(;(:mode=>:caseStudy));

#Ok version
#setMisc!(c;activeSafety=true,followPath=true,followDriver=false,PredictX0=false,FixedTp=false,NF=0,tp=4.0,tex=0.4,max_cpu_time=0.35,Ni=2,Nck=[10,10]);
#setWeights!(c;sr=0.05,path=10.0,driver=0.0)

#without path following
#setMisc!(c;activeSafety=true,followPath=false,followDriver=false,PredictX0=false,FixedTp=false,NF=0,tp=3.5,tex=0.4,max_cpu_time=0.4,Ni=2,Nck=[10,10]);
#setWeights!(c;sr=0.05,path=10.0,driver=5.0)

#4.23 working with path following
setMisc!(c;activeSafety=true,followPath=true,followDriver=false,PredictX0=false,FixedTp=false,NF=0,tp=3.5,tex=0.4,max_cpu_time=0.4,Ni=2,Nck=[10,10]);
setWeights!(c;sr=0.05,path=10.0,driver=5.0)
# 4.21
#setMisc!(c;activeSafety=true,followPath=false,followDriver=false,PredictX0=false,FixedTp=false,NF=0,tp=3.0,tex=0.5,max_cpu_time=0.45,Ni=3,Nck=[12,10,8]);
#setWeights!(c;sr=0.05,path=5.0,driver=0.0)

mdl,n,r,params,x,d=initializeSharedControl!(c)

# setup UDP port and UDP communication
x.s1 = UDPSocket();bind(x.s1,ip"73.113.149.41",36880); # change this ip to the server (or julia code)
x.s2 = UDPSocket();bind(x.s2,ip"73.113.149.41",12000); # change this ip to the server (or julia code)

global pa=params[1];
#global s=Settings(;reset=false,save=true,simulate=true,MPC=true,format=:png);
global s=Settings(;reset=false,save=true,evalConstraints=false,MPC=true);

r.eval_num=1;count=1;
# set infeasible_counter for active safety
x.infeasible_counter=100; x.infeasible_counter_max=4;  # initially the condition will be false
println("Getting Vehicle Position From MATLAB");
while(Bool(x.runJulia))
  if count!=r.eval_num
    #println("Running model for the: ",r.eval_num," time");
  end
  count=r.eval_num;

  getPlantData!(n,params,x,c,r);

  updateX0!(n,r,x.X0;(:userUpdate=>true));

  if c.m.activeSafety
    checkFeasibility!(mdl,d,n,c,x,params[1],r;feas_tol=0.005);
    if !x.feasible
      x.infeasible_counter=1;
    else
      x.infeasible_counter+=1;
    end
  end

#  if !c.m.activeSafety || !x.feasible || (x.infeasible_counter < x.infeasible_counter_max)
  #println("Running Optimization");
  sharedControl!(mdl,n,r,s,params,x);   # rerun optimization
#  end
  sendOptData!(n,r,x)
end
close(x.s1)
close(x.s2)
mdl=0;n=0;r=0;params=0;c=0;x=0;d=0;pa=0;s=0;

#include("postProcess.jl")


text:
using Soss

p = @model begin
    x ~ Normal()
    y ~ Normal(x,1) |> iid(5)
end;

q = @model λ begin
    μ = λ.x.Ex
    σ = λ.x.Ex²-λ.x.Ex^2
    x ~ Normal(μ, σ)
end;

y = rand(p()).y;

λ = ( _ℓ = -20
    , x = (Ex=0.0, Ex²=1.0)
)

function logweightedmean(ℓ1,x1,ℓ2,x2)
    if ℓ1 < ℓ2
        return logweightedmean(ℓ2,x2,ℓ1,x1)
    elseif isinf(ℓ2)
        return ℓ1
    else
        w = exp(ℓ2 - ℓ1)
        return (x1 + w*x2)/(1+w)
    end
end

importanceSample(p(), q(λ=λ), (y=y,))

function f(λ, n=10)
    for j in 1:n
        ℓx = importanceSample(p(), q(λ=λ), (y=y,))
        val = ℓx.val
        ℓ = ℓx.ℓ

        x = val.x
        x² = x*x
        x³ = x*x²
        x⁴ = x*x³
        λ = @set λ.x.Ex = logweightedmean(λ._ℓ+1, λ.x.Ex, ℓ, x)
        λ = @set λ.x.Ex² = logweightedmean(λ._ℓ+1, λ.x.Ex², ℓ, x²)
        λ = @set λ.x.Ex³ = logweightedmean(λ._ℓ+1, λ.x.Ex³, ℓ, x³)
        λ = @set λ.x.Ex⁴ = logweightedmean(λ._ℓ+1, λ.x.Ex⁴, ℓ, x⁴)

        λ = @set λ._ℓ = logaddexp(λ._ℓ, ℓ)
    end
    return λ
end

λ = f(λ,1000)


text:
##
println("hello world")
##
function test(x)
    x + 3.0
end

##
test(123)

##
function greeks( α , β )
    α ^2 - 2 * β
end

greeks(12.2,44.0)
##

f(x) = 2*x + 3

println(f(4))
##
using Plots
x = 1:10; y = rand(10); # These are the plotting data
plot(x,y, label="my label")

##


text:
# Gaussian Processes in Julia
# based on http://mrmartin.net/?p=223

#####################################################################
# 1. Univariate Normal
#####################################################################

# in order to generate a univariate random number, use any real numbers mu
# and sigma to define the distribution.
mu = 6.2          # mean
sigma_squared = 2 # variance

# randn returns a random number from the standard normal distribution, which
# is a univariate normal with mu = 0, sigma = 1, i.e. N(0, 1):
standard_normal_random_number = randn()

# which we can then convert to arbitrary univariate gaussian normal, N(mu, sigma):
normal_random_number = mu + sqrt(sigma_squared) * standard_normal_random_number

# We can generate many of them:
normal_random_numbers = mu + sqrt(sigma_squared) * randn(100000)

# and plot what is going on:
using Gadfly
using DataFrames
df = DataFrame(nums = normal_random_numbers);
draw(PNG("hist_normal_rand_nums.png", 50cm, 36cm),
  plot(df, x = "nums", Geom.histogram(bincount=50)));

#####################################################################
# 2. Multivariate Normal
#####################################################################

# The same trick works for the multivariate normal. Mu has to be a vector,
# and Sigma a symetric semidefinite matrix
Mu = [3, 0]

# a matrix is symetric iff M(a,b) = M(b,a), a positive semidefinite iff
# x*M*x'>=0 for any vector x.
Sigma =  [ 1.0000   -0.9195;  -0.9195    1.0000]

# to find A such that A'*A = Sigma, we find the eigendecompositon of Sigma:
# V'*D*V = Sigma
D, V = eig(Sigma)
A = V .* (D .^ 0.5);

# But a quicker way if Sigma has many dimensions is to use the Cholesky
# decomposition instead:
# C = chol(Sigma)
# which then gives C' * C == Sigma

# The standard multivariate normal is just a series of independently drawn
# univariate normal random numbers
standard_random_vector = randn(2, 1)

# this is converted to an arbitrary multivariate normal by:
normal_random_vector = Mu .+ A * standard_random_vector

# so now we can generate many of them, and plot what's going on:
standard_random_vector = randn(2, 1000)
normal_random_vector = Mu .+ A * standard_random_vector

draw(PNG("scatter_std_gaussian.png", 50cm, 36cm),
  plot(x = standard_random_vector[1,:], y = standard_random_vector[2,:]));

draw(PNG("scatter_normal_gaussian.png", 50cm, 36cm),
  plot(x = normal_random_vector[1,:], y = normal_random_vector[2,:]));

#####################################################################
# 3. A Kernel Matrix demonstration Brownian motion
#####################################################################
n=40

# The brownian motion kernel is defined through its very simple inverse
inverse = 2 .* eye(n)
inverse = inverse - (triu(ones(n,n),1) - triu(ones(n,n),2))'
inverse = inverse - (triu(ones(n,n),1) - triu(ones(n,n),2))

# subplot(1,3,1);imagesc(inverse);title('Brownian motion K^{-1}')
# subplot(1,3,2);imagesc(inv(inverse));title('Brownian motion K')
#draw some samples
#r=randn(n,10)
#subplot(1,3,3);plot(inverse\r);title('10 samples')

#####################################################################
# 4. Gaussian Process
#####################################################################

# one way of looking at a Gaussian Process is as a Multivariate Normal with
# an infinite number of dimensions. However, in order to model
# relationships between points, we construct the covariance matrix with a
# function that defines the value of the matrix for any pair of real numbers:
sigma_f = 1.1251  # parameter of the squared exponential kernel
l =  0.90441      # parameter of the squared exponential kernel

# This is one of many popular kernel functions, the squared exponential
# kernel. It favors smooth functions.
kernel_function(x, x2) = sigma_f^2 * exp((x-x2)^2 / (-2*l^2))

# we can also define an error function, which models the observation noise
sigma_n = 0.1  # known noise on observed data
# this is just iid gaussian noise with mean 0 and variance sigma_n^2s
error_function(x, x2) = sigma_n^2 * (x==x2)

# kernel functions can be added together. Here, we add the error kernel to
# the squared exponential kernel:
k(x, x2) = kernel_function(x,x2) + error_function(x,x2)

# We can find the mean and the variance of the GP at each point
# mean = zeros()
n = int(abs(-2 - 1)/0.01)
prediction_x = linspace(-2, 1, n)

# now, we would like to sample from a Gaussian Process defined by this
# kernel. First, for the subset of points we are interested to plot, we
# construct the kernel matrix (using our kernel function)
function kernel_matrix_from_func(xs, kernel_func)
  n = length(xs)
  K = zeros(n, n)
  for(i in 1:n)
    for(j in i:n)
      K[i,j] = K[j,i] = k(xs[i], xs[j])
    end
  end
  return K
end

K = kernel_matrix_from_func(prediction_x, k)

# We will use K as our covariance matrix, so lets decompose it
function decompose(M)
  D, V = eig(K)
  V .* (D .^ 0.5);
end
A = decompose(K)

# and then we can sample:
gaussian_process_sample = A * randn(n, 7)

# Lets generalize this:
function gp_sample(kernelMatrix, numSamples)
  decompose(kernelMatrix) * randn(size(kernelMatrix, 1), numSamples)
end


text:
export ExtensionalSet, IntensionalSet, settuple, orderedpair, cardinality

abstract type MathematicalSet <: AbstractExpression end
istree(s::MathematicalSet) = false # TODO: refactor structure conditionals to make more sense

struct ExtensionalSet{T} <: MathematicalSet
    elements::Set{T}
end
ExtensionalSet(elements::Vector{T}) where {T} = ExtensionalSet(Set{T}(elements))
ExtensionalSet(elements::Tuple{Vararg{T}}) where {T} = ExtensionalSet(Set{T}([elements...]))
elements(es::ExtensionalSet) = es.elements
Base.length(es::ExtensionalSet) = length(elements(es))
Base.isempty(es::ExtensionalSet) = Base.isempty(elements(es))
Base.hash(es::ExtensionalSet, h::UInt) = hash(elements(es), h)
Base.:(==)(es1::ExtensionalSet, es2::ExtensionalSet) = elements(es1) == elements(es2)
cardinality(es::ExtensionalSet) = length(es)

# expression methods
variables(es::ExtensionalSet) = length(es) > 0 ? reduce(∪, variables.(elements(es))) : Set{LogicalSymbol}()
operations(es::ExtensionalSet) = length(es) > 0 ? reduce(∪, operations.(elements(es))) : Set{LogicalOperation}()

function Base.show(io::IO, es::ExtensionalSet)
    if isempty(es)
        print(io, "∅")
        return
    end

    print(io, "{")
    i = 1
    truncated_elements = truncate([elements(es)...], 5)
    for el ∈ truncated_elements
        print(io, el)
        if i < cardinality(es)
            print(io, ", ")
        end
        i += 1
    end
    if length(truncated_elements) < cardinality(es)
        missed_elements = cardinality(es) - length(truncated_elements)
        print(io, "...$(missed_elements) more")
    end
    print(io, "}")
end

orderedpair(a, b) = ExtensionalSet([ExtensionalSet([a]), ExtensionalSet([a, b])])
function settuple(a...)
    a = [a...]
    if length(a) == 0
        return ∅
    end
    x = settuple(a[1:end-1]...)
    ExtensionalSet([ExtensionalSet([x]), ExtensionalSet([x, a[end]])])
end

struct IntensionalSet <: MathematicalSet
    transform::AbstractExpression
    rule::AbstractExpression
end
IntensionalSet(transform::AbstractExpression, rules::Vector{LogicalExpression}) = IntensionalSet(transform, reduce(∧, rules))
IntensionalSet(transform::AbstractExpression, rules::Set{LogicalExpression}) = IntensionalSet(transform, [rules...])
IntensionalSet(transform::AbstractExpression, rules::Tuple{Vararg{LogicalExpression}}) = IntensionalSet(transform, [rules...])
transform(is::IntensionalSet) = is.transform
rule(is::IntensionalSet) = is.rule

# expression methods
variables(is::IntensionalSet) = variables(rule(is))
operations(is::IntensionalSet) = operations(rule(is))

function Base.show(io::IO, is::IntensionalSet)
    print(io, "{")
    print(io, transform(is))
    print(io, " | ")
    print(io, rule(is))
    print(io, "}")
end


text:
using Gen
using GenVariableElimination
import Random

function make_data_set(n)
    Random.seed!(1)
    prob_outlier = 0.5
    true_inlier_noise = 0.5
    true_outlier_noise = 5.0
    true_slope = -1
    true_intercept = 2
    xs = collect(range(-5, stop=5, length=n))
    ys = Float64[]
    for (i, x) in enumerate(xs)
        if rand() < prob_outlier
            y = true_slope * x + true_intercept + randn() * true_inlier_noise
        else
            y = true_slope * x + true_intercept + randn() * true_outlier_noise
        end
        push!(ys, y)
    end
    (xs, ys)
end

@gen (static) function datum(x, inlier_std, outlier_std, slope, intercept)
    is_outlier ~ bernoulli(0.5)
    std = is_outlier ? outlier_std : inlier_std
    y ~ normal(x * slope + intercept, std)
    return y
end

@gen (static) function model(xs::Vector{Float64})
    n = length(xs)
    inlier_log_std ~ normal(0, 2)
    outlier_log_std ~ normal(0, 2)
    inlier_std = exp(inlier_log_std)
    outlier_std = exp(outlier_log_std)
    slope ~ normal(0, 2)
    intercept ~ normal(0, 2)
    data ~ Map(datum)(
        xs,
        fill(inlier_std, n),
        fill(outlier_std, n),
        fill(slope, n),
        fill(intercept, n))
end

@gen (static) function slope_proposal(trace)
    slope ~ normal(trace[:slope], 0.5)
end

@gen (static) function intercept_proposal(trace)
    intercept ~ normal(trace[:intercept], 0.5)
end

@gen (static) function inlier_std_proposal(trace)
    inlier_log_std ~ normal(trace[:inlier_log_std], 0.5)
end

@gen (static) function outlier_std_proposal(trace)
    outlier_log_std ~ normal(trace[:outlier_log_std], 0.5)
end

@load_generated_functions()

using PyCall
graphviz = pyimport("graphviz")
function addr_to_name(addr)
    name = replace(replace(replace(replace(string(addr), " => " => "-"), "(" => ""), ")" => ""), ":" => "")
    println("addr: $addr, name: $name")
    return name
end

# for debugging..
function mh_log_acceptance_ratio(trace, proposal, proposal_args)
    forward_trace = simulate(proposal, (trace, proposal_args...))
    println("FORWARD TRACE")
    display(get_choices(forward_trace))
    println(get_score(forward_trace))
    (new_trace, weight, _, discard) = update(trace, get_args(trace), map((_) -> NoChange(), get_args(trace)), get_choices(forward_trace))
    (backward_trace, _) = generate(proposal, (new_trace, proposal_args...))
    return (new_trace, weight, get_score(backward_trace), get_score(forward_trace))
end

function do_inference(xs, ys, num_iters)
    observations = choicemap()
    for (i, y) in enumerate(ys)
        observations[:data => i => :y] = y
    end

    # initial trace
    display(observations)
    println(xs)
    (trace, _) = generate(model, (xs,), observations)

    addresses = [:data => i => :is_outlier for i in 1:length(xs)]

    # generate sampler for is_outlier indciator variables
    is_outlier_proposal = generate_backwards_sampler_fixed_structure(trace, addresses)

    scores = Vector{Float64}(undef, num_iters)
    for i=1:num_iters

        # steps on the parameters
        for j=1:5
            (trace, _) = metropolis_hastings(trace, slope_proposal, ())
            (trace, _) = metropolis_hastings(trace, intercept_proposal, ())
            (trace, _) = metropolis_hastings(trace, inlier_std_proposal, ())
            (trace, _) = metropolis_hastings(trace, outlier_std_proposal, ())
        end

        # block gibbs step on the is_outlier variables (using variable elimination)
        @time trace, acc = metropolis_hastings(trace, is_outlier_proposal, ())
        println("acc: $acc") # should accept always, except in the beginning where near-deterministic proposals arise

        score = get_score(trace)
        scores[i] = score

        # print
        slope = trace[:slope]
        intercept = trace[:intercept]
        inlier_std = exp(trace[:inlier_log_std])
        outlier_std = exp(trace[:outlier_log_std])
        println("score: $score, slope: $slope, intercept: $intercept, inlier_std: $inlier_std, outlier_std: $outlier_std")
    end
    return scores
end

(xs, ys) = make_data_set(30)
do_inference(xs, ys, 10)
@time do_inference(xs, ys, 50)


text:
using Compat.Base64

function create_event()
    data_base64 = "[encoded data, 10300 characters]"

    headers = Dict(
        "content-type"      => "application/json",
        "http_minor"        => "1",
        "Keep-Alive"        => "1",
        "User-Agent"        => "GitHub-Hookshot/333881f",
        "X-Github-Delivery" => "42e04680-7d9d-11e5-8e29-6553a49e3e25",
        "X-Forwarded-For"   => "192.30.252.34",
        "Content-Length"    => "7724",
        "Host"              => "2a204a87.ngrok.io",
        "X-Github-Event"    => "commit_comment",
        "http_major"        => "1",
        "Accept"            => "*/*",
        "X-Hub-Signature"   => "sha1=494685b443b34bb0241a14f82604a0600b18f996",
    )

    return HTTP.Request("GET", "http://www.github.com", headers, base64decode(data_base64))
end


text:
module TestTMapreduce
using Test
using KissThreading

struct FreeAbelianSemigroup{T}
    coords::Dict{T, Int}
end

const FASG{T} = FreeAbelianSemigroup{T}

function Base.:*(x::FASG, y::FASG)
    FASG(merge(+, x.coords, y.coords))
end
Base.:(==)(x::FASG, y::FASG) = x.coords == y.coords
pure(x) = FASG(Dict(x=>1))

@testset "tmapreduce" begin
    for setup in [
        (f=identity, op=+, src=1:10,            init=0),
        (f=identity, op=+, src=1:10,            init=21),
        (f=identity, op=+, src=Int[],           init=0),
        (f=identity, op=+, src=Int[],           init=21),
        (f=x->2x,    op=*, src=collect(1:10),   init=21),
        (f=x->2x,    op=*, src=rand(Int, 10^5), init=rand(Int)),
        (f=pure,     op=*, src=randn(4),        init=pure(42.0)),
        (f=pure,     op=*, src=randn(10^4),     init=pure(42.0)),
       ]
        res_base = @inferred mapreduce(setup.f, setup.op, setup.src, init=setup.init)
        res_kiss = @inferred tmapreduce(setup.f, setup.op, setup.src, init=setup.init)
        @test res_base == res_kiss
    end
    @test tmapreduce(pure, *, 1:1, init=pure(-1)) == pure(-1) * pure(1)
    @test tmapreduce(pure, *, 1:1, init=pure(-1)) == mapreduce(pure, *, 1:1, init=pure(-1))

    # multiple src arrays
    @test 1:6 == @inferred tmapreduce(vcat, vcat, [1,3,5], [2,4,6], init=[])
    for setup in [
            (f=+, op=*, src=(1:10, 2:2:20), init=rand(Int)),
            (f=*, op=+, src=[rand(Int, 100) for _ in 1:5], init=rand(-10:10)),
            (f=pure∘tuple, op=*, src=[rand(-1:2, 100) for _ in 1:5], init=pure((1,2,3,4,5)))
       ]
        res_base = @inferred reduce(setup.op, map(setup.f, setup.src...), init=setup.init)
        res_kiss = @inferred tmapreduce(setup.f, setup.op, setup.src..., init=setup.init)
        @test res_base == res_kiss
    end

    # incompatible input size
    tmapreduce(+,*,randn(10), randn(10), init=0.0)
    @test_throws ArgumentError tmapreduce(+,*,randn(10), randn(11), init=0.0)
end

end#module


text:
# functions that wrap around other pieces to do higher level tasks, e.g.,
#     perform cross-validation of prediction scores on a datasets
#     optimisation of a parameter
#

"""
    test_cv_1()

 Calculates k-fold cross-validated scores for predictions based on an update rule and ratings.

## Arguments
* `rule::UpdateRule`: the update rule to test
* `input_competitions::DataFrame`: a DataFrame in format specified in IO.jl
* `k::Int = 10`: doing k-fold cross-validation
* `n_samples::Int = 1000`: number of re-samples of the underlying dataset to use (see SampleIterate.jl)
* `batch_size::Int = 30`: size of batches in re-sampled iteration (see SampleIterate.jl)
* `srule::ScoringRule = ScoreBrier(;normalise=true)`: scoring schema to use (see ScoringRules/)

## Outputs
* S = length k vector of the cross-validation scores

```
"""
function cross_validate(rule::UpdateRule, input_competitions::DataFrame;
                        k::Int = 10, # cross validation parameters
                        n_samples::Int = 1000, # iteration parameters
                        batch_size::Int = 30, # iteration parameters
                        srule::ScoringRule = ScoreBrier(;normalise=true) )
    player_list = merge(+, countmap( input_competitions[PlayerA] ), countmap( input_competitions[PlayerB] )  )
    input_ratings = RatingsList( player_list ) # a blank set of starter ratings
    sampledrule = UpdateSampleIterate(; rule=rule, batch_size = batch_size, n_samples = n_samples)

    # k-fold cross-validation
    # separate data into training and test data, at random
    indexes = Random.shuffle( 1:size(input_competitions,1) )
    n = Int( k*floor(length(indexes)/k) )
    n_block = Int(n/k)
    indexes = indexes[1:n] # randomly drop a small amount of left-over data :(

    S = zeros(Float64, k)
    for i=1:k
        test_indexes = indexes[ ((i-1)*n_block)+1 : i*n_block ]
        training_indexes = setdiff( indexes, test_indexes )

        # "train" by calculating ratings from training data
        training_data = input_competitions[ training_indexes, : ]
        out_ratings = update_ratings(rule, input_ratings, training_data)

        # score on the test data
        test_data = input_competitions[ test_indexes, : ]
        S[i] = score_ratings( srule, test_data, rule, out_ratings )
    end

    return S
end

function cross_validate(rule::UpdateRule, file::String;
                        k::Int = 10,
                        n_samples::Int = 1000,
                        batch_size::Int = 30,
                        srule::ScoringRule = ScoreBrier(;normalise=true) )
    (input_competitions,  player_list) = read_results( file )

    return cross_validate(rule, input_competitions; k=k, n_samples=n_samples, batch_size=batch_size, srule=srule)
end

function optimise(rule::UpdateRule, input_competitions::DataFrame;
                  k::Int = 10, # cross-validation parameters
                  n_samples::Int = 1000, # iteration parameters
                  batch_size::Int = 30, # iteration parameters
                  srule::ScoringRule = ScoreBrier(;normalise=true),
                  optimize_this # parameter that we need to optimize
                  )
    # really should be a separate cross-validation process for choosing "hyperparameters"
    #   need data partitioned into training, test, and validation datasets

    # line search does the search, but its kind of ugly to call
    #   need a nicer strategy to set up the calls

    # or have a separate function to optimise each parameter?, but then how many and how specific?
end


text:
# Welcome to your new Toolips server!
using Main.Toolips

PUBLIC = "../public"
IP = "127.0.0.1"
PORT = 8000
function main()
        # Essentials
    routes = make_routes()
    server_template = ServerTemplate(IP, PORT, routes)
        # Fun stuff (examples !, you should probably delete these.)
    delayed = Route("/delay", fn(delay))
    suicide = Route("/suicide", fn(suicide_fn))
    arguments = Route("/args", fn(args))
    server_template.add(delayed)
    server_template.add(suicide)
    server_template.add(arguments)
    global TLSERVER = server_template.start()
    return(TLSERVER)
end

# Routes
function make_routes()
        # Pages
    four04 = html("
404, Page not found!
")
    index = html("
Hello world!

Not so exciting, is it?
     well, it is a work in progress :p.
")
        # Routes
    routes = []
    homeroute = Route("/", index)
    four04route = Route("404", four04)
    push!(routes, homeroute)
    push!(routes, four04route)
    routes
end

suicide_fn = http -> stop!(TLSERVER)

args = http -> string(getargs(http))

function delay(http::Any)
        for character in "Hello World!"
            write(http, string(character))
            sleep(1)
        end
end

main()


text:
using DirectConvolution
using Documenter

DocMeta.setdocmeta!(DirectConvolution, :DocTestSetup, :(using DirectConvolution); recursive=true)
doctest(DirectConvolution, fix=true)

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
# File for plotting Figure 4

using JLD2
using DataFrames
using Gadfly

import Cairo, Fontconfig

tmp_dict = load("data/tree_depth.jld2")
df_tmp = tmp_dict["df_tmp"]

df_plot = stack(df_tmp)
rename!(df_plot, "variable" => "Approach")

df_plot.Approach .= replace.(df_plot.Approach, "unweighted" => "CART")

p = plot(Scale.color_discrete_manual("#264653", "#E7C15F"), Scale.x_discrete,
Theme(bar_spacing=2mm))

push!(p, layer(df_plot, x=:n, y=:value, color=:Approach, Geom.bar(position=:dodge)))
push!(p, Guide.YLabel("Tree Depth"))
push!(p, Guide.XLabel("Observations n"))

draw(PNG("figures/graphs/tree_depth.png", 20cm, 12cm, dpi=300), p)

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
# All square roots are periodic when written as continued fractions and can be
# written in the form:
#
# √N = a0 + 1 / (a1 + 1 / (a2 + 1 / (a3 + ...)))
#
# For example, let us consider √23:
# √23 = 4 + √23 − 4 = 4 + 1 / (1 / (√23 − 4)) = 4 + 1 / (1 + (√23 − 3) / 7)
#
# If we continue we would get the following expansion:
# √23 = 4 + 1 / (1 + 1 / (3 + 1 / (1 + 1 / (8 + ...))))
#
# The process can be summarised as follows:
#
# a0 = 4, 1 / (23√ - 4) = (√23 + 4) / 7 = 1 + (√23 - 3) / 7
# a1 = 1, 7 / (23√ - 3) = 7 * (23√+3)14 = 3 + (√23 - 3) / 2
# a2 = 3, 2 / (23√ - 3) = 2 * (23√+3)14 = 1 + (√23 - 4) / 7
# a3 = 1, 7 / (23√ - 4) = 7 * (23√+4)7 = 8 + √23 - 4
# a4 = 8, 1 / (23√ - 4) = (√23 + 4) / 7 = 1 + (√23 - 3) / 7
# a5 = 1, 7 / (23√ - 3) = 7 * (√23 + 3) / 14 = 3 + (√23 - 3) / 2
# a6 = 3, 2 / (23√ - 3) = 2 * (√23 + 3) / 14 = 1 + (√23 - 4) / 7
# a7 = 1, 7 / (23√ - 4) = 7 * (√23 + 4) / 7 = 8 + √23 - 4
#
# It can be seen that the sequence is repeating. For conciseness, we use the
# notation √23 = [4;(1,3,1,8)], to indicate that the block (1,3,1,8) repeats
# indefinitely.
#
# The first ten continued fraction representations of (irrational) square roots
# are:
#
# √2 = [1;(2)], period=1
# √3 = [1;(1,2)], period=2
# √5 = [2;(4)], period=1
# √6 = [2;(2,4)], period=2
# √7 = [2;(1,1,1,4)], period=4
# √8 = [2;(1,4)], period=2
# √10 = [3;(6)], period=1
# √11 = [3;(3,6)], period=2
# √12 = [3;(2,6)], period=2
# √13 = [3;(1,1,1,1,6)], period=5
#
# Exactly four continued fractions, for N ≤ 13, have an odd period.
#
# How many continued fractions for N ≤ 10000 have an odd period?

using ProjectEulerSolutions

# A continued fraction period ends when the a,m,d tuple matches a previously
# seen one.  Just maintain the list of tuples until one matches, measuring the
# length of the continued fraction period to find the odd numbered ones.
function p064solution(n_max::Integer=10)::Integer

    odd_count = 0
    for n in 2:n_max
        a0 = floor(sqrt(n))
        if sqrt(n) == a0
            continue
        end
        a0 = Integer(a0)
        a, m, d = a0, 0, 1
        s = Array{Tuple{Integer, Integer, Integer}, 1}()
        while true
            m = d * a - m
            d = div(n - m * m, d)
            a = div(a0 + m, d)
            key = (a, m, d)
            if length(s) > 0 && key == s[1]
                break
            else
                push!(s, (a, m, d))
            end
        end
        if length(s) % 2 == 1
            odd_count += 1
        end
    end
    return odd_count
end

p064 = Problems.Problem(p064solution)

Problems.benchmark(p064, 10_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# eql_root.jl

export eql_root

using Test: @test, @test_throws, @inferred

"""
    x = eql_root(a,b,c)
    Numerically stable method for computing the positive root
    of the quadratic polynomial `-ax^2 - 2bx + c, a >= 0`.
    Assumes solvable equations; will throw otherwise.

in
- `a` : The negative of the `x^2` term. Must be positive.
- `b` : Half the negative of the `x` term.
- `c` : The constant term.

out
- `x` : The positive root that satisfies `0 = -ax^2 - 2bx + c`.
"""
function eql_root(a::Real, b::Real, c::Real)
    (a < 0) && throw(DomainError(a, "a must be enonnegative"))

    if a == 0
        return c / b / 2 # trivially solve case where a == 0
    end

    det = sqrt(b^2 + a * c) # determinant / 2

    # handle positive and negative b appropriately:
    return (b > 0) ? c / (det + b) : (det - b) / a
end

function eql_root(a, b, c)
    (size(a) != size(b) || size(b) != size(c)) && throw(DimensionMismatch("all arguments must share dimensions"))

    eql_root.(a, b, c)
end

"""
    eql_root(:test)
self test
"""
function eql_root(x::Symbol)
    x != :test && throw("Must use :test as input for testing")
    tests = [1 1 -1;
             0 1 2;
             4 5 0;
             1 0 9; # 3 both times
             #1 0 -9; # -3 both times. omitted because out of prog scope
             4 4 5] # (2x+5)(2x-1), so 1/2.
    predicted = [-1, 1, 0, 3, 1/2]
    results = eql_root(tests[:,1], tests[:,2], tests[:,3])

    @test results == predicted
    @test eql_root(1, 1, -1) == -1
    return true
end


text:
### A Pluto.jl notebook ###
# v0.17.0

using Markdown
using InteractiveUtils

# ╔═╡ 39da52d0-3e5b-4dad-9de9-6a8b0337f124
begin
	# Run once if dependencies not installed
	using Pkg
	Pkg.add("Plots")
end

# ╔═╡ 4345a821-21fa-4a84-91a9-c3c539db2725
begin
	using Plots
end

# ╔═╡ 88843947-6008-41cc-b204-3f822e542ae3
md"""
# Chapter 1
## Section 2
"""

# ╔═╡ 98cdb6dc-3a40-11ec-1eb2-f7e085a3b6e6
function trapezoid(f, n::Integer, a::T, b::T)::T where T <: Real
	dx = (b - a) / n
	if a < b
		a, b = a + dx, b - dx
		sign = 1.
	else
		a, b = b + dx, a - dx
		sign = -1
	end
	integral = 0
	h = 2. * dx
	for x in a:h:b
		integral += 0.5 * dx * (f(x + dx) + 2 * f(x) + f(x - dx))
	end
	integral
end

# ╔═╡ e62b3b53-60c8-4450-978c-7558398fa085
trapezoid(sin, 256, 0., float(pi))

# ╔═╡ 7aa7729a-9c6c-457b-a271-8cd0775a65bf
function simpson(f, n::Integer, a::T, b::T)::T where T <: Real
	dx = (b - a) / n
	if a < b
		a, b = a + dx, b - dx
		sign = 1.
	else
		a, b = b + dx, a - dx
		sign = -1
	end
	integral = 0
	h = 2. * dx
	for x in a:h:b
		integral +=  (dx / 3) * (f(x + dx) + 4 * f(x) + f(x - dx))
	end
	integral
end

# ╔═╡ c1449e54-2bfa-4fce-a7b4-227451327b48
function bode(f, n::Integer, a::T, b::T)::T where T <: Real
	dx = (b - a) / n
	h = 4. * dx
	if a < b
		a, b = a + h / 2., b - h / 2.
		sign = 1.
	else
		a, b = b + h / 2., a - h / 2.
		sign = -1
	end
	integral = 0
	for x in a:h:b
		integral +=  (2 * dx / 45) * (
			7 * f(x + 2. * dx) + 32. * f(x + dx)
			+  12. * f(x) +
			7 * f(x - 2. * dx) + 32. * f(x - dx))
	end
	integral
end

# ╔═╡ 1dea0b53-fc17-43f3-b1df-495fff0529e9
begin
	# Exercice 1.2
	x = collect(2:8) # Weird result betwin 7 ~ 14
	func = [trapezoid, simpson, bode]
	y = Array{Float64}(undef, (length(func), length(x)))
	for (k, f) in enumerate(func)
		for i in 1:length(x)
			y[k, i] = abs(2. - f(sin, 1 << x[i], 0., float(pi)))
		end
	end
	plot(x, y', label=["trap" "simp" "bode"])
end

# ╔═╡ 76f93ada-aa5a-4b6b-abed-f006dbe6f89a
# Exercice 1.3
md"""

The function: $t^{-2/3} (1 - t) ^{-1/3)}$ have 2 singularity 0 and 1

We will need ton transform is integrale to be able to compute it

First we spliting the integral into two part

$$\int_0^1 t^{-2/3} (1 - t) ^{-1/3)} = \int_0^{0.5} t^{-2/3} (1 - t) ^{-1/3dt + \int_{0.5}^1 t^{-2/3} (1 - t) ^{-1/3}dt$$

Now we have isolated the singularities. Now let change variables

"""

# ╔═╡ 46f5c4bb-8f93-4c32-81ea-717a05232aec
begin

end

# ╔═╡ Cell order:
# ╟─88843947-6008-41cc-b204-3f822e542ae3
# ╠═39da52d0-3e5b-4dad-9de9-6a8b0337f124
# ╠═4345a821-21fa-4a84-91a9-c3c539db2725
# ╠═98cdb6dc-3a40-11ec-1eb2-f7e085a3b6e6
# ╠═e62b3b53-60c8-4450-978c-7558398fa085
# ╠═7aa7729a-9c6c-457b-a271-8cd0775a65bf
# ╠═c1449e54-2bfa-4fce-a7b4-227451327b48
# ╠═1dea0b53-fc17-43f3-b1df-495fff0529e9
# ╠═76f93ada-aa5a-4b6b-abed-f006dbe6f89a
# ╠═46f5c4bb-8f93-4c32-81ea-717a05232aec

max_stars_repo_name: zazbone/ComputationalPhysics

max_issues_repo_name: zazbone/ComputationalPhysics

max_forks_repo_name: zazbone/ComputationalPhysics


text:
"Internal! this should hold all PyTest specific args gathered from the invocation"
runner_args = Dict{AbstractString, Any}()

"Helper function that injects invocation args to all @pytest calls"
function set_parsed_args!(pa)
 global runner_args
 runner_args = pa
end