text:
module _Arctan

using Brainstorm._DataStructure: PreAllocatedStack
using Brainstorm._Math._GMP: set!, add!, mul!, div!, pow!

const guard_terms = 5
const guard_digits = 5

preallocatedstack(size) = PreAllocatedStack(0, [(big(0), big(0), big(0)) for i in 1:size])

arctan(digits::Integer, reciprocal_x::Integer) = begin
    p, q, o = _arctan(digits, false, reciprocal_x)

    set!(o, 10)
    pow!(o, digits) # one

    mul!(p, o)
    div!(p, q)
end

arctanh(digits::Integer, reciprocal_x::Integer) = begin
    p, q, o = _arctan(digits, true, reciprocal_x)

    set!(o, 10)
    pow!(o, digits) # one

    mul!(p, o)
    div!(p, q)
end

arctansum(digits::Integer, rx, coeff) = begin
    p, q, o = _arctansum(digits, false, rx, coeff)

    set!(o, 10)
    pow!(o, digits) # one

    mul!(p, o)
    div!(p, q)
end

arctanhsum(digits::Integer, rx, coeff) = begin
    p, q, o = _arctansum(digits, true, rx, coeff)

    set!(o, 10)
    pow!(o, digits) # one

    mul!(p, o)
    div!(p, q)
end

_arctan(digits, hyperbolic, rx) = begin
    n = terms(digits, rx)

    depth = 1
    while (1 << depth) <= n; depth += 1; end
    s = preallocatedstack(depth + 1)

    _arctan(hyperbolic, rx, n, s)
end

_arctansum(digits, hyperbolic, rx, f) = begin
    n = [terms(digits + guard_digits, x) for x in rx]

    depth = 1
    maxn = maximum(n)
    while (1 << depth) <= maxn; depth += 1; end
    s = preallocatedstack(depth + 2)

    for i in 1:length(rx)
        p, q, _ = _arctan(hyperbolic, rx[i], n[i], s)
        f[i] != 1 && mul!(p, f[i])

        if i > 1
            pop!(s)
            sp, sq = s[1]

            # sp = sp * q + p * sq
            mul!(p, sq)
            mul!(sp, q)
            add!(sp, p)

            # sq = sq * q
            mul!(sq, q)
        end
        push!(s)
    end
    pop!(s)

    p, q, t = s[1]
    p, q, t
end

# How many terms to compute
#
#         x^(2k+1) > 10^digits
# => (2k+1) log(x) > digits * log(10)
# =>        2k + 1 > digits * log(10) / log(x)
# =>            2k > digits * log(10) / log(x)
# =>             k > [digits * log(10) / log(x)] / 2
#
terms(digits::T, rx) where {T<:Integer} = trunc(T, (digits / 2) * log(rx, 10) + guard_terms)

_arctan(hyperbolic, rx, n, s) = begin
    # calculate P(0, N), Q(0, N) and T(0, N)
    binarysplitting!(s, hyperbolic, rx^2, zero(n), n)

    complete!(s, rx)
end

binarysplitting!(s, hyperbolic, sqrx, a, b) = begin
    # directly compute P(a, a + 1), Q(a, a + 1) and T(a, a + 1)
    (b - a) == 1 && (compute!(s, hyperbolic, sqrx, a, b); return)

    # recursively compute P(a, b), Q(a, b) and T(a, b)

    # m is the midpoint of a and b
    m = (a + b) >>> 1

    # recursively calculate P(a, m), Q(a, m) and T(a, m)
    binarysplitting!(s, hyperbolic, sqrx, a, m)

    # recursively calculate P(m, b), Q(m, b) and T(m, b)
    push!(s)
    binarysplitting!(s, hyperbolic, sqrx, m, b)
    pop!(s)

    # now combine
    combine!(s)
    return
end

compute!(s, hyperbolic, sqrx, a, b) = begin
    p, q, t = s[1]

    y = 2a + 3

    set!(p, hyperbolic || isodd(a) ? 1 : -1)
    set!(q, y)
    mul!(q, sqrx)
    set!(t, y)
end

combine!(s) = begin
    p₁, q₁, t₁ = s[1]
    p₂, q₂, t₂ = s[2]

    # p = q₂ * p₁ + t₁ * p₂
    mul!(p₁, q₂)
    mul!(p₂, t₁)
    add!(p₁, p₂)

    # q = q₁ * q₂
    mul!(q₁, q₂)

    # t = t₁ * t₂
    mul!(t₁, t₂)
end

complete!(s, rx) = begin
    p, q, t = s[1]

    add!(p, q)
    mul!(q, rx)

    p, q, t
end

end


text:
@deprecate lmm(f::Formula, fr::AbstractDataFrame) LinearMixedModel(f, fr)
@deprecate glmm(f::Formula, fr::AbstractDataFrame, d::Distribution) GeneralizedLinearMixedModel(f, fr, d)
@deprecate glmm(f::Formula, fr::AbstractDataFrame, d::Distribution, l::Link) GeneralizedLinearMixedModel(f, fr, d, l)

ext: jl


text:
using LinearAlgebra

"""
Truncated SVD to find the matrix of rank k that is the closest to A

# Arguments
- `A::Array{Float64,2}`
- `k::Array{Float64,1}`

"""
function approx_matrix_svd(A,k)
    (m,n) = size(A)
    F = svd(A)
    A_SVD = F.U[1:m,1:k] * Diagonal(F.S)[1:k,1:k] * F.V'[1:k,1:n]
    return A_SVD
end


text:
using Statistics
using Logging
using Printf
using DataDeps
using JLD2

using Oceananigans
using Oceananigans.Units

using Oceananigans.Diagnostics: accurate_cell_advection_timescale

Logging.global_logger(OceananigansLogger())

#####
##### Progress monitor
#####

mutable struct Progress
    interval_start_time :: Float64
end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9

    @info @sprintf("Time: %s, iteration: %d, extrema(h): (min=%.2e, max=%.2e), CFL: %.2e, wall time: %s",
                prettytime(sim.model.clock.time),
                sim.model.clock.iteration,
                minimum(sim.model.tracers.h),
                maximum(sim.model.tracers.h),
                sim.parameters.cfl(sim.model),
                prettytime(wall_time))

    p.interval_start_time = time_ns()

    return nothing
end

#####
##### Script starts here
#####

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

dd = DataDep("cubed_sphere_32_grid",
    "Conformal cubed sphere grid with 32×32 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cubed_sphere_32_grid.jld2",
    "b1dafe4f9142c59a2166458a2def743cd45b20a4ed3a1ae84ad3a530e1eff538" # sha256sum
)

DataDeps.register(dd)

cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"

# Central (λ, φ) for each face of the cubed sphere
central_longitude = (0, 90,  0, 180, -90,   0)
central_latitude  = (0,  0, 90,   0,   0, -90)

"""
    cubed_sphere_tracer_advection(; face_number, α)

Run a tracer advection experiment that initializes a cosine bell on face `face_number`
and advects it around the sphere over 12 days. `α` is the angle between the axis of
solid body rotation and the polar axis (degrees).
"""
function cubed_sphere_tracer_advection(; face_number, α)

    grid = ConformalCubedSphereGrid(cs32_filepath, Nz=1, z=(-1, 0))

    ## Prescribed velocities and initial condition according to Williamson et al. (1992) §3.1

    period = 12days  # Time to make a full rotation (s)
    R = grid.faces[1].radius  # radius of the sphere (m)
    u₀ = 2π*R / period  # advecting velocity (m/s)

    # U(λ, φ, z) = u₀ * (cosd(φ) * cosd(α) + sind(φ) * cosd(λ) * sind(α))
    # V(λ, φ, z) = - u₀ * sind(λ) * sind(α)

    Ψ(λ, φ, z) = - R * u₀ * (sind(φ) * cosd(α) - cosd(λ) * cosd(φ) * sind(α))

    Ψᶠᶠᶜ = Field(Face, Face,   Center, CPU(), grid)
    Uᶠᶜᶜ = Field(Face, Center, Center, CPU(), grid)
    Vᶜᶠᶜ = Field(Center, Face, Center, CPU(), grid)
    Wᶜᶜᶠ = Field(Center, Center, Face, CPU(), grid)  # So we can use CFL

    for (f, grid_face) in enumerate(grid.faces)
        for i in 1:grid_face.Nx+1, j in 1:grid_face.Ny+1
            Ψᶠᶠᶜ.data.faces[f][i, j, 1] = Ψ(grid_face.λᶠᶠᵃ[i, j], grid_face.φᶠᶠᵃ[i, j], 0)
        end
    end

    for (f, grid_face) in enumerate(grid.faces)
        Ψᶠᶠᶜ_face = Ψᶠᶠᶜ.data.faces[f]
        Uᶠᶜᶜ_face = Uᶠᶜᶜ.data.faces[f]
        Vᶜᶠᶜ_face = Vᶜᶠᶜ.data.faces[f]

        for i in 1:grid_face.Nx+1, j in 1:grid_face.Ny
            Uᶠᶜᶜ_face[i, j, 1] = (Ψᶠᶠᶜ_face[i, j, 1] - Ψᶠᶠᶜ_face[i, j+1, 1]) / grid.faces[f].Δyᶠᶜᵃ[i, j]
        end

        for i in 1:grid_face.Nx, j in 1:grid_face.Ny+1
            Vᶜᶠᶜ_face[i, j, 1] = (Ψᶠᶠᶜ_face[i+1, j, 1] - Ψᶠᶠᶜ_face[i, j, 1]) / grid.faces[f].Δxᶜᶠᵃ[i, j]
        end
    end

    ## Model setup

    model = HydrostaticFreeSurfaceModel(
                      grid = grid,
        momentum_advection = nothing,
                   tracers = :h,
                velocities = PrescribedVelocityFields(u=Uᶠᶜᶜ, v=Vᶜᶠᶜ, w=Wᶜᶜᶠ),
              free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1),
                  coriolis = nothing,
                   closure = nothing,
                  buoyancy = nothing
    )

    ## Cosine bell initial condition according to Williamson et al. (1992) §3.1

    h₀ = 1000
    λ₀ = central_longitude[face_number]
    φ₀ = central_longitude[face_number]

    # Great circle distance between (λ, φ) and the center of the cosine bell (λ₀, φ₀)
    # using the haversine formula
    r(λ, φ) = 2R * asin(√(sind((φ - φ₀) / 2)^2 + cosd(φ) * cosd(φ₀) * sind((λ - λ₀) / 2)^2))

    cosine_bell(λ, φ, z) = r(λ, φ) < R ? h₀/2 * (1 + cos(π * r(λ, φ) / R)) : 0

    Oceananigans.set!(model, h=cosine_bell)

    ## Simulation setup

    Δt = 10minutes

    cfl = CFL(Δt, accurate_cell_advection_timescale)

    simulation = Simulation(model,
                        Δt = Δt,
                 stop_time = 1period,
        iteration_interval = 1,
                  progress = Progress(time_ns()),
                parameters = (; cfl)
    )

    outputs = (u=Uᶠᶜᶜ, v=Vᶜᶠᶜ, h=model.tracers.h)

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, outputs,
            schedule = TimeInterval(1hour),
            filename = "cubed_sphere_tracer_advection_face$(face_number)_alpha$α",
            overwrite_existing = true)

    run!(simulation)

    return simulation
end

#####
##### Run all the experiments!
#####

include("animate_on_map_projection.jl")

function run_cubed_sphere_tracer_advection_validation()

    αs = (0, 45, 90)

    for f in 1:6, α in αs
        cubed_sphere_tracer_advection(face_number=f, α=α)
    end

    projections = [
        ccrs.NearsidePerspective(central_longitude=0,   central_latitude=30),
        ccrs.NearsidePerspective(central_longitude=180, central_latitude=-30)
    ]

    for f in 1:6, α in αs
        animate_tracer_advection(face_number=f, α=α, projections=projections)
    end
end


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

"""
    associate_attribute_group(application, attribute_group)
    associate_attribute_group(application, attribute_group, params::Dict{String,<:Any})

Associates an attribute group with an application to augment the application's metadata
with the group's attributes. This feature enables applications to be described with
user-defined details that are machine-readable, such as third-party integrations.

# Arguments
- `application`: The name or ID of the application.
- `attribute_group`: The name or ID of the attribute group that holds the attributes to
  describe the application.

"""
function associate_attribute_group(
    application, attributeGroup; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "PUT",
        "/applications/$(application)/attribute-groups/$(attributeGroup)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function associate_attribute_group(
    application,
    attributeGroup,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "PUT",
        "/applications/$(application)/attribute-groups/$(attributeGroup)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    associate_resource(application, resource, resource_type)
    associate_resource(application, resource, resource_type, params::Dict{String,<:Any})

Associates a resource with an application. Both the resource and the application can be
specified either by ID or name.

# Arguments
- `application`: The name or ID of the application.
- `resource`: The name or ID of the resource of which the application will be associated.
- `resource_type`: The type of resource of which the application will be associated.

"""
function associate_resource(
    application, resource, resourceType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "PUT",
        "/applications/$(application)/resources/$(resourceType)/$(resource)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function associate_resource(
    application,
    resource,
    resourceType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "PUT",
        "/applications/$(application)/resources/$(resourceType)/$(resource)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Creates a new application that is the top-level node in a hierarchy of related cloud
resource abstractions.

# Arguments
- `client_token`: A unique identifier that you provide to ensure idempotency. If you retry
  a request that completed successfully using the same client token and the same parameters,
  the retry succeeds without performing any further actions. If you retry a successful
  request using the same client token, but one or more of the parameters are different, the
  retry fails.
- `name`: The name of the application. The name must be unique in the region in which you
  are creating the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: The description of the application.
- `"tags"`: Key-value pairs you can use to associate with the application.
"""
function create_application(
    clientToken, name; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "POST",
        "/applications",
        Dict{String,Any}("clientToken" => clientToken, "name" => name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_application(
    clientToken,
    name,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "POST",
        "/applications",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("clientToken" => clientToken, "name" => name),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_attribute_group(attributes, client_token, name)
    create_attribute_group(attributes, client_token, name, params::Dict{String,<:Any})

Creates a new attribute group as a container for user-defined attributes. This feature
enables users to have full control over their cloud application's metadata in a rich
machine-readable format to facilitate integration with automated workflows and third-party
tools.

# Arguments
- `attributes`: A JSON string in the form of nested key-value pairs that represent the
  attributes in the group and describes an application and its components.
- `client_token`: A unique identifier that you provide to ensure idempotency. If you retry
  a request that completed successfully using the same client token and the same parameters,
  the retry succeeds without performing any further actions. If you retry a successful
  request using the same client token, but one or more of the parameters are different, the
  retry fails.
- `name`: The name of the attribute group.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: The description of the attribute group that the user provides.
- `"tags"`: Key-value pairs you can use to associate with the attribute group.
"""
function create_attribute_group(
    attributes, clientToken, name; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "POST",
        "/attribute-groups",
        Dict{String,Any}(
            "attributes" => attributes, "clientToken" => clientToken, "name" => name
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_attribute_group(
    attributes,
    clientToken,
    name,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "POST",
        "/attribute-groups",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "attributes" => attributes, "clientToken" => clientToken, "name" => name
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Deletes an application that is specified either by its application ID or name. All
associated attribute groups and resources must be disassociated from it before deleting an
application.

# Arguments
- `application`: The name or ID of the application.

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

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

Deletes an attribute group, specified either by its attribute group ID or name.

# Arguments
- `attribute_group`: The name or ID of the attribute group that holds the attributes to
  describe the application.

"""
function delete_attribute_group(
    attributeGroup; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "DELETE",
        "/attribute-groups/$(attributeGroup)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_attribute_group(
    attributeGroup,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "DELETE",
        "/attribute-groups/$(attributeGroup)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    disassociate_attribute_group(application, attribute_group)
    disassociate_attribute_group(application, attribute_group, params::Dict{String,<:Any})

Disassociates an attribute group from an application to remove the extra attributes
contained in the attribute group from the application's metadata. This operation reverts
AssociateAttributeGroup.

# Arguments
- `application`: The name or ID of the application.
- `attribute_group`: The name or ID of the attribute group that holds the attributes to
  describe the application.

"""
function disassociate_attribute_group(
    application, attributeGroup; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "DELETE",
        "/applications/$(application)/attribute-groups/$(attributeGroup)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function disassociate_attribute_group(
    application,
    attributeGroup,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "DELETE",
        "/applications/$(application)/attribute-groups/$(attributeGroup)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    disassociate_resource(application, resource, resource_type)
    disassociate_resource(application, resource, resource_type, params::Dict{String,<:Any})

Disassociates a resource from application. Both the resource and the application can be
specified either by ID or name.

# Arguments
- `application`: The name or ID of the application.
- `resource`: The name or ID of the resource.
- `resource_type`: The type of the resource that is being disassociated.

"""
function disassociate_resource(
    application, resource, resourceType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "DELETE",
        "/applications/$(application)/resources/$(resourceType)/$(resource)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function disassociate_resource(
    application,
    resource,
    resourceType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "DELETE",
        "/applications/$(application)/resources/$(resourceType)/$(resource)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves metadata information about one of your applications. The application can be
specified either by its unique ID or by its name (which is unique within one account in one
region at a given point in time). Specify by ID in automated workflows if you want to make
sure that the exact same application is returned or a ResourceNotFoundException is thrown,
avoiding the ABA addressing problem.

# Arguments
- `application`: The name or ID of the application.

"""
function get_application(application; aws_config::AbstractAWSConfig=global_aws_config())
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_application(
    application,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_associated_resource(application, resource, resource_type)
    get_associated_resource(application, resource, resource_type, params::Dict{String,<:Any})

Gets the resource associated with the application.

# Arguments
- `application`: The name or ID of the application.
- `resource`: The name or ID of the resource associated with the application.
- `resource_type`: The type of resource associated with the application.

"""
function get_associated_resource(
    application, resource, resourceType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)/resources/$(resourceType)/$(resource)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_associated_resource(
    application,
    resource,
    resourceType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)/resources/$(resourceType)/$(resource)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves an attribute group, either by its name or its ID. The attribute group can be
specified either by its unique ID or by its name.

# Arguments
- `attribute_group`: The name or ID of the attribute group that holds the attributes to
  describe the application.

"""
function get_attribute_group(
    attributeGroup; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/attribute-groups/$(attributeGroup)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_attribute_group(
    attributeGroup,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "GET",
        "/attribute-groups/$(attributeGroup)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves a list of all of your applications. Results are paginated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The upper bound of the number of results to return (cannot exceed 25). If
  this parameter is omitted, it defaults to 25. This value is optional.
- `"nextToken"`: The token to use to get the next page of results after a previous API
  call.
"""
function list_applications(; aws_config::AbstractAWSConfig=global_aws_config())
    return service_catalog_appregistry(
        "GET", "/applications"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_applications(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/applications",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists all attribute groups that are associated with specified application. Results are
paginated.

# Arguments
- `application`: The name or ID of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The upper bound of the number of results to return (cannot exceed 25). If
  this parameter is omitted, it defaults to 25. This value is optional.
- `"nextToken"`: The token to use to get the next page of results after a previous API
  call.
"""
function list_associated_attribute_groups(
    application; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)/attribute-groups";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_associated_attribute_groups(
    application,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)/attribute-groups",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists all resources that are associated with specified application. Results are paginated.

# Arguments
- `application`: The name or ID of the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The upper bound of the number of results to return (cannot exceed 25). If
  this parameter is omitted, it defaults to 25. This value is optional.
- `"nextToken"`: The token to use to get the next page of results after a previous API
  call.
"""
function list_associated_resources(
    application; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)/resources";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_associated_resources(
    application,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "GET",
        "/applications/$(application)/resources",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists all attribute groups which you have access to. Results are paginated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The upper bound of the number of results to return (cannot exceed 25). If
  this parameter is omitted, it defaults to 25. This value is optional.
- `"nextToken"`: The token to use to get the next page of results after a previous API
  call.
"""
function list_attribute_groups(; aws_config::AbstractAWSConfig=global_aws_config())
    return service_catalog_appregistry(
        "GET", "/attribute-groups"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_attribute_groups(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/attribute-groups",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists all of the tags on the resource.

# Arguments
- `resource_arn`: The Amazon resource name (ARN) that specifies the resource.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    sync_resource(resource, resource_type)
    sync_resource(resource, resource_type, params::Dict{String,<:Any})

Syncs the resource with current AppRegistry records. Specifically, the resource’s
AppRegistry system tags sync with its associated application. We remove the resource's
AppRegistry system tags if it does not associate with the application. The caller must have
permissions to read and update the resource.

# Arguments
- `resource`: An entity you can work with and specify with a name or ID. Examples include
  an Amazon EC2 instance, an Amazon Web Services CloudFormation stack, or an Amazon S3 bucket.
- `resource_type`: The type of resource of which the application will be associated.

"""
function sync_resource(
    resource, resourceType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "POST",
        "/sync/$(resourceType)/$(resource)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function sync_resource(
    resource,
    resourceType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "POST",
        "/sync/$(resourceType)/$(resource)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Assigns one or more tags (key-value pairs) to the specified resource. Each tag consists of
a key and an optional value. If a tag with the same key is already associated with the
resource, this action updates its value. This operation returns an empty response if the
call was successful.

# Arguments
- `resource_arn`: The Amazon resource name (ARN) that specifies the resource.
- `tags`: The new or modified tags for the resource.

"""
function tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config())
    return service_catalog_appregistry(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tags" => tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    resourceArn,
    tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tags" => tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Removes tags from a resource. This operation returns an empty response if the call was
successful.

# Arguments
- `resource_arn`: The Amazon resource name (ARN) that specifies the resource.
- `tag_keys`: A list of the tag keys to remove from the specified resource.

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Updates an existing application with new attributes.

# Arguments
- `application`: The name or ID of the application that will be updated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: The new description of the application.
- `"name"`: The new name of the application. The name must be unique in the region in which
  you are updating the application.
"""
function update_application(application; aws_config::AbstractAWSConfig=global_aws_config())
    return service_catalog_appregistry(
        "PATCH",
        "/applications/$(application)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_application(
    application,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "PATCH",
        "/applications/$(application)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Updates an existing attribute group with new details.

# Arguments
- `attribute_group`: The name or ID of the attribute group that holds the attributes to
  describe the application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"attributes"`: A JSON string in the form of nested key-value pairs that represent the
  attributes in the group and describes an application and its components.
- `"description"`: The description of the attribute group that the user provides.
- `"name"`: The new name of the attribute group. The name must be unique in the region in
  which you are updating the attribute group.
"""
function update_attribute_group(
    attributeGroup; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_catalog_appregistry(
        "PATCH",
        "/attribute-groups/$(attributeGroup)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_attribute_group(
    attributeGroup,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_catalog_appregistry(
        "PATCH",
        "/attribute-groups/$(attributeGroup)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
module GeneRegulatoryNetwork
export pcregress

using LinearAlgebra, Statistics, Arpack, Clustering
# using CommunityDetection   # https://github.com/afternone/CommunityDetection.jl

function pcregress(X::AbstractMatrix{T}) where T<:Real
  n=size(X,2)
  A=1.0 .-Matrix(I,n,n)
  for k in 1:n
      y=X[:,k]
      𝒳=X[:,1:end.≠k]
      _,ϕ=Arpack.eigs(𝒳'𝒳,nev=3,which=:LM)
      s=𝒳*ϕ
      s ./=(norm.(s[:,i] for i=1:size(s,2)).^2)'
      b=sum(y.*s,dims=1)
      𝒷=ϕ*b'
      A[k,A[k,:].==1.0]=𝒷
  end
  return A
end

function affprop(A::AbstractMatrix{T}, k) where T<:AbstractFloat
  # check arguments
  n = size(A, 1)
  size(A, 2) == n || throw(ArgumentError("A must be a square matrix ($(size(A)) given)."))
  r=Clustering.affprop(A, k);
  return r
end

end


text:
"""
    module Simulate_TwoCollidingBalls

Simulates a ball that is first sliding, then rolling and then colliding
with another ball on a table.
"""
module Simulate_TwoRollingAndCollidingBalls

using Modia3D
import Modia3D.ModiaMath

vmatSolids   = Material(color="Red"       , transparency=0.0)    # material of solids
vmatTable    = Material(color="Green"     , transparency=0.1)    # material of table

LxTable  = 3.0
LyTable  = 0.5
LzTable  = 0.02
diameter = 0.06

solidBall  = Solid(Modia3D.SolidSphere(diameter), "BilliardBall", vmatSolids ; contactMaterial="BilliardBall")
solidTable = Solid(SolidBox(LxTable, LyTable, LzTable) , "DryWood", vmatTable; contactMaterial="BilliardTable")

@assembly TwoCollidingBalls() begin
  world = Object3D()
  table = Object3D(world, solidTable, fixed=true , r=[1.5, 0.0, -LzTable/2])
  ball1 = Object3D(world, solidBall , fixed=false, r=[0.2, 0.0, diameter/2], v_start=[3.0, 0.0, 0.0], visualizeFrame=true)
  ball2 = Object3D(world, solidBall , fixed=false, r=[1.5, 0.0, diameter/2], visualizeFrame=true)
end

gravField         = UniformGravityField(g=9.81, n=[0,0,-1])
twoCollidingBalls = TwoCollidingBalls(sceneOptions=SceneOptions(gravityField=gravField, visualizeFrames=false,
                                                                defaultFrameLength=0.1, enableContactDetection=true))

model = SimulationModel( twoCollidingBalls )

result = ModiaMath.simulate!(model; stopTime=1.0, tolerance=1e-8, log=false)

ModiaMath.plot(result, [ ("ball1.r[1]", "ball2.r[1]"),
                         ("ball1.v[1]", "ball2.v[1]"),
                         ("ball1.w[2]", "ball2.w[2]")])

println("... success of examples/collisions/Simulate_TwoCollidingBalls.jl!")
end


text:
# Only run coverage from linux Julia 0.6 build on travis.
get(ENV, "TRAVIS_OS_NAME", "")       == "linux"   || exit()
get(ENV, "TRAVIS_JULIA_VERSION", "") == "0.6" || exit()

Pkg.add("Coverage")
using Coverage

cd(joinpath(dirname(@__FILE__), "..")) do
    processed = process_folder()
    Coveralls.submit(processed)
    Codecov.submit(processed)
end


text:
# The implementations are inspired by MatrixMarket.jl
# https://github.com/JuliaSparse/MatrixMarket.jl
# The MatrixMarket.jl package is licensed under the MIT Expat License:
# Copyright (c) 2013: Viral B. Shah.

# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:

# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

"""julia
    MatrixDepot

Give access to a wealth of sample and test matrices and accompanying data.
A set of matrices is generated locally (with arguments controlling the special case).
Another set is loaded from one of the publicly accessible matrix collections
`SuiteSparse Matrix Collection` (formerly `University of Florida Matrix Collection`)
and the `Matrix Market Collection`.

Access is like

    using MatrixDepot

    A = matrixdepot("hilb", 10) # locally generated hilbert matrix dimensions (10,10)

    A = matrixdepot("HB/1138_bus")     # named matrix of the SuiteSparse Collection
    A = matrixdepot(sp(1))             # same matrix using numerical id
    A = matrixdepot("Harwell*/*/1138_bus") # matrix from the Matrix Market Collection

    md = mdopen("*/bfly")   # named matrix with some extra data
    A = md.A
    co = md.coord
    tx = md("Gname_10.txt")

    md = mdopen("gravity", 10, false) # localy generated example with rhs and solution
    A = md.A
    b = md.b
    x = md.x

###### commands:
    mdinfo, listdir, listgroups, matrixdepot, mdopen, listdata, mdlist,
    metasymbols, @addgroup, @modifygroup, @rmgroup.
###### selector patterns:
    strings, string-patterns (using "*", "?", "/", "**"), regular expressions: for names
    builtin(42), user(3,5), sp(10:11,6,2833), mm(1): to access by integer id
###### predicate patterns:
    isboolean, isinteger, isreal, iscomplex
    isgeneral, issymmetric, ishermitian, isskew
    isbuiltin, isuser, islocal, isremote, isloaded, isunloaded
    keyword(string expression), logical, hasdata(symbol), @pred(expression)

    see also: "logical" for logical combinations of all kinds of patterns.
"""
module MatrixDepot
using LinearAlgebra, SparseArrays, SuiteSparse, Serialization
using CodecZlib
import Base: show

export matrixdepot
export listnames, listdir, listdata, listgroups, mdlist, mdinfo, metasymbols, mdopen
export @addgroup, @rmgroup, @modifygroup

# exports for predicate functions in `logical.jl`
export builtin, user, sp, mm, logical
export isgeneral, issymmetric, isskew, ishermitian
export iscomplex, isreal, isinteger, isboolean
export isremote, islocal, isloaded, isunloaded, isbuiltin, isuser
export @pred, keyword, hasdata

# The following functions are re-used as predicate functions / logical operators
import Base: isreal, isinteger
import LinearAlgebra: issymmetric, ishermitian
import SparseArrays: issparse
import Base: &, |, *, ~

include("types.jl")         # common data type definitions
include("higham.jl")        # test matrices
include("regu.jl")          # regularization test problem
include("graph.jl")         # adjacency matrices for graphs
include("data.jl")          # global variables and matrix data
include("common.jl")        # main functions
include("logical.jl")       # operations on patterns and predicates
include("download.jl")      # download data from the UF and MM sparse matrix collection
include("datareader.jl")    # read matrix data from local storage
include("matrixmarket.jl")  # read matrix data from local storage
include("markdown.jl")      # construct MD objects

function init(;ignoredb::Bool=false)
    GROUP = "group.jl"
    GENERATOR = "generator.jl"

    if !isdir(DATA_DIR)
        mkpath(DATA_DIR)
    end

    if !isdir(MY_DEPOT_DIR)
        mkpath(MY_DEPOT_DIR)
        open(joinpath(MY_DEPOT_DIR, GROUP), "w") do f
            write(f, "usermatrixclass = Dict(\n);")
        end
        open(joinpath(MY_DEPOT_DIR, GENERATOR), "w") do f
            write(f, "# include your matrix generators below \n")
        end
        println("created dir $MY_DEPOT_DIR")
    end

    for file in readdir(MY_DEPOT_DIR)
        if endswith(file, ".jl") && file != GENERATOR
            println("include $file for user defined matrix generators")
            include(joinpath(MY_DEPOT_DIR, file))
        end
    end
    include(joinpath(MY_DEPOT_DIR, GENERATOR))
    println("verify download of index files...")
    downloadindices(MATRIX_DB, ignoredb=ignoredb)
    println("used remote site is $(uf_remote.params.indexurl)")
    nothing
end

# will be called automatically once after `using`, `import`, `require`.
function __init__()
    try init() catch ex; @warn "exception during initialization: '$ex'"
    end
end

end # end module


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Utilities for computing and printing IterationStats.

function max_primal_violation(
  problem::QuadraticProgrammingProblem,
  primal_vec::AbstractVector{Float64},
)
  return norm(compute_primal_residual(problem, primal_vec), Inf)
end

"""
Given quadratic program and a primal solution vector as input this function
computes a vector containing the individual violations of each of the primal
constraints (both equalities and inequalities), variable lower bounds and upper
bounds respectively.
"""
function compute_primal_residual(
  problem::QuadraticProgrammingProblem,
  primal_vec::AbstractVector{Float64},
)
  activities = problem.constraint_matrix * primal_vec

  if isempty(equality_range(problem))
    equality_violation = []
  else
    equality_violation =
      problem.right_hand_side[equality_range(problem)] -
      activities[equality_range(problem)]
  end

  if isempty(inequality_range(problem))
    inequality_violation = []
  else
    inequality_violation =
      max.(
        problem.right_hand_side[inequality_range(problem)] -
        activities[inequality_range(problem)],
        0.0,
      )
  end

  lower_bound_violation = max.(problem.variable_lower_bound - primal_vec, 0.0)
  upper_bound_violation = max.(primal_vec - problem.variable_upper_bound, 0.0)
  return [
    equality_violation
    inequality_violation
    lower_bound_violation
    upper_bound_violation
  ]
end

""" Computes the primal objective.
"""
function primal_obj(
  problem::QuadraticProgrammingProblem,
  primal_solution::AbstractVector{Float64},
)
  return problem.objective_constant +
         problem.objective_vector' * primal_solution +
         0.5 * (primal_solution' * problem.objective_matrix * primal_solution)
end

struct DualStats
  dual_objective::Float64
  dual_residual::Vector{Float64}
  reduced_costs::Vector{Float64}
end

"""
  reduced_costs_dual_objective_contribution(
    variable_lower_bound,
    variable_upper_bound,
    reduced_costs
  )
This function returns the contribution of the
reduced costs to the dual objective value.
"""
function reduced_costs_dual_objective_contribution(
  variable_lower_bound::Vector{Float64},
  variable_upper_bound::Vector{Float64},
  reduced_costs::Vector{Float64},
)
  dual_objective_contribution = 0.0
  for i in 1:length(variable_lower_bound)
    if reduced_costs[i] == 0.0
      continue
    elseif reduced_costs[i] > 0.0
      # A positive reduced cost is associated with a binding lower bound.
      bound_value = variable_lower_bound[i]
    else
      # A negative reduced cost is associated with a binding upper bound.
      bound_value = variable_upper_bound[i]
    end
    if !isfinite(bound_value)
      return -Inf
    else
      dual_objective_contribution += bound_value * reduced_costs[i]
    end
  end

  return dual_objective_contribution
end

"""
  compute_reduced_costs_from_primal_gradient(
    variable_lower_bound,
    variable_upper_bound,
    primal_gradient
  )
where primal_gradient is the gradient of the Lagrangian with respect to the
primal vector. This function returns a vector of reduced costs.
"""
function compute_reduced_costs_from_primal_gradient(
  variable_lower_bound::Vector{Float64},
  variable_upper_bound::Vector{Float64},
  primal_gradient::AbstractVector{Float64},
)
  primal_size = length(primal_gradient)
  reduced_costs = zeros(primal_size)
  for i in 1:primal_size
    if primal_gradient[i] > 0.0
      bound_value = variable_lower_bound[i]
    else
      bound_value = variable_upper_bound[i]
    end
    if isfinite(bound_value)
      reduced_costs[i] = primal_gradient[i]
    end
  end

  return reduced_costs
end

"""
Given a dual vector (i.e., multipliers on the linear constraints) and a primal
vector, returns a DualStats object.
The objective function is linearized about the given
primal solution. Positive values of the dual variables correspond to the lower
bound on the corresponding constraint activities. Reduced costs are computed
from scratch as implied by the primal and dual vectors.
"""
function compute_dual_stats(
  problem::QuadraticProgrammingProblem,
  primal_solution::AbstractVector{Float64},
  dual_solution::AbstractVector{Float64},
)
  objective_product = problem.objective_matrix * primal_solution
  primal_gradient =
    compute_primal_gradient(problem, primal_solution, dual_solution)
  reduced_costs = compute_reduced_costs_from_primal_gradient(
    problem.variable_lower_bound,
    problem.variable_upper_bound,
    primal_gradient,
  )

  # Duals on the inequalities must be nonnegative.
  if !isempty(inequality_range(problem))
    dual_residual = max.(-dual_solution[inequality_range(problem)], 0.0)
  else
    dual_residual = []
  end
  reduced_cost_violations = primal_gradient .- reduced_costs
  dual_residual = [dual_residual; reduced_cost_violations]

  # The dual objective excluding the reduced costs.
  # The quadratic term 0.5 x^T Q x in the objective is linearized as
  #   0.5 x^T Q x >= 0.5 x_0^T Q x_0 + x_0^T Q (x - x_0).
  # The -0.5 in dual_obj's initializer comes from the two x_0^T Q x_0 terms in
  # the above.
  base_dual_objective =
    problem.right_hand_side' * dual_solution + problem.objective_constant -
    0.5 * objective_product' * primal_solution

  dual_objective =
    base_dual_objective + reduced_costs_dual_objective_contribution(
      problem.variable_lower_bound,
      problem.variable_upper_bound,
      reduced_costs,
    )

  return DualStats(dual_objective, dual_residual, reduced_costs)
end

""" Computes the dual objective. It uses the bounds on the primal variables to
correct for any infeasibility in the dual, producing a valid lower bound.
Returns -Inf if not all variables have finite bounds.
"""
function corrected_dual_obj(
  problem::QuadraticProgrammingProblem,
  dual_stats::DualStats,
)
  if norm(dual_stats.dual_residual, Inf) == 0.0
    return dual_stats.dual_objective
  else
    return -Inf
  end
end

function corrected_dual_obj(
  problem::QuadraticProgrammingProblem,
  primal_solution::AbstractVector{Float64},
  dual_solution::AbstractVector{Float64},
)
  dual_stats = compute_dual_stats(problem, primal_solution, dual_solution)
  return corrected_dual_obj(problem, dual_stats)
end

"""
Returns a `ConvergenceInformation` object given a `QuadraticProgrammingProblem`
with given primal and dual vectors.
"""
function compute_convergence_information(
  problem::QuadraticProgrammingProblem,
  qp_cache::CachedQuadraticProgramInfo,
  primal_iterate::Vector{Float64},
  dual_iterate::Vector{Float64},
  eps_ratio::Float64,
  candidate_type::PointType,
)
  num_constraints, num_vars = size(problem.constraint_matrix)
  @assert length(primal_iterate) == num_vars
  @assert length(dual_iterate) == num_constraints

  convergence_info = ConvergenceInformation()

  primal_residual = compute_primal_residual(problem, primal_iterate)
  convergence_info.primal_objective = primal_obj(problem, primal_iterate)
  convergence_info.l_inf_primal_residual = norm(primal_residual, Inf)
  convergence_info.l2_primal_residual = norm(primal_residual, 2)
  convergence_info.relative_l_inf_primal_residual =
    convergence_info.l_inf_primal_residual /
    (eps_ratio + qp_cache.l_inf_norm_primal_right_hand_side)
  convergence_info.relative_l2_primal_residual =
    convergence_info.l2_primal_residual /
    (eps_ratio + qp_cache.l2_norm_primal_right_hand_side)
  convergence_info.l_inf_primal_variable = norm(primal_iterate, Inf)
  convergence_info.l2_primal_variable = norm(primal_iterate, 2)

  dual_stats = compute_dual_stats(problem, primal_iterate, dual_iterate)
  convergence_info.dual_objective = dual_stats.dual_objective
  convergence_info.l_inf_dual_residual = norm(dual_stats.dual_residual, Inf)
  convergence_info.l2_dual_residual = norm(dual_stats.dual_residual, 2)
  convergence_info.relative_l_inf_dual_residual =
    convergence_info.l_inf_dual_residual /
    (eps_ratio + qp_cache.l_inf_norm_primal_linear_objective)
  convergence_info.relative_l2_dual_residual =
    convergence_info.l2_dual_residual /
    (eps_ratio + qp_cache.l2_norm_primal_linear_objective)
  convergence_info.l_inf_dual_variable = norm(dual_iterate, Inf)
  convergence_info.l2_dual_variable = norm(dual_iterate, 2)

  convergence_info.corrected_dual_objective =
    corrected_dual_obj(problem, dual_stats)

  gap = abs(convergence_info.primal_objective - convergence_info.dual_objective)
  abs_obj =
    abs(convergence_info.primal_objective) +
    abs(convergence_info.dual_objective)
  convergence_info.relative_optimality_gap = gap / (eps_ratio + abs_obj)

  convergence_info.candidate_type = candidate_type

  return convergence_info
end

"""
Returns an `InfeasibilityInformation` object given a
`QuadraticProgrammingProblem` and estimates for the primal and dual rays.
The rays do not need to be pre-scaled to have Inf-norm equal to 1.0.
"""
function compute_infeasibility_information(
  problem::QuadraticProgrammingProblem,
  primal_ray_estimate::Vector{Float64},
  dual_ray_estimate::Vector{Float64},
  candidate_type::PointType,
)

  infeas_info = InfeasibilityInformation()

  primal_ray_inf_norm = norm(primal_ray_estimate, Inf)
  if !iszero(primal_ray_inf_norm)
    primal_ray_estimate /= primal_ray_inf_norm
  end

  homogeneous_primal = linear_programming_problem(
    [isfinite(l) ? 0.0 : -Inf for l in problem.variable_lower_bound],
    [isfinite(u) ? 0.0 : Inf for u in problem.variable_upper_bound],
    problem.objective_vector,
    0.0,  # objective_constant
    problem.constraint_matrix,
    zeros(length(problem.right_hand_side)),
    problem.num_equalities,
  )

  homogeneous_residual =
    compute_primal_residual(homogeneous_primal, primal_ray_estimate)
  infeas_info.max_primal_ray_infeasibility = norm(homogeneous_residual, Inf)
  infeas_info.primal_ray_linear_objective =
    problem.objective_vector' * primal_ray_estimate
  infeas_info.primal_ray_quadratic_norm =
    norm(problem.objective_matrix * primal_ray_estimate, Inf)

  homogeneous_dual = linear_programming_problem(
    problem.variable_lower_bound,
    problem.variable_upper_bound,
    zeros(length(problem.objective_vector)),
    0.0,  # objective_constant
    problem.constraint_matrix,
    problem.right_hand_side,
    problem.num_equalities,
  )

  homogeneous_dual_stats =
    compute_dual_stats(homogeneous_dual, primal_ray_estimate, dual_ray_estimate)

  scaling_factor = max(
    norm(dual_ray_estimate, Inf),
    norm(homogeneous_dual_stats.reduced_costs, Inf),
  )
  if !iszero(scaling_factor)
    infeas_info.max_dual_ray_infeasibility =
      norm(homogeneous_dual_stats.dual_residual, Inf) / scaling_factor
    infeas_info.dual_ray_objective =
      homogeneous_dual_stats.dual_objective / scaling_factor
  else
    infeas_info.max_dual_ray_infeasibility = 0.0
    infeas_info.dual_ray_objective = 0.0
  end

  infeas_info.candidate_type = candidate_type

  return infeas_info
end

"""
Returns a `IterationStats` object given a `QuadraticProgrammingProblem`,
primal and dual vectors, estimates for primal and dual rays, and other
iteration statistics corresponding to the fields in the struct.
"""
function compute_iteration_stats(
  problem::QuadraticProgrammingProblem,
  qp_cache::CachedQuadraticProgramInfo,
  primal_iterate::Vector{Float64},
  dual_iterate::Vector{Float64},
  primal_ray_estimate::Vector{Float64},
  dual_ray_estimate::Vector{Float64},
  iteration_number::Integer,
  cumulative_kkt_matrix_passes::Float64,
  cumulative_time_sec::Float64,
  eps_optimal_absolute::Float64,
  eps_optimal_relative::Float64,
  step_size::Float64,
  primal_weight::Float64,
  candidate_type::PointType,
)
  num_constraints, num_vars = size(problem.constraint_matrix)
  @assert length(primal_iterate) == num_vars
  @assert length(primal_ray_estimate) == num_vars
  @assert length(dual_iterate) == num_constraints
  @assert length(dual_ray_estimate) == num_constraints

  stats = IterationStats()
  stats.iteration_number = iteration_number
  stats.cumulative_kkt_matrix_passes = cumulative_kkt_matrix_passes
  stats.cumulative_time_sec = cumulative_time_sec

  stats.convergence_information = [
    compute_convergence_information(
      problem,
      qp_cache,
      primal_iterate,
      dual_iterate,
      eps_optimal_absolute / eps_optimal_relative,
      candidate_type,
    ),
  ]
  stats.infeasibility_information = [
    compute_infeasibility_information(
      problem,
      primal_ray_estimate,
      dual_ray_estimate,
      candidate_type,
    ),
  ]
  stats.step_size = step_size
  stats.primal_weight = primal_weight
  stats.method_specific_stats = Dict{AbstractString,Float64}()

  return stats
end

"""
This code computes the unscaled iteration stats.
The input iterates to this function have been scaled according to
scaled_problem.
"""
function evaluate_unscaled_iteration_stats(
  scaled_problem::ScaledQpProblem,
  qp_cache::CachedQuadraticProgramInfo,
  termination_criteria::TerminationCriteria,
  record_iteration_stats::Bool,
  primal_solution::AbstractVector{Float64},
  dual_solution::AbstractVector{Float64},
  iteration::Int64,
  cumulative_time::Float64,
  cumulative_kkt_passes::Float64,
  eps_optimal_absolute::Float64,
  eps_optimal_relative::Float64,
  step_size::Float64,
  primal_weight::Float64,
  candidate_type::PointType,
)
  # Unscale iterates.
  original_primal_solution::Vector{Float64} =
    primal_solution ./ scaled_problem.variable_rescaling
  original_dual_solution::Vector{Float64} =
    dual_solution ./ scaled_problem.constraint_rescaling

  return compute_iteration_stats(
    scaled_problem.original_qp,
    qp_cache,
    original_primal_solution,
    original_dual_solution,
    original_primal_solution,  # ray estimate
    original_dual_solution,  # ray estimate
    iteration - 1,
    cumulative_kkt_passes,
    cumulative_time,
    eps_optimal_absolute,
    eps_optimal_relative,
    step_size,
    primal_weight,
    candidate_type,
  )
end

"""
   print_to_screen_this_iteration(termination_reason, iteration, verbosity,
      termination_evaluation_frequency)

Decides if we should print iteration stats to screen this iteration.
"""
function print_to_screen_this_iteration(
  termination_reason::Union{TerminationReason,Bool},
  iteration::Int64,
  verbosity::Int64,
  termination_evaluation_frequency::Int32,
)
  if verbosity >= 2
    if termination_reason == false
      num_of_evaluations = (iteration - 1) / termination_evaluation_frequency
      if verbosity >= 9
        display_frequency = 1
      elseif verbosity >= 6
        display_frequency = 3
      elseif verbosity >= 5
        display_frequency = 10
      elseif verbosity >= 4
        display_frequency = 20
      elseif verbosity >= 3
        display_frequency = 50
      else
        return iteration == 1
      end
      # print_to_screen_this_iteration is true every
      # display_frequency * termination_evaluation_frequency iterations.
      return mod(num_of_evaluations, display_frequency) == 0
    else
      return true
    end
  else
    return false
  end
end

"""
  display_iteration_stats_heading(show_infeasibility)

The heading for the iteration stats table. If show_infeasibility is true then
an extended table is printed that includes infeasibility information. See
README.md for documentation on what each heading means.
"""
function display_iteration_stats_heading(show_infeasibility::Bool)
  Printf.@printf(
    "%s | %s | %s | %s |",
    rpad("runtime", 24),
    rpad("residuals", 26),
    rpad(" solution information", 26),
    rpad("relative residuals", 23)
  )
  if show_infeasibility
    Printf.@printf(" %s | %s |", rpad("primal ray", 27), rpad("dual ray", 18))
  end
  println("")
  Printf.@printf(
    "%s %s %s | %s %s  %s | %s %s %s | %s %s %s |",
    rpad("#iter", 7),
    rpad("#kkt", 8),
    rpad("seconds", 7),
    rpad("pr norm", 8),
    rpad("du norm", 8),
    rpad("gap", 7),
    rpad(" pr obj", 9),
    rpad("pr norm", 8),
    rpad("du norm", 7),
    rpad("rel pr", 7),
    rpad("rel du", 7),
    rpad("rel gap", 7)
  )
  if show_infeasibility
    Printf.@printf(
      " %s %s %s | %s %s |",
      rpad("pr norm", 9),
      rpad("linear", 8),
      rpad("qu norm", 8),
      rpad("du norm", 9),
      rpad("dual obj", 8)
    )
  end
  print("\n")
end

function display_iteration_stats_heading(verbosity::Int64)
  if verbosity >= 7
    display_iteration_stats_heading(true)
  elseif verbosity >= 2
    display_iteration_stats_heading(false)
  end
end

"""
Make sure that a float is of a constant length, irrespective if it is negative
or positive.
"""
function lpad_float(number::Float64)
  return lpad(Printf.@sprintf("%.1e", number), 8)
end

"""
Displays a row of the iteration stats table.
"""
function display_iteration_stats(
  stats::IterationStats,
  show_infeasibility::Bool,
)
  # TODO: Decide if we want to show information about all the
  # convergence_information entries.
  if length(stats.convergence_information) > 0
    Printf.@printf(
      "%s  %.1e  %.1e | %.1e  %.1e  %s | %s  %.1e  %.1e | %.1e %.1e %.1e |",
      rpad(string(stats.iteration_number), 6),
      stats.cumulative_kkt_matrix_passes,
      stats.cumulative_time_sec,
      stats.convergence_information[1].l2_primal_residual,
      stats.convergence_information[1].l2_dual_residual,
      lpad_float(
        stats.convergence_information[1].primal_objective -
        stats.convergence_information[1].dual_objective,
      ),
      lpad_float(stats.convergence_information[1].primal_objective),
      stats.convergence_information[1].l2_primal_variable,
      stats.convergence_information[1].l2_dual_variable,
      stats.convergence_information[1].relative_l2_primal_residual,
      stats.convergence_information[1].relative_l2_dual_residual,
      stats.convergence_information[1].relative_optimality_gap
    )
  else
    Printf.@printf(
      "%s  %.1e  %.1e",
      rpad(string(stats.iteration_number), 6),
      stats.cumulative_kkt_matrix_passes,
      stats.cumulative_time_sec
    )
  end

  if show_infeasibility
    # TODO: Decide if we want to show information about all the
    # entries.
    if length(stats.infeasibility_information) > 0
      Printf.@printf(
        " %.1e  %s  %.1e  | %.1e  %s  |",
        stats.infeasibility_information[1].max_primal_ray_infeasibility,
        lpad_float(
          stats.infeasibility_information[1].primal_ray_linear_objective,
        ),
        stats.infeasibility_information[1].primal_ray_quadratic_norm,
        stats.infeasibility_information[1].max_dual_ray_infeasibility,
        lpad_float(stats.infeasibility_information[1].dual_ray_objective)
      )
    end
  end

  print("\n")
end

function display_iteration_stats(stats::IterationStats, verbosity::Int64)
  if verbosity >= 7
    display_iteration_stats(stats, true)
  else
    display_iteration_stats(stats, false)
  end
end

function print_infinity_norms(convergence_info::ConvergenceInformation)
  print("l_inf: ")
  Printf.@printf(
    "primal_res = %.3e, dual_res = %.3e, primal_var = %.3e, dual_var = %.3e",
    convergence_info.l_inf_primal_residual,
    convergence_info.l_inf_dual_residual,
    convergence_info.l_inf_primal_variable,
    convergence_info.l_inf_dual_variable
  )
  println()
end


text:
## ELASTIC u ##
function assemble_cell!(Ke, fe, model::Model{dim,1,Nothing,Nothing,E,Nothing}, cell, cv, nu) where {dim,E}

    reinit!(cv, cell)
    # get stiffness tensor of the cell
    cell_id = cell.current_cellid.x
    Dᵉ = model.material_properties[cell_id].elasticity.Dᵉ
    bodyforces = model.body_forces[cell_id].components

    # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
    @inbounds for q_point in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, q_point)
        for i in 1:nu
            δϵ2D = shape_symmetric_gradient(cv, q_point, i)
            δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
            δu = shape_value(cv, q_point, i)
            fe[i] += (δu ⋅ bodyforces) * dΩ # sign is ok
            for j in 1:i
                Δϵ2D = shape_symmetric_gradient(cv, q_point, j)
                Δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,Δϵ2D))
                Ke[i, j] += δϵ ⊡ Dᵉ ⊡ Δϵ * dΩ
            end
        end
    end
    symmetrize_lower!(Ke)
    apply_Neumann_bc!(fe, model, cell, nu)
end

## ELASTIC u,p ##
function assemble_cell!(Ke, fe, model::Model{dim,2,Nothing,Nothing,E,Nothing}, cell, ɛdev, cvu, cvp, nu, np, u▄, p▄) where {dim,E}

    reinit!(cvu, cell)
    reinit!(cvp, cell)

    cell_id = cell.current_cellid.x
    r = model.material_properties[cell_id]
    bodyforces = model.body_forces[cell_id].components
    Gmod, Kmod = r.elasticity.G, r.elasticity.K # get elast moduli

    # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
    @inbounds for q_point in 1:getnquadpoints(cvu)
        dΩ = getdetJdV(cvu, q_point)

        for i in 1:nu
            ɛdev[i] = dev(shape_symmetric_gradient(cvu, q_point, i))
        end

        for i in 1:nu
            #divδu = shape_divergence(cvu, q_point, i)
            δu = shape_value(cvu, q_point, i)
            δɛ = shape_symmetric_gradient(cvu, q_point, i)
            fe[i] += (δu ⋅ bodyforces) * dΩ
            for j in 1:i
                # Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * Gmod * ɛdev[i] ⊡ ɛdev[j] * dΩ # initial one
                Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * Gmod * δϵ ⊡ ɛdev[j] * dΩ #TODO check thid or revert change
            end
        end

        for i in 1:np
            δp = shape_value(cvp, q_point, i)
            for j in 1:nu
                divδu = shape_divergence(cvu, q_point, j)
                Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
            end
            for j in 1:i
                p = shape_value(cvp, q_point, j)
                Ke[BlockIndex((p▄, p▄), (i, j))] += - 1/Kmod * δp * p  * dΩ
            end

        end
    end

    symmetrize_lower!(Ke)

    # We integrate the Neumann boundary using the facevalues.
    # We loop over all the faces in the cell, then check if the face
    # is in our `"traction"` faceset.
    apply_Neumann_bc!(fe, model, cell, nu)

end

## Functions for viscosity based drucker-prager plasticity :

function assemble_cell!(Ke, fe, model::Model{dim,2,Nothing,TV,TE,TP}, cell, ɛdev, cvu, cvp, nu, np, u▄, p▄, nodal_vars_prev) where {dim,TV,TE,TP<:ViscousDruckerPrager}

    reinit!(cvu, cell)
    reinit!(cvp, cell)

    cell_id = cell.current_cellid.x
    r = model.material_properties[cell_id]
    bodyforces = model.body_forces[cell_id].components
    Gmod, Kmod = r.elasticity.G, r.elasticity.K # get elast moduli
    η = r.viscosity.η
    Δt = model.clock.Δt
    u_prev = nodal_vars_prev
    states = model.material_state[cell_id]

    if r.elasticity.ν != 0.5
        # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
        @inbounds for q_point in 1:getnquadpoints(cvu)

            state = states[q_point]
            ϵ_prev = function_symmetric_gradient(cvu, q_point,u_prev[1:nu])

            #ϵ_prev = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ_prev_2D))
            dΩ = getdetJdV(cvu, q_point)

            for i in 1:nu
                ɛdev[i] = dev(shape_symmetric_gradient(cvu, q_point, i))
            end

            for i in 1:nu
                #divδu = shape_divergence(cvu, q_point, i)
                δu = shape_value(cvu, q_point, i)
                δN = shape_symmetric_gradient(cvu, q_point, i)
                fe[i] += (ɛdev[i] ⊡ ((2η/Δt)*ϵ_prev) + δu ⋅ bodyforces) * dΩ
                # fe[i] += (δN ⊡ ((2η/Δt)*ϵ_prev) + δu ⋅ bodyforces) * dΩ
                for j in 1:i
                    ϵ = shape_symmetric_gradient(cvu, q_point, j)
                    Ke[BlockIndex((u▄, u▄), (i, j))] +=  ɛdev[i] ⊡ (2Gmod * ɛdev[j] + (2η/Δt) * ϵ) * dΩ # initial one
                    # Ke[BlockIndex((u▄, u▄), (i, j))] +=  δN ⊡ (2Gmod * ɛdev[j] + (2η/Δt) * ϵ) * dΩ
                end
            end

            for i in 1:np
                δp = shape_value(cvp, q_point, i)
                for j in 1:nu
                    divδu = shape_divergence(cvu, q_point, j)
                    Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
                end
                for j in 1:i
                    p = shape_value(cvp, q_point, j)
                    Ke[BlockIndex((p▄, p▄), (i, j))] += - 1/Kmod * δp * p  * dΩ
                end

            end
        end
    else
        # code Sister pde
        Z = Gmod*Δt/(η + Gmod*Δt)
        @inbounds for q_point in 1:getnquadpoints(cvu)

            state = states[q_point]
            ϵ_prev_2D = function_symmetric_gradient(cvu, q_point,u_prev[1:nu])
            ϵ_prev = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ_prev_2D))
            σ_prev = state.σ

            if cell_id == 1
                println("ϵ_prev : ", ϵ_prev)
                println("σ_prev : ", σ_prev)
            end

            dΩ = getdetJdV(cvu, q_point)

            # ∇N = ɛdev
            # for i in 1:nu
            #     ∇N[i] = shape_symmetric_gradient(cvu, q_point, i)
            # end

            for i in 1:nu
                #divδu = shape_divergence(cvu, q_point, i)
                δu = shape_value(cvu, q_point, i)
                δN_2D = shape_gradient(cvu, q_point, i)
                δN = Tensor{2,3}((i,j)->get_3D_func(i,j,δN_2D))
                fe[i] += (δN ⊡ ( (2Z*η/Δt)*ϵ_prev - (1 - Z)*σ_prev ) + δu ⋅ bodyforces) * dΩ
                # fe[i] += (δN ⊡ ((2η/Δt)*ϵ_prev) + δu ⋅ bodyforces) * dΩ
                for j in 1:i # 1:i initialy with symmetrize_lower
                    ϵ_2D = shape_symmetric_gradient(cvu, q_point, j)
                    ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ_2D))
                    Ke[BlockIndex((u▄, u▄), (i, j))] +=  δN ⊡ ((2Z*η/Δt) * ϵ) * dΩ # initial one
                    # Ke[BlockIndex((u▄, u▄), (i, j))] +=  δN ⊡ (2Gmod * ɛdev[j] + (2η/Δt) * ϵ) * dΩ
                end
            end

            for i in 1:np
                δp = shape_value(cvp, q_point, i)
                for j in 1:nu
                    divδu = shape_divergence(cvu, q_point, j)
                    Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
                end
                # for j in 1:i
                #     p = shape_value(cvp, q_point, j)
                #     Ke[BlockIndex((p▄, p▄), (i, j))] += 0.0
                # end

            end
        end
    end

    symmetrize_lower!(Ke)

    # We integrate the Neumann boundary using the facevalues.
    # We loop over all the faces in the cell, then check if the face
    # is in our `"traction"` faceset.
    apply_Neumann_bc!(fe, model, cell, nu)

end

# TODO : put config out of cell loop, use DiffResults to get func eval + jacobian in a single pass.
function assemble_cell_AD!(re, Ke, model::Model{dim,2,TD,TV,TE,TP}, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev) where {dim,TD,TV<:Viscosity,TE,TP<:ViscousDruckerPrager}
    assemble_res_cell!(re, model, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev)
    f!(re,nodal_vars_el) = assemble_res_cell!(re, model, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev)
    re2 = similar(re)
    jacobian_config = ForwardDiff.JacobianConfig(f!, re2, nodal_vars_el)#, chunk = Chunk(nodal_vars))
    ForwardDiff.jacobian!(Ke, f!, re2, nodal_vars_el, jacobian_config, Val{true}())
end

## DAMAGED-ELASTO-VISCO-PLASTIC ##
# import Base.getindex
# Base.getindex(barray::PseudoBlockArray{Float64,1,Array{Float64,1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, blockindex::BlockIndex{1}) =
# getindex(barray::PseudoBlockArray{T,N,R,BS} where BS<:Tuple{Vararg{AbstractUnitRange{Int64},N}} where R<:AbstractArray{T,N}, blockindex::BlockIndex{N}) where {T, N}
function assemble_cell_elast!(re, Ke, model::Model{DIM,2,TD,V,E,P}, cell, cvu, cvD, nu, nD, u▄, D▄, ue, De, De_prev) where {DIM,TD,V,E,P}

    reinit!(cvu, cell)
    reinit!(cvD, cell)

    cell_id = cell.current_cellid.x
    r = model.material_properties[cell_id]
    bodyforces = model.body_forces[cell_id].components

    # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
    @inbounds for q_point in 1:getnquadpoints(cvu)

        # get qp state
        state = model.material_state[cell_id][q_point]

        # get Gauss differential term
        dΩ = getdetJdV(cvu, q_point)

        # get strain
        ϵ2D = function_symmetric_gradient(cvu, q_point, ue)
        ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ2D))

        # with homogeneous elasticity
        # r_nodamage = Rheology(nothing,r.viscosity,r.elasticity,r.plasticity) # TODO find a better way to dispatch on compute_stress_tangent
        # σ, C = compute_stress_tangent(ϵ, r_nodamage, state, model.clock, noplast = true)

        # with damaged elasticity :
        D_prev = function_value(cvD, q_point, De_prev)
        C = compute_damaged_stiffness_tensor(r,ϵ,D_prev)
        σ = C ⊡ ϵ

        if (cell_id == 1) & (q_point == 1)
            println("σ = ", σ)
            println("ϵ = ", ϵ)
        end
        # update temp state :
        state.temp_σ = σ
        state.temp_ϵ = ϵ
        for i in 1:nu
            #
            δϵ2D = shape_symmetric_gradient(cvu, q_point, i)
            δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
            δu = shape_value(cvu, q_point, i)
            # increment residual
            #re[BlockIndex(u▄, i)] += ((δϵ ⊡ σ) - (δu ⋅ bodyforces)) * dΩ
            re[i] += ((δϵ ⊡ σ) - (δu ⋅ bodyforces)) * dΩ
            δϵC = δϵ ⊡ C
            for j in 1:i
                Δϵ2D = shape_symmetric_gradient(cvu, q_point, j)
                Δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,Δϵ2D))
                Ke[BlockIndex((u▄, u▄), (i, j))] += δϵC ⊡ Δϵ * dΩ
            end
        end

        for i in 1:nD
            #re[BlockIndex(D▄, i)] += ((D - D_prev)/model.clock.Δt - dDdt) * dΩ
            re[i+nu] = 1e9 #* dΩ
            Ke[BlockIndex((D▄, D▄), (i, i))] = 1.0 #* dΩ
        end
    end
    symmetrize_lower!(Ke) #when j in 1:i
    # We integrate the Neumann boundary using the facevalues.
    # We loop over all the faces in the cell, then check if the face
    # is in our `"traction"` faceset.
    apply_Neumann_bc!(re, model, cell, nu ; inc_sign = -)

end

function assemble_res_cell!(re, model::Model{dim,2,TD,Nothing,TE,TP}, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev) where {dim,TD,TE,TP}

    reinit!(cvu, cell)
    reinit!(cvD, cell)

    cell_id = cell.current_cellid.x
    r = model.material_properties[cell_id]
    bodyforces = model.body_forces[cell_id].components

    ue = nodal_vars_el[1:nu]
    De = exp.(nodal_vars_el[nu+1:end])
    De_prev = exp.(nodal_vars_el_prev[nu+1:end])
    dDdt = 0.0

    fill!(re, 0)

    if (cell_id == 1) & !(eltype(De) <: ForwardDiff.Dual)
        println("De at cell id 1 : ")
        display(De)
    end
    # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
    @inbounds for q_point in 1:getnquadpoints(cvu)

        # get qp state
        # sol_type = eltype(nodal_vars)
        # if sol_type <: ForwardDiff.Dual
        #     state_dual = PlasticMaterialState{sol_type}()
        # else
        state = model.material_state[cell_id][q_point]
        # end

        # get Gauss differential term
        dΩ = getdetJdV(cvu, q_point)

        # get strain
        ϵ2D = function_symmetric_gradient(cvu, q_point, ue)
        ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ2D))

        # get damage variable
        D = function_value(cvD, q_point, De)
        D_prev = function_value(cvD, q_point, De_prev)

        ###### TEST
        if D <= r.damage.D₀
            println("cell_id = ",cell_id)
            println("ϵ = ",ϵ)
            println("D = ",D)
            println("all_D_qp = ", De)
        end
        ######
        if D <= 1 # TEST TODO put back D < 1
            # get stress
            σ = compute_σij(r,D,ϵ)

            # get damage growth rate
            KI = compute_KI(r,σ,D)

            ###### TEST
            if isnan(KI)
                println("KI is NaN")
                println("\tcell_id = ",cell_id) # 1
                println("\tϵ = ",ϵ)
                println("\tσ = ",σ) # All NaN !!!
                println("\tD = ",D) # OK
            end
            ######

            dDdt = compute_subcrit_damage_rate(r, KI, D)

            #### TEST
            if (cell_id == 1) & (q_point == 1) & !(eltype(σ) <: ForwardDiff.Dual)
                println("σ at cell_id & qp = 1 : ", σ)
                println("σ_prev at cell_id & qp = 1 : ", state.temp_σ)
                println("KI at cell_id & qp = 1 : ", KI)
                println("dDdt at cell_id & qp = 1 : ", dDdt)
            end
            if D <= r.damage.D₀
                println("σ = ",σ)
            end
            ####

            for i in 1:nu
                #
                δϵ2D = shape_symmetric_gradient(cvu, q_point, i)
                δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
                δu = shape_value(cvu, q_point, i)
                # increment residual
                #re[BlockIndex(u▄, i)] += ((δϵ ⊡ σ) - (δu ⋅ bodyforces)) * dΩ
                re[i] += ((δϵ ⊡ σ) - (δu ⋅ bodyforces)) * dΩ
            end

            for i in 1:nD
                #re[BlockIndex(D▄, i)] += ((D - D_prev)/model.clock.Δt - dDdt) * dΩ
                #(D <= D_prev) && (regularization_term = abs(D-D_prev))
                δD = shape_value(cvD, q_point, i)

                re[i+nu] += δD * ((D - D_prev)/model.clock.Δt - dDdt ) * dΩ #
            end

            # update qp state
            sol_type = eltype(nodal_vars_el)
            if !(sol_type <: ForwardDiff.Dual) # no state update when AD
                set_temp_state!(r,model.clock,state,σ,ϵ)
            end
        end

        # if (D >= 1) | (D_prev + dDdt*model.clock.Δt >= 1) # loss of cohesion, drucker-Prager rheology
        #     rp = Rheology(nothing,nothing,r.elasticity,r.plasticity)
        #
        #     sol_type = eltype(nodal_vars_el)
        #     if sol_type <: ForwardDiff.Dual
        #         state_dual = PlasticMaterialState{sol_type}() # initializes with zero values
        #         add_state!(state_dual,state) # add current state
        #         σ, _ = compute_stress_tangent(ϵ, rp, state_dual, model.clock, noplast = false)
        #     else
        #         σ, _ = compute_stress_tangent(ϵ, rp, state, model.clock, noplast = false)
        #     end
        #
        #     dDdt = 0.0
        #     for i in 1:nu
        #         δϵ2D = shape_symmetric_gradient(cvu, q_point, i)
        #         δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
        #         δu = shape_value(cvu, q_point, i)
        #         # increment residual
        #         #re[BlockIndex(u▄, i)] += ((δϵ ⊡ σ) - (δu ⋅ bodyforces)) * dΩ
        #         re[i] += ((δϵ ⊡ σ) - (δu ⋅ bodyforces)) * dΩ
        #     end
        #
        #     for i in 1:nD
        #         #re[BlockIndex(D▄, i)] += ((D - D_prev)/model.clock.Δt - dDdt) * dΩ
        #         #(D <= D_prev) && (regularization_term = abs(D-D_prev))
        #         δD = shape_value(cvD, q_point, i)
        #         re[i+nu] += δD * ((D - D_prev)/model.clock.Δt - dDdt ) * dΩ #
        #     end
        # end
    end

    # We integrate the Neumann boundary using the facevalues.
    # We loop over all the faces in the cell, then check if the face
    # is in our `"traction"` faceset.
    apply_Neumann_bc!(re, model, cell, nu ; inc_sign = -)

end

# TODO : put config out of cell loop, use DiffResults to get func eval + jacobian in a single pass.
function assemble_cell_AD!(re, Ke, model::Model{dim,2,TD,Nothing,TE,TP}, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev) where {dim,TD,TE,TP}
    assemble_res_cell!(re, model, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev)
    f!(re,nodal_vars_el) = assemble_res_cell!(re, model, cell, cvu, cvD, nu, nD, u▄, D▄, nodal_vars_el, nodal_vars_el_prev)
    re2 = similar(re)
    jacobian_config = ForwardDiff.JacobianConfig(f!, re2, nodal_vars_el)#, chunk = Chunk(nodal_vars))
    ForwardDiff.jacobian!(Ke, f!, re2, nodal_vars_el, jacobian_config, Val{true}())

    ### TEST ###
    #Ke = compute_jacobian(f!, re2 ; eps = 1e-8)
    # if (cell.current_cellid.x == 1)
    #     println("passé ici !!!")
    #     J = compute_jacobian(f!, nodal_vars_el ; eps = 1e-9)
    #     println("finitediff jacobian for cell_id = 1, iter 2:")
    #     display(J[1:2,:])
    #     println("difference with AD for cell_id = 1, iter 2 :")
    #     display(abs.(Array(Ke)[1:2,:].-J[1:2,:]))
    # end
    ############
    return Ke
end

###################################################
#### DAMAGE with displacement only formulation ####
###################################################

function assemble_cell!(Ke, re, model::Model{dim,1,TD,V,E,P}, cell, cv, n_basefuncs, ue, noplast = false) where {dim,TD,V,E,P}
    # unpack
    cell_id = cell.current_cellid.x
    r, states = model.material_properties[cell_id], model.material_state[cell_id]
    bodyforces = model.body_forces[cell_id].components
    #norm_bodyforces = norm(bodyforces)
    #unit_bodyforces = bodyforces./norm_bodyforces

    # reinit cellvalues
    reinit!(cv, cell)

    @inbounds for q_point in 1:getnquadpoints(cv)
        # For each integration point, compute stress and material stiffness
        ϵ2D = function_symmetric_gradient(cv, q_point, ue)
        # Total strain recomputed each time because of the newton correction
        ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ2D))
        # @timeit "tangent computation"
        σ, C = compute_stress_tangent(ϵ, r, states[q_point], model.clock, noplast = noplast)
        dΩ = getdetJdV(cv, q_point)
        for i in 1:n_basefuncs #TODO change loops order
            δϵ2D = shape_symmetric_gradient(cv, q_point, i)
            δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
            δu = shape_value(cv,q_point,i) # keep it 2D to fit body_forces dim
            re[i] += ((δϵ ⊡ σ) - (δu⋅bodyforces)) * dΩ # add internal force and external body forces to residual. ⊡ : double contraction, ⋅ : single contraction
            δϵC = δϵ ⊡ C
            for j in 1:i
                Δϵ2D = shape_symmetric_gradient(cv, q_point, j)
                Δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,Δϵ2D))
                Ke[i, j] += δϵC ⊡ Δϵ * dΩ
            end
        end
    end
    symmetrize_lower!(Ke) #when j in 1:i

    # Add traction as a negative contribution to the element residual `re`:
    apply_Neumann_bc!(re, model, cell, n_basefuncs ; inc_sign = -)
end

### individual assembly of res function and derivative K used for linesearch

function assemble_cell_res!(re, model::Model{dim,1,TD,V,E,P}, cell, cv, n_basefuncs, ue, noplast = false) where {dim,TD,V,E,P}
    # unpack
    cell_id = cell.current_cellid.x
    r, states = model.material_properties[cell_id], model.material_state[cell_id]
    bodyforces = model.body_forces[cell_id].components
    #norm_bodyforces = norm(bodyforces)
    #unit_bodyforces = bodyforces./norm_bodyforces

    # reinit cellvalues
    reinit!(cv, cell)

    @inbounds for q_point in 1:getnquadpoints(cv)
        # For each integration point, compute stress and material stiffness
        ϵ2D = function_symmetric_gradient(cv, q_point, ue)
        # Total strain recomputed each time because of the newton correction
        ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ2D))
        # @timeit "tangent computation"
        σ, C = compute_stress_tangent(ϵ, r, states[q_point], model.clock, noplast = noplast)
        dΩ = getdetJdV(cv, q_point)
        for i in 1:n_basefuncs #TODO change loops order
            δϵ2D = shape_symmetric_gradient(cv, q_point, i)
            δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
            δu = shape_value(cv,q_point,i) # keep it 2D to fit body_forces dim
            re[i] += ((δϵ ⊡ σ) - (δu⋅bodyforces)) * dΩ # add internal force and external body forces to residual. ⊡ : double contraction, ⋅ : single contraction
        end
    end

    # Add traction as a negative contribution to the element residual `re`:
    apply_Neumann_bc!(re, model, cell, n_basefuncs)
end

function assemble_cell_K!(Ke, model::Model{dim,1,TD,V,E,P}, cell, cv, n_basefuncs, ue, noplast = false) where {dim,TD,V,E,P}
    # unpack
    cell_id = cell.current_cellid.x
    r, states = model.material_properties[cell_id], model.material_state[cell_id]
    #norm_bodyforces = norm(bodyforces)
    #unit_bodyforces = bodyforces./norm_bodyforces

    # reinit cellvalues
    reinit!(cv, cell)

    @inbounds for q_point in 1:getnquadpoints(cv)
        # For each integration point, compute stress and material stiffness
        ϵ2D = function_symmetric_gradient(cv, q_point, ue)
        # Total strain recomputed each time because of the newton correction
        ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ2D))
        # @timeit "tangent computation"
        σ, C = compute_stress_tangent(ϵ, r, states[q_point], model.clock, noplast = noplast)
        dΩ = getdetJdV(cv, q_point)
        for i in 1:n_basefuncs #TODO change loops order
            δϵ2D = shape_symmetric_gradient(cv, q_point, i)
            δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,δϵ2D))
            for j in 1:i
                Δϵ2D = shape_symmetric_gradient(cv, q_point, j)
                Δϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,Δϵ2D))
                Ke[i, j] += δϵ ⊡ C ⊡ Δϵ * dΩ
            end
        end
    end
    symmetrize_lower!(Ke) #when j in 1:i
end


text:
module CGcoefficient

export CG, threeJ, sixJ, nineJ,
       dCG, d3j, d6j, d9j,
       SqrtRational,
       HalfInt,
       iphase,
       is_same_parity,
       check_jm,
       check_couple

include("SqrtRational.jl")
include("util.jl")
include("WignerSymbols.jl")

end # module CGcoefficient


text:
using Test

using JuliaRL
using Random

# include("utils.jl")

# function test_env()
#     rng = MersenneTwister(0)
#     @testset "MountainCar" begin
#         @test test_construction(MountainCar, 0.1, 0.04)
#         @test test_construction(MountainCar)
#         @test test_construction(MountainCar, rng)

#         mc = MountainCar(0.1, 0.04, true)
#         @test get_state(mc) == [0.7647058823529412, 0.7857142857142857]
#         rng = MersenneTwister(0)
#         mc = MountainCar(rng)
#         @test get_state(mc) == [0.2236475079774124, 0.0]

#         mc = MountainCar(true)
#         action = 1
#         step!(mc, action)
#         @test get_state(mc) == [0.7038235294117646, 0.475]
#         rng = MersenneTwister(0)
#         start!(mc; rng=rng)
#         @test get_state(mc) == [0.8374397105749485, 0.5]
#     end
# end

function runtests()
    # test_env()
end

runtests()

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct ConstDensity{T<:Real,B<:Union{VarVolumeBounds,Missing}} <: DistLikeDensity
    value::LogDVal{T}
    bounds::B
end

ConstDensity(value::LogDVal) = ConstDensity(value, missing)

ConstDensity(::typeof(normalize), bounds::Union{VarVolumeBounds,Missing} = missing) =
    ConstDensity(LogDVal(-log_volume(spatialvolume(bounds))), bounds)

Base.@deprecate ConstDensity(::typeof(one)) ConstDensity(LogDVal(0))
Base.@deprecate ConstDensity(::typeof(one), ::Missing) ConstDensity(LogDVal(0), missing)
Base.@deprecate ConstDensity(::typeof(one), bounds::VarVolumeBounds) ConstDensity(LogDVal(0), bounds)

Base.convert(::Type{ConstDensity}, value::LogDVal) = ConstDensity(value)
Base.convert(::Type{AbstractDensity}, value::LogDVal) = convert(ConstDensity, value)

@inline eval_logval_unchecked(density::ConstDensity, v::Any) = logvalof(density.value)

var_bounds(density::ConstDensity) = density.bounds

ValueShapes.varshape(density::ConstDensity{T,<:VarVolumeBounds}) where T = ArrayShape{Real}(totalndof(density.bounds))

ValueShapes.varshape(density::ConstDensity{T,Missing}) where T = missing

Distributions.sampler(density::ConstDensity) = spatialvolume(var_bounds(density))

function Statistics.cov(density::ConstDensity{T,<:HyperRectBounds}) where T
    vol = spatialvolume(var_bounds(density))
    #flat_var = (vol.hi - vol.lo).^2 / 12
    flat_var = var.(Uniform.(vol.lo, vol.hi))
    Matrix(PDiagMat(flat_var))
end