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

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

 Use to associate an AWS CodeCommit repository or a repostory managed by AWS CodeStar
Connections with Amazon CodeGuru Reviewer. When you associate a repository, CodeGuru
Reviewer reviews source code changes in the repository's pull requests and provides
automatic recommendations. You can view recommendations using the CodeGuru Reviewer
console. For more information, see Recommendations in Amazon CodeGuru Reviewer in the
Amazon CodeGuru Reviewer User Guide.  If you associate a CodeCommit or S3 repository, it
must be in the same AWS Region and AWS account where its CodeGuru Reviewer code reviews are
configured. Bitbucket and GitHub Enterprise Server repositories are managed by AWS CodeStar
Connections to connect to CodeGuru Reviewer. For more information, see Associate a
repository in the Amazon CodeGuru Reviewer User Guide.    You cannot use the CodeGuru
Reviewer SDK or the AWS CLI to associate a GitHub repository with Amazon CodeGuru Reviewer.
To associate a GitHub repository, use the console. For more information, see Getting
started with CodeGuru Reviewer in the CodeGuru Reviewer User Guide.

# Arguments
- `repository`: The repository to associate.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientRequestToken"`: Amazon CodeGuru Reviewer uses this value to prevent the
  accidental creation of duplicate repository associations if there are failures and retries.
- `"KMSKeyDetails"`: A KMSKeyDetails object that contains:   The encryption option for this
  repository association. It is either owned by AWS Key Management Service (KMS)
  (AWS_OWNED_CMK) or customer managed (CUSTOMER_MANAGED_CMK).   The ID of the AWS KMS key
  that is associated with this respository association.
- `"Tags"`:  An array of key-value pairs used to tag an associated repository. A tag is a
  custom attribute label with two parts:    A tag key (for example, CostCenter, Environment,
  Project, or Secret). Tag keys are case sensitive.   An optional field known as a tag value
  (for example, 111122223333, Production, or a team name). Omitting the tag value is the same
  as using an empty string. Like tag keys, tag values are case sensitive.
"""
associate_repository(Repository; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("POST", "/associations", Dict{String, Any}("Repository"=>Repository, "ClientRequestToken"=>string(uuid4())); aws_config=aws_config)
associate_repository(Repository, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("POST", "/associations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Repository"=>Repository, "ClientRequestToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    create_code_review(name, repository_association_arn, type)
    create_code_review(name, repository_association_arn, type, params::Dict{String,<:Any})

 Use to create a code review with a  CodeReviewType  of RepositoryAnalysis. This type of
code review analyzes all code under a specified branch in an associated repository.
PullRequest code reviews are automatically triggered by a pull request.

# Arguments
- `name`:  The name of the code review. The name of each code review in your AWS account
  must be unique.
- `repository_association_arn`:  The Amazon Resource Name (ARN) of the
  RepositoryAssociation  object. You can retrieve this ARN by calling
  ListRepositoryAssociations .   A code review can only be created on an associated
  repository. This is the ARN of the associated repository.
- `type`:  The type of code review to create. This is specified using a  CodeReviewType
  object. You can create a code review only of type RepositoryAnalysis.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientRequestToken"`:  Amazon CodeGuru Reviewer uses this value to prevent the
  accidental creation of duplicate code reviews if there are failures and retries.
"""
create_code_review(Name, RepositoryAssociationArn, Type; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("POST", "/codereviews", Dict{String, Any}("Name"=>Name, "RepositoryAssociationArn"=>RepositoryAssociationArn, "Type"=>Type, "ClientRequestToken"=>string(uuid4())); aws_config=aws_config)
create_code_review(Name, RepositoryAssociationArn, Type, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("POST", "/codereviews", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name, "RepositoryAssociationArn"=>RepositoryAssociationArn, "Type"=>Type, "ClientRequestToken"=>string(uuid4())), params)); aws_config=aws_config)

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

 Returns the metadata associated with the code review along with its status.

# Arguments
- `code_review_arn`: The Amazon Resource Name (ARN) of the  CodeReview  object.

"""
describe_code_review(CodeReviewArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/codereviews/$(CodeReviewArn)"; aws_config=aws_config)
describe_code_review(CodeReviewArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/codereviews/$(CodeReviewArn)", params; aws_config=aws_config)

"""
    describe_recommendation_feedback(code_review_arn, recommendation_id)
    describe_recommendation_feedback(code_review_arn, recommendation_id, params::Dict{String,<:Any})

 Describes the customer feedback for a CodeGuru Reviewer recommendation.

# Arguments
- `code_review_arn`: The Amazon Resource Name (ARN) of the  CodeReview  object.
- `recommendation_id`:  The recommendation ID that can be used to track the provided
  recommendations and then to collect the feedback.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"UserId"`:  Optional parameter to describe the feedback for a given user. If this is not
  supplied, it defaults to the user making the request.   The UserId is an IAM principal that
  can be specified as an AWS account ID or an Amazon Resource Name (ARN). For more
  information, see  Specifying a Principal in the AWS Identity and Access Management User
  Guide.
"""
describe_recommendation_feedback(CodeReviewArn, RecommendationId; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/feedback/$(CodeReviewArn)", Dict{String, Any}("RecommendationId"=>RecommendationId); aws_config=aws_config)
describe_recommendation_feedback(CodeReviewArn, RecommendationId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/feedback/$(CodeReviewArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RecommendationId"=>RecommendationId), params)); aws_config=aws_config)

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

 Returns a  RepositoryAssociation  object that contains information about the requested
repository association.

# Arguments
- `association_arn`:  The Amazon Resource Name (ARN) of the  RepositoryAssociation  object.
  You can retrieve this ARN by calling  ListRepositoryAssociations .

"""
describe_repository_association(AssociationArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/associations/$(AssociationArn)"; aws_config=aws_config)
describe_repository_association(AssociationArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/associations/$(AssociationArn)", params; aws_config=aws_config)

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

Removes the association between Amazon CodeGuru Reviewer and a repository.

# Arguments
- `association_arn`:  The Amazon Resource Name (ARN) of the  RepositoryAssociation  object.
  You can retrieve this ARN by calling  ListRepositoryAssociations .

"""
disassociate_repository(AssociationArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("DELETE", "/associations/$(AssociationArn)"; aws_config=aws_config)
disassociate_repository(AssociationArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("DELETE", "/associations/$(AssociationArn)", params; aws_config=aws_config)

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

 Lists all the code reviews that the customer has created in the past 90 days.

# Arguments
- `type`:  The type of code reviews to list in the response.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:  The maximum number of results that are returned per call. The default is
  100.
- `"NextToken"`:  If nextToken is returned, there are more results available. The value of
  nextToken is a unique pagination token for each page. Make the call again using the
  returned token to retrieve the next page. Keep all other arguments unchanged.
- `"ProviderTypes"`:  List of provider types for filtering that needs to be applied before
  displaying the result. For example, providerTypes=[GitHub] lists code reviews from GitHub.
- `"RepositoryNames"`:  List of repository names for filtering that needs to be applied
  before displaying the result.
- `"States"`:  List of states for filtering that needs to be applied before displaying the
  result. For example, states=[Pending] lists code reviews in the Pending state.  The valid
  code review states are:    Completed: The code review is complete.     Pending: The code
  review started and has not completed or failed.     Failed: The code review failed.
  Deleting: The code review is being deleted.
"""
list_code_reviews(Type; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/codereviews", Dict{String, Any}("Type"=>Type); aws_config=aws_config)
list_code_reviews(Type, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/codereviews", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Type"=>Type), params)); aws_config=aws_config)

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

 Returns a list of  RecommendationFeedbackSummary  objects that contain customer
recommendation feedback for all CodeGuru Reviewer users.

# Arguments
- `code_review_arn`: The Amazon Resource Name (ARN) of the  CodeReview  object.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:  The maximum number of results that are returned per call. The default is
  100.
- `"NextToken"`:  If nextToken is returned, there are more results available. The value of
  nextToken is a unique pagination token for each page. Make the call again using the
  returned token to retrieve the next page. Keep all other arguments unchanged.
- `"RecommendationIds"`:  Used to query the recommendation feedback for a given
  recommendation.
- `"UserIds"`:  An AWS user's account ID or Amazon Resource Name (ARN). Use this ID to
  query the recommendation feedback for a code review from that user.   The UserId is an IAM
  principal that can be specified as an AWS account ID or an Amazon Resource Name (ARN). For
  more information, see  Specifying a Principal in the AWS Identity and Access Management
  User Guide.
"""
list_recommendation_feedback(CodeReviewArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/feedback/$(CodeReviewArn)/RecommendationFeedback"; aws_config=aws_config)
list_recommendation_feedback(CodeReviewArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/feedback/$(CodeReviewArn)/RecommendationFeedback", params; aws_config=aws_config)

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

 Returns the list of all recommendations for a completed code review.

# Arguments
- `code_review_arn`: The Amazon Resource Name (ARN) of the  CodeReview  object.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:  The maximum number of results that are returned per call. The default is
  100.
- `"NextToken"`:  Pagination token.
"""
list_recommendations(CodeReviewArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/codereviews/$(CodeReviewArn)/Recommendations"; aws_config=aws_config)
list_recommendations(CodeReviewArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/codereviews/$(CodeReviewArn)/Recommendations", params; aws_config=aws_config)

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

 Returns a list of  RepositoryAssociationSummary  objects that contain summary information
about a repository association. You can filter the returned list by  ProviderType ,  Name ,
 State , and  Owner .

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of repository association results returned by
  ListRepositoryAssociations in paginated output. When this parameter is used,
  ListRepositoryAssociations only returns maxResults results in a single page with a
  nextToken response element. The remaining results of the initial request can be seen by
  sending another ListRepositoryAssociations request with the returned nextToken value. This
  value can be between 1 and 100. If this parameter is not used, ListRepositoryAssociations
  returns up to 100 results and a nextToken value if applicable.
- `"Name"`: List of repository names to use as a filter.
- `"NextToken"`: The nextToken value returned from a previous paginated
  ListRepositoryAssociations request where maxResults was used and the results exceeded the
  value of that parameter. Pagination continues from the end of the previous results that
  returned the nextToken value.   Treat this token as an opaque identifier that is only used
  to retrieve the next items in a list and not for other programmatic purposes.
- `"Owner"`: List of owners to use as a filter. For AWS CodeCommit, it is the name of the
  CodeCommit account that was used to associate the repository. For other repository source
  providers, such as Bitbucket and GitHub Enterprise Server, this is name of the account that
  was used to associate the repository.
- `"ProviderType"`: List of provider types to use as a filter.
- `"State"`: List of repository association states to use as a filter. The valid repository
  association states are:    Associated: The repository association is complete.
  Associating: CodeGuru Reviewer is:     Setting up pull request notifications. This is
  required for pull requests to trigger a CodeGuru Reviewer review.    If your repository
  ProviderType is GitHub, GitHub Enterprise Server, or Bitbucket, CodeGuru Reviewer creates
  webhooks in your repository to trigger CodeGuru Reviewer reviews. If you delete these
  webhooks, reviews of code in your repository cannot be triggered.      Setting up source
  code access. This is required for CodeGuru Reviewer to securely clone code in your
  repository.       Failed: The repository failed to associate or disassociate.
  Disassociating: CodeGuru Reviewer is removing the repository's pull request notifications
  and source code access.     Disassociated: CodeGuru Reviewer successfully disassociated the
  repository. You can create a new association with this repository if you want to review
  source code in it later. You can control access to code reviews created in an associated
  repository with tags after it has been disassociated. For more information, see Using tags
  to control access to associated repositories in the Amazon CodeGuru Reviewer User Guide.
"""
list_repository_associations(; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/associations"; aws_config=aws_config)
list_repository_associations(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/associations", params; aws_config=aws_config)

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

Returns the list of tags associated with an associated repository resource.

# Arguments
- `resource_arn`:  The Amazon Resource Name (ARN) of the  RepositoryAssociation  object.
  You can retrieve this ARN by calling  ListRepositoryAssociations .

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    put_recommendation_feedback(code_review_arn, reactions, recommendation_id)
    put_recommendation_feedback(code_review_arn, reactions, recommendation_id, params::Dict{String,<:Any})

 Stores customer feedback for a CodeGuru Reviewer recommendation. When this API is called
again with different reactions the previous feedback is overwritten.

# Arguments
- `code_review_arn`: The Amazon Resource Name (ARN) of the  CodeReview  object.
- `reactions`:  List for storing reactions. Reactions are utf-8 text code for emojis. If
  you send an empty list it clears all your feedback.
- `recommendation_id`:  The recommendation ID that can be used to track the provided
  recommendations and then to collect the feedback.

"""
put_recommendation_feedback(CodeReviewArn, Reactions, RecommendationId; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("PUT", "/feedback", Dict{String, Any}("CodeReviewArn"=>CodeReviewArn, "Reactions"=>Reactions, "RecommendationId"=>RecommendationId); aws_config=aws_config)
put_recommendation_feedback(CodeReviewArn, Reactions, RecommendationId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("PUT", "/feedback", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CodeReviewArn"=>CodeReviewArn, "Reactions"=>Reactions, "RecommendationId"=>RecommendationId), params)); aws_config=aws_config)

"""
    tag_resource(tags, resource_arn)
    tag_resource(tags, resource_arn, params::Dict{String,<:Any})

Adds one or more tags to an associated repository.

# Arguments
- `tags`:  An array of key-value pairs used to tag an associated repository. A tag is a
  custom attribute label with two parts:    A tag key (for example, CostCenter, Environment,
  Project, or Secret). Tag keys are case sensitive.   An optional field known as a tag value
  (for example, 111122223333, Production, or a team name). Omitting the tag value is the same
  as using an empty string. Like tag keys, tag values are case sensitive.
- `resource_arn`:  The Amazon Resource Name (ARN) of the  RepositoryAssociation  object.
  You can retrieve this ARN by calling  ListRepositoryAssociations .

"""
tag_resource(Tags, resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("POST", "/tags/$(resourceArn)", Dict{String, Any}("Tags"=>Tags); aws_config=aws_config)
tag_resource(Tags, resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Tags"=>Tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes a tag from an associated repository.

# Arguments
- `resource_arn`:  The Amazon Resource Name (ARN) of the  RepositoryAssociation  object.
  You can retrieve this ARN by calling  ListRepositoryAssociations .
- `tag_keys`: A list of the keys for each tag you want to remove from an associated
  repository.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codeguru_reviewer("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)


text:
# Based upon alg 1 in
# A unified sparse optimization framework to learn parsimonious physics-informed models from data
# by K Champion et. al.

mutable struct SR3{U,T} <: AbstractOptimizer
    λ::U
    ν::U
    R::T
end

"""
    SR3(λ, ν, R)
    SR3(λ = 1e-1, ν = 1.0)

`SR3` is an optimizer framework introduced [by Zheng et. al., 2018](https://ieeexplore.ieee.org/document/8573778) and used within
[Champion et. al., 2019](https://arxiv.org/abs/1906.10612). `SR3` contains a sparsification parameter `λ`, a relaxation `ν`,
and a corresponding penalty function `R`, which should be taken from [ProximalOperators.jl](https://github.com/kul-forbes/ProximalOperators.jl).

# Examples
```julia
opt = SR3()
opt = SR3(1e-2)
opt = SR3(1e-3, 1.0)
```
"""
function SR3(λ = 1e-1, ν = 1.0)
    R = NormL1
    return SR3(λ, ν, R)
end

function set_threshold!(opt::SR3, threshold)
    opt.λ = threshold^2*opt.ν /2
    return
end

get_threshold(opt::SR3) = sqrt(2*opt.λ/opt.ν)

init(o::SR3, A::AbstractArray, Y::AbstractArray) =  A \ Y
init!(X::AbstractArray, o::SR3, A::AbstractArray, Y::AbstractArray) =  ldiv!(X, qr(A, Val(true)), Y)

function fit!(X::AbstractArray, A::AbstractArray, Y::AbstractArray, opt::SR3; maxiter::Int64 = 1, convergence_error::T = eps()) where T <: Real
    f = opt.R(get_threshold(opt))

    n, m = size(A)
    W = copy(X)

    # Init matrices
    P = inv(A'*A+I(m)/(opt.ν))
    X̂ = A'*Y

    w_i = similar(W)
    w_i .= W
    iters = 0

    for i in 1:maxiter
        iters += 1
        # Solve ridge regression
        X .= P*(X̂.+W./(opt.ν))
        # Add proximal iteration
        prox!(W, f, X, opt.ν*opt.λ)

        if norm(w_i - W, 2)/opt.ν < convergence_error
            break
        else
            w_i .= W
        end

    end

    X[abs.(X) .< get_threshold(opt)] .= zero(eltype(X))
    return iters
end

ext: jl


text: println("z")

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Sockets

# IP addresses:
addr = ip"185.43.124.6"    # this uses the @ip_str macro
typeof(addr)               # IPv4

connect("julialang.org", 80)
# TCPSocket(Base.Libc.WindowsRawSocket(0x00000000000003ec) open, 0 bytes waiting)
getaddrinfo("julialang.org")
# ip"151.101.130.49"

# Installing URIs
# (@v1.8) pkg> add URIs, HTTP
using URIs, HTTP

# URIs:
req = HTTP.Request(
"GET",
"https://www.google.com/search?q=mammoth",
["Content-Type" => "text/plain"],
"Hi there!"
)

uri = URI(HTTP.uri(req))  # URI("https://www.google.com/search?q=mammoth")
uri.scheme # "https"
uri.host   # "www.google.com"
uri.path   # "/search"
uri.query  # "q=mammoth"
queryparams(uri) # Dict{String, String} with 1 entry: "q" => "mammoth"


text:
module TestProductOfSets

using Test

import MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

MOI.Utilities.@mix_of_scalar_sets(
    _ScalarSets,
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T},
)

function test_scalar_set_index()
    sets = _ScalarSets{Float64}()
    @test MOI.Utilities.set_index(sets, MOI.EqualTo{Float64}) == 1
    @test MOI.Utilities.set_index(sets, MOI.GreaterThan{Float64}) == 2
    @test MOI.Utilities.set_index(sets, MOI.Interval{Float64}) == 4
    @test MOI.Utilities.set_index(sets, MOI.LessThan{Float64}) == 3
    @test MOI.Utilities.set_index(sets, MOI.ZeroOne) === nothing
end

function test_scalar_set_types()
    sets = _ScalarSets{Float64}()
    @test MOI.Utilities.set_types(sets) == [
        MOI.EqualTo{Float64},
        MOI.GreaterThan{Float64},
        MOI.LessThan{Float64},
        MOI.Interval{Float64},
    ]
end

function test_scalar_basic()
    sets = _ScalarSets{Float64}()
    ci = MOI.ConstraintIndex{
        MOI.ScalarAffineFunction{Float64},
        MOI.EqualTo{Float64},
    }(
        12345,
    )
    @test !MOI.is_valid(sets, ci)
    i = MOI.Utilities.set_index(sets, MOI.EqualTo{Float64})
    ci_value = MOI.Utilities.add_set(sets, i)
    ci = MOI.ConstraintIndex{
        MOI.ScalarAffineFunction{Float64},
        MOI.EqualTo{Float64},
    }(
        ci_value,
    )
    @test MOI.is_valid(sets, ci)
    @test MOI.Utilities.rows(sets, ci) == ci.value
end

function test_scalar_dimension()
    sets = _ScalarSets{Float64}()
    @test MOI.dimension(sets) == 0
    MOI.Utilities.add_set(sets, 1)
    @test MOI.dimension(sets) == 1
    MOI.Utilities.add_set(sets, 1)
    @test MOI.dimension(sets) == 2
    MOI.Utilities.add_set(sets, 2)
    @test MOI.dimension(sets) == 3
end

function test_scalar_empty()
    sets = _ScalarSets{Float64}()
    @test MOI.is_empty(sets)
    MOI.Utilities.add_set(sets, 1)
    @test !MOI.is_empty(sets)
    MOI.empty!(sets)
    @test MOI.is_empty(sets)
end

function test_scalar_ConstraintTypesPresent()
    sets = _ScalarSets{Float64}()
    @test MOI.get(sets, MOI.ListOfConstraintTypesPresent()) == []
    MOI.Utilities.add_set(sets, 1)
    @test MOI.get(sets, MOI.ListOfConstraintTypesPresent()) ==
          [(MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})]
end

function test_scalar_NumberOfConstraints()
    sets = _ScalarSets{Float64}()
    MOI.Utilities.add_set(sets, 1)
    MOI.Utilities.add_set(sets, 1)
    MOI.Utilities.add_set(sets, 2)
    MOI.Utilities.final_touch(sets)
    @test MOI.get(
        sets,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.ZeroOne}(),
    ) == 0
    for (x, S) in zip([2, 1, 0, 0], MOI.Utilities.set_types(sets))
        @test MOI.get(
            sets,
            MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},S}(),
        ) == x
    end
    return
end

function test_scalar_ListOfConstraintIndices()
    sets = _ScalarSets{Float64}()
    MOI.Utilities.add_set(sets, 2)
    MOI.Utilities.add_set(sets, 4)
    MOI.Utilities.add_set(sets, 1)
    for (x, S) in zip([[3], [1], [], [2]], MOI.Utilities.set_types(sets))
        ci = MOI.get(
            sets,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64},S}(),
        )
        @test ci == MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}.(x)
    end
    return
end

MOI.Utilities.@product_of_sets(
    _VectorSets,
    MOI.Nonpositives,
    MOI.Nonnegatives,
    MOI.EqualTo{T},
)

function test_vector_set_index()
    sets = _VectorSets{Float64}()
    @test MOI.Utilities.set_index(sets, MOI.Nonpositives) == 1
    @test MOI.Utilities.set_index(sets, MOI.Nonnegatives) == 2
    @test MOI.Utilities.set_index(sets, MOI.EqualTo{Float64}) == 3
    @test MOI.Utilities.set_index(sets, MOI.LessThan{Float64}) === nothing
end

function test_vector_set_types()
    sets = _VectorSets{Float64}()
    @test MOI.Utilities.set_types(sets) ==
          [MOI.Nonpositives, MOI.Nonnegatives, MOI.EqualTo{Float64}]
end

function test_vector_basic()
    sets = _VectorSets{Float64}()
    nonneg_i = MOI.Utilities.add_set(
        sets,
        MOI.Utilities.set_index(sets, MOI.Nonnegatives),
        2,
    )
    equalto_i = MOI.Utilities.add_set(
        sets,
        MOI.Utilities.set_index(sets, MOI.EqualTo{Float64}),
    )
    MOI.Utilities.final_touch(sets)
    VAF = MOI.VectorAffineFunction{Float64}
    @test !MOI.is_valid(sets, MOI.ConstraintIndex{VAF,MOI.Nonnegatives}(12345))
    nonneg_ci = MOI.ConstraintIndex{VAF,MOI.Nonnegatives}(nonneg_i)
    @test MOI.is_valid(sets, nonneg_ci)
    @test MOI.Utilities.rows(sets, nonneg_ci) == 1:2
    SAF = MOI.ScalarAffineFunction{Float64}
    equalto_ci = MOI.ConstraintIndex{SAF,MOI.EqualTo{Float64}}(equalto_i)
    @test MOI.is_valid(sets, equalto_ci)
    @test MOI.Utilities.rows(sets, equalto_ci) == 3
end

function test_vector_dimension()
    sets = _VectorSets{Float64}()
    MOI.Utilities.final_touch(sets)
    @test MOI.dimension(sets) == 0

    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 1, 1)
    MOI.Utilities.final_touch(sets)
    @test MOI.dimension(sets) == 1

    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 1, 1)
    MOI.Utilities.add_set(sets, 1, 2)
    MOI.Utilities.final_touch(sets)
    @test MOI.dimension(sets) == 3

    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 1, 1)
    MOI.Utilities.add_set(sets, 1, 2)
    MOI.Utilities.add_set(sets, 2, 3)
    MOI.Utilities.final_touch(sets)
    @test MOI.dimension(sets) == 6

    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 1, 1)
    MOI.Utilities.add_set(sets, 1, 2)
    MOI.Utilities.add_set(sets, 2, 3)
    MOI.Utilities.add_set(sets, 3)
    MOI.Utilities.final_touch(sets)
    @test MOI.dimension(sets) == 7
end

function test_vector_empty()
    sets = _VectorSets{Float64}()
    @test MOI.is_empty(sets)
    MOI.Utilities.add_set(sets, 1, 1)
    @test !MOI.is_empty(sets)
    MOI.empty!(sets)
    @test MOI.is_empty(sets)
end

function test_vector_ConstraintTypesPresent()
    sets = _VectorSets{Float64}()
    MOI.Utilities.final_touch(sets)
    @test MOI.get(sets, MOI.ListOfConstraintTypesPresent()) == []
    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 1, 1)
    MOI.Utilities.add_set(sets, 3)
    MOI.Utilities.final_touch(sets)
    @test MOI.get(sets, MOI.ListOfConstraintTypesPresent()) == [
        (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
    ]
end

function test_vector_NumberOfConstraints()
    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 1, 2)
    MOI.Utilities.add_set(sets, 1, 2)
    MOI.Utilities.add_set(sets, 2, 2)
    MOI.Utilities.add_set(sets, 3)
    MOI.Utilities.final_touch(sets)
    @test MOI.get(
        sets,
        MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},MOI.Zeros}(),
    ) == 0
    for (x, S) in zip([2, 1], MOI.Utilities.set_types(sets)[1:2])
        @test MOI.get(
            sets,
            MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},S}(),
        ) == x
    end
    @test MOI.get(
        sets,
        MOI.NumberOfConstraints{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
    ) == 1
    return
end

function test_vector_ListOfConstraintIndices()
    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 2, 2)
    MOI.Utilities.add_set(sets, 1, 4)
    MOI.Utilities.add_set(sets, 2, 3)
    MOI.Utilities.add_set(sets, 3)
    MOI.Utilities.final_touch(sets)
    VAF = MOI.VectorAffineFunction{Float64}
    @test MOI.get(sets, MOI.ListOfConstraintIndices{VAF,MOI.Zeros}()) ==
          MOI.ConstraintIndex{VAF,MOI.Zeros}[]
    for (x, S) in zip([[1], [1, 3]], MOI.Utilities.set_types(sets)[1:2])
        ci = MOI.get(sets, MOI.ListOfConstraintIndices{VAF,S}())
        @test ci == MOI.ConstraintIndex{VAF,S}.(x)
    end
    F, S = MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}
    @test MOI.get(sets, MOI.ListOfConstraintIndices{F,S}()) ==
          [MOI.ConstraintIndex{F,S}(1)]
    return
end

"""
    test_vector_ListOfConstraintIndices2()

Test a more complicated sequence of dimensions to check the `_UnevenIterator`
works appropriately.
"""
function test_vector_ListOfConstraintIndices2()
    sets = _VectorSets{Float64}()
    MOI.Utilities.add_set(sets, 2, 2)
    MOI.Utilities.add_set(sets, 2, 3)
    MOI.Utilities.add_set(sets, 2, 2)
    MOI.Utilities.add_set(sets, 2, 4)
    MOI.Utilities.final_touch(sets)
    S = MOI.Utilities.set_types(sets)[2]
    VAF = MOI.VectorAffineFunction{Float64}
    indices = MOI.get(sets, MOI.ListOfConstraintIndices{VAF,S}())
    @test indices == MOI.ConstraintIndex{VAF,S}.([1, 3, 6, 8])
end

end

TestProductOfSets.runtests()


text:
module OpenCV

include("typestructs.jl")
include("Vec.jl")
const dtypes = Union{UInt8, Int8, UInt16, Int16, Int32, Float32, Float64}
const size_t = UInt

using CxxWrap
@wrapmodule(joinpath(@__DIR__,"lib","libopencv_julia"), :cv_wrap)
function __init__()
    @initcxx
end
const Scalar = Union{Tuple{Number}, Tuple{Number, Number}, Tuple{Number, Number, Number}, NTuple{4, Number}}

include("Mat.jl")

const InputArray = Union{AbstractArray{T, 3} where {T <: dtypes}, CxxMat}

include("mat_conversion.jl")
include("types_conversion.jl")

function cpp_to_julia(var)
    return var
end
function julia_to_cpp(var)
    return var
end

# TODO: Vector to array conversion. There's some problems with the copy constructor.

# function julia_to_cpp(var::Array{T, 1}) where {T}
#     ret = CxxWrap.StdLib.StdVector{T}()
#     for x in var
#         push!(ret, julia_to_cpp(x)) # When converting an array keep expected type as final type.
#     end
#     return ret
# end

# function cpp_to_julia(var::CxxWrap.StdLib.StdVector{T}) where {T}
#     ret = Array{T, 1}()
#     for x in var
#         push!(ret, cpp_to_julia(x))
#     end
#     return ret
# end

function cpp_to_julia(var::Tuple)
    ret_arr = Array{Any, 1}()
    for it in var
        push!(ret_arr, cpp_to_julia(it))
    end
    return tuple(ret_arr...)
end
# using StaticArrays

include("cv_wrap.jl")
end


text:
using StatefulIterators
using Base.Test

function suite_123(make_iter; explicit=true, bits=true)

    print(make_iter(), " ... ")

    s = StatefulIterator(make_iter())
    @test read(s) == 1
    @test read(s) == 2
    c = explicit ? copy(s) : StatefulIterator(make_iter())
    @test read(s) == 3
    @test eof(s)

    if explicit
        @test !eof(c)
        @test read(c) == 3
        @test eof(c)
    else
        for x in (copy, peek, available)
            @test_throws ErrorException x(c)
        end
    end

    s = StatefulIterator(make_iter())
    if bits
        @test read(s, UInt8) == 0x01
        # TODO - document this behaviour (restarts on instance boundary)
        @test read(s, UInt8, 8) == [0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]
    else
        @test_throws Exception read(s, UInt8)
        @test_throws Exception read(s, UInt8, 1)
    end

    @test collect(take(StatefulIterator(make_iter()), 2)) == [1, 2]
    @test read(StatefulIterator(make_iter()), 2) == [1, 2]
    if bits
        @test read(StatefulIterator(make_iter()), UInt8, 2) == [0x01, 0x00]
    end

    s = StatefulIterator(make_iter())
    if explicit
        @test peek(s) == 1
        if bits
            @test peek(s, UInt8) == 0x01
            @test peek(s, UInt8, 2) == [0x01, 0x00]
        end
        @test read(s) == 1
        p = position(s)
        @test read(s) == 2
        seek(s, p)
        @test read(s) == 2
        seekend(s)
        @test done(s)
        seek(s, p)
        @test read(s) == 2
        @test available(s) == 1
        if bits
            @test available(s, UInt8) == sizeof(Int) / sizeof(UInt8)
        end
        seekstart(s)
        @test available(s) == 3
        if bits
            @test available(s, UInt8) == 3 * sizeof(Int) / sizeof(UInt8)
        end
        skip(s, 1)
        @test read(s) == 2
    else
        for x in (peek, position, available)
            @test_throws Exception x(s)
        end
        for x in (seek,)
            @test_throws Exception x(s, 1)
        end
        seekend(s)
        @test done(s)
    end

    @test size(read(StatefulIterator(make_iter()), 1, 2)) == (1, 2)
    @test size(read(StatefulIterator(make_iter()), 2, 1)) == (2, 1)

    s = StatefulIterator(make_iter())
    if bits
        a = zeros(Int, 2)
    else
        a = Any[0, 0]
    end
    read!(s, a)
    @test a == [1, 2]

    if bits
        s = StatefulIterator(make_iter())
        a = zeros(UInt8, 2)
        read!(s, a)
        @test a == [0x01, 0x00]
        s = StatefulIterator(make_iter())
        a = readbytes(s, 2)
        @test a == [0x01, 0x00]
        s = StatefulIterator(make_iter())
        a = zeros(UInt32, 3)
        read!(s, a)
        @test a == [0x00000001,0x00000000,0x00000002]
        s = StatefulIterator(make_iter())
        a = readbytes(s, 10)
        @test a == [0x01,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x02,0x00]
    end

    for n in 2:5
        s = StatefulIterator(make_iter())
        if n <= 3
            skip(s, n)
            @test done(s) == (n == 3)
        else
            @test_throws EOFError skip(s, n)
        end
    end

    println("ok")
end

suite_123(() -> [1,2,3])
suite_123(() -> 1:3)
suite_123(() -> 1:1:3)
suite_123(() -> Task(() -> (for i in 1:3; produce(i); end)); explicit=false, bits=false)
suite_123(() -> Any[1,2,3]; bits=false)

# chunking to larger bitlengths
@test read(StatefulIterator([0x01, 0x00, 0x00, 0x00]), UInt32) == one(UInt32)

# strange step size
@test available(StatefulIterator(1:2:7)) == length(collect(1:2:7))
for iter in ([1,2,3], 1:3, 1:2:7, 1:1:3, 1:2:6, 1:2:7, 2:2:3, 2:2:4)
    s = StatefulIterator(iter)
    skip(s, 1)
    @test collect(s) == collect(iter)[2:end]
    seekend(s)
    @test done(s)
end

# bug in GoCL
e = zeros(UInt8, 7)
e[1:4] = map(UInt8, collect("goxp"))
e[5] = 0x00
e[6:7] = reinterpret(UInt8, UInt16[1000])
s = StatefulIterator(e)
@test read(s, 4) == map(UInt8, collect("goxp"))
@test read(s) == 0x00
@test read(s, UInt16) == 1000

# text only
s = StatefulIterator("hello\nworld\n")
@test readline(s) == "hello\n"

function suite_speed(make_iter, target)
    s = make_iter()
    println(typeof(s))
    n = sum(s)
    @test n == target
    s = make_iter()
    @time sum(s)
end

suite_speed(() -> ones(Int, Int(1e6)), 1e6)
suite_speed(() -> StatefulIterator(ones(Int, Int(1e6))), 1e6)


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using Random
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
	using ViewerGL
	GL = ViewerGL
	include("../examples/3d/views.jl")
end

function build_copFE(FV::Lar.Cells, EV::Lar.Cells)
	copFE = Lar.u_coboundary_1(FV, EV) # unsigned
	faceedges = [findnz(copFE[f,:])[1] for f=1:size(copFE,1)]

	f_edgepairs = Array{Array{Int64,1}}[]
	for f=1:size(copFE,1)
		edgepairs = Array{Int64,1}[]
		for v in FV[f]
			push!(edgepairs, [e for e in faceedges[f] if v in EV[e]])
		end
		push!(f_edgepairs, edgepairs)
	end
	copFE.=0
	for f=1:size(copFE,1)
		facepairs = sort(sort.(f_edgepairs[f]))
		# setting a reference
		copFE[f, facepairs[1][1]] = 1
		for (e1,e2) in facepairs
			@assert copFE[f, e1] != 0
			if copFE[f, e2] == 0
				copFE[f, e2] =
					-((EV[e1][1]==EV[e2][1])||(EV[e1][2]==EV[e2][2])) * copFE[f,e1]+
					((EV[e1][1]==EV[e2][2])||(EV[e1][2]==EV[e2][1])) * copFE[f,e1]
			end
		end
	end
	return copFE
end

Random.seed!(1234);

scaling = 1.2
no_sphere = 10
expl = 2.0
pt1, pt2 = 5, 10
#pt1, pt2 = 10, 10

carry = Array{Tuple{Array{Float64,2},Array{Array{Int64,1},1}}}(undef, no_sphere)
for i = 1 : no_sphere
	carry[i] = Lar.struct2lar(Lar.Struct([
		Lar.t(rand(0.5:1e-2:5, 3)...),
		Lar.sphere(rand(1:1e-2:3) * scaling)([pt1, pt2])
	]))
end
V, FV = Lar.struct2lar(Lar.Struct(carry))
FV = sort(sort.(FV))

model = CAGD.Model(V);
function getFaceEdges(fV::Array{Int,1})
    return [fV[1], fV[2]], [fV[1], fV[3]], [fV[2], fV[3]]
end
EV = unique([(map(f -> getFaceEdges(f), FV)...)...])
FE = [[3*i+1, 3*i+2, 3*i+3] for i = 0 : length(FV)-1]
CAGD.addModelCells!(model, 1, EV, signed = true)
CAGD.addModelCells!(model, 2, build_copFE(FV, EV))
model = CAGD.mergeModelVertices(model, signed_merge=true)

if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model)
FC, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model)  end


text:
#
# Function to set bounds with optional user input
#

function setbounds(vars,l,u)
  # Throw error if the user set a bound to a variable which is not in list for this function
  for field in fieldnames(typeof(l))
    if (getfield(l,field) != nothing) || (getfield(u,field) !=nothing)
      found = false
      for var in vars
        if field == var.field
          found = true
          if var.type == Nothing
            error(" Bounds to intercept $(var.field) are not supported. For a constant use $(var.field) = const.")
          end
        end
      end
      if ! found
        error(" A bound was set to variable '$field', but '$field' is not a variable of the current fit function.")
      end
    end
  end
  # Total number variables
  n = 0
  for var in vars
    n += var.dim
  end
  lower = Vector{Float64}(undef,n)
  upper = Vector{Float64}(undef,n)
  idim = 1
  for var in vars
    ltmp = check_lower_bound_input(var,getfield(l,var.field))
    utmp = check_upper_bound_input(var,getfield(u,var.field))
    @. lower[idim:idim+var.dim-1] = ltmp
    @. upper[idim:idim+var.dim-1] = utmp
    idim += var.dim
  end
  i = 0
  for var in vars
    i += 1
    if lower[i] > upper[i]
      error(" Error in bounds. Lower bound of '$(var.field)' greater than upper bound. ")
    end
  end
  return lower, upper
end

function check_lower_bound_input(var,value)
  if value == nothing
    return [ -Inf for _ in 1:var.dim ]
  end
  if ! (value isa var.type )
    error("Lower bound of $(var.field) must be of type $(var.type), got $(typeof(value))")
  end
  if (value isa Vector) && (length(value) != var.dim)
    error("Lower bound of $(var.field) must be of dimension $(var.dim), got $(length(value))")
  end
  return value
end

function check_upper_bound_input(var,value)
  if value == nothing
    return [ +Inf for _ in 1:var.dim ]
  end
  if ! (value isa var.type)
    error("Upper bound of $(var.field) must be of type $(var.type), got $(typeof(value))")
  end
  if (value isa Vector) && (length(value) != var.dim)
    error("Upper bound of $(var.field) must be of dimension $(var.dim), got $(length(value))")
  end
  return value
end


text:
module Objects

export Object, Object1D, Object2D,
       Curve, BCurve, NCurve,
       Surface, BSurface, NSurface,
       evalpoint,
       evalderiv1, evalderiv2, evalderiv3,
       evaltangent, evalnormal,
       isocurve

# Abstract Basis types
# ========================================================================

abstract type Object end
abstract type Object1D <: Object end
abstract type Object2D <: Object end

mutable struct ObjectFunction1D{O<:Object1D}
    object::O
    index::Int
    deriv::Int

    function ObjectFunction1D{O}(object, index, deriv) where {O<:Object1D}
        if !(1 <= index <= length(object.basis)) throw(BoundsError()) end
        if !(0 <= deriv <= nderivs(object.basis)) throw(ArgumentError("Differentiation order not supported")) end
        new(object, index, deriv)
    end

    ObjectFunction1D(object, index) = ObjectFunction1D(object, index, 0)
end

mutable struct ObjectFunction2D{O<:Object2D}
    object::O
    indexU::Int
    indexV::Int
    deriv::Int

    function ObjectFunction2D{O}(object, indexU, indexV, deriv) where {O<:Object1D}
        if !(1 <= indexU <= size(object.basis, 1)) throw(BoundsError()) end
        if !(1 <= indexU <= size(object.basis, 2)) throw(BoundsError()) end
        if !(0 <= deriv <= nderivs(object.basis)) throw(ArgumentError("Differentiation order not supported")) end
        new(object, indexU, indexV, deriv)
    end

    ObjectFunction2D(object, indexU, indexV) = ObjectFunction2D(object, indexU, indexV, 0)
end

(o::Object)(u::T...) where {T<:Real} = evalpoint(o, u...)

include("Curves.jl")
include("Surfaces.jl")

end  # module Objects


text:
using Thermodynamics
using Test

function test_peos(eos::Thermodynamics.parameterizedeos, δ, τ)
    print("Testing module parameterized_eos:\n")
    print("δ = "*string(δ))
    print("\nτ = "*string(τ))
    print("\nϕ₀ = "*string(ϕ₀(eos, δ, τ)))
    print("\nϕ₀_δ = "*string(ϕ₀_δ(eos, δ, τ)))
    print("\nϕ₀_δδ = "*string(ϕ₀_δδ(eos, δ, τ)))
    print("\nϕ₀_τ = "*string(ϕ₀_τ(eos, δ, τ)))
    print("\nϕ₀_ττ = "*string(ϕ₀_ττ(eos, δ, τ)))
    print("\nϕ₀_δτ = "*string(ϕ₀_δτ(eos, δ, τ)))
    print("\nϕr = "*string(ϕr(eos, δ, τ)))
    print("\nϕr_δ = "*string(ϕr_δ(eos, δ, τ)))
    print("\nϕr_δδ = "*string(ϕr_δδ(eos, δ, τ)))
    print("\nϕr_τ = "*string(ϕr_τ(eos, δ, τ)))
    print("\nϕr_ττ = "*string(ϕr_ττ(eos, δ, τ)))
    print("\nϕr_δτ = "*string(ϕr_δτ(eos, δ, τ)))
    print("\n\nNew test conditions:")
    ρ = 467.5483
    T = 305
    @printf "\nT = %f K, ρ = %f kg/m³" T ρ

    p_ = 7.5e6
    print("\np = "*string(p_))
    p_ = p(eos, ρ, T)
    print("\np = "*string(p_))
    print("\nSolving in reverse, with p = "string(p_)*" and T="*string(T))
    #print(ρ_from_p(eos, p_, 300))
    ρ_from_p(eos, p_, T)
end


text:
module PlanarVisibility

import Base: length, push!, getindex
using LinearAlgebra

using DataStructures
using GeoInterface
using LightGraphs

include("points.jl")
include("intersect.jl")

#
# Types
#

"Environment of obstacles."
struct Environment
    # TODO: make types parametric?
    # TODO: support linestrings?
    polygons::Array{Polygon}
end

"Visibility graph."
struct VisibilityGraph
    env::Environment               # initial environment
    pointset::PointSet             # all points (including environment)
    graph::SimpleGraph{Int64}      # edges of visibility relation
end

VisibilityGraph(env::Environment) = construct_graph(env)

#
# Utilities
#

"Extract unique points from given environment."
function extract_points(env::Environment)
    points = PointSet()
    for polygon in env.polygons
        push!(points, polygon)
    end
    return points
end

"Extract edges from given environment."
function extract_edges(env::Environment, points::PointSet)
    # start with disconnected graph
    n = length(points)
    graph = SimpleGraph(n)

    # iterate over all polygons
    for polygon in env.polygons
        # iterate over outer and inner line strings (coordinates)
        for str in coordinates(polygon)
            # iterate over pairs of positions, with wrap-around
            for (pos1, pos2) in zip(str, [str[2:end]; str[1:1]])
                pos1 == pos2 && continue # skip with explicit close
                add_edge!(graph, points[pos1], points[pos2])
            end
        end
    end

    return graph
end

"Extract connected components (via polygons) from environment."
function extract_components(env::Environment, points::PointSet)
    # start with disjoint set forest
    components = IntDisjointSets(length(points))

    # for each polygon, merge all points
    for polygon in env.polygons
        linestrings = coordinates(polygon)
        length(linestrings) == 1 || error("nonconvex polygons not implemented!")
        positions = linestrings[1]
        first = points[positions[1]] # get index
        for other in positions[2:end]
            union!(components, first, points[other])
        end
    end

    return components
end

"Angle from origin to point rel. to horiz. line going east."
function angle_to_east(origin::Point, point::Point)::Float64
    dx, dy = coordinates(point) - coordinates(origin)
    angle = asin(abs(dy) / norm([dx, dy]))
    if dy > 0.0 # top
        if dx >= 0.0 # right
            return angle
        else # dx < 0, left
            return π - angle
        end
    elseif dy < 0.0
        if dx >= 0.0 # right
            return 2π - angle
        else # dx < 0, left
            return π + angle
        end
    else # dy == 0.0
        if dx >= 0.0
            return 0.0
        else # dx < 0.0
            return π
        end
    end
end

"Find permutation to iterate thru points in counterclockwise order."
function sortperm_ccw(points::Vector{Point}, origin::Point)::Vector{Int}
    by = p -> angle_to_east(origin, p)
    return sortperm(points, by=by)
end

"Sorted edge for unique representation"
function sorted(edge::Edge)
    if edge.src < edge.dst
        return edge
    else
        return reverse(edge)
    end
end

#
# Constructing visibility graphs
#

"Find all visible points from given origin."
function visible_points(points::PointSet, envgraph::SimpleGraph,
                        polycomps::IntDisjointSets, origin::Int64)
    # point at origin
    orig = points[origin]

    # starting with east-looking horizontal ray
    xmax = maximum([xcoord(coordinates(p)) for p in points.points])
    ray = Point([xmax + 1.0, ycoord(coordinates(orig))])

    # currently intersecting edges
    open_edges = Set{Edge}()
    for e in edges(envgraph)
        # skip non-intersecting edges
        intersect_segments(orig, ray, points[e.src], points[e.dst]) || continue
        # skip those edges where the endpoint is on the ray
        !collinear_intersect(orig, points[e.src], ray) || continue
        !collinear_intersect(orig, points[e.dst], ray) || continue
        push!(open_edges, sorted(e))
    end

    # counterclockwise sweep to consider all candidates
    perm = sortperm_ccw(points.points, orig)
    result = Int[]
    for index in perm
        index == origin && continue
        cand = points[index]

        # open edges: rm incident, clockwise edges
        for neighbor in neighbors(envgraph, index)
            edge = sorted(Edge(index, neighbor))
            edge in open_edges || continue
            if orientation(orig, cand, points[neighbor]) == CLOCKWISE
                delete!(open_edges, edge)
            end
        end

        visible = true
        for e in open_edges
            if intersect_segments(orig, cand, points[e.src], points[e.dst])
                visible = false
                break
            end
        end

        # check if line-of-sight is inside of polygon
        if visible && !has_edge(envgraph, origin, index)
            if in_same_set(polycomps, origin, index)
                visible = false
            end
        end
        # TODO: deal with non-convex polygons

        if visible
            push!(result, index)
        end

        # open edges: add incident, counterclockwise edges
        for neighbor in neighbors(envgraph, index)
            if orientation(orig, cand, points[neighbor]) == COUNTERCLOCKWISE
                push!(open_edges, sorted(Edge(index, neighbor)))
            end
        end
    end

    return result
end

function construct_graph(env::Environment)
    # extract points from environment
    points = extract_points(env)

    # extract edges from environment
    envgraph = extract_edges(env, points)

    # extract connected componenst from environment polygons
    polycomps = extract_components(env, points)

    # build graph from disconnected points
    n = length(points)
    graph = SimpleGraph(n)

    # iterate over all points and find out what can be seen from there
    for i in 1:n
        for j in visible_points(points, envgraph, polycomps, i)
            add_edge!(graph, i, j)
        end
    end

    # assemble result
    return VisibilityGraph(env, points, graph)
end

#
# Extending visibility graphs
#

end


text:
## Loading packages.
using MicroHH
using Tullio
using Statistics
using PyPlot
pygui(:qt5)

## Solve the Taylor Green vortex
function solve_taylorgreen(itot)
    # Loading settings.
    d_grid = Dict{String, Any}()
    d_grid["itot"] = itot
    d_grid["jtot"] = 1
    d_grid["ktot"] = itot ÷ 2

    d_grid["xsize"] = 1.
    d_grid["ysize"] = 1.
    d_grid["zsize"] = 0.5

    zsize = d_grid["zsize"]; ktot = d_grid["ktot"]
    dz = zsize / ktot
    d_grid["z"] = range(0.5*dz, step=dz, length=ktot) |> collect

    d_fields = Dict{String, Any}()
    d_fields["visc"] = (8 * pi^2 * 100)^(-1)
    d_fields["alpha"] = 0.

    d_boundary = Dict{String, Any}()
    d_boundary["mom_bot_type"] = "Neumann"
    d_boundary["mom_top_type"] = "Neumann"
    d_boundary["s_bot_type"] = "Neumann"
    d_boundary["s_top_type"] = "Neumann"

    d_timeloop = Dict{String, Any}()
    d_timeloop["start_time"] = 0.
    d_timeloop["end_time"] = 1.
    d_timeloop["save_time"] = 100.
    d_timeloop["check_time"] = 0.1
    d_timeloop["dt"] = 0.005 * (128 / itot)

    d_multidomain = Dict{String, Any}()
    d_multidomain["enable_nudge"] = false

    settings = [ Dict(
        "grid" => d_grid,
        "fields" => d_fields,
        "boundary" => d_boundary,
        "timeloop" => d_timeloop,
        "multidomain" => d_multidomain) ]

    # Initialize the model.
    n_domains = 1
    m = Model("taylorgreen", n_domains, settings, Float64)

    # Create the initials fields.
    f = m.fields[1]; g = m.grid[1]
    u = @view f.u[g.is:g.ie, g.js:g.je, g.ks:g.ke]
    w = @view f.w[g.is:g.ie, g.js:g.je, g.ks:g.keh]
    p = @view f.p[g.is:g.ie, g.js:g.je, g.ks:g.ke]
    x = @view g.x[g.is:g.ie]
    xh = @view g.xh[g.is:g.ie]
    z = @view g.z[g.ks:g.ke]
    zh = @view g.zh[g.ks:g.keh]
    @tullio u[i, j, k] = sin(2pi*xh[i])*cos(2pi*z[k])
    @tullio w[i, j, k] = -cos(2pi*x[i])*sin(2pi*zh[k])

    f_t = exp(-8. * pi^2 * d_fields["visc"] * d_timeloop["end_time"])
    u_ref = u[:, :, :] * f_t
    w_ref = w[:, :, :] * f_t
    p_ref = zeros(size(u_ref))
    @tullio p_ref[i, j, k] = (1/4 * (cos(4pi*x[i]) + cos(4pi*z[k])) - 1/4) * f_t^2

    # Run the model
    in_progress = prepare_model!(m)
    while in_progress
        in_progress = step_model!(m)
    end

    dx = 1. / itot
    u_err = sum(dx^2 * abs.(u .- u_ref))
    w_err = sum(dx^2 * abs.(w .- w_ref))
    p_err = sum(dx^2 * abs.(p .- p_ref))

    return u_err, w_err, p_err
end

## Solve the case at multiple resolutions.
u016, w016, p016 = solve_taylorgreen(16)
u032, w032, p032 = solve_taylorgreen(32)
u064, w064, p064 = solve_taylorgreen(64)
u128, w128, p128 = solve_taylorgreen(128)
u256, w256, p256 = solve_taylorgreen(256)

dxs = [ 1/16, 1/32, 1/64, 1/128, 1/256 ]
u_errs = [ u016, u032, u064, u128, u256 ]
w_errs = [ w016, w032, w064, w128, w256 ]
p_errs = [ p016, p032, p064, p128, p256 ]

## Plot the results
figure()
loglog(dxs, u_errs, "C0o-", label="u")
loglog(dxs, w_errs, "C1o-", label="w")
loglog(dxs, p_errs, "C2o-", label="p")
loglog(dxs, 0.6*dxs.^2, "k:", label="2nd")
legend(loc=0, frameon=false)
xlabel(L"$\Delta$ x")
ylabel("Error")
display(gcf())
show()


text:
using Oceananigans.Fields: validate_indices, Reduction
using Oceananigans.AbstractOperations: AbstractOperation
using Oceananigans.Grids: default_indices

restrict_to_interior(::Colon, loc, topo, N) = interior_indices(loc, topo, N)
restrict_to_interior(::Colon, ::Type{Nothing}, topo, N) = UnitRange(1, 1)
restrict_to_interior(index::UnitRange, ::Type{Nothing}, topo, N) = UnitRange(1, 1)

function restrict_to_interior(index::UnitRange, loc, topo, N)
    from = max(first(index), 1)
    to = min(last(index), last(interior_indices(loc, topo, N)))
    return UnitRange(from, to)
end

#####
##### Function output fallback
#####

function construct_output(output, grid, indices, with_halos)
    if !(indices isa typeof(default_indices(3)))
        output_type = output isa Function ? "Function" : ""
        @warn "Cannot slice $output_type $output with $indices: output will be unsliced."
    end

    return output
end

#####
##### Support for Field, Reduction, and AbstractOperation outputs
#####

function output_indices(output::Union{AbstractField, Reduction}, grid, indices, with_halos)
    indices = validate_indices(indices, location(output), grid)

    if !with_halos # Maybe chop those indices
        loc = location(output)
        topo = topology(grid)
        indices = restrict_to_interior.(indices, loc, topo, size(grid))
    end

    return indices
end

function construct_output(user_output::Union{AbstractField, Reduction}, grid, user_indices, with_halos)
    indices = output_indices(user_output, grid, user_indices, with_halos)
    return construct_output(user_output, indices)
end

construct_output(user_output::Field, indices) = view(user_output, indices...)
construct_output(user_output::Reduction, indices) = Field(user_output; indices)
construct_output(user_output::AbstractOperation, indices) = Field(user_output; indices)

#####
##### Time-averaging
#####

function construct_output(averaged_output::WindowedTimeAverage{<:Field}, grid, indices, with_halos)
    output = construct_output(averaged_output.operand, grid, indices, with_halos)
    return WindowedTimeAverage(output; schedule=averaged_output.schedule)
end


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

for typ in atomictypes
    lt = llvmtypes[typ]
    ilt = llvmtypes[inttype(typ)]
    rt = "$lt, $lt*"
    irt = "$ilt, $ilt*"
    @eval getindex(x::Atomic{$typ}) =
        llvmcall($"""
                 %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                 %rv = load atomic $rt %ptr acquire, align $(gc_alignment(typ))
                 ret $lt %rv
                 """, $typ, Tuple{Ptr{$typ}}, unsafe_convert(Ptr{$typ}, x))
    @eval setindex!(x::Atomic{$typ}, v::$typ) =
        llvmcall($"""
                 %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                 store atomic $lt %1, $lt* %ptr release, align $(gc_alignment(typ))
                 ret void
                 """, Cvoid, Tuple{Ptr{$typ}, $typ}, unsafe_convert(Ptr{$typ}, x), v)

    # Note: atomic_cas! succeeded (i.e. it stored "new") if and only if the result is "cmp"
    if typ <: Integer
        @eval atomic_cas!(x::Atomic{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                     %rs = cmpxchg $lt* %ptr, $lt %1, $lt %2 acq_rel acquire
                     %rv = extractvalue { $lt, i1 } %rs, 0
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     unsafe_convert(Ptr{$typ}, x), cmp, new)
    else
        @eval atomic_cas!(x::Atomic{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                     %icmp = bitcast $lt %1 to $ilt
                     %inew = bitcast $lt %2 to $ilt
                     %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                     %irv = extractvalue { $ilt, i1 } %irs, 0
                     %rv = bitcast $ilt %irv to $lt
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     unsafe_convert(Ptr{$typ}, x), cmp, new)
    end

    arithmetic_ops = [:add, :sub]
    for rmwop in [arithmetic_ops..., :xchg, :and, :nand, :or, :xor, :max, :min]
        rmw = string(rmwop)
        fn = Symbol("atomic_", rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        if rmwop in arithmetic_ops && !(typ <: ArithmeticTypes) continue end
        if typ <: Integer
            @eval $fn(x::Atomic{$typ}, v::$typ) =
                llvmcall($"""
                         %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                         %rv = atomicrmw $rmw $lt* %ptr, $lt %1 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, unsafe_convert(Ptr{$typ}, x), v)
        else
            rmwop === :xchg || continue
            @eval $fn(x::Atomic{$typ}, v::$typ) =
                llvmcall($"""
                         %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         %irv = atomicrmw $rmw $ilt* %iptr, $ilt %ival acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, unsafe_convert(Ptr{$typ}, x), v)
        end
    end
end


text: global const A=1


text:
#!/usr/bin/env julia

using Celeste
using Celeste: SDSSIO
using Gasp
using Gasp: grank
using Base.Threads

include("../bin/binutil.jl")

function run_io_simulation(args::Vector{String})
    if length(args) != 3
        println("""
            Usage:
              simulate_io.jl <strategy> <rcf_nsrcs_file> <boxes_file>
            """)
        exit(-1)
    end
    strategy, all_rcfs, all_rcf_nsrcs = decide_strategy(args[1], args[2])

    # parse the specified box file
    boxfile = args[3]
    boxes = Vector{BoundingBox}()
    f = open(boxfile)
    for ln in eachline(f)
        lp = split(ln, '\t')
        if length(lp) < 4
            Log.one_message("ERROR: malformed line in box file:\n> $ln ")
            continue
        end

        ss = split(lp[4], ' ')
        ramin = parse(Float64, ss[1])
        ramax = parse(Float64, ss[2])
        decmin = parse(Float64, ss[3])
        decmax = parse(Float64, ss[4])
        bb = BoundingBox(ramin, ramax, decmin, decmax)
        push!(boxes, bb)
    end
    close(f)

    field_extents = Celeste.ParallelRun.load_field_extents(strategy)

    dt, _ = Dtree(length(boxes), 0.25)
    ni, (ci, li) = initwork(dt)
    rundt = runtree(dt)
    Celeste.Log.info("dtree: initial: $(ni) ($(ci) to $(li))")
    l = SpinLock()
    iol = SpinLock()
    function do_work()
        try
            if rundt && threadid() == nthreads()
                Celeste.Log.info("dtree: running tree")
                while runtree(dt)
                    Gasp.cpu_pause()
                    ccall(:jl_gc_safepoint, Void, ())
                end
                return
            end

            while true
                box_idx = 0
                lock(l)
                while true
                    li == 0 && break
                    if ci > li
                        ni, (ci, li) = try getwork(dt)
                        catch exc
                            Celeste.Log.exception(exc)
                            break
                        end
                        if li == 0
                            Celeste.Log.info("dtree: out of work")
                        else
                            Celeste.Log.info("dtree: $(ni) work items ($(ci) to $(li))")
                        end
                    else
                        box_idx = ci
                        ci += 1
                        break
                    end
                end
                unlock(l)
                box_idx == 0 && return

                box = boxes[box_idx]
                rcfs = Celeste.ParallelRun.get_overlapping_fields(box, field_extents)
                tic()
                for rcf in rcfs
                    # Do this one RCF at a time, we don't have the memory to keep all images
                    states = SDSSIO.preload_rcfs(strategy, [rcf])
                    SDSSIO.read_photoobj(strategy, rcf, states[], drop_quickly = true)
                    SDSSIO.load_field_images(strategy, [rcf], states, #= drop_quickly = =# true)
                end
                bytes, mb = Base.prettyprint_getunits(Sys.maxrss(), length(Base._mem_units), Int64(1024))
                lock(iol)
                Celeste.Log.message("Loaded box $box_idx ($(box.ramin) $(box.ramax) $(box.decmin) $(box.decmax)) with $(length(rcfs)) RCFs in $(toq())s, maxrss is $bytes $(Base._mem_units[mb])")
                unlock(iol)
                # Force full GC to drop all the memory we allocated above
            end
        catch exc
            Celeste.Log.exception(exc)
        end
    end
    ccall(:jl_threading_run, Void, (Any,), Core.svec(do_work))
end
run_io_simulation(ARGS)


text:
#------------------------------------------------------#
# Build all constraints of the LaplacianOptModel here  #
#------------------------------------------------------#

function constraint_build_W_var_matrix(lom::LaplacianOptModel)
    n = lom.data["num_nodes"]
    ew_matrix = lom.data["edge_weights"]

    # Diagonal entries
    JuMP.@constraint(lom.model, [i=1:n], lom.variables[:W_var][i,i] == sum(ew_matrix[i,:] .* lom.variables[:z_var][i,:]) - lom.variables[:γ_var]*(n-1)/n)

    # Off-diagonal entries
    JuMP.@constraint(lom.model, [i=1:(n-1), j=(i+1):n], lom.variables[:W_var][i,j] == -ew_matrix[i,j] * lom.variables[:z_var][i,j] + lom.variables[:γ_var]/n)

    return
end

function constraint_topology_no_self_loops(lom::LaplacianOptModel)
    n = lom.data["num_nodes"]

    JuMP.@constraint(lom.model, [i=1:n], lom.variables[:z_var][i,i] == 0)

    return
end

function constraint_topology_vertex_cutset(lom::LaplacianOptModel)
    n = lom.data["num_nodes"]

    JuMP.@constraint(lom.model, [i=1:n], sum(lom.variables[:z_var][i,j] for j in setdiff((1:n), i)) >= 1)

    return
end

function constraint_topology_total_edges(lom::LaplacianOptModel)
    n = lom.data["num_nodes"]

    JuMP.@constraint(lom.model, sum(lom.variables[:z_var][i,j] for i=1:(n-1), j=(i+1):n) == (n-1))

    return
end

function constraint_lazycallback_wrapper(lom::LaplacianOptModel)

    function constraint_lazycallback_cuts(cb_cuts)

        status = JuMP.callback_node_status(cb_cuts, lom.model)

        # if status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL

        #     if lom.data["optimizer"] != "glpk"
        #         Memento.info(_LOGGER, "Callback status: solution is fractional")
        #     end

        if status == MOI.CALLBACK_NODE_STATUS_UNKNOWN
            Memento.error(_LOGGER, "Callback status: unknown - solution may not be integer feasible")

        else status == MOI.CALLBACK_NODE_STATUS_INTEGER

            if lom.data["eigen_cuts_full"]
                W_val = JuMP.callback_value.(Ref(cb_cuts), lom.variables[:W_var])

                LOpt.constraint_eigen_cuts_on_full_matrix(W_val, cb_cuts, lom)
            end

            if lom.data["soc_linearized_cuts"]
                W_val = JuMP.callback_value.(Ref(cb_cuts), lom.variables[:W_var])

                LOpt.constraint_soc_linearized_cuts_on_2minors(W_val, cb_cuts, lom)
            end

            if lom.data["topology_flow_cuts"]
                z_val = abs.(JuMP.callback_value.(Ref(cb_cuts), lom.variables[:z_var]))

                LOpt.get_rounded_zeros_and_ones!(z_val, lom.data["tol_zero"])

                LOpt.constraint_topology_flow_cuts(z_val, cb_cuts, lom)
            end

        end

    end

    if lom.data["lazy_callback_status"]
        MOI.set(lom.model, MOI.LazyConstraintCallback(), constraint_lazycallback_cuts)
    end

end

function constraint_eigen_cuts_on_full_matrix(W_val::Matrix{Float64}, cb_cuts, lom::LaplacianOptModel)

    if typeof(LA.eigvals(W_val)) != Vector{Float64}
        Memento.error(_LOGGER, "PSD matrix (W) cannot have complex eigenvalues")
    end

    if LA.eigmin(W_val) <= -lom.data["tol_psd"]

        if LA.eigvals(W_val)[1] == LA.eigmin(W_val)

            violated_eigen_vec = LA.eigvecs(W_val)[:,1]

            LOpt.get_rounded_zeros_and_ones!(violated_eigen_vec, lom.data["tol_zero"])

            con = JuMP.@build_constraint(violated_eigen_vec' * lom.variables[:W_var] * violated_eigen_vec >= 0)

            MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con)

            if lom.data["lazycuts_logging"]
                Memento.info(_LOGGER, "Polyhedral relaxation cuts: full-sized eigen")
            end

        else

            Memento.warn(_LOGGER, "Eigen cut corresponding to the negative eigenvalue could not be added")
        end
    end

    return
end

function constraint_soc_linearized_cuts_on_2minors(W_val::Matrix{Float64}, cb_cuts, lom::LaplacianOptModel)

    num_nodes = lom.data["num_nodes"]

    for i = 1:(num_nodes-1)
        for j = (i+1):num_nodes

            W_val_minor = W_val[i,j]^2 - W_val[i,i] * W_val[j,j]

            if  W_val_minor >= lom.data["tol_zero"]

                coefficient_vec = [-(W_val[i,j])^2
                                   2 * W_val[i,j] * W_val[i,i]]

                variable_vec = [lom.variables[:W_var][i,i], lom.variables[:W_var][i,j]]

                # Taylor's approximation of rotated SOC constraint (x^2 <= y*z)
                if !isapprox(W_val[i,i], 0, atol=1E-6)
                    con = JuMP.@build_constraint(coefficient_vec[1]*variable_vec[1] + coefficient_vec[2]*variable_vec[2]
                                                    <= W_val[i,i]^2 * lom.variables[:W_var][j,j])

                    MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con)
                end

                if lom.data["lazycuts_logging"]
                    Memento.info(_LOGGER, "Polyhedral relaxation cuts: 2x2 minor linearized SOC")
                end

            end

        end
    end

    return
end

function constraint_topology_flow_cuts(z_val::Matrix{Float64}, cb_cuts, lom::LaplacianOptModel)

    cc_lazy = Graphs.connected_components(Graphs.SimpleGraph(abs.(z_val)))

    if length(cc_lazy) == 2

        con = JuMP.@build_constraint(sum(lom.variables[:z_var][i,j] for i in cc_lazy[1], j in cc_lazy[2]) >= 1)

        MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con)

        if lom.data["lazycuts_logging"]
            Memento.info(_LOGGER, "Polyhedral relaxation cuts: graph connectivity (#cc = $(length(cc_lazy)))")
        end

    elseif length(cc_lazy) == 3

        con1 = JuMP.@build_constraint(sum(lom.variables[:z_var][i,j] for i in cc_lazy[1], j in cc_lazy[2]) >= 1)
        con2 = JuMP.@build_constraint(sum(lom.variables[:z_var][i,j] for i in cc_lazy[2], j in cc_lazy[3]) >= 1)

        MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con1)
        MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con2)

        if lom.data["lazycuts_logging"]
            Memento.info(_LOGGER, "Polyhedral relaxation cuts: graph connectivity (#cc = $(length(cc_lazy)))")
        end

    elseif length(cc_lazy) >= 4

        Memento.info(_LOGGER, "Polyhedral relaxation cuts: flow cuts not added for integer solutions with $(length(cc_lazy)) connected components")

    end

    return
end

function constraint_topology_multi_commodity_flow(lom::LaplacianOptModel)

    n = lom.data["num_nodes"]

    # Flow balance on source node
    JuMP.@constraint(lom.model, [k=1:(n-1)], sum(lom.variables[:flow_var][1,2:n,k]) - sum(lom.variables[:flow_var][2:n,1,k]) == 1)

    # Flow balance on target node
    JuMP.@constraint(lom.model, [k=1:(n-1)], sum(lom.variables[:flow_var][:,k+1,k]) - sum(lom.variables[:flow_var][k+1,:,k]) == 1) #(k-th commodity corresponds to k+1th node)

    # Flow balance on remaining nodes
    JuMP.@constraint(lom.model, [k=1:(n-1), v=(2:n); k!=v-1], sum(lom.variables[:flow_var][:,v,k]) - sum(lom.variables[:flow_var][v,:,k]) == 0)

    # Flow <= capacity
    JuMP.@constraint(lom.model, [k=1:(n-1), i=1:n, j=1:n, kdash=1:(n-1)], lom.variables[:flow_var][i,j,k] + lom.variables[:flow_var][j,i,kdash]' <= lom.variables[:z_var][i,j])

    return
end


text:
module FindInterval

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export find_interval

function find_interval(m)

    # Find max and min part of the guillotine cut to determine its height

    # flbx

    # MAX
    @variable(m, max_cut_flbx_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flbx_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flbx_left_ij][i, j] >=
        m[:guillotine_flbx][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flbx_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flbx_left_i[i=1:length(boxes)-1],
        m[:max_cut_flbx_left_i][i] .>=
        [m[:max_cut_flbx_left_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )

    @constraint(m, c_max_cut_flbx_left_ii[i=1:length(boxes)-1],
        m[:max_cut_flbx_left_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    @variable(m, max_cut_flbx_j_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flbx_j_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flbx_j_plus_left_ij][i, j] >=
        m[:guillotine_flbx_j_plus][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flbx_j_plus_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flbx_j_plus_left_i[i=1:length(boxes)-1],
        m[:max_cut_flbx_j_plus_left_i][i] .>=
        [m[:max_cut_flbx_j_plus_left_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )

    @constraint(m, c_max_cut_flbx_j_plus_left_ii[i=1:length(boxes)-1],
        m[:max_cut_flbx_j_plus_left_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    # MIN using two possibilities
    @variable(m, min_cut_flbx_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flbx_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flbx_left_ij][i, j] .>=
        [m[:guillotine_flbx][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flbx][i, j]) * m[:flbz_front_left_bot_z][i]]

    )
    @variable(m, min_cut_flbx_j_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flbx_j_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flbx_j_plus_left_ij][i, j] .>=
        [m[:guillotine_flbx_j_plus][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flbx_j_plus][i, j]) * m[:flbz_front_left_bot_z][i]]

    )

    #Find the lowest box of the j other boxes
    @variable(m, min_cut_flbx_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flbx_left_i[i=1:length(boxes)-1],
        [m[:min_cut_flbx_left_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flbx_left_i][i]
    )

    @variable(m, min_cut_flbx_j_plus_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flbx_j_plus_left_i[i=1:length(boxes)-1],
        [m[:min_cut_flbx_j_plus_left_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flbx_j_plus_left_i][i]
    )

    #Need to ensure that the min_guillotine cut must be lower than or equal to z_i
    @constraint(m, c_min_cut_flbx_left_ii[i=1:length(boxes)-1],
        m[:flbz_front_left_bot_z][i] .>=
        [m[:min_cut_flbx_left_i][i], m[:min_cut_flbx_j_plus_left_i][i]]
    )

    #the difference between max and min is the cut
    @expression(m, cut_flbx_left_i[i=1:length(boxes)-1],
        m[:max_cut_flbx_left_i][i] - m[:min_cut_flbx_left_i][i]
    )
    @expression(m, cut_flbx_j_plus_left_i[i=1:length(boxes)-1],
        m[:max_cut_flbx_j_plus_left_i][i] - m[:min_cut_flbx_j_plus_left_i][i]
    )

    """
    @variable(m, cut_left[i=1:length(boxes)] >= 0)
    @constraint(m, c_cut_left[i=1:length(boxes)],
        m[:cut_left][i] .>=
        [m[:cut_flbx_left_i][i], m[:cut_flbx_j_plus_left_i][i]]
    )
    """

    # flbx + cover
    @variable(m, max_cut_flbx_ij_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flbx_ij_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flbx_ij_plus_right_ij][i, j] >=
        m[:guillotine_flbx_ij_plus][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flbx_ij_plus_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flbx_ij_plus_right_i[i=1:length(boxes)-1],
        m[:max_cut_flbx_ij_plus_right_i][i] .>=
        [m[:max_cut_flbx_ij_plus_right_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )
    @constraint(m, c_max_cut_flbx_ij_plus_right_ii[i=1:length(boxes)-1],
        m[:max_cut_flbx_ij_plus_right_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    @variable(m, max_cut_flbx_i_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flbx_i_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flbx_i_plus_right_ij][i, j] >=
        m[:guillotine_flbx_i_plus][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flbx_i_plus_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flbx_i_plus_right_i[i=1:length(boxes)-1],
        m[:max_cut_flbx_i_plus_right_i][i] .>=
        [m[:max_cut_flbx_i_plus_right_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )
    @constraint(m, c_max_cut_flbx_i_plus_right_ii[i=1:length(boxes)-1],
        m[:max_cut_flbx_i_plus_right_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    # MIN using two possibilities
    @variable(m, min_cut_flbx_i_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flbx_i_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flbx_i_plus_right_ij][i, j] .>=
        [m[:guillotine_flbx_i_plus][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flbx_i_plus][i, j]) * m[:flbz_front_left_bot_z][i]]

    )
    @variable(m, min_cut_flbx_ij_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flbx_ij_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flbx_ij_plus_right_ij][i, j] .>=
        [m[:guillotine_flbx_ij_plus][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flbx_ij_plus][i, j]) * m[:flbz_front_left_bot_z][i]]

    )

    #Find the lowest box of the j other boxes
    @variable(m, min_cut_flbx_i_plus_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flbx_i_plus_right_i[i=1:length(boxes)-1],
        [m[:min_cut_flbx_i_plus_right_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flbx_i_plus_right_i][i]
    )

    @variable(m, min_cut_flbx_ij_plus_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flbx_ij_plus_right_i[i=1:length(boxes)-1],
        [m[:min_cut_flbx_ij_plus_right_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flbx_ij_plus_right_i][i]
    )

    #Need to ensure that the min_guillotine cut must be lower than or equal to z_i
    @constraint(m, c_min_cut_flbx_cover_right_ii[i=1:length(boxes)-1],
        m[:flbz_front_left_bot_z][i] .>=
        [m[:min_cut_flbx_i_plus_right_i][i], m[:min_cut_flbx_ij_plus_right_i][i]]
    )

    #the difference between max and min is the cut
    @expression(m, cut_flbx_i_plus_right[i=1:length(boxes)-1],
        m[:max_cut_flbx_i_plus_right_i][i] - m[:min_cut_flbx_i_plus_right_i][i]
    )
    @expression(m, cut_flbx_ij_plus_right[i=1:length(boxes)-1],
        m[:max_cut_flbx_ij_plus_right_i][i] - m[:min_cut_flbx_ij_plus_right_i][i]
    )

    """
    @variable(m, cut_right[i=1:length(boxes)] >= 0)
    @constraint(m, c_cut_right[i=1:length(boxes)],
        m[:cut_right][i] .>=
        [m[:cut_flbx_i_plus_right][i], m[:cut_flbx_ij_plus_right][i]]
    )
    """

    # ===== Y AXIS =====

    # flby
    @variable(m, max_cut_flby_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flby_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flby_right_ij][i, j] >=
        m[:guillotine_flby][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flby_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flby_right_i[i=1:length(boxes)-1],
        m[:max_cut_flby_right_i][i] .>=
        [m[:max_cut_flby_right_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )

    @constraint(m, c_max_cut_flby_right_ii[i=1:length(boxes)-1],
        m[:max_cut_flby_right_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    @variable(m, max_cut_flby_j_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flby_j_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flby_j_plus_right_ij][i, j] >=
        m[:guillotine_flby_j_plus][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flby_j_plus_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flby_j_plus_right_i[i=1:length(boxes)-1],
        m[:max_cut_flby_j_plus_right_i][i] .>=
        [m[:max_cut_flby_j_plus_right_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )

    @constraint(m, c_max_cut_flby_j_plus_right_ii[i=1:length(boxes)-1],
        m[:max_cut_flby_j_plus_right_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    # MIN using two possibilities
    @variable(m, min_cut_flby_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flby_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flby_right_ij][i, j] .>=
        [m[:guillotine_flby][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flby][i, j]) * m[:flbz_front_left_bot_z][i]]

    )
    @variable(m, min_cut_flby_j_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flby_j_plus_right_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flby_j_plus_right_ij][i, j] .>=
        [m[:guillotine_flby_j_plus][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flby_j_plus][i, j]) * m[:flbz_front_left_bot_z][i]]

    )

    #Find the lowest box of the j other boxes
    @variable(m, min_cut_flby_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flby_right_i[i=1:length(boxes)-1],
        [m[:min_cut_flby_right_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flby_right_i][i]
    )

    @variable(m, min_cut_flby_j_plus_right_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flby_j_plus_right_i[i=1:length(boxes)-1],
        [m[:min_cut_flby_j_plus_right_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flby_j_plus_right_i][i]
    )

    #Need to ensure that the min_guillotine cut must be lower than or equal to z_i
    @constraint(m, c_min_cut_flby_cover_right_ii[i=1:length(boxes)-1],
        m[:flbz_front_left_bot_z][i] .>=
        [m[:min_cut_flby_right_i][i], m[:min_cut_flby_j_plus_right_i][i]]
    )

    #the difference between max and min is the cut
    @expression(m, cut_flby_right[i=1:length(boxes)-1],
        m[:max_cut_flby_right_i][i] - m[:min_cut_flby_right_i][i]
    )
    @expression(m, cut_flby_j_plus_right[i=1:length(boxes)-1],
        m[:max_cut_flby_j_plus_right_i][i] - m[:min_cut_flby_j_plus_right_i][i]
    )

    """
    @variable(m, cut_right[i=1:length(boxes)] >= 0)
    @constraint(m, c_cut_right[i=1:length(boxes)],
        m[:cut_right][i] .>=
        [m[:cut_flby_right][i], m[:cut_flby_j_plus_right][i]]
    )
    """

    # flby + cover
    @variable(m, max_cut_flby_ij_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flby_ij_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flby_ij_plus_left_ij][i, j] >=
        m[:guillotine_flby_ij_plus][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flby_ij_plus_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flby_ij_plus_left_i[i=1:length(boxes)-1],
        m[:max_cut_flby_ij_plus_left_i][i] .>=
        [m[:max_cut_flby_ij_plus_left_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )
    @constraint(m, c_max_cut_flby_ij_plus_left_ii[i=1:length(boxes)-1],
        m[:max_cut_flby_ij_plus_left_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    @variable(m, max_cut_flby_i_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_max_cut_flby_i_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:max_cut_flby_i_plus_left_ij][i, j] >=
        m[:guillotine_flby_i_plus][i, j] * (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j])
    )

    @variable(m, max_cut_flby_i_plus_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_max_cut_flby_i_plus_left_i[i=1:length(boxes)-1],
        m[:max_cut_flby_i_plus_left_i][i] .>=
        [m[:max_cut_flby_i_plus_left_ij][i, j] for j in 1:length(boxes)-1 if i != j]
    )
    @constraint(m, c_max_cut_flby_i_plus_left_ii[i=1:length(boxes)-1],
        m[:max_cut_flby_i_plus_left_i][i] >=
        m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]
    )

    @variable(m, min_cut_flby_i_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flby_i_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flby_i_plus_left_ij][i, j] .>=
        [m[:guillotine_flby_i_plus][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flby_i_plus][i, j]) * m[:flbz_front_left_bot_z][i]]

    )
    @variable(m, min_cut_flby_ij_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j] >= 0)
    @constraint(m, c_min_cut_flby_ij_plus_left_ij[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:min_cut_flby_ij_plus_left_ij][i, j] .>=
        [m[:guillotine_flby_ij_plus][i, j] * m[:flbz_front_left_bot_z][j], (1 - m[:guillotine_flby_ij_plus][i, j]) * m[:flbz_front_left_bot_z][i]]

    )

    #Find the lowest box of the j other boxes
    @variable(m, min_cut_flby_i_plus_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flby_i_plus_left_i[i=1:length(boxes)-1],
        [m[:min_cut_flby_i_plus_left_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flby_i_plus_left_i][i]
    )

    @variable(m, min_cut_flby_ij_plus_left_i[i=1:length(boxes)-1] >= 0)
    @constraint(m, c_min_cut_flby_ij_plus_left_i[i=1:length(boxes)-1],
        [m[:min_cut_flby_ij_plus_left_ij][i, j] for j in 1:length(boxes)-1 if i != j] .>=
        m[:min_cut_flby_ij_plus_left_i][i]
    )

    #Need to ensure that the min_guillotine cut must be lower than or equal to z_i
    @constraint(m, c_min_cut_flby_cover_left_ii[i=1:length(boxes)-1],
        m[:flbz_front_left_bot_z][i] .>=
        [m[:min_cut_flby_i_plus_left_i][i], m[:min_cut_flby_ij_plus_left_i][i]]
    )

    #the difference between max and min is the cut
    @expression(m, cut_flby_i_plus_left[i=1:length(boxes)-1],
        m[:max_cut_flby_i_plus_left_i][i] - m[:min_cut_flby_i_plus_left_i][i]
    )
    @expression(m, cut_flby_ij_plus_left[i=1:length(boxes)-1],
        m[:max_cut_flby_ij_plus_left_i][i] - m[:min_cut_flby_ij_plus_left_i][i]
    )

    """
    @variable(m, cut_left[i=1:length(boxes)] >= 0)
    @constraint(m, c_cut_left[i=1:length(boxes)],
        m[:cut_left][i] .>=
        [m[:cut_flby_i_plus_left][i], m[:cut_flby_ij_plus_left][i]]
    )

    @variable(m, max_cut >= 0)
    @constraints(m, begin
        m[:max_cut] .>= [m[:cut_left][i] for i in 1:length(boxes)]
        m[:max_cut] .>= [m[:cut_right][i] for i in 1:length(boxes)]
        m[:max_cut] .>= [m[:cut_left][i] for i in 1:length(boxes)]
        m[:max_cut] .>= [m[:cut_right][i] for i in 1:length(boxes)]
    end)
    """

    @expression(m, guillotine_interval_objective[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        #xleft
        0.5 * m[:max_cut_flbx_left_ij][i, j] +
        0.5 * m[:max_cut_flbx_left_i][i] +
        0.5 * m[:max_cut_flbx_j_plus_left_ij][i, j] +
        0.5 * m[:max_cut_flbx_j_plus_left_i][i] +
        10 * m[:min_cut_flbx_left_ij][i, j] +
        10 * m[:min_cut_flbx_j_plus_left_ij][i, j] +
        8 * (H - m[:min_cut_flbx_left_i][i]) +
        8 * (H - m[:min_cut_flbx_j_plus_left_i][i]) +
        #0.00000001 * m[:cut_left][i] +

        #xright
        0.5 * m[:max_cut_flbx_ij_plus_right_ij][i, j] +
        0.5 * m[:max_cut_flbx_ij_plus_right_i][i] +
        0.5 * m[:max_cut_flbx_i_plus_right_ij][i, j] +
        0.5 * m[:max_cut_flbx_i_plus_right_i][i] +
        10 * m[:min_cut_flbx_i_plus_right_ij][i, j] +
        10 * m[:min_cut_flbx_ij_plus_right_ij][i, j] +
        8 * (H - m[:min_cut_flbx_i_plus_right_i][i]) +
        8 * (H - m[:min_cut_flbx_ij_plus_right_i][i]) +
        #0.00000001 * m[:cut_right][i] +

        #yleft
        0.5 * m[:max_cut_flby_right_ij][i, j] +
        0.5 * m[:max_cut_flby_right_i][i] +
        0.5 * m[:max_cut_flby_j_plus_right_ij][i, j] +
        0.5 * m[:max_cut_flby_j_plus_right_i][i] +
        10 * m[:min_cut_flby_right_ij][i, j] +
        10 * m[:min_cut_flby_j_plus_right_ij][i, j] +
        8 * (H - m[:min_cut_flby_right_i][i]) +
        8 * (H - m[:min_cut_flby_j_plus_right_i][i]) +
        #0.00000001 * m[:cut_right][i] +

        #yleft
        0.5 * m[:max_cut_flby_ij_plus_left_ij][i, j] +
        0.5 * m[:max_cut_flby_ij_plus_left_i][i] +
        0.5 * m[:max_cut_flby_i_plus_left_ij][i, j] +
        0.5 * m[:max_cut_flby_i_plus_left_i][i] +
        10 * m[:min_cut_flby_i_plus_left_ij][i, j] +
        10 * m[:min_cut_flby_ij_plus_left_ij][i, j] +
        8 * (H - m[:min_cut_flby_i_plus_left_i][i]) +
        8 * (H - m[:min_cut_flby_ij_plus_left_i][i])
        #0.00000001 * m[:cut_left][i]

        #0.00001 * m[:max_cut]
    )
    return m
end

end


text:
using LinkedLists

function compare(array, list)
    # Iterate over items in container.
    for (ridx, d) in enumerate(list)
        @assert(d==array[ridx])
    end
end

function compare_iterator(array, list)
    # Iterate over indices of container, which are nodes.
    for (ridx, iter) in enumerate(eachindex(list))
        @assert(getindex(list, iter)==array[ridx])
    end
end

function walkthrough()
    l=SLinkedList{Int}()
    @assert(isempty(l))
    @assert(length(l)==0)
    vals=[3,2,7,7,9]

    prepend!(l, vals)
    println(l)
    println(vals)
    compare(vals, l)
    compare_iterator(vals, l)

    @assert(!isempty(l))
    show(l)

    pushfirst!(vals, 8)
    pushfirst!(l, 8)
    compare(vals, l)

    @assert(length(l)==6)

    @assert(indextoposition(findall(x -> x∈[2,3], l) ,l) == findall(x -> x∈[2,3], vals))
    @assert(indextoposition(indexin([2,3,12], l), l) == indexin([2,3,12], vals))

    @assert(unique(l)==unique(vals))

    @assert(getindex(l, lastindex(l))==9)

    @assert(!in(13, l))
    @assert(in(2, l))
    @assert(2 in l)
    @assert(eltype(l)==Int)
    @assert(first(l)==8)
    @assert(last(l)==9)
    show(l)
    push!(l, 11)
    println(l)
    println(vals)
    @assert(pop!(l)==11)
    @assert(popfirst!(l)==8)
    popfirst!(vals)

    two=findfirst(x -> x==2, l)
    insert!(l, two, 6)
    insert!(vals, 3, 6)
    compare(vals, l)
    deleteat!(l, two)
    deleteat!(vals, 2)
    compare(vals, l)
    splice!(l, findfirst(x -> x==9, l))
    deleteat!(vals, 5)
    compare(vals, l)
    splice!(l, findfirst(x -> x==6, l), 12)
    vals[2]=12
    compare(vals, l)
    append!(l, [23, 24])
    append!(vals, [23, 24])
    compare(vals, l)
    prepend!(l, [17, 16])
    prepend!(vals, [17, 16])
    println(l)
    println(vals)
    compare(vals, l)
end

walkthrough()


text:
# General application of plasticity rule on network

function plasticity_update!(t_now::Real,ntw::RecurrentNetwork)
  for pop in ntw.populations
    post = pop.state
    for (conn,pre) in zip(pop.connections,pop.pre_states)
      for plast in conn.plasticities
        plasticity_update!(t_now,pre,conn,post,plast)
      end
    end
  end
  return nothing
end

# bounds are structures too!
# because!

abstract type PlasticityBounds end
struct PlasticityBoundsNonnegative <: PlasticityBounds end
struct PlasticityBoundsLowHigh <: PlasticityBounds
  low::Float64
  high::Float64
end
function (plb::PlasticityBoundsNonnegative)(w::R,Δw::R) where R
  ret = w+Δw
  return max(eps(100.0),ret)
end
function (plb::PlasticityBoundsLowHigh)(w::R,Δw::R) where R
  ret = w+Δw
  return min(plb.high,max(plb.low,ret))
end

# utility function (ALWAYS follow post<-pre , never pre,post)
@inline function find_which_spiked(t_spike::Real,pspost::PopulationState,pspre::PopulationState)
  return find_which_spiked(t_spike,pspost.trains,pspre.trains)
end
function find_which_spiked(t_spike::R,trains_post::Vector{Vector{R}},trains_pre::Vector{Vector{R}}) where R
  is_t_now(train) = isempty(train) ? false : last(train) == t_spike
  k_post::Int64 = something(findfirst(is_t_now,trains_post),0)
  k_pre::Int64  = something(findfirst(is_t_now,trains_pre),0)
  return (k_post,k_pre)
end
function find_which_spiked(label_spike::Symbol,idx_spike::Int64,
    pspost::AbstractPopulationState,pspre::AbstractPopulationState)
  k_post = label_spike == pspost.label ? idx_spike : 0
  k_pre = label_spike == pspre.label ? idx_spike : 0
  return (k_post,k_pre)
end

# two different signatures, but unified base plasticity update

# signature for ExpKernel: plasticity_update!(tfire,label_fire,neufire,post,conn,pre,plast)
@inline function plasticity_update!(tfire::Real,label_spike::Symbol,
    neufire::Integer,
    pspre::AbstractPopulationState,conn::Connection,
    pspost::AbstractPopulationState,
    plast::PlasticityRule)
  k_post,k_pre = find_which_spiked(label_spike,neufire,pspost,pspre)
  return plasticity_update!(tfire,k_post,k_pre,pspre,conn,pspost,plast)
end

# signature for general Hawkes
@inline function plasticity_update!(t_spike::Real,
    pspost::PopulationState,conn::Connection,pspre::PopulationState,
    plast::PlasticityRule)
  k_post,k_pre = find_which_spiked(t_spike,pspost,pspre)
  return plasticity_update!(t_spike,k_post,k_pre,pspre,conn,pspost,plast)
end

# Pairwise plasticity
struct PairSTDP <: PlasticityRule
  Aplus::Float64
  Aminus::Float64
  pre_trace::Trace
  post_trace::Trace
  bounds::PlasticityBounds
  function PairSTDP(τplus,τminus,Aplus,Aminus,n_post,n_pre;
       plasticity_bounds=PlasticityBoundsNonnegative())
    pre_trace = Trace(τplus,n_pre)
    post_trace = Trace(τminus,n_post)
    new(Aplus,Aminus,pre_trace,post_trace,plasticity_bounds)
  end
end
function reset!(pl::PairSTDP)
  reset!(pl.pre_trace)
  reset!(pl.post_trace)
  return nothing
end

function plasticity_update!(t_spike::Real,k_post_spike::Integer,k_pre_spike::Integer,
    ::AbstractPopulationState,conn::Connection,::AbstractPopulationState,plast::PairSTDP)
  if iszero(k_pre_spike) && iszero(k_post_spike)
    return nothing
  end
  weights = conn.weights
  npost,npre = size(weights)
  # update all pre and post traces to t_now
  propagate!(t_spike,plast.pre_trace)
  propagate!(t_spike,plast.post_trace)
  # update the plasticity trace variables, if needed
  update_now!(plast.pre_trace,k_pre_spike)
  update_now!(plast.post_trace,k_post_spike)
  if !iszero(k_pre_spike)
    # k is presynaptic: go along rows of k_pre column
    for i in 1:npost
      wik = weights[i,k_pre_spike]
      if wik > 0
        Δw = plast.post_trace.val[i]*plast.Aminus
        weights[i,k_pre_spike] =  plast.bounds(wik,Δw)
      end
    end
  end
  if !iszero(k_post_spike)
    # k is postsynaptic: go along columns of k_post row
    for j in 1:npre
      wkj = weights[k_post_spike,j]
      if wkj > 0
        Δw = plast.pre_trace.val[j]*plast.Aplus
        weights[k_post_spike,j] =  plast.bounds(wkj,Δw)
      end
    end
  end
  return nothing
end

### Triplets rule
struct PlasticityTriplets <: PlasticityRule
  A2plus::Float64
  A3plus::Float64
  A2minus::Float64
  A3minus::Float64
  o1::Trace # pOst tau_minus
  o2::Trace # pOst tau_y
  r1::Trace # pRe  tau_plus
  r2::Trace # pRe  tau_y
  bounds::PlasticityBounds
end
function PlasticityTriplets(τplus::R,τminus::R,τx::R,τy::R,
    A2plus::R,A3plus::R,A2minus::R,A3minus::R,n_post::I,n_pre::I;
      plasticity_bounds=PlasticityBoundsNonnegative()) where {R,I}
  o1 = Trace(τminus,n_post,ForPlasticity())
  o2 = Trace(τy,n_post,ForPlasticity())
  r1 = Trace(τplus,n_pre,ForPlasticity())
  r2 = Trace(τx,n_pre,ForPlasticity())
  PlasticityTriplets(A2plus,A3plus,A2minus,A3minus,
    o1,o2,r1,r2,plasticity_bounds)
end
function reset!(pl::PlasticityTriplets)
  reset!.((pl.o1,pl.o2,pl.r1,pl.r2))
  return nothing
end

function plasticity_update!(t_spike::Real,k_post_spike::Integer,k_pre_spike::Integer,
    ::AbstractPopulationState,conn::Connection,::AbstractPopulationState,
    plast::PlasticityTriplets)
  if iszero(k_pre_spike) && iszero(k_post_spike)
    return nothing
  end
  # update all pre and post traces to t_now
  propagate!.(t_spike,(plast.o1,plast.o2,plast.r1,plast.r2))
  # update the plasticity trace variables for o1 and r1
  if !iszero(k_post_spike)
    update_now!(plast.o1,k_post_spike)
  end
  if !iszero(k_pre_spike)
    update_now!(plast.r1,k_pre_spike)
  end
  # update synapses
  weights=conn.weights
  npost,npre = size(weights)
  if !iszero(k_pre_spike)
    # k is presynaptic: move vertically along k_pre column
    for i_post in 1:npost
      wik = weights[i_post,k_pre_spike]
      if wik > 0
        Δw = -plast.o1.val[i_post]*(plast.A2minus+plast.A3minus*plast.r2.val[k_pre_spike])
        weights[i_post,k_pre_spike] =  plast.bounds(wik,Δw)
      end
    end
  end
  if !iszero(k_post_spike)
    # k is postsynaptic: move horizontally along k_post row
    for j_pre in 1:npre
      wkj = weights[k_post_spike,j_pre]
      if wkj > 0
        Δw = plast.r1.val[j_pre]*(plast.A2plus+plast.A3plus*plast.o2.val[k_post_spike])
        weights[k_post_spike,j_pre] =  plast.bounds(wkj,Δw)
      end
    end
  end
  # update the plasticity trace variables for o2 and r2
  if !iszero(k_post_spike)
    update_now!(plast.o2,k_post_spike)
  end
  if !iszero(k_pre_spike)
    update_now!(plast.r2,k_pre_spike)
  end
  return nothing
end

# Inhibitory stabilization, Vogels-Sprekeler 2011

# Warning : this formulation likely assumes all positive weights !

struct PlasticityInhibitory <: PlasticityRule
  τ::Float64
  η::Float64
  α::Float64
  o::Trace # pOst
  r::Trace # pRe
  bounds::PlasticityBounds
  function PlasticityInhibitory(τ,η,n_post,n_pre;
      r_target=5.0,
      plasticity_bounds::PlasticityBounds=PlasticityBoundsNonnegative())
    α = 2*r_target*τ
    o = Trace(τ,n_post,ForPlasticity())
    r = Trace(τ,n_pre,ForPlasticity())
    new(τ,η,α,o,r,plasticity_bounds)
  end
end
function reset!(pl::PlasticityInhibitory)
  reset!.((pl.o,pl.r))
  return nothing
end

function plasticity_update!(t_spike::Real,k_post_spike::Integer,k_pre_spike::Integer,
    ::AbstractPopulationState,conn::Connection,::AbstractPopulationState,
    plast::PlasticityInhibitory)
  if iszero(k_pre_spike) && iszero(k_post_spike)
    return nothing
  end
  weights = conn.weights
  # update all pre and post traces to t_now
  propagate!.(t_spike,(plast.o,plast.r))
  # spike update, pre or post
  if !iszero(k_post_spike)
    update_now!(plast.o,k_post_spike)
  end
  if !iszero(k_pre_spike)
    update_now!(plast.r,k_pre_spike)
  end
  if !iszero(k_pre_spike)
    # k is presynaptic: move vertically along k_pre column
    for i_post in 1:npost
      wik = weights[i_post,k_pre_spike]
      if !iszero(wik)
        Δw = plast.η*(plast.o[i_post]-plast.α)
        weights[i_post,k_pre_spike] = plast.bounds(wik,Δw)
      end
    end
  end
  if !iszero(k_post_spike)
    # k is presynaptic: move vertically along k_pre column
    for j_pre in 1:npre
      wkj = weights[k_post_spike,j_pre]
      if !iszero(wkj)
        Δw = plast.η*plast.r[j_pre]
        weights[k_post_spike,j_pre] = plast.bounds(wkj,Δw)
      end
    end
  end
  return nothing
end