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

const Chars = Union{AbstractChar,Tuple{Vararg{<:AbstractChar}},AbstractVector{<:AbstractChar},Set{<:AbstractChar}}

# starts with and ends with predicates

"""
    startswith(s::AbstractString, prefix::AbstractString)

Return `true` if `s` starts with `prefix`. If `prefix` is a vector or set
of characters, test whether the first character of `s` belongs to that set.

See also [`endswith`](@ref).

# Examples
```jldoctest
julia> startswith("JuliaLang", "Julia")
true
```
"""
function startswith(a::AbstractString, b::AbstractString)
    a, b = Iterators.Stateful(a), Iterators.Stateful(b)
    all(splat(==), zip(a, b)) && isempty(b)
end
startswith(str::AbstractString, chars::Chars) = !isempty(str) && first(str)::AbstractChar in chars

"""
    endswith(s::AbstractString, suffix::AbstractString)

Return `true` if `s` ends with `suffix`. If `suffix` is a vector or set of
characters, test whether the last character of `s` belongs to that set.

See also [`startswith`](@ref).

# Examples
```jldoctest
julia> endswith("Sunday", "day")
true
```
"""
function endswith(a::AbstractString, b::AbstractString)
    a = Iterators.Stateful(Iterators.reverse(a))
    b = Iterators.Stateful(Iterators.reverse(b))
    all(splat(==), zip(a, b)) && isempty(b)
end
endswith(str::AbstractString, chars::Chars) = !isempty(str) && last(str) in chars

function startswith(a::Union{String, SubString{String}},
                    b::Union{String, SubString{String}})
    cub = ncodeunits(b)
    if ncodeunits(a) < cub
        false
    elseif _memcmp(a, b, sizeof(b)) == 0
        nextind(a, cub) == cub + 1
    else
        false
    end
end

function endswith(a::Union{String, SubString{String}},
                  b::Union{String, SubString{String}})
    cub = ncodeunits(b)
    astart = ncodeunits(a) - ncodeunits(b) + 1
    if astart < 1
        false
    elseif GC.@preserve(a, _memcmp(pointer(a, astart), b, sizeof(b))) == 0
        thisind(a, astart) == astart
    else
        false
    end
end

"""
    contains(haystack::AbstractString, needle)

Return `true` if `haystack` contains `needle`.
This is the same as `occursin(needle, haystack)`, but is provided for consistency with
`startswith(haystack, needle)` and `endswith(haystack, needle)`.

# Examples
```jldoctest
julia> contains("JuliaLang is pretty cool!", "Julia")
true

julia> contains("JuliaLang is pretty cool!", 'a')
true

julia> contains("aba", r"a.a")
true

julia> contains("abba", r"a.a")
false
```

!!! compat "Julia 1.5"
    The `contains` function requires at least Julia 1.5.
"""
contains(haystack::AbstractString, needle) = occursin(needle, haystack)

"""
    endswith(suffix)

Create a function that checks whether its argument ends with `suffix`, i.e.
a function equivalent to `y -> endswith(y, suffix)`.

The returned function is of type `Base.Fix2{typeof(endswith)}`, which can be
used to implement specialized methods.

!!! compat "Julia 1.5"
    The single argument `endswith(suffix)` requires at least Julia 1.5.

# Examples
```jldoctest
julia> endswith_julia = endswith("Julia");

julia> endswith_julia("Julia")
true

julia> endswith_julia("JuliaLang")
false
```
"""
endswith(s) = Base.Fix2(endswith, s)

"""
    startswith(prefix)

Create a function that checks whether its argument starts with `prefix`, i.e.
a function equivalent to `y -> startswith(y, prefix)`.

The returned function is of type `Base.Fix2{typeof(startswith)}`, which can be
used to implement specialized methods.

!!! compat "Julia 1.5"
    The single argument `startswith(prefix)` requires at least Julia 1.5.

# Examples
```jldoctest
julia> startswith_julia = startswith("Julia");

julia> startswith_julia("Julia")
true

julia> startswith_julia("NotJulia")
false
```
"""
startswith(s) = Base.Fix2(startswith, s)

"""
    contains(needle)

Create a function that checks whether its argument contains `needle`, i.e.
a function equivalent to `haystack -> contains(haystack, needle)`.

The returned function is of type `Base.Fix2{typeof(contains)}`, which can be
used to implement specialized methods.
"""
contains(needle) = Base.Fix2(contains, needle)

"""
    chop(s::AbstractString; head::Integer = 0, tail::Integer = 1)

Remove the first `head` and the last `tail` characters from `s`.
The call `chop(s)` removes the last character from `s`.
If it is requested to remove more characters than `length(s)`
then an empty string is returned.

# Examples
```jldoctest
julia> a = "March"
"March"

julia> chop(a)
"Marc"

julia> chop(a, head = 1, tail = 2)
"ar"

julia> chop(a, head = 5, tail = 5)
""
```
"""
function chop(s::AbstractString; head::Integer = 0, tail::Integer = 1)
    if isempty(s)
        return SubString(s)
    end
    SubString(s, nextind(s, firstindex(s), head), prevind(s, lastindex(s), tail))
end

# TODO: optimization for the default case based on
# chop(s::AbstractString) = SubString(s, firstindex(s), prevind(s, lastindex(s)))

"""
    chomp(s::AbstractString) -> SubString

Remove a single trailing newline from a string.

# Examples
```jldoctest
julia> chomp("Hello\\n")
"Hello"
```
"""
function chomp(s::AbstractString)
    i = lastindex(s)
    (i < 1 || s[i] != '\n') && (return SubString(s, 1, i))
    j = prevind(s,i)
    (j < 1 || s[j] != '\r') && (return SubString(s, 1, j))
    return SubString(s, 1, prevind(s,j))
end
function chomp(s::String)
    i = lastindex(s)
    if i < 1 || codeunit(s,i) != 0x0a
        return @inbounds SubString(s, 1, i)
    elseif i < 2 || codeunit(s,i-1) != 0x0d
        return @inbounds SubString(s, 1, prevind(s, i))
    else
        return @inbounds SubString(s, 1, prevind(s, i-1))
    end
end

"""
    lstrip([pred=isspace,] str::AbstractString) -> SubString
    lstrip(str::AbstractString, chars) -> SubString

Remove leading characters from `str`, either those specified by `chars` or those for
which the function `pred` returns `true`.

The default behaviour is to remove leading whitespace and delimiters: see
[`isspace`](@ref) for precise details.

The optional `chars` argument specifies which characters to remove: it can be a single
character, or a vector or set of characters.

# Examples
```jldoctest
julia> a = lpad("March", 20)
"               March"

julia> lstrip(a)
"March"
```
"""
function lstrip(f, s::AbstractString)
    e = lastindex(s)
    for (i::Int, c::AbstractChar) in pairs(s)
        !f(c) && return @inbounds SubString(s, i, e)
    end
    SubString(s, e+1, e)
end
lstrip(s::AbstractString) = lstrip(isspace, s)
lstrip(s::AbstractString, chars::Chars) = lstrip(in(chars), s)

"""
    rstrip([pred=isspace,] str::AbstractString) -> SubString
    rstrip(str::AbstractString, chars) -> SubString

Remove trailing characters from `str`, either those specified by `chars` or those for
which the function `pred` returns `true`.

The default behaviour is to remove trailing whitespace and delimiters: see
[`isspace`](@ref) for precise details.

The optional `chars` argument specifies which characters to remove: it can be a single
character, or a vector or set of characters.

# Examples
```jldoctest
julia> a = rpad("March", 20)
"March               "

julia> rstrip(a)
"March"
```
"""
function rstrip(f, s::AbstractString)
    for (i, c) in Iterators.reverse(pairs(s))
        f(c::AbstractChar) || return @inbounds SubString(s, 1, i::Int)
    end
    SubString(s, 1, 0)
end
rstrip(s::AbstractString) = rstrip(isspace, s)
rstrip(s::AbstractString, chars::Chars) = rstrip(in(chars), s)

"""
    strip([pred=isspace,] str::AbstractString) -> SubString
    strip(str::AbstractString, chars) -> SubString

Remove leading and trailing characters from `str`, either those specified by `chars` or
those for which the function `pred` returns `true`.

The default behaviour is to remove leading whitespace and delimiters: see
[`isspace`](@ref) for precise details.

The optional `chars` argument specifies which characters to remove: it can be a single
character, vector or set of characters.

!!! compat "Julia 1.2"
    The method which accepts a predicate function requires Julia 1.2 or later.

# Examples
```jldoctest
julia> strip("{3, 5}\\n", ['{', '}', '\\n'])
"3, 5"
```
"""
strip(s::AbstractString) = lstrip(rstrip(s))
strip(s::AbstractString, chars::Chars) = lstrip(rstrip(s, chars), chars)
strip(f, s::AbstractString) = lstrip(f, rstrip(f, s))

## string padding functions ##

"""
    lpad(s, n::Integer, p::Union{AbstractChar,AbstractString}=' ') -> String

Stringify `s` and pad the resulting string on the left with `p` to make it `n`
characters (code points) long. If `s` is already `n` characters long, an equal
string is returned. Pad with spaces by default.

# Examples
```jldoctest
julia> lpad("March", 10)
"     March"
```
"""
lpad(s, n::Integer, p::Union{AbstractChar,AbstractString}=' ') = lpad(string(s)::AbstractString, n, string(p))

function lpad(
    s::Union{AbstractChar,AbstractString},
    n::Integer,
    p::Union{AbstractChar,AbstractString}=' ',
) :: String
    n = Int(n)::Int
    m = signed(n) - Int(length(s))::Int
    m ≤ 0 && return string(s)
    l = length(p)
    q, r = divrem(m, l)
    r == 0 ? string(p^q, s) : string(p^q, first(p, r), s)
end

"""
    rpad(s, n::Integer, p::Union{AbstractChar,AbstractString}=' ') -> String

Stringify `s` and pad the resulting string on the right with `p` to make it `n`
characters (code points) long. If `s` is already `n` characters long, an equal
string is returned. Pad with spaces by default.

# Examples
```jldoctest
julia> rpad("March", 20)
"March               "
```
"""
rpad(s, n::Integer, p::Union{AbstractChar,AbstractString}=' ') = rpad(string(s)::AbstractString, n, string(p))

function rpad(
    s::Union{AbstractChar,AbstractString},
    n::Integer,
    p::Union{AbstractChar,AbstractString}=' ',
) :: String
    n = Int(n)::Int
    m = signed(n) - Int(length(s))::Int
    m ≤ 0 && return string(s)
    l = length(p)
    q, r = divrem(m, l)
    r == 0 ? string(s, p^q) : string(s, p^q, first(p, r))
end

"""
    split(str::AbstractString, dlm; limit::Integer=0, keepempty::Bool=true)
    split(str::AbstractString; limit::Integer=0, keepempty::Bool=false)

Split `str` into an array of substrings on occurrences of the delimiter(s) `dlm`.  `dlm`
can be any of the formats allowed by [`findnext`](@ref)'s first argument (i.e. as a
string, regular expression or a function), or as a single character or collection of
characters.

If `dlm` is omitted, it defaults to [`isspace`](@ref).

The optional keyword arguments are:
 - `limit`: the maximum size of the result. `limit=0` implies no maximum (default)
 - `keepempty`: whether empty fields should be kept in the result. Default is `false` without
   a `dlm` argument, `true` with a `dlm` argument.

See also [`rsplit`](@ref).

# Examples
```jldoctest
julia> a = "Ma.rch"
"Ma.rch"

julia> split(a, ".")
2-element Vector{SubString{String}}:
 "Ma"
 "rch"
```
"""
function split end

function split(str::T, splitter;
               limit::Integer=0, keepempty::Bool=true) where {T<:AbstractString}
    _split(str, splitter, limit, keepempty, T <: SubString ? T[] : SubString{T}[])
end
function split(str::T, splitter::Union{Tuple{Vararg{<:AbstractChar}},AbstractVector{<:AbstractChar},Set{<:AbstractChar}};
               limit::Integer=0, keepempty::Bool=true) where {T<:AbstractString}
    _split(str, in(splitter), limit, keepempty, T <: SubString ? T[] : SubString{T}[])
end
function split(str::T, splitter::AbstractChar;
               limit::Integer=0, keepempty::Bool=true) where {T<:AbstractString}
    _split(str, isequal(splitter), limit, keepempty, T <: SubString ? T[] : SubString{T}[])
end

function _split(str::AbstractString, splitter::F, limit::Integer, keepempty::Bool, strs::Vector) where F
    # Forcing specialization on `splitter` improves performance (roughly 30% decrease in runtime)
    # and prevents a major invalidation risk (1550 MethodInstances)
    i = 1 # firstindex(str)
    n = lastindex(str)::Int
    r = findfirst(splitter,str)::Union{Nothing,Int,UnitRange{Int}}
    if !isnothing(r)
        j, k = first(r), nextind(str,last(r))::Int
        while 0 < j <= n && length(strs) != limit-1
            if i < k
                if keepempty || i < j
                    push!(strs, @inbounds SubString(str,i,prevind(str,j)::Int))
                end
                i = k
            end
            (k <= j) && (k = nextind(str,j)::Int)
            r = findnext(splitter,str,k)::Union{Nothing,Int,UnitRange{Int}}
            isnothing(r) && break
            j, k = first(r), nextind(str,last(r))::Int
        end
    end
    if keepempty || i <= ncodeunits(str)::Int
        push!(strs, @inbounds SubString(str,i))
    end
    return strs
end

# a bit oddball, but standard behavior in Perl, Ruby & Python:
split(str::AbstractString;
      limit::Integer=0, keepempty::Bool=false) =
    split(str, isspace; limit=limit, keepempty=keepempty)

"""
    rsplit(s::AbstractString; limit::Integer=0, keepempty::Bool=false)
    rsplit(s::AbstractString, chars; limit::Integer=0, keepempty::Bool=true)

Similar to [`split`](@ref), but starting from the end of the string.

# Examples
```jldoctest
julia> a = "M.a.r.c.h"
"M.a.r.c.h"

julia> rsplit(a, ".")
5-element Vector{SubString{String}}:
 "M"
 "a"
 "r"
 "c"
 "h"

julia> rsplit(a, "."; limit=1)
1-element Vector{SubString{String}}:
 "M.a.r.c.h"

julia> rsplit(a, "."; limit=2)
2-element Vector{SubString{String}}:
 "M.a.r.c"
 "h"
```
"""
function rsplit end

function rsplit(str::T, splitter;
                limit::Integer=0, keepempty::Bool=true) where {T<:AbstractString}
    _rsplit(str, splitter, limit, keepempty, T <: SubString ? T[] : SubString{T}[])
end
function rsplit(str::T, splitter::Union{Tuple{Vararg{<:AbstractChar}},AbstractVector{<:AbstractChar},Set{<:AbstractChar}};
                limit::Integer=0, keepempty::Bool=true) where {T<:AbstractString}
    _rsplit(str, in(splitter), limit, keepempty, T <: SubString ? T[] : SubString{T}[])
end
function rsplit(str::T, splitter::AbstractChar;
                limit::Integer=0, keepempty::Bool=true) where {T<:AbstractString}
    _rsplit(str, isequal(splitter), limit, keepempty, T <: SubString ? T[] : SubString{T}[])
end

function _rsplit(str::AbstractString, splitter, limit::Integer, keepempty::Bool, strs::Array)
    n = lastindex(str)::Int
    r = something(findlast(splitter, str)::Union{Nothing,Int,UnitRange{Int}}, 0)
    j, k = first(r), last(r)
    while j > 0 && k > 0 && length(strs) != limit-1
        (keepempty || k < n) && pushfirst!(strs, @inbounds SubString(str,nextind(str,k)::Int,n))
        n = prevind(str, j)::Int
        r = something(findprev(splitter,str,n)::Union{Nothing,Int,UnitRange{Int}}, 0)
        j, k = first(r), last(r)
    end
    (keepempty || n > 0) && pushfirst!(strs, SubString(str,1,n))
    return strs
end
rsplit(str::AbstractString;
      limit::Integer=0, keepempty::Bool=false) =
    rsplit(str, isspace; limit=limit, keepempty=keepempty)

_replace(io, repl, str, r, pattern) = print(io, repl)
_replace(io, repl::Function, str, r, pattern) =
    print(io, repl(SubString(str, first(r), last(r))))
_replace(io, repl::Function, str, r, pattern::Function) =
    print(io, repl(str[first(r)]))

replace(str::String, pat_repl::Pair{<:AbstractChar}; count::Integer=typemax(Int)) =
    replace(str, isequal(first(pat_repl)) => last(pat_repl); count=count)

replace(str::String, pat_repl::Pair{<:Union{Tuple{Vararg{<:AbstractChar}},
                                            AbstractVector{<:AbstractChar},Set{<:AbstractChar}}};
        count::Integer=typemax(Int)) =
    replace(str, in(first(pat_repl)) => last(pat_repl), count=count)

_pat_replacer(x) = x
_free_pat_replacer(x) = nothing

function replace(str::String, pat_repl::Pair; count::Integer=typemax(Int))
    pattern, repl = pat_repl
    count == 0 && return str
    count < 0 && throw(DomainError(count, "`count` must be non-negative."))
    n = 1
    e = lastindex(str)
    i = a = firstindex(str)
    pattern = _pat_replacer(pattern)
    r = something(findnext(pattern,str,i), 0)
    j, k = first(r), last(r)
    if j == 0
        _free_pat_replacer(pattern)
        return str
    end
    out = IOBuffer(sizehint=floor(Int, 1.2sizeof(str)))
    while j != 0
        if i == a || i <= k
            GC.@preserve str unsafe_write(out, pointer(str, i), UInt(j-i))
            _replace(out, repl, str, r, pattern)
        end
        if k < j
            i = j
            j > e && break
            k = nextind(str, j)
        else
            i = k = nextind(str, k)
        end
        r = something(findnext(pattern,str,k), 0)
        r === 0:-1 || n == count && break
        j, k = first(r), last(r)
        n += 1
    end
    _free_pat_replacer(pattern)
    write(out, SubString(str,i))
    String(take!(out))
end

"""
    replace(s::AbstractString, pat=>r; [count::Integer])

Search for the given pattern `pat` in `s`, and replace each occurrence with `r`.
If `count` is provided, replace at most `count` occurrences.
`pat` may be a single character, a vector or a set of characters, a string,
or a regular expression.
If `r` is a function, each occurrence is replaced with `r(s)`
where `s` is the matched substring (when `pat` is a `AbstractPattern` or `AbstractString`) or
character (when `pat` is an `AbstractChar` or a collection of `AbstractChar`).
If `pat` is a regular expression and `r` is a [`SubstitutionString`](@ref), then capture group
references in `r` are replaced with the corresponding matched text.
To remove instances of `pat` from `string`, set `r` to the empty `String` (`""`).

# Examples
```jldoctest
julia> replace("Python is a programming language.", "Python" => "Julia")
"Julia is a programming language."

julia> replace("The quick foxes run quickly.", "quick" => "slow", count=1)
"The slow foxes run quickly."

julia> replace("The quick foxes run quickly.", "quick" => "", count=1)
"The  foxes run quickly."

julia> replace("The quick foxes run quickly.", r"fox(es)?" => s"bus\\1")
"The quick buses run quickly."
```
"""
replace(s::AbstractString, pat_f::Pair; count=typemax(Int)) =
    replace(String(s), pat_f, count=count)

# TODO: allow transform as the first argument to replace?

# hex <-> bytes conversion

"""
    hex2bytes(s::Union{AbstractString,AbstractVector{UInt8}})

Given a string or array `s` of ASCII codes for a sequence of hexadecimal digits, returns a
`Vector{UInt8}` of bytes  corresponding to the binary representation: each successive pair
of hexadecimal digits in `s` gives the value of one byte in the return vector.

The length of `s` must be even, and the returned array has half of the length of `s`.
See also [`hex2bytes!`](@ref) for an in-place version, and [`bytes2hex`](@ref) for the inverse.

# Examples
```jldoctest
julia> s = string(12345, base = 16)
"3039"

julia> hex2bytes(s)
2-element Vector{UInt8}:
 0x30
 0x39

julia> a = b"01abEF"
6-element Base.CodeUnits{UInt8, String}:
 0x30
 0x31
 0x61
 0x62
 0x45
 0x46

julia> hex2bytes(a)
3-element Vector{UInt8}:
 0x01
 0xab
 0xef
```
"""
function hex2bytes end

hex2bytes(s::AbstractString) = hex2bytes(String(s))
hex2bytes(s::Union{String,AbstractVector{UInt8}}) = hex2bytes!(Vector{UInt8}(undef, length(s) >> 1), s)

_firstbyteidx(s::String) = 1
_firstbyteidx(s::AbstractVector{UInt8}) = first(eachindex(s))
_lastbyteidx(s::String) = sizeof(s)
_lastbyteidx(s::AbstractVector{UInt8}) = lastindex(s)

"""
    hex2bytes!(d::AbstractVector{UInt8}, s::Union{String,AbstractVector{UInt8}})

Convert an array `s` of bytes representing a hexadecimal string to its binary
representation, similar to [`hex2bytes`](@ref) except that the output is written in-place
in `d`.   The length of `s` must be exactly twice the length of `d`.
"""
function hex2bytes!(d::AbstractVector{UInt8}, s::Union{String,AbstractVector{UInt8}})
    if 2length(d) != sizeof(s)
        isodd(sizeof(s)) && throw(ArgumentError("input hex array must have even length"))
        throw(ArgumentError("output array must be half length of input array"))
    end
    j = first(eachindex(d)) - 1
    for i = _firstbyteidx(s):2:_lastbyteidx(s)
        @inbounds d[j += 1] = number_from_hex(_nthbyte(s,i)) << 4 + number_from_hex(_nthbyte(s,i+1))
    end
    return d
end

@inline number_from_hex(c) =
    (UInt8('0') <= c <= UInt8('9')) ? c - UInt8('0') :
    (UInt8('A') <= c <= UInt8('F')) ? c - (UInt8('A') - 0x0a) :
    (UInt8('a') <= c <= UInt8('f')) ? c - (UInt8('a') - 0x0a) :
    throw(ArgumentError("byte is not an ASCII hexadecimal digit"))

"""
    bytes2hex(a::AbstractArray{UInt8}) -> String
    bytes2hex(io::IO, a::AbstractArray{UInt8})

Convert an array `a` of bytes to its hexadecimal string representation, either
returning a `String` via `bytes2hex(a)` or writing the string to an `io` stream
via `bytes2hex(io, a)`.  The hexadecimal characters are all lowercase.

# Examples
```jldoctest
julia> a = string(12345, base = 16)
"3039"

julia> b = hex2bytes(a)
2-element Vector{UInt8}:
 0x30
 0x39

julia> bytes2hex(b)
"3039"
```
"""
function bytes2hex end

function bytes2hex(a::Union{NTuple{<:Any, UInt8}, AbstractArray{UInt8}})
    b = Base.StringVector(2*length(a))
    @inbounds for (i, x) in enumerate(a)
        b[2i - 1] = hex_chars[1 + x >> 4]
        b[2i    ] = hex_chars[1 + x & 0xf]
    end
    return String(b)
end

function bytes2hex(io::IO, a::Union{NTuple{<:Any, UInt8}, AbstractArray{UInt8}})
    for x in a
        print(io, Char(hex_chars[1 + x >> 4]), Char(hex_chars[1 + x & 0xf]))
    end
end

# check for pure ASCII-ness
function ascii(s::String)
    for i in 1:sizeof(s)
        @inbounds codeunit(s, i) < 0x80 || __throw_invalid_ascii(s, i)
    end
    return s
end
@noinline __throw_invalid_ascii(s::String, i::Int) = throw(ArgumentError("invalid ASCII at index $i in $(repr(s))"))

"""
    ascii(s::AbstractString)

Convert a string to `String` type and check that it contains only ASCII data, otherwise
throwing an `ArgumentError` indicating the position of the first non-ASCII byte.

# Examples
```jldoctest
julia> ascii("abcdeγfgh")
ERROR: ArgumentError: invalid ASCII at index 6 in "abcdeγfgh"
Stacktrace:
[...]

julia> ascii("abcdefgh")
"abcdefgh"
```
"""
ascii(x::AbstractString) = ascii(String(x))

Base.rest(s::Union{String,SubString{String}}, i=1) = SubString(s, i)
function Base.rest(s::AbstractString, st...)
    io = IOBuffer()
    for c in Iterators.rest(s, st...)
        print(io, c)
    end
    return String(take!(io))
end

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
int_rules_1_1_1_7 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.1.7*P(x)*(a+b*x)^m*(c+d*x)^n*(e+f*x)^p*(g+h*x)^q =#
    @apply_utils Antiderivative(~Px * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (Coeff(~Px, ~x, ~n - 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) + Antiderivative((~Px - Coeff(~Px, ~x, ~n - 1) * (~x) ^ (~n - 1)) * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~Px, ~x) && (IGtQ(~p, 1) && (IGtQ(~n, 1) && (NeQ(Coeff(~Px, ~x, ~n - 1), 0) && (NeQ(~Px, Coeff(~Px, ~x, ~n - 1) * (~x) ^ (~n - 1)) && Not(MatchQ(~Px, ~(Qx') * (~c + ~(d') * (~x) ^ ~m) ^ ~q <-- FreeQ([c, d], ~x) && (PolyQ(Qx, ~x) && (IGtQ(q, 1) && (IGtQ(m, 1) && (NeQ(Coeff(Qx * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x, m - 1), 0) && GtQ(m * q, ~n * ~p))))))))))))
    @apply_utils Antiderivative(~Px * (~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (Coeff(~Px, ~x, (~n - ~m) - 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) + Antiderivative((~Px - Coeff(~Px, ~x, (~n - ~m) - 1) * (~x) ^ ((~n - ~m) - 1)) * (~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (PolyQ(~Px, ~x) && (IGtQ(~p, 1) && (IGtQ(~n - ~m, 0) && NeQ(Coeff(~Px, ~x, (~n - ~m) - 1), 0))))
    @apply_utils Antiderivative(~(u') * (~x) ^ ~(m') * (~(a') * (~x) ^ ~(p') + ~(b') * (~x) ^ ~(q')) ^ ~(n'), ~x) => Antiderivative(~u * (~x) ^ (~m + ~n * ~p) * (~a + ~b * (~x) ^ (~q - ~p)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~p, ~q], ~x) && (IntegerQ(~n) && PosQ(~q - ~p))
    @apply_utils Antiderivative(~(u') * (~x) ^ ~(m') * (~(a') * (~x) ^ ~(p') + ~(b') * (~x) ^ ~(q') + ~(c') * (~x) ^ ~(r')) ^ ~(n'), ~x) => Antiderivative(~u * (~x) ^ (~m + ~n * ~p) * (~a + ~b * (~x) ^ (~q - ~p) + ~c * (~x) ^ (~r - ~p)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~p, ~q, ~r], ~x) && (IntegerQ(~n) && (PosQ(~q - ~p) && PosQ(~r - ~p)))
    @apply_utils Antiderivative(~(u') * (~Px) ^ ~(p') * (~Qx) ^ ~(q'), ~x) => Antiderivative(~u * PolynomialQuotient(~Px, ~Qx, ~x) ^ ~p * (~Qx) ^ (~p + ~q), ~x) <-- FreeQ(~q, ~x) && (PolyQ(~Px, ~x) && (PolyQ(~Qx, ~x) && (EqQ(PolynomialRemainder(~Px, ~Qx, ~x), 0) && (IntegerQ(~p) && LtQ(~p * ~q, 0)))))
    @apply_utils Antiderivative(~Pp / ~Qq, ~x) => With([p = Expon(~Pp, ~x), q = Expon(~Qq, ~x)], (Coeff(~Pp, ~x, p) * log(RemoveContent(~Qq, ~x))) / (q * Coeff(~Qq, ~x, q)) <-- EqQ(p, q - 1) && EqQ(~Pp, Simplify((Coeff(~Pp, ~x, p) / (q * Coeff(~Qq, ~x, q))) * D(~Qq, ~x)))) <-- PolyQ(~Pp, ~x) && PolyQ(~Qq, ~x)
    @apply_utils Antiderivative(~Pp * (~Qq) ^ ~(m'), ~x) => With([p = Expon(~Pp, ~x), q = Expon(~Qq, ~x)], (Coeff(~Pp, ~x, p) * (~x) ^ ((p - q) + 1) * (~Qq) ^ (~m + 1)) / ((p + ~m * q + 1) * Coeff(~Qq, ~x, q)) <-- NeQ(p + ~m * q + 1, 0) && EqQ((p + ~m * q + 1) * Coeff(~Qq, ~x, q) * ~Pp, Coeff(~Pp, ~x, p) * (~x) ^ (p - q) * (((p - q) + 1) * ~Qq + (~m + 1) * ~x * D(~Qq, ~x)))) <-- FreeQ(~m, ~x) && (PolyQ(~Pp, ~x) && (PolyQ(~Qq, ~x) && NeQ(~m, -1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a1 + ~(b1') * (~x) ^ ~(n')) ^ ~p * (~a2 + ~(b2') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~a1 + ~b1 * (~x) ^ ~n) ^ (~p + 1) * (~a2 + ~b2 * (~x) ^ ~n) ^ (~p + 1)) / (2 * ~b1 * ~b2 * ~n * (~p + 1)) <-- FreeQ([~a1, ~b1, ~a2, ~b2, ~m, ~n, ~p], ~x) && (EqQ(~a2 * ~b1 + ~a1 * ~b2, 0) && (EqQ((~m - 2 * ~n) + 1, 0) && NeQ(~p, -1)))
    @apply_utils Antiderivative(~Pp * (~Qq) ^ ~(m') * (~Rr) ^ ~(n'), ~x) => With([p = Expon(~Pp, ~x), q = Expon(~Qq, ~x), r = Expon(~Rr, ~x)], (Coeff(~Pp, ~x, p) * (~x) ^ (((p - q) - r) + 1) * (~Qq) ^ (~m + 1) * (~Rr) ^ (~n + 1)) / ((p + ~m * q + ~n * r + 1) * Coeff(~Qq, ~x, q) * Coeff(~Rr, ~x, r)) <-- NeQ(p + ~m * q + ~n * r + 1, 0) && EqQ((p + ~m * q + ~n * r + 1) * Coeff(~Qq, ~x, q) * Coeff(~Rr, ~x, r) * ~Pp, Coeff(~Pp, ~x, p) * (~x) ^ ((p - q) - r) * ((((p - q) - r) + 1) * ~Qq * ~Rr + (~m + 1) * ~x * ~Rr * D(~Qq, ~x) + (~n + 1) * ~x * ~Qq * D(~Rr, ~x)))) <-- FreeQ([~m, ~n], ~x) && (PolyQ(~Pp, ~x) && (PolyQ(~Qq, ~x) && (PolyQ(~Rr, ~x) && (NeQ(~m, -1) && NeQ(~n, -1)))))
    @apply_utils Antiderivative(~Qr * (~(a') + ~(b') * (~Pq) ^ ~(n')) ^ ~(p'), ~x) => With([q = Expon(~Pq, ~x), r = Expon(~Qr, ~x)], (Coeff(~Qr, ~x, r) / (q * Coeff(~Pq, ~x, q))) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x, ~Pq) <-- EqQ(r, q - 1) && EqQ(Coeff(~Qr, ~x, r) * D(~Pq, ~x), q * Coeff(~Pq, ~x, q) * ~Qr)) <-- FreeQ([~a, ~b, ~n, ~p], ~x) && (PolyQ(~Pq, ~x) && PolyQ(~Qr, ~x))
    @apply_utils Antiderivative(~Qr * (~(a') + ~(b') * (~Pq) ^ ~(n') + ~(c') * (~Pq) ^ ~(n2')) ^ ~(p'), ~x) => Module([q = Expon(~Pq, ~x), r = Expon(~Qr, ~x)], (Coeff(~Qr, ~x, r) / (q * Coeff(~Pq, ~x, q))) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~Pq) <-- EqQ(r, q - 1) && EqQ(Coeff(~Qr, ~x, r) * D(~Pq, ~x), q * Coeff(~Pq, ~x, q) * ~Qr)) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && PolyQ(~Qr, ~x)))
    @apply_utils Antiderivative(~(u') * (~(a') * (~x) ^ ~(p') + ~(b') * (~x) ^ ~(q')) ^ ~(n'), ~x) => Antiderivative(~u * (~x) ^ (~n * ~p) * (~a + ~b * (~x) ^ (~q - ~p)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~p, ~q], ~x) && (IntegerQ(~n) && PosQ(~q - ~p))
    @apply_utils Antiderivative(~(u') * (~(a') * (~x) ^ ~(p') + ~(b') * (~x) ^ ~(q') + ~(c') * (~x) ^ ~(r')) ^ ~(n'), ~x) => Antiderivative(~u * (~x) ^ (~n * ~p) * (~a + ~b * (~x) ^ (~q - ~p) + ~c * (~x) ^ (~r - ~p)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~p, ~q, ~r], ~x) && (IntegerQ(~n) && (PosQ(~q - ~p) && PosQ(~r - ~p)))
    #= Antiderivative(sqrt((~a')+(~b')*(~x))*((~A')+(~B')*(~x))/(sqrt((~c')+(~d')*(~x))*sqrt((~e')+(~f')*(~x) )*sqrt((~g')+(~h')*(~x))),~x) :=  B*sqrt(a+b*x)*sqrt(e+f*x)*sqrt(g+h*x)/(f*h*sqrt(c+d*x)) - B*(b*g-a*h)/(2*f*h)*Antiderivative(sqrt(e+f*x)/(sqrt(a+b*x)*sqrt(c+d*x)*sqrt(g+ h*x)),x) + B*(d*e-c*f)*(d*g-c*h)/(2*d*f*h)*Antiderivative(sqrt(a+b*x)/((c+d*x)^(3/2)*sqrt( e+f*x)*sqrt(g+h*x)),x) <-- FreeQ([a,b,c,d,e,f,g,h,A,B],x) &&  EqQ(2*A*d*f-B*(d*e+c*f),0) =#
    @apply_utils Antiderivative((sqrt(~(a') + ~(b') * ~x) * (~(A') + ~(B') * ~x)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((~b * ~B * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~d * ~f * ~h * sqrt(~a + ~b * ~x)) - ((~B * (~b * ~g - ~a * ~h)) / (2 * ~f * ~h)) * Antiderivative(sqrt(~e + ~f * ~x) / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~g + ~h * ~x)), ~x)) + ((~B * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) / (2 * ~d * ~f * ~h)) * Antiderivative(sqrt(~c + ~d * ~x) / ((~a + ~b * ~x) ^ (3 / 2) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B], ~x) && EqQ(2 * ~A * ~d * ~f - ~B * (~d * ~e + ~c * ~f), 0)
    #= Antiderivative(sqrt((~a')+(~b')*(~x))*((~A')+(~B')*(~x))/(sqrt((~c')+(~d')*(~x))*sqrt((~e')+(~f')*(~x) )*sqrt((~g')+(~h')*(~x))),~x) :=  (2*A*d*f-B*(d*e+c*f))/(2*d*f)*Antiderivative(sqrt(a+b*x)/(sqrt(c+d*x)*sqrt(e+f*x) *sqrt(g+h*x)),x) + B/(2*d*f)*Antiderivative((sqrt(a+b*x)*(d*e+c*f+2*d*f*x))/(sqrt(c+d*x)*sqrt(e+f* x)*sqrt(g+h*x)),x) <-- FreeQ([a,b,c,d,e,f,g,h,A,B],x) &&  NeQ(2*A*d*f-B*(d*e+c*f),0) =#
    @apply_utils Antiderivative((sqrt(~(a') + ~(b') * ~x) * (~(A') + ~(B') * ~x)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (((~B * sqrt(~a + ~b * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~f * ~h * sqrt(~c + ~d * ~x)) + ((~B * (~d * ~e - ~c * ~f) * (~d * ~g - ~c * ~h)) / (2 * ~d * ~f * ~h)) * Antiderivative(sqrt(~a + ~b * ~x) / ((~c + ~d * ~x) ^ (3 / 2) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x)) - ((~B * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) / (2 * ~b * ~f * ~h)) * Antiderivative(1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x)) + ((2 * ~A * ~b * ~d * ~f * ~h + ~B * (~a * ~d * ~f * ~h - ~b * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h))) / (2 * ~b * ~d * ~f * ~h)) * Antiderivative(sqrt(~a + ~b * ~x) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B], ~x) && NeQ(2 * ~A * ~d * ~f - ~B * (~d * ~e + ~c * ~f), 0)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~(m') * (~(A') + ~(B') * ~x)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (1 / (~d * ~f * ~h * (2 * ~m + 3))) * Antiderivative(((~a + ~b * ~x) ^ (~m - 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(~a * ~A * ~d * ~f * ~h * (2 * ~m + 3) + (~A * ~b + ~a * ~B) * ~d * ~f * ~h * (2 * ~m + 3) * ~x + ~b * ~B * ~d * ~f * ~h * (2 * ~m + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B], ~x) && (IntegerQ(2 * ~m) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * ~x) / (sqrt(~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((~A * ~b - ~a * ~B) / ~b) * Antiderivative(1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + (~B / ~b) * Antiderivative(sqrt(~a + ~b * ~x) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B], ~x)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * (~(A') + ~(B') * ~x)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((~A * (~b) ^ 2 - ~a * ~b * ~B) * (~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) - (1 / (2 * (~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(((~A * ((2 * (~a) ^ 2 * ~d * ~f * ~h * (~m + 1) - 2 * ~a * ~b * (~m + 1) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) + (~b) ^ 2 * (2 * ~m + 3) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h)) - ~b * ~B * (~a * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) + 2 * ~b * ~c * ~e * ~g * (~m + 1))) - 2 * ((~A * ~b - ~a * ~B) * (~a * ~d * ~f * ~h * (~m + 1) - ~b * (~m + 2) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h))) * ~x) + ~d * ~f * ~h * (2 * ~m + 5) * (~A * (~b) ^ 2 - ~a * ~b * ~B) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B], ~x) && (IntegerQ(2 * ~m) && LtQ(~m, -1))
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~(m') * (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (2 * ~C * (~a + ~b * ~x) ^ ~m * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~d * ~f * ~h * (2 * ~m + 3)) + (1 / (~d * ~f * ~h * (2 * ~m + 3))) * Antiderivative(((~a + ~b * ~x) ^ (~m - 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((~a * ~A * ~d * ~f * ~h * (2 * ~m + 3) - ~C * (~a * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) + 2 * ~b * ~c * ~e * ~g * ~m)) + ((~A * ~b + ~a * ~B) * ~d * ~f * ~h * (2 * ~m + 3) - ~C * (2 * ~a * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h) + ~b * (2 * ~m + 1) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h))) * ~x + (~b * ~B * ~d * ~f * ~h * (2 * ~m + 3) + 2 * ~C * (~a * ~d * ~f * ~h * ~m - ~b * (~m + 1) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h))) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B, ~C], ~x) && (IntegerQ(2 * ~m) && GtQ(~m, 0))
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~(m') * (~(A') + ~(C') * (~x) ^ 2)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (2 * ~C * (~a + ~b * ~x) ^ ~m * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~d * ~f * ~h * (2 * ~m + 3)) + (1 / (~d * ~f * ~h * (2 * ~m + 3))) * Antiderivative(((~a + ~b * ~x) ^ (~m - 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((~a * ~A * ~d * ~f * ~h * (2 * ~m + 3) - ~C * (~a * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) + 2 * ~b * ~c * ~e * ~g * ~m)) + (~A * ~b * ~d * ~f * ~h * (2 * ~m + 3) - ~C * (2 * ~a * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h) + ~b * (2 * ~m + 1) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h))) * ~x + 2 * ~C * (~a * ~d * ~f * ~h * ~m - ~b * (~m + 1) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~C], ~x) && (IntegerQ(2 * ~m) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) / (sqrt(~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (~C * sqrt(~a + ~b * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~b * ~f * ~h * sqrt(~c + ~d * ~x)) + ((~C * (~d * ~e - ~c * ~f) * (~d * ~g - ~c * ~h)) / (2 * ~b * ~d * ~f * ~h)) * Antiderivative(sqrt(~a + ~b * ~x) / ((~c + ~d * ~x) ^ (3 / 2) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + (1 / (2 * ~b * ~d * ~f * ~h)) * Antiderivative((1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((2 * ~A * ~b * ~d * ~f * ~h - ~C * (~b * ~d * ~e * ~g + ~a * ~c * ~f * ~h)) + (2 * ~b * ~B * ~d * ~f * ~h - ~C * (~a * ~d * ~f * ~h + ~b * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h))) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B, ~C], ~x)
    @apply_utils Antiderivative((~(A') + ~(C') * (~x) ^ 2) / (sqrt(~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (~C * sqrt(~a + ~b * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~b * ~f * ~h * sqrt(~c + ~d * ~x)) + ((~C * (~d * ~e - ~c * ~f) * (~d * ~g - ~c * ~h)) / (2 * ~b * ~d * ~f * ~h)) * Antiderivative(sqrt(~a + ~b * ~x) / ((~c + ~d * ~x) ^ (3 / 2) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + (1 / (2 * ~b * ~d * ~f * ~h)) * Antiderivative((1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((2 * ~A * ~b * ~d * ~f * ~h - ~C * (~b * ~d * ~e * ~g + ~a * ~c * ~f * ~h)) - ~C * (~a * ~d * ~f * ~h + ~b * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~C], ~x)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C) * (~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) - (1 / (2 * (~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(((~A * ((2 * (~a) ^ 2 * ~d * ~f * ~h * (~m + 1) - 2 * ~a * ~b * (~m + 1) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) + (~b) ^ 2 * (2 * ~m + 3) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h)) - (~b * ~B - ~a * ~C) * (~a * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) + 2 * ~b * ~c * ~e * ~g * (~m + 1))) - 2 * ((~A * ~b - ~a * ~B) * (~a * ~d * ~f * ~h * (~m + 1) - ~b * (~m + 2) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) - ~C * (((~a) ^ 2 * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h) - (~b) ^ 2 * ~c * ~e * ~g * (~m + 1)) + ~a * ~b * (~m + 1) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h))) * ~x) + ~d * ~f * ~h * (2 * ~m + 5) * ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B, ~C], ~x) && (IntegerQ(2 * ~m) && LtQ(~m, -1))
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * (~(A') + ~(C') * (~x) ^ 2)) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((~A * (~b) ^ 2 + (~a) ^ 2 * ~C) * (~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) - (1 / (2 * (~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(((~A * ((2 * (~a) ^ 2 * ~d * ~f * ~h * (~m + 1) - 2 * ~a * ~b * (~m + 1) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) + (~b) ^ 2 * (2 * ~m + 3) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h)) + ~a * ~C * (~a * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) + 2 * ~b * ~c * ~e * ~g * (~m + 1))) - 2 * (~A * ~b * (~a * ~d * ~f * ~h * (~m + 1) - ~b * (~m + 2) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) - ~C * (((~a) ^ 2 * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h) - (~b) ^ 2 * ~c * ~e * ~g * (~m + 1)) + ~a * ~b * (~m + 1) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h))) * ~x) + ~d * ~f * ~h * (2 * ~m + 5) * (~A * (~b) ^ 2 + (~a) ^ 2 * ~C) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~C], ~x) && (IntegerQ(2 * ~m) && LtQ(~m, -1))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p') * (~(g') + ~(h') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p, ~q], ~x) && (PolyQ(~Px, ~x) && IntegersQ(~m, ~n))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p') * (~(g') + ~(h') * ~x) ^ ~(q'), ~x) => PolynomialRemainder(~Px, ~a + ~b * ~x, ~x) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x) + Antiderivative(PolynomialQuotient(~Px, ~a + ~b * ~x, ~x) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p, ~q], ~x) && (PolyQ(~Px, ~x) && EqQ(~m, -1))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p') * (~(g') + ~(h') * ~x) ^ ~(q'), ~x) => PolynomialRemainder(~Px, ~a + ~b * ~x, ~x) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x) + Antiderivative(PolynomialQuotient(~Px, ~a + ~b * ~x, ~x) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p, ~q], ~x) && PolyQ(~Px, ~x)
end

ext: jl

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.7 P(x) (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl


text:
module TestDBI
    using DBI
    using DBDSQLite

    # User database
    path = Pkg.dir("DBDSQLite", "test", "db", "users.sqlite3")
    run(`touch $path`)
    db = connect(SQLite3, path)

    stmt = prepare(
        db,
        "CREATE TABLE users (id INT NOT NULL, name VARCHAR(255))"
    )
    @assert executed(stmt) == 0
    execute(stmt)
    @assert executed(stmt) == 1
    finish(stmt)

    try
        stmt = prepare(
            db,
            "CREATE TABLE users (id INT NOT NULL, name VARCHAR(255))"
        )
    end
    errcode(db)
    errstring(db)

    stmt = prepare(db, "INSERT INTO users VALUES (1, 'Jeff Bezanson')")
    execute(stmt)
    finish(stmt)

    stmt = prepare(db, "INSERT INTO users VALUES (2, 'Viral Shah')")
    execute(stmt)
    finish(stmt)

    run(db, "INSERT INTO users VALUES (3, 'Stefan Karpinski')")

    stmt = prepare(db, "INSERT INTO users VALUES (?, ?)")
    execute(stmt, {4, "Jameson Nash"})
    execute(stmt, {5, "Keno Fisher"})
    finish(stmt)

    stmt = prepare(db, "SELECT * FROM users")
    execute(stmt)
    row = fetchrow(stmt)
    row = fetchrow(stmt)
    row = fetchrow(stmt)
    row = fetchrow(stmt)
    row = fetchrow(stmt)
    row = fetchrow(stmt)
    finish(stmt)

    stmt = prepare(db, "SELECT * FROM users")
    execute(stmt)
    rows = fetchall(stmt)
    finish(stmt)

    stmt = prepare(db, "SELECT * FROM users")
    execute(stmt)
    rows = fetchdf(stmt)
    finish(stmt)

    rows = select(db, "SELECT * FROM users")

    tabledata = tableinfo(db, "users")

    columndata = columninfo(db, "users", "id")
    columndata = columninfo(db, "users", "name")

    stmt = prepare(db, "DROP TABLE users")
    execute(stmt)
    finish(stmt)

    disconnect(db)

    rm(path)

    # China OK database
    path = Pkg.dir("DBDSQLite", "test", "db", "chinook.sqlite3")
    db = connect(SQLite3, path)

    stmt = prepare(db, "SELECT * FROM Employee")
    execute(stmt)
    df = fetchdf(stmt)
    finish(stmt)

    df = select(
        db,
        "SELECT * FROM sqlite_master WHERE type = 'table' ORDER BY name"
    )
    df = select(db, "SELECT * FROM Album")
    df = select(
        db,
        """
        SELECT a.*, b.AlbumId
        FROM Artist a
        LEFT OUTER JOIN Album b ON b.ArtistId = a.ArtistId
        ORDER BY name
        """
    )

    disconnect(db)
end


text:
using Church

dp(concentration::Real, base_measure::Function) = begin
    sticks = Mem(i::Int -> beta(1., concentration))
    atoms  = Mem(i::Int -> base_measure())
    loop(i::Int) =
        @If(bernoulli(sticks[i]), atoms[i], loop(i+1))
    d = () -> loop(1)
end

dp_mixture(concentration::Real, base_measure::Function, parameter::Function) = begin
    dp_ = dp(concentration, base_measure)
    () -> parameter(dp_())
end

d = dp(1., normal)
ds = [d() for i = 1:20]

f = dp_mixture(1., normal, x -> normal(x, 0.1))
fs = [f() for i = 1:20]

n_data = 10
n_components = 20
indicies = [categorical(n_components) for i = 1:n_data]
components = Mem((i::Int) -> normal())
data = [components[indicies[i]] for i = 1:n_data]


text:
module Core

const API_SIGN = "API\0"
const HEADTYPE = UInt32    # sizeof(HEADTYPE) == 4 bytes
const MAX_LEN =  0xffffff

isascii(m) = all(x -> x < 0x80, m) # ASCII

function write_one(socket, buf, api_sign=false)

  msg = take!(buf)

  len = length(msg)

  @assert isascii(msg)
  @assert len ≤ MAX_LEN

  api_sign ? write(socket, API_SIGN, hton(HEADTYPE(len)), msg) :
             write(socket,           hton(HEADTYPE(len)), msg)
end

function read_one(socket)

  len = ntoh(read(socket, HEADTYPE))

  @assert len ≤ MAX_LEN

  read(socket, len)
end

end


text:
mutable struct DoubleMacaulayMatrix{C}
    # pivot -> row
    leftpivs::Vector{Vector{Int}}
    leftcoeffs::Vector{Vector{C}}

    # idx -> row
    rightrows::Vector{Vector{Int}}
    rightcoeffs::Vector{Vector{C}}

    # pivot -> idx
    pivot2idx::Vector{Int}

    #
    lefthash2col::Dict{Int, Int}
    leftcol2hash::Vector{Int}

    righthash2col::Dict{Int, Int}
    rightcol2hash::Vector{Int}

    # size of leftpivs
    nlsize::Int

    # load of rightrows
    nrrows::Int
    # size of rightrows
    nrsize::Int

    # left cols filled
    nlcols::Int
    # right cols filled
    nrcols::Int
end

function initialize_double_matrix(basis::Basis{C}) where {C<:Coeff}
    n = length(basis.gens)
    leftrows   = Vector{Vector{Int}}(undef, n)
    leftcoeffs = Vector{Vector{C}}(undef, n)

    pivot2idx = Vector{Int}(undef, n)

    m = length(basis.gens)
    rightrows  = Vector{Vector{Int}}(undef, m)
    rightcoeffs= Vector{Vector{C}}(undef, m)

    lsz = n
    lefthash2col = Dict{Int, Int}()
    leftcol2hash = Vector{Int}(undef, lsz)

    rsz = m
    righthash2col = Dict{Int, Int}()
    rightcol2hash = Vector{Int}(undef, rsz)

    DoubleMacaulayMatrix(leftrows, leftcoeffs, rightrows,
                rightcoeffs, pivot2idx, lefthash2col, leftcol2hash,
                righthash2col, rightcol2hash, n, 0, m, 0, 0)
end

function convert_to_double_dense_row(matrix, monom, vector::Basis{C}, ht) where {C<:Coeff}

    if matrix.nrcols >= length(matrix.rightcol2hash)
        resize!(matrix.rightcol2hash, 2*length(matrix.rightcol2hash))
    end
    matrix.nrcols += 1
    matrix.rightcol2hash[matrix.nrcols] = monom
    matrix.righthash2col[monom] = matrix.nrcols

    rightrow = zeros(C, matrix.nrcols)
    rightrow[end] = one(C)

    exps, coeffs = vector.gens[1], vector.coeffs[1]
    for i in 1:length(exps)
        if !haskey(matrix.lefthash2col, exps[i])
            if matrix.nlcols >= length(matrix.leftcol2hash)
                resize!(matrix.leftcol2hash, 2*length(matrix.leftcol2hash))
            end
            matrix.nlcols += 1
            matrix.lefthash2col[exps[i]] = matrix.nlcols
            matrix.leftcol2hash[matrix.nlcols] = exps[i]
        end
    end

    leftrow = zeros(C, matrix.nlcols)
    for i in 1:length(exps)
        leftrow[matrix.lefthash2col[exps[i]]] = coeffs[i]
    end

    leftrow, rightrow
end

# reduces row by mul*cfs modulo ch at indices positions
#
# Finite field magic specialization
function reduce_by_pivot_simultaneous!(leftrow, leftexps, leftcfs::Vector{CoeffFF},
                rightrow, rightexps, rightcfs, magic)

    # mul = -densecoeffs[i]
    # actually.. not bad!
    mul = (magic.divisor - leftrow[leftexps[1]]) % magic

    @inbounds for j in 1:length(leftexps)
        idx = leftexps[j]
        leftrow[idx] = (leftrow[idx] + mul*leftcfs[j]) % magic
    end

    @inbounds for j in 1:length(rightexps)
        idx = rightexps[j]
        rightrow[idx] = (rightrow[idx] + mul*rightcfs[j]) % magic
    end

    mul
end

#
# Finite field magic specialization
function normalize_double_row_sparse!(leftcfs::Vector{CoeffFF}, rightcfs, magic)
    pinv = invmod(leftcfs[1], magic.divisor) % magic
    @inbounds for i in 2:length(leftcfs)
        # row[i] *= pinv
        leftcfs[i] = (leftcfs[i] * pinv) % magic
    end
    @inbounds leftcfs[1] = one(leftcfs[1])

    @inbounds for i in 1:length(rightcfs)
        # row[i] *= pinv
        rightcfs[i] = (rightcfs[i] * pinv) % magic
    end
end

#
# Finite field magic specialization
function normalize_double_row_sparse!(leftcfs::Vector{CoeffQQ}, rightcfs, magic)
    pinv = inv(leftcfs[1])
    @inbounds for i in 2:length(leftcfs)
        # row[i] *= pinv
        leftcfs[i] = leftcfs[i] * pinv
    end
    @inbounds leftcfs[1] = one(leftcfs[1])

    @inbounds for i in 1:length(rightcfs)
        rightcfs[i] = rightcfs[i] * pinv
    end
end

# reduces row by mul*cfs modulo ch at indices positions
#
# Rational field specialization
function reduce_by_pivot_simultaneous!(leftrow, leftexps, leftcfs::Vector{CoeffQQ},
                rightrow, rightexps, rightcfs, magic)

    # mul = -densecoeffs[i]
    # actually.. not bad!

    mul = -leftrow[leftexps[1]]

    @inbounds for j in 1:length(leftexps)
        idx = leftexps[j]
        leftrow[idx] = leftrow[idx] + mul*leftcfs[j]
    end

    @inbounds for j in 1:length(rightexps)
        idx = rightexps[j]
        rightrow[idx] = rightrow[idx] + mul*rightcfs[j]
    end

    mul
end

function reduce_double_dense_row_by_known_pivots_sparse!(
            matrix::DoubleMacaulayMatrix{C},
            leftrow, rightrow, magic) where {C}

    leftrows  = matrix.leftpivs
    rightrows = matrix.rightrows

    pivot2idx = matrix.pivot2idx

    # new row nonzero elements count
    k = 0
    uzero = C(0)

    # new pivot index
    np = -1

    if debug()
        @warn "in reduce" matrix.nlcols matrix.nrcols matrix.leftpivs
        @warn "hmm" leftrow
    end

    for i in 1:matrix.nlcols

        # if row element zero - no reduction
        @inbounds if leftrow[i] == uzero
            continue
        end

        # TODO: check this first?
        if !isassigned(leftrows, i)
            if np == -1
                np = i
            end
            k += 1
            continue
        end

        # exponents of reducer row at column i
        leftexps = leftrows[i]
        leftcfs  = matrix.leftcoeffs[i]

        # here map pivot --> when added
        rightexps = rightrows[pivot2idx[i]]
        rightcfs  = matrix.rightcoeffs[pivot2idx[i]]

        mul = reduce_by_pivot_simultaneous!(leftrow, leftexps, leftcfs,
                        rightrow, rightexps, rightcfs, magic)

    end

    return k == 0, np, k
end

function extract_sparse_row(row)
    newrow, newcfs, k = extract_sparse_row(row, 1, length(row))
    resize!(newrow, k)
    resize!(newcfs, k)
    newrow, newcfs, k
end

function extract_sparse_row(row::Vector{C}, np, k) where {C}
    newrow = Vector{Int}(undef, k)
    newcfs = Vector{C}(undef, k)

    # store new row in sparse format
    # where k - number of structural nonzeros in new reduced row, k > 0
    j = 1
    @inbounds for i in np:length(row) # from new pivot
        @inbounds if row[i] != 0
            newrow[j] = i
            newcfs[j] = row[i]
            j += 1
        end
    end

    newrow, newcfs, j - 1
end

function linear_relation!(
            matrix::DoubleMacaulayMatrix,
            monom::Int, vector::Basis{C},
            ht) where {C<:Coeff}

    magic = select_divisor(vector.coeffs, vector.ch)

    leftrow, rightrow = convert_to_double_dense_row(matrix, monom, vector, ht)

    if debug()
        @warn "start"
        println(monom)
        println(vector.gens, " ", vector.coeffs)
        println(leftrow)
        println(rightrow)
        println(matrix)
    end

    reduced, np, k = reduce_double_dense_row_by_known_pivots_sparse!(matrix, leftrow, rightrow, magic)

    if debug()
        @warn "reduced"
        println(reduced, " ", np, " ", k)
        println(leftrow)
        println(rightrow)
    end

    if reduced
        # pass
    else
        lexps, lcoeffs, _ = extract_sparse_row(leftrow, np, k)
        rexps, rcoeffs, _ = extract_sparse_row(rightrow)

        normalize_double_row_sparse!(lcoeffs, rcoeffs, magic)

        if debug()
            @warn "extracted"
            println(lexps, " ", lcoeffs)
            println(rexps, " ", rcoeffs)
        end

        while np >= matrix.nlsize
            matrix.nlsize *= 2
            resize!(matrix.leftpivs, matrix.nlsize)
            resize!(matrix.leftcoeffs, matrix.nlsize)
            resize!(matrix.pivot2idx, matrix.nlsize)
        end
        matrix.leftpivs[np] = lexps
        matrix.leftcoeffs[np] = lcoeffs

        matrix.nrrows += 1
        matrix.pivot2idx[np] = matrix.nrrows

        if matrix.nrrows >= matrix.nrsize
            matrix.nrsize *= 2
            resize!(matrix.rightrows, matrix.nrsize)
            resize!(matrix.rightcoeffs, matrix.nrsize)
        end
        matrix.rightrows[matrix.nrrows] = rexps
        matrix.rightcoeffs[matrix.nrrows] = rcoeffs
    end

    return reduced, rightrow
end


text:
using WAVI, Plots
"""
Produce a plot of the melt rate in MISMIP for specified melt rate parametrization (c.f. figure 4 in Favier 2019 10.5194/gmd-12-2255-2019)
***Options***
    - PICO_nboxP_zQ    : PICO with P boxes, ambient temperature taken at Qm depth for P = 10, 8, 5, 2 and Q = 500, 700
                        (e.g. "PICO_nbox2_z700")
    - PME              : Plume model emulator (using Lazeroms 2018 algorithm for grounding line and basal slope -- doi:10.5194/tc-12-49-2018)
    - MISMIP_1r        : Melt rate according to the MISMIP+ ice 1r experiment (doi: 10.5194/tc-14-2283-2020) scaled to match the mean melt rate
    - QuadL            : Quadratic formulation of melting with local dependency on thermal driving
    - QuadNL           : Quadratic formulation of melting with non-local dependency on thermal driving
"""

melt_rate_model = "QuadL"

function Favier2019_4km_init(melt_model)
# Grid and boundary conditions
nx = 160
ny = 20
nσ = 4
x0 = 0.0
y0 = -40000.0
dx = 4000.0
dy = 4000.0
h_mask = trues(nx, ny)
u_iszero = falses(nx + 1, ny); u_iszero[1,:] .= true
v_iszero = falses(nx, ny + 1); v_iszero[:,1] .= true; v_iszero[:,end] .= true
grid = Grid(nx=nx,
            ny=ny,
            nσ=nσ,
            x0=x0,
            y0=y0,
            dx=dx,
            dy=dy,
            h_mask=h_mask,
            u_iszero=u_iszero,
            v_iszero=v_iszero)

# Bed
bed = WAVI.mismip_plus_bed # function definition

# Inputing thickness profile, reading from binary file
#fname = "examples\\Favier2019_melt_params\\data\\MISMIP_ice0_2km_SteadyThickness.bin";
#fname = joinpath(dirname(@__FILE__), "data",  "WAVI_ice0_4km_thick_interpolated.bin")
fname = joinpath(dirname(@__FILE__), "data",  "WAVI_ice0_4km_thickness.bin")
h = Array{Float64,2}(undef, nx, ny)
read!(fname, h)
h = ntoh.(h)

# make the model
initial_conditions = InitialConditions(initial_thickness=h) # set thickness
model = Model(grid=grid,
            bed_elevation=bed,
            melt_rate = melt_model,
            initial_conditions=initial_conditions,
            solver_params=SolverParams(maxiter_picard=1)
            ) ;

# update the configuration
update_state!(model)

#contour plot melt rate
m = deepcopy(model.fields.gh.basal_melt)
m[model.fields.gh.grounded_fraction .== 1.] .= NaN
msat = deepcopy(m)
msat[msat .> 50] .= 50
#msat[zb .> -300] .= NaN;
x =  model.grid.xxh[:,1]; y = model.grid.yyh[1,:];
plt = heatmap(x ./ 1e3,y / 1e3, msat',
                fill=true,
                linewidth=0,
                colorbar=true,
                colorbar_title="melt rate (m/yr)",
                framestyle=:box)

xlims!((420, 640))
xlabel!("x (km)")
ylabel!("y (km)")

return model, plt
end

#Need these for PICO
nx = 160; ny = 20;
ice_front_mask = zeros(nx,ny);
ice_front_mask[end,:] .= 1;

if melt_rate_model == "PICO_nbox10_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 =1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.94e-5,
                nbox = 10);

elseif melt_rate_model == "PICO_nbox8_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 =1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.91e-5,
                nbox = 8);

elseif melt_rate_model == "PICO_nbox5_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 = 1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.87e-5,
                nbox = 5);

elseif melt_rate_model == "PICO_nbox2_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 = 1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.85e-5,
                nbox = 2);

elseif melt_rate_model == "PME"
    melt_model = PlumeEmulator(α=1.49)

elseif melt_rate_model == "MISMIP_1r"
    melt_model = MISMIPMeltRateOne(α = 0.184)

elseif melt_rate_model == "QuadL"
    melt_model = QuadraticMeltRate(γT = 0.745*1e-3)
    #melt_model = QuadraticMeltRate(γT = 0.745*1e-3, melt_partial_cell = true)

elseif melt_rate_model == "QuadNL"
    melt_model = QuadraticMeltRate(γT = 0.85*1e-3, flocal = false)
    #melt_model = QuadraticMeltRate(γT = 0.85*1e-3, flocal = false, melt_partial_cell = true)

else
    throw(ArgumentError("Specified melt rate model not found"))
end

model, plt = Favier2019_4km_init(melt_model);

m = deepcopy(model.fields.gh.basal_melt);
zb = model.fields.gh.b .* (model.fields.gh.grounded_fraction .== 1) + -  918.0 / 1028.0  .* model.fields.gh.h .* (model.fields.gh.grounded_fraction .< 1)
idx = (model.fields.gh.grounded_fraction .== 0)
mean_melt = sum(m[idx])./length(m[idx])
println("mean melt rate for shelf points is $mean_melt")


text: fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using Distributed
using Dates
import REPL
using Printf: @sprintf

if isdefined(Base, :Experimental)
  Stub = Base.Experimental
else
  Stub = Base
end

# The following is taken from https://github.com/JuliaLang/julia/blob/master/test/runtests.jl\
# Part of Julia. License is MIT: https://julialang.org/license

addprocs(n_procs)

const max_worker_rss = if haskey(ENV, "JULIA_TEST_MAXRSS_MB")
    parse(Int, ENV["JULIA_TEST_MAXRSS_MB"]) * 2^20
else
    typemax(Csize_t)
end

limited_worker_rss = max_worker_rss != typemax(Csize_t)

exit_on_error = true
n = 1
skipped = 0
seed = rand(UInt128)
running_under_rr() = false

@everywhere include("testdefs.jl")
@everywhere using Hecke
@everywhere using RandomExtensions

if short_test
  @everywhere short_test = true
else
  @everywhere short_test = false
end

if long_test
  @everywhere long_test = true
else
  @everywhere long_test = false
end

@everywhere include("setup.jl")

if with_gap
  @everywhere push!(Base.LOAD_PATH, "@v#.#")
  @everywhere using GAP
end

testgroupheader = "Test"
workerheader = "(Worker)"
name_align    = maximum([textwidth(testgroupheader) + textwidth(" ") + textwidth(workerheader); map(x -> textwidth(x) + 3 + ndigits(nworkers()), tests)])
elapsed_align = textwidth("Time (s)")
gc_align      = textwidth("GC (s)")
percent_align = textwidth("GC %")
alloc_align   = textwidth("Alloc (MB)")
rss_align     = textwidth("RSS (MB)")
printstyled(testgroupheader, color=:white)
printstyled(lpad(workerheader, name_align - textwidth(testgroupheader) + 1), " | ", color=:white)
printstyled("Time (s) | GC (s) | GC % | Alloc (MB) | RSS (MB)\n", color=:white)
results = []
print_lock = stdout isa Base.LibuvStream ? stdout.lock : ReentrantLock()
if stderr isa Base.LibuvStream
  stderr.lock = print_lock
end

function print_testworker_stats(test, wrkr, resp)
  @nospecialize resp
  lock(print_lock)
  try
    printstyled(test, color=:white)
    printstyled(lpad("($wrkr)", name_align - textwidth(test) + 1, " "), " | ", color=:white)
    time_str = @sprintf("%7.2f",resp[2])
    printstyled(lpad(time_str, elapsed_align, " "), " | ", color=:white)
    gc_str = @sprintf("%5.2f", resp[5].total_time / 10^9)
    printstyled(lpad(gc_str, gc_align, " "), " | ", color=:white)

    # since there may be quite a few digits in the percentage,
    # the left-padding here is less to make sure everything fits
    percent_str = @sprintf("%4.1f", 100 * resp[5].total_time / (10^9 * resp[2]))
    printstyled(lpad(percent_str, percent_align, " "), " | ", color=:white)
    alloc_str = @sprintf("%5.2f", resp[3] / 2^20)
    printstyled(lpad(alloc_str, alloc_align, " "), " | ", color=:white)
    rss_str = @sprintf("%5.2f", resp[6] / 2^20)
    printstyled(lpad(rss_str, rss_align, " "), "\n", color=:white)
  finally
    unlock(print_lock)
  end
  nothing
end

global print_testworker_started = (name, wrkr)->begin
  pid = running_under_rr() ? remotecall_fetch(getpid, wrkr) : 0
  at = lpad("($wrkr)", name_align - textwidth(name) + 1, " ")
  lock(print_lock)
  try
    printstyled(name, at, " |", " "^elapsed_align,
                "started at $(now())",
                (pid > 0 ? " on pid $pid" : ""),
                "\n", color=:white)
  finally
    unlock(print_lock)
  end
  nothing
end

function print_testworker_errored(name, wrkr, @nospecialize(e))
  lock(print_lock)
  try
    printstyled(name, color=:red)
    printstyled(lpad("($wrkr)", name_align - textwidth(name) + 1, " "), " |",
                " "^elapsed_align, " failed at $(now())\n", color=:red)
    if isa(e, Test.TestSetException)
      for t in e.errors_and_fails
        show(t)
        println()
      end
    elseif e !== nothing
      Base.showerror(stdout, e)
    end
    println()
  finally
    unlock(print_lock)
  end
end

all_tests = tests

local stdin_monitor
all_tasks = Task[]
try
  # Monitor stdin and kill this task on ^C
  # but don't do this on Windows, because it may deadlock in the kernel
  running_tests = Dict{String, DateTime}()
  if !Sys.iswindows() && isa(stdin, Base.TTY)
    t = current_task()
    stdin_monitor = @async begin
      term = REPL.Terminals.TTYTerminal("xterm", stdin, stdout, stderr)
      try
        REPL.Terminals.raw!(term, true)
        while true
          c = read(term, Char)
          if c == '\x3'
            Base.throwto(t, InterruptException())
            break
          elseif c == '?'
            println("Currently running: ")
            tests = sort(collect(running_tests), by=x->x[2])
            foreach(tests) do (test, date)
              println(test, " (running for ", round(now()-date, Minute), ")")
            end
          end
        end
      catch e
        isa(e, InterruptException) || rethrow()
      finally
        REPL.Terminals.raw!(term, false)
      end
    end
  end
  @Stub.sync begin
    for p in workers()
      @async begin
        push!(all_tasks, current_task())
        while length(tests) > 0
          test = popfirst!(tests)
          running_tests[test] = now()
          wrkr = p
          resp = try
            remotecall_fetch(runtests, wrkr, test, test_path(test); seed=seed, isolate = false)
          catch e
            isa(e, InterruptException) && return
            Any[CapturedException(e, catch_backtrace())]
          end
          delete!(running_tests, test)
          push!(results, (test, resp))
          if length(resp) == 1
            print_testworker_errored(test, wrkr, exit_on_error ? nothing : resp[1])
            if exit_on_error
              skipped = length(tests)
              empty!(tests)
            elseif n > 1
              # the worker encountered some failure, recycle it
              # so future tests get a fresh environment
              rmprocs(wrkr, waitfor=30)
              p = addprocs_with_testenv(1)[1]
              remotecall_fetch(include, p, "testdefs.jl")
              if use_revise
                Distributed.remotecall_eval(Main, p, revise_init_expr)
              end
            end
          else
            print_testworker_stats(test, wrkr, resp)
            if resp[end] > max_worker_rss
              # the worker has reached the max-rss limit, recycle it
              # so future tests start with a smaller working set
              if n > 1
                rmprocs(wrkr, waitfor=30)
                p = addprocs_with_testenv(1)[1]
                remotecall_fetch(include, p, "testdefs.jl")
                if use_revise
                  Distributed.remotecall_eval(Main, p, revise_init_expr)
                end
              else # single process testing
                error("Halting tests. Memory limit reached : $resp > $max_worker_rss")
              end
            end
          end
        end
        if p != 1
          # Free up memory =)
          rmprocs(p, waitfor=30)
        end
      end
    end
  end
catch e
  isa(e, InterruptException) || rethrow()
  # If the test suite was merely interrupted, still print the
  # summary, which can be useful to diagnose what's going on
  foreach(task -> begin
            istaskstarted(task) || return
            istaskdone(task) && return
            try
              schedule(task, InterruptException(); error=true)
            catch ex
              @error "InterruptException" exception=ex,catch_backtrace()
            end
          end, all_tasks)
  foreach(wait, all_tasks)
finally
  if @isdefined stdin_monitor
    schedule(stdin_monitor, InterruptException(); error=true)
  end
end

#=
`   Construct a testset on the master node which will hold results from all the
test files run on workers and on node1. The loop goes through the results,
inserting them as children of the overall testset if they are testsets,
handling errors otherwise.
Since the workers don't return information about passing/broken tests, only
errors or failures, those Result types get passed `nothing` for their test
expressions (and expected/received result in the case of Broken).
If a test failed, returning a `RemoteException`, the error is displayed and
the overall testset has a child testset inserted, with the (empty) Passes
and Brokens from the worker and the full information about all errors and
failures encountered running the tests. This information will be displayed
as a summary at the end of the test run.
If a test failed, returning an `Exception` that is not a `RemoteException`,
it is likely the julia process running the test has encountered some kind
of internal error, such as a segfault.  The entire testset is marked as
Errored, and execution continues until the summary at the end of the test
run, where the test file is printed out as the "failed expression".
=#
Test.TESTSET_PRINT_ENABLE[] = false
o_ts = Test.DefaultTestSet("Overall")
Test.push_testset(o_ts)
completed_tests = Set{String}()
for (testname, (resp,)) in results
  push!(completed_tests, testname)
  if isa(resp, Test.DefaultTestSet)
    Test.push_testset(resp)
    Test.record(o_ts, resp)
    Test.pop_testset()
  elseif isa(resp, Test.TestSetException)
    fake = Test.DefaultTestSet(testname)
    for i in 1:resp.pass
      Test.record(fake, VERSION < v"1.7.0-DEV.1196" ? Test.Pass(:test, nothing, nothing, nothing) : Test.Pass(:test, nothing, nothing, nothing, LineNumberNode(@__LINE__, @__FILE__)))
    end
    for i in 1:resp.broken
      Test.record(fake, Test.Broken(:test, nothing))
    end
    for t in resp.errors_and_fails
      Test.record(fake, t)
    end
    Test.push_testset(fake)
    Test.record(o_ts, fake)
    Test.pop_testset()
  else
    if !isa(resp, Exception)
      resp = ErrorException(string("Unknown result type : ", typeof(resp)))
    end
    # If this test raised an exception that is not a remote testset exception,
    # i.e. not a RemoteException capturing a TestSetException that means
    # the test runner itself had some problem, so we may have hit a segfault,
    # deserialization errors or something similar.  Record this testset as Errored.
    fake = Test.DefaultTestSet(testname)
    Test.record(fake, Test.Error(:nontest_error, testname, nothing, Any[(resp, [])], LineNumberNode(1)))
    Test.push_testset(fake)
    Test.record(o_ts, fake)
    Test.pop_testset()
  end
end
for test in all_tests
  (test in completed_tests) && continue
  fake = Test.DefaultTestSet(test)
  Test.record(fake, Test.Error(:test_interrupted, test, nothing, [("skipped", [])], LineNumberNode(1)))
  Test.push_testset(fake)
  Test.record(o_ts, fake)
  Test.pop_testset()
end
Test.TESTSET_PRINT_ENABLE[] = true
println()
Test.print_test_results(o_ts, 1)
if !o_ts.anynonpass
  println("    \033[32;1mSUCCESS\033[0m")
else
  println("    \033[31;1mFAILURE\033[0m\n")
  skipped > 0 &&
  println("$skipped test", skipped > 1 ? "s were" : " was", " skipped due to failure.")
  println("The global RNG seed was 0x$(string(seed, base = 16)).\n")
  Test.print_test_errors(o_ts)
  throw(Test.FallbackTestSetException("Test run finished with errors"))
end


text:
#######################################################
### ------------------ DEFINITIONS -------------------
#######################################################

abstract type Dimension end
abstract type TwoD   <: Dimension end
abstract type ThreeD <: Dimension end

struct Square <: TwoD end
struct Circle <: TwoD end
struct Cube   <: ThreeD end
struct Sphere <: ThreeD end

struct Atom{T <: Dimension}
    shape::T
    r::Matrix
    N::Int64
    sizes::Any
end

get_dimension(atom::Atom{T}) where T <: TwoD = 2
get_dimension(atom::Atom{T}) where T <: ThreeD = 3

abstract type Pump end
abstract type PlaneWave <: Pump end
abstract type Gaussian  <: Pump end

struct PlaneWave2D <: PlaneWave
    direction::Vector
end
struct PlaneWave3D <: PlaneWave
    direction::Vector
end
struct Gaussian2D  <: Gaussian
    w₀::Float64
end
struct Gaussian3D  <: Gaussian
    w₀::Float64
end

struct Laser{T <: Pump}
    pump::T
    s::Float64
    Δ::Float64
end
struct Lasers{T <: Pump}
    pump::Vector{T}
    s::Vector{Float64}
    Δ::Vector{Float64}
end

abstract type Physics end
abstract type Linear    <: Physics end
abstract type NonLinear <: Physics end

struct Scalar    <: Linear end
struct Vectorial <: Linear end
struct MeanField <: NonLinear end
struct BBGKY     <: NonLinear end

struct LinearOptics{T <: Linear}
    physic::T
    atoms::Atom
    laser::Laser
    kernelFunction::Function
    spectrum::Dict
    data::Dict
end

struct NonLinearOptics{T <: NonLinear}
    physic::T
    atoms::Atom
    laser::Laser
    excitations::Dict
    data::Dict
end

#######################################################
### ------------------ CONSTRUCTORS -------------------
#######################################################

# """
#     b₀_of(atoms::Cube)

# Formula : `ρ^2*N/( (4π/3 * (3/(16π))^3))^(1/3)  )`
# """
# b₀_of(atoms::Cube) = (ρ_of(atoms) .^ 2 * atoms.N / ((4π / 3) * (3 / (16π))^3))^(1 / 3)
# """
#     ρ_of(atoms::Cube)
# Formula : `N/( kL^3 )`
# """
# ρ_of(atoms::Cube) = atoms.N / atoms.kL^3

# """
#     b₀_of(atoms::Cube)

# Formula : `ρ^2*N/( (4π/3 * (3/(16π))^3))^(1/3)  )`
# """
# b₀_of(atoms::Sphere) = (ρ_of(atoms) .^ 2 * atoms.N / ((4π / 3) * (3 / (16π))^3))^(1 / 3)
# ρ_of(atoms::Sphere) = atoms.N / VolumeSphere(atoms.kR)

# ### --------------- LASER ---------------
# """
#     PlaneWave_3D(direction=:z, s=1e-6, Δ=0)

# `direction` can be [:x, :y, :z]
# """
# function PlaneWave_3D(direction=:z, s=1e-6, Δ=0)
#     if direction == :x
#         matrix_slace = 1
#     elseif direction == :y
#         matrix_slace = 2
#     elseif direction == :z
#         matrix_slace = 3
#     else
#         @error("No support for this direction yet")
#     end
#     return PlaneWave(matrix_slace, s, Δ)
# end
# """
#     PlaneWave_2D(direction=:z, s=1e-6, Δ=0)

# `direction` can be [:x, :y]
# """
# function PlaneWave_2D(direction=:x, s=1e-6, Δ=0)
#     if direction == :x
#         matrix_slace = 1
#     elseif direction == :y
#         matrix_slace = 2
#     else
#         @error("No support for this direction yet")
#     end
#     return PlaneWave(matrix_slace, s, Δ)
# end
# """


text:
#jl
"""
Struct Material
Used for solid Materials in Pipes. So far includes just name, thermal
conductivity, and surface roughness.
"""
struct Material
    name::String
    thermal_conductivity::Function
    roughness::Real
end
struct Pipe{T<:Real}
    length::T
    diameter::T
    radius::T
    flowarea::T
    thickness::T
    outer_diameter::T
    outer_radius::T
    totalarea::T
    roughness::T
    elevation_change::T
    therm_cond::Function
end

const valid_Pipe_diams = [0.125, 0.25, 0.375, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5,
                    3.0, 3.5, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 16.0]

function Pipe(l, d, t, ϵ, z, kfunc)
    # Assumes a circular Pipe of diameter d
    return Pipe(l,d,d/2,pi*(d/2)^2,t,d+2t,(d+2t)/2,pi*((d+2t)/2)^2,ϵ,z,kfunc)
end

mutable struct PipeProps{T<:Real}
    dᵢ::T
    dₒ::T
    thickness::T
    rᵢ::T
    rₒ::T
    units::String
end

function PipeProps(d, t, units="in")
  di = d - 2t
  ro = d/2
  ri = di/2
  return PipeProps(di, d, t, ri, ro, units)
end

const conversions = Dict(["cm","m"] => 1e-2,
                   ["m","cm"] => 100,
                   ["in","cm"] => 2.54,
                   ["cm","in"] => 1/2.54,
                   ["in","m"] => 2.54e-2,
                   ["m","in"] => 1/2.54e-2,
                   ["m","m"] => 1.0,
                   ["cm","cm"] => 1.0,
                   ["in","in"] => 1.0)

function convert(p::PipeProps, newUnits::String)
    key_ = [p.units, newUnits]
    conversion_factor = conversions[key_]
    p.dᵢ *= conversion_factor
    p.dₒ *= conversion_factor
    p.rᵢ *= conversion_factor
    p.rₒ *= conversion_factor
    p.units = newUnits
end

const schedule40_props = [PipeProps(0.405, 0.068), PipeProps(0.540, 0.088),
                    PipeProps(0.675, 0.091), PipeProps(0.840, 0.109),
                    PipeProps(1.050, 0.113), PipeProps(1.315, 0.133),
                    PipeProps(1.660, 0.140), PipeProps(1.900, 0.145),
                    PipeProps(2.375, 0.154), PipeProps(2.875, 0.203),
                    PipeProps(3.500, 0.216), PipeProps(4.000, 0.226),
                    PipeProps(4.500, 0.237), PipeProps(5.563, 0.258),
                    PipeProps(6.625, 0.280), PipeProps(8.625, 0.322),
                    PipeProps(10.750, 0.365), PipeProps(12.750, 0.406),
                    PipeProps(16.000, 0.500)]
# Note: These outer diameters and wall thicknesses were obtained at this web
# site: http://www.engineeringtoolbox.com/steel-dimensions-d_43.html

for props in schedule40_props
    convert(props, "m")
end

const props40_by_nominal_diameter = Dict(zip(valid_Pipe_diams, schedule40_props))

const schedule_props = Dict("40" => props40_by_nominal_diameter,
                           40 => props40_by_nominal_diameter)

function Pipe(schedule_info::Tuple, mat::Material, l, Δz)
  # schedule_info is a Tuple: the schedule number comes first and the nominal
  # diameter comes second
  props = schedule_props[schedule_info[1]][schedule_info[2]]
  return Pipe(l,props.dᵢ,props.dᵢ/2,pi*(props.rᵢ)^2,props.thickness,props.dₒ,
            props.rₒ,pi*(props.rₒ)^2,mat.roughness,Δz,mat.thermal_conductivity)
end

function Pipe(schedule_info::Tuple, mat::String, l, Δz)
  return Pipe(schedule_info, Pipe_Materials[mat], l, Δz)
end

function colebrook(Re, p::Pipe, f_g)
    if f_g < 0
        return f_g * 1e5
    else
        return 10^(-1 / (2 * √(f_g))) - (p.roughness / (3.7 * p.diameter) + 2.51 / (Re * √(f_g)))
    end
end

function getf(Re, p::Pipe)
    print("Sanity check: Re = ",Re)
    function dummysolve!(f_guess, fvec)
        fvec[1] = colebrook(Re, p, f_guess[1])
    end
    return NLsolve.nlsolve(dummysolve!, [0.001]).zero[1]

end

laminarf(Re, p::Pipe) = 64.0/Re

function swameejain(Re, p::Pipe)
    return 0.25/(log10(p.roughness/(3.7*p.diameter)+ 5.74/(Re^  0.9)))^2
end

function getf2(Re, p::Pipe)
    if Re > 2000
        return swameejain(Re, p)
    else
        return laminarf(Re, p)
    end
end

function stainless_steel_thermal_conductivity(T)
    # Source: http://www.mace.manchester.ac.uk/project/research/structures/strucfire/MaterialInFire/Steel/StainlessSteel/thermalProperties.htm
    return 14.6 + 1.27e-2 *(T-273.15)
end

function copper_thermal_conductivity(T)
    # Source: http://www-ferp.ucsd.edu/LIB/PROPS/PANOS/cu.html
    return 14.6 + 1.27e-2 *(T-273.15)
end

const ss_roughness_by_finish = Dict("2D" => 1E-6,
                              "2B" => 0.5E-6,
                              "2R" => 0.2E-6,
                              "BA" => 0.2E-6,
                              "2BB" => 0.1E-6)
# These surface roughness values are in m and were obtained from the following
# corporate site: http://www.outokumpu.com/en/products-properties/
# more-stainless/stainless-steel-surface-finishes/cold-rolled-finishes/Pages/default.aspx
# Note that the given values are the maximum roughness values in the range on
# the website, except for finish 2BB, which didn't have a roughness value
# so I took the lowest value for finish 2B since the compnay says this:
# "Due to the fact that surface roughness of the finish 2BB is lower than that
# of 2B, some modifications on lubrication during forming might be needed."

const ss_mat = Material("Stainless Steel", stainless_steel_thermal_conductivity,
                  ss_roughness_by_finish["2B"])

const cu_roughness = 0.03e-3
# Copper roughness.
# Source: http://www.pressure-drop.com/Online-Calculator/rauh.html

const cu_mat = Material("Copper", copper_thermal_conductivity, cu_roughness)

const pipematerials = Dict("Copper" => cu_mat, "Stainless Steel"=>ss_mat)


text:
"""
This file contains all the functions that required to extract region information
from a given textfile repository, and from there extraction of relevant regional
data from given datasets.

There are three major types of functions contained in this .jl module:
1)  extraction of region information and attributes
2)  extraction of data from a given dataset based upon the given region
    information and list of attributes
3)  transformation of longitude coordinate systems [-180,180] <=> [0,360]

"""

# Region Information and Attributes

function regionload()
    @debug "$(Dates.now()) - Loading information on possible regions ..."
    return readdlm(joinpath(@__DIR__,"regions.txt"),',',comments=true);
end

function regioninfodisplay(regioninfo)
    @info "$(Dates.now()) - The following regions are offered in the ClimateEasy.jl"
    for ii = 1 : size(regioninfo,1); @info "$(Dates.now()) - $(ii)) $(regioninfo[ii,7])" end
end

# Find Regions Bounds

function regionbounds(reg::AbstractString)
    reginfo = regionload(); regions = reginfo[:,1]; regid = (regions .== reg);
    N,S,E,W = reginfo[regid,[3,5,6,4]];
    @debug "$(Dates.now()) - The bounds of the region are, in [N,S,E,W] format, [$(N),$(S),$(E),$(W)]."
    return [N,S,E,W]
end

function regionbounds(reg::AbstractString,reginfo::AbstractArray)
    regions = reginfo[:,1]; regid = (regions .== reg)[1];
    N,S,E,W = reginfo[regid,[3,5,6,4]];
    @debug "$(Dates.now()) - The bounds of the region are, in [N,S,E,W] format, [$(N),$(S),$(E),$(W)]."
    return [N,S,E,W]
end

function regionbounds(regID::Int64)
    reginfo = regionload();
    N,S,E,W = reginfo[regID,[3,5,6,4]];
    @debug "$(Dates.now()) - The bounds of the region are, in [N,S,E,W] format, [$(N),$(S),$(E),$(W)]."
    return [N,S,E,W]
end

function regionbounds(regID::Int64,reginfo::AbstractArray)
    N,S,E,W = reginfo[regID,[3,5,6,4]];
    @debug "$(Dates.now()) - The bounds of the region are, in [N,S,E,W] format, [$(N),$(S),$(E),$(W)]."
    return [N,S,E,W]
end

# Find Short Region Name

function regionshortname(regID::Int64)
    reginfo = regionload(); return reginfo[regID,1];
end

function regionshortname(regID::Int64,reginfo::AbstractArray)
    return reginfo[regID,1];
end

# Find Full Region Name

function regionfullname(reg::AbstractString)
    reginfo = regionload(); regions = reginfo[:,1]; regid = (regions .== reg);
    return reginfo[regid,7][1];
end

function regionfullname(reg::AbstractString,reginfo::AbstractArray)
    regions = reginfo[:,1]; regid = (regions .== reg);
    return reginfo[regid,7][1];
end

function regionfullname(regID::Int64)
    reginfo = regionload(); return reginfo[regID,7];
end

function regionfullname(regID::Int64,reginfo::AbstractArray)
    return reginfo[regID,7];
end

# Find Region Parent

function regionparent(reg::AbstractString)
    reginfo = regionload(); regions = reginfo[:,1]; regid = (regions .== reg);
    return reginfo[regid,2][1];
end

function regionparent(reg::AbstractString,reginfo::AbstractArray)
    regions = reginfo[:,1]; regid = (regions .== reg);
    return reginfo[regid,2][1];
end

function regionparent(regID::Int64)
    reginfo = regionload(); return reginfo[regID,2];
end

function regionparent(regID::Int64,reginfo::AbstractArray)
    return reginfo[regID,2];
end

# Find if the Region is Global

function regionisglobe(reg::AbstractString)
    reginfo = regionload(); regions = reginfo[:,1]; regid = (regions .== reg);
    if regid == 1; return true; else; return false end
end

function regionisglobe(reg::AbstractString,reginfo::AbstractArray)
    regions = reginfo[:,1]; regid = (regions .== reg);
    if regid == 1; return true; else; return false end
end

function regionisglobe(regID::Int64)
    if regID == 1; return true; else; return false end
end

# Find if Point / Grid is in specified Regions

function ispointinregion(plon::AbstractFloat,plat::AbstractFloat,lon::Array,lat::Array)

    minlon = minimum(lon); maxlon = maximum(lon);
    minlat = minimum(lat); maxlat = maximum(lat);

    if     plon > maxlon; plon = plon - 360;
    elseif plon < minlon; plon = plon + 360;
    end

    if (plon > minlon) && (plon < maxlon) &&
         (plat > minlat) && (plat < maxlat)
          return true
    else; return false
    end

end

function ispointinregion(plon::AbstractFloat,plat::AbstractFloat,reg)

    N,S,E,W = regionbounds(reg); lon = [E,W]; lat = [N,S];

    minlon = minimum(lon); maxlon = maximum(lon);
    minlat = minimum(lat); maxlat = maximum(lat);

    if     plon > maxlon; plon = plon - 360;
    elseif plon < minlon; plon = plon + 360;
    end

    if (plon > minlon) && (plon < maxlon) &&
         (plat > minlat) && (plat < maxlat)
          return true
    else; return false
    end

end

function ispointinregion(pcoord::AbstractArray,lon::Array,lat::Array)

    plon,plat = pcoord;
    minlon = minimum(lon); maxlon = maximum(lon);
    minlat = minimum(lat); maxlat = maximum(lat);

    if     plon > maxlon; plon = plon - 360;
    elseif plon < minlon; plon = plon + 360;
    end

    if (plon > minlon) && (plon < maxlon) &&
         (plat > minlat) && (plat < maxlat)
          return true
    else; return false
    end

end

function ispointinregion(pcoord::Array{Any,2},reg)

    plon,plat = pcoord; N,S,E,W = regionbounds(reg); lon = [E,W]; lat = [N,S];

    minlon = minimum(lon); maxlon = maximum(lon);
    minlat = minimum(lat); maxlat = maximum(lat);

    if     plon > maxlon; plon = plon - 360;
    elseif plon < minlon; plon = plon + 360;
    end

    if (plon > minlon) && (plon < maxlon) &&
         (plat > minlat) && (plat < maxlat)
          return true
    else; return false
    end

end

function isgridinregion(bounds::Array,reg)

    N,S,E,W = bounds;
    rN,rS,rE,rW = regionbounds(reg); lon = [rE,rW]; lat = [rN,rS];

    minlon = minimum(lon); maxlon = maximum(lon);
    minlat = minimum(lat); maxlat = maximum(lat);

    if     E > maxlon; E = from0360to180(E);
    elseif E < minlon; E = from180to0360(E);
    end
    if     W > maxlon; W = from0360to180(W);
    elseif W < minlon; W = from180to0360(W);
    end

    if all([W,E] > minlon) && all([W,E] < maxlon) &&
        all([N,S] > minlat) && all([N,S] < maxlat)
          return true
    else; return false
    end

end

function isgridinregion(bounds::Array{Any,2},lon::Array,lat::Array)

    N,S,E,W = bounds;

    minlon = minimum(lon); maxlon = maximum(lon);
    if     E > maxlon; E = from0360to180(E);
    elseif E < minlon; E = from180to0360(E);
    end
    if     W > maxlon; W = from0360to180(W);
    elseif W < minlon; W = from180to0360(W);
    end

    if all([W,E] > minlon) && all([W,E] < maxlon) &&
        all([N,S] > minlat) && all([N,S] < maxlat)
          return true
    else; return false
    end

end

# Find Index of given position in Region
# Assumes that we points/grid are definitely in the region.

function regionpoint(plon::AbstractFloat,plat::AbstractFloat,lon::Array,lat::Array)

    minlon = minimum(lon); maxlon = maximum(lon);
    if     plon > maxlon; plon = from0360to180(plon);
    elseif plon < minlon; plon = from180to0360(plon);
    end

    @info "$(Dates.now()) - Finding grid points in data closest to requested location ..."
    ilon = argmin(abs.(lon.-plon)); ilat = argmin(abs.(lat.-plat));
    return [ilon,ilat]

end

function regionpoint(pcoord::Array{Any,2},lon::Array,lat::Array)

    plon,plat = pcoord; minlon = minimum(lon); maxlon = maximum(lon);
    if     plon > maxlon; plon = from0360to180(plon);
    elseif plon < minlon; plon = from180to0360(plon);
    end

    @info "$(Dates.now()) - Finding grid points in data closest to requested location ..."
    ilon = argmin(abs.(lon.-plon)); ilat = argmin(abs.(lat.-plat));
    return [ilon,ilat]

end

function regiongrid(bounds::Array,lon::Array,lat::Array)

    N,S,E,W = bounds;

    minlon = minimum(lon); maxlon = maximum(lon);
    if     E > maxlon; E = from0360to180(E);
    elseif E < minlon; E = from180to0360(E);
    end
    if     W > maxlon; W = from0360to180(W);
    elseif W < minlon; W = from180to0360(W);
    end

    @info "$(Dates.now()) - Finding indices of data matching given boundaries ..."
    iE = argmin(abs.(lon.-E)); iW = argmin(abs.(lon.-W));
    iN = argmin(abs.(lat.-N)); iS = argmin(abs.(lat.-S));
    return [iN,iS,iE,iW]

end

# Tranformation of Coordinates

function from180to0360(lon)

    @info "$(Dates.now()) - Longitude of point given in [-180,180] but data range is [0,360].  Adjusting coordinates of point to match."
    lon = lon + 360;

end

function from0360to180(lon)

    @info "$(Dates.now()) - Longitude of point given in [0,360] but data range is [-180,180].  Adjusting coordinates of point to match."
    lon = lon - 360;

end

# Data Extraction

function regionpermute(data)
    irow = size(data,1); icol = size(data,2); nd = ndims(data)
    if irow < icol
        @info "$(Dates.now()) - Number of rows smaller than number of columns, indicating that data is in (lat,lon) rather than (lon,lat) format.  Permuting to (lon,lat) formatting."
        if     nd == 2; data = transpose(data)
        elseif nd > 2;  data = permutedims(data,vcat(2,1,convert(Array,3:nd)))
        end
    end
    return data
end

function regionextract(data,coord,ndim)

    if     ndim == 2; data = data[coord[1],coord[2]];
    elseif ndim == 3; data = data[coord[1],coord[2],:];
    elseif ndim == 4; data = data[coord[1],coord[2],:,:];
    end

end

function regionextractpoint(data,plon,plat,lon::Array,lat::Array)

    icoord = regionpoint(plon,plat,lon,lat); ndim = ndims(data)

    @info "$(Dates.now()) - Extracting data from coordinates (lon=$(plon),lat=$(plat)) from global datasets ..."
    pdata = regionextract(data,icoord,ndim)

    return pdata

end

function regionextractgrid(data,reg,lon::Array,lat::Array)

    data = regionpermute(data);

    @info "$(Dates.now()) - Determining indices of longitude and latitude boundaries in parent dataset ..."
    bounds = regionbounds(reg); nlon = length(lon); ndim = ndims(data);
    igrid  = regiongrid(bounds,lon,lat);
    iN = igrid[1]; iS = igrid[2]; iE = igrid[3]; iW = igrid[4];

    @debug "$(Dates.now()) - Creating vector of latitude indices to extract ..."
    if     iN < iS; iNS = iN : iS
    elseif iS < iN; iNS = iS : iN
    end

    @debug "$(Dates.now()) - Creating vector of longitude indices to extract ..."
    if     iW < iE; iWE = iW : iE
    elseif iW > iE; iWE = 1 : (iE + nlon - iW); ilon = vcat(iW:nlon,1:(iW-1));

        @info "$(Dates.now()) - West indice larger than East indice.  Reshaping of longitude vector required."
        lon[1:(iW-1)] = lon[1:(iW-1)] .+ 360; lon = lon[ilon];

        @info "$(Dates.now()) - Reshaping data matrix to match longitude vector ..."
        if     ndim == 2; data = data[ilon,:];
        elseif ndim == 3; data = data[ilon,:,:];
        elseif ndim == 4; data = data[ilon,:,:,:];
        end
    end

    @info "$(Dates.now()) - Extracting data for the $(regionfullname(reg)) region from global datasets ..."
    rdata = regionextract(data,[iWE,iNS],ndim)
    rgrid = [lon[iWE],lat[iNS]]

    return rdata,rgrid

end


text:
import Base.+
import Base.-
import Base.*
import Base./
import Base.zero
import Base.iszero
import Base.one
import Base.inv
import Base.isnan
using Images
using ColorTypes

function zero(a::Tuple{Float64, Float64, Float64, Float64}
              )::Tuple{Float64, Float64, Float64, Float64}
    return (0.0,0.0,0.0,0.0)
end
function iszero((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
              )::Bool
    return x1==0.0 && x2==0.0 && x3==0.0 && x4==0.0
end

function one((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
             )::Tuple{Float64, Float64, Float64, Float64}
    return (1.0,0.0,0.0,0.0)
end

function +((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64},
           (y1,y2,y3,y4)::Tuple{Float64, Float64, Float64, Float64}
           )::Tuple{Float64, Float64, Float64, Float64}
    return (x1+y1,x2+y2,x3+y3,x4+y4)
end

function -((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64},
           (y1,y2,y3,y4)::Tuple{Float64, Float64, Float64, Float64}
           )::Tuple{Float64, Float64, Float64, Float64}
    return (x1-y1,x2-y2,x3-y3, x4-y4)
end
function -((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
           )::Tuple{Float64, Float64, Float64, Float64}
    return (-x1,-x2,-x3,-x4)
end

function *((x1,x2,x3, x4)::Tuple{Float64, Float64, Float64, Float64}, y::Float64
           )::Tuple{Float64, Float64, Float64, Float64}
    return (x1*y,x2*y,x3*y,x4*y)
end
function *(y::Float64, (x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
           )::Tuple{Float64, Float64, Float64, Float64}
    return (x1*y,x2*y,x3*y,x4*y)
end

# \cdot tab
function ⋅((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64},
          (y1,y2,y3,y4)::Tuple{Float64, Float64, Float64, Float64})::Float64
    return x1*y1+x2*y2+x3*y3+x4*y4
end

function norm(x::Tuple{Float64, Float64, Float64, Float64})::Float64
    return sqrt(x ⋅ x)
end
function normalize(a::Tuple{Float64, Float64, Float64, Float64}
                   )::Tuple{Float64, Float64, Float64, Float64}
    n=norm(a)
    if n==0.0
        return a
    end
    return (1.0/n)*a
end

function conj((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
             )::Tuple{Float64, Float64, Float64, Float64}
    return (x1,-x2,-x3,-x4)
end

function inv((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
             )::Tuple{Float64, Float64, Float64, Float64}
    return (1.0/(x1*x1+x2*x2+x3*x3+x4*x4))*(x1,-x2,-x3,-x4)
end

function isnan((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64}
               )::Bool
    return isnan(x1) || isnan(x2) || isnan(x3) || isnan(x4)
end

function *((x1,x2,x3,x4)::Tuple{Float64, Float64, Float64, Float64},
           (y1,y2,y3,y4)::Tuple{Float64, Float64, Float64, Float64}
           )::Tuple{Float64, Float64, Float64, Float64}
    return (x1*y1-x2*y2-x3*y3-x4*y4,
            x2*y1+x1*y2+x3*y4-x4*y3,
            x3*y1+x1*y3-x2*y4+x4*y2,
            x1*y4+x4*y1+x2*y3-x3*y2)
end

function /(a::Tuple{Float64, Float64, Float64, Float64},
           b::Tuple{Float64, Float64, Float64, Float64}
           )::Tuple{Float64, Float64, Float64, Float64}
    return a*inv(b)
end

function initPalette(;colorScheme::Int64=0,
                     colorRepetitions::Int64=1)::Tuple{Vector{RGB},Int64}
    colorstepsOneColor=256
    colorsteps=6*colorRepetitions*colorstepsOneColor
    gray = 1.0/convert(Float64,colorstepsOneColor)
    colors=Array{RGB}(UndefInitializer(),colorsteps)
    for ii in 1:colorstepsOneColor
        scaledgray=gray*ii
        red=RGB(scaledgray,0.0,0.7*scaledgray)
        green=RGB(0.0,scaledgray,0.8*scaledgray)
        blue=RGB(0.4*scaledgray,0.0,scaledgray)
        if colorScheme == 1
            color1=blue
            color2=red
            color3=green
        elseif colorScheme == 2
            color1=green
            color2=blue
            color3=red
        elseif colorScheme == 3
            color1=blue
            color2=green
            color3=red
        elseif colorScheme == 4
            color1=green
            color2=red
            color3=blue
        elseif colorScheme == 5
            color1=red
            color2=blue
            color3=green
        else
            color1=red
            color2=green
            color3=blue
        end
        for jj in 0:colorRepetitions-1
            colors[jj*6*colorstepsOneColor+ii]=color2
            colors[(jj*6+2)*colorstepsOneColor-(ii-1)]=color2
            colors[(jj*6+2)*colorstepsOneColor+ii]=color1
            colors[(jj*6+4)*colorstepsOneColor-(ii-1)]=color1
            colors[(jj*6+4)*colorstepsOneColor+ii]=color3
            colors[(jj*6+6)*colorstepsOneColor-(ii-1)]=color3
        end
    end
    return (colors,colorsteps)
end

function myimage((x,y,z,u)::Tuple{Float64, Float64, Float64, Float64},
                 radius::Float64,limit::Float64,size::Int64;
                 turnIt::Tuple{Float64, Float64, Float64, Float64}=(1.0,0.0,0.0,0.0),
                 colorScheme::Int64=0,
                 colorFactor::Int64=1,
                 colorOffset::Int64=0,
                 colorRepetitions::Int64=1)::Matrix{RGB}
    image=Matrix{RGB}(UndefInitializer(),size,size)
    step = radius*2.0/convert(Float64,size)
    (colors,colorsteps) = initPalette(colorScheme=colorScheme,colorRepetitions=colorRepetitions)
    black=RGB(0.0,0.0,0.0)
    turnItNorm=normalize(turnIt)
    xpos = x-radius
    colorLimit=div(colorsteps-colorOffset,colorFactor)
    for i in 1:size
        ypos = y-radius
        for j in 1:size
            n=1
            c=(xpos,ypos,z,u)*turnItNorm
            v=zero(c)
            w=zero(c)
            while true
                if norm(v)+norm(w)>=limit
                    image[i,j] = colors[colorOffset+n*colorFactor]
                    break
                end
                if n>colorLimit-1
                    image[i,j] = black
                    break
                end
                n += 1
                vtemp = v
                vinv = inv(v)
                if isnan(vinv)
                    vinv = zero(vinv)
                end
                winv = inv(w)
                if isnan(winv)
                    winv = zero(winv)
                end
                s=v*v+winv
                t=w*w-vinv
                v = s * s * (1.0/2273.0)  + w + c
                w = t * t * t * (1.0/3709.0) + vtemp + c
            end
            ypos += step
        end
        xpos += step
    end
    return image
end

function mydraw(fn::String,
                a::Tuple{Float64, Float64, Float64, Float64},
                radius::Float64,limit::Float64,size::Int64;
                turnIt::Tuple{Float64, Float64, Float64, Float64}=(1.0,0.0,0.0,0.0),
                colorScheme::Int64=0,
                colorFactor::Int64=1,
                colorOffset::Int64=0,
                colorRepetitions::Int64=1)
    image=myimage(a,radius,limit,size,
                  turnIt=turnIt,
                  colorScheme=colorScheme,
                  colorFactor=colorFactor,
                  colorOffset=colorOffset,
                  colorRepetitions=colorRepetitions)
    save(fn,image)
end