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

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

Describes the status of the StartReportCreation operation.  You can call this operation
only from the organization's management account and from the us-east-1 Region.

"""
function describe_report_creation(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "DescribeReportCreation"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_report_creation(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "DescribeReportCreation",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns a table that shows counts of resources that are noncompliant with their tag
policies. For more information on tag policies, see Tag Policies in the AWS Organizations
User Guide.  You can call this operation only from the organization's management account
and from the us-east-1 Region. This operation supports pagination, where the response can
be sent in multiple pages. You should check the PaginationToken response parameter to
determine if there are additional results available to return. Repeat the query, passing
the PaginationToken response parameter value as an input to the next request until you
recieve a null value. A null value for PaginationToken indicates that there are no more
results waiting to be returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"GroupBy"`: Specifies a list of attributes to group the counts of noncompliant resources
  by. If supplied, the counts are sorted by those attributes.
- `"MaxResults"`: Specifies the maximum number of results to be returned in each page. A
  query can return fewer than this maximum, even if there are more results still to return.
  You should always check the PaginationToken response value to see if there are more
  results. You can specify a minimum of 1 and a maximum value of 100.
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
- `"RegionFilters"`: Specifies a list of AWS Regions to limit the output by. If you use
  this parameter, the count of returned noncompliant resources includes only resources in the
  specified Regions.
- `"ResourceTypeFilters"`: Specifies that you want the response to include information for
  only resources of the specified types. The format of each resource type is
  service[:resourceType]. For example, specifying a resource type of ec2 returns all Amazon
  EC2 resources (which includes EC2 instances). Specifying a resource type of ec2:instance
  returns only EC2 instances.  The string for each service name and resource type is the same
  as that embedded in a resource's Amazon Resource Name (ARN). Consult the AWS General
  Reference for the following:   For a list of service name strings, see AWS Service
  Namespaces.   For resource type strings, see Example ARNs.   For more information about
  ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.   You can specify
  multiple resource types by using a comma separated array. The array can include up to 100
  items. Note that the length constraint requirement applies to each resource type filter.
- `"TagKeyFilters"`: Specifies that you want the response to include information for only
  resources that have tags with the specified tag keys. If you use this parameter, the count
  of returned noncompliant resources includes only resources that have the specified tag keys.
- `"TargetIdFilters"`: Specifies target identifiers (usually, specific account IDs) to
  limit the output by. If you use this parameter, the count of returned noncompliant
  resources includes only resources with the specified target IDs.
"""
function get_compliance_summary(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetComplianceSummary"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_compliance_summary(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetComplianceSummary",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns all the tagged or previously tagged resources that are located in the specified
Region for the AWS account. Depending on what information you want returned, you can also
specify the following:    Filters that specify what tags and resource types you want
returned. The response includes all tags that are associated with the requested resources.
 Information about compliance with the account's effective tag policy. For more information
on tag policies, see Tag Policies in the AWS Organizations User Guide.    This operation
supports pagination, where the response can be sent in multiple pages. You should check the
PaginationToken response parameter to determine if there are additional results available
to return. Repeat the query, passing the PaginationToken response parameter value as an
input to the next request until you recieve a null value. A null value for PaginationToken
indicates that there are no more results waiting to be returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ExcludeCompliantResources"`: Specifies whether to exclude resources that are compliant
  with the tag policy. Set this to true if you are interested in retrieving information on
  noncompliant resources only. You can use this parameter only if the
  IncludeComplianceDetails parameter is also set to true.
- `"IncludeComplianceDetails"`: Specifies whether to include details regarding the
  compliance with the effective tag policy. Set this to true to determine whether resources
  are compliant with the tag policy and to get details.
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
- `"ResourceARNList"`: Specifies a list of ARNs of resources for which you want to retrieve
  tag data. You can't specify both this parameter and any of the pagination parameters
  (ResourcesPerPage, TagsPerPage, PaginationToken) in the same request. If you specify both,
  you get an Invalid Parameter exception. If a resource specified by this parameter doesn't
  exist, it doesn't generate an error; it simply isn't included in the response. An ARN
  (Amazon Resource Name) uniquely identifies a resource. For more information, see Amazon
  Resource Names (ARNs) and AWS Service Namespaces in the AWS General Reference.
- `"ResourceTypeFilters"`: Specifies the resource types that you want included in the
  response. The format of each resource type is service[:resourceType]. For example,
  specifying a resource type of ec2 returns all Amazon EC2 resources (which includes EC2
  instances). Specifying a resource type of ec2:instance returns only EC2 instances.  The
  string for each service name and resource type is the same as that embedded in a resource's
  Amazon Resource Name (ARN). Consult the AWS General Reference for the following: For more
  information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces. You
  can specify multiple resource types by using an array. The array can include up to 100
  items. Note that the length constraint requirement applies to each resource type filter.
- `"ResourcesPerPage"`: Specifies the maximum number of results to be returned in each
  page. A query can return fewer than this maximum, even if there are more results still to
  return. You should always check the PaginationToken response value to see if there are more
  results. You can specify a minimum of 1 and a maximum value of 100.
- `"TagFilters"`: Specifies a list of TagFilters (keys and values) to restrict the output
  to only those resources that have the specified tag and, if included, the specified value.
  Each TagFilter must contain a key with values optional. A request can include up to 50
  keys, and each key can include up to 20 values.  Note the following when deciding how to
  use TagFilters:   If you don't specify a TagFilter, the response includes all resources
  that are currently tagged or ever had a tag. Resources that currently don't have tags are
  shown with an empty tag set, like this: \"Tags\": [].   If you specify more than one filter
  in a single request, the response returns only those resources that satisfy all filters.
  If you specify a filter that contains more than one value for a key, the response returns
  resources that match any of the specified values for that key.   If you don't specify any
  values for a key, the response returns resources that are tagged with that key and any or
  no value. For example, for the following filters: filter1= {keyA,{value1}},
  filter2={keyB,{value2,value3,value4}}, filter3= {keyC}:    GetResources({filter1}) returns
  resources tagged with key1=value1     GetResources({filter2}) returns resources tagged with
  key2=value2 or key2=value3 or key2=value4     GetResources({filter3}) returns resources
  tagged with any tag with the key key3, and with any or no value
  GetResources({filter1,filter2,filter3}) returns resources tagged with (key1=value1) and
  (key2=value2 or key2=value3 or key2=value4) and (key3, any or no value)
- `"TagsPerPage"`: AWS recommends using ResourcesPerPage instead of this parameter. A limit
  that restricts the number of tags (key and value pairs) returned by GetResources in
  paginated output. A resource with no tags is counted as having one tag (one key and value
  pair).  GetResources does not split a resource and its associated tags across pages. If the
  specified TagsPerPage would cause such a break, a PaginationToken is returned in place of
  the affected resource and its tags. Use that token in another request to get the remaining
  data. For example, if you specify a TagsPerPage of 100 and the account has 22 resources
  with 10 tags each (meaning that each resource has 10 key and value pairs), the output will
  consist of three pages. The first page displays the first 10 resources, each with its 10
  tags. The second page displays the next 10 resources, each with its 10 tags. The third page
  displays the remaining 2 resources, each with its 10 tags. You can set TagsPerPage to a
  minimum of 100 items up to a maximum of 500 items.
"""
function get_resources(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetResources"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_resources(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetResources", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

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

Returns all tag keys currently in use in the specified Region for the calling AWS account.
This operation supports pagination, where the response can be sent in multiple pages. You
should check the PaginationToken response parameter to determine if there are additional
results available to return. Repeat the query, passing the PaginationToken response
parameter value as an input to the next request until you recieve a null value. A null
value for PaginationToken indicates that there are no more results waiting to be returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
"""
function get_tag_keys(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetTagKeys"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_tag_keys(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetTagKeys", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    get_tag_values(key)
    get_tag_values(key, params::Dict{String,<:Any})

Returns all tag values for the specified key that are used in the specified AWS Region for
the calling AWS account. This operation supports pagination, where the response can be sent
in multiple pages. You should check the PaginationToken response parameter to determine if
there are additional results available to return. Repeat the query, passing the
PaginationToken response parameter value as an input to the next request until you recieve
a null value. A null value for PaginationToken indicates that there are no more results
waiting to be returned.

# Arguments
- `key`: Specifies the tag key for which you want to list all existing values that are
  currently used in the specified AWS Region for the calling AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
"""
function get_tag_values(Key; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetTagValues",
        Dict{String,Any}("Key" => Key);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_tag_values(
    Key, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetTagValues",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Key" => Key), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_report_creation(s3_bucket)
    start_report_creation(s3_bucket, params::Dict{String,<:Any})

Generates a report that lists all tagged resources in the accounts across your organization
and tells whether each resource is compliant with the effective tag policy. Compliance data
is refreshed daily. The report is generated asynchronously. The generated report is saved
to the following location:
s3://example-bucket/AwsTagPolicies/o-exampleorgid/YYYY-MM-ddTHH:mm:ssZ/report.csv  You can
call this operation only from the organization's management account and from the us-east-1
Region.

# Arguments
- `s3_bucket`: The name of the Amazon S3 bucket where the report will be stored; for
  example:  awsexamplebucket  For more information on S3 bucket requirements, including an
  example bucket policy, see the example S3 bucket policy on this page.

"""
function start_report_creation(S3Bucket; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "StartReportCreation",
        Dict{String,Any}("S3Bucket" => S3Bucket);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_report_creation(
    S3Bucket,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return resource_groups_tagging_api(
        "StartReportCreation",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("S3Bucket" => S3Bucket), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resources(resource_arnlist, tags)
    tag_resources(resource_arnlist, tags, params::Dict{String,<:Any})

Applies one or more tags to the specified resources. Note the following:   Not all
resources can have tags. For a list of services with resources that support tagging using
this operation, see Services that support the Resource Groups Tagging API.   Each resource
can have up to 50 tags. For other limits, see Tag Naming and Usage Conventions in the AWS
General Reference.    You can only tag resources that are located in the specified AWS
Region for the AWS account.   To add tags to a resource, you need the necessary permissions
for the service that the resource belongs to as well as permissions for adding tags. For
more information, see the documentation for each service.    Do not store personally
identifiable information (PII) or other confidential or sensitive information in tags. We
use tags to provide you with billing and administration services. Tags are not intended to
be used for private or sensitive data.

# Arguments
- `resource_arnlist`: Specifies the list of ARNs of the resources that you want to apply
  tags to. An ARN (Amazon Resource Name) uniquely identifies a resource. For more
  information, see Amazon Resource Names (ARNs) and AWS Service Namespaces in the AWS General
  Reference.
- `tags`: Specifies a list of tags that you want to add to the specified resources. A tag
  consists of a key and a value that you define.

"""
function tag_resources(
    ResourceARNList, Tags; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "TagResources",
        Dict{String,Any}("ResourceARNList" => ResourceARNList, "Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resources(
    ResourceARNList,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return resource_groups_tagging_api(
        "TagResources",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("ResourceARNList" => ResourceARNList, "Tags" => Tags),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resources(resource_arnlist, tag_keys)
    untag_resources(resource_arnlist, tag_keys, params::Dict{String,<:Any})

Removes the specified tags from the specified resources. When you specify a tag key, the
action removes both that key and its associated value. The operation succeeds even if you
attempt to remove tags from a resource that were already removed. Note the following:   To
remove tags from a resource, you need the necessary permissions for the service that the
resource belongs to as well as permissions for removing tags. For more information, see the
documentation for the service whose resource you want to untag.   You can only tag
resources that are located in the specified AWS Region for the calling AWS account.

# Arguments
- `resource_arnlist`: Specifies a list of ARNs of the resources that you want to remove
  tags from. An ARN (Amazon Resource Name) uniquely identifies a resource. For more
  information, see Amazon Resource Names (ARNs) and AWS Service Namespaces in the AWS General
  Reference.
- `tag_keys`: Specifies a list of tag keys that you want to remove from the specified
  resources.

"""
function untag_resources(
    ResourceARNList, TagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "UntagResources",
        Dict{String,Any}("ResourceARNList" => ResourceARNList, "TagKeys" => TagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resources(
    ResourceARNList,
    TagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return resource_groups_tagging_api(
        "UntagResources",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ResourceARNList" => ResourceARNList, "TagKeys" => TagKeys
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
"""
    CompositeGP{Targs} <: AbstractGP

A GP derived from other GPs via an affine transformation. Specification given by `args`.
You should generally _not_ construct this object manually.
"""
struct CompositeGP{Targs} <: SthenoAbstractGP
    args::Targs
    n::Int
    gpc::GPC
    function CompositeGP{Targs}(args::Targs, gpc::GPC) where {Targs}
        gp = new{Targs}(args, next_index(gpc), gpc)
        gpc.n += 1
        return gp
    end
end
CompositeGP(args::Targs, gpc::GPC) where {Targs} = CompositeGP{Targs}(args, gpc)

mean(f::CompositeGP, x::AbstractVector) = mean(f.args, x)

cov(f::CompositeGP, x::AbstractVector) = cov(f.args, x)
var(f::CompositeGP, x::AbstractVector) = var(f.args, x)

cov(f::CompositeGP, x::AbstractVector, x′::AbstractVector) = cov(f.args, x, x′)
var(f::CompositeGP, x::AbstractVector, x′::AbstractVector) = var(f.args, x, x′)

function cov(
    f::SthenoAbstractGP, f′::SthenoAbstractGP, x::AbstractVector, x′::AbstractVector,
)
    @assert f.gpc === f′.gpc
    if f.n === f′.n
        return cov(f.args, x, x′)
    elseif f isa WrappedGP && f.n > f′.n || f′ isa WrappedGP && f′.n > f.n
        return zeros(length(x), length(x′))
    elseif f.n >= f′.n
        return cov(f.args, f′, x, x′)
    else
        return cov(f, f′.args, x, x′)
    end
end

function var(
    f::SthenoAbstractGP, f′::SthenoAbstractGP, x::AbstractVector, x′::AbstractVector,
)
    @assert f.gpc === f′.gpc
    if f.n === f′.n
        return var(f.args, x, x′)
    elseif f isa WrappedGP && f.n > f′.n || f′ isa WrappedGP && f′.n > f.n
        return zeros(length(x))
    elseif f.n >= f′.n
        return var(f.args, f′, x, x′)
    else
        return var(f, f′.args, x, x′)
    end
end

ext: jl


text:
#======================================================================================
This is a set of functions for consistently tokenizing categorical values.

TODO for now there is no special handing of nulls

Obviously I'm no longer going to use these at some point...
======================================================================================#

"""
    tokenmap(A[, tokentype=Int64, sort_values=true])
    tokenmap(data, col[, tokentype=Int64; sort_values=true])

Tokenizes elements in the `Vector` or `NullableVector` `A` by creating a map between distinct
elements of `A` and tokens of type `tokentype`.  The mapping will be returned in the form
of a `Dict`.  If `sort_values`, the distinct values of `A` will be sorted before
tokenization (so that the map is order-preserving).

This can also be called on a column of a dataframe.

**TODO:** Note that currently this doesn't do anything special to handle nulls.
"""
function tokenmap{T,K}(A::Vector{T}, ::Type{K}=Int64;
                       sort_values::Bool=false)
    u = unique(A)
    sort_values && sort!(u)
    t = zero(K)
    mapping = Dict{T,K}();  sizehint!(mapping, length(u))
    for i ∈ 0:(length(u)-1)
        mapping[u[i+1]] = t
        t += one(K)
    end
    mapping
end

function tokenmap{T,K}(A::NullableVector{T}, ::Type{K}=Int64;
                       sort_values::Bool=false)
    A = convert(Vector{Nullable{T}}, A)
    tokenmap(A, K, sort_values=sort_values)
end

function tokenmap{K}(data::DataFrame, col::Symbol, ::Type{K}=Int64;
                     sort_values::Bool=false)
    tokenmap(data[col], K, sort_values=sort_values)
end
export tokenmap

"""
    tokenmaps(data, cols[, tokentype=Int64, sort_values=false])

Creates a token map for each column of the `DataFrame` `data` in `cols`, and returns
them in a `Dict`  of the form `col=>tokenmap` (see documentation for `tokenmap`).
"""
function tokenmaps{K}(data::DataFrame, cols::Vector{Symbol}, ::Type{K}=Int64;
                      sort_values::Bool=false)
    Dict(col=>tokenmap(data, col, K, sort_values=sort_values) for col ∈ cols)
end
export tokenmaps

"""
    tokenize(A, mapping)
    tokenize(A[, tokentype=Int64, sort_values=true])

Returns a vector which is the result of applying the map returned by `tokenmap` (see
documentation for that function) to elements of `A`.  This mapping can either be passed
explicitly as `mapping` or it can be generated if it is not passed.
"""
function tokenize{T,K}(A::Vector{T}, mapping::Dict{T,K})
    K[mapping[a] for a ∈ A]
end

function tokenize{T,K}(A::NullableVector{T}, mapping::Dict{Nullable{T},K})
    K[mapping[a] for a ∈ A]
end

function tokenize{T,K}(A::Union{Vector{T},NullableVector{T}}, ::Type{K}=Int64;
                       sort_values::Bool=false)
    mapping = tokenmap(A, K, sort_values=sort_values)
    tokenize(A, mapping)
end
export tokenize

"""
    detokenize(A, mapping[; is_inverse=false])

Inverts the tokenization of `A` (see `tokenmap` and `tokenize`).  The mapping provided is
assumed to be the forward mapping (from original objects to tokens) unless `is_inverse` in
which case it is assumed to be the inverse mapping (i.e. from tokens to original objects).
"""
function detokenize{T,K}(A::Union{Vector{T},NullableVector{T}}, mapping::Dict{K,T};
                         is_inverse::Bool=false)
    if !is_inverse
        inv_mapping = invert(mapping)
    else
        inv_mapping = mapping
    end
    T[inv_mapping[a] for a ∈ A]
end
export detokenize


text:
using QuNet
using Statistics

grid_size = 100
num_pairs = 50
num_trials = 100
max_paths = 4

failure_counts = []
for i in 1:num_trials
    println("num_trial $i")
    edge_costs = Dict("loss"=>0.05, "Z"=>0.05)
    net = GridNetwork(grid_size, grid_size, edge_costs=edge_costs)
    user_pairs = QuNet.make_user_pairs(net, 50)
    pathset, pur_paths, pathuse_count = QuNet.greedy_multi_path!(net, QuNet.purify, user_pairs, 4)
    push!(failure_counts, pathuse_count[1])
end
println(mean(failure_counts./num_pairs))
println(failure_counts)

# # Set up a Python backend before importing PyPlot
# using PyCall
# # Specify user GUI from options (:tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx)
# # pygui(:qt)
# using PyPlot
#
# # # Test PyPlot
# # # use x = linspace(0,2*pi,1000) in Julia 0.6
# # x = range(0; stop=2*pi, length=1000); y = sin.(3 * x + 4 * cos.(2 * x));
# # plot(x, y, color="red", linewidth=2.0, linestyle="--")
# # title("A sinusoidally modulated sinusoid")
# # # Display figure with Julia backend
# # display(gcf())
# # # Save as a Pdf
# # savefig("sinusoid_test.pdf")
# # #show()
#
# # Try heatmap
# # Install numpy
# using PyCall
# np = pyimport("numpy")
# kde = pyimport("scipy.stats.kde")
# data = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 3]], 200)
# Tdata = transpose(data)
# x = Tdata[1, :]
# y = Tdata[2, :]
#
# # # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
# # k = kde.gaussian_kde(Tdata)
# # xi, yi = np.mgrid[minimum(x):maximum(x):nbins*1j, minimum(y):maximum(y):nbins*1j]
# # zi = k(np.vstack([xi.flatten(), yi.flatten()]))
#
# fig, axs = plt.subplots(ncols=2, nrows=1, figsize=(21, 5))
# # Everything starts with a Scatterplot
# axs[0].set_title('Scatterplot')
# axs[0].plot(x, y, 'ko')
#
# # # plot a density
# # axes[1].set_title('Calculate Gaussian KDE')
# # axes[1].pcolormesh(xi, yi, zi.reshape(xi.shape), shading='auto', cmap=plt.cm.BuGn_r)
# #
# display(gcf())

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
import Random
Random.seed!(1337)

println("\n")

sp = sptenrand([10, 10, 10], 0.01)
dense = Array(sp)

println("... Testing sparse and dense tensor product equivalence.")
@test ttt(sp,sp) == ttt(dense,dense)

println("... Testing sparse tensor decomposition")
@test cp_als(reduce(ttt,SparseTensor.([[0, 0, 1],[0, 1, 0],[1, 0, 0]])),1) == ktensor([[0 0 1],[0 1 0],[1 0 0]].|>permutedims .|> x->x .|> Float64)

println("... Testing sparse tensor broadcast")
@test SparseTensor([0 0 1]) .* SparseTensor([1 0 0]) == [0 0 0]


text:
function bounded_sum(x::Int64,y::Int64,bound::Int64)
	return x + y < bound ? x + y : bound
end

function bounded_diff(x::Int64,y::Int64,bound::Int64)
	return x - y > bound ? x - y : bound
end

function random_array_fixed_sum(sum_value::Int64,array_size::Int64)
	if sum_value == 0
		return [ 0 for i in 1:array_size ]
	end

	total = sum_value
	array = []
	for i in 1:array_size-1
		push!(array,rand(1:total - (array_size - i) ))
		total -= array[i]
	end
	push!(array,total)

	return array
end

function node_distance(p1::Tuple,p2::Tuple)
	return round(evaluate(Euclidean(),collect(p1),collect(p2)),digits=3)
end

function random_elements_from_array(S::Array,num::Int)
	if length(S) <= num
		return S
	end

	V = copy(S)
	chosed = []
	while length(chosed) < num
		e = rand(S)
		setdiff!(S,[e])
		push!(chosed,e)
	end
	return chosed
end


text:
"""
    LanczosBasis(n) -> LanczosBasis(n,identity,one,L2(n))
    LanczosBasis(n,a,b,p::InnerProduct)

Orthogonal basis for `span{b(x), a(x)*b(x), ..., a(x)^(n-1)*b(x)}`.
"""
struct LanczosBasis{T,A,B,P} <: AbstractBasis
    length::Int
    a_fun::A
    b_fun::B
    inner_product::P

    # Recursion coefficients
    invnorms::Vector{T}
    dots1::Vector{T}
    dots2::Vector{T}
end

LanczosBasis(n) = LanczosBasis(n,identity,one,L2(n))
function LanczosBasis(n,a_fun,b_fun,p::InnerProduct)
    a = evaluate(a_fun,p)
    b = evaluate(b_fun,p)

    T = typeof(dot(a,b,p))
    invnorms = Vector{T}(undef,n)
    dots1 = Vector{T}(undef,max(0,n-1))
    dots2 = Vector{T}(undef,max(0,n-2))

    n == 0 && return LanczosBasis(n,a_fun,b_fun,p, invnorms,dots1,dots2)
    invnorms[1] = inv(norm(b,p))
    n == 1 && return LanczosBasis(n,a_fun,b_fun,p, invnorms,dots1,dots2)
    q1 = invnorms[1].*b

    q̃2 = a.*q1
    dots1[1] = dot(q̃2,q1,p)
    q̃2 = q̃2 .- dots1[1].*q1
    dots2_end = norm(q̃2,p)
    invnorms[2] = inv(dots2_end)
    q2 = invnorms[2].*q̃2

    for i = 3:n
        dots2[i-2] = dots2_end
        q̃3 = a.*q2
        dots1[i-1] = dot(q̃3,q2,p)
        q̃3 = q̃3 .- dots1[i-1].*q2 .- dots2[i-2].*q1
        dots2_end = norm(q̃3,p)
        invnorms[i] = inv(dots2_end)
        q3 = invnorms[i].*q̃3
        q1,q2 = q2,q3
    end
    return LanczosBasis(n,a_fun,b_fun,p, invnorms,dots1,dots2)
end

Base.length(B::LanczosBasis) = B.length

approxpoints(B::LanczosBasis) = evalpoints(B.inner_product)
approxtransform(B::LanczosBasis,f) = dot(B,f,w)

function iterate_basis(B::LanczosBasis, x)
    length(B) == 0 && return nothing
    q1 = B.invnorms[1]*B.b_fun(x)
    return q1,(2,B.a_fun(x),q1,q1)
end

function iterate_basis(B::LanczosBasis, x, (i,a,q2,q1))
    i > length(B) && return nothing
    if i == 2
        q2 = B.invnorms[i]*(a*q1 - B.dots1[i-1]*q1)
    else
        q1,q2 = q2, B.invnorms[i]*(a*q2 - B.dots1[i-1]*q2 - B.dots2[i-2]*q1)
    end
    return q2,(i+1,a,q2,q1)
end


text:
# Atmos RHS

#=
Ekman column:
    ∂_t ρ =  ∇ (μ ∇ ρ - w ρ)
    ∂_t ρθ =  ∇ (μ ∇ ρθ - w ρθ)
    ∂_t u =  ∇ (μ ∇ u - w u) + (v - v_g)
    ∂_t v =  ∇ (μ ∇ v - w v) - (u - u_g)
    ∂_t w =  ∇ (μ ∇ w - w w) - g - c_p θ ∂_z Π

where
    Π = (p/p_0)^{R/c_p}

top BCs are insulating and impenetrable:
    ∂_t T = 0
    u = u_g
    v = v_g
    w = 0.0
    ∂_t ρ = 0

and bottom BCs use bulk formulae for surface fluxes of heat and momentum:
    ∂_t ρθ = F₃ = -Ch ρ ||u|| (T_sfc - ρθ / ρ)
    ∂_t u  = F₁ = -Cd u ||u||
    ∂_t v  = F₂ = -Cd v ||u||
    w = 0.0
    ∂_t ρ = 0

We also use this model to accumulate fluxes it calculates
    ∂_t F_accum = -(F₁, F₂, F₃)
=#

function ∑tendencies_atm!(dY, Y, (parameters, T_sfc), t)

    UnPack.@unpack Ch, Cd, f, ν, ug, vg, C_p, MSLP, R_d, R_m, C_v, grav = parameters

    #@show(t, Y.x[3])
    (Yc, Yf, F_sfc) = Y.x
    (dYc, dYf, dF_sfc) = dY.x

    UnPack.@unpack ρ, u, v, ρθ = Yc
    UnPack.@unpack w = Yf
    dρ = dYc.ρ
    du = dYc.u
    dv = dYc.v
    dρθ = dYc.ρθ
    dw = dYf.w

    # Auxiliary calculations
    u_1 = parent(u)[1]
    v_1 = parent(v)[1]
    ρ_1 = parent(ρ)[1]
    ρθ_1 = parent(ρθ)[1]
    u_wind = sqrt(u_1^2 + v_1^2)

    # surface flux calculations
    surface_flux_ρθ =
        -calculate_sfc_fluxes_energy(
            DryBulkFormulaWithRadiation(),
            parameters,
            T_sfc[1],
            parent(ρθ)[1] / parent(ρ)[1],
            u_1,
            v_1,
            ρ_1,
            t,
        ) ./ C_p
    surface_flux_u = -Cd * u_1 * sqrt(u_1^2 + v_1^2)
    surface_flux_v = -Cd * v_1 * sqrt(u_1^2 + v_1^2)

    # accumulate in the required right units
    @inbounds begin
        dY.x[3][1] = -ρ_1 * surface_flux_u  #
        dY.x[3][2] = -ρ_1 * surface_flux_v  #
        dY.x[3][3] = -C_p * surface_flux_ρθ # W / m^2
    end

    # @inbounds begin
    #     dY.x[3][1] = - 10.0
    #     dY.x[3][2] = - 1.0
    #     dY.x[3][3] = - 1.0
    # end

    # Density tendency (located at cell centers)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))

    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    @. dρ = gradf2c(-w * If(ρ)) # Eq. 4.11

    # Potential temperature tendency (located at cell centers)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(surface_flux_ρθ), top = Operators.SetValue(0.0)) # Eq. 4.20, 4.21

    @. dρθ = gradf2c(-w * If(ρθ) + ν * gradc2f(ρθ / ρ)) # Eq. 4.12

    # u velocity tendency (located at cell centers)
    gradc2f = Operators.GradientC2F(top = Operators.SetValue(ug)) # Eq. 4.18
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(Cd * u_wind * u_1)) # Eq. 4.16

    A = Operators.AdvectionC2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    @. du = gradf2c(ν * gradc2f(u)) + f * (v - vg) - A(w, u) # Eq. 4.8

    # v velocity (centers)
    gradc2f = Operators.GradientC2F(top = Operators.SetValue(vg)) # Eq. 4.18
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(Cd * u_wind * v_1)) # Eq. 4.16

    A = Operators.AdvectionC2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    @. dv = gradf2c(ν * gradc2f(v)) - f * (u - ug) - A(w, v) # Eq. 4.9

    # w velocity (faces)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))

    B = Operators.SetBoundaryOperator(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    println(R_d ./ C_p)
    println(ρθ)
    println(t)
    Π(ρθ) = C_p .* (R_d .* ρθ ./ MSLP) .^ (R_m ./ C_v)

    @. dw = B(-(If(ρθ / ρ) * gradc2f(Π(ρθ))) - grav + gradc2f(ν * gradf2c(w)) - w * If(gradf2c(w))) # Eq. 4.10 # this makes everything unstable... use new ClimaAtmos rhs!

    return dY
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
#!/usr/bin/env julia

using Luxor, Colors

using Test

using Random
Random.seed!(42)

function spiral_logo_eps()
    gsave()
    scale(.3, .3)
    r = 200
    setcolor("gray")
    for i in 0:pi/8:2pi
        gsave()
        translate(r * cos(i), r * sin(i))
        rotate(i)
        julialogo()
        grestore()
    end
    grestore()
end

function expandingspiral_eps()
    gsave()
    scale(.3, .3)
    r = 200
    for i in pi:pi/12:6pi
        gsave()
        translate(i/3 * r * cos(i), i/3 * r * sin(i))
        scale(0.8, 0.8)
        rotate(i)
        julialogo()
        grestore()
    end
    grestore()
end

function dropshadow_eps()
    steps=20
    # white-gray ramp
    gramp = range(colorant"white", stop=colorant"gray60", length=steps)
    gsave()
    r = 200
    setopacity(0.1)
    for i in 1:steps
        sethue(gramp[i])
        translate(-0.6, -0.5)
        julialogo(color=false)
    end
    julialogo()
    grestore()
end

function colorgrid_eps()
    #cols = colormap("RdBu", 5; mid=0.5, logscale=false)
    #cols = sequential_palette(rand(10:360), 5, b=0.1)
    cols = distinguishable_colors(25)
    gsave()
    c = 0
    for row in 100:100:500
        for column in 100:100:500
            gsave()
            setcolor(color(cols[c+=1]))
            translate(row, column)
            scale(0.3, 0.3)
            julialogo(color=false)
            grestore()
        end
    end
    grestore()
end

function draw_logo(fname)
    Drawing(1600, 1600, fname)
    origin()
    background("white")

    translate(-500, -200)
    spiral_logo_eps()

    translate(750, 0)
    expandingspiral_eps()

    translate(-1000, 500)
    dropshadow_eps()

    translate(700, -100)
    colorgrid_eps()

    @test finish() == true
    println("...finished test: output in $(fname)")
end

draw_logo("julia-logo-draw-eps.eps")


text:
using Base.Threads: @spawn
using Base.Iterators: drop
using ProgressMeter: @showprogress
using NamedTupleTools: namedtuple

read_parquet(path, cols::Vector{Symbol}; kwargs...) = read_parquet(path, String.(cols); kwargs...)

read_parquet(path; kwargs...) = read_parquet(path, String[]; kwargs...)

function read_parquet(path, cols::Vector{String}; multithreaded=true, verbose = false)
	"""function for reading parquet"""

	if multithreaded
		# use a bounded channel to limit
	    c1 = Channel{Bool}(Threads.nthreads())
	    atexit(()->close(c1))
	end

	nc = ncols(ParFile(path))

	colnames = [sch.name for sch in  drop(ParFile(path).schema.schema, 1)]

	if length(cols) == 0
		colnums = collect(1:nc)
	else
		colnums = [findfirst(==(c), colnames) for c in cols]
	end

	results = Vector{Any}(undef, length(colnums))

	filemetadata = metadata(path)

	if multithreaded
		@showprogress for (i, j) in enumerate(colnums)
			put!(c1, true)
			results[i] = @spawn begin
				res = read_column(path, filemetadata, j)
				take!(c1)
				res
			end
		end
	else
		@showprogress for (i, j) in enumerate(colnums)
			results[i] = read_column(path, filemetadata, j)
		end
	end

	symbol_col_names = collect(Symbol(col) for col in colnames[colnums])

	if multithreaded
		fnl_results = collect(fetch(result) for result in results)
		return namedtuple(symbol_col_names, fnl_results)
	else
		return namedtuple(symbol_col_names, results)
	end
end


text:
# Autogenerated wrapper script for isoband_jll for i686-w64-mingw32
export libisoband

JLLWrappers.@generate_wrapper_header("isoband")
JLLWrappers.@declare_library_product(libisoband, "libisoband.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libisoband,
        "bin\\libisoband.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    CWaylandClientCore

Wrapper around libwayland-client.

Implemented differently:
- naming and calling convention - Julia-lised the functions to have shorter names and be distinguished by different object overloads. Eg. `wl_display_disconnect(struct wl_display *display)` -> `disconnect(display::Ptr{WlDisplay})`
- array marshalling functions - their library implementations require C unions and are more trouble then they're worth. Simply made them splat the array and call the non-array versions.

Not implemented:
- log handler
"""
module CWaylandClientCore

abstract type WlDisplay end
abstract type WlProxy end
abstract type WlProxyReal <: WlProxy end
abstract type WlProxyWrapper <: WlProxy end
abstract type WlEventQueue end

# Display
"Connect to the default display."
function connect()
    ccall((:wl_display_connect, "libwayland-client"), Ptr{WlDisplay}, (Ptr{Cvoid},), C_NULL)
end
"Connect to a named display."
function connect(path::String)
    ccall((:wl_display_connect, "libwayland-client"), Ptr{WlDisplay}, (Cstring,), path)
end
"Connect to a display via a file descriptor."
function connect(fd::RawFD)
    ccall((:wl_display_connect_to_fd, "libwayland-client"), Ptr{WlDisplay}, (RawFD,), fd)
end
"Disconnect from `display`."
function disconnect!(display::Ptr{WlDisplay})
    ccall((:wl_display_disconnect, "libwayland-client"), Cvoid, (Ptr{WlDisplay},), display)
end
"Get `display`'s file descriptor."
function getfd(display::Ptr{WlDisplay})
    ccall((:wl_display_get_fd, "libwayland-client"), RawFD, (Ptr{WlDisplay},), display)
end
"Get a local error."
function geterror(display::Ptr{WlDisplay})
    ccall((:wl_display_get_error, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Get a remote error."
function geterror(display::Ptr{WlDisplay}, interface::Ref{WlInterface}, id::Integer)
    ccall((:wl_display_get_protocol_error, "libwayland-client"), UInt32,
    (Ptr{WlDisplay}, Ref{WlInterface}, UInt32),
    display, interface, id)
end
"Flush all requests to display."
function flush(display::Ptr{WlDisplay})
    ccall((:wl_display_flush, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Block until all requests have been processed by `display`."
function roundtrip(display::Ptr{WlDisplay})
    ccall((:wl_display_roundtrip, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Block until all requests on `queue` have been processed by `display`."
function roundtrip(display::Ptr{WlDisplay}, queue::Ptr{WlEventQueue})
    ccall((:wl_display_roundtrip_queue, "libwayland-client"), Cint,
    (Ptr{WlDisplay}, Ptr{WlEventQueue}),
    display, queue))
end
"Start preparing all events on `display` for reading."
function prepare_read(display::Ptr{WlDisplay})
    ccall((:wl_display_prepare_read, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Start preparing all events on `display`'s' `queue` for reading."
function prepare_read(display::Ptr{WlDisplay}, queue::Ptr{WlEventQueue})
    ccall((:wl_display_prepare_read_queue, "libwayland-client"), Cint,
    (Ptr{WlDisplay}, Ptr{WlEventQueue}),
    display, queue))
end
"Cancel the prepared read on `display`."
function cancel_read(display::Ptr{WlDisplay})
    ccall((:wl_display_cancel_read, "libwayland-client"), Cvoid, (Ptr{WlDisplay},), display)
end
"Read the prepared events on `display`."
function read_events(display::Ptr{WlDisplay})
    ccall((:wl_display_read_events, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Dispatch event handlers on all events for `display`. Block if there are none."
function dispatch(display::Ptr{WlDisplay})
    ccall((:wl_display_dispatch, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Dispatch event handlers on all events for `display`'s `queue`. Block if there are none."
function dispatch(display::Ptr{WlDisplay}, queue::Ptr{WlEventQueue})
    ccall((:wl_display_dispatch_queue, "libwayland-client"), Cint,
    (Ptr{WlDisplay}, Ptr{WlEventQueue}),
    display, queue)
end
"Dispatch event handlers on all events for `display`. Skip if there are none."
function dispatch_pending(display::Ptr{WlDisplay})
    ccall((:wl_display_dispatch_pending, "libwayland-client"), Cint, (Ptr{WlDisplay},), display)
end
"Dispatch event handlers on all events for `display`'s `queue`. Skip if there are none."
function dispatch_pending(display::Ptr{WlDisplay}, queue::Ptr{WlEventQueue})
    ccall((:wl_display_dispatch_queue_pending, "libwayland-client"), Cint,
    (Ptr{WlDisplay}, Ptr{WlEventQueue}),
    display, queue)
end
# Event queue
"Create a `WlEventQueue`."
function create_queue(display::Ptr{WlDisplay})
    ccall((:wl_event_queue_destroy, "libwayland-client"), Ptr{WlEventQueue}, (Ptr{WlDisplay},), display)
end
"Destroy a `WlEventQueue`."
function destroy!(queue::Ptr{WlEventQueue})
    ccall((:wl_event_queue_destroy, "libwayland-client"), Cvoid, (Ptr{WlEventQueue},), queue)
end
# Proxy
"Create a `WlProxy`."
function create_proxy(factory::Ptr{<: WlProxy}, interface::Ref{WlInterface})
    ret = ccall((:wl_proxy_create, "libwayland-client"), Ptr{WlProxy},
    (Ptr{WlProxy}, Ref{WlInterface}),
    factory, interface)
    if ret == C_NULL
        error("Failure creating proxy.")
    else
        return ret
    end
end
"Destroy a `WlProxy`."
function destroy!(proxy::Ptr{WlProxyReal})
    ccall((:wl_proxy_destroy, "libwayland-client"), Cvoid, (Ptr{WlProxy},), proxy)
end
"Get `proxy`'s version."
function getversion(proxy::Ptr{<: WlProxy})
    ccall((:wl_proxy_get_version, "libwayland-client"), UInt32, (Ptr{WlProxy},), proxy)
end
"Get `proxy`'s ID."
function getid(proxy::Ptr{<: WlProxy})
    ccall((:wl_proxy_get_id, "libwayland-client"), UInt32, (Ptr{WlProxy},), proxy)
end
"Get `proxy`'s class."
function getclass(proxy::Ptr{<: WlProxy})
    ccall((:wl_proxy_get_class, "libwayland-client"), Cstring, (Ptr{WlProxy},), proxy)
end
"Set `proxy`'s event queue."
function setqueue!(proxy::Ptr{<: WlProxy}, queue::Ptr{WlEventQueue})
    ccall((:wl_proxy_set_queue, "libwayland-client"), Cvoid, (Ptr{WlProxy},), proxy)
end
# Proxy wrapper
"Create a `WlProxyWrapper` for `proxy`."
function create_wrapper(proxy::Ptr{<: WlProxy})
    wrapper = ccall((:wl_proxy_create_wrapper, "libwayland-client"), Ptr{WlProxyWrapper}, (Ptr{<: WlProxy},), proxy)
    if wrapper == C_NULL
        error("Failed to create a wrapper.")
    else
        return wrapper
    end
end
"Destroy a `WlProxyWrapper`."
function destroy!(wrapper::Ptr{WlProxyWrapper})
    ccall((:wl_proxy_wrapper_destroy, "libwayland-client"), Cvoid, (Ptr{WlProxyWrapper},), wrapper)
end
# Marshalling messages
"Marshal a message."
function marshal(proxy::Ptr{<: WlProxy}, opcode::Integer, args...)
    ccall((:wl_proxy_marshal, "libwayland-client"), Cvoid,
    (Ptr{WlProxy}, UInt32, Any...),
    proxy, opcode, args...)
end
marshal(proxy::Ptr{<: WlProxy}, opcode::Integer, args::Vector) = marshal(proxy, opcode, args...)
"Marshal a message, which may be a constructor."
function marshal_constructor(proxy::Ptr{<: WlProxy}, opcode::Integer, interface::Ref{WlInterface}, args...)
    ccall((:wl_proxy_marshal_constructor, "libwayland-client"), Ptr{WlProxy},
    (Ptr{WlProxy}, UInt32, Ref{WlInterface}, Any...),
    proxy, opcode, interface, args...)
end
marshal_constructor(proxy::Ptr{<: WlProxy}, opcode::Integer, args::Vector, interface::Ref{WlInterface}) = marshal_constructor(proxy, opcode, interface, args...)
"Versioned [`marshal_constructor`](@ref)"
function marshal_constructor_versioned(proxy::Ptr{<: WlProxy}, opcode::Integer, interface::Ref{WlInterface}, version::Integer, args...)
    ccall((:wl_proxy_marshal_constructor_versioned, "libwayland-client"), Ptr{WlProxy},
    (Ptr{WlProxy}, UInt32, Ref{WlInterface}, UInt32, Any...),
    proxy, opcode, interface, version, args...)
end
marshal_constructor_versioned(proxy::Ptr{<: WlProxy}, opcode::Integer, args::Vector, interface::Ref{WlInterface}, version::Integer) = marshal_constructor_versioned(proxy, opcode, interface, args...)
# Listening
"Attach listener struct to `proxy`, overriding the previous."
function addlistener!(proxy::Ptr{WlProxyReal}, implementation::Ref, data::Ref)
    ccall((:wl_proxy_add_listener, "libwayland-client"), Cint,
    (Ptr{WlProxy}, Ref, Ref),
    proxy, implementation, data)
end
function adddispatcher!(proxy::Ptr{WlProxyReal}, dispatcher_func::CFunction, dispatcher_data::Ref, data::Ref)
    ccall((:wl_proxy_add_dispatcher, "libwayland-client"), Cint,
    (Ptr{WlProxy}, CFunction, Ref, Ref),
    proxy, dispatcher_func, dispatcher_data, data)
end
"Get the `proxy`'s listener."
function getlistener(proxy::Ptr{WlProxyReal})
    ccall((:wl_proxy_get_listener, "libwayland-client"), Ptr, (Ptr{WlProxy},), proxy)
end
# User data (the same data which is added to a proxy when adding listeners and is supplied to the listeners)
"Set `proxy`'s user data, overriding the previous."
function setuserdata!(proxy::Ptr{WlProxyReal}, data::Ref)
    ccall((:wl_proxy_set_user_data, "libwayland-client"), Cvoid,
    (Ptr{WlProxy}, Ref),
    proxy, data)
end
"Get `proxy`'s user data."
function getuserdata(proxy::Ptr{WlProxyReal})
    ccall((:wl_proxy_get_user_data, "libwayland-client"), Ref, (Ptr{WlProxy},), proxy)
end

end # module


text:
module BackwardsLinalg
import LinearAlgebra: svd, qr, lq, eigen
using LinearAlgebra, Requires

export svd, qr, lq, symeigen, rsvd

struct ZeroAdder end
Base.:+(a, zero::ZeroAdder) = a
Base.:+(zero::ZeroAdder, a) = a
Base.:-(a, zero::ZeroAdder) = a
Base.:-(zero::ZeroAdder, a) = -a
Base.:-(zero::ZeroAdder) = zero

include("qr.jl")
include("svd.jl")
include("rsvd.jl")
include("symeigen.jl")
include("zygote.jl")

function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("cudalib.jl")
end
end


text:
function distributeCandies(nums::Array{Int64})
    a = Set(nums)
    if length(a) > length(nums)/2
        length(nums) / 2
    else
        length(a)
    end
end

using Test
@test distributeCandies([1,1,2,2,3,3]) == 3
@test distributeCandies([1,1,2,3]) == 2

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : false
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 64 - 1
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 64 - 1
const nz          = haskey(ENV, "NZ"         ) ? parse(Int , ENV["NZ"]         ) : 64 - 1
###
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using Plots, Printf, Statistics, LinearAlgebra, MAT
# CPU functions
@views av_xza(A) = (A[1:end-1,:,1:end-1] .+ A[1:end-1,:,2:end] .+ A[2:end,:,1:end-1] .+ A[2:end,:,2:end]).*0.25
@views av_zi(A)  = (A[2:end-1,2:end-1,2:end-2] .+ A[2:end-1,2:end-1,3:end-1]).*0.5
@views av_za(A)  = (A[:,:,1:end-1] .+ A[:,:,2:end]).*0.5
@views inn(A)    =  A[2:end-1,2:end-1,2:end-1]

@parallel function smooth!(A2::Data.Array, A::Data.Array, fact::Data.Number)
    @inn(A2)   = @inn(A) + 1.0/6.1/fact*(@d2_xi(A) + @d2_yi(A) + @d2_zi(A))
    return
end

@parallel function compute_maxloc!(Musτ2::Data.Array, Musτ::Data.Array)
    @inn(Musτ2) = @maxloc(Musτ)
    return
end

macro Mu_eff() esc(:(1.0/(1.0/@all(Musτ) + 1.0/(G*dt)))) end
@parallel function compute_iter_params!(dτ_Rho::Data.Array, Gdτ::Data.Array, Musτ::Data.Array, Vpdτ::Data.Number, G::Data.Number, dt::Data.Number, Re::Data.Number, r::Data.Number, max_lxyz::Data.Number)
    @all(dτ_Rho) = Vpdτ*max_lxyz/Re/@Mu_eff()
    @all(Gdτ)    = Vpdτ^2/@all(dτ_Rho)/(r+2)
    return
end

@parallel_indices (ix,iy,iz) function assign_τ!(τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array,  τxz::Data.Array,  τyz::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τzz_o::Data.Array, τxy_o::Data.Array, τxz_o::Data.Array, τyz_o::Data.Array)
    if (ix<=size(τxx,1) && iy<=size(τxx,2) && iz<=size(τxx,3))  τxx_o[ix,iy,iz] = τxx[ix,iy,iz]  end
    if (ix<=size(τyy,1) && iy<=size(τyy,2) && iz<=size(τyy,3))  τyy_o[ix,iy,iz] = τyy[ix,iy,iz]  end
    if (ix<=size(τzz,1) && iy<=size(τzz,2) && iz<=size(τzz,3))  τzz_o[ix,iy,iz] = τzz[ix,iy,iz]  end
    if (ix<=size(τxy,1) && iy<=size(τxy,2) && iz<=size(τxy,3))  τxy_o[ix,iy,iz] = τxy[ix,iy,iz]  end
    if (ix<=size(τxz,1) && iy<=size(τxz,2) && iz<=size(τxz,3))  τxz_o[ix,iy,iz] = τxz[ix,iy,iz]  end
    if (ix<=size(τyz,1) && iy<=size(τyz,2) && iz<=size(τyz,3))  τyz_o[ix,iy,iz] = τyz[ix,iy,iz]  end
    return
end

macro inn_yz_Gdτ(ix,iy,iz) esc(:( Gdτ[$ix  ,$iy+1,$iz+1] )) end
macro inn_xz_Gdτ(ix,iy,iz) esc(:( Gdτ[$ix+1,$iy  ,$iz+1] )) end
macro inn_xy_Gdτ(ix,iy,iz) esc(:( Gdτ[$ix+1,$iy+1,$iz  ] )) end
macro inn_yz_Mus(ix,iy,iz) esc(:( Mus[$ix  ,$iy+1,$iz+1] )) end
macro inn_xz_Mus(ix,iy,iz) esc(:( Mus[$ix+1,$iy  ,$iz+1] )) end
macro inn_xy_Mus(ix,iy,iz) esc(:( Mus[$ix+1,$iy+1,$iz  ] )) end
macro av_xyi_Gdτ(ix,iy,iz) esc(:( (Gdτ[$ix  ,$iy  ,$iz+1] + Gdτ[$ix+1,$iy  ,$iz+1] + Gdτ[$ix  ,$iy+1,$iz+1] + Gdτ[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_xzi_Gdτ(ix,iy,iz) esc(:( (Gdτ[$ix  ,$iy+1,$iz  ] + Gdτ[$ix+1,$iy+1,$iz  ] + Gdτ[$ix  ,$iy+1,$iz+1] + Gdτ[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_yzi_Gdτ(ix,iy,iz) esc(:( (Gdτ[$ix+1,$iy  ,$iz  ] + Gdτ[$ix+1,$iy+1,$iz  ] + Gdτ[$ix+1,$iy  ,$iz+1] + Gdτ[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_xyi_Mus(ix,iy,iz) esc(:( (Mus[$ix  ,$iy  ,$iz+1] + Mus[$ix+1,$iy  ,$iz+1] + Mus[$ix  ,$iy+1,$iz+1] + Mus[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_xzi_Mus(ix,iy,iz) esc(:( (Mus[$ix  ,$iy+1,$iz  ] + Mus[$ix+1,$iy+1,$iz  ] + Mus[$ix  ,$iy+1,$iz+1] + Mus[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_yzi_Mus(ix,iy,iz) esc(:( (Mus[$ix+1,$iy  ,$iz  ] + Mus[$ix+1,$iy+1,$iz  ] + Mus[$ix+1,$iy  ,$iz+1] + Mus[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro inn_yz_Gr(ix,iy,iz) esc(:( @inn_yz_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro inn_xz_Gr(ix,iy,iz) esc(:( @inn_xz_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro inn_xy_Gr(ix,iy,iz) esc(:( @inn_xy_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro av_xyi_Gr(ix,iy,iz) esc(:( @av_xyi_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro av_xzi_Gr(ix,iy,iz) esc(:( @av_xzi_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro av_yzi_Gr(ix,iy,iz) esc(:( @av_yzi_Gdτ($ix,$iy,$iz)/(G*dt) )) end
@parallel_indices (ix,iy,iz) function compute_Pt_τ!(Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τzz_o::Data.Array, τxy_o::Data.Array, τxz_o::Data.Array, τyz_o::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Mus::Data.Array, Gdτ::Data.Array, r::Data.Number, G::Data.Number, dt::Data.Number, _dx::Data.Number, _dy::Data.Number, _dz::Data.Number)
    if (ix<=size(Pt,1)  && iy<=size(Pt,2)  && iz<=size(Pt,3))    Pt[ix,iy,iz] = Pt[ix,iy,iz] - r*Gdτ[ix,iy,iz]*(_dx*(Vx[ix+1,iy,iz] - Vx[ix,iy,iz]) + _dy*(Vy[ix,iy+1,iz] - Vy[ix,iy,iz]) + _dz*(Vz[ix,iy,iz+1] - Vz[ix,iy,iz]))  end
    if (ix<=size(τxx,1) && iy<=size(τxx,2) && iz<=size(τxx,3))  τxx[ix,iy,iz] = ( τxx[ix,iy,iz] + τxx_o[ix,iy,iz]*@inn_yz_Gr(ix,iy,iz) + 2.0*@inn_yz_Gdτ(ix,iy,iz)*(_dx*(Vx[ix+1,iy+1,iz+1] - Vx[ix,iy+1,iz+1])))/(1.0 + @inn_yz_Gdτ(ix,iy,iz)/@inn_yz_Mus(ix,iy,iz) + @inn_yz_Gr(ix,iy,iz))  end
    if (ix<=size(τyy,1) && iy<=size(τyy,2) && iz<=size(τyy,3))  τyy[ix,iy,iz] = ( τyy[ix,iy,iz] + τyy_o[ix,iy,iz]*@inn_xz_Gr(ix,iy,iz) + 2.0*@inn_xz_Gdτ(ix,iy,iz)*(_dy*(Vy[ix+1,iy+1,iz+1] - Vy[ix+1,iy,iz+1])))/(1.0 + @inn_xz_Gdτ(ix,iy,iz)/@inn_xz_Mus(ix,iy,iz) + @inn_xz_Gr(ix,iy,iz))  end
    if (ix<=size(τzz,1) && iy<=size(τzz,2) && iz<=size(τzz,3))  τzz[ix,iy,iz] = ( τzz[ix,iy,iz] + τzz_o[ix,iy,iz]*@inn_xy_Gr(ix,iy,iz) + 2.0*@inn_xy_Gdτ(ix,iy,iz)*(_dz*(Vz[ix+1,iy+1,iz+1] - Vz[ix+1,iy+1,iz])))/(1.0 + @inn_xy_Gdτ(ix,iy,iz)/@inn_xy_Mus(ix,iy,iz) + @inn_xy_Gr(ix,iy,iz))  end
    if (ix<=size(τxy,1) && iy<=size(τxy,2) && iz<=size(τxy,3))  τxy[ix,iy,iz] = ( τxy[ix,iy,iz] + τxy_o[ix,iy,iz]*@av_xyi_Gr(ix,iy,iz) + 2.0*@av_xyi_Gdτ(ix,iy,iz)*(0.5*(_dy*(Vx[ix+1,iy+1,iz+1] - Vx[ix+1,iy,iz+1]) + _dx*(Vy[ix+1,iy+1,iz+1] - Vy[ix,iy+1,iz+1])))) / (1.0 + @av_xyi_Gdτ(ix,iy,iz)/@av_xyi_Mus(ix,iy,iz) + @av_xyi_Gr(ix,iy,iz))  end
    if (ix<=size(τxz,1) && iy<=size(τxz,2) && iz<=size(τxz,3))  τxz[ix,iy,iz] = ( τxz[ix,iy,iz] + τxz_o[ix,iy,iz]*@av_xzi_Gr(ix,iy,iz) + 2.0*@av_xzi_Gdτ(ix,iy,iz)*(0.5*(_dz*(Vx[ix+1,iy+1,iz+1] - Vx[ix+1,iy+1,iz]) + _dx*(Vz[ix+1,iy+1,iz+1] - Vz[ix,iy+1,iz+1])))) / (1.0 + @av_xzi_Gdτ(ix,iy,iz)/@av_xzi_Mus(ix,iy,iz) + @av_xzi_Gr(ix,iy,iz))  end
    if (ix<=size(τyz,1) && iy<=size(τyz,2) && iz<=size(τyz,3))  τyz[ix,iy,iz] = ( τyz[ix,iy,iz] + τyz_o[ix,iy,iz]*@av_yzi_Gr(ix,iy,iz) + 2.0*@av_yzi_Gdτ(ix,iy,iz)*(0.5*(_dz*(Vy[ix+1,iy+1,iz+1] - Vy[ix+1,iy+1,iz]) + _dy*(Vz[ix+1,iy+1,iz+1] - Vz[ix+1,iy,iz+1])))) / (1.0 + @av_yzi_Gdτ(ix,iy,iz)/@av_yzi_Mus(ix,iy,iz) + @av_yzi_Gr(ix,iy,iz))  end
    return
end

macro av_xi_dτ_Rho(ix,iy,iz)  esc(:( (dτ_Rho[$ix  ,$iy+1,$iz+1] + dτ_Rho[$ix+1,$iy+1,$iz+1])*0.5 )) end
macro av_yi_dτ_Rho(ix,iy,iz)  esc(:( (dτ_Rho[$ix+1,$iy  ,$iz+1] + dτ_Rho[$ix+1,$iy+1,$iz+1])*0.5 )) end
macro av_zi_dτ_Rho(ix,iy,iz)  esc(:( (dτ_Rho[$ix+1,$iy+1,$iz  ] + dτ_Rho[$ix+1,$iy+1,$iz+1])*0.5 )) end
@parallel_indices (ix,iy,iz) function compute_V!(Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, dτ_Rho::Data.Array, _dx::Data.Number, _dy::Data.Number, _dz::Data.Number, nx_1, nx_2, ny_1, ny_2, nz_1, nz_2)
    if (ix<=nx_1 && iy<=ny_2 && iz<=nz_2)  Vx[ix+1,iy+1,iz+1] = Vx[ix+1,iy+1,iz+1] + (_dx*(τxx[ix+1,iy,iz] - τxx[ix,iy,iz]) + _dy*(τxy[ix,iy+1,iz] - τxy[ix,iy,iz]) + _dz*(τxz[ix,iy,iz+1] - τxz[ix,iy,iz]) - _dx*(Pt[ix+1,iy+1,iz+1] - Pt[ix,iy+1,iz+1]))*@av_xi_dτ_Rho(ix,iy,iz)  end
    if (ix<=nx_2 && iy<=ny_1 && iz<=nz_2)  Vy[ix+1,iy+1,iz+1] = Vy[ix+1,iy+1,iz+1] + (_dy*(τyy[ix,iy+1,iz] - τyy[ix,iy,iz]) + _dx*(τxy[ix+1,iy,iz] - τxy[ix,iy,iz]) + _dz*(τyz[ix,iy,iz+1] - τyz[ix,iy,iz]) - _dy*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy,iz+1]))*@av_yi_dτ_Rho(ix,iy,iz)  end
    if (ix<=nx_2 && iy<=ny_2 && iz<=nz_1)  Vz[ix+1,iy+1,iz+1] = Vz[ix+1,iy+1,iz+1] + (_dz*(τzz[ix,iy,iz+1] - τzz[ix,iy,iz]) + _dx*(τxz[ix+1,iy,iz] - τxz[ix,iy,iz]) + _dy*(τyz[ix,iy+1,iz] - τyz[ix,iy,iz]) - _dz*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy+1,iz]))*@av_zi_dτ_Rho(ix,iy,iz)  end
    return
end

@parallel_indices (ix,iy,iz) function compute_Res!(∇V::Data.Array, Rx::Data.Array, Ry::Data.Array, Rz::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, _dx::Data.Number, _dy::Data.Number, _dz::Data.Number)
    if (ix<=size(∇V,1) && iy<=size(∇V,2) && iz<=size(∇V,3))  ∇V[ix,iy,iz] = _dx*(Vx[ix+1,iy,iz] - Vx[ix,iy,iz]) + _dy*(Vy[ix,iy+1,iz] - Vy[ix,iy,iz]) + _dz*(Vz[ix,iy,iz+1] - Vz[ix,iy,iz]) end
    if (ix<=size(Rx,1) && iy<=size(Rx,2) && iz<=size(Rx,3))  Rx[ix,iy,iz] = _dx*(τxx[ix+1,iy,iz] - τxx[ix,iy,iz]) + _dy*(τxy[ix,iy+1,iz] - τxy[ix,iy,iz]) + _dz*(τxz[ix,iy,iz+1] - τxz[ix,iy,iz]) - _dx*(Pt[ix+1,iy+1,iz+1] - Pt[ix,iy+1,iz+1]) end
    if (ix<=size(Ry,1) && iy<=size(Ry,2) && iz<=size(Ry,3))  Ry[ix,iy,iz] = _dy*(τyy[ix,iy+1,iz] - τyy[ix,iy,iz]) + _dx*(τxy[ix+1,iy,iz] - τxy[ix,iy,iz]) + _dz*(τyz[ix,iy,iz+1] - τyz[ix,iy,iz]) - _dy*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy,iz+1]) end
    if (ix<=size(Rz,1) && iy<=size(Rz,2) && iz<=size(Rz,3))  Rz[ix,iy,iz] = _dz*(τzz[ix,iy,iz+1] - τzz[ix,iy,iz]) + _dx*(τxz[ix+1,iy,iz] - τxz[ix,iy,iz]) + _dy*(τyz[ix,iy+1,iz] - τyz[ix,iy,iz]) - _dz*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy+1,iz]) end
    return
end

@parallel_indices (iy,iz) function bc_x!(A::Data.Array)
    A[  1, iy,  iz] = A[    2,   iy,   iz]
    A[end, iy,  iz] = A[end-1,   iy,   iz]
    return
end

@parallel_indices (ix,iz) function bc_y!(A::Data.Array)
    A[ ix,  1,  iz] = A[   ix,    2,   iz]
    A[ ix,end,  iz] = A[   ix,end-1,   iz]
    return
end

@parallel_indices (ix,iy) function bc_z!(A::Data.Array)
    A[ ix,  iy,  1] = A[   ix,   iy,    2]
    A[ ix,  iy,end] = A[   ix,   iy,end-1]
    return
end

@views function Stokes3D_()
    # Physics
    lx, ly, lz = 10.0, 10.0, 10.0  # domain extends
    ξ          = 1.0               # Maxwell relaxation time
    μs0        = 1.0               # matrix viscosity
    μsi        = 1e-3              # inclusion viscosity
    G          = 1.0               # elastic shear modulus
    εbg        = 1.0               # background strain-rate
    dt         = μs0/(G*ξ)
    # Numerics
    nt         = 1#5                 # number of time steps
    iterMax    = 200#2e5               # maximum number of pseudo-transient iterations
    nout       = 2000               # error checking frequency
    Re         = 6.0π              # Reynolds number
    r          = 1.0               # Bulk to shear elastic modulus ratio
    CFL        = 0.8/sqrt(3)       # CFL number
    ε          = 1e-8              # nonlinear absolute tolerence
    # nx, ny, nz = 127, 127, 127     # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    # Derived numerics
    dx, dy, dz = lx/nx, ly/ny, lz/nz      # cell sizes
    max_lxyz   = max(lx,ly,lz)
    Vpdτ       = min(dx,dy,dz)*CFL
    _dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz
    xc, yc, zc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny), LinRange(dz/2, lz-dz/2, nz)
    # Array allocations
    Pt         = @zeros(nx  ,ny  ,nz  )
    dτ_Rho     = @zeros(nx  ,ny  ,nz  )
    Gdτ        = @zeros(nx  ,ny  ,nz  )
    ∇V         = @zeros(nx  ,ny  ,nz  )
    Vy         = @zeros(nx  ,ny+1,nz  )
    τxx        = @zeros(nx  ,ny-2,nz-2)
    τyy        = @zeros(nx-2,ny  ,nz-2)
    τzz        = @zeros(nx-2,ny-2,nz  )
    τxy        = @zeros(nx-1,ny-1,nz-2)
    τxz        = @zeros(nx-1,ny-2,nz-1)
    τyz        = @zeros(nx-2,ny-1,nz-1)
    τxx_o      = @zeros(nx  ,ny-2,nz-2)
    τyy_o      = @zeros(nx-2,ny  ,nz-2)
    τzz_o      = @zeros(nx-2,ny-2,nz  )
    τxy_o      = @zeros(nx-1,ny-1,nz-2)
    τxz_o      = @zeros(nx-1,ny-2,nz-1)
    τyz_o      = @zeros(nx-2,ny-1,nz-1)
    Rx         = @zeros(nx-1,ny-2,nz-2)
    Ry         = @zeros(nx-2,ny-1,nz-2)
    Rz         = @zeros(nx-2,ny-2,nz-1)
    Mus2       = @zeros(nx  ,ny  ,nz  )
    Musτ       = @zeros(nx  ,ny  ,nz  )
    # Initial conditions
    Rad2       =  zeros(nx  ,ny  ,nz  )
    Vx         =  zeros(nx+1,ny  ,nz  )
    Vz         =  zeros(nx  ,ny  ,nz+1)
    Rad2      .= [(((ix-1)*dx +0.5*dx -0.5*lx)^2 + ((iy-1)*dy +0.5*dy -0.5*ly)^2 + ((iz-1)*dz +0.5*dz -0.5*lz)^2) for ix=1:size(Rad2,1), iy=1:size(Rad2,2), iz=1:size(Rad2,3)]
    Vx         = Data.Array( -εbg.*[((ix-1)*dx -0.5*lx) for ix=1:size(Vx,1), iy=1:size(Vx,2), iz=1:size(Vx,3)] )
    Vz         = Data.Array(  εbg.*[((iz-1)*dz -0.5*lz) for ix=1:size(Vz,1), iy=1:size(Vz,2), iz=1:size(Vz,3)] )
    Mus        = μs0*ones(nx,ny,nz)
    Mus[Rad2.<1.0] .= μsi
    Mus        = Data.Array(Mus)
    Mus2      .= Mus
    for ism=1:10#15
        @parallel smooth!(Mus2, Mus, 1.0)
        Mus, Mus2 = Mus2, Mus
    end
    Musτ      .= Mus
    @parallel compute_maxloc!(Musτ, Mus)
    @parallel (1:size(Musτ,2), 1:size(Musτ,3)) bc_x!(Musτ)
    @parallel (1:size(Musτ,1), 1:size(Musτ,3)) bc_y!(Musτ)
    @parallel (1:size(Musτ,1), 1:size(Musτ,2)) bc_z!(Musτ)
    nx_1, nx_2, ny_1, ny_2, nz_1, nz_2 = nx-1, nx-2, ny-1, ny-2, nz-1, nz-2
    # Preparation of visualisation
    if do_viz || do_save_viz
        ENV["GKSwstype"]="nul"; if do_viz !ispath("../../figures") && mkdir("../../figures") end
        y_sl2, y_sl = Int(ceil((ny-2)/2)), Int(ceil(ny/2))
        xc, zc  = dx/2:dx:(lx-dx/2), dz/2:dz:(lz-dz/2)
        xv, zv  = 0:dx:lx, 0:dz:lz
    end
    # Time loop
    @parallel compute_iter_params!(dτ_Rho, Gdτ, Musτ, Vpdτ, G, dt, Re, r, max_lxyz)
    t=0.0; ittot=0; evo_t=[]; evo_τzz=[]; t_tic = 0.0
    for it = 1:nt
        err=2*ε; iter=0; err_evo1=[]; err_evo2=[]
        @parallel assign_τ!(τxx, τyy, τzz, τxy, τxz, τyz, τxx_o, τyy_o, τzz_o, τxy_o, τxz_o, τyz_o)
        # Pseudo-transient iteration
        while err > ε && iter <= iterMax
            if (it==1 && iter==11) t_tic = Base.time() end
            @parallel compute_Pt_τ!(Pt, τxx, τyy, τzz, τxy, τxz, τyz, τxx_o, τyy_o, τzz_o, τxy_o, τxz_o, τyz_o, Vx, Vy, Vz, Mus, Gdτ, r, G, dt, _dx, _dy, _dz)
            @parallel compute_V!(Vx, Vy, Vz, Pt, τxx, τyy, τzz, τxy, τxz, τyz, dτ_Rho, _dx, _dy, _dz, nx_1, nx_2, ny_1, ny_2, nz_1, nz_2)
            @parallel (1:size(Vy,2), 1:size(Vy,3)) bc_x!(Vy)
            @parallel (1:size(Vz,2), 1:size(Vz,3)) bc_x!(Vz)
            @parallel (1:size(Vx,1), 1:size(Vx,3)) bc_y!(Vx)
            @parallel (1:size(Vz,1), 1:size(Vz,3)) bc_y!(Vz)
            @parallel (1:size(Vx,1), 1:size(Vx,2)) bc_z!(Vx)
            @parallel (1:size(Vy,1), 1:size(Vy,2)) bc_z!(Vy)
            iter += 1
            if iter % nout == 0
                @parallel compute_Res!(∇V, Rx, Ry, Rz, Vx, Vy, Vz, Pt, τxx, τyy, τzz, τxy, τxz, τyz, _dx, _dy, _dz)
                Vmin, Vmax = minimum(Vx), maximum(Vx)
                Pmin, Pmax = minimum(Pt), maximum(Pt)
                norm_Rx    = norm(Rx)/(Pmax-Pmin)*lx/sqrt(length(Rx))
                norm_Ry    = norm(Ry)/(Pmax-Pmin)*lx/sqrt(length(Ry))
                norm_Rz    = norm(Rz)/(Pmax-Pmin)*lx/sqrt(length(Rz))
                norm_∇V    = norm(∇V)/(Vmax-Vmin)*lx/sqrt(length(∇V))
                # norm_Rx = norm(Rx)/length(Rx); norm_Ry = norm(Ry)/length(Ry); norm_Rz = norm(Rz)/length(Rz); norm_∇V = norm(∇V)/length(∇V)
                err = maximum([norm_Rx, norm_Ry, norm_Rz, norm_∇V])
                if isnan(err) error("NaN") end
                push!(err_evo1,maximum([norm_Rx, norm_Ry, norm_Rz, norm_∇V])); push!(err_evo2,iter)
                @printf("Step = %d, iter = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_Rz=%1.3e, norm_∇V=%1.3e] \n", it, iter, err, norm_Rx, norm_Ry, norm_Rz, norm_∇V)
            end
        end
        ittot += iter; t += dt
        # push!(evo_t, t); push!(evo_τzz, maximum(τzz))
    end
    # Performance
    wtime    = Base.time() - t_tic
    A_eff    = (10*2 + 6*1 + 3)/1e9*nx*ny*nz*sizeof(Data.Number) # Effective main memory access per iteration [GB] (Lower bound of required memory access: Te has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    wtime_it = wtime/(ittot-10)                           # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                             # Effective memory throughput [GB/s]
    @printf("Total iters = %d (%d steps), time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", ittot, nt, wtime, round(T_eff, sigdigits=3))
    # Visualisation
    if do_viz
        p1 = heatmap(xc, zc, Array(Pt)[:,y_sl,:]', aspect_ratio=1, xlims=extrema(xc), zlims=extrema(zc), c=:viridis, title="Pressure")
        p2 = heatmap(xc, zv, Array(Vz)[:,y_sl,:]', aspect_ratio=1, xlims=extrema(xc), zlims=extrema(zv), c=:viridis, title="Vz")
        p4 = heatmap(xc[2:end-1], xv[2:end-1], log10.(abs.(Array(Rz)[:,y_sl2,:]')), aspect_ratio=1,  xlims=extrema(xc[2:end-1]), zlims=extrema(zv[2:end-1]), c=:viridis, title="log10(Rz)")
        #p5 = plot(err_evo2,err_evo1, legend=false, xlabel="# iterations", ylabel="log10(error)", linewidth=2, markershape=:circle, markersize=3, labels="max(error)", yaxis=:log10)
        p3 = plot(evo_t, evo_τzz, legend=false, xlabel="time", ylabel="max(τzz)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
            #plot!(evo_t, 2.0.*εbg.*μs0.*(1.0.-exp.(.-evo_t.*G./μs0)), linewidth=2.0) # analytical solution
        plot(p1, p2, p4, p3)
        savefig("../../figures/Stokes_3D_ve_perf_$(nx).png")
    end
    if do_save
        !ispath("../../output") && mkdir("../../output")
        open("../../output/out_Stokes3D_ve_perf.txt","a") do io
            println(io, "1 $(nx) $(ny) $(nz) $(ittot) $(wtime) $(A_eff) $(wtime_it) $(T_eff)")
        end
    end
    if do_save_viz
        !ispath("../../out_visu") && mkdir("../../out_visu")
        matwrite("../../out_visu/Stokes_3D_ve_perf.mat", Dict("Pt_3D"=> Array(Pt), "Mus_3D"=> Array(Mus), "Txz_3D"=> Array(τxz), "Vz_3D"=> Array(Vz), "dx_3D"=> dx, "dy_3D"=> dy, "dz_3D"=> dz); compress = true)
    end
    return xc, yc, zc, Pt
end

if use_return
    xc, yc, zc, P = Stokes3D_();
else
    Stokes3D = begin Stokes3D_(); return; end
end


text:
#jl #! format: off
# # PTDF with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Sourabh Dalvi

# ## Introduction

# PowerSimulations.jl supports linear PTDF optimal power flow formulation. This example shows a
# single multi-period optimization of economic dispatch with a linearized DC-OPF representation of
# using PTDF power flow and how to extract duals values or locational marginal prices for energy.

# ## Dependencies
using SIIPExamples
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder
using DataFrames

# Since we'll be retrieving duals, we need a solver that returns duals values
# here we use Ipopt.
using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

# We can use the same RTS data and some of the initialization as in
# [OperationsProblem example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/01_operations_problems.ipynb)
sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

# Here, we want do define an economic dispatch (linear generation decisions) with
# linear DC-OPF using PTDF network representation.
# So, starting with the network, we can select from _almost_ any of the endpoints on this
# tree:
print_tree(PowerSimulations.PM.AbstractPowerModel)

# For now, let's just choose a standard PTDF formulation.
ed_template = template_economic_dispatch(network = StandardPTDFModel)

# Currently  energy budget data isn't stored in the RTS-GMLC dataset.
set_device_model!(ed_template, HydroEnergyReservoir, HydroDispatchRunOfRiver)

# Calculate the PTDF matrix.
PTDF_matrix = PTDF(sys)

# Now we can build a 4-hour economic dispatch / OPF problem with the RTS data.
# Here, we have to pass the keyword argument `constraint_duals` to OperationsProblem
# with the name of the constraint for which duals are required for them to be returned in the results.
problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 1,
    optimizer = solver,
    balance_slack_variables = true,
    constraint_duals = [
        :CopperPlateBalance,
        :network_flow__Line,
        :network_flow__TapTransformer,
    ],
    PTDF = PTDF_matrix,
)
build!(problem, output_dir = mktempdir())

# And solve the problem and collect the results
solve!(problem)

# Here we collect the dual values from the results for the `:CopperPlateBalance` and `:network_flow`
# constraints. In the case of PTDF network formulation we need to compute the final LMP for each bus in the system by
# subtracting the duals (μ) of `:network_flow` constraints multiplied by the PTDF matrix
# from the  dual (λ) of `:CopperPlateBalance` constraint.
res = ProblemResults(problem)
duals = get_duals(res)
λ = convert(Array, duals[:CopperPlateBalance][:, :CopperPlateBalance])
flow_duals = outerjoin(
    [duals[k] for k in [:network_flow__Line, :network_flow__TapTransformer]]...,
    on = :DateTime,
)
μ = Matrix(flow_duals[:, PTDF_matrix.axes[1]])

# Here we create Dict to store the calculate congestion component of the LMP which is a product of μ and the PTDF matrix.
buses = get_components(Bus, sys)
congestion_lmp = Dict()
for bus in buses
    congestion_lmp[get_name(bus)] = μ * PTDF_matrix[:, get_number(bus)]
end
congestion_lmp = DataFrame(congestion_lmp)

# Finally here we get the LMP for each node in a lossless DC-OPF using the PTDF formulation.
LMP = λ .- congestion_lmp


text:
export SSP

"""
    Atmostphere
All objects describing atmospheric condtions should be a subtype of this.
This includes objects: SSP
"""
abstract type Atmosphere end

# --------------- #
# - SSP - #
# --------------- #

"""
    SSP
Sound speed profile.
### Fields:
* `sound_speed_profile`  - x-depth, y-sound speed.
### Available Constructors:
* `SSP(;sound_speed_profile= Vec2(10.0f0, 10.0f0))`
* `SSP(sound_speed_profile)`
"""
struct SSP{T<:AbstractArray} <: Atmosphere
    sound_speed_profile::Vec2{T}
end

SSP(;sound_speed_profile= Vec2(10.0f0, 10.0f0)) =
    SSP(sound_speed_profile)


text:
using GaussQuadrature
using Base.Test

# write your own tests here
include("test_legendre.jl")


text:
function calculate_measure!(measure::ExpectedReturn{S,EstimatedCase}, w) where {S<:AmbiguitySet}
    return dot(mean(measure.ambiguity_set), w)
end

function calculate_measure!(measure::Variance{S,EstimatedCase}, w::Union{Vector{VariableRef},Real}) where {S<:AmbiguitySet}
    return dot(cov(measure.ambiguity_set) * w, w)
end

function calculate_measure!(measure::Variance{S,EstimatedCase}, w::Vector{AffExpr}) where {S<:AmbiguitySet}
    model = owner_model(w)

    # Cholesky decomposition of the covariance matrix
    Σ = PDMat(Symmetric(cov(measure.ambiguity_set)))
    sqrt_Σ = collect(Σ.chol.U)

    # Extra dimention to represent the portfolio variance
    @variable(model, risk);
    @constraint(model, [risk; 0.5; sqrt_Σ * w] in JuMP.RotatedSecondOrderCone())

    return risk
end

function calculate_measure!(measure::SqrtVariance{S,EstimatedCase}, w) where {S<:AmbiguitySet}
    model = owner_model(w)

    # Cholesky decomposition of the covariance matrix
    Σ = PDMat(Symmetric(cov(measure.ambiguity_set)))
    sqrt_Σ = collect(Σ.chol.U)

    # Extra dimention to represent the square root of the portfolio variance
    @variable(model, risk);
    @constraint(model, [risk; sqrt_Σ * w] in JuMP.SecondOrderCone())

    return risk
end

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PDEAssembles.jl/blob/master/LICENSE

if VERSION < v"1.0.0"
    push!(LOAD_PATH, joinpath(homedir(), ".julia", "dev"))
    using Base.Test
else
    using Test
end

using PDEAssembler
using PDEAssembler: get_unit_square, Poisson

field_elements, boundary_elements = get_unit_square()

field_problem = Problem(Poisson, "Poisson problem in 1x1 square", 1)
update!(field_elements, "source", 10.0)
update!(field_elements, "density", 6.0)
add_elements!(field_problem, field_elements)

boundary_problem = Problem(Poisson, "Poisson boundary", 1)
update!(boundary_elements, "fixed u", 0.0)
update!(boundary_elements, "density", 0.0)
add_elements!(boundary_problem, boundary_elements)

problems = (field_problem, boundary_problem)
M, C, K, f = get_global_matrices(problems, 0.0)
u = cholfact(Symmetric(K)) \ full(f)
X = first(problems)("geometry", 0.0)

N = length(u)
x = [X[i][1] for i=1:N]
y = [X[i][2] for i=1:N]

#using Plots
#surface(x, y, u)

@test isapprox(maximum(u), 0.7381696594268394)


text:
df(sem_fit) = df(sem_fit.model)

df(model::AbstractSemSingle) = n_dp(model) - n_par(model)

df(model::SemEnsemble) = n_dp(model) - n_par(model)

function n_dp(model::AbstractSemSingle)
    nman = n_man(model)
    ndp = 0.5(nman^2 + nman)
    if !isnothing(model.imply.μ)
        ndp += n_man(model)
    end
    return ndp
end

n_dp(model::SemEnsemble) = sum(n_dp.(model.sems))


text:
using MethodsFor
using Test

@testset "MethodsFor.jl" begin
    # Write your tests here.
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AliasPattern <: SwaggerModel
    phrase::Any # spec type: Union{ Nothing, String } # spec name: phrase
    variable::Any # spec type: Union{ Nothing, String } # spec name: variable
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function AliasPattern(;phrase=nothing, variable=nothing, type=nothing)
        o = new()
        validate_property(AliasPattern, Symbol("phrase"), phrase)
        setfield!(o, Symbol("phrase"), phrase)
        validate_property(AliasPattern, Symbol("variable"), variable)
        setfield!(o, Symbol("variable"), variable)
        validate_property(AliasPattern, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type AliasPattern

const _property_map_AliasPattern = Dict{Symbol,Symbol}(Symbol("phrase")=>Symbol("phrase"), Symbol("variable")=>Symbol("variable"), Symbol("type")=>Symbol("type"))
const _property_types_AliasPattern = Dict{Symbol,String}(Symbol("phrase")=>"String", Symbol("variable")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ AliasPattern }) = collect(keys(_property_map_AliasPattern))
Swagger.property_type(::Type{ AliasPattern }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AliasPattern[name]))}
Swagger.field_name(::Type{ AliasPattern }, property_name::Symbol) =  _property_map_AliasPattern[property_name]

const _allowed_AliasPattern_type = ["NotSpecified", "Extract"]

function check_required(o::AliasPattern)
    true
end

function validate_property(::Type{ AliasPattern }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "AliasPattern", :enum, val, _allowed_AliasPattern_type)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ dba34531-0a3d-4598-a4ba-52cf81203caf
using Pkg; Pkg.add(["Plots","LinearAlgebra"])

# ╔═╡ ed657f60-a5f3-11eb-2974-8714af0f16dc
using LinearAlgebra,Plots

# ╔═╡ 97916012-f5e5-4957-99c3-4f91d7b86b9d
function ansatz(r1,r2,α)
    norm_r1 = norm(r1)
    norm_r2 = norm(r2)
    r12 = norm(r1-r2)
    wf = exp(-2*norm_r1)*exp(-2*norm_r2)*exp(r12/(2*(1+α*r12)))
    return wf
end;

# ╔═╡ 9a7c2034-7b47-4c70-b469-60d0e1bcf159
function prob_density(r1,r2,α)
    return ansatz(r1,r2,α)^2
end;

# ╔═╡ 9f80362c-5e59-4e73-b878-ddba6429974a
function E_local(r1,r2,α)
    norm_r1 = norm(r1)
    norm_r2 = norm(r2)
    r12 = norm(r1-r2)
    dot_product = dot(r1/norm_r1-r2/norm_r2,r1-r2)
    energy = -4+dot_product/(r12*(1+α*r12)^2)-1/(r12*(1+α*r12)^3)-1/(4*(1+α*r12)^4)+1/r12
    return energy
end;

# ╔═╡ 3c677d47-d02c-4328-9e12-7ad2ce540ade
function metropolis(N, α)

    L = 1
    r1 = rand(3)*2*L .- L
    r2 = rand(3)*2*L .- L #random number from -L to L
    E = 0
    E2 = 0
    Eln_average = 0
    ln_average = 0
    rejection_ratio = 0
    step = 0
    max_steps = -1

    for i in 1:N
        chose = rand()
        step = step + 1
        if chose < 0.5
            r1_trial = r1 + 0.5*(rand(3)*2*L .- L)
            r2_trial = r2
        else
            r2_trial = r2 + 0.5*(rand(3)*2*L .- L)
            r1_trial = r1
        end
        if prob_density(r1_trial,r2_trial,α) >= prob_density(r1,r2,α)
            r1 = r1_trial
            r2 = r2_trial
        else
            dummy = rand()
            if dummy < prob_density(r1_trial,r2_trial,α)/prob_density(r1,r2,α)
                r1 = r1_trial
                r2 = r2_trial
            else
                rejection_ratio += 1/N
            end
        end

        if step > max_steps
            E += E_local(r1,r2,α)/(N-max_steps)
            E2 += E_local(r1,r2,α)^2/(N-max_steps)
            r12 = norm(r1-r2)
            Eln_average += (E_local(r1,r2,α)*-r12^2/(2*(1+α*r12)^2))/(N-max_steps)
            ln_average += -r12^2/(2*(1+α*r12)^2)/(N-max_steps)
        end
    end
    return E, E2, Eln_average, ln_average, rejection_ratio
end;

# ╔═╡ 4ad41744-94c3-4624-b5b3-f5bf679a7957

function iterate_alpha()
	α = 0

	alpha_iterations = 15
	N_metropolis = 5000
	random_walkers = 200
	γ = 0.5

	energy_plot = []
	alpha_plot = []
	variance_plot = []

	for i in 1:alpha_iterations
	    E = 0
	    E2 = 0
	    dE_dalpha = 0
	    Eln = 0
	    ln = 0
	    rejection_ratio = 0

	    for j in 1:random_walkers
	        E_met, E2_met, Eln_met, ln_met, rejections_met = metropolis(N_metropolis,α)
	        E += E_met/random_walkers
	        E2 += E2_met/random_walkers
	        Eln += Eln_met/random_walkers
	        ln += ln_met/random_walkers
	        rejection_ratio += rejections_met/random_walkers
	    end

		α = α + 0.05

	    energy_plot = append!(energy_plot, E)
	    alpha_plot = append!(alpha_plot, α)
	    variance_plot = append!(variance_plot, E2-E^2)
	end
	return energy_plot,alpha_plot,variance_plot
end;

# ╔═╡ 82e68fc3-f949-43f8-8704-f05e19dec3aa
energy,alpha,variance = iterate_alpha();

# ╔═╡ 6f50f51d-1ec8-4031-af1a-0d6cd69171a0
plot(energy,
	label = false,
	title = "Energy expenctation Plot",
	linecolor = :blue,
	ylabel = "<E>",
	xlabel = "timestep"
)

# ╔═╡ d797e924-293c-496b-8ac8-f06b59b9ad89
plot(variance,
	title = "Variance in energy",
	label = false,
	linecolor = :blue,
	ylabel = "Var(E)",
	xlabel = "timestep"
)

# ╔═╡ 03294208-6397-4602-9157-58ca24154003
plot(alpha,
	title = "Variation Parameter",
	label = false,
	linecolor = :blue,
	ylabel = "α",
	xlabel = "timestep"
)

# ╔═╡ 53b76a30-942b-4191-982a-f06fb38bd52f
scatter(alpha,
	energy,
	title = "Energy vs α",
	xlabel = "α",
	ylabel = "<E>",
	label = false,
	markercolor = :blue
)

# ╔═╡ Cell order:
# ╠═dba34531-0a3d-4598-a4ba-52cf81203caf
# ╠═ed657f60-a5f3-11eb-2974-8714af0f16dc
# ╠═97916012-f5e5-4957-99c3-4f91d7b86b9d
# ╠═9a7c2034-7b47-4c70-b469-60d0e1bcf159
# ╠═9f80362c-5e59-4e73-b878-ddba6429974a
# ╠═3c677d47-d02c-4328-9e12-7ad2ce540ade
# ╠═4ad41744-94c3-4624-b5b3-f5bf679a7957
# ╠═82e68fc3-f949-43f8-8704-f05e19dec3aa
# ╠═6f50f51d-1ec8-4031-af1a-0d6cd69171a0
# ╠═d797e924-293c-496b-8ac8-f06b59b9ad89
# ╠═03294208-6397-4602-9157-58ca24154003
# ╠═53b76a30-942b-4191-982a-f06fb38bd52f

max_stars_repo_name: Ved-Mahajan/VMC

max_issues_repo_name: Ved-Mahajan/VMC

max_forks_repo_name: Ved-Mahajan/VMC


text:
module GaussianProcesses

export k, K, GP_Regressor

using LinearAlgebra, Distributions

# exponential quadratic kernel
k(a::Array, b::Array; h::Real, d::Function) = exp(-d(a,b)^2 / (2 * h^2))
k(a::Real, b::Real; h::Real, d::Function) = k([a], [b]; h=h, d=d)

# kernel matrix
K(A::Vector, B::Vector; h::Real, d::Function) = [k(a, b; h = h, d = d) for a in A, b in B]

# prior distribution: mean 0, covariance K
mutable struct GP_Regressor
	"MVN distribution of prior"
	prior::Distributions.MvNormal
	"MVN distribution of posterior"
	posterior::Distributions.MvNormal
	"equivalent smoothing kernel weights (rows of matrix)"
	equiv_weights::Array
	"input points for sampling and plotting"
	x_grid::Vector
	"sample observations"
	x_train::Vector
	"sample outcomes"
	y_train::Vector
	"standard deviation of iid Gaussian noise in observations"
	σ::Real
	"characteristic length scale"
	h::Real
	"distance function"
	d::Function

	# inner constructor: call this function to instantiate
	function GP_Regressor(X_train::Vector, y_train::Vector, X::Vector; σ::Real, h::Real, d::Function)

		# prior
		Kpred = K(X, X; h=h, d=d)  # covariance
		prior = MvNormal(Kpred)

		# posterior terms
		Kinv = inv(K(X_train, X_train; h=h, d=d) + σ^2 * I)
		Kcross = K(X, X_train; h=h, d=d)
		Kequiv = Kcross * Kinv  # equivalent smoothing kernel weights (each row)

		mean_posterior = Kequiv * y_train
		cov_posterior = Kpred .- Kequiv * Kcross'

		posterior = MvNormal(mean_posterior, Matrix(Hermitian(cov_posterior)))

		new(prior, posterior, Kequiv, X, X_train, y_train, σ, h, d)
	end
end

end  # module

max_stars_repo_name: tobanw/gaussian-process-regression

max_issues_repo_name: tobanw/gaussian-process-regression

max_forks_repo_name: tobanw/gaussian-process-regression


text:
abstract type AbstractNode end

struct Component{N} <: AbstractNode
    type::Symbol
    name::Symbol
    freeports::MVector{N,Bool}
    vertex::RefValue{Int}
    parameters::Vector{Num}
    state_vars::Vector{Num}
    default::Dict{Num,Any}
    equations::Vector{Equation}
    function Component{N}(m::Symbol, n::Symbol, np::Int, v::Int,
         p::Vector, x::Vector, d::Dict, eq::Vector) where N
        new(m, n, ones(MVector{np,Bool}), Ref(v), p, x, d, eq)
    end
end
function Component(type::Symbol, name::Symbol=type;
    numports::Int=1, vertex::Int=0, parameters::Vector=[], state_vars::Vector=[],
    default::Dict=Dict(), equations::Vector=[])
    Component{numports}(type, name, numports, vertex, parameters, state_vars, default, equations)
end

abstract type Junction <: AbstractNode end

struct EqualEffort <: Junction
    name::Symbol
    freeports::Vector{Bool}
    vertex::RefValue{Int}
    EqualEffort(; name::Symbol=Symbol("0"), v::Int=0) = new(name, [], Ref(v))
end

struct EqualFlow <: Junction
    name::Symbol
    freeports::Vector{Bool}
    weights::Vector{Int}
    vertex::RefValue{Int}
    EqualFlow(; name::Symbol=Symbol("1"), v::Int=0) = new(name, [], [], Ref(v))
end

struct Port
    node::AbstractNode
    index::Int
    function Port(node::AbstractNode, index)
        ports = freeports(node)
        any(ports) || error("Node $node has no free ports")
        ports[index] || error("Port $index in node $node is already connected")
        new(node, index)
    end
    Port(node::Junction, index) = new(node, index)
end
Port(node::AbstractNode) = Port(node, nextfreeport(node))

struct Bond <: lg.AbstractSimpleEdge{Int}
    src::Port
    dst::Port
end
function Bond(srcnode::AbstractNode, dstnode::AbstractNode)
    Bond(Port(srcnode), Port(dstnode))
end

struct BondGraph <: lg.AbstractGraph{Int64}
    type::Symbol
    name::Symbol
    nodes::Vector{T} where T <: AbstractNode
    bonds::Vector{Bond}
end
BondGraph(name::Symbol=:BG) = BondGraph(:BG, name, AbstractNode[], Bond[])

# New component
function new(type,name::Symbol=type;library=standard_library)
    d = library[type]
    p = collect(keys(d[:parameters]))
    if haskey(d,:state_vars)
        x = collect(keys(d[:state_vars]))
    else
        x = Vector{Num}([])
    end
    Component(type, name; numports=d[:numports], parameters=p, state_vars=x, equations=d[:equations])
end

# All models
Model = Union{AbstractNode,BondGraph}

# Component type
type(n) = n.type
type(n::EqualEffort) = Symbol("0")
type(n::EqualFlow) = Symbol("1")

# Vertex
vertex(n::AbstractNode) = n.vertex[]
set_vertex!(n::AbstractNode, v::Int) = n.vertex[] = v

# Ports
freeports(n::AbstractNode) = n.freeports
numports(n::AbstractNode) = length(n.freeports)
updateport!(n::AbstractNode, idx::Int) = freeports(n)[idx] = !freeports(n)[idx]

nextfreeport(n::AbstractNode) = findfirst(freeports(n))
function nextfreeport(j::Junction)
    push!(j.freeports,true)
    numports(j)
end
function nextfreeport(j::EqualFlow)
    push!(j.weights,0)
    push!(j.freeports,true)
    numports(j)
end

nextsrcport(n::AbstractNode) = nextfreeport(n)
function nextsrcport(n::EqualFlow)
    i = nextfreeport(n)
    n.weights[i] = -1
    i
end

nextdstport(n::AbstractNode) = nextfreeport(n)
function nextdstport(n::EqualFlow)
    i = nextfreeport(n)
    n.weights[i] = 1
    i
end

# Nodes in Bonds
srcnode(b::Bond) = b.src.node
dstnode(b::Bond) = b.dst.node
in(n::AbstractNode, b::Bond) = n == srcnode(b) || n == dstnode(b)
Base.iterate(b::Bond) = (b.src,true)
function Base.iterate(b::Bond,state)
    if state
        return (b.dst,false)
    else
        return nothing
    end
end

# Parameters
params(n::AbstractNode) = n.parameters
params(j::Junction) = Vector{Num}([])

# Components
components(bg::BondGraph) = bg.nodes

# State variables
state_vars(n::Component) = n.state_vars
state_vars(j::Junction) = Vector{Num}([])
function state_vars(m::BondGraph)
    i = 0
    states = Vector{Tuple{Model,Num}}([])
    for c in components(m), x in state_vars(c)
        push!(states,(c,x))
    end
    @variables x[1:length(states)](t)
    return OrderedDict(x[i] => y for (i,y) in enumerate(states))
end

# Set parameter values
function set_param!(n::Component,var,val)
    any(isequal.(params(n),var)) || error("Component does not have parameter.")
    default_value(n)[var] = val
end

# Set initial conditions
function set_initial_value!(n::Component,var,val)
    any(isequal.(state_vars(n),var)) || error("Component does not have state variable.")
    default_value(n)[var] = val
end

# Get default values
default_value(n::Component) = n.default
default_value(n::Component,v::Num) = default_value(n::Component)[v]
default_value(j::Junction) = Dict{Num,Any}()

# Bonds
bonds(m::BondGraph) = m.bonds

# Searching
getnodes(bg::BondGraph, m::Symbol) = filter(x -> type(x) == m, bg.nodes)
getnodes(bg::BondGraph, n::AbstractString) = filter(x -> x.name == Symbol(n), bg.nodes)

getbonds(bg::BondGraph, t::Tuple) = getbonds(bg, t[1], t[2])
getbonds(bg::BondGraph, n1::AbstractNode, n2::AbstractNode) = filter(b -> n1 in b && n2 in b, bg.bonds)

# I/O
show(io::IO, node::Component) = print(io, "$(node.type):$(node.name)")
show(io::IO, node::Junction) = print(io, "$(node.name)")
show(io::IO, port::Port) = print(io, "Port $(port.node) ($(port.index))")
show(io::IO, b::Bond) = print(io, "Bond $(srcnode(b)) ⇀ $(dstnode(b))")
show(io::IO, bg::BondGraph) = print(io, "BondGraph $(bg.type):$(bg.name) ($(lg.nv(bg)) Nodes, $(lg.ne(bg)) Bonds)")

# Comparisons
# These definitions will need to expand when equations etc. are added
==(n1::AbstractNode, n2::AbstractNode) = type(n1) == type(n2) && n1.name == n2.name

struct BondGraphNode <: AbstractNode
    bondgraph::BondGraph
    type::Symbol
    name::AbstractString
    freeports::Vector{Bool}
    vertex::RefValue{Int}
    function BondGraphNode(bg::BondGraph, type::Symbol=:BG, name::AbstractString=bg.name; vertex::Int=0)
        new(bg, type, name, Bool[], Ref(vertex))
    end
end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
# Forest-Ruth and Position Extended Forest-Ruth Like algorithm
# Leapfrog is 2nd order, FR and PEFRL are 4th order.

function fr_step!(r::AbstractVecOrMat{V3{T}}, v::AbstractVecOrMat{V3{T}},
         m::AbstractVector; G, dt) where T
    nplanets = length(m)
    a = zeros(V3{T}, nplanets)
    θ = 1/(2-2^(1/3))
    @inbounds begin
        for j=1:nplanets
            r[j] += θ/2*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += θ*dt*a[j]
            r[j] += (1-θ)/2*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += (1-2θ)*dt*a[j]
            r[j] += (1-θ)/2*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += θ*dt*a[j]
            r[j] += θ/2*dt*v[j]
        end
    end
    return r, v
end

function pefrl_step!(r::AbstractVecOrMat{V3{T}}, v::AbstractVecOrMat{V3{T}},
         m::AbstractVector; G, dt) where T
    nplanets = length(m)
    a = zeros(V3{T}, nplanets)
    ξ = +0.1786178958448091E+00
    λ = -0.2123418310626054E+00
    χ = -0.6626458266981849E-01
    @inbounds begin
        for j=1:nplanets
            r[j] = r[j]
            v[j] = v[j]
            r[j] += ξ*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += (1-2λ)/2*dt*a[j]
            r[j] += χ*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += λ*dt*a[j]
            r[j] += (1-2(ξ+χ))*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += λ*dt*a[j]
            r[j] += χ*dt*v[j]
        end
        update_acceleration!(a, r, m, G)
        for j=1:nplanets
            v[j] += (1-2λ)/2*dt*a[j]
            r[j] += ξ*dt*v[j]
        end
    end
    return r, v
end

function fr!(r::AbstractVecOrMat{V3{T}}, v::AbstractVecOrMat{V3{T}},
         m::AbstractVector; G=G_year_AU, n, dt) where T
    leapfrog_looper!(FR,r, v, m; G=G, n=n, dt=dt)
end

function pefrl!(r::AbstractVecOrMat{V3{T}}, v::AbstractVecOrMat{V3{T}},
         m::AbstractVector; G=G_year_AU, n, dt) where T
    leapfrog_looper!(PEFRL,r, v, m; G=G, n=n, dt=dt)
end

function fr(planets::AbstractVector{Body{T}}; G=G_year_AU, n, dt, keep_history=false) where {T}
    leapfrog(FR, planets; G=G, n=n, dt=dt, keep_history=keep_history)
end

function pefrl(planets::AbstractVector{Body{T}}; G=G_year_AU, n, dt, keep_history=false) where {T}
    leapfrog(PEFRL, planets; G=G, n=n, dt=dt, keep_history=keep_history)
end

steper(::Type{FR}) = fr_step!
steper(::Type{PEFRL}) = pefrl_step!


text:
"""Count the occurences of the 1-bit in bit position b across all elements of A."""
function bitcount(A::AbstractArray{T},b::Int) where {T<:Union{Integer,AbstractFloat}}
    N = sizeof(T)*8             # number of bits in T
    @boundscheck b <= N || throw(error("Can't count bit $b for $N-bit type $T."))
    UIntT = whichUInt(T)
    n = 0                       # counter
    shift = N-b                 # shift desired bit b
    mask = one(UIntT) << shift  # leftshift always pushes 0s
    for a in A                  # mask everything but b and shift
        n += (reinterpret(UIntT,a) & mask) >>> shift   # to have either 0x00 or 0x01
    end
    return n
end

"""Count the occurences of the 1-bit in every bit position b across all elements of A."""
function bitcount(A::AbstractArray{T}) where {T<:Union{Integer,AbstractFloat}}
    N = 8*sizeof(T)             # determine the size [bit] of elements in A
    n = fill(0,N)               # preallocate counters
    for b in 1:N                # loop over bit positions and count each
        n[b] = bitcount(A,b)
    end
    return n
end

"""Entropy [bit] for bitcount functions. Maximised to 1bit for random uniformly
distributed bits in A."""
function bitcount_entropy(A::AbstractArray)
    N = prod(size(A))
    ps = bitcount(A) / N
    e = [abs(entropy([p,1-p],2)) for p in ps]
    return e
end

"""Count pairs of bits across elements of A for bit position b therein.
Returns a 4-element array with counts for 00,01,10,11."""
function bitpaircount(A::AbstractArray{T},b::Int) where {T<:Union{Integer,AbstractFloat}}
    N = sizeof(T)*8             # number of bits in T
    @boundscheck b <= N || throw(error("Can't count bit $b for $N-bit type $T."))
    UIntT = whichUInt(T)

    n = [0,0,0,0]                   # counter for 00,01,10,11
    shift = N-b-1                   # shift to the 2nd last position either 0x00,0x10
    shift2 = N-b                    # shift to last position 0x0 or 0x1
    mask = one(UIntT) << shift2     # mask everything except b

    # a1 is bit from previous entry in A, a2 is the current
    # a1 is shifted to sit in the 2nd last position
    # a2 sits in the last position
    a1 = (reinterpret(UIntT,A[1]) & mask) >>> shift
    @inbounds for i in 2:length(A)
        a2 = (reinterpret(UIntT,A[i]) & mask) >>> shift2
        n[(a1 | a2)+1] += 1
        a1 = a2 << 1
    end
    return n
end

"""Count pairs of bits across elements of A for every bit position therein.
Returns a 4xn-array with counts for 00,01,10,11 in rows and every bit position
in columns."""
function bitpaircount(A::AbstractArray{T}) where {T<:Union{Integer,AbstractFloat}}
    N = 8*sizeof(T)         # number of bits in T
    n = fill(0,4,N)         # store the 4 possible pair counts for every bit position
    for b in 1:N
        n[:,b] = bitpaircount(A,b)
    end
    return n
end

"""Calculates the conditional probabilities of pairs 00,01,10,11 in the bit sequences
of array A. Returns a 4xn array with conditional probabilities p(nextbit=0|previousbit=0),
p(1|0),p(0|1),p(1|1) in rows and every bit position in columns. Returns NaN when
all bits are either 0,1 (in which case the conditional probability is not defined)."""
function bit_condprobability(A::Array{T}) where {T<:Union{Integer,AbstractFloat}}
    N = prod(size(A[2:end]))        # elements in array (A[1] is )
    n1 = bitcount(A[1:end-1])
    n0 = N.-n1
    npair = bitpaircount(A)
    pcond = similar(npair,Float64)
    pcond[1,:] = npair[1,:] ./ n0
    pcond[2,:] = npair[2,:] ./ n0
    pcond[3,:] = npair[3,:] ./ n1
    pcond[4,:] = npair[4,:] ./ n1
    return pcond
end

"""Calculates the conditional entropy for 00,01,10,11 for every bit position across
all elements of A."""
function bit_condentropy(   A::Array{T},
                            set_zero_insignificant::Bool=true,
                            confidence::Real=0.99) where {T<:Union{Integer,AbstractFloat}}

    pcond = bit_condprobability(A)

    # set NaN (occurs when occurences n=0) 0*-Inf = 0 here.
    pcond[isnan.(pcond)] .= 0

    # conditional entropy H given bit = 0 (H0) and bit = 1 (H1)
    H0 = [entropy([p00,p01],2) for (p00,p01) in zip(pcond[1,:],pcond[2,:])]
    H1 = [entropy([p10,p11],2) for (p10,p11) in zip(pcond[3,:],pcond[4,:])]

    # set insignficant to zero
    if set_zero_insignificant
        # get chance p for 1 (or 0) from binom distr
        p = binom_confidence(length(A),confidence)
        I₀ = 1 - entropy([p,1-p],2)         # free entropy of random [bit]
        H0[H0 .<= I₀] .= 0                  # set insignficant to zero
        H1[H1 .<= I₀] .= 0
    end

    return H0,H1
end

function bitinformation(A::AbstractArray,dims::Symbol)
    if dims == :all_dimensions
        N = prod(size(A))-1                     # elements in array
        n1 = bitcount(A)-bitcount([A[end]])     # occurences of bit = 1
        bi = fill(0.0,length(n1))
        permu = collect(1:ndims(A))                 # permutation
        perm0 = vcat(permu[2:end],permu[1])         # used to permute permu
        for idim in 1:ndims(A)
            Aperm = PermutedDimsArray(A,permu)
            permute!(permu,perm0)
            npair = bitpaircount(Aperm)
            bi .+= bitinformation(n1,npair,N)
        end
        # arithmetic mean of information across the dimensions
        return bi ./ ndims(A)
    else
        return bitinformation(A)
    end
end

function bitinformation(A::AbstractArray;
                        dims::Int=1,
                        set_zero_insignificant::Bool=true,
                        confidence::Real=0.99)

    N = prod(size(A))-1                         # elements in array
    n1 = bitcount(A)-bitcount([A[end]])     # occurences of bit = 1

    if dims > 1
        permu = collect(1:ndims(A))                 # permutation
        perm0 = vcat(permu[2:end],permu[1])         # used to permute permu
        for _ in 1:dims-1
            permute!(permu,perm0)
        end
        A = PermutedDimsArray(A,permu)
    end

    npair = bitpaircount(A)
    return bitinformation(n1,npair,N;set_zero_insignificant,confidence)
end

function bitinformation(n1::Array{Int,1},
                        npair::Array{Int,2},
                        N::Int;
                        set_zero_insignificant::Bool=true,
                        confidence::Real=0.99)

    nbits = length(n1)
    @assert size(npair) == (4,nbits)    "npair array must be of size 4xnbits"
    @assert maximum(n1) <= N            "N cannot be smaller than maximum(n1)"

    n0 = N.-n1                      # occurences of bit = 0
    q0 = n0/N                       # respective probabilities
    q1 = n1/N

    pcond = similar(npair,Float64)      # preallocate conditional probability

    pcond[1,:] = npair[1,:] ./ n0       # p(0|0) = n(00)/n(0)
    pcond[2,:] = npair[2,:] ./ n0       # p(1|0) = n(01)/n(0)
    pcond[3,:] = npair[3,:] ./ n1       # p(0|1) = n(10)/n(1)
    pcond[4,:] = npair[4,:] ./ n1       # p(1|1) = n(11)/n(1)

    # set NaN (occurs when occurences n=0) 0*-Inf = 0 here.
    pcond[isnan.(pcond)] .= 0

    # unconditional entropy
    H = [entropy([q0i,q1i],2) for (q0i,q1i) in zip(q0,q1)]

    # conditional entropy given bit = 0, bit = 1
    H0 = [entropy([p00,p01],2) for (p00,p01) in zip(pcond[1,:],pcond[2,:])]
    H1 = [entropy([p10,p11],2) for (p10,p11) in zip(pcond[3,:],pcond[4,:])]

    # Information content
    I = @. H - q0*H0 - q1*H1

    # remove information that is insignificantly different from a random 50/50
    if set_zero_insignificant
        p = binom_confidence(N,confidence)  # get chance p for 1 (or 0) from binom distr
        I₀ = 1 - entropy([p,1-p],2)         # free entropy of random [bit]
        I[I .<= I₀] .= 0                    # set insignficant to zero
    end

    return I
end