text:
#=
Single-point derivatives of scalar->scalar maps.
=#
function finite_difference_derivative(f, x::T, fdtype::Type{T1}=Val{:central},
    returntype::Type{T2}=eltype(x), f_x::Union{Nothing,T}=nothing) where {T<:Number,T1,T2}

    epsilon = compute_epsilon(fdtype, x)
    if fdtype==Val{:forward}
        return (f(x+epsilon) - f(x)) / epsilon
    elseif fdtype==Val{:central}
        return (f(x+epsilon) - f(x-epsilon)) / (2*epsilon)
    elseif fdtype==Val{:complex} && returntype<:Real
        return imag(f(x+im*epsilon)) / epsilon
    end
    fdtype_error(returntype)
end
# Single point derivative implementations end here.

#=
Multi-point implementations of scalar derivatives for efficiency.
=#
struct DerivativeCache{CacheType1, CacheType2, fdtype, returntype}
    fx      :: CacheType1
    epsilon :: CacheType2
end

function DerivativeCache(
    x          :: AbstractArray{<:Number},
    fx         :: Union{Nothing,AbstractArray{<:Number}} = nothing,
    epsilon    :: Union{Nothing,AbstractArray{<:Real}} = nothing,
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(x)) where {T1,T2}

    if fdtype==Val{:complex} && !(eltype(returntype)<:Real)
        fdtype_error(returntype)
    end

    if fdtype!=Val{:forward} && typeof(fx)!=Nothing
        @warn("Pre-computed function values are only useful for fdtype==Val{:forward}.")
        _fx = nothing
    else
        # more runtime sanity checks?
        _fx = fx
    end

    if typeof(epsilon)!=Nothing && typeof(x)<:StridedArray && typeof(fx)<:Union{Nothing,StridedArray} && 1==2
        @warn("StridedArrays don't benefit from pre-allocating epsilon.")
        _epsilon = nothing
    elseif typeof(epsilon)!=Nothing && fdtype==Val{:complex}
        @warn("Val{:complex} makes the epsilon array redundant.")
        _epsilon = nothing
    else
        if typeof(epsilon)==Nothing || eltype(epsilon)!=real(eltype(x))
          epsilon = fill(zero(real(eltype(x))), size(x))
        end
        _epsilon = epsilon
    end
    DerivativeCache{typeof(_fx),typeof(_epsilon),fdtype,returntype}(_fx,_epsilon)
end

#=
Compute the derivative df of a scalar-valued map f at a collection of points x.
=#
function finite_difference_derivative(
    f,
    x          :: AbstractArray{<:Number},
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(x),      # return type of f
    fx         :: Union{Nothing,AbstractArray{<:Number}} = nothing,
    epsilon    :: Union{Nothing,AbstractArray{<:Real}} = nothing) where {T1,T2}

    df = fill(zero(returntype), size(x))
    finite_difference_derivative!(df, f, x, fdtype, returntype, fx, epsilon)
end

function finite_difference_derivative!(
    df         :: AbstractArray{<:Number},
    f,
    x          :: AbstractArray{<:Number},
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(x),
    fx         :: Union{Nothing,AbstractArray{<:Number}} = nothing,
    epsilon    :: Union{Nothing,AbstractArray{<:Real}}   = nothing) where {T1,T2}

    cache = DerivativeCache(x, fx, epsilon, fdtype, returntype)
    finite_difference_derivative!(df, f, x, cache)
end

function finite_difference_derivative!(df::AbstractArray{<:Number}, f, x::AbstractArray{<:Number},
    cache::DerivativeCache{T1,T2,fdtype,returntype}) where {T1,T2,fdtype,returntype}

    fx, epsilon = cache.fx, cache.epsilon
    if typeof(epsilon) != Nothing
        epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))
        @. epsilon = compute_epsilon(fdtype, x, epsilon_factor)
    end
    if fdtype == Val{:forward}
        if typeof(fx) == Nothing
            @. df = (f(x+epsilon) - f(x)) / epsilon
        else
            @. df = (f(x+epsilon) - fx) / epsilon
        end
    elseif fdtype == Val{:central}
        @. df = (f(x+epsilon) - f(x-epsilon)) / (2 * epsilon)
    elseif fdtype == Val{:complex} && returntype<:Real
        epsilon_complex = eps(eltype(x))
        @. df = imag(f(x+im*epsilon_complex)) / epsilon_complex
    else
        fdtype_error(returntype)
    end
    df
end

#=
Optimized implementations for StridedArrays.
Essentially, the only difference between these and the AbstractArray case
is that here we can compute the epsilon one by one in local variables and avoid caching it.
=#
function finite_difference_derivative!(df::StridedArray, f, x::StridedArray,
    cache::DerivativeCache{T1,T2,fdtype,returntype}) where {T1,T2,fdtype,returntype}

    epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))
    if fdtype == Val{:forward}
        fx = cache.fx
        @inbounds for i ∈ eachindex(x)
            epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
            x_plus = x[i] + epsilon
            if typeof(fx) == Nothing
                df[i] = (f(x_plus) - f(x[i])) / epsilon
            else
                df[i] = (f(x_plus) - fx[i]) / epsilon
            end
        end
    elseif fdtype == Val{:central}
        @inbounds for i ∈ eachindex(x)
            epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
            epsilon_double_inv = one(typeof(epsilon)) / (2*epsilon)
            x_plus, x_minus = x[i]+epsilon, x[i]-epsilon
            df[i] = (f(x_plus) - f(x_minus)) * epsilon_double_inv
        end
    elseif fdtype == Val{:complex}
        epsilon_complex = eps(eltype(x))
        @inbounds for i ∈ eachindex(x)
            df[i] = imag(f(x[i]+im*epsilon_complex)) / epsilon_complex
        end
    else
        fdtype_error(returntype)
    end
    df
end


text:
fix_literal(id) = if id isa Symbol; Meta.quot(id) else id end
(~)(id) = fix_literal(id)

commutativity(op) = :( ($op)(a, b) == ($op)(b, a) ) |> Rule
right_associative(op) = :( ($op)(a, $(op)(b,c)) => ($op)($(op)(a,b), c) ) |> Rule
left_associative(op) = :( ($op)($(op)(a,b), c) => ($op)(a, $(op)(b,c)) ) |> Rule

associativity_left(op) = :( ($op)(a, $(op)(b,c)) => ($op)($(op)(a,b), c) ) |> Rule
associativity_right(op) = :(  ($op)($(op)(a,b), c) => ($op)(a, $(op)(b,c)) ) |> Rule

associativity(op) = [associativity_left(op), associativity_right(op)]

identity_left(op, id) = let id = ~id; :( ($op)($id, a) => a ) |> Rule end
identity_right(op, id) = let id = ~id; :( ($op)(a, $id) => a ) |> Rule end

inverse_left(op, id, invop) = let id = ~id; :( ($op)(($invop)(a), a) => $id ) |> Rule end
inverse_right(op, id, invop) = let id = ~id; :( ($op)(a, ($invop)(a)) => $id ) |> Rule end

# distributivity of two operations
# example: `@distrib (⋅) (⊕)`
function distrib_left(outop, inop)
	@assert Base.isbinaryoperator(outop)
	@assert Base.isbinaryoperator(inop)
	:( ($outop)(a, ($inop)(b,c)) == ($inop)(($outop)(a,b),($outop)(a,c)) ) |> Rule
end

function distrib_right(outop, inop)
	@assert Base.isbinaryoperator(outop)
	@assert Base.isbinaryoperator(inop)
	:( ($outop)(($inop)(a,b), c) == ($inop)(($outop)(a,c),($outop)(b,c)) ) |> Rule
end

function monoid(op, id)
	let id = ~id
		@assert Base.isbinaryoperator(op)
		[associativity_left(op), associativity_right(op),
		 identity_left(op, id), identity_right(op,id)]
	end
end
macro monoid(op, id) monoid(op, id) end

function commutative_monoid(op, id)
	let id = ~id;
		@assert Base.isbinaryoperator(op)
		[commutativity(op), associativity_left(op),
		associativity_right(op), identity_left(op, id)]
	end
end
macro commutative_monoid(op, id) commutative_monoid(op, id) end

# constructs a semantic theory about a an abelian group
# The definition of a group does not require that a ⋅ b = b ⋅ a
# for all elements a and b in G. If this additional condition holds,
# then the operation is said to be commutative, and the group is called an abelian group.
function commutative_group(op, id, invop)
	let id = ~id;
		@assert Base.isbinaryoperator(op)
		# @assert Base.isunaryoperator(invop)
		commutative_monoid(op, id) ∪ [inverse_right(op, id, invop)]
	end
end
abelian_group(op, id, invop) = commutative_group(op, id, invop)
macro commutative_group(op, id, invop) commutative_group(op, id, invop) end
macro abelian_group(op, id, invop) commutative_group(op, id, invop) end

distrib(outop, inop) = [
	distrib_left(outop, inop), distrib_right(outop, inop),
]

ext: jl


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

## core text I/O ##

"""
    print([io::IO], xs...)

Write to `io` (or to the default output stream [`stdout`](@ref)
if `io` is not given) a canonical (un-decorated) text representation
of values `xs` if there is one, otherwise call [`show`](@ref).
The representation used by `print` includes minimal formatting and tries to
avoid Julia-specific details.

Printing `nothing` is deprecated and will throw an error in the future.

# Examples
```jldoctest
julia> print("Hello World!")
Hello World!
julia> io = IOBuffer();

julia> print(io, "Hello", ' ', :World!)

julia> String(take!(io))
"Hello World!"
```
"""
function print(io::IO, x)
    lock(io)
    try
        show(io, x)
    finally
        unlock(io)
    end
    return nothing
end

function print(io::IO, xs...)
    lock(io)
    try
        for x in xs
            print(io, x)
        end
    finally
        unlock(io)
    end
    return nothing
end

"""
    println([io::IO], xs...)

Print (using [`print`](@ref)) `xs` followed by a newline.
If `io` is not supplied, prints to [`stdout`](@ref).

# Examples
```jldoctest
julia> println("Hello, world")
Hello, world

julia> io = IOBuffer();

julia> println(io, "Hello, world")

julia> String(take!(io))
"Hello, world\\n"
```
"""
println(io::IO, xs...) = print(io, xs..., '\n')

## conversion of general objects to strings ##

"""
    sprint(f::Function, args...; context=nothing, sizehint=0)

Call the given function with an I/O stream and the supplied extra arguments.
Everything written to this I/O stream is returned as a string.
`context` can be either an [`IOContext`](@ref) whose properties will be used,
or a `Pair` specifying a property and its value. `sizehint` suggests the capacity
of the buffer (in bytes).

The optional keyword argument `context` can be set to `:key=>value` pair
or an `IO` or [`IOContext`](@ref) object whose attributes are used for the I/O
stream passed to `f`.  The optional `sizehint` is a suggersted (in bytes)
to allocate for the buffer used to write the string.

# Examples
```jldoctest
julia> sprint(show, 66.66666; context=:compact => true)
"66.6667"

julia> sprint(showerror, BoundsError([1], 100))
"BoundsError: attempt to access 1-element Array{Int64,1} at index [100]"
```
"""
function sprint(f::Function, args...; context=nothing, sizehint::Integer=0)
    s = IOBuffer(sizehint=sizehint)
    if context !== nothing
        f(IOContext(s, context), args...)
    else
        f(s, args...)
    end
    String(resize!(s.data, s.size))
end

tostr_sizehint(x) = 0
tostr_sizehint(x::AbstractString) = lastindex(x)
tostr_sizehint(x::Float64) = 20
tostr_sizehint(x::Float32) = 12

function print_to_string(xs...; env=nothing)
    if isempty(xs)
        return ""
    end
    # specialized for performance reasons
    s = IOBuffer(sizehint=tostr_sizehint(xs[1]))
    if env !== nothing
        env_io = IOContext(s, env)
        for x in xs
            print(env_io, x)
        end
    else
        for x in xs
            print(s, x)
        end
    end
    String(resize!(s.data, s.size))
end

string_with_env(env, xs...) = print_to_string(xs...; env=env)

"""
    string(xs...)

Create a string from any values using the [`print`](@ref) function.

# Examples
```jldoctest
julia> string("a", 1, true)
"a1true"
```
"""
string(xs...) = print_to_string(xs...)

# note: print uses an encoding determined by `io` (defaults to UTF-8), whereas
#       write uses an encoding determined by `s` (UTF-8 for `String`)
print(io::IO, s::AbstractString) = for c in s; print(io, c); end
write(io::IO, s::AbstractString) = (len = 0; for c in s; len += write(io, c); end; len)
show(io::IO, s::AbstractString) = print_quoted(io, s)

# optimized methods to avoid iterating over chars
write(io::IO, s::Union{String,SubString{String}}) =
    GC.@preserve s unsafe_write(io, pointer(s), reinterpret(UInt, sizeof(s)))
print(io::IO, s::Union{String,SubString{String}}) = (write(io, s); nothing)

## printing literal quoted string data ##

# this is the inverse of print_unescaped_chars(io, s, "\\\")

function print_quoted_literal(io, s::AbstractString)
    print(io, '"')
    for c = s; c == '"' ? print(io, "\\\"") : print(io, c); end
    print(io, '"')
end

"""
    repr(x; context=nothing)

Create a string from any value using the [`show`](@ref) function.

The optional keyword argument `context` can be set to an `IO` or [`IOContext`](@ref)
object whose attributes are used for the I/O stream passed to `show`.

Note that `repr(x)` is usually similar to how the value of `x` would
be entered in Julia.  See also [`repr(MIME("text/plain"), x)`](@ref) to instead
return a "pretty-printed" version of `x` designed more for human consumption,
equivalent to the REPL display of `x`.

# Examples
```jldoctest
julia> repr(1)
"1"

julia> repr(zeros(3))
"[0.0, 0.0, 0.0]"

julia> repr(big(1/3))
"3.33333333333333314829616256247390992939472198486328125e-01"

julia> repr(big(1/3), context=:compact => true)
"3.33333e-01"

```
"""
repr(x; context=nothing) = sprint(show, x; context=context)

# IOBuffer views of a (byte)string:

"""
    IOBuffer(string::String)

Create a read-only `IOBuffer` on the data underlying the given string.

# Examples
```jldoctest
julia> io = IOBuffer("Haho");

julia> String(take!(io))
"Haho"

julia> String(take!(io))
"Haho"
```
"""
IOBuffer(str::String) = IOBuffer(unsafe_wrap(Vector{UInt8}, str))
IOBuffer(s::SubString{String}) = IOBuffer(view(unsafe_wrap(Vector{UInt8}, s.string), s.offset + 1 : s.offset + sizeof(s)))

# join is implemented using IO

"""
    join([io::IO,] strings, delim, [last])

Join an array of `strings` into a single string, inserting the given delimiter between
adjacent strings. If `last` is given, it will be used instead of `delim` between the last
two strings. If `io` is given, the result is written to `io` rather than returned as
as a `String`.  For example,

# Examples
```jldoctest
julia> join(["apples", "bananas", "pineapples"], ", ", " and ")
"apples, bananas and pineapples"
```

`strings` can be any iterable over elements `x` which are convertible to strings
via `print(io::IOBuffer, x)`. `strings` will be printed to `io`.
"""
function join(io::IO, strings, delim, last)
    first = true
    local prev
    for str in strings
        if @isdefined prev
            first ? (first = false) : print(io, delim)
            print(io, prev)
        end
        prev = str
    end
    if @isdefined prev
        first || print(io, last)
        print(io, prev)
    end
    nothing
end
function join(io::IO, strings, delim="")
    # Specialization of the above code when delim==last,
    # which lets us emit (compile) less code
    first = true
    for str in strings
        first ? (first = false) : print(io, delim)
        print(io, str)
    end
end

join(strings) = sprint(join, strings)
join(strings, delim) = sprint(join, strings, delim)
join(strings, delim, last) = sprint(join, strings, delim, last)

## string escaping & unescaping ##

need_full_hex(c::Union{Nothing, AbstractChar}) = c !== nothing && isxdigit(c)
escape_nul(c::Union{Nothing, AbstractChar}) =
    (c !== nothing && '0' <= c <= '7') ? "\\x00" : "\\0"

"""
    escape_string(str::AbstractString[, esc])::AbstractString
    escape_string(io, str::AbstractString[, esc::])::Nothing

General escaping of traditional C and Unicode escape sequences. The first form returns the
escaped string, the second prints the result to `io`.

Backslashes (`\\`) are escaped with a double-backslash (`"\\\\"`). Non-printable
characters are escaped either with their standard C escape codes, `"\\0"` for NUL (if
unambiguous), unicode code point (`"\\u"` prefix) or hex (`"\\x"` prefix).

The optional `esc` argument specifies any additional characters that should also be
escaped by a prepending backslash (`\"` is also escaped by default in the first form).

# Examples
```jldoctest
julia> escape_string("aaa\\nbbb")
"aaa\\\\nbbb"

julia> escape_string("\\xfe\\xff") # invalid utf-8
"\\\\xfe\\\\xff"

julia> escape_string(string('\\u2135','\\0')) # unambiguous
"ℵ\\\\0"

julia> escape_string(string('\\u2135','\\0','0')) # \\0 would be ambiguous
"ℵ\\\\x000"
```

## See also
[`unescape_string`](@ref) for the reverse operation.
"""
function escape_string(io::IO, s::AbstractString, esc="")
    a = Iterators.Stateful(s)
    for c in a
        if c in esc
            print(io, '\\', c)
        elseif isascii(c)
            c == '\0'          ? print(io, escape_nul(peek(a))) :
            c == '\e'          ? print(io, "\\e") :
            c == '\\'          ? print(io, "\\\\") :
            '\a' <= c <= '\r'  ? print(io, '\\', "abtnvfr"[Int(c)-6]) :
            isprint(c)         ? print(io, c) :
                                 print(io, "\\x", string(UInt32(c), base = 16, pad = 2))
        elseif !isoverlong(c) && !ismalformed(c)
            isprint(c)         ? print(io, c) :
            c <= '\x7f'        ? print(io, "\\x", string(UInt32(c), base = 16, pad = 2)) :
            c <= '\uffff'      ? print(io, "\\u", string(UInt32(c), base = 16, pad = need_full_hex(peek(a)) ? 4 : 2)) :
                                 print(io, "\\U", string(UInt32(c), base = 16, pad = need_full_hex(peek(a)) ? 8 : 4))
        else # malformed or overlong
            u = bswap(reinterpret(UInt32, c))
            while true
                print(io, "\\x", string(u % UInt8, base = 16, pad = 2))
                (u >>= 8) == 0 && break
            end
        end
    end
end

escape_string(s::AbstractString, esc=('\"',)) = sprint(escape_string, s, esc, sizehint=lastindex(s))

function print_quoted(io, s::AbstractString)
    print(io, '"')
    escape_string(io, s, ('\"','$')) #"# work around syntax highlighting problem
    print(io, '"')
end

# general unescaping of traditional C and Unicode escape sequences

# TODO: handle unescaping invalid UTF-8 sequences
"""
    unescape_string(str::AbstractString)::AbstractString
    unescape_string(io, str::AbstractString)::Nothing

General unescaping of traditional C and Unicode escape sequences. The first form returns
the escaped string, the second prints the result to `io`.

The following escape sequences are recognised:
 - Escaped backslash (`\\\\`)
 - Escaped double-quote (`\\\"`)
 - Standard C escape sequences (`\\a`, `\\b`, `\\t`, `\\n`, `\\v`, `\\f`, `\\r`, `\\e`)
 - Unicode code points (`\\u` or `\\U` prefixes with 1-4 trailing hex digits)
 - Hex bytes (`\\x` with 1-2 trailing hex digits)
 - Octal bytes (`\\` with 1-3 trailing octal digits)

# Examples
```jldoctest
julia> unescape_string("aaa\\\\nbbb") # C escape sequence
"aaa\\nbbb"

julia> unescape_string("\\\\u03c0") # unicode
"π"

julia> unescape_string("\\\\101") # octal
"A"
```

## See also
[`escape_string`](@ref).
"""
function unescape_string(io, s::AbstractString)
    a = Iterators.Stateful(s)
    for c in a
        if !isempty(a) && c == '\\'
            c = popfirst!(a)
            if c == 'x' || c == 'u' || c == 'U'
                n = k = 0
                m = c == 'x' ? 2 :
                    c == 'u' ? 4 : 8
                while (k += 1) <= m && !isempty(a)
                    nc = peek(a)
                    n = '0' <= nc <= '9' ? n<<4 + (nc-'0') :
                        'a' <= nc <= 'f' ? n<<4 + (nc-'a'+10) :
                        'A' <= nc <= 'F' ? n<<4 + (nc-'A'+10) : break
                    popfirst!(a)
                end
                if k == 1 || n > 0x10ffff
                    u = m == 4 ? 'u' : 'U'
                    throw(ArgumentError("invalid $(m == 2 ? "hex (\\x)" :
                                        "unicode (\\$u)") escape sequence"))
                end
                if m == 2 # \x escape sequence
                    write(io, UInt8(n))
                else
                    print(io, Char(n))
                end
            elseif '0' <= c <= '7'
                k = 1
                n = c-'0'
                while (k += 1) <= 3 && !isempty(a)
                    c = peek(a)
                    n = ('0' <= c <= '7') ? n<<3 + c-'0' : break
                    popfirst!(a)
                end
                if n > 255
                    throw(ArgumentError("octal escape sequence out of range"))
                end
                write(io, UInt8(n))
            else
                print(io, c == 'a' ? '\a' :
                          c == 'b' ? '\b' :
                          c == 't' ? '\t' :
                          c == 'n' ? '\n' :
                          c == 'v' ? '\v' :
                          c == 'f' ? '\f' :
                          c == 'r' ? '\r' :
                          c == 'e' ? '\e' :
                          (c == '\\' || c == '"') ? c :
                          throw(ArgumentError("invalid escape sequence \\$c")))
            end
        else
            print(io, c)
        end
    end
end
unescape_string(s::AbstractString) = sprint(unescape_string, s, sizehint=lastindex(s))

macro b_str(s)
    v = codeunits(unescape_string(s))
    QuoteNode(v)
end

"""
    @raw_str -> String

Create a raw string without interpolation and unescaping.
The exception is that quotation marks still must be escaped. Backslashes
escape both quotation marks and other backslashes, but only when a sequence
of backslashes precedes a quote character. Thus, 2n backslashes followed by
a quote encodes n backslashes and the end of the literal while 2n+1 backslashes
followed by a quote encodes n backslashes followed by a quote character.

# Examples
```jldoctest
julia> println(raw"\\ \$x")
\\ \$x

julia> println(raw"\\"")
"

julia> println(raw"\\\\\\"")
\\"

julia> println(raw"\\\\x \\\\\\"")
\\\\x \\"
```
"""
macro raw_str(s); s; end

## multiline strings ##

"""
    indentation(str::AbstractString; tabwidth=8) -> (Int, Bool)

Calculate the width of leading white space. Return the width and a flag to indicate
if the string is empty.

# Examples
```jldoctest
julia> Base.indentation("")
(0, true)

julia> Base.indentation("  a")
(2, false)

julia> Base.indentation("\\ta"; tabwidth=3)
(3, false)
```
"""
function indentation(str::AbstractString; tabwidth=8)
    count = 0
    for ch in str
        if ch == ' '
            count += 1
        elseif ch == '\t'
            count = div(count + tabwidth, tabwidth) * tabwidth
        else
            return count, false
        end
    end
    count, true
end

"""
    unindent(str::AbstractString, indent::Int; tabwidth=8)

Remove leading indentation from string.

# Examples
```jldoctest
julia> Base.unindent("   a\\n   b", 2)
" a\\n b"

julia> Base.unindent("\\ta\\n\\tb", 2, tabwidth=8)
"      a\\n      b"
```
"""
function unindent(str::AbstractString, indent::Int; tabwidth=8)
    indent == 0 && return str
    # Note: this loses the type of the original string
    buf = IOBuffer(sizehint=sizeof(str))
    cutting = true
    col = 0     # current column (0 based)
    for ch in str
        if cutting
            if ch == ' '
                col += 1
            elseif ch == '\t'
                col = div(col + tabwidth, tabwidth) * tabwidth
            elseif ch == '\n'
                # Now we need to output enough indentation
                for i = 1:col-indent
                    print(buf, ' ')
                end
                col = 0
                print(buf, '\n')
            else
                cutting = false
                # Now we need to output enough indentation to get to
                # correct place
                for i = 1:col-indent
                    print(buf, ' ')
                end
                col += 1
                print(buf, ch)
            end
        elseif ch == '\t'       # Handle internal tabs
            upd = div(col + tabwidth, tabwidth) * tabwidth
            # output the number of spaces that would have been seen
            # with original indentation
            for i = 1:(upd-col)
                print(buf, ' ')
            end
            col = upd
        elseif ch == '\n'
            cutting = true
            col = 0
            print(buf, '\n')
        else
            col += 1
            print(buf, ch)
        end
    end
    # If we were still "cutting" when we hit the end of the string,
    # we need to output the right number of spaces for the indentation
    if cutting
        for i = 1:col-indent
            print(buf, ' ')
        end
    end
    String(take!(buf))
end

function String(chars::AbstractVector{<:AbstractChar})
    sprint(sizehint=length(chars)) do io
        for c in chars
            print(io, c)
        end
    end
end

max_stars_repo_name: pochoi/julia

max_issues_repo_name: pochoi/julia

max_forks_repo_name: pochoi/julia


text:
################################################################################
#
#  The algorithm of Husert
#
################################################################################

function _issimilar_husert(A::fmpz_mat, B::fmpz_mat)
  QA = change_base_ring(FlintQQ, A)
  QB = change_base_ring(FlintQQ, B)
  fl, QC = _issimilar_husert_generic(QA, QB)
  return fl, change_base_ring(FlintZZ, QC)
end

# If successful, returns X such that X * A * X^-1 == B
function _issimilar_husert_generic(A, B)
  @assert issquarefree(minpoly(A))
  m = nrows(A)
  f = charpoly(A)
  fac = factor(f)
  mus = typeof(f)[]
  ns = Int[]
  ds = Int[]
  for (p, n) in fac
    push!(mus, p)
    push!(ns, n)
    push!(ds, degree(p))
  end

  s = length(mus)
  Ks = AnticNumberField[]
  vecsA = Matrix{nf_elem}(undef, m, sum(ns))
  vecsB = Matrix{nf_elem}(undef, m, sum(ns))
  k = 1
  for i in 1:s
    K, _ = NumberField(mus[i], "a", cached = false)
    EA = eigenspace(change_base_ring(K, A), gen(K), side = :left)
    EB = eigenspace(change_base_ring(K, B), gen(K), side = :left)
    push!(Ks, K)
    if nrows(EA) != nrows(EB)
      return false, zero_matrix(FlintQQ, 0, 0)
    end
    @assert nrows(EA) == nrows(EB)
    for j in 1:nrows(EA)
      for l in 1:m
        vecsA[l, k] = EA[j, l]
        vecsB[l, k] = EB[j, l]
      end
      k += 1
    end
  end

  absolute_basis_vec_A = Vector{Vector{fmpq}}()
  absolute_basis_vec_B = Vector{Vector{fmpq}}()
  absolute_basis_A = zero_matrix(FlintQQ, m, m)
  absolute_basis_B = zero_matrix(FlintQQ, m, m)
  for i in 1:m
    vA = _to_absolute_basis(view(vecsA, i, :), m, ns, Ks)
    vB = _to_absolute_basis(view(vecsB, i, :), m, ns, Ks)
    push!(absolute_basis_vec_A, vA)
    push!(absolute_basis_vec_B, vB)
    for j in 1:m
      absolute_basis_A[i, j] = vA[j]
      absolute_basis_B[i, j] = vB[j]
    end
  end

  # Now construct the colon ideal
  # First the Q-basis of \prod Mat(n_i, K_i)
  actions = Vector{fmpq_mat}()
  another_basis_of_actions = Vector{Vector{dense_matrix_type(nf_elem)}}()

  for i in 1:s
    ni = ns[i]
    for i1 in 1:ni
      for i2 in 1:ni
        for j in 1:degree(Ks[i])
          V = Vector{dense_matrix_type(nf_elem)}(undef, s)
          M = representation_matrix(gen(Ks[i])^(j - 1))
          z = zero_matrix(Ks[i], ni, ni)
          z[i1, i2] = gen(Ks[i])^(j-1)
          V[i] = z
          for j in 1:s
            if j != i
              V[j] = zero_matrix(Ks[j], ns[j], ns[j])
            end
          end
          push!(another_basis_of_actions, V)
          # I need to construct the block matrix, where M is at block (i1, i2)
          # and all the other ones are zero
          # The overall matrix must have size ni * degree(Ks[i])
          Z1 = zero_matrix(FlintQQ, (i1 - 1) * ds[i], ni * ds[i])
          Z2 = zero_matrix(FlintQQ, ds[i], (i2 - 1) * ds[i])
          Z3 = hcat(Z2, M)
          Z3 = hcat(Z3, zero_matrix(FlintQQ, ds[i], (ni - i2) * ds[i]))
          Z4 = zero_matrix(FlintQQ, (ni - i1) * ds[i], ni * ds[i])
          MM = vcat(Z1, Z3, Z4)
          blocklengthbefore = sum(Int[ns[o] * ds[o] for o in 1:(i - 1)])
          blocklengthafter = sum(Int[ns[o] * ds[o] for o in (i + 1):s])
          D = diagonal_matrix(
                zero_matrix(FlintQQ, blocklengthbefore, blocklengthbefore),
                MM,
                zero_matrix(FlintQQ, blocklengthafter, blocklengthafter))
          push!(actions, D)
        end
      end
    end
  end

  absolute_basis_B_inv = inv(absolute_basis_B)
  absolute_basis_A_inv = inv(absolute_basis_A)

  # So D is a Q-basis of Mat_n(K)

  # Now determine the colon ideal

  la = length(actions)

  M = zero_matrix(FlintQQ, la, 0)

  for a in absolute_basis_vec_A
    _M = zero_matrix(FlintQQ, la, m)
    for i in 1:length(actions)
      D = actions[i]
      v = matrix(FlintQQ, 1, m, a) * D * absolute_basis_B_inv
      for j in 1:m
        _M[i, j] = v[1, j]
      end
    end
    M = hcat(M, _M)
  end

  M = sub(hnf(FakeFmpqMat(transpose(M)), :upperright), 1:la, 1:la)
  N = inv(M)

  SS = N

  M = zero_matrix(FlintQQ, length(actions), 0)

  for a in absolute_basis_vec_B
    _M = zero_matrix(FlintQQ, length(actions), m)
    for i in 1:length(actions)
      D = actions[i]
      v = matrix(FlintQQ, 1, m, a) * D * absolute_basis_B_inv
      for j in 1:m
        _M[i, j] = v[j]
      end
    end
    M = hcat(M, _M)
  end
  M = sub(hnf(FakeFmpqMat(transpose(M)), :upperright), 1:la, 1:la)
  N = inv(M)

  bcolonb = N

  basis_of_colon_ideal = fmpq_mat[sum(SS[i, j] * actions[i]
                                      for i in 1:length(actions))
                                      for j in 1:nrows(SS)]

  basis_of_colon_idealAA = fmpq_mat[sum(bcolonb[i, j] * actions[i]
                                      for i in 1:length(actions))
                                      for j in 1:nrows(bcolonb)]

  other_basis_of_colon_ideal = Vector{dense_matrix_type(nf_elem)}[]
  other_basis_of_colon_idealAA = Vector{dense_matrix_type(nf_elem)}[]

  for j in 1:nrows(SS)
    ob = foldl((x, y) -> x .+ y,
               [SS[i, j] .* another_basis_of_actions[i] for i in 1:nrows(SS)])
    push!(other_basis_of_colon_ideal, ob)
  end

  for j in 1:nrows(bcolonb)
    ob = foldl((x, y) -> x .+ y,
               [bcolonb[i, j] .* another_basis_of_actions[i] for i in 1:nrows(bcolonb)])
    push!(other_basis_of_colon_idealAA, ob)
  end

  for a in absolute_basis_vec_A
    for i in 1:length(basis_of_colon_ideal)
      @assert denominator(matrix(FlintQQ, 1, m, a) *
                            basis_of_colon_ideal[i] * absolute_basis_B_inv) == 1
    end
  end

  A = _matrix_algebra(Ks, ns)

  idealAgens = elem_type(A)[]

  ordergens = elem_type(A)[]

  dec = decompose(A)

  for b in other_basis_of_colon_ideal
    z = zero(A)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      B, mB = dec[i]
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(B))
    end
    push!(idealAgens, z)
  end

  for b in other_basis_of_colon_idealAA
    z = zero(A)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      B, mB = dec[i]
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(B))
    end
    push!(ordergens, z)
  end

  ide = ideal_from_lattice_gens(A, idealAgens)
  ideO = Order(A, ordergens)

  @assert ideO == right_order(ide)

  #@assert right_order(ide) == ideO

  fl, y = _isprincipal(ide, ideO, :right)

  if fl
    dec = decompose(A)
    D = Vector{Generic.MatSpaceElem{nf_elem}}(undef, length(dec))
    for i in 1:length(dec)
      B, mB = dec[i]
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = BtoC(mB\y)::elem_type(C)
      D[i] = matrix(z)
    end
    DD = diagonal_matrix(fmpq_mat[_explode(D[i]) for i in 1:length(D)])
    return fl, transpose(absolute_basis_A * DD * absolute_basis_B_inv)
  else
    return false, zero_matrix(FlintQQ, 0, 0)
  end
end

function _explode(x::Generic.MatSpaceElem{nf_elem})
  K = base_ring(x)
  d = degree(K)
  n = nrows(x)
  m = ncols(x)
  z = zero_matrix(FlintQQ, 0, n * d)
  for i in 1:n
    zz = zero_matrix(FlintQQ, d, 0)
    for j in 1:m
      M = representation_matrix(x[i, j])
      zz = hcat(zz, M)
    end
    z = vcat(z, zz)
  end
  return z
end

Base.:(*)(x::fmpq, y::AbstractAlgebra.Generic.MatSpaceElem{nf_elem}) = base_ring(y)(x) * y

function _to_absolute_basis(v, m, ns, Ks)
  w = Vector{fmpq}(undef, m)
  k = 1
  for i in 1:sum(ns)
    vi = v[i]
    K = parent(vi)
    for j in 1:degree(K)
      w[k] = coeff(vi, j - 1)
      k += 1
    end
  end
  return w
end

function _to_relative_basis(v, m, ns, Ks)
  w = Vector{nf_elem}(undef, sum(ns))
  k = 1
  l = 1
  for i in 1:length(ns)
    a = gen(Ks[i])
    for j in 1:ns[i]
      w[k] = zero(Ks[i])
      for n in 1:degree(Ks[i])
        w[k] = w[k] + a^(n - 1) * v[l]
        l += 1
      end
      k += 1
    end
  end
  return w
end

function _matrix_algebra(Ks, ns)
  s = length(Ks)
  algs = AlgAss{fmpq}[]
  for i in 1:s
    A = matrix_algebra(Ks[i], ns[i])
    B, BtoA = AlgAss(A)
    C, CtoB = restrict_scalars(B, FlintQQ)
    C.isomorphic_full_matrix_algebra = (A, CtoB * BtoA)
    push!(algs, C)
  end
  return _my_direct_product(algs)
end


text:
using CausalDiscovery
using Test

# @testset "CausalDiscovery.jl" begin
#   # Write your own tests here.
#   include("CISC.jl")
#   include("model.jl")
# end

@testset "Autumn" begin
  # include("Autumn/sexpr.jl")
  # include("Autumn/transform.jl")
  include("Autumn/particles.jl")
  include("Autumn/compile.jl")
  #commented out because it will not work until support for cells is added
  # include("Autumn/best_effort_ranking.jl")

end


text:
#    Main Module of the  PosDefManifold Package for julia language
#    v 0.0.3 - last update 16th of May 2019
#
#    MIT License
#    Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#    https://sites.google.com/site/marcocongedo/home

# __precompile__()

module PosDefManifold

using LinearAlgebra, Statistics, Base.Threads

# Special instructions and variables
BLAS.set_num_threads(Sys.CPU_THREADS-Threads.nthreads())

# constants
const sqrt2=√2
const invsqrt2=1/sqrt2
const maxpos=1e15

# aliases
𝚺 = sum                 # tab-completition: \bfSigma
𝛍 = mean                # tab-completition: \bfmu
𝕄 = Matrix	            # tab-completion: \bbM
𝔻 = Diagonal	        # tab-completition: \bbD
ℍ = Hermitian          # tab completion \bbH
𝕃 = LowerTriangular    # tab completition \bbL

# types
MatrixVector = 𝕄Vector = Vector{𝕄}            # vector of Matrices
MatrixVector₂= 𝕄Vector₂= Vector{𝕄Vector}      # vector of vectors of Matrices

DiagonalVector = 𝔻Vector = Vector{𝔻}          # vector of Diagonal Matrices
DiagonalVector₂= 𝔻Vector₂= Vector{𝔻Vector}    # vector of vectors of Diagonal matrices

LowerTriangularVector = 𝕃Vector = Vector{𝕃}       # vector of LowerTriangular Matrices
LowerTriangularVector₂= 𝕃Vector₂= Vector{𝕃Vector} # vector of vectors of Diagonal matrices

HermitianVector = ℍVector = Vector{ℍ}           # vector of Hermitian matrices
HermitianVector₂= ℍVector₂= Vector{ℍVector}     # vector of vectors of Hermitian matrices

RealOrComplex=Union{Real, Complex}

AnyMatrix=Union{𝔻{T}, 𝕃{T}, ℍ{T}, 𝕄{T}} where T<:RealOrComplex
AnyMatrixVector=Union{𝕄Vector, 𝔻Vector, 𝕃Vector, ℍVector}
AnyMatrixVector₂=Union{𝕄Vector₂, 𝔻Vector₂, 𝕃Vector₂, ℍVector₂}

@enum Metric begin
    Euclidean    =  1  # distance: δ_e; mean: Arithmetic
    invEuclidean =  2  # distance: δ_i; mean: Harmonic
    ChoEuclidean =  3  # distance: δ_c; mean: Cholesky Euclidean
    logEuclidean =  4  # distance: δ_l; mean: Log Euclidean
    logCholesky  =  5  # distance: δ_c; mean: Log-Cholesky
    Fisher       =  6  # distance: δ_f; mean: Fisher (Cartan, Karcher, Pusz-Woronowicz,...)
    logdet0      =  7  # distance: δ_s; mean: LogDet (S, α, Bhattacharyya, Jensen,...)
    Jeffrey      =  8  # distance: δ_j; mean: Jeffrey (symmetrized Kullback-Leibler)
    VonNeumann   =  9  # distance: δ_v; mean: Not Availale
    Wasserstein  = 10 # distance: δ_w; mean: Wasserstein (Bures, Hellinger, ...)
    #...
end

import
    Statistics.mean,
    LinearAlgebra.tr

export
    # From this module

    #constants
    sqrt2,
    invsqrt2,
    minpos,
    maxpos,

    #aliases
    𝚺,
    𝛍,
    𝕄,
    𝔻,
    ℍ,
    𝕃,

    #types
    MatrixVector, 𝕄Vector,
    MatrixVector₂, 𝕄Vector₂,
    DiagonalVector, 𝔻Vector,
    DiagonalVector₂, 𝔻Vector₂,
    LowerTriangularVector, 𝕃Vector,
    LowerTriangularVector₂, 𝕃Vector₂,
    HermitianVector, ℍVector,
    HermitianVector₂, ℍVector₂,
    RealOrComplex,
    AnyMatrix,
    AnyMatrixVector,
    AnyMatrixVector₂,
    Metric,
        Euclidean,
        invEuclidean,
        ChoEuclidean,
        logEuclidean,
        logCholesky,
        Fisher,
        logdet0,
        Jeffrey,
        VonNeumann,
        Wasserstein,

    # from LinearAlgebra.jl
    typeofMatrix, typeofMat,
    dim,
    det1,
    tr1,
    normalizeCol!,
    colProd,
    sumOfSqr, ss,
    sumOfSqrDiag, ssd,
    colNorm,
    sumOfSqrTril, sst,
    tr,
    quadraticForm, qf,
    fidelity,
    fDiag, 𝑓𝔻,
    DiagOfProd, dop,
    mgs,
    evd,
    spectralFunctions,
    ispos,
    pow,
    invsqrt,
    sqr,
    powerIterations, powIter,
    choL,

    # from SignalProcessing.jl
    randChi², randχ²,
    randEigvals, randλ,
    randEigvalsMat, randΛ,
    randUnitaryMat, randOrthMat, randU,
    randPosDefMat, randP,
    regularize!,
    gram,
    trade,

    # from RiemannianGeometry.jl
    geodesic,
    distanceSqr, distance²,
    distance,
    distanceSqrMat, distance²Mat,
    distanceMat,
    laplacian,
    laplacianEigenMaps, laplacianEM,
    spectralEmbedding,
    mean,
    means,
    generalizedMean,
    geometricMean, gMean,
    logdet0Mean, ld0Mean,
    wasMean,
    powerMean,
    logMap,
    expMap,
    vecP,
    matP,
    procrustes,

    # from classification.jl
    softmax,

    # from Test.jl
    testall

include("linearAlgebra.jl")
include("signalProcessing.jl")
include("riemannianGeometry.jl")
include("classification.jl")
include("test.jl")

println("\n⭐ "," Welcome to the PosDefManifold package", " ⭐\n")
@info(" ")
println(" Your Machine ",gethostname()," (",Sys.MACHINE, ")")
println(" runs on kernel ",Sys.KERNEL," with word size ",Sys.WORD_SIZE,".")
println(" CPU  Threads: ",Sys.CPU_THREADS)
# Sys.BINDIR # julia bin directory
println(" Base.Threads: ", "$(Threads.nthreads())")
println(" BLAS Threads: ", "$(Sys.CPU_THREADS-Threads.nthreads())", "\n")

end # module


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")

using gridding
import PyPlot.plt
using LaTeXStrings

etas = linspace(-1.0, 1.0, 100)

alphas = Float64[0., 0.5, 1.0, 1.5, 2.0]

fig = plt.figure(figsize=(6,6))
ax = fig[:add_subplot](111)

for alpha in alphas
    label = @sprintf("%.1f", alpha)
    ax[:plot](etas, spheroid(etas, alpha), label=label)
end

ax[:set_xlabel](L"$\eta$")
ax[:legend](loc="upper right")
plt.savefig("../plots/spheroids.png")

fig = plt.figure(figsize=(6,6))
ax = fig[:add_subplot](111)

for alpha in alphas
    label = @sprintf("%.1f", alpha)
    ax[:plot](etas, corrfun(etas, alpha), label=label)
end

ax[:set_xlabel](L"$\eta$")
ax[:legend](loc="upper right")
plt.savefig("../plots/corrfuns.png")

fig = plt.figure(figsize=(6,6))
ax = fig[:add_subplot](111)

for alpha in alphas
    label = @sprintf("%.1f", alpha)
    ax[:plot](etas, gcffun(etas, alpha), label=label)
end

ax[:set_xlabel](L"$\eta$")
ax[:legend](loc="upper right")
plt.savefig("../plots/gcffuns.png")


text:
include("ghtype_tests.jl")
include("event_tests.jl")
include("read_only_api_tests.jl")

include("branches.jl")
include("owners.jl")
## include("comments.jl") ## these are blocking tests
## include("events.jl") ## these are blocking tests
include("statuses.jl")
include("contents.jl")
include("repos.jl")


text:
module StoreTest

using FactCheck

include(Pkg.dir("GenomeAnnotations", "src","Config.jl"))
include(Pkg.dir("GenomeAnnotations", "src", "Store.jl" ) )
local_store_path =  Pkg.dir("GenomeAnnotations", "testdata","store", "annotation")

SET_LOCAL_STORE_PATH(local_store_path)
setup_genome_annotations( )

# test ucscgb
include( Pkg.dir("GenomeAnnotations", "src","UCSCGB.jl") )
install_ucsc_annotation ( "mm10" )

# test annotation list now we have some annotations
a = current_annotation_files( )
println(a)
end


text:
"""
    set_coordinates(x, y)

Input x and y airfoil coordinates into XFOIL. Coordinates must start at the
trailing edge and loop counterclockwise.
"""
set_coordinates

"""
    set_coordinates_cs(x, y)

`set_coordinates` for complex step enabled version of XFOIL.
"""
set_coordinates_cs

# definition of set_coordinates
for (T, name, globals, library) in
    ((:Float64, :set_coordinates, :xfoilglobals, :libxfoil_light),
    (:ComplexF64, :set_coordinates_cs, :xfoilglobals_cs, :libxfoil_light_cs))

    @eval begin

        function $(name)(x, y)

            nx = length(x)

            @assert nx == length(y) "x and y-coordinate arrays must match in length"
            @assert nx <= IBX "x and y-coordinate array lengths must be less than $IBX"

            $(globals).nb[1] = nx
            $(globals).xb[1:nx] = x
            $(globals).yb[1:nx] = y

            ccall( (:xfoil_, $(library)), Nothing, ())

            return nothing
        end

    end

end

"""
    pane(; kwargs...)

Run XFOIL's PANE command (Repanel the input airfoil)

# Arguments
- `npan::Integer=140`: Number of panel nodes
- `cvpar::Float64=1.0`: Panel bunching parameter
- `cterat::Float64=0.15`: TE/LE panel density ratio
- `ctrrat::Float64=0.2`: Refined-area/LE panel density ratio
- `xsref1::Float64=1.0`: Top side refined area x/c limits
- `xsref2::Float64=1.0`:
- `xpref1::Float64=1.0`: Bottom side refined area x/c limits
- `xpref2::Float64=1.0`:
"""
pane

"""
    pane_cs(; kwargs...)

`pane` for complex step enabled version of XFOIL.
"""
pane_cs

# definition of pane
for (T, name, globals, library) in
    ((:Float64, :pane, :xfoilglobals, :libxfoil_light),
    (:ComplexF64, :pane_cs, :xfoilglobals_cs, :libxfoil_light_cs))

    @eval begin

        function $(name)(; npan=140, cvpar=1.0, cterat=0.15, ctrrat=0.2, xsref1=1.0,
            xsref2=1.0, xpref1=1.0, xpref2=1.0)

            $(globals).npan[1] = npan
            $(globals).cvpar[1] = cvpar
            $(globals).cterat[1] = cterat
            $(globals).ctrrat[1] = ctrrat
            $(globals).xsref1[1] = xsref1
            $(globals).xsref2[1] = xsref2
            $(globals).xpref1[1] = xpref1
            $(globals).xpref2[1] = xpref2

            ccall( (:pangen_, $(library)), Nothing, ())

            return nothing
        end

    end

end

"""
    solve_alpha(alpha, re; mach=0.0, iter=50, ncrit=9)

Compute the flow solution at the specified angle of attack. Returns cl, cd, cdp,
cm, and a convergence flag indicating whether the solution converged.

# Arguments:
 - `alpha`: Angle of attack (degrees)
 - `re`: Reynolds number
 - `mach`: Mach number
 - `iter`: Number of iterations
 - `ncrit`: turbulence level
"""
solve_alpha

"""
    solve_alpha_cs(alpha, re; mach=0.0, iter=50, ncrit=9)

`solve_alpha` for the complex step enabled version of XFOIL.
"""
solve_alpha_cs

# definition of solve_alpha
for (T, name, globals, library) in
    ((:Float64, :solve_alpha, :xfoilglobals, :libxfoil_light),
    (:ComplexF64, :solve_alpha_cs, :xfoilglobals_cs, :libxfoil_light_cs))

    @eval begin

        function $(name)(alpha, re; mach=0.0, iter=50, ncrit=9)

            $(globals).adeg[1] = alpha
            $(globals).reinf1[1] = re
            $(globals).minf1[1] = mach
            $(globals).itmax[1] = iter
            $(globals).acrit[1] = ncrit

            ccall((:oper_, $(library)), Nothing, ())

            cl = $(globals).cl[1]
            cd = $(globals).cd[1]
            cdp = $(globals).cdp[1]
            cm = $(globals).cm[1]
            lvconv = convert(Bool, $(globals).lvconv[1])

            return cl[1], cd[1], cdp[1], cm[1], lvconv
        end

    end

end

"""
    bldump()

Return boundary layer parameters `s`, `x`, `y`, `ue`, `dstar`, `theta`, and `cf`
"""
bldump

"""
    bldump_cs()

`bldump` for the complex step enabled version of XFOIL.
"""
bldump_cs

# definition of bldump
for (T, name, globals, library) in
    ((:Float64, :bldump, :xfoilglobals, :libxfoil_light),
    (:ComplexF64, :bldump_cs, :xfoilglobals_cs, :libxfoil_light_cs))

    @eval begin

        function $(name)()

            nelem = zeros(Int32, 1)
            s = zeros($T, IZX)
            x = zeros($T, IZX)
            y = zeros($T, IZX)
            ue = zeros($T, IZX)
            dstar = zeros($T, IZX)
            theta = zeros($T, IZX)
            cf = zeros($T, IZX)

            ccall( (:bldump_, $(library)), Nothing,
                (Ref{Int32}, Ref{$T}, Ref{$T}, Ref{$T}, Ref{$T},
                Ref{$T}, Ref{$T}, Ref{$T}),
                nelem, s, x, y, ue, dstar, theta, cf)

            return s[1:nelem[1]], x[1:nelem[1]], y[1:nelem[1]], ue[1:nelem[1]],
                dstar[1:nelem[1]], theta[1:nelem[1]], cf[1:nelem[1]]
        end

    end

end


text:
@testset "IDW" begin
  # basic problem
  geodata = georef((variable=[1.,0.,1.],), [25. 50. 75.;  25. 75. 50.])
  domain  = CartesianGrid(100,100)
  problem = EstimationProblem(geodata, domain, :variable)

  solver = IDW(:variable => (neighbors=3,))

  solution = solve(problem, solver)

  if visualtests
    @test_reference "data/idw.png" contourf(solution,size=(800,400))
  end

  # haversine distance
  geodata = georef((variable=[4.0,-1.0,3.0],), [50.0 100.0 200.0; -30.0 30.0 10.0])
  domain  = CartesianGrid((1.0, -89.0), (359.0, 89.0), dims=(200, 100))
  problem = EstimationProblem(geodata, domain, :variable)

  solver = IDW(:variable => (neighbors=3, distance=Haversine(1.0)))

  solution = solve(problem, solver)

  if visualtests
    @test_reference "data/idw-haversine.png" contourf(solution,size=(800,200))
  end
end


text:
import Base.Docs: Binding

# Keywords
# ––––––––

const builtins = ["abstract", "baremodule", "begin", "bitstype", "break",
                  "catch", "ccall", "const", "continue", "do", "else",
                  "elseif", "end", "export", "finally", "for", "function",
                  "global", "if", "immutable", "import", "importall", "let",
                  "local", "macro", "module", "quote", "return", "try", "type",
                  "typealias", "using", "while"]

const builtin_completions =
  [d(:text=>k, :type=>:keyword)
   for k in builtins]

# Module completions
# ––––––––––––––––––

const identifier_pattern = r"^@?[_\p{L}][\p{Xwd}!]*+$"

moduleusings(mod) = ccall(:jl_module_usings, Any, (Any,), mod)

filtervalid(names) = @>> names map(string) filter(x->ismatch(identifier_pattern, x))

accessible(mod::Module) =
  [names(mod, true, true);
   map(names, moduleusings(mod))...] |> unique |> filtervalid

Base.getindex(b::Binding) = getfield(b.mod, b.var)

completiontype(x) =
  isa(x, Module) ? "module" :
  isa(x, DataType) ? "type" :
  isa(x, Function) ? "λ" :
  "constant"

const meta_cache = Dict{Any,Any}()

function withmeta(completion::AString, mod::Module)
  isdefined(mod, Symbol(completion)) || return completion
  b = Binding(mod, Symbol(completion))
  mod = b.mod
  haskey(meta_cache, (mod, completion)) && return meta_cache[(mod, completion)]
  x = b[]
  c = d(:text => completion,
        :type => Base.startswith(completion, "@") ? "macro" : completiontype(x),
        :rightLabel => string(mod))
  c[:displayText] = signature(b)
  c[:description] = description(b)
  c
end

for name in map(string, names(Base))
  meta_cache[(Base, name)] = withmeta(name, Base)
end

withmeta(completions::Vector, mod::Module) =
  [withmeta(completion, mod) for completion in completions]

function namecompletions_(mod::Module, qualified = false)
  if !qualified
    [withmeta(filter!(x->!Base.isdeprecated(mod, Symbol(x)), accessible(mod)), mod); builtin_completions]
  else
    withmeta(filter!(x->!Base.isdeprecated(mod, Symbol(x)), filtervalid(names(mod, true))), mod)
  end
end

const namecompletions = memoize_debounce(namecompletions_)

const prefix_pattern = r"(@?[_\p{L}][\p{Xwd}!]*+\.?)+$|@$"

function prefix(line)
  match = Base.match(prefix_pattern, line)
  match == nothing && return String[]
  split(match.match, ".")
end

# Path Completions
# ––––––––––––––––

function funcprefix(line)
  m = match(r"\b([^\s]*)\(([^\)]*)$", line) # matches `foo(bar`
  m == nothing && return
  return m.captures[1], m.captures[2]
end

function pathprefix(line)
  m = match(r"\"?([\w\\/:\-.]*)$", line)
  m == nothing && return
  return m.match, m.captures[1]
end

function stringmeta(cs, prefix)
  map(c -> d(:text => c, :_prefix => prefix, :type => :file), cs)
end

function pathmeta(cs, path, prefix)
  stringmeta(map(c -> replace(joinpath(path, c), r"^\./", ""), cs), prefix)
end

function children(dir, ext = ""; depth = 0, out = String[], prefix = "")
  isdir(dir) || return []
  for f in readdir(dir)
    path = joinpath(dir, f)
    if isfile(path) && ext ≠ nothing && endswith(path, ext)
      push!(out, joinpath(prefix, f))
    elseif isdir(path)
      push!(out, joinpath(prefix, f, ""))
      if depth > 0
        children(path, ext, depth = depth-1, out = out, prefix = joinpath(prefix, f))
      else
      end
    end
  end
  return out
end

function includepath()
  path = Base.source_path()
  path == nothing && return pwd()
  return dirname(path)
end

function pathcompletions(line)
  m = funcprefix(line)
  m == nothing && return
  func, pre = m
  m = pathprefix(pre)
  m == nothing && return
  pre, path = m
  dir = dirname(path)
  if func == "include"
    pathmeta(children(joinpath(includepath(), dir), ".jl", depth = 2), dir, path)
  elseif func == "readcsv"
    pathmeta(children(joinpath(pwd(), dir), ".csv", depth = 2), dir, path)
  elseif func == "cd"
    pathmeta(children(joinpath(pwd(), dir), nothing), dir, path)
  elseif func == "open"
    pathmeta(children(joinpath(pwd(), dir)), dir, path)
  end
end

# Package Completions
# –––––––––––––––––––

packages(dir = Pkg.dir()) =
  @>> dir readdir filter(x->!ismatch(r"^\.|^METADATA$|^REQUIRE$", x))

all_packages() = packages(Pkg.dir("METADATA"))

required_packages() =
  @>> Pkg.dir("REQUIRE") readstring lines

unused_packages() = setdiff(all_packages(), required_packages())

pkgmeta(xs) = map(x -> d(:text=>x, :type=>"package"), xs)

function pkgcompletions(line)
  if ismatch(r"^using", line)
    return pkgmeta(packages())
  end
  m = funcprefix(line)
  m == nothing && return
  func, _ = m
  if func in ["Pkg.add", "Pkg.clone"]
    pkgmeta(all_packages())
  elseif func in ["Pkg.pin", "Pkg.checkout"]
    pkgmeta(packages())
  elseif func in ["Pkg.rm"]
    pkgmeta(required_packages())
  end
end

# Completions
# –––––––––––

const providers = [pkgcompletions, pathcompletions]

function completions(line::AString, mod::Module = Main; default = true)
  for provider in providers
    cs = provider(line)
    cs ≠ nothing && return cs
  end
  pre = prefix(line)
  mod = getthing(mod, pre[1:end-1])
  if isa(mod, Module)
    qualified = length(pre) > 1
    !(qualified || default) && return
    return namecompletions(mod, qualified)
  end
  return []
end

completions(mod::Module) = completions("", mod, default = true)


text:
test("IntPoint ↔ Tuple") do
	x=1; y=0
	xy=(x,y)
	ip_direct=IntPoint(x,y)
	ip_fromtp=IntPoint(xy)
	@test ip_direct==ip_fromtp

	@test easyClipper.IntPoint2Tuple(ip_direct)==xy
	@test all( easyClipper.IntPoint2Vector(ip_direct) .== [x,y] )
end
test("Tuple extrema") do
	pointsvec = [(0,0),(1,0),(0,2),(-1,1)]
	minpoint,maxpoint=easyClipper.tupleextrema(pointsvec)
	@test minpoint == (-1,0)
	@test maxpoint == (1,2)
end
test("Rescaling julia DataTypes to full-range Cint") do
	pointsvec = NTuple{2,Float64}[(0.,0.),(1.,2.)]
	minx,xwidth=easyClipper.min_width(pointsvec)
	@test minx == (0.,0.)
	@test xwidth == (1.,2.)
	cmax,o,t=easyClipper.maxCint_one_two(Float64)
	Cpointsvec = easyClipper._scale_j2c(pointsvec,minx,xwidth)
	@test Cpointsvec[1] == (-cmax,-cmax)
	@test Cpointsvec[2] == (cmax,cmax)

	pointsvec = easyClipper.vector2ntuple(2,Float32,[rand(Float32,2) for x in 1:50])
	minx,xwidth=easyClipper.min_width(pointsvec)
	if xwidth[1]>0 && xwidth[2]>0
		Cpointsvec = easyClipper._scale_j2c(pointsvec,minx,xwidth)
		minpoint,maxpoint=easyClipper.tupleextrema(Cpointsvec)
		@test minpoint == (-cmax,-cmax)
		@test maxpoint == (cmax,cmax)

		Rpointsvec = easyClipper._scale_c2j(Cpointsvec,minx,xwidth)
		@test all( all(isapprox.(x,y)) for (x,y) in zip(pointsvec,Rpointsvec) )
	end

	pointsvec = [(0x01,0x03),(0xd3,0x2f),(0x39,0x22),(0x08,0x34)]
	minx,xwidth=easyClipper.min_width(pointsvec)
	Cpointsvec = easyClipper._scale_j2c(pointsvec,minx,xwidth)
	minpoint,maxpoint=easyClipper.tupleextrema(Cpointsvec)
	@test minpoint == (-cmax,-cmax)
	@test maxpoint == (cmax,cmax)

	Rpointsvec = easyClipper._scale_c2j(Cpointsvec,minx,xwidth)
	@test all( all(isapprox.(x,y)) for (x,y) in zip(pointsvec,Rpointsvec) )

end


text:
# Convergence checks
#
# for single inputs
function check_convergence(
    μ::AbstractVector,
    u::AbstractVector,
    Kv::AbstractVector,
    cache::SinkhornConvergenceCache,
    atol::Real,
    rtol::Real,
)
    # unpack
    tmp = cache.tmp
    norm_μ = cache.norm_source

    # do not overwrite `Kv` but reuse it for computing `u` if not converged
    tmp .= u .* Kv
    norm_uKv = sum(abs, tmp)
    tmp .= abs.(μ .- tmp)
    norm_diff = sum(tmp)

    isconverged = norm_diff < max(atol, rtol * max(norm_μ, norm_uKv))

    return isconverged, norm_diff
end

# for batches
function check_convergence(
    μ::AbstractVecOrMat,
    u::AbstractMatrix,
    Kv::AbstractMatrix,
    cache::SinkhornBatchConvergenceCache,
    atol::Real,
    rtol::Real,
)
    # unpack
    tmp = cache.tmp
    tmp2 = cache.tmp2
    norm_μ = cache.norm_source
    norm_uKv = cache.norm_uKv
    norm_diff = cache.norm_diff
    isconverged = cache.isconverged

    # do not overwrite `Kv` but reuse it for computing `u` if not converged
    tmp .= u .* Kv
    tmp2 .= abs.(tmp)
    sum!(norm_uKv, tmp2)
    tmp .= abs.(μ .- tmp)
    sum!(norm_diff, tmp)

    # check stopping criterion
    @. isconverged = norm_diff < max(atol, rtol * max(norm_μ, norm_uKv))

    return all(isconverged), norm_diff
end

# Common solve! operation
function solve!(solver::Union{SinkhornSolver,SinkhornBarycenterSolver})
    # unpack solver
    atol = solver.atol
    rtol = solver.rtol
    maxiter = solver.maxiter
    check_convergence = solver.check_convergence
    cache = solver.cache
    convergence_cache = solver.convergence_cache

    isconverged = false
    to_check_step = check_convergence
    # initial step if needed
    init_step!(solver)
    for iter in 1:maxiter
        # computations before the Sinkhorn iteration (e.g., absorption step)
        prestep!(solver, iter)
        # perform Sinkhorn iteration
        step!(solver, iter)

        # check source marginal
        # always check convergence after the final iteration
        to_check_step -= 1
        if to_check_step == 0 || iter == maxiter
            # reset counter
            to_check_step = check_convergence

            isconverged, abserror = OptimalTransport.check_convergence(solver)
            @debug string(solver.alg) *
                   " (" *
                   string(iter) *
                   "/" *
                   string(maxiter) *
                   ": absolute error of source marginal = " *
                   string(maximum(abserror))

            if isconverged
                @debug "$(solver.alg) ($iter/$maxiter): converged"
                break
            end
        end
    end

    if !isconverged
        @warn "$(solver.alg) ($maxiter/$maxiter): not converged"
    end

    return nothing
end


text:
@doc """ query a gene's synonym including itself with a gene name
""" ->
function query_gene(genename::ASCIIString)
    if !isdefined(GeneMisc, :gene_id) || !isdefined(GeneMisc,:id_genes)
        load_index()
    end
    id_genes[gene_id[genename]]
end

@doc """ query synonym for genes
""" ->
function query_gene(genenames::Array{ASCIIString,1})
    Dict(zip(genenames, map(query_gene, genenames)))
end


text:
@testset "Test equality method" begin
    model = GLPK.Optimizer()
    x = MOI.add_variables(model, 2)

    con = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0; 1.0], x), 0.0),
            MOI.LessThan(1.0)
          )

    con_set = MOI.get(model, MOI.ConstraintSet(), con)

    @test is_equality(con_set) == false
end


text:
"
This function performs the true jump method for piecewise deterministic Markov processes. It takes the following arguments:

- **initial_state** : a `Vector` of `Int64`, representing the initial states of the system.
- **F** : a `Function` or a callable type, which itself takes three arguments; x, a `Vector` of `Int64` representing the states, parameters, a `Vector` of `Float64` representing the parameters of the system, and t, a `Float64` representing the time of the system.
- **nu** : a `Matrix` of `Int64`, representing the transitions of the system, organised by row.
- **parameters** : a `Vector` of `Float64` representing the parameters of the system.
- **tf** : the final simulation time (`Float64`).
"
function truejump(initial_state::AbstractVector{Int64},F::Base.Callable,nu::AbstractMatrix{Int64},parameters::AbstractVector{Float64},tf::Float64)
    # Args
    args = SSAArgs(initial_state,F,nu,parameters,tf,:tjm,true)
    # Set up time array
    ta = Vector{Float64}()
    t = 0.0
    push!(ta,t)
    # Set up initial x
    nstates = length(initial_state)
    x = initial_state'
    xa = copy(initial_state)
    # Number of propensity functions
    numpf = size(nu,1)
    # Main loop
    termination_status = "finaltime"
    nsteps = 0
    while t <= tf
        ds = rand(Exponential(1.0))
        f = (u)->(quadgk((u)->sum(F(x,parameters,u)),t,u)[1]-ds)
        newt = fzero(f,t)
        if newt>tf
          break
        end
        t=newt
        pf = F(x,parameters,t)
        # Update time
        sumpf = sum(pf)
        if sumpf == 0.0
            termination_status = "zeroprop"
            break
        end
        push!(ta,t)
        # Update event
        ev = pfsample(pf,sumpf,numpf)
        if x isa SVector
            @inbounds x[1] += nu[ev,:]
        else
            deltax = view(nu,ev,:)
            for i in 1:nstates
                @inbounds x[1,i] += deltax[i]
            end
        end
        for xx in x
            push!(xa,xx)
        end
        # update nsteps
        nsteps += 1
    end
    stats = SSAStats(termination_status,nsteps)
    xar = transpose(reshape(xa,length(x),nsteps+1))
    return SSAResult(ta,xar,stats,args)
end


text:
using Documenter, PressureDrop

isCI = get(ENV, "CI", nothing) == "true" #Travis populates this env variable by default

makedocs(
    clean = true,
    pages = [   "Overview" => "index.md",
                "Core functions" => "core.md", #includes types
                "Plotting" => "plotting.md",
                "Pressure & temperature correlations" => "correlations.md",
                "PVT properties" => "pvt.md",
                "Valve calculations" => "valves.md",
                "Utilities" => "utilities.md",
                "Extending" => "extending.md",
                "Similar tools" => "similartools.md"],
    sitename="PressureDrop.jl",
    format = Documenter.HTML(prettyurls = isCI)
)

#tag versions before pushing
if isCI
    deploydocs(repo = "github.com/jnoynaert/PressureDrop.jl.git")
end


text:
function inertiaforces!(node::Int64, params::Parameters)
  om2 = params.omega^2

  for i in 2:3
    for j in 1:6
      if j <= 3
        params.f[j, i, node] = params.f[j, i, node] +
          om2 * params.nodes[node].massinertia * params.del[j, i, node]
      elseif j <= 5
        params.f[j, i, node] =  params.f[j, i, node] +
          om2 * params.nodes[node].radialinertia * params.del[j, i, node]
      else
        params.f[j, i, node] =  params.f[j, i, node] +
          om2 * params.nodes[node].axialinertia * params.del[j, i, node]
      end
    end
  end
end


text:
# by @walra356
using CairoMakie
CairoMakie.activate!() #HIDE
let
    supertitle = "heatmap of camera chip exposure with x and y sections"
  # define transformation from pixel coordinates to physical coordinates
    edges(p,Δx,x0) = collect(p .* Δx) .-(x0+0.5Δx)
  # parameters camera chip (pixel coordinates)
    dimX = 40; dimY = 30                       # pixel dimensions
    (px, py) = (1:dimX, 1:dimY)                  # pixel positions
    (ix, iy) = (10, 3)                        # selected pixel
  # physical coordinates
    (Δx, Δy) = (4.0, 5,0)                       # pixel lattice periods (μm)
    (Ox, Oy) = (edges(ix,Δx,0), edges(iy,Δy,0)) # manual offset (μm)
    (x, y) = (edges(px,Δx,Ox), edges(py,Δy,Oy)) # pixel positions (μm)
  # create exposure matrix σ (model alternative for measurement file)
    valmax = 500.0                              # maximum value exposure matrix
    exposure(i, j , ix, iy, wx, wy) = 0.9valmax *
                        (exp(-(((i-ix)/wx)^2 + ((Δy/Δx)*(j-iy)/wy)^2)) + 0.1)
    σ = round.(Int, [exposure(i, j, ix, iy, 6, 6) for i=1:dimX, j=1:dimY])
    σ1 = σ[px,py]  # raw matrix
    sh = σ[:,iy]   # horizontal section (pixel row iy from bottom)
    sv = σ[ix,:]   # vertical section (pixel column ix from left)
    footnote = "Origin of coordinate system was shifted by ($Ox μm, $Oy μm)
                w.r.t. left-bottom edge of camera chip
                \nPixel lattice periods: (Δx = $Δx μm, Δy = $Δy μm);
                pixel position of light spot: (ix, y0) = ($(ix), $(iy))"

  # start actual plot
    theme = Theme(fontsize = 12, colormap = :gist_earth, resolution = (750,850))
    set_theme!(theme)
    fig = Figure()
  # collect attributes
    fsize = fig.scene.attributes.fontsize.val
    fresx = fig.scene.attributes.resolution.val[1]
    fresy = fig.scene.attributes.resolution.val[2]
    attr = (xlabelsize = 6fsize/5, ylabelsize = 6fsize/5, titlesize = 7fsize/5,
            xautolimitmargin = (0,0), yautolimitmargin = (0,0), )
    attr1 = (attr..., title = "physical scale",
        xticks = [x[1]-Δx/2, x[ix], 50, x[end]+Δx/2],
        yticks = [y[1]-Δy/2, y[iy], 50, y[end]+Δy/2],
        xlabel = "x (μm)", ylabel = "y (μm)",
        aspect = (Δx * size(σ1,1))/(Δy * size(σ1,2)), )
    attrc = (label = "exposure (counts)", ticksize=6fsize/5, tickalign = 1,
        width = 15, height = Relative(0.88), ) # tweaked
    attrh = (attr..., title = "horizontal section",
        xlabel = "x (μm)", ylabel = "exposure (counts)",
        xticks = [x[1]-Δx/2, x[ix], 50, x[end]+Δx/2],
        yticks = [0, 250, valmax],
        limits = ((x[1]-Δx/2, x[end]+Δx/2),(0, valmax)),
        height = fresy/5, aspect = 2.05, ) # tweaked
    attrv = (attr..., xlabel = "exposure (counts)",
        xticks = [0, 250, valmax],
        yticks = [y[1]-Δy/2, y[iy], 50, y[end]+Δy/2],
        xreversed = true, yaxisposition = :right,
        limits = ((0,valmax),(y[1]-Δy/2, y[end]+Δy/2)),
        width = fresy/5, aspect = 0.58, ) # tweaked
  # create axes and add contents
    ax1 = Axis(fig; attr1...)
    axh = Axis(fig; attrh...)
    axv = Axis(fig; attrv...)
    hm1 = heatmap!(ax1, x, y, σ1)
    axc = Colorbar(fig, hm1; attrc...)
    lfn = Label(fig, text = footnote, textsize = 6fsize/5)
    lst = Label(fig, text = supertitle, textsize = 2fsize)
    lvs = Label(fig, text = "vertical section", padding = (0, 5, 0, 0),
        rotation = pi/2, textsize = 7fsize/5)
    sch = scatter!(axh, x, sh, direction = :x, color = :red, markersize = 5, )
    scv = scatter!(axv, sv, y, color = :red, markersize = 5 )
    lnh = lines!(ax1, x , fill(0, length(x)), color = :white, linewidth = 1)
    lnv = lines!(ax1, fill(0, length(y)), y, color = :white, linewidth = 1)
  # create layout and show figure
    fig[1,2] = axh
    fig[2,1] = axv
    fig[2,1, Left()] = lvs
    fig[2,2] = ax1
    fig[2,3] = axc
    fig[3,:] = lfn
    fig[0,:] = lst
    fig
    save(joinpath(@__DIR__, "output", "heatmapScaleSections.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Base.LinAlg.LAPACK.getrf!
using Base.LinAlg.LAPACK.getrs!
using Base.LinAlg.BLAS.gemm!
using Base.LinAlg.BLAS.gemv!
using ArrayViews

abstract LinearSolver

abstract StaticCondensation <: LinearSolver
abstract DirectStaticCondensation <: StaticCondensation

type CholeskySolver <: DirectStaticCondensation

    nel::Int
    nb::Int
    nbslv::Int
    bwd::Int

<<<<<<< HEAD
    nbe::Int
    nie::Int

    A2::Array{Float64, 2}
    X::Array{Float64,3}
    iAii::Array{Float64,2}

    function CholekySolver(bmap::Array{Int,2}, nb, nbslv, Q)
        nel = size(bmap,2)
        bwd = bandwidth(bmap, nbslv)

        nbe = 4 + 4*(Q-2)
        nie = Q*Q-nbe

        A2 = zeros(bwd+1, nbslv)
        X = zeros(nie, nbe, nel)
        iAii = zeros(div(nie*(nie+1),2), nel)
        new(bmap, nel, nb, nbslv, bwd, nbe, nie, A2, X, iAii)
    end
=======
    neb::Int
    nei::Int

    A::Array{Float64,2}
    Aii::Array{Float64,2}
    X::Array{Float64,3}
>>>>>>> ef030efb9656085dd4f7e05f695cbdb6bdca1eb0
end

function sym2pckd!(A, AP, uplo='L')
    n = size(A,1)
    np = div(n*(n+1),2)

    #AP = zeros(eltype(A), np)

    if uplo=='U'
        count = 1
        for i = 1:n
            for j = 1:i
                AP[count] = A[j,i]
                count += 1
            end
        end
    else
        count = 1
        for i = 1:n
            for j = i:n
                AP[count] = A[i,j]
                count += 1
            end
        end
    end

    return AP
end

using ArrayViews
using XLAPACK
using Base.LinAlg.BLAS
function add_matrix(solver::CholeskySolver, bmap, e, Abb, Aib, Aii)

    Aiip = view(solver.Aii, :, e)
    sym2pckd!(Aii, Aiip, 'L')
    pptrf!('L', Aiip)
    X = view(solver.X, :, :, e)

    copy!(X, Aib)
    pptrs!('L', Aiip, X)

    gemm!('T', 'N', -1.0, Aib, X, 1.0, Abb)

    assemble(solver, Abb)

<<<<<<< HEAD
function add_matrix(solver::CholeskySolver, e, Abb, Aib, Aii)

    iAii = view(solver.iAii, :, e)
    X = view(solver.X, :, :, e)

    sym2pckd!(Aii, iAii, 'L')

    tmp,info = pptrf!('L', iAii)

    copy!(X, Aib)

    pptrs!('L', iAii, X)

    gemm!('T', 'N', -1.0, X, Aib, 1.0, Abb)

    assemble(solver, e, Abb)
end

function assemble(solver::ChokeskySolver, e, Ae)
=======
>>>>>>> ef030efb9656085dd4f7e05f695cbdb6bdca1eb0

end

function assemble(solver::ChokeskySolver, bmap, A)

    nei = solver.nei
    neb = solver.neb
    Ag = solver.A
    nbslv = solver.nbslv
    for i = 1:neb
        ig = bmap[i]
        if ig > nbslv
            continue
        end
        for k = 1:neb
            kg = bmap[k]
            if kg <= nbslv && kg >= ig
                if (kg >= ig)
                    Ag[kg-ig+1][ig] += A[k, i]
                end
            end
        end
    end

end

    bm = view(solver.bmap, :, 2)
    nbe = solver.nbe

    for i = 1:nbe
        ig = bm[i]
        if ig <= nblsv
            for j = 1:nbe
                jg = bm[j]
                if jg >= ig
                    solver.A2[ig, jg-ig+1] += A[j,i]
                end
            end
        end
    end

end

function trf(solver::CholeskySolver)

    pbtrf('L', solver.bwd, solver.A2)

end

function static_cond(solver::CholeskySolver, e, fbi)

    nbe = solver.nbe
    nbi = solver.nie

    fb = view(fbi, 1:nbe)
    fi = view(fbi, (nbe+1):(nbe+nie))

    gemv('T', -1.0, view(solver.X, :, :, e), fi, 1.0, fb)

end

function solve!(solver::Cholesky, xbi)

    bmap = solver.bmap
    nbslv = solver.nbslv
    bwd = solver.bwd

    xb = view(xbi, 1:nbslv)
    pbtrs!('L', solver.bwd, solver.A2, xb)

    #xbie = zeros(solver.nbe


text:
# Solving Ordinary Differential Equations I: Nonstiff Problems
# by Ernst Hairer, Gerhard Wanner, and Syvert P Norsett.
# III.5 Variable Step Size Multistep Methods: Formulae 5.9
function ϕ_and_ϕstar!(cache, du, k)
  @inbounds begin
    @unpack dts, ϕstar_nm1, ϕ_n, ϕstar_n, β = cache
    ξ = dt = dts[1]
    ξ0 = zero(dt)
    β[1] = one(dt)
    if typeof(cache) <: OrdinaryDiffEqMutableCache
      ϕ_n[1] .= du
      ϕstar_n[1] .= du
    else
      ϕ_n[1] = du
      ϕstar_n[1] = du
    end
    for i = 2:k
      ξ0 += dts[i]
      β[i] = β[i-1] * ξ/ξ0
      ξ += dts[i]
      if typeof(cache) <: OrdinaryDiffEqMutableCache
        @.. broadcast=false ϕ_n[i] = ϕ_n[i-1] - ϕstar_nm1[i-1]
        @.. broadcast=false ϕstar_n[i] = β[i] * ϕ_n[i]
      else
        ϕ_n[i] = ϕ_n[i-1] - ϕstar_nm1[i-1]
        ϕstar_n[i] = β[i] * ϕ_n[i]
      end
    end
  end # inbounds
end

function ϕ_and_ϕstar!(cache::Union{VCABMConstantCache,VCABMCache}, du, k)
  @inbounds begin
    @unpack dts, ϕstar_nm1, ϕ_n, ϕstar_n, β = cache
    ξ = dt = dts[1]
    ξ0 = zero(dt)
    β[1] = one(dt)
    if typeof(cache) <: OrdinaryDiffEqMutableCache
      ϕ_n[1] .= du
      ϕstar_n[1] .= du
    else
      ϕ_n[1] = du
      ϕstar_n[1] = du
    end
    for i = 2:k
      ξ0 += dts[i]
      β[i] = β[i-1] * ξ/ξ0
      ξ += dts[i]
      if typeof(cache) <: OrdinaryDiffEqMutableCache
        @.. broadcast=false ϕ_n[i] = ϕ_n[i-1] - ϕstar_nm1[i-1]
        @.. broadcast=false ϕstar_n[i] = β[i] * ϕ_n[i]
      else
        ϕ_n[i] = ϕ_n[i-1] - ϕstar_nm1[i-1]
        ϕstar_n[i] = β[i] * ϕ_n[i]
      end
    end
    cache.ξ = ξ
    cache.ξ0 = ξ0
  end # inbounds
end

function expand_ϕ_and_ϕstar!(cache, i)
  @unpack ξ, ξ0, β, dts, ϕstar_nm1, ϕ_n, ϕstar_n = cache
  ξ0 += dts[i]
  β[i] = β[i-1] * ξ/ξ0
  if typeof(cache) <: OrdinaryDiffEqMutableCache
    @.. broadcast=false ϕ_n[i] = ϕ_n[i-1] - ϕstar_nm1[i-1]
    @.. broadcast=false ϕstar_n[i] = β[i] * ϕ_n[i]
  else
    ϕ_n[i] = ϕ_n[i-1] - ϕstar_nm1[i-1]
    ϕstar_n[i] = β[i] * ϕ_n[i]
  end
end

function ϕ_np1!(cache, du_np1, k)
  @inbounds begin
    @unpack ϕ_np1, ϕstar_n = cache
    for i = 1:k
      if i != 1
        if typeof(cache) <: OrdinaryDiffEqMutableCache
          @.. broadcast=false ϕ_np1[i] = ϕ_np1[i-1] - ϕstar_n[i-1]
        else
          ϕ_np1[i] = ϕ_np1[i-1] - ϕstar_n[i-1]
        end
      else
        if typeof(cache) <: OrdinaryDiffEqMutableCache
          ϕ_np1[i] .= du_np1
        else
          ϕ_np1[i] = du_np1
        end
      end
    end
  end #inbounds
end

# Solving Ordinary Differential Equations I: Nonstiff Problems
# by Ernst Hairer, Gerhard Wanner, and Syvert P Norsett.
# III.5 Variable Step Size Multistep Methods: Formulae 5.9 & 5.10
# ------------------------------------------------------------
# Note that `g` is scaled by `dt` in here
function g_coefs!(cache, k)
  @inbounds begin
    @unpack dts,c,g = cache
    ξ = dt = dts[1]
    for i = 1:k
      if i > 2
        ξ += dts[i-1]
      end
      for q = 1:k-(i-1)
        if i > 2
          c[i,q] =  muladd(-dt/ξ, c[i-1,q+1], c[i-1,q])
        elseif i == 1
          c[i,q] = inv(q)
        elseif i == 2
          c[i,q] = inv(q*(q+1))
        end
      end # q
      g[i] = c[i,1] * dt
    end # i
  end # inbounds
end

# Coefficients for the implicit Adams methods
global const γstar = [1,-1/2,-1/12,-1/24,-19/720,-3/160,-863/60480,-275/24192,-33953/3628800,-0.00789255,-0.00678585,-0.00592406,-0.00523669,-0.0046775,-0.00421495,-0.0038269]


text:
# Run the model for the ensemble of solutions -
include("SolveBalances.jl")
include("Common.jl")
using PyPlot
using PyCall
@pyimport numpy as np

# Load the data -
# 1 mg/ml
# MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2e_C3a.txt",',')
# MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3c_C5a.txt",',')

# Control -
# MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2a_C3a.txt",',')
# MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3ai_C5a_original.txt",',')

# 0.1 mg/ml -
# MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2d_C3a.txt",',')
# MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3b_C5a_original.txt",',')

# 0.01 mg/ml -
# MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2c_C3a.txt",',')
# MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3aiii_C5a_original.txt",',')

# 0.001 mg/ml -
MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2b_C3a.txt",',')
MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3aii_C5a_original.txt",',')

# Load the ensemble files from disk -
pc_array_full = readdlm("./data/pc_array_O1_O2.dat")
ec_array_full = readdlm("./data/ec_array_O1_O2.dat")

# Find lowest error solution -
total_error = sum(ec_array_full[:,1:end],1)

# Which col is the min error?
max_index = indmax(total_error)

# min error parameter set -
max_error_parameter_set = pc_array_full[:,max_index]
number_of_parameters = length(max_error_parameter_set)

# Generate N = 100 new *random* parameter family -
number_of_parameter_sets = 100
coefficient_of_variation = 0.29
random_parameter_family = zeros(number_of_parameters,number_of_parameter_sets)
for parameter_set_index = 1:number_of_parameter_sets
  for parameter_index = 1:number_of_parameters
    mean_value = max_error_parameter_set[parameter_index]
    random_parameter_family[parameter_index,parameter_set_index] = rand_normal(mean_value,coefficient_of_variation*mean_value)
  end
end

# Setup time scale -
tStart = 0.0;
tStop = 25.0;
tStep = 0.1;
number_of_timesteps = 200
time_experimental = linspace(tStart,tStop,number_of_timesteps)

# initialize data array -
data_array = zeros(length(time_experimental),1)

data_dictionary = Dict()
initial_condition_array = zeros(19)
initial_condition_array[1] = 1.0      # 1 zymosan
initial_condition_array[2] = 1.9      # 2 C4
initial_condition_array[3] = 0.322    # 3 C2
initial_condition_array[8] = 7.57     # 4 C3
initial_condition_array[14] = 0.195   # 5 C5
initial_condition_array[17] = 2.23    # 6 Factor H
initial_condition_array[18] = 0.417   # 7 C4BP

initial_condition_array[19] = 0.55     # 19 C3a (ic initial)
initial_condition_array[15] = 0.006   # 19 C3a (ic initial)
data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

# Setup the data set -
aic_data_set = Set{Array{Float64,1}}()

# Run the simulation for this rank -
number_of_samples = 1
for index = 1:number_of_parameter_sets

  # Grab the parameter set from the cache -
  parameter_array = random_parameter_family[:,index]
  data_dictionary["PARAMETER_ARRAY"] = parameter_array

  # How many parameters?
  number_of_parameters = length(parameter_array)

  # Run the model -
  (T,X) = SolveBalances(tStart,tStop,tStep,data_dictionary)

  # Calculate the error -
  IC3a = np.interp(MEASURED_ARRAY_C3a[:,1],T,X[:,19])
  IC5a = np.interp(MEASURED_ARRAY_C5a[:,1],T,X[:,15])

  # Compute the error values -
  MC3a = 0.11002*MEASURED_ARRAY_C3a[:,2]
  MC5a = 0.00009615*MEASURED_ARRAY_C5a[:,2]

  # Error -
  error_c3a = (1/norm(MC3a))*sum((IC3a - MC3a).^2)
  error_c5a = (1/norm(MC5a))*sum((IC5a - MC5a).^2)

  # How many steps?
  number_of_steps = length(MC3a)

  # Compute the aic -
  error_array = zeros(2)
  error_array[1] = 2*number_of_parameters+number_of_steps*log(error_c3a)
  error_array[2] = 2*number_of_parameters+number_of_steps*log(error_c5a)

  # Add -
  push!(aic_data_set,error_array)
  @show (index,2*number_of_parameters+number_of_steps*log(error_c3a),2*number_of_parameters+number_of_steps*log(error_c5a))
end

# ok, so lets comput the mean and std -
mean_array_C3 = Float64[]
mean_array_C5 = Float64[]
while (length(aic_data_set)>0)

  data_array = pop!(aic_data_set)
  push!(mean_array_C3,data_array[1])
  push!(mean_array_C5,data_array[2])
end

mean_value_C3 = mean(mean_array_C3)
mean_value_C5 = mean(mean_array_C5)
std_value_C3 = std(mean_array_C3)
std_value_C5 = std(mean_array_C5)

results_array = [mean_value_C3 std_value_C3 ; mean_value_C5 std_value_C5]


text:
using OrdinaryDiffEq, DifferentialEquations
using Flux

tspan = (0.0f0, 10.0f0)

"""
Simple ball in one spatial dimension, for 2 dimensions total:
    v, a
"""
function ball_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    g = p
    Ft = U_func_time(t, u)
    Fs = U_func_space(u)

    du[1] = u[2]
    du[2] = -g
    du .+= Ft .+ Fs
end

"""
Simple uncontrolled ball in one spatial dimension, for 2 dimensions total:
    v, a
Note: non-inplace version
"""
function ball_system(u, p, t)
    g = p

    du = zeros(eltype(p), size(u))
    du[1] = u[2]
    du[2] = -g
    return du
end

"""
Similar to the falling ball system, but with deformations recorded in a third
variable. Comes with nonlinearities, and requires an extra parameter to set
the strength

In words, the three variables now are:
    Position; normal
    Velocity; affected by gravity, increases (to 0) if close to the ground,
            and increases again as deformation decreases back to 0
    Deformation; 0 until near the ground, where it deforms at a certain rate
            and then reforms at an ideally slower rate
"""
function deformed_ball_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    # g, deformation_rate, restoration_rate = p
    g, r, k, damping = p

    Ft = U_func_time(t)
    Fs = U_func_space(u)

    # May allow the ball to go slightly negative
    # r = 0.1; # Radius of ball
    Δ = max(0.0, r - u[1])
    damping_force = Δ>0.0 ? -damping*u[2] : 0.0

    du[1] = u[2] + Fs[1]
    du[2] = -g + k*Δ + damping_force + Fs[2] + Ft

    # ground_closeness = deformation_rate*σ(100.0*(0.01-u[1]));
    # du3 = ground_closeness*(-u[2]) - restoration_rate*u[3] + Fs[3]

    # du[1] = u[2] + Fs[1]
    # du[2] = -g - du3 + Fs[2] + Ft
    # du[3] = du3
    # du[2] = -g + restoration_rate*u[3] + Fs[2] + Ft
    # du[3] = ground_closeness - restoration_rate*u[3] + Fs[3]
end

# Parameters
p1 = 9.81
u1 = [10.0f0, 0.0f0]

# p2 = [9.81, 1000.0, 100.0]
p2 = [9.81, 0.1, 10000.0, 5.0]
u2 = [1.0f0, 0.0f0]

ts = range(tspan[1], tspan[2], length=1000)

function solve_ball_system(;U_func_time = U_func_time_trivial,
                           U_func_space = U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->
                ball_system(du, u, p, t,
                            U_func_time=U_func_time,
                            U_func_space=U_func_space),
                    u1, tspan, p1)
    sol = solve(prob, Tsit5(), saveat=ts);
    return sol
end

function solve_deformed_ball_system(U_func_time = U_func_time_trivial,
                                    U_func_space = U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->
                deformed_ball_system(du, u, p, t,
                            U_func_time=U_func_time,
                            U_func_space=U_func_space),
                    u2, tspan, p2)
    sol = solve(prob, Tsit5(), saveat=ts);
    return Array(sol)
end

# Also define true dynamics in NN syntax
core_dyn_true = Dense(Float32.([0.0 1.0; 0.0 0.0]),
                      Float32.([0.0;-p1]),
                      identity)

export solve_ball_system, ball_system,
        solve_deformed_ball_system, deformed_ball_system,
        core_dyn_true

## Old scratch
# Try to define using a jump system
# bounce_rate(u,p,t) = u[1]<0.0 ? 1.0 : 0.0
# bounce_effect!(integrator) = integrator.u[2] += 10.0
# bounce_jump = VariableRateJump(bounce_rate, bounce_effect!)
#
# prob = ODEProblem(ball_system, u0, tspan, p)
# jump_problem = JumpProblem(prob, Direct(), bounce_jump)
# sol = solve(jump_problem, Tsit5())

# SCRATCH: Bounce using callbacks
# From: http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-1:-Bouncing-Ball-1
# function condition(u,t,integrator) # Event when event_f(u,t) == 0
#   u[1]
# end
# function affect!(integrator) # The actual bounce
#   integrator.u[2] = -integrator.u[2]
# end
# # Note: by default the callback will cause extra saves
# cb = ContinuousCallback(condition,affect!,
#                         save_positions = (false,false))
# prob = ODEProblem(ball_system, u0, tspan, p)
# sol = solve(prob, Tsit5(), saveat=ts, callback=cb)
            #initialize_save=false)
# plot(sol)


text:
@info "@time using GaussLaw1D"
@time using GaussLaw1D
using Test

@testset "Smoketest $shape" for shape in [:parallel, :cylindrical, :spherical]
    prob = Problem(Constant(1.0nC/(mm^3)),
        limits=(1.0mm, 2.0mm),
        potential=(0.0V, 100.0V),
        shape=shape,
    )
    @info "@time solve $shape problem"
    sol = @time solve(prob)
end

@testset "Manufactured solution $shape" for shape in [:parallel, :cylindrical, :spherical]
    phi_true = r-> 1/2*(r/mm)^2*V
    E_true = r -> -(r/mm)*(V/mm)
    ε = GaussLaw1D.VACUUM_PERMITIVITY
    n_div_dims = if shape == :parallel
        1
    elseif shape == :cylindrical
        2
    elseif shape == :spherical
        3
    end
    rho_true = r -> -n_div_dims*ε * (V*mm^-2)

    r1 = 1.0mm
    r2 = 2.0mm
    limits = (1.0mm, 2.0mm)
    p = Problem(rho_true,
        limits=(r1,r2),
        potential=phi_true.(limits),
        shape=shape,
    )

    sol = solve(p)
    @test E_true.(p.radii) ≈ sol.E
    @test phi_true.(p.radii) ≈ sol.potential
    @test rho_true.(p.radii) ≈ sol.rho
end


text:
# This file was generated, do not modify it.

# In this tutorial, we are exploring the application of Ridge and Lasso

using MLJ
import RDatasets: dataset
using PrettyPrinting
MLJ.color_off() # hide
import Distributions
const D = Distributions

@load LinearRegressor pkg=MLJLinearModels
@load RidgeRegressor pkg=MLJLinearModels
@load LassoRegressor pkg=MLJLinearModels

hitters = dataset("ISLR", "Hitters")
@show size(hitters)
names(hitters) |> pprint

y, X = unpack(hitters, ==(:Salary), col->true);

no_miss = .!ismissing.(y)

# And keep only the corresponding features values.
y = collect(skipmissing(y))
X = X[no_miss, :]

# Let's now split our dataset into a train and test sets.
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=424);

using PyPlot
ioff() # hide

figure(figsize=(8,6))
plot(y, ls="none", marker="o")

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14), ylabel("Salary", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g1.svg")) # hide

figure(figsize=(8,6))
hist(y, bins=50, density=true)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Salary", fontsize=14); ylabel("Density", fontsize=14)

edfit = D.fit_mle(D.Exponential, y)
xx = range(minimum(y), 2500, length=100)
yy = pdf.(edfit, xx)
plot(xx, yy, lw=3, label="Exponential distribution fit")

legend(fontsize=12)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g2.svg")) # hide

Xc = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))
scitype(Xc)

model = @pipeline(Standardizer(),
                     OneHotEncoder(),
                     LinearRegressor())

pipe  = machine(model, Xc, y)
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

figure(figsize=(8,6))

res = ŷ .- y[test]
stem(res)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14); ylabel("Residual (ŷ - y)", fontsize=14)

ylim([-1300, 1000])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g3.svg")) # hide

figure(figsize=(8,6))
hist(res, bins=30, density=true, color="green")

xx = range(-1100, 1100, length=100)
ndfit = D.fit_mle(D.Normal, res)
lfit  = D.fit_mle(D.Laplace, res)

plot(xx, pdf.(ndfit, xx), lw=3, color="orange", label="Normal fit")
plot(xx, pdf.(lfit, xx), lw=3, color="magenta", label="Laplace fit")

legend(fontsize=12)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Residual (ŷ - y)", fontsize=14); ylabel("Density", fontsize=14)
xlim([-1100, 1100])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g4.svg")) # hide

pipe.model.linear_regressor = RidgeRegressor()
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

r  = range(model, :(linear_regressor.lambda), lower=1e-2, upper=100_000, scale=:log10)
tm = TunedModel(model=model, ranges=r, tuning=Grid(resolution=50),
                resampling=CV(nfolds=3, rng=4141), measure=rms)
mtm = machine(tm, Xc, y)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.linear_regressor.lambda, sigdigits=4)

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

figure(figsize=(8,6))

res = ŷ .- y[test]
stem(res)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14);
ylabel("Residual (ŷ - y)", fontsize=14)
xlim(1, length(res))

ylim([-1300, 1000])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g5.svg")) # hide

mtm.model.model.linear_regressor = LassoRegressor()
mtm.model.range = range(model, :(linear_regressor.lambda), lower=500, upper=100_000, scale=:log10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.linear_regressor.lambda, sigdigits=4)

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

coefs, intercept = fitted_params(mtm.fitresult).linear_regressor
@show coefs
@show intercept

coef_vals = [c[2] for c in coefs]
sum(coef_vals .≈ 0) / length(coefs)

figure(figsize=(8,6))
stem(coef_vals)

# name of the features including one-hot-encoded ones
all_names = [:AtBat, :Hits, :HmRun, :Runs, :RBI, :Walks, :Years,
             :CAtBat, :CHits, :CHmRun, :CRuns, :CRBI, :CWalks,
             :League__A, :League__N, :Div_E, :Div_W,
             :PutOuts, :Assists, :Errors, :NewLeague_A, :NewLeague_N]

idxshow = collect(1:length(coef_vals))[abs.(coef_vals) .> 10]
xticks(idxshow .- 1, all_names[idxshow], rotation=45, fontsize=12)
yticks(fontsize=12)
ylabel("Amplitude", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g6.svg")) # hide

@load ElasticNetRegressor pkg=MLJLinearModels

mtm.model.model.linear_regressor = ElasticNetRegressor()
mtm.model.range = [range(model, :(linear_regressor.lambda), lower=0.1, upper=100, scale=:log10),
                    range(model, :(linear_regressor.gamma),  lower=500, upper=10_000, scale=:log10)]
mtm.model.tuning = Grid(resolution=10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
@show round(best_mdl.linear_regressor.lambda, sigdigits=4)
@show round(best_mdl.linear_regressor.gamma, sigdigits=4)

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

PyPlot.close_figs() # hide


text:
@testset "agreement/confusion matrix tests" begin
    hard_label_pairs = zip([1, 1, 3, 1, 3, 1, 2, 1], [2, 2, 1, 1, 3, 2, 3, 1])
    c = confusion_matrix(3, hard_label_pairs)
    @test c == [
        2 3 0
        0 0 1
        1 0 1
    ]
    kappa, percent_agreement = cohens_kappa(3, hard_label_pairs)
    chance = Lighthouse._probability_of_chance_agreement(3, hard_label_pairs)
    @test chance == (5 * 3 + 1 * 3 + 2 * 2) / 8^2
    @test accuracy(c) == percent_agreement == 3 / 8
    @test kappa == (3 / 8 - chance) / (1 - chance)
    stats = binary_statistics(c, 3)
    @test stats.predicted_positives == 2
    @test stats.predicted_negatives == 6
    @test stats.actual_positives == 2
    @test stats.actual_negatives == 6
    @test stats.true_positives == 1
    @test stats.true_negatives == 2
    @test stats.false_positives == 1
    @test stats.false_negatives == 1
    @test stats.true_positive_rate == 0.5
    @test stats.true_negative_rate == 1 / 3
    @test stats.false_positive_rate == 1 / 6
    @test stats.false_negative_rate == 0.5
    @test stats.precision == 0.5

    labels = rand(StableRNG(42), 1:3, 100)
    hard_label_pairs = zip(labels, labels)
    c = confusion_matrix(3, hard_label_pairs)
    @test tr(c) == length(labels) == sum(c)
    kappa, percent_agreement = cohens_kappa(3, hard_label_pairs)
    @test accuracy(c) == percent_agreement == 1
    @test kappa == 1
    for i in 1:3
        stats = binary_statistics(c, i)
        @test stats.predicted_positives == stats.true_positives == stats.actual_positives
        @test stats.predicted_negatives == stats.true_negatives == stats.actual_negatives
        @test stats.false_positives == stats.false_negatives == 0
        @test stats.true_positive_rate == stats.true_negative_rate == 1
        @test stats.false_positive_rate == stats.false_negative_rate == 0
        @test stats.precision == 1
    end

    n, k = 1_000_000, 2
    rng = StableRNG(42)
    hard_label_pairs = zip(rand(rng, 1:k, n), rand(rng, 1:k, n))
    c = confusion_matrix(k, hard_label_pairs)
    kappa, percent_agreement = cohens_kappa(3, hard_label_pairs)
    @test percent_agreement == accuracy(c)
    @test isapprox(percent_agreement, 0.5; atol=0.02)
    @test isapprox(kappa, 0.0; atol=0.02)
    stats = binary_statistics(c, 1)
    @test isapprox(stats.predicted_positives, 500_000; atol=2000)
    @test isapprox(stats.predicted_negatives, 500_000; atol=2000)
    @test isapprox(stats.actual_positives, 500_000; atol=2000)
    @test isapprox(stats.actual_negatives, 500_000; atol=2000)
    @test isapprox(stats.true_positives, 250_000; atol=2000)
    @test isapprox(stats.true_negatives, 250_000; atol=2000)
    @test isapprox(stats.false_positives, 250_000; atol=2000)
    @test isapprox(stats.false_negatives, 250_000; atol=2000)
    @test isapprox(stats.true_positive_rate, 0.5; atol=0.02)
    @test isapprox(stats.true_negative_rate, 0.5; atol=0.02)
    @test isapprox(stats.false_positive_rate, 0.5; atol=0.02)
    @test isapprox(stats.false_negative_rate, 0.5; atol=0.02)
    @test isapprox(stats.precision, 0.5; atol=0.02)

    @test confusion_matrix(10, ()) == zeros(10, 10)
    @test all(ismissing, cohens_kappa(10, ()))
    @test ismissing(accuracy(zeros(10, 10)))
    stats = binary_statistics(zeros(10, 10), 1)
    @test stats.predicted_positives == 0
    @test stats.predicted_negatives == 0
    @test stats.actual_positives == 0
    @test stats.actual_negatives == 0
    @test stats.true_positives == 0
    @test stats.true_negatives == 0
    @test stats.false_positives == 0
    @test stats.false_negatives == 0
    @test stats.true_positive_rate == 1
    @test stats.true_negative_rate == 1
    @test stats.false_positive_rate == 0
    @test stats.false_negative_rate == 0
    @test ismissing(stats.precision)

    for p in 0:0.1:1
        @test Lighthouse._cohens_kappa(p, p) == 0
        if p > 0
            @test Lighthouse._cohens_kappa(p / 2, p) < 0
            @test Lighthouse._cohens_kappa(p, p / 2) > 0
        end
    end
end

@testset "`calibration_curve`" begin
    rng = StableRNG(42)
    probs = rand(rng, 1_000_000)
    bitmask = rand(rng, Bool, 1_000_000)
    bin_count = 12
    bins, fractions, totals, mean_squared_error = calibration_curve(probs, bitmask;
                                                                    bin_count=bin_count)
    @test bin_count == length(bins)
    @test first(first(bins)) == 0.0 && last(last(bins)) == 1.0
    @test all(!ismissing, fractions)
    @test all(!iszero, totals)
    @test all(isapprox.(fractions, 0.5; atol=0.02))
    @test all(isapprox.(totals, length(probs) / bin_count; atol=1000))
    @test sum(totals) == length(probs)
    @test isapprox(ceil(mean(fractions) * length(bitmask)), count(bitmask); atol=1)
    @test isapprox(mean_squared_error, inv(bin_count); atol=0.002)

    rng = StableRNG(42)
    probs = range(0.0, 1.0; length=1_000_000)
    bitmask = [rand(rng) <= p for p in probs]
    bin_count = 10
    bins, fractions, totals, mean_squared_error = calibration_curve(probs, bitmask;
                                                                    bin_count=bin_count)
    ideal = range(mean(first(bins)), mean(last(bins)); length=bin_count)
    @test bin_count == length(bins)
    @test first(first(bins)) == 0.0 && last(last(bins)) == 1.0
    @test all(!ismissing, fractions)
    @test all(!iszero, totals)
    @test all(isapprox.(fractions, ideal; atol=0.01))
    @test all(totals .== 1_000_000 / bin_count)
    @test isapprox(ceil(mean(fractions) * length(bitmask)), count(bitmask); atol=1)
    @test isapprox(mean_squared_error, 0.0; atol=0.00001)

    bitmask = reverse(bitmask)
    bins, fractions, totals, mean_squared_error = calibration_curve(probs, bitmask;
                                                                    bin_count=bin_count)
    @test bin_count == length(bins)
    @test first(first(bins)) == 0.0 && last(last(bins)) == 1.0
    @test all(!ismissing, fractions)
    @test all(!iszero, totals)
    @test all(isapprox.(fractions, reverse(ideal); atol=0.01))
    @test all(totals .== 1_000_000 / bin_count)
    @test isapprox(ceil(mean(fractions) * length(bitmask)), count(bitmask); atol=1)
    @test isapprox(mean_squared_error, 1 / 3; atol=0.01)
end


text:
using Statistics

using CUDA.CUDNN:
    cudnnNormalizationForward,
    cudnnNormalizationForward!,
    cudnnNormalizationForwardInference,
    cudnnNormalizationForwardTraining,
    cudnnNormalizationBackward,
    cudnnActivationDescriptor,
    cudnnNormMode_t,
        CUDNN_NORM_PER_ACTIVATION, # 0, bnScale, bnBias tensor dims are 1xCxHxWx.. (one value per CHW...-slice, normalized over N slice)
        CUDNN_NORM_PER_CHANNEL,    # 1, bnScale, bnBias tensor dims are 1xCx1x1 (one value per C-dim normalized over Nx1xHxW subtensors)
    cudnnNormOps_t,
        CUDNN_NORM_OPS_NORM,                # 0, /* do normalization only */
        CUDNN_NORM_OPS_NORM_ACTIVATION,     # 1, /* do Norm, then activation */
        CUDNN_NORM_OPS_NORM_ADD_ACTIVATION, # 2, /* do Norm, then elemWiseAdd, then activation */
    cudnnNormAlgo_t,
        CUDNN_NORM_ALGO_STANDARD, # 0
        CUDNN_NORM_ALGO_PERSIST,  # 1
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,      # 0
        CUDNN_ACTIVATION_RELU,         # 1
        CUDNN_ACTIVATION_TANH,         # 2
        CUDNN_ACTIVATION_CLIPPED_RELU, # 3
        CUDNN_ACTIVATION_ELU,          # 4
        CUDNN_ACTIVATION_IDENTITY,     # 5
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C  # 2, /* each image point is vector of element of C, vector length in data type */

@testset "cudnn/normalization" begin

    function normtest(
        x;

        training = false,

        # Inference parameters:
        z = nothing, # for residual addition to the result of the normalization operation, prior to the activation
        mode::cudnnNormMode_t = CUDNN_NORM_PER_CHANNEL, # Per-channel layer is based on the paper Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift, S. Ioffe, C. Szegedy, 2015.
        normOps::cudnnNormOps_t = CUDNN_NORM_OPS_NORM,  # Currently CUDNN_NORM_OPS_NORM_ACTIVATION and CUDNN_NORM_OPS_NORM_ADD_ACTIVATION are only supported in the NHWC layout (training,backward), not supported (inference)
        algo::cudnnNormAlgo_t = CUDNN_NORM_ALGO_STANDARD, # trigger the new semi-persistent NHWC kernel when CUDNN_NORM_ALGO_PERSIST
        alpha::Real = 1,
        beta::Real = 0,
        epsilon::Real = 1e-5, # Has to be >= 0. Should be the same in forward and backward functions.
        groupCnt::Integer = 1, # Place hold for future work, should be set to 1 now

        # Main argument defaults:
        format::cudnnTensorFormat_t = CUDNN_TENSOR_NCHW, # or NHWC
        _sdims = (mode == CUDNN_NORM_PER_CHANNEL    && format == CUDNN_TENSOR_NCHW ? (1,1,size(x,3),1) :
                  mode == CUDNN_NORM_PER_CHANNEL    && format == CUDNN_TENSOR_NHWC ? (size(x,1),1,1,1) :
                  mode == CUDNN_NORM_PER_ACTIVATION && format == CUDNN_TENSOR_NCHW ? (size(x)[1:3]...,1) :
                  mode == CUDNN_NORM_PER_ACTIVATION && format == CUDNN_TENSOR_NHWC ? (size(x)[1:3]...,1) :
                  error("Unknown mode $mode and format $format")),
        scale = fill!(similar(x, _sdims), 1),
        bias =  fill!(similar(x, _sdims), 0),
        xmean =  fill!(similar(x, _sdims), 0),
        xvar = fill!(similar(x, _sdims), 1),

        # Training-only parameters:
        exponentialAverageFactor::Real = 0.1,
        savedMean = nothing, # Optionally save intermediate results from the forward pass here - can be reused to speed up backward pass. NULL if unused.
        savedInvVariance = nothing,

        # Activation parameters:
        activationMode::cudnnActivationMode_t = CUDNN_ACTIVATION_IDENTITY,
        activationReluNanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
        activationCoef::Real = 1,
        activationDesc::Union{Nothing,cudnnActivationDescriptor} = (normOps == CUDNN_NORM_OPS_NORM ? nothing : cudnnActivationDescriptor(activationMode, activationReluNanOpt, Cdouble(activationCoef))),
    )
        if training
            dims = findall(size(xmean) .== 1)
            m = mean(x; dims)
            v = var(x; dims, mean=m, corrected=false)
            y = bias .+ scale .* (x .- m) ./ sqrt.(epsilon .+ v)
        else
            y = bias .+ scale .* (x .- xmean) ./ sqrt.(epsilon .+ xvar)
        end
        y0 = randn!(similar(x))
        y1 = alpha * y
        y2 = y1 + beta * y0
        @test y1 ≈ cudnnNormalizationForward(x, xmean, xvar, bias, scale; training, z, mode,
                                             normOps, algo, alpha, epsilon, groupCnt,
                                             format, exponentialAverageFactor, savedMean,
                                             savedInvVariance, activationDesc)
        @test y2 ≈ cudnnNormalizationForward!(copy(y0), x, xmean, xvar, bias, scale;
                                              training, z, mode, normOps, algo, alpha, beta,
                                              epsilon, groupCnt, format,
                                              exponentialAverageFactor, savedMean,
                                              savedInvVariance, activationDesc)
    end

    x, z, s = (CUDA.randn(x...) for x in ((5,4,3,2),(5,4,3,2),(1,1,3,1)))
    normtest(x)
    normtest(x; training = true)
    normtest(x; mode = CUDNN_NORM_PER_ACTIVATION)
    normtest(x; algo = CUDNN_NORM_ALGO_PERSIST)
    normtest(x; algo = CUDNN_NORM_ALGO_PERSIST, format = CUDNN_TENSOR_NHWC)
    normtest(x; alpha = 2)
    normtest(x; beta = 2)
    normtest(x; epsilon = 0)
    normtest(x; format = CUDNN_TENSOR_NHWC)
    normtest(x; scale = fill!(s, 2))
    normtest(x; bias  = fill!(s, 2))
    normtest(x; xmean  = fill!(s, 2))
    normtest(x; xvar = fill!(s, 2))
    normtest(x; exponentialAverageFactor = 0.01)
    normtest(x; savedMean = similar(s))
    normtest(x; savedInvVariance = similar(s))
    # cudnn-8.0.5: Currently, CUDNN_NORM_OPS_NORM_ACTIVATION and CUDNN_NORM_OPS_NORM_ADD_ACTIVATION are not supported in inference.
    #normtest(x; normOps = CUDNN_NORM_OPS_NORM_ACTIVATION, activationMode = CUDNN_ACTIVATION_RELU, format = CUDNN_TENSOR_NHWC)
    #normtest(x; normOps = CUDNN_NORM_OPS_NORM_ADD_ACTIVATION, activationMode = CUDNN_ACTIVATION_RELU, z, format = CUDNN_TENSOR_NHWC)
    #normtest(x; groupCnt = 2) # cudnn-8.0.5: Currently only groupCnt=1 is supported
end


text:
using FactCheck
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("Unit tests confirm correct algorithm implementation") do

    context("true things are true") do
        @fact 1 --> 1
    end

    context("errors are thrown when expected") do
        @fact_throws find(foo)
    end

    context("pending facts are pending") do
        @pending foo(12) --> "foo"
    end
end

max_stars_repo_name: milktrader/Tuxedo

max_issues_repo_name: milktrader/Tuxedo

max_forks_repo_name: milktrader/Tuxedo