text:
module WeightedApproximation

using ..Platforms, ..FrameFunInterface, ..ApproximationProblems, BasisFunctions,
    LinearAlgebra

import ..FrameFunInterface: AZ_A, AZ_Z, solver
import BasisFunctions: approximate

struct WeightedApproximationStyle <: ProblemStyle
    weight
end

function approximate(pstyle::WeightedApproximationStyle, samplingstyle::SamplingStyle, solverstyle::SolverStyle, fun, ap;
            verbose = false, options...)
    dict = dictionary(ap)

    verbose && println("Approximate: using problem style $pstyle")
    verbose && println("Approximate: using sampling style $samplingstyle")
	verbose && println("Approximate: using solver style $solverstyle")
    # Trigger computation of sampling parameter L first
    L = samplingparameter(samplingstyle, ap; verbose=verbose, options...)
    S = samplingoperator(samplingstyle, ap; verbose=verbose, options...)
	if verbose
		println()
		println("Approximate: sampling operator with size $(size(S)) is")
		println(S)
		println()
	end
    A = apply(S, dict; options...)
    B = apply(S, fun; options...)
    W = DiagonalOperator(apply(S, pstyle.weight; options...))

    C = solve(solverstyle, ap, W*A, W*B; problemstyle=pstyle, S=S, verbose=verbose, samplingstyle=samplingstyle, options...)
    F = Expansion(dictionary(ap), C)
    res = norm(W*(A*C-B))
    verbose && println("Approximate: ended with residual $res\n")
    F, A, B, C, S, L
end

export WFun
"""
Creates an approximate solution to `fun` using a weighed L2 measure, given by the `weight`.
"""
WFun(fun, weight, args...; options...) =
    Fun(fun, args...;
        problemstyle=WeightedApproximationStyle(weight), solverstyle=AZStyle(), samplingstyle=OversamplingStyle(),
        options...)

solver(pstyle::WeightedApproximationStyle, style::AZStyle, ap, A; options...) =
    solver(style, ap, A, AZ_Zt(pstyle, ap; options...); options...)

AZ_A(pstyle::WeightedApproximationStyle, ap; options...) =
    WAZ_W(pstyle, ap; options...)*AZ_A(DictionaryOperatorStyle(), ap; options...)

AZ_Z(pstyle::WeightedApproximationStyle, ap; weight_epsilon=1e-12, options...) =
    pinv(WAZ_W(pstyle, ap; options...), weight_epsilon)*AZ_Z(DictionaryOperatorStyle(), ap; options...)

WAZ_W(pstyle, ap::ApproximationProblem; samplingstyle=SamplingStyle(ap), options...) =
    DiagonalOperator(apply(samplingoperator(samplingstyle, ap; options...), pstyle.weight; options...))

end


text:
module Windows
using ..Util
import SpecialFunctions: besseli
using LinearAlgebra: Diagonal, SymTridiagonal, eigen!, mul!, rmul!
using FFTW

export  rect,
        hanning,
        hamming,
        tukey,
        cosine,
        lanczos,
        triang,
        bartlett,
        gaussian,
        bartlett_hann,
        blackman,
        kaiser,
        dpss,
        dpsseig

# this gets shared between the docstrings of all the window functions because
# they share these args.
const zerophase_docs = """
If `zerophase` is `false` (the default) the window is centered around index
`(n+1)/2`, which is commonly used for FIR filter design. These are often used
for FIR filter design, and usually odd-length. Note that for even-length windows
this will cause the window center to be located between samples.

If `zerophase` is `true` the window is centered around index 1 (with the
negative half wrapped to the end of the vector). Additionally this creates a
"periodic" window, which means that if there is no padding then the left and
right endpoints of the window wrap around to the same sample, so the window
length is the same as an `n+1`-length non-`zerophase` window. Alternatively you
can think of the continuous `zerophase` window being of width `n` and the
non-`zerophase` window being of length `n-1`. `zerophase` windows are often used
in FFT processing, and are usually even-length.
"""

function twoD_docs(arg=nothing)
    argstr = if arg !== nothing
        "`$arg`, "
    else
        ""
    end
    """
    Providing a `dims` `Tuple` rather than a single `n` constructs a 2D window.
    $argstr`padding` and `zerophase` can then be given either as a single value
    for both horizontal and vertical or a 2-tuple to specify them separately.
    """
end

"""
    padplot(plotstr)

Takes a multiline string and pre-pads it so that it shows up as
preformatted code when included in a docstring.
"""
function padplot(plotstr)
    minpad = typemax(Int)
    lines = split(plotstr, "\n")
    # count the minimum number of spaces preceeding any line
    for line in lines
        line == "" && continue
        minpad = min(minpad, length(match(r"^ *", line).match))
    end

    pad = " " ^ (4-minpad)
    join((string(pad, line) for line in lines), "\n")
end

include("winplots.jl")

"""
    makewindow(winfunc::Function, n::Integer, padding::Integer, zerophase::Bool)

Generate a discrete window vector of the given length with `padding` zeros.
`winfunc` should be a function giving the window value in the range of
[-0.5, 0.5]. The function is assumed to be 0 outside of this range.

$zerophase_docs

Example:

```
function hanning(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        0.5*(1+cos(2pi*x))
    end
end
```
"""
function makewindow(winfunc::Function, n::Integer, padding::Integer, zerophase::Bool)
    if n < 0
        throw(ArgumentError("`n` must be nonnegative"))
    end
    if padding < 0
        throw(ArgumentError("`padding` must be nonnegative"))
    end
    win = zeros(n+padding)
    if n == 1
        win[1] = winfunc(0.0)
    elseif zerophase
        # note that the endpoint of the window gets set in both lines. In the
        # unpadded case this will set the same index (which shouldn't make a
        # difference if the window is symmetric), but it's necessary for when
        # there's padding, which ends up in the center of the vector length
        # n÷2+1
        win[1:n÷2+1] .= winfunc.(range(0.0, stop=(n÷2)/n, length=n÷2+1))
        # length n÷2
        win[end-n÷2+1:end] .= winfunc.(range(-(n÷2)/n, stop=-1/n, length=n÷2))
    else
        win[1:n] .= winfunc.(range(-0.5, stop=0.5, length=n))
    end

    win
end

#
# Window functions
#

# equations rendered into Unicode with https://arthursonzogni.com/Diagon/

"""
$rect_winplot

    rect(n::Integer; padding::Integer=0, zerophase::Bool=false)
    rect(dims; padding=0, zerophase=false)

Rectangular window of length `n`, padded with `padding` zeros. This window is 1
within the window, and 0 outside of it.

$(twoD_docs())

$zerophase_docs
"""
function rect(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        1.0
    end
end

"""
$hanning_winplot

    hanning(n::Integer; padding::Integer=0, zerophase::Bool=false)
    hanning(dims; padding=0, zerophase=false)

Hanning window of length `n` with `padding` zeros. The Hanning (or Hann) window
is a raised-cosine window that reaches zero at the endpoints.

The window is defined by sampling the continuous function:

           1 + cos(2πx)
    w(x) = ──────────── = cos²(πx)
                2

in the range `[-0.5, 0.5]`

The `hanning` window satisfies the Constant Overlap-Add (COLA) property with an
hop of 0.5, which means that adding together a sequence of delayed windows with
50% overlap will result in a constant signal. This is useful when synthesizing
a signal from a number of overlapping frames (each with a roughly rectangular
window), to eliminate windowing amplitude modulation.

Note that the `hanning` window is the `cosine` window squared.

$(twoD_docs())

$zerophase_docs
"""
function hanning(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        0.5*(1+cos(2pi*x))
    end
end

"""
$hamming_winplot

    hamming(n::Integer; padding::Integer=0, zerophase::Bool=false)
    hamming(dims; padding=0, zerophase=false)

Hamming window of length `n` with `padding` zeros. The Hamming window does not
reach zero at the endpoints and so has a shallower frequency roll-off when
compared to the Hanning window, but is designed to cancel the first side-lobe.

The window is defined by sampling the continuous function:

    w(x) = 0.54 + 0.46*cos(2pi*x)

in the range `[-0.5, 0.5]`

$(twoD_docs())

$zerophase_docs
"""
function hamming(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        0.54 + 0.46*cos(2pi*x)
    end
end

"""
$tukey_winplot

    tukey(n::Integer, α::Real; padding::Integer=0, zerophase::Bool=false)
    tukey(dims, α; padding=0, zerophase=false)

Tukey window of length `n` with `padding` zeros. The Tukey window has a flat top
and reaches zero at the endpoints, with a sinusoidal transition area
parameterized by `α`. For `α == 0`, the window is equivalent to a rectangular
window. For `α == 1`, the window is a Hann window.

The window is defined by sampling the continuous function:

           ⎛              ⎛    ⎛    1 - α⎞⎞
           ⎜      1 + cos ⎜2πα ⎜x + ─────⎟⎟             1 - α
           ⎜              ⎝    ⎝      2  ⎠⎠         x ≤ ─────
           ⎜      ─────────────────────────               2
           ⎜                  2
           ⎜
    w(x) = ⎜      1                                 -α/2 < x ≤ α/2
           ⎜
           ⎜              ⎛    ⎛    1 - α⎞⎞
           ⎜      1 + cos ⎜2πα ⎜x - ─────⎟⎟             1 - α
           ⎜              ⎝    ⎝      2  ⎠⎠         x > ─────
           ⎜      ─────────────────────────               2
           ⎝                  2

in the range `[-0.5, 0.5]`

$(twoD_docs("α"))

$zerophase_docs
"""
function tukey(n::Integer, α::Real; padding::Integer=0, zerophase::Bool=false)
    # check that α is reasonable
    !(0 <= α <= 1) && error("α must be in the range 0 <= α <= 1.")

    # if α is less than machine precision, call it zero and return the
    # rectangular window for this length.  if we don't short circuit this
    # here, it will blow up below.
    abs(α) <= eps() && return rect(n; padding=padding, zerophase=zerophase)

    makewindow(n, padding, zerophase) do x
        if x <= -(1-α)/2
            0.5*(1 + cos(2pi/α*(x+(1-α)/2)))
        elseif x <= (1-α)/2
            1.0
        else
            0.5*(1 + cos(2pi/α*(x-(1-α)/2)))
        end
    end
end

"""
$cosine_winplot

    cosine(n::Integer; padding::Integer=0, zerophase::Bool=false)
    cosine(dims; padding=0, zerophase=false)

Cosine window of length `n` with `padding` zeros. The cosine window is the first
lobe of a cosine function (with the zero crossings at +/- π as endpoints). Also
called the sine window.

The window is defined by sampling the continuous function:

    w(x) = cos(πx)

in the range `[-0.5, 0.5]`

Note that the cosine window is the square root of the `hanning` window, so it is
sometimes used when you are applying the window twice, such as the analysis and
synthesis steps of an STFT.

$(twoD_docs())

$zerophase_docs
"""
function cosine(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        cos(pi*x)
    end
end

"""
$lanczos_winplot

    lanczos(n::Integer; padding::Integer=0, zerophase::Bool=false)
    lanczos(dims; padding=0, zerophase=false)

Lanczos window of length `n` with `padding` zeros. The Lanczos window is the
main lobe of a `sinc` function.

The window is defined by sampling the continuous function:

                      sin(2πx)
    w(x) = sinc(2x) = ────────
                         2πx

in the range `[-0.5, 0.5]`

$(twoD_docs())

$zerophase_docs
"""
function lanczos(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        sinc(2x)
    end
end

"""
$triang_winplot

    triang(n::Integer; padding::Integer=0, zerophase::Bool=false)
    triang(dims; padding=0, zerophase=false)

Triangular window of length `n` with `padding` zeros. The Triangular window does
not reach zero at the endpoints. For odd `n` the `triang` window is the center
`n` points of an `n+2`-point [`bartlett`](@ref) window (i.e. the samples just
outside the window would be zero). For even `n` the window slope is the same as
the `n-1` window but delayed by a half sample so the zero points would be 1/2
sample past the ends of the window.

The window is defined by sampling the continuous function:

            ⎛    2(n-1)
            ⎜1 - ────── abs(x)     n is even
            ⎜       n
    w(x) =  ⎜
            ⎜    2(n-1)
            ⎜1 - ────── abs(x)     n is odd
            ⎝     n+1

in the range `[-0.5, 0.5]`.

$(twoD_docs())

$zerophase_docs

When `zerophase` is `true` substitute `n+1` for `n` in the above window
expressions.
"""
function triang(n::Integer; padding::Integer=0, zerophase::Bool=false)
    # for the purpose of calculating the slope of the window, consider `n` to be
    # 1 larger to compensate for the fact that `zerophase` gives a periodic
    # window
    m = zerophase ? n+1 : n
    scale = iseven(m) ? 2(m-1)/m : 2(m-1)/(m+1)
    makewindow(n, padding, zerophase) do x
        1 - scale*abs(x)
    end
end

"""
$bartlett_winplot

    bartlett(n::Integer; padding::Integer=0, zerophase::Bool=false)
    bartlett(dims; padding=0, zerophase=false)

Bartlett window of length `n`. The Bartlett window is a triangular window that
reaches 0 at the endpoints. This is equivalent to convolving two rectangular
windows of length `(n-1)/2` and adding the zero endpoints. See `triang` for a
window that does not reach zero at the endpoints.

The window is defined by sampling the continuous function:

    1 - abs(2x)

in the range `[-0.5, 0.5]`

$(twoD_docs())

$zerophase_docs
"""
function bartlett(n::Integer; padding::Integer=0, zerophase::Bool=false)
    makewindow(n, padding, zerophase) do x
        1 - abs(2x)
    end
end

"""
$gaussian_winplot

    gaussian(n::Integer, σ::Real; padding::Integer=0, zerophase::Bool=false)
    gaussian(dims, σ; padding=0, zerophase=false)

Gives an n-sample gaussian window defined by sampling the function:

            ⎛        2⎞
            ⎜-1   ⎛x⎞ ⎟
            ⎜── ⋅ ⎜─⎟ ⎟
            ⎝ 2   ⎝σ⎠ ⎠
    w(x) = e

in the range `[-0.5,0.5]`. This means that for `σ=0.5` the endpoints of the
window will correspond to 1 standard deviation away from the center.

$(twoD_docs("σ"))

$zerophase_docs
"""
function gaussian(n::Integer, σ::Real; padding::Integer=0, zerophase::Bool=false)
    σ > 0.0 || error("σ must be positive")
    makewindow(n, padding, zerophase) do x
        exp(-0.5*(x/σ)^2)
    end
end

"""
$bartlett_hann_winplot

    bartlett_hann(n::Integer; padding::Integer=0, zerophase::Bool=false)
    bartlett_hann(dims; padding=0, zerophase=false)

Bartlett-Hann window of length `n` with `padding` zeros. The Bartlett-Hann
window is a weighted sum of the Bartlett and Hann windows.

The window is defined by sampling the continuous function:

    w(x) = 0.62 - 0.48*abs(x) + 0.38*cos(2π*x)

in the range `[-0.5, 0.5]`

$(twoD_docs())

$zerophase_docs
"""
function bartlett_hann(n::Integer; padding::Integer=0, zerophase::Bool=false)
    a0, a1, a2 = 0.62, 0.48, 0.38
    makewindow(n, padding, zerophase) do x
        a0 - a1*abs(x) + a2*cos(2pi*x)
    end
end

"""
$blackman_winplot

    blackman(n::Integer; padding::Integer=0, zerophase::Bool=false)
    blackman(dims; padding=0, zerophase=false)

Approximates the "Exact" Blackman window. This is the generalized Blackman
window with α = 0.16.

The window is defined by sampling the continuous function:

    w(x) = 0.42 + 0.5*cos(2π*x) + 0.08*cos(4π*x)

in the range `[-0.5, 0.5]`

$(twoD_docs())

$zerophase_docs
"""
function blackman(n::Integer; padding::Integer=0, zerophase::Bool=false)
    a0, a1, a2 = 0.42, 0.5, 0.08
    makewindow(n, padding, zerophase) do x
        a0 + a1*cos(2pi*x) + a2*cos(4pi*x)
    end
end

"""
$kaiser_winplot

    kaiser(n::Integer, α::Real; padding::Integer=0, zerophase::Bool=false)
    kaiser(dims, α; padding=0, zerophase=false)

Kaiser window of length `n` parameterized by `α`. The Kaiser window approximates
the DPSS window (given by `dpss`), using a simplified definition relying on a
Bessel function. Larger values for `α` give a wider main lobe but have lower
sidelobes. Typically `α` is set around 3.

The window is defined by sampling the continuous function:

             ⎛  ⎛   _________⎞⎞
             ⎜  ⎜  ╱        2⎟⎟
    w(x) = I₀⎝πα⎝╲╱ 1 - (2x) ⎠⎠
           ────────────────────
                   I₀(πα)

in the range `[-0.5, 0.5]`

Where I₀(⋅) is the zeroth-order modified Bessel function of the first kind.

$(twoD_docs("α"))

$zerophase_docs
"""
function kaiser(n::Integer, α::Real; padding::Integer=0, zerophase::Bool=false)
    pf = 1.0/besseli(0,pi*α)
    makewindow(n, padding, zerophase) do x
        pf*besseli(0, pi*α*(sqrt(1 - (2x)^2)))
    end
end

# Discrete prolate spheroid sequences (Slepian tapers)
#
# See Gruenbacher, D. M., & Hummels, D. R. (1994). A simple algorithm
# for generating discrete prolate spheroidal sequences. IEEE
# Transactions on Signal Processing, 42(11), 3276-3278.
"""
$dpss_winplot

    dpss(n::Integer, nw::Real, ntapers::Integer=iceil(2*nw)-1;
         padding::Integer=0, zerophase::Bool=false)

The first `ntapers` discrete prolate spheroid sequences (Slepian
tapers) as an `n` × `ntapers` matrix. The signs of the tapers
follow the convention that the first element of the skew-symmetric
(odd) tapers is positive. The time-bandwidth product is given by
`nw`.

The DPSS window maximizes the energy concentration in the main lobe.

$zerophase_docs
"""
function dpss(n::Integer, nw::Real, ntapers::Integer=ceil(Int, 2*nw)-1;
              padding::Integer=0, zerophase::Bool=false)
    if isodd(n) && zerophase
        throw(ArgumentError("`dpss` does not currently support odd-length zerophase windows"))
    end
    # TODO there's probably a more efficient way to compute the zero-phase
    # version of this
    if zerophase
        n += 1
    end
    0 < ntapers <= n || error("ntapers must be in interval (0, n]")
    0 <= nw < n/2 || error("nw must be in interval [0, n/2)")

    # Construct symmetric tridiagonal matrix
    v = cospi(2*nw/n)
    dv = Vector{Float64}(undef, n)
    ev = Vector{Float64}(undef, n - 1)
    @inbounds dv[1] = v * abs2((n - 1) / 2)
    @inbounds @simd for i = 1:(n-1)
        dv[i + 1] = v * abs2((n - 1) / 2 - i)
        ev[i] = 0.5 * (i * n - i^2)
    end
    mat = SymTridiagonal(dv, ev)

    # Get tapers
    eigvec = eigen!(mat, n-ntapers+1:n).vectors
    rv = reverse(eigvec, dims=2)::Matrix{Float64}

    # Slepian's convention; taper starts with a positive element
    sgn = ones(size(rv, 2))
    for i = 2:2:size(rv, 2)
        s = zero(eltype(rv))
        for j = 1:n
            s = sign(rv[j, i])
            s != 0 && break
        end
        @assert s != 0
        sgn[i] = s
    end
    rmul!(rv, Diagonal(sgn))

    if zerophase
        # remove the last element (should be equal to the first element)
        rv = rv[1:end-1, :]
    end

    if padding > 0
        rv = [rv; zeros(padding, ntapers)]
    end

    if zerophase
        # TODO: figure out how to handle odd-length zerophase windows. See the
        # tests for one approach, but upsampling by 2 doesn't work here becasue
        # the overall amplitude varies with `n`
        rv = ifftshift(rv, 1)
    end

    rv
end

# Eigenvalues of DPSS, following Percival & Walden p. 390, exercise 8.1
# See also implementation in MNE:
# https://github.com/mne-tools/mne-python/blob/d7082cf909ccab581667bc1f1ed3c23e6a24b567/mne/time_frequency/multitaper.py#L226
"""
    dpsseig(A, nw)

Eigenvalues of the DPSS matrix, representing the ratios of the
power within the main lobe to the total power (main and sidelobes).
`A` is the output of [`dpss`](@ref), and `nw` is the
time-bandwidth product provided to [`dpss`](@ref) as input.
"""
function dpsseig(A::Matrix{Float64}, nw::Real)
    0 <= nw < size(A, 1)/2 || error("nw must be in interval [0, n/2)")

    w = nw/size(A, 1)

    # Compute coefficients
    seq = Vector{Float64}(undef, size(A, 1))
    seq[1] = 1.0
    for i = 1:size(A, 1)-1
        seq[i+1] = 2 * sinc(2w*i)
    end

    q = Vector{Float64}(undef, size(A, 2))
    nfft = nextfastfft(2*size(A, 1)-1)

    tmp1 = Vector{Float64}(undef, nfft)
    tmp2 = Vector{Complex{Float64}}(undef, nfft >> 1 + 1)
    p1 = plan_rfft(tmp1)
    p2 = plan_brfft(tmp2, nfft)

    for i = 1:size(A, 2)
        fill!(tmp1, 0)
        copyto!(tmp1, 1, A, (i-1)*size(A, 1)+1, size(A, 1))
        mul!(tmp2, p1, tmp1)
        for j = 1:length(tmp2)
            @inbounds tmp2[j] = abs2(tmp2[j])
        end
        mul!(tmp1, p2, tmp2)

        eig = 0.0
        for j = 1:size(A, 1)
            eig += seq[j]*tmp1[j]
        end
        q[i] = 2w * eig / nfft
    end

    q
end

# tensor product window functions in 2-d, defined as w(x,y) = w(x) \otimes w(y)

# convert the given argument to a 2-Tuple so we can distribute it to the
# input windows when making a 2D window
argdup(arg::Tuple) = arg
argdup(arg::Real) = (arg, arg)
const IntegerOr2 = Union{Tuple{<:Integer, <:Integer}, Integer}
const RealOr2 = Union{Tuple{<:Real, <:Real}, Real}
const BoolOr2 = Union{Tuple{Bool, Bool}, Bool}

for func in (:rect, :hanning, :hamming, :cosine, :lanczos,
             :triang, :bartlett, :bartlett_hann, :blackman)
    @eval begin
        function $func(dims::Tuple; padding::IntegerOr2=0,
                                    zerophase::BoolOr2=false)
            length(dims) == 2 || throw(ArgumentError("`dims` must be length 2"))
            paddings = argdup(padding)
            zerophases = argdup(zerophase)
            w1 = $func(dims[1], padding=paddings[1], zerophase=zerophases[1])
            w2 = $func(dims[2], padding=paddings[2], zerophase=zerophases[2])
            w1 * w2'
        end
    end
end

for func in (:tukey, :gaussian, :kaiser)
    @eval begin
        function $func(dims::Tuple, arg::RealOr2;
                       padding::IntegerOr2=0, zerophase::BoolOr2=false)
            length(dims) == 2 || throw(ArgumentError("`dims` must be length 2"))
            args = argdup(arg)
            paddings = argdup(padding)
            zerophases = argdup(zerophase)
            w1 = $func(dims[1], args[1]; padding=paddings[1], zerophase=zerophases[1])
            w2 = $func(dims[2], args[2]; padding=paddings[2], zerophase=zerophases[2])
            w1 * w2'
        end
    end
end

end # end module definition

ext: jl


text:
module DataManipulation
  files = ["DataTransform"]
  for f in files
    include("$f.jl")
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
#### Simulations for "Post-error with attractor network" Berlemont & Nadal###
# Author: Kevin Berlemont
# 10 July 2018
# Shared privately with the reviewers.
########################################################
using Distributions,CSV,Query
using DataFrames
using Plots
gr()
srand(1234)

param = Dict("T_rest"=>500.0,
"mu0"=>30.0,
"noise_amp"=>0.02,
"I0E1"=>0.3255,
"I0E2"=>0.3255,
"JN11"=>0.2609,
"JN22"=>0.2609,
"JN12"=>0.0497,
"JN21"=>0.0497,
"JAext" => 0.00052,
"Tnmda" => 100.0,
"Tampa" => 2.0,
"gamma" => 0.641,
"Tstim"=>5000.0,
"coh"=>5.0,
"false_time"=>200.0,
"Threshold"=>20.0
)

fixed_param=Dict("a"=>270.0,"b"=>108.0,"d"=>0.1540)
#Parameters of the simulation for EUler methods
simu_param=Dict("dt"=>0.5,"time_wind"=>4.0,"slide_wind"=>4.0) #4=2/dt

using FileIO
using DataFrames
using CSV

Tr=500
ICD=0.047
println(Tr)

for i=1:50
    print(i)
    print("\n")
    for Tr in [200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500,1600,
        1700,1800,1900,2000,3000,4000,5000]
        for c in 0.0:1.0:20.0

        c_min=-c
        c_max = c
        number = 2

        distri = make_distributions(c_min,c_max,number,5000)

        param["T_rest"] = Tr
        If=0.0

        result= sim_mf_decisionmaking_pattern_fast(distri,ICD,0.0,
        param,fixed_param,simu_param)

        CSV.write("./I47/Coh$c-Final200Exp$i-T$Tr-Ir$Ir.csv",result)

    end

end


text:
push!(LOAD_PATH, "../src/")
push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("SummationByParts"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("MeshMovement"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/Utils"))

using PdePumiInterface
using SummationByParts
using ODLCommonTools
using ArrayViews
using Utils
using MPI
# using MeshMovement
using FFD
using Base.Test

# Include functions that will be used for performing tests
include("./test_parallel_FFD_2D_mesh.jl")
include("./test_parallel_FFD_3D_mesh.jl")

#MPI.Finalize()


text:
include("benchmark.jl")

println("####################################################")
println("Ingore the first set of benchmarks. They are run so as to make sure that everything is compiled once.")
println("####################################################")
generate_benchmark_file(file_name = tempname())

println("####################################################")
println("Everything has been compiled once. The following are the final benchmarks.")
println("####################################################")
generate_benchmark_file()


text:
using ReservoirComputing

train_data = ones(Int, 2, 10)
rule = 30
generations = 3
expansion_size = 10
permutations = 2

reca = RECA_discrete(train_data, rule, generations, expansion_size, permutations)

@test isequal(reca.train_data, train_data)
@test isequal(reca.rule, rule)
@test isequal(reca.generations, generations)
@test isequal(reca.expansion_size, expansion_size)
@test isequal(reca.permutations, permutations)


text:
@testset "1027.longest-arithmetic-subsequence.jl" begin
    @test longest_arith_seq_length([20, 1, 15, 3, 10, 5, 8]) == 4
    @test longest_arith_seq_length([3, 6, 9, 12]) == 4
    @test longest_arith_seq_length([9, 4, 7, 2, 10]) == 3
    @test longest_arith_seq_length([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) == 10
end


text:
	"""
	options = ContinuationPar(dsmin = 1e-4,...)

Returns a variable containing parameters to affect the `continuation` algorithm used to solve `F(x,p) = 0`.

# Arguments
- `dsmin, dsmax` are the minimum, maximum arclength allowed value. It controls the density of points in the computed branch of solutions.
- `ds = 0.01` is the initial arclength.
- `theta` is a parameter in the arclength constraint. It is very **important** to tune it. See the docs of [`continuation`](@ref). We quote them here "It should be tuned for the continuation to work properly especially in the case of large problems where the < x - x_0, dx_0 > component in the constraint equation might be favoured too much. Also, large thetas favour p as the corresponding term in N involves the term 1-theta."
- `pMin, pMax` allowed parameter range for `p`
- `maxSteps = 100` maximum number of continuation steps
- `newtonOptions::NewtonPar`: options for the Newton algorithm
- `saveToFile = false`: save to file. A name is automatically generated.
- `saveSolEveryStep::Int64 = 0` at which continuation steps do we save the current solution
- `plotEveryStep = 10`

## Handling eigen elements, their computation is triggered by the argument `detectBifurcation` (see below)
- `nev = 3` number of eigenvalues to be computed. It is automatically increased to have at least `nev` unstable eigenvalues. To be set for proper  bifurcation detection. See [Detection of bifurcation points](@ref) for more informations.
- `saveEigEveryStep = 1`	record eigen vectors every specified steps. **Important** for memory limited resource, *e.g.* GPU.
- `saveEigenvectors	= true`	**Important** for memory limited resource, *e.g.* GPU.

## Handling bifurcation detection
- `precisionStability = 1e-10` lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbits
- `detectFold = true` detect Fold bifurcations? It is a useful option although the detection of Fold is cheap. Indeed, it may happen that there is a lot of Fold points and this can saturate the memory in memory limited devices (e.g. on GPU)
- `detectBifurcation::Int` ∈ {0, 1, 2, 3} If set to 0, nothing is done. If set to 1, the eigen-elements are computed. If set to 2, the bifurcations points are detected during the continuation run, but not located precisely. If set to 3, a bisection algorithm is used to locate the bifurcations points (slower). The possibility to switch off detection is a useful option. Indeed, it may happen that there are a lot of bifurcation points and this can saturate the memory of memory limited devices (e.g. on GPU)
- `dsminBisection = 1e-16` dsmin for the bisection algorithm for locating bifurcation points
- `nInversion = 2` number of sign inversions in bisection algorithm
- `maxBisectionSteps = 15` maximum number of bisection steps
- `tolBisectionEigenvalue = 1e-16` tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps

## Handling `ds` adaptation (see [`continuation`](@ref) for more information)
- `a  = 0.5` aggressiveness factor. It is used to adapt `ds` in order to have a number of newton iterations per continuation step roughly constant. The higher `a` is, the larger the step size `ds` is changed at each continuation step.
- `thetaMin = 1.0e-3` minimum value of `theta`
- `doArcLengthScaling = false` trigger further adaptation of `theta`

## Handling event detection
- `detectEvent::Int` ∈ {0, 1, 2} If set to 0, nothing is done. If set to 1, the event locations are seek during the continuation run, but not located precisely. If set to 2, a bisection algorithm is used to locate the event (slower).
- `tolParamBisectionEvent = 1e-16` tolerance on parameter to locate event

## Misc
- `finDiffEps::T  = 1e-9` ε used in finite differences computations
- `detectLoop` [WORK IN PROGRESS] detect loops in the branch and stop the continuation

!!! tip "Mutating"
    For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package `Setfield.jl` to drastically simplify the mutation of different fields. See tutorials for more examples.
"""
@with_kw struct ContinuationPar{T, S <: AbstractLinearSolver, E <: AbstractEigenSolver}
	# parameters for arclength continuation
	dsmin::T	= 1e-3
	dsmax::T	= 1e-1
	@assert dsmax >= dsmin "You must provide a valid interval (ordered) for ds"
	ds::T		= 1e-2;		@assert dsmax >= abs(ds) >= dsmin
	@assert dsmin > 0 "The interval for ds must be positive"
	@assert dsmax > 0

	# parameters for scaling arclength step size
	theta::T					= 0.5 		# parameter in the dot product used for the extended system
	doArcLengthScaling::Bool  	= false
	gGoal::T					= 0.5
	gMax::T						= 0.8
	thetaMin::T					= 1.0e-3
	a::T						= 0.5  		# aggressiveness factor
	tangentFactorExponent::T 	= 1.5

	# parameters bound
	pMin::T	= -1.0
	pMax::T	=  1.0; 			@assert pMax >= pMin "You must provide a valid interval [pMin, pMax]"

	# maximum number of continuation steps
	maxSteps::Int64  = 100

	# Newton solver parameters
	finDiffEps::T  = 1e-9 					# constant for finite differences
	newtonOptions::NewtonPar{T, S, E} = NewtonPar()
	η::T = 150.								# parameter to estimate tangent at first point

	saveToFile::Bool = false 				# save to file?
	saveSolEveryStep::Int64 = 0				# what steps do we save the current solution

	# parameters for eigenvalues
	nev::Int64 = 3 							# number of eigenvalues
	saveEigEveryStep::Int64 = 1				# what steps do we keep the eigenvectors
	saveEigenvectors::Bool	= true			# useful options because if puts a high memory pressure

	plotEveryStep::Int64 = 10

	# handling bifurcation points
	precisionStability::T = 1e-10			# lower bound for stability of equilibria and periodic orbits
	detectFold::Bool = true 				# detect fold points?
	detectBifurcation::Int64 = 0			# detect other bifurcation points?
	dsminBisection::T = 1e-16				# dsmin for the bisection algorithm when locating bifurcation points
	nInversion::Int64 = 2					# number of sign inversions in bisection algorithm
	maxBisectionSteps::Int64 = 15			# maximum number of bisection steps
	tolBisectionEigenvalue::T = 1e-16 		# tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps. Must be small otherwise Shooting and friends will fail detecting bifurcations.

	# handling event detection
	detectEvent::Int64 = 0				# event location
	tolParamBisectionEvent::T = 1e-16	# tolerance on value of parameter

	@assert iseven(nInversion) "The option `nInversion` number must be even"
	@assert 0 <= detectBifurcation <= 3 "The option `detectBifurcation` must belong to {0,1,2,3}"
	@assert 0 <= detectEvent <= 2 "The option `detectEvent` must belong to {0,1,2}"
	@assert (detectBifurcation > 1 && detectEvent == 0) || (detectBifurcation <= 1 && detectEvent >= 0)  "One of these options must be disabled detectBifurcation = $detectBifurcation and detectEvent = $detectEvent"
	@assert tolBisectionEigenvalue >= 0 "The option `tolBisectionEigenvalue` must be positive"
	detectLoop::Bool = false				# detect if the branch loops
	@assert 0 <= theta <=1 "theta must belong to [0, 1]"
	@assert plotEveryStep > 0 "plotEveryStep must be positive. You can turn off plotting by passing plot = false to `continuation`"
	@assert ~(detectBifurcation > 1 && saveEigEveryStep > 1) "We must at least save all eigenvalues for detection of bifurcation points. Please use saveEigEveryStep = 1 or detectBifurcation = 1."
end

@inline computeEigenElements(cp::ContinuationPar) = cp.detectBifurcation > 0
@inline computeEigenvalues(cp::ContinuationPar) = cp.detectBifurcation > 0
@inline saveEigenvectors(cp::ContinuationPar) = computeEigenvalues(cp) * (cp.saveEigenvectors)

# clamp ds value
clampDs(ds, contparams::ContinuationPar) = sign(ds) * clamp(abs(ds), contparams.dsmin, contparams.dsmax)


text:
debug_level = get(ENV, "MODERNGL_DEBUGGING", "false") == "true"
open(joinpath(@__DIR__, "deps.jl"), "w") do io
    println(io, "const enable_opengl_debugging = $debug_level")
end


text:
# Problem 5 - Project Euler
# http://projecteuler.net/index.php?section=problems&id=5

function p005()
    ans = 1
    for i in 1:20
        ans = lcm(ans, i)
    end
    ans
end

println(p005())

max_stars_repo_name: emergent/ProjectEuler

max_issues_repo_name: emergent/ProjectEuler

max_forks_repo_name: emergent/ProjectEuler


text:
function write_vtp(dataset::VTKPolyData, filepath_no_ext::String)
    polydata = PyVTK(dataset)
    return write_vtp(polydata, filepath_no_ext)
end

function write_vtp(polydata::PyObject, filepath_no_ext::String)
    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetFileName(string(filepath_no_ext, ".vtp"))
    writer.SetDataModeToBinary()
    writer.SetDataModeToAppended()
    writer.SetCompressorTypeToZLib()
    writer.SetInputData(polydata)
    writer.Update()
    return string(filepath_no_ext, ".vtp")
end


text:
## ----------  HS GPA or other endowment percentile cutoff

"""
	$(SIGNATURES)

Switches governing admissions by cutoff rule.
"""
mutable struct AdmissionsCutoffSwitches{I1, F1 <: Real} <: AbstractAdmissionsSwitches{I1, F1}
    nColleges :: I1
    # Name of the variable that holds the individual percentiles
    # pctVar :: Symbol
    # Minimum HS GPA percentile required for each college; should be increasing
    minPctV :: Vector{F1}
    # Minimum probability of all college sets
    minCollSetProb :: F1
end

"""
	$(SIGNATURES)

Admissions are governed by a single indicator, such as a test score. Students can attend colleges for which they qualify in the sense that their indicator exceeds the college's cutoff value. Students may be allowed to attend other colleges with a fixed probability.

This means that there are only two probabilities. With a high probability, the student can attend colleges `1:n` where `n` is the best college for which the student qualifies. With a low probability, the student may draw each college set `1 : m` where `m != n`.
"""
struct AdmissionsCutoff{I1, F1 <: Real} <: AbstractAdmissionsRule{I1, F1}
    switches :: AdmissionsCutoffSwitches{I1, F1}
end

function Base.show(io :: IO, a :: AdmissionsCutoff)
    nc = n_colleges(a);
    print(io, typeof(a), " with cutoff percentiles",
        round.(min_percentiles(a), digits = 2));
end

StructLH.describe(a :: AdmissionsCutoffSwitches) = [
    "Admission rule"  " ";
    # "Cutoff rule based on"  "$(a.pctVar)";
    "Min percentile by college"  "$(a.minPctV)"
];

function validate_admissions(a :: AdmissionsCutoff{I1, F1}) where {I1, F1}
    isValid = true;
    if any(diff(min_percentiles(a)) .<= 0.0)
        @warn "Min percentiles should be increasing"
        isValid = false;
    end
    isValid  ||  println(a)
    return isValid
end

make_admissions(switches :: AdmissionsCutoffSwitches{I1, F1}) where {I1, F1} =
    AdmissionsCutoff(switches);

make_test_adm_cutoff_switches(nc) =
    AdmissionsCutoffSwitches(nc,
        collect(range(0.0, 0.8, length = nc)), 0.05);
make_test_admissions_cutoff(nc) =
    AdmissionsCutoff(make_test_adm_cutoff_switches(nc));

min_percentiles(a :: AdmissionsCutoff{I1, F1}) where {I1, F1} = a.switches.minPctV;
# percentile_var(switches :: AdmissionsCutoffSwitches{I1, F1}) where {I1, F1} =
#     switches.pctVar;

# Highest college for which a student qualifies
# Last GPA cutoff that is smaller than student's endowment percentile.
function highest_college(a :: AdmissionsCutoff{I1, F1}, endowPct :: F2) where
    {I1, F1, F2 <: Real}

    @assert (zero(F2) <= endowPct <= one(F2))  "Expecting percentiles between 0 and 1";
    qIdx = findlast(x -> x <= endowPct, min_percentiles(a));
    @check qIdx >= 1
    return qIdx
end

function highest_college(a :: AdmissionsCutoff{I1, F1}, endowPctV :: AbstractVector{F2}) where {I1, F1, F2 <: Real}

    highV = Vector{Int}(undef, length(endowPctV));
    for (j, endowPct) in enumerate(endowPctV)
        highV[j] = highest_college(a, endowPct);
    end
    return highV
end

# Probability that one student draws one college set
function prob_coll_set(a :: AdmissionsCutoff{I1, F1},
    admProbFct :: AF1,
    iSet :: Integer, endowPct :: F2) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1, F2 <: Real}

    @assert (zero(F2) <= endowPct <= one(F2))  "Expecting percentiles between 0 and 1";
    qIdx = highest_college(a, endowPct);
    if qIdx == iSet
        # Best college the student qualifies for
        nSets = n_college_sets(a);
        prob = one(F1) - min_coll_set_prob(a) * (nSets - 1);
        # nNot = nSets - qIdx;
        # prob = (1.0 - min_coll_set_prob(a) * nNot) / qIdx;
    else
        # Student does not qualify
        prob = min_coll_set_prob(a);
    end
    return prob
end

# Prob that a given student draws each college set
# If student qualifies for college 2 out of 5, the probabilities are:
#   (1 - minProb * 3) / 2  for first 2 (where he qualifies)
#   minProb for the last 3 where he does not qualify
function prob_coll_sets(a :: AdmissionsCutoff{I1, F1},
    admProbFct :: AF1,
    endowPct :: F2) where {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1, F2 <: Real}

    @assert (zero(F2) <= endowPct <= one(F2))  "Expecting percentiles between 0 and 1";
    nSets = n_college_sets(a);
    probV = [prob_coll_set(a, admProbFct, iSet, endowPct)  for iSet = 1 : nSets];
    @check sum(probV) ≈ 1
    return probV
end

# function prob_coll_sets(a :: AdmissionsCutoff{I1, F1}, hsGpaPctV :: AbstractVector{F2}) where {I1, F1, F2 <: Real}

#     nSets = n_college_sets(a);
#     n = length(hsGpaPctV);
#     prob_jsM = Matrix{F1}(undef, n, nSets);
#     for j = 1 : n
#         prob_jsM[j, :] = prob_coll_sets(a, hsGpaPctV[j]);
#     end

#     @assert all_at_least(prob_jsM, a.switches.minCollSetProb)
#     @assert all_at_most(prob_jsM, 1.0)
#     @check all(isapprox.(sum(prob_jsM, dims = 2), 1.0, atol = 1e-6))
#     return prob_jsM
# end

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


text:
function isisogram(s::AbstractString)

end


text:
module OutlierDetectionNeighbors
    using OutlierDetectionInterface
    using OutlierDetectionInterface:SCORE_UNSUPERVISED
    const OD = OutlierDetectionInterface

    import NearestNeighbors
    const NN = NearestNeighbors

    import Distances
    const DI = Distances

    include("utils.jl")
    include("models/abod.jl")
    include("models/cof.jl")
    include("models/dnn.jl")
    include("models/knn.jl")
    include("models/lof.jl")

    const UUID = "51249a0a-cb36-4849-8e04-30c7f8d311bb"
    const MODELS = [:ABODDetector,
                    :COFDetector,
                    :DNNDetector,
                    :KNNDetector,
                    :LOFDetector]

    for model in MODELS
        @eval begin
            OD.@default_frontend $model
            OD.@default_metadata $model $UUID
            export $model
        end
    end
end


text:
using PdfFraud
using JSON
using CSV

# get representations and their sha256 codes
X, sha256, page_nums = PdfFraud.load_pdf_data()

# files
indexf = joinpath(PdfFraud.proto_path, "index.json")
labelf = "/home/vit/vyzkum/bulletproof/pdf_experiment/gt-clustering-pdfs.txt"

# get the labels
labels, labeled_sha256, filenames = PdfFraud.get_labels(sha256, indexf, labelf);
labels_x = vcat(map(x -> repeat([x[1]], x[2]), zip(labels, page_nums))...)

######### there are some pdfs that share sha256 ##########
sum(map(x -> x != nothing, labels)) # some labeled data is not present in the dataset
shas = Array{Any,1}()
for fn in labeldata[!,:file]
	sha = sha256[map(x -> x == fn, filenames)]
	if length(sha) > 0
		push!(shas, sha[1])
	else
		println(fn)
		push!(shas, nothing)
	end
end

fn = "2016_06_14_Faktura.pdf"
filenames .== fn
is = index["original_name"] .== fn
index["original_name"][is]
sha = "0836db253a5b6cf36737a1e98806377f32b798e6e2d6d44fc99718fa65116cda"
labeldata[!,:file] .== fn
filenames[filenames .== fn]
fname = index["original_name"][map(y -> sha == y, index["sha256"])]


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using Test

for file in readdir(@__DIR__)
    if file in ["Utilities.jl"]
        continue
    end
    @testset "$(file)" begin
        include(file)
    end
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
function sum_homogeneous(k, bs::AbstractVector{Blade{T}}) where {T}
    kvec = zero(KVector{k}, T)
    for b ∈ bs
        add!(kvec, b)
    end
    kvec
end

sum_homogeneous(k, bs::Blade{T}...) where {T} = sum_homogeneous(k, collect(bs))

function (+)(bs::Blade{T}...) where {T}
    grades = unique(grade.(bs))
    if length(grades) == 1
        k = grades[1]
        k == 0 && return scalar(reduce(+, getproperty.(bs, :coef)))
        return sum_homogeneous(k, SVector{length(bs)}(bs))
    else
        mvec = zero(Multivector, T)
        for b ∈ bs
            add!(mvec, b)
        end
        return mvec
    end
end

@commutative (+)(b::Blade, kvec::KVector) = +(b, blades(kvec)...)
@commutative (+)(b::Blade, mv::Multivector) = +(b, blades(mv)...)

convert(T::Type{<:Blade}, b::Blade) = T(b.coef, b.index)
convert(T::Type{<:Blade}, x::Number) = scalar(x)

promote_rule(T1::Type{<:Blade}, T2::Type{<:Blade}) = Blade{promote_type(eltype(T1), eltype(T2))}
promote_rule(T1::Type{<:KVector}, T2::Type{<:KVector}) = Multivector
promote_rule(T1::Type{<:Multivector}, T2::Type{<:KVector}) = Multivector
# promote_rule(T1::Type{<:Blade}, T2::Type{<:Number}) = Blade{promote_type(eltype(T1), T2)}

(+)(kvecs::KVector{K}...) where {K} = KVector{K}(.+(getproperty.(kvecs, :coefs)...))
(+)(mvecs::Multivector...) = Multivector(.+(getproperty.(mvecs, :coefs)...))
(+)(x::GeometricAlgebraType, y::GeometricAlgebraType) = +(promote(x, y)...)
@type_commutative (+)(x::GeometricAlgebraType, y::Number) = +(promote(x, y)...)

sum(x::AbstractVector{<:GeometricAlgebraType}) = reduce(+, x)

(-)(x::Blade) = Blade(-x.coef, x.index)
(-)(x::Multivector) = Multivector(map(-, x.coefs))
(-)(x::KVector) = typeof(x)(map(-, x.coefs))
(-)(x::GeometricAlgebraType, y::GeometricAlgebraType) = x + (-y)

"""
    x ∧ y
Outer product of `x` with `y`.
"""
function ∧ end

"""
    x ⋅ y
Inner product of `x` with `y`. This product is in general non-associative, and is conventinally executed right to left in absence of parenthesis.
For example, `A ⋅ B ⋅ C == A ⋅ (B ⋅ C)`, and in most cases, `A ⋅ B ⋅ C ≠ (A ⋅ B) ⋅ C`.
"""
function ⋅ end

"""
    x ⦿ y
Scalar product between `x` and `y`.
"""
function ⦿ end

"""
    x ∨ y
Meet of `x` and `y`
"""
function ∨ end

"""
    lcontract(x, y)
Left contraction of `x` with `y`.
"""
function lcontract end

"""
    rcontract(x, y)
Right contraction of `x` with `y`.
"""
function rcontract end

function geom end

"""
Duality operator.
"""
function dual end

geom!(res, args...) = add!(res, geom(args...))

@generated op_result_type(op, T1::Type{<:KVector}, T2::Type{<:KVector}) = KVector{result_grade(op, grade(T1), grade(T2)), promote_type(vectype(T1), vectype(T2))}

(*)(x::Blade, y::Blade) = geom(x, y)

function (*)(x::GeometricAlgebraType, y::GeometricAlgebraType)
    res = zero(op_result_type(*, x, y))
    geom!(res, x, y)
    res
end

function geom!(res::Multivector, x::Blade, y::KVector)
    for b ∈ blades(y)
        geom!(res, x, b)
    end
end

function geom!(res::Multivector, x::KVector, y::Blade)
    for b ∈ blades(x)
        geom!(res, b, y)
    end
end

function geom!(res::Multivector, x::KVector, y::KVector)
    for b ∈ blades(x)
        geom!(res, b, y)
    end
end

function geom!(res::Multivector, x::Multivector, y)
    for kvec ∈ blades(x)
        geom!(res, kvec, y)
    end
end

function geom!(res::Multivector, x::Multivector, y::Multivector)
    for kvec ∈ blades(y)
        geom!(res, x, kvec)
    end
end

@commutative (*)(x::KVector, y::Number) = typeof(x)(x.coefs .* y, x.start)
@commutative (*)(x::Multivector, y::Number) = typeof(x)(x.coefs .* y)
@commutative (*)(x::Blade, y::Number) = Blade(x.coef * y, x.index)
(*)(xs::GeometricAlgebraType...) = reduce(*, xs)

(^)(x::GeometricAlgebraType, p::Integer) = prod(ntuple(i -> x, p))

@commutative (⋅)(x::GeometricAlgebraType, y::Number) = scalar(zero(promote_type(eltype(x), typeof(y))))
(⋅)(x::GeometricAlgebraType, y::Blade) = is_scalar(y) ? x ⋅ y.coef : grade_projection(x * y, Val(result_grade(⋅, grade(x), grade(y))))
(⋅)(x::Multivector, y::Blade) = sum(kx ⋅ y for kx ∈ kvectors(x))
(⋅)(x::Blade, y::GeometricAlgebraType) = is_scalar(x) ? x.coef ⋅ y : grade_projection(x * y, Val(result_grade(⋅, grade(x)), grade(x)))
(⋅)(x::Blade, y::Multivector) = sum(x ⋅ ky for ky ∈ kvectors(y))
(⋅)(x::Blade, y::Blade) = is_scalar(x) || is_scalar(y) ? scalar(zero(promote_type(eltype(x), eltype(y)))) : grade_projection(x * y, Val(result_grade(⋅, grade(x), grade(y))))
(⋅)(x::T, y::T) where {T<:Number} = scalar(zero(T))
(⋅)(x::Number, y::Number) = scalar(zero(promote_type(typeof(x), typeof(y))))

(⦿)(x::Number, y::Number) = x * y

@commutative (∧)(x::GeometricAlgebraType, y::Number) = x * y
(∧)(x::Number, y::Number) = scalar(x * y)

for op ∈ [:∧, :⋅, :⦿]

    # ⋅ is not associative
    fold = op == :⋅ ? foldr : reduce

    if op == :⦿
        @eval begin
            ($op)(x::GeometricAlgebraType, y::GeometricAlgebraType) = first(grade_projection(x * y, Val(0)).coefs)
            ($op)(x::Blade, y::Blade) = grade_projection(x * y, Val(0)).coef
        end
    else
        @eval ($op)(x::GeometricAlgebraType, y::GeometricAlgebraType) = grade_projection(x * y, Val(result_grade($op, grade(x), grade(y))))
    end

    @eval begin
        ($op)(x::Any, y::Any, z::Any...) = $fold($op, vcat(x, y, z...))
        ($op)(x::Multivector, y::Multivector) = sum($op(kx, ky) for kx ∈ kvectors(x), ky ∈ kvectors(y))
        ($op)(x::GeometricAlgebraType, y::Multivector) = sum($op(x, ky) for ky ∈ kvectors(y))
        ($op)(x::Multivector, y::GeometricAlgebraType) = sum($op(kx, y) for kx ∈ kvectors(x))
    end
end

permsign(i, j) =
    1 - 2 * parity(sortperm(SVector{length(i) + length(j),Int}(vcat(i, j))))

reverse_sign(grade) = (-1) ^ (grade * (grade - 1) ÷ 2)

reverse(x::GeometricAlgebraType) = x * reverse_sign(grade(x))

function grade_projection(x::Multivector, ::Val{K}) where {K}
    K > log2(length(x)) && return scalar(zero(eltype(x)))
    res = zero(KVector{K}, eltype(x))
    res.coefs .= x.coefs[indices(res)]
    res
end

grade_projection(x::KVector{K}, ::Val{K}) where {K} = x

grade_projection(x::KVector, ::Val{K}) where {K} = zero(KVector{K}, eltype(x))

magnitude2(x) = x ⦿ x

magnitude(x) = sqrt(abs(magnitude2(x)))

inv(x::GeometricAlgebraType) = x / magnitude2(x)

(/)(x::GeometricAlgebraType, y::GeometricAlgebraType) = x * inv(y)
@type_commutative (/)(x::Any, y::GeometricAlgebraType) = x * inv(y)

"""
Return the grade(s) that can be present in the result of an operation.
"""
result_grade(::typeof(⋅), r, s) = abs(r - s)
result_grade(::typeof(∧), r, s) = r + s
result_grade(::typeof(lcontract), r, s) = s - r
result_grade(::typeof(rcontract), r, s) = r - s
result_grade(::typeof(⦿), _, _) = 0
result_grade(::typeof(*), r, s) = result_grade(⋅, r, s):2:result_grade(∧, r, s)


text:
@inline function simulate!(_sol,p::model, x...;
    outputs=isempty(_sol) ? p.opts.outputs : (@views @inbounds _sol.results[end].opts.outputs),
    overwrite_sol::Bool=false,
    sol::Nothing=nothing,
    kw...)

    simulate(p, x...;
        sol = overwrite_sol ? deepcopy(_sol) : _sol,
        outputs = outputs,
        kw...)
end

@inline function get_run(
    I::current,
    V::voltage,
    P::power,
    t0, tf::Q, p::model, sol::solution) where {
    current  <: Union{Number,Symbol,Function,Nothing},
    voltage  <: Union{Number,Symbol,Function,Nothing},
    power    <: Union{Number,Symbol,Function,Nothing},
    Q,
    }

    if !( sum(!(method == Nothing) for method in (current,voltage,power)) == 1 )
        error("Cannot select more than one input")
    end

    if     current ≠ Nothing
        method, input, name = method_I(), I, "I"
    elseif voltage ≠ Nothing
        method, input, name = method_V(), V, "V"
    elseif power ≠ Nothing
        method, input, name = method_P(), P, "P"
    else
        error("Method not supported")
    end

    if input isa Function
        input = redefine_func(input)
    end

    value = [0.0]
    tf = (Q == Nothing ? 1e6 : (@inbounds Float64(tf[end])))

    run_type = run_determination(method, input)
    run      = run_type(input,value,method,t0,tf,name,run_info())

    return run
end

@inline run_determination(::AbstractMethod,::Any)      = run_constant
@inline run_determination(::AbstractMethod,::Function) = run_function

function custom_res!(p::model,x::T,func_RHS::Q,sol::solution;kw...) where {T<:Function,Q<:Function}
    func = redefine_func(x)
    p.θ[:_residual_val] = 0.0
    return (t,Y,YP,p) -> func_RHS(t,Y,YP,p) - func(t,Y,YP,p)
end
function custom_res!(p::model,x::T,func_RHS::Q,sol::solution;kw...) where {T<:Number,Q<:Function}
    p.θ[:_residual_val] = Float64(x)
    return func_RHS
end
function custom_res!(p::model,x::T,func_RHS::Q,sol::solution;hold_val::Number=0.0,kw...) where {T<:Symbol,Q<:Function}
    if check_is_hold(x,sol)
        p.θ[:_residual_val] = hold_val
    end
    return func_RHS
end

@inline within_bounds(run::AbstractRun) = run.info.flag == -1

@inline function initialize_simulation!(sol::model_struct, p::model{jac}, run::T, bounds::boundary_stop_conditions_immutable, opts::AbstractOptionsModel, res_I_guess=nothing) where {jac<:AbstractJacobian,method<:AbstractMethod,T<:AbstractRun{method,<:Any},model_struct<:Union{solution,Vector{Float64}}}
    if !haskey(p.funcs,run)
        get_method_funcs!(p,run)
        funcs = p.funcs(run)
    end
    funcs = p.funcs(run)
    θ_tot = funcs.J_full.θ_tot
    update_θ!(θ_tot, funcs.J_full.θ_keys, p)

    cache = p.cache

    initial_guess! = p.funcs.initial_guess!
    keep_Y = opts.var_keep.Y

    # update the θ_tot vector from the dict p.θ
    check_errors_parameters_runtime(p, opts)

    # if this is a new sol?
    new_run = model_struct == Array{Float64,1} || isempty(sol)

    ## initializing the states vector Y and time t
    Y0 = keep_Y ? zeros(Float64,length(cache.Y0)) : cache.Y0
    YP0 = cache.YP0
    if model_struct == Array{Float64,1}
        Y0 = copy(sol)
        sol = solution()
        SOC = calc_SOC(Y0, p)
    elseif new_run
        SOC = opts.SOC
        if opts.reinit
            initial_guess!(Y0, SOC, θ_tot, res_I_guess)
        end
    else # continue from previous simulation
        Y0 .= @inbounds keep_Y ? copy(sol.Y[end]) : sol.Y[end]
        SOC = @inbounds sol.SOC[end]
    end

    initial_current!(Y0,YP0,p,run,sol,res_I_guess)

    ## getting the DAE integrator function
    initialize_states!(p,Y0,YP0,run,opts,funcs,SOC)

    int = retrieve_integrator(run,p,funcs,Y0,YP0,opts,new_run)

    set_vars!(sol, p, Y0, YP0, int.t, run, opts, bounds; init_all=new_run, SOC=SOC)
    set_var!(sol.Y, Y0, new_run || keep_Y)

    check_simulation_stop!(sol, 0.0, Y0, YP0, run, p, bounds, opts)
    return int, funcs, sol
end

@inline function retrieve_integrator(run::T, p::model, funcs::Jac_and_res{<:Sundials.IDAIntegrator}, Y0, YP0, opts::AbstractOptionsModel, new_run::Bool) where {method<:AbstractMethod,T<:AbstractRun{method,<:Any}}
    """
    If the sol has previously been evaluated for a constant run simulation, you can reuse
    the integrator function with its cache instead of creating a new one
    """

    if isempty(funcs.int) || opts.calc_integrator
        int = create_integrator(run, p, funcs, Y0, YP0, opts)
    else
        # reuse the integrator cache
        int = @inbounds funcs.int[end]
        @inbounds int.p.value .= run.value

        mem = int.mem
        # reinitialize at t = 0 with new Y0/YP0 and tolerances
        Sundials.IDASStolerances(mem, opts.reltol, opts.abstol)
        Sundials.IDAReInit(mem, 0.0, Y0, YP0)
        int.t = 0.0
    end

    postfix_integrator!(int, run, opts, new_run)

    return int
end

@inline function estimate_steady_state(p::model, sol::solution, opts::AbstractOptionsModel = p.opts, run::AbstractRun = (@views @inbounds sol.results[end].run);
    itermax::Int64=100)
    funcs = p.funcs(sol)
    t = @inbounds sol.t[end]
    γ = Float64(1.0)
    Y  = @inbounds copy(sol.Y[end])
    @inbounds Y[p.ind.I[1]] = 0
    YP = p.cache.YP0 .= 0

    R_full = funcs.R_full
    J_full = funcs.J_full
    factor = J_full.factor

    res = p.cache.Y0
    update = similar(res) .= +Inf
    J = J_full.sp

    @inbounds for iter in 1:itermax
        R_full(res,t,Y,YP,p,run)
        J_full(t,Y,YP,γ,p,run)
        factorize!(factor, J)

        res[p.ind.I[1]] = 0
        update .= factor\res

        Y .-= update
        if norm(update) < opts.reltol
            return Y
        elseif iter == itermax
            error("Could not converge to a steady state in $itermax iterations.")
        end
    end
end

@inline function estimate_SOC(p::model, x...;kw...)
    Y = estimate_steady_state(p, x...; kw...)
    SOC = calc_SOC(Y,p)

    return SOC
end

@inline function create_integrator(run::T, p::model, funcs::Q, Y0, YP0, opts::AbstractOptionsModel) where {T<:AbstractRun,Q<:Jac_and_res}
    R_full = funcs.R_full
    J_full = funcs.J_full
    DAEfunc = DAEFunction(
        (res,YP,Y,run,t) -> R_full(res,t,Y,YP,p,run);
        jac = (J,YP,Y,run,γ,t) -> J_full(J,t,Y,YP,γ,p,run),
        jac_prototype = J_full.sp,
    )

    prob = DAEProblem(DAEfunc, YP0, Y0, (0.0, run.tf), run, differential_vars=p.cache.id)

    int = init(prob, Sundials.IDA(;linear_solver=:KLU, init_all=false), tstops=Float64[], abstol=opts.abstol, reltol=opts.reltol, save_everystep=false, save_start=false, verbose=false)

    if isempty(funcs.int)
        push!(funcs.int, int)
    end

    return int
end
@inline function postfix_integrator!(int::Sundials.IDAIntegrator, run::AbstractRun, opts::AbstractOptionsModel, new_run::Bool)
    tstops = int.opts.tstops.valtree
    empty!(tstops)

    append!(tstops,opts.tstops)
    append!(tstops,opts.tdiscon .- opts.reltol/2)
    push!(tstops,run.tf)
    if !new_run prepend!(tstops, 1.0) end

    sort!(tstops)

    # the sol can fail is tstops includes 0
    if (@inbounds tstops[1]) ≤ 0.0 deleteat!(tstops, 1:findfirst(tstops .≤ 0)) end
    return nothing
end

@inline function solve!(sol,int::R1,run::R2,p,bounds,opts::R3,funcs) where {R1<:Sundials.IDAIntegrator,R2<:AbstractRun,R3<:AbstractOptionsModel}
    keep_Y = opts.var_keep.Y
    Y  = int.u.v
    YP = int.du.v

    iter = 1
    status = within_bounds(run)
    @inbounds while status
        step!(int)
        iter += 1
        t = int.t

        set_vars!(sol, p, Y, YP, t, run, opts, bounds)

        check_simulation_stop!(sol, t, Y, YP, run, p, bounds, opts)

        status = check_solve(run, sol, int, p, bounds, opts, funcs, keep_Y, iter, Y, t)
    end

    run.info.iterations = iter
    return nothing
end

@inline function exit_simulation!(p::R1, sol::R2, run::R3, bounds::R4, int::R5, opts::R6; cancel_interp::Bool=false) where {R1<:model,R2<:solution,R3<:AbstractRun,R4<:boundary_stop_conditions_immutable,R5<:Sundials.IDAIntegrator,R6<:AbstractOptionsModel}
    # if a stop condition (besides t = tf) was reached
    if !cancel_interp
        if opts.interp_final && !(run.info.flag == 0) && int.t > 1
            interp_final_points!(p, sol, run, bounds, int, opts)
        else
            set_var!(sol.Y, opts.var_keep.Y ? copy(int.u.v) : int.u.v, opts.var_keep.Y)
        end
    end

    iterations_start = isempty(sol) ? 0 : (@inbounds @views sol.results[end].run_index[end])

    run_index = (1:run.info.iterations) .+ iterations_start

    tspan = (run.t0, @inbounds sol.t[end])

    results = run_results(
        run,
        tspan,
        run.info,
        run_index,
        int,
        opts,
        bounds,
        p.N,
        p.numerics,
        p,
    )

    push!(sol.results, results)

    return nothing
end

@views @inbounds @inline function interp_final_points!(p::R1, sol::R2, run::R3, bounds::R4, int::R5, opts::R6) where {R1<:model,R2<:solution,R3<:AbstractRun,R4<:boundary_stop_conditions_immutable,R5<:Sundials.IDAIntegrator,R6<:AbstractOptionsModel}
    if opts.var_keep.YP
        YP = length(sol.YP) > 1 ? bounds.prev.t_final_interp_frac.*(sol.YP[end] .- sol.YP[end-1]) .+ sol.YP[end-1] : sol.YP[end]
    else
        YP = Float64[]
    end

    t = bounds.prev.t_final_interp_frac*(int.t - int.tprev) + int.tprev

    set_var!(sol.Y,  bounds.prev.t_final_interp_frac.*(int.u.v .- sol.Y[end]) .+ sol.Y[end], opts.var_keep.Y)
    set_vars!(sol, p, sol.Y[end], YP, t, run, opts, bounds; modify! = set_var_last!)

    return nothing
end

@inline function save_start_init!(Y0::Vector{Float64}, run::AbstractRun, p::model, SOC::Float64)
    key = save_start_info(
        run.method,
        round(SOC, digits=4),
        round(value(run), digits=4),
    )

    save_start_dict = p.cache.save_start_dict

    key_exists = key ∈ keys(save_start_dict)
    if key_exists
        @inbounds Y0[(1:p.N.alg) .+ p.N.diff] .= save_start_dict[key]
    end

    return key, key_exists
end
@inline initialize_states!(p::model, Y0::T, YP0::T, run::AbstractRun, opts::AbstractOptionsModel, funcs::Jac_and_res; kw...) where {T<:Vector{Float64}} = newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg;kw...)
@inline function initialize_states!(p::model, Y0::T, YP0::T, run::AbstractRun, opts::AbstractOptionsModel, funcs::Jac_and_res, SOC::Number;kw...) where {T<:Vector{Float64}}
    if opts.save_start
        key, key_exists = save_start_init!(Y0, run, p, SOC)

        newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg;kw...)

        if !key_exists
            p.cache.save_start_dict[key] = @inbounds Y0[(1:p.N.alg) .+ p.N.diff]
        end
    else
        newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg;kw...)
    end

    return nothing
end

@inline factorize!(factor::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64},A::SparseMatrixCSC{Float64, Int64}) = LinearAlgebra.lu!(factor,A)
@inline factorize!(factor::KLUFactorization{Float64, Int64},A::SparseMatrixCSC{Float64, Int64}) = klu!(factor,A)
@inline function newtons_method!(p::model,Y::R1,YP::R1,run,opts::AbstractOptionsModel,R_alg::T1,R_diff::T2,J_alg::T3;
    itermax::Int64=100, t::Float64=0.0
    ) where {R1<:Vector{Float64},T1<:residual_combined,T2<:residual_combined,T3<:jacobian_combined}

    res    = p.cache.res
    Y_old  = p.cache.Y_alg
    Y_new  = @views @inbounds Y[p.N.diff+1:end]
    YP    .= 0.0
    J      = J_alg.sp
    γ      = 0.0
    factor = J_alg.factor

    # starting loop for Newton's method
    @inbounds for iter in 1:itermax
        # updating res, Y, and J
        R_alg(res,t,Y,YP,p,run)
        J_alg(t,Y,YP,γ,p,run)
        factorize!(factor, J)

        Y_old .= Y_new
        Y_new .-= factor\res
        if norm(Y_old .- Y_new) < opts.reltol_init # || maximum(abs, res) < opts.abstol_init
            break
        elseif iter == itermax
            error("Could not initialize DAE in $itermax iterations.")
        end
    end
    # calculate the differential equations for YP0
    R_diff(YP,t,Y,YP,p,run)

    return nothing
end

"""
Current
"""
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_constant{method,in},sol, res_I_guess) where {method<:method_I,in<:Number}
    input = run.input
    @inbounds run.value .= Y0[p.ind.I[1]] = input
    return nothing
end
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_constant{method,in},sol::solution, res_I_guess) where {method<:method_I,in<:Symbol}
    input = run.input
    if check_is_hold(input,sol)
        @inbounds run.value .= Y0[p.ind.I[1]] = calc_I((@views @inbounds sol.Y[end]), p)
    elseif input == :rest
        @inbounds run.value .= Y0[p.ind.I[1]] = 0.0
    else
        error("Unsupported input symbol.")
    end
    return nothing
end
@inline function initial_current!(Y0::Vector{Float64},YP0::Vector{Float64},p,run::run_function{method,func},sol, res_I_guess) where {method<:method_I,func<:Function}
    run.value .= Y0[p.ind.I[1]] = run.func(0.0,Y0,YP0,p)
    return nothing
end

"""
Power
"""
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_constant{method,in},sol, res_I_guess) where {method<:method_P,in<:Number}
    input = run.input
    @inbounds run.value .= Y0[p.ind.I[1]] = input/(calc_V(Y0,p)*p.θ[:I1C])
    return nothing
end
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_constant{method,in},sol::solution, res_I_guess) where {method<:method_P,in<:Symbol}
    input = run.input
    if check_is_hold(input,sol)
        @inbounds run.value .= Y0[p.ind.I[1]] = calc_P((@views @inbounds sol.Y[end]), p)
    elseif input == :rest
        @inbounds run.value .= Y0[p.ind.I[1]] = 0.0
    else
        error("Unsupported input symbol.")
    end
    return nothing
end
@inline function initial_current!(Y0::Vector{Float64},YP0::Vector{Float64},p,run::run_function{method,func},sol, res_I_guess) where {method<:method_P,func<:Function}
    run.value .= Y0[p.ind.I[1]] = run.func(0.0,Y0,YP0,p)/(calc_V(Y0,p)*p.θ[:I1C])
    return nothing
end

"""
Voltage and η_plating
"""
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_constant{method,in},sol::solution, res_I_guess) where {method<:method_V,in<:Number}
    input = run.input
    @inbounds run.value .= input
    if !isempty(sol)
        @inbounds Y0[p.ind.I[1]] = calc_I((@views @inbounds sol.Y[end]), p)
    else
        OCV = calc_V(Y0,p)
        @inbounds Y0[p.ind.I[1]] = input > OCV ? +1.0 : -1.0
    end
    return nothing
end
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_constant{method,in},sol::solution, res_I_guess) where {method<:method_V,in<:Symbol}
    input = run.input
    if check_is_hold(input,sol)
        Y = @views @inbounds sol.Y[end]
        @inbounds run.value .= calc_V(Y, p)
        @inbounds Y0[p.ind.I[1]] = calc_V(Y, p)
    else
        error("Unsupported input symbol.")
    end
    return nothing
end
@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_function{method,func},sol::solution, res_I_guess) where {method<:method_V,func<:Function}
    @inbounds run.value .= run.func(0.0,Y0,YP0,p)
    if !isempty(sol)
        @inbounds Y0[p.ind.I[1]] = calc_I((@views @inbounds sol.Y[end]), p)
    else
        # Arbitrary guess for the initial current.
        OCV = calc_V(Y0,p)
        @inbounds Y0[p.ind.I[1]] = value(run) > OCV ? +1.0 : -1.0
    end
    return nothing
end


text:
function _escape(repl_dict::AbstractDict)
    function escape(s::String)
        for (old, new) in repl_dict
            s = replace.(s, Ref(old => new))
        end
        return s
    end
end

escape_latex_dict = Dict("&" => "\\&", "%" => "\\%", "\$" => "\\\$", "#" => "\\#", "_" => "\\_", "{" => "\\{", "}" => "\\}")

function _escape(s::Symbol)
  if s == :ampersand
    _escape(Dict("&" => "\\&"))
  elseif s == :underscore
    _escape(Dict("_" => "\\_"))
  elseif s == :underscore2space
    _escape(Dict("_" => " "))
  elseif s == :latex
    _escape(escape_latex_dict)
  else
    @error "Please provide `:ampersand`, `:underscore`, `:underscore2space`, `:latex`, a `::Dict` or a custom function."
  end
end


text:
"""
    canonicaldomain([ctype::CanonicalType, ]d::Domain)

Return an associated canonical domain, if any, of the given domain.

For example, the canonical domain of an Interval `[a,b]` is the interval `[-1,1]`.

Optionally, a canonical type argument may specify an alternative canonical domain.
Canonical domains help with establishing equality between domains, with finding
maps between domains and with finding parameterizations.

If a domain implements a canonical domain, it should also implement
`mapfrom_canonical` and `mapto_canonical`.
"""
canonicaldomain(d) = d

"Does the domain have a canonical domain?"
hascanonicaldomain(d) = !(d === canonicaldomain(d))

identitymap(d) = IdentityMap{eltype(d)}(dimension(d))

"Return a map to a domain from its canonical domain."
mapfrom_canonical(d) = identitymap(d)
mapfrom_canonical(d, x) = mapfrom_canonical(d)(x)

"Return a map from the domain to its canonical domain."
mapto_canonical(d) = leftinverse(mapfrom_canonical(d))
mapto_canonical(d, x) = mapto_canonical(d)(x)

"Supertype of kinds of canonical domains."
abstract type CanonicalType end

canonicaldomain(ctype::CanonicalType, d) = d
hascanonicaldomain(ctype::CanonicalType, d) = !(d === canonicaldomain(ctype, d))

mapfrom_canonical(ctype::CanonicalType, d, x) = mapfrom_canonical(ctype, d)(x)
mapto_canonical(ctype::CanonicalType, d, x) = mapto_canonical(ctype, d)(x)

"A canonical domain that is equal but simpler (e.g. a 1-dimensional ball is an interval)."
struct Equal <: CanonicalType end

canonicaldomain(::Equal, d) = d
mapfrom_canonical(::Equal, d) = identitymap(d)
mapto_canonical(::Equal, d) = leftinverse(mapto_canonical(Equal(), d))

simplify(d) = canonicaldomain(Equal(), d)
simplifies(d) = hascanonicaldomain(Equal(), d)

"A canonical domain that is isomorphic but may have different element type."
struct Isomorphic <: CanonicalType end

canonicaldomain(::Isomorphic, d) = canonicaldomain(Equal(), d)
mapfrom_canonical(::Isomorphic, d) = mapfrom_canonical(Equal(), d)
mapto_canonical(::Isomorphic, d) = leftinverse(mapto_canonical(Isomorphic(), d))

canonicaldomain(::Isomorphic, d::Domain{SVector{1,T}}) where {T} =
    convert(Domain{T}, d)
mapfrom_canonical(::Isomorphic, d::Domain{SVector{1,T}}) where {T} =
    NumberToVector{T}()

canonicaldomain(::Isomorphic, d::Domain{NTuple{N,T}}) where {N,T} =
    convert(Domain{SVector{N,T}}, d)
mapfrom_canonical(::Isomorphic, d::Domain{NTuple{N,T}}) where {N,T} =
    VectorToTuple{N,T}()

"A parameter domain that can be mapped to the domain."
struct Parameterization <: CanonicalType end

canonicaldomain(ctype::Parameterization, d) =
    hascanonicaldomain(d) ? canonicaldomain(ctype, canonicaldomain(d)) : d
mapfrom_canonical(ctype::Parameterization, d) =
    hascanonicaldomain(d) ? mapfrom_canonical(d) ∘ mapfrom_canonical(ctype, canonicaldomain(d)) : mapfrom_canonical(d)
mapto_canonical(ctype::Parameterization, d) = leftinverse(mapfrom_canonical(ctype, d))

# We define some convenience functions:
"Return a parameter domain which supports a `parameterization`."
parameterdomain(d::Domain) = canonicaldomain(Parameterization(), d)

"Return a parameterization of the given domain."
parameterization(d::Domain) = mapfrom_canonical(Parameterization(), d)

"Does the domain have a parameterization?"
hasparameterization(d) = hascanonicaldomain(Parameterization(), d)

mapfrom_parameterdomain(d::Domain) = mapfrom_canonical(Parameterization(), d)
mapto_parameterdomain(d::Domain) = mapto_canonical(Parameterization(), d)
mapfrom_parameterdomain(d::Domain, x) = mapfrom_canonical(Parameterization(), d, x)
mapto_parameterdomain(d::Domain, x) = mapto_canonical(Parameterization(), d, x)

"Return a map from domain `d1` to domain `d2`."
mapto(d1, d2) = mapto1(d1, d2)
mapto(d1::D, d2::D) where {D} = d1 == d2 ? identitymap(d1) : mapto1(d1,d2)

# simplify the first argument
mapto1(d1, d2) =
    hasparameterization(d1) ? mapto(parameterdomain(d1), d2) ∘ mapto_parameterdomain(d1) : mapto2(d1, d2)
# simplify the second argument
mapto2(d1, d2) =
    hasparameterization(d2) ? mapfrom_parameterdomain(d2) ∘ mapto(d1, parameterdomain(d2)) : no_known_mapto(d1,d2)

no_known_mapto(d1, d2) = d1 == d2 ? identitymap(d1) : error("No map known between $(d1) and $(d2).")

## Equality for domains

## Note: here we generically define the equality of Domains, using the framework
# of canonical domains above. It has the benefit of automatically recognizing
# equality between some domains, but a side-effect is that the default
# (implicitly defined) equality for a concrete domain type is overruled.
# To be safe, concrete domains should specialize `==`.
==(d1::Domain, d2::Domain) = isequal1(d1, d2)
# simplify the first argument
isequal1(d1, d2) = simplifies(d1) ? simplify(d1)==d2 : isequal2(d1, d2)
# simplify the second argument
isequal2(d1, d2) = simplifies(d2) ? d1==simplify(d2) : d1===d2


text:
"""
_compute_mutliplier: Additional function to estimate Lagrange multiplier of the problems
    (guarantee if LICQ holds)

`_compute_mutliplier(pb :: AbstractNLPModel, x :: T, gx :: T, cx :: T, Jx :: MT; active_prec_c :: Float64 = 1e-6, active_prec_b :: Float64 = 1e-6)`
"""
function _compute_mutliplier(pb    :: AbstractNLPModel,
                             x     :: T,
                             gx    :: T,
                             cx    :: T,
                             Jx    :: MT;
                             active_prec_c :: AbstractFloat = 1e-6,
                			       active_prec_b :: AbstractFloat = 1e-6
                            ) where {MT, T}

  n  = length(x)
  nc = cx == _init_field(typeof(cx)) ? 0 : length(cx)

  #active res_bounds
  Ib = findall(x->(norm(x) <= active_prec_b),
	  		      min(abs.(x - pb.meta.lvar),
				      abs.(x - pb.meta.uvar)))
  if nc != 0
    #active constraints
    Ic = findall(x->(norm(x) <= active_prec_c),
                     min(abs.(cx-pb.meta.ucon),
                     abs.(cx-pb.meta.lcon)))

    Jc = hcat(Matrix(1.0I, n, n)[:,Ib], Jx'[:,Ic])
  else
    Ic = []
    Jc = hcat(Matrix(1.0I, n, n)[:,Ib])
  end

  mu, lambda = zeros(n), zeros(nc)
  if (Ib != []) || (Ic  != [])
    l = Jc \ (- gx)
    mu[Ib], lambda[Ic] = l[1:length(Ib)], l[length(Ib)+1:length(l)]
  end

  return mu, lambda
end


text:
#######################################################################
# Utilities for parsing Julia files
#######################################################################

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Parsing Julia files
#
#   Based on
#   https://discourse.julialang.org/t/parsing-a-julia-file/32622
#
#=
parsefile(file) = parse(join(["quote", readstring(file), "end"], ";"))
parsecode(code::String)::Vector =
    # https://discourse.julialang.org/t/parsing-a-julia-file/32622
    filter(x->!(x isa LineNumberNode),
           Meta.parse(join(["quote", code, "end"], ";")).args[1].args)
=#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# String → AST
# Parses `text` as Julia code
parseJuliaCode(text :: String) =
    Meta.parse(join(["quote", text, "end"], "\n"))

# String → AST
# Parses file `filePath` as Julia code
parseJuliaFile(filePath :: String) =
    parseJuliaCode(read(filePath, String))

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
using Diderot: Arc, Node, RestrictLowDistance
using Diderot.IndexFund

const IF = IndexFund

function make_instance(number)
    weights = [0.4, 0.2, 0.2, 0.1, 0.1]
    similarity = [
        1.0    0.1   -0.1    0.08   0.51
        0.1    1.0    0.07  -0.61   0.55
       -0.1    0.07   1.0    0.27   0.46
        0.08  -0.61   0.27   1.0   -0.3
        0.51   0.55   0.46  -0.3    1.0
    ]
    return IF.Instance(weights, similarity, number)
end

@testset "model" begin
    instance = make_instance(2)

    @test Diderot.domain_type(instance) == Int
    @test Diderot.value_type(instance) == Float64
    @test Diderot.initial_state(instance)== BitSet()
    @test length(instance) == 5

    @test Diderot.transitions(instance, BitSet(), 1) ==
        Dict(Arc(BitSet(), 1, 0.4 * 1.0)    => BitSet(1),
             Arc(BitSet(), 2, 0.4 * 0.1)    => BitSet(2),
             Arc(BitSet(), 3, 0.4 * (-0.1)) => BitSet(3),
             Arc(BitSet(), 4, 0.4 * 0.08)   => BitSet(4),
             Arc(BitSet(), 5, 0.4 * 0.51)   => BitSet(5))
    @test Diderot.transitions(instance, BitSet(1), 1) ==
        Dict(Arc(BitSet(1), 1, 0.4 * 1.0)    => BitSet(1),
             Arc(BitSet(1), 2, 0.4 * 0.1)    => BitSet(1:2),
             Arc(BitSet(1), 3, 0.4 * (-0.1)) => BitSet([1,3]),
             Arc(BitSet(1), 4, 0.4 * 0.08)   => BitSet([1,4]),
             Arc(BitSet(1), 5, 0.4 * 0.51)   => BitSet([1,5]))
    @test Diderot.transitions(instance, BitSet(1:2), 1) ==
        Dict(Arc(BitSet(1:2), 1, 0.4 * 1.0)    => BitSet(1:2))
end

@testset "exact diagram" begin
    N = Node{BitSet,Int,Float64}

    instance = make_instance(2)
    diagram = Diderot.Diagram(instance)
    Diderot.top_down!(diagram, instance)

    @test length(diagram.layers) == 6

    @test length(diagram.layers[1]) == 1
    @test diagram.layers[1][BitSet()] == N(0.0, nothing)

    @test length(diagram.layers[2]) == 5
    @test diagram.layers[2][BitSet(1)] ≈ N(0.4, Arc(BitSet(), 1, 0.4))
    @test diagram.layers[2][BitSet(2)] ≈ N(0.04, Arc(BitSet(), 2, 0.04))
    # ...

    @test length(diagram.layers[end]) == 1

    solution = Diderot.longest_path(diagram)
    @test solution.decisions == [1, 5, 5, 1, 5]
    @test solution.objective ≈ 0.71
end

@testset "restriction" begin
    instance = make_instance(2)

    @testset "width 1" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=RestrictLowDistance(1))
        @test length(diagram.layers) == 6
        @test all(l -> length(l) == 1, diagram.layers)
        solution = Diderot.longest_path(diagram)
        @test solution.objective <= 0.72
    end

    @testset "width 2" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=RestrictLowDistance(2))
        @test length(diagram.layers) == 6
        @test all(l -> length(l) <= 2, diagram.layers)
        solution = Diderot.longest_path(diagram)
        @test solution.objective <= 0.72
    end
end

@testset "relaxation" begin
    instance = make_instance(2)

    @testset "width 1" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=IF.relax(1))
        @test length(diagram.layers) == 6
        @test all(l -> length(l) == 1, diagram.layers)
        solution = Diderot.longest_path(diagram)
        @test solution.objective >= 0.7
    end

    @testset "width 2" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=IF.relax(2))
        @test length(diagram.layers) == 6
        @test all(l -> length(l) <= 2, diagram.layers)
        solution = Diderot.longest_path(diagram)
        @test solution.objective >= 0.7
    end
end

@testset "branch and bound" begin
    instance = make_instance(2)

    # need at least width 5 for exact cutsets to branch on
    @testset "width 5" begin
        solution = Diderot.branch_and_bound(
            instance, restrict=RestrictLowDistance(5), relax=IF.relax(5))
        @test solution.decisions == [1, 5, 5, 1, 5]
        @test solution.objective ≈ 0.71
    end

    @testset "width 10" begin
        solution = Diderot.branch_and_bound(
            instance, restrict=RestrictLowDistance(10), relax=IF.relax(10))
        @test solution.decisions == [1, 5, 5, 1, 5]
        @test solution.objective ≈ 0.71
    end
end


text:
module ReactiveAmqp

include("includes.jl")

export AmqpConnectionDef, compose_auth, amqp_conn_define!
export source!, queue!, execute_queues!
export map_ack, safe_map_ack!, sink!
export init_config!

end


text:
using Documenter, DiscretePredictors

makedocs(
    modules = [DiscretePredictors],
    doctest = true,
    format = Documenter.HTML(prettyurls=!("local" in ARGS)),
    sitename = "DiscretePredictors.jl",
    authors = "Vishnu Raj",
    linkcheck = false,
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "man/examples.md",
            # "man/contributing.md",
        ],
        "Library" => Any[
            "High Level Interface" => "lib/interface.md",
            "Predictors" => "lib/predictors.md"
        ]
    ],
    # Use clean URLs, unless built as a "local" build
    # html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo    = "github.com/v-i-s-h/DiscretePredictors.jl.git",
    target  = "build",
    # julia   = "0.6",
    # osname  = 'linux',
    deps    = nothing,
    make    = nothing,
)


text:
#reload("MetabolicFluxAnalysis.jl")
using MetabolicFluxAnalysis
using CoordinateTransformations

# current tests..
fluxi_1_f = open("testdata/test_01.fluxi")
fluxi_1_s = readstring(fluxi_1_f)
fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

# now get and then parametrize the flux polytope
fp = MetabolicFluxAnalysis.get_flux_polyhedron( fluxi1.net , fluxi1.linear_constraints["ConstraintsExchange"] , max_val=200 )

# compute transform to make full-dimensional:
( fp_fd , at_fd ) = MetabolicFluxAnalysis.find_minimal_affine_subspace(fp)
fp_fd_a = MetabolicFluxAnalysis.remove_redundant_constraints( fp_fd )
( fp_rounded , at_rounded ) = MetabolicFluxAnalysis.round_polytope(fp_fd_a)

at_opt2flux = CoordinateTransformations.compose(at_fd , at_rounded)

# test chebyshev comp.
x_int_cheby = MetabolicFluxAnalysis.solve_chebyshev(fp_fd.G,fp_fd.h[:])

@test all( fp_fd.G * x_int_cheby .< fp_fd.h )
print( fp_fd.h - fp_fd.G * x_int_cheby  )

## test with toy polytopes..
d_test = 100
g_test = [eye(d_test);-eye(d_test)]; h_test = 10*rand(2*d_test,1);
a_test = rand(2,d_test)       ; b_test = 0.1 * rand(2,1);
f_test = MetabolicFluxAnalysis.Polyhedron(g_test,h_test,a_test,b_test)
rfd_parametrization = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( f_test ; remove_red_constr=1)
@assert MetabolicFluxAnalysis.isinside( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt[:,:] )

@test MetabolicFluxAnalysis.isinside( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt[:,:] )

# generate some unifrom samples over polytope:
xx = MetabolicFluxAnalysis.sample_unif_hr( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt , 1000 ,1000 );

xx_raw = MetabolicFluxAnalysis.eval_opt2raw(rfd_parametrization,xx)

@test MetabolicFluxAnalysis.checksamples( rfd_parametrization.p_opt , xx     ; eq_tol=1e-4)
@test MetabolicFluxAnalysis.checksamples( rfd_parametrization.p_raw , xx_raw ; eq_tol=1e-4)

##

# some more basic tests
# check that things behave correctly for unbounded polyhedra
using MetabolicFluxAnalysis
bb_01 = MetabolicFluxAnalysis.bounding_box([-10.;-10;-10],[10.;10;20])
unbounded_polyhedron = MetabolicFluxAnalysis.Polyhedron( bb_01.G[1:1,:],bb_01.h[1:1,:],bb_01.A[:,:],bb_01.b[:,:] )
lpm_result = MetabolicFluxAnalysis.solve_lp_max_v(unbounded_polyhedron,[0.;1;0])

#(mve_f,mve_E) =  MetabolicFluxAnalysis.mve_run(fp_fd_a.G,fp_fd_a.h[:])

# test rounding transform:
cube_g = [1.0*eye(4);-eye(4)]
cube_h = 1.0*collect(1:8)

p_c = MetabolicFluxAnalysis.Polyhedron(cube_g,cube_h[:,:])

# remove_redundant_constraintss
param1= fp_fd_a.G
param2 = fp_fd_a.h[:]
param3 = [79.514 ;52.79]
(mve_f,mve_E) =  MetabolicFluxAnalysis.mve_solver( param1,param2,param3 )

(mve_f,mve_E) =  MetabolicFluxAnalysis.mve_run(cube_g,cube_h)

# transform somehow:

p_fixed = MetabolicFluxAnalysis.Polyhedron(fp_fd,MetabolicFluxAnalysis.bounding_box([-1000.;-1000],[1000.;1000]) )

(mve_f,mve_E) = MetabolicFluxAnalysis.mve_run( p_fixed.G , p_fixed.h[:] )
(mve_f,mve_E) = MetabolicFluxAnalysis.mve_run( fp_fd.G , fp_fd.h[:] )

# compute rounding transform:
(mve_f,mve_E) = MetabolicFluxAnalysis.mve_run( fp_fd.G , fp_fd.h[:] )


text:
function add_born_sources!(issp::Int64, pac::PFdtdc, pass::Vector{PFdtdss}, pap::PFdtdp)

	δx24I=pac.δx24I; δz24I=pac.δz24I;
	δxI=pac.δxI; δzI=pac.δzI;
	δt=pac.δt
	δtI=pac.δtI
	δmodtt=pac.δmodtt; modttI=pac.modttI;
	δmodrrvx=pac.δmodrrvx; δmodrrvz=pac.δmodrrvz

	nx=pac.nx; nz=pac.nz

	born_svalue_stack=pap.born_svalue_stack
	p=pap.p; pp=pap.pp; ppp=pap.ppp;
	dpdx=pap.dpdx; dpdz=pap.dpdz;

	# secondary sources for Born modeling
	# adding born sources from pressure(:,:,1) to pressure(:,:,2)
	# upto until [it-2]
	# lambdaI scatterrer source term at [it-1]
	# p is at [it], pp is at [it-1], ppp is at [it-2]
	# dpdx is at [it-1] and dpdz is at [it-1]
	# modrrvx scatterrer source term at [it-1]
	# modrrvz scatterrer source term at [it-1]

	born_stacktt!(born_svalue_stack,p,pp,ppp,δmodtt,nx,nz,δtI,δt)
	born_stackrrvx!(born_svalue_stack,dpdx,δmodrrvx,nx,nz,δx24I,δt)
	born_stackrrvz!(born_svalue_stack,dpdz,δmodrrvz,nx,nz,δz24I,δt)
	born_stack!(p,born_svalue_stack,modttI,nx,nz,δt)
end
@inbounds @fastmath function born_stacktt!(born_svalue_stack,p,pp,ppp,δmodtt,nx,nz,δtI,δt)
	pppw=ppp[1]
	ppw=pp[1]
	pw=p[1]
	for ix=1:nx
		@simd for iz=1:nz
			born_svalue_stack[iz,ix] +=
			δt * ((-1.0 * (pppw[iz, ix, 1] + pw[iz, ix,  1]
		   	- 2.0 * ppw[iz, ix,  1]) * δmodtt[iz, ix] * δtI * δtI))
		end
	end
end
@inbounds @fastmath function born_stackrrvx!(born_svalue_stack,dpdx,δmodrrvx,nx,nz,δx24I,δt)
	dpdxw=dpdx[1]
	for ix=3:nx-1
		@simd for iz=1:nz
			born_svalue_stack[iz,ix] +=
				δt * ((27.e0*dpdxw[iz,ix,1] * δmodrrvx[iz,ix] -27.e0*dpdxw[iz,ix-1,1] * δmodrrvx[iz,ix-1] -dpdxw[iz,ix+1,1] * δmodrrvx[iz,ix+1] +dpdxw[iz,ix-2,1] * δmodrrvx[iz,ix-2] ) * (δx24I))
		end
	end

end
@inbounds @fastmath function born_stackrrvz!(born_svalue_stack,dpdz,δmodrrvz,nx,nz,δz24I,δt)
	dpdzw=dpdz[1]
	for ix=1:nx
		@simd for iz=3:nz-1
			born_svalue_stack[iz,ix] +=
				δt * ((27.e0*dpdzw[iz,ix,1] * δmodrrvz[iz,ix] -27.e0*dpdzw[iz-1,ix,1] * δmodrrvz[iz-1,ix] -dpdzw[iz+1,ix,1] * δmodrrvz[iz+1,ix] +dpdzw[iz-2,ix,1] * δmodrrvz[iz-2,ix] ) * (δz24I))

		end
	end
end
@inbounds @fastmath function born_stack!(p,born_svalue_stack,modttI,nx,nz,δt)
	pw2=p[2]
	for ix=1:nx
		@simd for iz=1:nz
			pw2[iz,ix,1] += born_svalue_stack[iz,ix] * modttI[iz,ix] * δt  #* δxI * δzI
		end
	end
end


text:
function ReadCSV(fname)
    inDim = readdlm(fname,',')
    (rowN,colN) = size(inDim)
    row = [];
    col = [];
    val = [];
    if colN >2
        for x = 2:colN
            for y = 2:rowN
                push!(row, inDim[y,1])
                push!(col, inDim[1,x])
                push!(val, inDim[y,x])
            end
        end
        return Assoc(row,col,val,(+))
    else #2D associative?
    #Currently Ignoring
        return
    end
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
import MacroTools
using TraceCalls: apply_code!, revert_code!, update_code_revertible, RevertibleCodeUpdate,
    CodeUpdate, source

# a = parse_file("docstring.jl")[1] # TODO: see comment in that file
# @test TraceCalls.strip_docstring(a).head == :function

################################################################################
# RevertibleCodeUpdate

# That's a hacky way of reusing the test file in ClobberingReload
push!(LOAD_PATH, joinpath(@__DIR__, "ClobberingReloadTests"))
using AA

include("incl.jl")

counter = fill(0)
function add_counter(fdef)
    di = MacroTools.splitdef(fdef)
    di[:body] = quote $counter .+= 1; $(di[:body]) end
    MacroTools.combinedef(di)
end
upd_high = update_code_revertible(AA.high) do code
    add_counter(code)
end
upd_module = update_code_revertible(AA) do code
    if TraceCalls.is_function_definition(code) add_counter(code) end
end
upd_include = update_code_revertible("incl.jl") do code
    if TraceCalls.is_function_definition(code) add_counter(code) end
end
@test AA.high(1) == 10
@test counter[] == 0
upd_high() do
    @test AA.high(1) == 10
    @test AA.high(1.0) == 2
end
@test AA.high(1) == 10
@test AA.bar(1.0) == 2
@test counter[] == 2 # only the two calls within `upd` increase the counter
upd_module() do
    @test AA.high(1) == 10
    @test AA.high(1.0) == 2
end
@test counter[] == 5 # three calls, since `bar` also becomes counting
upd_include() do
    @test apple() == :orange
    @test plant() == "happy"
    dog()
    inlined()
end
apple()
@test counter[] == 9
upd_include() do
    # The number of counter hits here depends on the SimpleTraits macroexpansion.
    # Right now, each call has four hits - September'17
    f(1.0)
    f(0.0)
end
@test counter[] == 17
## upd_include isn't written to test generated functions. Could be fixed, but I'm testing
## them already with `@trace` in runtests.jl
# @test generated_mouse("hey") == "hey DataType"
# upd_include() do
#     @test generated_mouse("hey") == "hey DataType"
# end
# @test generated_mouse("hey") == "hey DataType"
# @test counter[] == 18

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

source(Base.which) # check that it works with Base
source(Base.vcat)  # nested where (as of July '17)


text:
function set_initQbase(xmax,ymax,restart_file,init_rho,init_u,init_v,init_p,specific_heat_ratio,out_file_front,out_ext,out_dir,restart_num)
    Qbase=[]
    cellxmax = xmax - 1
    cellymax = ymax - 1

    restart_check=0
    try Qbase=setup_restart_value(cellxmax,cellymax,out_dir,restart_file)
        println("Restart "*restart_file)
        restart_check=2
    catch
        restart_check=1
    end

    if restart_check == 1
        Qbase=setup_init_value(cellxmax,cellymax,init_rho,init_u,init_v,init_p)
        println("Start Initial condition")
        restart_num=0
        output_result(0,Qbase,cellxmax,cellymax,specific_heat_ratio,out_file_front,out_ext,out_dir)
    end

    return Qbase,cellxmax,cellymax,restart_num
end

function setup_init_value(cellxmax,cellymax,init_rho,init_u,init_v,init_p)
    Qbase=zeros(cellxmax, cellymax, 4)
    for i in 1:cellxmax
        for j in 1:cellymax
            Qbase[i,j,1]=init_rho
            Qbase[i,j,2]=init_u
            Qbase[i,j,3]=init_v
            Qbase[i,j,4]=init_p
        end
    end
    return Qbase
end

function setup_restart_value(cellxmax,cellymax,out_dir,restart_file)
    Qbase=zeros(cellxmax,cellymax,4)

    skipnum=1
    fff=[]
    open("result/"*restart_file, "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    k=1
    for i in 2:cellxmax-1
        for j in 2:cellymax-1
            temp = split(fff[k+skipnum]," ")
            for k in 1:4
                Qbase[i,j,k]=parse(Float64,temp[k])
            end
            k = k+1
        end
    end
    return Qbase
end

function base_to_conservative(Qbase,cellxmax,cellymax,specific_heat_ratio)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """
    Qcon = zeros(cellxmax,cellymax,4)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            Qcon[i,j,1]=Qbase[i,j,1]
            Qcon[i,j,2]=Qbase[i,j,1]*Qbase[i,j,2]
            Qcon[i,j,3]=Qbase[i,j,1]*Qbase[i,j,3]
            Qcon[i,j,4]=Qbase[i,j,4]/(specific_heat_ratio-1)+Qbase[i,j,1]*(Qbase[i,j,2]^2+Qbase[i,j,3]^2)/2
        end
    end
    return Qcon
end

function conservative_to_base(Qcon,cellxmax,cellymax,specific_heat_ratio)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """
    Qbase = zeros(cellxmax,cellymax,4)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            Qbase[i,j,1]=Qcon[i,j,1]
            Qbase[i,j,2]=Qcon[i,j,2]/Qcon[i,j,1]
            Qbase[i,j,3]=Qcon[i,j,3]/Qcon[i,j,1]
            Qbase[i,j,4]=(Qcon[i,j,4]-Qcon[i,j,1]*(Qbase[i,j,2]^2+Qbase[i,j,3]^2)/2)*(specific_heat_ratio-1)
        end
    end
    return Qbase
end


text:
__precompile__(true)
"""
StrAPI package

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module StrAPI

using ModuleInterfaceTools
using ModuleInterfaceTools: m_eval, _stdout, _stderr

const NEW_ITERATE = VERSION >= v"0.7.0-DEV.5127"

const MaybeSub{T} = Union{T, SubString{T}} where {T<:AbstractString}

const CodeUnitTypes = Union{UInt8, UInt16, UInt32}

symstr(s...)   = Symbol(string(s...))
quotesym(s...) = Expr(:quote, symstr(s...))

@api public! found, find_result, basetype, charset, encoding, cse, codepoints

@api public StringError

@api develop NEW_ITERATE, CodeUnitTypes, CodePoints, MaybeSub, symstr, quotesym,
             _stdout, "@preserve"

@api base convert, getindex, length, map, collect, hash, sizeof, size, strides,
          pointer, unsafe_load, string, read, write, reverse,
          nextind, prevind, typemin, typemax, rem, size, ndims, first, last, eltype,
          isless, -, +, *, ^, cmp, promote_rule, one, repeat, filter,
          print, show, isimmutable, chop, chomp, replace, ascii, uppercase, lowercase,
          lstrip, rstrip, strip, lpad, rpad, split, rsplit, join, IOBuffer,
          containsnul, unsafe_convert, cconvert

# Conditionally import or export names that are only in v0.6 or in master
@api base! codeunit, codeunits, ncodeunits, codepoint, thisind, firstindex, lastindex

@static NEW_ITERATE ? (@api base iterate) : (@api base start, next, done)

@static if V6_COMPAT
    include("compat.jl")
else # !V6_COMPAT
    import Base.GC: @preserve

    function find end
    function ind2chr end
    function chr2ind end

    # Handle changes in array allocation
    create_vector(T, len)  = Vector{T}(undef, len)

    # Add new short name for deprecated hex function
    outhex(v, p=1) = string(v, base=16, pad=p)

    get_iobuffer(siz) = IOBuffer(sizehint=siz)

    const utf8crc         = Base._crc32c
    const is_lowercase    = islowercase
    const is_uppercase    = isuppercase
    const lowercase_first = lowercasefirst
    const uppercase_first = uppercasefirst

    using Base: unsafe_crc32c, Fix2

    # Location of some methods moved from Base.UTF8proc to Base.Unicode
    const UC = Base.Unicode
    const Unicode = UC

    import Base.CodeUnits

    @api base IteratorSize

    const is_letter = isletter

    pwc(c, io, str) = printstyled(io, str; color = c)

    const graphemes = UC.graphemes

end # !V6_COMPAT

@api base isequal, ==, in

pwc(c, l) = pwc(c, _stdout(), l)

pr_ul(l)     = pwc(:underline, l)
pr_ul(io, l) = pwc(:underline, io, l)

@api develop! pwc, pr_ul

const str_next = @static NEW_ITERATE ? iterate : next
str_done(str::AbstractString, i::Integer) = i > ncodeunits(str)

@api develop! str_next, str_done
@api develop unsafe_crc32c, Fix2, CodeUnits
@api public! is_lowercase, is_uppercase, lowercase_first, uppercase_first

function found end
function find_result end

"""Get the base type (of CodeUnitTypes) of a character or aligned/swapped type"""
function basetype end

"""Get the character set used by a string or character type"""
function charset end

"""Get the character set used by a string type"""
function encoding end

"""Get the character set / encoding used by a string type"""
function cse end

function _write end
function _print end
function _isvalid end
function _lowercase end
function _uppercase end
function _titlecase end
@api develop! _write, _print, _isvalid, _lowercase, _uppercase, _titlecase

const curmod = @static V6_COMPAT ? current_module() : @__MODULE__

include("errors.jl")
include("traits.jl")
include("codepoints.jl")
include("uni.jl")

@api modules Uni, StrErrors, Unicode

@api develop boundserr, strerror, nulerr, neginderr, codepoint_error,
             argerror, ascii_err, ncharerr, repeaterr, index_error

# Possibly import functions, give new names with underscores

# Todo: Should probably have a @api function for importing/defining renamed functions
const namlst = Symbol[]

for (pref, lst) in
    ((0,
      ((:textwidth,      :text_width),
       (:occursin,       :occurs_in),
       (:startswith,     :starts_with),
       (:endswith,       :ends_with))),

     (1,
      ((:xdigit, :hex_digit),
       (:cntrl,  :control),
       (:punct,  :punctuation),
       (:print,  :printable))),

     (2, (:ascii, :digit, :space, :numeric,
          :valid, :defined, :assigned, :empty,
          :latin, :bmp, :unicode))
     ), nam in lst

    oldname, newname =
        (pref == 0 ? nam : pref == 1
         ? (symstr("is", nam[1]), symstr("is_", nam[2]))
         : (symstr("is", nam), symstr("is_", nam)))

    m_eval(curmod,
           (isdefined(Base, oldname)
            ? Expr(:const, Expr(:(=), newname, oldname))
            : Expr(:function, newname)))

    push!(namlst, newname)
end
@eval @api public! $(namlst...)

# Handle renames where function was deprecated

# Todo: have function for defining and making public
function is_alphabetic end
function is_alphanumeric end
function is_graphic end
@api public! is_alphabetic, is_alphanumeric, is_graphic, is_letter

# import and add new names from UTF8proc/Unicode

const is_grapheme_break  = UC.isgraphemebreak
const is_grapheme_break! = UC.isgraphemebreak!
const category_code      = UC.category_code
const category_abbrev    = UC.category_abbrev
const category_string    = UC.category_string

@api public! is_grapheme_break, is_grapheme_break!, category_code, category_abbrev, category_string

const fnd = find
@api public! fnd, find

@api develop create_vector, outhex, get_iobuffer
@api develop! utf8crc, ind2chr, chr2ind

# Operations for find/search operations

abstract type FindOp end

struct First <: FindOp end
struct Last  <: FindOp end
struct Next  <: FindOp end
struct Prev  <: FindOp end
struct Each  <: FindOp end
struct All   <: FindOp end

abstract type Direction <: FindOp end

struct Fwd   <: Direction end
struct Rev   <: Direction end

@api public FindOp, Direction, Fwd, Rev, First, Last, Next, Prev, Each, All

@api freeze

end # module StrAPI


text:
@deprecate fadj(g) adjacency_list(g)
@deprecate fadj(g, v) out_neighbors(g, v)
@deprecate badj(g) adjacency_list(g, :in)
@deprecate badj(g, v) in_neighbors(g, v)
@deprecate induced_subgraph subgraph
@deprecate indegree(g, v) in_degree(g, v)
@deprecate outdegree(g, v) out_degree(g, v)

@deprecate radius(vec) minimum(vec)
@deprecate diameter(vec) maximum(vec)

@deprecate adjacency_spectrum(g)  eigvals(Matrix(adjacency_matrix(g)))
@deprecate laplacian_spectrum(g)  eigvals(Matrix(laplacian_matrix(g)))

@deprecate is_cyclic(g) has_cycles(g)


text:
import ForwardDiff

using CCBlade: residual, Outputs, Rotor, Section, OperatingPoint

struct PartialsWrt{TF}
    phi::TF

    r::TF
    chord::TF
    theta::TF

    Vx::TF
    Vy::TF
    rho::TF
    mu::TF
    asound::TF

    Rhub::TF
    Rtip::TF
    pitch::TF
    precone::TF

end

# Need this for the mapslices call in output_partials.
PartialsWrt(x::AbstractArray) = PartialsWrt(x...)

# Convenience function to access fields within an array of structs, from Andrew
# Ning.
Base.getproperty(obj::Array{PartialsWrt{TF}, N}, sym::Symbol) where {TF, N} = getfield.(obj, sym)

function pack_inputs(phi, rotor, section, op)
    r = section.r
    chord = section.chord
    theta = section.theta
    Rhub = rotor.Rhub
    Rtip = rotor.Rtip
    pitch = rotor.pitch
    precone = rotor.precone
    Vx = op.Vx
    Vy = op.Vy
    rho = op.rho
    mu = op.mu
    asound = op.asound
    x = [phi, r, chord, theta, Vx, Vy, rho, mu, asound, Rhub, Rtip, pitch, precone]

    return x
end

function unpack_inputs(x, af, B, turbine)
    phi = x[1]
    r = x[2]
    chord = x[3]
    theta = x[4]
    Vx, Vy, rho, mu, asound = x[5], x[6], x[7], x[8], x[9]
    Rhub, Rtip, pitch, precone = x[10], x[11], x[12], x[13]
    rotor = Rotor(Rhub, Rtip, B, turbine, pitch, precone)
    section = Section(r, chord, theta, af)
    op = OperatingPoint(Vx, Vy, rho, mu, asound)
    return phi, rotor, section, op
end

function residual_partials(phi, rotor, section, op)
    af = section.af
    B = rotor.B
    turbine = rotor.turbine

    # Get a version of the residual function that's compatible with ForwardDiff.
    function R(inputs)
        _phi, _rotor, _section, _op = unpack_inputs(inputs, af, B, turbine)
        res, out = residual(_phi, _rotor, _section, _op)
        return res
    end

    # Do it.
    x = pack_inputs(phi, rotor, section, op)
    return PartialsWrt(ForwardDiff.gradient(R, x))

end

function output_partials(phi, rotor, section, op)
    af = section.af
    B = rotor.B
    turbine = rotor.turbine

    # Get a version of the output function that's compatible with ForwardDiff.
    function R(inputs)
        _phi, _rotor, _section, _op = unpack_inputs(inputs, af, B, turbine)
        res, out = residual(_phi, _rotor, _section, _op)
        return [getfield(out, i) for i in fieldnames(typeof(out))]
    end

    # Do it.
    x = pack_inputs(phi, rotor, section, op)
    return Outputs(mapslices(PartialsWrt, ForwardDiff.jacobian(R, x), dims=2)...)
end