text:
# --------------module TensorBasic-----------
# huangrzh@icloud.com
#---------------------------------------------------

module TensorBasic

# package code goes here
export tensor_product,tensor_permute,tensor_reshape

# Sveral tensor product
# By Huang Ruizhen, 2016/7/8
function tensor_product(A, aindex, B, bindex)
  a_length = length( aindex );
  b_length = length( bindex );

  size_a = size(A);
  size_b = size(B);
  comab = intersect(aindex,bindex);
  if length(comab)==0
    cindex = [aindex;bindex];
    A = reshape(A,length(A));
    B = reshape(B,1,length(B));
    return reshape(A*B,size_a...,size_b...),cindex;
  else
    com_in_a = [findfirst(aindex,xa) for xa in comab];
    com_in_b = [findfirst(bindex,xb) for xb in comab];

    if size_a[com_in_a] != size_b[com_in_b]
      error("Error: Dimentions don't match!");
    end

    diff_in_a = setdiff(collect(1:1:a_length), com_in_a);
    diff_in_b = setdiff(collect(1:1:b_length), com_in_b);
    cindex = [ aindex[diff_in_a]; bindex[diff_in_b] ];

    # mutiply
    if [ com_in_a; diff_in_a ] != collect(1:1:a_length)
        A = permutedims( A, [ com_in_a; diff_in_a ] );
    end
    if [ com_in_b;diff_in_b ] != collect(1:1:b_length)
        B = permutedims( B, [ com_in_b; diff_in_b ] );
    end

    sda = prod(size_a[diff_in_a]);
    sc = prod(size_a[com_in_a]);
    sdb = prod(size_b[diff_in_b]);
    A = reshape(A,sc,sda);
    B = reshape(B,sc,sdb);

    C = A.' * B ;
    return reshape(C,size_a[diff_in_a]..., size_b[diff_in_b]...),cindex;
  end
end

function tensor_product(cindex, A, aindex, B, bindex)
  (C,indc) = tensor_product(A,aindex,B,bindex);
  return tensor_permute(C,indc,cindex);
end

# Tensor permute
# By Huang Ruizhen, 2016.10.15
function tensor_permute(T, oldidx, newidx)
  exchidx = indexin(newidx, oldidx);
  return permutedims(T,exchidx);
end

# hrz, 2016/10/15
function tensor_reshape(T, indx_all, inds)
  size_t = size(T);
  dim_t = ndims(T);
  dim_all = length(indx_all);
  T_ = T;
  if dim_all > dim_t
    size_t = tuple(size_t..., ones(Int64, dim_all-dim_t)...);
    T_ = reshape(T, size_t);
  end

  ind_new = vcat(inds...);
  if indx_all != ind_new
    T_ = tensor_permute(T_, indx_all, ind_new);
  end
  size_t = size(T_);
  n_inds = length(inds);
  dims_i = zeros(Int64, n_inds);
  for i in 1:1:n_inds
    dims_i[i] = prod(size_t[indexin(inds[i],ind_new)]);
  end
  return reshape(T_, dims_i...);
end

end # module


text:
using Samsara
using Test

include("test_Samsara.jl")
include("test_circularsys.jl")
# include("./methods_for_testing_samsara_for_circularsys.jl") # TODO Implementer dette som en test-fil. No er det mest for å teste og utprøve medan eg godter meg over kor bra detteer..
include("./test_linked_node.jl")
include("./test_linked_gate.jl")

ext: jl

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
"""
The sparsity pattern.

- `I`: Input index
- `J`: Ouput index

`(i, j)` means the `j`th element of the output depends on
the `i`th element of the input. Therefore `length(I) == length(J)`
"""
struct Sparsity
    m::Int
    n::Int
    I::Vector{Int} # Input
    J::Vector{Int} # Output
end

SparseArrays.sparse(s::Sparsity) = sparse(s.I, s.J, true, s.m, s.n)

Sparsity(m, n) = Sparsity(m, n, Int[], Int[])

function Base.push!(S::Sparsity, i::Int, j::Int)
    push!(S.I, i)
    push!(S.J, j)
end

# Tags:
struct Input end
struct Output end

struct ProvinanceSet{T}
    set::T # Set, Array, Int, Tuple, anything!
end

# note: this is not strictly set union, just some efficient way of concating
Base.union(p::ProvinanceSet{<:Tuple},
           q::ProvinanceSet{<:Integer}) = ProvinanceSet((p.set..., q.set,))
Base.union(p::ProvinanceSet{<:Integer},
           q::ProvinanceSet{<:Tuple}) = ProvinanceSet((p.set, q.set...,))
Base.union(p::ProvinanceSet{<:Integer},
           q::ProvinanceSet{<:Integer}) = ProvinanceSet((p.set, q.set,))
Base.union(p::ProvinanceSet{<:Tuple},
           q::ProvinanceSet{<:Tuple}) = ProvinanceSet((p.set..., q.set...,))
Base.union(p::ProvinanceSet,
           q::ProvinanceSet) = ProvinanceSet(union(p.set, q.set))
Base.union(p::ProvinanceSet,
           q::ProvinanceSet,
           rs::ProvinanceSet...) = union(union(p, q), rs...)
Base.union(p::ProvinanceSet) = p

function Base.push!(S::Sparsity, i::Int, js::ProvinanceSet)
    for j in js.set
        push!(S, i, j)
    end
end

Cassette.@context SparsityContext

const TagType = Union{Input, Output, ProvinanceSet}
Cassette.metadatatype(::Type{<:SparsityContext}, ::DataType) = TagType

metatype(x, ctx) = hasmetadata(x, ctx) && istagged(x, ctx) && typeof(metadata(x, ctx))
function ismetatype(x, ctx, T)
    hasmetadata(x, ctx) && istagged(x, ctx) && (metadata(x, ctx) isa T)
end

# Dummy type when you getindex
struct Tainted end

# getindex on the input
@inline function Cassette.overdub(ctx::SparsityContext,
                                  f::typeof(getindex),
                                  X::Tagged,
                                  idx::Int...)
    if ismetatype(X, ctx, Input)
        i = LinearIndices(untag(X, ctx))[idx...]
        val = Cassette.fallback(ctx, f, X, idx...)
        tag(val, ctx, ProvinanceSet(i))
    else
        Cassette.recurse(ctx, f, X, idx...)
    end
end

# setindex! on the output
@inline function Cassette.overdub(ctx::SparsityContext,
                                  f::typeof(setindex!),
                                  Y::Tagged,
                                  val::Tagged,
                                  idx::Int...)
    S, path = ctx.metadata
    if ismetatype(Y, ctx, Output)
        set = metadata(val, ctx)
        if set isa ProvinanceSet
            i = LinearIndices(untag(Y, ctx))[idx...]
            push!(S, i, set)
        end
        Cassette.fallback(ctx, f, Y, val, idx...)
    else
        Cassette.recurse(ctx, f, Y, val, idx...)
    end
end

function get_provinance(ctx, arg::Tagged)
    if metadata(arg, ctx) isa ProvinanceSet
        metadata(arg, ctx)
    else
        ProvinanceSet(())
    end
end

get_provinance(ctx, arg) = ProvinanceSet(())

# Any function acting on a value tagged with ProvinanceSet
function _overdub_union_provinance(::Val{eval}, ctx::SparsityContext, f, args...) where {eval}
    idxs = findall(x->ismetatype(x, ctx, ProvinanceSet), args)
    if isempty(idxs)
        Cassette.fallback(ctx, f, args...)
    else
        provinance = union(map(arg->get_provinance(ctx, arg), args[idxs])...)
        if eval
            val = Cassette.fallback(ctx, f, args...)
            tag(val, ctx, provinance)
        else
            tag(Tainted(), ctx, provinance)
        end
    end
end

@inline function Cassette.overdub(ctx::SparsityContext, f, args...)
    haspsets = any(x->ismetatype(x, ctx, ProvinanceSet), args)
    hasinput = any(x->ismetatype(x, ctx, Input), args)
    if haspsets && !hasinput # && !canrecurse(ctx, f, args...)
        _overdub_union_provinance(Val{true}(), ctx, f, args...)
    else
        Cassette.recurse(ctx, f, args...)
    end
end

#=
# Examples:
#
using UnicodePlots

sspy(s::Sparsity) = spy(sparse(s))

julia> sparsity!([0,0,0], [23,53,83]) do Y, X
           Y[:] .= X
           Y == X
       end
(true, Sparsity([1, 2, 3], [1, 2, 3]))

julia> sparsity!([0,0,0], [23,53,83]) do Y, X
           for i=1:3
               for j=i:3
                   Y[j] += X[i]
               end
           end; Y
       end
([23, 76, 159], Sparsity(3, 3, [1, 2, 3, 2, 3, 3], [1, 1, 1, 2, 2, 3]))

julia> sspy(ans[2])
     Sparsity Pattern
     ┌─────┐
   1 │⠀⠄⠀⠀⠀│ > 0
   3 │⠀⠅⠨⠠⠀│ < 0
     └─────┘
     1     3
     nz = 6

julia> sparsity!(f, zeros(Int, 3,3), [23,53,83])
([23, 53, 83], Sparsity(9, 3, [2, 5, 8], [1, 2, 3]))

julia> sspy(ans[2])
     Sparsity Pattern
     ┌─────┐
   1 │⠀⠄⠀⠀⠀│ > 0
     │⠀⠀⠠⠀⠀│ < 0
   9 │⠀⠀⠀⠐⠀│
     └─────┘
     1     3
     nz = 3
=#


text:
# ustring.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files ustring.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
   Unicode string handling functions

   These C API functions provide general Unicode string handling.

   Some functions are equivalent in name, signature, and behavior to the ANSI C <string.h>
   functions. (For example, they do not check for bad arguments like NULL string pointers.)
   In some cases, only the thread-safe variant of such a function is implemented here
   (see u_strtok_r()).

   Other functions provide more Unicode-specific functionality like locale-specific
   upper/lower-casing and string comparison in code point order.

   ICU uses 16-bit Unicode (UTF-16) in the form of arrays of UChar code units.
   UTF-16 encodes each Unicode code point with either one or two UChar code units.
   (This is the default form of Unicode, and a forward-compatible extension of the original,
   fixed-width form that was known as UCS-2. UTF-16 superseded UCS-2 with Unicode 2.0
   in 1996.)

   Some APIs accept a 32-bit UChar32 value for a single code point.

   ICU also handles 16-bit Unicode text with unpaired surrogates.
   Such text is not well-formed UTF-16.
   Code-point-related functions treat unpaired surrogates as surrogate code points,
   i.e., as separate units.

   Although UTF-16 is a variable-width encoding form (like some legacy multi-byte encodings),
   it is much more efficient even for random access because the code unit values
   for single-unit characters vs. lead units vs. trail units are completely disjoint.
   This means that it is easy to determine character (code point) boundaries from
   random offsets in the string.

   Unicode (UTF-16) string processing is optimized for the single-unit case.
   Although it is important to support supplementary characters
   (which use pairs of lead/trail code units called "surrogates"),
   their occurrence is rare. Almost all characters in modern use require only
   a single UChar code unit (i.e., their code point values are <=0xffff).

   For more details see the User Guide Strings chapter (http://icu-project.org/userguide/strings.html).
   For a discussion of the handling of unpaired surrogates see also
   Jitterbug 2145 and its icu mailing list proposal on 2002-sep-18.
"""
module ustring end

macro libstr(s)     ; _libicu(s, iculib,     "u_str")     ; end

_tolower(dest::Ptr{UInt16}, destsiz, src::Ptr{UInt16}, srclen, err) =
    ccall(@libstr(ToLower), Cint,
          (Ptr{UChar}, Cint, Ptr{UChar}, Cint, Ptr{UInt8}, Ptr{UErrorCode}),
          dest, destsiz, src, srclen, locale[], err)

_toupper(dest::Ptr{UInt16}, destsiz, src::Ptr{UInt16}, srclen, err) =
    ccall(@libstr(ToUpper), Cint,
          (Ptr{UChar}, Cint, Ptr{UChar}, Cint, Ptr{UInt8}, Ptr{UErrorCode}),
          dest, destsiz, src, srclen, locale[], err)

_foldcase(dest::Ptr{UInt16}, destsiz, src::Ptr{UInt16}, srclen, err) =
    ccall(@libstr(FoldCase), Cint,
          (Ptr{UChar}, Cint, Ptr{UChar}, Cint, Cint, Ptr{UErrorCode}),
          dest, destsiz, src, srclen, 0, err)

_totitle(dest::Ptr{UInt16}, destsiz, src::Ptr{UInt16}, srclen, breakiter, err) =
    ccall(@libstr(ToTitle), Cint,
          (Ptr{UChar}, Cint, Ptr{UChar}, Cint, Ptr{Cvoid}, Ptr{UInt8}, Ptr{UErrorCode}),
          dest, destsiz, src, srclen, breakiter, locale[], err)

"""
   Case-folds the characters in a string.

   Case-folding is locale-independent and not context-sensitive, but there is an option for whether to
   include or exclude mappings for dotted I and dotless i that are marked with 'T' in CaseFolding.txt.

   The result may be longer or shorter than the original.
   The source string and the destination buffer are allowed to overlap.

   Arguments:
   src       The original string
   options   Either U_FOLD_CASE_DEFAULT or U_FOLD_CASE_EXCLUDE_SPECIAL_I

   Returns:  Case-folded string
"""
function foldcase end

for f in (:tolower, :toupper, :foldcase)
    uf = Symbol(string('_',f))
    @eval begin
        function ($f)(src::T) where {T<:WordStrings}
            @preserve src begin
                dstlen = srclen = ncodeunits(src)
                dest, pnt = _allocate(UInt16, dstlen)
                err = Ref{UErrorCode}(0)
                siz = ($uf)(pnt, dstlen, pointer(src), srclen, err)
                # Retry with large enough buffer if got buffer overflow
                if err[] == U_BUFFER_OVERFLOW_ERROR
                    err[] = 0
                    dstlen = siz
                    dest, pnt = _allocate(UInt16, dstlen)
                    siz = ($uf)(pnt, dstlen, pointer(src), srclen, err)
                end
                FAILURE(err[]) && error("ICU: internal error $(err[]) failed to map case")
                siz != dstlen ? Str(cse(T), dest[1:siz<<1]) : Str(cse(T), dest)
            end
        end
    end
end

"""
   Titlecase a string.

   Casing is locale-dependent and context-sensitive.
   Titlecasing uses a break iterator to find the first characters of words that are to be titlecased.
   It titlecases those characters and lowercases all others.

   The titlecase break iterator can be provided to customize for arbitrary styles, using rules and
   dictionaries beyond the standard iterators.  It may be more efficient to always provide an iterator
   to avoid opening and closing one for each string.
   The standard titlecase iterator for the root locale implements the algorithm of Unicode TR 21.

   This function uses only the setText(), first() and next() methods of the provided break iterator.

   The result may be longer or shorter than the original.

   Arguments:
   src       The original string
   titleIter A break iterator to find the first characters of words that are to be titlecased.
             If none is provided, then a standard titlecase break iterator is opened.
   locale    The locale to consider, or "" for the root locale (not passed -> default locale)

   Returns: Title-cased string
"""
function totitle end

function totitle(src::T, bi) where {T<:WordStrings}
    dstlen = srclen = ncodeunits(src)
    @preserve src begin
        dest, pnt = _allocate(UInt16, dstlen)
        err = Ref{UErrorCode}(0)
        siz = _totitle(pnt, dstlen, pointer(src), srclen, bi, err)
        # Retry with large enough buffer if got buffer overflow
        if err[] == U_BUFFER_OVERFLOW_ERROR
            err[] = 0
            dstlen = siz
            dest, pnt = _allocate(UInt16, siz)
            siz = _totitle(pnt, dstlen, pointer(src), srclen, bi, err)
        end
        FAILURE(err[]) && error("ICU: internal error $(err[]) failed to map case")
        siz != dstlen ? Str(cse(T), dest[1:siz<<1]) : Str(cse(T), dest)
    end
end
totitle(str::T) where {T<:WordStrings} = totitle(str, get_break_iterator())


text:
const utf8_offset = [
    0x00000000, 0x00003080,
    0x000e2080, 0x03c82080,
    0xfa082080, 0x82082080,
]

const utf8_trailing = [
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,
]

function endof(s::MutableUTF8String)
    d = s.data
    i = length(d)
    i == 0 && return i
    while !Base.is_utf8_start(d[i])
        i -= 1
    end
    i
end
length(s::MutableUTF8String) = int(ccall(:u8_strlen, Csize_t, (Ptr{Uint8},), s.data))

next(s::MutableUTF8String, i::Int) = (s[i], i+1+utf8_trailing[s.data[i]+1])

isvalid(s::MutableUTF8String, i::Integer) = (1 <= i <= endof(s.data)) && Base.is_utf8_start(s.data[i])

# common base method for utf8 strings would avoid duplication of getindex
function getindex(s::MutableUTF8String, i::Int)
    d = s.data
    b = d[i]
    if !Base.is_utf8_start(b)
        j = i-1
        while 0 < j && !Base.is_utf8_start(d[j])
            j -= 1
        end
        if 0 < j && i <= j+utf8_trailing[d[j]+1] <= length(d)
            # b is a continuation byte of a valid UTF-8 character
            error("invalid UTF-8 character index")
        end
        return '\ufffd'
    end
    trailing = utf8_trailing[b+1]
    if length(d) < i + trailing
        return '\ufffd'
    end
    c::Uint32 = 0
    for j = 1:trailing+1
        c <<= 6
        c += d[i]
        i += 1
    end
    c -= utf8_offset[trailing+1]
    char(c)
end

function getindex(s::MutableUTF8String, r::UnitRange{Int})
    isempty(r) && return empty_utf8
    i, j = first(r), last(r)
    d = s.data
    if !Base.is_utf8_start(d[i])
        i = nextind(s,i)
    end
    if j > endof(s)
        throw(BoundsError())
    end
    j = nextind(s,j)-1
    UTF8String(d[i:j])
end

function setindex!(s::MutableUTF8String, x, i0::Real)
    d = s.data
    !Base.is_utf8_start(d[i0]) && error("invalid UTF-8 character index")
    endi = i0+utf8_trailing[d[i0]+1]
    setindex!(s, string(x), i0:endi)
end
function setindex!(s::MutableUTF8String, r::ByteString, I::UnitRange)
    st, ed = first(I), last(I)
    ed = nextind(s, ed)-1
    _splice!(s.data, I, r.data)
    r
end
function setindex!(s::MutableUTF8String, c::Char, I::UnitRange)
    st, ed = first(I), last(I)
    ed = nextind(s, ed)-1
    if isascii(c)
        s.data[st:ed] = uint8(c)
    else
        iob = IOBuffer()
        print(iob, c)
        nd = nextind(s, st)
        while nd <= ed
            print(iob, c)
            nd = nextind(s, nd)
        end
        _splice!(s.data, st:ed, takebuf_array(iob))
    end
    c
end

search(s::MutableUTF8String, c::Char, i::Integer) = search(utf8(s), c, i)
rsearch(s::MutableUTF8String, c::Char, i::Integer) = rsearch(utf8(s), c, i)

function string(a::MutableString...)
    data = Array(Uint8,0)
    for d in a
        append!(data,d.data)
    end
    MutableUTF8String(data)
end

function reverse!(s::MutableUTF8String)
    indata = copy(s.data)
    if ccall(:u8_reverse, Cint, (Ptr{Uint8}, Ptr{Uint8}, Csize_t), s.data, indata, length(indata)) == 1
        error("invalid UTF-8 data")
    end
    nothing
end

lcfirst!(s::MutableUTF8String) = (s[1] = lowercase(s[1]); nothing)
ucfirst!(s::MutableUTF8String) = (s[1] = uppercase(s[1]); nothing)

function uppercase!(s::MutableUTF8String)
    (length(s) == 0) && return
    st = 1
    iob = IOBuffer()
    d = s.data
    while !done(s, st)
        c, ed = next(s, st)
        seek(iob, 0)
        print(iob, uppercase(c))
        p = position(iob)
        _splice!(d, st:(ed-1), iob.data, p)
        st += p
    end
end

function lowercase!(s::MutableUTF8String)
    (length(s) == 0) && return
    st = 1
    iob = IOBuffer()
    d = s.data
    while !done(s, st)
        c, ed = next(s, st)
        print(iob, lowercase(c))
        p = position(iob)
        _splice!(d, st:(ed-1), iob.data, p)
        seek(iob, 0)
        st += p
    end
end

function replace!{T}(str::MutableUTF8String, pattern::T, repl::Union(Function,ByteString,Char), limit::Integer=0)
    n = 1
    es = endof(str)
    r = search(str,pattern,1)
    j, k = first(r), last(r)
    k = nextind(str,k) - 1
    while j != 0
        rv::Union(ByteString,Char) = isa(repl, Function) ? repl(SubString(str, j, k)) : repl
        if isa(rv, ByteString)
            _splice!(str.data, j:k, rv.data)
            lrv = length(rv.data)
            lr = k-j+1
            if lrv != lr
                k = j+lrv
                es = es - lr + lrv
            end
        else
            str[r] = rv
        end
        ((k+1) > es) && break
        r = search(str, pattern, k+1)
        j, k = first(r), last(r)
        k = nextind(str,k) - 1
        n == limit && break
        n += 1
    end
end

function replace!(str::MutableUTF8String, re::Regex, repl::Union(Function,ByteString,Char), limit::Integer=0)
    regex = Base.compile(re).regex
    extra = re.extra
    n = length(str.data)
    offset = int32(0)
    opts = re.options & Base.PCRE.EXECUTE_MASK
    opts_nonempty = opts | Base.PCRE.ANCHORED | Base.PCRE.NOTEMPTY_ATSTART
    prevempty = false
    ovec = Array(Int32, 3)
    n_repl = 1
    while true
        result = ccall((:pcre_exec, :libpcre), Int32,
                       (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Int32,
                       Int32, Int32, Ptr{Int32}, Int32),
                       regex, extra, str.data, n,
                       offset, prevempty ? opts_nonempty : opts, ovec, 3)

        if result < 0
            if prevempty && offset < n
                offset = int32(nextind(str, offset + 1) - 1)
                prevempty = false
                continue
            else
                break
            end
        end

        rv::Union(ByteString,Char) = isa(repl, Function) ? repl(SubString(str, ovec[1]+1, ovec[2])) : repl
        r = (ovec[1]+1):(nextind(str, ovec[2])-1)
        if isa(rv, ByteString)
            _splice!(str.data, r, rv.data)

            lrv = length(rv)
            lr = length(r)
            if lrv != lr
                n = n - lr + lrv
            end

            prevempty = offset == ovec[2]
            offset = ovec[1]+1+lrv
        else
            str[r] = rv
            prevempty = offset == ovec[2]
            offset = ovec[2]
        end
        n_repl == limit && break
        n_repl += 1
    end
    nothing
end

function map!(f, s::MutableUTF8String)
    st = 1
    d = s.data
    while !done(s, st)
        s[st] = f(s[st])
        st = st + 1 + utf8_trailing[d[st] + 1]
    end
end


text:
module ProximalAlgorithms

using ProximalCore
using ProximalCore: prox, prox!, gradient, gradient!

const RealOrComplex{R} = Union{R,Complex{R}}
const Maybe{T} = Union{T,Nothing}

# various utilities

include("utilities/ad.jl")
include("utilities/fb_tools.jl")
include("utilities/iteration_tools.jl")

# acceleration utilities

include("accel/traits.jl")
include("accel/lbfgs.jl")
include("accel/anderson.jl")
include("accel/nesterov.jl")
include("accel/broyden.jl")
include("accel/noaccel.jl")

# algorithm interface

struct IterativeAlgorithm{IteratorType, H, S, D, K}
    maxit::Int
    stop::H
    solution::S
    verbose::Bool
    freq::Int
    display::D
    kwargs::K
end

"""
    IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, display, kwargs...)

Wrapper for an iterator type `T`, adding termination and verbosity options on top of it.

This is a conveniency constructor to allow for "partial" instantiation of an iterator of type `T`.
The resulting "algorithm" object `alg` can be called on a set of keyword arguments, which will be merged
to `kwargs` and passed on to `T` to construct an iterator which will be looped over.
Specifically, if an algorithm is constructed as

    alg = IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, display, kwargs...)

then calling it with

    alg(; more_kwargs...)

will internally loop over an iterator constructed as

    T(; alg.kwargs..., more_kwargs...)

# Note
This constructor is not meant to be used directly: instead, algorithm-specific constructors
should be defined on top of it and exposed to the user, that set appropriate default functions
for `stop`, `solution`, `display`.

# Arguments
* `T::Type`: iterator type to use
* `maxit::Int`: maximum number of iteration
* `stop::Function`: termination condition, `stop(::T, state)` should return `true` when to stop the iteration
* `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
* `verbose::Bool`: whether the algorithm state should be displayed
* `freq::Int`: every how many iterations to display the algorithm state
* `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
* `kwargs...`: keyword arguments to pass on to `T` when constructing the iterator
"""
IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, display, kwargs...) =
    IterativeAlgorithm{T, typeof(stop), typeof(solution), typeof(display), typeof(kwargs)}(
        maxit, stop, solution, verbose, freq, display, kwargs
    )

function (alg::IterativeAlgorithm{IteratorType})(; kwargs...) where IteratorType
    iter = IteratorType(; alg.kwargs..., kwargs...)
    for (k, state) in enumerate(iter)
        if k >= alg.maxit || alg.stop(iter, state)
            alg.verbose && alg.display(k, iter, state)
            return (alg.solution(iter, state), k)
        end
        alg.verbose && mod(k, alg.freq) == 0 && alg.display(k, iter, state)
    end
end

# algorithm implementations

include("algorithms/forward_backward.jl")
include("algorithms/fast_forward_backward.jl")
include("algorithms/zerofpr.jl")
include("algorithms/panoc.jl")
include("algorithms/douglas_rachford.jl")
include("algorithms/drls.jl")
include("algorithms/primal_dual.jl")
include("algorithms/davis_yin.jl")
include("algorithms/li_lin.jl")
include("algorithms/sfista.jl")
include("algorithms/panocplus.jl")

end # module


text:
using DataFrames

"""
    accuracy(target, output; normalize = true)

If `normalize` is `true`, the fraction of correctly classified observations is returned.  Otherwise
the total number is returned.

Adapted from https://github.com/JuliaML/MLMetrics.jl/blob/c03435006ab64f84376835ad588d9cbf546da506/src/classification/multiclass.jl#L194.
"""
function accuracy{T}(target::AbstractVector{T},
                     output::AbstractVector{T};
                     normalize = true)
    @assert length(target) == length(output)
    correct = 0

    @inbounds for i in eachindex(target, output)
        correct += target[i] == output[i]
    end

    normalize ? Float64(correct / length(target)) : Float64(correct)
end

"""
    kfolds(n::Integer, [k = 5]) -> Tuple

Compute the train/validation assignments for `k` repartitions of `n` observations, and return them
in the form of two vectors. The first vector contains the index-vectors for the training subsets,
and the second vector the index-vectors for the validation subsets respectively.

From https://github.com/JuliaML/MLDataPattern.jl/blob/d9fcafe1cdd7230e249d3ae85163eef088e3d659/src/folds.jl#L226.
"""
function kfolds(n::Integer, k::Integer = 5)
    2 <= k <= n || throw(ArgumentError("n must be positive and k must to be within 2:$(max(2,n))"))
    # Compute the size of each fold. This is important because
    # in general the number of total observations might not be
    # divideable by k. In such cases it is custom that the remaining
    # observations are divided among the folds. Thus some folds
    # have one more observation than others.
    sizes = fill(floor(Int, n / k), k)
    for i = 1:(n % k)
        sizes[i] = sizes[i] + 1
    end

    # Compute start offset for each fold
    offsets = cumsum(sizes) .- sizes .+ 1
    # Compute the validation indices using the offsets and sizes
    val_indices = map((o, s) -> (o:o+s-1), offsets, sizes)
    # The train indices are then the indicies not in validation
    train_indices = map(idx -> setdiff(1:n, idx), val_indices)
    # We return a tuple of arrays
    train_indices, val_indices
end

function create_fitness{N, C}(data, Vn::Type{Val{N}}, Vc::Type{Val{C}};
                              size_factor::Float64 = 0.0, depth_factor::Float64 = 0.0)
    # it is more efficient to iterate along columns
    features = convert(Matrix{Float64}, data[1:end-1]) |> transpose
    classes = Classification{N, C}.(getfield.(data[end], :level))

    function fitness(tree::DecisionTree{N, C})
        correct = sum(eachindex(classes)) do i
            classes[i] == decide(view(features, :, i), tree)
        end

        penalty = size_factor * treesize(tree) + depth_factor * treedepth(tree)
        correct - penalty
    end

    function accuracy(tree::DecisionTree{N, C})
        correct = sum(eachindex(classes)) do i
            classes[i] == decide(view(features, :, i), tree)
        end

        correct / length(classes)
    end

    fitness, accuracy
end

max_stars_repo_name: phipsgabler/soft-computing-final

max_issues_repo_name: phipsgabler/soft-computing-final

max_forks_repo_name: phipsgabler/soft-computing-final


text:
using Scientific_Programming_in_Julia
using EcosystemCore
using Random
Random.seed!(2)

n_grass  = 200
max_size = 10

n_sheep         = 10
Δenergy_sheep   = 0.2
energy_sheep    = 4.0
sheep_reproduce = 0.8
sheep_foodprob  = 0.6

n_wolves       = 2
Δenergy_wolf   = 8.0
energy_wolf    = 10.0
wolf_reproduce = 0.1
wolf_foodprob  = 0.2

# all animals can mate - regardless of sex
EcosystemCore.mates(a::Animal{S},b::Animal{S}) where S<:Species = true
EcosystemCore.mates(a::Agent, b::Agent) = false

gs = [Grass(id,max_size) for id in 1:n_grass]
ss = [Sheep(id,energy_sheep,Δenergy_sheep,sheep_reproduce, sheep_foodprob) for id in (n_grass+1):(n_grass+n_sheep)]
ws = [Wolf(id,energy_wolf,Δenergy_wolf,wolf_reproduce, wolf_foodprob) for id in (n_grass+n_sheep+1):(n_grass+n_sheep+n_wolves)]
as = vcat(gs,ss,ws)

w = World(as)

counts = Dict(n=>[c] for (n,c) in agent_count(w))
function _save(w::World)
    for (n,c) in agent_count(w)
        push!(counts[n],c)
    end
end

logcb = every_nth(w->(@info w), 1)
savecb = every_nth(_save, 1)
cbs = [logcb, savecb]

simulate!(w, 200, callbacks=cbs)

using Plots
plt = plot()
for (n,c) in counts
    plot!(plt, c, label="$n", lw=2)
end
savefig(plt, "docs/src/lecture_02/pred-prey.png")
# savefig(plt, "docs/src/lecture_03/pred-prey.png")
display(plt)


text:
import Yao.Registers: tracedist

export QuGAN, psi, toy_qugan, QuGANGo!
"""
Quantum GAN.

Reference:
    Benedetti, M., Grant, E., Wossnig, L., & Severini, S. (2018). Adversarial quantum circuit learning for pure state approximation, 1–14.
"""
struct QuGAN{N} <: QCOptProblem
    target::DefaultRegister
    generator::MatrixBlock{N}
    discriminator::MatrixBlock
    reg0::DefaultRegister
    witness_op::MatrixBlock
    circuit::AbstractBlock
    gdiffs
    ddiffs

    function QuGAN(target::DefaultRegister, gen::MatrixBlock, dis::MatrixBlock)
        N = nqubits(target)
        c = sequence(gen, addbit(1), dis)
        witness_op = put(N+1, (N+1)=>P0)
        gdiffs = collect(gen, AbstractDiff)
        ddiffs = collect(dis, AbstractDiff)
        new{N}(target, gen, dis, zero_state(N), witness_op, c, gdiffs, ddiffs)
    end
end

# INTERFACES
circuit(qg::QuGAN) = qg.circuit
diff_blocks(qg::QuGAN) = [qg.gdiffs...; qg.ddiffs...]
loss(qg::QuGAN) = p0t(qg) - p0g(qg)
function gradient(qg::QuGAN)
    ggrad_g = opdiff.(()->psi_discgen(qg), qg.gdiffs, Ref(qg.witness_op))
    dgrad_g = opdiff.(()->psi_discgen(qg), qg.ddiffs, Ref(qg.witness_op))
    dgrad_t = opdiff.(()->psi_disctarget(qg), qg.ddiffs, Ref(qg.witness_op))
    [-ggrad_g; dgrad_t - dgrad_g]
end

"""probability to get evidense qubit 0 on generation set."""
p0g(qg::QuGAN) = expect(qg.witness_op, psi_discgen(qg)) |> real
"""probability to get evidense qubit 0 on target set."""
p0t(qg::QuGAN) = expect(qg.witness_op, psi_disctarget(qg)) |> real
"""generated wave function"""
psi(qg::QuGAN) = copy(qg.reg0) |> qg.generator
"""input |> generator |> discriminator"""
psi_discgen(qg::QuGAN) = copy(qg.reg0) |> qg.circuit
"""target |> discriminator"""
psi_disctarget(qg::QuGAN) = copy(qg.target) |> qg.circuit[2:end]
"""tracedistance between target and generated wave function"""
tracedist(qg::QuGAN) = tracedist(qg.target, psi(qg))[]

"""
    toy_qugan(target::DefaultRegister, depth_gen::Int, depth_disc::Int) -> QuGAN

Construct a toy qugan.
"""
function toy_qugan(target::DefaultRegister, depth_gen::Int, depth_disc::Int)
    n = nqubits(target)
    generator = dispatch!(random_diff_circuit(n, depth_gen, pair_ring(n)), :random) |> autodiff(:QC)
    discriminator = dispatch!(random_diff_circuit(n+1, depth_disc, pair_ring(n+1)), :random) |> autodiff(:QC)
    return QuGAN(target, generator, discriminator)
end

# Optimization
"""
    QuGANGo!{QT<:QuGAN, OT} <: QCOptGo!{QT}

Iterative training of quantum generative optimization problem,
QT is the type of quantum optimization problem,
OT is the optimizer/learning_rate parameter type.
"""
struct QuGANGo!{QT<:QuGAN, OT} <: QCOptGo!{QT}
    qop::QT
    goptim::OT
    doptim::OT
    niter::Int
end
QuGANGo!(qg::QuGAN, glr, dlr) = QuGANGo!(qg, glr, dlr, typemax(Int))
Base.length(qgg::QuGANGo!) = qgg.niter

function Base.iterate(qgg::QuGANGo!{<:Any, <:Real}, state=1)
    state > qgg.niter && return nothing

    qg = qgg.qop
    ng = length(qg.gdiffs)
    grad = gradient(qg)
    dispatch!(+, qg.generator, -grad[1:ng]*qgg.goptim)
    dispatch!(-, qg.discriminator, -grad[ng+1:end]*qgg.doptim)
    Dict("step"=>state,"gradient"=>grad), state+1
end

function Base.iterate(qgg::QuGANGo!{<:Any, <:Adam}, state=(1, parameters(qgg.qop.generator), parameters(qgg.qop.discriminator)))
    state[1] > qgg.niter && return nothing

    qg = qgg.qop
    ng = length(qg.gdiffs)
    grad = gradient(qg)
    dispatch!(qg.generator, update!(state[2], grad[1:ng], qgg.goptim))
    dispatch!(qg.discriminator, update!(state[3], -grad[ng+1:end], qgg.doptim))
    Dict("step"=>state[1], "gradient"=>grad), (state[1]+1, state[2], state[3])
end


text:
"""
# `module Data`

Provides methods to read files containing simulation data. Primarily this is
intended to load `.xyz` files and convert them to JuLIP-compatible data:
```
data = IPFitting.Data.load_data("mydata.xyz")
```
where `mydata.xyz` contains multiple configurations, will read in those
configurations, then convert them into a `JuLIP.Atoms` object, extract
energies, forces, etc store each configuration in a `Dat` and return them
as a `Vector{Dat}`. This can then be stored using `JLD` or `JLD2`.

If `d::Dat`, then one can call `energy(d), forces(d)`, etc to extract the
loaded information.
"""
module Data

using JuLIP, ProgressMeter, FileIO, Printf, StatsBase, PrettyTables, DataFrames
using IPFitting: Dat, vec_obs, devec_obs, observation, hasobservation
using IPFitting.DataTypes
using StringDistances
import JuLIP: Atoms, energy, forces, virial
import Base: length, Dict
using PyCall

import ASE     # use ASE since this will have already figured out how to
               # load `ase` without problems
using ASE: ASEAtoms
#ase_io = ASE.ase_io

const ase_io = PyNULL()

function __init__()
   copy!(ase_io, pyimport_conda("ase.io", "ase", "rmg"))
end

export configtype, weight, load_data, truncate_string

Atoms(d::Dat) = d.at
length(d::Dat) = length(d.at)
configtype(d::Dat) = d.configtype
energy(d::Dat) = haskey(d.D, ENERGY) ? devec_obs(Val(:E), d.D[ENERGY]) : nothing
forces(d::Dat) = haskey(d.D, FORCES) ? devec_obs(Val(:F), d.D[FORCES]) : nothing
virial(d::Dat) = haskey(d.D, VIRIAL) ? devec_obs(Val(:V), d.D[VIRIAL]) : nothing

function read_energy(atpy, energy_key)
   for key in keys(atpy.info)
      if lowercase(key) == lowercase(energy_key)
         return atpy.info[key]
      end
   end
   return nothing
end

function read_forces(atpy, force_key)
   for key in keys(atpy.arrays)
      if lowercase(key) == lowercase(force_key)
         return atpy.arrays[key]' |> vecs
      end
   end
   return nothing
end

function read_virial(atpy, virial_key)
   for key in keys(atpy.info)
      if lowercase(key) == lowercase(virial_key)
         return JMat(atpy.info[key]...)
      end
   end
   # for key in keys(atpy.info)
   #    if compare(lowercase(key), "stress", Levenshtein()) > 0.8
   #       @info("No virial key found: converting stress to virial!")
   #       return -1.0 * JMat(atpy.info[key]...) * atpy.get_volume()
   #    end
   # end
   return nothing
end

function read_configtype(atpy)
   for key in keys(atpy.info)
      if lowercase(key) == "config_type"
         return atpy.info[key]
      end
   end
   return nothing
end

function truncate_string(s, n)
   if length(s) <= n
      return s
   end
   n1 = n2 = (n-2) ÷ 2
   if n1 + n2 + 2 < n
      n1 += 1
   end
   return "$(s[1:n1])..$(s[end-n2+1:end])"
end

function count_scalars(obs)
   if obs == nothing
      return 0
   else
      return length(hcat(obs...))
   end
end

function keys_info(key_dict)
   s = "keys found: "
   for key in sort(collect(keys(key_dict)))
      n = key_dict[key]
      s *= "\"$(key)\" [$n], "
   end
   return s[1:end-2]
end

read_Atoms(atpy) = Atoms(ASEAtoms( atpy ))

"""
```
function read_xyz(fname; verbose=true, index = ":",
                         include = nothing, exclude = nothing)
```
Loads the atoms objects contained in an xyz file, attempts to read the
DFT data stored inside and returns a `Vector{Dat}`.
"""
function read_xyz(fname; energy_key = "dft_energy", force_key = "dft_force", virial_key = "dft_virial", auto=false, verbose=true, index = ":",
                         include = nothing, exclude = nothing)
   if verbose
      println("Reading in $fname ...")
   end
   at_list = ase_io.read(fname, index=index)
   #at_list = ase_io_read(fname, index=index)
   data = Vector{Dat}(undef, length(at_list))
   idx = 0
   if verbose
      println("Processing data ...")
      # tic()
   end

   info_dict = countmap(collect(Iterators.flatten([collect(keys(at.info)) for at in at_list])))
   array_dict = countmap(collect(Iterators.flatten([collect(keys(at.arrays)) for at in at_list])))

   if auto
      energy_key = findmax([[compare(energy_key, key, Levenshtein()), key] for key in keys(info_dict)])[1][2]
      force_key = findmax([[compare(force_key, key, Levenshtein()), key] for key in keys(array_dict)])[1][2]
      virial_key = findmax([[compare(virial_key, key, Levenshtein()), key] for key in keys(info_dict)])[1][2]
   end

   @info("Keys used: E => \"$(energy_key)\", F => \"$(force_key)\", V => \"$(virial_key)\"")

   if verbose
      for (key,val) in Dict("Info " => info_dict, "Array " => array_dict)
         @info(key * keys_info(val))
      end
   end

   ct_dict = Dict()
   ct_datf = DataFrame()#A = Any[], B = Any[], C = Any[], D = Any[], E = Any[], F = Any[])

   dt = verbose ? 1.0 : 0.0
   @showprogress dt for (i,atpy) in enumerate(at_list)

      # get the config type and decide whether to keep or skip this config
      ct = read_configtype(atpy)
      if ct == nothing
         verbose && @warn("$(idx+1) has no config_type")
         ct = "nothing"
      end
      if exclude != nothing
         if ct in exclude
            continue
         end
      end
      if include != nothing
         if !(ct in include)
            continue
         end
      end

      idx += 1
      at = read_Atoms(atpy)
      E = read_energy(atpy, energy_key)
      F = read_forces(atpy, force_key)
      V = read_virial(atpy, virial_key)

      E_c = count_scalars(E)
      F_c = count_scalars(F)
      V_c = count_scalars(V)

      if ct in names(ct_datf)
         ct_datf[!, Symbol(ct)] .+= [1, length(atpy), E_c, F_c, V_c]
      else
         ct_datf[!, Symbol(ct)] = [1, length(atpy), E_c, F_c, V_c]
      end

      EFV = ""
      (E != nothing) && (EFV *= "E")
      (F != nothing) && (EFV *= "F")
      (V != nothing) && (EFV *= "V")

      data[idx] = Dat( at,
                       ct;
                       E = E, F = F, V = V )
   end
   # verbose && toc()

   ct_datf_trans = DataFrame([[names(ct_datf)]; collect.(eachrow(ct_datf))], [:column; Symbol.(axes(ct_datf, 1))]) #apparently only "simple" way to take a dataframe transpose......
   rename!(ct_datf_trans, Symbol.(["config_type", "#cfgs", "#envs", "#E", "#F", "#V"]))

   totals = sum.(eachcol(ct_datf_trans[:, 2:end]))
   missings = [0, 0, abs(totals[1] - totals[3]), abs(3*totals[2] - totals[4]), abs(9*totals[1] - totals[5])]
   push!(ct_datf_trans, vcat(["total", totals]...))
   push!(ct_datf_trans, vcat(["missing", missings]...))

   pretty_table(ct_datf_trans, crop = :horizontal, body_hlines = [length(ct_datf_trans[!, 1])-2])
   #
   # totals = [0,0,0,0,0]
   #
   # print("┏━━━━━━━━━━━━━━┳━━━━━━━┳━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┓\n")
   # print("┃ config type  ┃ #cfgs ┃ #envs  ┃      #E      ┃      #F      ┃      #V      ┃\n")
   # print("┠──────────────╂───────╂────────╂──────────────╂──────────────╂──────────────┨\n")
   # for ct in sort(collect(keys(ct_dict)))
   #    s = @sprintf("┃ %12s ┃ %5s ┃ %6s ┃ %6s [%4s]┃ %6s [%4s]┃ %6s [%4s]┃\n",
   #       truncate_string(ct, 12), ct_dict[ct][1], ct_dict[ct][2],
   #          ct_dict[ct][3], abs(ct_dict[ct][3] - ct_dict[ct][1]),
   #          ct_dict[ct][4], abs(ct_dict[ct][4] - 3*ct_dict[ct][2]),
   #          ct_dict[ct][5], abs(ct_dict[ct][5] - 9*ct_dict[ct][1]))
   #    print(s)
   #    for i in 1:5
   #       totals[i] += ct_dict[ct][i]
   #    end
   # end
   #
   # print("┠──────────────╂───────╂────────╂──────────────╂──────────────╂──────────────┨\n")
   # s = @sprintf("┃    totals    ┃ %5s ┃ %6s ┃ %12s ┃ %12s ┃ %12s ┃\n",
   #          totals[1], totals[2], totals[3], totals[4], totals[5])
   # print(s)
   # print("┗━━━━━━━━━━━━━━┻━━━━━━━┻━━━━━━━━┻━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━┛\n")

   return data[1:idx]
end

function load_data(fname; kwargs...)
   if fname[end-3:end] == ".xyz"
      return read_xyz(fname; kwargs...)
   elseif fname[end-3:end] == ".jld" || fname[end-4:end] == ".jld2"
      return load(fname)
   end
   error("unknown file format $(fname[end-3:end])")
end

end


text:
module KannanBioinformaticsBISC195

export normalizeDNA,
       composition,
       gc_content,
       complement,
       reverse_complement,
       parse_fasta

# # uncomment the following line if you intend to use BioSequences types
# using BioSequences

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("invalid base $base")
    end
    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

function composition(seq)
    sequence = normalizeDNA(seq) # make uppercase string, check invalid bases
    a = c = g = t = n = 0 # sets all 4 variables to `0`
    for base in sequence
        a = count(==('A'), sequence)
        c = count(==('C'), sequence)
        g = count(==('G'), sequence)
        t = count(==('T'), sequence)
        n = count(==('N'), sequence)
    end
    composition = Dict('A' => a, 'G' => g, 'T' => t, 'N' => n, 'C' => c)
    return composition
end

function gc_content(seq)
    comp = composition(seq)
    gs = comp['G']
    cs = comp['C']
    seqlength = length(seq)
    return (gs+cs) / (seqlength)
end

function complement(base::Char)
    complements = Dict('A' => 'T',
                       'T' => 'A',
                       'G' => 'C',
                       'C' => 'G',
                       'N' => 'N')

    !(base in keys(complements)) && error("Invalid base $base")
    return complements[base]
end

function complement(seq::AbstractString)
    seq = normalizeDNA(seq)
    seq = collect(seq)
    complement_seq = []
    for base in eachindex(seq)
        new_seq = complement(seq[base])
        push!(complement_seq, new_seq)
    end
    complement_seq = join(complement_seq)
    return complement_seq
end

function reverse_complement(seq)
    seq = normalizeDNA(seq)
    reverse_sequence = reverse(seq)
    sequence = collect(reverse_sequence)
    reverse_complement = map(complement, sequence)
    reverse_complement = join(reverse_complement)
end

function parse_fasta(path)
    header = String[]
    sequence = []
    seqonly = []
    for line in eachline(path)
        if startswith(line, '>')
            if length(header) > length(sequence)
                seqonlynew = join(seqonly)
                push!(sequence, seqonlynew)
                seqonly = []
            end
            headerline = split(line, '>')
            headerline = headerline[2]
            push!(header, headerline)
        elseif line == ""
            continue
        else
            #if line is not a header or an empty string, push it to seqonly
            #keep pushing these lines until you reach another header or empty string line
            #then join all of these lines together and push the sequence to sequence
            line = normalizeDNA(line)
            push!(seqonly, line)
        end
    end
    seqonlynew = join(seqonly)
    push!(sequence, seqonlynew)
    return(header, sequence)
end

end # module KannanBioinformaticsBISC195


text:
#EDAData: PSF Reader Interface
#-------------------------------------------------------------------------------

#==Main data structures
===============================================================================#
mutable struct PSFReader
	reader::PSFReaderLib.DataReader
	x::Vector
end

#==Open/read/close functions
===============================================================================#
function Base.open(::Type{PSFReader}, path::String)
	local x
	reader = PSFReaderLib._open(path)
	try
		x = PSFReaderLib.readsweep(reader)
	catch
		x = Nothing[] #No sweep
	end
	return PSFReader(reader, x)
end

function Base.read(r::PSFReader, signame::String)
#	if typeof(r.x) <: Vector{Nothing}
		#Very hacky... Figure out something better
#		return PSFReaderLib.readscalar(r.reader, signame)
#	end
	y = read(r.reader, signame)
	return DataF1(r.x, y)
end

Base.close(r::PSFReader) = return #No close interface provided by library

#==Accessors
===============================================================================#
Base.names(r::PSFReader) = names(r.reader)

#==Generate friendly show functions
===============================================================================#
function Base.show(io::IO, r::PSFReader)
	print(io, "PSFReader($(r.reader.filepath))")
end

#==High-level interface
===============================================================================#
open_psf(path::String) = open(PSFReader, path) #Use PSF reader from this module.
open_psf(fn::Function, path::String) = open(fn, PSFReader, path) #For do/end method

#Last line


text:
using Plots

Plots.gr()

x= -π:0.01:π
y1=sin.(x)
y2=cos.(x)
p=plot(
    x,
    hcat(y1,y2),
    title="sin curve",
    label=["sin","cos"],
    lw=3
)

xlabel!(p,"x")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
#__precompile__()

module RcImages
import ..RchannelExperiments: plot_res
using ..RchannelExperiments

using Images, ImageView, ImageSegmentation, ImageMagick
using CoordinateTransformations, OffsetArrays
using ProgressMeter
import PyPlot
import FileIO
import JLD
# const JLD = JLD2
const P = PyPlot
import PyCall
PyCall.@pyimport matplotlib.animation as anim
using StatsBase
using Parameters
import DSP
import Interpolations
const IP=Interpolations

export ExpImgs, ExpImgsResults, channel_width, prep_img, get_time_series

"""
All parameters needed to turn pixels into meters, assuming
an image without distortion.

TODO: this should eventually be in Setup
"""
@with_kw struct Pixel2Meter @deftype Float64
    h1 # distance surface to pupil of camera
    h2 # thickness of water-pool or perspex
    h3 # depth from bottom of pool to r-channel center
    n1 = 1.0
    n2 = [1.339, 1.48][1] # [water, perspex]
    n3 = 1.31 # ice
    ppm::Float64 # pixel (of full resolution) per meter at surface
end

"""
This holds all data needed to process the automatic pictures
from one Experiment.
"""
@with_kw struct ExpImgs @deftype Int
    # Directory holding the pictures
    dir::String
    # time to be subtracted from EXIF:DateTimeOriginal to get correct time
    time_correction::Dates.Period=Dates.Second(0)
    experiment_start::Dates.DateTime # start-time of experiment.  Plotting will be relative to this time
    p1::Tuple{Int,Int}
    p2::Tuple{Int,Int}
    halfheight_crop = 1200
    thin_num = 4
    color_loc::Tuple{Float64,Float64} = (0.5,0.5) # location of where to pic the color (in original image, coords 0..1)
    ns::Vector{Int}
    minhalfwidth_orig  # minimal half-width of R-channel in pixels of original image
    algo::Symbol=[:both, :thresh, :edge][1]
    gauss_w = 1
    quant::Float64=0.6
    gap = 19
    median_filter_region::Tuple{Int,Int} = (3,7) # (top-bottom, left-right)
    p2m::Pixel2Meter
    extras::Dict{Symbol}=Dict{Symbol,Any}()
end
function Base.size(exi::ExpImgs)
    @unpack halfheight_crop, thin_num, p1, p2 = exi
    return fld1(2*halfheight_crop+1,thin_num), fld1(round(Int, sqrt((p2[1]-p1[1])^2 + (p2[2]-p1[2])^2)) +1, thin_num)
end

function image_files(exi::ExpImgs, nss=exi.ns)
    @unpack dir = exi
    ["$dir/$f" for f in readdir(dir) if lowercase(splitext(f)[2])==".jpg"][nss]
end
"""
Hold the result from the processing
"""
@with_kw struct ExpImgsResults
    exi::ExpImgs
    time::Vector{Dates.Time} # time of day of image
    t::Vector{Float64}     # time since start of experiment (s)
    center_dist::Vector{Int} # distance from original image-center to center line of R-channel
    ts::Matrix{Int} # tops
    bs::Matrix{Int} # bottoms
    ts_dist::Matrix{Float64}
    bs_dist::Matrix{Float64}
    # products
    dia_mean::Vector{Float64} # mean diameter
    dia_std::Vector{Float64} # std of diameter
    dia_quant::Matrix{Float64} # 10%, 50%, 90% quantiles
    scalloping::Vector{Float64} # scalloping factor
    tb_cor::Vector{Float64}
    # images
    thumbs::Vector{Any} # thumbnail images
    imgs::Vector{Any} # full images
    extras::Dict{Symbol}=Dict{Symbol,Any}()
end
Base.length(res::ExpImgsResults) = length(res.time)

include("optics.jl")
include("helpers.jl")
include("preparation.jl")
include("width.jl")
include("plotting.jl")
end # module


text:
@system Germination(Stage, Planting, Emergence) begin
    germinateable(planted) ~ flag
    #HACK: assume germination occurs at 50% of emergence (2021-06-14: SHK, KDY)
    germinated(emergence, emerged) => begin
        (emergence >= 0.5) || emerged
    end ~ flag
    germinating(germinateable & !germinated) ~ flag
end


text:
# ProximalOperators.jl - library of nonsmooth functions and associated proximal mappings

__precompile__()

module ProximalOperators

typealias RealOrComplex Union{Real, Complex}
typealias HermOrSym{T,S} Union{Hermitian{T,S}, Symmetric{T,S}}

export prox, prox!

export ProximableFunction

export IndAffine, IndHalfspace,
       IndBallLinf, IndBallL0, IndBallL1, IndBallL2, IndBallRank,
       IndBox, IndNonnegative, IndNonpositive,
       IndExpPrimal, IndExpDual, IndPSD, IndSOC, IndRotatedSOC,
       IndFree,
       IndPoint, IndZero,
       IndSimplex,
       IndSphereL2, IndBinary,
       SumPositive, HingeLoss, HuberLoss,
       LogBarrier,
       LeastSquares,
       Maximum,
       NormL0, NormL1, NormL2, NormL21, NormLinf, NuclearNorm, SqrNormL2, ElasticNet,
       FirmThreshold,
       DistL2, SqrDistL2,
       Zero

export Conjugate,
       Postcompose,
       Precompose,
       SlicedSeparableSum,
       SeparableSum,
       Tilt,
       Regularize

abstract ProximableFunction
abstract ProximableConvex <: ProximableFunction
abstract IndicatorConvex <: ProximableConvex
abstract IndicatorConvexCone <: IndicatorConvex
abstract IndicatorNonconvex <: ProximableFunction

include("utilities/symmetricpacked.jl")

include("calculus/conjugate.jl")
include("calculus/postcompose.jl")
include("calculus/precompose.jl")
include("calculus/separableSum.jl")
include("calculus/slicedSeparableSum.jl")
include("calculus/tilt.jl")
include("calculus/regularize.jl")

include("functions/distL2.jl")
include("functions/elasticNet.jl")
include("functions/firmThreshold.jl")
include("functions/logBarrier.jl")
include("functions/normL2.jl")
include("functions/normL1.jl")
include("functions/normL21.jl")
include("functions/normL0.jl")
include("functions/nuclearNorm.jl")
include("functions/hingeLoss.jl")
include("functions/huberLoss.jl")
include("functions/indAffine.jl")
include("functions/indBallL0.jl")
include("functions/indBallL1.jl")
include("functions/indBallL2.jl")
include("functions/indBallRank.jl")
include("functions/indBinary.jl")
include("functions/indBox.jl")
include("functions/indNonnegative.jl")
include("functions/indNonpositive.jl")
include("functions/indZero.jl")
include("functions/indExp.jl")
include("functions/indFree.jl")
include("functions/indPoint.jl")
include("functions/indPSD.jl")
include("functions/indSimplex.jl")
include("functions/indSOC.jl")
include("functions/indSphereL2.jl")
include("functions/indHalfspace.jl")
include("functions/sqrDistL2.jl")
include("functions/sqrNormL2.jl")
include("functions/sumLargest.jl")
include("functions/sumPositive.jl")
include("functions/maximum.jl")
include("functions/normLinf.jl")
include("functions/leastSquares.jl")

function Base.show(io::IO, f::ProximableFunction)
  println(io, "description : ", fun_name(f))
  println(io, "domain      : ", fun_dom(f))
  println(io, "expression  : ", fun_expr(f))
  print(  io, "parameters  : ", fun_params(f))
end

fun_name(  f) = "n/a"
fun_dom(   f) = "n/a"
fun_expr(  f) = "n/a"
fun_params(f) = "n/a"

is_prox_accurate(f::ProximableFunction) = true

"""
  prox(f::ProximableFunction, x::AbstractArray, γ::Real=1.0)

Computes the proximal point of `x` with respect to function `f`
and parameter `γ > 0`, that is

  y = argmin_z { f(z) + 1/(2γ)||z-x||² }

and returns `y` and `f(y)`.
"""

function prox(f::ProximableFunction, x::AbstractArray, gamma::Real=1.0)
  y = similar(x)
  fy = prox!(f, x, y, gamma)
  return y, fy
end

"""
  prox!(f::ProximableFunction, x::AbstractArray, y::AbstractArray, γ::Real=1.0)

Computes the proximal point of `x` with respect to function `f`
and parameter `γ > 0`, and writes the result in `y`. Returns `f(y)`.
"""

function prox!(f::ProximableFunction, x::AbstractArray, y::AbstractArray, gamma::Real=1.0)
  throw(MethodException(
    "prox! is not implemented for f::", typeof(f),
    ", x::", typeof(x), ", y::", typeof(y), ", gamma::", typeof(gamma)
  ))
end

"""
  prox!(f::ProximableFunction, x::AbstractArray, γ::Real=1.0)

Computes the proximal point of `x` with respect to function `f`
and parameter `γ > 0` *in place*, that is

  x ← argmin_z { f(z) + 1/(2γ)||z-x||² }

and returns `f(x)`.
"""

prox!(f::ProximableFunction, x::AbstractArray, gamma::Real=1.0) = prox!(f, x, x, gamma)

end


text:
__precompile__()

module cme257demo

using Distributions.Exponential

export hello, demofun, exprand

function hello()
  print("hello from cme257demo!")
end

function demofun(x::Int, y::Int)
  x + y
end

function exprand()
  d = Exponential()
  rand(d)
end

end # module


text:
@testset "IO" begin
    @testset "FileRoller" begin
        levels = Dict(
            "not_set" => 0,
            "debug" => 10,
            "info" => 20,
            "warn" => 30,
            "error" => 40,
            "fubar" => 50
        )
        roller_prefix = tempname()
        @info("Path to roller_prefix: $roller_prefix")
        handler = DefaultHandler(FileRoller(roller_prefix; max_sz=1024))

        logger = Logger(
            "IO.FileRoller",
            Dict("Roller" => handler),
            "debug",
            levels,
            DefaultRecord,
            true
        )

        for i = 1:100
            log(logger, "debug", "some-msg")
            log(logger, "info", "some-msg")
            log(logger, "warn", "some-msg")
            log(logger, "error", "some-msg")
            log(logger, "fubar", "some-msg")
        end

        @test ispath(dirname(roller_prefix))
        @test isfile(handler.io.filepath)
        @test dirname(roller_prefix) == dirname(handler.io.filepath)
        @test occursin(roller_prefix, handler.io.filepath)
        @test handler.io.filepath != roller_prefix
        @test handler.io.filepath != "$(roller_prefix).001"
    end
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^-1*arcsec(a*x^5), x) == :(-1//5*arcsec(a*x^5)*log(1+exp(2*im*arcsec(a*x^5)))+1//10*im*arcsec(a*x^5)^2+1//10*im*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a*x^5))))
@test integrate(x^3*arcsec(x^(1/2)), x) == :(-3//20*(-1+x)^5//2+-1//4*(-1+x)^(1/2)+-1//4*(-1+x)^3//2+-1//28*(-1+x)^7//2+1//4*x^4*arcsec(x^(1/2)))
@test integrate(x^2*arcsec(x^(1/2)), x) == :(-2//9*(-1+x)^3//2+-1//3*(-1+x)^(1/2)+-1//15*(-1+x)^5//2+1//3*x^3*arcsec(x^(1/2)))
@test integrate(x*arcsec(x^(1/2)), x) == :(-1//2*(-1+x)^(1/2)+-1//6*(-1+x)^3//2+(1/2)*x^2*arcsec(x^(1/2)))
@test integrate(arcsec(x^(1/2)), x) == :(-1*(-1+x)^(1/2)+x*arcsec(x^(1/2)))
@test integrate(x^-1*arcsec(x^(1/2)), x) == :(im*arcsec(x^(1/2))^2+im*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(x^(1/2))))+-2*arcsec(x^(1/2))*log(1+exp(2*im*arcsec(x^(1/2)))))
@test integrate(x^-2*arcsec(x^(1/2)), x) == :((1/2)*arctan((-1+x)^(1/2))+(1/2)*x^-1*(-1+x)^(1/2)+-1*x^-1*arcsec(x^(1/2)))
@test integrate(x^-3*arcsec(x^(1/2)), x) == :(3//16*arctan((-1+x)^(1/2))+-1//2*x^-2*arcsec(x^(1/2))+1//8*x^-2*(-1+x)^(1/2)+3//16*x^-1*(-1+x)^(1/2))
@test integrate(x^-4*arcsec(x^(1/2)), x) == :(5//48*arctan((-1+x)^(1/2))+-1//3*x^-3*arcsec(x^(1/2))+1//18*x^-3*(-1+x)^(1/2)+5//48*x^-1*(-1+x)^(1/2)+5//72*x^-2*(-1+x)^(1/2))
@test integrate(x^2*arcsec(a*x^-1), x) == :(-1//3*a^3*(1+-1*a^-2*x^2)^(1/2)+1//3*x^3*arccos(x*a^-1)+1//9*a^3*(1+-1*a^-2*x^2)^3//2)
@test integrate(x*arcsec(a*x^-1), x) == :((1/2)*x^2*arccos(x*a^-1)+1//4*a^2*arcsin(x*a^-1)+-1//4*a*x*(1+-1*a^-2*x^2)^(1/2))
@test integrate(arcsec(a*x^-1), x) == :(x*arccos(x*a^-1)+-1*a*(1+-1*a^-2*x^2)^(1/2))
@test integrate(x^-1*arcsec(a*x^-1), x) == :(arccos(x*a^-1)*log(1+exp(2*im*arccos(x*a^-1)))+-1//2*im*arccos(x*a^-1)^2+-1//2*im*Polylogarithms.polylog(2,-1*exp(2*im*arccos(x*a^-1))))
@test integrate(x^-2*arcsec(a*x^-1), x) == :(a^-1*arctanh((1+-1*a^-2*x^2)^(1/2))+-1*x^-1*arccos(x*a^-1))
@test integrate(x^-3*arcsec(a*x^-1), x) == :(-1//2*x^-2*arccos(x*a^-1)+(1/2)*a^-1*x^-1*(1+-1*a^-2*x^2)^(1/2))
@test integrate(x^-4*arcsec(a*x^-1), x) == :(-1//3*x^-3*arccos(x*a^-1)+1//6*a^-3*arctanh((1+-1*a^-2*x^2)^(1/2))+1//6*a^-1*x^-2*(1+-1*a^-2*x^2)^(1/2))
@test integrate(x^-1*arcsec(a*x^n), x) == :((1/2)*im*n^-1*arcsec(a*x^n)^2+(1/2)*im*n^-1*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a*x^n)))+-1*n^-1*arcsec(a*x^n)*log(1+exp(2*im*arcsec(a*x^n))))
@test integrate(x^4*arcsec(a+b*x), x) == :(1//5*x^5*arcsec(a+b*x)+-1//40*b^-5*(3+40*a^2+40*a^4)*arctanh((1+-1*(a+b*x)^-2)^(1/2))+1//5*a^5*b^-5*arcsec(a+b*x)+-1//20*b^-2*x^3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+-1//120*b^-5*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*(9+58*a^2)+1//30*a*b^-5*(1+-1*(a+b*x)^-2)^(1/2)*(20+53*a^2)*(a+b*x)+11//60*a*b^-3*x^2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x^3*arcsec(a+b*x), x) == :(1//4*x^4*arcsec(a+b*x)+-1//4*a^4*b^-4*arcsec(a+b*x)+(1/2)*a*b^-4*(1+2*a^2)*arctanh((1+-1*(a+b*x)^-2)^(1/2))+-1//12*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(2+17*a^2)*(a+b*x)+-1//12*b^-2*x^2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+1//3*a*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2)
@test integrate(x^2*arcsec(a+b*x), x) == :(1//3*x^3*arcsec(a+b*x)+-1//6*b^-3*(1+6*a^2)*arctanh((1+-1*(a+b*x)^-2)^(1/2))+1//3*a^3*b^-3*arcsec(a+b*x)+-1//6*x*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+5//6*a*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x*arcsec(a+b*x), x) == :((1/2)*x^2*arcsec(a+b*x)+a*b^-2*arctanh((1+-1*(a+b*x)^-2)^(1/2))+-1//2*a^2*b^-2*arcsec(a+b*x)+-1//2*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(arcsec(a+b*x), x) == :(-1*b^-1*arctanh((1+-1*(a+b*x)^-2)^(1/2))+b^-1*(a+b*x)*arcsec(a+b*x))
@test integrate(x^-1*arcsec(a+b*x), x) == :(arcsec(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+arcsec(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+(1/2)*im*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a+b*x)))+-1*im*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-1*im*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-1*arcsec(a+b*x)*log(1+exp(2*im*arcsec(a+b*x))))
@test integrate(x^-2*arcsec(a+b*x), x) == :(-1*x^-1*arcsec(a+b*x)+-1*b*a^-1*arcsec(a+b*x)+2*b*a^-1*(1+-1*a^2)^-1//2*arctan((1+a)^(1/2)*(1+-1a)^-1//2*tan((1/2)*arcsec(a+b*x))))
@test integrate(x^-3*arcsec(a+b*x), x) == :(-1//2*x^-2*arcsec(a+b*x)+(1/2)*a^-2*b^2*arcsec(a+b*x)+-1*a^-2*b^2*(1+-1*a^2)^-3//2*(1+-2*a^2)*arctan((1+a)^(1/2)*(1+-1a)^-1//2*tan((1/2)*arcsec(a+b*x)))+(1/2)*b*a^-1*x^-1*(1+-1*a^2)^-1*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x^-4*arcsec(a+b*x), x) == :(-1//3*x^-3*arcsec(a+b*x)+-1//3*a^-3*b^3*arcsec(a+b*x)+1//3*a^-3*b^3*(1+-1*a^2)^-5//2*(2+-5*a^2+6*a^4)*arctan((1+a)^(1/2)*(1+-1a)^-1//2*tan((1/2)*arcsec(a+b*x)))+1//6*b*a^-1*x^-2*(1+-1*a^2)^-1*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+-1//6*a^-2*b^2*x^-1*(1+-1*a^2)^-2*(1+-1*(a+b*x)^-2)^(1/2)*(2+-5*a^2)*(a+b*x))
@test integrate(x^3*arcsec(a+b*x)^2, x) == :(1//3*b^-4*log(a+b*x)+1//4*x^4*arcsec(a+b*x)^2+1//12*b^-4*(a+b*x)^2+-1*a*x*b^-3+3*a^2*b^-4*log(a+b*x)+-1//4*a^4*b^-4*arcsec(a+b*x)^2+im*a*b^-4*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))+-1*im*a*b^-4*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))+-2*im*a^3*b^-4*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))+2*im*a^3*b^-4*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))+-1//3*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arcsec(a+b*x)+-1//6*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^3*arcsec(a+b*x)+a*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*arcsec(a+b*x)+-4*im*a^3*b^-4*arcsec(a+b*x)*arctan(exp(im*arcsec(a+b*x)))+-3*a^2*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arcsec(a+b*x)+-2*im*a*b^-4*arcsec(a+b*x)*arctan(exp(im*arcsec(a+b*x))))
@test integrate(x^2*arcsec(a+b*x)^2, x) == :(1//3*x*b^-2+1//3*x^3*arcsec(a+b*x)^2+-2*a*b^-3*log(a+b*x)+-1//3*im*b^-3*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))+1//3*im*b^-3*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))+1//3*a^3*b^-3*arcsec(a+b*x)^2+-2*im*a^2*b^-3*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))+2*im*a^2*b^-3*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))+-1//3*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*arcsec(a+b*x)+2//3*im*b^-3*arcsec(a+b*x)*arctan(exp(im*arcsec(a+b*x)))+2*a*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arcsec(a+b*x)+4*im*a^2*b^-3*arcsec(a+b*x)*arctan(exp(im*arcsec(a+b*x))))
@test integrate(x*arcsec(a+b*x)^2, x) == :(b^-2*log(a+b*x)+(1/2)*x^2*arcsec(a+b*x)^2+-1//2*a^2*b^-2*arcsec(a+b*x)^2+-1*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arcsec(a+b*x)+-2*im*a*b^-2*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))+2*im*a*b^-2*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))+-4*im*a*b^-2*arcsec(a+b*x)*arctan(exp(im*arcsec(a+b*x))))
@test integrate(arcsec(a+b*x)^2, x) == :(b^-1*arcsec(a+b*x)^2*(a+b*x)+-2*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))+2*im*b^-1*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))+4*im*b^-1*arcsec(a+b*x)*arctan(exp(im*arcsec(a+b*x))))
@test integrate(x^-1*arcsec(a+b*x)^2, x) == :(2*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+2*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-1//2*Polylogarithms.polylog(3,-1*exp(2*im*arcsec(a+b*x)))+arcsec(a+b*x)^2*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+arcsec(a+b*x)^2*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-1*arcsec(a+b*x)^2*log(1+exp(2*im*arcsec(a+b*x)))+im*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a+b*x)))*arcsec(a+b*x)+-2*im*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+-2*im*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))*arcsec(a+b*x))
@test integrate(x^-2*arcsec(a+b*x)^2, x) == :(-1*x^-1*arcsec(a+b*x)^2+-1*b*a^-1*arcsec(a+b*x)^2+-2*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+2*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-2*im*b*a^-1*(1+-1*a^2)^-1//2*arcsec(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+2*im*b*a^-1*(1+-1*a^2)^-1//2*arcsec(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x))))
@test integrate(x^2*arcsec(a+b*x)^3, x) == :(b^-3*Polylogarithms.polylog(3,-1*im*exp(im*arcsec(a+b*x)))+-1*b^-3*Polylogarithms.polylog(3,im*exp(im*arcsec(a+b*x)))+-1*b^-3*arctanh((1+-1*(a+b*x)^-2)^(1/2))+1//3*x^3*arcsec(a+b*x)^3+b^-3*(a+b*x)*arcsec(a+b*x)+-6*a^2*b^-3*Polylogarithms.polylog(3,im*exp(im*arcsec(a+b*x)))+6*a^2*b^-3*Polylogarithms.polylog(3,-1*im*exp(im*arcsec(a+b*x)))+1//3*a^3*b^-3*arcsec(a+b*x)^3+im*b^-3*arcsec(a+b*x)^2*arctan(exp(im*arcsec(a+b*x)))+im*b^-3*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+-1*im*b^-3*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+-3*im*a*b^-3*arcsec(a+b*x)^2+-3*im*a*b^-3*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a+b*x)))+6*a*b^-3*arcsec(a+b*x)*log(1+exp(2*im*arcsec(a+b*x)))+-1//2*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*arcsec(a+b*x)^2+-6*im*a^2*b^-3*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+3*a*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*arcsec(a+b*x)^2*(a+b*x)+6*im*a^2*b^-3*arcsec(a+b*x)^2*arctan(exp(im*arcsec(a+b*x)))+6*im*a^2*b^-3*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x))
@test integrate(x*arcsec(a+b*x)^3, x) == :((1/2)*x^2*arcsec(a+b*x)^3+-6*a*b^-2*Polylogarithms.polylog(3,-1*im*exp(im*arcsec(a+b*x)))+-3*b^-2*arcsec(a+b*x)*log(1+exp(2*im*arcsec(a+b*x)))+6*a*b^-2*Polylogarithms.polylog(3,im*exp(im*arcsec(a+b*x)))+-1//2*a^2*b^-2*arcsec(a+b*x)^3+3//2*im*b^-2*arcsec(a+b*x)^2+3//2*im*b^-2*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a+b*x)))+-1//2*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*arcsec(a+b*x)^2*(3a+3*b*x)+-6*im*a*b^-2*arcsec(a+b*x)^2*arctan(exp(im*arcsec(a+b*x)))+-6*im*a*b^-2*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+6*im*a*b^-2*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x))
@test integrate(arcsec(a+b*x)^3, x) == :(-6*b^-1*Polylogarithms.polylog(3,im*exp(im*arcsec(a+b*x)))+6*b^-1*Polylogarithms.polylog(3,-1*im*exp(im*arcsec(a+b*x)))+b^-1*arcsec(a+b*x)^3*(a+b*x)+-6*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+6*im*b^-1*arcsec(a+b*x)^2*arctan(exp(im*arcsec(a+b*x)))+6*im*b^-1*Polylogarithms.polylog(2,im*exp(im*arcsec(a+b*x)))*arcsec(a+b*x))
@test integrate(x^-1*arcsec(a+b*x)^3, x) == :(arcsec(a+b*x)^3*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+arcsec(a+b*x)^3*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-1*arcsec(a+b*x)^3*log(1+exp(2*im*arcsec(a+b*x)))+6*im*Polylogarithms.polylog(4,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+6*im*Polylogarithms.polylog(4,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+6*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+6*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+-3//2*Polylogarithms.polylog(3,-1*exp(2*im*arcsec(a+b*x)))*arcsec(a+b*x)+-3//4*im*Polylogarithms.polylog(4,-1*exp(2*im*arcsec(a+b*x)))+-3*im*arcsec(a+b*x)^2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-3*im*arcsec(a+b*x)^2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+3//2*im*arcsec(a+b*x)^2*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a+b*x))))
@test integrate(x^-2*arcsec(a+b*x)^3, x) == :(-1*x^-1*arcsec(a+b*x)^3+-1*b*a^-1*arcsec(a+b*x)^3+-6*im*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+-6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+6*im*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))*arcsec(a+b*x)+-3*im*b*a^-1*(1+-1*a^2)^-1//2*arcsec(a+b*x)^2*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x)))+3*im*b*a^-1*(1+-1*a^2)^-1//2*arcsec(a+b*x)^2*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arcsec(a+b*x))))
@test integrate(x*(a+b*arcsec(c+d*x^2)), x) == :((1/2)*a*x^2+-1//2*b*d^-1*arctanh((1+-1*(c+d*x^2)^-2)^(1/2))+(1/2)*b*d^-1*(c+d*x^2)*arcsec(c+d*x^2))
@test integrate(x^2*(a+b*arcsec(c+d*x^3)), x) == :(1//3*a*x^3+-1//3*b*d^-1*arctanh((1+-1*(c+d*x^3)^-2)^(1/2))+1//3*b*d^-1*(c+d*x^3)*arcsec(c+d*x^3))
@test integrate(x^3*(a+b*arcsec(c+d*x^4)), x) == :(1//4*a*x^4+-1//4*b*d^-1*arctanh((1+-1*(c+d*x^4)^-2)^(1/2))+1//4*b*d^-1*(c+d*x^4)*arcsec(c+d*x^4))
@test integrate(x^(-1+n)*arcsec(a+b*x^n), x) == :(-1*b^-1*n^-1*arctanh((1+-1*(a+b*x^n)^-2)^(1/2))+b^-1*n^-1*(a+b*x^n)*arcsec(a+b*x^n))
@test integrate(arcsec(c*exp(a+b*x)), x) == :((1/2)*im*b^-1*arcsec(c*exp(a+b*x))^2+(1/2)*im*b^-1*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(c*exp(a+b*x))))+-1*b^-1*arcsec(c*exp(a+b*x))*log(1+exp(2*im*arcsec(c*exp(a+b*x)))))
@test integrate(x^-2*exp(arcsec(a*x)), x) == :(-1//2*x^-1*exp(arcsec(a*x))+(1/2)*a*(1+-1*a^-2*x^-2)^(1/2)*exp(arcsec(a*x)))
@test integrate(x^-3*exp(arcsec(a*x)), x) == :(-1//5*a^2*cos(2*arcsec(a*x))*exp(arcsec(a*x))+1//10*a^2*exp(arcsec(a*x))*sin(2*arcsec(a*x)))
@test integrate(x^-4*exp(arcsec(a*x)), x) == :(-3//40*a^3*cos(3*arcsec(a*x))*exp(arcsec(a*x))+-1//8*a^2*x^-1*exp(arcsec(a*x))+1//8*a^3*(1+-1*a^-2*x^-2)^(1/2)*exp(arcsec(a*x))+1//40*a^3*exp(arcsec(a*x))*sin(3*arcsec(a*x)))
@test integrate((d*x+a*d*b^-1)^-1*arcsec(a+b*x), x) == :((1/2)*im*d^-1*arcsec(a+b*x)^2+(1/2)*im*d^-1*Polylogarithms.polylog(2,-1*exp(2*im*arcsec(a+b*x)))+-1*d^-1*arcsec(a+b*x)*log(1+exp(2*im*arcsec(a+b*x))))

max_stars_repo_path: test/5 Inverse trig functions/5.5 Inverse secant/5.5.2 Inverse secant functions.jl


text:
using Gaugefields
using Wilsonloop
using LinearAlgebra

function calculate_topological_charge_plaq(U::Array{T,1},temp_UμνTA,temps) where T
    UμνTA = temp_UμνTA
    numofloops = calc_UμνTA!(UμνTA,"plaq",U,temps)
    Q = calc_Q(UμνTA,numofloops,U)
    return Q
end

function calculate_topological_charge_clover(U::Array{T,1},temp_UμνTA,temps) where T
    UμνTA = temp_UμνTA
    numofloops = calc_UμνTA!(UμνTA,"clover",U,temps)
    Q = calc_Q(UμνTA,numofloops,U)
    return Q
end

function calculate_topological_charge_improved(U::Array{T,1},temp_UμνTA,Qclover,temps) where T
    UμνTA = temp_UμνTA
    numofloops = calc_UμνTA!(UμνTA,"rect",U,temps)
    Qrect = 2*calc_Q(UμνTA,numofloops,U)
    c1 = -1/12
    c0 = 5/3
    Q = c0*Qclover + c1*Qrect
    return Q
end

function calc_Q(UμνTA,numofloops,U::Array{<: AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
    Q = 0.0
    if Dim == 4
        ε(μ,ν,ρ,σ) = epsilon_tensor(μ,ν,ρ,σ)
    else
        error("Dimension $Dim is not supported")
    end
    for μ=1:Dim
        for ν=1:Dim
            if ν == μ
                continue
            end
            Uμν = UμνTA[μ,ν]
            for ρ =1:Dim
                for σ=1:Dim
                    if ρ == σ
                        continue
                    end
                    Uρσ = UμνTA[ρ,σ]
                    s = tr(Uμν,Uρσ)
                    Q += ε(μ,ν,ρ,σ)*s/numofloops^2
                end
            end
        end
    end

    return -real(Q)/(32*(π^2))
end

function calc_UμνTA!(temp_UμνTA,name::String,U::Array{<: AbstractGaugefields{NC,Dim},1},temps) where {NC,Dim}
    loops_μν,numofloops = calc_loopset_μν_name(name,Dim)
    calc_UμνTA!(temp_UμνTA,loops_μν,U,temps)
    return numofloops
end

function calc_UμνTA!(UμνTA,loops_μν,U::Array{<: AbstractGaugefields{NC,Dim},1},temps) where {NC,Dim}
    for μ=1:Dim
        for ν=1:Dim
            if ν == μ
                continue
            end
            evaluate_gaugelinks!(temps[1],loops_μν[μ,ν],U,temps[2:3])
            Traceless_antihermitian!(UμνTA[μ,ν],temps[1])
        end
    end
    return
end

function calc_loopset_μν_name(name,Dim)
    loops_μν= Array{Vector{Wilsonline{Dim}},2}(undef,Dim,Dim)
    if name == "plaq"
        numofloops = 1
        for μ=1:Dim
            for ν=1:Dim
                loops_μν[μ,ν] = Wilsonline{Dim}[]
                if ν == μ
                    continue
                end
                plaq = make_plaq(μ,ν,Dim=Dim)
                push!(loops_μν[μ,ν],plaq)
            end
        end
    elseif name == "clover"
        numofloops = 4
        for μ=1:Dim
            for ν=1:Dim
                loops_μν[μ,ν] = Wilsonline{Dim}[]
                if ν == μ
                    continue
                end
                loops_μν[μ,ν] = make_cloverloops(μ,ν,Dim=Dim)
            end
        end
    elseif name == "rect"
        numofloops = 8
        for μ=1:4
            for ν=1:4
                if ν == μ
                    continue
                end
                loops = Wilsonline{Dim}[]
                loop_righttop = Wilsonline([(μ,2),(ν,1),(μ,-2),(ν,-1)])
                loop_lefttop = Wilsonline([(ν,1),(μ,-2),(ν,-1),(μ,2)])
                loop_rightbottom = Wilsonline([(ν,-1),(μ,2),(ν,1),(μ,-2)])
                loop_leftbottom= Wilsonline([(μ,-2),(ν,-1),(μ,2),(ν,1)])
                push!(loops,loop_righttop)
                push!(loops,loop_lefttop)
                push!(loops,loop_rightbottom)
                push!(loops,loop_leftbottom)

                loop_righttop = Wilsonline([(μ,1),(ν,2),(μ,-1),(ν,-2)])
                loop_lefttop = Wilsonline([(ν,2),(μ,-1),(ν,-2),(μ,1)])
                loop_rightbottom = Wilsonline([(ν,-2),(μ,1),(ν,2),(μ,-1)])
                loop_leftbottom= Wilsonline([(μ,-1),(ν,-2),(μ,1),(ν,2)])
                push!(loops,loop_righttop)
                push!(loops,loop_lefttop)
                push!(loops,loop_rightbottom)
                push!(loops,loop_leftbottom)

                loops_μν[μ,ν] = loops
            end
        end
    else
        error("$name is not supported")
    end
    return loops_μν,numofloops
end

function make_cloverloops(μ,ν;Dim=4)
    loops = Wilsonline{Dim}[]
    loop_righttop = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)])
    loop_lefttop = Wilsonline([(ν,1),(μ,-1),(ν,-1),(μ,1)])
    loop_rightbottom = Wilsonline([(ν,-1),(μ,1),(ν,1),(μ,-1)])
    loop_leftbottom= Wilsonline([(μ,-1),(ν,-1),(μ,1),(ν,1)])
    push!(loops,loop_righttop)
    push!(loops,loop_lefttop)
    push!(loops,loop_rightbottom)
    push!(loops,loop_leftbottom)
    return loops
end

#topological charge
function epsilon_tensor(mu::Int,nu::Int,rho::Int,sigma::Int)
    sign=1 # (3) 1710.09474 extended epsilon tensor
    if mu < 0
        sign*=-1
        mu=-mu
    end
    if nu < 0
        sign*=-1
        nu=-nu
    end
    if rho < 0
        sign*=-1
        rh=-rho
    end
    if sigma < 0
        sign*=-1
        sigma=-sigma
    end
    epsilon = zeros(Int,4,4,4,4)
    epsilon[ 1, 2, 3, 4 ] = 1
    epsilon[ 1, 2, 4, 3 ] = -1
    epsilon[ 1, 3, 2, 4 ] = -1
    epsilon[ 1, 3, 4, 2 ] = 1
    epsilon[ 1, 4, 2, 3 ] = 1
    epsilon[ 1, 4, 3, 2 ] = -1
    epsilon[ 2, 1, 3, 4 ] = -1
    epsilon[ 2, 1, 4, 3 ] = 1
    epsilon[ 2, 3, 1, 4 ] = 1
    epsilon[ 2, 3, 4, 1 ] = -1
    epsilon[ 2, 4, 1, 3 ] = -1
    epsilon[ 2, 4, 3, 1 ] = 1
    epsilon[ 3, 1, 2, 4 ] = 1
    epsilon[ 3, 1, 4, 2 ] = -1
    epsilon[ 3, 2, 1, 4 ] = -1
    epsilon[ 3, 2, 4, 1 ] = 1
    epsilon[ 3, 4, 1, 2 ] = 1
    epsilon[ 3, 4, 2, 1 ] = -1
    epsilon[ 4, 1, 2, 3 ] = -1
    epsilon[ 4, 1, 3, 2 ] = 1
    epsilon[ 4, 2, 1, 3 ] = 1
    epsilon[ 4, 2, 3, 1 ] = -1
    epsilon[ 4, 3, 1, 2 ] = -1
    epsilon[ 4, 3, 2, 1 ] = 1
    return epsilon[mu,nu,rho,sigma]*sign
end

function test()
    NX = 8
    NY = 8
    NZ = 8
    NT = 8
    Nwing = 0
    NC = 3

    Dim = 4

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")
    filename="./conf_08080808.ildg"

    ildg = ILDG(filename)
    i = 1
    L = [NX,NY,NZ,NT]
    load_gaugefield!(U,i,ildg,L,NC)

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    β = 5.7
    temp_UμνTA= Matrix{typeof(U[1])}(undef,Dim,Dim)
    for μ=1:Dim
        for ν=1:Dim
            temp_UμνTA[μ,ν] = similar(U[1])
        end
    end

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println(" plaq_t = $plaq_t")

    dt = 0.01

    g = Gradientflow(U,eps = dt)

    for itrj=1:10
        flow!(U,g)
        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        Qplaq = calculate_topological_charge_plaq(U,temp_UμνTA,[temp1,temp2,temp3])
        Qclover = calculate_topological_charge_clover(U,temp_UμνTA,[temp1,temp2,temp3])
        Qimproved= calculate_topological_charge_improved(U,temp_UμνTA,Qclover,[temp1,temp2,temp3])
        println("$itrj $dt $plaq_t $Qplaq $Qclover $Qimproved # itrj dt plaq Qplaq Qclover Qimproved ")
    end

end
test()


text:
function Fract_KelvinVoigt_G_reduce()
    dt = 0.001
    # fractional modulus is NaN not inf at 0.0 so start at dt
    t = collect(dt:dt:10.0)

    fm_model = Fract_KelvinVoigt.Ga(t, [1.0, 1.0, 1.0, 0.0])
    kelvinvoigt = KelvinVoigt.Ga(t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i]), eachindex(t))
end
@test Fract_KelvinVoigt_G_reduce()

function FractS_KelvinVoigt_G_reduce()
    dt = 0.001
    # fractional modulus is NaN not inf at 0.0 so start at dt
    t = collect(dt:dt:10.0)

    fms_model = FractS_KelvinVoigt.Ga(t, [1.0, 1.0, 1.0])
    kelvinvoigt = KelvinVoigt.Ga(t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i]), eachindex(t))
end
@test FractS_KelvinVoigt_G_reduce()

function FractD_KelvinVoigt_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = FractD_KelvinVoigt.Ga(t, [1.0, 1.0, 0.0])
    kelvinvoigt = KelvinVoigt.Ga(t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i]), eachindex(t))
end
@test FractD_KelvinVoigt_G_reduce()

function Fract_KelvinVoigt_J_reduce(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fm_model = Fract_KelvinVoigt.Ja(t, [1.0, 0.99999, 1.0, 0.00001])
    kelvinvoigt = KelvinVoigt.Ja(t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i], atol=tol), eachindex(t))
end
@test Fract_KelvinVoigt_J_reduce(tol)

function FractS_KelvinVoigt_J_reduce(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fms_model = FractS_KelvinVoigt.Ja(t, [1.0, 0.99999, 1.0])
    kelvinvoigt = KelvinVoigt.Ja(t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i], atol=tol), eachindex(t))
end
@test FractS_KelvinVoigt_J_reduce(tol)

function FractD_KelvinVoigt_J_reduce(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = FractD_KelvinVoigt.Ja(t, [1.0, 1.0, 0.00001])
    kelvinvoigt = KelvinVoigt.Ja(t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i], atol=tol), eachindex(t))
end
@test FractD_KelvinVoigt_J_reduce(tol)

function Fract_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = Fract_KelvinVoigt.Gpa(chirp, [1.0, 1.0, 1.0, 0.0])
    kelvinvoigt = KelvinVoigt.Gpa(chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test Fract_KelvinVoigt_Gp_reduce()

function FractS_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = FractS_KelvinVoigt.Gpa(chirp, [1.0, 1.0, 1.0])
    kelvinvoigt = KelvinVoigt.Gpa(chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractS_KelvinVoigt_Gp_reduce()

function FractD_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = FractD_KelvinVoigt.Gpa(chirp, [1.0, 1.0, 0.0])
    kelvinvoigt = KelvinVoigt.Gpa(chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractD_KelvinVoigt_Gp_reduce()

function Fract_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = Fract_KelvinVoigt.Gppa(chirp, [1.0, 1.0, 1.0, 0.0])
    kelvinvoigt = KelvinVoigt.Gppa(chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test Fract_KelvinVoigt_Gp_reduce()

function FractS_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = FractS_KelvinVoigt.Gppa(chirp, [1.0, 1.0, 1.0])
    kelvinvoigt = KelvinVoigt.Gppa(chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractS_KelvinVoigt_Gp_reduce()

function FractD_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = FractD_KelvinVoigt.Gppa(chirp, [1.0, 1.0, 0.0])
    kelvinvoigt = KelvinVoigt.Gppa(chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractD_KelvinVoigt_Gp_reduce()


text:
# Autogenerated wrapper script for raxml_ng_jll for x86_64-linux-gnu
export raxml_ng

JLLWrappers.@generate_wrapper_header("raxml_ng")
JLLWrappers.@declare_executable_product(raxml_ng)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        raxml_ng,
        "bin/raxml-ng",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SIIPExamples
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations

using PowerSystemCaseBuilder
using Dates
using DataFrames
using TimeSeries

using Cbc #solver

include(
    joinpath(SIIPExamples.TEST_DIR, SIIPExamples.PSY_EX_FOLDER, "04_parse_tabulardata.jl"),
)

MultiDay = collect(
    DateTime("2020-01-01T00:00:00"):Hour(1):(DateTime("2020-01-01T00:00:00") + Hour(8783)),
);

remove_time_series!(sys, DeterministicSingleTimeSeries)

for gen in get_components(ThermalGen, sys)
    varcost = get_operation_cost(gen)
    data = TimeArray(MultiDay, repeat([get_cost(get_variable(varcost))], 8784))
    _time_series = SingleTimeSeries("variable_cost", data)
    add_time_series!(sys, gen, _time_series)
    #set_variable_cost!(sys, gen, _time_series)
end

horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval)

uc_template = template_unit_commitment(network = CopperPlatePowerModel)

set_device_model!(uc_template, ThermalMultiStart, ThermalMultiStartUnitCommitment)

solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

problem = OperationsProblem(
    uc_template,
    sys,
    horizon = 4,
    optimizer = solver,
    balance_slack_variables = true,
)
build!(problem, output_dir = mktempdir())

solve!(problem)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#=
utilities for Pardiso
only works with Float64

TODO add to docs following example:
```julia
import Hypatia

ENV["OMP_NUM_THREADS"] = length(Sys.cpu_info())
import Pardiso

syssolver = Hypatia.Solvers.NaiveElimSparseSystemSolver{Float64}(fact_cache = Hypatia.PardisoNonSymCache())
solver = Hypatia.Solvers.Solver{Float64}(syssolver)

include("test/native.jl")
nonnegative1(Float64, solver = solver)
```
=#

mutable struct PardisoNonSymCache{Float64} <: SparseNonSymCache{Float64}
    analyzed::Bool
    pardiso::Pardiso.MKLPardisoSolver
    function PardisoNonSymCache{Float64}()
        cache = new{Float64}()
        cache.analyzed = false
        cache.pardiso = Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(cache.pardiso, Pardiso.REAL_NONSYM)
        return cache
    end
end

mutable struct PardisoSymCache{Float64} <: SparseSymCache{Float64}
    analyzed::Bool
    pardiso::Pardiso.MKLPardisoSolver
    function PardisoSymCache{Float64}()
        cache = new{Float64}()
        cache.analyzed = false
        cache.pardiso = Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(cache.pardiso, Pardiso.REAL_SYM_INDEF)
        return cache
    end
end

PardisoSparseCache = Union{PardisoSymCache{Float64}, PardisoNonSymCache{Float64}}
int_type(::PardisoSparseCache) = Int32

function update_fact(
    cache::PardisoSparseCache{Float64},
    A::SparseMatrixCSC{Float64, Int32},
    )
    pardiso = cache.pardiso

    if !cache.analyzed
        Pardiso.pardisoinit(pardiso)
        # solve transposed problem (Pardiso accepts CSR matrices)
        Pardiso.fix_iparm!(pardiso, :N)
        Pardiso.set_phase!(pardiso, Pardiso.ANALYSIS)
        Pardiso.pardiso(pardiso, A, Float64[])
        cache.analyzed = true
    end

    Pardiso.set_phase!(pardiso, Pardiso.NUM_FACT)
    Pardiso.pardiso(pardiso, A, Float64[])

    return
end

function inv_prod(
    cache::PardisoSparseCache{Float64},
    x::Vector{Float64},
    A::SparseMatrixCSC{Float64, Int32},
    b::Vector{Float64},
    )
    pardiso = cache.pardiso

    Pardiso.set_phase!(pardiso, Pardiso.SOLVE_ITERATIVE_REFINE)
    Pardiso.pardiso(pardiso, x, A, b)

    return x
end

function free_memory(cache::PardisoSparseCache{Float64})
    Pardiso.set_phase!(cache.pardiso, Pardiso.RELEASE_ALL)
    Pardiso.pardiso(cache.pardiso)
    return
end


text:
using Test
using CliOptions

include("testutils.jl")

@testset "ParseResult" begin
    let result = CliOptions.ParseResult(Dict("foo" => "FOO", "bar" => "BAR"))

        @testset "show(::ParseResult)" begin
            @test repr(result) == "CliOptions.ParseResult(:bar,:foo)"
        end

        @testset "getindex; Symbol" begin
            @test result[:foo] == "FOO"
            @test result[:bar] == "BAR"
        end

        @testset "getindex; String" begin
            @test result["foo"] == "FOO"
            @test result["bar"] == "BAR"
        end

        @testset "propertynames" begin
            props = propertynames(result)
            @test sort(props) == [:bar, :foo]
            props = propertynames(result; private = true)
            @test sort(props) == [:_argvals, :_defaults, :_errors, :_resolved, :bar, :foo]
        end

        @testset "getproperty" begin
            @test getproperty(result, :_defaults) isa Dict
            @test getproperty(result, :_argvals) isa Dict
            @test getproperty(result, :_resolved) isa Dict
            @test getproperty(result, :_errors) isa Vector
            @test getproperty(result, :foo) == "FOO"
            @test getproperty(result, :bar) == "BAR"
        end
    end

    @testset "update_defaults(); \"$(v[1])\", $(v[2])" for v in [
        (Dict{String,Any}("n" => 9, "c" => 7, "f" => true, "p" => "bar"),
            "",
            (9, 7, true, "bar")),
        (Dict{String,Any}("n" => 9, "c" => 7, "f" => false, "p" => "bar"),
            "-n 4 -c -f foo",
            (4, 8, true, "foo")),
        (Dict{String,Any}(),
            "-n 4 -c -f foo",
            (4, 1, true, "foo")),
        (Dict{String,Any}(),
            "",
            (0, 0, false, "qux")),
    ]
        defaults, args, expected = v
        spec = CliOptionSpec(
            Option(Int8, "-n"; default = 0),
            CounterOption("-c"),
            FlagOption("-f"),
            Positional("p", default = "qux"),
            onerror = error,
        )
        options = parse_args(spec, split(args))
        options = update_defaults(options, defaults)
        @test options.n == expected[1]
        @test options.c == expected[2]
        @test options.f == expected[3]
        @test options.p == expected[4]
    end
end


text:
module BatchRenaming

export rename!

function match_order(ary, rg)
    medias = filter(contains(rg), ary)
    media_nums = map(x->parse(Int, match(rg, x).captures[1]), medias)
    order = sortperm(media_nums)
    medias[order]
end

function rename!(r1, r2; presuffix="", dryrun=true)
    ary = readdir()

    ary1 = match_order(ary, r1)
    ary2 = match_order(ary, r2)

    @assert length(ary1) == length(ary2)

    N1 = maximum(length, ary1)
    N2 = maximum(length, ary2)
    suf = string(presuffix, splitext(first(ary2))[end])
    if dryrun
        @show suf
        println("File matching:")
        for i in eachindex(ary1)
            println(rpad(ary1[i], N1," "), " <---- ", ary2[i])
        end
        println("Renaming:")
        for i in eachindex(ary1)
            new_name = string(first(splitext(ary1[i])), suf)
            println(rpad(ary2[i], N2," "), " ----> ", new_name)
        end
    else
        for i in eachindex(ary1)
            new_name = string(first(splitext(ary1[i])), suf)
            mv(ary2[i],new_name; force=false)
        end
    end
end

end


text:
using Base.Broadcast
import Base.Broadcast: broadcast!, _broadcast!, broadcast_t
using Base.Broadcast: map_newindexer
using Base: @propagate_inbounds, @pure, RefValue

@inline function const_kernel(state, A, op, len)
    idx = linear_index(state)
    @inbounds if idx <= len
        A[idx] = op()
    end
    return
end
@inline function const_kernel2(state, A, x, len)
    idx = linear_index(state)
    @inbounds if idx <= len
        A[idx] = x
    end
    return
end

function broadcast!(f, A::GPUArray)
    gpu_call(const_kernel, A, (A, f, UInt32(length(A))))
    A
end
function broadcast!(f::typeof(identity), A::GPUArray, val::Number)
    valconv = convert(eltype(A), val)
    gpu_call(const_kernel2, A, (A, valconv, UInt32(length(A))))
    A
end
@inline function broadcast_t(f, T::Type{Bool}, shape, it, A::GPUArrays.GPUArray, Bs::Vararg{Any,N}) where N
    C = similar(A, T, shape)
    keeps, Idefaults = map_newindexer(shape, A, Bs)
    _broadcast!(f, C, keeps, Idefaults, A, Bs, Val{N}, it)
    return C
end
@inline function broadcast_t(f, T::Type{Bool}, shape, it, A::GPUArrays.GPUArray, B::GPUArrays.GPUArray, Bs::Vararg{Any,N}) where N
    C = similar(A, T, shape)
    Bs = (B, Bs...)
    keeps, Idefaults = map_newindexer(shape, A, Bs)
    _broadcast!(f, C, keeps, Idefaults, A, Bs, Val{N}, it)
    return C
end

@inline function broadcast_t(
        f, ::Type{T}, shape, iter, A::GPUArray, Bs::Vararg{Any,N}
    ) where {N, T}
    C = similar(A, T, shape)
    keeps, Idefaults = map_newindexer(shape, A, Bs)
    _broadcast!(f, C, keeps, Idefaults, A, Bs, Val{N}, iter)
    return C
end
@inline function broadcast_t(
        f, ::Type{T}, shape, iter, A::GPUArray, B::GPUArray, rest::Vararg{Any,N}
    ) where {N, T}
    C = similar(A, T, shape)
    Bs = (B, rest...)
    keeps, Idefaults = map_newindexer(shape, A, Bs)
    _broadcast!(f, C, keeps, Idefaults, A, Bs, Val{N}, iter)
    return C
end

@inline function broadcast_t(
        f, T, shape, iter, A::Any, B::GPUArray, rest::Vararg{Any, N}
    ) where N
    C = similar(B, T, shape)
    Bs = (B, rest...)
    keeps, Idefaults = map_newindexer(shape, A, Bs)
    _broadcast!(f, C, keeps, Idefaults, A, Bs, Val{N}, iter)
    return C
end
function broadcast_t(f::Any, ::Type{Any}, ::Any, ::Any, A::GPUArrays.GPUArray, args::Vararg{Any,N}) where N
    error("Return type couldn't be inferred for broadcast. Func: $f, $(typeof(A)), $args")
end

deref(x) = x
deref(x::RefValue) = (x[],) # RefValue doesn't work with CUDAnative so we use Tuple

function _broadcast!(
        func, out::GPUArray,
        keeps::K, Idefaults::ID,
        A::AT, Bs::BT, ::Type{Val{N}}, unused2 # we don't need those arguments
    ) where {N, K, ID, AT, BT}

    shape = UInt32.(size(out))
    args = (A, Bs...)
    descriptor_tuple = ntuple(length(args)) do i
        BInfo(args[i], keeps[i], Idefaults[i])
    end
    gpu_call(broadcast_kernel!, out, (func, out, shape, descriptor_tuple, (A, deref.(Bs)...)))
    out
end

function Base.foreach(func, over::GPUArray, Bs...)
    shape = UInt32.(size(over))
    keeps, Idefaults = map_newindexer(shape, over, Bs)
    args = (over, Bs...)
    descriptor_tuple = ntuple(length(args)) do i
        BInfo(args[i], keeps[i], Idefaults[i])
    end
    gpu_call(foreach_kernel, over, (func, shape, descriptor_tuple, (over, deref.(Bs)...)))
    return
end

arg_length(x::Tuple) = (UInt32(length(x)),)
arg_length(x::AbstractArray) = UInt32.(size(x))
arg_length(x) = () # Scalar

struct BInfo{Typ, N}
    size::NTuple{N, UInt32}
    keep::NTuple{N, UInt32}
    idefault::NTuple{N, UInt32}
end

function BInfo(val, keep, idefault)
    N = length(keep)
    typ = Broadcast.containertype(val)
    BInfo{typ, N}(arg_length(val), UInt32.(keep), UInt32.(idefault))
end

@propagate_inbounds @inline function _broadcast_getindex(
        ::BInfo{Array}, A, I
    )
    A[I]
end
@propagate_inbounds @inline function _broadcast_getindex(
        ::BInfo{Tuple}, A, I
    )
    A[I]
end
@propagate_inbounds @inline function _broadcast_getindex(
        ::BInfo{Array}, A::Ref, I
    )
    A[]
end

@inline _broadcast_getindex(any, A, I) = A

@inline function broadcast_kernel!(state, func, out, shape, descriptor, args)
    ilin = @linearidx(out, state)
    @inbounds out[ilin] = apply_broadcast(ilin, func, shape, descriptor, args)
    return
end
function foreach_kernel(state, func, shape, descriptor, args)
    ilin = @linearidx(args[1], state)
    apply_broadcast(ilin, func, shape, descriptor, args)
    return
end

function mapidx_kernel(state, f, A, args)
    ilin = @linearidx(A, state)
    f(ilin, A, args...)
    return
end

for N = 0:15
    nargs = N + 1
    inner_expr = []
    valargs = []
    for i = 1:N
        val_i = Symbol("val_", i); I_i = Symbol("I_", i);
        desi = Symbol("deref_", i)
        inner = quote
            # destructure the keeps and As tuples
            $desi = descriptor[$i]
            # reverse-broadcast the indices
            $I_i = newindex(
                I, ilin,
                $desi.keep,
                $desi.idefault,
                $desi.size
            )
            # extract array values
            @inbounds $val_i = _broadcast_getindex($desi, args[$i], $I_i)
        end
        push!(inner_expr, inner)
        push!(valargs, val_i)
    end
    @eval begin
        @inline function apply_broadcast(ilin, func, shape, descriptor, args::NTuple{$N, Any})
            # this will hopefully get dead code removed,
            # if only arrays with linear index are involved, because I should be unused in that case
            I = gpu_ind2sub(shape, ilin)
            $(inner_expr...)
            # call the function and store the result
            func($(valargs...))
        end
    end
end

function mapidx(f, A::GPUArray, args::NTuple{N, Any}) where N
    gpu_call(mapidx_kernel, A, (f, A, args))
end

# don't do anything for empty tuples
@pure newindex(I, ilin, keep::Tuple{}, Idefault::Tuple{}, size::Tuple{}) = UInt32(1)

# optimize for 1D arrays
@pure function newindex(I::NTuple{1}, ilin, keep::NTuple{1}, Idefault, size)
    if Bool(keep[1])
        return ilin
    else
        return Idefault[1]
    end
end

# differently shaped arrays
@generated function newindex(I, ilin::T, keep::NTuple{N}, Idefault, size) where {N, T}
    exprs = Expr(:tuple)
    for i = 1:N
        push!(exprs.args, :(T(Bool(keep[$i]) ? T(I[$i]) : T(Idefault[$i]))))
    end
    :(Base.@_inline_meta; gpu_sub2ind(size, $exprs))
end


text:
module ScientificTypes

# Dependencies
using ScientificTypesBase
using Tables
using CategoricalArrays
using ColorTypes
using PersistenceDiagramsBase
using PrettyTables
using Dates

import StatisticalTraits: info

# re-exports from ScientificTypes
export Scientific, Found, Unknown, Known, Finite, Infinite,
    OrderedFactor, Multiclass, Count, Continuous, Textual,
    Binary, ColorImage, GrayImage, Image, Table,
    ScientificTimeType, ScientificDate, ScientificDateTime,
    ScientificTime
export scitype, scitype_union, elscitype, nonmissing, trait

# re-export from StatisticalTraits
export info

# exports
export coerce, coerce!, autotype, schema

# -------------------------------------------------------------
# Abbreviations

const ST   = ScientificTypesBase
const Arr  = AbstractArray
const CArr = CategoricalArray
const Cat  = CategoricalValue

# Indicate the convention, see init.jl where it is set.
struct DefaultConvention <: Convention end

include("init.jl")

# -------------------------------------------------------------
# Includes

include("coerce.jl")
include("schema.jl")
include("autotype.jl")

include("convention/scitype.jl")
include("convention/coerce.jl")

end # module


text:
using QuantumESPRESSOBase
using Test

@testset "QuantumESPRESSOBase.jl" begin
    include("Inputs.jl")
end


text:
# This file was generated, do not modify it. # hide
using MLJ
using PrettyPrinting
MLJ.color_off() # hide
X, y = @load_iris
DecisionTreeClassifier = @load DecisionTreeClassifier pkg=DecisionTree

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials