text:
__precompile__(true)
"""
Support for Character Sets, Encodings, and Character Set Encodings

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module CharSetEncodings

using ModuleInterfaceTools

@api extend! StrAPI

# Define symbols used for characters, codesets, codepoints

const cse_info =
    ((:Binary,   UInt8),		# really, no character set at all, not text
     (:Text1,    UInt8),        	# Unknown character set, 1 byte
     (:ASCII,    UInt8),		# (7-bit subset of Unicode)
     (:Latin,    UInt8),		# ISO-8859-1 (8-bit subset of Unicode)
     (:_Latin,   UInt8),		# Latin subset of Unicode (0-0xff)
     (:UTF8,     UInt8,  :UTF32),	# Validated UTF-8
     (:RawUTF8,  UInt8,  :UniPlus),	# Unvalidated UTF-8
     (:Text2,    UInt16),		# Unknown character set, 2 byte
     (:UCS2,     UInt16),		# BMP (16-bit subset of Unicode)
     (:_UCS2,    UInt16),		# UCS-2 Subset of Unicode (0-0xd7ff, 0xe000-0xffff)
     (:UTF16,    UInt16, :UTF32),	# Validated UTF-16
     (:RawUTF16, UInt16, :UniPlus),	# Unvalidated UTF-16
     (:Text4,    UInt32),		# Unknown character set, 4 byte
     (:UTF32,    UInt32),		# corresponding to codepoints (0-0xd7ff, 0xe000-0x10fff)
     (:_UTF32,   UInt32))		# Full validated UTF-32

@api develop cse_info

include("charsets.jl")
include("encodings.jl")
include("cses.jl")
include("traits.jl")

@api freeze

end # module CharSetEncodings


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    fit_mle(gev, bm)

Fit generalized extreme value distribution `gev` to block maxima
`bm` with constrained maximum likelihood estimation.
"""
function fit_mle(::Type{GeneralizedExtremeValue}, bm::BlockMaxima)
  # retrieve maxima values
  x = collect(bm)
  n = length(x)

  # MLE for case ξ ≠ 0
  mle = Model(with_optimizer(Ipopt.Optimizer, print_level=0, sb="yes"))
  @variable(mle, μ, start=0.0)
  @variable(mle, σ, start=1.0)
  @variable(mle, ξ, start=0.1)
  @NLobjective(mle, Max,
    -n*log(σ)
    -(1 + 1/ξ)*sum(log(1 + ξ*(x[i]-μ)/σ) for i in 1:n)
    -sum((1 + ξ*(x[i]-μ)/σ)^(-1/ξ) for i in 1:n)
  )
  @NLconstraint(mle, [i=1:n], 1 + ξ*(x[i]-μ)/σ ≥ 1e-6)
  @NLconstraint(mle, abs(ξ) ≥ 1e-6)
  @constraint(mle, σ ≥ 1e-6)

  # MLE for case ξ = 0
  mle₀ = Model(with_optimizer(Ipopt.Optimizer, print_level=0, sb="yes"))
  @variable(mle₀, μ₀, start=0.0)
  @variable(mle₀, σ₀, start=1.0)
  @NLobjective(mle₀, Max,
    -n*log(σ₀)
    -sum((x[i]-μ₀)/σ₀ for i in 1:n)
    -sum(exp(-(x[i]-μ₀)/σ₀) for i in 1:n)
  )
  @constraint(mle₀, σ₀ ≥ 1e-6)

  # attempt to solve both cases
  optimize!(mle)
  optimize!(mle₀)

  # retrieve solver status
  status  = termination_status(mle)
  status₀ = termination_status(mle₀)

  # acceptable statuses
  OK = (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)

  if status ∈ OK && status₀ ∈ OK
    # choose the maximum amongst the two
    if objective_value(mle) > objective_value(mle₀)
      GeneralizedExtremeValue(value(μ), value(σ), value(ξ))
    else
      GeneralizedExtremeValue(value(μ₀), value(σ₀), 0.)
    end
  elseif status ∈ OK
    GeneralizedExtremeValue(value(μ), value(σ), value(ξ))
  elseif status₀ ∈ OK
    GeneralizedExtremeValue(value(μ₀), value(σ₀), 0.)
  else
    error("could not fit distribution to maxima")
  end
end

"""
    fit_mle(gp, pm)

Fit generalized Pareto distribution `gp` to peak over threshold
maxima `pm` with constrained maximum likelihood estimation.
"""
function fit_mle(::Type{GeneralizedPareto}, pm::PeakOverThreshold)
  # retrieve maxima values
  x = collect(pm)
  n = length(x)

  # compute excess
  y = x .- pm.u

  # MLE for case ξ ≠ 0
  mle = Model(with_optimizer(Ipopt.Optimizer, print_level=0, sb="yes"))
  @variable(mle, σ, start=1.0)
  @variable(mle, ξ, start=0.1)
  @NLobjective(mle, Max,
    -n*log(σ)
    -(1 + 1/ξ)*sum(log(1 + ξ*y[i]/σ) for i in 1:n)
  )
  @NLconstraint(mle, [i=1:n], 1 + ξ*y[i]/σ ≥ 1e-6)
  @NLconstraint(mle, abs(ξ) ≥ 1e-6)
  @constraint(mle, σ ≥ 1e-6)

  # MLE for case ξ = 0
  mle₀ = Model(with_optimizer(Ipopt.Optimizer, print_level=0, sb="yes"))
  @variable(mle₀, σ₀, start=1.0)
  @NLobjective(mle₀, Max,
    -n*log(σ₀)
    -sum(y[i] for i in 1:n)/σ₀
  )
  @constraint(mle₀, σ₀ ≥ 1e-6)

  # attempt to solve both cases
  optimize!(mle)
  optimize!(mle₀)

  # retrieve solver status
  status  = termination_status(mle)
  status₀ = termination_status(mle₀)

  # acceptable statuses
  OK = (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)

  if status ∈ OK && status₀ ∈ OK
    # choose the maximum amongst the two
    if objective_value(mle) > objective_value(mle₀)
      GeneralizedPareto(0.0, value(σ), value(ξ))
    else
      GeneralizedPareto(0.0, value(σ₀), 0.0)
    end
  elseif status ∈ OK
    GeneralizedPareto(0.0, value(σ), value(ξ))
  elseif status₀ ∈ OK
    GeneralizedPareto(0.0, value(σ₀), 0.0)
  else
    error("could not fit distribution to maxima")
  end
end

ext: jl


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

"""
    create_application(author, description, name)
    create_application(author, description, name, params::Dict{String,<:Any})

Creates an application, optionally including an AWS SAM file to create the first
application version in the same call.

# Arguments
- `author`: The name of the author publishing the app.Minimum length=1. Maximum
  length=127.Pattern \"^[a-z0-9](([a-z0-9]|-(?!-))*[a-z0-9])?\";
- `description`: The description of the application.Minimum length=1. Maximum length=256
- `name`: The name of the application that you want to publish.Minimum length=1. Maximum
  length=140Pattern: \"[a-zA-Z0-9-]+\";

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"homePageUrl"`: A URL with more information about the application, for example the
  location of your GitHub repository for the application.
- `"labels"`: Labels to improve discovery of apps in search results.Minimum length=1.
  Maximum length=127. Maximum number of labels: 10Pattern: \"^[a-zA-Z0-9+-_:/@]+\";
- `"licenseBody"`: A local text file that contains the license of the app that matches the
  spdxLicenseID value of your application.
 The file has the format
  file://<path>/<filename>.Maximum size 5 MBYou can specify only one of licenseBody and
  licenseUrl; otherwise, an error results.
- `"licenseUrl"`: A link to the S3 object that contains the license of the app that matches
  the spdxLicenseID value of your application.Maximum size 5 MBYou can specify only one of
  licenseBody and licenseUrl; otherwise, an error results.
- `"readmeBody"`: A local text readme file in Markdown language that contains a more
  detailed description of the application and how it works.
 The file has the format
  file://<path>/<filename>.Maximum size 5 MBYou can specify only one of readmeBody and
  readmeUrl; otherwise, an error results.
- `"readmeUrl"`: A link to the S3 object in Markdown language that contains a more detailed
  description of the application and how it works.Maximum size 5 MBYou can specify only one
  of readmeBody and readmeUrl; otherwise, an error results.
- `"semanticVersion"`: The semantic version of the application:
 https://semver.org/
- `"sourceCodeArchiveUrl"`: A link to the S3 object that contains the ZIP archive of the
  source code for this version of your application.Maximum size 50 MB
- `"sourceCodeUrl"`: A link to a public repository for the source code of your application,
  for example the URL of a specific GitHub commit.
- `"spdxLicenseId"`: A valid identifier from https://spdx.org/licenses/.
- `"templateBody"`: The local raw packaged AWS SAM template file of your application.
 The
  file has the format file://<path>/<filename>.You can specify only one of templateBody
  and templateUrl; otherwise an error results.
- `"templateUrl"`: A link to the S3 object containing the packaged AWS SAM template of your
  application.You can specify only one of templateBody and templateUrl; otherwise an error
  results.
"""
function create_application(
    author, description, name; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "POST",
        "/applications",
        Dict{String,Any}("author" => author, "description" => description, "name" => name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_application(
    author,
    description,
    name,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "POST",
        "/applications",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "author" => author, "description" => description, "name" => name
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_application_version(application_id, semantic_version)
    create_application_version(application_id, semantic_version, params::Dict{String,<:Any})

Creates an application version.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.
- `semantic_version`: The semantic version of the new version.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"sourceCodeArchiveUrl"`: A link to the S3 object that contains the ZIP archive of the
  source code for this version of your application.Maximum size 50 MB
- `"sourceCodeUrl"`: A link to a public repository for the source code of your application,
  for example the URL of a specific GitHub commit.
- `"templateBody"`: The raw packaged AWS SAM template of your application.
- `"templateUrl"`: A link to the packaged AWS SAM template of your application.
"""
function create_application_version(
    applicationId, semanticVersion; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "PUT",
        "/applications/$(applicationId)/versions/$(semanticVersion)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_application_version(
    applicationId,
    semanticVersion,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "PUT",
        "/applications/$(applicationId)/versions/$(semanticVersion)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_cloud_formation_change_set(application_id, stack_name)
    create_cloud_formation_change_set(application_id, stack_name, params::Dict{String,<:Any})

Creates an AWS CloudFormation change set for the given application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.
- `stack_name`: This property corresponds to the parameter of the same name for the AWS
  CloudFormation CreateChangeSet
  API.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"capabilities"`: A list of values that you must specify before you can deploy certain
  applications.
 Some applications might include resources that can affect permissions in
  your AWS
 account, for example, by creating new AWS Identity and Access Management (IAM)
  users.
 For those applications, you must explicitly acknowledge their capabilities by
  specifying this parameter.The only valid values are CAPABILITY_IAM, CAPABILITY_NAMED_IAM,
  CAPABILITY_RESOURCE_POLICY, and CAPABILITY_AUTO_EXPAND.The following resources require you
  to specify CAPABILITY_IAM or
 CAPABILITY_NAMED_IAM:
 AWS::IAM::Group,
  AWS::IAM::InstanceProfile,
 AWS::IAM::Policy, and
 AWS::IAM::Role.
 If the application
  contains IAM resources, you can specify either CAPABILITY_IAM
 or CAPABILITY_NAMED_IAM. If
  the application contains IAM resources
 with custom names, you must specify
  CAPABILITY_NAMED_IAM.The following resources require you to specify
  CAPABILITY_RESOURCE_POLICY:
 AWS::Lambda::Permission,
 AWS::IAM:Policy,
  AWS::ApplicationAutoScaling::ScalingPolicy,
 AWS::S3::BucketPolicy,
 AWS::SQS::QueuePolicy,
  and
 AWS::SNS:TopicPolicy.Applications that contain one or more nested applications require
  you to specify
 CAPABILITY_AUTO_EXPAND.If your application template contains any of the
  above resources, we recommend that you review
 all permissions associated with the
  application before deploying. If you don't specify
 this parameter for an application that
  requires capabilities, the call will fail.
- `"changeSetName"`: This property corresponds to the parameter of the same name for the
  AWS CloudFormation CreateChangeSet
  API.
- `"clientToken"`: This property corresponds to the parameter of the same name for the AWS
  CloudFormation CreateChangeSet
  API.
- `"description"`: This property corresponds to the parameter of the same name for the AWS
  CloudFormation CreateChangeSet
  API.
- `"notificationArns"`: This property corresponds to the parameter of the same name for the
  AWS CloudFormation CreateChangeSet
  API.
- `"parameterOverrides"`: A list of parameter values for the parameters of the application.
- `"resourceTypes"`: This property corresponds to the parameter of the same name for the
  AWS CloudFormation CreateChangeSet
  API.
- `"rollbackConfiguration"`: This property corresponds to the parameter of the same name
  for the AWS CloudFormation CreateChangeSet
  API.
- `"semanticVersion"`: The semantic version of the application:
 https://semver.org/
- `"tags"`: This property corresponds to the parameter of the same name for the AWS
  CloudFormation CreateChangeSet
  API.
- `"templateId"`: The UUID returned by CreateCloudFormationTemplate.Pattern:
  [0-9a-fA-F]{8}-[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12}
"""
function create_cloud_formation_change_set(
    applicationId, stackName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "POST",
        "/applications/$(applicationId)/changesets",
        Dict{String,Any}("stackName" => stackName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_cloud_formation_change_set(
    applicationId,
    stackName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "POST",
        "/applications/$(applicationId)/changesets",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("stackName" => stackName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Creates an AWS CloudFormation template.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"semanticVersion"`: The semantic version of the application:
 https://semver.org/
"""
function create_cloud_formation_template(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "POST",
        "/applications/$(applicationId)/templates";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_cloud_formation_template(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "POST",
        "/applications/$(applicationId)/templates",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Deletes the specified application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

"""
function delete_application(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "DELETE",
        "/applications/$(applicationId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_application(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "DELETE",
        "/applications/$(applicationId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Gets the specified application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"semanticVersion"`: The semantic version of the application to get.
"""
function get_application(applicationId; aws_config::AbstractAWSConfig=global_aws_config())
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_application(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the policy for the application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

"""
function get_application_policy(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/policy";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_application_policy(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/policy",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_cloud_formation_template(application_id, template_id)
    get_cloud_formation_template(application_id, template_id, params::Dict{String,<:Any})

Gets the specified AWS CloudFormation template.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.
- `template_id`: The UUID returned by CreateCloudFormationTemplate.Pattern:
  [0-9a-fA-F]{8}-[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12}

"""
function get_cloud_formation_template(
    applicationId, templateId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/templates/$(templateId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_cloud_formation_template(
    applicationId,
    templateId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/templates/$(templateId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the list of applications nested in the containing application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxItems"`: The total number of items to return.
- `"nextToken"`: A token to specify where to start paginating.
- `"semanticVersion"`: The semantic version of the application to get.
"""
function list_application_dependencies(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/dependencies";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_application_dependencies(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/dependencies",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists versions for the specified application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxItems"`: The total number of items to return.
- `"nextToken"`: A token to specify where to start paginating.
"""
function list_application_versions(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/versions";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_application_versions(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "GET",
        "/applications/$(applicationId)/versions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists applications owned by the requester.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxItems"`: The total number of items to return.
- `"nextToken"`: A token to specify where to start paginating.
"""
function list_applications(; aws_config::AbstractAWSConfig=global_aws_config())
    return serverlessapplicationrepository(
        "GET", "/applications"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_applications(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "GET",
        "/applications",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_application_policy(application_id, statements)
    put_application_policy(application_id, statements, params::Dict{String,<:Any})

Sets the permission policy for an application. For the list of actions supported for this
operation, see
 Application
 Permissions
 .

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.
- `statements`: An array of policy statements applied to the application.

"""
function put_application_policy(
    applicationId, statements; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "PUT",
        "/applications/$(applicationId)/policy",
        Dict{String,Any}("statements" => statements);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_application_policy(
    applicationId,
    statements,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "PUT",
        "/applications/$(applicationId)/policy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("statements" => statements), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    unshare_application(application_id, organization_id)
    unshare_application(application_id, organization_id, params::Dict{String,<:Any})

Unshares an application from an AWS Organization.This operation can be called only from the
organization's master account.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.
- `organization_id`: The AWS Organization ID to unshare the application from.

"""
function unshare_application(
    applicationId, organizationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "POST",
        "/applications/$(applicationId)/unshare",
        Dict{String,Any}("organizationId" => organizationId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function unshare_application(
    applicationId,
    organizationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "POST",
        "/applications/$(applicationId)/unshare",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("organizationId" => organizationId), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Updates the specified application.

# Arguments
- `application_id`: The Amazon Resource Name (ARN) of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"author"`: The name of the author publishing the app.Minimum length=1. Maximum
  length=127.Pattern \"^[a-z0-9](([a-z0-9]|-(?!-))*[a-z0-9])?\";
- `"description"`: The description of the application.Minimum length=1. Maximum length=256
- `"homePageUrl"`: A URL with more information about the application, for example the
  location of your GitHub repository for the application.
- `"labels"`: Labels to improve discovery of apps in search results.Minimum length=1.
  Maximum length=127. Maximum number of labels: 10Pattern: \"^[a-zA-Z0-9+-_:/@]+\";
- `"readmeBody"`: A text readme file in Markdown language that contains a more detailed
  description of the application and how it works.Maximum size 5 MB
- `"readmeUrl"`: A link to the readme file in Markdown language that contains a more
  detailed description of the application and how it works.Maximum size 5 MB
"""
function update_application(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return serverlessapplicationrepository(
        "PATCH",
        "/applications/$(applicationId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_application(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return serverlessapplicationrepository(
        "PATCH",
        "/applications/$(applicationId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ba49db2e-cb57-11eb-0b6d-95ae95c35719
begin
	using Pkg
	Pkg.activate("..")
	using RiceBPH
	using InteractiveDynamics
	using CairoMakie
end

# ╔═╡ 03c3bbc0-f91d-4e91-b2c7-dbcd28eceaed
using PlutoUI

# ╔═╡ ecd1ccd8-fc7f-4577-aec9-f6ac0442ee87
using DelimitedFiles

# ╔═╡ f834aec9-676e-4161-b275-6488f9a5ffd2
using Statistics

# ╔═╡ 8e1f5590-b7b6-4b2c-a3b1-1ec573e25a1f
Model = RiceBPH.Model

# ╔═╡ 93a151a3-a978-44a4-b93e-b3f3916a202a
function ac(model)
    return function (agent)
        if isnan(model.food[agent.pos...])
            return RGBf0(0, 0, 0)
        elseif agent.age < model.age_init
            "#FF0000"
        else
           	"#0000FF"
        end
    end
end

# ╔═╡ ccfcd584-f319-44cd-9220-e0f2af9fe736
const heatkwargs=(
			nan_color=RGBAf0(1, 1, 0.270588),
			colormap=[RGBAf0(0.270588, 1, 0.270588, i) for i in 0:0.01:1],
			colorrange=(0, 1),
		)

# ╔═╡ 47b279c9-19fd-4095-87f4-c477f10c4d47
function plotkwargs(model)
	(
		ac=ac(model),
		am=Model.am,
		heatarray=Model.heatarray,
		resolution=(1200,1200),
		heatkwargs=heatkwargs
	)
end

# ╔═╡ aa6ada26-a723-4198-b5c9-78e8481bc46c
base = joinpath("/tmp", "snapshot")

# ╔═╡ 686d5e42-6834-4290-9bc8-14497d399ce7
@bind seed PlutoUI.Slider(1:1000, show_value=true)

# ╔═╡ 76bc81a3-ec4f-4763-9f2d-ec31a8fa1f68
function init_model()
	model = Model.init_model(;
		envmap="../src/maps/019-1x2.csv",
		nb_bph_init=200,
		init_position=:corner,
		pr_killed=0.05,
		seed = seed
	)
end

# ╔═╡ 9b197d84-be93-4929-92a5-ad350e405435
names = let model = init_model()
	rm(base, recursive=true, force=true)
	mkpath(base)
	kw = plotkwargs(model)
	map(0:250:2880) do i
		fig, _ = abm_plot(model; kw...)
		Model.run!(model, Model.agent_step!, Model.model_step!, 300)
		name = joinpath(base, "$i.png")
		save(name, fig)
		name
	end
end

# ╔═╡ 47ed7f95-97af-4e0c-9e5e-1c1ea48e51b8
@bind plt PlutoUI.Select(reverse(names), default=names[end-1])

# ╔═╡ b9d81098-7130-4159-8f83-029d28cec1cc
LocalResource(plt)

# ╔═╡ 31fd3ba1-e3ec-4e9a-a79c-4aca3a0d33f0
collect(0:250:2880) |> clipboard

# ╔═╡ 36ec99d4-7d51-4c7d-bab3-270b2cf14eb7
map1x2 = readdlm("../src/maps/019-1x2.csv")

# ╔═╡ 48d8ee85-7b31-43f6-a492-e1858a625bfc
heatmapkwargs = (zauto=true, transpose=true, showscale=false)

# ╔═╡ bd8e91af-d793-42a1-863a-81ce019b06f4
p = CairoMakie.heatmap(transpose(map1x2); resolution=(1200,1200), heatkwargs...)

# ╔═╡ 18a1bc29-8202-4a42-8ce0-811d88ddc8bc
save("/tmp/map1x2.png", p, resolution=(1200, 1200))

# ╔═╡ 21843312-1ead-4982-9f72-1302f96001e4
LocalResource("/tmp/map1x2.png")

# ╔═╡ a4cfa06f-ce3d-4931-89d2-86db58b26333
adf = let model = init_model()
	adf, _ = Model.run!(model, Model.agent_step!, Model.model_step!, 2880, adata=Model.adata)
	adf
end

# ╔═╡ 2672bbdc-a564-44a5-8fc4-f7208243b969
let f = lines(adf.step / 24, adf.count_bph)
	save("/tmp/bph-count-before.png", f)
	f
end

# ╔═╡ 7abc52f1-1a81-4bc8-ac2c-adf383f72a6e
begin
	function ma(X, k)
		pad = zeros(eltype(X), k ÷ 2)
		[pad; [mean(X[i:(i+k)]) for i in 1:length(X) - k]; pad]
	end
	function proc(X; thres=0)
		X = ma(X, 24 * 7 ÷ 2)
		X = (X .- mean(X)) ./ std(X)
		@. X ≥ thres
	end
end

# ╔═╡ 0565bd26-81af-4740-9d9c-fc878fdde641
let f = lines(adf.step / 24, proc(adf.count_bph))
	save("/tmp/bph-count-after.png", f)
	f
end

# ╔═╡ e995b7fa-f2e7-4e85-8e54-b741057652d1
let f = Figure()
	a = Axis(f[1,1])
	b = Axis(f[1,2])
	CairoMakie.lines!(a, adf.step / 24, adf.count_bph)
	lines!(b, adf.step / 24, proc(adf.count_bph))
	f
end

# ╔═╡ Cell order:
# ╠═ba49db2e-cb57-11eb-0b6d-95ae95c35719
# ╠═8e1f5590-b7b6-4b2c-a3b1-1ec573e25a1f
# ╠═76bc81a3-ec4f-4763-9f2d-ec31a8fa1f68
# ╠═93a151a3-a978-44a4-b93e-b3f3916a202a
# ╠═ccfcd584-f319-44cd-9220-e0f2af9fe736
# ╠═47b279c9-19fd-4095-87f4-c477f10c4d47
# ╠═aa6ada26-a723-4198-b5c9-78e8481bc46c
# ╠═9b197d84-be93-4929-92a5-ad350e405435
# ╠═03c3bbc0-f91d-4e91-b2c7-dbcd28eceaed
# ╟─686d5e42-6834-4290-9bc8-14497d399ce7
# ╠═47ed7f95-97af-4e0c-9e5e-1c1ea48e51b8
# ╠═b9d81098-7130-4159-8f83-029d28cec1cc
# ╠═31fd3ba1-e3ec-4e9a-a79c-4aca3a0d33f0
# ╠═ecd1ccd8-fc7f-4577-aec9-f6ac0442ee87
# ╠═36ec99d4-7d51-4c7d-bab3-270b2cf14eb7
# ╠═48d8ee85-7b31-43f6-a492-e1858a625bfc
# ╠═bd8e91af-d793-42a1-863a-81ce019b06f4
# ╠═18a1bc29-8202-4a42-8ce0-811d88ddc8bc
# ╠═21843312-1ead-4982-9f72-1302f96001e4
# ╠═a4cfa06f-ce3d-4931-89d2-86db58b26333
# ╠═f834aec9-676e-4161-b275-6488f9a5ffd2
# ╠═2672bbdc-a564-44a5-8fc4-f7208243b969
# ╠═0565bd26-81af-4740-9d9c-fc878fdde641
# ╠═e995b7fa-f2e7-4e85-8e54-b741057652d1
# ╠═7abc52f1-1a81-4bc8-ac2c-adf383f72a6e

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
# TODO: make this the definition for siteind
"""
    firstsiteind(M::MPO, j::Int; plev = 0, kwargs...)

    siteind(M::MPO, j::Int; plev = 0, kwargs...)

Return the first site Index found of the MPO.

By default, it returns the first site Index with prime
level of 0. You can choose different filters, like tags,
with the `kwargs`.
"""
function firstsiteind(M::Union{MPS, MPO}, j::Int; plev::Int = 0, kwargs...)
  N = length(M)
  (N==1) && return firstind(M[1]; plev = plev, kwargs...)

  if j == 1
    si = uniqueind(M[j], M[j+1]; plev = plev, kwargs...)
  elseif j == N
    si = uniqueind(M[j], M[j-1]; plev = plev, kwargs...)
  else
    si = uniqueind(M[j], M[j-1], M[j+1]; plev = plev, kwargs...)
  end
  return si
end

# TODO: make this alias?
#const firstsiteind = siteind

# TODO: rename to allsiteinds?
"""
    siteinds(M::MPO, j::Int; kwargs...)

Return the site Indices found of the MPO.
"""
function ITensors.siteinds(M::Union{MPS, MPO}, j::Int; kwargs...)
  N = length(M)
  (N==1) && return inds(M[1]; kwargs...)

  if j == 1
    si = uniqueinds(M[j], M[j+1]; kwargs...)
  elseif j == N
    si = uniqueinds(M[j], M[j-1]; kwargs...)
  else
    si = uniqueinds(M[j], M[j-1], M[j+1]; kwargs...)
  end
  return si
end

# TODO: rename to allsiteinds?
"""
    siteinds(M::MPO; kwargs...)

Get a Vector of the site indices of M.
"""
ITensors.siteinds(M::MPO; kwargs...) =
  [siteinds(M, j; kwargs...) for j in 1:length(M)]

"""
    firstsiteinds(M::Union{MPS, MPO}; kwargs...)

Get a Vector of the first site Index found on each site of M.
"""
firstsiteinds(M::Union{MPS, MPO}; kwargs...) =
  [firstsiteind(M, j; kwargs...) for j in 1:length(M)]

function findsiteinds(ψ::Union{MPS, MPO},
                      inds::Vector)
  return [findfirstsiteinds(ψ, inds[n]) for n in 1:length(ψ)]
end

function findcommonsiteinds(o::ITensor,
                            ψ::Union{MPS, MPO})
  ns = Int[]
  for (n, ψn) in enumerate(ψ)
    if hascommoninds(o, ψn)
      push!(ns, n)
    end
  end
  ns
end

"""
    movesite(::Union{MPS, MPO}, n1::Int, n2::Int)

Create a new MPS where the site at `n1` is moved to `n2`
"""
function movesite(ψ::Union{MPS, MPO},
                  n1::Int, n2::Int;
                  orthocenter::Int = n2,
                  kwargs...)
  n1 == n2 && return copy(ψ)
  ψ = orthogonalize(ψ, n2)
  r = n1:n2-1
  ortho = "left"
  if n1 > n2
    r = reverse(n2:n1-1)
    ortho = "right"
  end
  for n in r
    ψ = swapbondsites(ψ, n; ortho = ortho, kwargs...)
  end
  ψ = orthogonalize(ψ, orthocenter)
  return ψ
end

function movesite(ns::Vector{Int},
                  n1::Int, n2::Int)
  n1 == n2 && return copy(ns)
  r = n1:n2-1
  if n1 > n2
    r = reverse(n2:n1-1)
  end
  for n in r
    #ns = swapbondsites(ns, n)
    ns = replace(ns, n => n+1, n+1 => n)
  end
  return ns
end

# TODO: make a permutesites(::MPS/MPO, perm)
# function that that a permutation of the sites
# p(1:N) for N sites
function movesites(ψ::Union{MPS, MPO},
                   ns, ns′; kwargs...)
  ψ = copy(ψ)
  N = length(ns)
  @assert N == length(ns′)
  p = sortperm(ns′)
  ns = ns[p]
  ns′ = ns′[p]
  ns = collect(ns)
  for i in 1:N
    ψ = movesite(ψ, ns[i], ns′[i]; kwargs...)
    ns = movesite(ns, ns[i], ns′[i])
  end
  return ψ
end

"""
    swapbondsites(ψ::MPO, b::Int; kwargs...)

Swap the sites `b` and `b+1`.
"""
function ITensors.swapbondsites(ψ::MPO, b::Int; kwargs...)
  ortho = get(kwargs, :ortho, "right")
  if ortho == "left"
    orthocenter = b+1
  elseif ortho == "right"
    orthocenter = b
  end
  if ITensors.leftlim(ψ) < b - 1
    ψ = orthogonalize(ψ, b)
  elseif ITensors.rightlim(ψ) > b + 2
    ψ = orthogonalize(ψ, b + 1)
  end
  #kwargs = setindex(values(kwargs), [2, 1], :parm)
  ψ = replacesites(ψ, ψ[b] * ψ[b+1];
                   firstsite = b, lastsite = b+1,
                   orthocenter = orthocenter,
                   perm = [2, 1],
                   kwargs...)
  return ψ
end

"""
    MPS(::ITensor, sites)

Construct an MPS from an ITensor by decomposing it site
by site.
"""
function MPS(A::ITensor, sites;
             firstlinkind::Union{Nothing, Index} = nothing,
             lastlinkind::Union{Nothing, Index} = nothing,
             orthocenter::Int = length(sites),
             kwargs...)
  N = length(sites)
  @assert order(A) == N + !isnothing(firstlinkind) + !isnothing(lastlinkind)
  @assert hasinds(A, sites)
  @assert isnothing(firstlinkind) || hasind(A, firstlinkind)
  @assert isnothing(lastlinkind) || hasind(A, lastlinkind)

  # TODO: generalize to other orthocenters.
  # To minimize work, may need to loop from
  # 1:orthocenter and reverse(orthocenter:N)
  @assert orthocenter == N

  ψ = Vector{ITensor}(undef, N)
  Ã = A
  l = firstlinkind
  for n in 1:N-1
    s = sites[n]
    Lis = isnothing(l) ? (s,) : (l, s)
    L, R = factorize(Ã, Lis; kwargs...)
    l = commonind(L, R)
    ψ[n] = L
    Ã = R
  end
  ψ[N] = Ã
  return MPS(ψ)
end

"""
    MPO(::ITensor, sites)

Construct an MPO from an ITensor by decomposing it site
by site.
"""
function MPO(A::ITensor, sites;
             firstlinkind::Union{Nothing, Index} = nothing,
             lastlinkind::Union{Nothing, Index} = nothing,
             orthocenter::Int = length(sites),
             kwargs...)
  N = length(sites)
  @assert order(A) == 2*N + !isnothing(firstlinkind) + !isnothing(lastlinkind)
  @assert hasinds(A, sites)
  @assert hasinds(A, prime.(sites))
  @assert isnothing(firstlinkind) || hasind(A, firstlinkind)
  @assert isnothing(lastlinkind) || hasind(A, lastlinkind)

  # TODO: generalize to other orthocenters.
  # To minimize work, may need to loop from
  # 1:orthocenter and reverse(orthocenter:N)
  @assert orthocenter == N

  ψ = Vector{ITensor}(undef, N)
  Ã = A
  l = firstlinkind
  for n in 1:N-1
    s = sites[n]
    Lis = isnothing(l) ? (s, s') : (l, s, s')
    L, R = factorize(Ã, Lis; kwargs...)
    l = commonind(L, R)
    ψ[n] = L
    Ã = R
  end
  ψ[N] = Ã
  M = MPO(ψ)
  ITensors.setleftlim!(M, orthocenter-1)
  ITensors.setrightlim!(M, orthocenter+1)
  return M
end

"""
    setindex!(::Union{MPS, MPO}, ::Union{MPS, MPO},
              r::UnitRange{Int64})

Sets a contiguous range of MPS/MPO tensors
"""
function Base.setindex!(ψ::MPST, ϕ::MPST,
                        r::UnitRange{Int64}) where {MPST <: Union{MPS, MPO}}
  @assert length(r) == length(ϕ)
  # TODO: accept r::Union{AbstractRange{Int}, Vector{Int}}
  # if r isa AbstractRange
  #   @assert step(r) = 1
  # else
  #   all(==(1), diff(r))
  # end
  llim = ITensors.leftlim(ψ)
  rlim = ITensors.rightlim(ψ)
  for (j, n) in enumerate(r)
    ψ[n] = ϕ[j]
  end
  if llim + 1 ≥ r[1]
    ITensors.setleftlim!(ψ, ITensors.leftlim(ϕ) + r[1] - 1)
  end
  if rlim - 1 ≤ r[end]
    ITensors.setrightlim!(ψ, ITensors.rightlim(ϕ) + r[1] - 1)
  end
  return ψ
end

# TODO: add a version of replacesites that
# determines `firstsite` and `lastsite` from common
# site indices of ψ and A
"""
    replacesites(ψ::Union{MPS, MPO},
                 A::ITensor;
                 firstsite::Int,
                 lastsite::Int,
                 orthocenter::Int = firstsite,
                 perm = nothing,
                 kwargs...)

Replace the sites between `firstsite` and `lastsite` with tensors made
from decomposing `A` into an MPS.

The MPS must be orthogonalized such that
```
firstsite ≤ ITensors.orthocenter(ψ) ≤ lastsite
```

Choose the new orthogonality center with `orthocenter`, which
should be within `firstsite:lastsite`.

Optionally, permute the order of the sites with `perm`.
"""
function replacesites(ψ::MPST,
                      A::ITensor;
                      firstsite::Int,
                      lastsite::Int,
                      orthocenter::Int = lastsite,
                      perm = nothing,
                      kwargs...) where {MPST <: Union{MPS, MPO}}
  # Replace the sites of ITensor ψ
  # with the tensor A, splitting up A
  # into MPS tensors
  ψ = copy(ψ)

  @assert firstsite ≤ ITensors.orthocenter(ψ) ≤ lastsite
  @assert firstsite ≤ ITensors.leftlim(ψ) + 1
  @assert ITensors.rightlim(ψ) - 1 ≤ lastsite

  # TODO: allow orthocenter outside of this
  # range, and orthogonalize/truncate as needed
  @assert firstsite ≤ orthocenter ≤ lastsite

  # Check that A has the proper common
  # indices with ψ
  l = linkind(ψ, firstsite-1)
  r = linkind(ψ, lastsite)
  # TODO: replace with commonsiteinds(::MPST, ::ITensor)?
  sites = collect(Iterators.flatten([siteinds(ψ, j) for j in firstsite:lastsite]))
  indsA = filter(x -> !isnothing(x), [l, sites..., r])

  @assert hassameinds(A, indsA)

  # For MPO case, restrict to 0 prime level
  sites = filter(hasplev(0), sites)

  if !isnothing(perm)
    sites = sites[[perm...]]
  end

  ψA = MPST(A, sites; firstlinkind = l,
                      lastlinkind = r,
                      kwargs...)
  #@assert prod(ψA) ≈ A

  ψ[firstsite:lastsite] = ψA

  return ψ
end

"""
    MPO([::Type{ElT} = Float64}, ]sites, ops::Vector{String})

Make an MPO with pairs of sites `s[i]` and `s[i]'`
and operators `ops` on each site.
"""
function ITensors.MPO(::Type{ElT},
             sites::Vector{<:Index},
             ops::Vector{String}) where {ElT <: Number}
  N = length(sites)
  its = Vector{ITensor}(undef, N)
  links = Vector{Index}(undef, N)
  for ii ∈ eachindex(sites)
    si = sites[ii]
    d = dim(si)
    spin_op = op(sites, ops[ii], ii)
    if hasqns(si)
      links[ii] = Index([QN() => 1], "Link,n=$ii")
    else
      links[ii] = Index(1, "Link,n=$ii")
    end
    local this_it
    if ii == 1
      this_it = emptyITensor(ElT, links[ii], si', dag(si))
      for jj in 1:d, jjp in 1:d
        this_it[links[ii](1), si[jj], si'[jjp]] = spin_op[si[jj], si'[jjp]]
      end
    elseif ii == N
      this_it = emptyITensor(ElT, dag(links[ii-1]), si', dag(si))
      for jj in 1:d, jjp in 1:d
        this_it[links[ii-1](1), si[jj], si'[jjp]] = spin_op[si[jj], si'[jjp]]
      end
    else
      this_it = emptyITensor(ElT, dag(links[ii-1]), links[ii], si', dag(si))
      for jj in 1:d, jjp in 1:d
        this_it[links[ii-1](1),
                links[ii](1),
                si[jj],
                si'[jjp]] = spin_op[si[jj], si'[jjp]]
      end
    end
    its[ii] = this_it
  end
  MPO(its)
end


text:
#Function for getting solution to radial Schrodinger equation using numerov algorithm
#23/07/2017

function numerovfunc(atom,nn,ll,jj)

#Get parameters
ZZ, spin, eneigval, params = GetAtomParams(atom,nn,ll,jj)
alpha_c = params[1]
a_1 = params[2]
a_2 = params[3]
a_3 = params[4]
a_4 = params[5]
r_c = params[6]

#Create vector to integrate over
hh = 0.01 #Choose this to be small so that O(hh^6) may be ignored
h2 = hh^2 #For efficiency
r_i = hh/2 #Define starting point
r_o = 2*nn*(nn+15)  #Define end point
x_i = log(r_i) #Change of variables from Callagher, 2005
x_o = log(r_o)
xx = collect(x_i:hh:x_o)
rr = exp(xx) #For efficiency

#Set up Schrodinger equation from Gallagher (2005), note Pritchard (2012)

LdotS = (jj*(jj+1)-ll*(ll+1)-spin*(spin+1))/2 #Spin-orbit coupling

spinorbitpotential = finestrucconst^2*LdotS./(2*rr.^3) #Fine structure splitting

radialcharge = 1 + (ZZ-1)*exp(-a_1*rr) - rr.*(a_3+a_4*rr).*exp(-a_2*rr) #Effective nuclear charge

coulombpotential = -radialcharge./rr - (alpha_c./(2*(rr.^4))).*(1-exp(-((rr/r_c).^6))) #Corrected Coulomb potential

totalpotential = spinorbitpotential+coulombpotential #total potential

cenfugterm = (ll+1/2)^2 #centrifugal term from Schrodinger equation

#Apply Numerov algorithm
G_x = 2*exp(2*xx).*(totalpotential-eneigval) + cenfugterm #Coeficient in differential equation

T_x = h2*G_x/12 #For efficiency

Ysoln_vec = zeros(size(xx)) #Generate a placeholder vector for solutions

Ysoln_vec[end] = -1e-10 #The function must go to zero at infinity

Ysoln_vec[end-1] = -2e-10

#To perform the iteration with convenient indexing the vectors xx, G_x and soln_vec are flipped

fYsoln = flipdim(Ysoln_vec,1)
fT_x = flipdim(T_x,1)

for ii = 3:length(xx)
        fYsoln[ii] = ((2+10*fT_x[ii-1])*fYsoln[ii-1] -  (1-fT_x[ii-2])*fYsoln[ii-2])/(1-fT_x[ii])
end

Ysoln_vec = flipdim(fYsoln,1) #Return solutions to proper ordering

#Normalise (method adapted from Gallagher, 2005)
normconstvec = zeros(length(rr)-1)
for mm = 2:length(rr)
    deltar = rr[mm]-rr[mm-1]
    normconstvec[mm-1] = (Ysoln_vec[mm]^2)*rr[mm]*deltar
end

normconst = sqrt(sum(normconstvec))
normY_sol = Ysoln_vec/normconst

#Saving and loading seems to hurt performance
#filename = filenamemake(atom,nn,ll,jj)

#save(filename,"rr",rr,"normY_sol",normY_sol)

return normY_sol, rr
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
# General structure for reachability methods

# Performs layer-by-layer propagation
# It is called by all solvers under reachability
# TODO: also called by ReluVal and FastLin, so move to general utils (or to network.jl)
function forward_network(solver, nnet::Network, input)
    reach = input
    for layer in nnet.layers
        reach = forward_layer(solver, layer, reach)
    end
    return reach
end

# Checks whether the reachable set belongs to the output constraint
# It is called by all solvers under reachability
# Note vertices_list is not defined for HPolytope: to be defined
function check_inclusion(reach::Vector{<:LazySet}, output)
    for poly in reach
        issubset(poly, output) || return ReachabilityResult(:violated, reach)
    end
    return ReachabilityResult(:holds, reach)
end

function check_inclusion(reach::P, output) where P<:LazySet
    return ReachabilityResult(issubset(reach, output) ? :holds : :violated, [reach])
end

# return a vector so that append! is consistent with the relu forward_partition
forward_partition(act::Id, input) = [input]

function forward_partition(act::ReLU, input)
    N = dim(input)
    N > 30 && @warn "Got dim(X) == $N in `forward_partition`. Expecting 2ᴺ = $(2^big(N)) output sets."

    output = HPolytope{Float64}[]

    for h in 0:(big"2"^N)-1
        P = Diagonal(1.0.*digits(h, base = 2, pad = N))
        orthant = HPolytope(Matrix(I - 2.0P), zeros(N))
        S = intersection(input, orthant)
        if !isempty(S)
            push!(output, linear_map(P, S))
        end
    end
    return output
end

###################### For ReluVal and Neurify ######################

struct SymbolicInterval{F<:AbstractPolytope}
    Low::Matrix{Float64}
    Up::Matrix{Float64}
    domain::F
end

# Data to be passed during forward_layer
struct SymbolicIntervalGradient{F<:AbstractPolytope, N<:Real}
    sym::SymbolicInterval{F}
    LΛ::Vector{Vector{N}}
    UΛ::Vector{Vector{N}}
end
# Data to be passed during forward_layer
const SymbolicIntervalMask = SymbolicIntervalGradient{<:Hyperrectangle, Bool}

function _init_symbolic_grad_general(domain, N)
    n = dim(domain)
    I = Matrix{Float64}(LinearAlgebra.I(n))
    Z = zeros(n)
    symbolic_input = SymbolicInterval([I Z], [I Z], domain)
    symbolic_mask = SymbolicIntervalGradient(symbolic_input,
                                             Vector{Vector{N}}(),
                                             Vector{Vector{N}}())
end
function init_symbolic_grad(domain)
    VF = Vector{HalfSpace{Float64, Vector{Float64}}}
    domain = HPolytope(VF(constraints_list(domain)))
    _init_symbolic_grad_general(domain, Float64)
end
function init_symbolic_mask(interval)
    _init_symbolic_grad_general(interval, Bool)
end

domain(sym::SymbolicInterval) = sym.domain
domain(grad::SymbolicIntervalGradient) = domain(grad.sym)
_sym(sym::SymbolicInterval) = sym
_sym(grad::SymbolicIntervalGradient) = grad.sym

upper(sym::SymbolicInterval) = AffineMap(@view(sym.Up[:, 1:end-1]), sym.domain, @view(sym.Up[:, end]))
lower(sym::SymbolicInterval) = AffineMap(@view(sym.Low[:, 1:end-1]), sym.domain, @view(sym.Low[:, end]))
upper(grad::SymbolicIntervalGradient) = upper(grad.sym)
lower(grad::SymbolicIntervalGradient) = lower(grad.sym)

upper_bound(a::AbstractVector, set::LazySet) = a'σ(a, set)
lower_bound(a::AbstractVector, set::LazySet) = a'σ(-a, set) # ≡ -ρ(-a, set)
bounds(a::AbstractVector, set::LazySet) = (a'σ(-a, set), a'σ(a, set))  # (lower, upper)

upper_bound(S::LazySet, j::Integer) = upper_bound(Arrays.SingleEntryVector(j, dim(S), 1.0), S)
lower_bound(S::LazySet, j::Integer) = lower_bound(Arrays.SingleEntryVector(j, dim(S), 1.0), S)
bounds(S::LazySet, j::Integer) = (lower_bound(S, j), upper_bound(S, j))

const _SymIntOrGrad = Union{SymbolicInterval, SymbolicIntervalGradient}
LazySets.dim(sym::_SymIntOrGrad) = size(_sym(sym).Up, 1)
LazySets.high(sym::_SymIntOrGrad) = [upper_bound(upper(sym), j) for j in 1:dim(sym)]
LazySets.low(sym::_SymIntOrGrad) = [lower_bound(lower(sym), j) for j in 1:dim(sym)]
# radius of the symbolic interval in the direction of the
# jth generating vector. This is not the axis aligned radius,
# or the bounding radius, but rather a radius with respect to
# a node in the network. Equivalent to the upper-upper
# bound minus the lower-lower bound
function radius(sym::_SymIntOrGrad, j::Integer)
    upper_bound(upper(sym), j) - lower_bound(lower(sym), j)
end


text:
using MyPkg1
using Test

# @test f1(2,1)==6

@testset "MyPkg1.jl" begin
    @test f1(2,1)==6
    @test f1(3,2)==10
end


text:
# This file was generated, do not modify it. # hide
# hideall

using CSV, DataFrames, Latexify, Formatting
E = CSV.read("emojis.tsv.gz", DataFrame, delim='\t')
cclist = unique(E.country_code)
sort!(cclist)

voclist = Dict()

for g in groupby(E, :country_code)
    D = Dict{String7,Float64}()
    voclist[g.country_code[1]] = D
    for row in eachrow(g)
         D[row.emoji] = row.freq
    end
end

k = 20
# plot()
table = Matrix{String}(undef, length(cclist) + 1, k+1)
header = ["region"]
for i in 1:k
    push!(header, string(i))
end

table[1, :] .= header
for (i, cc) in enumerate(cclist)
    voc = collect(voclist[cc])
    sort!(voc, by=p -> p[end], rev=true)
    resize!(voc, k)
    # display("text/markdown", "$cc: " * join(first.(voc), ", "))
    table[i+1, 1] = cc
    table[i+1, 2:end] .= first.(voc)
end
# display(latexify(table))

#plot!(title="raw frequencies per region", yscale=:log10, ylabel="frequency", xlabel="emoji rank", legend=:outertopright, fmt=:png) |> display
println(latexify(table, latex=false, env=:mdtable))


text:
using MarkdownStyles
using Documenter

DocMeta.setdocmeta!(MarkdownStyles, :DocTestSetup, :(using MarkdownStyles); recursive=true)

makedocs(;
    modules=[MarkdownStyles],
    authors="Felipe S. S. Schneider <schneider.felipe@posgrad.ufsc.br> and contributors",
    repo="https://github.com/schneiderfelipe/MarkdownStyles.jl/blob/{commit}{path}#{line}",
    sitename="MarkdownStyles.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://schneiderfelipe.github.io/MarkdownStyles.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/schneiderfelipe/MarkdownStyles.jl",
)


text:
function vhem_step_E(base::H3M{Z}, reduced::H3M{Z}, τ::Integer, N::Integer) where Z
    ## Expectations
    # logη[i,j][β,ρ][m,l]
    logη = Matrix{Matrix{Matrix{Float64}}}(undef, length(base.M), length(reduced.M))

    # lhmm[i,j]: expected log-likelihood
    lhmm = zeros(length(base.M), length(reduced.M))

    ## Aggregated summaries
    ## Coviello14, p. 15, eqn. 23

    # νagg1[i,j][ρ]: expected number of times that
    # HMM j starts from ρ, when modeling sequences
    # generated by HMM i.
    νagg1 = Matrix{Vector{Float64}}(undef, length(base.M), length(reduced.M))

    # νagg[i,j][ρ,β]: expected number of times that
    # HMM j is in state ρ when HMM i is in state β,
    # when both HMMs are modeling sequences generated
    # by HMM i.
    νagg = Matrix{Matrix{Float64}}(undef, length(base.M), length(reduced.M))

    # ξagg[i,j][ρ,ρ']: expected number of transitions
    # from state ρ to ρ' of the HMM j when modeling
    # sequences generated by HMM i.
    ξagg = Matrix{Matrix{Float64}}(undef, length(base.M), length(reduced.M))

    ## Optimal assignments probabilities
    logz = zeros(length(base.M), length(reduced.M))

    # TODO: "Summary statistics" struct

    logωj = log.(reduced.ω)

    Threads.@threads for i in OneTo(length(base))
        Mi = base.M[i]

        # Initial probabilities (π in the paper)
        # Transition matrices (a or A in the paper)
        logai = log.(Mi.a)
        logAi = log.(Mi.A)

        logzacc = LogSumExpAcc()

        for (j, Mj) in enumerate(reduced.M)
            Ki, Kj = size(Mi, 1), size(Mj, 1)

            ## Summary statistics
            logν = zeros(τ, Kj, Ki)
            logξ = zeros(τ, Kj, Kj, Ki)

            ## Aggregate summaries
            νagg1[i,j] = zeros(Kj)
            νagg[i,j]  = zeros(Kj, Ki)
            ξagg[i,j]  = zeros(Kj, Kj)

            logη[i,j], _, (lhmm[i,j], logϕ, logϕ1) = loglikelihood_va(Mi, Mj, τ)

            for ρ in OneTo(Kj), β in OneTo(Ki)
                logν[1,ρ,β]     = logai[β] + logϕ1[β,ρ]
                νagg1[i,j][ρ]  += exp(logν[1,ρ,β])
                νagg[i,j][ρ,β] += exp(logν[1,ρ,β])
            end

            logtmps = zeros(Kj, Ki)

            for t in 2:τ
                for β in OneTo(Ki)
                    for ρp in OneTo(Kj)
                        acc = LogSumExpAcc()
                        for βp in OneTo(Ki)
                            add!(acc, logν[t-1,ρp,βp] + logAi[βp,β])
                        end
                        logtmps[ρp,β] = sum(acc)
                    end
                end

                for ρ in OneTo(Kj)
                    for β in OneTo(Ki)
                        acc = LogSumExpAcc()
                        for ρp in OneTo(Kj)
                            logξ[t,ρp,ρ,β] = logtmps[ρp,β] + logϕ[t,β,ρp,ρ]
                            ξagg[i,j][ρp,ρ] += exp(logξ[t,ρp,ρ,β])
                            add!(acc, logξ[t,ρp,ρ,β])
                        end
                        logν[t,ρ,β] = sum(acc)
                        νagg[i,j][ρ,β] += exp(logν[t,ρ,β])
                    end
                end
            end

            # Compute optimal assignment probabilities
            logz[i,j] = logωj[j] + (N * base.ω[i] * lhmm[i,j])
            add!(logzacc, logz[i,j])
        end

        logz[i,:] .-= sum(logzacc)
    end

    logz, logη, lhmm, νagg, νagg1, ξagg
end

function vhem_step(base::H3M{Z}, reduced::H3M{Z}, τ::Integer, N::Integer) where Z
    ## E step
    logz, logη, lhmm, νagg, νagg1, ξagg = vhem_step_E(base, reduced, τ, N)

    # TODO: Use logz below instead ?
    z = exp.(logz)

    ## M step

    # 1. Compute H3M weights
    newω = zeros(length(reduced))
    normω = 0.0

    for j in OneTo(length(reduced)), i in OneTo(length(base))
        newω[j] += z[i,j]
        normω += z[i,j]
    end

    for j in OneTo(length(reduced))
        newω[j] /= normω
    end

    # 2. Compute H3M models
    newM = Vector{Z}(undef, length(reduced))

    # NOTE: Make sure that no variable is capture inside the loop!
    # That's why we have `normω` above and `norm` below.
    Threads.@threads for j in OneTo(length(reduced))
        Mj, ωj = reduced.M[j], reduced.ω[j]
        Kj = size(Mj, 1)

        # 2.a Compute initial probabilities
        newa = zeros(Kj)
        norm = 0.0

        for ρ in OneTo(Kj)
            for i in OneTo(length(base))
                newa[ρ] += z[i,j] * base.ω[i] * νagg1[i,j][ρ]
            end
            norm += newa[ρ]
        end

        newa /= norm

        # 2.b Compute transition matrix
        newA = zeros(Kj, Kj)

        for ρ in OneTo(Kj)
            norm = 0.0
            for ρp in OneTo(Kj)
                for i in OneTo(length(base))
                    newA[ρ,ρp] += z[i,j] * base.ω[i] * ξagg[i,j][ρ,ρp]
                end
                norm += newA[ρ,ρp]
            end

            for ρp in OneTo(Kj)
                newA[ρ,ρp] /= norm
            end
        end

        # 2.c Compute observation mixtures
        newB = Vector{UnivariateDistribution}(undef, Kj)

        for (ρ, Mjρ) in enumerate(Mj.B)
            norm = 0.0
            newc = zeros(ncomponents(Mjρ))
            newd = Vector{Normal}(undef, length(Mj.B[ρ].components))

            for (l, Mjρl) in enumerate(Mj.B[ρ].components)
                newc[l] = Ω(base, j, ρ, z, νagg) do i, β, m
                    exp(logη[i,j][β,ρ][m,l])
                end

                newμ = Ω(base, j, ρ, z, νagg) do i, β, m
                    exp(logη[i,j][β,ρ][m,l]) * base.M[i].B[β].components[m].μ
                end

                newσ2 = Ω(base, j, ρ, z, νagg) do i, β, m
                    exp(logη[i,j][β,ρ][m,l]) * (base.M[i].B[β].components[m].σ^2 + (base.M[i].B[β].components[m].μ - Mjρl.μ)^2)
                end

                newμ  /= newc[l]
                newσ2 /= newc[l]
                norm  += newc[l]

                newd[l] = Normal(newμ, sqrt(newσ2))
            end

            newc /= norm
            newB[ρ] = MixtureModel(newd, newc)
        end

        # 2.d Build HMM
        newM[j] = HMM(newa, newA, newB)
    end

    H3M(newM, newω), lhmm, z
end

# TODO: Cleanup (perform computations in the M-step loop instead ?)
function Ω(f, b::H3M, j::Integer, ρ::Integer, z::AbstractMatrix, νagg)
    tot = 0.0
    for (i, ωi) in enumerate(b.ω)
        s = 0.0
        for β in 1:size(b.M[i],1)
            ss = 0.0
            for (m, c) in enumerate(b.M[i].B[β].prior.p)
                ss += c * f(i, β, m)
            end
            s += νagg[i,j][ρ,β] * ss
        end
        tot += z[i,j] * s
    end
    tot
end


text:
module LibSodium

using libsodium_jll

const OUTPUT_DIR = abspath(@__DIR__, "..", "gen", "out")

include(joinpath(OUTPUT_DIR, "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t
include(joinpath(OUTPUT_DIR, "common.jl"))
include(joinpath(OUTPUT_DIR, "api.jl"))

foreach(names(@__MODULE__, all=true)) do s
    if any([startswith(string(s), prefix) for prefix in ["sodium_", "crypto_"]])
        @eval export $s
    end
end

end


text:
@testset "▶ /object                     " begin
    f = Figure()
    text("blah", (0,0))
    @test gca().objects[1].anchor == (0.0,0.0)
    @test gca().objects[1].text == "blah"

    # vline, hline
    clo()
    hline(0.5; lw=0.1)
    vline(0.2; col="red")
    @test gca().objects[1] isa G.StraightLine2D
    @test gca().objects[1].anchor == 0.5
    @test gca().objects[1].horiz == true
    @test gca().objects[1].linestyle.lwidth == 0.1
    @test gca().objects[2].anchor == 0.2
    @test gca().objects[2].linestyle.color == c"red"
    @test gca().objects[2].horiz == false

    # box
    clo()
    box((0.5, 0.2), (0.0,0.0); col="red", fill="blue")
    @test gca().objects[1] isa G.Box2D
    @test gca().objects[1].anchor == (0.0,0.0)
    @test gca().objects[1].size == (0.5,0.2)
    @test gca().objects[1].fillstyle.fill == c"blue"
    @test gca().objects[1].linestyle.color == c"red"
    @test gca().objects[1].nobox == true

    # colorbar
    clo()
    colorbar(0, 1, [c"blue", c"red", c"green"])
    @test gca().objects[1] isa G.Colorbar
    @test gca().objects[1].zmin == 0.0
    @test gca().objects[1].zmax == 1.0
    @test gca().objects[1].cmap == [c"blue", c"red", c"green"]
    @test gca().objects[1].ticks.places == collect(range(0,stop=1,length=5)[2:end-1])
    @test isnothing(gca().objects[1].size)
    @test gca().objects[1].pixels == 100
    @test gca().objects[1].nobox == true
    @test gca().objects[1].position == "right"
    @test gca().objects[1].offset == (0.3,0.0)

    # clf()
    # X = randn(2,2)
    # h = heatmap(X)
    # colorbar(h)
    # @test gca().objects[1].zmin == minimum(X)
    # @test gca().objects[1].zmax == maximum(X)
end

@testset "▶ apply_gle/object            " begin
    f = Figure()

    begin # TEXT
        clf()
        text("blah", (0,0))
        s = G.assemble_figure(f, debug=true)
        isin(s, "\ngsave")
        isin(s, "\nset just cc")
        isin(s, "\namove xg(0.0) yg(0.0)")
        isin(s, "\nwrite \"blah\"")
        isin(s, "\ngrestore")
    end

    begin # STRAIGHTLINE
        clf()
        plot([0, 1],[0, 1])
        hline(0.1; col="red")
        vline(0.5; lw=0.5)
        s = G.assemble_figure(f, debug=true)
        isin(s, "d1 line color") # plot
        # plot horizontal line
        isin(s, "gsave")
        isin(s, "set color rgba(1.0,0.0,0.0,1.0)")
        isin(s, "amove xg(xgmin) yg(0.1)")
        isin(s, "aline xg(xgmax) yg(0.1)")
        isin(s, "grestore")
        isin(s, "set lwidth 0.5")
        isin(s, "amove xg(0.5) yg(ygmin)")
        isin(s, "aline xg(0.5) yg(ygmax)")
    end

    begin # BOX2D
        clf()
        plot([0,1],[0,1])
        box((0.4,0.6),(0.3,0.7);fill="red",alpha=0.1,col="blue",nobox=false)
        s = G.assemble_figure(f, debug=true)
        isin(s, "gsave")
        isin(s, "set color rgba(0.0,0.0,1.0,1.0)") # blue
        isin(s, "box xg(0.4)-xg(0) yg(0.6)-yg(0) fill rgba(1.0,0.0,0.0,0.102)")
    end

    begin # COLORBAR
        clf()

        X = [ 0.112125   1.83719    0.610405    -0.17842     0.117574;
             0.340722   0.382532   1.0889      -0.0145842  -0.293014;
             1.39024    0.140773   1.28426      0.79504     0.172258;
             0.519013   0.471506  -0.904603     1.11415     0.359543;
             0.813583  -0.495029  -0.00136247  -0.807024    1.04971 ]

        h = heatmap(X)

        # default colorbar on the right side
        colorbar(h)
        s = G.assemble_figure(f, debug=true)
        for λ ∈ ["amove xg(xgmax)+0.3 yg(ygmin)+0.0",
                 "begin box name cmap nobox",
                 "colormap \"y\" 0 1 0 1 1 100 0.25 abs(yg(ygmax)-yg(ygmin)) palette",
                 "end box",
                 # ticks lines
                 "amove xg(xgmax)+0.3+0.25 yg(ygmin)+0.0+abs(yg(ygmax)-yg(ygmin))*0.25",
                 "rline 0.25/3 0",
                 "amove xg(xgmax)+0.3+0.25 yg(ygmin)+0.0+abs(yg(ygmax)-yg(ygmin))*0.5",
                 "amove xg(xgmax)+0.3+0.25 yg(ygmin)+0.0+abs(yg(ygmax)-yg(ygmin))*0.749",
                 # # ticks writing
                 "set just lc"
                 ]
            isin(s, λ)
        end
    end
end


text:
# Load packages
using Dojo
using Plots
using Random
using MeshCat

# Open visualizer
vis= Visualizer()
open(vis)

# Include new files
include( "../utils.jl")

mech = get_mechanism(:box, timestep=0.05, g=-9.81, friction_coefficient=0.2, radius=0.00, side=0.50, mode=:box);
initialize!(mech, :box, x=[0,-1,1.], v=[0,2,1.], ω=[2,5,10.])
storage = simulate!(mech, 5.0, record=true,
    opts=SolverOptions(btol=1e-6, rtol=1e-6, verbose=false))
visualize(mech, storage, vis=vis, show_contact=false)

################################################################################
# Generate & Save Dataset
################################################################################
init_kwargs = Dict(:xlims => [[0,0,0.2], [1,1,0.4]],
				   :vlims => [-2ones(3), [2,2,-1.]],
				   :ωlims => [-6ones(3), 6ones(3)])
mech_kwargs = Dict(:friction_coefficient => 0.1, :radius => 0.0, :side => 0.5)
generate_dataset(:box, H=0.40, N=35,
	opts=SolverOptions(btol=3e-4, rtol=3e-4),
	init_kwargs=init_kwargs,
	mech_kwargs=mech_kwargs,
	sleep_ratio=0.01,
	show_contact=false)

################################################################################
# Load Dataset
################################################################################
params0, trajs0, pairs0 = open_dataset(:box; N = 35, mech_kwargs...)

data0 = params0[:data]

################################################################################
# Optimization Objective: Evaluation & Gradient
################################################################################
# @benchmark clean_loss(:box, pairs0, data0, opts=SolverOptions(btol=3e-4, rtol=3e-4))
# @profiler clean_loss(:box, pairs0, data0, opts=SolverOptions(btol=3e-4, rtol=3e-4))
clean_loss(:box, pairs0, data0, opts=SolverOptions(btol=3e-4, rtol=3e-4))

[clean_loss(:box, pairs0, data0 + [i;zeros(55)], opts=SolverOptions(btol=3e-4, rtol=3e-4))
	for i in Vector(-0.10:0.01:0.1)]
[clean_loss(:box, pairs0, data0 + [0;i;zeros(54)], opts=SolverOptions(btol=3e-4, rtol=3e-4))
	for i in Vector(-0.10:0.01:0.1)]

plot(hcat([p[1][1:3] for p in pairs0]...)')
plot(hcat([p[1][4:6] for p in pairs0]...)')
plot(hcat([p[1][7:10] for p in pairs0]...)')
plot(hcat([p[1][11:13] for p in pairs0]...)')

################################################################################
# Optimization Algorithm: Quasi Newton:
# We learn a single coefficient of friction and a 8 contact locations [x,y,z] -> 25 params in total
################################################################################

include("../quasi_newton.jl")
function d2data(d)
	friction_coefficient = d[1]
	data = [friction_coefficient; 0;0;0; +d[2:4];
			friction_coefficient; 0;0;0; +d[5:7];
			friction_coefficient; 0;0;0; +d[8:10];
			friction_coefficient; 0;0;0; +d[11:13];
			friction_coefficient; 0;0;0; +d[14:16];
			friction_coefficient; 0;0;0; +d[17:19];
			friction_coefficient; 0;0;0; +d[20:22];
			friction_coefficient; 0;0;0; +d[23:25];
			]
	return data
end
∇d2data = ForwardDiff.jacobian(d -> d2data(d), zeros(25))

d0 = [0.40,
	+0.50, +0.50, -0.50,
	+0.50, -0.50, -0.50,
	-0.50, +0.50, -0.50,
	-0.50, -0.50, -0.50,
	+0.50, +0.50, +0.50,
	+0.50, -0.50, +0.50,
	-0.50, +0.50, +0.50,
	-0.50, -0.50, +0.50]
lower = [0.00,
	+0.05, +0.05, -1.00,
	+0.05, -1.00, -1.00,
	-1.00, +0.05, -1.00,
	-1.00, -1.00, -1.00,
	+0.05, +0.05, +0.05,
	+0.05, -1.00, +0.05,
	-1.00, +0.05, +0.05,
	-1.00, -1.00, +0.05]
upper = [0.80,
	+1.00, +1.00, -0.05,
	+1.00, -0.05, -0.05,
	-0.05, +1.00, -0.05,
	-0.05, -0.05, -0.05,
	+1.00, +1.00, +1.00,
	+1.00, -0.05, +1.00,
	-0.05, +1.00, +1.00,
	-0.05, -0.05, +1.00]

function f0(d; rot=0)
	return clean_loss(:box, pairs0, d2data(d), n_sample=35, rot=rot, opts=SolverOptions(btol=3e-4, rtol=3e-4))[1]
end

function fgH0(d; rot=0)
	f, g, H = clean_loss(:box, pairs0, d2data(d), n_sample=35, rot=rot, opts=SolverOptions(btol=3e-4, rtol=3e-4))
	return f, ∇d2data' * g, ∇d2data' * H * ∇d2data
end
dsol = quasi_newton_solve(f0, fgH0, d0, iter=50, gtol=1e-8, ftol=1e-6,
	lower=lower, upper=upper, reg=1e-9)

d0
f0(d0)
f0(dsol)
norm(data0 - d2data(dsol), Inf)

# We can learn the coefficient of friction and the side dimenson of the cube
# form 35*0.40 seconds of recording. We use the simulator to evaluate the loss
# and its gradients by differentiating through the simulator. With gradient
# information we can use L-BFGS

################################################################################
# Visualization
################################################################################

include("../quasi_newton.jl")
function d2data(d)
	data = [d[1]; 0;0;0; +d[2:4];
			d[5]; 0;0;0; +d[6:8];
			d[9]; 0;0;0; +d[10:12];
			d[13]; 0;0;0; +d[14:16];
			d[17]; 0;0;0; +d[18:20];
			d[21]; 0;0;0; +d[22:24];
			d[25]; 0;0;0; +d[26:28];
			d[29]; 0;0;0; +d[30:32];
			]
	return data
end
∇d2data = ForwardDiff.jacobian(d -> d2data(d), zeros(32))

d0 = [
	0.40, +0.50, +0.50, -0.50,
	0.40, +0.50, -0.50, -0.50,
	0.40, -0.50, +0.50, -0.50,
	0.40, -0.50, -0.50, -0.50,
	0.40, +0.50, +0.50, +0.50,
	0.40, +0.50, -0.50, +0.50,
	0.40, -0.50, +0.50, +0.50,
	0.40, -0.50, -0.50, +0.50]
lower = [
	0.00, +0.05, +0.05, -1.00,
	0.00, +0.05, -1.00, -1.00,
	0.00, -1.00, +0.05, -1.00,
	0.00, -1.00, -1.00, -1.00,
	0.00, +0.05, +0.05, +0.05,
	0.00, +0.05, -1.00, +0.05,
	0.00, -1.00, +0.05, +0.05,
	0.00, -1.00, -1.00, +0.05]
upper = [
	0.80, +1.00, +1.00, -0.05,
	0.80, +1.00, -0.05, -0.05,
	0.80, -0.05, +1.00, -0.05,
	0.80, -0.05, -0.05, -0.05,
	0.80, +1.00, +1.00, +1.00,
	0.80, +1.00, -0.05, +1.00,
	0.80, -0.05, +1.00, +1.00,
	0.80, -0.05, -0.05, +1.00]

function f0(d; rot=0)
	return clean_loss(:box, pairs0, d2data(d), n_sample=35, rot=rot, opts=SolverOptions(btol=3e-4, rtol=3e-4))[1]
end

function fgH0(d; rot=0)
	f, g, H = clean_loss(:box, pairs0, d2data(d), n_sample=35, rot=rot, opts=SolverOptions(btol=3e-4, rtol=3e-4))
	return f, ∇d2data' * g, ∇d2data' * H * ∇d2data
end

dsol, Dsol = quasi_newton_solve(f0, fgH0, d0, iter=50, gtol=1e-8, ftol=1e-6,
	lower=lower, upper=upper, reg=1e-9)

d0
f0(d0)
f0(dsol)
norm(data0 - d2data(dsol), Inf)

# We can learn the coefficient of friction and the side dimenson of the cube
# form 35*0.40 seconds of recording. We use the simulator to evaluate the loss
# and its gradients by differentiating through the simulator. With gradient
# information we can use L-BFGS

################################################################################
# Visualization
################################################################################

# Visualize Dataset with solution
datasol = ∇d2data * dsol
for i = 1:8
	datasol[(i-1)*7 + 4] = 0.05
end
set_simulator_data!(mech, datasol)
for traj in trajs0
	visualize(mech, traj, vis=vis, show_contact = true)
	sleep(0.5)
end

# Visualize Progress made during 'learning'
for (i,dsol) in enumerate(Dsol[1:6])
	datasol = ∇d2data * dsol
	for i = 1:8
		datasol[(i-1)*7 + 4] = 0.05
	end
	set_simulator_data!(mech, datasol)
	visualize(mech, trajs0[i], vis=vis, show_contact = true)
	sleep(1.5)
end


text:
macro import_layout_huge()
    :(
        using ComplexVisual:
            cv_destroy,
            CV_Layout, CV_Abstract2DLayout, CV_2DLayoutWrapper,
            CV_2DContainer, CV_2DLayout,
            CV_Framed2DLayout, CV_MinimalFramed2DLayout,
            CV_2DLayoutPosition, CV_2DLayoutCanvas, cv_create_context,
            cv_get_seen_boxes, cv_width, cv_height,
            cv_translate, cv_add_canvas!, cv_add_rectangle!, cv_add_padding!,
            cv_ensure_size!,
            cv_canvas_for_layout, cv_anchor, cv_global2local, cv_local2global,
            cv_pixel2math,
            CV_StateLayout, cv_get_state_counter,
            CV_SceneSetupChain, CV_2DScene
    )
end

import Base:show

# {{{ the layouts, composition-concept, macro for getter-methods

"""
A layout is able to position other objects (e.g. canvas, rectangles)
relative to already positioned objects.

Typically a layout "grows" (in the size) with new objects. Even additional
information may be stored in a layout, see `CV_2DLayoutWrapper`.

At some time a layout gets "framed", see `CV_Framed2DLayout`. At this time
such a layout has a `can_layout`  and `cc_can_layout` field.

At the ende a layout becomes a "scene", see `CV_2DScene`. At this time
there are some callback-functions for interaction available as fields.
"""
abstract type CV_Layout         end

show(io::IO, l::CV_Layout) = cv_show_impl(io, l)
show(io::IO, m::MIME{Symbol("text/plain")}, l::CV_Layout) =
    cv_show_impl(io, m, l)

"""
A 2D layout.

see also `CV_Layout`.
"""
abstract type CV_Abstract2DLayout   <: CV_Layout end

"""
That's the root for a composition "hierarchy" of immutable structs. There
are "getter"-methods to access the fields.

Why a composition "hierarchy"?

Sometimes some additional informations (often `CV_2DLayoutPosition`,
`CV_Math2DCanvas` or `CV_Std2DCanvas`) have to be stored with the layout.
In order to allow for *inferable* types (i.e. no `Vector{...}` is possible)
and to store such objects in immutable structs the composition concept is used.
If some additional informations has to be stored a "wrapper" is used and
the "original" layout is saved inside in a field called `parent_layout`.

Field access is done with (type-inferable) "getter"-methods.
"""
abstract type CV_2DLayoutWrapper <: CV_Abstract2DLayout end

"""
A Layout which supports the fields `can_layout` and `cc_can_layout`.
At this time the size of the layout cannot grow anymore (because the canvas
is already allocated).
"""
abstract type CV_Framed2DLayout <: CV_2DLayoutWrapper
    # can_layout      :: CV_2DLayoutCanvas
    # cc_can_layout
end

"""
internal macro to (semi-)automatically create the "getter"-methods
for 2D layout (wrappers).

For the owner `getproperty` is skipped. `getfield` gets called directly.
"""
macro layout_composition_getter(field, owner_type)
    func_sym = Symbol("cv_get_", field)
    func_quote = Meta.quot(func_sym)
    field_sym = Meta.quot(Symbol(field))
    return esc(quote
        if !isdefined(ComplexVisual, $func_quote)
            $func_sym(al::CV_2DLayoutWrapper) = $func_sym(al.parent_layout)
        end
        $func_sym(layout::$owner_type) = getfield(layout, $field_sym)
    end)
end

cv_destroy(l::CV_2DLayoutWrapper) = cv_destroy(l.parent_layout)

# }}}

struct CV_2DLayout  <: CV_Abstract2DLayout         # {{{
    seen_boxes      :: CV_RectangleStore{Int32}

    function CV_2DLayout()
        return new(CV_RectangleStore(Int32))
    end
end
cv_destroy(layout::CV_2DLayout) = nothing
@layout_composition_getter(seen_boxes, CV_2DLayout)
# }}}

struct CV_MinimalFramed2DLayout{parentT,
            canT <: CV_2DCanvas, cccT} <: CV_Framed2DLayout # {{{
    parent_layout   :: parentT
    can_layout      :: canT
    cc_can_layout   :: cccT
end

@layout_composition_getter(can_layout,    CV_MinimalFramed2DLayout)
@layout_composition_getter(cc_can_layout, CV_MinimalFramed2DLayout)

function cv_destroy(layout::CV_MinimalFramed2DLayout)
    cv_destroy(layout.cc_can_layout)
    cv_destroy(layout.can_layout)
    cv_destroy(layout.parent_layout)
    return nothing
end
# }}}

"""
A container with a `bounding_box` and `user_box` (pixel-)coordinates.

For the `bounding_box` the bottom left corner is at `(0, 0)`.

The `user_box` coordinates may be shifted but must not be scaled.
"""
abstract type  CV_2DContainer   <: CV_2DCanvas end

"""
`CV_2DCanvas` with size and trafo adapted to `CV_2DLayout`.
"""
struct CV_2DLayoutCanvas{afT} <: CV_2DContainer  # {{{
    surface      :: Cairo.CairoSurfaceImage{UInt32}
    pixel_width  :: Int32
    pixel_height :: Int32
    bounding_box :: CV_Rectangle{Int32} # zero-based
    user_box     :: CV_Rectangle{Int32} # user-coordinates (result of layout)
                                        # typically nonzero-based
    anchor_func  :: afT
    function CV_2DLayoutCanvas(user_box::CV_Rectangle{Int32},
            anchor_func=(can, name) -> cv_anchor(can.bounding_box, name))
        width, height = cv_width(user_box), cv_height(user_box)
        surface = cv_create_cairo_image_surface(width, height)
        self = new{typeof(anchor_func)}(
            surface, width, height,
            CV_Rectangle(height, Int32(0), Int32(0), width),
            user_box, anchor_func)
        return self
    end
end

function cv_anchor(can::CV_2DLayoutCanvas, anchor_name::Symbol)
    return can.anchor_func(can, anchor_name) :: Tuple{Int32, Int32}
end

function cv_create_context(canvas::CV_2DContainer; prepare::Bool=true,
        fill_with::Union{CV_ContextStyle, Nothing}=cv_white)
    con = CV_2DCanvasContext(canvas)
    if prepare
        ctx = con.ctx
        reset_transform(ctx)

        if fill_with isa CV_ContextStyle
            set_operator(ctx, Cairo.OPERATOR_SOURCE)
            cv_prepare(con, fill_with)
            rectangle(ctx, 0, 0, canvas.pixel_width, canvas.pixel_height)
            fill(ctx)
            set_operator(ctx, Cairo.OPERATOR_OVER)
        end

        ubox = canvas.user_box
        translate(ctx, -ubox.left, -ubox.bottom)
    end
    return con
end

"""
"finalize" layout construction by building a canvas, on which all
added objects are visible.
"""
function cv_canvas_for_layout(layout::CV_Abstract2DLayout,
        anchor_func=(can, name) -> cv_anchor(can.bounding_box, name))
    bb = cv_get_seen_boxes(layout).bounding_box
    if bb.empty
        cv_error("Bounding box for Layout is empty.")
    end
    canvas = CV_2DLayoutCanvas(bb, anchor_func)
    return canvas
end
# }}}

"""
Destination Position/Size (a rectangle) in the layout. It's
for an optional canvas or for an optional "drawing callback".

This struct is callable: If called (with a CV_2DCanvasContext) then
(a) if a canvas was given, then this canvas is copied to its target
    position (rectangle) with style "style".
(b) if a drawing_callback was given, then the style is used and the
    drawing-callback is called.
"""
struct CV_2DLayoutPosition{canT<:Union{CV_2DCanvas, Nothing},
                           dcbT, styleT<:CV_ContextStyle}  # {{{
    rectangle   :: CV_Rectangle{Int32}
    canvas      :: canT
    drawing_cb  :: dcbT
    style       :: styleT
end

show(io::IO, p::CV_2DLayoutPosition) = cv_show_impl(io, p)
show(io::IO, m::MIME{Symbol("text/plain")}, p::CV_2DLayoutPosition) =
    cv_show_impl(io, m, p)

"""
if canvas was given:
copy canvas to its target position (given in rectangle) with style "style".

otherweise call: drawing_cb.

See `CV_2DLayoutPosition`
"""
function (cl::CV_2DLayoutPosition)(cc::CV_2DCanvasContext)
    ctx = cc.ctx
    if !(cl.canvas isa Nothing)
        cv_prepare(cc, cl.style)
        rect = cl.rectangle
        set_source_surface(ctx, cl.canvas.surface, rect.left, rect.bottom)
        rectangle(ctx, rect.left, rect.bottom, cv_width(rect), cv_height(rect))
        fill(ctx)
    end
    if !(cl.drawing_cb isa Nothing)
        cv_prepare(cc, cl.style)
        cl.drawing_cb(cc, cl)
    end
    return nothing
end

# }}}

"""
use 2DLayoutPosition to transform a global pixel position `(gx, gy)` to
local/relative pixels w.r.t. the positon's coordinates.
"""
function cv_global2local(canvas::CV_2DContainer,
                        cl::CV_2DLayoutPosition, gx::Integer, gy::Integer)
    ubox = canvas.user_box
    rect = cl.rectangle
    ux, uy = Int32(gx) + ubox.left, Int32(gy) + ubox.bottom
    return (ux - rect.left, uy - rect.bottom)
end

"""
use 2DLayoutPosition to transform a local pixel position `(lx, ly)`
w.r.t. the positon's coordinates to global pixels.

This is the opposite of `cv_global2local`.
"""
function cv_local2global(canvas::CV_2DContainer,
                         cl::CV_2DLayoutPosition, lx::Integer, ly::Integer)
    ubox = canvas.user_box
    rect = cl.rectangle
    ux, uy = Int32(lx) + rect.left, Int32(ly) + rect.bottom
    return (ux - ubox.left, uy - ubox.bottom)
end

"""
use 2DLayoutPosition (for a `CV_Math2DCanvas`) to convert a global
pixel position `(gx, gy)` coordinates in math units.
"""
function cv_pixel2math(canvas::CV_2DContainer,
                       cl::CV_2DLayoutPosition{canT,styleT},
                       gx::Integer, gy::Integer) where {styleT,
                                                canT<:CV_Math2DCanvas}
    lx, ly = cv_global2local(canvas, cl, gx, gy)
    return cv_pixel2math(cl.canvas, lx, ly)
end

"""
Return anchor-point for a `CV_2DLayoutPosition`.

If there is position belongs to a canvas then its `cv_anchor` function
is used. Otherwise the anchor function of rectangle (of the position)
is used.
"""
function cv_anchor(lres::CV_2DLayoutPosition, anchor_name::Symbol)
    rect = lres.rectangle
    if lres.canvas isa Nothing
        # Fallback to anchors of the rectangle
        return cv_anchor(rect, anchor_name)
    else
        # support all anchors of the canvas
        inner = cv_anchor(lres.canvas, anchor_name)
        return (rect.left + inner[1], rect.bottom + inner[2])
    end
end

"""
width and/or height of `CV_2DLayoutPosition` is given by its rectangle.
"""
cv_width(lpos::CV_2DLayoutPosition) = cv_width(lpos.rectangle)
cv_height(lpos::CV_2DLayoutPosition) = cv_height(lpos.rectangle)

"""
build a new anchor using the x-coordinate of the `anchor1` of object `obj1`
and the y-coordinate of the `anchor2` of `obj2`.
"""
function cv_anchor(obj1, anchor1::Symbol, obj2, anchor2::Symbol)
    return (cv_anchor(obj1, anchor1)[1], cv_anchor(obj2, anchor2)[2])
end

"""
Move a position (typically a anchor) by `deltax` and `deltay`.
"""
function cv_translate(pos::Tuple{N, N}, deltax, deltay) where {N<:Number}
    return tuple(N(pos[1] + deltax), N(pos[2] + deltay))
end

const cv_add_canvas_default_style = cv_op_over

"""
add a canvas to the layout scheme. `where` describes the position
to place the canvas and `anchor` describes which part of the canvas
is located at the position `where`.

With `style` one can set the mode how the canvas is copied to the
target surface.
"""
function cv_add_canvas!(layout::CV_Abstract2DLayout,
                        canvas::CV_2DCanvas,
                        anchor::Tuple{Integer, Integer},
                        where::Tuple{Integer, Integer};
                        style::CV_ContextStyle=cv_add_canvas_default_style)
    where = (where[1] - anchor[1], where[2] - anchor[2])
    rect = cv_rect_blwh(Int32, where[2], where[1],
                               canvas.pixel_width, canvas.pixel_height)
    cv_add_rectangle!(cv_get_seen_boxes(layout), rect)
    return CV_2DLayoutPosition(rect, canvas, nothing, style)
end

const cv_add_rectangle_default_style = cv_op_over

"""
add a rectangle of given width and height to the layout scheme. `where`
describes the position to place the rectangle and `anchor` describes which
part of the rectangle is located at the position `where`.
"""
function cv_add_rectangle!(layout::CV_Abstract2DLayout,
                           width::Integer, height::Integer,
                           anchor::Tuple{Integer, Integer},
                           where::Tuple{Integer, Integer},
                           drawing_cb=nothing,
                           style=cv_add_rectangle_default_style)
    where = (where[1] - anchor[1], where[2] - anchor[2])
    rect = cv_rect_blwh(Int32, where[2], where[1], width, height)
    return cv_add_rectangle!(layout, rect, drawing_cb, style)
end

"""
add a rectangle (absolute coordinates are given in the rectangle) to the
layout scheme. The result is a `CV_2DLayoutPosition` which (when called),
calls the given drawing_callback. The drawing_callback
has the form:  `drawing_cb(cc::CV_2DCanvasContext, pos::CV_2DLayoutPosition)`
"""
function cv_add_rectangle!(layout::CV_Abstract2DLayout,
                           rect::CV_Rectangle{Int32}, drawing_cb=nothing,
                           style=cv_op_over)
    cv_add_rectangle!(cv_get_seen_boxes(layout), rect)
    return CV_2DLayoutPosition(rect, nothing, drawing_cb, style)
end

"""
add padding (with the help of cv_add_rectangle!) to the current
target layout area.
"""
function cv_add_padding!(layout::CV_Abstract2DLayout,
                         north::Integer, east::Integer=north,
                         south::Integer=north, west::Integer=east) # {{{
    o = one(Int32)
    bb = cv_get_seen_boxes(layout).bounding_box
    style = cv_op_source
    if north > 0
        cv_add_rectangle!(layout,
            CV_Rectangle(bb.bottom, bb.left,
                         bb.bottom - Int32(north), bb.left + o), nothing, style)
    end
    if south > 0
        cv_add_rectangle!(layout,
            CV_Rectangle(bb.top + Int32(south), bb.left, bb.top, bb.left + o),
            nothing, style)
    end
    if east > 0
        cv_add_rectangle!(layout,
            CV_Rectangle(bb.top, bb.right, bb.top - o, bb.right + Int32(east)),
            nothing, style)
    end
    if west > 0
        cv_add_rectangle!(layout,
            CV_Rectangle(bb.top, bb.left - Int32(west), bb.top - o, bb.left),
            nothing, style)
    end
    return nothing
end # }}}

function cv_ensure_size!(layout::CV_Abstract2DLayout,
                         minwidth::Integer, minheight::Integer) # {{{
    minwidth, minheight = Int32(minwidth), Int32(minheight)

    bb = cv_get_seen_boxes(layout).bounding_box
    missing = minwidth - cv_width(bb)
    if missing > 0
        half = cv_half(missing)
        rest = missing - half

        cv_add_rectangle!(layout,
            CV_Rectangle(bb.top, bb.left - half, bb.bottom, bb.right + rest))
    end

    bb = cv_get_seen_boxes(layout).bounding_box
    missing = minheight - cv_height(bb)
    if missing >0
        half = cv_half(missing)
        rest = missing - half
        cv_add_rectangle!(layout,
            CV_Rectangle(bb.top + half, bb.left, bb.bottom - rest, bb.right))
    end
    return nothing
end # }}}

"""
Layout with degree of freedoms: state_counter
"""
struct CV_StateLayout{parentT, maxV} <: CV_2DLayoutWrapper # {{{
    parent_layout   :: parentT
    state_counter   :: CV_CyclicValue{maxV}
end

@layout_composition_getter(state_counter, CV_StateLayout)

# }}}

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
include("organ.jl")

include("leaf.jl")
include("sheath.jl")
include("nodalunit.jl")
include("root.jl")


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

"""
    Generator(f, iter)

Given a function `f` and an iterator `iter`, construct an iterator that yields
the values of `f` applied to the elements of `iter`.
The syntax `f(x) for x in iter [if cond(x)::Bool]` is syntax for constructing an instance of this
type. The `[if cond(x)::Bool]` expression is optional and acts as a "guard", effectively
filtering out values where the condition is false.

```jldoctest
julia> g = (abs2(x) for x in 1:5 if x != 3);

julia> for x in g
           println(x)
       end
1
4
16
25

julia> collect(g)
4-element Array{Int64,1}:
  1
  4
 16
 25
```
"""
struct Generator{I,F}
    f::F
    iter::I
end

Generator(f, I1, I2, Is...) = Generator(a->f(a...), zip(I1, I2, Is...))

Generator(::Type{T}, iter::I) where {T,I} = Generator{I,Type{T}}(T, iter)

start(g::Generator) = (@_inline_meta; start(g.iter))
done(g::Generator, s) = (@_inline_meta; done(g.iter, s))
function next(g::Generator, s)
    @_inline_meta
    v, s2 = next(g.iter, s)
    g.f(v), s2
end

## iterator traits

abstract type IteratorSize end
struct SizeUnknown <: IteratorSize end
struct HasLength <: IteratorSize end
struct HasShape <: IteratorSize end
struct IsInfinite <: IteratorSize end

"""
    iteratorsize(itertype::Type) -> IteratorSize

Given the type of an iterator, returns one of the following values:

* `SizeUnknown()` if the length (number of elements) cannot be determined in advance.
* `HasLength()` if there is a fixed, finite length.
* `HasShape()` if there is a known length plus a notion of multidimensional shape (as for an array).
   In this case the [`size`](@ref) function is valid for the iterator.
* `IsInfinite()` if the iterator yields values forever.

The default value (for iterators that do not define this function) is `HasLength()`.
This means that most iterators are assumed to implement [`length`](@ref).

This trait is generally used to select between algorithms that pre-allocate space for their
result, and algorithms that resize their result incrementally.

```jldoctest
julia> Base.iteratorsize(1:5)
Base.HasShape()

julia> Base.iteratorsize((2,3))
Base.HasLength()
```
"""
iteratorsize(x) = iteratorsize(typeof(x))
iteratorsize(::Type) = HasLength()  # HasLength is the default

abstract type IteratorEltype end
struct EltypeUnknown <: IteratorEltype end
struct HasEltype <: IteratorEltype end

"""
    iteratoreltype(itertype::Type) -> IteratorEltype

Given the type of an iterator, returns one of the following values:

* `EltypeUnknown()` if the type of elements yielded by the iterator is not known in advance.
* `HasEltype()` if the element type is known, and [`eltype`](@ref) would return a meaningful value.

`HasEltype()` is the default, since iterators are assumed to implement [`eltype`](@ref).

This trait is generally used to select between algorithms that pre-allocate a specific
type of result, and algorithms that pick a result type based on the types of yielded
values.

```jldoctest
julia> Base.iteratoreltype(1:5)
Base.HasEltype()
```
"""
iteratoreltype(x) = iteratoreltype(typeof(x))
iteratoreltype(::Type) = HasEltype()  # HasEltype is the default

iteratorsize(::Type{<:AbstractArray}) = HasShape()
iteratorsize(::Type{Generator{I,F}}) where {I,F} = iteratorsize(I)
length(g::Generator) = length(g.iter)
size(g::Generator) = size(g.iter)
indices(g::Generator) = indices(g.iter)
ndims(g::Generator) = ndims(g.iter)

iteratoreltype(::Type{Generator{I,T}}) where {I,T} = EltypeUnknown()

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia