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

"""
    create_experiment_template(actions, client_token, description, role_arn, stop_conditions)
    create_experiment_template(actions, client_token, description, role_arn, stop_conditions, params::Dict{String,<:Any})

Creates an experiment template.  To create a template, specify the following information:
  Targets: A target can be a specific resource in your AWS environment, or one or more
resources that match criteria that you specify, for example, resources that have specific
tags.    Actions: The actions to carry out on the target. You can specify multiple actions,
the duration of each action, and when to start each action during an experiment.    Stop
conditions: If a stop condition is triggered while an experiment is running, the experiment
is automatically stopped. You can define a stop condition as a CloudWatch alarm.   For more
information, see the AWS Fault Injection Simulator User Guide.

# Arguments
- `actions`: The actions for the experiment.
- `client_token`: Unique, case-sensitive identifier that you provide to ensure the
  idempotency of the request.
- `description`: A description for the experiment template. Can contain up to 64 letters
  (A-Z and a-z).
- `role_arn`: The Amazon Resource Name (ARN) of an IAM role that grants the AWS FIS service
  permission to perform service actions on your behalf.
- `stop_conditions`: The stop conditions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: The tags to apply to the experiment template.
- `"targets"`: The targets for the experiment.
"""
function create_experiment_template(
    actions,
    clientToken,
    description,
    roleArn,
    stopConditions;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return fis(
        "POST",
        "/experimentTemplates",
        Dict{String,Any}(
            "actions" => actions,
            "clientToken" => clientToken,
            "description" => description,
            "roleArn" => roleArn,
            "stopConditions" => stopConditions,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_experiment_template(
    actions,
    clientToken,
    description,
    roleArn,
    stopConditions,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return fis(
        "POST",
        "/experimentTemplates",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "actions" => actions,
                    "clientToken" => clientToken,
                    "description" => description,
                    "roleArn" => roleArn,
                    "stopConditions" => stopConditions,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Deletes the specified experiment template.

# Arguments
- `id`: The ID of the experiment template.

"""
function delete_experiment_template(id; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "DELETE",
        "/experimentTemplates/$(id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_experiment_template(
    id, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "DELETE",
        "/experimentTemplates/$(id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Gets information about the specified AWS FIS action.

# Arguments
- `id`: The ID of the action.

"""
function get_action(id; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "GET", "/actions/$(id)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_action(
    id, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "GET",
        "/actions/$(id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Gets information about the specified experiment.

# Arguments
- `id`: The ID of the experiment.

"""
function get_experiment(id; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "GET", "/experiments/$(id)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_experiment(
    id, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "GET",
        "/experiments/$(id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Gets information about the specified experiment template.

# Arguments
- `id`: The ID of the experiment template.

"""
function get_experiment_template(id; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "GET",
        "/experimentTemplates/$(id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_experiment_template(
    id, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "GET",
        "/experimentTemplates/$(id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the available AWS FIS actions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, make another call with the returned nextToken value.
- `"nextToken"`: The token for the next page of results.
"""
function list_actions(; aws_config::AbstractAWSConfig=global_aws_config())
    return fis("GET", "/actions"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
end
function list_actions(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "GET", "/actions", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

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

Lists your experiment templates.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, make another call with the returned nextToken value.
- `"nextToken"`: The token for the next page of results.
"""
function list_experiment_templates(; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "GET",
        "/experimentTemplates";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_experiment_templates(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "GET",
        "/experimentTemplates",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists your experiments.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, make another call with the returned nextToken value.
- `"nextToken"`: The token for the next page of results.
"""
function list_experiments(; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "GET", "/experiments"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_experiments(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "GET",
        "/experiments",
        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 the tags for the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "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 fis(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_experiment(client_token, experiment_template_id)
    start_experiment(client_token, experiment_template_id, params::Dict{String,<:Any})

Starts running an experiment from the specified experiment template.

# Arguments
- `client_token`: Unique, case-sensitive identifier that you provide to ensure the
  idempotency of the request.
- `experiment_template_id`: The ID of the experiment template.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: The tags to apply to the experiment.
"""
function start_experiment(
    clientToken, experimentTemplateId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "POST",
        "/experiments",
        Dict{String,Any}(
            "clientToken" => clientToken, "experimentTemplateId" => experimentTemplateId
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_experiment(
    clientToken,
    experimentTemplateId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return fis(
        "POST",
        "/experiments",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "clientToken" => clientToken,
                    "experimentTemplateId" => experimentTemplateId,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Stops the specified experiment.

# Arguments
- `id`: The ID of the experiment.

"""
function stop_experiment(id; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "DELETE",
        "/experiments/$(id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function stop_experiment(
    id, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "DELETE",
        "/experiments/$(id)",
        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})

Applies the specified tags to the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tags`: The tags for the resource.

"""
function tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "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 fis(
        "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)
    untag_resource(resource_arn, params::Dict{String,<:Any})

Removes the specified tags from the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tagKeys"`: The tag keys to remove.
"""
function untag_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "DELETE",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return fis(
        "DELETE",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Updates the specified experiment template.

# Arguments
- `id`: The ID of the experiment template.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"actions"`: The actions for the experiment.
- `"description"`: A description for the template.
- `"roleArn"`: The Amazon Resource Name (ARN) of an IAM role that grants the AWS FIS
  service permission to perform service actions on your behalf.
- `"stopConditions"`: The stop conditions for the experiment.
- `"targets"`: The targets for the experiment.
"""
function update_experiment_template(id; aws_config::AbstractAWSConfig=global_aws_config())
    return fis(
        "PATCH",
        "/experimentTemplates/$(id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_experiment_template(
    id, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return fis(
        "PATCH",
        "/experimentTemplates/$(id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
# NOTE: at the moment these are ad hoc implementations for each group that
#       we decide we need. eventually we can hopefully simplify and merge
#       many of these codes
#       The first version was written before NamedTupleTools (or before I
#        knew about it -> this will simplify the code here a bit)

using NamedTupleTools: delete, merge, namedtuple

abstract type SymmetryGroup end

"""
`struct NoSym <: SymmetryGroup ` : no symmetrisation other than
permutation symmetry already baked into the AA basis.
"""

"""
`struct NoSym <: SymmetryGroup` : no symmetry beyond the standard
permutation symmetry. This is currently not used, but could be incorporated
to provide a more streamlined experience for the user.
"""
struct NoSym <: SymmetryGroup
end

# write_dict(G::O3) =
#       Dict("__id__" => "ACE_O3",
#            "lsym" => lsym(G),
#            "msym" => msym(G) )

# read_dict(::Val{:ACE_O3}, D::Dict) =
#       O3(Symbol(D["lsym"]), Symbol(D["msym"]))

is_refbasisfcn(G::NoSym, AA) = true

get_sym_spec(G::NoSym, bb) = bb

function coupling_coeffs(symgrp::NoSym, bb, rotc::Rot3DCoeffs{T, TP}) where {T, TP}
   return [ one(TP) ], [bb]
end

# -----------------------  O3 SYMMETRY

# this is a prototype implemenation; eventually (asap!) we need to allow
# rotation of multiple features at once, e.g., spin-orbit coupling!

"""
`struct O3 <: SymmetryGroup` : this is the default symmetry group; describing
the action of a single O3 group on the basis.

Standard Usage:
```julia
O3()
```
will create an `O3{:l, :m}` instance, i.e. the group will expect the symbols
`:l, :m` in the relevant 1p basis.

But if the Ylm component of the 1p basis uses different symbols then one can
tell `O3` this via `O3(lsym, rsym)`. E.g. if the variable w.r.t. which we
symmetrize is a spin `s` then we might call it `O3(:ls, :ms)`. The main thing
to remember is that the symbols in the Ylm basis and in the O3 basis must
match.
"""
struct O3{LSYM, MSYM} <: SymmetryGroup
end

O3(lsym::Symbol = :l, msym::Symbol = :m) = O3{lsym, msym}()

lsym(G::O3{LSYM, MSYM}) where {LSYM, MSYM} = LSYM
getl(G::O3, b::NamedTuple) = b[lsym(G)]

msym(G::O3{LSYM, MSYM}) where {LSYM, MSYM} = MSYM
getm(G::O3, b::NamedTuple) = b[msym(G)]

write_dict(G::O3) =
      Dict("__id__" => "ACE_O3",
           "lsym" => lsym(G),
           "msym" => msym(G) )

read_dict(::Val{:ACE_O3}, D::Dict) =
      O3(Symbol(D["lsym"]), Symbol(D["msym"]))

is_refbasisfcn(G::O3, AA) = all( bi[msym(G)] == 0 for bi in AA )

get_sym_spec(G::O3, bb) = delete.(bb, (msym(G),))

function coupling_coeffs(symgrp::O3, bb, rotc::Rot3DCoeffs{T, TP}) where {T, TP}
   # bb = [ b1, b2, b3, ... ]
   # bi = (μ = ..., n = ..., l = ..., m = ...)
   #    (μ, n) -> n; only the l and m are used in the angular basis
   if length(bb) == 0
      return coco_init(rotc.phi), [bb,]
   end
   # convert to a format that the Rotations3D implementation can understand
   # this utility function splits the bb = (b1, b2 ...) with each
   # b1 = (μ = ..., n = ..., l = ..., m = ...) into
   #    l, and a new n = (μ, n)
   ll, nn = _b2llnn(symgrp, bb)
   # now we can call the coupling coefficient construiction!!
   U, Ms = rpe_basis(rotc, nn, ll)

   # but now we need to convert the m spec back to complete basis function
   # specifications (provided by sending in a prototype b = bb[1])
   rpebs = [ _nnllmm2b(symgrp, bb[1], nn, ll, mm) for mm in Ms ]

   return U, rpebs
end

function rpe_basis(A::Rot3DCoeffs,
						 nn::SVector{N, TN},
						 ll::SVector{N, Int}) where {N, TN}
	Ure, Mre = Rotations3D.re_basis(A, ll)
	G = _gramian(nn, ll, Ure, Mre)
   S = svd(G)
   rk = rank(Diagonal(S.S); rtol =  1e-7)
	Urpe = S.U[:, 1:rk]'
	return Diagonal(sqrt.(S.S[1:rk])) * Urpe * Ure, Mre
end

function _gramian(nn, ll, Ure, Mre)
   N = length(nn)
   nre = size(Ure, 1)
   G = zeros(Complex{Float64}, nre, nre)
   for σ in permutations(1:N)
      if (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mre), (iU2, mm2) in enumerate(Mre)
         if mm1[σ] == mm2
            for i1 = 1:nre, i2 = 1:nre
               G[i1, i2] += coco_dot(Ure[i1, iU1], Ure[i2, iU2])
            end
         end
      end
   end
   return G
end

# TODO: replace all this awful code with NamedTupleTools

_nnllmm2b(G, b, nn, ll, mm) = [ _nlm2b(G, b, n, l, m) for (n, l, m) in zip(nn, ll, mm) ]

@generated function _nlm2b(G::O3{LSYM, MSYM}, b::NamedTuple{ALLKEYS},
                           n::NamedTuple{NKEYS},
                           l, m) where {LSYM, MSYM, ALLKEYS, NKEYS}
   code =
      ( "( _b = (" * prod("$(k) = n.$(k), " for k in NKEYS)
               * "$(LSYM) = l, $(MSYM) = m ); "
         *
        " b = (" * prod("$(k) = _b.$(k), " for k in ALLKEYS) * ") )" )
   :( $(Meta.parse(code)) )
end

function _b2llnn(G::O3, bb)
   @assert all( iszero(b[msym(G)]) for b in bb )
   ll = SVector( [b[lsym(G)] for b in bb]... )
   nn = SVector( [_all_but_lm(G, b) for b in bb]... )
   return ll, nn
end

"""
return a NamedTuple containing all values in b except those corresponding
to l and m keys

TODO: get rid of this ridiculousness and replace with NamedTupleTools methods
"""
@generated function _all_but_lm(G::O3{LSYM, MSYM},
                                b::NamedTuple{NAMES}) where {LSYM, MSYM, NAMES}
   code = "n = ("
   for k in NAMES
      if !(k in (LSYM, MSYM))
         code *= "$(k) = b.$(k), "
      end
   end
   code *= ")"
   quote
      $(Meta.parse(code))
      n
   end
end

# -------------- O3 ⊗ O3

# this is a preliminary implementation; eventually we may want a more
# general description composition of arbitrary isometry combinations

@doc raw"""
`struct O3O3 <: SymmetryGroup` : This type implements the ``O(3) \otimes O(3)`` symmetry
group. This is useful when a particle has two euclidean vector attributes, say
``{\bm r}`` and ``{\bm s}`` and the action of the group on the pair is
```math
   (Q_r, Q_s)[ (\boldsymbol{r}, \boldsymbol{s}) ] = (Q_r \boldsymbol{r}, Q_s \boldsymbol{s})
```
A canA canonical application is magnetism: it is known that spin-orbit coupling
is a very weak effect. By ignoring it, i.e., letting positions and spins rotate
independently of one another, one makes a small modelling error. This leads
precisely to the ``O(3) \otimes O(3)`` symmetry.

To construct this group, use
```julia
symgrp = O3(:lr, :mr) ⊗ O3(:ls, :ms)
```
or replace those symbols with the appropriate symbols used to specify the
corresponding `Ylm1pbasis` objects.
"""
struct O3O3{LSYM1, MSYM1, LSYM2, MSYM2} <: SymmetryGroup
   G1::O3{LSYM1, MSYM1}
   G2::O3{LSYM2, MSYM2}
end

import Base: kron

function kron(G1::O3, G2::O3)
   @assert lsym(G1) != lsym(G2)
   @assert msym(G1) != msym(G2)
   return O3O3(G1, G2)
end

⊗(G1::O3, G2::O3) = kron(G1, G2)
export ⊗

write_dict(G::O3O3) =
      Dict("__id__" => "ACE_O3O3",
           "G1" => write_dict(G.G1),
           "G2" => write_dict(G.G2) )

read_dict(::Val{:ACE_O3O3}, D::Dict) =
      read_dict(D["G1"]) ⊗ read_dict(D["G2"])

is_refbasisfcn(G::O3O3, AA) = all( bi[msym(grp)] == 0
                                   for bi in AA, grp in (G.G1, G.G2) )

get_sym_spec(G::O3O3, bb) = delete.(bb, Ref( (msym(G.G1), msym(G.G2)) ))

function coupling_coeffs(symgrp::O3O3, bb, rotc::Rot3DCoeffs)
   # bb = [ b1, b2, b3, ... ]
   # bi = (μ = ..., n = ..., l1 = ..., m1 = ..., l2 = ..., m2 = ...)
   #    (μ, n, ...) -> n; only the l and m are used in the angular basis
   if length(bb) == 0
      return coco_init(rotc.phi), [bb,]
   end

   # the prototype namedtuple describing a single 1p basis fcn
   PROTOTUPLE = prototype(bb[1])
   NU = length(bb)

   # convert to (nn, ll, mm) format for Rotations3D
   ll1, ll2, nn, ll12 = _b2llnn(symgrp, bb)
   # ... and construct the coupling coefficients for the individual subgroups
   U1, M1 = Rotations3D.re_basis(rotc, ll1)
   U2, M2 = Rotations3D.re_basis(rotc, ll2)

   nU1, nM1 = size(U1)
   nU2, nM2 = size(U2)
   @assert nM1 == length(M1)
   @assert nM2 == length(M2)
   UT = promote_type(eltype(U1), eltype(U2))

   # there is admissible combination:
   if nU1 == 0 || nU2 == 0
      return UT[], SVector{NU, PROTOTUPLE}[]
   end

   # now combine them into the effective coupling coeffs and combined Ms
   # each column Ure[:, i] corresponds to one rotation-invariant basis fcn
   M1M2TUPLE = namedtuple(msym(symgrp.G1), msym(symgrp.G2))
   Mre = Vector{typeof(Vector(M1M2TUPLE.(M1[1], M2[1])))}(undef, nM1 * nM2)
   Ure = zeros( UT, (nU1 * nU2, nM1 * nM2) )
   jdx = 0
   for j1 = 1:nM1, j2 = 1:nM2
      jdx += 1
      Mre[jdx] = Vector(M1M2TUPLE.(M1[j1], M2[j2]))
      idx = 0
      for i1 = 1:nU1, i2 = 1:nU2
         idx += 1
         Ure[idx, jdx] = U1[i1, j1] * U2[i2, j2]
      end
   end

   # insert another reduction step
   # in my tests this never reduces the size, but I haven't had the time
   # to actually prove it isn't needed, so we will keep it for now
   Gre = [ sum(coco_dot.(Ure[i1, :], Ure[i2, :])) for i1 = 1:size(Ure, 1), i2 = 1:size(Ure, 1) ]
   Sre = svd(Gre)
   rk = rank(Diagonal(Sre.S); rtol =  1e-7)
   Ure = Sre.U[:, 1:rk]' * Ure

   # now symmetrize w.r.t. permutations
   G = _gramian(nn, ll12, Ure, Mre)
   S = svd(G)
   rk = rank(Diagonal(S.S); rtol =  1e-7)
   Urpe = S.U[:, 1:rk]'
   U = Diagonal(sqrt.(S.S[1:rk])) * Urpe * Ure

   # reconstruct the basis function specifications
   rpebs = [ PROTOTUPLE.(merge.(nn, ll12, mm12)) for mm12 in Mre ]

   return U, rpebs
end

function _b2llnn(G::O3O3{L1, M1, L2, M2}, bb) where {L1, M1, L2, M2}
   N = length(bb)
   @assert all( iszero(b[M]) for b in bb, M in (M1, M2) )
   ll1 = ntuple(i -> bb[i][L1], N)  |> SVector
   ll2 = ntuple(i -> bb[i][L2], N)  |> SVector
   nn  = ntuple(i -> delete(bb[i], (L1, M1, L2, M2)), N)
   ll12 = ntuple(i -> select(bb[i], (L1, L2)), N)
   return ll1, ll2, nn, ll12
end

ext: jl


text:
import Einsum
import Benchmarks

let
	"""
	Common fate bechmark proposed by @faroit
	"""
	function commonfate(A,H,C)
	  Einsum.@einsum P[a,b,f,t,c] := A[a,b,f,j]*H[t,j]*C[c,j]
	  return nothing
	end

	function assess_benchmark(dim=30, r=10)
		A = randn(dim,dim,dim,r)
	    H = randn(dim,r)
	    C = randn(dim,r)

	    result = Benchmarks.@benchmark commonfate(A,H,C)
	    stat = Benchmarks.SummaryStatistics(result)

	    # convert ns -> seconds
	    t = stat.elapsed_time_center / 1_000_000_000.0

	    return (dim,t)
	end

	assess_benchmark()
end


text:
"""
    load_gid_mesh(filename) -> mesh
"""
function load_gid_mesh(fn::AbstractString)
    open(fn, "r") do io
        load_gid_mesh(io)
    end
end

"""
    load_gid_mesh(stream) ->mesh
"""
function load_gid_mesh(file)

    lineCount = countlines(file);
    seekstart(file);
    readline(file); # MESH dimension 3 ElemType Triangle Nnode 3
    readline(file); # Coordinates

    # Read the vertices
    vertices = Vector{Pt{3,Float64}}(undef,lineCount)
    vertexCount = 0;
    while !eof(file)
            line = split(readline(file));
            size(line, 1) == 4 || break; # End Coordinates
            vertexCount += 1;
            #vertices[vertexCount] = Pt{3,Float64}(map(Meta.parse, line[2:4]));
            vertices[vertexCount] = Pt{3,Float64}(parse.(Float64,line[2:4]))
    end

    readline(file); # "\n"
    readline(file); # Elements

    # Read the triangles
    triangles = Vector{Pt{3,Int}}(undef,lineCount-vertexCount);
    triangleCount = 0;
    while !eof(file)
            line = split(readline(file));
            size(line, 1) == 4 || break; # End Elements
            triangleCount += 1;
            triangles[triangleCount] = Pt{3,Int}(map(Meta.parse, line[2:4]));
    end

    Mesh(vertices[1:vertexCount], triangles[1:triangleCount])
end


text:
export eikonal_const_vel

"""
`eikonal_const_vel(src, z, x, [y, ] velocity)`

Compute traveltime tables between `src` locations and subsurface model locations given by `z`, `x`, and optionally `y`. It assumes `velocity` is constant.

Parameters
----------
* `src` : `(ns, n)`, `AbstractMatrix{<:Real}`

The `(z_s, x_s [, x_y])` coordinates for every source. In 2D, `n = 2` and in 3D `n = 3`.

* `z` : `(nz,)`, `AbstractVector{<:Real}`

Model depth `z` coordinates.

* `x` : `(nx,)`, `AbstractVector{<:Real}`

Model horizontal `x` coordinates.

* `y` : `(ny,)`, `AbstractVector{<:Real}`, optional

Model horizontal `y` coordinates.

* `velocity` : `Real`

Wave speed

Returns
-------

* `trav` : `(nz, nx, [ny, ] ns)`, `AbstractArray{<:Real, M}`

Traveltime between each src and locations given by `z`, `x` and optionally `y`:

```math
t(z, x, y, s) = \\frac{\\sqrt{(z-s_z)^2 + (x-s_x)^2 + (y-s_y)^2}}{v}
```
"""
function eikonal_const_vel(src::AbstractMatrix{<:Real},
                           z::AbstractVector{<:Real},
                           x::AbstractVector{<:Real},
                           velocity::Real)
    if velocity < 0
        @warn("Velocity cannot not be negative")
        throw(DomainError())
    end
    size(src, 2) == 2 || throw(DimensionMismatch("Second dimension of `src` must match number of model arguments"))

    ns = size(src, 1)
    trav = zeros(length(z), length(x), ns)
    Z = [i for i in z, j in x]
    X = [j for i in z, j in x]
    for is in 1:ns
        src_z, src_x = src[is, :]
        trav[:,:,is] = @. sqrt.((Z - src_z)^2 + (X - src_x)^2)/velocity
    end
    return trav
end
function eikonal_const_vel(src::AbstractMatrix{<:Real},
                           z::AbstractVector{<:Real},
                           x::AbstractVector{<:Real},
                           y::AbstractVector{<:Real},
                           velocity::Real)
    if velocity < 0
        @warn("Velocity cannot not be negative")
        throw(DomainError())
    end
    size(src, 2) == 3 || throw(DimensionMismatch("Second dimension of `src` must match number of model arguments"))

    ns = size(src, 1)
    trav = zeros(length(z), length(x), length(y), ns)
    Z = [i for i in z, j in x, k in y]
    X = [j for i in z, j in x, k in y]
    Y = [k for i in z, j in x, k in y]
    for is in 1:ns
        src_z, src_x, src_y = src[is, :]
        trav[:,:,:,is] = @. sqrt((Z - src_z)^2 + (X - src_x)^2 + (Y - src_y)^2)/velocity
    end
    return trav
end


text:
#========================================================================================#
"""
	IN504Boids

Module Boids: Reynolds' Boids model of flocking in birds.

Author: Niall Palfreyman, 27/04/22
"""
module Boids

export demo									                # Externally available names
using Agents, LinearAlgebra, InteractiveDynamics, GLMakie   # Required packages

#-----------------------------------------------------------------------------------------
# Module definitions:

#-----------------------------------------------------------------------------------------
"""
	Boid

To construct a new agent implementation, use the @agent macro and add all necessary extra
attributes. Again, id, pos and vel are automatically inserted.
"""
@agent Boid ContinuousAgent{2} begin
    speed::Float64						# Speed of boid
    cohering::Float64					# Extent to which boids want to stay together
    proximity::Float64					# Uncomfortable proximity distance
    separating::Float64					# Extent to which boids move out of proximity
    aligning::Float64					# Extent to which boids align their flight
    sight::Float64						# Distance within which I notice other boids
end

#-----------------------------------------------------------------------------------------
const boid_polygon = Polygon(Point2f[(-0.5, -0.5), (1, 0), (-0.5, 0.5)])
"""
	boid_marker(boid)

The graphical marker for boids
"""
function boid_marker(b::Boid)
    φ = atan(b.vel[2], b.vel[1])
    scale(rotate2D(boid_polygon, φ), 2)
end

#-----------------------------------------------------------------------------------------
# Module methods:

#-----------------------------------------------------------------------------------------
"""
	boidworld()

Create the Boids model.
"""
function boidWorld(;
	n_boids				= 100,
	speed				= 1.0,
	cohering			= 0.25,
	proximity			= 4.0,
	separating			= 0.25,
	aligning			= 0.01,
	sight				= 5.0,
	extent				= (100, 100),
	spacing				= sight / 1.5,
)
	world = ABM( Boid, ContinuousSpace(extent,spacing), scheduler = Schedulers.randomly)
	for _ in 1:n_boids
		add_agent!( world,
			Tuple(2*rand(world.rng, 2) .- 1),				# Random facing direction
			speed,
			cohering,
			proximity,
			separating,
			aligning,
			sight,
		)
	end
	return world
end

#-----------------------------------------------------------------------------------------
"""
    agent_step!( boid, world)

Define how a boid moves within the Boid world, depending on its propensity for cohering,
separating and aligning.
"""
function agent_step!( boid, world)
	idNbrs = nearby_ids(boid, world, boid.sight)		# All neighbour ids within my sight

	N = 0												# Number of neighbours
	heading = separating = cohering = (0.0, 0.0)		# Initialise flight parameters
	for id in idNbrs
		# Calculate neighbour's position and bearing:
		N += 1											# Increment number of neighbours
		position = world[id].pos						# Calculate neighbour's position
		bearing = position .- boid.pos					# Calculate neighbour's bearing from here

		cohering = cohering .+ bearing					# 'cohering' accumulates bearing of neighbours
		if edistance(boid.pos, position, world) < boid.proximity
			# 'separating' accumulates recoil from uncomfortably near neighbours:
			separating = separating .- bearing
		end
		heading = heading .+ world[id].vel				# 'align' accumulates heading of neighbours
	end

	# Calculate average values from accumulations:
	N				= max(N,1)								# Allow for no neighbours in sight
	cohering		= cohering .* boid.cohering ./ N		# Cohere along average bearning
	separating		= separating .* boid.separating ./ N	# Separate along average bearing
	heading			= heading .* boid.aligning ./ N			# Align with average heading

	# Compute velocity based on the three parameters cohering, separating, aligning:
	boid.vel = (boid.vel .+ cohering .+ separating .+ heading) ./ 2
	boid.vel = boid.vel ./ norm(boid.vel)

	# Move boid according to new velocity and speed
	move_agent!(boid, world, boid.speed)
end

#-----------------------------------------------------------------------------------------
function demo()
    bworld = boidWorld()							# Create a Boid world ...

    abmvideo(										# ... then make a pretty video of it.
        "boids.mp4", bworld, agent_step!;
        am = boid_marker,
        framerate = 20,
        frames = 100,
        title = "Boids"
    )
end

end		# ... of module Boids


text:
# validation.jl: contains all functions enabling a comparison with the exact polynomials
# About
#   Author       - Peter Opsomer (Peter.Opsomer@cs.kuleuven.be)
#   History      - Created October 2014, last edit October 2015

# R_JACOBI Recurrence coefficients for monic Jacobi polynomials.
#    ab=R_JACOBI(n,a,b) generates the first n recurrence
#    coefficients for monic Jacobi polynomials with parameters
#    a and b. These are orthogonal on [-1,1] relative to the
#    weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
#    are stored in the first column, the n beta-coefficients in
#    the second column, of the nx2 array ab. The call ab=
#    R_JACOBI(n,a) is the same as ab=R_JACOBI(n,a,a) and
#    ab=R_JACOBI(n) the same as ab=R_JACOBI(n,0,0).
#    Supplied by Dirk Laurie, 6-22-1998; edited by Walter Gautschi, 4-4-2002.
# About
#   Translated from https://www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
function r_jacobi(N,a = 0,b = a)
if((N<=0)|(a<=-1)|(b<=-1))
	error("parameter(s) out of range")
end
nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1 return [nu mu]; end
N=N-1; n=linspace(1,N,N); nab=2*n+a+b;
A=[nu; (b^2-a^2)./(nab.*(nab+2))];
n=2:N; nab=nab[n];
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
return [A'; [mu B1 B']]';
end

# GAUSS Gauss quadrature rule.
#    Given a weight function w encoded by the nx2 array ab of the
#    first n recurrence coefficients for the associated orthogonal
#    polynomials, the first column of ab containing the n alpha-
#    coefficients and the second column the n beta-coefficients,
#    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
#    the n-point Gauss quadrature rule for the weight function w.
#    The nodes, in increasing order, are stored in the first
#    column, the n corresponding weights in the second column, of
#    the nx2 array xw.
#
# About
#   Translated from https://www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
function gauss(N,ab)
N0 = size(ab,1);
if N0 < N
	error("input array ab too short")
end
J=zeros(N,N);
for n=1:N
	J[n,n] = ab[n,1];
end
for n=2:N
  J[n,n-1] = sqrt(ab[n,2]);
  J[n-1,n] = J[n,n-1];
end
(D,V)=eig(J);
I=sortperm(D); #Actually already sorted..
V=V[:,I];
return (D, ab[1,2]*V[1,:]'.^2)
end

# Real orthonormal polynomials using the recurrence relation and the OPQ-library
# of Walter Gautschi. These polynomials exhibit degrading accuracy from n about
# 128 to about a relative error of 1e-8 around n=2^12, which can be checked by
# evaluating exactPolys(0,0,@(x) 1,4096): then P(1,n+1) should be sqrt((2*n+1)/2).
# Input
#   alpha,beta,h - Parts of the weight w(x) = (1-x).^alpha.*(1+x).^beta.*h(x)
#   N            - Maximum degree
# Output
#   P            - P(x,n) giving orthonormal poly of degree n at point x
#                  as a function, ||pi(x,n)||_w = 1
#   gammaP       - gammaP(n+1)*2^n = \gamma_n = is the leading order coefficient
#		   of the orthonormal polynomial of degree n. Multiplication by
#		   2^(-n) to avoid overflow of saved value gammaP(n+1)
#   alphaP,betaP - alphaP(n+1) & betaP(n+1) are the recurrence coefficients of
#		   the monic orthogonal polynomial of degree n
# Note: This uses another algorithm than the Matlab version.
function exactPolys(alpha,beta,h,N)
# Use alpha and beta to construct quadrature rules for computing the
# real alphaP & betaP. ab(:,1) != alphaP because it doesn't use h(x)
ab = r_jacobi(3*N,alpha,beta);
(z, w) = gauss(3*N,ab);
z = (1+0im)*z
w = (1+0im)*w
function inwprod(f,g)
	sum(w.*h(z+0im).*f(z+0im).*conj(g(z+0im) ) );
end
alphaP = zeros(N+1,1);
betaP = zeros(N+1,1);
gammaP = zeros(N+2,1);
gammaP[1] = 1/sqrt(inwprod(x -> ones(x), x -> ones(x) ) );
alphaP[1] = inwprod(x -> x, x -> ones(x) )*gammaP[1]^2;
function mon(x,n,alphaP,betaP)
	pin = [(1+1im)*zeros(length(x),1) (1+0im)*ones(length(x),1)];
	for j=1:n
	    pin = [pin[:,2] (x-alphaP[j]+0im).*pin[:,2]-betaP[j].*pin[:,1]];
	end
	return pin[:,2]
end
gammaP[2] = 1/sqrt(inwprod(x -> mon(x,1,alphaP,betaP), x -> mon(x,1,alphaP,betaP) ) );
start = time();
prevToc = start;
for k = 1:N
        if time()-prevToc > 1.0
		prevToc = time()
		print("OPQ iteration ", k, ", ", round((k-1)/N,3), "% done, estim sec left = ", (prevToc-start)*(N-k+1.0)/(k-1.0), '\n');
        end
	alphaP[k+1] = inwprod(x -> x.*mon(x,k,alphaP,betaP), x -> mon(x,k,alphaP,betaP) )*gammaP[k+1]^2;
        betaP[k+1] = (gammaP[k]/gammaP[k+1] )^2;
        gammaP[k+2] = 1/sqrt(inwprod(x -> mon(x,k+1,alphaP,betaP), x -> mon(x,k+1,alphaP,betaP) ) );
end
function P(x,n) return (mon(x,n,alphaP,betaP)*gammaP[n+1])[1]; end
return (P,gammaP,alphaP,betaP)
end

# Plot the convergence results with LS convergence in loglog, in 1 figure.
# Input
#   r     - Actual values
#   a     - Approximation
#   t     - Title of the plot
#   shift - Degree shift with respect to 2.0.^(1:maxP2)
function plotConv(r,a,t,shift=0)
(maxP2,maxOrder) = size(a);
data =  abs(a - repmat(r,1,maxOrder) )./abs(repmat(r,1,maxOrder) );
colors = "bgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmyk";
ns = 2.0.^(1:maxP2)+shift;
for i = 1:maxOrder
	idx = maxP2;
	while data[idx,i] < 0.1^11
		idx = idx-1;
	end
	if i == 1
	        loglog(ns,data[:,i],string("o",colors[1]),label="1 term")
	else
	        loglog(ns,data[:,i],string("o",colors[i]),label=string(i," terms") )
	end
	loglog(ns,data[idx,i]./((ns./ns[idx]).^i), string(colors[i],"-") );
end
legend(loc=3);
xlabel("n");
ylabel("Relative error");
title(t);
show()
end


text:
# acc plotting
using PGFPlots
using QHull
using JLD2
using LazySets
using LinearAlgebra

include("models/problems.jl")
include("OverApprox/src/OA_relational_util.jl")
include("MIP/src/overt_to_mip.jl")
include("MIP/src/mip_utils.jl")
include("models/acc/acc.jl")

JLD2.@load "examples/jmlr/data/acc/acc_reachability_data_1step_55.jld2"
one_step_state_sets = symbolic_state_sets
one_step_meas_sets = symbolic_meas_sets

JLD2.@load "examples/jmlr/data/acc/acc_reachability_data_55.jld2"

# monte carlo simulate
mc_state_sets, xvec, x0, mc_meas_sets, yvec, y0 = monte_carlo_simulate(query, input_set, n_sim=1000000)

# clean up sets
init, reachable_state_sets = clean_up_sets(concrete_state_sets, symbolic_state_sets, concretization_intervals)

reachable_meas_sets = clean_up_meas_sets(concrete_meas_sets, symbolic_meas_sets, concretization_intervals)

# plot styles
concrete_line_style = "solid, black, mark=none"
symbolic_line_style = "solid, red, very thick, mark=none"
mc_line_style = "solid, blue, ultra thick, mark=none"
############################################################
# First plot: rigorous plot of meas w/ unsafe=10
############################################################
fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$y\$", xlabel="timesteps", title="Adaptive Cruise Control, Measurement Sets")
conc_meas_sets = get_interval_subsets(one_step_meas_sets, [1])
sym_meas_sets = get_interval_subsets(reachable_meas_sets, [1])
mc_meas_sets = get_interval_subsets(mc_meas_sets, [1])
push!(fig, PGFPlots.Plots.Linear([-1., length(mc_meas_sets) + 1.], [10., 10.], style="red, dashed", mark="none", legendentry="Unsafe Set"))
for t in 1:query.ntime
    timestep = Float64[t,t]
    if t == query.ntime
        push!(fig, PGFPlots.Plots.Linear(timestep, conc_meas_sets[t], style=concrete_line_style, legendentry="Concrete Sets"))
        push!(fig, PGFPlots.Plots.Linear(timestep, sym_meas_sets[t], style=symbolic_line_style, legendentry="OVERT Hybrid Symbolic Sets"))
        push!(fig, PGFPlots.Plots.Linear(timestep, mc_meas_sets[t], style=mc_line_style, legendentry="Monte Carlo Simulations"))
    else
        push!(fig, PGFPlots.Plots.Linear(timestep, conc_meas_sets[t], style=concrete_line_style))
        push!(fig, PGFPlots.Plots.Linear(timestep, sym_meas_sets[t], style=symbolic_line_style))
        push!(fig, PGFPlots.Plots.Linear(timestep, mc_meas_sets[t], style=mc_line_style))
    end
end
fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

save("examples/jmlr/plots/acc/acc_meas.tex", fig)
save("examples/jmlr/plots/acc/acc_meas.pdf", fig)

############################################################
### Second plot: car positions SYMBOLIC ONLY
############################################################
# in tex, replace with: #ab0000
define_color("lead_color", 0xffc2ca)
input_lead_style = "dashed, lead_color, mark=none"
sym_lead_style = "solid, lead_color, ultra thick, mark=none"
# ego in tex, replace with #d84a43
define_color("ego_color", 0x800000)
input_ego_style = "dashed, ego_color, mark=none"
sym_ego_style = "solid, ego_color, ultra thick, mark=none"
# d safe, in tex replace with #ff8080
define_color("dsafe_color", 0xc76666)
input_dsafe_style = "dashed, dsafe_color, mark=none"
sym_dsafe_style = "densely dotted, dsafe_color, very thick, mark=none"

fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="Distance Along Roadway (m)", xlabel="timesteps", title="Adaptive Cruise Control, State Sets")
dim = [1] # x_lead
sym_lead_sets = get_interval_subsets(reachable_state_sets, dim)
dim = [4] # x_ego
sym_ego_sets = get_interval_subsets(reachable_state_sets, dim)
dim = [5] # v_ego
sym_dsafe_sets = get_interval_subsets(reachable_state_sets, dim) .* 1.4 .+ Ref([10., 10.])
# construct dsafe sets extending from behind lead car position
sym_dsafe_behind_ego = [[minimum(sym_lead_sets[i]) - maximum(sym_dsafe_sets[i]), minimum(sym_lead_sets[i])] for i in 1:length(sym_dsafe_sets)]

# plot init sets
dims=[1] # lead
lead_input = [low(input_set)[dims]..., high(input_set)[dims]...]
push!(fig, PGFPlots.Plots.Linear([0.,0.], lead_input, style=input_lead_style, legendentry="Lead Position Input Set"))
##############
dims=[5] # v_ego
d_safe_i = [low(input_set)[dims]..., high(input_set)[dims]...].*1.4 .+ 10.
dsafe_input = [minimum(lead_input) - maximum(d_safe_i), minimum(lead_input)]
push!(fig, PGFPlots.Plots.Linear([0.,0.], dsafe_input, style=input_dsafe_style, legendentry="Minimum Safe Distance Input Set"))
##############
dims=[4] # ego
ego_input = [low(input_set)[dims]..., high(input_set)[dims]...]
push!(fig, PGFPlots.Plots.Linear([0.,0.], ego_input, style=input_ego_style, legendentry="Ego Position Input Set"))

for t in 1:query.ntime
    timestep = Float64[t,t]
    if t == query.ntime
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_lead_sets[t], style=sym_lead_style, legendentry="Lead Position, OVERT Hybrid Symbolic"))
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_dsafe_behind_ego[t], style=sym_dsafe_style, legendentry="Minimum Safe Distance, OVERT Hybrid Symbolic"))
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_ego_sets[t], style=sym_ego_style, legendentry="Ego Position, OVERT Hybrid Symbolic"))
    else
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_lead_sets[t], style=sym_lead_style))
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_dsafe_behind_ego[t], style=sym_dsafe_style))
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_ego_sets[t], style=sym_ego_style))
    end
end
fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

save("examples/jmlr/plots/acc/acc.tex", fig)
save("examples/jmlr/plots/acc/acc.pdf", fig)

############################################################
### Third plot: 2D plots of either velocity or accel
############################################################
concrete_set_style = "solid, black, mark=none, fill=black"
symbolic_set_style = "solid, ego_color, mark=none, fill=ego_color"
mc_set_style = "solid, blue, mark=none, fill=blue"
concrete_set_trans = concrete_set_style*", fill opacity=0.5"
symbolic_set_trans = symbolic_set_style*", fill opacity=0.5"
mc_set_trans = mc_set_style*", fill opacity=0.5"

dims=[2,5] # velocity
fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$v_{lead}~(m/s)\$", xlabel="\$v_{ego}~(m/s)\$", title="ACC Velocity Reachable Sets")

# plot init set in both concrete and hybrid colors
inputx, inputy = get_rectangle(input_set, dims)
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=concrete_set_style, legendentry="Concrete Sets"))
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=symbolic_set_style, legendentry="OVERT Hybrid Symbolic Sets"))
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=mc_set_style, legendentry="Monte Carlo Simulations"))

for t in 1:query.ntime
    if t == query.ntime
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=concrete_set_style))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=symbolic_set_style))
        # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style_transparent))
        # plot mc_sim points
        points = xvec[:, t, dims]
        #plot hull
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_set_style))
    else
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=concrete_set_trans))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=symbolic_set_trans))
        # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style_transparent))
        # plot mc simulations hull
        points = xvec[:, t, dims]
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_set_style))
    end
end

fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/acc/acc_x$(dims[1])$(dims[2]).tex", fig)
PGFPlots.save("examples/jmlr/plots/acc/acc_x$(dims[1])$(dims[2]).pdf", fig)


text:
function fit!(
    loss::AbstractLoss,
    triplets,
    X::Embedding;
    verbose::Bool = true,
    print_every::Int = 10,
    max_iterations::Int64 = 1000
)

    @assert max_iterations >= 10 "Iterations must be at least 10"
    # @assert maximum([maximum(t) for t in triplets]) == nitems(X)

    R = fill(Inf, max_iterations) # risk
    ∇R = zeros(Float64, size(X))  # gradient

    tolerance = 1e-7 # convergence tolerance
    η = 1.0          # learning rate
    best_X = X       # best embedding found so far

    # Perform main iterations
    i = 1   # Number of iterations (so far)
    lr_increments = 0   # Number of learning rate increments
    pviolations = 0.0 # Percentage of violations

    while i < max_iterations && lr_increments < 5
        i += 1

        # Calculate gradient descent and risk
        R[i], ∇R = gradient(loss, triplets, X)
        if any(isnan.(∇R))
            break
        end

        # Update the embedding according to the gradient
        X.X = X.X .- (η / length(triplets) * nitems(X)) * ∇R

        if R[i-1] < minimum(R)
            best_X = X
        end

        # Update learning rate
        if R[i-1] > R[i] + tolerance
            lr_increments = 0
            η *= 1.01
        else
            lr_increments += 1
            η *= 0.5
        end

        # Print out progress
        if verbose && (i % print_every == 0)
            # If we have a large number of triplets, computing the number of violations
            # can be costly. Therefore, we only perform this operation every 10 iterations
            # If more details are needed, you can set the environment variable
            # JULIA_DEBUG=TripletEmbeddings before starting Julia.
            pviolations = percent_violations(triplets, X)
            @info @sprintf "loss = %s, iteration = %d, risk = %.6f, misclassifications = %.2f%%, ||∇ₓR|| = %.6f" typeof(loss) i R[i]/length(triplets) 100 * pviolations norm(∇R/length(triplets))
        end

        @debug "Iteration = $i, risk = $(R[i]), lr_increments = $lr_increments" X.X percent_violations(triplets, X)
    end

    if !verbose
        pviolations = percent_violations(triplets, X)
    end

    X = best_X
    return pviolations
end

function percent_violations(triplets::Triplets, X::Embedding)
    D = pairwise(SqEuclidean(), X.X, dims=2)

    nviolations = 0

    for t in triplets
        nviolations += D[t[:i], t[:j]] > D[t[:i], t[:k]]
    end

    return nviolations/length(triplets)
end

function percent_violations(triplets::LabeledTriplets, X::Embedding)
    D = pairwise(SqEuclidean(), X.X, dims=2)

    nviolations = 0

    for t in triplets
        nviolations += (D[t[:i], t[:j]] > D[t[:i], t[:k]] && t[:y] == -1) || ((D[t[:i], t[:j]] < D[t[:i], t[:k]] && t[:y] == +1))
    end

    return nviolations/length(triplets)
end


text:
####
#	We really need partial evaluation !!!
####

function partialeval(m::IdentityModel, ds::ArrayNode, skipnode)
    ds === skipnode && return(m, skipnode, true)
    (m, skipnode, false)
end

function partialeval(m::ArrayModel, ds::ArrayNode, skipnode)
    ds === skipnode && return(m, skipnode, true)
    (identity_model(), m(ds), false)
end

# (m::BagModel)(x::WeightedBagNode{<: AbstractNode}) = m.bm(m.a(m.im(x.data), x.bags, x.weights))

function partialeval(m::BagModel, ds::BagNode, skipnode)
    ds === skipnode && return(m, skipnode, true)
    im, ids, keep = partialeval(m.im, ds.data, skipnode)
    if keep
        return (BagModel(im, m.a,  m.bm), BagNode(ids, ds.bags, ds.metadata), true)
    end
    (identity_model(), m.bm(m.a(ids, ds.bags)), false)
end

function partialeval(m::BagModel, ds::WeightedBagNode, skipnode)
    ds === skipnode && return(m, skipnode, true)
    im, ids, keep = partialeval(m.im, ds.data, skipnode)
    if keep
        return (BagModel(im, m.a,  m.bm), WeightedBagNode(ids, ds.bags, ds.weights, ds.metadata),  true)
    end
    (identity_model(), m.bm(m.a(ids, ds.bags, ds.weights)), false)
end

function Mill.partialeval(m::LazyModel{N}, ds::LazyNode{N}, skipnode) where {N}
    ds === skipnode && return(m, skipnode, true)
    (identity_model(), m(ds), false)
end

function partialeval(m::ProductModel{M}, ds::ProductNode{P}, newnode) where {P<:NamedTuple, M<:NamedTuple}
    ks = keys(m.ms)
    mods = map(ks) do k
        partialeval(m.ms[k], ds.data[k], newnode)
    end
    ms = map(f -> f[1], mods)
    dd = map(f -> f[2], mods)
    if any(f[3] for f in mods)
        return (ProductModel((;zip(ks, ms)...), m.m), ProductNode((;zip(ks, dd)...), ds.metadata), true)
    end
    (identity_model(), m.m(vcat(dd...)), false)
end

function partialeval(m::ProductModel{M}, ds::ProductNode{P}, newnode) where {P<:Tuple, M<:Tuple}
    mods = map(1:length(m.ms)) do k
        partialeval(m.ms[k], ds.data[k], newnode)
    end
    ms = map(f -> f[1], mods)
    dd = map(f -> f[2], mods)
    if any(f[3] for f in mods)
        return (ProductModel(tuple(ms...), m.m), ProductNode(tuple(dd...), ds.metadata), true)
    end
    (identity_model(), m.m(vcat(dd...)), false)
end


text:
using GR
using LaTeXStrings

formula = L"- \frac{{\hbar ^2}}{{2m}}\frac{{\partial ^2 \psi (x,t)}}{{\partial x^2 }} + U(x)\psi (x,t) = i\hbar \frac{{\partial \psi (x,t)}}{{\partial t}}"

selntran(0)

settextfontprec(232, 3)
settextalign(2, 3)
setbordercolorind(6)

mathtex(0.5, 0.5, formula)

tbx, tby = inqmathtex(0.5, 0.5, formula)
Δx = (tbx[2] - tbx[1]) / √2
Δy = (tby[3] - tby[1]) / √2
path([0.5+Δx, Δx, 2π], [0.5, Δy, 0], "MAs");

updatews()


text:
using AsyPlots
using Test

# 2D
Pen()
Arrow()
Arrow3()
Arrow(5)
Arrow(name="MidArrow",size=3)
Point(1,2,label="hey",pen=Pen(linewidth=3))
Path([0 0; 1 1];arrow=Arrow())
Circle((0,0),1)
Polygon([im,1,0])

Point(1,2,3,label="hey",pen=Pen(linewidth=3))
Path([0 0 0; 1 1 1])
Polygon([0 0 0; 1 1 1; 0 0 1])

z = randn(5,5)
S = Surface(z)

P = Point(1,2,color="DarkGreen")
Q = Path([0 0; 0 1; 1 1; 0 2],spline=true,label="Label(\"A\",Relative(0.8))")
C = Circle((0,0),1.0,linewidth=3,opacity=0.2,color="Purple")

X = collect(range(0,stop=1,length=100))
Y1 = cumsum(randn(100))
Y2 = cumsum(randn(100))
plot(X,Y1)
plot(X,[Y1,Y2])
plot(Y1)
plot(Y1,linewidth=3,opacity=0.4)
plot([x->x^n for n=1:4],0,1)
plot(x->sin(x),0,2π)
plot(t->cos(t)^3,(0,π/2))
plot([t->sin(t),t->cos(t)^3],(0,π/2))

heatmap(randn(20,20))


text:
module OpenEXR

export load_exr, save_exr

using FileIO, Colors

module C
using OpenEXR_jll
using Colors
const ImfHalf = Float16
const ImfRgba = RGBA{ImfHalf}
include("OpenEXR_common.jl")
include("OpenEXR_api.jl")
const IMF_WRITE_RGBA = IMF_WRITE_RGB + IMF_WRITE_A
end  # module C

const MAGIC = Cint(C.IMF_MAGIC)

@enum Compression::Cint begin
    NO_COMPRESSION = C.IMF_NO_COMPRESSION
    RLE_COMPRESSION = C.IMF_RLE_COMPRESSION
    ZIPS_COMPRESSION = C.IMF_ZIPS_COMPRESSION
    ZIP_COMPRESSION = C.IMF_ZIP_COMPRESSION
    PIZ_COMPRESSION = C.IMF_PIZ_COMPRESSION
    PXR24_COMPRESSION = C.IMF_PXR24_COMPRESSION
    B44_COMPRESSION = C.IMF_B44_COMPRESSION
    B44A_COMPRESSION = C.IMF_B44A_COMPRESSION
    DWAA_COMPRESSION = C.IMF_DWAA_COMPRESSION
    DWAB_COMPRESSION = C.IMF_DWAB_COMPRESSION
end

@enum RgbaChannels::Cint begin
    WRITE_R = C.IMF_WRITE_R
    WRITE_G = C.IMF_WRITE_G
    WRITE_B = C.IMF_WRITE_B
    WRITE_A = C.IMF_WRITE_A
    WRITE_Y = C.IMF_WRITE_Y
    WRITE_C = C.IMF_WRITE_C
    WRITE_RGB = C.IMF_WRITE_RGB
    WRITE_RGBA = C.IMF_WRITE_RGBA
    WRITE_YC = C.IMF_WRITE_YC
    WRITE_YA = C.IMF_WRITE_YA
    WRITE_YCA = C.IMF_WRITE_YCA
end

function check(ret)
    ret == typeof(ret)(0) && error(unsafe_string(C.ImfErrorMessage()))
end

"""
    load_exr(filename)::(Array{RGBA{Float16},2}, RgbaChannels)

Returns the image data contained in `filename` along with flags representing the on-disk
storage format.
"""
function load_exr(filename)
    infile = C.ImfOpenInputFile(filename)  # open the file
    check(infile)
    chans = RgbaChannels(C.ImfInputChannels(infile))
    try
        # get the header
        hdr = C.ImfInputHeader(infile)

        # read its data window
        xmin = Ref{Cint}()
        ymin = Ref{Cint}()
        xmax = Ref{Cint}()
        ymax = Ref{Cint}()
        C.ImfHeaderDataWindow(hdr, xmin, ymin, xmax, ymax)

        # compute the window size
        width = xmax[] - xmin[] + 1
        height = ymax[] - ymin[] + 1

        # allocate space for the result and get its strides
        data = Array{C.ImfRgba,2}(undef, height, width)
        (xstride, ystride) = strides(data)

        # get the pointer to the data, shifting it according to the expected window
        dataptr =
            Base.unsafe_convert(Ptr{C.ImfRgba}, data) - xmin[] * xstride - ymin[] * ystride

        # copy the data
        check(C.ImfInputSetFrameBuffer(infile, dataptr, ystride, xstride))
        check(C.ImfInputReadPixels(infile, ymin[], ymax[]))

        # return the loaded raster along with the channels
        return (data, chans)
    finally
        check(C.ImfCloseInputFile(infile))
    end
end

"""
    save_exr(filename, image[, compression[, channels]])

Save `image` as an OpenEXR file in `filename`, storing the data in a format
indicated by `channels` using the `compression` algorithm.
"""
function save_exr(
    filename,
    image::AbstractArray{C.ImfRgba,2},
    compression::Compression = ZIP_COMPRESSION,
    channels::RgbaChannels = WRITE_RGBA,
)
    # get the size of the data
    (height, width) = size(image)

    # create a new header
    hdr = C.ImfNewHeader()
    check(hdr)
    try
        # set the compression
        C.ImfHeaderSetCompression(hdr, compression)

        # set the correct window sizes
        C.ImfHeaderSetDataWindow(hdr, 0, 0, width - 1, height - 1)
        C.ImfHeaderSetDisplayWindow(hdr, 0, 0, width - 1, height - 1)

        # open the output file
        outfile = C.ImfOpenOutputFile(filename, hdr, channels)
        check(outfile)
        try
            # get the strides and a pointer to the raster
            (xstride, ystride) = strides(image)
            dataptr = Base.unsafe_convert(Ptr{C.ImfRgba}, image)

            # copy the data
            check(C.ImfOutputSetFrameBuffer(outfile, dataptr, ystride, xstride))
            check(C.ImfOutputWritePixels(outfile, height))
        finally
            check(C.ImfCloseOutputFile(outfile))
        end
    finally
        C.ImfDeleteHeader(hdr)
    end
    nothing
end

function save_exr(
    filename,
    image::AbstractArray{T,2},
    compression::Compression = ZIP_COMPRESSION,
    channels::RgbaChannels = WRITE_RGBA,
) where {T}
    save_exr(filename, (c -> convert(C.ImfRgba, c)).(image), compression, channels)
end

"""
    load(filename)::Array{[RGB|RGBA|Gray|GrayA]{Float16},2}

Returns the image data contained in `filename`.
"""
function load(filename::AbstractString)
    (rgba, chans) = load_exr(filename)
    if chans == WRITE_YA
        return (c -> convert(GrayA{Float16}, c)).(rgba)
    elseif chans == WRITE_Y
        return (c -> convert(Gray{Float16}, c)).(rgba)
    elseif UInt32(chans) & UInt32(WRITE_A) == 0x00
        return (c -> convert(RGB{Float16}, c)).(rgba)
    else
        return rgba
    end
end

"""
    save(filename, image[, compression])

Save `image` as an OpenEXR file in `filename` using the `compression` algorithm.
"""
function save(
    filename::AbstractString,
    image::AbstractArray{T,2},
    compression::Compression = ZIP_COMPRESSION,
) where {T<:Transparent3}
    save_exr(filename, image, compression, WRITE_RGBA)
end

function save(
    filename::AbstractString,
    image::AbstractArray{T,2},
    compression::Compression = ZIP_COMPRESSION,
) where {T<:Color3}
    save_exr(filename, image, compression, WRITE_RGB)
end

function save(
    filename::AbstractString,
    image::AbstractArray{T,2},
    compression::Compression = ZIP_COMPRESSION,
) where {T<:TransparentGray}
    save_exr(filename, image, compression, WRITE_YA)
end

function save(
    filename::AbstractString,
    image::AbstractArray{T,2},
    compression::Compression = ZIP_COMPRESSION,
) where {T<:AbstractGray}
    save_exr(filename, image, compression, WRITE_Y)
end

# FileIO interface
load(f::File{DataFormat{:EXR}}, args...) = load(f.filename, args...)
save(f::File{DataFormat{:EXR}}, args...) = save(f.filename, args...)

end  # module OpenEXR


text:
#!/usr/bin/env julia

module VCF

using GZip
using JLD
using ProgressMeter

#Eat the second argument for non-GZip IOstreams
Base.position(b::Base.AbstractIOBuffer, ::Bool) = position(b)

@doc "Initial guess for how many nonzero entries will be read in" ->
const INITIAL_SIZE_GUESS = 1_000_000

@doc "Fudge factor used to reallocate storage vectors" ->
const SLACK_FACTOR = 1.25

@doc "Update progress bar each time this many lines are read" ->
const UPDATE_LINE_INTERVAL = 500

@doc """
Read VCF files into a [0,1,2]-genotype matrix

Implementation note:

    This is a very simplistic VCF file reader that skips a lot of the data
    stored in the file.

    All header lines are ignored.

    Only field 2 (the chromosomal position) and fields 10+ (the variants) are
    read in.

    The phases of the variants are ignored.

Input:

    filename: Name of VCF file. Can be compressed in GZip format

Output:

    Is: Row indexes of variants
    Js: Column indexes of variants
    Vs: Number of variants at position I for sample J

Side effects:

    Also displays progress meter

Reference:

    Based on the VCF Version 4.2 Specification at

    http://samtools.github.io/hts-specs/VCFv4.2.pdf
""" ->
function Base.read(filename)
    Is = Int64[]
    Js = Int64[]
    Vs = Int8[]
    sizehint!(Is, INITIAL_SIZE_GUESS)
    sizehint!(Js, INITIAL_SIZE_GUESS)
    sizehint!(Vs, INITIAL_SIZE_GUESS)

    csize = filesize(filename)
    p = Progress(csize, 2, "", 50) #Progress bar

    mode = :startline
    lineno = pos = fieldidx = skiplines = 0
    thisentry = 0
    rowid = colid = 1
    tmpbuf = IOBuffer(19)
    n = 0

    #Now actually try to read file
    if csize==0
        warn("File $filename is empty")
        @goto done
    end

    stream = if endswith(filename, ".gz")
        GZip.open(filename)
    else
        IOBuffer(Mmap.mmap(filename))
    end

    #Iterate through file
    while true
        c = read(stream, UInt8)
        n += 1
        if mode == :startline
            lineno += 1
            fieldidx = 1
            mode = c == '#' ? :skipline : :seekpos #Skip comment

        elseif mode == :skipline
            if c == '\n'
                skiplines += 1
                mode = :startline
            end

        elseif mode == :seekpos
            startline = false

            if c == '\t'
                fieldidx += 1

                if fieldidx==2
                    mode = :readpos
                    seek(tmpbuf, 0)
                end
            end

        elseif mode == :readpos
            if c == '\t'
                pos = parse(Int64, takebuf_array(tmpbuf))
                fieldidx = 3
                mode = :seekdata
            else
                write(tmpbuf, c)
            end

        elseif mode == :seekdata
            if c == '\t'
                fieldidx += 1
                if fieldidx == 10
                    mode = :readdata
                    colid = thisentry = 0
                end
            end

        elseif mode == :readdata
            if c == '0' && position(tmpbuf)==0
                continue

            elseif '0' ≤ c ≤ '9' || c == '.'
                #Horrible hack to force tmpbuf pointer forward without storing a value
                tmpbuf.ptr += 1
                #write(tmpbuf, c)

            elseif c == '\n'
                rowid += 1
                mode = :savedata

            elseif (c == '\t' || c == '\r')
                #Reset IO buffer
                if position(tmpbuf) > 0
                    thisentry += 1
                    seek(tmpbuf, 0)
                end
                colid += 1

                mode = :savedata

            elseif (c == '|' || c == '\\') #Ploidy divider
                #Reset IO buffer
                if position(tmpbuf) > 0
                    thisentry += 1
                    seek(tmpbuf, 0)
                end

            elseif c==':' #Extra genome-level data I don't understand
                mode = :skipfield

            else
                error("Unknown char $(Char(c)) (0x$(hex(c))) on line $lineno, file offset $n")
            end

        elseif mode == :skipfield
            if c == '\t'
                mode = :readdata
            elseif c == '\n'
                mode = :startline
            end
        end

        if mode == :savedata
            if lineno % UPDATE_LINE_INTERVAL == 0 #Update progress bar
                coffset = position(stream, true)
                update!(p, coffset)
            end

            if thisentry != 0#position(tmpbuf) > 0
                push!(Is, pos)
                push!(Js, colid)
                push!(Vs, thisentry) #parse(Int8, takebuf_array(tmpbuf)))
                thisentry = 0

                if length(Vs) == INITIAL_SIZE_GUESS
                    #More data than initially guessed
                    #Estimate memory consumption and resize in-memory storage
                    coffset = position(stream, true)
                    estnnz = ceil(Int, length(Vs)*csize/coffset)
                    info("Number of nonzero entries read so far:", length(Vs))
                    info("Estimated number of nonzero entries:", estnnz)
                    estnnz = ceil(Int, SLACK_FACTOR*estnnz)
                    sizehint!(Is, estnnz)
                    sizehint!(Js, estnnz)
                    sizehint!(Vs, estnnz)
                end
            end

            mode = c=='\n' ? :startline : :readdata
        end

        eof(stream) && break
    end

    @label done

    finish!(p)
    println("Actual number of header lines:", skiplines)
    println("Actual number of record lines:", lineno-skiplines)
    println("Actual number of nonzeros:", length(Vs))

    if length(Vs) > 0
        nrows = maximum(Is)
        ncols = maximum(Js)
        println("Number of record rows: ", rowid)
        println("Number of logical rows: ", nrows)
        println("Number of logical cols: ", ncols)
        println("Average density: ", length(Vs)/(Int64(nrows)*Int64(ncols)))
    end

    try
        close(stream)
    end

    Is, Js, Vs
end

function Base.parse{T<:Integer}(::Type{T}, s::AbstractArray{UInt8})
    n = zero(T)
    @fastmath for c in s
        n::T = 10n + (c - T(0x30))
        if n<0 || !('0'<=c<='9')
            error(string(Char(c),':',s,':',n))
        end
    end
    n
end

end #module VCF


text:
# A keyword argument of a function or an operator.

mutable struct KeywordClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    name::Symbol

    KeywordClause(;
             over = nothing,
             name::Union{Symbol, AbstractString}) =
        new(over, Symbol(name))
end

KeywordClause(name; over = nothing) =
    KeywordClause(over = over, name = name)

KeywordClause(name, over) =
    KeywordClause(over = over, name = name)

"""
    KW(; over = nothing, name)
    KW(name; over = nothing)
    KW(over, name)

A keyword argument of a function or an operator.

# Examples

```jldoctest
julia> c = FUN(:SUBSTRING, :zip, KW(:FROM, 1), KW(:FOR, 3));

julia> print(render(c))
SUBSTRING("zip" FROM 1 FOR 3)
```

```jldoctest
julia> c = OP(:BETWEEN, :year_of_birth, 2000, KW(:AND, 2010));

julia> print(render(c))
("year_of_birth" BETWEEN 2000 AND 2010)
```
"""
KW(args...; kws...) =
    KeywordClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(KW), pats::Vector{Any}) =
    dissect(scr, KeywordClause, pats)

function PrettyPrinting.quoteof(c::KeywordClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(KW), quoteof(c.name))
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::KeywordClause, c′) =
    KeywordClause(over = rebase(c.over, c′), name = c.name)


text:
using FractionalDiffEq, Plots, SpecialFunctions

α = 0.5

function testfun(t, y)
    return 40320/gamma(9-α)*t^(8-α) - 3*gamma(5+α/2)/gamma(5-α/2)*t^(4-α/2) + 9/4*gamma(α+1) + (3/2*t^α/2-t^4)^3-y^1.5
end

prob = FODEProblem(testfun, α, 0, 1, 0.001)

result = solve(prob, PECE())

tfun(t) = t^8-3t^4.25+9/4*t^0.5

tspan = collect(0:0.001:1)

plot(tspan, result)

plot!(tspan, tfun, lw=2)


text:
include("../aimajulia.jl");

using Test;

using Main.aimajulia;

#The following CSP tests are from the aima-python doctest

@test (solution(depth_first_graph_search(aimajulia.australia_csp))...,) == (("NSW","B"),("Q","G"),("NT","B"),("T","B"),("V","G"),("SA","R"),("WA","G"));

d = CSPDict(ConstantFunctionDict(42));

@test d["life"] == 42;

@test (!(typeof(backtracking_search(aimajulia.australia_csp)) <: Nothing) == true);

@test (!(typeof(backtracking_search(aimajulia.australia_csp,
                                    select_unassigned_variable=minimum_remaining_values)) <: Nothing) == true);

@test (!(typeof(backtracking_search(aimajulia.australia_csp,
                                    order_domain_values=least_constraining_values)) <: Nothing) == true);

@test (!(typeof(backtracking_search(aimajulia.australia_csp,
                                    select_unassigned_variable=minimum_remaining_values,
                                    order_domain_values=least_constraining_values)) <: Nothing) == true);

@test (!(typeof(backtracking_search(aimajulia.australia_csp, inference=forward_checking)) <: Nothing) == true);

@test (!(typeof(backtracking_search(aimajulia.australia_csp,
                                    inference=maintain_arc_consistency)) <: Nothing) == true);

@test (!(typeof(backtracking_search(aimajulia.australia_csp,
                                    select_unassigned_variable=minimum_remaining_values,
                                    order_domain_values=least_constraining_values,
                                    inference=maintain_arc_consistency)) <: Nothing) == true);

topological_sorted_nodes, parent_dict = topological_sort(aimajulia.australia_csp, "NT");

@test topological_sorted_nodes == Any["NT","SA","Q","NSW","V","WA"];

@test haskey(parent_dict, "NT") == false;

@test parent_dict == Dict{Any,Any}(Pair("NSW","Q"), Pair("Q","SA"), Pair("V","NSW"), Pair("SA","NT"), Pair("WA","SA"));

@test length(backtracking_search(NQueensCSP(8))) == 8;

e = SudokuCSP(aimajulia.easy_sudoku_grid);

@test display(e, infer_assignment(e)) == ". . 3 | . 2 . | 6 . .\n9 . . | 3 . 5 | . . 1\n. . 1 | 8 . 6 | 4 . .\n------+-------+------\n. . 8 | 1 . 2 | 9 . .\n7 . . | . . . | . . 8\n. . 6 | 7 . 8 | 2 . .\n------+-------+------\n. . 2 | 6 . 9 | 5 . .\n8 . . | 2 . 3 | . . 9\n. . 5 | . 1 . | 3 . .";

AC3(e);

@test display(e, infer_assignment(e)) == "4 8 3 | 9 2 1 | 6 5 7\n9 6 7 | 3 4 5 | 8 2 1\n2 5 1 | 8 7 6 | 4 9 3\n------+-------+------\n5 4 8 | 1 3 2 | 9 7 6\n7 2 9 | 5 6 4 | 1 3 8\n1 3 6 | 7 9 8 | 2 4 5\n------+-------+------\n3 7 2 | 6 8 9 | 5 1 4\n8 1 4 | 2 5 3 | 7 6 9\n6 9 5 | 4 1 7 | 3 8 2";

@test !(typeof(backtracking_search(SudokuCSP(aimajulia.harder_sudoku_grid), select_unassigned_variable=minimum_remaining_values, inference=forward_checking)) <: Nothing)

@test solve_zebra(ZebraCSP(), backtracking_search) == (5,
                                                        1,
                                                        75472,
                                                        Dict{Any,Any}(Pair{Any,Any}("Tea",2),
                                                                    Pair{Any,Any}("Red",3),
                                                                    Pair{Any,Any}("Kools",1),
                                                                    Pair{Any,Any}("Green",5),
                                                                    Pair{Any,Any}("Horse",2),
                                                                    Pair{Any,Any}("Zebra",5),
                                                                    Pair{Any,Any}("OJ",4),
                                                                    Pair{Any,Any}("Milk",3),
                                                                    Pair{Any,Any}("Coffee",5),
                                                                    Pair{Any,Any}("Ukranian",2),
                                                                    Pair{Any,Any}("Japanese",5),
                                                                    Pair{Any,Any}("Snails",3),
                                                                    Pair{Any,Any}("Spaniard",4),
                                                                    Pair{Any,Any}("Water",1),
                                                                    Pair{Any,Any}("Winston",3),
                                                                    Pair{Any,Any}("Norwegian",1),
                                                                    Pair{Any,Any}("Fox",1),
                                                                    Pair{Any,Any}("Dog",4),
                                                                    Pair{Any,Any}("Ivory",4),
                                                                    Pair{Any,Any}("Englishman",3),
                                                                    Pair{Any,Any}("Yellow",1),
                                                                    Pair{Any,Any}("LuckyStrike",4),
                                                                    Pair{Any,Any}("Parliaments",5),
                                                                    Pair{Any,Any}("Blue",2),
                                                                    Pair{Any,Any}("Chesterfields",2)));

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork


text:
function SGP4(arg0::TLE, arg1::AttitudeProvider, arg2::jdouble)
    return SGP4((TLE, AttitudeProvider, jdouble), arg0, arg1, arg2)
end

function SGP4(arg0::TLE, arg1::AttitudeProvider, arg2::jdouble, arg3::Frame)
    return SGP4((TLE, AttitudeProvider, jdouble, Frame), arg0, arg1, arg2, arg3)
end


text:
# Plot from MPB-data (::AbstractFourierLattice) -------------------------------------------
using .PyPlot

function plot_lattice_from_mpbparams(filepath::String; kwargs...)
    Rs, flat, isoval, _ = lattice_from_mpbparams(filepath)
    plot(flat, Rs; isoval=isoval, kwargs...) # via Crystalline.jl overload of PyPlot.jl's `plot`
end


text:
include("gmm.jl")
import Random, Dates
using RCall
Random.seed!(0);

ShowTime() = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")

# Generate Data
N = 10000
m = [3., -1., 2.]
s = [.1, .05, .15]
w = [.5, .3, .2]
K = length(m)
y = [begin
       k = rand(Categorical(w))
       rand(Normal(m[k], s[k]))
     end for i in 1:N]

metrics = Dict(m => Float64[] for m in (:elbo, :ll, :lp, :lq))

# TRAIN_VERSION = 1
TRAIN_VERSION = 2
if TRAIN_VERSION == 1
  Random.seed!(6); # VERSION I
else
  Random.seed!(0); # VERSION II
  # Random.seed!(100); # VERSION II
end
mp = State(K)
loss(y) = -compute_elbo(mp, y, N, metrics) / N
ps = ADVI.vparams(mp)

opt = ADAM(1e-1)
# opt = RMSProp(1e-1)
minibatch_size = 500
niters = 1000

if TRAIN_VERSION == 1
  # VERSION I
  @time Flux.train!(loss, ps, [(y[sample(1:N, minibatch_size, replace=false)], )
                               for t in 1:niters], opt)
else
  @time for t in 1:niters
    idx = Distributions.sample(1:N, minibatch_size, replace=false)

    # Flux.train!(loss, ps, [(y[idx], )], opt)
    # OR
    gs = Tracker.gradient(() -> loss(y[idx]), ps)
    Flux.Tracker.update!(opt, ps, gs)
    #= looking at each grad
    for k in keys(gs.grads)
      println(gs.grads[k])
    end
    =#

    if t % 100 == 0
      println("$(ShowTime()) | $(t)/$(niters) ")
    end
  end
end

println(metrics[:elbo][1:100:end] / N)
R"plot"(metrics[:elbo]/N, typ="l", xlab="", ylab="")

samps = [rsample(mp)[2] for i in 1:100]
m_post = cat([s.m for s in samps]..., dims=2)
s_post = cat([s.s for s in samps]..., dims=2)
w_post = vcat([s.w for s in samps]...)

println("m mean: $(mean(m_post, dims=2)) | true: $m")
println("s mean: $(mean(s_post, dims=2)) | true: $s")
println("w mean: $(mean(w_post, dims=1)) | true: $w")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
include("../src/Kiltro.jl")
using .Kiltro

#======================== MESH ===========================#
# Mesh input to the problem, this case a hexaedra of l=1.0.
points = [0.0 0.0 0.0; 1.0 0.0 0.0; 1.0 1.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.5; 1.0 0.0 0.5; 1.0 1.0 0.5; 0.0 1.0 0.5; 0.0 0.0 1.0; 1.0 0.0 1.0; 1.0 1.0 1.0; 0.0 1.0 1.0]
elements = [1 2 3 4 5 6 7 8; 5 6 7 8 9 10 11 12]
faces = [1 2 3 4; 1 5 6 2; 2 6 7 3; 3 7 8 4; 1 4 8 5; 5 9 10 6; 6 10 11 7; 7 11 12 8; 5 8 12 9; 9 12 11 10]

nnode = size(points,1)
ndim = size(points,2)

DirichletBoundary1 = falses(nnode)
DirichletBoundary2 = falses(nnode)
DirichletBoundary3 = falses(nnode)
DirichletBoundary4 = falses(nnode)
DirichletBoundary5 = falses(nnode)
DirichletBoundary1[faces[1,:]] .= true
DirichletBoundary2[faces[2,:]] .= true
DirichletBoundary3[faces[5,:]] .= true
DirichletBoundary4[faces[6,:]] .= true
DirichletBoundary5[faces[9,:]] .= true
NeumannBoundary = falses(nnode)
NeumannBoundary[faces[10,:]] .= true

#====================== MATERIAL =========================#
#material parameters for neo-Hooke hyperelastic model
mu = 15.0e-3
kappa = 100.0*mu

#===================== FORMULATION =======================#
#sets the function space for the displacement problem
nvar = 3
x, w = QuadratureRule(2,-1.0,1.0)
AllGauss, Basis, Jm = FunctionSpace(x,w)

#================= BOUNDARY CONDITIONS ===================#
#sets the boundary conditions, with symmetries (dirichlet) and the external force (neumann)
dirichlet_flags = zeros(nnode,3).+NaN
dirichlet_flags[DirichletBoundary1,3] .= 0.0
dirichlet_flags[DirichletBoundary2,2] .= 0.0
dirichlet_flags[DirichletBoundary3,1] .= 0.0
dirichlet_flags[DirichletBoundary4,2] .= 0.0
dirichlet_flags[DirichletBoundary5,1] .= 0.0

neumann_flags = zeros(nnode,3).+NaN
neumann_flags[NeumannBoundary,3] .= 0.005

#======================= SOLVER ==========================#
TotalDisp = Solve(AllGauss,Basis,Jm,points,elements,mu,kappa,dirichlet_flags,neumann_flags,ndim,nvar)

#====================== SOLUTION =========================#

println(TotalDisp)

max_stars_repo_name: jdlaubrie/Kiltro

max_issues_repo_name: jdlaubrie/Kiltro

max_forks_repo_name: jdlaubrie/Kiltro


text:
cd("D:\\OneDrive\\Current\\meanVar\\Code\\hox_krox")
include("bp_sharpening_experiment_function.jl")
cd("D:\\OneDrive\\Backups\\ArchivedProjects\\UCI\\MeanVar Independance\\arxdata\\hox_krox_data")
using JLD

###################### Analysis Consts

const x = X[1,31:61]-50
const y = Y[:,1]
const boundary_xs = Vector{Float64}(y)

function calculate_boundary_stds(u)
  boundary_stds = Vector{Float64}(length(u))
  for ti in eachindex(u)
    gₕ = u[ti][:,31:61,7]
    for i in 1:size(gₕ,1)
      cur_g = gₕ[i,:] .< 0.5
      idx = findfirst(cur_g)
      boundary_xs[i] = x[idx]
    end
    boundary_stds[ti] = sqrt(var(boundary_xs))
  end
  boundary_stds
end

function calculate_boundary_stds2(u)
  boundary_stds = Vector{Float64}(length(u))
  for ti in eachindex(u)
    bad_distances = Vector{Float64}(0)
    gₕ = u[ti][:,31:61,7]
    numy = size(gₕ,1)
    numx = size(gₕ,2)
    for i in 1:numy
      cur_g = gₕ[i,:] .< 0.5
      idx = findfirst(cur_g)
      idx2 = findlast(.!(cur_g))
      mid_idx = (idx+idx2)/2
      midpoint = (x[idx] + x[idx2])/2
      bad_left_idxs  = (1:round(Int,mid_idx,RoundDown))[cur_g[1:round(Int,mid_idx,RoundDown)]]
      bad_right_idxs = (round(Int,mid_idx,RoundUp):numx)[.!cur_g[round(Int,mid_idx,RoundUp):numx]]
      bad_left_distances = x[bad_left_idxs] - midpoint
      bad_right_distances = x[bad_right_idxs] - midpoint
      append!(bad_distances,[bad_left_distances;bad_right_distances])
    end
    boundary_stds[ti] = 2std(bad_distances)/numy
  end
  boundary_stds
end

function calculate_boundary_stds3(u,reduction=mean)
  boundary_stds = Vector{Float64}(length(u))
  for ti in eachindex(u)
    gₕ = u[ti][:,31:61,7]
    for i in 1:size(gₕ,1)
      cur_g = gₕ[i,:] .< 0.5
      idx = findfirst(cur_g)
      idx2 = findlast(.!(cur_g))
      boundary_xs[i] = max(x[idx2] - x[idx],0)
    end
    boundary_stds[ti] = reduction(boundary_xs)
  end
  boundary_stds
end

function calculate_boundary_stds4(u,cell_length=10)
  boundary_stds = Vector{Float64}(length(u))
  for ti in eachindex(u)
    gₕ = u[ti][:,31:61,7]
    for i in 1:size(gₕ,1)
      cur_g = gₕ[i,:] .< 0.5
      idx = findfirst(cur_g)
      idx2 = findlast(.!(cur_g))
      boundary_xs[i] = max(x[idx2] - x[idx],0)
    end
    boundary_stds[ti] = sum(boundary_xs .> cell_length)
  end
  boundary_stds
end

effective_σRA(u) = var(u[end][:,46,4])

####################### Load the data

using JLD
data1 = load("wildtype/bp_100.0_10.jld")
data2 = load("bp_knocks/bp_5.0_7.jld")
data3 = load("cyp_knocks/cyp_33000.0_1.jld")
t = data1["t"]
u1 = data1["u"]
t1 = data1["t"]
u2 = data2["u"]
t2 = data2["t"]
u3 = data3["u"]
t3 = data3["t"]

####################### Generate Plots

## Boundary Sharpening Diagrams

gₕ = u1[2][:,31:61,7]
gₖ = u1[2][:,31:61,8]

PlotUtils.default_cgrad(default = :sequential, sequential = :redsblues)

p11 = heatmap(x,y,gₕ,cbar =  false,
              ylabel="Wildtype",clims=(0.0,2.0))

              p11 = heatmap(x,y,gₕ,cbar =  true,
                            ylabel="Wildtype",clims=(0.0,2.0))

gₕ = u1[19][:,31:61,7] # 1/2 hour
gₖ = u1[19][:,31:61,8]

p12 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u1[37][:,31:61,7] # 1 hour
gₖ = u1[37][:,31:61,8]

p13 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u1[2*37-1][:,31:61,7] # 2 hours
gₖ = u1[2*37-1][:,31:61,8]

p14 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

p1 = plot(p11,p12,p13,p14,layout=grid(1,4))

gₕ = u2[2][:,31:61,7]
gₖ = u2[2][:,31:61,8]

p21 = heatmap(x,y,gₕ,cbar =  false, ylabel="BP Knockdown",clims=(0.0,2.0))

gₕ = u2[19][:,31:61,7] # 1/2 hour
gₖ = u2[19][:,31:61,8]

p22 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u2[37][:,31:61,7] # 1 hour
gₖ = u2[37][:,31:61,8]

p23 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u2[2*37-1][:,31:61,7] # 2 hours
gₖ = u2[2*37-1][:,31:61,8]

p24 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u3[2][:,31:61,7]
gₖ = u3[2][:,31:61,8]

p31 = heatmap(x,y,gₕ,cbar =  false, ylabel="Cyp Knockdown",clims=(0.0,2.0))

gₕ = u3[19][:,31:61,7] # 1/2 hour
gₖ = u3[19][:,31:61,8]

p32 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u3[37][:,31:61,7] # 1 hour
gₖ = u3[37][:,31:61,8]

p33 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

gₕ = u3[2*37-1][:,31:61,7] # 2 hours
gₖ = u3[2*37-1][:,31:61,8]

p34 = heatmap(x,y,gₕ,cbar =  false,clims=(0.0,2.0))

p = plot(p11,p12,p13,p14,p21,p22,p23,p24,p31,p32,p33,p34,layout=grid(3,4))

####################### Calculate boundary variances

#### Frequency

boundary_stds_bp100 = Matrix{Float64}(length(t),10)
for i in 1:10
  data = load("wildtype/bp_100.0_$(i).jld"); u = data["u"]
  t = data["t"]
  boundary_stds = calculate_boundary_stds4(u)
  boundary_stds_bp100[:,i] = boundary_stds
end
p2 = plot(t[1:73]/60/60 .+ 10,mean(boundary_stds_bp100,2)[1:73],lw=3,
        label="Wildtype",xlabel="Time (hpf)",
        ylabel="# Displaced\nCells",
        yaxis = [0,Inf])
boundary_stds_bp5 = Matrix{Float64}(length(t),10)
for i in 1:10
  data = load("bp_knocks/bp_5.0_$(i).jld"); u = data["u"]
  t = data["t"]
  boundary_stds = calculate_boundary_stds4(u)
  boundary_stds_bp5[:,i] = boundary_stds
end
plot!(p2,t[1:73]/60/60 .+ 10,mean(boundary_stds_bp5,2)[1:73],
      lw=3,label="BP Knockdown")

boundary_stds_cyp = Matrix{Float64}(length(t),10)
for i in 1:10
  data = load("cyp_knocks/cyp_33000.0_$(i).jld"); u = data["u"]
  t = data["t"]
  boundary_stds = calculate_boundary_stds4(u)
  boundary_stds_cyp[:,i] = boundary_stds
end
plot!(p2,t[1:73]/60/60 .+ 10,mean(boundary_stds_cyp,2)[1:73],
      lw=3,label="Cyp Knockdown")

### Maximum Length

boundary_stds_bp100 = Matrix{Float64}(length(t),10)
for i in 1:10
  data = load("wildtype/bp_100.0_$(i).jld"); u = data["u"]
  t = data["t"]
  boundary_stds = calculate_boundary_stds3(u,maximum)
  boundary_stds_bp100[:,i] = boundary_stds
end

p3 = plot(t[1:73]/60/60 .+ 10,mean(boundary_stds_bp100,2)[1:73]./10,
          lw=3,legend=false,
          xlabel="Time (hpf)", ylabel="Max Boundary\nLength\n(cell lengths)",
          yaxis = [0,12])

boundary_stds_bp5 = Matrix{Float64}(length(t),10)
for i in 1:10
  data = load("bp_knocks/bp_5.0_$(i).jld"); u = data["u"]
  t = data["t"]
  boundary_stds = calculate_boundary_stds3(u,maximum)
  boundary_stds_bp5[:,i] = boundary_stds
end
plot!(p3,t[1:73]/60/60 .+ 10,mean(boundary_stds_bp5,2)[1:73]./10,lw=3)

boundary_stds_cyp = Matrix{Float64}(length(t),10)
for i in 1:10
  data = load("cyp_knocks/cyp_33000.0_$(i).jld"); u = data["u"]
  t = data["t"]
  boundary_stds = calculate_boundary_stds3(u,maximum)
  boundary_stds_cyp[:,i] = boundary_stds
end
plot!(p3,t[1:73]/60/60 .+ 10,mean(boundary_stds_cyp,2)[1:73]./10,lw=3)

l = @layout [  a{0.75h}
             [b c]]

plot(p,p2,p3,layout=l,size=(750,500))

################################################################################

## Noise Success

Vbps = 5.0:10.0:105.0
ϵₕₖs = 0.1:0.025:0.3
num_sims = length(Vbps)*length(ϵₕₖs)
σRAs = Vector{Float64}(num_sims)
σGenes = Vector{Float64}(num_sims)
boundary_sharpness = Vector{Float64}(num_sims)
mean_boundary_sharpness = Vector{Float64}(num_sims)
number_boundary_over = Vector{Float64}(num_sims)
iter = 0
t[73]./60./60 == 2 #hpf
for Vbp in Vbps, ϵₕₖ in ϵₕₖs
  data = load("noise_noise_plot/bp_$(Int(Vbp))_epshk_$(ϵₕₖ).jld")
  u = data["u"]
  iter += 1
  σRAs[iter] = effective_σRA(u)
  boundary_sharpness[iter] = calculate_boundary_stds3(u,maximum)[73]
  mean_boundary_sharpness[iter] = calculate_boundary_stds3(u,mean)[73]
  number_boundary_over[iter] = calculate_boundary_stds4(u)[73]
  σGenes[iter] = ϵₕₖ
end

sharpened = boundary_sharpness.< 30

p1 = scatter(σGenes[sharpened],σRAs[sharpened],
        yscale=:log10,markersize=10,markercolor=:green,
        title="Maximum Displacement < 3 Cell Lengths",
        xlabel="Gene Regulatory Noise",
        ylabel="Effective RA Signal Noise",
        labels="Successful Shapening Event")
scatter!(p1,σGenes[!sharpened],σRAs[!sharpened],yscale=:log10,markersize=10,
        yscale=:log10,markersize=10,markercolor=:red,
        labels="Unsuccessful Shapening Event")

mean_sharpened = mean_boundary_sharpness .< 5

p2 = scatter(σGenes[mean_sharpened],σRAs[mean_sharpened],
        yscale=:log10,markersize=10,markercolor=:green,
        title="Mean Displacement < 1/2 Cell Lengths",
        xlabel="Gene Regulatory Noise",
        ylabel="Effective RA Signal Noise",
        labels="Successful Shapening Event")
scatter!(p2,σGenes[!mean_sharpened],σRAs[!mean_sharpened],
        yscale=:log10,markersize=10,markercolor=:red,
        labels="Unsuccessful Shapening Event")

number_sharpened = number_boundary_over .<= 3

plot(p1,p2)

p5 = scatter(σGenes[number_sharpened],σRAs[number_sharpened],
        yscale=:log10,markersize=10,markercolor=:green,
        xlabel="Gene Regulatory Noise",
        ylabel="Effective RA Signal Noise",
        labels="Successful Shapening Event")
scatter!(σGenes[!number_sharpened],σRAs[!number_sharpened],
         yscale=:log10,markersize=10,markercolor=:red,
         labels="Unsuccessful Shapening Event")

################################################################################

## Cyp Shift

RAR_wt = u1[end][:,31:61,4]
RAR_cyp = u3[end][:,31:61,4]
plot(X[1,31:61],mean(RAR_wt,1)',lw=3,labels="Wildtype",xlabel="Anterior-Posterior (um)",ylabel="Average [RAR] (uM)")
plot!(X[1,31:61],mean(RAR_cyp,1)',lw=3,labels="Cyp Knockdown")

################################################################################

## Supplement: Sharpening Index

boundary_stds_bp100 = Matrix{Float64}(length(t),10)
for i in 1:10
 data = load("wildtype/bp_100.0_$(i).jld"); u = data["u"]
 t = data["t"]
 boundary_stds = calculate_boundary_stds2(u)
 boundary_stds_bp100[:,i] = boundary_stds
end
plot(t,mean(boundary_stds_bp100,2),lw=3,label="Wildtype",xlabel="Time", ylabel="Sharpness Index")
boundary_stds_bp5 = Matrix{Float64}(length(t),10)
for i in 1:10
 data = load("bp_knocks/bp_5.0_$(i).jld"); u = data["u"]
 t = data["t"]
 boundary_stds = calculate_boundary_stds2(u)
 boundary_stds_bp5[:,i] = boundary_stds
end
plot!(t,mean(boundary_stds_bp5,2),lw=3,label="BP Knockdown")
boundary_stds_cyp = Matrix{Float64}(length(t),10)
for i in 1:10
 data = load("cyp_knocks/cyp_33000.0_$(i).jld"); u = data["u"]
 t = data["t"]
 boundary_stds = calculate_boundary_stds2(u)
 boundary_stds_cyp[:,i] = boundary_stds
end
plot!(t,mean(boundary_stds_cyp,2),lw=3,label="Cyp Knockdown")

####################### Gradient Plots

p1 = surface(X[:,31:61],Y[:,31:61],u1[end][:,31:61,2],title="[RAin]")
p2 = surface(X[:,31:61],Y[:,31:61],u1[end][:,31:61,4],zlims=[0.0,2.1],clims=(0.0,2.1),title="RAR")
p3 = surface(X[:,31:61],Y[:,31:61],u2[end][:,31:61,2],title="[RAin], BP Knockdown")
p4 = surface(X[:,31:61],Y[:,31:61],u2[end][:,31:61,4],zlims=[0.0,2.1],title="[RAR], BP Knockdown")

plot(p1,p2,p3,p4,cbar=false)

###################### Gradient Heatmaps

p1 = heatmap(x,y,u1[end][:,31:61,2],title="RAin",cbar=false)
p2 = heatmap(x,y,u1[end][:,31:61,4],zlims=[0.0,2.1],cbar=false,clims=(0.0,2.1),title="RAR")

p3 = heatmap(x,y,u2[end][:,31:61,2],cbar=false,title="RAin")
p4 = heatmap(x,y,u2[end][:,31:61,4],cbar=false,zlims=[0.0,2.1],title="RAR")

plot(p1,p2,p3,p4,cbar=false)

max_stars_repo_name: ChrisRackauckas/RAMeanVar

max_issues_repo_name: ChrisRackauckas/RAMeanVar

max_forks_repo_name: ChrisRackauckas/RAMeanVar


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

include("model_Resource.jl")
include("model_RouteFilterListResult.jl")
include("model_RouteFilterPropertiesFormat.jl")
include("model_RouteFilterRuleListResult.jl")
include("model_RouteFilterRulePropertiesFormat.jl")
include("model_SubResource.jl")
include("model_ExpressRouteCircuitPeering.jl")
include("model_PatchRouteFilter.jl")
include("model_PatchRouteFilterRule.jl")
include("model_RouteFilter.jl")
include("model_RouteFilterRule.jl")


text:
"""
Solution data structure definition
"""
mutable struct Solution

    arc_selection::Dict{Tuple{Int, Int}, Float64} # selection rate of the arcs
    inst::Instance # corresponding instance
    score::Float64 # objective value

    """
    Solution object builder
    """
    function Solution(inst::Instance)

        this = new()

        # initialization of the containers which will be filled dynamically
        this.arc_selection = Dict{Tuple{Int, Int}, Float64}()

        this.inst = inst
        return this

    end

end

"""
Prints the solution
"""
function pprint(sol::Solution)

    ln0()
    ln0("Score: $(sol.score)")
    ln0()
    ln0("Selected arcs:")
    for a in sol.inst.feasible_arcs
        if sol.arc_selection[a] >= 1e-8
            ln0("$(a[1]) $(a[2]) => $(sol.arc_selection[a])")
        end
    end

end

"""
Exports the solutions as a .sol file
"""
function export_to_file(sol::Solution)

    open("solutions/$(sol.inst.name)_cost=$(round(Int, sol.score)).sol", "w") do out
        redirect_stdout(out) do
            println(round(Int, sol.score))
            println()
            for a in sol.inst.feasible_arcs
                if sol.arc_selection[a] >= 1e-8
                    println("$(a[1]) $(a[2])")
                end
            end
        end
    end

end

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
using BlockArrays

struct PartiallyImplementedBlockVector <: AbstractBlockArray{Float64,1} end

@testset "partially implemented block array" begin
    # the error thrown before was incorrect
    A = PartiallyImplementedBlockVector()
    try
        getblock(A, 1)
    catch err
        @test err isa ErrorException && err.msg == "getblock for PartiallyImplementedBlockVector is not implemented"
    end
    try
        getblock!(zeros(5), A, Block(1))
    catch err
        @test err isa ErrorException && err.msg == "getblock! for PartiallyImplementedBlockVector is not implemented"
    end
    try
        BlockArrays.setblock!(A, zeros(5), Block(1))
    catch err
        @test err isa ErrorException && err.msg == "setblock! for PartiallyImplementedBlockVector is not implemented"
    end
    try
        BlockArrays.blocksize(A, 2)
    catch err
        @test err isa ErrorException && err.msg == "blocksizes for PartiallyImplementedBlockVector is not implemented"
    end
end

@testset "Array block interface" begin
    @test blocksizes(1) == BlockArrays.BlockSizes{0}()
    @test 1[Block()] == 1

    A = randn(5)
    @test blocksizes(A) == BlockArrays.BlockSizes([5])
    @test A[Block(1)] == A
    view(A,Block(1))[1] = 2
    @test A[1] == 2
    @test_throws BoundsError A[Block(2)]

    A = randn(5,5)
    @test A[Block(1,1)] == A
end

@testset "Triangular/Symmetric/Hermitian block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (1:4))
    A .= randn.() .+ randn.().*im

    @test UpperTriangular(A)[Block(2,2)] == UpperTriangular(A[2:3,2:3])
    @test UpperTriangular(A)[Block(2,3)] == A[2:3,4:6]
    @test UpperTriangular(A)[Block(3,2)] == zeros(3,2)
    @test Symmetric(A)[Block(2,2)] == Symmetric(A[2:3,2:3])
    @test Symmetric(A)[Block(2,3)] == A[2:3,4:6]
    @test Symmetric(A)[Block(3,2)] == transpose(A[2:3,4:6])
    @test Hermitian(A)[Block(2,2)] == Hermitian(A[2:3,2:3])
    @test Hermitian(A)[Block(2,3)] == A[2:3,4:6]
    @test Hermitian(A)[Block(3,2)] == A[2:3,4:6]'
end

@testset "Adjoint/Transpose block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (2:5))
    A .= randn.() .+ randn.().*im

    @test blocksizes(A') == BlockArrays.BlockSizes(2:5, 1:4)
    @test blocksizes(Transpose(A)) == BlockArrays.BlockSizes(2:5, 1:4)

    @test A'[Block(2,2)] == A[Block(2,2)]' == A[2:3,3:5]'
    @test transpose(A)[Block(2,2)] == transpose(A[2:3,3:5])
    @test A'[Block(2,3)] == A[Block(3,2)]'
    @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

    @test BlockArray(A') == A'
    @test BlockArray(transpose(A)) == transpose(A)
end


text:
construct_device!(canonical::Canonical, sys::PSY.System,
                  model::DeviceModel{B, Br},
                  ::Type{CopperPlatePowerModel};
                  kwargs...) where {B<:PSY.DCBranch,
                                    Br<:AbstractBranchFormulation} = nothing

construct_device!(canonical::Canonical, sys::PSY.System,
                  model::DeviceModel{B, Br},
                  ::Type{CopperPlatePowerModel};
                  kwargs...) where {B<:PSY.ACBranch,
                                    Br<:AbstractBranchFormulation} = nothing

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{B, Br},
                           ::Type{S};
                           kwargs...) where {B<:PSY.Branch,
                                             Br<:AbstractBranchFormulation,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(B, sys)

    if validate_available_devices(devices,B)
        return
    end

    branch_rate_bounds!(canonical, devices, Br, S)

    branch_rate_constraint!(canonical, devices, Br, S)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{PSY.MonitoredLine, FlowMonitoredLine},
                           ::Type{S};
                           kwargs...) where {S<:PM.AbstractPowerModel}

    devices = PSY.get_components(PSY.MonitoredLine, sys)

    if validate_available_devices(devices, PSY.MonitoredLine)
        return
    end

    branch_rate_bounds!(canonical,
                        devices,
                        model.formulation,
                        S)

    branch_rate_constraint!(canonical,
                        devices,
                        model.formulation,
                        S)

    branch_flow_constraint!(canonical,
                        devices,
                        model.formulation,
                        S)

    return

end

 construct_device!(canonical::Canonical, sys::PSY.System,
                   model::DeviceModel{B, Br},
                   ::Type{S};
                   kwargs...) where {B<:PSY.Branch,
                                     Br<:Union{Type{StaticLineUnbounded},
                                               Type{StaticTransformerUnbounded}},
                                     S<:PM.AbstractPowerModel} = nothing

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{B, Br},
                           ::Type{S};
                           kwargs...) where {Br<:AbstractBranchFormulation,
                                             B<:PSY.DCBranch,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(B, sys)

    if validate_available_devices(devices, B)
        return
    end

    branch_rate_constraint!(canonical, devices, Br, S)

    return

end