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

## object-oriented Regex interface ##

include("pcre.jl")

const DEFAULT_COMPILER_OPTS = PCRE.UTF | PCRE.NO_UTF_CHECK | PCRE.ALT_BSUX | PCRE.UCP
const DEFAULT_MATCH_OPTS = PCRE.NO_UTF_CHECK

"""
    An abstract type representing any sort of pattern matching expression (typically a regular
    expression).
    `AbstractPattern` objects can be used to match strings with [`match`](@ref).
"""
abstract type AbstractPattern end

"""
    Regex(pattern[, flags])

A type representing a regular expression. `Regex` objects can be used to match strings
with [`match`](@ref).

`Regex` objects can be created using the [`@r_str`](@ref) string macro. The
`Regex(pattern[, flags])` constructor is usually used if the `pattern` string needs
to be interpolated. See the documentation of the string macro for details on flags.
"""
mutable struct Regex <: AbstractPattern
    pattern::String
    compile_options::UInt32
    match_options::UInt32
    regex::Ptr{Cvoid}

    function Regex(pattern::AbstractString, compile_options::Integer,
                   match_options::Integer)
        pattern = String(pattern)
        compile_options = UInt32(compile_options)
        match_options = UInt32(match_options)
        if (compile_options & ~PCRE.COMPILE_MASK) != 0
            throw(ArgumentError("invalid regex compile options: $compile_options"))
        end
        if (match_options & ~PCRE.EXECUTE_MASK) !=0
            throw(ArgumentError("invalid regex match options: $match_options"))
        end
        re = compile(new(pattern, compile_options, match_options, C_NULL))
        finalizer(re) do re
            re.regex == C_NULL || PCRE.free_re(re.regex)
        end
        re
    end
end

function Regex(pattern::AbstractString, flags::AbstractString)
    options = DEFAULT_COMPILER_OPTS
    for f in flags
        if f == 'a'
            options &= ~PCRE.UCP
        else
            options |= f=='i' ? PCRE.CASELESS  :
                       f=='m' ? PCRE.MULTILINE :
                       f=='s' ? PCRE.DOTALL    :
                       f=='x' ? PCRE.EXTENDED  :
                       throw(ArgumentError("unknown regex flag: $f"))
        end
    end
    Regex(pattern, options, DEFAULT_MATCH_OPTS)
end
Regex(pattern::AbstractString) = Regex(pattern, DEFAULT_COMPILER_OPTS, DEFAULT_MATCH_OPTS)

function compile(regex::Regex)
    if regex.regex == C_NULL
        if PCRE.PCRE_COMPILE_LOCK === nothing
            regex.regex = PCRE.compile(regex.pattern, regex.compile_options)
            PCRE.jit_compile(regex.regex)
        else
            l = PCRE.PCRE_COMPILE_LOCK::Threads.SpinLock
            lock(l)
            try
                if regex.regex == C_NULL
                    regex.regex = PCRE.compile(regex.pattern, regex.compile_options)
                    PCRE.jit_compile(regex.regex)
                end
            finally
                unlock(l)
            end
        end
    end
    regex
end

"""
    @r_str -> Regex

Construct a regex, such as `r"^[a-z]*\$"`, without interpolation and unescaping (except for
quotation mark `"` which still has to be escaped). The regex also accepts one or more flags,
listed after the ending quote, to change its behaviour:

- `i` enables case-insensitive matching
- `m` treats the `^` and `\$` tokens as matching the start and end of individual lines, as
  opposed to the whole string.
- `s` allows the `.` modifier to match newlines.
- `x` enables "comment mode": whitespace is enabled except when escaped with `\\`, and `#`
  is treated as starting a comment.
- `a` disables `UCP` mode (enables ASCII mode). By default `\\B`, `\\b`, `\\D`, `\\d`, `\\S`,
  `\\s`, `\\W`, `\\w`, etc. match based on Unicode character properties. With this option,
  these sequences only match ASCII characters.

See [`Regex`](@ref) if interpolation is needed.

# Examples
```jldoctest
julia> match(r"a+.*b+.*?d\$"ism, "Goodbye,\\nOh, angry,\\nBad world\\n")
RegexMatch("angry,\\nBad world")
```
This regex has the first three flags enabled.
"""
macro r_str(pattern, flags...) Regex(pattern, flags...) end

function show(io::IO, re::Regex)
    imsxa = PCRE.CASELESS|PCRE.MULTILINE|PCRE.DOTALL|PCRE.EXTENDED|PCRE.UCP
    opts = re.compile_options
    if (opts & ~imsxa) == (DEFAULT_COMPILER_OPTS & ~imsxa)
        print(io, "r\"")
        escape_raw_string(io, re.pattern)
        print(io, "\"")
        if (opts & PCRE.CASELESS ) != 0; print(io, 'i'); end
        if (opts & PCRE.MULTILINE) != 0; print(io, 'm'); end
        if (opts & PCRE.DOTALL   ) != 0; print(io, 's'); end
        if (opts & PCRE.EXTENDED ) != 0; print(io, 'x'); end
        if (opts & PCRE.UCP      ) == 0; print(io, 'a'); end
    else
        print(io, "Regex(")
        show(io, re.pattern)
        print(io, ',')
        show(io, opts)
        print(io, ')')
    end
end

"""
   `AbstractMatch` objects are used to represent information about matches found in a string
   using an `AbstractPattern`.
"""
abstract type AbstractMatch end

struct RegexMatch <: AbstractMatch
    match::SubString{String}
    captures::Vector{Union{Nothing,SubString{String}}}
    offset::Int
    offsets::Vector{Int}
    regex::Regex
end

function keys(m::RegexMatch)
    idx_to_capture_name = PCRE.capture_names(m.regex.regex)
    return map(eachindex(m.captures)) do i
        # If the capture group is named, return it's name, else return it's index
        get(idx_to_capture_name, i, i)
    end
end

function show(io::IO, m::RegexMatch)
    print(io, "RegexMatch(")
    show(io, m.match)
    capture_keys = keys(m)
    if !isempty(capture_keys)
        print(io, ", ")
        for (i, capture_name) in enumerate(capture_keys)
            print(io, capture_name, "=")
            show(io, m.captures[i])
            if i < length(m)
                print(io, ", ")
            end
        end
    end
    print(io, ")")
end

# Capture group extraction
getindex(m::RegexMatch, idx::Integer) = m.captures[idx]
function getindex(m::RegexMatch, name::Symbol)
    idx = PCRE.substring_number_from_name(m.regex.regex, name)
    idx <= 0 && error("no capture group named $name found in regex")
    m[idx]
end
getindex(m::RegexMatch, name::AbstractString) = m[Symbol(name)]

haskey(m::RegexMatch, idx::Integer) = idx in eachindex(m.captures)
function haskey(m::RegexMatch, name::Symbol)
    idx = PCRE.substring_number_from_name(m.regex.regex, name)
    return idx > 0
end
haskey(m::RegexMatch, name::AbstractString) = haskey(m, Symbol(name))

iterate(m::RegexMatch, args...) = iterate(m.captures, args...)
length(m::RegexMatch) = length(m.captures)
eltype(m::RegexMatch) = eltype(m.captures)

function occursin(r::Regex, s::AbstractString; offset::Integer=0)
    compile(r)
    return PCRE.exec_r(r.regex, String(s), offset, r.match_options)
end

function occursin(r::Regex, s::SubString; offset::Integer=0)
    compile(r)
    return PCRE.exec_r(r.regex, s, offset, r.match_options)
end

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

Return `true` if `s` starts with the regex pattern, `prefix`.

!!! note
    `startswith` does not compile the anchoring into the regular
    expression, but instead passes the anchoring as
    `match_option` to PCRE. If compile time is amortized,
    `occursin(r"^...", s)` is faster than `startswith(s, r"...")`.

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

!!! compat "Julia 1.2"
    This method requires at least Julia 1.2.

# Examples
```jldoctest
julia> startswith("JuliaLang", r"Julia|Romeo")
true
```
"""
function startswith(s::AbstractString, r::Regex)
    compile(r)
    return PCRE.exec_r(r.regex, String(s), 0, r.match_options | PCRE.ANCHORED)
end

function startswith(s::SubString, r::Regex)
    compile(r)
    return PCRE.exec_r(r.regex, s, 0, r.match_options | PCRE.ANCHORED)
end

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

Return `true` if `s` ends with the regex pattern, `suffix`.

!!! note
    `endswith` does not compile the anchoring into the regular
    expression, but instead passes the anchoring as
    `match_option` to PCRE. If compile time is amortized,
    `occursin(r"...\$", s)` is faster than `endswith(s, r"...")`.

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

!!! compat "Julia 1.2"
    This method requires at least Julia 1.2.

# Examples
```jldoctest
julia> endswith("JuliaLang", r"Lang|Roberts")
true
```
"""
function endswith(s::AbstractString, r::Regex)
    compile(r)
    return PCRE.exec_r(r.regex, String(s), 0, r.match_options | PCRE.ENDANCHORED)
end

function endswith(s::SubString, r::Regex)
    compile(r)
    return PCRE.exec_r(r.regex, s, 0, r.match_options | PCRE.ENDANCHORED)
end

"""
    match(r::Regex, s::AbstractString[, idx::Integer[, addopts]])

Search for the first match of the regular expression `r` in `s` and return a `RegexMatch`
object containing the match, or nothing if the match failed. The matching substring can be
retrieved by accessing `m.match` and the captured sequences can be retrieved by accessing
`m.captures` The optional `idx` argument specifies an index at which to start the search.

# Examples
```jldoctest
julia> rx = r"a(.)a"
r"a(.)a"

julia> m = match(rx, "cabac")
RegexMatch("aba", 1="b")

julia> m.captures
1-element Vector{Union{Nothing, SubString{String}}}:
 "b"

julia> m.match
"aba"

julia> match(rx, "cabac", 3) === nothing
true
```
"""
function match end

function match(re::Regex, str::Union{SubString{String}, String}, idx::Integer,
               add_opts::UInt32=UInt32(0))
    compile(re)
    opts = re.match_options | add_opts
    matched, data = PCRE.exec_r_data(re.regex, str, idx-1, opts)
    if !matched
        PCRE.free_match_data(data)
        return nothing
    end
    n = div(PCRE.ovec_length(data), 2) - 1
    p = PCRE.ovec_ptr(data)
    mat = SubString(str, unsafe_load(p, 1)+1, prevind(str, unsafe_load(p, 2)+1))
    cap = Union{Nothing,SubString{String}}[unsafe_load(p,2i+1) == PCRE.UNSET ? nothing :
                                        SubString(str, unsafe_load(p,2i+1)+1,
                                                  prevind(str, unsafe_load(p,2i+2)+1)) for i=1:n]
    off = Int[ unsafe_load(p,2i+1)+1 for i=1:n ]
    result = RegexMatch(mat, cap, unsafe_load(p,1)+1, off, re)
    PCRE.free_match_data(data)
    return result
end

match(r::Regex, s::AbstractString) = match(r, s, firstindex(s))
match(r::Regex, s::AbstractString, i::Integer) = throw(ArgumentError(
    "regex matching is only available for the String type; use String(s) to convert"
))

findnext(re::Regex, str::Union{String,SubString}, idx::Integer) = _findnext_re(re, str, idx, C_NULL)

# TODO: return only start index and update deprecation
function _findnext_re(re::Regex, str::Union{String,SubString}, idx::Integer, match_data::Ptr{Cvoid})
    if idx > nextind(str,lastindex(str))
        throw(BoundsError())
    end
    opts = re.match_options
    compile(re)
    alloc = match_data == C_NULL
    if alloc
        matched, data = PCRE.exec_r_data(re.regex, str, idx-1, opts)
    else
        matched = PCRE.exec(re.regex, str, idx-1, opts, match_data)
        data = match_data
    end
    if matched
        p = PCRE.ovec_ptr(data)
        ans = (Int(unsafe_load(p,1))+1):prevind(str,Int(unsafe_load(p,2))+1)
    else
        ans = nothing
    end
    alloc && PCRE.free_match_data(data)
    return ans
end
findnext(r::Regex, s::AbstractString, idx::Integer) = throw(ArgumentError(
    "regex search is only available for the String type; use String(s) to convert"
))
findfirst(r::Regex, s::AbstractString) = findnext(r,s,firstindex(s))

"""
    findall(
        pattern::Union{AbstractString,AbstractPattern},
        string::AbstractString;
        overlap::Bool = false,
    )

Return a `Vector{UnitRange{Int}}` of all the matches for `pattern` in `string`.
Each element of the returned vector is a range of indices where the
matching sequence is found, like the return value of [`findnext`](@ref).

If `overlap=true`, the matching sequences are allowed to overlap indices in the
original string, otherwise they must be from disjoint character ranges.

# Examples
```jldoctest
julia> findall("a", "apple")
1-element Vector{UnitRange{Int64}}:
 1:1

julia> findall("nana", "banana")
1-element Vector{UnitRange{Int64}}:
 3:6

julia> findall("a", "banana")
3-element Vector{UnitRange{Int64}}:
 2:2
 4:4
 6:6
```

!!! compat "Julia 1.3"
     This method requires at least Julia 1.3.
"""
function findall(t::Union{AbstractString,AbstractPattern}, s::AbstractString; overlap::Bool=false)
    found = UnitRange{Int}[]
    i, e = firstindex(s), lastindex(s)
    while true
        r = findnext(t, s, i)
        isnothing(r) && break
        push!(found, r)
        j = overlap || isempty(r) ? first(r) : last(r)
        j > e && break
        @inbounds i = nextind(s, j)
    end
    return found
end

"""
    findall(c::AbstractChar, s::AbstractString)

Return a vector `I` of the indices of `s` where `s[i] == c`. If there are no such
elements in `s`, return an empty array.

# Examples
```jldoctest
julia> findall('a', "batman")
2-element Vector{Int64}:
 2
 5
```

!!! compat "Julia 1.7"
     This method requires at least Julia 1.7.
"""
findall(c::AbstractChar, s::AbstractString) = findall(isequal(c),s)

"""
    count(
        pattern::Union{AbstractChar,AbstractString,AbstractPattern},
        string::AbstractString;
        overlap::Bool = false,
    )

Return the number of matches for `pattern` in `string`. This is equivalent to
calling `length(findall(pattern, string))` but more efficient.

If `overlap=true`, the matching sequences are allowed to overlap indices in the
original string, otherwise they must be from disjoint character ranges.

!!! compat "Julia 1.3"
     This method requires at least Julia 1.3.

!!! compat "Julia 1.7"
      Using a character as the pattern requires at least Julia 1.7.
"""
function count(t::Union{AbstractChar,AbstractString,AbstractPattern}, s::AbstractString; overlap::Bool=false)
    n = 0
    i, e = firstindex(s), lastindex(s)
    while true
        r = findnext(t, s, i)
        isnothing(r) && break
        n += 1
        j = overlap || isempty(r) ? first(r) : last(r)
        j > e && break
        @inbounds i = nextind(s, j)
    end
    return n
end

"""
    SubstitutionString(substr)

Stores the given string `substr` as a `SubstitutionString`, for use in regular expression
substitutions. Most commonly constructed using the [`@s_str`](@ref) macro.

```jldoctest
julia> SubstitutionString("Hello \\\\g<name>, it's \\\\1")
s"Hello \\g<name>, it's \\1"

julia> subst = s"Hello \\g<name>, it's \\1"
s"Hello \\g<name>, it's \\1"

julia> typeof(subst)
SubstitutionString{String}

```

"""
struct SubstitutionString{T<:AbstractString} <: AbstractString
    string::T
end

ncodeunits(s::SubstitutionString) = ncodeunits(s.string)::Int
codeunit(s::SubstitutionString) = codeunit(s.string)::CodeunitType
codeunit(s::SubstitutionString, i::Integer) = codeunit(s.string, i)::Union{UInt8, UInt16, UInt32}
isvalid(s::SubstitutionString, i::Integer) = isvalid(s.string, i)::Bool
iterate(s::SubstitutionString, i::Integer...) = iterate(s.string, i...)::Union{Nothing,Tuple{AbstractChar,Int}}

function show(io::IO, s::SubstitutionString)
    print(io, "s\"")
    escape_raw_string(io, s.string)
    print(io, "\"")
end

"""
    @s_str -> SubstitutionString

Construct a substitution string, used for regular expression substitutions.  Within the
string, sequences of the form `\\N` refer to the Nth capture group in the regex, and
`\\g<groupname>` refers to a named capture group with name `groupname`.

```jldoctest
julia> msg = "#Hello# from Julia";

julia> replace(msg, r"#(.+)# from (?<from>\\w+)" => s"FROM: \\g<from>; MESSAGE: \\1")
"FROM: Julia; MESSAGE: Hello"
```
"""
macro s_str(string) SubstitutionString(string) end

# replacement

struct RegexAndMatchData
    re::Regex
    match_data::Ptr{Cvoid}
    RegexAndMatchData(re::Regex) = (compile(re); new(re, PCRE.create_match_data(re.regex)))
end

findnext(pat::RegexAndMatchData, str, i) = _findnext_re(pat.re, str, i, pat.match_data)

_pat_replacer(r::Regex) = RegexAndMatchData(r)

_free_pat_replacer(r::RegexAndMatchData) = PCRE.free_match_data(r.match_data)

replace_err(repl) = error("Bad replacement string: $repl")

function _write_capture(io, re::RegexAndMatchData, group)
    len = PCRE.substring_length_bynumber(re.match_data, group)
    ensureroom(io, len+1)
    PCRE.substring_copy_bynumber(re.match_data, group,
        pointer(io.data, io.ptr), len+1)
    io.ptr += len
    io.size = max(io.size, io.ptr - 1)
end

const SUB_CHAR = '\\'
const GROUP_CHAR = 'g'
const KEEP_ESC = [SUB_CHAR, GROUP_CHAR, '0':'9'...]

function _replace(io, repl_s::SubstitutionString, str, r, re::RegexAndMatchData)
    LBRACKET = '<'
    RBRACKET = '>'
    repl = unescape_string(repl_s.string, KEEP_ESC)
    i = firstindex(repl)
    e = lastindex(repl)
    while i <= e
        if repl[i] == SUB_CHAR
            next_i = nextind(repl, i)
            next_i > e && replace_err(repl)
            if repl[next_i] == SUB_CHAR
                write(io, SUB_CHAR)
                i = nextind(repl, next_i)
            elseif isdigit(repl[next_i])
                group = parse(Int, repl[next_i])
                i = nextind(repl, next_i)
                while i <= e
                    if isdigit(repl[i])
                        group = 10group + parse(Int, repl[i])
                        i = nextind(repl, i)
                    else
                        break
                    end
                end
                _write_capture(io, re, group)
            elseif repl[next_i] == GROUP_CHAR
                i = nextind(repl, next_i)
                if i > e || repl[i] != LBRACKET
                    replace_err(repl)
                end
                i = nextind(repl, i)
                i > e && replace_err(repl)
                groupstart = i
                while repl[i] != RBRACKET
                    i = nextind(repl, i)
                    i > e && replace_err(repl)
                end
                #  TODO: avoid this allocation
                groupname = SubString(repl, groupstart, prevind(repl, i))
                if all(isdigit, groupname)
                    _write_capture(io, re, parse(Int, groupname))
                else
                    group = PCRE.substring_number_from_name(re.re.regex, groupname)
                    group < 0 && replace_err("Group $groupname not found in regex $(re.re)")
                    _write_capture(io, re, group)
                end
                i = nextind(repl, i)
            else
                replace_err(repl)
            end
        else
            write(io, repl[i])
            i = nextind(repl, i)
        end
    end
end

struct RegexMatchIterator
    regex::Regex
    string::String
    overlap::Bool

    function RegexMatchIterator(regex::Regex, string::AbstractString, ovr::Bool=false)
        new(regex, string, ovr)
    end
end
compile(itr::RegexMatchIterator) = (compile(itr.regex); itr)
eltype(::Type{RegexMatchIterator}) = RegexMatch
IteratorSize(::Type{RegexMatchIterator}) = SizeUnknown()

function iterate(itr::RegexMatchIterator, (offset,prevempty)=(1,false))
    opts_nonempty = UInt32(PCRE.ANCHORED | PCRE.NOTEMPTY_ATSTART)
    while true
        mat = match(itr.regex, itr.string, offset,
                    prevempty ? opts_nonempty : UInt32(0))

        if mat === nothing
            if prevempty && offset <= sizeof(itr.string)
                offset = nextind(itr.string, offset)
                prevempty = false
                continue
            else
                break
            end
        else
            if itr.overlap
                if !isempty(mat.match)
                    offset = nextind(itr.string, mat.offset)
                else
                    offset = mat.offset
                end
            else
                offset = mat.offset + ncodeunits(mat.match)
            end
            return (mat, (offset, isempty(mat.match)))
        end
    end
    nothing
end

"""
    eachmatch(r::Regex, s::AbstractString; overlap::Bool=false)

Search for all matches of the regular expression `r` in `s` and return an iterator over the
matches. If `overlap` is `true`, the matching sequences are allowed to overlap indices in the
original string, otherwise they must be from distinct character ranges.

# Examples
```jldoctest
julia> rx = r"a.a"
r"a.a"

julia> m = eachmatch(rx, "a1a2a3a")
Base.RegexMatchIterator(r"a.a", "a1a2a3a", false)

julia> collect(m)
2-element Vector{RegexMatch}:
 RegexMatch("a1a")
 RegexMatch("a3a")

julia> collect(eachmatch(rx, "a1a2a3a", overlap = true))
3-element Vector{RegexMatch}:
 RegexMatch("a1a")
 RegexMatch("a2a")
 RegexMatch("a3a")
```
"""
eachmatch(re::Regex, str::AbstractString; overlap = false) =
    RegexMatchIterator(re, str, overlap)

## comparison ##

function ==(a::Regex, b::Regex)
    a.pattern == b.pattern && a.compile_options == b.compile_options && a.match_options == b.match_options
end

## hash ##
const hashre_seed = UInt === UInt64 ? 0x67e195eb8555e72d : 0xe32373e4
function hash(r::Regex, h::UInt)
    h += hashre_seed
    h = hash(r.pattern, h)
    h = hash(r.compile_options, h)
    h = hash(r.match_options, h)
end

## String operations ##

"""
    *(s::Regex, t::Union{Regex,AbstractString,AbstractChar}) -> Regex
    *(s::Union{Regex,AbstractString,AbstractChar}, t::Regex) -> Regex

Concatenate regexes, strings and/or characters, producing a [`Regex`](@ref).
String and character arguments must be matched exactly in the resulting regex,
meaning that the contained characters are devoid of any special meaning
(they are quoted with "\\Q" and "\\E").

!!! compat "Julia 1.3"
    This method requires at least Julia 1.3.

# Examples
```jldoctest
julia> match(r"Hello|Good bye" * ' ' * "world", "Hello world")
RegexMatch("Hello world")

julia> r = r"a|b" * "c|d"
r"(?:a|b)\\Qc|d\\E"

julia> match(r, "ac") == nothing
true

julia> match(r, "ac|d")
RegexMatch("ac|d")
```
"""
function *(r1::Union{Regex,AbstractString,AbstractChar}, rs::Union{Regex,AbstractString,AbstractChar}...)
    mask = PCRE.CASELESS | PCRE.MULTILINE | PCRE.DOTALL | PCRE.EXTENDED # imsx
    match_opts   = nothing # all args must agree on this
    compile_opts = nothing # all args must agree on this
    shared = mask
    for r in (r1, rs...)
        r isa Regex || continue
        if match_opts === nothing
            match_opts = r.match_options
            compile_opts = r.compile_options & ~mask
        else
            r.match_options == match_opts &&
                r.compile_options & ~mask == compile_opts ||
                throw(ArgumentError("cannot multiply regexes: incompatible options"))
        end
        shared &= r.compile_options
    end
    unshared = mask & ~shared
    Regex(string(wrap_string(r1, unshared), wrap_string.(rs, Ref(unshared))...), compile_opts | shared, match_opts)
end

*(r::Regex) = r # avoids wrapping r in a useless subpattern

wrap_string(r::Regex, unshared::UInt32) = string("(?", regex_opts_str(r.compile_options & unshared), ':', r.pattern, ')')
# if s contains raw"\E", split '\' and 'E' within two distinct \Q...\E groups:
wrap_string(s::AbstractString, ::UInt32) =  string("\\Q", replace(s, raw"\E" => raw"\\E\QE"), "\\E")
wrap_string(s::AbstractChar, ::UInt32) = string("\\Q", s, "\\E")

regex_opts_str(opts) = (isassigned(_regex_opts_str) ? _regex_opts_str[] : init_regex())[opts]

# UInt32 to String mapping for some compile options
const _regex_opts_str = Ref{ImmutableDict{UInt32,String}}()

@noinline init_regex() = _regex_opts_str[] = foldl(0:15, init=ImmutableDict{UInt32,String}()) do d, o
    opt = UInt32(0)
    str = ""
    if o & 1 != 0
        opt |= PCRE.CASELESS
        str *= 'i'
    end
    if o & 2 != 0
        opt |= PCRE.MULTILINE
        str *= 'm'
    end
    if o & 4 != 0
        opt |= PCRE.DOTALL
        str *= 's'
    end
    if o & 8 != 0
        opt |= PCRE.EXTENDED
        str *= 'x'
    end
    ImmutableDict(d, opt => str)
end

"""
    ^(s::Regex, n::Integer)

Repeat a regex `n` times.

!!! compat "Julia 1.3"
    This method requires at least Julia 1.3.

# Examples
```jldoctest
julia> r"Test "^2
r"(?:Test ){2}"

julia> match(r"Test "^2, "Test Test ")
RegexMatch("Test Test ")
```
"""
^(r::Regex, i::Integer) = Regex(string("(?:", r.pattern, "){$i}"), r.compile_options, r.match_options)

max_stars_repo_name: ranocha/julia

max_issues_repo_name: ranocha/julia

max_forks_repo_name: ranocha/julia


text:
## Functions to compute the reduced shape

# for reductions that expand 0 dims to 1
reduced_dims(a::AbstractArray, region) = reduced_dims(size(a), region)

# for reductions that keep 0 dims as 0
reduced_dims0(a::AbstractArray, region) = reduced_dims0(size(a), region)
reduced_dims{N}(siz::NTuple{N,Int}, d::Int, rd::Int) = (d == 1 ? tuple(rd, siz[d+1:N]...) :
                                                        d == N ? tuple(siz[1:N-1]..., rd) :
                                                        1 < d < N ? tuple(siz[1:d-1]..., rd, siz[d+1:N]...) :
                                                        siz)::typeof(siz)
reduced_dims{N}(siz::NTuple{N,Int}, d::Int) = reduced_dims(siz, d, 1)
reduced_dims0{N}(siz::NTuple{N,Int}, d::Int) = 1 <= d <= N ? reduced_dims(siz, d, (siz[d] == 0 ? 0 : 1)) : siz

function reduced_dims{N}(siz::NTuple{N,Int}, region)
    rsiz = [siz...]
    for i in region
        if 1 <= i <= N
            rsiz[i] = 1
        end
    end
    tuple(rsiz...)::typeof(siz)
end

function reduced_dims0{N}(siz::NTuple{N,Int}, region)
    rsiz = [siz...]
    for i in region
        if i <= i <= N
            rsiz[i] = (rsiz[i] == 0 ? 0 : 1)
        end
    end
    tuple(rsiz...)::typeof(siz)
end

function regionsize(a, region)
    s = 1
    for d in region
        s *= size(a,d)
    end
    s
end

###### Generic reduction functions #####

## initialization

for (Op, initfun) in ((:AddFun, :zero), (:MulFun, :one), (:MaxFun, :typemin), (:MinFun, :typemax))
    @eval initarray!{T}(a::AbstractArray{T}, ::$(Op), init::Bool) = (init && fill!(a, $(initfun)(T)); a)
end

for (Op, initval) in ((:AndFun, true), (:OrFun, false))
    @eval initarray!(a::AbstractArray, ::$(Op), init::Bool) = (init && fill!(a, $initval); a)
end

reducedim_initarray{R}(A::AbstractArray, region, v0, ::Type{R}) = fill!(similar(A,R,reduced_dims(A,region)), v0)
reducedim_initarray{T}(A::AbstractArray, region, v0::T) = reducedim_initarray(A, region, v0, T)

reducedim_initarray0{R}(A::AbstractArray, region, v0, ::Type{R}) = fill!(similar(A,R,reduced_dims0(A,region)), v0)
reducedim_initarray0{T}(A::AbstractArray, region, v0::T) = reducedim_initarray0(A, region, v0, T)

# TODO: better way to handle reducedim initialization
#
# The current scheme is basically following Steven G. Johnson's original implementation
#
function reducedim_init{T}(f, op::AddFun, A::AbstractArray{T}, region)
    if method_exists(zero, (Type{T},))
        x = evaluate(f, zero(T))
        z = zero(x) + zero(x)
        Tr = typeof(z) == typeof(x) && !isbits(T) ? T : typeof(z)
    else
        z = zero(sum(f, A))
        Tr = typeof(z)
    end
    return reducedim_initarray(A, region, z, Tr)
end

function reducedim_init{T}(f, op::MulFun, A::AbstractArray{T}, region)
    if method_exists(zero, (Type{T},))
        x = evaluate(f, zero(T))
        z = one(x) * one(x)
        Tr = typeof(z) == typeof(x) && !isbits(T) ? T : typeof(z)
    else
        z = one(prod(f, A))
        Tr = typeof(z)
    end
    return reducedim_initarray(A, region, z, Tr)
end

reducedim_init{T}(f, op::MaxFun, A::AbstractArray{T}, region) = reducedim_initarray0(A, region, typemin(evaluate(f, zero(T))))
reducedim_init{T}(f, op::MinFun, A::AbstractArray{T}, region) = reducedim_initarray0(A, region, typemax(evaluate(f, zero(T))))
reducedim_init{T}(f::Union(AbsFun,Abs2Fun), op::MaxFun, A::AbstractArray{T}, region) =
    reducedim_initarray(A, region, zero(evaluate(f, zero(T))))

reducedim_init(f, op::AndFun, A::AbstractArray, region) = reducedim_initarray(A, region, true)
reducedim_init(f, op::OrFun, A::AbstractArray, region) = reducedim_initarray(A, region, false)

# specialize to make initialization more efficient for common cases

typealias CommonReduceResult Union(Uint64,Uint128,Int64,Int128,Float32,Float64,Complex64,Complex128)

for (IT, RT) in ((:CommonReduceResult, :T), (:SmallSigned, :Int), (:SmallUnsigned, :Uint))
    @eval begin
        reducedim_init{T<:$IT}(f::Union(IdFun,AbsFun,Abs2Fun), op::AddFun, A::AbstractArray{T}, region) =
            reducedim_initarray(A, region, zero($RT))
        reducedim_init{T<:$IT}(f::Union(IdFun,AbsFun,Abs2Fun), op::MulFun, A::AbstractArray{T}, region) =
            reducedim_initarray(A, region, one($RT))
    end
end
reducedim_init(f::Union(IdFun,AbsFun,Abs2Fun), op::AddFun, A::AbstractArray{Bool}, region) =
    reducedim_initarray(A, region, 0)

## generic (map)reduction

has_fast_linear_indexing(a::AbstractArray) = false
has_fast_linear_indexing(a::Array) = true

function check_reducdims(R, A)
    # Check whether R has compatible dimensions w.r.t. A for reduction
    #
    # It returns an integer value value (useful for choosing implementation)
    # - If it reduces only along leading dimensions, e.g. sum(A, 1) or sum(A, (1, 2)),
    #   it returns the length of the leading slice. For the two examples above,
    #   it will be size(A, 1) or size(A, 1) * size(A, 2).
    # - Otherwise, e.g. sum(A, 2) or sum(A, (1, 3)), it returns 0.
    #
    lsiz = 1
    had_nonreduc = false
    for i = 1:ndims(A)
        sRi = size(R, i)
        sAi = size(A, i)
        if sRi == 1
            if sAi > 1
                if had_nonreduc
                    lsiz = 0  # to reduce along i, but some previous dimensions were non-reducing
                else
                    lsiz *= sAi  # if lsiz was set to zero, it will stay to be zero
                end
            end
        else
            sRi == sAi ||
                throw(DimensionMismatch("Reduction on array of size $(size(A)) with output of size $(size(R))"))
            had_nonreduc = true
        end
    end
    return lsiz
end

@ngenerate N typeof(R) function _mapreducedim!{T,N}(f, op, R::AbstractArray, A::AbstractArray{T,N})
    lsiz = check_reducdims(R, A)
    isempty(A) && return R
    @nextract N sizeR d->size(R,d)
    sizA1 = size(A, 1)

    if has_fast_linear_indexing(A) && lsiz > 16
        # use mapreduce_impl, which is probably better tuned to achieve higher performance
        nslices = div(length(A), lsiz)
        ibase = 0
        for i = 1:nslices
            @inbounds R[i] = mapreduce_impl(f, op, A, ibase+1, ibase+lsiz)
            ibase += lsiz
        end
    elseif size(R, 1) == 1 && sizA1 > 1
        # keep the accumulator as a local variable when reducing along the first dimension
        @nloops N i d->(d>1? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin
            @inbounds r = (@nref N R j)
            for i_1 = 1:sizA1
                @inbounds v = evaluate(f, (@nref N A i))
                r = evaluate(op, r, v)
            end
            @inbounds (@nref N R j) = r
        end
    else
        # general implementation
        @nloops N i A d->(j_d = sizeR_d==1 ? 1 : i_d) begin
            @inbounds v = evaluate(f, (@nref N A i))
            @inbounds (@nref N R j) = evaluate(op, (@nref N R j), v)
        end
    end
    return R
end

mapreducedim!(f, op, R::AbstractArray, A::AbstractArray) = _mapreducedim!(f, op, R, A)

function mapreducedim!(f::Function, op, R::AbstractArray, A::AbstractArray)
    is(op, +) ? _mapreducedim!(f, AddFun(), R, A) :
    is(op, *) ? _mapreducedim!(f, MulFun(), R, A) :
    is(op, &) ? _mapreducedim!(f, AndFun(), R, A) :
    is(op, |) ? _mapreducedim!(f, OrFun(), R, A) :
    _mapreducedim!(f, op, R, A)
end

reducedim!{RT}(op, R::AbstractArray{RT}, A::AbstractArray) = mapreducedim!(IdFun(), op, R, A, zero(RT))

mapreducedim(f, op, A::AbstractArray, region, v0) = mapreducedim!(f, op, reducedim_initarray(A, region, v0), A)
mapreducedim{T}(f, op, A::AbstractArray{T}, region) = mapreducedim!(f, op, reducedim_init(f, op, A, region), A)

reducedim(op, A::AbstractArray, region, v0) = mapreducedim(IdFun(), op, A, region, v0)
reducedim(op, A::AbstractArray, region) = mapreducedim(IdFun(), op, A, region)

##### Specific reduction functions #####

for (fname, Op) in [(:sum, :AddFun), (:prod, :MulFun),
                    (:maximum, :MaxFun), (:minimum, :MinFun),
                    (:all, :AndFun), (:any, :OrFun)]

    fname! = symbol(string(fname, '!'))
    @eval begin
        $(fname!)(f::Union(Function,Func{1}), r::AbstractArray, A::AbstractArray; init::Bool=true) =
            mapreducedim!(f, $(Op)(), initarray!(r, $(Op)(), init), A)
        $(fname!)(r::AbstractArray, A::AbstractArray; init::Bool=true) = $(fname!)(IdFun(), r, A; init=init)

        $(fname)(f::Union(Function,Func{1}), A::AbstractArray, region) =
            mapreducedim(f, $(Op)(), A, region)
        $(fname)(A::AbstractArray, region) = $(fname)(IdFun(), A, region)
    end
end

for (fname, fbase, Fun) in [(:sumabs, :sum, :AbsFun),
                            (:sumabs2, :sum, :Abs2Fun),
                            (:maxabs, :maximum, :AbsFun),
                            (:minabs, :minimum, :AbsFun)]
    fname! = symbol(string(fname, '!'))
    fbase! = symbol(string(fbase, '!'))
    @eval begin
        $(fname!)(r::AbstractArray, A::AbstractArray; init::Bool=true) =
            $(fbase!)($(Fun)(), r, A; init=init)
        $(fname)(A::AbstractArray, region) = $(fbase)($(Fun)(), A, region)
    end
end

##### findmin & findmax #####

# Generate the body for a reduction function reduce!(f, Rval, Rind, A), using a comparison operator f
# Rind contains the index of A from which Rval was taken
function gen_findreduction_body(N, f::Function)
    F = Expr(:quote, f)
    quote
        (isempty(Rval) || isempty(A)) && return Rval, Rind
        for i = 1:$N
            (size(Rval, i) == size(A, i) || size(Rval, i) == 1) || throw(DimensionMismatch("Find-reduction on array of size $(size(A)) with output of size $(size(Rval))"))
            size(Rval, i) == size(Rind, i) || throw(DimensionMismatch("Find-reduction: outputs must be of the same size"))
        end
        @nexprs $N d->(sizeR_d = size(Rval,d))
        # If we're reducing along dimension 1, for efficiency we can make use of a temporary.
        # Otherwise, keep the result in Rval/Rind so that we traverse A in storage order.
        k = 0
        @inbounds if size(Rval, 1) < size(A, 1)
            @nloops $N i d->(d>1? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin
                tmpRv = (@nref $N Rval j)
                tmpRi = (@nref $N Rind j)
                for i_1 = 1:size(A,1)
                    k += 1
                    tmpAv = (@nref $N A i)
                    if ($F)(tmpAv, tmpRv)
                        tmpRv = tmpAv
                        tmpRi = k
                    end
                end
                (@nref $N Rval j) = tmpRv
                (@nref $N Rind j) = tmpRi
            end
        else
            @nloops $N i A d->(j_d = sizeR_d==1 ? 1 : i_d) begin
                k += 1
                tmpAv = (@nref $N A i)
                if ($F)(tmpAv, (@nref $N Rval j))
                    (@nref $N Rval j) = tmpAv
                    (@nref $N Rind j) = k
                end
            end
        end
        Rval, Rind
    end
end

eval(ngenerate(:N, :(typeof((Rval,Rind))), :(_findmin!{T,N}(Rval::AbstractArray, Rind::AbstractArray, A::AbstractArray{T,N})), N->gen_findreduction_body(N, <)))
findmin!{R}(rval::AbstractArray{R}, rind::AbstractArray, A::AbstractArray; init::Bool=true) = _findmin!(initarray!(rval, typemax(R), init), rind, A)
findmin{T}(A::AbstractArray{T}, region) =
    isempty(A) ? (similar(A,reduced_dims0(A,region)), zeros(Int,reduced_dims0(A,region))) :
                  _findmin!(reducedim_initarray0(A, region, typemax(T)), zeros(Int,reduced_dims0(A,region)), A)

eval(ngenerate(:N, :(typeof((Rval,Rind))), :(_findmax!{T,N}(Rval::AbstractArray, Rind::AbstractArray, A::AbstractArray{T,N})), N->gen_findreduction_body(N, >)))
findmax!{R}(rval::AbstractArray{R}, rind::AbstractArray, A::AbstractArray; init::Bool=true) = _findmax!(initarray!(rval, typemin(R), init), rind, A)
findmax{T}(A::AbstractArray{T}, region) =
    isempty(A) ? (similar(A,reduced_dims0(A,region)), zeros(Int,reduced_dims0(A,region))) :
                  _findmax!(reducedim_initarray0(A, region, typemin(T)), zeros(Int,reduced_dims0(A,region)), A)

ext: jl

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Test
using CompilerPluginTools

module TestIntrinsic
using CompilerPluginTools

@intrinsic_stub device main(gate)
@intrinsic_stub device gate(gate, loc::Int)
@intrinsic_stub device ctrl(gate, loc::Int, ctrl::Int)
end

@test TestIntrinsic.isintrinsic(TestIntrinsic.main)
@test TestIntrinsic.isintrinsic(TestIntrinsic.ctrl)
@test TestIntrinsic.isintrinsic(TestIntrinsic.gate)

@test_throws IntrinsicError TestIntrinsic.main(:gate)
@test_throws IntrinsicError TestIntrinsic.gate(:gate, 1)
@test_throws IntrinsicError TestIntrinsic.ctrl(:gate, 1, 2)


text:
function Array{NamedD{Srcs},1}(grid::StepRangeLen, ageom::AGeom, fields::Vector{Symbol})
	return [NamedD(grid,Srcs(ageom[i].ns),fields) for i in 1:length(ageom)]
end

"""
A mutable type that bundles multi-component source wavelets.
```julia
srcwav=SrcWav(tgrid, ageom, [:P, :Vx])

```
Here, we initialized wavelets, for `:P` and `:Vx` fields, in time domain for sources and supersources in `ageom`.

## Indexing
* `srcwav[i]` : wavelets for ith supersource of all fields
* `srcwav[i][:P]` : extract field `:P`
* `srcwav[i][:ns]` : number of sources (same as in `ageom`)
* `srcwav.grid` : returns `tgrid`
As mutable objects in Julia are like containers that might hold different values over time, `srcwav` can be modified.
"""
SrcWav=Array{NamedD{Srcs},1}

function issimilar(ageom::AGeom, srcwav::SrcWav)
	test=[]
	push!(test, length(ageom)==length(srcwav))
	nss=length(ageom)
	push!(test, [ageom[i].ns for i in 1:nss]==[srcwav[i][:n] for i in 1:nss])
	return all(test)
end
issimilar(srcwav::SrcWav, ageom::AGeom)=issimilar(ageom, srcwav)

#=

#"""
#Allocate `SrcWav` with zeros depending on the acquisition ageometry.
#"""
#function SrcWav_zeros(ageom::AGeom,  fields::Vector{Symbol}, tgrid::StepRangeLen)
#	wavsrc = [zeros(length(tgrid),ageom.ns[iss]) for iss=1:ageom.nss, ifield=1:length(fields)]
#	return SrcWav(ageom.nss, ageom.ns, fields, wavsrc, deepcopy(tgrid))
#end

"""
Return minimum, maximum and peak frequencies of `SrcWav`
"""
#function freqs(src::SrcWav)
#	freqmin = minimum([Utils.findfreq(src.wav[i,j][:,:],src.tgrid,attrib=:min) for i in 1:src.nss, j in 1:length(src.fields)])
#	freqmax = maximum([Utils.findfreq(src.wav[i,j][:,:],src.tgrid,attrib=:max) for i in 1:src.nss, j in 1:length(src.fields)])
#	freqpeak = mean([Utils.findfreq(src.wav[i,j][:,:],src.tgrid,attrib=:peak) for i in 1:src.nss, j in 1:length(src.fields)])
#	return freqmin, freqmax, freqpeak
#end

"""
Constructor for `SrcWav` data type.
Uses same source wavelet, i.e., `wav` for all sources and supersources

# Arguments

* `nss::Int64` : number of supersources
* `ns::Int64` : number of sources
* `fields::Vector{Symbol}` : number of fields the sources are exciting
* `wav::Array{Float64}` : a source wavelet that is used for all sources and supersources
* `tgrid` : time grid for the wavelet
"""

#=
"""
Constructor of `SrcWav`, which is typical for a input model such that
the model has `nλ` wavelengths.

# Arguments

* `nss::Int64` : number of supersources
* `ns::Int64` : number of source per each supersource
* `fields::Vector{Symbol}` :

# Keyword Arguments

* `mod::Medium` :
* `nλ::Int64=10` : number of wavelengths in the mod
* `wav_func::Function=(fqdom, tgrid)->Utils.Wavelets.ricker(fqdom,tgrid)` : which wavelet to generate, see Utils.Wavelets.jl
* `tmaxfrac::Float64=1.0` : by default the maximum modelling time is computed using the average velocity and the diagonal distance of the model, use this fraction to increase of reduce the maximum time
"""
function update!(s::Datat,
		fields::Vector{Symbol},
		mod::Medium,
		nλ::Int64,
		wav_func::Function=(fqdom, tgrid)->Utils.Wavelets.ricker(fqdom,tgrid),
		tmaxfrac::Float64=1.0
		)

	x=mod.mgrid[2]; z=mod.mgrid[1]
	# dominant wavelength using mod dimensions
	λdom=mean([(abs(x[end]-x[1])), (abs(z[end]-z[1]))])/real(nλ)
	# average P velocity
	vavg=mean(mod[:vp])

	fqdom = vavg/λdom

	# maximum distance the wave travels
	d = sqrt((x[1]-x[end]).^2+(z[1]-z[end]).^2)

	# use two-way maximum distance to get tmax
	tmax=2.0*d*inv(vavg)*tmaxfrac

	# choose sampling interval to obey max freq of source wavelet
	δmin = minimum([step(mod.mgrid[2]), step(mod.mgrid[1])])
	vpmax = mod.bounds[:vp][2]
	δt=0.5*δmin/vpmax

	# check if δt is reasonable
	#(δt > 0.1/fqdom) : error("decrease spatial sampling or nλ")

	tgrid=range(0.0, stop=tmax, step=δt)
	wav = wav_func(fqdom, tgrid);

	wavsrc = [hcat([wav_func(fqdom, tgrid) for is in 1:ns]...) for iss=1:nss, ifield=1:length(fields)]
	src=SrcWav(nss, fill(ns, nss), fields, wavsrc, deepcopy(tgrid))
	print(src)
	# choose ricker waveletes of fdom
	return src
end
=#

#=

"""
Generate band-limited random source signals
"""
function SrcWav_fixed_random(nss::Int64, ns::Int64, fields::Vector{Symbol};
			  distvec=[Normal() for iss in 1:nss],
			  sparsepvec=[1. for iss in 1:nss],
			  fmin::Float64=0.0,
			  fmax::Float64=0.0,
			  tgrid::StepRangeLen=nothing,
			  tmaxfrac::Float64=1.0)
	wavsrc = [repeat(zeros(length(tgrid)),inner=(1,ns)) for iss=1:nss, ifield=1:length(fields)]
	for ifield in 1:length(fields), iss in 1:nss, is in 1:ns
		wavsrc[iss, ifield][:,is] = Utils.get_tapered_random_tmax_signal(tgrid, fmin=fmin, fmax=fmax, tmaxfrac=tmaxfrac, dist=distvec[iss],
								 sparsep=sparsepvec[iss])
	end
	src=SrcWav(nss, fill(ns, nss), fields, wavsrc, deepcopy(tgrid))
	print(src)
	return src
end

"""
Function that returns SrcWav after time reversal
"""
function SrcWav_tr(src::SrcWav)
	return SrcWav(src.nss,src.ns,src.fields,
	    [reverse(src.wav[i,j],dims=1) for i in 1:src.nss, j in 1:length(src.fields)],deepcopy(src.tgrid))
end

"""
Pad `SrcWav`
tgrids should be same in all SrcWav
"""
function SrcWav_uspos(src::Vector{SrcWav}, acq::Vector{AGeom})
	np = length(src) == length(acq) ? length(src) : error("unequal sizez")

	# unique source positions
	nus=AGeom_get(acq,:nus)
	uspos=AGeom_get(acq,:uspos)
	# all zeros for all unique positions
	wavout = [[zeros(src[ip].length(tgrid),nus) for iss=1:src[ip].nss, ifield=1:length(src[ip].fields)] for ip=1:np]
	# fill source wavelets when necessary
	for ip=1:np
		for ifield=1:length(src[ip].fields), iss=1:src[ip].nss, is=1:acq[ip].ns[iss]
			is0=findall([[uspos[1][i]-acq[ip].sz[iss][is],
		      uspos[2][i]-acq[ip].sx[iss][is]] == [0., 0.,] for i in 1:nus])

			wavout[ip][iss, ifield][:,is0] = src[ip].wav[iss, ifield][:,is]
		end
	end
	# output src
	return [SrcWav(src[ip].nss,fill(nus, src[ip].nss),
	     src[ip].fields,wavout[ip],src[ip].tgrid) for ip=1:np]
end

"""
return a vector of the order

"""
function SrcWav_getvec(src::Vector{SrcWav}, field::Symbol)
	np = length(src);
	vect = [getfield(src[ip],field) for ip=1:np]
	return vec(hcat(hcat(vect...)...))
end

=#
=#


text:
export Xcorr
#TODO make more efficient

"""
`Xcorr([domainType=Float64::Type,] dim_in::Tuple, h::AbstractVector)`

`Xcorr(x::AbstractVector, h::AbstractVector)`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns the cross correlation between `x` and `h`. Uses `xcross`.

"""
struct Xcorr{T,H <:AbstractVector{T}} <: LinearOperator
	dim_in::Tuple{Int}
	h::H
end

# Constructors
function Xcorr(DomainType::Type, DomainDim::NTuple{N,Int},  h::H) where {H<:AbstractVector, N}
	eltype(h) != DomainType && error("eltype(h) is $(eltype(h)), should be $(DomainType)")
	N != 1 && error("Xcorr treats only SISO, check Filt and MIMOFilt for MIMO")
	Xcorr{DomainType,H}(DomainDim,h)
end
Xcorr(x::H, h::H) where {H} = Xcorr(eltype(x), size(x), h)

# Mappings

function mul!(y::H,A::Xcorr{T,H},b::H) where {T,H}
	y .= xcorr(b,A.h)
end

function mul!(y::H,L::AdjointOperator{Xcorr{T,H}},b::H) where {T,H}
    A = L.A
	l =floor(Int64,size(A,1)[1]/2)
	idx = l+1:l+length(y)
	y .= conv(b,A.h)[idx]
end

# Properties

domainType(L::Xcorr{T}) where {T} = T
codomainType(L::Xcorr{T}) where {T} = T

#TODO find out a way to verify this,
is_full_row_rank(L::Xcorr)    = true
is_full_column_rank(L::Xcorr) = true

size(L::Xcorr) = ( 2*max(L.dim_in[1], length(L.h))-1, ), L.dim_in

fun_name(A::Xcorr)  = "◎"


text:
"""
# PolyFields.jl
Welcome to PolyFields.jl!
"""
module PolyFields

#==============================================================================#
# Dependencies
#==============================================================================#

import Statistics: mean, std
import LinearAlgebra: norm, dot, det, mul!, inv!, lu!
import Dates

using Printf
using JSON
using Parameters
using RecipesBase

using Random
using StaticArrays
using Interpolations
using FFTW

#==============================================================================#
# Exports
#==============================================================================#

# Types
export AbstractSystem, AbstractSpecies, AbstractConstraint, AbstractPropagator, AbstractInteraction, AbstractFieldUpdater
export FieldGrid, PseudoSpectralSolver, FFTHolder, UnitCell
export Monomer, Point, Homopolymer, Diblock, Multiblock
export FloryInteraction, EdwardsInteraction
export Ensemble, Canonical, Grand, Compressibility
export FieldSystem, SCFTOptions, SCFTResults
export DomainUpdater, EulerUpdater, AndersonUpdater, MomentumUpdater

# Methods
export make_basis, basis_matrix, basis_dimensions
export add_monomer!, add_species!, add_cell!, add_interaction!, set_interaction!, validate
export density!, muphi!, potentials!, potential!, uniform_fields!, residuals!
export free_energy, free_energy_bulk, scfstress, scfstress!, field_error, stress_error
export step!, scft!, fieldinit!, savefields, loadfields, scft!
export energy, energy_bulk, monomer_fraction, chain_grid
export updater_with_system

#==============================================================================#
# Load files
#==============================================================================#

# General
include("types.jl")
include("utils.jl")

# Cell and system
include("system/fft.jl")
include("system/mde.jl")
include("system/monomer.jl")
include("system/basis.jl")
include("system/cell.jl")
include("system/compress.jl")
include("system/system.jl")

# Species architectures
include("species/base.jl")
include("species/point.jl")
include("species/multiblock.jl")

# # Interaction classes
include("interactions/base.jl")
include("interactions/flory.jl")
include("interactions/edwards.jl")

# # Field updaters
include("updaters/base.jl")
include("updaters/euler.jl")
include("updaters/anderson.jl")
include("updaters/momentum.jl")
include("updaters/domain.jl")

include("fileio.jl")
include("simulate.jl")
# # include("viz.jl")

end


text:
const LOCAL_REPO_NAME = "REPO"

function get_project_deps(
    api::Forge,
    ci::CIService,
    repo::Union{GitHub.Repo,GitLab.Project};
    options::Options,
    subdir::AbstractString,
)
    mktempdir() do f
        url_with_auth = get_url_with_auth(api, ci, repo)
        local_path = joinpath(f, LOCAL_REPO_NAME)
        @mock git_clone(url_with_auth, local_path)

        @mock cd(local_path) do
            master_branch = @mock git_get_master_branch(options.master_branch)
            @mock git_checkout(master_branch)
        end

        # Get all the compat dependencies from the local Project.toml file
        project_file = @mock joinpath(local_path, subdir, "Project.toml")
        deps = get_project_deps(project_file; include_jll=options.include_jll)

        return deps
    end
end

function get_project_deps(project_file::AbstractString; include_jll::Bool=false)
    project_deps = Set{DepInfo}()
    project = TOML.parsefile(project_file)

    if haskey(project, "deps")
        deps = project["deps"]
        add_compat_section!(project)
        compat = project["compat"]

        for dep in deps
            name = dep[1]
            uuid = UUIDs.UUID(dep[2])

            # Ignore STDLIB packages and JLL ones if flag set
            if !Pkg.Types.is_stdlib(uuid) &&
                (!endswith(lowercase(strip(name)), "_jll") || include_jll)
                package = Package(name, uuid)
                compat_entry = DepInfo(package)
                dep_entry = convert(String, strip(get(compat, name, "")))

                if !isempty(dep_entry)
                    compat_entry.version_spec = semver_spec(dep_entry)
                    compat_entry.version_verbatim = dep_entry
                end

                push!(project_deps, compat_entry)
            end
        end
    end

    return project_deps
end

function clone_all_registries(f::Function, registry_list::Vector{Pkg.RegistrySpec})
    registry_temp_dirs = Vector{String}()

    for registry in registry_list
        tmp_dir = @mock mktempdir(; cleanup=true)
        local_registry_path = joinpath(tmp_dir, registry.name)
        push!(registry_temp_dirs, local_registry_path)
        @mock git_clone(registry.url, local_registry_path)
    end

    f(registry_temp_dirs)

    for tmp_dir in registry_temp_dirs
        @mock rm(tmp_dir; force=true, recursive=true)
    end
end

function get_latest_version!(deps, registries; options::Options)
    for registry in registries
        registry_toml_path = joinpath(registry, "Registry.toml")
        registry_toml = TOML.parsefile(registry_toml_path)
        packages = registry_toml["packages"]

        for dep in deps
            uuid = string(dep.package.uuid)

            if uuid in keys(packages)
                versions_toml_path = joinpath(
                    registry, packages[uuid]["path"], "Versions.toml"
                )
                versions_dict = TOML.parsefile(versions_toml_path)
                include_yanked = options.include_yanked
                exclude_yanked = !include_yanked
                if exclude_yanked
                    filter!(pair -> !get(pair[2], "yanked", false), versions_dict)
                end
                versions = VersionNumber.(collect(keys(versions_dict)))

                max_version = maximum(versions)
                dep.latest_version = _max(dep.latest_version, max_version)
            end
        end
    end
end

function get_existing_registries!(deps::Set{DepInfo}, depot::String; options::Options)
    registries = readdir(joinpath(depot, "registries"); join=true)
    get_latest_version!(deps, registries; options)

    return deps
end

function get_latest_version_from_registries!(
    deps::Set{DepInfo}, registry_list::Vector{Pkg.RegistrySpec}; options::Options
)
    @mock clone_all_registries(registry_list) do registry_temp_dirs
        get_latest_version!(deps, registry_temp_dirs; options)
    end

    return deps
end


text:
export gpu,cpu, NodalCorr

mutable struct NodalCorr <: NoiseData
    n::Int64
    id::Array{String,1}                 # id
    name::Array{String,1}               # name
    loc::Array{<:InstrumentPosition,1}  # loc
    fs::Float64                         # fs
    src::String                         # src
    t::Array{Array{Int64,2},1}          # time
    corr::AbstractArray{Float32, 2}     # actual data

    function NodalCorr(
        n        ::Int64,
        id       ::Array{String,1},
        name     ::Array{String,1},
        loc      ::Array{<:InstrumentPosition,1},
        fs       ::Float64,
        src      ::String,
        t        ::Array{Array{Int64,2},1} ,
        corr     ::AbstractArray{<:AbstractFloat,2},
        )

        return new(n, id, name, loc, fs, src, t, corr)
    end
end

NodalCorr(;
          n        ::Int64                       = 0,
          id       ::Array{String,1}             = Array{String}(undef,0),
          name     ::Array{String,1}             = Array{String}(undef,0),
          loc      ::Array{InstrumentPosition,1} = Array{InstrumentPosition}(undef,0),
          fs       ::Float64                     = zero(Float64),
          src      ::String                      = "",
          t        ::Array{Array{Int64,2},1}     = Array{Array{Int64,2},1}(undef,0),
          corr     ::AbstractArray{<:AbstractFloat,2} = Array{Float32,2}(undef, 0, 2),
          ) = CorrData(n, id, name, loc, fs, src, t, corr)

# functions for adapting NodalData to GPU
cpu(m::NodalData) = adapt(Array,m)
gpu(x::NodalData) = use_cuda[] ? cu(x) : x

# function for adapting NodalCorr to GPU
cpu(m::NodalCorr) = adapt(Array,m)
gpu(x::NodalCorr) = use_cuda[] ? cu(x) : x

Adapt.adapt_structure(to, N::NodalData) = NodalData(
    N.n,
    N.ox,
    N.oy,
    N.oz,
    N.info,
    N.id,
    N.name,
    N.loc,
    N.fs,
    N.gain,
    N.resp,
    N.units,
    N.src,
    N.misc,
    N.notes,
    N.t,
    adapt(to,N.data),
)

Adapt.adapt_structure(to, N::NodalCorr) = NodalCorr(
    N.n,
    N.id,
    N.name,
    N.loc,
    N.fs,
    N.gain,
    N.resp,
    N.units,
    N.src,
    N.t,
    adapt(to,N.corr),
)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 59afc220-9195-11eb-2384-29ae6f4aa263
md"""
### REPL 101

- `]` enters the package manager
- `;` enters the terminal
- `?` gives you the documentation of a function
"""

# ╔═╡ e5c840d8-8978-11eb-37fa-77fb05a3d821
md"""
# Multi-dimensional `Array`s
_continuation_ ([read more](https://docs.julialang.org/en/v1/manual/arrays/))

There are a bunch of basic functions one has to know to effectively work with arrays.

##### Constructing arrays

```julia
A = zeros(T, sizes...)
B = ones(T, sizes...)
C = rand(T, sizes...)
D = copy(A)
E = reshape(A, new_sizes...)
```

##### Working with arrays
```julia
eltype(A)
length(A)
size(A)
eachindex(A) # iterator for visiting each position in A
```

##### Concatenating arrays
Arrays can be concatenated with the `;` syntax
```julia
A = rand(3)
B = zeros(3)
C = [A; B]
```

or with the `cat` functions
```julia
C = cat(A,B,dims=1) # What would happens if `dims=2`
```

### Linear indexing
When exactly one index `i` is provided, that index no longer represents a location in a particular dimension of the array but the `i`th element using the **column-major order**
```julia
A = [2 6; 4 7; 3 1]

A[5] == vec(A)[5] == 7
```

- Question: How to understand the matrix-creation notation (for example with `A`) from the concatenation operator `;`?

### Array `@views`

**Slicing operations like x[1:2] create a copy by default in Julia**

A “view” is a data structure that acts like an array, but the underlying data is actually part of another array (reference!).

```julia
A = [1 2; 3 4]

b = view(A, :, 1)
b = @views A[:,1]

b[2] = 99
A
```

Note: Sometimes it's not faster to use `view`s!
"""

# ╔═╡ 8636e926-8b38-11eb-3f98-f58135f3d02e
md"""
### Broadcasting and loop fusions
It is common to have "vectorized" versions of functions, which simply `map` a function `f(x)` to each element of an array `a`. Unlike some languages, in `Julia` this is **not required** for performance: it's just a convenient for-loop.

This is achieved through the `dot` syntax
```julia
a = [0 π/2; -π/2 π/6]
sin.(a)
```

Due to historical (and parsing) reasons, the `dot` syntax for infix operators in on the left
```julia
a = [1.0, 2.0, 3.0]
a .^ 3 # NOT the power of a vector
```

On complicated expression with several `dot` calls, the operation is fused together (there will be a single loop)
```julia
b = 2 .* a.^2 .+ sin.(a)
```

The `@.` macro can be used to convert all function / operator calls in an expression to a `dot`-call
```julia
b = @. 2 * a^2 + sin(a)
```

and a `$` inserted to bypass the dot
```julia
@. sqrt(abs($sort(x))) # equivalent to sqrt.(abs.(sort(x)))
```

Naturally calls like
```julia
sin.(sort(cos.(X)))
```
can't be completely fused together.

**Singleton (size=1) and missing dimensions are expanded to match the extents of the other arguments by virtually repeating the value.**
```julia
a = rand(2,2)
b = zeros(2,2,3)
a .+ b
```

Dot calls are just syntatic sugar for `broadcast(f, As...)` so you can extend _broadcasting_ for custom types.
"""

# ╔═╡ 540516aa-8b38-11eb-26f8-d31b72022689
md"""
#### Exercise:
  - Convince yourself that the loop is indeed fused by `@time`ing a complex `dot`ed expression vs the expression terms seperately computed. Run the code once beforehand to avoid timing the compilation time!
"""

# ╔═╡ 53f03512-8b38-11eb-031b-55d81ef5aeac
md"""
##### Solution:

```julia
a = rand(100_000)

f(a) = @. a ^ 2 + sin(2a) / (a ^ 3)

function g(a)
	b = a .^ 2
	c = sin.(2a)
	d = a .^ 3
	e = c ./ d
	f = b + e
	return f
end

# Make sure they are doing the same thing
@assert f(a) == g(a)

@time f(a);
@time g(a);
```
"""

# ╔═╡ 533f5a58-8b38-11eb-1c75-61000cf8394f
md"""

# `using LinearAlgebra`

"High-level" mathematical functions to operate on (multi-)dimensional `Array`s.

Assume that when you call methods from [`LinearAlgebra`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) that libraries like OpenBLAS (default) will be called for standard types such as `Float64`.

If your life depends on OpenBLAS-like operations its worth to check other implementations such as Intel's MKL or even pure Julia's like `Octavian.jl` (especially if your matrices are more interesting than dense ones).

### `LinearAlgebra` vs `GenericLinearAlgebra`

Generic programming allied with multiple dispatch allows one to share types with generic algorithms. An example of this is `GenericLinearAlgebra`, which implements some of `LinearAlgebra` functions in pure Julia.
"""

# ╔═╡ e6cc166a-8ed2-11eb-23be-2d657eb00795
md"""
#### Exercise: Quantum Ising

- Consider a finite ``1``D spin-chain Hamiltonian with ``N`` sites coupled to a magnetic field `` H_N = -\sum^N_i \sigma^z_i \otimes \sigma^z_{i+1} - h \sum^N_i \sigma^x_i`` with ``\sigma^z_i \otimes \sigma^z_{i+1} := ... \otimes \mathbb{1}_{i-1} \otimes \sigma^z_i \otimes \sigma^z_{i+1} \otimes \mathbb{1}_{i+1}  \otimes ...``

  - Define the identity ``\mathbb{1}_2`` and Pauli matrices.
  - Take the Kronecker products using the function `kron` at all sites `i` and sum them.
  - Diagonalise the Hamiltonian `using LinearAlgebra`'s `eigen` function
    - Which states (columns of the eigenvectors matrix) have the lowest energy for ``h=0`` and ``h \gg 1``?
"""

# ╔═╡ 247b2cf8-8ed3-11eb-19ee-fd0f17cec519
#= md"""
##### Solution:

```julia
function H(;N,h)
	id = [1 0; 0 1]
	σˣ = [0 1; 1 0]
    σᶻ = [1 0; 0 -1]

	# final Hamiltonian
	hamiltonian = zeros(Float64, 2^N, 2^N)

	# kinetic term: sum only until (N-1) because there's no site to the right of N
	for i in 1:(N-1)

		# For each site contribution, σᶻ will be at position k and k+1
		mat(k) = k == i || k == (i + 1) ? σᶻ : id

		# kronecker products
		out = mat(1)
		for j in 2:N
			out = kron(out, mat(j))
		end

		# Add contribution to the Hamiltonian
		hamiltonian -= out
	end

	# magnetic term: can take sum everywhere
	for i in 1:N

		# For each site contribution, σᶜ will be at position k
		mat(k) = k == i ? σˣ : id

		# kronecker products
		out = mat(1)
		for j in 2:N
			out = kron(out, mat(j))
		end

		# Add contribution to the Hamiltonian
		hamiltonian -= h .* out
	end
	return hamiltonian
end
```

With a function that can generate the Hamiltonian generally, it's easy to check for the eigenstates.

Since the eigenvalues are sorted by ascending value, we just need to look at the first terms

```julia
using LinearAlgebra

h1 = H(N=8, h=0.0);
h2 = H(N=8, h=100);

vals1, vecs1 = eigen(h1);
vals2, vecs2 = eigen(h2);

@show vecs1[:,1]; # Vectors all pointing in the same direction
@show vecs2[:,1]; # Vectors aligned with σˣ (equal super-position of eigenstates)
```
""" =#

# ╔═╡ 9ba8170a-7b50-11eb-0aac-adc493e1f386
md"""
# A first taste of functional programming

**Functional programming** a programming paradigm where programs are constructed by applying and composing functions.

Function definitions are trees of expressions that map values to other values.
Unlike imperative programming, in which a sequence of imperative statements which updates the running state of the program.

There are 2 main conceptes from FP which have clear benefits for scientific code
- Higher-order functions
- Pure functions

## Higher order functions

In Julia the functions are already first-class citizens (can pass them as arguments, return them and assign them to variable names).

A higher-order function is a function that
- takes other functions as arguments
and/or
- returns a function as result

### Function composition and piping
Many times in scientific code it's required that functions are chained together.

This can be quite of eye sore
```julia
sqrt(abs(sum([1,2,3])))
```

Function composition `\circ`
```julia
(sqrt ∘ abs ∘ sum)([1,2,3])
```

and piping
```julia
sum([1,2,3]) |> abs |> sqrt
```
alleviate this problem for _unary_ (single-argument) functions.
"""

# ╔═╡ c9862c2a-897e-11eb-3592-ada39cc0637b
md"""
#### Exercise: Becoming a pastry chef
(Re)create some syntatic sugar such as
- a function composition operator `∘` for _unary_ functions
- a reverse pipe `<|`

Note: infix operators such as `∘` need to be wrapped around `()` in method definitions for parsing reasons.
"""

# ╔═╡ 44b9eb9c-8980-11eb-17b8-1937532d3a28
md"""
##### Solution:
```julia
(∘)(f::Function, g::Function) = (x...) -> f(g(x...))

(<|)(f::Function, args...) = f(args...)
```
"""

# ╔═╡ af38a46c-8981-11eb-26de-35f9e487423c
md"""
One of the most glaring differences between functional and imperative programming is that functional avoids _side effects_, which are needed in imperative style to control the state of the program.

Here are examples of higher-order functions that take us closer to the functional heaven.

### map

`map(f, [a1,a2,...]) = [f(a1), f(a2), ...]`

`map(g, [a1,a2,...], [b1,b2,...]) = [g(a1,b1), g(a2,b2), ...]`

Maps a function over the elements of a container and collect the results

```julia
# data
v1 = [1, 2, 3]
v2 = [10, 20, 30]

# imperative
out = zero(v1)
for i in eachindex(v1) # equivalent to 1:length(v1)
	out[i] = v1[i] + v2[i]
end

# functional
out = map(+, v1, v2)
```

##### `map` vs `broadcast`
  - Broadcast only handles containers with th "shapes" `M×N×...` (i.e., a size and dimensionality) while map is more general (unknown length iterator)
  - Map requires all arguments to have the same length (and hence cannot combine arrays and scalars)

### foreach
`foreach(f, [a1,a2,...]) = f(a1); f(a2); ...; nothing`

Map a function over the elements of a container but without collecting the results

```julia
# data
v = [1, 2, 3]

# imperative
for i in eachindex(v)
	println(v[i])
end

# functional
foreach(println, v)
```

### foldl and foldr
`foldl(op, [a1, a2, a3, ...]) = op(op(op(a1, a2), a3), ...)`

`foldr(op, [a1, a2, a3, ...]) = op(a1, op(a2, op(a3, op(...))))`

Folds are ubiquitous in functional and scientific programming. They are a class of functions that process some data structure and return a value. Basically a left- and right- associative reduce.

```julia
# data
v = [1, 2, 3, 4]

# binary function
op = *

# imperative
out = v[1]
for i in 2:length(v)
	out = op(out, v[i])
end

# functional
out = foldl(op, v)
```

Question: why is there a distinction between left and right folds?

For most cases these binary functions are not associative!

We implicitely have left- and right- associativeness hard coded in our brain. What's `0 - 1 - 2 - 3 - 4` ?
- left-associative `((((0 - 1) - 2) - 3) - 4) == -10`
- right-associative `(0 - (1 - (2 - (3 - 4)))) == 2`

An initial value cal also be passed as a _keyword_
```julia
foldl(=>, 1:4) == ((1 => 2) => 3) => 4
foldl(=>, 1:4; init=0) == (((0 => 1) => 2) => 3) => 4
```
"""

# ╔═╡ afc5f13c-8ee2-11eb-2dfe-0dcbc74c120e
md"""
#### Exercise: folding left and right
- Use a folding operator to find the `min`imum element in a container
  - Define your own `min(a,b)` function or use Julia's
"""

# ╔═╡ bded83c2-8ee2-11eb-3b00-a7681ba476dc
md"""
##### Solution:

```julia
my_min(x,y) = x < y ? x : y
min_fold(c) = foldl(my_min, c)
```
"""

# ╔═╡ 0d6a34a4-8ee3-11eb-08ae-2763ca230fd5
md"""
### mapreduce
`mapreduce(f, op, [a1, a2, a3, ...]) = op(op(op(f(a1), f(a2)), f(a3)), ...)`

Applies a function `f` to each element of the container and `reduces` the result using a binary function `op`.

```julia
# data
v = [1, 2, 3, 4]

# operations
f = sin
op = *

# imperative
out = 1.0
for i in eachindex(v)
	out = op(out, f(v[i]))
end

# functional
out = mapreduce(f, op, v; init=1.0)
```

Unlike a `fold`, a `reduce` operation **assumes associativity** with the binary operations. If this cannot be guaranteed, a `mapfoldl` or `mapfoldr` can be used instead.

Question: Where could a `reduce` operation be preferred (computationally!) over a fold?

It doesn't end here: check also `filter`, `reduce` and possible multiple chains with `Transducers.jl`.
"""

# ╔═╡ 59e7aaf0-8cac-11eb-0ea8-3fb87d36ad90
md"""
#### Exercise: Wilson chains
- Consider a finite ``1``D Wilson chain with ``N`` sites `` H_N = -\sum^N_i \alpha^{-i} \sigma^z_i \otimes \sigma^z_{i+1}`` for ``\alpha \ge 1``
  - Build the Hamiltonian in a functional way. Suggestion:
    - Use `mapfold` or `fold` to take the Kronecker products

- Numerically renormalize the problem by iterating
  - Diagonalise ``H_N =: U_N D_N U^\dagger_N`` `using LinearAlgebra`'s `eigen` function
  - Truncate ``H_N \rightarrow \tilde{H}_N`` by taking only the half lowest eigenvalues
  - Coupling ``\tilde{H}_N`` to a next site in the chain ``H_{N+1} = \tilde{D}_N \otimes \mathbb{1}_2 - \alpha^{-(N+1)} \tilde{\sigma}^z_N \otimes \sigma^z_{N+1}`` where ``\tilde{\sigma}^z_N = \tilde{U}_N \sigma^z_N \tilde{U}_N^†``

Tip: A diagonal matrix can be created out of the vector of eigenvalues with the `diagm` function and the adjoint of a matrix is achieved with the `adjoint` function or `'` operator after the matrix.
"""

# ╔═╡ 45c45594-8ca4-11eb-1db9-6187459bf306
#= md"""
##### Solution:

```julia
using LinearAlgebra

const id = [1 0; 0 1]
const σᶻ = [1 0; 0 -1]

# Wilson chain Hamiltonian
function H(N, α)
	T(i) = -α^(-i) * mapreduce(j -> j == i || j == i + 1 ? σᶻ : id, kron, 1:N)
	return mapreduce(T, +, 1:(N-1))
end

# Couple H to the chain site `n`
function add_chain(H, n, α)
	# truncates half of the highest eigenvalues
	λ̃, Ũ = let
		λ, U = eigen(H)
		N₂ = length(λ) ÷ 2
		λ[1:N₂], U[1:N₂, 1:N₂]
	end

	σᶻₙ = let
		n = ceil(Int, log2(length(λ̃))) # number of sites in H̃
		mapreduce(i -> i == n ? σᶻ : id, kron, 1:n)
	end

	# Rotate σᶻ → σ̃ᶻₙ and add chain
	return kron(diagm(λ̃), id) - α^(-(n+1)) * kron(Ũ * σᶻₙ * Ũ', σᶻ)
end

# Numerical renormalization group
rg(;N,α,steps=0) = foldl((h, i) -> add_chain(h, N+i, α), 1:steps; init=H(N, α))
```
""" =#

# ╔═╡ 4bee5324-8eeb-11eb-326b-efba53215f34
md"""
### Pure functions
_Pure thoughts_

Another great tool of functional programming we should steal for scientifc programming is the concept of pure functions. These functions are very close to mathematical functions.

- The function return values are identical for identical arguments
- The application of the function has no side effects
  - No mutation of non-local variables or mutable reference arguments

Examples of _indecent thoughts_:

Mutation of non local variables ⚰⚰⚰
```julia
f() = x[1] += 1

x = [1]
@show x

f()
@show x
```

Different return values with identical arguments 🪦🪦🪦
```julia
f() = x

x = 1
@show f()

x = 2
@show f()
```

##### Why do I care?
- Becase you're not a sociopath
- If the result of a pure expression is not used, it can be removed without affecting other expressions.
- Pure functions have no side-effects and you can intelectually refeer to them as being _referentially transparent_, just like mathematical functions
- Since they only depend on their arguments, different function calls can't interfer with each other (great for parallel programming!)
- No side effects means your compiler can, theoretically, safely apply
- Unit tests are valid and can be injected anywhere
"""

# ╔═╡ 28f63950-8eef-11eb-0481-55fad44619d5
md"""
### Meta-discussion: mutable vs immutable/pure algorithms

Immutability doesn't really exist: immutability implies time-independence... and there's nothing really stopping time (at least until the heat-death of the universe).

The very process of storing information (that is ordering bits) requires mutation.
But we can achieve immutability at least syntatically.

##### Tips to scientific code right, denying mutation and and promoting good hygiene

- Use `let` blocks to reduce global scope pollution
    - Global variables are **very** prone to be mutated since they don't have to be passed as an argument explicitely

- Pure thoughts: decompose programs into (pure) functions:
    - Break software into chunks to fit into the most limited memory: human memory.

- Give functions and variables meaningful names
    - Use `Pluto` notebooks to prototype

- Use tuples / structs to avoid repetition
    - `a1 = 1, a2 = 2` becomes `as = (1, 2)`

- Be defensive
    - Add `@assert`s to ensure validity of your inputs / results
    - Generate unit tests for your functions: these are as important as the problem you are ultimately solving

- Be smart by not oversmarting yourself:
    - avoid _premature optimisation_: write clear and concise code and only think about optimisations after unit testing
    - avoid _premature pessimisation_: take a chill pill and sketch on paper the data structures / algorithm design before writing any code
    - require of your code the same standards you require others' calculations / experiments / general care in life

Read more on good Scientific Practises
- [1](https://swcarpentry.github.io/good-enough-practices-in-scientific-computing/)
- [2](https://arxiv.org/pdf/1210.0530v3.pdf)
- [3](https://blog.higher-order.com/blog/2009/04/27/a-critique-of-impure-reason/)
"""

# ╔═╡ 8ac22e72-8978-11eb-1674-e1b95403e215
md"""
# Performance
We will step out of the beautiful pure world of functional programming and dive deep into the messy world of algorithm optimisation.

While some topics will be related to avoiding common pitfalls when writing in Julia, other universal topics on algorithm optimisation will also be covered.

## Profiling
In order to know where and what to optimise we need tools to diagnose time spent, memory allocations and possibly how machine code is being generated.

- For quick-and-dirty diagnostics (time and memory tracking) we can use the `@time` macro behind function calls
- For an accurate measure of the latter preferer `using BenchmarkTools` and the `@btime` macro (this will give you the lower bound of `@benchmark`, which is what you want to measure!)
- For more serious profiling tools, consider reading the [manual section](https://docs.julialang.org/en/v1/manual/profile/)

The code can be inspected at several stages with the macros
- The AST after parsing: `@macroexpand`
- The AST after type inference and some optimizations: `@code_typed` (prefer `@code_warntype`)
- The LLVM and assembly: `@code_llvm`, `@code_native`

Read more on [introspection](https://docs.julialang.org/en/v1/devdocs/reflection/)
"""

# ╔═╡ 8d9f2162-8f07-11eb-1c1d-b994b88811e1
md"""
#### Exercise: Not all "created" equal
Benchmark the different ways of mapping a function to a container and realise that computationally not all operations are equivalent
- `map`
- `broadcast`
- list comprehension
- explicit `for`-loop
"""

# ╔═╡ bff46124-8f08-11eb-1140-d78145e2fe05
#= md"""

```julia
using BenchmarkTools

fmap(x) = map(x -> 2x, x)

fdot(x) = 2 .* x

fcomprehension(x) = [2x for x in x]

function floop(x)
    y = similar(x) # prefer similar to zero as it only has to allocate (and not zero)
    for i in eachindex(x)
        y[i] = 2*x[i]
    end
    return y
end

x = rand(10_000)

# The $ interpolates values into the expression avoiding the global-variable problem

@btime fmap($x);
# > 7.977 μs (2 allocations: 78.20 KiB)
@btime fcomprehension($x);
# > 7.890 μs (2 allocations: 78.20 KiB)
@btime fdot($x);
# > 8.238 μs (2 allocations: 78.20 KiB)
@btime floop($x);
# > 9.495 μs (2 allocations: 78.20 KiB)
```
""" =#

# ╔═╡ fdfddce2-919d-11eb-1b48-5d31c8d4b250
md"""
### Global (`isa Any`) variables: Electric Boogaloo
_It's not over yet_

[1st performance tip](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables) on the Julia's official documentation: **Avoid global variables**

A global variable might have its value, and therefore its type, change at any point. This makes it difficult for the compiler to optimize code using global variables.
- Variables should be local
- Or passed as arguments to functions (this way the code will be specialized for the input types)

```julia
x = rand(100_000);

function sum_global()
	s = 0.0
	for i in x
	   s += i
	end
	return s
end

function sum_arg(x)
	s = 0.0
	for i in x
	   s += i
	end
	return s
end

@time sum_global()
# 0.016230 seconds (399.49 k allocations: 7.622 MiB)

@time sum_arg(x)
# 0.000123 seconds (1 allocation: 16 bytes)
```

This unfortunate dichotomy between the functions is due to one thing: Julia can't **specialize** the function `sum_global` since it cannot guarantee the type of `x`.
While for `sum_arg`, it will specialize it for the type of `x` that is encountered (a complex `x` would trigger compilation).

**Type inference did not require the type to be annotated!**

`@code_warntype` can detect these `isa Any` variables or the also related type instabilities!
"""

# ╔═╡ 8a97bcb8-91a3-11eb-37b4-4f9004e01e0b
md"""
#### Exercise: Purging type instabilities / "untyped" containers
- Remember that an array with abstract type (e.g. `Any` or `Real`) will end up being an array of pointers and can't be operated on efficiently. `abstract` types can also prevent Julia from triggering the appropriate function specialization. Compare the `@code_native` outputs from a simple function operating on the fields of
```julia
mutable struct MyAmbiguousType
   a::AbstractFloat # less ambiguous than `Any` but still ambiguous
end

mutable struct MyType{T<:AbstractFloat}
   a::T
end
```

- Can you spot the type instability? Fix it and consequently `@btime` and `@code_native` the new function vs the old one.
```julia
function f()
	x = 1
	for i = 1:10
		x = x/2
	end
	return x
end
```
"""

# ╔═╡ 0ae35bf0-91a6-11eb-1a75-235aaf78d346
#= md"""
##### Solution:

```julia
add1(m::MyType) = m.a + 1

@code_native add1(MyType{Float64}(0.3))
@code_native add1(MyType{AbstractFloat}(0.3))
```
""" =#

# ╔═╡ 700b1240-9197-11eb-0cfd-1bb3d7393084
md"""
## Understanding timecales
There's a rather famous table comparing computer with human timecales

| Action (3GHz) | Average latency | Human time |
|:---------- | ---------- |:------------:|
| 1 clock-cyle   | 0.3 ns | 1 s |
| L1 cache    | 0.9 ns  | 3 s |
| L2 cache    | 2.8 ns  | 9 s |
| L3 cache    | 12.9 ns  | 43 s |
| RAM    | 70-100 ns  | 3.5 to 5.5 min |
| SSD/IO | 7-150 μs | 2h to 2 days |
| Reboot | 30-60s | 1000-2000 years |

In a single cycle a photon can only travel <10 cms to the RAM better not be too far away (imagine an electron).

The slowness of RAM can be mitigated by transfering data into the caches. **BUT** when the CPU requests data from the RAM, it is checked if it's in the cache. If not, there will be a _cache miss_ and the program will stall until the data is fetched from the RAM.

Naturally, to reduce these cache misses, consider
- _Temporal locality_: If you need to access a piece of memory multiple times, do it close in time
- _Spatial locality_: Access memory which is close to each other (since the CPU fetches chunks of data at a time)

As a corollary
- Use little memory
- Access data sequentially (since the CPU can prefetch data that you may need)

It doesn't end here: read more on [alignment issues](https://biojulia.net/post/hardware/#alignment) [(show-in-class)](https://juliasimd.github.io/LoopVectorization.jl/stable/examples/matrix_multiplication/)

### Allocations
Memory allocation can be a significant bottleneck in critical operations.
Dynamic languages such as Julia usually employ a _garbage collector_ to allocate and deallocate objects in the RAM for us
```julia
a = [1,2,3] # allocation
a = nothing # the previous value of `a` is now garbage (since in this case no other variables are pointing at it) and shall be collected automatically
```
- Allocation and deallocation create _overhead_ which can be very costly
- More allocations results in more memory usage which results in more cache misses

The 3 most encountered problems where this can be fixed are

- When updating some value _inplace_
```julia
N = 1000
a = rand(10,10);

function f1(a, N)
	x = zero(a)
	for i in 1:N
		x += i * a # remember that updating operators such as `+=` reassign `x`
	end
	x
end

@btime f1($a, $N)
# 436.185 μs (2001 allocations: 1.71 MiB)
```

In this case use the **inplace (broadcasted) assignment**
```julia
function f2(a, N)
	x = zero(a)
	for i in 1:N
		x .+= i * a # "add element-wise to `x`
	end
	x
end

@btime f2($a, $N)
# 354.372 μs (1001 allocations: 875.88 KiB)
```

- When running a compution whose output memory can be recycled (see the exercise!)

- When taking slices (see the exercise!)

#### Single instruction, multiple data: `@simd`
CPUs operate on data present in registers inside the CPU, which are meant to hold small fixed size slots, like floats (see the `r`s in `@code_native`).
Since this is a major bottleneck, modern CPUs have bigger registries (instead of 64-bit, 256+), which allow a **S**ingle **I**nstruction operate on **M**ultiple **D**ata.

Show the difference of
```julia
using StaticArrays
a = @SVector Int32[1,2,3,4,5,6,7,8]
code_native(+, (typeof(a), typeof(a)), debuginfo=:none)

a = @SVector Int64[1,2,3,4,5,6,7,8]
code_native(+, (typeof(a), typeof(a)), debuginfo=:none)

a = @SVector Int64[1,2,3,4]
code_native(+, (typeof(a), typeof(a)), debuginfo=:none)
```

- SIMD needs uninterrupted iteration of fixed length
- Bound-checking causes branching so deactivate it with `@inbounds for i in ...` (This is also may be desirable for critical non-`@simd` loops)
  - Avoid branching at all costs. Even with _branch prediction_, a misprediction (common for random braches) will cost several CPU cycles.
- SIMD needs associative operations (since the loop will be reordered)
Since float addition is **not** associative, automatic `@simd` is not "automatically" on for, e.g., float addition
```julia
@show 0.1 + (0.2 + 0.3)
@show (0.1 + 0.2) + 0.3
```

Actually, IEEE 754 float arithmetic is tricky! Consider a number that is undefined or unrepresentable, a `NaN`
```julia
@show NaN == NaN # not even reflexive!
```

[Reference](https://biojulia.net/post/hardware/#simd)
"""

# ╔═╡ a3346d6a-919a-11eb-13e9-59293be49654
md"""
#### Exercise: The 3 ecology Rs: RECYCLE, REUSE and REDUCE

- Avoid extra allocations by _recycling_
  - Rewrite `loopinc` using an in-place version of `xinc` –> `xinc!`
  - Compare the performance of the new `loop` function
  - Note: The point is not to optimise away the `ret[2]` part so keep it!
```julia
xinc(x) = [x, x+1, x+2]
function loop()
	y = 0
	for i = 1:10^7
		ret = xinc(i)
		y += ret[2]
	end
	return y
end
```

- Avoid extra allocations by _reusing_
  - Create a `(50, 100_000)` `rand`om matrix `A` and a vector `x` with size `(100_000,)`
  - Compare the performance of ``\sum_i \sum^{80000}_{j=1}A_{ij}x_j`` using `slice`s and `@view`s
  - Note that performance from `@view` is just a [rule of thumb](https://docs.julialang.org/en/v1/manual/performance-tips/#Copying-data-is-not-always-bad).

- Optimise away by _reducing_ (and everything else)
```julia
function work()
	A = zeros(N,N)
	for i in 1:N
        for j in 1:N
            val = mod(v[i],256);
            A[i,j] = B[i,j]*(sin(val)*sin(val)-cos(val)*cos(val))
        end
	end
	return A
end
```

given the parameters
```julia
N = 4_000
B = [float(i-j) for i in 1:N, j in 1:N]
v = [i for i in 1:N]
```

Pro-tip: Compare different implementations of `work!` using the `≈` (`\approx`) operator, since the `==` may be too strict given the shenanigans we encountered with float arithmetics.
"""

# ╔═╡ fd81c0c0-91fd-11eb-1ddb-1546f52626af
#= md"""
##### Solution
This excellent exercise was adapted from [here](https://github.com/crstnbr/JuliaCologne21/blob/master/Day2/exercise_solutions/solution_optimization2.ipynb)

```julia
using Test
using BenchmarkTools

# Parameters
N = 1_000
B = [float(i-j) for i in 1:N, j in 1:N]
v = [i for i in 1:N]

# Base
function work()
	A = zeros(N,N)
	for i in 1:N
        for j in 1:N
            val = mod(v[i],256);
            A[i,j] = B[i,j]*(sin(val)*sin(val)-cos(val)*cos(val))
        end
	end
	return A
end

@btime work()
# 670.186 ms (11469981 allocations: 197.94 MiB)
```

- Purge globals!
```julia
@code_warntype work()

function work1(B, v, N)
	A = zeros(N,N)
	for i in 1:N
        for j in 1:N
            val = mod(v[i],256);
            A[i,j] = B[i,j]*(sin(val)*sin(val)-cos(val)*cos(val))
        end
	end
	return A
end

@code_warntype work(B, v, N)

@test work() ≈ work1(B, v, N)
@btime work1($B, $v, $N)
# 88.320 ms (2 allocations: 7.63 MiB)
```

- Analytic optimisations (these are the best)
```julia
@testset "My trig identities" begin
	x = rand()
	@test 1-2*cos(x)*cos(x) ≈ sin(x)*sin(x)-cos(x)*cos(x)
	@test -cos(2*x) ≈ sin(x)*sin(x)-cos(x)*cos(x)
end

function work2(B, v, N)
    A = zeros(N,N)
    for i in 1:N
        for j in 1:N
            val = mod(v[i],256);
            A[i,j] = B[i,j]*(-cos(2*val));
        end
    end
	return A
end

@test work() ≈ work2(B, v, N)
@btime work2($B, $v, $N)
# 37.520 ms (2 allocations: 7.63 MiB)
```

- Pull-out `val` computation
```julia
function work3(B, v, N)
    A = zeros(N,N)
    for i in 1:N
		val = -cos(2*mod(v[i],256))
        for j in 1:N
            A[i,j] = B[i,j]*val;
        end
    end
	return A
end

function work4(B, v, N)
	val = [-cos(2*mod(x,256)) for x in v]

    A = zeros(N,N)
    for i in 1:N
        for j in 1:N
            A[i,j] = B[i,j]*val[i];
        end
    end
	return A
end

@test work() ≈ work3(B, v, N)
@test work() ≈ work4(B, v, N)

@btime work3($B, $v, $N)
# 20.755 ms (2 allocations: 7.63 MiB)

@btime work4($B, $v, $N)
# 20.262 ms (3 allocations: 7.64 MiB)
```

- Switch order of the loops: favour data locality!
```julia
function work5(B, v, N)
	val = [-cos(2*mod(x,256)) for x in v]

    A = zeros(N,N)
    for j in 1:N
        for i in 1:N
            A[i,j] = B[i,j]*val[i];
        end
    end
	return A
end

@test work() ≈ work5(B, v, N)

@btime work5($B, $v, $N)
# 4.947 ms (3 allocations: 7.64 MiB)
```

- `@inbounds` and `@simd`
```julia
function work6(B, v, N)
	val = [-cos(2*mod(x,256)) for x in v]

    A = zeros(N,N)
    for j in 1:N
        for i in 1:N
            @inbounds A[i,j] = B[i,j] * val[i];
        end
    end
	return A
end

@test work() ≈ work6(B, v, N)

@btime work6($B, $v, $N)
# 3.298 ms (3 allocations: 7.64 MiB)
```

- Broadcast it for beauty points
```julia
work7(B, v) = return B .* [-cos(2*mod(x,256)) for x in v]

@test work() ≈ work7(B, v)

@btime work7($B, $v, $N)
# 2.019 ms (3 allocations: 7.64 MiB)
```
""" =#

# ╔═╡ 4eb2573e-8998-11eb-2274-379a03bed49c
md"""
# Iteration Utilities
For more examples see [here](https://docs.julialang.org/en/v1/base/iterators/)

### zip
Run multiple iterators at the same time, until any of them is exhausted
```julia
a = [1,2,3]
b = (10,20,30)

for z in zip(a,b)
	println(z) # (1,10) ... (2, 20) ... (3, 30)
end
```

Question: What is the mathematical operation equivalent of zipping?

### enumerate
An iterator that yields `(i, x)` where `i` is a counter starting at `1`, and `x` is the `i`-th value from the given iterator
```julia
a = [10, 20, 30]

for (i, aᵢ) in enumerate(a)
	println("The $i-th entry of a is $aᵢ)
end
```
"""

# ╔═╡ b58a861c-7b50-11eb-286d-c5a5dd03429f
#= md"""
# Fixed points

A fixed point of a function is an element of the function's domain that is mapped to itself by the function.

$ x,\;f(x),\;f \circ f(x),\;f \circ f \circ f(x),\;... \rightarrow x^*$

Let's write a `Julia` function that finds the fixed point of some function `f`

```julia

# Identity operator
const id = x -> x

# Ideally... Why will it fail?
function ideal_fixed_point(f::Function)
    fix(x) = f ∘ fix(x)
    return fix
end

bad_g = fixed_point(id)
# @show bad_g(3) why will it fail??
```

```julia
# Need a stopping criterion
function fixed_point(f::Function)
	function fix(xᵢ)
		xᵢ₊₁ = f(xᵢ)
		if xᵢ₊₁ == xᵢ
			return xᵢ
		else
			return fix(xᵢ₊₁)
		end
	end
	return fix
end

good_g = fixed_point(id)
@show good_g(3);

# TODO: Write this function as a 1-liner!
fixed_point(f::Function) = nothing
```

```julia
# The equality condition is too strict for computations! Change it for an approximation (\approx)
approx_fixed_point(f::Function) = fix(xᵢ) = (xᵢ₊₁ = f(xᵢ)) ≈ xᵢ ? xᵢ : fix(xᵢ₊₁)

g = approx_fixed_point(id)
@show g(3);
```

```julia
# Add another stopping criteria to avoid getting stuck
approx_stop_fixed_point(f::Function; maxᵢ=10^4) = fix(xᵢ, i=1) = (xᵢ₊₁ = f(xᵢ)) ≈ xᵢ || i == maxᵢ ? xᵢ : fix(xᵢ₊₁, i+1)

# Kondo poor-man's scaling equations
kondo_rg_eqs(x; dk=-10^-3) = x .+ dk * [-2 * x[2]^2, -2 * x[1] * x[2]]
kondo_rg = approx_stop_fixed_point(kondo_rg_eqs)
```

```julia
using Plots

dj = 0.06

# ps = hcat([[[jz, j₊] ; kondo_rg([jz, j₊])] for j₊ in 0.0:dj:1.0 for jz in -dj-j₊:-dj:-1.0]...)

plot(ps[[1,3],:], ps[[2,4],:];
	xlims = (-1,dj),
	ylims=(-dj,1),
	framestyle = :origin,
	legend=false)
```

```julia
# Another example
logistic_map(r) = approx_stop_fixed_point(x -> r * x * (1-x))
```
""" =#

# ╔═╡ Cell order:
# ╟─59afc220-9195-11eb-2384-29ae6f4aa263
# ╟─e5c840d8-8978-11eb-37fa-77fb05a3d821
# ╟─8636e926-8b38-11eb-3f98-f58135f3d02e
# ╟─540516aa-8b38-11eb-26f8-d31b72022689
# ╟─53f03512-8b38-11eb-031b-55d81ef5aeac
# ╟─533f5a58-8b38-11eb-1c75-61000cf8394f
# ╟─e6cc166a-8ed2-11eb-23be-2d657eb00795
# ╟─247b2cf8-8ed3-11eb-19ee-fd0f17cec519
# ╟─9ba8170a-7b50-11eb-0aac-adc493e1f386
# ╟─c9862c2a-897e-11eb-3592-ada39cc0637b
# ╟─44b9eb9c-8980-11eb-17b8-1937532d3a28
# ╟─af38a46c-8981-11eb-26de-35f9e487423c
# ╟─afc5f13c-8ee2-11eb-2dfe-0dcbc74c120e
# ╟─bded83c2-8ee2-11eb-3b00-a7681ba476dc
# ╟─0d6a34a4-8ee3-11eb-08ae-2763ca230fd5
# ╟─59e7aaf0-8cac-11eb-0ea8-3fb87d36ad90
# ╟─45c45594-8ca4-11eb-1db9-6187459bf306
# ╟─4bee5324-8eeb-11eb-326b-efba53215f34
# ╟─28f63950-8eef-11eb-0481-55fad44619d5
# ╟─8ac22e72-8978-11eb-1674-e1b95403e215
# ╟─8d9f2162-8f07-11eb-1c1d-b994b88811e1
# ╟─bff46124-8f08-11eb-1140-d78145e2fe05
# ╟─fdfddce2-919d-11eb-1b48-5d31c8d4b250
# ╟─8a97bcb8-91a3-11eb-37b4-4f9004e01e0b
# ╟─0ae35bf0-91a6-11eb-1a75-235aaf78d346
# ╟─700b1240-9197-11eb-0cfd-1bb3d7393084
# ╟─a3346d6a-919a-11eb-13e9-59293be49654
# ╟─fd81c0c0-91fd-11eb-1ddb-1546f52626af
# ╟─4eb2573e-8998-11eb-2274-379a03bed49c
# ╟─b58a861c-7b50-11eb-286d-c5a5dd03429f

max_stars_repo_name: Fromeworld/JuliaSeminar

max_issues_repo_name: Fromeworld/JuliaSeminar

max_forks_repo_name: Fromeworld/JuliaSeminar


text:
module StaticLint
using SymbolServer
using CSTParser
using CSTParser: isidentifier
using CSTParser: EXPR, PUNCTUATION, IDENTIFIER, KEYWORD, OPERATOR
using CSTParser: Call, UnaryOpCall, BinaryOpCall, WhereOpCall, Import, Using, Export, TopLevel, ModuleH, BareModule, Quote, Quotenode, MacroName, MacroCall, Macro, x_Str, FileH, Parameters, FunctionDef
using CSTParser: setparent!, kindof, valof, typof, parentof

const noname = EXPR(CSTParser.NoHead, nothing, 0, 0, nothing, CSTParser.NoKind, false, nothing, nothing)
baremodule CoreTypes # Convenience
using ..SymbolServer
const DataType = SymbolServer.stdlibs["Core"].vals["DataType"]
const Function = SymbolServer.stdlibs["Core"].vals["Function"]
const Module = SymbolServer.stdlibs["Core"].vals["Module"]
const String = SymbolServer.stdlibs["Core"].vals["String"]
const Int = SymbolServer.stdlibs["Core"].vals["Int"]
const Float64 = SymbolServer.stdlibs["Core"].vals["Float64"]

end

include("bindings.jl")
include("scope.jl")

mutable struct Meta
    binding::Union{Nothing,Binding}
    scope::Union{Nothing,Scope}
    ref::Union{Nothing,Binding,SymbolServer.SymStore}
    error
end
Meta() = Meta(nothing, nothing, nothing, nothing)

function Base.show(io::IO, m::Meta)
    m.binding !== nothing && show(io, m.binding)
    m.ref !== nothing && printstyled(io, " * ", color = :red)
    m.scope !== nothing && printstyled(io, " new scope", color = :green)
    m.error !== nothing && printstyled(io, " lint ", color = :red)
end
hasmeta(x::EXPR) = x.meta isa Meta
hasbinding(m::Meta) = m.binding isa Binding
hasref(m::Meta) = m.ref !== nothing
hasscope(m::Meta) = m.scope isa Scope
scopeof(m::Meta) = m.scope
bindingof(m::Meta) = m.binding

mutable struct State{T}
    file::T
    targetfile::Union{Nothing,T}
    included_files::Vector{String}
    scope::Scope
    delayed::Bool
    urefs::Vector{EXPR}
    server
end

function (state::State)(x::EXPR)
    delayed = state.delayed # store states

    resolve_import(x, state)
    if typof(x) === Export # Allow delayed resolution
        state.delayed = true
    end
    mark_bindings!(x, state)
    add_binding(x, state)
    mark_globals(x, state)
    handle_macro(x, state)
    s0 = scopes(x, state)
    _resolve_ref(x, state)
    followinclude(x, state)

    traverse(x, state)

    # return to previous states
    state.scope != s0 && (state.scope = s0)
    state.delayed = delayed
    return state.scope
end

"""
    traverse(x, state)

Iterates across the child nodes of an EXPR in execution order (rather than
storage order) calling `state` on each node.
"""
function traverse(x::EXPR, state)
    if typof(x) === CSTParser.BinaryOpCall && (CSTParser.is_assignment(x) && !CSTParser.is_func_call(x.args[1]) || typof(x.args[2]) === CSTParser.Tokens.DECLARATION) && !(CSTParser.is_assignment(x) && typof(x.args[1]) === CSTParser.Curly)
        state(x.args[3])
        state(x.args[2])
        state(x.args[1])
    elseif typof(x) === CSTParser.WhereOpCall
        @inbounds for i = 3:length(x.args)
            state(x.args[i])
        end
        state(x.args[1])
        state(x.args[2])
    elseif typof(x) === CSTParser.Generator
        @inbounds for i = 2:length(x.args)
            state(x.args[i])
        end
        state(x.args[1])
    elseif typof(x) === CSTParser.Flatten && x.args !== nothing && length(x.args) === 1 && x.args[1].args !== nothing && length(x.args[1]) >= 3 && length(x.args[1].args[1]) >= 3
        for i = 3:length(x.args[1].args[1].args)
            state(x.args[1].args[1].args[i])
        end
        for i = 3:length(x.args[1].args)
            state(x.args[1].args[i])
        end
        state(x.args[1].args[1].args[1])
    elseif x.args !== nothing
        @inbounds for i in 1:length(x.args)
            state(x.args[i])
        end
    end
end

"""
    followinclude(x, state)

Checks whether the arguments of a call to `include` can be resolved to a path.
If successful it checks whether a file with that path is loaded on the server
or a file exists on the disc that can be loaded.
If this is successful it traverses the code associated with the loaded file.

"""
function followinclude(x, state::State)
    if typof(x) === Call && typof(x.args[1]) === IDENTIFIER && valof(x.args[1]) == "include"
        path = get_path(x, state)
        if isempty(path)
        elseif hasfile(state.server, path)
        elseif canloadfile(state.server, path)
            loadfile(state.server, path)
        elseif hasfile(state.server, joinpath(dirname(getpath(state.file)), path))
            path = joinpath(dirname(getpath(state.file)), path)
        elseif canloadfile(state.server, joinpath(dirname(getpath(state.file)), path))
            path = joinpath(dirname(getpath(state.file)), path)
            loadfile(state.server, path)
        elseif !isempty((basepath = _is_in_basedir(getpath(state.file)); basepath))
            # Special handling for include method used within Base
            path = joinpath(basepath, path)
            if hasfile(state.server, path)
                # skip
            elseif canloadfile(state.server, path)
                loadfile(state.server, path)
            else
                path = ""
            end
        else
            path = ""
        end
        if !isempty(path)
            if path in state.included_files
                seterror!(x, IncludeLoop)
                return
            end
            oldfile = state.file
            state.file = getfile(state.server, path)
            push!(state.included_files, getpath(state.file))
            setroot(state.file, getroot(oldfile))
            setscope!(getcst(state.file), nothing)
            state(getcst(state.file))
            state.file = oldfile
            pop!(state.included_files)
        elseif !is_in_fexpr(x, CSTParser.defines_function)
            seterror!(x, MissingFile)
        end
    end
end

include("server.jl")
include("imports.jl")
include("references.jl")
include("macros.jl")
include("linting/checks.jl")
include("type_inf.jl")
include("utils.jl")
end


text:
using LinearAlgebra
using FermiCG
using Printf
using Arpack
using Test
using OrderedCollections

@testset "EST" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0, 0, 0.1]))
    push!(atoms,Atom(2,"H",[0, 1,-1]))
    push!(atoms,Atom(3,"H",[0, 1, 1]))
    push!(atoms,Atom(4,"H",[0, 2, 0]))
    push!(atoms,Atom(5,"H",[0, 4, 0]))
    push!(atoms,Atom(6,"H",[0, 5,-1]))
    push!(atoms,Atom(7,"H",[0, 5, 1]))
    push!(atoms,Atom(8,"H",[0, 6, 0]))
    #basis = "6-31g"
    basis = "sto-3g"

    na = 4
    nb = 4

    mol     = Molecule(0,1,atoms,basis)
    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,na,nb)
    # @printf(" FCI Energy: %12.8f\n", e_fci)

    FermiCG.pyscf_write_molden(mol,mf.mo_coeff,filename="scf.molden")

    C = mf.mo_coeff
    rdm_mf = C[:,1:2] * C[:,1:2]'
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Build Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    clusters    = [(1:4),(5:8)]
    init_fspace = [(2,2),(2,2)]

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    rdm1 = zeros(size(ints.h1))
    rdm1a = rdm_mf*.5
    rdm1b = rdm_mf*.5

    rdm1a = rdm1
    rdm1b = rdm1

    display(rdm1a)
    display(rdm1b)

    #e_cmf, D1,D2,temp,temp2 = FermiCG.cmf_ci(ints, clusters, init_fspace, rdm1a, verbose=0)
    #Da = D1
    #Db = D1

    e_cmf, U, Da, Db  = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1, rdm1,
                                       max_iter_oo=40, verbose=0, gconv=1e-6, method="bfgs")
    ints = FermiCG.orbital_rotation(ints,U)

    rdm1a = Da
    rdm1b = Db

    display(Da)
    println()
    display(Db)

    #for ci in clusters
    #    ints_i = subset(ints, ci.orb_list, rdm1a, rdm1b)
    #	print(ints_i.h1)
    #end

    cb_cmf = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=1, max_roots=5,rdm1a=Da,rdm1b=Db,init_fspace=init_fspace)
    cb_est = FermiCG.compute_cluster_est_basis(ints, clusters, rdm1a, rdm1b, thresh_schmidt=1e-4, init_fspace=init_fspace)

    clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters)
    cluster_ops = FermiCG.compute_cluster_ops(cb_est, ints);

    ci_vector = FermiCG.ClusteredState(clusters)
    FermiCG.expand_to_full_space!(ci_vector, cb_est, na, nb)
    display(ci_vector,thresh=-1)

    ovlp = OrderedDict()
    for ci in clusters
    	for (key,value) in cb_est[ci.idx]
            println(key)
            #display(cb_est[ci.idx][key])
            ovlp[ci.idx,key] = cb_cmf[ci.idx].basis[key]'*cb_est[ci.idx].basis[key]
        end
    end

    H = FermiCG.build_full_H(ci_vector, cluster_ops, clustered_ham)
    #display(size(H))
    #display(H)
    #println()

    display(ci_vector,root=1)
    e,v = Arpack.eigs(H, nev = 8, which=:SR)
    for ei in e
        @printf(" Energy: %18.12f\n",real(ei))
    end

    ref_e = [-10.18725871
    -10.15805362
    -10.14789239
    -10.12277607
    -10.11982911
    -10.11553898
    -10.03259372
    -10.02109892]

    @test isapprox(e, ref_e, atol=1e-6)

    S_1_42 =[0.616641;
	  0.766495;
	  0.027335;
	  0.176888;
	  0.013518]
    @test isapprox(abs.(ovlp[1,(4,2)]), abs.(S_1_42), atol=1e-4)

    S_1_41 = [ 0.935067 -0.307226;
          0.208181  0.478093;
         -0.272371 -0.777974;
         -0.090131  0.267947]
    @test isapprox(abs.(ovlp[1,(4,1)]), abs.(S_1_41), atol=1e-4)

    S_1_33 = [-0.32085  -0.        0.851868;
         -0.        0.616641  0.      ;
          0.       -0.766495 -0.      ;
         -0.711982 -0.       -0.110697;
         -0.137815  0.       -0.329198]
    @test isapprox(abs.(ovlp[1,(3,3)]), abs.(S_1_33), atol=1e-4)

    S_1_23 = [ 0.674829  0.229042 -0.228784  0.656271 -0.061039  0.       -0.024558 -0. ;
         -0.635893  0.420242  0.073237  0.548461  0.306119 -0.       -0.023685 -0.      ;
         -0.       -0.        0.        0.       -0.       -0.935067 -0.       -0.307226;
          0.036266 -0.099898  0.164285  0.035816 -0.0381   -0.       -0.387746 -0.      ;
          0.226091  0.792939  0.421423 -0.369864 -0.06511  -0.       -0.002342  0.      ]
    @test isapprox(abs.(ovlp[1,(2,3)]) , abs.(S_1_23), atol=1e-4)
    S_2_42 = [-0.032297;
          0.981557;
         -0.012785;
         -0.187085;
         -0.004896]
    @test isapprox(abs.(ovlp[2,(4,2)]) , abs.(S_2_42), atol=1e-4)
    S_2_41 = [ 0.944586  0.281996;
          0.011955 -0.030939;
         -0.318901  0.908279;
          0.07692   0.30749 ;]
    @test isapprox(abs.(ovlp[2,(4,1)]) , abs.(S_2_41), atol=1e-4)
    S_2_24 = [ 0.032297;
         -0.981557;
         -0.012785;
          0.187085;
         -0.004896]
    @test isapprox(abs.(ovlp[2,(2,4)]) , abs.(S_2_24), atol=1e-4)

end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG