text:
module BasicModelInterface

"""
    initialize(model, [config_file])

Perform startup tasks for the model.

Perform all tasks that take place before entering the model's time
loop, including opening files and initializing the model state. Model
inputs are read from a text-based configuration file, specified by
`config_file`.

# Arguments
config_file : str, optional
    The path to the model configuration file.

# Notes
Models should be refactored, if necessary, to use a
configuration file. CSDMS does not impose any constraint on
how configuration files are formatted, although YAML is
recommended. A template of a model's configuration file
with placeholder values is used by the BMI.
"""
function initialize end

"""
    update(model)

Advance model state by one time step.

Perform all tasks that take place within one pass through the model's
time loop. This typically includes incrementing all of the model's
state variables. If the model's state variables don't change in time,
then they can be computed by the :func:`initialize` method and this
method can return with no action.
"""
function update end

"""
    update_until(model, time)

Advance model state until the given time.

The given `time` must be a model time later than the current model time.
"""
function update_until end

"""
    finalize(model)

Perform tear-down tasks for the model.

Perform all tasks that take place after exiting the model's time
loop. This typically includes deallocating memory, closing files and
printing reports.
"""
function finalize end

"""
    get_component_name(model)

Name of the component.
"""
function get_component_name end

"""
    get_input_item_count(model)

Count of a model's input variables.
"""
function get_input_item_count end

"""
    get_output_item_count(model)

Count of a model's output variables.
"""
function get_output_item_count end

"""
    get_input_var_names(model)

List of a model's input variables.

Input variable names must be CSDMS Standard Names, also known
as *long variable names*.

# Notes
Standard Names enable the CSDMS framework to determine whether
an input variable in one model is equivalent to, or compatible
with, an output variable in another model. This allows the
framework to automatically connect components.

Standard Names do not have to be used within the model.
"""
function get_input_var_names end

"""
    get_output_var_names(model)

List of a model's output variables.

Output variable names must be CSDMS Standard Names, also known
as *long variable names*.
"""
function get_output_var_names end

"""
    get_var_grid(model, name)

Get grid identifier integer for the given variable.

The `name` can be an input or output variable name, a CSDMS Standard Name.
"""
function get_var_grid end

"""
    get_var_type(model, name)

Get data type of the given variable.

The `name` can be an input or output variable name, a CSDMS Standard Name.
"""
function get_var_type end

"""
    get_var_units(model, name)

Get units of the given variable.

Standard unit names, in lower case, should be used, such as
``meters`` or ``seconds``. Standard abbreviations, like ``m`` for
meters, are also supported. For variables with compound units,
each unit name is separated by a single space, with exponents
other than 1 placed immediately after the name, as in ``m s-1``
for velocity, ``W m-2`` for an energy flux, or ``km2`` for an
area.

The `name` can be an input or output variable name, a CSDMS Standard Name.

# Notes
CSDMS uses the [UDUNITS](http://www.unidata.ucar.edu/software/udunits)
standard from Unidata.
"""
function get_var_units end

"""
    get_var_itemsize(model, name)

Get memory use for each array element in bytes.

The `name` can be an input or output variable name, a CSDMS Standard Name.
"""
function get_var_itemsize end

"""
    get_var_nbytes(model, name)

Get size, in bytes, of the given variable.

The `name` can be an input or output variable name, a CSDMS Standard Name.
"""
function get_var_nbytes end

"""
    get_var_location(model, name)

Get the grid element type that the a given variable is defined on.

The grid topology can be composed of *nodes*, *edges*, and *faces*.

- node: A point that has a coordinate pair or triplet: the most
    basic element of the topology.
- edge: A line or curve bounded by two nodes.
- face: A plane or surface enclosed by a set of edges. In a 2D
    horizontal application one may consider the word "polygon",
    but in the hierarchy of elements the word "face" is most common.

The `name` can be an input or output variable name, a CSDMS Standard Name.

# Returns
The grid location on which the variable is defined. Must be one of
`"node"`, `"edge"`, or `"face"`.

# Notes
CSDMS uses the [ugrid conventions](http://ugrid-conventions.github.io/ugrid-conventions)
to define unstructured grids.
"""
function get_var_location end

"""
    get_current_time(model)

Current time of the model.
"""
function get_current_time end

"""
    get_start_time(model)

Start time of the model.
"""
function get_start_time end

"""
    get_end_time(model)

End time of the model.
"""
function get_end_time end

"""
    get_time_units(model)

Time units of the model; e.g., `days` or `s`.

# Notes
CSDMS uses the [UDUNITS](http://www.unidata.ucar.edu/software/udunits)
standard from Unidata.
"""
function get_time_units end

"""
    get_time_step(model)

Current time step of the model.
"""
function get_time_step end

"""
    get_value(model, name, dest)

Get a copy of values of the given variable.

This is a getter for the model, used to access the model's
current state. It returns a *copy* of a model variable, with
the return type, size and rank dependent on the variable.

# Arguments
- `name`: An input or output variable name, a CSDMS Standard Name.
- `dest`: An array into which to place the values.

# Returns
The same array that was passed as an input buffer.
"""
function get_value end

"""
    get_value_ptr(model, name)

Get a reference to values of the given variable.

This is a getter for the model, used to access the model's
current state. It returns a reference to a model variable,
with the return type, size and rank dependent on the variable.

The `name` can be an input or output variable name, a CSDMS Standard Name.

Returns a reference to a model variable.
"""
function get_value_ptr end

"""
    get_value_at_indices(model, name, dest, inds)

Get values at particular indices.

# Arguments
- `name`: An input or output variable name, a CSDMS Standard Name.
- `dest`: An array into which to place the values.
- `inds`: The indices into the variable array.

# Returns
array_like
    Value of the model variable at the given location.
"""
function get_value_at_indices end

"""
    set_value(model, name, value)

Specify a new value for a model variable.

This is the setter for the model, used to change the model's
current state. It accepts, through *value*, a new value for a
model variable, with the type, size and rank of *value*
dependent on the variable.

# Arguments
- `name`: An input or output variable name, a CSDMS Standard Name.
- `value`: The new value for the specified variable.
"""
function set_value end

"""
    set_value_at_indices(model, name, inds, value)

Specify a new value for a model variable at particular indices.

# Arguments
- `name`: An input or output variable name, a CSDMS Standard Name.
- `inds`: The indices into the variable array.
- `value`: The new value for the specified variable.
"""
function set_value_at_indices end

# Grid information

"""
    get_grid_rank(model, grid)

Get number of dimensions of the computational grid.
"""
function get_grid_rank end

"""
    get_grid_size(model, grid)

Get the total number of elements in the computational grid.
"""
function get_grid_size end

"""
    get_grid_type(model, grid)

Get the grid type as a string.
"""
function get_grid_type end

# Uniform rectilinear

"""
    get_grid_shape(model, grid)

Get dimensions of the computational grid.
"""
function get_grid_shape end

"""
    get_grid_spacing(model, grid, spacing)

Get distance between nodes of the computational grid.

# Arguments
- `grid`: A grid identifier.
- `spacing`: An array to hold the spacing between grid rows and columns.

Returns the filled `spacing` array.
"""
function get_grid_spacing end

"""
    get_grid_origin(model, grid, origin)

Get coordinates for the lower-left corner of the computational grid.

# Arguments
- `grid`: A grid identifier.
- `origin`: An array to hold the coordinates of the lower-left corner of
    the grid.

Returns the filled `origin` array.
"""
function get_grid_origin end

# Non-uniform rectilinear, curvilinear

"""
    get_grid_x(model, grid, x)

Get coordinates of grid nodes in the x direction.

# Arguments
- `grid`: A grid identifier.
- `x`: An array to hold the x-coordinates of the grid nodes.

Returns the filled `x` array.
"""
function get_grid_x end

"""
    get_grid_y(model, grid, y)

Get coordinates of grid nodes in the y direction.

# Arguments
- `grid`: A grid identifier.
- `y`: An array to hold the y-coordinates of the grid nodes.

Returns the filled `y` array.
"""
function get_grid_y end

"""
    get_grid_z(model, grid, z)

Get coordinates of grid nodes in the z direction.

# Arguments
- `grid`: A grid identifier.
- `z`: An array to hold the z-coordinates of the grid nodes.

Returns the filled `z` array.
"""
function get_grid_z end

"""
    get_grid_node_count(model, grid)

Get the number of nodes in the grid.
"""
function get_grid_node_count end

"""
    get_grid_edge_count(model, grid)

Get the number of edges in the grid.
"""
function get_grid_edge_count end

"""
    get_grid_face_count(model, grid)

Get the number of faces in the grid.
"""
function get_grid_face_count end

"""
    get_grid_edge_nodes(model, grid, edge_nodes)

Get the edge-node connectivity.

# Arguments
- `grid`: A grid identifier.
- `edge_nodes`: An array of integers to place the edge-node connectivity.
    For each edge, connectivity is given as node at edge tail,
    followed by node at edge head. Therefore this array must be twice
    the number of nodes long.

Returns the filled `edge_nodes` array.
"""
function get_grid_edge_nodes end

"""
    get_grid_face_edges(model, grid, face_edges)

Get the face-edge connectivity.

# Arguments
- `grid`: A grid identifier.
- `face_edges`: An array of integers in which to place the face-edge connectivity.

Returns the filled `face_edges` array.
"""
function get_grid_face_edges end

"""
    get_grid_face_nodes(model, grid, face_nodes)

Get the face-node connectivity.

# Arguments
- `grid`: A grid identifier.
- `face_nodes`: An array of integers in which to place the face-node connectivity.
    For each face, the nodes (listed in a counter-clockwise direction) that form the
    boundary of the face.

Returns the filled `face_nodes` array.
"""
function get_grid_face_nodes end

"""
    get_grid_nodes_per_face(model, grid, nodes_per_face)

Get the number of nodes for each face.

# Arguments
- `grid`: A grid identifier.
- `nodes_per_face`: An array in which to place the number of edges per face.

Returns the filled `nodes_per_face` array.
"""
function get_grid_nodes_per_face end

end # module


text:
# ------------------------------------------
#   ACE Basis

import ACE1.PairPotentials: PolyPairBasis
import ACE1.Transforms: PolyTransform, MultiTransform

export basis_params, degree_params, transform_params

"""
`basis_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct one of the basis.
All parameters are passed as keyword argument and the kind of
parameters required depend on "type".

## ACE (RPI) basis
Returns a dictionary containing the complete set of parameters
required to construct an ACE basis (`RPIBasis`). All parameters
are passed as keyword argument. If no default is given then
the argument is required.

### Parameters
* `type = "ace"`
* `species` : single species or list of species (mandatory)
* `N` : correlation order, positive integer (mandatory)
* `maxdeg` : maximum degree, positive real number (note the precise
notion of degree is specified by further parameters) (mandatory)
* `r0 = 2.5` : rough estimate for nearest neighbour distance
* `radial = radial_basis_params(; r0 = r0)` : one-particle basis
parameters; cf `?basis_params` of type "radial" for details
* `transform = transform_params(; r0 = r0)` : distance transform
parameters; cf `?transform_params()` for details
* `degree = degree_params()` : class of sparse polynomial degree
to select the basis; see `?degree_params` for details

## Pair basis
Returns a dictionary containing the complete set of parameters
required to construct an pair basis (`PolyPairBasis`). All
parameters are passed as keyword argument.

### Parameters
* `type = "pair"`
* `species` : single species or list of species (mandatory)
* `maxdeg` : maximum degree, positive real number (note the precise
notion of degree is specified by further parameters) (mandatory)
* `r0 = 2.5` : rough estimate for nearest neighbour distance
* `rcut = 5.0`: outer cutoff, Å
* `rin = 0.0`: inner cutoff, Å
* `pcut = 2`: outer cutoff parameter;
      * `pcut=2`: function and first derivative go to zero at  the outer cutoff
      * `pcut=1`: function forced to go through zero at the outer cutoff
      * `pcut=0`: no constraint at the outer cutoff
* `pin = 0`: inner cutoff parameter
      * `pin=2`: function and first derivative go to zero at  the inner cutoff
      * `pin=1`: function forced to go through zero at the inner cutoff
      * `pin=0`: no constraint at the inner cutoff
* `transform = transform_params(; r0 = r0)` : distance transform
parameters; cf `?transform_params()` for details

## Radial basis of ACE
Returns a dictionary containing the complete set of parameters
required to construct radial basis for ACE. All parameters are
passed as keyword argument.

### Parameters
* `type = "radial"`
* `r0 = 2.5` : rough estimate for nearest neighbour distance
* `rcut = 5.0`: outer cutoff, Å
* `rin = 0.5 * r0`: inner cutoff, Å
* `pcut = 2`: outer cutoff parameter;
      * `pcut=2`: function and first derivative go to zero at  the outer cutoff
      * `pcut=1`: function forced to go through zero at the outer cutoff
      * `pcut=0`: no constraint at the outer cutoff
* `pin = 2`: inner cutoff parameter
      * `pin=2`: function and first derivative go to zero at  the inner cutoff
      * `pin=1`: function forced to go through zero at the inner cutoff
      * `pin=0`: no constraint at the inner cutoff
"""
function basis_params(;
      type = nothing,
      kwargs...)
      @assert haskey(_bases, type) "type $(type) not found among available types of basis ($(keys(_bases)))"
      return _bases[type][2](; kwargs...)
end

function generate_basis(params::Dict)
      @assert params["type"] != "radial"
      params = copy(params)
      basis_constructor = _bases[params["type"]][1]
      delete!(params, "type")
      return basis_constructor(params)
end

# ------------------------------------------
#  ace basis

"""
`ace_basis_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct an ACE basis (`RPIBasis`).
All parameters are passed as keyword argument. If no default is given then
the argument is required.

### Parameters
* `species` : single species or list of species (mandatory)
* `N` : correlation order, positive integer (mandatory)
* `maxdeg` : maximum degree, positive real number (note the precise notion of
degree is specified by further parameters) (mandatory)
* `r0 = 2.5` : rough estimate for nearest neighbour distance
* `radial = radial_basis_params(; r0 = r0)` : one-particle basis parameters;
cf `?radial_basis_params` for details
* `transform = transform_params(; r0 = r0)` : distance transform parameters;
cf `?transform_params()` for details
* `degree = degree_params()` : class of sparse polynomial degree to select
the basis; see `?degree_params` for details
"""
function ace_basis_params(;
      species = nothing,
      N::Integer = nothing,
      maxdeg = nothing,
      r0 = 2.5,
      radial = radial_basis_params(; r0 = r0),
      transform = transform_params(; r0 = r0),
      degree = degree_params(),
      type = "ace"
   )

   @assert !isnothing(species) "`species` is mandatory for `ace_basis_params`"
   @assert isinteger(N) "correlation order `N` must be a positive integer"
   @assert N > 0 "correlation order `N` must be a positive integer"
   @assert isreal(maxdeg) "Maximum polynomial degree `maxdeg` must be a real positive number"
   @assert maxdeg > 0 "Maximum polynomial degree `maxdeg` must be a real positive number"
   @assert isreal(r0) "`r0` must be a real positive number "
   @assert r0 > 0 "`r0` must be a real positive number "
   @assert type == "ace" "`type` must be set to \"ace\" for `ace_basis_params`"

   if !haskey(radial, "type")
      radial = convert(Dict{String, Any}, radial)
      radial["type"] = "radial"
   end

   return Dict(
         "type" => "ace",
         "species" => _species_to_params(species),
         "N" => N,
         "maxdeg" => maxdeg,
         "radial" => radial,
         "transform" => transform,
         "degree" => degree
         )
end

"""Returns ACE1.Utils.rpi_basis """
function generate_ace_basis(params::Dict)
   species = _params_to_species(params["species"])
   trans = generate_transform(params["transform"])
   D = generate_degree(params["degree"])
   maxdeg = params["maxdeg"]
   radial = generate_radial_basis(params["radial"], D, maxdeg, species, trans)
   return ACE1.Utils.rpi_basis(;
            species = species,
            N = params["N"],
            trans = trans,
            D = D,
            maxdeg = maxdeg,
            rbasis = radial,
         )
end

# ------------------------------------------
#  pair basis

"""
`pair_basis_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct an pair basis (`PolyPairBasis`).
All parameters are passed as keyword argument.

### Parameters
* `species` : single species or list of species (mandatory)
* `maxdeg` : maximum degree, positive real number (note the precise notion of degree
is specified by further parameters) (mandatory)
* `r0 = 2.5` : rough estimate for nearest neighbour distance
* `rcut = 5.0`: outer cutoff, Å
* `rin = 0.0`: inner cutoff, Å
* `pcut = 2`: outer cutoff parameter;
      * `pcut=2`: function and first derivative go to zero at  the outer cutoff
      * `pcut=1`: function forced to go through zero at the outer cutoff
      * `pcut=0`: no constraint at the outer cutoff
* `pin = 0`: inner cutoff parameter
      * `pin=2`: function and first derivative go to zero at  the inner cutoff
      * `pin=1`: function forced to go through zero at the inner cutoff
      * `pin=0`: no constraint at the inner cutoff
* `transform = transform_params(; r0 = r0)` : distance transform parameters;
cf `?transform_params()` for details
"""
function pair_basis_params(;
      species = nothing,
      maxdeg = nothing,
      r0 = 2.5,
      rcut = 5.0,
      rin = 0.0,
      pcut = 2,
      pin = 0,
      transform = transform_params(; r0=r0),
      type = "pair",
      )

      @assert !isnothing(species) "`species` is mandatory for `ace_basis_params`"
      @assert isreal(maxdeg) "Maximum polynomial degree `maxdeg` must be a real positive number"
      @assert maxdeg > 0 "Maximum polynomial degree `maxdeg` must be a real positive number"
      @assert isreal(r0) "`r0` must be a real positive number "
      @assert r0 > 0 "`r0` must be a real positive number "
      @assert type == "pair" "`type` must be set to \"pair\" for `pair_basis_params`"

      return Dict(
            "type" => "pair",
            "species" => _species_to_params(species),
            "maxdeg" => maxdeg,
            "rcut" => rcut,
            "rin" => rin,
            "pcut" => pcut,
            "pin" => pin,
            "transform" => transform)
end

""" Returns PolyPairBasis """
function generate_pair_basis(params::Dict)
      species = _params_to_species(params["species"])
      trans = generate_transform(params["transform"])
      rad_basis = transformed_jacobi(params["maxdeg"], trans, params)
      return PolyPairBasis(rad_basis, species)

end

# ------------------------------------------
#  rad_basis

"""
`radial_basis_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct radial basis for ACE.
All parameters are passed as keyword argument.

### Parameters
* `r0 = 2.5` : rough estimate for nearest neighbour distance
* `rcut = 5.0`: outer cutoff, Å
* `rin = 0.5 * r0`: inner cutoff, Å
* `pcut = 2`: outer cutoff parameter;
      * `pcut=2`: function and first derivative go to zero at  the outer cutoff
      * `pcut=1`: function forced to go through zero at the outer cutoff
      * `pcut=0`: no constraint at the outer cutoff
* `pin = 2`: inner cutoff parameter
      * `pin=2`: function and first derivative go to zero at  the inner cutoff
      * `pin=1`: function forced to go through zero at the inner cutoff
      * `pin=0`: no constraint at the inner cutoff
"""
function radial_basis_params(;
      r0 = 2.5,
      rcut = 5.0,
      rin = 0.5 * r0,
      pcut = 2,
      pin = 2,
      type = "radial")

      @assert isreal(r0) "`r0` must be a real positive number "
      @assert r0 > 0 "`r0` must be a real positive number "
      @assert type == "radial" "`type` must be set to \"radial\" for `radial_basis_params`"

      return Dict(
            "type" => "radial",
            "rcut" => rcut,
            "rin" => rin,
            "pcut" => pcut,
            "pin" => pin )
end

function generate_radial_basis(params::Dict, D, maxdeg, species, trans)
   maxn = ACE1.RPI.get_maxn(D, maxdeg, species)
   return transformed_jacobi(maxn, trans, params)
end

# ------------------------------------------
#  basis helper functions

_bases = Dict("pair" => (generate_pair_basis, pair_basis_params),
              "ace" => (generate_ace_basis, ace_basis_params),
              "radial" => (nothing, radial_basis_params))

transformed_jacobi(maxn::Integer, trans::MultiTransform, params::Dict) =
      OrthPolys.transformed_jacobi(maxn, trans; pcut = params["pcut"], pin = params["pin"])

transformed_jacobi(maxn::Integer, trans::PolyTransform, params::Dict) =
      OrthPolys.transformed_jacobi(maxn, trans, params["rcut"], params["rin"];
                         pcut = params["pcut"], pin = params["pin"])

# ------------------------------------------
#  degree

"""
`degree_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct a specification
for polynomial degree. All parameters are passed as keyword argument
and the kind of parameters required depend on "type".

## `SparsePSHDegree`
Returns a dictionary containing the complete set of parameters required
to construct `ACE1.RPI.SparsePSHDegree`. See `?SparsePSHDegree`.

### Parameters
* `type = "sparse"`
* `wL = 1.5`
* `csp = 1.0`
* `chc = 0.0`
* `chc = 0.0`
* `ahc = 0.0`
* `bhc = 0.0`
* `p = 1.0`

## `SparsePSHDegreeM`
Returns a dictionary containing the complete set of parameters required
to construct `ACE1.RPI.SparsePSHDegree`. Also see `?SparsePSHDegreeM`.

NB `maxdeg` of ACE basis (`RPIBasis`) has to be set to `1.0`.

### Parameters
* `Dd` : Dictionary specifying max degrees (mandatory)
* `Dn = Dict("default" => 1.0)` : Dictionary specifying weights for degree
of radial basis functions (n)
* `Dl = Dict("default" => 1.5)` : Dictionary specifying weights for degree
of angular basis functions (l)

Each dictionary should have a "default" entry. In addition, different degrees
or weights can be specified for each correlation order and/or correlation
order-species combination. For example

```
"Dd" => Dict(
      "default" => 10,
      3 => 9,
      (4, "C") => 8,
      (4, "H") => 0)
```

in combination with N=4 and maxdeg=1.0, will set maximum polyonmial degree on
N=1 and N=2 functions to 10, to 9 for N=3 functions and will only allow
N=4 basis functions on carbon atoms, up to polynomial degree 8.

"""
function degree_params(;
      type = "sparse",
      kwargs...)
      @assert haskey(_degrees, type)
      return _degrees[type][2](; kwargs...)
end

function generate_degree(params::Dict)
      @assert haskey(_degrees, params["type"])
      # we ignore p for `SparsePSHDegree`, for now
      if params["type"] == "sparse" && haskey(params, "p")
            delete!(params, "p")
      end
      degree_measure = _degrees[params["type"]][1]
      kwargs = Dict([Symbol(key) => val for (key, val) in params]...)
      delete!(kwargs, :type)
      return degree_measure(; kwargs...)
end

"""
`sparse_degree_params(; kwargs...)`: returns a dictionary containing the
complete set of parameters required to construct `ACE1.RPI.SparsePSHDegree`.
See `?SparsePSHDegree`.

### Parameters
* `wL = 1.5`
* `csp = 1.0`
* `chc = 0.0`
* `chc = 0.0`
* `ahc = 0.0`
* `bhc = 0.0`
* `p = 1.0`

NB `p = 1` is current ignored, but we put it in so we can experiment later
with `p = 2`, `p = inf`.
"""
function sparse_degree_params(;
      wL::Real = 1.5,
      csp::Real = 1.0,
      chc::Real = 0.0,
      ahc::Real = 0.0,
      bhc::Real = 0.0,
      p::Real = 1.0 )

   @assert wL > 0
   @assert csp >= 0 && chc >= 0
   @assert csp > 0 || chc > 0
   @assert ahc >= 0 && bhc >= 0
   @assert p == 1
   return Dict( "type" => "sparse",
                "wL" => wL,
                "csp" => csp,
                "chc" => chc,
                "ahc" => ahc,
                "bhc" => bhc,
                "p" => p )
end

"""
`sparse_degree_M_params(;kwargs...)`: Returns a dictionary containing the
complete set of parameters required to construct `ACE1.RPI.SparsePSHDegree`.
Also see `?SparsePSHDegreeM`.

NB `maxdeg` of ACE basis (`RPIBasis`) has to be set to `1.0`.

### Parameters
* `Dd` : Dictionary specifying max degrees (mandatory)
* `Dn = Dict("default" => 1.0)` : Dictionary specifying weights for degree
of radial basis functions (n)
* `Dl = Dict("default" => 1.5)` : Dictionary specifying weights for degree
of angular basis functions (l)

Each dictionary should have a "default" entry. In addition, different degrees
or weights can be specified for each correlation order and/or correlation
order-species combination. For example

```
"Dd" => Dict(
      "default" => 10,
      3 => 9,
      (4, "C") => 8,
      (4, "H") => 0)
```

in combination with N=4 and maxdeg=1.0, will set maximum polyonmial degree on
N=1 and N=2 functions to 10, to 9 for N=3 functions and will only allow
N=4 basis functions on carbon atoms, up to polynomial degree 8.
"""
function sparse_degree_M_params(;
      Dd::Dict = nothing,
      Dn::Dict = Dict("default" => 1.0),
      Dl::Dict = Dict("default" => 1.5))

      @assert !isnothing(Dd) "`Dd`` is a mandatory."

      return Dict(
            "type" => "sparseM",
            "Dd" => _AtomicNumber_to_params(Dd),
            "Dn" => _AtomicNumber_to_params(Dn),
            "Dl" => _AtomicNumber_to_params(Dl))
end

SparsePSHDegreeM(; Dn::Dict, Dl::Dict, Dd::Dict) =
      ACE1.RPI.SparsePSHDegreeM(_params_to_AtomicNumber(Dn),
                                _params_to_AtomicNumber(Dl),
                                _params_to_AtomicNumber(Dd))

_degrees = Dict(
      "sparse" => (ACE1.RPI.SparsePSHDegree, sparse_degree_params),
      "sparseM" => (SparsePSHDegreeM, sparse_degree_M_params))

# ------------------------------------------
#  transform
#  this is a little more interesting since there are quite a
#  few options.

"""
`transform_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct one of the transforms.
All parameters are passed as keyword argument and the kind of
parameters required depend on "type".

## Polynomial transform
Returns a dictionary containing the complete set of parameters required
to construct `ACE1.Transforms.PolyTransform``. All parameters are passed
 as keyword argument. Also see `?PolyTransform`

Implements the distance transform

```math
   x(r) = \\Big(\\frac{1 + r_0}{1 + r}\\Big)^p
```

### Parameters
* `type = "polynomial"`
* `p = 2`
* `r0 = 2.5`

## Multitransform
Returns a dictionary containing the complete set of parameters required
to construct `ACE.Transform.multitransform`. All parameters are passed
as keyword argument.

### Parameters
* `transforms` : dictionary specifying transforms for each species pair. Can be
given per-pair (i.e. only for `("element1", "element2")` and not for
`("element2", "element1")`) or can be different for `("element1", "element2")` and
`("element2", "element1")`. For example
```
transforms = Dict(
      ("C", "C") => Dict("type"=> "polynomial"),
      ("C", "H") => Dict("type"=> "polynomial"),
      ("H", "H") => Dict("type" => "polynomial"))
```

* `rin`, `rcut`: values for inner and outer cutoffs, alternative to `cutoffs`
* `cutoffs` : dictionary specifying inner and outer cutoffs for each element pair
(either symmetrically or non-symmetrically). Alternative to `rin` & `rcut`.
For example
```
cutoffs => Dict(
      ("C", "C") => (1.1, 4.5),
      ("C", "H") => (0.9, 4.5),
      ("H", "H") => (1.23, 4.5)),
```

## identity
`IdTransform_params(;)` : returns `Dict("type" => "identity")`,
needed to construct `ACE1.Transforms.IdTransform`.
"""
function transform_params(;
      type = "polynomial",
      kwargs...)
   @assert haskey(_transforms, type)
   return _transforms[type][2](; kwargs...)
end

"""
`PolyTransform_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct `ACE1.Transforms.PolyTransform``.
All parameters are passed as keyword argument. Also see `?PolyTransform`

Implements the distance transform

```math
   x(r) = \\Big(\\frac{1 + r_0}{1 + r}\\Big)^p
```

### Parameters
* `p = 2`
* `r0 = 2.5`
"""
function PolyTransform_params(; p = 2, r0 = 2.5)
   @assert isreal(p) "`p`` must be a real positive number"
   @assert p > 0 "`p`` must be a real positive number"
   @assert isreal(r0) "`r0` must be a real positive number "
   @assert r0 > 0 "`r0` must be a real positive number "

   return Dict("type" => "polynomial",
               "p" => p,
               "r0" => r0)
end

"""
`multitransform_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to construct `ACE.Transform.multitransform`.
All parameters are passed as keyword argument.

### Parameters
* `transforms` : dictionary specifying transforms for each species pair. Can be
given per-pair (i.e. only for `("element1", "element2")` and not for
`("element2", "element1")`) or can be different for `("element1", "element2")` and
`("element2", "element1")`. For example
```
transforms = Dict(
      ("C", "C") => Dict("type"=> "polynomial"),
      ("C", "H") => Dict("type"=> "polynomial"),
      ("H", "H") => Dict("type" => "polynomial"))
```

* `rin`, `rcut`: values for inner and outer cutoffs, alternative to `cutoffs`
* `cutoffs` : dictionary specifying inner and outer cutoffs for each element pair
(either symmetrically or non-symmetrically). Alternative to `rin` & `rcut`.
For example
```
cutoffs => Dict(
      ("C", "C") => (1.1, 4.5),
      ("C", "H") => (0.9, 4.5),
      ("H", "H") => (1.23, 4.5)),
```
"""
function multitransform_params(;
      transforms = nothing,
      rin = nothing,
      rcut = nothing,
      cutoffs=nothing)

      @assert !isnothing(transforms) "`transforms` must be specified."
      @assert (!isnothing(rin) && !isnothing(rcut)) || !isnothing(cutoffs) "Either `rin` & `rcut` or `cutoffs` must be given."

      return Dict("type" => "multitransform",
                  "transforms" => _species_to_params(transforms),
                  "rin" => rin,
                  "rcut" => rcut,
                  "cutoffs" => _species_to_params(cutoffs))
end

function generate_multitransform(; transforms, rin = nothing, rcut = nothing, cutoffs = nothing)
      transforms = _params_to_species(transforms)
      transforms = Dict(key => generate_transform(params) for (key, params) in transforms)
      cutoffs = _params_to_species(cutoffs)
      return ACE1.Transforms.multitransform(transforms, rin = rin, rcut = rcut, cutoffs = cutoffs)
end

"""
`IdTransform_params(;)` : returns `Dict("type" => "identity")`,
needed to construct `ACE1.Transforms.IdTransform`.
"""
IdTransform_params(;) = Dict("type" => "identity")

function generate_transform(params::Dict)
   @assert haskey(_transforms, params["type"])
   TTransform = _transforms[params["type"]][1]
   kwargs = Dict([Symbol(key) => val for (key, val) in params]...)
   delete!(kwargs, :type)
   return TTransform(; kwargs...)
end

# In this dictionary we "register" all the transforms for which we have
# supplied an interface. At the moment I've done it just for one of them
# others can introduce more. The key is a string the specifies the key
# user supplies for the `type` parameter. The value is a tuple containing
# the corresponding transform type and function that generates the defaul
# parameters
_transforms = Dict(
      "polynomial" => (ACE1.Transforms.PolyTransform, PolyTransform_params),
      "multitransform" => (generate_multitransform, multitransform_params),
      "identity" => (ACE1.Transforms.IdTransform, IdTransform_params)
)

# ------------------------------------------
# helper functions

# -- Symbol to string
_species_to_params(species::Union{Symbol, AbstractString}) =
      [ string(species), ]

_species_to_params(species::Union{Tuple, AbstractArray}) =
      collect( string.(species) )

# accept tuples of Symbol or String for dictionary key
# values can be anything
_species_to_params(dict::Dict{Tuple{Tsym, Tsym}, Tval}) where Tsym <: Union{Symbol, AbstractString} where Tval <: Any =
      Dict(Tuple(_species_to_params(key)) => val for (key, val) in dict)

_species_to_params(dict::Nothing) = nothing

# -- String to Symbol
_params_to_species(species::Union{AbstractArray{T}, Tuple{T, T}}) where T <: AbstractString  =
      Symbol.(species)

_params_to_species(dict::Dict{Tuple{Tsym, Tsym}, Tval}) where Tsym <: AbstractString where Tval <: Any =
      Dict(Tuple(_params_to_species(d)) => val for (d, val) in dict)

_params_to_species(dict::Nothing) = nothing

function _AtomicNumber_to_params(dict)
      new_dict = Dict()
      for (key, val) in dict
            if typeof(key) <: Tuple
                 key = Tuple(typeof(entry) <: AtomicNumber ? string(chemical_symbol(entry)) : entry for entry in key)
            end
            new_dict[key] = val
      end
      return new_dict
end

function _params_to_AtomicNumber(dict)
      new_dict = Dict()
      for (key, val) in dict
            if typeof(key) <: Tuple
                  key = Tuple(typeof(entry) <: AbstractString && length(entry) == 1 ?
                              AtomicNumber(Symbol(entry)) : entry for entry in key)
            end
            new_dict[key] = val
      end
      return new_dict
end

ext: jl


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

"""
A QuadraticProgrammingProblem struct specifies a quadratic programming problem
with the following format:

```
minimize 1/2 x' * objective_matrix * x + objective_vector' * x
          + objective_constant

s.t. constraint_matrix[1:num_equalities, :] * x =
     right_hand_side[1:num_equalities]

     constraint_matrix[(num_equalities + 1):end, :] * x >=
     right_hand_side[(num_equalities + 1):end]

     variable_lower_bound <= x <= variable_upper_bound
```
The variable_lower_bound may contain `-Inf` elements and variable_upper_bound
may contain `Inf` elements when the corresponding variable bound is not present.
"""
mutable struct QuadraticProgrammingProblem
  """
  The vector of variable lower bounds.
  """
  variable_lower_bound::Vector{Float64}

  """
  The vector of variable upper bounds.
  """
  variable_upper_bound::Vector{Float64}

  """
  The symmetric and positive semidefinite matrix that defines the quadratic
  term in the objective.
  """
  objective_matrix::SparseMatrixCSC{Float64,Int64}

  """
  The linear coefficients of the objective function.
  """
  objective_vector::Vector{Float64}

  """
  The constant term of the objective function.
  """
  objective_constant::Float64

  """
  The matrix of coefficients in the linear constraints.
  """
  constraint_matrix::SparseMatrixCSC{Float64,Int64}

  """
  The vector of right-hand side values in the linear constraints.
  """
  right_hand_side::Vector{Float64}

  """
  The number of equalities in the problem. This value splits the rows of the
  constraint_matrix between the equality and inequality parts.
  """
  num_equalities::Int64
end

"""
Estimates the variable and constraint hardness for a first-order method given
a solution vector. By hardness we mean roughly how much does the constraint or
variable contribute to the worst-case runtime bound. This is useful for
debugging situations where the method is performing poorly.
"""
function print_variable_and_constraint_hardness(
  problem::QuadraticProgrammingProblem,
  primal_solution::AbstractVector{Float64},
  dual_solution::AbstractVector{Float64},
)

  row_norms = get_row_l2_norms(problem.constraint_matrix)
  constraint_hardness = row_norms .* abs.(dual_solution)
  print("Constraint hardness: ")
  Printf.@printf(
    "median_hardness=%f, mean_hardness=%f, quantile_99=%f, hardest=%f\n",
    median(constraint_hardness),
    mean(constraint_hardness),
    quantile(constraint_hardness, 0.99),
    maximum(constraint_hardness)
  )

  col_norms = get_col_l2_norms(problem.constraint_matrix)
  variable_hardness = col_norms .* abs.(primal_solution)
  print("Variable hardness: ")
  Printf.@printf(
    "median_hardness=%f, mean_hardness=%f, quantile_99=%f, hardest=%f\n",
    median(variable_hardness),
    mean(variable_hardness),
    quantile(variable_hardness, 0.99),
    maximum(variable_hardness)
  )
end

function get_row_l2_norms(matrix::SparseMatrixCSC{Float64,Int64})
  row_norm_squared = zeros(size(matrix, 1))
  nzval = nonzeros(matrix)
  rowval = rowvals(matrix)
  for i in 1:length(nzval)
    row_norm_squared[rowval[i]] += nzval[i]^2
  end

  return sqrt.(row_norm_squared)
end

function get_col_l2_norms(matrix::SparseMatrixCSC{Float64,Int64})
  col_norms = zeros(size(matrix, 2))
  for j in 1:size(matrix, 2)
    col_norms[j] = norm(nonzeros(matrix[:, j]), 2)
  end
  return col_norms
end

function get_row_l_inf_norms(matrix::SparseMatrixCSC{Float64,Int64})
  row_norm = zeros(size(matrix, 1))
  nzval = nonzeros(matrix)
  rowval = rowvals(matrix)
  for i in 1:length(nzval)
    row_norm[rowval[i]] = max(abs(nzval[i]), row_norm[rowval[i]])
  end

  return row_norm
end

function get_col_l_inf_norms(matrix::SparseMatrixCSC{Float64,Int64})
  col_norms = zeros(size(matrix, 2))
  for j in 1:size(matrix, 2)
    col_norms[j] = norm(nonzeros(matrix[:, j]), Inf)
    typeof(matrix[:, j])
  end
  return col_norms
end

"""
  print_problem_details(qp)

This is primarily useful for detecting when a problem is poorly conditioned and
needs rescaling.
"""
function print_problem_details(qp::QuadraticProgrammingProblem)
  println(
    "  There are ",
    size(qp.constraint_matrix, 2),
    " variables, ",
    size(qp.constraint_matrix, 1),
    " constraints (including ",
    qp.num_equalities,
    " equalities) and ",
    SparseArrays.nnz(qp.constraint_matrix),
    " nonzero coefficients.",
  )

  print("  Absolute value of nonzero constraint matrix elements: ")
  Printf.@printf(
    "largest=%f, smallest=%f, avg=%f\n",
    maximum(abs, nonzeros(qp.constraint_matrix)),
    minimum(abs, nonzeros(qp.constraint_matrix)),
    sum(abs, nonzeros(qp.constraint_matrix)) /
    length(nonzeros(qp.constraint_matrix))
  )

  col_norms = get_col_l_inf_norms(qp.constraint_matrix)
  row_norms = get_row_l_inf_norms(qp.constraint_matrix)

  print("  Constraint matrix, infinity norm: ")
  Printf.@printf(
    "max_col=%f, min_col=%f, max_row=%f, min_row=%f\n",
    maximum(col_norms),
    minimum(col_norms),
    maximum(row_norms),
    minimum(row_norms)
  )

  if length(nonzeros(qp.objective_matrix)) > 0
    print("  Absolute value of objective matrix elements: ")
    Printf.@printf(
      "largest=%f, smallest=%f, avg=%f\n",
      maximum(abs, nonzeros(qp.objective_matrix)),
      minimum(abs, nonzeros(qp.objective_matrix)),
      sum(abs, nonzeros(qp.objective_matrix)) /
      length(nonzeros(qp.constraint_matrix))
    )
  end

  print("  Absolute value of objective vector elements: ")
  Printf.@printf(
    "largest=%f, smallest=%f, avg=%f\n",
    maximum(abs, qp.objective_vector),
    minimum(abs, qp.objective_vector),
    sum(abs, qp.objective_vector) / length(qp.objective_vector)
  )

  print("  Absolute value of rhs vector elements: ")
  Printf.@printf(
    "largest=%f, smallest=%f, avg=%f\n",
    maximum(abs, qp.right_hand_side),
    minimum(abs, qp.right_hand_side),
    sum(abs, qp.right_hand_side) / length(qp.right_hand_side)
  )

  bound_gaps = qp.variable_upper_bound - qp.variable_lower_bound
  finite_bound_gaps = bound_gaps[isfinite.(bound_gaps)]

  print("  Gap between upper and lower bounds: ")
  Printf.@printf(
    "#finite=%i of %i, largest=%f, smallest=%f, avg=%f\n",
    length(finite_bound_gaps),
    length(bound_gaps),
    length(finite_bound_gaps) > 0 ? maximum(finite_bound_gaps) : NaN,
    length(finite_bound_gaps) > 0 ? minimum(finite_bound_gaps) : NaN,
    length(finite_bound_gaps) > 0 ?
    sum(finite_bound_gaps) / length(finite_bound_gaps) : NaN
  )
end

"""
Creates a QuadraticProgrammingProblem with the objective matrix set to zero.

# Arguments
- `variable_lower_bound::Vector{Float64}`: The lower bound of variables.
- `variable_upper_bound::Vector{Float64}`: The upper bound of variables.
- `objective_vector::Vector{Float64}`: The linear coefficients of
  objective function.
- `objective_constant::Float64`: The constant of objective function
- `variable_upper_bound::Vector{Float64}`: The constant of objective
  function.
- `constraint_matrix::Union{SparseMatrixCSC{Float64,Int64},Array{Float64,2}}`:
  The matrix of coefficients in the linear constraints.
- `right_hand_side::Vector{Float64}`: The vector of right-hand side values in
  the linear constraints.
- `num_equalities::Int64`: The number of equalities in the problem. This value
  splits the rows of the constraint_matrix between the equality and inequality
  parts.

"""
function linear_programming_problem(
  variable_lower_bound::Vector{Float64},
  variable_upper_bound::Vector{Float64},
  objective_vector::Vector{Float64},
  objective_constant::Float64,
  constraint_matrix::Union{SparseMatrixCSC{Float64,Int64},Array{Float64,2}},
  right_hand_side::Vector{Float64},
  num_equalities::Int64,
)

  num_variables = length(variable_lower_bound)
  objective_matrix = spzeros(num_variables, num_variables)
  return QuadraticProgrammingProblem(
    variable_lower_bound,
    variable_upper_bound,
    objective_matrix,
    objective_vector,
    objective_constant,
    sparse(constraint_matrix),
    right_hand_side,
    num_equalities,
  )
end

"""
Returns true if and only if the objective matrix is zero.
"""
function is_linear_programming_problem(problem::QuadraticProgrammingProblem)
  return nnz(problem.objective_matrix) == 0
end

"""
A ScaledQpProblem struct specifies a original quadratic programming problem,
a scaled quadratic programming problem, and the scaling vector, which requires
to satisfy the condition that:

orginal_qp = unscale_problem(scaled_qp, constraint_rescaling,variable_rescaling)
"""
mutable struct ScaledQpProblem
  original_qp::QuadraticProgrammingProblem
  scaled_qp::QuadraticProgrammingProblem
  constraint_rescaling::Vector{Float64}
  variable_rescaling::Vector{Float64}
end

equality_range(problem::QuadraticProgrammingProblem) = 1:problem.num_equalities

function inequality_range(problem::QuadraticProgrammingProblem)
  return (problem.num_equalities+1):size(problem.constraint_matrix, 1)
end


text:
module PandorasBoxFlagToggleSwitch

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/flagToggleSwitch" FlagToggleSwitch(x::Integer, y::Integer, onlyOff::Bool=false, onlyOn::Bool=false, silent::Bool=false, flag::String="")

const placements = Ahorn.PlacementDict(
    "Flag Toggle Switch (Activate) (Pandora's Box)" => Ahorn.EntityPlacement(
        FlagToggleSwitch,
        "point",
        Dict{String, Any}(
            "onlyOn" => true
        )
    ),
    "Flag Toggle Switch (Deactivate) (Pandora's Box)" => Ahorn.EntityPlacement(
        FlagToggleSwitch,
        "point",
        Dict{String, Any}(
            "onlyOff" => true
        )
    ),
    "Flag Toggle Switch (Pandora's Box)" => Ahorn.EntityPlacement(
        FlagToggleSwitch
    ),
)

function switchSprite(entity::FlagToggleSwitch)
    onlyOff = get(entity.data, "onlyOff", false)
    onlyOn = get(entity.data, "onlyOn", false)

    if onlyOff
        return "objects/pandorasBox/flagToggleSwitch/switch13"

    elseif onlyOn
        return "objects/pandorasBox/flagToggleSwitch/switch15"

    else
        return "objects/pandorasBox/flagToggleSwitch/switch01"
    end
end

function Ahorn.selection(entity::FlagToggleSwitch)
    x, y = Ahorn.position(entity)
    sprite = switchSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleSwitch, room::Maple.Room)
    sprite = switchSprite(entity)

    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
path = joinpath(dirname(@__FILE__), "..", "..")
push!(Base.LOAD_PATH, path)
using ProxSDP, JuMP

# Create a JuMP model using ProxSDP as the solver
model = Model(with_optimizer(ProxSDP.Optimizer, log_verbose=true, tol_primal = 1e-5, tol_dual = 1e-5))

@variable(model, X[1:2,1:2], PSD)

x = X[1,1]
y = X[2,2]

@constraint(model, ub_x, x <= 2)
@constraint(model, ub_y, y <= 30)
@constraint(model, con, 1x + 5y <= 3)

# ProxSDP supports maximization or minimization
# of linear functions
@objective(model, Max, 5x + 3 * y)

# Then we can solve the model
JuMP.optimize!(model)

# And ask for results!
JuMP.objective_value(model)

JuMP.value(x)

JuMP.value(y)


text:
"Tableau for symplectic Euler-A method"
function getTableauSymplecticEulerA()
    a_q = [[0.0 0.0]
           [1.0 0.0]]
    b_q = [1.0, 0.0]
    c_q = [0.0, 1.0]

    a_p = [[0.0 0.0]
           [0.0 1.0]]
    b_p = [0.0, 1.0]
    c_p = [0.0, 1.0]

    o = 1

    TableauEPRK(:symplectic_euler_a, o, a_q, a_p, b_q, b_p, c_q, c_p)
end

"Tableau for symplectic Euler-B method"
function getTableauSymplecticEulerB()
    a_q = [[0.0 0.0]
           [0.0 1.0]]
    b_q = [0.0, 1.0]
    c_q = [0.0, 1.0]

    a_p = [[0.0 0.0]
           [1.0 0.0]]
    b_p = [1.0, 0.0]
    c_p = [0.0, 1.0]

    o = 1

    TableauEPRK(:symplectic_euler_b, o, a_q, a_p, b_q, b_p, c_q, c_p)
end