text:
#!/usr/bin/julia

import Base.print;
import Base.+;

#
## Boolean type
#
mutable struct SBool
    value::Bool;

    function SBool(b::Bool)
        this = new();
        this.value = b;
        return this;
    end
end

#
## String type
#
mutable struct SString
    value::AbstractString;

    function SString(s::AbstractString)
        this = new();
        this.value = s;
        return this;
    end
end

#
## String methods
#
function say(s::SString)
    println(s.value);
    SBool(true);
end

function print(s::SString)
    print(s.value);
    SBool(true);
end

function +(a::SString, b::SString)
    SString(a.value * b.value);
end

#
## Expression executor
#
function execute_expr(statement)
    if haskey(statement, :self)
        self = statement[:self];

        if isa(self, Dict)
            self = execute(self);
        end

        if haskey(statement, :call)
            for call in statement[:call]
                method = Symbol(call[:method]);

                if haskey(call, :arg)

                    # Evaluate and collect the arguments
                    args = Any[];
                    for arg in call[:arg]
                        if isa(arg, Dict)
                            push!(args, execute_expr(arg))
                        else
                            push!(args, arg);
                        end
                    end

                    # Call method with arguments
                    self = eval(:($method($self, $args...)));
                else
                    # Call method without arguments
                    self = eval(:($method($self)));
                end
            end
        end

        return(self);
    end
end

#
## Parse-tree executor
#
function execute(structure)
    results = Any[];
    for statement in structure[:main]
        push!(results, execute_expr(statement));
    end
    results[end];
end

#
## The AST
#
structure = Dict{Symbol, Any}(
    :main => [
        Dict{Symbol, Any}(
            :call => [Dict{Symbol, Any}(:method => "print")],
            :self => Dict{Symbol, Any}(
                :main => [
                    Dict{Symbol, Any}(
                        :call => [Dict{Any, Any}(:method => "+", :arg => [Dict{Symbol, Any}(:self => SString("llo"))])],
                        :self => SString("he"),
                    )
                ],
            )
        ),
        Dict{Symbol, Any}(
            :call => [Dict{Symbol, Any}(:method => "say")],
            :self => SString(" world!");
        ),
    ]
);

#
## Begin execution
#
execute(structure);


text:
using Test, AxisKeys, NamedDims
using Statistics, OffsetArrays, Tables, UniqueVectors, LazyStack

@testset "wrapdims with $outside outside" for outside in [:NamedDimsArray, :KeyedArray]

    if outside == :NamedDimsArray
        AxisKeys.nameouter() = true
    else
        AxisKeys.nameouter() = false
    end

    # include("_basic.jl")

    include("_functions.jl")

    # include("_fast.jl")

    # include("_packages.jl")

end
@testset "fast findfirst & findall" begin

    include("_findrange.jl")

end
@testset "ambiguities" begin

    @test isempty(detect_unbound_args(AxisKeys))
    alist = detect_ambiguities(AxisKeys, Base)
    @test_broken isempty(alist)

end

ext: jl


text:
function legendre_pos_m(lmax, x, Plm)
    """ Associated Legendre polynomials only with positive M. """

    a = (1. - x^2)
    b = sqrt(a)
    Plm[1] = 1.
    if lmax > 0
        Plm[2] = x
        Plm[3] = - b
        if lmax > 1
            Plm[4] = 0.5 * (3. * x^2 - 1.)
            Plm[5] = -3. * x * b
            Plm[6] =  3. * a

            if lmax > 2
                Plm[7] = 0.5 * (5. * x^3 - 3. * x)
                Plm[8] = -3. / 2. * (5. * x^2 - 1. ) * b
                Plm[9] =  15. * x * a
                Plm[10] =  -15. * a * b
            end
        end
    end

    i_lminus = 4
    i_l = 7
    for l=4:lmax
        i_lplus = l*(l+1) ÷ 2 + 1

        Plm[i_lplus] = ((2*l-1) * x * Plm[i_l] - (l-1) * Plm[i_lminus]) / l

        for m=1:l
            @inbounds Plm[i_lplus+m] = ((l+1-m) * x * Plm[i_lplus+m-1] - (l+m-1) * Plm[i_l+m-1]) / b
        end

        i_lminus = i_l
        i_l = i_lplus

    end
    return nothing
end

function fact_quotient(l, m)
    if 2m > 21
        fact = factorial(big(2m))
    else
        fact = factorial(2m)
    end
    return 1/(binomial(l+m, l-m) * fact)
end
function get_ylm(lmax)
    ylm = Array{Float64}(undef, lmax+1, lmax+1)
    for l=0:lmax for m=0:lmax
        if m == 0
            ylm[l+1, m+1] = sqrt((2l +1) / 4pi)
        else
            # abs needed as very small numbers could become negative
            ylm[l+1, m+1] = (-1)^m * sqrt(abs((2l + 1) / 2pi * fact_quotient(l,  m)))
        end

    end end
    return ylm
end

function get_sincos(lmax, phi, sincos)
    for m=1:lmax
        sincos[1, m]= sin(m * phi)
        sincos[2, m]= cos(m * phi)
    end
    return nothing
end

function spherical_harmonics_series(lmax, x, y, z, Plm, Ylm, ylm, sincos)
    theta = atan(sqrt(x^2 + y^2), z)
    phi = atan(y, x)
    legendre_pos_m(lmax, cos(theta), Plm)
    get_sincos(lmax, phi, sincos)
    c = 1
    for l=0:lmax
        i_l = (l+1)^2 - l

        @inbounds Ylm[i_l] = ylm[l+1, 1] * Plm[c]
        c += 1
        for m=1:l
            @inbounds Ylm[i_l-m] = ylm[l+1, m+1] * Plm[c] * sincos[1, m]
            @inbounds Ylm[i_l+m] = ylm[l+1, m+1] * Plm[c] * sincos[2, m]
            c+=1
        end
    end
    return nothing
end

function sum_environmental_to_structural2(d, r_cut, cutoff_width, lmax, nmax, shift_3d)
    n = size(d, 1)
    d = reshape(d, n, 3)
    n_spherical = (lmax+1)^2
    g = zeros( n_spherical * (nmax+1))
    r_cut_start = r_cut - cutoff_width

    Plm = Array{Float64}(undef, lmax*(lmax+1) ÷ 2 + lmax+1)
    Ylm = Array{Float64}(undef, (lmax+1)^2)
    ylm = get_ylm(lmax)
    sincos = Array{Float64}(undef, 2, lmax)
    polynomials = ones(Float64, nmax+1)

    for i=1:n
        @inbounds d1 = d[i, 1]
        @inbounds d2 = d[i, 2]
        @inbounds d3 = d[i, 3]
        @inbounds r1 = d1 - shift_3d[1]
        @inbounds r2 = d2 - shift_3d[2]
        @inbounds r3 = d3 - shift_3d[3]

        r_radial = sqrt(r1^2 + r2^2 + r3^2)
        spherical_harmonics_series(lmax, r1, r2, r3, Plm, Ylm, ylm, sincos)

       # polynomials
        @inbounds polynomials[1] = 1
        for k=1:nmax
            @inbounds polynomials[k+1] = polynomials[k] * r_radial
        end

        c = 1
        if d1 > r_cut_start || d2 > r_cut_start
            if d1 > r_cut_start && d2 > r_cut_start
                f_cut = 0.25 * (cos(pi / cutoff_width * (d1 - r_cut_start)) + 1.) * (cos(pi / cutoff_width * (d2 - r_cut_start)) + 1.)
            elseif d1 > r_cut_start
                f_cut = 0.5 * (cos(pi / cutoff_width * (d1 - r_cut_start)) + 1.)
            else
                f_cut = 0.5 * (cos(pi / cutoff_width * (d2 - r_cut_start)) + 1.)
            end

            for j=1:n_spherical for k=1:nmax+1
                @inbounds g[c] += Ylm[j] * polynomials[k] * f_cut
                c+=1
            end end
        else
            for j=1:n_spherical for k=1:nmax+1
                @inbounds g[c] += Ylm[j] *  polynomials[k]
                c+=1
            end end
        end
    end
    return g
end

function sum_environmental_to_structural2(d, r_cut, cutoff_width, r_centers, sigmas)
    n = size(d, 1)
    d = reshape(d, n, 3)
    s = size(sigmas, 1)
    g = zeros( s)
    r_cut_start = r_cut - cutoff_width

    for i=1:n
        @inbounds d1 = d[i, 1]
        @inbounds d2 = d[i, 2]
        @inbounds d3 = d[i, 3]

        c = 1
        if d1 > r_cut_start || d2 > r_cut_start
            if d1 > r_cut_start && d2 > r_cut_start
                f_cut = 0.25 * (cos(pi / cutoff_width * (d1 - r_cut_start)) + 1.) * (cos(pi / cutoff_width * (d2 - r_cut_start)) + 1.)
            elseif d1 > r_cut_start
                f_cut = 0.5 * (cos(pi / cutoff_width * (d1 - r_cut_start)) + 1.)
            else
                f_cut = 0.5 * (cos(pi / cutoff_width * (d2 - r_cut_start)) + 1.)
            end

            for j=1:s
                @inbounds g[j] += exp(-0.5 * ((d1-r_centers[j, 1])^2 + (d2-r_centers[j, 2])^2 + (d3-r_centers[j, 3])^2) / sigmas[j]^2) * f_cut
            end
        else

            for j=1:s
                @inbounds g[j] += exp(-0.5 * ((d1-r_centers[j, 1])^2 + (d2-r_centers[j, 2])^2 + (d3-r_centers[j, 3])^2) / sigmas[j]^2)

            end
        end
    end
    return g
end


text:
using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg = 3, element_type = Tet(),
             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
             volume_integral = VolumeIntegralWeakForm())

equations = CompressibleEulerEquations3D(1.4)
initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test

# example where we tag two separate boundary segments of the mesh
top_boundary(x, tol=50*eps()) = abs(x[2] - 1) < tol
rest_of_boundary(x, tol=50*eps()) = !top_boundary(x, tol)
is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)
vertex_coordinates, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, 4)
mesh = VertexMappedMesh(vertex_coordinates, EToV, dg, is_on_boundary = is_on_boundary)

boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = (; :top => boundary_condition_convergence_test,
                        :rest => boundary_condition_convergence_test)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                    source_terms = source_terms,
                                    boundary_conditions = boundary_conditions)

tspan = (0.0, 0.1)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
export LinearCongruentialGenerator

@doc raw"""
    LinearCongruentialGenerator(a, c, m, seed)

A linear congruential generator.

A linear congruential generator (LCG) is a primitive psueudo-random number
generator (PRNG), producing samples that approximate a ``\text{Uniform}(0, 1)``
distribution. Starting from a seed ``X_0``, a sequence is generated according
to the recurence relation,

```math
X_{n+1} = (aX_n + c) \mod m
```

which are then turned into samples ``Z_n = \frac{X_n}{m} \in [0, 1)``.

Although the generator is favoured for esoteric purposes due to its simplicity
and speed, in practice it is highly flawed and can result in samples that have
minimal random structure when viewed through an appropriate projection (as in
[Entacher](https://dl.acm.org/doi/10.1145/272991.273009)). The generator is
also highly sensitive to the choice of ``a``, ``c``, and ``m``.

# Arguments
* `a::Integer`: the "multiplier"; must be between 0 and m, exclusive.
* `c::Integer`: the "increment"; must be between 0 and m-1, inclusive.
* `m::Integer`: the "modulus"; must be non-negative
    proposal distribution.
* `seed::Integer`: the seed value; must be between 0 and m-1 inclusive

# Examples
```julia
# Generating 10 uniform random variables
# Parameters taken from "Numerical Recipes"
a = 1664525
c = 1013904223
m = 2^32
seed = 1729
s = LinearCongruentialGenerator(a, c, m, seed)
samples = sample(s, 10)
```

```julia
# Batch sampling
s = LinearCongruentialGenerator(a, c, m, seed)
samples1 = sample(s, 10)
samples2 = sample(s, 10)
# samples1 ≠ samples2 in general
```
"""
struct LinearCongruentialGenerator <: Sampler
    a::Integer
    c::Integer
    m::Integer
    seed::Integer
    state::Vector{Int}
    function LinearCongruentialGenerator(
        a::Integer,
        c::Integer,
        m::Integer,
        seed::Integer
    )
        # Validation
        if !(0 < m)
            error("m must be positive")
        elseif !(0 < a < m)
            error("a must be between 1 and m, exclusive")
        elseif !(0 <= c < m)
            error("c must be between 0 and m-1, inclusive")
        elseif !(0 ≤ seed < m)
            error("seed must be between 0 and m-1, inclusive")
        end
        state = [seed]
        new(a, c, m, seed, state)
    end
end

@doc raw"""
    sample(s::LinearCongruentialGenerator, N)

Generate samples using a linear congruential generator.

# Arguments
* `s::LinearCongruentialGenerator`: a linear congruential generator.
* `N::Integer`: the number of samples to generate.

# Returns
* `samples::Array{Float64}`: a 2-D array in which each column is a sample.
"""
function sample(s::LinearCongruentialGenerator, N::Integer)
    samples = Array{Float64}(undef, 1, N)
    for i in 1:N
        s.state[1] = (s.a * s.state[1] + s.c) % s.m
        samples[i] = s.state[1] / s.m
    end

    return samples
end


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

name = "libint_compiler"
version = v"2.7.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/evaleev/libint.git", "3bf3a07b58650fe2ed4cd3dc6517d741562e1249")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libint/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-shared=yes --enable-static=no --enable-1body=1 --enable-eri=1
make -j`nproc` export
cp libint-2.7.0-beta.3.tgz ../
cd ..
tar -xzvf libint-2.7.0-beta.3.tgz
cd libint-2.7.0-beta.3
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DLIBINT2_BUILD_SHARED_AND_STATIC_LIBS=ON ..
cmake --build . --target install -- -j${nproc}
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc)
]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"))
    Dependency(PackageSpec(name="Eigen_jll", uuid="bc6bbf8a-a594-5541-9c57-10b0d0312c70"))
    Dependency(PackageSpec(name="boost_jll", uuid="28df3c45-c428-5900-9ff8-a3135698ca75"))
]

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


text:
# ---
# title: 657. Robot Return to Origin
# id: problem657
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/robot-return-to-origin/description/>
# hidden: true
# ---
#
# There is a robot starting at position (0, 0), the origin, on a 2D plane. Given
# a sequence of its moves, judge if this robot **ends up at (0, 0)** after it
# completes its moves.
#
# The move sequence is represented by a string, and the character moves[i]
# represents its ith move. Valid moves are R (right), L (left), U (up), and D
# (down). If the robot returns to the origin after it finishes all of its moves,
# return true. Otherwise, return false.
#
# **Note** : The way that the robot is "facing" is irrelevant. "R" will always
# make the robot move to the right once, "L" will always make it move left, etc.
# Also, assume that the magnitude of the robot's movement is the same for each
# move.
#
#
#
# **Example 1:**
#
#
#
#     Input: moves = "UD"
#     Output: true
#     **Explanation** : The robot moves up once, and then down once. All moves have the same magnitude, so it ended up at the origin where it started. Therefore, we return true.
#
#
# **Example 2:**
#
#
#
#     Input: moves = "LL"
#     Output: false
#     **Explanation** : The robot moves left twice. It ends up two "moves" to the left of the origin. We return false because it is not at the origin at the end of its moves.
#
#
# **Example 3:**
#
#
#
#     Input: moves = "RRDD"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: moves = "LDRRLRUULR"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= moves.length <= 2 * 104`
#   * `moves` only contains the characters `'U'`, `'D'`, `'L'` and `'R'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
println("Running tests at $(pwd())")

using Test, OpenStreetMapX

const m = OpenStreetMapX.get_map_data("data/reno_east3.osm",use_cache=false);

using Random
Random.seed!(0);
const pointA = point_to_nodes(generate_point_in_bounds(m), m)
const pointB = point_to_nodes(generate_point_in_bounds(m), m)

sr1, shortest_distance1, shortest_time1 = OpenStreetMapX.shortest_route(m, pointA, pointB)
@assert length(sr1) > 0

ENV["PLOTS_TEST"] = "true" #copied from Plots.jl tests
ENV["GKSwstype"] = "100" #copied from Plots.jl tests
ENV["MPLBACKEND"]="agg" # no GUI - copied from PyPlot.jl tests

using OpenStreetMapXPlot
println("OpenStreetMapXPlot is located at $(pathof(OpenStreetMapXPlot))")

import Plots
Plots.gr()
Plots.default(show=false, reuse=true)

const trk = [LLA(m.bounds.min_y,m.bounds.min_x,0.0), LLA(m.bounds.max_y,m.bounds.max_x,0.0)]

@testset "Plots backend " begin
	p = OpenStreetMapXPlot.plotmap(m)
	@test typeof(p) == Plots.Plot{Plots.GRBackend}
	@test addroute!(p,m,sr1;route_color="red") == p
	@test plot_nodes!(p,m,[sr1[1],sr1[end]],start_numbering_from=nothing,fontsize=13,color="pink") == p
        @test plot_nodes_as_symbols!(p, m, sr1, symbols=["*","+","#"],colors=["red","green","blue"],fontsize=13) == p
        @test addroute!(p,m,trk;route_color="red") == p
end

@testset "PyPlot backend" begin
	p=OpenStreetMapXPlot.plotmap(m,use_plain_pyplot=true)
	@test p == :osm_use_pyplot
	@test addroute!(p,m,sr1;route_color="red") == p
	@test plot_nodes!(p,m,[sr1[1],sr1[end]],start_numbering_from=nothing,fontsize=13,color="pink") == p
        @test plot_nodes_as_symbols!(p, m, sr1, symbols=["*","+","#"],colors=["red","green","blue"],fontsize=13) == p
	@test addroute!(p,m,trk;route_color="red") == p
end


text:
"""
    PoissonBinomial(p)

A *Poisson-binomial distribution* describes the number of successes in a sequence of independent trials, wherein each trial has a different success rate.
It is parameterized by a vector `p` (of length ``K``), where ``K`` is the total number of trials and `p[i]` corresponds to the probability of success of the `i`th trial.

```math
P(X = k) = \\sum\\limits_{A\\in F_k} \\prod\\limits_{i\\in A} p[i] \\prod\\limits_{j\\in A^c} (1-p[j]), \\quad \\text{ for } k = 0,1,2,\\ldots,K,
```

where ``F_k`` is the set of all subsets of ``k`` integers that can be selected from ``\\{1,2,3,...,K\\}``.

```julia
PoissonBinomial(p)   # Poisson Binomial distribution with success rate vector p

params(d)            # Get the parameters, i.e. (p,)
succprob(d)          # Get the vector of success rates, i.e. p
failprob(d)          # Get the vector of failure rates, i.e. 1-p
```

External links:

* [Poisson-binomial distribution on Wikipedia](http://en.wikipedia.org/wiki/Poisson_binomial_distribution)

"""
mutable struct PoissonBinomial{T<:Real,P<:AbstractVector{T}} <: DiscreteUnivariateDistribution
    p::P
    pmf::Union{Nothing,Vector{T}} # lazy computation of the probability mass function

    function PoissonBinomial{T}(p::AbstractVector{T}; check_args::Bool=true) where {T <: Real}
        check_args && @check_args(PoissonBinomial, all(x -> zero(x) <= x <= one(x), p))
        return new{T,typeof(p)}(p, nothing)
    end
end

function PoissonBinomial(p::AbstractVector{T}; check_args::Bool=true) where {T<:Real}
    return PoissonBinomial{T}(p; check_args=check_args)
end

function Base.getproperty(d::PoissonBinomial, x::Symbol)
    if x === :pmf
        z = getfield(d, :pmf)
        if z === nothing
            y = poissonbinomial_pdf(d.p)
            isprobvec(y) || error("probability mass function is not normalized")
            setfield!(d, :pmf, y)
            return y
        else
            return z
        end
    else
        return getfield(d, x)
    end
end

@distr_support PoissonBinomial 0 length(d.p)

#### Conversions

function PoissonBinomial(::Type{PoissonBinomial{T}}, p::AbstractVector{S}) where {T, S}
    return PoissonBinomial(AbstractVector{T}(p))
end
function PoissonBinomial(::Type{PoissonBinomial{T}}, d::PoissonBinomial{S}) where {T, S}
    return PoissonBinomial(AbstractVector{T}(d.p), check_args=false)
end

#### Parameters

ntrials(d::PoissonBinomial) = length(d.p)
succprob(d::PoissonBinomial) = d.p
failprob(d::PoissonBinomial{T}) where {T} = one(T) .- d.p

params(d::PoissonBinomial) = (d.p,)
partype(::PoissonBinomial{T}) where {T} = T

#### Properties

mean(d::PoissonBinomial) = sum(succprob(d))
var(d::PoissonBinomial) = sum(p * (1 - p) for p in succprob(d))

function skewness(d::PoissonBinomial{T}) where {T}
    v = zero(T)
    s = zero(T)
    p,  = params(d)
    for i in eachindex(p)
        v += p[i] * (one(T) - p[i])
        s += p[i] * (one(T) - p[i]) * (one(T) - T(2) * p[i])
    end
    return s / sqrt(v) / v
end

function kurtosis(d::PoissonBinomial{T}) where {T}
    v = zero(T)
    s = zero(T)
    p,  = params(d)
    for i in eachindex(p)
        v += p[i] * (one(T) - p[i])
        s += p[i] * (one(T) - p[i]) * (one(T) - T(6) * (one(T) - p[i]) * p[i])
    end
    s / v / v
end

entropy(d::PoissonBinomial) = entropy(d.pmf)
median(d::PoissonBinomial) = median(Categorical(d.pmf)) - 1
mode(d::PoissonBinomial) = argmax(d.pmf) - 1
modes(d::PoissonBinomial) = modes(DiscreteNonParametric(support(d), d.pmf))

#### Evaluation

quantile(d::PoissonBinomial, x::Float64) = quantile(Categorical(d.pmf), x) - 1

function mgf(d::PoissonBinomial, t::Real)
    expm1_t = expm1(t)
    mapreduce(*, succprob(d)) do p
        1 + p * expm1_t
    end
end

function cf(d::PoissonBinomial, t::Real)
    cis_t = cis(t)
    mapreduce(*, succprob(d)) do p
        1 - p + p * cis_t
    end
end

pdf(d::PoissonBinomial, k::Real) = insupport(d, k) ? d.pmf[Int(k+1)] : zero(eltype(d.pmf))
logpdf(d::PoissonBinomial, k::Real) = log(pdf(d, k))

cdf(d::PoissonBinomial, k::Int) = integerunitrange_cdf(d, k)

# leads to numerically more accurate results
for f in (:ccdf, :logcdf, :logccdf)
    @eval begin
        $f(d::PoissonBinomial, k::Real) = $(Symbol(f, :_int))(d, k)
        $f(d::PoissonBinomial, k::Int) = $(Symbol(:integerunitrange_, f))(d, k)
    end
end

# Computes the pdf of a poisson-binomial random variable using
# simple, fast recursive formula
#
#      Marlin A. Thomas & Audrey E. Taub (1982)
#      Calculating binomial probabilities when the trial probabilities are unequal,
#      Journal of Statistical Computation and Simulation, 14:2, 125-131, DOI: 10.1080/00949658208810534
#
function poissonbinomial_pdf(p)
    S = zeros(eltype(p), length(p) + 1)
    S[1] = 1
    @inbounds for (col, p_col) in enumerate(p)
        q_col = 1 - p_col
        for row in col:(-1):1
            S[row + 1] = q_col * S[row + 1] + p_col * S[row]
        end
        S[1] *= q_col
    end
    return S
end

# Computes the pdf of a poisson-binomial random variable using
# fast fourier transform
#
#     Hong, Y. (2013).
#     On computing the distribution function for the Poisson binomial
#     distribution. Computational Statistics and Data Analysis, 59, 41–51.
#
function poissonbinomial_pdf_fft(p::AbstractArray{T}) where {T <: Real}
    n = length(p)
    ω = 2 * one(T) / (n + 1)

    x = Vector{Complex{T}}(undef, n+1)
    lmax = ceil(Int, n/2)
    x[1] = one(T)/(n + 1)
    for l=1:lmax
        logz = zero(T)
        argz = zero(T)
        for j=1:n
            zjl = 1 - p[j] + p[j] * cospi(ω*l) + im * p[j] * sinpi(ω * l)
            logz += log(abs(zjl))
            argz += atan(imag(zjl), real(zjl))
        end
        dl = exp(logz)
        x[l + 1] = dl * cos(argz) / (n + 1) + dl * sin(argz) * im / (n + 1)
        if n + 1 - l > l
            x[n + 1 - l + 1] = conj(x[l + 1])
        end
    end
    [max(0, real(xi)) for xi in _dft(x)]
end

# A simple implementation of a DFT to avoid introducing a dependency
# on an external FFT package just for this one distribution
function _dft(x::Vector{T}) where T
    n = length(x)
    y = zeros(complex(float(T)), n)
    @inbounds for j = 0:n-1, k = 0:n-1
        y[k+1] += x[j+1] * cis(-π * float(T)(2 * mod(j * k, n)) / n)
    end
    return y
end

#### Sampling

sampler(d::PoissonBinomial) = PoissBinAliasSampler(d)

## ChainRules definitions

# Compute matrix of partial derivatives [∂P(X=j-1)/∂pᵢ]_{i=1,…,n; j=1,…,n+1}
#
# This implementation uses the same dynamic programming "trick" as for the computation of
# the primals.
#
# Reference (for the primal):
#
#      Marlin A. Thomas & Audrey E. Taub (1982)
#      Calculating binomial probabilities when the trial probabilities are unequal,
#      Journal of Statistical Computation and Simulation, 14:2, 125-131, DOI: 10.1080/00949658208810534
function poissonbinomial_pdf_partialderivatives(p::AbstractVector{<:Real})
    n = length(p)
    A = zeros(eltype(p), n, n + 1)
    @inbounds for j in 1:n
        A[j, end] = 1
    end
    @inbounds for (i, pi) in enumerate(p)
        qi = 1 - pi
        for k in (n - i + 1):n
            kp1 = k + 1
            for j in 1:(i - 1)
                A[j, k] = pi * A[j, k] + qi * A[j, kp1]
            end
            for j in (i+1):n
                A[j, k] = pi * A[j, k] + qi * A[j, kp1]
            end
        end
        for j in 1:(i-1)
            A[j, end] *= pi
        end
        for j in (i+1):n
            A[j, end] *= pi
        end
    end
    @inbounds for j in 1:n, i in 1:n
        A[i, j] -= A[i, j+1]
    end
    return A
end

for f in (:poissonbinomial_pdf, :poissonbinomial_pdf_fft)
    pullback = Symbol(f, :_pullback)
    @eval begin
        function ChainRulesCore.frule(
            (_, Δp)::Tuple{<:Any,<:AbstractVector{<:Real}}, ::typeof($f), p::AbstractVector{<:Real}
        )
            y = $f(p)
            A = poissonbinomial_pdf_partialderivatives(p)
            return y, A' * Δp
        end
        function ChainRulesCore.rrule(::typeof($f), p::AbstractVector{<:Real})
            y = $f(p)
            A = poissonbinomial_pdf_partialderivatives(p)
            function $pullback(Δy)
                p̄ = ChainRulesCore.InplaceableThunk(
                    Δ -> LinearAlgebra.mul!(Δ, A, Δy, true, true),
                    ChainRulesCore.@thunk(A * Δy),
                )
                return ChainRulesCore.NoTangent(), p̄
            end
            return y, $pullback
        end
    end
end


text:
# some Intervals functions that make my life easier

containsfirst(i::Interval) = first(inclusivity(i))
containslast(i::Interval) = last(inclusivity(i))

function Base.convert(::Type{Interval{T}}, i::Interval{N}) where {T,N}
    ifirst = convert(T, first(i))
    ilast = convert(T, last(i))
    return Interval(ifirst, ilast, containsfirst(i), containslast(i) )
end

Base.promote_rule(::Type{Interval{T}}, ::Type{Interval{S}}) where {T,S} = Interval{promote_type(T,S)}

function Base.intersect(a::AbstractInterval{T}, b::AbstractInterval{S}) where {T,S}
    a,b = promote(a,b)
    return intersect(a,b)
end


text:
# -------- #
# Families #
# -------- #

abstract BasisFamily
immutable Cheb <: BasisFamily end
immutable Lin <: BasisFamily end
immutable Spline <: BasisFamily end

# ---------- #
# Parameters #
# ---------- #

abstract BasisParams

type ChebParams <: BasisParams
    n::Int
    a::Float64
    b::Float64
end

function Base.show(io::IO, p::ChebParams)
    m = string("Chebyshev interpoland parameters with ",
               "$(p.n) basis functions from $(p.a), $(p.b)")
    print(io, m)
end

type SplineParams <: BasisParams
    breaks::Vector{Float64}
    evennum::Int
    k::Int
end

# constructor to take a, b, n and form linspace for breaks
SplineParams(n::Int, a::Real, b::Real, k::Int=3) =
    SplineParams(collect(linspace(a, b, n)), 0, k)

function Base.show(io::IO, p::SplineParams)
    m = string("$(p.k) order spline interpoland parameters from ",
               "$(p.breaks[1]), $(p.breaks[end])")
    print(io, m)
end

type LinParams <: BasisParams
    breaks::Vector{Float64}
    evennum::Int
end

# constructor to take a, b, n and form linspace for breaks
LinParams(n::Int, a::Real, b::Real) =
    LinParams(collect(linspace(a, b, n)), 0)

function Base.show(io::IO, p::LinParams)
    m = string("Piecewise linear interpoland parameters ",
               "from $(p.breaks[1]), $(p.breaks[end])")
    print(io, m)
end

for (T, TP) in [(Cheb, ChebParams), (Lin, LinParams), (Spline, SplineParams)]
    @eval _param(::$(T)) = $TP
    @eval _param(::Type{$(T)}) = $TP
end

# params of same type are equal if all their fields are equal
=={T<:BasisParams}(p1::T, p2::T) =
    all(map(nm->getfield(p1, nm) == getfield(p2, nm), fieldnames(p1)))::Bool

# Bases of different dimension can't be equal
=={T1<:BasisParams,T2<:BasisParams}(::T1, ::T2) = false

# ---------- #
# Basis Type #
# ---------- #

type Basis{N, BF<:BasisFamily, BP<:BasisParams}
    basistype::Vector{BF}  # Basis family
    n::Vector{Int}         # number of points and/or basis functions
    a::Vector{Float64}     # lower bound of domain
    b::Vector{Float64}     # upper bound of domain
    params::Vector{BP}     # params to construct basis
end

function Base.show{N}(io::IO, b::Basis{N})
    m = """
    $N dimensional Basis on the hypercube formed by $(b.a) × $(b.b).
    Basis families are $(join(map(x->string(typeof(x)), b.basistype), " × "))
    """
    print(io, m)
end

# constructor that takes all arguments and ensures each has N elemnets
function Basis{BF<:BasisFamily, BP<:BasisParams}(basistype::Vector{BF},
                                                 n::Vector{Int},
                                                 a::Vector{Float64},
                                                 b::Vector{Float64},
                                                 params::Vector{BP})
    N = length(basistype)
    @assert all(map(length, Any[n, a, b, params]) .== N)
    Basis{N,BF,BP}(basistype, n, a, b, params)
end

# univariate basis. Helper to wrap all args in arrays
function Basis{BF<:BasisFamily, BP<:BasisParams}(bt::BF, n::Int, a::Float64,
                                                 b::Float64, params::BP)
    Basis{1,BF,BP}([bt], [n], [a], [b], [params])
end

# constructor to allow `Basis(Spline, breaks, evennum, k)` instead of just
# `Basis(Spline(), breaks, evennum, k)`
Basis{BT<:BasisFamily}(::Type{BT}, args...) = Basis(BT(), args...)

# combining basis -- fundef-esque method
function Basis(b1::Basis, bs::Basis...)
    Nb = length(bs)
    basistype = vcat(b1.basistype, [bs[i].basistype for i=1:Nb]...)
    n = vcat(b1.n, [bs[i].n for i=1:Nb]...)
    a = vcat(b1.a, [bs[i].a for i=1:Nb]...)
    b = vcat(b1.b, [bs[i].b for i=1:Nb]...)
    params = vcat(b1.params, [bs[i].params for i=1:Nb]...)
    N = length(n)

    # determine if BF and BP type are concrete or use fall-back abstracts
    BF1, BP1 = map(x->typeof(x[1]), Any[b1.basistype, b1.params])
    BF = all(map(x->isa(x, BF1), basistype[2:end])) ? BF1 : BasisFamily
    BP = all(map(x->isa(x, BP1), params[2:end])) ? BP1 : BasisParams
    Basis{N,BF,BP}(basistype, n, a, b, params)::Basis{N,BF,BP}
end

# construct basis out of multiple Params (type assertion for stability)
Basis(p::BasisParams, ps::BasisParams...) =
    Basis(Basis(p), Basis[Basis(p) for p in ps]...)::Basis{length(ps) + 1}

# fundefn type method
Basis{T<:BasisFamily}(bt::T, n::Vector, a::Vector, b::Vector) =
    Basis(map(_param(T), n, a, b)...)::Basis{length(n)}

# special method for Spline that adds `k` argument
Basis(::Spline, n::Vector, a::Vector, b::Vector, k::Vector=ones(Int, length(n))) =
    Basis(map(SplineParams, n, a, b, k)...)::Basis{length(n)}

# ----------------- #
# Basis API methods #
# ----------------- #

# separating Basis -- just re construct it from the nth set of params
function Base.getindex{N}(basis::Basis{N}, n::Int)
    n < 0 || n > N && error("n must be between 1 and $N")
    Basis(basis.params[n])::Basis{1}
end

# AbstractArray like methods for Basis
Base.ndims{N}(::Basis{N}) = N
Base.length(b::Basis) = prod(b.n)
Base.size(b::Basis, i::Int) = length(b[i])  # uses method on previous line
Base.size{N}(b::Basis{N}) = tuple(b.n...)::NTuple{N, Int64}

# Bases of different dimension can't be equal
=={N,M}(::Basis{N}, ::Basis{M}) = false

# basis are equal if all fields of the basis are equal
=={N}(b1::Basis{N}, b2::Basis{N}) =
    all(map(nm->getfield(b1, nm) == getfield(b2, nm), fieldnames(b1)))::Bool

function nodes(b::Basis{1})
    x = nodes(b.params[1])
    (x, Vector{Float64}[x])
end

function nodes{N}(b::Basis{N})  # funnode method
    xcoord = Vector{Float64}[nodes(b.params[j]) for j in 1:N]
    x = gridmake(xcoord...)::Matrix{Float64}
    return x, xcoord
end


text:
module util

# Write your package code here.
struct Cleaner
pattern::Regex
exclude::Union{Array{Regex,1},Missing}
dir::String
clean_symlinks::bool
force::bool
Cleaner(pattern,exclude=missing,dir=pwd(),clean_symlinks=false,force=true) = new(pattern,exclude,dir,clean_symlinks,force)
end

function (c::Cleaner)()
    dir = c.dir
    if !isdirpath(dir)
        return throw(ArgumentError("$dir is not a valid directory"))
    end
    for (root, dirs, files) = walkdir(dir,follow_symlinks=follow_symlinks)
        for file in files
            name = joinpath(root,file)
            m = match(c.pattern,name)
            if !isnothing(m)
                kill = true
                if !ismissing(c.exclude)
                    for k in c.exclude
                    if !isnothing(match(k,name))
                        kill = false
                        break
                    end
                end
                end
                if kill == true
                rm(name)
                end
            end
        end

    end
end
export Cleaner
end


text:
for func in [:clblasiSamax, :clblasiDamax, :clblasiCamax, :clblasiZamax]
    @eval @blasfun $func(N::Csize_t,
                              iMax::cl.CL_mem, offiMax::Csize_t,
                              X::cl.CL_mem, offx::Csize_t, incx::Cint,
                              scratch_buff::cl.CL_mem,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(N::Csize_t,
                               iMax::cl.CL_mem, offiMax::Csize_t,
                               X::cl.CL_mem, offx::Csize_t, incx::Cint,
                               scratch_buff::cl.CL_mem)
end


text:
using Test
using Forecast
using Random
using Logging

const tests = [
    "acf",
    "ar",
    "arsim",
    "boxcox",
    "ccf",
    "d",
    "datasets",
    "forecast_ar",
    "hma",
    "loess",
    "p",
    "pacf",
    "plots",
    "sma",
    "stl",
    "summarize",
    "utils",
    "utils_datetime",
]

printstyled("\nTest Summary List:\n", color=:underline)

Random.seed!(36)
Base.disable_logging(Base.CoreLogging.Error) # disable warnings

for t in tests
    @testset "Test $t" begin
        Random.seed!(36)
        include("$t.jl")
        println()
    end
end

println()


text:
function solve()
    n = parse(Int, readline())
    as = zeros(Int, n)
    bs = zeros(Int, n)
    for i=1:n
        a, b = [parse(Int, x) for x in split(readline())]
        as[i] = a
        bs[i] = b
    end
    sort!(as)
    sort!(bs)
    if n % 2 == 0
        m = n ÷ 2
        ma = (as[m] + as[m+1])
        mb = (bs[m] + bs[m+1])
        mb - ma + 1
    else
        m = (n+1) ÷ 2
        ma = as[m]
        mb = bs[m]
        mb - ma + 1
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
## main function
"""
    run_ed(descriptor, tasks...; runmode)
    run_ed(descriptor, tasks, runmode)

Main entry to the exact diagonalization functionality.

# Parameters
 - `descriptor`: A [`EDDataDescriptor`](@ref). This struct describes what should be diagonalized and how.
 - `tasks`: A vector of [`EDTask`](@ref)s. What to compute from the diagonalized Hamiltonians.
 - `runmode`: How to parallelize the task. See [`RunMode`](@ref).

# Returns
A vector of [`EDDerivedData`](@ref)s each corresponding to a task. Note that saving the data is on the caller!
"""
run_ed(descriptor, tasks...; runmode) = run_ed(descriptor, tasks, runmode)

function run_ed(descriptor, tasks, runmode)
    diagtype = descriptor.diagtype
    setup(runmode)
    array_type = array_initializer(runmode)

    model = initialize_model!(descriptor.model, array_type) # -> _compute_interactions
    task_array_init = (type, sizes...) -> array_type(type, sizes..., ed_size(diagtype), parameter_dims(model)...)
    initialize!.(tasks, task_array_init, spectral_size(diagtype, model), )
    parameter_chunks = split_workload(model, num_workers(runmode))
    start_computation(runmode, parameter_chunks) do chunk
        # call _compute_core! for each chunk
        _compute_core!(tasks, diagtype, model, chunk)
    end
    return assemble.(tasks, Ref(descriptor))
end

function _compute_core!(tasks, diagtype, model, parameter_chunk)
    tasks = initialize_local.(tasks)
    do_parameters(model, parameter_chunk) do parameter_index, H
        try
            diagtime = 0.0
            tasktime = 0.0
            start = time()
            diagonalize!(H, diagtype) do diag_index, eigvals, eigvecs
                diagtime += time() - start
                start = time()
                compute_task!.(tasks, Ref(eigvals), Ref(eigvecs), Ref(diag_index), Ref(parameter_index))
                tasktime += time() - start
                start = time()
            end
            logmsg("Diagonalization took TDIAG=$(round(diagtime; digits=2))s and tasks took TTASK=$(round(tasktime; digits=2))s.")
        catch e;
            logmsg("Error during diagonalization occured for index $parameter_index : $e")
            display(stacktrace(catch_backtrace()))
            rethrow(e)
            failed_task!.(tasks, :, Ref(parameter_index))
        end
    end
end

# function _compute_parallel!(tasks, edd, posdata)
#     N = edd.system_size
#     nshots = edd.shots
#     nρs = length(edd.ρs)
#     for task in tasks
#         initialize!(task, edd, _sharedarray_constructor)
#     end

#     interactions = SharedArray{Float64}((N,N,nshots,nρs))
#     _compute_interactions!(interactions, edd, posdata)

#     workloads = _workload_chunks(nshots*nρs, length(procs(interactions)))

#     @sync for (i, p) in enumerate(procs(interactions))
#         @async remotecall_wait(_compute_parallel_job!, p,
#             tasks, interactions, workloads[i], edd)
#     end
# end

# function _compute_parallel_job!(tasks, interactions, workload, desc)
#     # build operators and initial state here
#     # -> is it worth to just to this once and use SharedArrays?
#     # probably not much of a difference as the opjects here are quite small in RAM
#     spin_ops = real.(symmetrize_operator.(op_list(σx/2, desc.system_size), Ref(desc.basis)))
#     field_operator = sum(spin_ops)

#     logmsg("Range: $(workload) on #$(indexpids(interactions))")
#     _compute_core!(tasks, interactions, workload, desc.fields, field_operator, desc.basis, desc.diagtype)
# end

# function _compute_threaded!(tasks, edd, posdata)
#     N = edd.system_size
#     nshots = edd.shots
#     nρs = length(edd.ρs)
#     for task in tasks
#         initialize!(task, edd, _array_constructor)
#     end

#     spin_ops = real.(symmetrize_operator.(op_list(σx/2, N), Ref(edd.basis)))
#     field_operator = sum(spin_ops)

#     interactions = FArray{4}(undef, (N,N,nshots,nρs))
#     _compute_interactions!(interactions, edd, posdata)

#     # use one thread less as the main thread should not be used for work (I think)
#     workers = max(1, Threads.nthreads()-1)
#     workloads = _workload_chunks(nshots*nρs, workers)
#     @sync for i in 1:workers
#         Threads.@spawn _compute_core!(tasks, interactions, workloads[i], edd.fields, field_operator, edd.basis, edd.diagtype)
#     end
# end

# ## TODO better design? less copy-paste?
# function _compute_core!(tasks, interactions, workload, field_values, field_operator, symmetry, diagtype::Full)
#     tasks = initialize_local.(tasks)
#     nshots = size(interactions, 3)
#     matrix = zeros(eltype(field_operator), size(field_operator)) # preallocate, but DENSE
#     for index in workload
#         i, shot = _flat_to_indices(index, nshots)
#         J = @view interactions[:,:, shot, i]
#         model = real.(symmetrize_operator(xxzmodel(J, -0.73), symmetry))

#         for (k, h) in enumerate(field_values)
#             copyto!(matrix, model + h*field_operator) # this also converts from sparse to dense!
#             try
#                 evals, evecs = eigen!(Hermitian(matrix))
#                 compute_task!.(tasks, i, shot, k, Ref(evals), Ref(evecs))
#             catch e;
#                 logmsg("Error occured for #field=$k shot=$shot #rho=$i: $e")
#                 display(stacktrace(catch_backtrace()))
#                 failed_task!.(tasks, i, shot, k)
#                 continue
#             end
#         end
#         logmsg(@sprintf("Done %03i - #rho =%2i - %03i/%03i", index, i, shot, nshots))
#     end
# end

# function _compute_core!(tasks, interactions, workload, field_values, field_operator, symmetry, diagtype::Sparse)
#     tasks = initialize_local.(tasks)
#     nshots = size(interactions, 3)
#     nev = diagtype.count
#     sigma = diagtype.σ
#     for index in workload
#         i, shot = _flat_to_indices(index, nshots)
#         J = @view interactions[:,:, shot, i]
#         model = real.(symmetrize_operator(xxzmodel(J, -0.73), symmetry))

#         for (k, h) in enumerate(field_values)
#             try
#                 evals, evecs, nconv, niter, _, _ = eigs(Hermitian(model + h*field_operator); nev, sigma, which=:LM, tol=0, check=2)
#                 evals, evecs = LinearAlgebra.sorteig!(evals, evecs)
#                 logmsg("converged=",nconv,"/",nev," | niter=",niter)
#                 compute_task!.(tasks, i, shot, k, Ref(evals), Ref(evecs))
#             catch e;
#                 logmsg("Error occured for #field=$k shot=$shot #rho=$i: $e")
#                 display(stacktrace(catch_backtrace()))
#                 failed_task!.(tasks, i, shot, k)
#                 continue
#             end
#         end
#         logmsg(@sprintf("Done %03i - #rho =%2i - %03i/%03i", index, i, shot, nshots))
#     end
# end

# run_ed(desc::EDDataDescriptor, tasks::EDTask...) = run_ed(desc, tasks)

# function run_ed(desc::EDDataDescriptor, tasks::Vector{EDTask})
#     posdata = load_positions(desc.geometry, desc.dimension, desc.system_size, desc.pathdata)
#     return run_ed(desc, posdata, tasks)
# end

# run_ed(desc::EDDataDescriptor, posdata::PositionData, tasks::EDTask...) = run_ed(desc, posdata, tasks)

# function run_ed(desc::EDDataDescriptor, posdata::PositionData, tasks::Vector{<:EDTask})
#     # decide on mode of operation
#     logmsg("ToDo: rho=$(desc.ρs)")
#     logmsg("with $(desc.shots) realizations and  $(length(desc.fields)) field values")
#     wcount = length(workers())
#     if wcount > 1
#         logmsg("Running ED on $(wcount) PROCESSES")
#         _compute_parallel!(tasks, desc, posdata)
#     else
#         logmsg("Running ED on $(Threads.nthreads()) THREADS")
#         _compute_threaded!(tasks, desc, posdata)
#     end
#     return assemble.(tasks, Ref(desc))
# end


text:
module AdS5_3_1

using Jecco
using LinearAlgebra
using SparseArrays: SparseMatrixCSC
using OrdinaryDiffEq

import Base.Threads.@threads
import Base.Threads.@spawn

# abstract types and structs used throughout
include("types.jl")

include("system.jl")

include("potential.jl")

include("initial_data.jl")

# Kreiss-Oliger dissipation and filtering
include("dissipation.jl")

# time marching orders
include("setup_rhs.jl")

# nested system
include("equation_inner_coeff.jl")
include("equation_outer_coeff.jl")
include("inner_to_outer.jl")
include("solve_nested.jl")

# evolution equations
include("compute_boundary_t.jl")
include("compute_gauge_t.jl")
include("compute_bulkevolved_t.jl")

# finding the Apparent Horizon
include("equations_AH.jl")
include("find_AH.jl")

# input/output
include("IO.jl")

include("diagnostics.jl")

# run the model
include("run.jl")

# expressions for the Vacuum Expectation Values
include("VEVs.jl")

# random utilities useful for post-processing
include("utils.jl")
include("design_new_data.jl")
include("recipes.jl")

# always set the number of BLAS threads to 1 upon loading the module. by default
# it uses a bunch of them and we don't want that since they trample over each
# other when solving the nested systems equations. it's much better to thread
# over the loop. see also the discussion here:
# https://github.com/JuliaLang/julia/issues/33409
#
# this saves us the tedious task of always setting OMP_NUM_THREADS=1 before
# launching julia.
function __init__()
    LinearAlgebra.BLAS.set_num_threads(1)
    nothing
end

export SpecCartGrid3D
export Potential, ZeroPotential, Phi8Potential
export BlackBrane, BlackBranePert, BlackBrane_xi1, PhiGaussian_u
export Advect_xi, ConstantAH, AHF
export AffineNull, EvolTest0
export NoDiag, DiagAH
export BulkEvolved, BulkConstrained, Boundary, Gauge, Bulk
export System, SystemPartition
export BulkEvolvedPartition, BulkConstrainedPartition, BulkDerivPartition
export Integration, InOut, run_model
export get_energy, get_Jx, get_Jy, get_px, get_py, get_pz, get_pxy, get_Ophi
export BoundaryTimeSeries, XiTimeSeries, BulkTimeSeries, VEVTimeSeries, ConstrainedTimeSeries
export convert_to_mathematica

end


text:
"""
This type convert a bilevel problem as an MPCCModel.

Bilevel Program:
````
min_{x,y} f1(x,y)
s.t.  lcon    <= c(x)   <= ucon
      lvar    <= x      <= uvar
      y in argmin_{y} f2(x,y)
                 s.t. ll_lcon <= g(x,y) <= ll_ucon (l1,  l2)
                      ll_lvar <= y      <= ll_uvar (mu1, mu2)
````
MPCC reformulation of the Bilevel Program:
       min_{x,y,l1,l2,mu1,mu2} f1(x,y)
       s.t.  lvar <= x      <= uvar
             lcon <= c(x)   <= ucon
             0    <= ∇_y f2(x,y) - ∇_y g(x,y)' * (l1 - l2) - (mu1 - mu2) <= 0
             G(x):=(  g(x,y)  - ll_lcon)
                      ll_ucon - g(x,y) )
                      ll_uvar - y
                      y - ll_lvar )
             H(x):=( l1, l2, mu1, mu2 )
             0 <= G(x) _|_ H(x) >= 0
````
Gradient, Jacobian and Hessian are computed using automatic differentiation

set no_third_derivative as true, to ignore the 3rd derivative of g.
"""
mutable struct BPMPCCModel <: AbstractMPCCModel
  meta     :: MPCCModelMeta
  counters :: MPCCCounters

  # Functions
  f1 :: Function
  f2 :: Function
  c  :: Function
  g  :: Function

  ll_lcon :: AbstractVector
  ll_ucon :: AbstractVector
  ll_lvar :: AbstractVector
  ll_uvar :: AbstractVector

  no_third_derivative :: Bool

  m      :: Int #size of the lower level dimension
  n      :: Int #size of the upper level dimension + m
  llncon :: Int #number of constraints in the lower level problem
end

#m0 is the initial guess for the lower leverl problem
function BPMPCCModel(f1, x0::AbstractVector, m0::AbstractVector;
                     y0::AbstractVector = eltype(x0)[],
                     lvar::AbstractVector = eltype(x0)[],
                     uvar::AbstractVector = eltype(x0)[],
                     lcon::AbstractVector = eltype(x0)[],
                     ucon::AbstractVector = eltype(x0)[],
                     c = (args...)->throw(NotImplementedError("c")),
                     g = (args...)->throw(NotImplementedError("g")),
                     f2 = (args...)->throw(NotImplementedError("f2")),
                     ll_lcon::AbstractVector = eltype(x0)[],
                     ll_ucon::AbstractVector = eltype(x0)[],
                     ll_lvar::AbstractVector = eltype(x0)[],
                     ll_uvar::AbstractVector = eltype(x0)[],
                     name::String = "Generic",
                     lin::AbstractVector{<: Integer}=Int[],
                     no_third_derivative::Bool=false)

  m, n, llncon = length(ll_lvar), length(x0), length(ll_ucon)

  nvar = n + 3 * m + 2 * llncon
  length(lvar) == 0 && (lvar = -Inf*ones(n))
  length(uvar) == 0 && (uvar =  Inf*ones(n))
  lvar = vcat(lvar, -Inf*ones(3 * m + 2 * llncon))
  uvar = vcat(uvar,  Inf*ones(3 * m + 2 * llncon))

  ncon = maximum([length(lcon); length(ucon); length(y0)]) + length(m0)
  lcon = vcat(lcon, zeros(m))
  ucon = vcat(ucon, zeros(m))

  ncc  = 2 * m + 2 * llncon
  lccG = zeros(ncc)
  lccH = zeros(ncc)
  ini0 = vcat(x0, m0, zeros(2 * m + 2 * llncon))

  meta = MPCCModelMeta(nvar, x0 = ini0, lvar = lvar, uvar = uvar,
                                        ncon = ncon, lcon = lcon,
                                        ucon = ucon, ncc = ncc,
                                        lccG = lccG, lccH = lccH)

  return BPMPCCModel(meta, MPCCCounters(), f1, f2, c, g,
                     ll_lcon, ll_ucon, ll_lvar, ll_uvar,
                     no_third_derivative, m, n, llncon)
end

"""
_var_bp(:: BPMPCCModel, :: AbstractVector):
return the details of the different variables in the MPCC corresponding
to the bilevel program.
return: x, y, l1, l2, m1, m2
"""
function _var_bp(nlp :: BPMPCCModel, x :: AbstractVector)

 if length(x) != nlp.meta.nvar throw("Error: wrong vector size") end

 m, n, ncon = nlp.m, nlp.n, nlp.llncon

 x0, y0 = x[1:n], x[n+1:n+m]
 l1, l2 = x[n+m+1:n+m+ncon], x[n+m+ncon+1:n+m+2*ncon]
 m1, m2 = x[n+m+2*ncon+1:n+2*m+2*ncon], x[n+2*m+2*ncon+1:n+3*m+2*ncon]

 return x0, y0, l1, l2, m1, m2
end

"""
JLag: gives the gradient w.r.t. y of the Lagrangian of the lower level problem.

JLag(:: BPMPCCModel, :: Function, :: Function, :: AbstractVector)
"∇_y f2(x,y) - ∇_y g(x,y)' * v - mu"
"""
function JLag(nlp :: BPMPCCModel,
              f   :: Function,
              c   :: Function,
              x   :: AbstractVector)

    Jtv = Array{Union{Nothing, Number}}(nothing, nlp.m)

    x0, y0, l1, l2, m1, m2 = _var_bp(nlp, x)

    Jtv = ForwardDiff.gradient(f, x)[nlp.n+1:nlp.n+nlp.m] - (m1 - m2)
    increment!(nlp, :neval_grad)

    if nlp.llncon > 0
        Jtv += - ForwardDiff.jacobian(c, x)[:,nlp.n+1:nlp.n+nlp.m]' * (l1-l2)
        increment!(nlp, :neval_jac)
    end

    return Jtv
end

"""
HLag: gives the hessian w.r.t. y of the Lagrangian of the lower level problem.

HLag(:: BPMPCCModel, :: Function, :: Function, :: AbstractVector)
"∇^2_y f2(x,y) - sum_i ∇^2_y g_i(x,y)' * v_i"
"""
function HLag(nlp :: BPMPCCModel,
              f   :: Function,
              c   :: Function,
              x   :: AbstractVector)
    Jtv = zeros(nlp.n+nlp.m, nlp.n+nlp.m)
    x0, y0, l1, l2, m1, m2 = _var_bp(nlp, x)
    lag(x) = nlp.llncon > 0 ? f(x) - dot(c(x), l1-l2) : f(x)
    Jtv = ForwardDiff.hessian(lag, x)
    Jtvy = Jtv[nlp.n+1:nlp.n+nlp.m,nlp.n+1:nlp.n+nlp.m]
    Jtvx = Jtv[nlp.n+1:nlp.n+nlp.m,1:nlp.n]
    return Jtvy, Jtvx
end

function obj(nlp :: BPMPCCModel, x :: AbstractVector)
  increment!(nlp, :neval_obj)
  return nlp.f1(x)
end

function grad!(nlp :: BPMPCCModel, x :: AbstractVector, g :: AbstractVector)
  increment!(nlp, :neval_grad)
  ForwardDiff.gradient!(view(g, 1:length(x)), nlp.f1, x)
  return g
end

function cons_nl!(nlp :: BPMPCCModel, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_cons)
  Jl = JLag(nlp, nlp.f1, nlp.g, x)
  c[1:nlp.meta.ncon] .= vcat(nlp.c(x), Jl)
  return c
end

function consG!(nlp :: BPMPCCModel, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_consG)
  gxy = nlp.g(x)
  y   = x[nlp.n+1:nlp.n+nlp.m]
  c[1:nlp.meta.ncc] .= vcat( gxy - nlp.ll_lcon,
                            -gxy + nlp.ll_ucon,
                             y - nlp.ll_lvar,
                            -y + nlp.ll_uvar)
  return c
end

function consH!(nlp :: BPMPCCModel, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_consH)
  c[1:nlp.meta.ncc] .= x[nlp.n+nlp.m+1:nlp.meta.nvar]
  return c
end

function jac_nl_structure!(nlp  :: BPMPCCModel,
                           rows :: AbstractVector{<: Integer},
                           cols :: AbstractVector{<: Integer})
  m, n = nlp.meta.ncon, nlp.meta.nvar
  I = ((i,j) for i = 1:m, j = 1:n)
  rows .= getindex.(I, 1)[:]
  cols .= getindex.(I, 2)[:]
  return rows, cols
end

function jac_nl_coord!(nlp  :: BPMPCCModel,
                       x    :: AbstractVector,
                       vals :: AbstractVector)
  increment!(nlp, :neval_jac)
  m, n, llncon = nlp.m, nlp.n, nlp.llncon

  Jx  = ForwardDiff.jacobian(nlp.c, x)
  Jg  = ForwardDiff.jacobian(nlp.g, x)
  Jyy, Jyx = HLag(nlp, nlp.f2, nlp.g, x)
  Jyl = llncon >0 ? vcat(Jg[:,n+1:n+m], -Jg[:,n+1:n+m]) : zeros(0, m)
  Jym = hcat(diagm(0 => ones(m)),-diagm(0 => ones(m)))
  Jl  = hcat(Jyx, Jyy, Jyl', Jym)
  J   = vcat(Jx, Jl)

  vals .= J[:]
  return vals
end

function jnlprod!(nlp :: BPMPCCModel,
                  x   :: AbstractVector,
                  v   :: AbstractVector,
                  Jv  :: AbstractVector)
  increment!(nlp, :neval_jprod)
  #Jv[1:nlp.meta.ncon] = ForwardDiff.derivative(t -> nlp.c(x + t * v), 0)
  Jv .= jac_nl(nlp, x) * v
  return Jv
end

function jnltprod!(nlp :: BPMPCCModel,
                   x   :: AbstractVector,
                   v   :: AbstractVector,
                   Jtv :: AbstractVector)
  increment!(nlp, :neval_jtprod)
  #Jtv[1:nlp.meta.nvar] = ForwardDiff.gradient(x -> dot(nlp.c(x), v), x)
  Jtv .= jac_nl(nlp, x)' * v
  return Jtv
end

function jacG_structure!(nlp  :: BPMPCCModel,
                         rows :: AbstractVector{<: Integer},
                         cols :: AbstractVector{<: Integer})
  m, n = nlp.meta.ncc, nlp.meta.nvar
  I = ((i,j) for i = 1:m, j = 1:n)
   rows[1 : n*m] .= getindex.(I, 1)[:]
   cols[1 : n*m] .= getindex.(I, 2)[:]
  return rows, cols
end

function jacG_coord!(nlp  :: BPMPCCModel,
                     x    :: AbstractVector,
                     vals :: AbstractVector)
  increment!(nlp, :neval_jacG)
  ncc, nvar = nlp.meta.ncc, nlp.meta.nvar

  Jg  = ForwardDiff.jacobian(nlp.g, x)
  Idy = hcat(zeros(nlp.m,nlp.n),diagm(0 => ones(nlp.m)),zeros(nlp.m,ncc))
  Jx  = vcat(Jg, -Jg, Idy, - Idy)

  vals[1 : nvar * ncc] .= Jx[:]
  return vals
end

function jGprod!(nlp :: BPMPCCModel,
                 x   :: AbstractVector,
                 v   :: AbstractVector,
                 Jv  :: AbstractVector)
  increment!(nlp, :neval_jGprod)
  #Jv[1:nlp.meta.ncc] = ForwardDiff.derivative(t -> nlp.G(x + t * v), 0)
  Jv[1:nlp.meta.ncc] .= jacG(nlp, x) * v
  return Jv
end

function jGtprod!(nlp :: BPMPCCModel,
                  x   :: AbstractVector,
                  v   :: AbstractVector,
                  Jtv :: AbstractVector)
  increment!(nlp, :neval_jGtprod)
  #Jtv[1:nlp.meta.nvar] = ForwardDiff.gradient(x -> dot(nlp.G(x), v), x)
  Jtv[1:nlp.meta.nvar] .= jacG(nlp, x)' * v
  return Jtv
end

function jacH_structure!(nlp  :: BPMPCCModel,
                         rows :: AbstractVector{<: Integer},
                         cols :: AbstractVector{<: Integer})
  m, n = nlp.meta.ncc, nlp.meta.nvar
  I = ((i,j) for i = 1:m, j = 1:n)
  rows[1 : n*m] .= getindex.(I, 1)[:]
  cols[1 : n*m] .= getindex.(I, 2)[:]
  return rows, cols
end

function jacH_coord!(nlp  :: BPMPCCModel,
                     x    :: AbstractVector,
                     vals :: AbstractVector)
  m, n = nlp.meta.ncc, nlp.meta.nvar
  increment!(nlp, :neval_jacH)
  Jx = hcat(zeros(m, nlp.n + nlp.m), diagm(0 => ones(m)))

  vals[1 : n*m] .= Jx[:]
  return vals
end

function jHprod!(nlp :: BPMPCCModel,
                 x   :: AbstractVector,
                 v   :: AbstractVector,
                 Jv  :: AbstractVector)
  increment!(nlp, :neval_jHprod)
  #Jv[1:nlp.meta.ncc] = ForwardDiff.derivative(t -> nlp.G(x + t * v), 0)
  Jv[1:nlp.meta.ncc] .= jacH(nlp, x) * v
  return Jv
end

function jHtprod!(nlp :: BPMPCCModel,
                  x   :: AbstractVector,
                  v   :: AbstractVector,
                  Jtv :: AbstractVector)
  increment!(nlp, :neval_jHtprod)
  #Jtv[1:nlp.meta.nvar] = ForwardDiff.gradient(x -> dot(nlp.G(x), v), x)
  Jtv[1:nlp.meta.nvar] .= jacH(nlp, x)' * v
  return Jtv
end

function hess_structure!(nlp :: BPMPCCModel,
                        rows :: AbstractVector{<: Integer},
                        cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows[1 : nlp.meta.nnzh] .= getindex.(I, 1)
  cols[1 : nlp.meta.nnzh] .= getindex.(I, 2)
  return rows, cols
end

function hess_coord!(nlp        :: BPMPCCModel,
                     x          :: AbstractVector,
                     vals       :: AbstractVector;
                     obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  ℓ(x) = obj_weight * nlp.f1(x)
  Hx = ForwardDiff.hessian(ℓ, x)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hess_coord!(nlp        :: BPMPCCModel,
                     x          :: AbstractVector,
                     y          :: AbstractVector,
                     vals       :: AbstractVector;
                     obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  ncon, ncc, m, nt = nlp.meta.ncon, nlp.meta.ncc, nlp.m, nlp.no_third_derivative
  llncon = nlp.llncon

  temp(x) = nt ? obj_weight * nlp.f1(x) + dot(nlp.c(x),y[1:ncon - m]) : obj_weight * nlp.f1(x) + dot(cons_nl(nlp,x), y[1:ncon])
  ℓ(x) = llncon > 0 ? temp(x) + dot(nlp.g(x), y[ncon+1:ncon+llncon]) - dot(nlp.g(x), y[ncon+llncon+1:ncon+2*llncon]) : temp(x)
  Hx = ForwardDiff.hessian(ℓ, x)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hessG(nlp :: BPMPCCModel, x :: AbstractVector, y :: AbstractVector)
  increment!(nlp, :neval_hessG)
  lln = nlp.llncon
  ℓ(x) = lln > 0 ? dot(nlp.g(x), y[1:lln]) - dot(nlp.g(x), y[lln+1:2*lln]) : 0
  Hx = ForwardDiff.hessian(ℓ, x)
  return tril(Hx)
end

function hessG_structure!(nlp :: BPMPCCModel,
                         rows :: AbstractVector{<: Integer},
                         cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows[1 : nlp.meta.nnzh] .= getindex.(I, 1)
  cols[1 : nlp.meta.nnzh] .= getindex.(I, 2)
  return rows, cols
end

function hessG_coord!(nlp  :: BPMPCCModel,
                      x    :: AbstractVector,
                      y    :: AbstractVector,
                      vals :: AbstractVector)
  increment!(nlp, :neval_hessG)
  lln = nlp.llncon
  ℓ(x) = lln > 0 ? dot(nlp.g(x), y[1:lln]) - dot(nlp.g(x), y[lln+1:2*lln]) : 0
  Hx = ForwardDiff.hessian(ℓ, x)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hessH(nlp :: BPMPCCModel, x :: AbstractVector, y :: AbstractVector)
  increment!(nlp, :neval_hessH)
  Hx = zeros(nlp.meta.nvar, nlp.meta.nvar)
  return tril(Hx)
end

function hessH_structure!(nlp  :: BPMPCCModel,
                          rows :: AbstractVector{<: Integer},
                          cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows .= getindex.(I, 1)
  cols .= getindex.(I, 2)
  return rows, cols
end

function hessH_coord!(nlp  :: BPMPCCModel,
                      x    :: AbstractVector,
                      y    :: AbstractVector,
                      vals :: AbstractVector)
  increment!(nlp, :neval_hess)
  ℓ(x) = dot(nlp.H(x), y)
  Hx = zeros(nlp.meta.nvar, nlp.meta.nvar)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hprod!(nlp        :: BPMPCCModel,
                x          :: AbstractVector,
                v          :: AbstractVector,
                Hv         :: AbstractVector;
                obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  ℓ(x) = obj_weight * nlp.f1(x)
  Hv .= ForwardDiff.derivative(t -> ForwardDiff.gradient(ℓ, x + t * v), 0)
  return Hv
end

function hprod!(nlp        :: BPMPCCModel,
                x          :: AbstractVector,
                y          :: AbstractVector,
                v          :: AbstractVector,
                Hv         :: AbstractVector;
                obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  ncon, ncc, m, nt = nlp.meta.ncon, nlp.meta.ncc, nlp.m, nlp.no_third_derivative
  llncon = nlp.llncon

  temp(x) = nt ? obj_weight * nlp.f1(x) + dot(nlp.c(x),y[1:ncon - m]) : obj_weight * nlp.f1(x) + dot(cons_nl(nlp,x), y[1:ncon])
  ℓ(x) = llncon > 0 ? temp(x) + dot(nlp.g(x), y[ncon+1:ncon+llncon]) - dot(nlp.g(x), y[ncon+llncon+1:ncon+2*llncon]) : temp(x)
  Hv .= ForwardDiff.derivative(t -> ForwardDiff.gradient(ℓ, x + t * v), 0)
  return Hv
end

function hGprod!(nlp :: BPMPCCModel,
                 x   :: AbstractVector,
                 y   :: AbstractVector,
                 v   :: AbstractVector,
                 Hv  :: AbstractVector)
  increment!(nlp, :neval_hGprod)
  lln = nlp.llncon
  ℓ(x) = lln > 0 ? dot(nlp.g(x), y[1:lln]) - dot(nlp.g(x), y[lln+1:2*lln]) : 0
  Hv .= ForwardDiff.derivative(t -> ForwardDiff.gradient(ℓ, x + t * v), 0)
  return Hv
end

function hHprod!(nlp :: BPMPCCModel,
                 x   :: AbstractVector,
                 y   :: AbstractVector,
                 v   :: AbstractVector,
                 Hv  :: AbstractVector)
  increment!(nlp, :neval_hHprod)
  Hv .= zeros(nlp.meta.nvar)
  return Hv
end


text:
using .CUDAKernels, .CUDAKernels.CUDA
import ModelParameters.Flatten

export CuGPU

"""
    CuGPU <: GPU

    CuGPU()
    CuGPU{threads_per_block}()

```julia
ruleset = Ruleset(rule; proc=CuGPU())
# or
output = sim!(output, rule; proc=CuGPU())
```
"""
struct CuGPU{X} <: GPU end
CuGPU() = CuGPU{32}()

# CUDA setup

kernel_setup(::CuGPU{N}) where N = CUDAKernels.CUDADevice(), (N, N)

# _proc_setup
# Convert all arrays in SimData to CuArrays
@noinline function _proc_setup(::CuGPU, simdata::AbstractSimData)
    Adapt.adapt(CuArray, simdata)
end

_copyto_output!(outgrid, grid::GridData, proc::GPU) = copyto!(outgrid, gridview(grid))

# Thread-safe CUDA atomic ops

for (f, op) in atomic_ops
    atomic_f = Symbol(:atomic_, f)
    @eval begin
        function ($f)(d::WritableGridData{<:Any,R}, ::CuGPU, x, I...) where R
            A = parent(dest(d))
            i = Base._to_linear_index(A, (I .+ R)...)
            (CUDA.$atomic_f)(pointer(A, i), x)
        end
    end
end


text:
abstract NonlinearSolver{T}

solve!(s::NonlinearSolver) = error("solve! not implemented for $(typeof(s))")

function solve!{T}(s::NonlinearSolver{T}, x₀::Vector{T})
    setInitialConditions!(s, x₀)
    solve!(s)
end

# default parameters for nonlinear solvers
const DEFAULT_atol=2eps()
const DEFAULT_rtol=2eps()
const DEFAULT_stol=2eps()
const DEFAULT_nmax=100
const DEFAULT_ϵ=sqrt(eps())

abstract NonlinearFunctionParameters{T}

function function_stages!{T}(y::Vector{T}, b::Vector{T}, params::NonlinearFunctionParameters{T})
    error("No function_stages implemented for this integrator.")
end

immutable NonlinearSolverParameters{T}
    nmax::Int   # maximum number of iterations

    atol::T     # absolute tolerance
    rtol::T     # relative tolerance
    stol::T     # successive tolerance

    atol²::T
    rtol²::T
    stol²::T

    function NonlinearSolverParameters(nmax, atol, rtol, stol)
        @assert nmax > 0
        @assert atol > 0
        @assert rtol > 0
        @assert stol > 0

        new(nmax, atol, rtol, stol, atol^2, rtol^2, stol^2)
    end
end

type NonlinearSolverStatus{T}
    i::Int      # iteration number
    r₀::T       # initial residual (absolute)
    rₐ::T       # residual (absolute)
    rᵣ::T       # residual (relative)
    rₛ::T       # residual (successive)

    NonlinearSolverStatus() = new(0, 0., 0., 0., 0.)
end

Base.show(io::IO, status::NonlinearSolverStatus) = print(io,
                        status.i, ", ", status.rₐ,", ",
                        status.rᵣ,", ", status.rₛ)

function solverStatusOK(status::NonlinearSolverStatus, params::NonlinearSolverParameters)
    return (status.rₐ < params.atol² ||
            status.rᵣ < params.rtol  ||
            status.rₛ < params.stol²)&&
           status.i  ≤ params.nmax
end

function getLinearSolver(T, n, linear_solver)
    if linear_solver == nothing
        linear_solver = LUSolverLAPACK(zeros(T, n, n), zeros(T, n))
    else
        @assert typeof(linear_solver) <: LinearSolver{T}
        @assert n == linear_solver.n
    end
    return linear_solver
end

function residual_absolute{T}(x::Vector{T})
    local r::T = 0.
    for xᵢ in x
        r = max(r, xᵢ*xᵢ)
    end
    r
end

function residual_relative{T}(δx::Vector{T}, x::Vector{T})
    @assert length(x) == length(δx)
    local r::T = 0.
    for i in 1:length(x)
        r = max(r, abs(δx[i] / x[i]))
    end
    r
end


text:
using Walker2014Sampling
using Documenter

DocMeta.setdocmeta!(Walker2014Sampling, :DocTestSetup, :(using Walker2014Sampling); recursive=true)

makedocs(;
    modules=[Walker2014Sampling],
    authors="Iván Gutiérrez <ivangutierrez1988@gmail.com> and contributors",
    repo="https://github.com/igutierrezm/Walker2014Sampling.jl/blob/{commit}{path}#{line}",
    sitename="Walker2014Sampling.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://igutierrezm.github.io/Walker2014Sampling.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/igutierrezm/Walker2014Sampling.jl",
)


text:
@testset "display_obj.jl" begin
    ds = DS.DisplayStack()

    panel = DS.Panel(" 我是標題 ", [10, 30], [3, 5])
    label1 = DS.Label("我是字串1", [5, 10])
    push!(ds, :panel=>panel)
    push!(ds, :label=>label1)
    render(ds)
    @test T.read_stream(T.out_stream) ==
        "\e[34;1m\e[3;5H" *
        "\e[s┌ 我是標題 ──────────────────┐\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s└────────────────────────────┘\e[u\e[1B" *
        "\e[0m" *
        "\e[32m\e[5;10H我是字串1\e[0m" *
        "\e[24;1H"
    @test ds[:panel] == panel
    @test ds[:label] == label1

    label2 = DS.Label("我是字串2", [6, 11])
    push!(ds, label2)
    render(ds)
    @test T.read_stream(T.out_stream) ==
        "\e[34;1m\e[3;5H" *
        "\e[s┌ 我是標題 ──────────────────┐\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s│                            │\e[u\e[1B" *
        "\e[s└────────────────────────────┘\e[u\e[1B" *
        "\e[0m" *
        "\e[32m\e[5;10H我是字串1\e[0m" *
        "\e[32m\e[6;11H我是字串2\e[0m" *
        "\e[24;1H"
end


text:
module similarity_scorer

export createScorer, computeSimilarity, Scorer, updateIndex!, MlScorer, PositionPosting, BasicScorer, ProposedScorer

using util
using aggregation_strategy
using Method
using DataStructures

abstract type Scorer end

struct PositionPosting
    report::Report
    positions::Vector{Int16}
end

struct Posting
    report::Report
    freq::Float64
end

@inline function get_dfcache(method, doc_freq::DocFreq)
    if doc_freq.nDocs != doc_freq.previous_n_docs
        doc_freq.cache = computeDocFreq(method, doc_freq)
        doc_freq.previous_n_docs = doc_freq.nDocs
    end

    return doc_freq.cache
end

#################################################
# Basic
#################################################

struct BasicScorer{M <: MethodType} <: Scorer
    aggType::AggregationType
    method::M
end

function updateIndex!(scorer::BasicScorer, query::Report, doc_freq::DocFreq)
    nothing
end

function computeSimilarity(scorer::BasicScorer, query::Report, candidates::Vector{Report}, docFreq::DocFreq)
    results = fill(-9999999.99, length(candidates))
    df_vec = get_dfcache(scorer.method, docFreq)

    if length(query.stacks) == 0
        @info "Query $(query.id) is empty"
        return results
    end

    Threads.@threads for report_idx = 1:length(candidates)
        candidate = candidates[report_idx]

        if length(candidate.stacks) != 0
            agg = Aggregation(scorer.aggType, length(query.stacks), length(candidate.stacks))

            for (q_stack_idx, q_stack) in enumerate(query.stacks)
                for (c_stack_idx, c_stack) in enumerate(candidate.stacks)
                    sim = similarity(scorer.method, q_stack, c_stack, df_vec)
                    update_score_matrix!(agg, q_stack_idx, c_stack_idx, sim)
                end
            end

            results[report_idx] = aggregate(agg)
        end
    end

    return results
end

#################################################
# TraceSim Scorer - Reduce the algorithm runtime
#################################################

struct TraceSimScorer <: Scorer
    aggType::AggregationType
    method::TraceSim
end

function updateIndex!(scorer::TraceSimScorer, query::Report, doc_freq::DocFreq)
    nothing
end

"""
For each stack, create a vector that the indices of the sorted vector.
"""
function generate_sorted_indices(report::Report)
    if length(report.stats.sorted_idxs) > 0
        return report.stats.sorted_idxs
    end

    vec = report.stats.sorted_idxs

    for stack in report.stacks
        push!(vec, sortperm(stack))
    end

    return vec
end

function computeSimilarity(scorer::TraceSimScorer, query::Report, candidates::Vector{Report}, docFreq::DocFreq)
    results = fill(-9999999.99, length(candidates))
    df_vec = get_dfcache(scorer.method, docFreq)

    if length(query.stacks) == 0
        @info "Query $(query.id) is empty"
        return results
    end

    q_weights  = Vector{Vector{Float64}}()

    for q_stack in query.stacks
        push!(q_weights, compute_frame_weights(scorer.method, q_stack, df_vec))
    end

    q_sorted_idxs = generate_sorted_indices(query)

    Threads.@threads for report_idx = 1:length(candidates)
        candidate = candidates[report_idx]

        if length(candidate.stacks) != 0
            agg = Aggregation(scorer.aggType, length(query.stacks), length(candidate.stacks))
            c_sorted_idxs  = generate_sorted_indices(candidate)

            for (q_stack_idx, q_stack) in enumerate(query.stacks)
                for (c_stack_idx, c_stack) in enumerate(candidate.stacks)
                    sim = similarity(scorer.method, q_stack, c_stack, q_weights[q_stack_idx], q_sorted_idxs[q_stack_idx], c_sorted_idxs[c_stack_idx], df_vec)
                    update_score_matrix!(agg, q_stack_idx, c_stack_idx, sim)
                end
            end

            results[report_idx] = aggregate(agg)
        end
    end

    return results
end

#################################################
# Quick Align without Inverted Index
#################################################

struct FaSTScorer <: Scorer
    method::FaST
    aggType::AggregationType
end

function get_data_structure(scorer::FaSTScorer, report::Report)
    return aggregatePositionByFuncKeepDim(report)
end

function updateIndex!(scorer::FaSTScorer, new_report::Report, doc_freq::DocFreq)
    new_report.v_positional_bow = get_data_structure(scorer, new_report)
end

function computeSimilarity(scorer::FaSTScorer, query::Report, candidates::Vector{Report}, docFreq::DocFreq)
    similarities = fill(-9999999.99, length(candidates))
    df_vec = get_dfcache(scorer.method, docFreq)

    if length(query.stacks) == 0
        @info "Query $(query.id) is empty"
        return similarities
    end

    q_bow = get_data_structure(scorer, query)

    Threads.@threads for report_idx = 1:length(candidates)
        candidate = candidates[report_idx]

        if length(candidate.stacks) != 0
            agg = Aggregation(scorer.aggType, length(q_bow), length(candidate.v_positional_bow))

            for (q_stack_idx, q_pos_bow) in enumerate(q_bow)
                for (c_stack_idx, c_pos_bow) in enumerate(candidate.v_positional_bow)
                    sim = similarity(scorer.method, q_pos_bow, c_pos_bow, df_vec)
                    update_score_matrix!(agg, q_stack_idx, c_stack_idx, sim)
                end
            end

            similarities[report_idx] = aggregate(agg)
        end
    end

    return similarities
end

#################################################
# Information Retrieval
#################################################

struct PlainIRScorer{M <: MethodType} <: Scorer
    method::M
    aggType::AggregationType
    ngram::Int32
end

function transform_BOW_dim(method, report::Report, ngram, vocab)::Vector{Vector{Pair{UInt32, Float64}}}
    v_bow = Vector{Vector{Pair{UInt32, Float64}}}(undef, length(report.stacks))

    for (idx, stack) in enumerate(report.stacks)
        bow = Dict{UInt32, Float64}()

        for (pos, func_id) in enumerate(stack)
            if haskey(bow, func_id)
                bow[func_id] += 1.0
            else
                bow[func_id] = 1.0
            end
        end

        for ngram_tkn in generate_ngrams(report.stacks, ngram)
            ngram_id = vocab.vocab[ngram_tkn]
            if haskey(bow, ngram_id)
                bow[ngram_id] += 1.0
            else
                bow[ngram_id] = 1.0
            end
        end

        v_bow[idx] = sort!(collect(bow))
    end

    return v_bow
end

function generate_bow_vector_dim!(scorer::PlainIRScorer, report::Report, vocab::Vocab, df_vec::Vector{Float64})
    bow = transform_BOW_dim(scorer.method, report, scorer.ngram, vocab)
    report.v_bow = bow
    report.stats.doc_lens = [length(stack) for stack in report.stacks ]

    return report.v_bow
end

function updateIndex!(scorer::PlainIRScorer, new_report::Report, doc_freq::DocFreq)
    generate_bow_vector_dim!(scorer, new_report, doc_freq.vocab, get_dfcache(scorer.method, doc_freq))
end

function computeSimilarity(scorer::PlainIRScorer, query::Report, candidates::Vector{Report}, docFreq::DocFreq)
    q_v_bow = generate_bow_vector_dim!(scorer, query, docFreq.vocab, get_dfcache(scorer.method, docFreq))
    similarities = fill(-9999999.99, length(candidates))
    df_vec = get_dfcache(scorer.method, docFreq)

    if length(query.stacks) == 0
        @info "Query $(query.id) is empty"
        return similarities
    end

    Threads.@threads for report_idx = 1:length(candidates)
        candidate = candidates[report_idx]

        if length(candidate.stacks) != 0
            agg = Aggregation(scorer.aggType, length(q_v_bow), length(candidate.v_bow))
            for (q_stack_idx, q_bow) in enumerate(q_v_bow)
                for (c_stack_idx, c_bow) in enumerate(candidate.v_bow)
                    sim = similarity(scorer.method, q_bow, c_bow, df_vec, candidate.stats.doc_lens[c_stack_idx])
                    update_score_matrix!(agg, q_stack_idx, c_stack_idx, sim,)
                end
            end

            similarities[report_idx] = aggregate(agg)
        end
    end

    return similarities
end

#################################################
# Scorer that used our method to compute the similarity
#################################################

mutable struct ProposedScorer{M <: MethodType} <: Scorer
    method::M
    index::Dict{UInt32,Vector{PositionPosting}}
    is_sum_weights::Bool
end

function compute_weight_sum!( new_report::Report, scorer::ProposedScorer, function2positions::Vector{Pair{UInt32, Vector{Int16}}}, df_vec::Vector{Float64})
    new_report.stats.weight_sum != -Inf && return

    weight_sum = 0.0

    for (func_id, positions) in function2positions
        @inbounds df = df_vec[func_id]

        for pos in positions
            w = frame_weight(scorer.method, pos, df)
            weight_sum += w
        end
    end

    new_report.stats.weight_sum = weight_sum
end

function updateIndex!(scorer::ProposedScorer, new_report::Report, doc_freq::DocFreq)
    index = scorer.index
    df_vec = get_dfcache(scorer.method, doc_freq)

    function2positions = aggregatePositionByFunction(new_report)

    for (func_id, positions) in function2positions
        postingList = haskey(index, func_id) ? index[func_id] : Vector{PositionPosting}()

        if length(postingList) == 0
            index[func_id] = postingList
        end

        push!(postingList, PositionPosting(new_report, positions))
    end

    if scorer.is_sum_weights
        compute_weight_sum!(new_report, scorer, function2positions, df_vec)
    end

end

function computeSimilarity(scorer::ProposedScorer, query::Report, candidates::Vector{Report}, docFreq::DocFreq)
    index = scorer.index
    method = scorer.method
    q_bow = aggregatePositionByFunction(query)
    df_vec =  get_dfcache(scorer.method, docFreq)

    candidate2idx = Dict{Int64,Int64}()
    sizehint!(candidate2idx, length(candidates))

    for (cand_idx, cand) in  enumerate(candidates)
        candidate2idx[cand.id] = cand_idx
    end

    if scorer.is_sum_weights
        compute_weight_sum!(query, scorer, q_bow, df_vec)
    end

    stub = Vector{Posting}()
    similarities = zeros(length(candidates))

    for token_idx = 1:length(q_bow)
        func_id, q_positions = q_bow[token_idx]

        postingList = get(index, func_id, stub)
        ndf = df_vec[func_id]

        length(postingList) == 0 && continue

        Threads.@threads for posting in postingList
            candidate_idx = get(candidate2idx, posting.report.id, -1)
            candidate_idx < 0 && continue

            c_positions = posting.positions
            score = 0.0
            j = 1
            min_len_positions = length(q_positions) < length(c_positions) ? length(q_positions) + 1 : length(c_positions) + 1

            @inbounds while j != min_len_positions
                c_gap_score = gap_score(method, c_positions[j], ndf)
                q_gap_score = gap_score(method, q_positions[j], ndf)

                score -= (q_gap_score + c_gap_score)
                m = match_score(method, q_positions[j], c_positions[j], ndf, -q_gap_score, -c_gap_score)
                score += m
                # println("#Cand:$(candidate_idx)\tfunc: $(func_id)\tc_pos: $(c_positions[j]), q_pos: $(q_positions[j]), ndf: $(ndf)")
                # println("===> c_gap_score: $(c_gap_score), q_gap_score: $(q_gap_score), m: $(m)")

                j += 1
            end

            similarities[candidate_idx] += score
        end
    end

    return compute_final_similarity!(similarities, method, query, candidates, df_vec)
end

#########################################################
# Information retrieval techniques
########################################################
mutable struct IRScorer{M <: IRMethod} <: Scorer
    method::M
    index::Dict{UInt32,Vector{Posting}}
    doc_sum::Float64
    n_docs::Float64
    ngram::Int32
end

function get_posting_list!(index, func_id)
    postingList = haskey(index, func_id) ? index[func_id] : Vector{Posting}()

    if length(postingList) == 0
        index[func_id] = postingList
    end

    return postingList
end

function transform_BOW(method, report::Report, ngram, vocab)::Vector{Pair{UInt32, Float64}}
    bow = Dict{UInt32, Float64}()

    for stack in report.stacks
        for (pos, func_id) in enumerate(stack)
            if haskey(bow, func_id)
                bow[func_id] += 1.0
            else
                bow[func_id] = 1.0
            end
        end
    end

    for ngram_tkn in generate_ngrams(report.stacks, ngram)
        ngram_id = vocab.vocab[ngram_tkn]
        if haskey(bow, ngram_id)
            bow[ngram_id] += 1.0
        else
            bow[ngram_id] = 1.0
        end
    end

    return collect(bow)
end

function generate_bow_vector!(scorer::IRScorer, report::Report, vocab::Vocab, df_vec::Vector{Float64})
    bow = transform_BOW(scorer.method, report, scorer.ngram, vocab)
    report.bow = normalize_vec(scorer.method, bow, df_vec)

    return report.bow
end

function updateIndex!(scorer::IRScorer, new_report::Report, doc_freq::DocFreq)
    index = scorer.index
    weighted_bow = generate_bow_vector!(scorer, new_report, doc_freq.vocab, get_dfcache(scorer.method, doc_freq))
    new_report.bow = weighted_bow

    for idx in eachindex(weighted_bow)
        func_id, term_weight = weighted_bow[idx]
        postingList = get_posting_list!(index, func_id)

        push!(postingList, Posting(new_report, term_weight))
    end

    doc_len = set_doc_len!(new_report)
    scorer.doc_sum += doc_len
    scorer.n_docs += 1
end

function computeSimilarity(scorer::IRScorer, query::Report, candidates::Vector{Report}, docFreq::DocFreq)
    index = scorer.index
    method = scorer.method
    q_bow = generate_bow_vector!(scorer, query, docFreq.vocab, get_dfcache(scorer.method, docFreq))

    candidate2idx = Dict{Int64,Int64}(map(k -> (k[2].id, k[1]), enumerate(candidates)))
    idf_vec =  get_dfcache(scorer.method, docFreq)
    doc_avg = scorer.doc_sum/scorer.n_docs

    # number_of_match = zeros(length(candidates))
    # lengths = [ cand.stats.doc_len for cand in candidates]

    stub = Vector{Posting}()
    similarities = zeros(length(candidates))

    for token_idx = 1:length(q_bow)
    # for token_idx = 1:length(q_bow)
        func_id, q_freq = q_bow[token_idx]

        posting_list = get(index, func_id, stub)
        idf = idf_vec[func_id]

        length(posting_list) == 0 && continue

        Threads.@threads for posting in posting_list
            candidate_idx = get(candidate2idx, posting.report.id, -1)
            candidate_idx < 0 && continue

            # println("#Cand:$(posting.report.id)\tfunc: $(func_id)\t q_freq=$(q_freq)\tc_freq=$(posting.freq)")
            # println("\t$(compute_term_sim(method, idf, q_freq, posting.freq, posting.report.stats.doc_len , doc_avg))")
            similarities[candidate_idx] += compute_term_sim(method, idf, q_freq, posting.freq, posting.report.stats.doc_len , doc_avg)
            # number_of_match[candidate_idx]+=1.0
        end
    end

    # mean_match = sum([n/l for (n,l) in zip(number_of_match,lengths)])/ length(candidates)
    # println("Mean $(mean_match)")

    return compute_final_similarity!(similarities, method, query, candidates, idf_vec)
end

#########################################################
# Create Scorer function
########################################################

@inline function createScorer(methodName::String, args, upt_doc_freq)::Scorer
    if methodName == "trace_sim"
        @info("TraceSim")
        aggStrategy = parseAggStrategy(args[:aggregate])
        method = TraceSim(args[:alpha], args[:beta], args[:gamma],
            get(args, :sigmoid, false),
            get(args, :b, 0.0),
            get(args, :idf, false),
            get(args, :reciprocal_func, true),
            get(args, :no_norm, false),
            get(args, :sum_match, false),
            )

        return TraceSimScorer(aggStrategy, method)
    elseif methodName == "pdm"
        @info("PDM")

        aggStrategy = parseAggStrategy(args[:aggregate])
        method = PDM(args[:c], args[:o])

        return BasicScorer{PDM}(aggStrategy, method)
    elseif methodName == "prefix_match"
        @info("PrefixMatch")
        aggStrategy = parseAggStrategy(args[:aggregate])
        method = PrefixMatch()

        return BasicScorer{PrefixMatch}(aggStrategy, method)
    elseif methodName == "fast"
        @info("FaST")

        method = FaST(args[:alpha], args[:beta], args[:gamma], upt_doc_freq)
        return ProposedScorer{FaST}(method, Dict{UInt32,MutableLinkedList{PositionPosting}}(), !upt_doc_freq)
    elseif methodName == "plain_fast"
        @info("FaST Without Inverted Index")
        aggStrategy = parseAggStrategy(args[:aggregate])
        method = FaST(args[:alpha], args[:beta], args[:gamma], upt_doc_freq)

        return FaSTScorer(method, aggStrategy)
    elseif methodName == "plain_tf_idf"
        @info("TF-IDF Without Inverted Index")
        aggStrategy = parseAggStrategy(args[:aggregate])

        method = TFIDF(get(args, :tf_scheme, "sqrt"), get(args, :idf_scheme, "smooth"),  false, get(args, :add_query_freq, false))

        return PlainIRScorer{TFIDF}(method, aggStrategy, get(args, :ngram, 1))
    elseif methodName == "plain_cosine_tf_idf"
        @info("Cosine Without Inverted Index")
        aggStrategy = parseAggStrategy(args[:aggregate])
        method = CosineTFIDF(get(args, :tf_scheme, "sqrt"), get(args, :idf_scheme, "smooth"),  !upt_doc_freq)

        return PlainIRScorer{CosineTFIDF}(method, aggStrategy, get(args, :ngram, 1))
    elseif methodName == "tf_idf"
        method = TFIDF(get(args, :tf_scheme, "sqrt"), get(args, :idf_scheme, "smooth"),  false, get(args, :add_query_freq, false))
        return IRScorer{TFIDF}(method, Dict{UInt32,Vector{Posting}}(), 0, 0, get(args, :ngram, 1))
    elseif methodName == "bm25"
        method = BM25(arg[:k1], args[:b], 0.0)
        return IRScorer{BM25}(method, Dict{UInt32,Vector{Posting}}(), 0, 0, get(args, :ngram, 1))
    elseif methodName == "cosine_tf_idf"
        @info("CosineTFIDF")
        method = CosineTFIDF(get(args, :tf_scheme, "sqrt"), get(args, :idf_scheme, "smooth"), !upt_doc_freq)
        return IRScorer{CosineTFIDF}(method, Dict{UInt32,Vector{Posting}}(), 0, 0, get(args, :ngram, 1))
    elseif methodName == "fast_pdm"
        @info("Fast PDM")
        method = PDM(args[:c], args[:o])
        return ProposedScorer{PDM}(method, Dict{UInt32,MutableLinkedList{PositionPosting}}(), false)
    elseif methodName == "fast_trace_sim"
        @info("Fast TraceSim")
        method = TraceSim(args[:alpha], args[:beta], args[:gamma],
            get(args, :sigmoid, false),
            get(args, :b, 0.0),
            get(args, :idf, false),
            get(args, :reciprocal_func, true),
            get(args, :no_norm, false)
            )

        return ProposedScorer{TraceSim}(method, Dict{UInt32,MutableLinkedList{PositionPosting}}(), false)
    end
end

end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
"""
   magsq!(magusq::Nodes{Primal/Dual},u::Edges{Primal/Dual})

Calculate the in-placed squared magnitude of vector grid data `u`, placing the result on
the cell centers in `magusq`.
"""
function magsq!(umagsq::Nodes{C},u::Edges{C}) where {C <: CellType}

  usq = u∘u
  usq_nodes = Nodes(C,u)

  grid_interpolate!(usq_nodes,usq.u)
  umagsq .= usq_nodes

  grid_interpolate!(usq_nodes,usq.v)
  umagsq .+= usq_nodes

  return umagsq
end

"""
   magsq(u::Edges{Primal/Dual}) -> Nodes{Primal/Dual}

Calculate the squared magnitude of vector grid data `u`, placing the result on
the cell centers.
"""
magsq(u::Edges{C}) where {C<:CellType} = magsq!(Nodes(C,u),u)

"""
    mag!(magu::Nodes{Primal/Dual},u::Edges{Primal/Dual})

Calculate the in-placed magnitude of vector grid data `u`, placing the result on
the cell centers in `magu`.
"""
function mag!(magu::Nodes{C},u::Edges{C}) where {C <: CellType}
    magsq!(magu,u)
    @. magu = sqrt(magu)
end

"""
   mag(u::Edges{Primal/Dual}) -> Nodes{Primal/Dual}

Calculate the magnitude of vector grid data `u`, placing the result on
the cell centers.
"""
mag(u::Edges{C}) where {C<:CellType} = mag!(Nodes(C,u),u)

"""
    directional_derivative!(out,f,q)

Compute the directional derivative of `f` in the direction of `q`, ``q\\cdot\\nabla f``,
and put the result into `out`. Note that the result is not scaled by any grid spacing.
"""
function directional_derivative!(out::Edges{C},f::Edges{C},q::Edges{C}) where {C <: CellType}

    qdf = EdgeGradient(C,f)
    df = EdgeGradient(C,f)
    qt = EdgeGradient(C,f)

    out .= 0.0
    product!(qdf,transpose(grid_interpolate!(qt,q)),grad!(df,f))
    grid_interpolate!(out,qdf)

    #=
    # some temp storage
    q_dual = Nodes(othertype(C),f)
    q_primal = Nodes(C,f)
    outx = XEdges(C,f)
    outy = YEdges(C,f)

    # compute the gradient of f
    df = EdgeGradient(C,f)
    grad!(df,f)

    # compute qy*(dfx/dy) after first transforming q.v to dual nodes and then transforming product to primal x-edges
    grid_interpolate!(outx,grid_interpolate!(q_dual,q.v) ∘ df.dudy)
    out.u .= outx

    # compute qx*(dfx/dx) after first transforming q.u to primal nodes
    # and then interpolating product to primal x-edges
    grid_interpolate!(outx,grid_interpolate!(q_primal,q.u) ∘ df.dudx)
    out.u .+= outx

    # compute qx*(dfy/dx) after first transforming qx to dual nodes and then transforming product to primal y-edges
    grid_interpolate!(outy,grid_interpolate!(q_dual,q.u) ∘ df.dvdx)
    out.v .= outy

    # compute qy*(dfy/dy) after first transforming qy to primal nodes
    # and then interpolating product to primal y-edges
    grid_interpolate!(outy,grid_interpolate!(q_primal,q.v) ∘ df.dvdy)
    out.v .+= outy
    =#

    return out

end

"""
    directional_derivative_conserve!(out,f,q)

Compute the conservative form of the directional derivative of `f` in the direction
of `q`, ``\\nabla\\cdot(qf)``, and put the result
into `out`. Note that the result is not scaled by any grid spacing. This form is only
appropriate if `q` is divergence free.
"""
function directional_derivative_conserve!(out::Edges{C},f::Edges{C},q::Edges{C}) where {C <: CellType}

    # some temp storage
    q_dual = Nodes(othertype(C),f)
    f_dual = Nodes(othertype(C),f)

    qf = q ∘ f
    qfnode = NodePair(C,f)

    # interpolate both fx and qx to primal nodes, multiply element by element, assign to NodePair.u
    grid_interpolate!(qfnode.u,qf.u)

    # interpolate both f.u and q.v to dual nodes, multiply element by element, assign to NodePair.v
    qfnode.v .= grid_interpolate!(q_dual,q.v) .* grid_interpolate!(f_dual,f.u)

    # take divergence and assign to out.u
    divergence!(out.u,qfnode)

    qfnode = NodePair(othertype(C),f)

    # interpolate both fy and qx to dual nodes, multiply element by element, assign to NodePair.u
    qfnode.u .= grid_interpolate!(q_dual,q.u) .* grid_interpolate!(f_dual,f.v)

    # interpolate both fy and qy to primal nodes, multiply element by element, assign to NodePair.v
    grid_interpolate!(qfnode.v,qf.v)

    # take divergence and assign to out.v
    divergence!(out.v,qfnode)

    return out

end

"""
    curl_cross!(out,a,b)

Compute the curl of the cross product of `a` and `b`,
``\\nabla\\times(a\\times b)``,
and put the result into `out`. Note that the result is not scaled by any grid spacing.
"""
function curl_cross!(out::Edges{C},a::Edges{C},b::Edges{C}) where {C <: CellType}

    # some temp storage
    bx_dual = Nodes(othertype(C),a)
    by_dual = Nodes(othertype(C),a)

    ax_dual = Nodes(othertype(C),a)
    ay_dual = Nodes(othertype(C),a)

    acrossb = Nodes(othertype(C),a)

    grid_interpolate!((ax_dual,ay_dual),a)
    grid_interpolate!((bx_dual,by_dual),b)

    acrossb .= ax_dual.*by_dual - ay_dual.*bx_dual
    curl!(out,acrossb)

    return out

end

"""
    convective_derivative!(out,q)

Compute the convective derivative of `q` in the form
``u\\cdot\\nabla u``
and put the result into `out`. Note that the result is not scaled by any grid spacing.
"""
convective_derivative!(out::Edges{C},q::Edges{C}) where {C<:CellType} = directional_derivative!(out,q,q)

"""
    convective_derivative_rot!(out,q)

Compute the rotational form of the convective derivative of `q` in the form
``\\frac{1}{2}\\nabla|u|^2-u\\times(\\nabla\\times u)``
and put the result into `out`. Note that the result is not scaled by any grid spacing.
"""
function convective_derivative_rot!(out::Edges{C},q::Edges{C}) where {C<:CellType}

    # some temp storage
    qx_primal = Nodes(C,q)
    qy_primal = Nodes(C,q)
    ω = curl(q)

    grid_interpolate!((qx_primal,qy_primal),q)
    grad!(out,qx_primal ∘ qx_primal + qy_primal ∘ qy_primal)
    out .*= 0.5

    q_dual = Nodes(othertype(C),q)
    qcrossω = Edges(C,q)
    grid_interpolate!(qcrossω.u,grid_interpolate!(q_dual, q.v) ∘ ω)
    grid_interpolate!(qcrossω.v,grid_interpolate!(q_dual,-q.u) ∘ ω)

    out .-= qcrossω

    return out

end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
@testset "704.binary-search.jl" begin
    @test search704([-1,0,3,5,9,12], 9) == 5
    @test search704([-1,0,3,5,9,12], 2) == 0
end


text:
using Libdl

check_llvm() = isone(length(filter(lib->occursin("LLVM\b", basename(lib)), Libdl.dllist())))
@static if Sys.ARCH === :x86_64 || Sys.ARCH === :i686
    if Base.libllvm_version >= v"8" && check_llvm()
        include("build_x86.jl")
    else
        include("build_x86_cpuidonly.jl")
    end
else
    include("build_generic.jl")
end

open(joinpath(@__DIR__, "..", "src", "cpu_info.jl"), "w") do f
    write(f, cpu_info_string)
end


text:
using MultivariateStats,LinearAlgebra,Flux.Data.MNIST,Measures,Plots
pyplot()

imgs, labels   = MNIST.images(), MNIST.labels()
x = hcat([vcat(float.(im)...) for im in imgs]...)
pca = fit(PCA, x; maxoutdim=2)
M = projection(pca)

function compareDigits(dA,dB)
    imA, imB = imgs[labels .== dA], imgs[labels .== dB]
    xA = hcat([vcat(float.(im)...) for im in imA]...)
    xB = hcat([vcat(float.(im)...) for im in imB]...)
    zA, zB = M'*xA, M'*xB
    default(ms=0.8, msw=0, xlims=(-5,12.5), ylims=(-7.5,7.5),
            legend = :topright, xlabel="PC 1", ylabel="PC 2")
    scatter(zA[1,:],zA[2,:], c=:red,  label="Digit $(dA)")
    scatter!(zB[1,:],zB[2,:], c=:blue, label="Digit $(dB)")
end

plots = []
for k in 1:5
    push!(plots,compareDigits(2k-2,2k-1))
end
plot(plots...,size = (800, 500), margin = 5mm)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
if(Pkg.installed("PyPlot") == nothing)
    Pkg.add("PyPlot")
end

using PyPlot
close("all")

function plotfunction(PDATA::PAINTER.PAINTER_Data,OIDATA::PAINTER.PAINTER_Input)
    nx = OIDATA.nx
    nw = OIDATA.nw
    wvl = OIDATA.wvl
    FOV = OIDATA.FOV
    x = PDATA.x
    w = PDATA.w .> 0.

    indpix = linspace(-(FOV / 2), (FOV / 2), nx)
    pos = round.(Int, [1, round(Int,nx / 4), round(Int,nx / 2), round(Int,nx * 3 / 4), nx])

    count_y = 0
    count_x = 0
    SubRow  = 6
    SubColumn = 5

    if nw == 30

    for n in 1:nw
        subplot(SubColumn, SubRow, n)
        imshow(x[:, :, n] .* max.(0, w[:, :, n]), origin = "lower")
        titlestring = @sprintf("%2.4f µm", wvl[n] * 1e6)
        title(titlestring)
        xticks([])
        yticks([])
        if( n == (nw + 1 - SubRow + count_x) )
            xticks(collect(pos - 1), round.(Int,indpix[pos] * 100000) / 100)
            xlabel("FOV (mas)")
            count_x += 1
        end
        if(n == (1 + count_y * SubRow))
            yticks(collect(pos - 1), round.(Int,indpix[pos] * 100000) / 100)
            ylabel("FOV (mas)")
            count_y += 1
        end
    end

    else
        NW = ceil(Int, nw/30)
        m = 0
        for n in 1:NW:nw
            m+=1
            subplot(SubColumn, SubRow, m)
            imshow(x[:, :, n] .* max(0, w[:, :, n]), origin = "lower")
            titlestring = @sprintf("%2.4f µm", wvl[n] * 1e6)
            title(titlestring)
            xticks([])
            yticks([])
        if( m == (round(Int,nw/NW)  + 1 - SubRow + count_x) )
            xticks(collect(pos - 1), round(Int,indpix[pos] * 100000) / 100)
            xlabel("FOV (mas)")
            count_x += 1
        end
        if( m == (1 + count_y * SubRow))
            yticks(collect(pos - 1), round(Int,indpix[pos] * 100000) / 100)
            ylabel("FOV (mas)")
            count_y += 1
        end
        end
    end

end

    PlotFct = plotfunction
    FOV = 0.01
    indwvl = 1:6:180

    nx = 64

    eps1 = 1e-4
    eps2 = 1e-4

    admm = true
    aff = true     # plot is enabled

    nbitermax = 1
    savepath = "data.jld"
    dptype = "all"
    dpprm = 5
    Folder = ""
    # first estimate is flat normalized (ones everywhere)
    xinit3D = PAINTER.mask(nx, round(Int, nx/2 ) )
    # initialize algorithm and run admm

    rho_y = 10.
    rho_spat = 4.
    rho_ps = rho_spat
    rho_spec = .5
    lambda_L1 = 5*1e-3
    nbitermax = 1

    OIDATA, PDATA = PAINTER.painter(nbitermax = nbitermax, nx = nx, rho_y = rho_y,
                            rho_spat = rho_spat, lambda_L1 = lambda_L1,
                            rho_spec = rho_spec, rho_ps = rho_ps,
                            eps1 = eps1, eps2 = eps2, FOV = FOV, indwvl = indwvl, admm=admm,
                            PlotFct = PlotFct, aff = aff, dptype = dptype,
                            dpprm = dpprm, Folder = Folder, xinit3D = xinit3D)

# save data struture in .jld files
    println("save data")
    PAINTER.paintersave(savepath,PDATA,OIDATA)

# load data struture in .jld files
    println("load data")
    PDATA, OIDATA = PAINTER.painterload(savepath)

# Warm start of the algorithm to Iteration 100
    nbitermax = 100
    println("PAINTER: Warm Start")
    OIDATA, PDATA = PAINTER.painter(OIDATA,PDATA,nbitermax,true, PlotFct = PlotFct)

    nothing


text:
module Isaac

include("testsets.jl")

# linsearch  method
include("linesearch.jl")

# specialised Krylov subspace methods using finite-differences
include("fdkrylov.jl")

# the Newton-Krylov solver based on C T Kelley's nsoli.m
include("nsoli.jl")

# # the main transformed Newton solver
# include("solvermain.jl")

# # interfaces to Optim.jl and NLsolve.jl
# include("interfaces.jl")

end