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

"""
    search(string::AbstractString, chars::Chars, [start::Integer])

Search for the first occurrence of the given characters within the given string. The second
argument may be a single character, a vector or a set of characters, a string, or a regular
expression (though regular expressions are only allowed on contiguous strings, such as ASCII
or UTF-8 strings). The third argument optionally specifies a starting index. The return
value is a range of indexes where the matching sequence is found, such that `s[search(s,x)] == x`:

`search(string, "substring")` = `start:end` such that `string[start:end] == "substring"`, or
`0:-1` if unmatched.

`search(string, 'c')` = `index` such that `string[index] == 'c'`, or `0` if unmatched.

```jldoctest
julia> search("Hello to the world", "z")
0:-1

julia> search("JuliaLang","Julia")
1:5
```
"""
function search(s::AbstractString, c::Chars, i::Integer)
    if isempty(c)
        return 1 <= i <= nextind(s,endof(s)) ? i :
               throw(BoundsError(s, i))
    end
    if i < 1 || i > nextind(s,endof(s))
        throw(BoundsError(s, i))
    end
    while !done(s,i)
        d, j = next(s,i)
        if d in c
            return i
        end
        i = j
    end
    return 0
end
search(s::AbstractString, c::Chars) = search(s,c,start(s))

in(c::Char, s::AbstractString) = (search(s,c)!=0)

function _searchindex(s, t, i)
    if isempty(t)
        return 1 <= i <= nextind(s,endof(s)) ? i :
               throw(BoundsError(s, i))
    end
    t1, j2 = next(t,start(t))
    while true
        i = search(s,t1,i)
        if i == 0 return 0 end
        c, ii = next(s,i)
        j = j2; k = ii
        matched = true
        while !done(t,j)
            if done(s,k)
                matched = false
                break
            end
            c, k = next(s,k)
            d, j = next(t,j)
            if c != d
                matched = false
                break
            end
        end
        if matched
            return i
        end
        i = ii
    end
end

function _search_bloom_mask(c)
    UInt64(1) << (c & 63)
end

function _searchindex(s::Array, t::Array, i)
    n = length(t)
    m = length(s)

    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return search(s, t[1], i)
    end

    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = t[end]
    for j in 1:n
        bloom_mask |= _search_bloom_mask(t[j])
        if t[j] == tlast && j < n
            skip = n - j - 1
        end
    end

    i -= 1
    while i <= w
        if s[i+n] == tlast
            # check candidate
            j = 0
            while j < n - 1
                if s[i+j+1] != t[j+1]
                    break
                end
                j += 1
            end

            # match found
            if j == n - 1
                return i+1
            end

            # no match, try to rule out the next character
            if i < w && bloom_mask & _search_bloom_mask(s[i+n+1]) == 0
                i += n
            else
                i += skip
            end
        elseif i < w
            if bloom_mask & _search_bloom_mask(s[i+n+1]) == 0
                i += n
            end
        end
        i += 1
    end

    0
end

searchindex(s::ByteArray, t::ByteArray, i) = _searchindex(s,t,i)

"""
    searchindex(s::AbstractString, substring, [start::Integer])

Similar to [`search`](:func:`search`), but return only the start index at which
the substring is found, or `0` if it is not.
"""
searchindex(s::AbstractString, t::AbstractString, i::Integer) = _searchindex(s,t,i)
searchindex(s::AbstractString, t::AbstractString) = searchindex(s,t,start(s))
searchindex(s::AbstractString, c::Char, i::Integer) = _searchindex(s,c,i)
searchindex(s::AbstractString, c::Char) = searchindex(s,c,start(s))

function searchindex(s::String, t::String, i::Integer=1)
    # Check for fast case of a single byte
    # (for multi-byte UTF-8 sequences, use searchindex on byte arrays instead)
    if endof(t) == 1
        search(s, t[1], i)
    else
        searchindex(s.data, t.data, i)
    end
end

function search(s::ByteArray, t::ByteArray, i)
    idx = searchindex(s,t,i)
    if isempty(t)
        idx:idx-1
    else
        idx:(idx > 0 ? idx + endof(t) - 1 : -1)
    end
end

function search(s::AbstractString, t::AbstractString, i::Integer=start(s))
    idx = searchindex(s,t,i)
    if isempty(t)
        idx:idx-1
    else
        idx:(idx > 0 ? idx + endof(t) - 1 : -1)
    end
end

function rsearch(s::AbstractString, c::Chars)
    j = search(RevString(s), c)
    j == 0 && return 0
    endof(s)-j+1
end

"""
    rsearch(s::AbstractString, chars::Chars, [start::Integer])

Similar to [`search`](:func:`search`), but returning the last occurrence of the given characters within the
given string, searching in reverse from `start`.

```jldoctest
julia> rsearch("aaabbb","b")
6:6
```
"""
function rsearch(s::AbstractString, c::Chars, i::Integer)
    e = endof(s)
    j = search(RevString(s), c, e-i+1)
    j == 0 && return 0
    e-j+1
end

function _rsearchindex(s, t, i)
    if isempty(t)
        return 1 <= i <= nextind(s,endof(s)) ? i :
               throw(BoundsError(s, i))
    end
    t = RevString(t)
    rs = RevString(s)
    l = endof(s)
    t1, j2 = next(t,start(t))
    while true
        i = rsearch(s,t1,i)
        if i == 0 return 0 end
        c, ii = next(rs,l-i+1)
        j = j2; k = ii
        matched = true
        while !done(t,j)
            if done(rs,k)
                matched = false
                break
            end
            c, k = next(rs,k)
            d, j = next(t,j)
            if c != d
                matched = false
                break
            end
        end
        if matched
            return nextind(s,l-k+1)
        end
        i = l-ii+1
    end
end

function _rsearchindex(s::Array, t::Array, k)
    n = length(t)
    m = length(s)

    if n == 0
        return 0 <= k <= m ? max(k, 1) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return rsearch(s, t[1], k)
    end

    w = m - n
    if w < 0 || k <= 0
        return 0
    end

    bloom_mask = UInt64(0)
    skip = n - 1
    tfirst = t[1]
    for j in n:-1:1
        bloom_mask |= _search_bloom_mask(t[j])
        if t[j] == tfirst && j > 1
            skip = j - 2
        end
    end

    i = min(k - n + 1, w + 1)
    while i > 0
        if s[i] == tfirst
            # check candidate
            j = 1
            while j < n
                if s[i+j] != t[j+1]
                    break
                end
                j += 1
            end

            # match found
            if j == n
                return i
            end

            # no match, try to rule out the next character
            if i > 1 && bloom_mask & _search_bloom_mask(s[i-1]) == 0
                i -= n
            else
                i -= skip
            end
        elseif i > 1
            if bloom_mask & _search_bloom_mask(s[i-1]) == 0
                i -= n
            end
        end
        i -= 1
    end

    0
end

rsearchindex(s::ByteArray,t::ByteArray,i) = _rsearchindex(s,t,i)

"""
    rsearchindex(s::AbstractString, substring, [start::Integer])

Similar to [`rsearch`](:func:`rsearch`), but return only the start index at which the substring is found, or `0` if it is not.
"""
rsearchindex(s::AbstractString, t::AbstractString, i::Integer) = _rsearchindex(s,t,i)
rsearchindex(s::AbstractString, t::AbstractString) = (isempty(s) && isempty(t)) ? 1 : rsearchindex(s,t,endof(s))

function rsearchindex(s::String, t::String)
    # Check for fast case of a single byte
    # (for multi-byte UTF-8 sequences, use rsearchindex instead)
    if endof(t) == 1
        rsearch(s, t[1])
    else
        _rsearchindex(s.data, t.data, length(s.data))
    end
end

function rsearchindex(s::String, t::String, i::Integer)
    # Check for fast case of a single byte
    # (for multi-byte UTF-8 sequences, use rsearchindex instead)
    if endof(t) == 1
        rsearch(s, t[1], i)
    elseif endof(t) != 0
        _rsearchindex(s.data, t.data, nextind(s, i)-1)
    elseif i > sizeof(s)
        return 0
    elseif i == 0
        return 1
    else
        return i
    end
end

function rsearch(s::ByteArray, t::ByteArray, i::Integer)
    idx = rsearchindex(s,t,i)
    if isempty(t)
        idx:idx-1
    else
        idx:(idx > 0 ? idx + endof(t) - 1 : -1)
    end
end

function rsearch(s::AbstractString, t::AbstractString, i::Integer=endof(s))
    idx = rsearchindex(s,t,i)
    if isempty(t)
        idx:idx-1
    else
        idx:(idx > 0 ? idx + endof(t) - 1 : -1)
    end
end

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

Determine whether the second argument is a substring of the first.

```jldoctest
julia> contains("JuliaLang is pretty cool!", "Julia")
true
```
"""
contains(haystack::AbstractString, needle::AbstractString) = searchindex(haystack,needle)!=0

in(::AbstractString, ::AbstractString) = error("use contains(x,y) for string containment")

# ByteArray optimizations

# find the index of the first occurrence of a value in a byte array

function search(a::ByteArray, b::Union{Int8,UInt8}, i::Integer)
    if i < 1
        throw(BoundsError(a, i))
    end
    n = length(a)
    if i > n
        return i == n+1 ? 0 : throw(BoundsError(a, i))
    end
    p = pointer(a)
    q = ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p+i-1, b, n-i+1)
    q == C_NULL ? 0 : Int(q-p+1)
end
function search(a::ByteArray, b::Char, i::Integer)
    if isascii(b)
        search(a,UInt8(b),i)
    else
        search(a,string(b).data,i).start
    end
end
search(a::ByteArray, b::Union{Int8,UInt8,Char}) = search(a,b,1)

function rsearch(a::ByteArray, b::Union{Int8,UInt8}, i::Integer)
    if i < 1
        return i == 0 ? 0 : throw(BoundsError(a, i))
    end
    n = length(a)
    if i > n
        return i == n+1 ? 0 : throw(BoundsError(a, i))
    end
    p = pointer(a)
    q = ccall(:memrchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p, b, i)
    q == C_NULL ? 0 : Int(q-p+1)
end
function rsearch(a::ByteArray, b::Char, i::Integer)
    if isascii(b)
        rsearch(a,UInt8(b),i)
    else
        rsearch(a,string(b).data,i).start
    end
end
rsearch(a::ByteArray, b::Union{Int8,UInt8,Char}) = rsearch(a,b,length(a))

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
export safe_add_physical_group

"""
    safe_add_physical_group(name::String, entities::Vector{NTuple{2,Int32}})

A combined `gmsh.model.add_physical_group` and `gmsh.model.set_physical_name` that
accounts for physical groups that already exist with name `name`.

Returns the tag of the physical group.
"""
function safe_add_physical_group(name::String, entities::Vector{NTuple{2,I}}) where {I<:Integer}
    dims = unique!(getindex.(entities, 1))
    for dim in dims
        current_groups = gmsh.model.get_physical_groups(dim)
        current_names = [gmsh.model.get_physical_name(dim, grp[2]) for grp in current_groups]
        if name ∈ current_names
            id = findfirst(x->x == name, current_names)
            tag = current_groups[id][2]
            tags = gmsh.model.get_entities_for_physical_group(dim, tag)
            append!(entities, [(dim, t) for t in tags])
            unique!(entities)
            gmsh.model.remove_physical_groups([(dim, tag)])
            p = gmsh.model.add_physical_group(dim,
                                             [ent[2] for ent in filter(e->e[1] == dim, entities)],
                                             tag
                                             )
        else
            p = gmsh.model.add_physical_group(dim,
                                             [ent[2] for ent in filter(e->e[1] == dim, entities)]
                                             )
        end
        gmsh.model.set_physical_name(dim, p, name)
    end
    return nothing
end

ext: jl

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: budgets
using AWS.Compat
using AWS.UUIDs

"""
    create_budget(account_id, budget)
    create_budget(account_id, budget, params::Dict{String,<:Any})

Creates a budget and, if included, notifications and subscribers.   Only one of BudgetLimit
or PlannedBudgetLimits can be present in the syntax at one time. Use the syntax that
matches your case. The Request Syntax section shows the BudgetLimit syntax. For
PlannedBudgetLimits, see the Examples section.

# Arguments
- `account_id`: The accountId that is associated with the budget.
- `budget`: The budget object that you want to create.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"NotificationsWithSubscribers"`: A notification that you want to associate with a
  budget. A budget can have up to five notifications, and each notification can have one SNS
  subscriber and up to 10 email subscribers. If you include notifications and subscribers in
  your CreateBudget call, AWS creates the notifications and subscribers for you.
"""
create_budget(AccountId, Budget; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateBudget", Dict{String, Any}("AccountId"=>AccountId, "Budget"=>Budget); aws_config=aws_config)
create_budget(AccountId, Budget, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateBudget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "Budget"=>Budget), params)); aws_config=aws_config)

"""
    create_budget_action(account_id, action_threshold, action_type, approval_model, budget_name, definition, execution_role_arn, notification_type, subscribers)
    create_budget_action(account_id, action_threshold, action_type, approval_model, budget_name, definition, execution_role_arn, notification_type, subscribers, params::Dict{String,<:Any})

 Creates a budget action.

# Arguments
- `account_id`:
- `action_threshold`:
- `action_type`:  The type of action. This defines the type of tasks that can be carried
  out by this action. This field also determines the format for definition.
- `approval_model`:  This specifies if the action needs manual or automatic approval.
- `budget_name`:
- `definition`:
- `execution_role_arn`:  The role passed for action execution and reversion. Roles and
  actions must be in the same account.
- `notification_type`:
- `subscribers`:

"""
create_budget_action(AccountId, ActionThreshold, ActionType, ApprovalModel, BudgetName, Definition, ExecutionRoleArn, NotificationType, Subscribers; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateBudgetAction", Dict{String, Any}("AccountId"=>AccountId, "ActionThreshold"=>ActionThreshold, "ActionType"=>ActionType, "ApprovalModel"=>ApprovalModel, "BudgetName"=>BudgetName, "Definition"=>Definition, "ExecutionRoleArn"=>ExecutionRoleArn, "NotificationType"=>NotificationType, "Subscribers"=>Subscribers); aws_config=aws_config)
create_budget_action(AccountId, ActionThreshold, ActionType, ApprovalModel, BudgetName, Definition, ExecutionRoleArn, NotificationType, Subscribers, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateBudgetAction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "ActionThreshold"=>ActionThreshold, "ActionType"=>ActionType, "ApprovalModel"=>ApprovalModel, "BudgetName"=>BudgetName, "Definition"=>Definition, "ExecutionRoleArn"=>ExecutionRoleArn, "NotificationType"=>NotificationType, "Subscribers"=>Subscribers), params)); aws_config=aws_config)

"""
    create_notification(account_id, budget_name, notification, subscribers)
    create_notification(account_id, budget_name, notification, subscribers, params::Dict{String,<:Any})

Creates a notification. You must create the budget before you create the associated
notification.

# Arguments
- `account_id`: The accountId that is associated with the budget that you want to create a
  notification for.
- `budget_name`: The name of the budget that you want AWS to notify you about. Budget names
  must be unique within an account.
- `notification`: The notification that you want to create.
- `subscribers`: A list of subscribers that you want to associate with the notification.
  Each notification can have one SNS subscriber and up to 10 email subscribers.

"""
create_notification(AccountId, BudgetName, Notification, Subscribers; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateNotification", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification, "Subscribers"=>Subscribers); aws_config=aws_config)
create_notification(AccountId, BudgetName, Notification, Subscribers, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateNotification", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification, "Subscribers"=>Subscribers), params)); aws_config=aws_config)

"""
    create_subscriber(account_id, budget_name, notification, subscriber)
    create_subscriber(account_id, budget_name, notification, subscriber, params::Dict{String,<:Any})

Creates a subscriber. You must create the associated budget and notification before you
create the subscriber.

# Arguments
- `account_id`: The accountId that is associated with the budget that you want to create a
  subscriber for.
- `budget_name`: The name of the budget that you want to subscribe to. Budget names must be
  unique within an account.
- `notification`: The notification that you want to create a subscriber for.
- `subscriber`: The subscriber that you want to associate with a budget notification.

"""
create_subscriber(AccountId, BudgetName, Notification, Subscriber; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateSubscriber", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification, "Subscriber"=>Subscriber); aws_config=aws_config)
create_subscriber(AccountId, BudgetName, Notification, Subscriber, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("CreateSubscriber", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification, "Subscriber"=>Subscriber), params)); aws_config=aws_config)

"""
    delete_budget(account_id, budget_name)
    delete_budget(account_id, budget_name, params::Dict{String,<:Any})

Deletes a budget. You can delete your budget at any time.  Deleting a budget also deletes
the notifications and subscribers that are associated with that budget.

# Arguments
- `account_id`: The accountId that is associated with the budget that you want to delete.
- `budget_name`: The name of the budget that you want to delete.

"""
delete_budget(AccountId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteBudget", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName); aws_config=aws_config)
delete_budget(AccountId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteBudget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    delete_budget_action(account_id, action_id, budget_name)
    delete_budget_action(account_id, action_id, budget_name, params::Dict{String,<:Any})

 Deletes a budget action.

# Arguments
- `account_id`:
- `action_id`:  A system-generated universally unique identifier (UUID) for the action.
- `budget_name`:

"""
delete_budget_action(AccountId, ActionId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteBudgetAction", Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName); aws_config=aws_config)
delete_budget_action(AccountId, ActionId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteBudgetAction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    delete_notification(account_id, budget_name, notification)
    delete_notification(account_id, budget_name, notification, params::Dict{String,<:Any})

Deletes a notification.  Deleting a notification also deletes the subscribers that are
associated with the notification.

# Arguments
- `account_id`: The accountId that is associated with the budget whose notification you
  want to delete.
- `budget_name`: The name of the budget whose notification you want to delete.
- `notification`: The notification that you want to delete.

"""
delete_notification(AccountId, BudgetName, Notification; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteNotification", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification); aws_config=aws_config)
delete_notification(AccountId, BudgetName, Notification, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteNotification", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification), params)); aws_config=aws_config)

"""
    delete_subscriber(account_id, budget_name, notification, subscriber)
    delete_subscriber(account_id, budget_name, notification, subscriber, params::Dict{String,<:Any})

Deletes a subscriber.  Deleting the last subscriber to a notification also deletes the
notification.

# Arguments
- `account_id`: The accountId that is associated with the budget whose subscriber you want
  to delete.
- `budget_name`: The name of the budget whose subscriber you want to delete.
- `notification`: The notification whose subscriber you want to delete.
- `subscriber`: The subscriber that you want to delete.

"""
delete_subscriber(AccountId, BudgetName, Notification, Subscriber; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteSubscriber", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification, "Subscriber"=>Subscriber); aws_config=aws_config)
delete_subscriber(AccountId, BudgetName, Notification, Subscriber, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DeleteSubscriber", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification, "Subscriber"=>Subscriber), params)); aws_config=aws_config)

"""
    describe_budget(account_id, budget_name)
    describe_budget(account_id, budget_name, params::Dict{String,<:Any})

Describes a budget.  The Request Syntax section shows the BudgetLimit syntax. For
PlannedBudgetLimits, see the Examples section.

# Arguments
- `account_id`: The accountId that is associated with the budget that you want a
  description of.
- `budget_name`: The name of the budget that you want a description of.

"""
describe_budget(AccountId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudget", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName); aws_config=aws_config)
describe_budget(AccountId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    describe_budget_action(account_id, action_id, budget_name)
    describe_budget_action(account_id, action_id, budget_name, params::Dict{String,<:Any})

 Describes a budget action detail.

# Arguments
- `account_id`:
- `action_id`:  A system-generated universally unique identifier (UUID) for the action.
- `budget_name`:

"""
describe_budget_action(AccountId, ActionId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetAction", Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName); aws_config=aws_config)
describe_budget_action(AccountId, ActionId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetAction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    describe_budget_action_histories(account_id, action_id, budget_name)
    describe_budget_action_histories(account_id, action_id, budget_name, params::Dict{String,<:Any})

 Describes a budget action history detail.

# Arguments
- `account_id`:
- `action_id`:  A system-generated universally unique identifier (UUID) for the action.
- `budget_name`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
- `"TimePeriod"`:
"""
describe_budget_action_histories(AccountId, ActionId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetActionHistories", Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName); aws_config=aws_config)
describe_budget_action_histories(AccountId, ActionId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetActionHistories", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    describe_budget_actions_for_account(account_id)
    describe_budget_actions_for_account(account_id, params::Dict{String,<:Any})

 Describes all of the budget actions for an account.

# Arguments
- `account_id`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
"""
describe_budget_actions_for_account(AccountId; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetActionsForAccount", Dict{String, Any}("AccountId"=>AccountId); aws_config=aws_config)
describe_budget_actions_for_account(AccountId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetActionsForAccount", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId), params)); aws_config=aws_config)

"""
    describe_budget_actions_for_budget(account_id, budget_name)
    describe_budget_actions_for_budget(account_id, budget_name, params::Dict{String,<:Any})

 Describes all of the budget actions for a budget.

# Arguments
- `account_id`:
- `budget_name`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
"""
describe_budget_actions_for_budget(AccountId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetActionsForBudget", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName); aws_config=aws_config)
describe_budget_actions_for_budget(AccountId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetActionsForBudget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    describe_budget_performance_history(account_id, budget_name)
    describe_budget_performance_history(account_id, budget_name, params::Dict{String,<:Any})

Describes the history for DAILY, MONTHLY, and QUARTERLY budgets. Budget history isn't
available for ANNUAL budgets.

# Arguments
- `account_id`:
- `budget_name`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
- `"TimePeriod"`: Retrieves how often the budget went into an ALARM state for the specified
  time period.
"""
describe_budget_performance_history(AccountId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetPerformanceHistory", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName); aws_config=aws_config)
describe_budget_performance_history(AccountId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgetPerformanceHistory", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    describe_budgets(account_id)
    describe_budgets(account_id, params::Dict{String,<:Any})

Lists the budgets that are associated with an account.  The Request Syntax section shows
the BudgetLimit syntax. For PlannedBudgetLimits, see the Examples section.

# Arguments
- `account_id`: The accountId that is associated with the budgets that you want
  descriptions of.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: An optional integer that represents how many entries a paginated response
  contains. The maximum is 100.
- `"NextToken"`: The pagination token that you include in your request to indicate the next
  set of results that you want to retrieve.
"""
describe_budgets(AccountId; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgets", Dict{String, Any}("AccountId"=>AccountId); aws_config=aws_config)
describe_budgets(AccountId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeBudgets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId), params)); aws_config=aws_config)

"""
    describe_notifications_for_budget(account_id, budget_name)
    describe_notifications_for_budget(account_id, budget_name, params::Dict{String,<:Any})

Lists the notifications that are associated with a budget.

# Arguments
- `account_id`: The accountId that is associated with the budget whose notifications you
  want descriptions of.
- `budget_name`: The name of the budget whose notifications you want descriptions of.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: An optional integer that represents how many entries a paginated response
  contains. The maximum is 100.
- `"NextToken"`: The pagination token that you include in your request to indicate the next
  set of results that you want to retrieve.
"""
describe_notifications_for_budget(AccountId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeNotificationsForBudget", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName); aws_config=aws_config)
describe_notifications_for_budget(AccountId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeNotificationsForBudget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    describe_subscribers_for_notification(account_id, budget_name, notification)
    describe_subscribers_for_notification(account_id, budget_name, notification, params::Dict{String,<:Any})

Lists the subscribers that are associated with a notification.

# Arguments
- `account_id`: The accountId that is associated with the budget whose subscribers you want
  descriptions of.
- `budget_name`: The name of the budget whose subscribers you want descriptions of.
- `notification`: The notification whose subscribers you want to list.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: An optional integer that represents how many entries a paginated response
  contains. The maximum is 100.
- `"NextToken"`: The pagination token that you include in your request to indicate the next
  set of results that you want to retrieve.
"""
describe_subscribers_for_notification(AccountId, BudgetName, Notification; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeSubscribersForNotification", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification); aws_config=aws_config)
describe_subscribers_for_notification(AccountId, BudgetName, Notification, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("DescribeSubscribersForNotification", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "Notification"=>Notification), params)); aws_config=aws_config)

"""
    execute_budget_action(account_id, action_id, budget_name, execution_type)
    execute_budget_action(account_id, action_id, budget_name, execution_type, params::Dict{String,<:Any})

 Executes a budget action.

# Arguments
- `account_id`:
- `action_id`:  A system-generated universally unique identifier (UUID) for the action.
- `budget_name`:
- `execution_type`:  The type of execution.

"""
execute_budget_action(AccountId, ActionId, BudgetName, ExecutionType; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("ExecuteBudgetAction", Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName, "ExecutionType"=>ExecutionType); aws_config=aws_config)
execute_budget_action(AccountId, ActionId, BudgetName, ExecutionType, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("ExecuteBudgetAction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName, "ExecutionType"=>ExecutionType), params)); aws_config=aws_config)

"""
    update_budget(account_id, new_budget)
    update_budget(account_id, new_budget, params::Dict{String,<:Any})

Updates a budget. You can change every part of a budget except for the budgetName and the
calculatedSpend. When you modify a budget, the calculatedSpend drops to zero until AWS has
new usage data to use for forecasting.  Only one of BudgetLimit or PlannedBudgetLimits can
be present in the syntax at one time. Use the syntax that matches your case. The Request
Syntax section shows the BudgetLimit syntax. For PlannedBudgetLimits, see the Examples
section.

# Arguments
- `account_id`: The accountId that is associated with the budget that you want to update.
- `new_budget`: The budget that you want to update your budget to.

"""
update_budget(AccountId, NewBudget; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateBudget", Dict{String, Any}("AccountId"=>AccountId, "NewBudget"=>NewBudget); aws_config=aws_config)
update_budget(AccountId, NewBudget, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateBudget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "NewBudget"=>NewBudget), params)); aws_config=aws_config)

"""
    update_budget_action(account_id, action_id, budget_name)
    update_budget_action(account_id, action_id, budget_name, params::Dict{String,<:Any})

 Updates a budget action.

# Arguments
- `account_id`:
- `action_id`:  A system-generated universally unique identifier (UUID) for the action.
- `budget_name`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ActionThreshold"`:
- `"ApprovalModel"`:  This specifies if the action needs manual or automatic approval.
- `"Definition"`:
- `"ExecutionRoleArn"`:  The role passed for action execution and reversion. Roles and
  actions must be in the same account.
- `"NotificationType"`:
- `"Subscribers"`:
"""
update_budget_action(AccountId, ActionId, BudgetName; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateBudgetAction", Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName); aws_config=aws_config)
update_budget_action(AccountId, ActionId, BudgetName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateBudgetAction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "ActionId"=>ActionId, "BudgetName"=>BudgetName), params)); aws_config=aws_config)

"""
    update_notification(account_id, budget_name, new_notification, old_notification)
    update_notification(account_id, budget_name, new_notification, old_notification, params::Dict{String,<:Any})

Updates a notification.

# Arguments
- `account_id`: The accountId that is associated with the budget whose notification you
  want to update.
- `budget_name`: The name of the budget whose notification you want to update.
- `new_notification`: The updated notification to be associated with a budget.
- `old_notification`: The previous notification that is associated with a budget.

"""
update_notification(AccountId, BudgetName, NewNotification, OldNotification; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateNotification", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "NewNotification"=>NewNotification, "OldNotification"=>OldNotification); aws_config=aws_config)
update_notification(AccountId, BudgetName, NewNotification, OldNotification, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateNotification", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "NewNotification"=>NewNotification, "OldNotification"=>OldNotification), params)); aws_config=aws_config)

"""
    update_subscriber(account_id, budget_name, new_subscriber, notification, old_subscriber)
    update_subscriber(account_id, budget_name, new_subscriber, notification, old_subscriber, params::Dict{String,<:Any})

Updates a subscriber.

# Arguments
- `account_id`: The accountId that is associated with the budget whose subscriber you want
  to update.
- `budget_name`: The name of the budget whose subscriber you want to update.
- `new_subscriber`: The updated subscriber that is associated with a budget notification.
- `notification`: The notification whose subscriber you want to update.
- `old_subscriber`: The previous subscriber that is associated with a budget notification.

"""
update_subscriber(AccountId, BudgetName, NewSubscriber, Notification, OldSubscriber; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateSubscriber", Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "NewSubscriber"=>NewSubscriber, "Notification"=>Notification, "OldSubscriber"=>OldSubscriber); aws_config=aws_config)
update_subscriber(AccountId, BudgetName, NewSubscriber, Notification, OldSubscriber, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = budgets("UpdateSubscriber", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "BudgetName"=>BudgetName, "NewSubscriber"=>NewSubscriber, "Notification"=>Notification, "OldSubscriber"=>OldSubscriber), params)); aws_config=aws_config)


text:
# Extrapolate correlation functions on three dimensions

## KLUDGE: We only store correlation functions for 3D image here, but
## still have to be subtype of AbstractArray to be able to tell
## dimensionality and element type of an image.

struct ExtrapolatedData{T, N} <: AbstractArray{T, N}
    periodic   :: Bool
    corrdata   :: Dict{AbstractTracker{T}, Directional.CorrelationData}
    corrlen    :: Int
    directions :: Vector{Symbol}
    shape      :: NTuple{N, Int}
end

# AnnealingAPI.jl Interface
AnnealingAPI.tracked_data(x :: ExtrapolatedData) = x.corrdata |> keys
AnnealingAPI.tracked_length(x :: ExtrapolatedData) = x.corrlen
AnnealingAPI.tracked_directions(x :: ExtrapolatedData) = x.directions

# Array interface
Base.size(x :: ExtrapolatedData) = x.shape
Base.getindex(x :: ExtrapolatedData{T}, idx :: Vararg{Int}) where T = zero(T)

# Other standard functions
Base.show(io :: IO, x :: ExtrapolatedData) = print(io, "$(size(x)) $(typeof(x))")
Base.show(io :: IO, :: MIME"text/plain", x :: ExtrapolatedData) = show(io, x)

function new_shape(shape :: NTuple{N, Int}, dimensions :: Int) where N
    mindim = min(N, dimensions)
    minelt = minimum(shape)
    return (shape[1:mindim]..., Tuple(minelt for n in 1:(dimensions - N))...)
end

function extrapolate_vector(data           :: Vector{Float64},
                            data_grid_step :: Float64,
                            extr_grid_step :: Float64)
    grid = range(0; length = length(data), step = data_grid_step)
    itp  = interpolate((grid,), data, Gridded(Linear()))
    ext  = extrapolate(itp, Line())

    new_grid = range(0; length = length(data), step = extr_grid_step)
    return [ext(x) for x in new_grid]
end

function extrapolate_data(data       :: Directional.CorrelationData,
                          directions :: Vector{Symbol})
    # TODO: use real extrapolation here
    corrlen = length(data)
    unit_length = Directional.unit_length

    local function statistics_gen(dict)
        Iterators.map(directions) do direction
            if direction ∈ data
                direction => dict[direction]
            else
                exdata = (extrapolate_vector(dict[dir], unit_length(dir), unit_length(direction))
                          for dir in Directional.directions(data))
                direction => reduce(.+, exdata)
            end
        end
    end

    total   = statistics_gen(data.total)
    success = statistics_gen(data.success)

    return Directional.CorrelationData(directions, Dict(success), Dict(total))
end

"""
    ExtrapolatedData(tracker, dimensions, directions)

Extrapolate correlation functions to other dimensions and/or
directions. Correlation functions are defined for objects of
`ExtrapolatedData` type as they are defined for
`CorrelationTracker`s.

# Example
```jldoctest mylabel
julia> begin
       a = rand(MersenneTwister(1), 0:1, (20, 20));
       tracker = CorrelationTracker(a; periodic = true);
       Directional.s2(tracker, 0)
       end
┌────────┬────────┐
│      x │      y │
├────────┼────────┤
│ 0.4975 │ 0.4975 │
│  0.225 │ 0.2325 │
│  0.265 │   0.25 │
│   0.24 │   0.25 │
│  0.235 │  0.245 │
│ 0.2525 │ 0.2275 │
│ 0.2525 │ 0.2675 │
│ 0.2425 │ 0.2475 │
│  0.245 │  0.235 │
│ 0.2575 │ 0.2475 │
└────────┴────────┘
```

```jdoctest mylabel
julia> begin
       extra = ExtrapolatedData(tracker, 3, [:x, :y, :z, :xyz])
       Directional.s2(extra, 0)
       end
┌────────┬────────┬─────────┬──────────┐
│      x │      y │       z │      xyz │
├────────┼────────┼─────────┼──────────┤
│ 0.4975 │ 0.4975 │  0.4975 │   0.4975 │
│  0.225 │ 0.2325 │ 0.22875 │ 0.249796 │
│  0.265 │   0.25 │  0.2575 │ 0.242679 │
│   0.24 │   0.25 │   0.245 │ 0.243923 │
│  0.235 │  0.245 │    0.24 │ 0.246077 │
│ 0.2525 │ 0.2275 │    0.24 │ 0.248253 │
│ 0.2525 │ 0.2675 │    0.26 │ 0.269904 │
│ 0.2425 │ 0.2475 │   0.245 │ 0.291554 │
│  0.245 │  0.235 │    0.24 │ 0.313205 │
│ 0.2575 │ 0.2475 │  0.2525 │ 0.334856 │
└────────┴────────┴─────────┴──────────┘
```
"""
function ExtrapolatedData(tracker    :: CorrelationTracker{T, N},
                          dimensions :: Int,
                          directions :: Vector{Symbol}) where {T, N}
    shape = new_shape(size(tracker), dimensions)
    corrdata = CorrdataDict{T}(tracker => extrapolate_data(data, directions)
                               for (tracker, data) in tracker.corrdata)
    return ExtrapolatedData{T, dimensions}(tracker.periodic,
                                           corrdata,
                                           tracker.corrlen,
                                           directions,
                                           shape)
end

# TODO: Move to CorrelationFunctions.jl maybe?
function join_data(data1 :: Directional.CorrelationData,
                   data2 :: Directional.CorrelationData)
    @assert Directional.directions(data1) == Directional.directions(data2)
    directions = Directional.directions(data1)

    success = Dict(dir => data1.success[dir] + data2.success[dir]
                   for dir in directions)
    total   = Dict(dir => data1.total[dir] + data2.total[dir]
                   for dir in directions)
    return Directional.CorrelationData(directions, success, total)
end

"""
    ExtrapolatedData(extrapolated_data1, extrapolated_data2)

Return `ExtrapolatedData` object which contains correlation statistics
from both `extrapolated_data1` and `extrapolated_data2`. The two input
objects must share the same properties such as the shape of underlying
data, tracked functions, directions and so on.
"""
function ExtrapolatedData(data1 :: ExtrapolatedData{T, N},
                          data2 :: ExtrapolatedData{T, N}) where {T, N}
    if data1.periodic             != data2.periodic        ||
        size(data1)               != size(data2)           ||
        tracked_data(data1)       != tracked_data(data2)   ||
        tracked_length(data1)     != tracked_length(data2) ||
        tracked_directions(data1) != tracked_directions(data2)
        error("Two ExtrapolationData objects must share the same properties to be joined")
    end

    joined_data = Dict(tracker => join_data(tracker(data1), tracker(data2))
                       for tracker in tracked_data(data1))

    return ExtrapolatedData(data1.periodic,
                            joined_data,
                            tracked_length(data1),
                            tracked_directions(data1),
                            size(data1))
end

Directional.l2(tracker :: ExtrapolatedData{T}, phase) where T =
    tracker.corrdata[L2Tracker{T}(phase)]

Directional.s2(tracker :: ExtrapolatedData{T}, phase) where T =
    tracker.corrdata[S2Tracker{T}(phase)]

Directional.surfsurf(tracker :: ExtrapolatedData{T}, phase) where T =
    tracker.corrdata[SSTracker{T}(phase)]

Directional.surfvoid(tracker :: ExtrapolatedData{T}, phase) where T =
    tracker.corrdata[SVTracker{T}(phase)]


text:
using Revise
using ADCME
using ADCMEKit
using NNFEM
using PyPlot
using ProgressMeter
using MAT
using ADCMEKit
using Clustering

close("all")
include("load_domain_function.jl")

NT = 100
Δt = 30/NT
domain = load_crack_domain()
# domain = load_crack_domain(option="viscoelasticity")

# visualize_mesh(domain)
H = domain.elements[1].mat[1].H
EBC_func = nothing
FBC_func = nothing
body_func = nothing
globaldata = GlobalData(missing, missing, missing, missing, domain.neqs, EBC_func, FBC_func, body_func)
# globaldata, domain = LinearStaticSolver(globaldata, domain)
Fext = compute_external_force(globaldata, domain)
d = LinearStaticSolver(globaldata, domain, domain.state, H, Fext)
sess = Session(); init(sess)
domain.state = run(sess, d)

# error()
#=
figure()
subplot(211)
visualize_scalar_on_scoped_body(domain.state[1:domain.nnodes], zeros(size(domain.state)...), domain)
subplot(212)
visualize_scalar_on_scoped_body(domain.state[domain.nnodes+1:end], zeros(size(domain.state)...), domain)
=#

strain = getStrain(domain)
stress = getStress(domain)
state = domain.state

# error()

domain = load_crack_domain(option="mixed")

ts = GeneralizedAlphaSolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globaldata, ts)
Fext = compute_external_force(domain, globaldata, ts)

d0 = state
v0 = zeros(2domain.nnodes)
a0  = zeros(2domain.nnodes)
# σ0 = zeros(getNGauss(domain),3)
# ϵ0 = zeros(getNGauss(domain),3)
σ0 = stress
ϵ0 = strain

gnodes = getGaussPoints(domain)

using Random; Random.seed!(233)

kr = kmeans(gnodes', 20)
A = kr.assignments
# for i = 1:20
#   scatter(gnodes[A .== i,1], gnodes[A .== i,2])
# end
# gca().invert_yaxis()

vs = Variable(1.5*ones(20))

ETA = constant(zeros(getNGauss(domain)))
for i = 1:20
  mask = zeros(getNGauss(domain))
  mask[A.==i] .= 1.0
  global ETA += vs[i] * mask
end
ETA *= 1e10

μ = zeros(getNGauss(domain))
λ = zeros(getNGauss(domain))
η = zeros(getNGauss(domain))
k = 0
for i = 1:domain.neles
  e = domain.elements[i]
  for mat in e.mat
    global k += 1
    μ[k] = mat.μ
    λ[k] = mat.λ
    η[k] = mat.η
  end
end

# η = 1e10 + X2 * gnodes[:,1] + X3 * gnodes[:,2]
# η = @. 1e10 *( 2 + abs(gnodes[:,1] - maximum(gnodes[:,1])/2)/maximum(gnodes[:,1]) - (gnodes[:,2])/maximum(gnodes[:,2]))
# η = 1e10 * ( 2 .- (gnodes[:,2])/maximum(gnodes[:,2]))

# η = 1e10 *Variable( 1.5 ) * ones(getNGauss(domain))

# η = 1e10 *( Variable(1.5) + Variable(1.5) * (gnodes[:,2])/maximum(gnodes[:,2]))
η = ETA

d, v, a, σ, ϵ = ViscoelasticitySolver(
  globaldata, domain, d0, v0, a0, σ0, ϵ0, Δt, NT, μ, λ, η, Fext, ubd, abd
)

d_ = matread("data/viscoelasticity_linear2.mat")["d"]

loss = sum((d-d_)^2)
sess = Session(); init(sess)

@info run(sess, loss)
BFGS!(sess, loss, 30)
# lineview(sess, pl, loss, [20.0], [log(1e10)])

# d_ = run(sess, d)
# matwrite("data/viscoelasticity_linear2.mat", Dict("d"=>d_))

figure()
η_ = run(sess, η[1:9:end])
visualize_scalar_on_scoped_body(η_, zeros(domain.nnodes*2), domain)

figure()
η2 = 1e10 *( 2 .- (gnodes[:,2])/maximum(gnodes[:,2]))
η2 = η2[1:9:end]
visualize_scalar_on_scoped_body(η2, zeros(domain.nnodes*2), domain)

#=
d_, σ_ = run(sess, [d, σ])
domain.history["state"] = []
domain.history["stress"] = []
for i = 1:NT+1
  push!(domain.history["state"], d_[i,:])
  push!(domain.history["stress"], σ_[i,:,:])
end

###################### check surface ########################
idx = []
k = 0
pts = getGaussPoints(domain)
xcoord = Float64[]
for e in domain.elements
    global k
    if sum(map(x->x≈0.0, e.coords[:,2]))≠2
        k += length(e.weights)
        continue
    end
    l = argmin(pts[k+1:k+length(e.weights), 2])
    push!(idx, k + l)
    push!(xcoord, pts[k+l, 1])
    k += length(e.weights)
end
ii = sortperm(xcoord)
xcoord = xcoord[ii]
idx = idx[ii]s

σ_surface = σ_[:, idx, :]
s = zeros(size(σ_surface, 1), size(σ_surface, 2))
for i = 1:size(σ_surface, 1)
    for j = 1:size(σ_surface, 2)
        s[i,j] = NNFEM.postprocess_stress(σ_surface[i,j,:], "vonMises")
    end
end

figure(figsize=(8,3))
semilogy(s[1, :], "-", label="x0")
semilogy(s[NT+1, :], "--", label="x")
legend()
title("Stress")
savefig("Stress_test.png")

figure(figsize=(8,3))
semilogy(s[:, 13], "-")
title("Stress")

slip_idx = findall(domain.nodes[:,2] .≈ 0.0)
coords = domain.nodes[slip_idx,1]
ii = sortperm(coords)
coords = coords[ii]
x_ = d_[:, slip_idx[ii]]
y_ = d_[:, slip_idx[ii] .+ domain.nnodes]
# plot(coords, x_)

x0 = d0[slip_idx[ii]]
y0 = d0[slip_idx[ii] .+ domain.nnodes]

figure(figsize=(8,3))
subplot(211)
plot((x_' .- x0)[:,1:10:end])
subplot(212)
plot(x0, label="x")
# plot(x_[1, :] .- x0, "-", label="x0")
# plot(x_[NT+1, :] .- x0, "--", label="x")
legend()

figure(figsize=(8,3))
subplot(211)
S = 1:20:NT
for i = 1:length(S)
  plot(-(y_' .- y0)[:,S[i]],"C$i",  label="$i")
end
legend()

subplot(212)
S = 1:20:NT
for i = 1:length(S)
  plot(-y_'[:,S[i]],"C$i",  label="$i")
end
plot(-y0, label="y")
# plot(y_[1, :] .- y0, "-", label="x0")
# plot(y_[NT+1, :] .- y0, "--", label="x")
# plot(y0, ":", label="static")
legend()
savefig("disp.png")

figure()
plot(y_[:, 10])
# plot(y_[1, :] .- y0, "-", label="x0")
# plot(y_[NT+1, :] .- y0, "--", label="x")
# plot(y0, ":", label="static")
legend()
=#

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
module MLInterpret

using Random, Statistics, Combinatorics, Distributed
using ProgressMeter, PyCall, PyCallUtils, PandasLite
using PyCall: python
using PandasLite: PandasWrapped
import StatsBase: sample, Weights
using Pkg.GitTools: clone

export interpret, fit_surrogate, dai_interpret, sbrl_interpret, dnn_interpret

hasproperty(pyo, s) = !ispynull(pyo) && PyCall.hasproperty(pyo, s)

bash(cmd) = run(`bash -c "$cmd"`)

const pdfcat = joinpath(@__DIR__, "pdfcat")

macro savefig(dst, ex)
    dst, ex = esc(dst), esc(ex)
    quote
        @imports matplotlib.pyplot as plt
        plt.matplotlib.use("agg")
        plt.close(plt.gcf())
        $ex; mkpath(dirname($dst))
        plt.savefig($dst, bbox_inches = "tight")
    end
end

macro savefig(dir, n, ex)
    quote
        dst = joinpath($dir, lpad($n, 2, "0") * ".pdf")
        @savefig dst $ex
    end |> esc
end

macro indir(dir, ex)
    cwd = gensym()
    quote
        mkpath($dir); $cwd = pwd(); cd($dir)
        try $ex finally cd($cwd) end
    end |> esc
end

macro redirect(src, ex)
    src = src == :devnull ? "/dev/null" : src
    quote
        io = open($(esc(src)), "a")
        o, e = stdout, stderr
        redirect_stdout(io)
        redirect_stderr(io)
        try
            $(esc(ex)); sleep(0.01)
        finally
            flush(io); close(io)
            redirect_stdout(o)
            redirect_stderr(e)
        end
    end
end

istree(model) = occursin(r"gbm|booster|tree|forest"i, string(model)) && hasproperty(model, :predict)

function interpret(model, X, y; nskater = 50000, nsurro = 500000, nshap = 50000, nshap2 = 3000)
    X.columns = X.columns.astype("str")
    X = X.fillna(X.mean(axis = 0).fillna(0))
    tree = istree(model) ? model : fit_surrogate(X, y)
    pyo =  hasproperty(model, :predict) ? model : tree
    @indir "mli" begin
        surrogate_interpret(pyo.predict, X, y, nsample = nsurro)
        cols = shap_interpret(tree, X, nsample = nshap)
        pdp_interpret(pyo, X, nsample = nskater, cols = cols)
        shap2_interpret(X, y, nsample = nshap2, cols = cols)
        # skater_interpret(pyo.predict, X, y, nsample = nskater)
    end
end

interpret(X, y; ka...) = interpret(PyNULL(), X, y; ka...)

function sample(x::PandasWrapped, nsample)
    nsample >= length(x) && return x
    x = x.sample(nsample, random_state = 12345)
    x.reset_index(inplace = true, drop = true)
    return x
end

sample(X::PandasWrapped, y, nsample) = sample(X, nsample), sample(y, nsample)

function fit_surrogate(X, y; max_depth = 8)
    @from unidecode imports unidecode
    cols = [filter(!isspace, unidecode(c)) for c in X.columns]
    X = DataFrame(X, columns = cols)
    # TODO: tuning parameters
    @from lightgbm imports LGBMRegressor
    @from sklearn.metrics imports r2_score
    surrogate = LGBMRegressor(
        n_jobs = 20, max_bin = 15, boosting_type = "rf",
        subsample = 0.632, colsample_bytree = 0.632,
        subsample_freq = 1, num_leaves = 2^max_depth - 1
    )
    surrogate.fit(X, y, eval_set = [(X, y)])
    ŷ = surrogate.predict(X)
    println("r2: ", r2_score(y, ŷ))
    return surrogate
end

function sbrl_interpret(X, y; nsample = length(y), sample_method = "sample", cols = X.columns)
    @from skater.core.global_interpretation.interpretable_models.brlc imports BRLC
    @from skater.core.global_interpretation.interpretable_models.bigdatabrlc imports BigDataBRLC
    if sample_method == "sample"
        sbrl_model = BRLC(iterations = 30000, n_chains = 10, max_rule_len = 4, drop_features = true)
        X, y = sample(X[cols], y, nsample)
    elseif sample_method == "surrogate"
        pct = min(nsample / length(y), 1)
        sbrl_model = BigDataBRLC(iterations = 30000, n_chains = 10, drop_features = true, sub_sample_percentage = pct)
    end
    @indir "mli" begin
        sbrl_model.fit(X[cols], Float32.(y .> quantile(y, 0.7)))
        @redirect "sbrl.txt" sbrl_model.print_model()
        # sbrl_model.save_model("sbrl.pkl")
        bash("rm tdata_*")
    end
end

function surrogate_interpret(model, X, y; nsample = length(X))
    X, y = sample(X, y, nsample)
    @from skater.model imports InMemoryModel
    @from skater.core.global_interpretation.tree_surrogate imports TreeSurrogate
    skater_model = InMemoryModel(model, feature_names = X.columns, examples = X.iloc[1:2])
    for n in broadcast(^, 2, [3, 4, 5])
        surrogate_explainer = TreeSurrogate(skater_model, max_leaf_nodes = 2n)
        surrogate_explainer.fit(X, y, use_oracle = false)
        surrogate_explainer.plot_global_decisions(file_name = "surrogate_tree-$n.png")
    end
end

function skater_interpret(model, X, y; nsample = length(X), ntop = 20)
    X, y = sample(X, y, nsample)
    @from skater.core.explanations imports Interpretation
    @from skater.model imports InMemoryModel
    interpreter = Interpretation(X, y)
    skater_model = InMemoryModel(model, feature_names = X.columns, examples = X.iloc[1:2])
    sr = interpreter.feature_importance.feature_importance(skater_model, n_samples = 20000)
    sr.sort_values(ascending = false, inplace = true)
    sr.to_csv("perturb_feaimpt.csv", header = false, encoding = "gbk")
    ppd = interpreter.partial_dependence.plot_partial_dependence
    cols = sr.head(ntop).index.tolist()
    @showprogress for (n, c) in enumerate(cols)
        @savefig "pdp" n ppd([c], skater_model, with_variance = true, grid_resolution = 10, sample = false)
    end
    bash("$python $pdfcat -o pdp.pdf pdp/*.pdf && rm -rf pdp")
end

function pdp_interpret(model, X; nsample = length(X), ntop = 7, cols = X.columns)
    X = sample(X, nsample)
    @from pdpbox.info_plots imports actual_plot, actual_plot_interact
    @from warnings imports filterwarnings
    filterwarnings("ignore"; :module => "pdpbox")
    plot_params = Dict("font_family" => "sans-serif")
    for (n, c) in enumerate(cols)
        @savefig "actual" n actual_plot(model, X, c, c, predict_kwds = Dict(), plot_params = plot_params)
    end
    bash("$python $pdfcat -o actual.pdf actual/*.pdf && rm -rf actual")
    for (n, cc) in enumerate(combinations(cols[1:min(ntop, end)], 2))
        @savefig "actual2" n actual_plot_interact(model, X, cc, cc, predict_kwds = Dict(), plot_params = plot_params)
    end
    bash("$python $pdfcat -o actual2.pdf actual2/*.pdf && rm -rf actual2")
end

function pyrun(str)
    println(str)
    cmd = Cmd(["python", "-c", str])
    faketime = joinpath(DEPOT_PATH[1], "faketime/lib/libfaketimeMT.so.1")
    withenv("LD_PRELOAD" => faketime, "FAKETIME" => "-10year") do
        lines = readlines(cmd)
        println(join(lines, '\n'))
        lines[end]
    end
end

function dai_interpret(X, y; nsample = length(X))
    whl="http://localhost:12345/static/h2oai_client-1.7.0-py3-none-any.whl"
    run(pipeline(`$python -m pip install $whl`, stdout = devnull))
    X, y = sample(X, y, nsample)
    target = Series(y, name = "target")
    pred = Series(y, name = "pred")
    df = pd.concat([X, target, pred], axis = 1)
    df.to_parquet("pred.parquet")
    desc = """
    from h2oai_client import Client
    h2oai = Client(address='http://127.0.0.1:12345', username='username', password='password')
    pred = h2oai.upload_dataset_sync('pred.parquet')
    mli = h2oai.run_interpretation_sync('', pred.key, 'target', prediction_col='pred')
    print(mli.description)
    """ |> pyrun
    rm("pred.parquet")
    touch(desc * ".mli")
end

function shap_interpret(model, X; nsample = length(X), ntop = 20)
    X = sample(X, nsample)
    @from shap imports TreeExplainer, summary_plot, dependence_plot, force_plot, save_html
    explainer = TreeExplainer(model)
    shap_vals = explainer.shap_values(X)
    cols = shap_topfeas(shap_vals, X, ntop)
    @savefig "summary/dot.pdf" summary_plot(shap_vals, X)
    @savefig "summary/bar.pdf" summary_plot(shap_vals, X, plot_type = "bar")
    for (n, c) in enumerate(cols)
        @savefig "depend" n dependence_plot(c, shap_vals, X)
    end
    bash("$python $pdfcat -o shap.pdf summary/*.pdf depend/*.pdf && rm -rf summary depend")
    i = randsubseq(1:length(X), min(1, 5000 / length(X)))
    j = [findfirst(c .== X.columns) for c in cols]
    html = force_plot(explainer.expected_value, shap_vals[i, j], X.iloc[i, j])
    save_html("force_plot.html", html)
    return cols
end

function shap2_interpret(X, y; nsample = length(X), ntop = 7, cols = X.columns)
    model = fit_surrogate(X[cols], y)
    X, y = sample(X[cols], y, nsample)
    @from shap imports TreeExplainer, summary_plot, dependence_plot
    explainer = TreeExplainer(model)
    shap_vals = explainer.shap_interaction_values(X)
    @savefig "summary/dot.pdf" summary_plot(shap_vals, X)
    @savefig "summary/bar.pdf" summary_plot_bar2(shap_vals, X)
    for (n, cc) in enumerate(combinations(cols[1:min(ntop, end)], 2))
        @savefig "depend" n dependence_plot(cc, shap_vals, X)
    end
    bash("$python $pdfcat -o shap2.pdf summary/*.pdf depend/*.pdf && rm -rf summary depend")
end

function shap_topfeas(shap_vals, X, ntop = size(X, 2))
    feaimpt = vec(sum(abs, shap_vals, dims = 1))
    is = sortperm(feaimpt, rev = true)[1:min(ntop, end)]
    collect(X.columns[is])
end

function summary_plot_bar2(shap_vals, X)
    M = dropdims(sum(abs, shap_vals, dims = 1), dims = 1)
    for i in axes(M, 1) M[i, i] = 0 end
    feaimpt = vec(sum(M, dims = 2))
    is = sortperm(feaimpt, rev = true)
    M, cols = M[is, is], X.columns[is]
    @imports matplotlib.pyplot as plt
    plt.imshow(M)
    plt.yticks(axes(M, 1) .- 1, cols)
    plt.xticks(axes(M, 2) .- 1, cols)
    plt.gca().xaxis.tick_top()
end

function dnn_interpret(a...; ka...)
    id, = addprocs(1)
    @eval @everywhere using MLI
    remotecall_fetch(_dnn_interpret, id, a...; ka...)
    rmprocs(id)
end

function _dnn_interpret(h5::String, a...; ka...)
    @from keras.models imports load_model
    model = load_model(h5, compile = false)
    _dnn_interpret(model, a...; ka...)
end

function _dnn_interpret(model, X; alg = "deeplift", width = 20, nsample = 5000, nviz = 100, recep = 0, cols = [])
    X = sample_clip(X, nsample, recep)
    if model.output.shape.ndims == 3
        @from keras.layers imports Lambda
        @from keras.models imports Model
        output = Lambda(py"lambda x: x[:, -1, :]")(model.output)
        model = Model(model.inputs, output)
    end
    Xviz = sample_rank(X, vec(model.predict(X)), nviz)
    @from shap imports DeepExplainer, GradientExplainer, image_plot
    class = alg == "deeplift" ? DeepExplainer : GradientExplainer
    shap_vals, = class(model, X).shap_values(Xviz)
    @imports matplotlib.pyplot as pl
    @from shap.plots imports colors
    @from matplotlib.backends.backend_pdf imports PdfPages
    @imports numpy as np
    max_val = np.nanpercentile(abs.(shap_vals), 99.9)
    xmax = np.nanpercentile(X, 99.9)
    xmin = np.nanpercentile(X, 0.1)
    smax = np.nanpercentile(abs.(shap_vals), 99.9)
    pdf = PdfPages("dnn.pdf")
    @showprogress for n in 1:size(Xviz, 1)
        x = permutedims(Xviz[n, :, :])
        sv = permutedims(shap_vals[n, :, :])
        figsize = [3 * (size(x, 1) + 1), 2.5 * (size(x, 2) + 1)]
        if figsize[1] > width
            @. figsize = width * figsize / figsize[1]
        end
        fig, axes = pl.subplots(ncols = 2, sharey = true, figsize = figsize)
        axes[1].imshow(x, cmap = pl.get_cmap("gray"), vmin = xmin, vmax = xmax)
        if !isempty(cols)
            axes[1].set_yticks(0:size(x, 1))
            axes[1].set_yticklabels(cols)
            labelsize = 5 * figsize[1] / size(x, 1)
            axes[1].tick_params(axis = "y", labelsize = labelsize, length = 0)
            axes[1].get_xaxis().set_visible(false)
        else
            axes[1].set_axis_off()
        end
        axes[2].imshow(x, cmap = pl.get_cmap("gray"), alpha = 0.15)
        axes[2].imshow(sv, cmap = colors.red_transparent_blue, vmin = -smax, vmax = smax)
        axes[2].set_axis_off()
        fig.tight_layout()
        pdf.savefig(fig, bbox_inches = "tight")
        pl.close(fig)
    end
    pdf.close()
end

function sample_clip(X, nsample, recep)
    N, T, F = size(X)
    recep = recep == 0 ? T : recep
    nsample = min(nsample, N * T ÷ recep)
    Xs = zeros(Float32, nsample, recep, F)
    for ns in 1:nsample
        n, t = rand(1:N), rand(1:(T - recep + 1))
        Xs[ns, :, :] = X[n, t:(t + recep - 1), :]
    end
    return Xs
end

function sample_rank(X, y, nsample)
    nsample = min(nsample, length(y))
    perm = sortperm(y, by = abs, rev = true)
    wv = Weights(1 ./ axes(y, 1)[invperm(perm)])
    ns = sample(1:length(y), wv, nsample, replace = false)
    return X[ns, :, :]
end

const deps = joinpath(@__DIR__, "../deps")

function install_brl()
    skater = mktempdir()
    clone("https://github.com/AStupidBear/Skater.git", skater)
    run(`sed -i "s|sudo python|$python|g" $skater/setup.sh`)
    run(`sed -i "s|.*install python.*||g" $skater/setup.sh`)
    run(`sed -i "s|apt-get install|apt-get install -y|g" $skater/setup.sh`)
    run(`bash -c "cd $skater && $python setup.py install --ostype=linux-ubuntu --rl=True"`)
end

end # module


text:
# ---
# title: 324. Wiggle Sort II
# id: problem324
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort
# link: <https://leetcode.com/problems/wiggle-sort-ii/description/>
# hidden: true
# ---
#
# Given an unsorted array `nums`, reorder it such that `nums[0] < nums[1] >
# nums[2] < nums[3]...`.
#
# **Example 1:**
#
#
#
#     Input:nums = [1, 5, 1, 1, 6, 4]
#     Output: One possible answer is [1, 4, 1, 5, 1, 6].
#
# **Example 2:**
#
#
#
#     Input:nums = [1, 3, 2, 2, 3, 1]
#     Output: One possible answer is [2, 3, 1, 3, 1, 2].
#
# **Note:**
# You may assume all input has valid answer.
#
# **Follow Up:**
# Can you do it in O(n) time and/or in-place with O(1) extra space?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file investigates 2D hiddden interpolation mean square error(MSE) with respect to inteprolation point number.

using FractalTools
using Makie

# Construct interpolation data
f(x, y) = [
    x^2 + y^2 + 1,
    x^2 - y^2,
]
vtx = [
    BigFloat.([0.0, 0.0]),
    BigFloat.([1.0, 0.0]),
    BigFloat.([0.5, 1.0])
    ]
freevar = 0.001
npts    = 50 : 5 : 150
ntpts   = 2 * npts[end]

# Construct test data
tpts = getdata(vtx, ntpts)

# Compute errors
mse = map(npts) do npt
    @info npt
    # Construct interpolation data
    pts = getdata(f, vtx, npt)

    # Construct interpolant
    interp = interpolate(pts, HInterp2D(fill(freevar, 2, 2)))

    # Compute error
    fvals = getindex.(map(pnt -> f(pnt...), tpts), 1)
    ivals = getindex.(map(pnt -> interp(pnt...), tpts), 1)
    sum((fvals - ivals).^2) / length(tpts)
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Number of Points", ylabel="MSE", title="2D Hidden Interpolation MSE")
stem!(ax, npts, mse, color=:black)
save(joinpath(@__DIR__, "hinterp2d_mse.png"), fig)
display(fig)


text:
module CarterBoyerLindquist

using Accessors
using Parameters
using DocStringExtensions
using StaticArrays

using SciMLBase
using DiffEqCallbacks

import GeodesicBase:
    AbstractMetricParams,
    inner_radius,
    AbstractGeodesicPoint,
    get_endpoint,
    geodesic_point_type,
    unpack_solution,
    SciMLBase

import GeodesicTracer:
    DiscreteCallback,
    terminate!,
    integrator_problem,
    metric_callback,
    create_callback_set,
    constrain,
    alpha_beta_to_vel

abstract type AbstractFirstOrderMetricParams{T} <: AbstractMetricParams{T} end

include("carter-method-bl-impl.jl")
include("carter-method-bl-interface.jl")
include("implementation.jl")
include("callbacks.jl")

"""
$(TYPEDSIGNATURES)

Radius of marginally stable orbit.

From Bardeen et al. (1972) eq. (2.21):

```math
r_\\text{ms} = M \\left\\{ 3 + Z_2 \\pm \\sqrt{(3 - Z_1)(3 + Z_1 + 2 Z_2)} \\right\\}.
```

The choice of ``\\pm`` is chosen by the sign of ``a``.
"""
rms(M, a, ±) = M * (3 + Z₂(M, a) ± √((3 - Z₁(M, a)) * (3 + Z₁(M, a) + 2 * Z₂(M, a))))
rms(M, a) = a > 0.0 ? rms(M, a, -) : rms(M, a, +)

end # module


text:
#=
Computes quantities related to the Gaussian linear state space model

    x_{t+1} = A x_t + C w_{t+1}

        y_t = G x_t

The shocks {w_t} are iid and N(0, I)

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-07-28

References
----------

TODO: Come back and update to match `LinearStateSpace` type from py side
TODO: Add docstrings

http://quant-econ.net/jl/linear_models.html

=#

"""
A type that describes the Gaussian Linear State Space Model
of the form:

    x_{t+1} = A x_t + C w_{t+1}

        y_t = G x_t

where {w_t} and {v_t} are independent and standard normal with dimensions
k and l respectively.  The initial conditions are mu_0 and Sigma_0 for x_0
~ N(mu_0, Sigma_0).  When Sigma_0=0, the draw of x_0 is exactly mu_0.

#### Fields

- `A::Matrix` Part of the state transition equation.  It should be `n x n`
- `C::Matrix` Part of the state transition equation.  It should be `n x m`
- `G::Matrix` Part of the observation equation.  It should be `k x n`
- `k::Int` Dimension
- `n::Int` Dimension
- `m::Int` Dimension
- `mu_0::Vector` This is the mean of initial draw and is of length `n`
- `Sigma_0::Matrix` This is the variance of the initial draw and is `n x n` and
                    also should be positive definite and symmetric

"""
type LSS{TSampler<:MVNSampler}
    A::Matrix
    C::Matrix
    G::Matrix
    k::Int
    n::Int
    m::Int
    mu_0::Vector
    Sigma_0::Matrix
    dist::TSampler
end

function LSS(A::ScalarOrArray, C::ScalarOrArray, G::ScalarOrArray,
             mu_0::ScalarOrArray,
             Sigma_0::Matrix=zeros(size(G, 2), size(G, 2)))
    k = size(G, 1)
    n = size(G, 2)
    m = size(C, 2)

    # coerce shapes
    A = reshape(vcat(A), n, n)
    C = reshape(vcat(C), n, m)
    G = reshape(vcat(G), k, n)

    mu_0 = reshape([mu_0;], n)

    dist = MVNSampler(mu_0,Sigma_0)
    LSS(A, C, G, k, n, m, mu_0, Sigma_0, dist)
end

# make kwarg version
function LSS(A::ScalarOrArray, C::ScalarOrArray, G::ScalarOrArray;
             mu_0::Vector=zeros(size(G, 2)),
             Sigma_0::Matrix=zeros(size(G, 2), size(G, 2)))
    return LSS(A, C, G, mu_0, Sigma_0)
end

function simulate(lss::LSS, ts_length=100)
    x = Array{Float64}(lss.n, ts_length)
    x[:, 1] = rand(lss.dist)
    w = randn(lss.m, ts_length - 1)
    for t=1:ts_length-1
        x[:, t+1] = lss.A * x[:, t] .+ lss.C * w[:, t]
    end
    y = lss.G * x

    return x, y
end

"""
Simulate num_reps observations of x_T and y_T given x_0 ~ N(mu_0, Sigma_0).

#### Arguments

- `lss::LSS` An instance of the Gaussian linear state space model.
- `t::Int = 10` The period that we want to replicate values for.
- `num_reps::Int = 100` The number of replications we want

#### Returns

- `x::Matrix` An n x num_reps matrix, where the j-th column is the j_th
              observation of x_T
- `y::Matrix` An k x num_reps matrix, where the j-th column is the j_th
              observation of y_T

"""
function replicate(lss::LSS, t::Integer, num_reps::Integer=100)
    x = Array{Float64}(lss.n, num_reps)
    for j=1:num_reps
        x_t, _ = simulate(lss, t+1)
        x[:, j] = x_t[:, end]
    end

    y = lss.G * x
    return x, y
end

replicate(lss::LSS; t::Integer=10, num_reps::Integer=100) =
    replicate(lss, t, num_reps)

immutable LSSMoments
    lss::LSS
end

Base.start(L::LSSMoments) = (copy(L.lss.mu_0), copy(L.lss.Sigma_0))
Base.done(L::LSSMoments, x) = false
function Base.next(L::LSSMoments, moms)
    A, C, G = L.lss.A, L.lss.C, L.lss.G
    mu_x, Sigma_x = moms

    mu_y, Sigma_y = G * mu_x, G * Sigma_x * G'

    # Update moments of x
    mu_x2 = A * mu_x
    Sigma_x2 = A * Sigma_x * A' + C * C'

    (mu_x, mu_y, Sigma_x, Sigma_y), (mu_x2, Sigma_x2)
end

"""
Create an iterator to calculate the population mean and
variance-convariance matrix for both x_t and y_t, starting at
the initial condition (self.mu_0, self.Sigma_0).  Each iteration
produces a 4-tuple of items (mu_x, mu_y, Sigma_x, Sigma_y) for
the next period.

#### Arguments

- `lss::LSS` An instance of the Gaussian linear state space model

"""
moment_sequence(lss::LSS) = LSSMoments(lss)

"""
Compute the moments of the stationary distributions of x_t and
y_t if possible.  Computation is by iteration, starting from the
initial conditions lss.mu_0 and lss.Sigma_0

#### Arguments

- `lss::LSS` An instance of the Guassian linear state space model
- `;max_iter::Int = 200` The maximum number of iterations allowed
- `;tol::Float64 = 1e-5` The tolerance level one wishes to achieve

#### Returns

- `mu_x::Vector` Represents the stationary mean of x_t
- `mu_y::Vector`Represents the stationary mean of y_t
- `Sigma_x::Matrix` Represents the var-cov matrix
- `Sigma_y::Matrix` Represents the var-cov matrix

"""
function stationary_distributions(lss::LSS; max_iter=200, tol=1e-5)
    # Initialize iteration
    m = moment_sequence(lss)
    mu_x, mu_y, Sigma_x, Sigma_y = first(m)

    i = 0
    err = tol + 1.0

    for (mu_x1, mu_y, Sigma_x1, Sigma_y) in m
        i > max_iter && error("Convergence failed after $i iterations")
        i += 1
        err_mu = maximum(abs, mu_x1 - mu_x)
        err_Sigma = maximum(abs, Sigma_x1 - Sigma_x)
        err = max(err_Sigma, err_mu)
        mu_x, Sigma_x = mu_x1, Sigma_x1

        if err < tol && i > 1
            break
        end
    end

    return mu_x, mu_y, Sigma_x, Sigma_y
end

function geometric_sums(lss::LSS, bet, x_t)
    I = eye(lss.n)
    S_x = (I - bet .* lss.A) \ x_t
    S_y = lss.G * S_x
    return S_x, S_y
end


text:
using IndexableBitVectors
using FactCheck

srand(12345)

function test_access{T}(::Type{T})
    b = convert(T, Bool[])
    @fact_throws b[0]
    @fact_throws b[1]

    b = convert(T, [true])
    @fact_throws b[0]
    @fact b[1] => 1
    @fact_throws b[2]

    b = convert(T, [false])
    @fact_throws b[0]
    @fact b[1] => 0
    @fact_throws b[2]

    b = convert(T, [false, false, true, true])
    @fact_throws b[0]
    @fact b[1] => 0
    @fact b[2] => 0
    @fact b[3] => 1
    @fact b[4] => 1
    @fact_throws b[5]

    b = convert(T, trues(1024))
    @fact b[1] => true
    @fact b[end] => true
    @fact all([b[i] for i in 1:1024]) => true

    b = convert(T, falses(1024))
    @fact b[1] => false
    @fact b[end] => false
    @fact all([b[i] for i in 1:1024]) => false

    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024; @fact b[i] => bitv[i]; end
end

function test_rank{T}(::Type{T})
    b = convert(T, Bool[])
    @fact rank0(b, 0) => 0
    @fact_throws rank0(b, 1)
    @fact rank1(b, 0) => 0
    @fact_throws rank1(b, 1)

    b = convert(T, [false])
    @fact rank0(b, 0) => 0
    @fact rank0(b, 1) => 1
    @fact_throws rank0(b, 2)
    @fact rank1(b, 0) => 0
    @fact rank1(b, 1) => 0
    @fact_throws rank1(b, 2)

    b = convert(T, [true])
    @fact rank0(b, 0) => 0
    @fact rank0(b, 1) => 0
    @fact_throws rank0(b, 2)
    @fact rank1(b, 0) => 0
    @fact rank1(b, 1) => 1
    @fact_throws rank1(b, 2)

    b = convert(T, [false, false, true, true])
    # rank0
    @fact rank0(b, 0) => 0
    @fact rank0(b, 1) => 1
    @fact rank0(b, 2) => 2
    @fact rank0(b, 3) => 2
    @fact rank0(b, 4) => 2
    @fact rank(0, b, 0) => 0
    @fact rank(0, b, 4) => 2
    @fact_throws rank(0, b, 5) "out of bound"
    # rank1
    @fact rank1(b, 0) => 0
    @fact rank1(b, 1) => 0
    @fact rank1(b, 2) => 0
    @fact rank1(b, 3) => 1
    @fact rank1(b, 4) => 2
    @fact rank(1, b, 0) => 0
    @fact rank(1, b, 4) => 2
    @fact_throws rank(1, b, 5) "out of bound"

    b = convert(T, trues(1024))
    for i in 1:1024; @fact rank0(b, i) => 0; end
    for i in 1:1024; @fact rank1(b, i) => i; end

    b = convert(T, falses(1024))
    for i in 1:1024; @fact rank0(b, i) => i; end
    for i in 1:1024; @fact rank1(b, i) => 0; end

    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024; @fact rank0(b, i) => rank0(bitv, i); end
    for i in 1:1024; @fact rank1(b, i) => rank1(bitv, i); end
end

function test_select{T}(::Type{T})
    b = convert(T, Bool[])
    @fact select0(b, 0) => 0
    @fact select1(b, 0) => 0

    b = convert(T, [false])
    @fact select0(b, 0) => 0
    @fact select0(b, 1) => 1
    @fact select0(b, 2) => 0
    @fact select1(b, 0) => 0
    @fact select1(b, 1) => 0
    @fact select1(b, 2) => 0

    b = convert(T, [true])
    @fact select0(b, 0) => 0
    @fact select0(b, 1) => 0
    @fact select0(b, 2) => 0
    @fact select1(b, 0) => 0
    @fact select1(b, 1) => 1
    @fact select1(b, 2) => 0

    b = convert(T, [false, false, true, true])
    # select0
    @fact select0(b, 0) => 0
    @fact select0(b, 1) => 1
    @fact select0(b, 2) => 2
    @fact select0(b, 3) => 0
    # select1
    @fact select1(b, 0) => 0
    @fact select1(b, 1) => 3
    @fact select1(b, 2) => 4
    @fact select1(b, 3) => 0

    b = convert(T, trues(1024))
    for i in 1:1024; @fact select0(b, i) => 0; end
    for i in 1:1024; @fact select1(b, i) => i; end

    b = convert(T, falses(1024))
    for i in 1:1024; @fact select0(b, i) => i; end
    for i in 1:1024; @fact select1(b, i) => 0; end

    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024; @fact select0(b, i) => select0(bitv, i); end
    for i in 1:1024; @fact select1(b, i) => select1(bitv, i); end
end

facts("BitVector") do
    context("access") do
        test_access(BitVector)
    end
    context("rank") do
        test_rank(BitVector)
    end
    context("select") do
        test_select(BitVector)
    end
end

facts("CompactBitVector") do
    context("access") do
        test_access(CompactBitVector)
    end
    context("rank") do
        test_rank(CompactBitVector)
    end
    context("select") do
        test_select(CompactBitVector)
    end
end

facts("SucVector") do
    context("access") do
        test_access(SucVector)
    end
    context("rank") do
        test_rank(SucVector)
    end
    context("select") do
        test_select(SucVector)
    end
end

facts("CSucVector") do
    context("access") do
        test_access(SucVector)
    end
    context("rank") do
        test_rank(CSucVector)
    end
    context("select") do
        test_select(SucVector)
    end
end

facts("RRR") do
    context("access") do
        test_access(RRR)
    end
    context("rank") do
        test_rank(RRR)
    end
    context("select") do
        test_select(RRR)
    end
end

facts("RRRNP") do
    context("access") do
        test_access(RRRNP)
    end
    context("rank") do
        test_rank(RRRNP)
    end
    context("select") do
        test_select(RRRNP)
    end
end


text:
using DeconvOptim
using Test
using FFTW, Noise, Statistics, Zygote
using Random

 # fix seed for reproducibility
Random.seed!(42)

@testset "Utils" begin
    include("utils.jl")
end

include("mappings.jl")

include("forward_models.jl")

include("lossfunctions.jl")

 # testing is rather hard, but include at least some basic testing
include("regularizer.jl")

include("main.jl")


text:
const TEST_FLOAT32 = false
const TEST_DOUBLE64 = false
const TEST_STATIC_SIZED = false

using Manifolds
using ManifoldsBase
using ManifoldsBase: number_of_coordinates
using FunManifolds

using LinearAlgebra
using Distributions
using DoubleFloats
using ForwardDiff
using Random
using ReverseDiff
using StaticArrays
using Test

using Manifolds: find_eps

# local methods
function Manifolds.find_eps(f1::Function, fs::Function...)
    return find_eps(f1(0.0), map(g -> g(0.0), fs)...)
end
function Manifolds.find_eps(
    f1::FunManifolds.VectorizedFunction,
    fs::FunManifolds.VectorizedFunction...,
)
    return find_eps(f1(0.0), map(g -> g(0.0), fs)...)
end

using Manifolds: get_iterator, _write, _read

function discretize(M::DiscretizedCurves, f)
    rep_size = representation_size(M.manifold)
    p = similar(f(M.grid[1]), rep_size..., length(M.grid))
    for i in get_iterator(M)
        copyto!(_write(M, rep_size, p, i), f(M.grid[i]))
    end
    return p
end


text:
module Lookup
using Combinatorics
using ..PokerCard

export get_lexographically_next_bit_sequence, flush_lookup, unsuited_lookup

# ref lookup.py in deuces

"""
Number of Distinct Hand Values:
Straight Flush   10
Four of a Kind   156      [(13 choose 2) * (2 choose 1)]
Full Houses      156      [(13 choose 2) * (2 choose 1)]
Flush            1277     [(13 choose 5) - 10 straight flushes]
Straight         10
Three of a Kind  858      [(13 choose 3) * (3 choose 1)]
Two Pair         858      [(13 choose 3) * (3 choose 2)]
One Pair         2860     [(13 choose 4) * (4 choose 1)]
High Card      + 1277     [(13 choose 5) - 10 straights]
-------------------------
TOTAL            7462
Here we create a lookup table which maps:
    5 card hand's unique prime product => rank in range [1, 7462]
Examples:
* Royal flush (best hand possible)          => 1
* 7-5-4-3-2 unsuited (worst hand possible)  => 7462
"""

###############################
# SOME CONSTANTS
###############################
MAX_STRAIGHT_FLUSH  = 10
MAX_FOUR_OF_A_KIND  = 166
MAX_FULL_HOUSE      = 322
MAX_FLUSH           = 1599
MAX_STRAIGHT        = 1609
MAX_THREE_OF_A_KIND = 2467
MAX_TWO_PAIR        = 3325
MAX_PAIR            = 6185
MAX_HIGH_CARD       = 7462
MAX_TO_RANK_CLASS = Dict(
    MAX_STRAIGHT_FLUSH => 1,
    MAX_FOUR_OF_A_KIND => 2,
    MAX_FULL_HOUSE => 3,
    MAX_FLUSH => 4,
    MAX_STRAIGHT => 5,
    MAX_THREE_OF_A_KIND => 6,
    MAX_TWO_PAIR => 7,
    MAX_PAIR => 8,
    MAX_HIGH_CARD => 9
)
RANK_CLASS_TO_STRING = Dict(
    1 => "Straight Flush",
    2 => "Four of a Kind",
    3 => "Full House",
    4 => "Flush",
    5 => "Straight",
    6 => "Three of a Kind",
    7 => "Two Pair",
    8 => "Pair",
    9 => "High Card"
)
flush_lookup = Dict()
unsuited_lookup = Dict()
###############################
# BITHACK FUNCTION
###############################
function get_lexographically_next_bit_sequence(bits)
    """
    Bit hack from here:
    http://www-graphics.stanford.edu/~seander/bithacks.html#NextBitPermutation
    Generator even does this in poker order rank
    so no need to sort when done! Perfect.
    """
    bits = bits%UInt
    t = (bits | (bits -1)) + 1
    a = (t & -t)
    b = UInt(a) / UInt((bits & -bits))
    c = (UInt(b) >> UInt(1)) - 1
    next = UInt(t) | UInt(c)
    return (next)
end
###############################
# STRAIGHTS AND HIGHCARDS
###############################
function straight_and_highcards(straights, highcards)
    """
    Unique five card sets. Straights and highcards.
    Reuses bit sequences from flush calculations.
    """
    ### STRAIGHTS ###
    # since straights are weaker than flushes, we start iterating from there
    rank = MAX_FLUSH + 1
    for s in straights
        # iterate through the straights
        prime_product = PokerCard.prime_product_from_rankbits(s)
        unsuited_lookup[prime_product] = rank
        rank += 1
    end

    ### HIGH CARDS ###
    # high cards are weaker than pairs, we'll iterate after the weakest pair
    rank = MAX_PAIR
    for h in highcards
        prime_product = PokerCard.prime_product_from_rankbits(h)
        unsuited_lookup[prime_product] = rank
        rank += 1
    end
end # straight_and_highcards
###############################
# FLUSHES
###############################
function flushes()
    """
    Straight flushes and flushes.
    Lookup is done on 13 bit integer (2^13 > 7462):
    xxxbbbbb bbbbbbbb => integer hand index
    """
    flushes = []
    global straight_flushes = [
        7936, # int('0b1111100000000', 2), # royal flush
        3968, # int('0b111110000000', 2),
        1984, # int('0b11111000000', 2),
        992, # int('0b1111100000', 2),
        496, # int('0b111110000', 2),
        248, # int('0b11111000', 2),
        124, # int('0b1111100', 2),
        62, # int('0b111110', 2),
        31, # int('0b11111', 2),
        4111 # int('0b1000000001111', 2) # 5 high flush
    ]
    # start the current bit sequence

    # 1277 = number of high cards (ie no pairs = all unique ranks)
    # 1277 + len(str_flushes) is number of hands with all cards unique rank

    for i in 1:(1277 + length(straight_flushes) -1 )
        if i==1
            v = parse(UInt, "11111"; base=2)
            global current_bit_sequence = get_lexographically_next_bit_sequence(v)
        else
            global current_bit_sequence = get_lexographically_next_bit_sequence(current_bit_sequence)
        end

        notSF = true
        for sf in straight_flushes
            # if current_bit_sequence XOR sf == 0
            # then we have a sf and should not add
            if current_bit_sequence ⊻ sf == false
                notSF = false
            end
        end
        if notSF
            push!(flushes, current_bit_sequence)
        end
    end
    # we started from the lowest straight pattern, now we want to start ranking from
    # the most powerful hands, so we reverse
    reverse!(flushes)
    # now add to the lookup map:
    # start with straight flushes and the rank of 1
    # since theyit is the best hand in poker
    # rank 1 = Royal Flush!
    rank  = 1
    for sf in straight_flushes
        prime_product = PokerCard.prime_product_from_rankbits(sf)
        flush_lookup[prime_product] = rank
        rank += 1
    end
    # we start the counting for flushes starting from  max full house, which
    # is the worst rank that a full house can have (2,2,2,3,3)
    rank = MAX_FULL_HOUSE + 1
    for f in flushes
        prime_product = PokerCard.prime_product_from_rankbits(f)
        flush_lookup[prime_product] = rank
        rank += 1
    end
    @assert length(flush_lookup) == 1287 # 1277 flushes + 10 straight FLUSHES

    # make sure we account for the overlaps
    straight_and_highcards(straight_flushes, flushes)
end # flushes
###############################
# MULTIPLES
###############################
function multiples()
    """
    Pair, Two Pair, Three of a Kind, Full House, and 4 of a Kind.
    """
    backwards_ranks = reverse!(collect(1:13))

    ### FOUR OF A KIND
    rank = MAX_STRAIGHT_FLUSH + 1
    for i in backwards_ranks
        kickers = copy(backwards_ranks)
        filter!(e -> e ≠ i, kickers)
        for k in kickers
            product = (PokerCard.PRIMES[i] ^ 4) * PokerCard.PRIMES[k]
            unsuited_lookup[product] = rank
            rank +=1
        end
    end

    ### FULL HOUSE
    rank = MAX_FOUR_OF_A_KIND + 1
    for i in backwards_ranks
        pairranks = copy(backwards_ranks)
        filter!(e -> e ≠ i, pairranks)
        for pr in pairranks
            product = (PokerCard.PRIMES[i] ^ 3) * PokerCard.PRIMES[pr] ^ 2
            unsuited_lookup[product] = rank
            rank += 1
        end
    end

    ### THREE OF A KIND
    rank = MAX_STRAIGHT + 1
    for r in backwards_ranks
        kickers = copy(backwards_ranks)
        filter!(e -> e ≠ r, kickers)
        kickers_combos = collect(combinations(kickers, 2))
        for kickers in kickers_combos
            c1, c2 = kickers
            product = PokerCard.PRIMES[r]^3 * PokerCard.PRIMES[c1] * PokerCard.PRIMES[c2]
            unsuited_lookup[product] = rank
            rank += 1
        end

    end
    ### TWO PAIR
    rank = MAX_THREE_OF_A_KIND
    tpgen = collect(combinations(backwards_ranks, 2))
    for tp in tpgen
        pair1, pair2 = tp
        kickers = copy(backwards_ranks)
        filter!(e -> e ∉ tp, kickers)
        for kicker in kickers
            product = (PokerCard.PRIMES[pair1]^2) * (PokerCard.PRIMES[pair2])^2 * (PokerCard.PRIMES[kicker])
            unsuited_lookup[product] = rank
            rank += 1
        end
    end
    ### PAIRS
    rank = MAX_TWO_PAIR + 1
    for pairrank in backwards_ranks
        kickers = copy(backwards_ranks)
        filter!(e -> e ≠ pairrank, kickers)
        kgen = collect(combinations(kickers, 3))
        for kickers in kgen
            k1, k2, k3 = kickers
            product = PokerCard.PRIMES[pairrank]^2 * PokerCard.PRIMES[k1] * PokerCard.PRIMES[k2] * PokerCard.PRIMES[k3]
            unsuited_lookup[product] = rank
            rank += 1
        end
    end
end # multiples

function write_table_to_disk(table, filepath)
    open(filepath, "w") do io
        for val in table
            prime_prod = val[1]
            rank = val[2]
            write(io, "$prime_prod, $rank \n")
        end
    end
end

flushes()
multiples()

# write the table to disk
# write_table_to_disk(flush_lookup, "Data/LOOKUP_TABLE_FLUSH.csv")
# write_table_to_disk(unsuited_lookup, "Data/LOOKUP_TABLE_UNSUITED.csv")

end # Lookup module

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
module Tricks

using Base: rewrap_unionall, unwrap_unionall, uncompressed_ast
using Base: CodeInfo

export static_hasmethod, static_methods, compat_hasmethod, static_fieldnames, static_fieldcount, static_fieldtypes

# This is used to create the CodeInfo returned by static_hasmethod.
_hasmethod_false(@nospecialize(f), @nospecialize(t)) = false
_hasmethod_true(@nospecialize(f), @nospecialize(t)) = true

"""
    static_hasmethod(f, type_tuple::Type{<:Tuple)

Like `hasmethod` but runs at compile-time (and does not accept a worldage argument).
"""
@generated function static_hasmethod(@nospecialize(f), @nospecialize(t::Type{T}),) where {T<:Tuple}
    # The signature type:
    world = typemax(UInt)
    method_insts = Core.Compiler.method_instances(f.instance, T, world)

    ftype = Tuple{f, fieldtypes(T)...}
    covering_method_insts = [mi for mi in method_insts if ftype <: mi.def.sig]

    method_doesnot_exist = isempty(covering_method_insts)
    ret_func = method_doesnot_exist ? _hasmethod_false : _hasmethod_true
    ci_orig = uncompressed_ast(typeof(ret_func).name.mt.defs.func)
    ci = ccall(:jl_copy_code_info, Ref{CodeInfo}, (Any,), ci_orig)

    # Now we add the edges so if a method is defined this recompiles
    if method_doesnot_exist
        # No method so attach to method table
        mt = f.name.mt
        typ = rewrap_unionall(Tuple{f, unwrap_unionall(T).parameters...}, T)
        ci.edges = Core.Compiler.vect(mt, typ)
    else  # method exists, attach edges to all instances
        ci.edges = covering_method_insts
    end
    return ci
end

function create_codeinfo_with_returnvalue(argnames, spnames, sp, value)
    expr = Expr(:lambda,
        argnames,
        Expr(Symbol("scope-block"),
            Expr(:block,
                Expr(:return, value),
            )
        )
    )
    if spnames !== nothing
        expr = Expr(Symbol("with-static-parameters"), expr, spnames...)
    end
    ci = ccall(:jl_expand, Any, (Any, Any), expr, @__MODULE__)
    return ci
end

"""
    static_methods(f, [type_tuple::Type{<:Tuple])
    static_methods(@nospecialize(f)) = _static_methods(Main, f, Tuple{Vararg{Any}})

Like `methods` but runs at compile-time (and does not accept a worldage argument).
"""
static_methods(@nospecialize(f)) = static_methods(f, Tuple{Vararg{Any}})
@generated function static_methods(@nospecialize(f) , @nospecialize(_T::Type{T})) where {T <: Tuple}
    list_of_methods = methods(f.instance, T)
    ci = create_codeinfo_with_returnvalue([Symbol("#self#"), :f, :_T], [:T], (:T,), :($list_of_methods))

    # Now we add the edges so if a method is defined this recompiles
    mt = f.name.mt
    ci.edges = Core.Compiler.vect(mt, Tuple{Vararg{Any}})
    return ci
end

@static if VERSION < v"1.3"
    const compat_hasmethod = hasmethod
else
    const compat_hasmethod = static_hasmethod
end

Base.@pure static_fieldnames(t::Type) = Base.fieldnames(t)
Base.@pure static_fieldtypes(t::Type) = Base.fieldtypes(t)
Base.@pure static_fieldcount(t::Type) = Base.fieldcount(t)

end  # module


text:
require("fcsc_types.jl");
require("fcsc_learning_multi.jl");
require("fcsc_coding_fastmulti.jl");
require("fcsc_fastmulti.jl");
require("fcsc_simple.jl");
require("normalize.jl");
require("helpers.jl");
require("lagrangian.jl");

d = zeros(Float64, 10, 10, 2);
d[1:3, 1:3, 1] = [1 2 3;2 3 1;3 2 1];
d[1:3, 1:3, 2] = rotr90([1 2 3;2 3 1;3 2 1]);

X1 = [0 0 0 0 1 1 0 0 0 0;
      0 0 0 0 1 1 0 0 0 0;
      0 0 0 0 1 1 0 0 0 0;
      1 1 1 1 1 1 1 1 1 1;
      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;
      0 0 0 0 1 1 0 0 0 0;
      0 0 0 0 1 1 0 0 0 0;
      0 0 0 0 1 1 0 0 0 0];
X2 = eye(10)+rotr90(eye(10));
X1 = float64(X1); X2 = float64(X2);

opts = FCSCOpts(Float64);
opts.max_iter = 10;
opts.FIXED_D = true; opts.Dinit = d;
(nil,z_seq) = fcsc_simple(cat(3, X1, X2), 2, 5, opts);
(nil,z_bat) = fcsc_fastmulti(cat(3, X1, X2), 2, 5, opts);

@printf("\n%e\n", norm(z_seq[:]-z_bat[:]))

opts = FCSCOpts(Float64);
opts.max_iter = 10;
opts.DEBUG = true;
@time (s1,) = fcsc_simple(cat(3, X1, X2), 2, 5, opts);
@time (s2,) = fcsc_fastmulti(cat(3, X1, X2), 2, 5, opts);

@printf("\n%e\n", norm(s1[:]-s2[:]))

max_stars_repo_name: bkong/FCSC


text:
"""
    regpower(X,Y,m)

adalah fungsi untuk mencari persamaan pangkat `y=Ax^m` dari variabel dependen `Y` dan
independen `X` dengan pangkat `m`.

# Example
```jldoctest
julia> tk = [0.2,0.4,0.6,0.8,1.0];

julia> dk = [0.1960,0.7850,1.7665,3.1405,4.9075];

julia> A = regpower(tk,dk,2)
4.907303370786516
```
return nilai koefisien `A`.
"""
function regpower(X,Y,m)
  sumxy = (X.^m)'*Y
  sumx2 = (X.^m)'*(X.^m)
  A = sumxy/sumx2
end


text:
elapsedtime(adf::AbstractDataFrame) = DataFrame(Average = mean(adf[:Time]))

function compare(fs::Vector, replications::Integer)
	n = length(fs)

	# Force JIT compilation
	for i in 1:n
		fs[i]()
	end

	times = Array(Float64, n * replications)
	indices = Array(Int, n * replications)
	for i in 1:n
		f = fs[i]
		f()
		for itr in 1:replications
			index = (i - 1) * replications + itr
			times[index] = @elapsed f()
			indices[index] = i
		end
	end

	df = DataFrame(Any[times, indices], [:Time, :Function])
	df = by(df, :Function, elapsedtime)
	df[:Relative] = df[:Average] / minimum(df[:Average])
	df[:Function] = DataArray(UTF8String, size(df, 1))
	for i in 1:size(df, 1)
		df[i, :Function] = string(fs[i])
	end
	df[:Replications] = replications

	return df
end

compare(replications::Integer, fs::Function...) = compare([fs...], replications)


text:
#=
An instance of `ILMSystem` holds the base cache, any problem-specific extra cache,
and containers for physical parameters, boundary condition data/functions, forcing data/functions,
a timestep function (which returns the timestep) and motion data/functions.

If the system needs to be regenerated because of surface motion, only the two
caches will need to be regenerated.
=#

Base.length(sys::ILMSystem) = length(sys.base_cache)

"""
    construct_system(prob::AbstractILMProblem) -> ILMSystem

Return a system of type of type `ILMSystem` from the given problem
instance `prob`.
"""
@inline construct_system(prob::AbstractILMProblem) = __init(prob)

"""
    update_system(sysold::ILMSystem,body::Body/BodyList)

From an existing system `sysold`, return a new system based on
the Body or BodyList `body`.
"""
function update_system(sysold::ILMSystem,body::Union{Body,BodyList})
    prob = regenerate_problem(sysold,body)
    return __init(prob)
end

"""
    update_system!(sys::ILMSystem,u,sysold,t)

From an existing system `sysold` at time `t`, return a new system `sys`
in place, based on the solution vector `u`, whose body state
information will be used to replace the body information and
subsequent operators in `sysold`.
"""
function update_system!(sys::ILMSystem,u,sysold::ILMSystem,t)
    @unpack base_cache, motions = sysold
    @unpack bl = base_cache
    x = aux_state(u)
    bodies = deepcopy(bl)
    update_body!(bodies,x,motions)
    sysnew = update_system(sysold,bodies)
    sys.base_cache = sysnew.base_cache
    sys.extra_cache = sysnew.extra_cache
    return sys
end

"""
    __init(prob::AbstractILMProblem)

Initialize `ILMSystem` with the given problem `prob` specification.
Depending on the type of problem, this sets up a base cache of scalar or
vector type, as well as an optional extra cache
"""
function __init(prob::AbstractILMProblem{DT,ST,DTP}) where {DT,ST,DTP}
    @unpack g, bodies, phys_params, bc, forcing, timestep_func, motions = prob

    base_cache = _construct_base_cache(bodies,g,DT,ST,DTP,prob,Val(length(bodies)))

    extra_cache = prob_cache(prob,base_cache)

    return ILMSystem{_static_surfaces(motions),typeof(prob),length(base_cache),typeof(phys_params),typeof(bc),typeof(forcing),
                    typeof(timestep_func),typeof(motions),typeof(base_cache),typeof(extra_cache)}(
              phys_params,bc,forcing,timestep_func,motions,base_cache,extra_cache)

end

@inline _construct_base_cache(bodies,g,DT,ST,DTP,::PT,::Val{0}) where {PT <: AbstractScalarILMProblem} =
              SurfaceScalarCache(g,ddftype=DT,scaling=ST,dtype=DTP)

@inline _construct_base_cache(bodies,g,DT,ST,DTP,::PT,::Val{N}) where {PT <: AbstractScalarILMProblem,N} =
              SurfaceScalarCache(bodies,g,ddftype=DT,scaling=ST,dtype=DTP)

@inline _construct_base_cache(bodies,g,DT,ST,DTP,::PT,::Val{0}) where {PT <: AbstractVectorILMProblem} =
              SurfaceVectorCache(g,ddftype=DT,scaling=ST,dtype=DTP)

@inline _construct_base_cache(bodies,g,DT,ST,DTP,::PT,::Val{N}) where {PT <: AbstractVectorILMProblem,N} =
              SurfaceVectorCache(bodies,g,ddftype=DT,scaling=ST,dtype=DTP)

_static_surfaces(::Nothing) = true
_static_surfaces(::Any) = false

# Extend surface_velocity! and allow for null motions
function RigidBodyTools.surface_velocity!(u::AbstractVector,v::AbstractVector,sys::ILMSystem,t)
  @unpack base_cache, motions = sys
  @unpack bl = base_cache
  fill!(u,0.0)
  fill!(v,0.0)
  if !isnothing(motions)
    surface_velocity!(u,v,bl,motions,t)
  end
end

RigidBodyTools.surface_velocity!(u::VectorData,sys::ILMSystem,t) = surface_velocity!(u.u,u.v,sys,t)

# Create the basic solve function, to be extended
function solve(prob::AbstractILMProblem,sys::ILMSystem) end

## Extend functions on `BasicILMCache` type to `ILMSystem`
for f in [:zeros_surface,:zeros_surfacescalar,
          :zeros_grid,:zeros_gridcurl,:zeros_griddiv,:zeros_gridgrad,:zeros_gridgradcurl,
          :similar_surface,:similar_surfacescalar,
          :similar_grid,:similar_gridcurl,:similar_griddiv,:similar_gridgrad,:similar_gridgradcurl,
          :ones_surface,:ones_surfacescalar,
          :ones_grid,:ones_gridgrad,:ones_gridcurl,:ones_griddiv,:ones_gridgradcurl,
          :x_grid,:y_grid,:x_gridcurl,:y_gridcurl,:x_griddiv,:y_griddiv,
          :normals,:areas,:points,:arcs,
          :create_nRTRn,:create_GLinvD,:create_CLinvCT,:create_CL2invCT,
          :create_RTLinvR,:create_GLinvD_symm,
          :create_GLinvD_cross,:create_surface_filter]
   @eval $f(sys::ILMSystem) = $f(sys.base_cache)
end

for f in [:regularize!, :interpolate!, :regularize_normal!,
          :normal_interpolate!,
          :regularize_normal_cross!,:normal_cross_interpolate!,
          :regularize_normal_dot!,:normal_dot_interpolate!,
          :regularize_normal_symm!,:normal_interpolate_symm!,
          :surface_curl!,:surface_divergence!,:surface_grad!,
          :surface_divergence_symm!,:surface_grad_symm!,
          :surface_curl_cross!,:surface_divergence_cross!,:surface_grad_cross!]
  @eval $f(a,b,sys::ILMSystem) = $f(a,b,sys.base_cache)
end

for f in [:norm,:integrate]
   @eval $f(a,sys::ILMSystem) = $f(a,sys.base_cache)
   @eval $f(a,sys::ILMSystem,i::Int) = $f(a,sys.base_cache,i)
end

for f in [:view]
  @eval $f(a::PointData,sys::ILMSystem,i::Int) = $f(a,sys.base_cache,i)
end

for f in [:dot,:copyto!]
   @eval $f(a,b,sys::ILMSystem) = $f(a,b,sys.base_cache)
   @eval $f(a,b,sys::ILMSystem,i::Int) = $f(a,b,sys.base_cache,i)
end

for f in [:RegularizationMatrix,:InterpolationMatrix]
  @eval CartesianGrids.$f(sys::ILMSystem,args...) = $f(sys.base_cache,args...)
end