text:
using Test
using WavePropBase
using StaticArrays
using WavePropBase.Geometry
using WavePropBase.Integration
using LinearAlgebra
using Random
Random.seed!(1)

@testset "Reference segment" begin
    qstd = Fejer(20)
    d     = domain(qstd)
    # simple test on smooth integrand
    for shand in [IMT(), Kress(), KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        @test isapprox(integrate(cos,q),sin(1),rtol=1e-2)
    end
    # non-smooth integrand, singularity at 0
    f     = (x) -> log(abs(x))
    Ie    = -1
    for shand in [IMT(), Kress(),KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        Ia    = integrate(f,q)
        @test isapprox(Ia,Ie,rtol=1e-4)
        # check that the `naive` integration woudl have failed the test
        Istd    = integrate(f,qstd)
        @test !isapprox(Istd,Ie,rtol=1e-4)
    end
    # test that KressP() can handle singularity at both endpoints, and Kress() cannot
    f     = (x) -> log(x) + log(1-x)
    Ie    = -2
    q     = SingularQuadratureRule(qstd,KressP())
    Ia    = integrate(f,q)
    @test isapprox(Ia,Ie,rtol=1e-5)
    q     = SingularQuadratureRule(qstd,Kress())
    Ia    = integrate(f,q)
    @test !isapprox(Ia,Ie,rtol=1e-5)
    # non-smooth integrand with singularity inside
    s     = 0.1
    f     = (x) -> log(abs(x-s))
    Ie    = -1.32508297337145
    # FIXME: for the time being, usign IMT with a singular quadrature rule for
    # interior points yields a rounding problem because the change of variables
    # accumulates points very quickly around zero. That is OK when the
    # singularity is at zero, but shifting such nodes to another point s ∈ (0,1]
    # has to be done with care since these can get rounded to `s`, which in turn
    # yields a singular value for the integrand. Note that the same can happen
    # with other singularity handler types, but IMT is worse due to the
    # exponential nature of the change of variables.
    for shand in [Kress(),KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        x,w  = q(s) # singular nodes and weights which depend on s
        Ia    = integrate(f,x,w)
        @test isapprox(Ia,Ie,rtol=1e-5)
        # check that the `naive` integration woudl have failed the test
        Istd    = integrate(f,qstd)
        @test !isapprox(Istd,Ie,rtol=1e-5)
    end
end

@testset "Duffy" begin
    q1d   = Fejer(5)
    qstd  = TensorProductQuadrature(q1d,q1d)
    duffy = Duffy()
    qsin  = SingularQuadratureRule(qstd,duffy)
    x,w  = qsin()
    # regular kernel
    k    = x -> cos(x[1])*sin(x[2])
    Ie   = 1/2*(sin(1)-cos(1))
    Ia   = integrate(k,qsin)
    @test isapprox(Ie,Ia,rtol=1e-4)
    # singular kernel at right vertex
    s    = SVector(1,0)
    k    = x -> 1/norm(x-s)
    Ia   = integrate(k,qsin)
    Ie   = acosh(sqrt(2))
    @test isapprox(Ie,Ia,rtol=1e-5)
    # singular kernel at wrong vertex (should not be accurate, see test below)
    s    = SVector(0,1)
    k    = x -> 1/norm(x-s)
    Ie   = acosh(sqrt(2))
    Ia   = integrate(k,qsin)
    @test !isapprox(Ie,Ia,rtol=1e-5)
end

@testset "2d Kress" begin
    q1d   = Fejer{10}()
    qstd  = TensorProductQuadrature(q1d,q1d)
    s1d   = Kress(order=2)
    sing_handler = TensorProductSingularityHandler(s1d,s1d)
    qsin = SingularQuadratureRule(qstd,sing_handler)
    x,w  = qsin()
    # regular kernel
    k = (x) -> cos(x[1])
    Ie   = sin(1)
    Isin = sum(k.(x).*w)
    @test isapprox(Ie,Isin,rtol=1e-6)
    # singular kernel at origin
    k  = (x) -> 1/norm(x)
    Ie = 1/2*log(17 + 12*sqrt(2)) # Wolfram alpha is your friend ;-)
    Isin = sum(k.(x).*w)
    @test isapprox(Ie,Isin,rtol=1e-2)
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
#= ==========================================================================================
n-th dimensional interpolation
========================================================================================== =#

# Method for a variable number of inputs, which are the variables themselves
# it returns a vector of constants

function interpolate(z, variables::Array{Float64,1}...)
    m = Array{Int64}(undef, length(variables));
    for i in 1:length(variables)
        m[i] = length(unique(variables[i]));
    end
    Vt = 1;
    for i in 1:length(variables)
        x = unique(variables[i]);
        V = Array{Float64}(undef, m[i], m[i]);
        for j in 0:m[i] - 1
            for k in 1:m[i]
                V[k, j + 1] = x[k]^j; # Si truena fue esto
            end
        end
        Vt = kron(Vt, inv(V))
    end
    return Vt * z, m
end

# Method that takes a data array as an input
# it returns a vector of constants

function interpotale(data::Array{Float64,2})
    m = Array{Int64}(undef, size(data, 2) - 1);
    for i in 1:size(data, 2) - 1
        m[i] = length(unique(data[:,i]));
    end
    Vt = 1;
    for i in 1:size(data, 2) - 1
        x = unique(data[:,i]);
        V = Array{Float64}(undef, m[i], m[i]);
        for j in 0:m[i] - 1
            for k in 1:m[i]
                V[k, j + 1] = x[k]^j; # Si truena fue esto
            end
        end
        Vt = kron(Vt, inv(V))
    end
    return Vt * data[:,end], m
end

# Function that writes the polynomial in LaTeX
function pol2tex(c::Array{Float64,1}, m::Array{Int64,1})
    M = 1;
    for i in 1:length(m)
        M *= m[i];
    end
    c = round.(c, sigdigits=5);
    for i in 1:length(c)
        if abs(c[i]) < 1e-10
            c[i] = 0;
        end
    end

    function foranid(n::Int64, i::Array{Int64,1})
        if n == length(m) + 1
            if c[k] != 0
                if abs(c[k]) != 1
                    model[k] = string(abs(c[k]));
                else
                    model[k] = "";
                end
                for j in 1:n - 1
                    if i[j] != 0
                        if i[j] == 1
                            model[k] *= "x_" * "{" * string(j) * "}";
                        else
                            model[k] *= "x_" * "{" * string(j) * "}^{" * string(i[j]) * "}";
                        end
                    end
                end
            end
            k += 1;
        else
            for q in 0:m[n] - 1
                model = foranid(n + 1, [i;q]);
            end
        end
        return model
    end
    k = 1;
    model = Array{String}(undef, M, 1);
    for q in 0:m[1] - 1
        model = foranid(2, [q]);
    end

    first = false;
    if c[1] != 0
        modelo = string(c[1]);
        first = true;
    end

    for i in 2:M
        if first == false
            if c[i] > 0
                modelo = model[i];
            elseif c[i] <  0
                modelo = " -" + model[i];
            end
            first = true;
        else
            if c[i] > 0
                modelo *= " + " * model[i];
            elseif c[i] <  0
                modelo *= " - " * model[i];
            end
        end
    end
    for i = 10:-1:1
        modelo = replace(modelo, "e-" * string(i) => "\times 10^{-" * string(i) * "}");
    end
    return modelo
end

# Function that substitutes the variable names of the LaTeX string into sth desired by the user
function sustex(modelo::String, vars::Array{String,1})
    for i = length(vars):-1:1
        modelo = replace(modelo, "x_{" * string(i) * "}" => vars[i]);
    end
    return modelo
end

#= # Función del polinomio
function polin(c::Array{Float64,1}, m::Array{Int64,1}, x::Array{Float64,1})

    function foranid(n::Int64, i::Array{Int64,1})
        if n == length(m) + 1
            term = c[k];
            for j in 1:n - 1
                term *= x[j]^i[j]
            end
            z += term;
            k += 1;
        else
            for q in 0:m[n] - 1
                z = foranid(n + 1, [i;q]);
            end
        end
        return z
    end
    k, z = 1, 0; # initialize
    for q in 0:m[1] - 1
        z = foranid(2, [q]);
    end
    return z;
end =#

ext: jl


text:
using Turing
using StatsFuns

using Random
using Test

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@turing_testset "is.jl" begin
    function reference(n :: Int)
        logweights = zeros(Float64, n)
        samples = Array{Dict{Symbol,Any}}(undef, n)
        for i = 1:n
            samples[i] = reference()
            logweights[i] = samples[i][:logweight]
        end
        logevidence = logsumexp(logweights) - log(n)
        results = Dict{Symbol,Any}()
        results[:lp] = logevidence
        results[:lp] = logweights
        results[:samples] = samples
        return results
    end

    function reference()
        x = rand(Normal(4,5))
        y = rand(Normal(x,1))
        log_lik = logpdf(Normal(x,2), 3) + logpdf(Normal(y,2), 1.5)
        d = Dict()
        d[:logweight] = log_lik
        d[:a] = x
        d[:b] = y
        return d
    end

    let n = 10
      @model normal() = begin
        a ~ Normal(4,5)
        3 ~ Normal(a,2)
        b ~ Normal(a,1)
        1.5 ~ Normal(b,2)
        a, b
      end

      alg = IS()
      seed = 0

      _f = normal();
      for i=1:100
        Random.seed!(seed)
        # Move the rng twice - equivalent to sampling from prior
        rand(Normal(4,5))
        rand(Normal(0,1))
        exact = reference(n)
        Random.seed!(seed)
        tested = sample(_f, alg, n)
        t_vals = get(tested, [:a, :b, :lp])
        for i = 1:n
            @test exact[:samples][i][:a] == t_vals.a[i]
            @test exact[:samples][i][:b] == t_vals.b[i]
            @test exact[:lp][i] == t_vals.lp[i]
        end
      end
    end

    @turing_testset "logevidence" begin
        Random.seed!(100)

        @model function test()
            a ~ Normal(0, 1)
            x ~ Bernoulli(1)
            b ~ Gamma(2, 3)
            1 ~ Bernoulli(x / 2)
            c ~ Beta()
            0 ~ Bernoulli(x / 2)
            x
        end

        chains = sample(test(), IS(), 10000)

        @test all(isone, chains[:x])
        @test chains.logevidence ≈ - 2 * log(2)
    end
end


text:
include("../coin-or-common.jl")

name = "Cbc"
version = Cbc_version

# Collection of sources required to build CbcBuilder
sources = [
    GitSource("https://github.com/coin-or/Cbc.git", Cbc_gitsha),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Cbc*

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la
update_configure_scripts

# old and custom autoconf
sed -i s/elf64ppc/elf64lppc/ configure

mkdir build
cd build/

export CPPFLAGS="${CPPFLAGS} -DNDEBUG -I${includedir} -I${includedir}/coin"
if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L$prefix/bin"
elif [[ ${target} == *linux* ]]; then
    export LDFLAGS="-ldl -lrt"
fi

../configure \
    --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --with-pic \
    --disable-pkg-config \
    --disable-debug \
    --enable-shared \
    lt_cv_deplibs_check_method=pass_all \
    --with-asl-lib="-lasl" \
    --with-blas-lib="-lopenblas" \
    --with-lapack-lib="-lopenblas" \
    --with-clp-lib="-lClp" \
    --with-coindepend-lib="-lCgl -lOsiClp -lClp -lOsi -lCoinUtils" \
    --enable-cbc-parallel

make -j${nproc}
make install
"""

# The products that we will ensure are always built
products = [
    LibraryProduct("libCbc", :libCbc),
    LibraryProduct("libCbcSolver", :libcbcsolver),
    ExecutableProduct("cbc", :cbc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("ASL_jll", ASL_version),
    Dependency("Cgl_jll", Cgl_version),
    Dependency("Clp_jll", Clp_version),
    Dependency("Osi_jll", Osi_version),
    Dependency("CoinUtils_jll", CoinUtils_version),
    Dependency("OpenBLAS32_jll", OpenBLAS32_version),
    Dependency("CompilerSupportLibraries_jll")
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The minimum distortion problem

# This example arises from computational geometry, in particular the problem of
# embedding a general finite metric space into a euclidean space.
#
# It is known that the 4-point metric space defined by the star graph:
#
#     x
#      \\
#       x — x
#      /
#     x
#
# where distances are computed by length of the shortest path between vertices,
# cannot be exactly embedded into a euclidean space of any dimension.
#
# Here we will formulate and solve an SDP to compute the best possible embedding,
# that is, the embedding f() that minimizes the distortion c such that
#
#     (1 / c) * D(a, b) ≤ ||f(a) - f(b)|| ≤ D(a, b)
#
# for all points (a, b), where D(a, b) is the distance in the metric space.
#
# Any embedding can be characterized by its Gram matrix Q, which is PSD, and
#
#     ||f(a) - f(b)||^2 = Q[a, a] + Q[b, b] - 2 * Q[a, b]
#
# We can therefore constrain
#
#     D[i, j]^2 ≤ Q[i, i] + Q[j, j] - 2 * Q[i, j] ≤ c^2 * D[i, j]^2
#
# and minimize c^2, which gives us the SDP formulation below.
#
# For more detail, see "Lectures on discrete geometry" by J. Matoušek, Springer,
# 2002, pp. 378-379.

using JuMP
import SCS
import Test

function example_min_distortion()
    model = Model(SCS.Optimizer)
    set_silent(model)
    D = [
        0.0 1.0 1.0 1.0;
        1.0 0.0 2.0 2.0;
        1.0 2.0 0.0 2.0;
        1.0 2.0 2.0 0.0
    ]
    @variable(model, c² >= 1.0)
    @variable(model, Q[1:4, 1:4], PSD)
    for i in 1:4
        for j in (i + 1):4
            @constraint(model, D[i, j]^2 <= Q[i, i] + Q[j, j] - 2 * Q[i, j])
            @constraint(model, Q[i, i] + Q[j, j] - 2 * Q[i, j] <= c² * D[i, j]^2)
        end
    end
    @objective(model, Min, c²)
    optimize!(model)
    Test.@test termination_status(model) == MOI.OPTIMAL
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT
    Test.@test objective_value(model) ≈ 4/3 atol = 1e-4
    return
end

example_min_distortion()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

import Base.Threads: @sync, @spawn, nthreads, threadid
using Random
using Test

const nx = 64
const np = 1000000

# +
function serial_deposition( steps )

    xmin, xmax = -6, 6
    Lx = xmax - xmin
    rho = zeros(Float64, nx)
    rng = MersenneTwister(123)
    xp = randn(rng, np);
    wp = similar(xp)
    fill!(wp, 1.0 / np)
    ntid = nthreads()

    for step = 1:steps

        fill!(rho, 0.0)

        for i in 1:np
            x_norm = (xp[i]-xmin) / Lx
            ip = trunc(Int,  x_norm * nx)+1
            rho[ip] += wp[i]
        end

    end

    rho

end

# +
function parallel_deposition( steps )

    xmin, xmax = -6, 6
    Lx = xmax - xmin
    rho = zeros(Float64, nx)
    rng = MersenneTwister(123)
    xp = randn(rng, np);
    wp = similar(xp)
    fill!(wp, 1.0 / np)
    ntid = nthreads()
    rho_local = [zero(rho) for _ in 1:ntid]
    chunks = Iterators.partition(1:np, np÷ntid)

    for step in 1:steps

        @sync for chunk in chunks
            @spawn begin
                tid = threadid()
                fill!(rho_local[tid], 0.0)
                for i in chunk
                    x_norm = (xp[i]-xmin) / Lx
                    ip = trunc(Int,  x_norm * nx)+1
                    rho_local[tid][ip] += wp[i]
                end
            end
        end

        rho .= reduce(+,rho_local)

     end

     rho

end
# -

@show nthreads()

serial_deposition(1) # trigger compilation
@time rho = serial_deposition(100)
@test sum(rho) ≈ 1.0

parallel_deposition(1) # trigger compilation
@time rho = parallel_deposition(100)
@test sum(rho) ≈ 1.0


text:
function loadData(path::String)
    datasets = DataFrame[] # Init empty vector
    if ispath(path)
        if isfile(path)
            push!(datasets, CSV.read(path))
            return datasets
        else
            files = readdir(path) # grab all files in the given folder
            for f in files
                push!(datasets, CSV.read(path * "/" * f))
            end
            return datasets
        end
    else
        error("Enter a valid file or folder")
    end
end

function arData(datasets::Vector{DataFrame})
    return Array.(datasets)
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
using Test, BenchmarkTools
include("../src/SpeedyWeather.jl")
using .SpeedyWeather

# original speedy T30 with 96x48 grid
P = Parameters(NF=Float64;trunc=30,nlat=48,nlon=96)
G = GeoSpectral(P)
B = Boundaries(P,G)

geopot_surf_spectral = B.geopot_surf
geopot_surf_grid = gridded(geopot_surf_spectral,G)
geopot_surf_fourier = fourier(geopot_surf_grid,G)

@btime fourier_inverse(fourier($geopot_surf_grid,$G),$G);
@btime legendre_inverse(legendre($geopot_surf_fourier,$G),$G);
@btime gridded(spectral($geopot_surf_grid,$G),$G);


text:
# Chebyshev Polynomials of first and second kind
@register0 Chebyshev AbstractCCOP0
export Chebyshev
"""
   Chebyshev{<:Number}(coeffs::AbstractVector, var=:x)

Chebyshev polynomial of the first kind.

Construct a polynomial from its coefficients `a`, lowest order first, optionally in
terms of the given variable `x`. `x` can be a character, symbol, or string.

# Examples

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> Chebyshev([1, 0, 3, 4])
Chebyshev(1⋅T₀(x) + 3⋅T₂(x) + 4⋅T₃(x))

julia> Chebyshev([1, 2, 3, 0], :s)
Chebyshev(1⋅T₀(s) + 2⋅T₁(s) + 3⋅T₂(s))

julia> one(Chebyshev)
Chebyshev(1.0⋅T₀(x))
```

!!! note
    This is copied from the `ChebyshevT` example from the `Polynomials` package by Miles Lucas.

!!! note
    The sample chapter available online of [Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme](https://archive.siam.org/books/ot99/OT99SampleChapter.pdf) gives a very nice overview of these polynomials.
"""
Chebyshev

basis_symbol(::Type{<:Chebyshev})  = "T"
Polynomials.domain(::Type{<:Chebyshev}) = Polynomials.Interval{Open, Open}(-1, 1)
weight_function(::Type{<: Chebyshev}) = x -> one(x)/sqrt(one(x) - x^2)
generating_function(::Type{<: Chebyshev}) =  (t,x) -> (1-t*x)/(1-2*t*x - t^2)
function classical_hypergeometric(P::Type{<:Chebyshev}, n, x) where {α}
    as = (-n,n)
    bs = (one(eltype(P))/2, )
    pFq(as, bs, (1-x)/2)
end

function eval_basis(::Type{Chebyshev}, n, x)
    if abs(x) <= 1
        cos(n * acos(x))
    elseif x > 1
        cosh(n*acosh(x))
    else
        (-1)^n * cosh(n*acosh(-x))
    end
end

abcde(::Type{<:Chebyshev}) = NamedTuple{(:a,:b,:c,:d,:e)}((-1, 0, 1, -1, 0))

k0(P::Type{<:Chebyshev}) = one(eltype(P))
k1k0(P::Type{<:Chebyshev}, n::Int)  = iszero(n) ? one(eltype(P)) : 2*one(eltype(P))

norm2(P::Type{<:Chebyshev}, n::Int) = (iszero(n) ? 1 : 1/2) * pi
ω₁₀(P::Type{<:Chebyshev}, n::Int) = iszero(n) ? one(eltype(P))/sqrt(2) : one(eltype(P)) # √(norm2(n+1)/norm2(n)

leading_term(P::Type{<:Chebyshev}, n::Int)    = iszero(n) ? one(eltype(P)) : (2*one(eltype(P)))^(n-1)

# directly adding these gives a large (20x) speed up in polynomial evaluation
An(P::Type{<:Chebyshev}, n::Int) = iszero(n) ? one(eltype(P)) : 2*one(eltype(P))
Bn(P::Type{<:Chebyshev}, n::Int) = zero(eltype(P))
Cn(P::Type{<:Chebyshev}, n::Int) = one(eltype(P))

#
C̃n(P::Type{<:Chebyshev}, ::Val{1}) = one(eltype(P))/2
ĉ̃n(P::Type{<:Chebyshev}, ::Val{0}) = one(eltype(P))/4
ĉ̃n(P::Type{<:Chebyshev}, ::Val{1}) = Inf
γ̃n(P::Type{<:Chebyshev}, n::Int) = (n==1) ? one(eltype(P))/2 : n*one(eltype(P))/4

function ⊗(::Type{<:Chebyshev}, p1::Chebyshev{T,X}, p2::Chebyshev{S,Y}) where {T,X,S,Y}

    isconstant(p1) &&  return p2 * p1[0]
    isconstant(p2) &&  return p1 * p2[0]
    assert_same_variable(X,Y)

    R = promote_type(T,S)
    z1 = _c_to_z(convert(Vector{R}, p1.coeffs))
    z2 = _c_to_z(convert(Vector{R}, p2.coeffs))
    prod = Polynomials.fastconv(z1, z2)
    ret = Chebyshev(_z_to_c(prod), X)
    return truncate(ret)

end

## Defining p' and ∫dp directly speeds things up, and works around an issue with ĉ
function Polynomials.derivative(p::Chebyshev{T,X,N}, order::Int = 1) where {T,X,N}
    order < 0 && throw(ArgumentError("Order of derivative must be non-negative"))
    R  = eltype(one(T)/1)
    P = Chebyshev{R,X}
    order == 0 && return P(coeffs(p))
    hasnan(p) && return P([NaN])
    order > length(p) && return zero(P)

    q =  P(copy(coeffs(p)))
    n = length(p)
    der = Vector{R}(undef, n)

    for j in n:-1:2
        der[j] = 2j * q[j]
        q[j - 2] += j * q[j] / (j - 2)
    end
    if n > 1
        der[2] = 4q[2]
    end
    der[1] = q[1]

    dp = P(der)
    return order > 1 ?  derivative(dp, order - 1) : dp

end

function Polynomials.integrate(p::P, C::S) where {T, X, P<:Chebyshev{T,X}, S <: Number}
    R = promote_type(eltype(one(T) / 1), S)
    if hasnan(p) || isnan(C)
        return ⟒{P}{R,X}([NaN])
    end
    n = length(p)
    if n == 1
        return ⟒{P}{R,X}([C, p[0]])
    end
    a2 = Vector{R}(undef, n + 1)
    a2[1] = zero(R)
    a2[2] = p[0]
    a2[3] = p[1] / 4
    @inbounds for i in 2:n - 1
        a2[i + 2] = p[i] / (2 * (i + 1))
        a2[i] -= p[i] / (2 * (i - 1))
    end
    a2[1] = C - sum(a2[1+i] for i in 2:2:n)

    return ⟒(P){R,X}(a2)
end

function Polynomials.companion(p::Chebyshev{T}) where {T}
    d = length(p) - 1
    d < 1 && throw(ArgumentError("Series must have degree greater than 1"))
    d == 1 && return diagm(0 => [-p[0] / p[1]])
    R = eltype(one(T) / one(T))

    scl = vcat(1.0, fill(R(√0.5), d - 1))

    diag = vcat(√0.5, fill(R(0.5), d - 2))
    comp = diagm(1 => diag,
                 -1 => diag)
    monics = p.coeffs ./ p.coeffs[end]
    comp[:, end] .-= monics[1:d] .* scl ./ scl[end] ./ 2
    return R.(comp)
end

function Base.divrem(num::Chebyshev{T,X}, den::Chebyshev{S,Y}) where {T,X,S,Y}
    X != Y && throw(ArgumentError("Polynomials must have same variable"))
    n = length(num) - 1
    m = length(den) - 1

    R = typeof(one(T) / one(S))
    P = Chebyshev{R,X}

    if n < m
        return zero(P), convert(P, num)
    elseif m == 0
        den[0] ≈ 0 && throw(DivideError())
        return num ./ den[end], zero(P)
    end

    znum = _c_to_z(num.coeffs)
    zden = _c_to_z(den.coeffs)
    quo, rem = _z_division(znum, zden)
    q_coeff = _z_to_c(quo)
    r_coeff = _z_to_c(rem)
    return P(q_coeff), P(r_coeff)
end

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

# nodes/weights for fitting a polynomial using Lagrange polynomial
# cf. https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
function lagrange_barycentric_nodes_weights(::Type{<: Chebyshev}, n::Int)

    xs = [cospi((j+1/2)/(n+1)) for j in 0:n]
    ws = [(-1)^j*sinpi((j+1/2)/(n+1)) for j in 0:n]  # XX one loop

    xs, ws
end

# # noded for integrating against the weight function
# function gauss_nodes_weights(::Type{<:Chebyshev}, n::Int)
#     xs = cos.(pi/2n * (2*(n:-1:1).-1))
#     ws = pi/n * ones(n)
#     xs, ws
# end

## cf. fastgaussquadrature
#gauss_nodes_weights(p::Type{P}, n) where {P <: Chebyshev} =
#    FastGaussQuadrature.gausschebyshev(n)

##
## fitting coefficients
cks(val::Val{:interpolating},::Type{Chebyshev}, f, n::Int) = cks(val, Chebyshev{Float64}, f, n)
cks(val::Val{:lsq},::Type{Chebyshev}, f, n::Int) = cks(val, Chebyshev{Float64}, f, n)
cks(val::Val{:series},::Type{Chebyshev}, f) = cks(val, Chebyshev{Float64}, f)

# march forward to compute ck; c0 needs division by 1/2, as used
# ck = 2/(n+1) ∑ f(xⱼ) ⋅ T_j(xⱼ)
# The xs are the same as used by `lagrange_barycentric_nodes_weights`, but
# those weights refer to a different basis.
function cks(::Val{:interpolating}, P::Type{Chebyshev{T}}, f, n::Int) where {T}
    R = float(T)
    k = R(0):n
    ks = (k .+ 1/2)/(n+1)
    xks = cospi.(ks)
    fs = f.(xks)
    a = ones(R,n+1)
    b = copy(xks)
    cks = zeros(R, n+1)
    cks[0+1] = dot(a, fs)/2
    cks[1+1] = dot(b, fs)
    for r = 2:n
        for j in 1:n+1
            a[j], b[j] = b[j], muladd(xks[j], 2b[j],- a[j])
        end
        cks[r+1] = dot(b, fs)
    end
    (2/(n+1))*cks
end

## use discrete cosine  transformation to compute the ck.
## (slower than rrecursion, though maybe a bit more accurate, as implemented)
## https://archive.siam.org/books/ot99/OT99SampleChapter.pdf
# dct(P::Type{Chebyshev{T}}, f, n::Int) where  {T} =  [_dct(T,f,k,n) for k in 0:n]
# dct(::Type{Chebyshev}, f, n::Int) = dct(Chebyshev{Float64}, f,n)
# function _dct(T, f, k::Int, n::Int)
#     tot = zero(float(T))
#     for j in 0:n
#         θⱼ = (j+1/2)/(n+1)
#         tot += f(cospi(θⱼ)) * cospi(k * θⱼ)
#     end
#     ck = (2/(n+1)) * tot
#     ck
# end

# march forward to compute ck
# https://archive.siam.org/books/ot99/OT99SampleChapter.pdf
# ck ≈ 2/n ∑'' f(xⱼ) ⋅ T_k(xⱼ)
# xⱼ = cos(jπ/n)
function cks(::Val{:lsq}, P::Type{Chebyshev{T}}, f, n::Int) where {T}
    R = float(T)
    ks = (R(0):n)/n
    xks = cospi.(ks)
    fs = f.(xks)
    fs[1],fs[end] = fs[1]/2, fs[end]/2  # handle end point ''s

    a = ones(R,n+1)
    b = copy(xks)
    out = zeros(R, n+1)

    out[0+1] = dot(a, fs)/2
    out[1+1] = dot(b, fs)

    for r = 2:n
        for j in 1:n+1
            a[j], b[j] = b[j], muladd(xks[j], 2b[j],- a[j])
        end
        out[r+1] = dot(b, fs)
    end
    (2/n)*out
end

## Fit to a series
## use a heuristic to identify `n`
# ref: https://archive.siam.org/books/ot99/OT99SampleChapter.pdf
# Chebyshev interpolation can be interpreted as an approximation to Chebyshev series
# (or vice versa), provided that the coefficients decay fast and the discretization is accurate. In
# other words, Chebyshev series can be a good approximation to near minimax approximations
# (Chebyshev), which in turn are close to minimax approximations.
function cks(::Val{:series}, P::Type{Chebyshev{T}}, f) where {T}

    n = 3
    thresh = 2^n * 10 * eps(T)

    p = fit(Val(:interpolating), P, f, 2^n)

    while n < 10
        if maximum(abs, coeffs(p)[end-2^(n-1):end]) <= thresh
            break
        end
        p = fit(Val(:interpolating), P, f, 2^n)
        n += 1
        thresh *= 2
    end

    coeffs(chop(p, atol= thresh))

end

# ## Chebyshev interpolation of the second kind
# ## uses -1, 1 and roots of U_{n-1},
# ## (slower than recursion)
# ## used discrete cosine  transformation to compute the ck:
# ## https://archive.siam.org/books/ot99/OT99SampleChapter.pdf
# dct1(P::Type{ChebyshevU{T}}, f, n::Int) where  {T} =  [_dct1(T,f,k,n) for k in 0:n]
# dct1(::Type{Chebyshev},f, n::Int) = dct1(Chebyshev{Float64}, f, n)

# function _dct1(T, f, k::Int, n::Int)

#     tot = zero(float(T))
#     # j=0, n separate due to ∑''
#     j=0
#     θⱼ = j/n
#     tot += f(cospi(θⱼ)) * cospi(k*θⱼ)/2

#     j = n
#     θⱼ = j/n
#     tot += f(cospi(θⱼ)) * cospi(k*θⱼ)/2

#     for j in 1:n-1
#         θⱼ = j/n
#         tot += f(cospi(θⱼ)) * cospi(k*θⱼ)
#     end

#     ck = (2/n) * tot
#     ck

# end

##
## --------------------------------------------------
##
@register0 MonicChebyshev AbstractCCOP0
export MonicChebyshev
ϟ(::Type{<:MonicChebyshev}) = Chebyshev
ϟ(::Type{<:MonicChebyshev{T}}) where {T} = Chebyshev{T}
@register_monic(MonicChebyshev)
C̃n(P::Type{<:MonicChebyshev}, ::Val{1}) = one(eltype(P))

##

@register0 OrthonormalChebyshev AbstractCCOP0
export OrthonormalChebyshev
ϟ(::Type{<:OrthonormalChebyshev}) = Chebyshev
ϟ(::Type{<:OrthonormalChebyshev{T}}) where {T} = Chebyshev{T}
@register_orthonormal(OrthonormalChebyshev)

##
##  --------------------------------------------------
##
## Contrib
#=
zseries =#

function _c_to_z(cs::AbstractVector{T}) where {T}
    n = length(cs)
    U = typeof(one(T) / 2)
    zs = zeros(U, 2n - 1)
    zs[n:end] = cs ./ 2
    return zs .+ reverse(zs)
end

function _z_to_c(z::AbstractVector{T}) where {T}
    n = (length(z) + 1) ÷ 2
    cs = z[n:end]
    cs[2:n] *= 2
    return cs
end

function _z_division(z1::AbstractVector{T}, z2::AbstractVector{S}) where {T,S}
    R = eltype(one(T) / one(S))
    length(z1)
    length(z2)
    if length(z2) == 1
        z1 ./= z2
        return z1, zero(R)
    elseif length(z1) < length(z2)
        return zero(R), R.(z1)
    end
    dlen = length(z1) - length(z2)
    scl = z2[1]
    z2 ./= scl
    quo = Vector{R}(undef, dlen + 1)
    i = 1
    j = dlen + 1
    while i < j
        r = z1[i]
        quo[i] = z1[i]
        quo[end - i + 1] = r
        tmp = r .* z2
        z1[i:i + length(z2) - 1] .-= tmp
        z1[j:j + length(z2) - 1] .-= tmp
        i += 1
        j -= 1
    end

    r = z1[i]
    quo[i] = r
    tmp = r * z2
    z1[i:i + length(z2) - 1] .-= tmp
    quo ./= scl
    rem = z1[i + 1:i - 2 + length(z2)]
    return quo, rem
end

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

@register0 ChebyshevU AbstractCCOP0
export ChebyshevU
"""
    ChebyshevU{T}

Implements the
[Chebyshev](https://en.wikipedia.org/wiki/Chebyshev_polynomials)
polynomials of the second kind. These have weight function
`w(x) = sqrt(1-x^2)` over the domain `[-1,1]`.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> p = ChebyshevU([1,2,3])
ChebyshevU(1⋅U₀(x) + 2⋅U₁(x) + 3⋅U₂(x))

julia> convert(Polynomial, p)
Polynomials.Polynomial(-2.0 + 4.0*x + 12.0*x^2)

julia> derivative(p)
ChebyshevU(4.0⋅U₀(x) + 12.0⋅U₁(x))

julia> roots(p)
2-element Vector{Float64}:
 -0.6076252185107651
  0.27429188517743175
```
"""
ChebyshevU

basis_symbol(::Type{<:ChebyshevU})  = "U"
weight_function(::Type{<:ChebyshevU})  = x -> sqrt(one(x) - x^2)
generating_function(::Type{<: ChebyshevU}) = (t, x) -> 1 /  (1 - 2t*x + t^2)
function classical_hypergeometric(P::Type{<:ChebyshevU}, n, x) where {α}
    as = (-n,n+2)
    bs = ((3one(eltype(P)))/2, )
    (n+1)*pFq(as, bs, (1-x)/2)
end

Polynomials.domain(::Type{<:ChebyshevU}) = Polynomials.Interval(-1, 1)

abcde(::Type{<:ChebyshevU}) = NamedTuple{(:a,:b,:c,:d,:e)}((-1, 0, 1, -3, 0))

kn(P::Type{<:ChebyshevU}, n::Int) = (2 * one(eltype(P)))^n
k1k0(P::Type{<:ChebyshevU}, n::Int)  = 2 * one(eltype(P))
k1k_1(P::Type{<:ChebyshevU}, n::Int)  = 4 * one(eltype(P))

norm2(P::Type{<:ChebyshevU}, n::Int) = pi/2
ω₁₀(P::Type{<:ChebyshevU}, n::Int) = one(eltype(P))

# directly adding these gives a 5x speed up in polynomial evaluation
An(::Type{<:ChebyshevU}, n::Int) = 2
Bn(::Type{<:ChebyshevU}, n::Int) = 0
Cn(::Type{<:ChebyshevU}, n::Int) = 1
# work around cancellation
ĉ̃n(P::Type{<:ChebyshevU}, n::Int)  = -one(eltype(P)) /(4n+4)

function ⊗(::Type{<:ChebyshevU}, p::ChebyshevU{T,X}, q::ChebyshevU{S,Y}) where {T,X,S,Y}
    isconstant(p) &&  return q * p[0]
    isconstant(q) &&  return p * q[0]
    assert_same_variable(X,Y)

    M, N = degree(p), degree(q)
    R = promote_type(T, S)
    out = zeros(R, M + N + 1)

    # use Um * Un = sum(U_{m-n+2k} k in 0:n) (wikipedia page)
    for i in 0:M
        pᵢ = p[i]
        for j in 0:N
            qⱼ = q[j]
            pᵢqⱼ = pᵢ*qⱼ
            for k in max(i,j)-min(i,j):2:i+j
                out[k+1] += pᵢqⱼ
            end
        end
    end

    ChebyshevU{R,X}(out)
end

## σ_P = σ_Q:e
Base.convert(::Type{Q}, p::P) where {Q <: Chebyshev, P<: ChebyshevU} = _convert_cop(Q,p)
Base.convert(::Type{Q}, p::P) where {Q <: ChebyshevU, P<: Chebyshev} = _convert_cop(Q,p)

# function Polynomials.integrate(p::ChebyshevU{T}, C::Number=0) where {T}

#     # int Un dx = T_n+1/(n+1)
#     R = eltype(one(T)/1)
#     n = length(p)
#     ts = zeros(R, n+1)
#     for i in 1:n
#         ts[i+1] = p[i-1]/(i)
#     end

#     q = Chebyshev(ts, p.var)
#     q = q - q(0) + R(C)

#     return convert(ChebyshevU, q)
# end

##
## --------------------------------------------------
##
## fitting

## return xs ws
function lagrange_barycentric_nodes_weights(P::Type{<:ChebyshevU}, n::Int)
    xs = cospi.((0:n)/n)
    ws = ones(eltype(P), n+1)
    ws[1] /= 2
    ws[end] /= 2
    for i in  2:2:n+1
        ws[i] *= -1
    end
    xs, ws
end

@register0 MonicChebyshevU AbstractCCOP0
export MonicChebyshevU
ϟ(::Type{<:MonicChebyshevU}) = ChebyshevU
ϟ(::Type{<:MonicChebyshevU{T}}) where {T} = ChebyshevU{T}
@register_monic(MonicChebyshevU)

@register0 OrthonormalChebyshevU AbstractCCOP0
export OrthonormalChebyshevU
ϟ(::Type{<:OrthonormalChebyshevU}) = ChebyshevU
ϟ(::Type{<:OrthonormalChebyshevU{T}}) where {T} = ChebyshevU{T}
@register_orthonormal(OrthonormalChebyshevU)


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	configure_clp(solver_settings_path::String)

Reads user-specified solver settings from clp\_settings.yml in the directory specified by the string solver\_settings\_path.

Returns a MathOptInterface OptimizerWithAttributes Clp optimizer instance to be used in the GenX.generate_model() method.

The Clp optimizer instance is configured with the following default parameters if a user-specified parameter for each respective field is not provided:

 - PrimalTolerance = 1e-7 (Primal feasibility tolerance)
 - DualTolerance = 1e-7 (Dual feasibility tolerance)
 - DualObjectiveLimit = 1e308 (When using dual simplex (where the objective is monotonically changing), terminate when the objective exceeds this limit)
 - MaximumIterations = 2147483647 (Terminate after performing this number of simplex iterations)
 - MaximumSeconds = -1.0	(Terminate after this many seconds have passed. A negative value means no time limit)
 - LogLevel = 1 (Set to 1, 2, 3, or 4 for increasing output. Set to 0 to disable output)
 - PresolveType = 0 (Set to 1 to disable presolve)
 - SolveType = 5 (Solution method: dual simplex (0), primal simplex (1), sprint (2), barrier with crossover (3), barrier without crossover (4), automatic (5))
 - InfeasibleReturn = 0 (Set to 1 to return as soon as the problem is found to be infeasible (by default, an infeasibility proof is computed as well))
 - Scaling = 3 (0 0ff, 1 equilibrium, 2 geometric, 3 auto, 4 dynamic (later))
 - Perturbation = 100 (switch on perturbation (50), automatic (100), don't try perturbing (102))

"""
function configure_clp(solver_settings_path::String)

    solver_settings = YAML.load(open(solver_settings_path))

    # Optional solver parameters ############################################
    # Primal feasibility tolerance
    PrimalTolerance = 1e-7
    if (haskey(solver_settings, "Feasib_Tol"))
        PrimalTolerance = solver_settings["Feasib_Tol"]
    end

    # Dual feasibility tolerance
    DualTolerance = 1e-7
    if (haskey(solver_settings, "Feasib_Tol "))
        DualTolerance = solver_settings["Feasib_Tol"]
    end

    # When using dual simplex (where the objective is monotonically changing), terminate when the objective exceeds this limit
    DualObjectiveLimit = 1e308
    if (haskey(solver_settings, "DualObjectiveLimit"))
        DualObjectiveLimit = solver_settings["DualObjectiveLimit"]
    end

    # Terminate after performing this number of simplex iterations
    MaximumIterations = 2147483647
    if (haskey(solver_settings, "MaximumIterations"))
        MaximumIterations = solver_settings["MaximumIterations"]
    end

    # Terminate after this many seconds have passed. A negative value means no time limit
    MaximumSeconds = -1.0
    if (haskey(solver_settings, "TimeLimit"))
        MaximumSeconds = solver_settings["TimeLimit"]
    end

    # Set to 1, 2, 3, or 4 for increasing output. Set to 0 to disable output
    LogLevel = 1
    if (haskey(solver_settings, "LogLevel"))
        LogLevel = solver_settings["LogLevel"]
    end

    # Set to 1 to disable presolve
    PresolveType = 0
    if (haskey(solver_settings, "Pre_Solve"))
        PresolveType = solver_settings["Pre_Solve"]
    end

    # Solution method: dual simplex (0), primal simplex (1), sprint (2), barrier with crossover (3), barrier without crossover (4), automatic (5)
    SolveType = 5
    if (haskey(solver_settings, "Method"))
        SolveType = solver_settings["Method"]
    end

    # Set to 1 to return as soon as the problem is found to be infeasible (by default, an infeasibility proof is computed as well)
    InfeasibleReturn = 0
    if (haskey(solver_settings, "InfeasibleReturn"))
        InfeasibleReturn = solver_settings["InfeasibleReturn"]
    end

    # 0 off, 1 equilibrium, 2 geometric, 3 auto, 4 dynamic(later)
    Scaling = 3
    if (haskey(solver_settings, "Scaling"))
        Scaling = solver_settings["Scaling"]
    end

    # switch on perturbation (50), automatic (100), don't try perturbing (102)
    Perturbation = 100
    if (haskey(solver_settings, "Perturbation"))
        Perturbation = solver_settings["Perturbation"]
    end
    ########################################################################

    OPTIMIZER = optimizer_with_attributes(
        Clp.Optimizer,
        "PrimalTolerance" => PrimalTolerance,
        "DualTolerance" => DualTolerance,
        "DualObjectiveLimit" => DualObjectiveLimit,
        "MaximumIterations" => MaximumIterations,
        "MaximumSeconds" => MaximumSeconds,
        "LogLevel" => LogLevel,
        "PresolveType" => PresolveType,
        "SolveType" => SolveType,
        "InfeasibleReturn" => InfeasibleReturn,
        "Scaling" => Scaling,
        "Perturbation" => Perturbation,
    )

    return OPTIMIZER
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
const START = 0
const MOV = 1
const LEAVE = 2
const WAITING = 3

mutable struct Node
	id::Int
	status::Int
	home::Array
	position::Array
	current_cell::Array
	seen::Array
	waiting_time::Int
	waited_time::Int
	starting_position::Array
	waiting_upper_bound::Int
	node_seen::Set

	function Node(id::Int,status::Int,home,position::Array,cells_num::Array,upper_bound::Int)
		seen = zeros(Int,cells_num[1],cells_num[2])
		new(id,status,home,position,home,seen,0,0,position,upper_bound,Set())
	end

end

function move(node::Node,canvas::Canvas,alpha::Float64,k::Float64,seconds::Float64)
	if node.status == START  || node.status == LEAVE
		choose_destination(node,canvas,alpha,k)
		procede(node,canvas,seconds)
		return
	elseif node.status == WAITING
		wait_state(node,seconds)
		return
	end
	#Moving state
	procede(node,canvas,seconds)
end

#TODO: What the hell means the distance speed descrition?
function procede(node::Node,canvas::Canvas,seconds::Float64)
	node.status = MOV
	# I'm using a fixed value, but it should be constructed
	speed = 1.5
	distance = get_home_cell(node,canvas).position - node.starting_position
	angle = atan(distance[2],distance[1])
	node.position = node.position + [cos(angle),sin(angle)] * speed * seconds # Make the node procede

	if is_inside_home_cell(node,canvas) || has_surpassed_cell(node,angle,canvas)
		fix_surpassed_cell(node,angle,canvas)
		wait_state(node,seconds)
	end

	if get_cell_position(canvas,node.position) != node.current_cell
		delete_node_from_cell(node,canvas,node.current_cell)
		node.current_cell = get_cell_position(canvas,node.position)
		add_node_to_cell(node,canvas,node.current_cell)
	end
end

#If the angle between the starting position and the actual position is changed, we have surpassed the home cell
function has_surpassed_cell(node,angle,canvas)
	distance = get_home_cell(node,canvas).position - node.position
	return angle != atan(distance[2],distance[1])
end

function fix_surpassed_cell(node,angle,canvas)
	if has_surpassed_cell(node,angle,canvas)
		node.position = get_home_cell(node,canvas).position
	end
end

function choose_destination(node,canvas,alpha,k)
	p = -1;
	current_cell = node.home

	for i in 1:length(canvas.cells)
		for j in 1:length(canvas.cells[i])
			cell = canvas.cells[i][j]
			dist = distance(node,cell,k)
			p_dist = alpha * dist + (1 - alpha) * node.seen[i,j] # Node weight
			if p < p_dist
				p = p_dist
				current_cell = [i,j]
			end
		end
	end
	node.home = current_cell
	node.starting_position = node.position
end

#Distance formula from the paper
function distance(node,cell,k)
	position_distance = round(evaluate(Euclidean(),node.position,cell.position),digits=4)
	return 1 / (((1 + k) * position_distance) ^ 2);
end

function is_inside_home_cell(node,canvas::Canvas)
	home = get_home_cell(node,canvas)
	xl = home.position[1] - (home.l / 2)
	xr = home.position[1] + (home.l / 2)
	yt = home.position[2] + (home.l / 2)
	yb = home.position[2] - (home.l / 2)
	r = Rect(xl,yb,xr,yt)
	return inside((node.position[1],node.position[2]),r)
end

function wait_state(node,seconds)
	if node.status == WAITING
		node.waited_time += seconds
		if node.waiting_time <= node.waited_time
			node.status = LEAVE
		end
	else
		node.status = WAITING
		node.waiting_time = min(round(rand(Levy(1,2))),node.waiting_upper_bound)
		node.waited_time = 0
	end
end

#TODO: This should generate float points
function generate_nodes(canvas::Canvas,n::Int,waiting_upper_bound::Int)
	nodes = []
	cells_num = [ length(canvas.cells),length(canvas.cells[1]) ]
	for i in 1:n
		p = generate_point(canvas)
		node = Node(i,0,get_cell_position(canvas,p),p,cells_num,waiting_upper_bound)
		add_node_to_home_cell(node,canvas)
		push!(nodes,node)
	end
	return nodes
end

function get_home_cell(node::Node,canvas::Canvas)
	return canvas.cells[node.home[1]][node.home[2]]
end

function add_node_to_home_cell(node,canvas)
	push!(get_home_cell(node,canvas).nodes,node)
end

function add_node_to_cell(node,canvas,cell)
	push!(canvas.cells[cell[1]][cell[2]].nodes,node)
end

function delete_node_from_cell(node,canvas,cell)
	delete!(canvas.cells[cell[1]][cell[2]].nodes,node)
end

function update_seen(node::Node,canvas::Canvas,event_writer::EventWriter,time::DateTime)
	node.seen[node.home[1],node.home[2]] = length(get_home_cell(node,canvas).nodes) - 1

	for near_node in get_home_cell(node,canvas).nodes
		if near_node.id == node.id
			continue
		end
		if define_event_and_add(event_writer,node,near_node,time)
			push!(near_node.node_seen,node)
		end
	end

	near_cells = get_near_cells(canvas,node.position)
	for cell in near_cells
		node.seen[cell.x,cell.y] = 0

		for near_node in cell.nodes
			dist = evaluate(Euclidean(),node.position,near_node.position)
			if dist <= canvas.r
				if define_event_and_add(event_writer,node,near_node,time)
					push!(near_node.node_seen,node)
				end
				node.seen[cell.x,cell.y] += 1
			end
		end
	end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module CoreDocs

import ..esc, ..push!, ..getindex, ..unsafe_load, ..Csize_t

function doc!(source::LineNumberNode, mod::Module, str, ex)
    push!(DOCS, (mod, ex, str, source.file, source.line))
end
const DOCS = Array{Any, 1}()

isexpr(x, h) = isa(x, Expr) && x.head === h

lazy_iterpolate(s::AbstractString) = Expr(:call, Core.svec, s)
lazy_iterpolate(x) = isexpr(x, :string) ? Expr(:call, Core.svec, x.args...) : x

function docm(source::LineNumberNode, mod::Module, str, x)
    out = esc(Expr(:call, doc!, QuoteNode(source), mod, lazy_iterpolate(str), Expr(:quote, x)))
    isexpr(x, :module) ? Expr(:toplevel, out, esc(x)) :
    isexpr(x, :call) ? out : Expr(:block, esc(x), out)
end
docm(source::LineNumberNode, mod::Module, x) =
    isexpr(x, :->) ? docm(source, mod, x.args[1], x.args[2].args[2]) : error("invalid '@doc'.")

end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
module Kraken
using Dates

function trans(r)
    return (1103515245 * r + 12435) & 0xFFFFFFFF
end

function crop(r)
    return trunc(Int64, r/256) % 32768
end

function gen_one(index=missing, date=missing, seed_set='a', tld_set_nr=3, tmp=true)
    domain = ""
    if index === missing
        index = trunc(Int64, rand()*10000)
    end
    if date === missing
        date = Dates.Date(Dates.now())
    end
    tld_sets = Dict{Int8,Any}(1=>["com", "net", "tv", "cc"],
                            2=>["dyndns.org", "yi.org", "dynserv.com", "mooo.com"],
                            3=>missing)
    seeds = Dict{Char,Dict}('a'=>Dict{String,Int64}("ex"=>24938314, "nex"=>24938315),
                            'b'=>Dict{String,Int64}("ex"=>1600000, "nex"=>1600001))
    tld = tld_sets[tld_set_nr]

    domain_nr = trunc(Int64, index/2)
    if tmp
        r = 3*domain_nr + seeds[seed_set]["ex"]
    else
        r = 3*domain_nr + seeds[seed_set]["nex"]
    end
    discards = trunc(Int64, (time()-1207000000)/604800) + 2
    if domain_nr % 9 < 8
        if domain_nr % 9 >= 6
            discards -= 1
        end
        for _ in 1:discards
            r = crop(trans(r))
        end
    end

    rands = [0,0,0]
    for i in 1:3
        r = trans(r)
        rands[i] = crop(r)
    end
    len_domain = (rands[1]*rands[2]+rands[3])%6 + 7
    for _ in 1:len_domain
        r = trans(r)
        ch = crop(r)%26 + 97
        domain *= Char(ch)
    end
    if tld !== missing
        domain *= "."*tld[domain_nr%4+1]
    end
    return domain
end

function gen_many(cnt=100)
    res = []
    for i in 1:cnt
        res = [res..., gen_one(i*2,missing,rand(['a','b']),3,rand(Bool))]
    end
    return res
end

end  # module Kraken


text:
module FactCheck2TestSets
using MacroTools, ParserCombinator

# package code goes here
include("utils.jl")
include("jlparse.jl")
include("migration.jl")
include("parsercombinator.jl")
include("comments.jl")
include("sugar.jl")

end # module


text:
# code for Video 5.1:
julia> using RDatasets
julia> iris = dataset("datasets", "iris")
julia> features = convert(Array, iris[:, 1:4]) # measurements or independent variables
# 1) sepal_length, 2) sepal_width, 3) petal_length, 4) petal_width
julia> labels = convert(Array, iris[:, 5])     # targets or dependent variable, here the species name

# general application cycle of an algorithm:
train = # Array of length nrow(iris) and random true / false values
# then the training data set can be constructed as:
features[train]
labels[train]
# and the testing data set is then:
test = !train
features[test]
labels[test]

train  = Array(Bool, nrow(iris)) # this has the right dimension, but contains only false values

julia> using Distributions
julia> rand(Bernoulli(0.75), nrow(iris))
150-element Array{Int64,1}:
 1
 1
 0
 1
 1
 1
 1
 1
 1
 1
 1
 1
 ⋮

julia> train = rand(Bernoulli(0.75), nrow(iris)) .== 1
150-element BitArray{1}:
  true
  true
  true
  true
  true
  true
  true
 false
  true
 false
  true

# training and testing a model:
model = build_model(labels[train], features[train, :])
model = prune_model(model, 0.7)
print_tree(model))
predictions = apply_model(model, features[test,:])

### MLBase:
julia> Pkg.add("MLBase")
julia> using MLBase
## Data preprocessing:
julia> species = ["setosa", "virginica", "versicolor"]
3-element Array{ASCIIString,1}:
 "setosa"
 "virginica"
 "versicolor"

julia> arrs = repeach(species, [2,1,5])
8-element Array{ASCIIString,1}:
 "setosa"
 "setosa"
 "virginica"
 "versicolor"
 "versicolor"
 "versicolor"
 "versicolor"
 "versicolor"

julia> lm = labelmap(species)
LabelMap (with 3 labels):
[1] setosa
[2] virginica
[3] versicolor

julia> labelencode(lm, arrs)
8-element Array{Int64,1}:
 1
 1
 2
 3
 3
 3
 3
 3

julia> groupindices(lm, arrs)
3-element Array{Array{Int64,1},1}:
 [1,2]
 [3]
 [4,5,6,7,8]

# Performance evaluation:
julia> gt = ["setosa", "versicolor", "virginica", "setosa", "virginica"]
5-element Array{ASCIIString,1}:
 "setosa"
 "versicolor"
 "virginica"
 "setosa"
 "virginica"

julia> pred = ["setosa", "virginica", "setosa", "setosa", "virginica"]
5-element Array{ASCIIString,1}:
 "setosa"       # correct predition
 "virginica"    # error predition
 "setosa"       # error predition
 "setosa"       # correct predition
 "virginica"    # correct predition

julia> gtl = labelencode(lm, gt)
5-element Array{Int64,1}:
 1
 3
 2
 1
 2

julia> predl = labelencode(lm, pred)
5-element Array{Int64,1}:
 1
 2
 1
 1
 2

julia> correctrate(gtl, predl)
0.6

julia> errorrate(gtl, predl)
0.4

julia> confusmat(3, gtl, predl)
3x3 Array{Int64,2}:
 2  0  0
 1  1  0
 0  1  0

 ### Receiver Operating Characteristics (ROC)
julia> r = roc(gtl, predl)
MLBase.ROCNums{Int64}
  p = 5     # positive in ground-truth
  n = 0     # negative in ground-truth
  tp = 3    # true positives
  tn = 0    # true negatives
  fp = 0    # false positives
  fn = 0    # false negatives

true_positive_rate(r) # 0.6
true_negative_rate(r) # NaN
precision(r)          # 1.0
recall(r)             # 0.6

# for a real example: see section 2 and code file 5_2.jl

### Cross validation schemes
# for real examples: see section 2 and code file 5_2.j

max_stars_repo_path: Section 5/5_1.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 5/5_1.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 5/5_1.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
#!/usr/bin/env julia
using Lumberjack
using ArgParse

VERSION = "0.2.2"
function parse_commandline()
    s = ArgParseSettings()
    s.epilog = "MentaLiST -- The MLST pipeline developed by the PathOGiST research group. https://github.com/WGS-TB/MentaLiST\n" *
      "To cite: Feijao P, Yao H, Fornika D, Gardy J, Hsiao W, Chauve C, Chindelevitch L. 10/01/2018. Microbial Genomics 4(2): doi:10.1099/mgen.0.000146\n"
    s.version = "MentaLiST $VERSION."
    s.allow_ambiguous_opts= true # to allow -1 and -2.
    s.preformatted_epilog = true
    @add_arg_table s begin
      "call"
        help = "MLST caller, given a sample and a k-mer database."
        action = :command
      "build_db"
        help = "Build a MLST k-mer database, given a list of FASTA files."
        action = :command
      "list_pubmlst"
        help = "List all available MLST schemes from www.pubmlst.org. "
        action = :command
      "download_pubmlst"
        help = "Dowload a MLST scheme from pubmlst and build a MLST k-mer database."
        action = :command
      "list_cgmlst"
        help = "List all available cgMLST schemes from www.cgmlst.org."
        action = :command
      "download_cgmlst"
          help = "Dowload a MLST scheme from cgmlst.org and build a MLST k-mer database."
          action = :command
      "download_enterobase"
          help = "Dowload a MLST scheme from Enterobase (enterobase.warwick.ac.uk) and build a MLST k-mer database."
          action = :command
      "-v", "--version"
          action = :show_version
          help = "show version information and exit"
      end

    # Calling MLST options:
    s["call"].preformatted_epilog = true
    s["call"].epilog = "MentaLiST MLST calling function. Calls alleles on a given MLST database.\n" *
      "You can create a custom DB with 'create_db' or other MentaLiST functions that download schemes from pubmlst, cgmlst.org or Enterobase.\n\n" *
      "Examples:\n" *
      "mentalist call -o my_sample.mlst --db my_scheme.db -1 sample_1.fastq.gz -2 sample_2.fastq.gz # one paired-end sample.\n" *
      "mentalist call -o all_samples.mlst --db my_scheme.db -1 *.fastq.gz -2 *.fastq.gz # multiple paired-end samples."

    @add_arg_table s["call"] begin
        "-o"
          help = "Output file with MLST call"
          arg_type = String
          required = true
        "--db"
          help = "Kmer database"
          required = true
          arg_type = String
        "-t", "--mutation_threshold"
          help = "Maximum number of mutations when looking for novel alleles."
          arg_type = Int
          default = 6
        "--kt"
          help = "Minimum # of times a kmer is seen to be considered present in the sample (solid)."
          arg_type = Int
          default = 10
        "--output_votes"
          help = "Outputs the results for the original voting algorithm."
          action = :store_true
        "--output_special"
          help = "Outputs a FASTA file with the alleles from 'special cases' such as incomplete coverage, novel, and multiple alleles. "
          action = :store_true
        "-i", "--sample_input_file"
          help = "Input TXT file for multiple samples. First column has the sample name, second the FASTQ file. Repeat the sample name for samples with more than one file (paired reads, f.i.)"
        "-1"
          nargs = '*'
          help = "FastQ input files, one per sample, forward reads (or unpaired reads)."
          arg_type = String
        "-2"
          nargs = '*'
          help = "FastQ input files, one per sample, reverse reads."
          arg_type = String
    end

    ## Common option for all db building commands:
    s_db = ArgParseSettings()
    @add_arg_table s_db begin
      "--db"
        help = "Output file (kmer database)"
        arg_type = String
        required = true
      "-k"
        help = "Kmer size"
        required = true
        arg_type = Int8
      "--threads"
        arg_type = Int
        default = 2
        help = "Number of threads used in parallel."
    end

    # Build DB from FASTA, options:
    import_settings(s["build_db"], s_db)
    @add_arg_table s["build_db"] begin
        "-f", "--fasta_files"
            nargs = '+'
            arg_type = String
            help = "Fasta files with the MLST scheme"
            required = true
        "-p", "--profile"
            arg_type = String
            help = "Profile file for known genotypes."
    end

    # Listing functions, common options:
    s_list = ArgParseSettings()
    @add_arg_table s_list begin
      "-p", "--prefix"
      help = "Only list schemes where the species name starts with this prefix."
      arg_type = String
    end

    # List pubmlst
    import_settings(s["list_pubmlst"], s_list)

    # List cgmlst
    import_settings(s["list_cgmlst"], s_list)

    # Common options for MLST download functions:
    s_db_download = ArgParseSettings()
    import_settings(s_db_download, s_db) # import build_db common options
    @add_arg_table s_db_download begin
      "-o", "--output"
        help = "Output folder for the scheme Fasta files."
        arg_type = String
        required = true
      "-s", "--scheme"
        help = "Species name or scheme ID."
        arg_type = String
        required = true
    end

    # Download pubmlst:
    import_settings(s["download_pubmlst"], s_db_download) # import common options

    # Download cgmlst:
    import_settings(s["download_cgmlst"], s_db_download)  # import common options

    # Download enterobase:
    import_settings(s["download_enterobase"], s_db_download) # import common options
    # add (and override) specific enterobase options:
    s["download_enterobase"].error_on_conflict = false  # do not error-out when trying to override an option
    @add_arg_table s["download_enterobase"] begin
      "-s", "--scheme"
        help = "Letter identifying which scheme: (S)almonella, (Y)ersinia, or (E)scherichia/Shigella."
        arg_type = String
        required = true
      "-t", "--type"
        help = "Choose the type: 'cg' or 'wg' for cgMLST or wgMLST scheme, respectively."
        arg_type = String
        required = true
    end

    return parse_args(s)
end

# Check if files exist:
function check_files(files)
  dont_exist = [file for file in files if !isfile(file)]
  if length(dont_exist) > 0
    exit_error("The following input file(s) could not be found: $(join(dont_exist,',')).")
  end
end

### Error FUNCTIONS
function exit_error(msg)
  try
    Lumberjack.error(msg)
  catch
    println("exiting ...")
  end
  exit(-1)
end

#### Main COMMAND functions:
function call_mlst(args)
  run_calling_pipeline(args)
end

function list_pubmlst(args)
  list_pubmlst_schema(args["prefix"])
end

function download_pubmlst(args)
  loci_files, profile_file = download_pubmlst_scheme(args["scheme"], args["output"])
  info("Building the k-mer database ...")
  args["fasta_files"] = loci_files
  args["profile"] = profile_file
  build_db(args)
end

function list_cgmlst(args)
  list_cgmlst_schema(args["prefix"])
end

function download_cgmlst(args)
  loci_files = download_cgmlst_scheme(args["scheme"], args["output"])
  info("Building the k-mer database ...")
  args["fasta_files"] = loci_files
  args["profile"] = nothing
  build_db(args)
end

function download_enterobase(args)
  loci_files = download_enterobase_scheme(args["scheme"], args["type"], args["output"])
  info("Building the k-mer database ...")
  args["fasta_files"] = loci_files
  args["profile"] = nothing
  build_db(args)
end

function build_db(args, version=VERSION)
  # check if files exist:
  check_files(args["fasta_files"])
  # get arguments and call the kmer db builder for each locus:
  k::Int8 = args["k"]
  db_file = args["db"]
  profile = args["profile"]
  info("Opening FASTA files ... ")
  results, loci = kmer_class_for_each_locus(k, args["fasta_files"])
  # Combine results:
  info("Combining results for each locus ...")
  kmer_classification = combine_loci_classification(k, results, loci)

  for (index, fasta_file) in enumerate(args["fasta_files"])
    scheme_fasta_directory = pop!(split(dirname(fasta_file), "/"))
    if scheme_fasta_directory == ""
      scheme_fasta_directory = "."
    end
    args["fasta_files"][index] = scheme_fasta_directory * "/" * basename(fasta_file)
  end
  info("Saving DB ...")
  save_db(k, kmer_classification, loci, db_file, profile, args, version)
  info("Done!")
end

##### Main function: just calls the appropriate commands, with arguments:

args = parse_commandline()
# determine command:
if args["%COMMAND%"] == "call"
  include("calling_functions.jl")
  call_mlst(args["call"])
elseif args["%COMMAND%"] == "build_db"
  addprocs(args["build_db"]["threads"])
  include("build_db_functions.jl")
  build_db(args["build_db"])

elseif args["%COMMAND%"] == "list_pubmlst"
  include("mlst_download_functions.jl")
  list_pubmlst(args["list_pubmlst"])

elseif args["%COMMAND%"] == "download_pubmlst"
  include("mlst_download_functions.jl")
  addprocs(args["download_pubmlst"]["threads"])
  include("build_db_functions.jl")
  download_pubmlst(args["download_pubmlst"])

elseif args["%COMMAND%"] == "list_cgmlst"
  include("mlst_download_functions.jl")
  list_cgmlst(args["list_cgmlst"])

elseif args["%COMMAND%"] == "download_cgmlst"
  include("mlst_download_functions.jl")
  addprocs(args["download_cgmlst"]["threads"])
  include("build_db_functions.jl")
  download_cgmlst(args["download_cgmlst"])

elseif args["%COMMAND%"] == "download_enterobase"
  addprocs(args["download_enterobase"]["threads"])
  include("mlst_download_functions.jl")
  include("build_db_functions.jl")
  download_enterobase(args["download_enterobase"])

end

max_stars_repo_name: matnguyen/MentaLiST

max_issues_repo_name: matnguyen/MentaLiST

max_forks_repo_name: matnguyen/MentaLiST


text:
module Queue

using ...SimpleTasks.Types

export pop_message

"""
    pop_message(queue::QueueService)

Pop a message of the queue.
"""
function pop_message(queue::QueueService)
    error("pop_message for $queue is not implemented")
end

"""
    push_message(queue::QueueService)

Push a message onto the queue
"""
function push_message(queue::QueueService; attributes::Dict{AbstractString,
        Union{AbstractString, Number, Array{UInt8, 1}}} = MessageAtributeType[],
        messsage_body::AbstractString = "")
    error("push_message for $queue is not implemented")
end

end # module Queue

max_forks_repo_name: seung-lab/SimpleTasks


text:
# Density at a given lattice site
density(fᵢ::DenseVector) = sum(fᵢ)
function density(fᵢ::DenseArray)
    shape = size(fᵢ)[2:end]
    reshape([sum(fᵢ[:, i]) for i = 1:prod(shape)], shape)
end
density(sim::Simulation) = density(sim.populations)
# Momentum at given site
function momentum!(μ_out::DenseVector, fᵢ::DenseVector, cᵢ::DenseMatrix)
    @assert length(μ_out) == size(cᵢ, 1) && length(fᵢ) == size(cᵢ, 2)
    μ_out[:] = 0
    for i in 1:length(μ_out)
        for j in 1:length(fᵢ)
            if cᵢ[i, j] == 1
                μ_out[i] += fᵢ[j]
            elseif cᵢ[i, j] == -1
                μ_out[i] -= fᵢ[j]
            end
        end
    end
    μ_out
end
momentum(fᵢ::DenseVector, cᵢ::DenseMatrix) = cᵢ * fᵢ
# Velocities at a given lattice site
velocity(μ::DenseVector, ρ::Number) = μ / ρ
velocity(fᵢ::DenseVector, cᵢ::DenseMatrix, ρ::Number) = velocity(momentum(fᵢ, cᵢ), ρ)
velocity(fᵢ::DenseVector, cᵢ::DenseMatrix) = velocity(fᵢ, cᵢ, density(fᵢ))
function velocity!(ν_out::DenseVector, μ::DenseVector, ρ::Number)
    ν_out[:] = μ
    scale!(ν_out, 1 / ρ)
end

# Momentum and velocity functions that take the full grid, or a simulation object
for name in [:momentum, :velocity]
    @eval begin
        function $name(fᵢ::DenseArray, cᵢ::DenseMatrix)
            const shape = tuple(size(cᵢ, 1), size(fᵢ)[2:end]...)
            result = zeros(eltype(fᵢ), shape)
            for i = 1:prod(shape[2:end])
                result[1:shape[1], i] = $name(fᵢ[:, i], cᵢ)
            end
            result
        end
        $name(sim::Simulation) = $name(sim.populations, sim.lattice.celerities)
    end
end

#= Computes the equilibrium particle distributions $f^{eq}$:

    momentum: Macroscopic momentum μ at current lattice site
    celerities: d by n matrix of celerities ē for the lattice, with d the dimensionality of the
        lattice and n the number of particle distributions.
    weights: Weights associated with each celerity
    ρ: Density

    $f^{eq}$ = weights .* [ρ + 3ē⋅μ + \frac{9}{2ρ} (ē⋅μ)² - \frac{3}{2ρ} μ⋅μ]$
=#
function equilibrium{T, I}(
        ρ::T, momentum::DenseVector{T}, celerities::DenseMatrix{I}, weights::DenseVector{T})
    # computes momentum projected on each particle celerity first
    @assert length(momentum) == size(celerities, 1)
    @assert length(weights) == size(celerities, 2)
    μ_on_ē = celerities.'momentum
    weights .* (
        ρ
        + 3μ_on_ē
        + 9/(2ρ) * (μ_on_ē .* μ_on_ē)
        - 3/(2ρ) * dot(momentum, momentum)
    )
end
equilibrium{T, I}(lattice::Lattice{T, I}, fᵢ::DenseVector{T}) =
    equilibrium(density(fᵢ), momentum(fᵢ, lattice.celerities), lattice.celerities, lattice.weights)

equilibrium{T}(lattice::Lattice, ρ::T, momentum::DenseVector{T}) =
    equilibrium(ρ, momentum, lattice.celerities, lattice.weights)
equilibrium{T}(lattice::Symbol, ρ::T, momentum::DenseVector{T}) =
    equilibrium(getfield(LB, lattice), ρ, momentum)
# Computes equilibrium without memory-allocation
# Incoming must have the right size
function equilibrium!{T, I}(
        feq::DenseVector{T}, ρ::T, momentum::DenseVector{T}, celerities::DenseMatrix{I},
        weights::DenseVector{T})
    @assert length(feq) == size(celerities, 2)
    @assert length(momentum) == size(celerities, 1)
    @assert length(weights) == size(celerities, 2)
    const inv_tworho = 1 / (2ρ)
    # Interestingly enough, celerities[:, i] is fairly slower than celerities[j:j+d], where j is
    # incremented at each step of the loop
    const inc = length(momentum)
    const d = inc - 1
    j = 1
    for i in 1:length(feq)
        const μ_on_ē = 3dot(celerities[j:j+d], momentum)
        j += inc
        feq[i] = weights[i] * (
            ρ + μ_on_ē + (μ_on_ē * μ_on_ē - 3dot(momentum, momentum)) * inv_tworho)
    end
end


text:
# global constant e0 and mass2 is expected
function interactionDynamic(config, qd, τIn, τOut)
    para = config.para

    dτ = abs(τOut - τIn)

    kDiQ = sqrt(dot(qd, qd))
    vd = 4π * e0^2 / (kDiQ^2 + mass2)
    if kDiQ <= para.qgrid.grid[1]
        q = para.qgrid.grid[1] + 1.0e-6
        wd = vd * Grid.linear2D(para.dW0, para.qgrid, para.τgrid, q, dτ)
        # the current interpolation vanishes at q=0, which needs to be corrected!
    else
        wd = vd * Grid.linear2D(para.dW0, para.qgrid, para.τgrid, kDiQ, dτ) # dynamic interaction, don't forget the singular factor vq
    end

    return vd / β, wd
end

function vertexDynamic(config, qd, qe, τIn, τOut)
    vd, wd = interactionDynamic(config, qd, τIn, τOut)
    ve, we = interactionDynamic(config, qe, τIn, τOut)

    return -vd, -wd, ve, we
end


text:
@testset "Solver: Random QP (LU)" begin
    """
        minimize   x' P x + q' x
        subject to    x >= 0
    """

    n = 10

    _P = rand(n)
    P = Diagonal(_P)
    q = rand(n)
    θ = [_P; q]
    z = ones(2 * n)
    r = zeros(2 * n)
    rz = zeros(2 * n, 2 * n)
    rθ = zeros(2 * n, 2 * n)
    κ = [1.0]

    # residual
    function _r!(r, z, θ, κ)
        x = z[1:10]
        y = z[11:20]
        P = θ[1:10]
        q = θ[11:20]

        r[1:10] = 2.0 * Diagonal(P) * x + q - y
        r[11:20] = Diagonal(x) * y .- κ[1]
        nothing
    end

    @variables r_sym[1:20]
    @variables z_sym[1:20]
    @variables θ_sym[1:20]
    @variables κ_sym[1:1]

    parallel = Symbolics.SerialForm()
    _r!(r_sym, z_sym, θ_sym, κ_sym)
    rf! = eval(Symbolics.build_function(r_sym, z_sym, θ_sym, κ_sym,
        parallel = parallel)[2])
    rz_exp = Symbolics.jacobian(r_sym, z_sym)
    rθ_exp = Symbolics.jacobian(r_sym, θ_sym)
    rz_sp = similar(rz_exp, Float64)
    rθ_sp = similar(rθ_exp, Float64)
    rzf! = eval(Symbolics.build_function(rz_exp, z_sym, θ_sym,
        parallel = parallel)[2])
    rθf! = eval(Symbolics.build_function(rθ_exp, z_sym, θ_sym,
        parallel = parallel)[2])

    # options
    opts = RoboDojo.InteriorPointOptions(diff_sol = true, verbose=false)
    idx = RoboDojo.IndicesOptimization(
            2n, 2n,
            [collect(1:n), collect(n .+ (1:n))],
            [collect(1:n), collect(n .+ (1:n))],
            Vector{Vector{Int}}[],
            Vector{Vector{Int}}[],
            collect(1:n),
            collect(n .+ (1:n)),
            collect(1:0),
            Vector{Int}[],
            collect(n .+ (1:n)))

    # solver
    ip = RoboDojo.interior_point(z, θ,
        idx=idx,
        r! = rf!, rz! = rzf!, rθ! = rθf!,
        rz = rz_sp,
        rθ = rθ_sp,
        opts = opts)

    # solve
    status = RoboDojo.interior_point_solve!(ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test all(ip.z[idx.ortz[1]] .> 0.0)
    @test all(ip.z[idx.ortz[2]] .> 0.0)
    @test ip.κ[1] < opts.κ_tol
    @test norm(ip.δz, 1) != 0.0
end


text:
using Test
using Plots

opts = TrajectoryOptimization.SolverOptions()
opts.verbose = true
opts.cost_tolerance = 1e-8
opts.cost_tolerance_intermediate = 1e-8
opts.constraint_tolerance = 1e-8
opts.constraint_tolerance_intermediate = sqrt(opts.constraint_tolerance)
opts.gradient_norm_tolerance = 1e-12
opts.gradient_norm_tolerance_intermediate = 1e-12
opts.use_gradient_aula = false
opts.active_constraint_tolerance = 0.0
opts.penalty_scaling = 10.0
opts.penalty_initial = 1.0
opts.constraint_decrease_ratio = .25
opts.iterations = 500
opts.iterations_outerloop = 10
opts.iterations_innerloop = 300
opts.outer_loop_update_type = :default
opts.use_gradient_aula = false

model, obj = TrajectoryOptimization.Dynamics.double_integrator
u_min = -0.1
u_max = 0.1
obj_con = TrajectoryOptimization.ConstrainedObjective(obj, tf=20.0, u_min=u_min, u_max=u_max)#, x_min=x_min, x_max=x_max)

integrator = :rk4
dt = 0.05
solver = TrajectoryOptimization.Solver(model,obj_con,integration=integrator,dt=dt,opts=opts)
get_num_constraints(solver)
get_num_terminal_constraints(solver)
U = zeros(solver.model.m, solver.N-1)
results, stats = TrajectoryOptimization.solve(solver,U)
solver.opts.square_root = true
results_sqrt, stats_sqrt = TrajectoryOptimization.solve(solver,U)

# S_sqrt = [zeros(nn,nn) for k = 1:N]
cond_normal = zeros(N)
cond_sqrt = zeros(N)

for k = 1:N
    # S_sqrt[k] = results2.S[k]'*results2.S[k]
    cond_normal[k] = cond(results.S[k])
    cond_sqrt[k] = cond(results_sqrt.S[k])
end

plot(cond_normal)
plot!(cond_sqrt)
plot(to_array(results.U)',label="")
plot(to_array(results.X)',label="")

plot(log.(stats["gradient_norm"]))

plot(stats["max_condition_number_S"])
plot!(stats_sqrt["max_condition_number_S"])

plot(stats["max_condition_number"])
plot!(stats_sqrt["max_condition_number"])


text:
"""
Type: OneDAtX

Methods: update!, reinit!, copy

A structure designed to track line search information from one iteration to
another. Given f : ℜⁿ → ℜ, define h(θ) = f(x + θ*d) where x and d
are vectors of same dimension and θ is a scalar, more specifically the step size.

Tracked data can include:
 - x             : the current step size
 - fx [opt]      : h(θ) at the current iteration
 - gx [opt]      : h'(θ)
 - f₀ [opt]      : h(0)
 - g₀ [opt]      : h'(0)
 - d [opt]       : search direction
 - res [opt]     : residual
 - current_time  : the time at which the line search algorithm started.
 - current_score : the score at which the line search algorithm started.

Constructors:
 `OneDAtX(:: T, :: S; fx :: T = _init_field(T), gx :: T = _init_field(T), f₀ :: T = _init_field(T), g₀ :: T = _init_field(T), current_time :: Float64 = NaN) where {S, T <: Number}`

 `OneDAtX(:: T; fx :: T = _init_field(T), gx :: T = _init_field(T), f₀ :: T = _init_field(T), g₀ :: T = _init_field(T), current_time :: Float64 = NaN, current_score :: T = _init_field(T))  where T <: Number`

Note:
 - By default, unknown entries are set using `_init_field`.
 - By default the type of `current_score` is `eltype(x)` and cannot be changed once the State is created.
   To have a vectorized `current_score` of length n, use `OneDAtX(x, Array{eltype(x),1}(undef, n))`.
"""
mutable struct OneDAtX{S, T <: Number} <: AbstractState{S, T}
  x::T
  fx::T  # h(θ)
  gx::T  # h'(θ)
  f₀::T  # h(0)
  g₀::T  # h'(0)

  d::T
  res::T

  current_time::Float64
  current_score::S

  function OneDAtX(
    t::T,
    current_score::S;
    fx::T = _init_field(T),
    gx::T = _init_field(T),
    f₀::T = _init_field(T),
    g₀::T = _init_field(T),
    d::T = _init_field(T),
    res::T = _init_field(T),
    current_time::Float64 = NaN,
  ) where {S, T <: Number}
    return new{S, T}(t, fx, gx, f₀, g₀, d, res, current_time, current_score)
  end
end

function OneDAtX(
  t::T;
  fx::T = _init_field(T),
  gx::T = _init_field(T),
  f₀::T = _init_field(T),
  g₀::T = _init_field(T),
  d::T = _init_field(T),
  res::T = _init_field(T),
  current_time::Float64 = NaN,
  current_score::T = _init_field(T),
) where {T <: Number}
  return OneDAtX(
    t,
    current_score,
    fx = fx,
    gx = gx,
    f₀ = f₀,
    g₀ = g₀,
    d = d,
    res = res,
    current_time = current_time,
  )
end


text:
# This file exemplifies the workflow from data input to optimization result generation
#QUESTION using ClustForOpt_priv.col in module Main conflicts with an existing identifier., using ClustForOpt_priv.cols in module Main conflicts with an existing identifier.

using ClustForOpt
#using ClustForOpt_priv
#using Gurobi

# load data
ts_input_data, = load_timeseries_data("DAM", "GER";K=365, T=24) #DAM

 # run clustering
using Gurobi
env = Gurobi.Env()
clust_res_ar = []
for i=1:2
  push!(clust_res_ar, run_clust(ts_input_data;method="kmedoids_exact",representation="medoid",iterations=5,n_init=1,gurobi_env=env)) # default k-means
end

 # optimization

opt_res = run_opt("battery",clust_res_ar[2].best_results)
 #opt_res = run_opt("gas_turbine",clust_res.best_results[5])

 ###
 # run optimization for all k=1:9
opt_res_all = []
obj=[]
for i=1:2
  push!(opt_res_all,run_opt("battery", clust_res_ar[i].best_results))
  push!(obj,opt_res_all[i].obj)
end
 # run reference case
opt_res_full = run_opt("battery",ts_input_data)
 #using PyPlot
 # figure()
 # plot(obj/opt_res_full.obj)


text:
using FillArrays
export discreteswitchedsystem, DiscreteSwitchedLinearSystem, StateDepDiscreteSwitchedLinearSystem, ConstrainedDiscreteSwitchedLinearSystem

const DiscreteSwitchedLinearSystem = HybridSystem{OneStateAutomaton, <:DiscreteIdentitySystem, <:LinearDiscreteSystem, AutonomousSwitching}
const StateDepDiscreteSwitchedLinearSystem = HybridSystem{<:AbstractAutomaton, <:ConstrainedDiscreteIdentitySystem, <:ConstrainedLinearDiscreteSystem, AutonomousSwitching}
const ConstrainedDiscreteSwitchedLinearSystem = HybridSystem{<:AbstractAutomaton, <:DiscreteIdentitySystem, <:LinearDiscreteSystem, AutonomousSwitching}

"""
    discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix})

Creates the discrete switched linear system defined by
```math
x_{k+1} = A_{\\sigma_k} x_k, \\sigma_k = 1, \\ldots, m.
```
where `m` is the length of `A`.

    discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix}, S::AbstractVector)

Creates the state dependent discrete switched linear system defined by
```math
x_{k+1} = A_{\\sigma_k} x_k, \\sigma_k = 1, \\ldots, m, x_k \\in S[\\sigma_k].
```
where `m` is the length of `A` and `S`.

    discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix}, G::AbstractAutomaton)

Creates the constrained discrete switched linear system defined by
```math
x_{k+1} = A_{\\sigma_k} x_k,
```
where ``\\sigma_1, \\ldots, \\sigma_k`` is a valid sequence of events of the automaton `G`.

    discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix}, G::AbstractAutomaton, S::AbstractVector)

Creates the state-dependent constrained discrete switched linear system defined by
```math
x_{k+1} = A_{\\sigma_k} x_k, x_k \\in S[q_k]
```
where ``q_0, \\sigma_1, q_1, \\ldots, q_{k-1}, \\sigma_k, q_k`` is a valid sequence of events of the automaton `G` with intermediate states ``q_0, \\ldots, q_k``.
"""
function discreteswitchedsystem end
function discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix}, G::AbstractAutomaton=OneStateAutomaton(length(A)); kws...)
    n = nstates(G)
    modes = DiscreteIdentitySystem.(map(s -> _getstatedim(A, G, s), states(G)))
    rm = LinearDiscreteSystem.(A)
    sw = Fill(AutonomousSwitching(), n)
    HybridSystem(G, modes, rm, sw, Dict{Symbol, Any}(kws))
end
function discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix}, G::AbstractAutomaton, S::AbstractVector; kws...)
    n = nstates(G)
    modes = ConstrainedDiscreteIdentitySystem.(map(s -> _getstatedim(A, G, s), states(G)), S)
    guards = _guards(A, G, S)
    rm = ConstrainedLinearDiscreteSystem.(A, guards)
    sw = Fill(AutonomousSwitching(), n)
    HybridSystem(G, modes, rm, sw, Dict{Symbol, Any}(kws))
end
function discreteswitchedsystem(A::AbstractVector{<:AbstractMatrix}, S::AbstractVector; kws...)
    G = LightAutomaton(length(A))
    for σ in 1:length(A)
        for σnext in 1:length(A)
            add_transition!(G, σ, σnext, σ)
        end
    end
    discreteswitchedsystem(A, G, S; kws...)
end

_guards(A, G, S::Fill) = Fill(S[1], length(A))
function _guard(G::AbstractAutomaton, S::AbstractVector, σ)
    s = nothing
    for t in transitions(G)
        if symbol(G, t) == σ
            st = S[source(G, t)]
            if s === nothing
                s = st
            else
                @assert s === st
            end
        end
    end
    @assert s !== nothing
    return s
end
_guards(A, G, S::AbstractVector) = map(σ -> _guard(G, S, σ), 1:length(A))
function _getstatedim(A, G, s)
    d = -1
    function _setcheck(T, i)
        for t in T
            n = size(A[symbol(G, t)], i)
            if d == -1
                d = n
            else
                @assert d == n
            end
        end
    end
    _setcheck(out_transitions(G, s), 2)
    _setcheck(in_transitions(G, s), 1)
    # If d == -1, that means there is no ingoing nor outgoing transition.
    # In that case, using a zero dimension won't hurt.
    d == -1 ? 0 : d
end


text:
for (fn, name) in [(:o, "Open"), (:h, "High"), (:l, "Low"), (:c, "Close"), (:v, "Volume")]
    @eval begin
        function $fn(ta::TimeArray)
            ta[$name]
        end
    end
end


text:
module ConjugateGradient

using LazyAlgebra
using LinearAlgebra
using Printf

"""
    ConjugateGradient.reason(code) -> str

yields a textual description of the code returned by the conjugate gradient
method.

""" reason
const NOT_POSITIVE_DEFINITE = -1
const TOO_MANY_ITERATIONS   =  0
const F_TEST_SATISFIED      =  1
const G_TEST_SATISFIED      =  2
const X_TEST_SATISFIED      =  3
reason(code::Integer) = (
    code == NOT_POSITIVE_DEFINITE ? "LHS operator is not positive definite" :
    code == TOO_MANY_ITERATIONS   ? "too many iterations" :
    code == F_TEST_SATISFIED      ? "function reduction test satisfied" :
    code == G_TEST_SATISFIED      ? "gradient test satisfied" :
    code == X_TEST_SATISFIED      ? "variables change test satisfied" :
    "unknown conjugate gradient result")

struct Context{V}
    p::V # current search direction
    q::V # q = A*p
    r::V # current residuals: r = A*x - b
    z::V # preconditioned residuals: z = M*r
    maxiter::Int
    restart::Int
    fatol::Float64 # absolute tolerance for function reduction
    frtol::Float64 # relative tolerance for function reduction
    gatol::Float64 # absolute tolerance for gradient norm
    grtol::Float64 # relative tolerance for gradient norm
    xatol::Float64 # absolute tolerance for norm of variables change
    xrtol::Float64 # absolute tolerance for norm of variables change
    function Context{V}(p::V, q::V, r::V, z::V=r;
                        maxiter::Integer = typemax(Int),
                        restart::Integer = min(50, length(p)),
                        ftol::NTuple{2,Real} = (0.0, 1e-8),
                        gtol::NTuple{2,Real} = (0.0, 1e-5),
                        xtol::NTuple{2,Real} = (0.0, 1e-6)) where {V}
        maxiter ≥ 0 ||
            bad_argument("bad maximum number of iterations",
                         " (maxiter = ", maxiter, ")")
        restart ≥ 1 ||
            bad_argument("bad number of iterations for restarting",
                         " (restart = ", restart,")")
        ftol[1] ≥ 0 ||
            bad_argument("bad function reduction absolute tolerance",
                         " (ftol[1] = ", ftol[1], ")")
        0 ≤ ftol[2] < 1 ||
            bad_argument("bad function reduction relative tolerance",
                         " (ftol[2] = ", ftol[2], ")")
        gtol[1] ≥ 0 ||
            bad_argument("bad gradient absolute tolerance",
                         " (gtol[1] = ", gtol[1], ")")
        0 ≤ gtol[2] < 1 ||
            bad_argument("bad gradient relative tolerance",
                         " (gtol[2] = ", gtol[2], ")")
        xtol[1] ≥ 0 ||
            bad_argument("bad variables change absolute tolerance",
                         " (xtol[1] = ", xtol[1], ")")
        0 ≤ xtol[2] < 1 ||
            bad_argument("bad variables change relative tolerance",
                         " (xtol[2] = ", xtol[1], ")")
        return new{V}(p, q, r, z, maxiter, restart, ftol[1], ftol[2],
                      gtol[1], gtol[2], xtol[1], xtol[2])
    end
end

function Base.show(io::IO, ctx::Context{V}) where {V}
    print(io, "ConjugateGradient.Context(")
    print(io, "x::$V; # array of size $(size(ctx.r))\n")
    print(io, "    precond = $(ctx.z !== ctx.r),\n")
    print(io, "    maxiter = $(ctx.maxiter),\n")
    print(io, "    restart = $(ctx.restart),\n")
    print(io, "    ftol = ($(ctx.fatol), $(ctx.frtol)),\n")
    print(io, "    gtol = ($(ctx.gatol), $(ctx.grtol)),\n")
    print(io, "    xtol = ($(ctx.xatol), $(ctx.xrtol)))")
end

"""
    ConjugateGradient.Context(x) -> ctx

yields a structure with all parameters and storage for temporary variables
needed for running the linear conjugate gradient algorithm.  Argument `x`
specifies the variables of the problem and is used to allocate temporary
variables by calling `LazyAlgebra.vcreate(x)`. The returned context holds no
references on `x`.  Algorithm parameters are specified by the following
keywords:

* `precond` is to specify whether to allocate temporary variables to store the
  preconditioned residuals.  By default, `precond = false`.  If false, only the
  un-preconditioned version of the algorithm can be run.

* `maxiter` is to specify the maximum number of iterations to perform which is
  practically unlimited by default.

* `restart` is to specify the number of consecutive iterations before
  restarting the conjugate gradient recurrence.  Restarting the algorithm is to
  cope with the accumulation of rounding errors.  By default, `restart =
  min(50,length(x)+1)`.  Set `restart` to a value less or equal zero or greater
  than `maxiter` if you do not want that any restarts ever occur.

* `ftol = (fatol,frtol)` is to specify the absolute and relative tolerances for
  the function reduction.  By default, `ftol = (0.0,1e-8)`.

* `gtol = (gatol,grtol)` is to specify the absolute and relative tolerances for
  stopping the algorithm based on the gradient of the objective function.
  Convergence occurs when the Mahalanobis norm of the residuals (which is that
  of the gradient of the associated objective function) is less or equal the
  largest of `gatol` and `grtol` times the Mahalanobis norm of the initial
  residuals.  By default, `gtol = (0.0,1e-5)`.

* `xtol = (xatol,xrtol)` is to specify the absolute and relative tolerances for
  the change in variables.  By default, `xtol = (0.0,1e-6)`.

Just recall the constructor with an instance of `ConjugateGradient.Context`
to re-use the same temporary variables (if possible) but different parameters:

    ConjugateGradient.Context(ctx; precond=…, maxiter=…, ...)

"""
function Context(x::V; precond::Bool = false, kwds...) where {V}
    p = vcreate(x)
    q = vcreate(x)
    r = vcreate(x)
    z = (precond ? vcreate(x) : r)
    return Context{V}(p, q, r, z; kwds...)
end

function Context(ctx::Context{V};
                 precond::Bool = false,
                 maxiter::Integer = ctx.maxiter,
                 restart::Integer = ctx.restart,
                 ftol::NTuple{2,Real} = (ctx.fatol, ctx.frtol),
                 gtol::NTuple{2,Real} = (ctx.gatol, ctx.grtol),
                 xtol::NTuple{2,Real} = (ctx.xatol, ctx.xrtol)) where {V}
    return Context{V}(ctx.p,
                      ctx.q,
                      ctx.r,
                      (precond && ctx.z === ctx.r ? vcreate(ctx.r) : ctx.z);
                      maxiter = maxiter,
                      restart = restart,
                      ftol = ftol,
                      gtol = gtol,
                      xtol = xtol)
end

"""
    ConjugateGradient.solve!(x, A, b[, M], ctx[, io=devnull]) -> status

runs the (preconditioned) linear conjugate gradient algorithm to solve the
system of equations `A*x = b` in `x` and according to the settings in `ctx`.

Argument `x` stores the initial solution on entry and the estimated solution on
return.

Argument `A` implements the *left-hand-side (LHS) matrix* of the equations.  It
is used as `LazyAlgebra.vmul!(dst,A,src)` to store in `dst` the result of
applying `A` to `src` and where `src` and `dst` are similar to arguments `x`
and `b`.  If none of these is suitable, the method `LazyAlgebra.vmul!` can be
extended.  Note that, as `A` and `M` must be symmetric, it may be faster to
apply their adjoint.

Argument `b` is the *right-hand-side (RHS) vector* of the equations.  It is
left unchanged.

Argument `ctx` is a [`ConjugateGradient.Context`](@ref) structure storing all
temporary variables and parameters of the algorithm.  This argument is reusable
and is required to avoid any additional allocations.

Optional argument `M` is a preconditioner.  If `M` is unspecified or if `M` is
`LazyAlgebra.Identity()`, the unpreconditioned version of the algorithm is run.
The preconditioner can be specified in various forms (as for the LHS operator
`A`).

Optional argument `io` is to specify an `IO` instance to which print various
information at each iterations.  Nothing is printed if `io` is `devnull` which
is the default.

## Convergence criteria

Provided `A` be positive definite, the solution `x` of the equations `A*x = b`
is unique and is also the minimum of the following convex quadratic objective
function:

    f(x) = (1/2)*x'*A*x - b'*x + ϵ

where `ϵ` is an arbitrary constant.  The gradient of this objective function
is:

    ∇f(x) = A*x - b

hence solving `A*x = b` for `x` yields the minimum of `f(x)`.  The variations
of `f(x)` between successive iterations, the norm of the gradient `∇f(x)` or
the norm of the variation of variables `x` may be used to decide the
convergence of the algorithm (see keywords `ftol`, `gtol` and `xtol` below).

Let `x_{k}`, `f_{k} = f(x_{k})` and `∇f_{k} = ∇f(x_{k})` denote the variables,
the objective function and its gradient at iteration `k`.  The argument `x`
gives the initial variables `x_{0}`.  Starting with `k = 0`, the different
possibilities for the convergence of the algorithm are listed below.

* The convergence in the function reduction between succesive iterations occurs
  at iteration `k ≥ 1` if:

  ```
  f_{k-1} - f_{k} ≤ max(fatol, frtol*max_{k' ≤ k}(f_{k'-1} - f_{k'}))
  ```

* The convergence in the gradient norm occurs at iteration `k ≥ 0` if:

  ```
  ‖∇f_{k}‖_M ≤ max(gatol, grtol*‖∇f_{0}‖_M)
  ```

  where `‖u‖_M = sqrt(u'*M*u)` is the Mahalanobis norm of `u` with precision
  matrix `M` which is equal to the usual Euclidean norm of `u` if ; no
  preconditioner is used or if `M` is the identity.

* The convergence in the variables occurs at iteration `k ≥ 1` if:

  ```
  ‖x_{k} - x_{k-1}‖ ≤ max(xatol, xrtol*‖x_{k}‖)
  ```

In the conjugate gradient algorithm, the objective function is always reduced
at each iteration, but be aware that the gradient and the change of variables
norms are not always reduced at each iteration.

## Returned Status

The returned value `status` is negative to indicate an error, zero to indicate
that the maximum number of iterations have been reached, positive to indicate
that the method has converged.  More specfically, `status` is one of:

- `ConjugateGradient.NOT_POSITIVE_DEFINITE` if the left-hand-side matrix `A` is
  found to be not positive definite;

- `ConjugateGradient.TOO_MANY_ITERATIONS` if the maximum number of iterations
  have been reached;

- `ConjugateGradient.F_TEST_SATISFIED` if convergence occured because the
  function reduction satisfies the criterion specified by `ftol`;

- `ConjugateGradient.G_TEST_SATISFIED` if convergence occured because the
  gradient norm satisfies the criterion specified by `gtol`;

- `ConjugateGradient.X_TEST_SATISFIED` if convergence occured because the norm
  of the variation of variables satisfies the criterion specified by `xtol`.

Method [`ConjugateGradient.reason`](@ref) may be called to get a textual
explanation about the returned status.

"""
function solve!(x::V, A, b::V, ctx::Context{V},
                io::IO = devnull) where {V}
    # Run the unpreconditioned version of the algorithm.
    return solve!(x, A, b, Identity(), ctx, io)
end

function solve!(x::V, A, b::V, M, ctx::Context{V},
                io::IO = devnull) where {V}
    # Get workspace variables.
    precond = !isa(M, Identity)
    verbose = !isa(io, Base.DevNull)
    p, q, r = ctx.p, ctx.q, ctx.r
    z = (precond ? ctx.z : ctx.r)
    if precond && z === r
        error("workspace variables Z must be different from R with ",
              "a preconditioner")
    end

    # Enforce types of some variables (FIXME: This will not work for BigFloat).
    local rho::Float64, oldrho::Float64
    local alpha::Float64, gamma::Float64
    local psi::Float64, psimax::Float64
    local gtest::Float64

    # Initialize local variables.
    rho = 0.0
    psi = 0.0
    psimax = 0.0
    xtest = (ctx.xatol > 0 || ctx.xrtol > 0)
    if verbose
        t0 = time()
    end

    # Conjugate gradient iterations.
    k = 0
    while true
        # Is this the initial or a restarted iteration?
        restart = (k == 0 || ctx.restart > 0 && rem(k, ctx.restart) == 0)

        # Compute residuals and their squared norm.
        if restart
            # Compute residuals.
            if k > 0 || vnorm2(x) != 0
                # Compute r = b - A*x.
                vcombine!(r, 1, b, -1, vmul!(r, A, x))
            else
                # Spare applying A since x = 0.
                vcopy!(r, b)
            end
        else
            # Update residuals.
            vupdate!(r, -alpha, q)
        end
        if precond
            # Apply preconditioner.
            vmul!(z, M, r) # z = M*r
        end
        oldrho = rho
        rho = vdot(r, z) # rho = ‖r‖_M^2
        if k == 0
            gtest = tolerance(ctx.gatol, ctx.gatol, sqrt(rho))
        end
        if verbose
            t = (time() - t0)*1E3 # elapsed time in ms
            if precond
                if k == 0
                    print(io,
                          "# Iter.   Time (ms)     Δf(x)       ‖∇f(x)‖  ",
                          "   ‖∇f(x)‖_M\n",
                          "# --------------------------------------------",
                          "-------------\n")
                end
                @printf(io, "%7d %11.3f %12.4e %12.4e %12.4e\n",
                        k, t, psi, vnorm2(r), sqrt(rho))
            else
                if k == 0
                    print(io,
                          "# Iter.   Time (ms)     Δf(x)       ‖∇f(x)‖\n",
                          "# --------------------------------------------\n")
                end
                @printf(io, "%7d %11.3f %12.4e %12.4e\n",
                        k, t, psi, sqrt(rho))
            end
        end
        if sqrt(rho) ≤ gtest
            # Normal convergence in the gradient norm.
            verbose && println(io, "# Convergence in the gradient norm.")
            return G_TEST_SATISFIED
        end
        if k ≥ ctx.maxiter
            verbose && println(io, "# Too many iteration(s).")
            return TOO_MANY_ITERATIONS
        end

        # Compute search direction.
        if restart
            # Restarting or first iteration.
            vcopy!(p, z)
        else
            # Apply recurrence.
            beta = rho/oldrho
            vcombine!(p, 1, z, beta, p)
        end

        # Compute optimal step size.
        vmul!(q, A, p)
        gamma = vdot(p, q)
        if !(gamma > 0)
            verbose && println(io, "# Operator is not positive definite.")
            return NOT_POSITIVE_DEFINITE
        end
        alpha = rho/gamma

        # Update variables and check for convergence.
        vupdate!(x, +alpha, p)
        psi = alpha*rho/2  # psi = f(x_{k}) - f(x_{k+1}) ≥ 0
        psimax = max(psi, psimax)
        if psi ≤ tolerance(ctx.fatol, ctx.frtol, psimax)
            # Normal convergence in the function reduction.
            verbose && println(io, "# Convergence in the function reduction.")
            return F_TEST_SATISFIED
        end
        if xtest && alpha*vnorm2(p) ≤ tolerance(ctx.xatol, ctx.xrtol, x)
            # Normal convergence in the variables.
            verbose && println(io, "# Convergence in the variables.")
            return X_TEST_SATISFIED
        end

        # Increment iteration number.
        k += 1
    end
end

"""

Given absolute and relative tolerances `atol` and `rtol` (both finite and
nonnegative), the calls:

    tolerance(atol, rtol, val) -> max(0, atol, rtol*abs(val))
    tolerance(atol, rtol, arr) -> max(0, atol, rtol*vnorm2(arr))

yield the tolerances for the scalar `val` or for the array `arr`.  The result
is nonnegative.

If `rtol ≤ 0`, the computation of `vnorm2(arr)` is not performed.

"""
tolerance(atol::Real, rtol::Real, val::Real) =
    tolerance(promote(atol, rtol, val)...)

tolerance(atol::T, rtol::T, val::T) where {T<:Real} =
    max(zero(T), atol, rtol*abs(val))

function tolerance(atol::Real, rtol::Real,
                   arr::AbstractArray{T}) where {T<:AbstractFloat}
    # NOTE: This is to ensure type stability.
    val = (rtol > 0 ? vnorm2(arr)::T : zero(T))
    tolerance(atol, rtol, val)
end

"""
    bad_argument(args...)

throws an `ArgumentError` exception with `args...` converted into a string
as error message.

"""
bad_argument(msg::AbstractString) = throw(ArgumentError(msg))
@noinline bad_argument(args...) = bad_argument(string(args...))

end # module


text:
using JUDI.TimeModeling, PyPlot

## Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) + 0.4f0
v0 = ones(Float32,n) + 0.4f0
v[:,Int(round(end/2)):end] = 3f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

## Set up receiver geometry
nxrec = 120
xrec = linspace(50f0, 1150f0, nxrec)
yrec = 0f0
zrec = linspace(50f0, 50f0, nxrec)

# receiver sampling and recording time
timeR = 1000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(600f0)
ysrc = convertToCell(0f0)
zsrc = convertToCell(20f0)

# source sampling and number of time steps
timeS = 1000f0  # ms
dtS = 4f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.01f0     # MHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

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

# Write shots as segy files to disk
opt = Options(save_data_to_disk=false, file_path=pwd(), file_name="observed_shot", optimal_checkpointing=true)

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*Ps', q)

# Pr*F*Ps'*q
d_obs = Pr*F*Ps'*q
q_ad = Ps*F*Pr'*d_obs

# time_modeling(model, src_geometry, src_data, rec_geometry, rec_data, dm, nsrc, op, mode)

# Data in, wavefields out
u = time_modeling(model, q.geometry, q.data, nothing, nothing, nothing, 1, 'F', 1, Options())
v = time_modeling(model, nothing, nothing, d_obs.geometry, d_obs.data, nothing, 1, 'F', -1, Options())

# Wavefield in, data out
d_q = time_modeling(model, nothing, v.data, recGeometry, nothing, nothing, 1, 'F', 1, Options())
#qad_q = time_modeling(model, q.geometry, nothing, nothing, u, nothing, 1, 'F', 1, Options())

# Ps*F'*Pr'*d
#qad  = time_modeling(model, q.geometry, [], recGeometry, dobs.data, [], 1, 'F', -1)

# F*Ps'*q
#wf = time_modeling(mode, q.geometry, q.data, [], [], [], 1, 'F', 1)


text:
function LazyLoadedGeoMagneticFields(arg0::DataProvidersManager)
    return LazyLoadedGeoMagneticFields((DataProvidersManager,), arg0)
end

function get_field(obj::LazyLoadedGeoMagneticFields, arg0::GeoMagneticFieldFactory_FieldModel, arg1::jdouble)
    return jcall(obj, "getField", GeoMagneticField, (GeoMagneticFieldFactory_FieldModel, jdouble), arg0, arg1)
end

function get_igrf(obj::LazyLoadedGeoMagneticFields, arg0::jdouble)
    return jcall(obj, "getIGRF", GeoMagneticField, (jdouble,), arg0)
end

function get_wmm(obj::LazyLoadedGeoMagneticFields, arg0::jdouble)
    return jcall(obj, "getWMM", GeoMagneticField, (jdouble,), arg0)
end


text:
using NBInclude
using Base.Test

try
    using POMDPs
    POMDPs.add("DiscreteValueIteration")
    POMDPs.add("MCTS")

    using DiscreteValueIteration
    using MCTS

    @nbinclude("../examples/GridWorld.ipynb")
catch ex
    @warn "The Grid World tutorial notebook failed with the following error:"
    showerror(STDERR, ex)
    @warn "This failure is being ignored."
end

try
    @nbinclude("../examples/Tiger.ipynb")
catch ex
    @warn "The Tiger tutorial notebook failed with the following error:"
    showerror(STDERR, ex)
    @warn "This failure is being ignored."
end

try
    @nbinclude("../examples/rl-tuto/reinforcement_learning_tutorial.ipynb")
catch ex
    @warn "The RL tutorial notebook failed with the following error:"
    showerror(STDERR, ex)
    @warn "This failure is being ignored."
end


text:
"""
    karger_min_cut(g)

Perform [Karger Minimum Cut](https://en.wikipedia.org/wiki/Karger%27s_algorithm)
to find the minimum cut of graph `g` with some probability of success.
A cut is a partition of `vertices(g)` into two non-empty sets.
The size of a cut is the number of edges crossing the two non-empty sets.

### Implementation Notes
The cut is represented by an integer array.
If `cut[v] == 1` then `v` is in the first non-empty set.
If `cut[v] == 2` then `v` is in the second non-empty set.
`cut[1] = 1`.

If |V| < 2 then `cut[v] = 0` for all `v`.

### Performance
Runtime: O(|E|)
Memory: O(|E|)
"""
function karger_min_cut(g::AbstractGraph{T}) where T <: Integer

    nvg = nv(g)
    nvg < 2 && return zeros(Int, nvg)
    nvg == 2 && return [1, 2]

    connected_vs = IntDisjointSets(nvg)
    num_components = nvg

    for e in shuffle(collect(edges(g)))
        s = src(e)
        d = dst(e)
        if !in_same_set(connected_vs, s, d)
            union!(connected_vs, s, d)
            num_components -= one(T)
            (num_components <= 2) && break
        end
    end

    return [(in_same_set(connected_vs, one(T), v) ? 1 : 2) for v in vertices(g)]
end

"""
    karger_cut_cost(g, cut)

Find the number of crossing edges in the a cut of graph `g` where the cut is represented
by the integer array, `cut`.
"""
karger_cut_cost(g::AbstractGraph{T}, cut::Vector{<:Integer}) where T <: Integer =
count((e::Edge{T})->cut[src(e)] != cut[dst(e)], edges(g))

"""
    karger_cut_edges(g, cut)

Find the crossing edges the a cut of graph `g` where the cut is represented
by the integer array, `cut`.
"""
karger_cut_edges(g::AbstractGraph{T}, cut::Vector{<:Integer}) where T <: Integer =
[e for e in edges(g) if cut[src(e)] != cut[dst(e)]]