text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

const calculus = [
    :df,
    :int,
    :limit,
    :logb,
    :solve,
    :pf,
    :structr,
    :coeff,
    :coeffn,
    :part,
    :factorize,
    :remainder,
    :resultant,
    :interpol,
    :deg,
    :lcof,
    :lpower,
    :lterm,
    :reduct,
    :totaldeg,
    :pochhammer,
    :fibonaccip,
    :remfac,
    :gcd,
    :lcm,
    :compact,
    :changevar,
    :odesolve,
]

const cnan = [
    :clear,
    :matrix,
    :operator,
    :listargp,
    :infix,
    :precedence,
    :depend,
    :nodepend,
    :realvalued,
    :notrealvalued,
    :set,
    :unset,
    :mkid,
    :even,
    :odd,
    :linear,
    :noncom,
    :symmetric,
    :antisymmetric,
    :order,
    :korder,
    :on,
    :off,
    :define,
    :select,
    :nospur,
    :spur,
    :vector,
    :index,
    :reminx,
    :g,
    :mass,
    :mshell,
]

const alg = [
    :sum,
    :prod,
    :max,
    :min,
    :map,
    :eps,
]

const iops = [
    :+,
    :-,
    :*,
    :/,
    :^,
    ://
]

const cmat = [
    #:mateigen,
    #:cofactor
]

Expr(:block,[:(@inline $i(r...)=Base.$i(r...)) for i ∈ [alg;iops[6:end]]]...) |> eval
#Expr(:toplevel,[:(import Base: $i) for i ∈ [alg;iops]]...) |> eval
:(export $([calculus;cnan;alg;iops;cmat]...)) |> eval
#:(export $(Symbol.("@",[calculus;alg;iops])...)) |> eval

for fun in [calculus;alg;iops]
    @eval begin
        $(Reduce.parsegen(fun,:args))
        $(Reduce.unfoldgen(fun,:args))
        #=macro $fun(expr,s...)
            :($$(QuoteNode(fun))($(esc(expr)),$(esc(s))...))
        end=#
    end
end

for fun in [calculus;alg]
    @eval begin
        function $fun(expr::String,s...;be=0)
            convert(String, $fun(RExpr(expr),s...;be=be))
        end
    end
end

for fun in cnan
    @eval begin
        $fun(r::RExpr...) = string($(string(fun)),"(",join(string.(r),","),")") |> rcall |> RExpr
        $fun(r...) = $fun(RExpr.(r)...) |> parse
    end
end

for fun in cmat
    @eval begin
        $(Reduce.parsegen(fun,:args))
        function $fun(expr::Union{Array{Any,2},Expr,Symbol},s...)
            $fun(RExpr(expr),RExpr.(s)...)
        end
        function $fun(expr::Union{Array{T,2},Expr,Symbol},s...) where T <: ExprSymbol
            $fun(RExpr(expr),RExpr.(s)...)
        end
    end
end

for fun in iops
    @eval begin
        function $fun(a::Union{<:Number,Expr,Symbol},r::RExpr,s...)
            $fun(RExpr(a),r,RExpr.(s)...) |> parse
        end
        function $fun(a::T,b::ExprSymbol,s...) where T <: Number
            $fun(RExpr(a),RExpr(b),RExpr.(s)...) |> parse
        end
        function $fun(a::A,b::B,c::ExprSymbol,s...) where {A<:Number,B<:Number}
            $fun(RExpr(a),RExpr(b),RExpr(c),RExpr.(s)...) |> parse
        end
    end
end

const MatExpr = Union{Array{Any,2},Array{Expr,2},Array{Symbol,2},Expr,Symbol,<:Number}
const Mat = Union{Vector,Adjoint,Array{Any,2},Array{Expr,2},Array{Symbol,2}}
const MatOnly = Union{Array{Any,2},Array{Expr,2},Array{Symbol,2}}
const ESN = Union{Expr,Symbol,<:Number}

^(expr::Array{Any,2},s::Integer) = ^(RExpr(expr),s) |> parse |> mat
^(expr::Array{T,2},s::Integer) where T <: ExprSymbol = ^(RExpr(expr),s) |> parse |> mat
^(expr::RExpr,s::Integer) = ^(expr,RExpr(s))
^(a::T,s::ExprSymbol) where T <: Number = ^(RExpr(a),RExpr(s)) |> parse

*(a::T,s::S) where T <: MatOnly where S <: ESN = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: ESN where S <: MatOnly = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: MatOnly where S <: MatOnly = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: Vector where S <: Adjoint = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: Vector where S <: MatExpr = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: Adjoint where S <: Vector = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: Adjoint where S <: MatExpr = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: MatExpr where S <: Vector = *(RExpr(a),RExpr(s)) |> parse |> mat
*(a::T,s::S) where T <: MatExpr where S <: Adjoint = *(RExpr(a),RExpr(s)) |> parse |> mat

for o in [:+,:-]
    @eval begin
        function $o(a::T,s::S) where T <: MatOnly where S <: ESN
            $o(RExpr(a),RExpr(s)*RExpr(ones(size(a)))) |> parse |> mat
        end
        function $o(a::T,s::S) where T <: ESN where S <: MatOnly
            $o(RExpr(a)*RExpr(ones(size(s))),RExpr(s)) |> parse |> mat
        end
        $o(a::T,s::S) where T <: MatOnly where S <: MatOnly = $o(RExpr(a),RExpr(s)) |> parse |> mat
        $o(a::T,s::S) where T <: Vector where S <: Adjoint = $o(RExpr(a),RExpr(s)) |> parse |> mat
        $o(a::T,s::S) where T <: Vector where S <: MatOnly = $o(RExpr(a),RExpr(s)) |> parse |> mat
        function $o(a::T,s::ESN) where T <: Vector
            $o(RExpr(a),RExpr(s)*RExpr(ones(size(a)))) |> parse |> mat
        end
        $o(a::T,s::S) where T <: Adjoint where S <: Vector = $o(RExpr(a),RExpr(s)) |> parse |> mat
        $o(a::T,s::S) where T <: Adjoint where S <: MatOnly = $o(RExpr(a),RExpr(s)) |> parse |> mat
        function $o(a::T,s::ESN) where T <: Adjoint
            $o(RExpr(a),RExpr(s)*RExpr(ones(size(a)))) |> parse |> mat
        end
        $o(a::T,s::S) where T <: MatOnly where S <: Vector = $o(RExpr(a),RExpr(s)) |> parse |> mat
        function $o(a::ESN,s::S) where S <: Vector
            $o(RExpr(a)*RExpr(ones(size(s))),RExpr(s)) |> parse |> mat
        end
        $o(a::T,s::S) where T <: MatOnly where S <: Adjoint = $o(RExpr(a),RExpr(s)) |> parse |> mat
        function $o(a::ESN,s::S) where S <: Adjoint
            $o(RExpr(a)*RExpr(ones(size(s))),RExpr(s)) |> parse |> mat
        end
    end
end

for o in [:/,://]
    @eval begin
        $o(a::MatOnly,b::ESN) = $o(RExpr(a),RExpr(b)) |> parse |> mat
    end
end
//(expr,b::T) where T <: AbstractFloat = //(RExpr(expr),RExpr(b)) |> parse
//(a::T,expr) where T <: AbstractFloat = //(RExpr(a),RExpr(expr)) |> parse
//(expr::ExprSymbol,b::T) where T <: AbstractFloat = //(RExpr(expr),RExpr(b)) |> parse
//(a::T,expr::ExprSymbol) where T <: AbstractFloat = //(RExpr(a),RExpr(expr)) |> parse
function //(a::T,b::T) where T <: AbstractFloat
    isnan(a) | isnan(b) | (isinf(a) & isinf(b)) && return NaN
    return //(RExpr(a),RExpr(b)) |> parse |> eval
end

export inv, \

inv(r::RExpr) = AbstractTensors.:^(r,-1)
inv(r::T) where T <: MatExpr = AbstractTensors.:^(r,-1)

\(a,b) = Base.:\(a,b)
\(a::T,s::S) where T <: MatExpr where S <: Vector = (RExpr(a)^-1)*RExpr(s) |> parse |> mat
\(a::T,s::S) where T <: Vector where S <: MatExpr = (RExpr(a)^-1)*RExpr(s) |> parse |> mat

function solve(a::Array{T,1},s::Array{Symbol,1}) where T <: Any
    out = solve(list(a),list(s))
    return T <: String ? out.str : parse(out)
end
solve(a::T,s::Symbol) where T <: Vector = solve(a,[s])
solve(a::Expr,s::Array{Symbol,1}) = solve(a.head == :block ? a.args : [a],s)
solve(a::Expr,s::Symbol) = solve(a,[s])
solve(a::T,s::S) where T <: Tuple where S <: Tuple = solve(RExpr(a),RExpr(s)) |> parse
solve(a::T,s::Symbol) where T <: Tuple = solve(a,(s,))

odesolve(r...) = odesolve(RExpr.(r)...) |> parse

changevar(r...) = changevar(RExpr.(r)...) |> parse
map(a,b) = map(RExpr(a),RExpr(b)) |> parse
map(a::S,b::T) where S <: ExprSymbol where T <: MatOnly = map(RExpr(a),RExpr(b)) |> parse |> mat
select(a,b)  = select(RExpr(a),RExpr(b)) |> parse

order(::Nothing) = order(R"nil") |> parse
korder(::Nothing) = korder(R"nil") |> parse

export ∑, ∏, sub

const (∑, ∏) = (sum, prod)

@doc """
    sum(expr,k,lolim,uplim)

This implements the Gosper algorithm for the summation of series. The operator `sum` returns the indefinite or definite summation of a given expression, used with the syntax:
```
SUM(EXPR:expression, K:kernel, [LOLIM:expression [, UPLIM:expression]])
```
If there is no closed form solution, these operators return the input unchanged. `UPLIM` and `LOLIM` are optional parameters specifying the lower limit and upper limit of the summation. If `UPLIM` is not supplied, the upper limit is taken as `K` (the summation variable itself).

For example:
```Julia
Algebra.sum(:(n^3),:n)
Algebra.sum(:(a+k*r),:k,0,:(n-1))
Algebra.sum(:(1/((p+(k-1)*q)*(p+k*q))),:k,1,:(n+1))
```
""" Reduce.Algebra.sum

@doc """
    prod(expr,k,lolim,uplim)

The operator `prod` returns the product of the given expression, used with the syntax:
```
PROD(EXPR:expression, K:kernel, [LOLIM:expression [, UPLIM:expression]])
```
If there is no closed form solution, these operators return the input unchanged. `UPLIM` and `LOLIM` are optional parameters specifying the lower limit and upper limit of the product. If `UPLIM` is not supplied, the upper limit is taken as `K` (the product variable itself).

For example:
```Julia
Algebra.prod(:(k/(k-2)),:k)
```
""" Reduce.Algebra.prod

"""
    sub(::Union{Dict,Pair},expr)

Make variable substitutions using Reduce's native sub command. Syntax:
```
R"(⟨substitution_list⟩,⟨EXPRN1:algebraic⟩)"
```
where `⟨substitution_list⟩` is a list of one or more equations of the form
```
⟨VAR:kernel⟩ = ⟨EXPRN:algebraic⟩
```
or a kernel that evaluates to such a list.

The `sub` operator gives the algebraic result of replacing every occurrence of the variable `var` in the expression `EXPRN1` by the expression `EXPRN`. Specifically, `EXPRN1` is first evaluated using all available rules. Next the substitutions are made, and finally the substituted expression is reevaluated. When more than one variable occurs in the substitution list, the substitution is performed by recursively walking down the tree representing `EXPRN1`, and replacing every `VAR` found by the appropriate `EXPRN`. The `EXPRN` are not themselves searched for any occurrences of the various `VAR`s. The trivial case `sub`(EXPRN1)` returns the algebraic value of `EXPRN1`.

*Examples:*
```Julia
julia> Algebra.sub((:(x=a+y),:(y=y+1)),:(x^2+y^2))
:((a + y) ^ 2 + (y + 1) ^ 2)
```
and with `@rcall s = (x=a+y,y=y+1)`,
```
julia> Algebra.sub(:s,:(x^2+y^2))
:((a + y) ^ 2 + (y + 1) ^ 2)
```
Note that the global assignments `R"x:=a+y"`, etc., do not take place.

`EXPRN1` can be any valid algebraic expression whose type is such that a substitution process is defined for it (e.g., scalar expressions, lists and matrices). An error will occur if an expression of an invalid type for substitution occurs either in `EXPRN` or `EXPRN1`.
"""
sub(syme::String,expr::RExpr) = "sub($syme,$expr)" |> rcall |> RExpr
sub(syme::String,expr::T) where T = convert(T,sub(syme,RExpr(expr)))
sub(s::Dict{String,String},expr) = sub(Reduce._syme(s),expr)
sub(s::Dict{<:Any,<:Any},expr) = sub(Dict([=>(string.(RExpr.([b[1],b[2]]))...) for b ∈ collect(s)]...),expr)
sub(s::Pair{<:Any,<:Any},expr) = sub(Dict(s),expr)
sub(s::Array{<:Pair{<:Any,<:Any},1},expr) = sub(Dict(s...),expr)
sub(s::Symbol,expr) = sub(string(s),expr)
sub(s::T,expr) where T <: Tuple = sub(string(list(s)),expr)
sub(s::Expr,expr) = sub(string(list(s)),expr)

@doc """
    on(::Symbol...)

Takes a list of switch names as argument and turns them on.
""" Reduce.Algebra.on

@doc """
    off(::Symbol...)

Takes a list of switch names as argument and turns them off.
""" Reduce.Algebra.off

@doc """
    max(r...)

`max` can take an arbitrary number of expressions as their arguments. If all arguments evaluate to numerical values, the maximum of the argument list is returned. If any argument is non-numeric, an appropriately reduced expression is returned. For example:

```Julia
julia> Algebra.max(2,-3,4,5)
5

julia> Algebra.max(:a,2,3)
:(max(3, a))
```

`max` of an empty list returns 0.
""" Reduce.Algebra.max

@doc """
    min(r...)

`min` can take an arbitrary number of expressions as their arguments. If all arguments evaluate to numerical values, the minimum of the argument list is returned. If any argument is non-numeric, an appropriately reduced expression is returned. For example:

```Julia
julia> Algebra.min(2,-2)
-2

julia> Algebra.min(:x)
:x
```

`min` of an empty list returns 0.
""" Reduce.Algebra.min

@doc """
    df(r...)

The operator `df` is used to represent partial differentiation with respect to one or more variables. It is used with the syntax:
```Julia
R"df(⟨EXPRN:algebraic⟩[,⟨VAR:kernel⟩ <,⟨NUM:integer⟩ >])"
```
The first argument is the expression to be differentiated. The remaining arguments specify the differentiation variables and the number of times they are applied.

The number `num` may be omitted if it is `1`. For example,
```Julia
reduce> df(y,x)

reduce> df(y,x,2)

reduce> df(y,x1,2,x2,x3,2)
```
The evaluation of `df(y,x)` proceeds as follows: first, the values of `y` and `x` are found. Let us assume that `x` has no assigned value, so its value is `x`. Each term or other part of the value of `y` that contains the variable `x` is differentiated by the standard rules. If `z` is another variable, not `x` itself, then its derivative with respect to `x` is taken to be `0`, unless `z` has previously been declared to `depend` on `x`, in which case the derivative is reported as the symbol `df(z,x)`.
""" Reduce.Algebra.df

@doc """
    int(r...)

`int` is an operator in REDUCE for indefinite integration using a combination of the Risch-Norman algorithm and pattern matching. It is used with the syntax:
```Julia
R"int(⟨EXPRN:algebraic⟩,⟨VAR:kernel⟩)"
```
This will return correctly the indefinite integral for expressions comprising polynomials, log functions, exponential functions and tan and atan. The arbitrary constant is not represented. If the integral cannot be done in closed terms, it returns a formal integral for the answer in one of two ways:

1. It returns the input, `int(…,…)` unchanged.
2. It returns an expression involving `int`s of some other functions (sometimes more complicated than the original one, unfortunately).

Rational functions can be integrated when the denominator is factorizable by the program. In addition it will attempt to integrate expressions involving error functions, dilogarithms and other trigonometric expressions. In these cases it might not always succeed in finding the solution, even if one exists.

*Examples:*
```Julia
julia> Algebra.int(:(log(x)),:x)
:((log(x) - 1) * x)

julia> Algebra.int(:(e^x),:x)
:(e ^ x)
```
The program checks that the second argument is a variable and gives an error if it is not.

*Note:* If the `int` operator is called with 4 arguments, REDUCE will implicitly call the definite integration package (DEFINT) and this package will interpret the third and fourth arguments as the lower and upper limit of integration, respectively. For details, consult the documentation on the DEFINT package.
""" Reduce.Algebra.int

@doc """
    pochhammer(a,k)

The Pochhammer notation ``(a)_k`` (also called Pochhammer’s symbol) is supported by the binary operator `pochhammer(a,k)`. For a non-negative integer `k`, it is defined as ([http://dlmf.nist.gov/5.2.iii](http://dlmf.nist.gov/5.2.iii))

``(a)_0	= 1,``

``(a)_k	= a(a + 1)(a + 2)⋅⋅⋅(a + k - 1).``

For ``a ⁄= 0,±1,±2,…``, this is equivalent to

``(a)k	= \\frac{\\Gamma (a+-k-)}{\\Gamma (a)}``

With `rounded` off, this expression is evaluated numerically if `a` and `k` are both integral, and otherwise may be simplified where appropriate. The simplification rules are based upon algorithms supplied by Wolfram Koepf.
""" Reduce.Algebra.pochhammer

@doc """
    pf(expr,var)

`R"pf(⟨exp⟩,⟨var⟩)"` transforms the expression `⟨exp⟩` into a list of partial fractions with respect to the main variable, `⟨var⟩`. `pf` does a complete partial fraction decomposition, and as the algorithms used are fairly unsophisticated (factorization and the extended Euclidean algorithm), the code may be unacceptably slow in complicated cases.
""" Reduce.Algebra.pf

@doc """
    solve(r...)

`solve` is an operator for solving one or more simultaneous algebraic equations. It is used with the syntax:
```Julia
R"SOLVE(⟨EXPRN:algebraic⟩[,⟨VAR:kernel⟩∣,⟨VARLIST:list of kernels⟩])"
```
`exprn` is of the form `⟨expression⟩` or `{⟨expression1⟩,⟨expression2⟩, …}`. Each expression is an algebraic equation, or is the difference of the two sides of the equation. The second argument is either a kernel or a list of kernels representing the unknowns in the system. This argument may be omitted if the number of distinct, non-constant, top-level kernels equals the number of unknowns, in which case these kernels are presumed to be the unknowns.

For one equation, `solve` recursively uses factorization and decomposition, together with the known inverses of `log`, `sin`, `cos`, `^`, `acos`, `asin`, and linear, quadratic, cubic, quartic, or binomial factors. Solutions of equations built with exponentials or logarithms are often expressed in terms of Lambert’s `W` function. This function is (partially) implemented in the special functions package.

Linear equations are solved by the multi-step elimination method due to Bareiss, unless the switch `cramer` is on, in which case Cramer’s method is used. The Bareiss method is usually more efficient unless the system is large and dense.

Non-linear equations are solved using the Groebner basis package (chapter 16.28). Users should note that this can be quite a time consuming process.

*Examples:*
```Julia
Algebra.solve(:(log(sin(x+3))^5 == 8),:x)
Algebra.solve(:(a*log(sin(x+3))^5 - b), :(sin(x+3)))
Algebra.solve((:(a*x+y==3),:(y=-2)),(:x,:y))
```
`solve` returns a list of solutions. If there is one unknown, each solution is an equation for the unknown. If a complete solution was found, the unknown will appear by itself on the left-hand side of the equation. On the other hand, if the solve package could not find a solution, the “solution” will be an equation for the unknown in terms of the operator `root_of`. If there are several unknowns, each solution will be a list of equations for the unknowns. For example,
```Julia
julia> Algebra.solve(:(x^2==1),:x)
(:(x = 1), :(x = -1))

julia> Algebra.solve(:(x^7-x^6+x^2==1),:x)
(:(x = root_of(x_ ^ 6 + x_ + 1, x_, tag_1)), :(x = 1))

julia> Algebra.solve((:(x+3y==7),:(y-x==1)),(:x,:y))
(:(x = 1), :(y = 2))
```
The `tag` argument is used to uniquely identify those particular solutions.
```
""" Reduce.Algebra.solve

@doc """
    even(r...)

An operator can be declared to be even in its first argument by the declarations `even`. Expressions involving an operator declared in this manner are transformed if the first argument contains a minus sign. Any other arguments are not affected. For example, the declaration
```Julia
julia> Algebra.even(:f1)
```
means that
```
        f1(-a)    ->    f1(a)
        f1(-a,-b) ->    f1(a,-b)
```
""" Reduce.Algebra.even

@doc """
    odd(r...)

An operator can be declared to be odd in its first argument by the declarations `odd`. Expressions involving an operator declared in this manner are transformed if the first argument contains a minus sign. Any other arguments are not affected. In addition, if say `f` is declared odd, then `f(0)` is replaced by zero unless `f` is also declared non zero by the declaration `nonzero`. For example, the declarations
```Julia
julia> Algebra.odd(:f2)
```
means that
```
        f2(-a)    ->   -f2(a)
        f2(0)     ->    0
```
To inhibit the last transformation, say `nonzero(:f2)`.
""" Reduce.Algebra.odd

@doc """
    linear(r...)

An operator can be declared to be linear in its first argument over powers of its second argument. If an operator `f` is so declared, `f` of any sum is broken up into sums of `f`s, and any factors that are not powers of the variable are taken outside. This means that `f` must have (at least) two arguments. In addition, the second argument must be an identifier (or more generally a kernel), not an expression.

*Example:* If `f` were declared linear, then
```
f(a*x^5+b*x+c,x) ->  f(x^5,x)*a + f(x,x)*b + f(1,x)*c
```
More precisely, not only will the variable and its powers remain within the scope of the `f` operator, but so will any variable and its powers that had been declared to `depend` on the prescribed variable; and so would any expression that contains that variable or a dependent variable on any level, e.g. `cos(sin(x))`.

To declare operators `f` and `g` to be linear operators, use:
```Julia
julia> Algebra.linear(:f,:g)
```
The analysis is done of the first argument with respect to the second; any other arguments are ignored. It uses the following rules of evaluation:
```
f(0) 		-> 0
f(-y,x) 	-> -f(y,x)
f(y+z,x) 	-> f(y,x)+f(z,x)
f(y*z,x) 	-> z*f(y,x)   	if z does not depend on x
f(y/z,x) 	-> f(y,x)/z	if z does not depend on x
```
""" Reduce.Algebra.linear

@doc """
    noncom(r...)

An operator can be declared to be non-commutative under multiplication by the declaration `noncom`.

*Example:* After the declaration
```Julia
julia> Algebra.noncom(:u,:v);
```
the expressions `u(x)*u(y)-u(y)*u(x)` and `u(x)*v(y)-v(y)*u(x)` will remain unchanged on simplification, and in particular will not simplify to zero.

Note that it is the operator (`u` and `v` in the above example) and not the variable that has the non-commutative property.
""" Reduce.Algebra.noncom

@doc """
    symmetric(r...)

An operator can be declared to be symmetric with respect to its arguments by the declaration `symmetric`. For example
```Julia
julia> Algebra.symmetric(:u,:v);
```
means that any expression involving the top level operators `u` or `v` will have its arguments reordered to conform to the internal order used by REDUCE. The user can change this order for kernels by the command `korder`.
For example, `u(x,v(1,2))` would become `u(v(2,1),x)`, since numbers are ordered in decreasing order, and expressions are ordered in decreasing order of complexity.
""" Reduce.Algebra.symmetric

@doc """
    antisymmetric(r...)

the declaration `antisymmetric` declares an operator antisymmetric. For example,
```Julia
julia> Algebra.antisymmetric(:l,:m);
```
means that any expression involving the top level operators `l` or `m` will have its arguments reordered to conform to the internal order of the system, and the sign of the expression changed if there are an odd number of argument interchanges necessary to bring about the new order.

For example, `l(x,m(1,2))` would become `-l(-m(2,1),x)` since one interchange occurs with each operator. An expression like `l(x,x)` would also be replaced by `0`.
""" Reduce.Algebra.antisymmetric

@doc """
    operator(r...)

The user may add new prefix operators to the system by using the declaration `operator`. For example:
```Julia
julia> Algebra.operator(:h,:g1,:arctan)
```
adds the prefix operators `h`, `g1` and `arctan` to the system.

This allows symbols like `h(w)`, `h(x,y,z)`, `g1(p+q)`, `arctan(u/v)` to be used in expressions, but no meaning or properties of the operator are implied. The same operator symbol can be used equally well as a 0-, 1-, 2-, 3-, etc.-place operator.

To give a meaning to an operator symbol, or express some of its properties, `let` statements can be used, or the operator can be given a definition as a procedure.
""" Reduce.Algebra.operator

@doc """
    depend(r...)

There are several facilities in REDUCE, such as the differentiation operator and the linear operator facility, that can utilize knowledge of the dependency between various variables, or kernels. Such dependency may be expressed by the command `depend`. This takes an arbitrary number of arguments and sets up a dependency of the first argument on the remaining arguments. For example,
```Julia
julia> Algebra.depend(:x,:y,:z)
```
says that `x` is dependent on both `y` and `z`.
```Julia
julia> Algebra.depend(:z,:(cos(x)),:y)
```
says that `z` is dependent on `cos(x)` and `y`.
""" Reduce.Algebra.depend

@doc """
    nodepend(r...)

Dependencies introduced by `depend` can be removed by `nodepend`. The arguments of this are the same as for `depend`. For example, given the above dependencies,
```Julia
julia> Algebra.nodepend(:z,:(cos(x)))
```
says that `z` is no longer dependent on `cos(x)`, although it remains dependent on `y`.
""" Reduce.Algebra.nodepend

@doc """
    set(a,b)

In some cases, it is desirable to perform an assignment in which both the left- and right-hand sides of an assignment are evaluated. In this case, the `set` statement can be used with the syntax:
```Julia
R"set(⟨expression⟩,⟨expression⟩)"
```
For example, the statements
```
        j := 23;
        set(mkid(a,j),x);
```
assigns the value `x` to `a23`.
""" Reduce.Algebra.set

@doc """
    unset(r)

To remove a value from such a variable, the `unset` statement can be used with the syntax:
```
R"unset(⟨expression⟩)"
```
For example, the statement
```
        j := 23;
        unset(mkid(a,j));
```
clears the value of `a23`.
""" Reduce.Algebra.unset

@doc """
    fibonacci(n,x)

Fibonacci Polynomials are computed by the binary operator `fibonaccip`. `fibonaccip(n,x)` returns the `n`th Fibonacci polynomial in the variable `x`. If `n` is a positive or negative integer, it will be evaluated following the definition:

\$F_0(x) = 0; F_1(x) = 1; F_n(x) = xF_{n-1}(x) + F_{n-2}(x)\$
""" Reduce.Algebra.fibonaccip

@doc """
    mkid(u,v)

In many applications, it is useful to create a set of identifiers for naming objects in a consistent manner. In most cases, it is sufficient to create such names from two components. The operator `mkid` is provided for this purpose. Its syntax is:
```Julia
R"mkid(U:id,V:id|non-negative integer)"
```
for example
```Julia
julia> Algebra.mkid(:a,3)
:a3

julia> Algebra.mkid(:apple,:s)
:apples
```
while `mkid(:(a+b),2)` gives an error.
""" Reduce.Algebra.mkid

@doc """
    infix(r...)

Users can add new infix operators by using the declarations `infix` and `precedence`. For example,
```Julia
julia> Algebra.infix(:mm)
```
The declaration `infix(:mm)` would allow one to use the symbol `mm` as an infix operator:
`R"a mm b"` instead of `R"mm(a,b)"`.
""" Reduce.Algebra.infix

@doc """
    precedence(a,b)

Users can add new infix operators by using the declarations `infix` and `precedence`. For example,
```Julia
julia> Algebra.precedence(:mm,:-)
```
The declaration `precedence(:mm,:-)` says that `mm` should be inserted into the infix operator precedence list just after the `-` operator. This gives it higher precedence than `-` and lower precedence than `*` . Thus `R"a - b mm c - d"` means `R"a - (b mm c) - d"`, while `R"a * b mm c * d"` means `R"(a * b) mm (c * d)"`.
""" Reduce.Algebra.precedence

@doc """
    order(r...)

The declaration `order` may be used to order variables on output. The syntax is:
```Julia
julia> Algebra.order(v1,...vn)
```
where the `vi` are kernels. Thus,
```Julia
julia> Algebra.order(:x,:y,:z)
```
orders `x` ahead of `y`, `y` ahead of `z` and all three ahead of other variables not given an order. `order(nothing)` resets the output order to the system default. The order of variables may be changed by further calls of `order`, but then the reordered variables would have an order lower than those in earlier `order` calls. Thus,
```Julia
julia> Algebra.order(:x,:y,:z)

julia> Algebra.order(:y,:x)
```
would order `z` ahead of `y` and `x`. The default ordering is usually alphabetic.
""" Reduce.Algebra.order

@doc """
    factor(r...)

This declaration takes a list of identifiers or kernels as argument. `factor` is not a factoring command (use `factorize` or the `factor` switch for this purpose); rather it is a separation command. All terms involving fixed powers of the declared expressions are printed as a product of the fixed powers and a sum of the rest of the terms.

For example, after the declaration
```Julia
julia> Algebra.factor(:x)
```
the polynomial \$(x + y + 1)^2\$ will be printed as
```
         2                  2
        x  + 2*x*(y + 1) + y  + 2*y + 1
```
All expressions involving a given prefix operator may also be factored by putting the operator name in the list of factored identifiers. For example:
```Julia
julia> Algebra.factor(:x,:cos,:(sin(x))
```
causes all powers of `x` and `sin(x)` and all functions of `cos` to be factored.
""" Reduce.Algebra.factor

@doc """
    remfac(r...)

The declaration `remfac(v1,...,vn)` removes the factoring flag from the expressions `v1` through `vn`.
""" Reduce.Algebra.remfac

@doc """
    korder(r...)

The internal ordering of variables (more specifically kernels) can have a significant effect on the space and time associated with a calculation. In its default state, REDUCE uses a specific order for this which may vary between sessions. However, it is possible for the user to change this internal order by means of the declaration `korder`. The syntax for this is:
```Julia
julia> Algebra.korder(v1,...,vn)
```
where the `vi` are kernels. With this declaration, the `vi` are ordered internally ahead of any other kernels in the system. `v1` has the highest order, `v2` the next highest, and so on. A further call of `korder` replaces a previous one. `korder(nothing)` resets the internal order to the system default.
""" Reduce.Algebra.korder

@doc """
    realvalued(r...)

The declaration `realvalued` may be used to restrict variables to the real numbers. The syntax is:
```Julia
	Algebra.realvalued(v1,...vn)
```
For such variables the operator `impart` gives the result zero. Thus, with
```Julia
julia> Algebra.realvalued(:x,:y)
```
the expression `impart(x+sin(y))` is evaluated as zero. You may also declare an operator as real valued with the meaning, that this operator maps real arguments always to real values. Example:
```Julia
julia> Algebra.operator(:h); Algebra.realvalued(:h,:x)

julia> Algebra.impart(:(h(x)))
0

julia> Algebra.impart(:(h(w)))
:(impart(h(w)))
```
Such declarations are not needed for the standard elementary functions.
""" Reduce.Algebra.realvalued

@doc """
    notrealvalued(r...)

To remove the `realvalued` propery from a variable or an operator use the declaration `notrealvalued` with the syntax:
```Julia
julia> Algebra.notrealvalued(v1,...vn)
```
""" Reduce.Algebra.notrealvalued

@doc """
    factorize(r...)

It is also possible to factorize a given expression explicitly. The operator `factorize` that invokes this facility is used with the syntax
```Julia
R"factorize(EXPRN:polynomial[,INTEXP:prime integer])"
```
the optional argument of which will be described later. Thus to find and display all factors of the cyclotomic polynomial ``x^{105} - 1``, one could write:
```
julia> Algebra.factorize(:(x^105-1))
```
The result is a list of factor,exponent pairs. In the above example, there is no overall numerical factor in the result, so the results will consist only of polynomials in x. The number of such polynomials can be found by using the operator `length`. If there is a numerical factor, as in factorizing ``12x^2 - 12``, that factor will appear as the first member of the result. It will however not be factored further. Prime factors of such numbers can be found, using a probabilistic algorithm, by turning on the switch `ifactor`. For example,
```
julia> Algebra.on(:ifactor); Algebra.factorize(:(12x^2-12))
```
would result in the output
```
((2, 2), (3, 1), (:(x ^ 2 + 1), 1), (:(x + 1), 1), (:(x - 1), 1))
```
If the first argument of `factorize` is an integer, it will be decomposed into its prime components, whether or not `ifactor` is on.

Note that the `ifactor` switch only affects the result of `factorize`. It has no effect if the `factor` switch is also on.
""" Reduce.Algebra.factorize

@doc """
    remainder(a,b)

This operator is used with the syntax
```Julia
R"REMAINDER(EXPRN1:polynomial,EXPRN2:polynomial)"
```
It returns the remainder when `EXPRN1` is divided by `EXPRN2`. This is the true remainder based on the internal ordering of the variables, and not the pseudo-remainder. The pseudo-remainder and in general pseudo-division of polynomials can be calculated after loading the `polydiv` package. Please refer to the documentation of this package for details.

*Examples:*
```Julia
julia> Algebra.on(:exp); Algebra.remainder(:((x+y)*(x+2*y)),:(x+3*y))
:(2 * y ^ 2)

julia> Algebra.remainder(:(2*x+y),2)
:y
```
""" Reduce.Algebra.remainder

@doc """
    resultant(a,b,var)

This is used with the syntax
```Julia
R"resultant(EXPRN1:polynomial,EXPRN2:polynomial,VAR:kernel)"
```
It computes the resultant of the two given polynomials with respect to the given variable, the coefficients of the polynomials can be taken from any domain. The result can be identified as the determinant of a Sylvester matrix, but can often also be thought of informally as the result obtained when the given variable is eliminated between the two input polynomials. If the two input polynomials have a non-trivial GCD their resultant vanishes.

The switch `bezout` controls the computation of the resultants. It is off by default. In this case a subresultant algorithm is used. If the switch Bezout is turned on, the resultant is computed via the Bezout Matrix. However, in the latter case, only polynomial coefficients are permitted.
""" Reduce.Algebra.resultant

@doc """
    interpol(val,var,mp)

Syntax:
```Julia
R"interpol(⟨values⟩,⟨variable⟩,metapoints)"
```
where `⟨values⟩` and `⟨points⟩` are lists of equal length and `<variable>` is an algebraic expression (preferably a kernel).

`interpol` generates an interpolation polynomial ``f`` in the given variable of degree `length(⟨values⟩)-1`. The unique polynomial ``f`` is defined by the property that for corresponding elements ``v`` of `⟨values⟩` and ``p`` of `⟨points⟩` the relation ``f(p) = v`` holds.

The Aitken-Neville interpolation algorithm is used which guarantees a stable result even with rounded numbers and an ill-conditioned problem.
""" Reduce.Algebra.interpol

@doc """
    deg(p,var)

This operator is used with the syntax
```Julia
R"deg(EXPRN:polynomial,VAR:kernel)"
```
It returns the leading degree of the polynomial `EXPRN` in the variable `VAR`. If `VAR` does not occur as a variable in `EXPRN`, 0 is returned.

*Examples:*
```Julia
julia> Algebra.on(:exp)

julia> Algebra.deg(:((a+b)*(c+2*d)^2),:a)
1

julia> Algebra.deg(:((a+b)*(c+2*d)^2),:d)
2

julia> Algebra.deg(:((a+b)*(c+2*d)^2),:e)
0
```
Note also that if `ratarg` is on,
```
        deg((a+b)^3/a,a)       ->  3
```
since in this case, the denominator `a` is considered part of the coefficients of the numerator in `a`. With `ratarg` off, however, an error would result in this case.
""" Reduce.Algebra.deg

@doc """
    lcof(expr,var)

`lcof` is used with the syntax
```Julia
R"lcof(EXPRN:polynomial,VAR:kernel)"
```
It returns the leading coefficient of the polynomial `EXPRN` in the variable `VAR`. If `VAR` does not occur as a variable in `EXPRN`, `EXPRN` is returned.

*Examples:*
```Julia
julia> Algebra.on(:exp)

julia> Algebra.lcof(:((a+b)*(c+2*d)^2),:a)
:(c ^ 2 + 4 * c * d + 4 * d ^ 2)

julia> Algebra.lcof(:((a+b)*(c+2*d)^2),:d)
:(4 * (a + b))

julia> Algebra.lcof(:((a+b)*(c+2*d)),:e)
:(a * c + 2 * a * d + b * c + 2 * b * d)
```
""" Reduce.Algebra.lcof

@doc """
    lpower(exprn,var)

Syntax:
```Julia
R"lpower(EXPRN:polynomial,VAR:kernel)"
```
`lpower` returns the leading power of `EXPRN` with respect to `VAR`. If `EXPRN` does not depend on `VAR`, 1 is returned.
*Examples:*
```Julia
julia> Algebra.on(:exp)

julia> Algebra.lpower(:((a+b)*(c+2*d)^2),:a)
:a

julia> Algebra.lpower(:((a+b)*(c+2*d)^2),:d)
:(d ^ 2)

julia> Algebra.lpower(:((a+b)*(c+2*d)),:e)
1
```
""" Reduce.Algebra.lpower

@doc """
    lterm(exprn,var)

Syntax:
```Julia
R"lterm(EXPRN:polynomial,VAR:kernel)"
```
`lterm` returns the leading term of `EXPRN` with respect to `VAR`. If `EXPRN` does not depend on `VAR`, `EXPRN` is returned.

*Examples:*
```Julia
julia> Algebra.on(:exp)

julia> Algebra.lterm(:((a+b)*(c+2*d)^2),:a)
:(a * (c ^ 2 + 4 * c * d + 4 * d ^ 2))

julia> Algebra.lterm(:((a+b)*(c+2*d)^2),:d)
:(4 * d ^ 2 * (a + b))

julia> Algebra.lterm(:((a+b)*(c+2*d)),:e)
:(a * c + 2 * a * d + b * c + 2 * b * d)
```
""" Reduce.Algebra.lterm

@doc """
    reduct(exprn,var)

Syntax:
```Julia
R"reduct(EXPRN:polynomial,VAR:kernel)"
```
Returns the reductum of `EXPRN` with respect to `VAR` (i.e., the part of `EXPRN` left after the leading term is removed). If `EXPRN` does not depend on the variable `VAR`, 0 is returned.

*Examples:*
```Julia
julia> Algebra.on(:exp)

julia> Algebra.reduct(:((a+b)*(c+2*d)),:a)
:(b * (c + 2d))

julia> Algebra.reduct(:((a+b)*(c+2*d)),:d)
:(c * (a + b))

julia> Algebra.reduct(:((a+b)*(c+2*d)),:e)
0
```
""" Reduce.Algebra.reduct

@doc """
    totaldeg(expr,var)

Syntax:
```Julia
julia> Algebra.totaldeg(:(a*x^2+b*x+c), :x)
2

julia> Algebra.totaldeg(:(a*x^2+b*x+c), (:a,:b,:c))
1

julia> Algebra.totaldeg(:(a*x^2+b*x+c), (:x, :a))
3

julia> Algebra.totaldeg(:(a*x^2+b*x+c), (:x,:b))
2

julia> Algebra.totaldeg(:(a*x^2+b*x+c), (:p,:q,:r))
0
```
`totaldeg(u, kernlist)` finds the total degree of the polynomial `u` in the variables in `kernlist`. If `kernlist` is not a list it is treated as a simple single variable. The denominator of `u` is ignored, and "degree" here does not pay attention to fractional powers. Mentions of a kernel within the argument to any operator or function (eg `sin`, `cos`, `log`, `sqrt`) are ignored. Really `u` is expected to be just a polynomial.
""" Reduce.Algebra.totaldeg

@doc """
    clear(r...)

The user may remove all assignments and substitution rules from any expression by the command `clear`, in the form
```
R"cilear ⟨expression⟩,…,⟨expression⟩ = ⟨terminator⟩"
```
e.g.
```Julia
julia> Algebra.clear(:x,:(h(x,y)))
```
Because of their *instant evaluation* property, array and matrix elements cannot be cleared with `clear`. For example, if `a` is an array, you must say
```Julia
R"a(3) := 0"
```
rather than
```Julia
R"clear a(3)"
```
to “clear” element `a(3)`.

On the other hand, a whole array (or matrix) `a` can be cleared by the command `clear(:a)`. This means much more than resetting to 0 all the elements of `a`. The fact that `a` is an array, and what its dimensions are, are forgotten, so `a` can be redefined as another type of object, for example an operator.

If you need to clear a variable whose name must be computed, see the `unset` statement.
""" Reduce.Algebra.clear

"""
    mateign(exprn,id)

Syntax:
```
        MATEIGEN(EXPRN:matrix_expression,ID):list.
```
`mateigen` calculates the eigenvalue equation and the corresponding eigenvectors of a matrix, using the variable `id` to denote the eigenvalue. A square free decomposition of the characteristic polynomial is carried out. The result is a list of lists of 3 elements, where the first element is a square free factor of the characteristic polynomial, the second its multiplicity and the third the corresponding eigenvector (as an *n* by *1* matrix). If the square free decomposition was successful, the product of the first elements in the lists is the minimal polynomial. In the case of degeneracy, several eigenvectors can exist for the same eigenvalue, which manifests itself in the appearance of more than one arbitrary variable in the eigenvector. To extract the various parts of the result use the operations defined on lists.

*Example:* The command
```
        mateigen(mat((2,-1,1),(0,1,1),(-1,1,1)),eta);
```
gives the output
```
        {{ETA - 1,2,

          [ARBCOMPLEX(1)]
          [             ]
          [ARBCOMPLEX(1)]
          [             ]
          [      0      ]

          },

         {ETA - 2,1,

          [      0      ]
          [             ]
          [ARBCOMPLEX(2)]
          [             ]
          [ARBCOMPLEX(2)]

          }}
```
""" Reduce.Algebra.mateign

@doc """
    cofactor(exprn,row,column)

Syntax:
```
  cofactor(EXPRN:matrix_expression,ROW:integer,COLUMN:integer):algebraic
```
The operator `cofactor` returns the cofactor of the element in row `ROW` and column `COLUMN` of the matrix `EXPRN`. Errors occur if `ROW` or `COLUMN` do not simplify to integer expressions or if `EXPRN` is not square.
""" Reduce.Algebra.cofactor

@doc """
    limit(exprn,var,limpint)

Syntax
```
LIMIT(⟨EXPRN:algebraic⟩,⟨VAR:kernel⟩,⟨LIMPOINT:algebraic⟩) : algebraic
```
This is the standard way of calling limit, applying all of the methods. The result is the limit of `EXPRN` as `VAR` approaches `LIMPOINT`.
""" Reduce.Algebra.limit

@doc """
    compact(exprn,list)

COMPACT is a package of functions for the reduction of a polynomial in the presence of side relations. `compact` applies the side relations to the polynomial so that an equivalent expression results with as few terms as possible. For example, the evaluation of
```Julia
Algebra.compact(:(s*(1-sin(x^2))+c*(1-cos(x^2))+sin(x^2)+cos(x^2)),
    (:(cos(x^2)+sin(x^2)=1),))
```
yields the result
```Julia
:(sin(x ^ 2) * c + 1 + cos(x ^ 2) * s)
```
The switch `trcompact` can be used to trace the operation.
""" Reduce.Algebra.compact

@doc """
    define(r...)

The command `define` allows a user to supply a new name for any identifier or replace it by any well-formed expression. Its argument is a list of expressions of the form
```
⟨identifier⟩ = 	⟨number⟩∣⟨identifier⟩∣⟨operator⟩∣
		⟨reserved word⟩∣⟨expression⟩
```
*Example:*
```Julia
Algebra.define(:(x==y+z))
```
""" Reduce.Algebra.define

@doc """
    changevar(depvars,newvars,eqlist,diffeq)

The operator `changevar` does a variable transformation in a set of differential equations. Syntax:
```
changevar(⟨depvars⟩,⟨newvars⟩,⟨eqlist⟩,⟨diffeq⟩)
```
`⟨diffeq⟩` is either a single differential equation or a list of differential equations, `⟨depvars⟩` are the dependent variables to be substituted, `⟨newvars⟩` are the new depend variables, and `⟨eqlist⟩` is a list of equations of the form `⟨depvar⟩=⟨expression⟩` where `⟨expression⟩` is some function in the new dependent variables.

The three lists `⟨depvars⟩`, `⟨newvars⟩`, and `⟨eqlist⟩` must be of the same length. If there is only one variable to be substituted, then it can be given instead of the list. The same applies to the list of differential equations, i.e., the following two commands are equivalent
```Julia
Algebra.operator(:u)
Algebra.changevar(:u,:y,:(x==e^y),:(df(u(x),x) - log(x)))
Algebra.changevar((:u,),(:y,),(:(x=e^y),),(:(df(u(x),x) - log(x)),))
```
except for one difference: the first command returns the transformed differential equation, the second one a list with a single element.
""" Reduce.Algebra.changevar

@doc """
    map(fnc,obj)

The `map` operator applies a uniform evaluation pattern to all members of a composite structure: a matrix, a list, or the arguments of an operator expression. The evaluation pattern can be a unary procedure, an operator, or an algebraic expression with one free variable.

It is used with the syntax:
```
   map(FNC:function,OBJ:object)
```
Here `OBJ` is a list, a matrix or an operator expression. `FNC` can be one of the following:

1. the name of an operator with a single argument: the operator is evaluated once with each element of `OBJ` as its single argument;
2. an algebraic expression with exactly one free variable, i.e. a variable preceded by the tilde symbol. The expression is evaluated for each element of `OBJ`, with the element substituted for the free variable;
3. a replacement rule of the form `var => rep` where `var` is a variable (a kernel without a subscript) and `rep` is an expression that contains `var`. The replacement expression `rep` is evaluated for each element of `OBJ` with the element substituted for `var`. The variable `var` may be optionally preceded by a tilde.

The rule form for `FNC` is needed when more than one free variable occurs.

*Examples:*
```Julia
julia> Algebra.map(:abs, (1,-2,:a,:(-a)))
(1, 2, :(abs(a)), :(abs(a)))

julia> Algebra.map(:(int(~w,x)), [:(x^2) :(x^5); :(x^4) :(x^5)])
2×2 Array{Any,2}:
 :(x ^ 3 // 3)  :(x ^ 6 // 6)
 :(x ^ 5 // 5)  :(x ^ 6 // 6)

 julia> Algebra.map(:(~w*6), :(x^2/3 == y^3/2 -1))
:(2 * x ^ 2 = 3 * (y ^ 3 - 2))
```
You can use `map` in nested expressions. However, you cannot apply `map` to a non-composite object, e.g. an identifier or a number.
""" Reduce.Algebra.map

@doc """
    select(fnc,lst)

The `select` operator extracts from a list, or from the arguments of an n–ary operator, elements corresponding to a boolean predicate. It is used with the syntax:
```
select(⟨FNC:function⟩,⟨LST:list⟩)
```
`FNC` can be one of the following forms:

1. the name of an operator with a single argument: the operator is evaluated once on each element of `LST`;
2. an algebraic expression with exactly one free variable, i.e. a variable preceded by the tilde symbol. The expression is evaluated for each element of `⟨LST⟩`, with the element substituted for the free variable;
3. a replacement rule of the form `⟨var⟩ => ⟨rep⟩` where `⟨var⟩` is a variable (a kernel without subscript) and `⟨rep⟩` is an expression that contains `⟨var⟩`. `⟨rep⟩` is evaluated for each element of `LST` with the element substituted for `⟨var⟩`. `⟨var⟩` may be optionally preceded by a tilde.

The rule form for `FNC` is needed when more than one free variable occurs.

The result of evaluating `FNC` is interpreted as a boolean value corresponding to the conventions of REDUCE. These values are composed with the leading operator of the input expression.

*Examples:*
```Julia
julia> Algebra.select(:(~w>0), (1,-1,2,-3,3))
(1, 2, 3)
```
```
    select(evenp deg(~w,y),part((x+y)^5,0):=list)
           -> {X^5 ,10*X^3*Y^2 ,5*X*Y^4}
    select(evenp deg(~w,x),2x^2+3x^3+4x^4) -> 4X^4 + 2X^2
```
""" Reduce.Algebra.select

@doc """
    g(id,exprn...)

Syntax:
```
        G(ID:identifier[,EXPRN:vector_expression])
                :gamma_matrix_expression.
```
`g` is an n-ary operator used to denote a product of γ matrices contracted with Lorentz four-vectors. Gamma matrices are associated with fermion lines in a Feynman diagram. If more than one such line occurs, then a different set of γ matrices (operating in independent spin spaces) is required to represent each line. To facilitate this, the first argument of `g` is a line identification identifier (not a number) used to distinguish different lines.
""" Reduce.Algebra.g

@doc """
    eps(exprn1,...,exprn4)

Syntax:
```
         EPS(EXPRN1:vector_expression,...,EXPRN4:vector_exp):vector_exp.
```
The operator `eps` has four arguments, and is used only to denote the completely antisymmetric tensor of order 4 and its contraction with Lorentz four-vectors. Thus

\$ϵ_{ijkl} = \\begin{cases} +1 & \\text{if }i,j,k,l\\text{ is an even permutation of 0,1,2,3} \\\\ - 1  & \\text{if }i,j,k,l\\text{ is an odd permutation of 0,1,2,3} \\\\ 0 & \\text{otherwise} \\end{cases}.\$

A contraction of the form \$ϵ_{ijμν}p_μq_ν\$ may be written as `eps(i,j,p,q)`, with `i` and `j` flagged as indices, and so on.
""" Reduce.Algebra.eps


text:
module GcodeParser

"""
    stripComments(line::String)::String

Return a copy of string `line` with stripped comments inside parentheses and all characters after a semicolon.

This function also removes whitespace as it it not needed for further parsing.

# Examples
```julia-repl
julia> stripComments("G92 (G10(aaa)))) ((comment)G) Z0.2 ; this is a comment")
"G92Z0.2"
```
"""
function stripComments(line::String)::String
    re1 = r"\(.*\)";    # Remove anything inside the outer parentheses
    re2 = r"[^:]\;.*";  # Remove anything after a semi-colon to the end of the line, including preceding spaces

    line = replace(line,  re1 => s"");
    line = replace(line,  re2 => s"");
    line = filter(x -> !isspace(x), line) # Remove whitespace

    return line;
end

"""
    parseLine(line::String, returnPair::Bool = true)::Array{Union{String,Pair{String,String}},1}

Parse a single line of g-code and return an array of `Pair{String,String}` or an array of `String` containing the parsed commands.

The first command usually defines what to do (ie. `G01` - linear interpolation) and following commands are the arguments (ie. `X 14.312`);

# Examples
```julia-repl
julia> parseLine("G10 X5.Y3. E6.")
4-element Array{Union{Pair{String,String}, String},1}:
 "G" => "10"
 "X" => "5."
 "Y" => "3."
 "E" => "6."
```

Return array of strings
```julia-repl
julia> parseLine("G10 X5.Y3. E6.", false)
4-element Array{Union{Pair{String,String}, String},1}:
 "G10"
 "X5."
 "Y3."
 "E6."
```
"""
function parseLine(line::String, returnPair::Bool = true)::Array{Union{String,Pair{String,String}},1}
    line = stripComments(line);

    # Match commands
    gcode_regex = r"/(%.*)|({.*)|((?:\$\$)|(?:\$[a-zA-Z0-9#]*))|([a-zA-Z][0-9\+\-\.]+)|(\*[0-9]+)/igm";

    # array of matched strings
    matches = collect(String(m.match) for m in eachmatch(gcode_regex, line));

    if returnPair
        return collect(first(m, 1) => last(m, length(m) - 1) for m in matches);
    end

    return matches;
end

function parseFile(path::String, callbacks::Dict{String,Function}, dataObject)
    open(path) do f
        line = 1
        while !eof(f)
            x = readline(f);
            cmds = parseLine(x);
            if length(cmds) == 0 continue; end;

            letter = cmds[1].first;
            number = cmds[1].second;

            command = "$letter$number";
            if haskey(callbacks, command)
                if dataObject === nothing
                    callbacks[command](cmds);
                else
                    callbacks[command](cmds, dataObject);
                end
            end

            line += 1;
        end
    end
end

export parseLine;
export parseFile;

end # module

ext: jl


text:
# Autogenerated wrapper script for libdeflate_jll for x86_64-apple-darwin
export libdeflate

JLLWrappers.@generate_wrapper_header("libdeflate")
JLLWrappers.@declare_library_product(libdeflate, "@rpath/libdeflate.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libdeflate,
        "lib/libdeflate.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type Entry{T} end

mutable struct DiagonalEntry{T,N,N²} <: Entry{T}
    D::SMatrix{N,N,T,N²}
    Dinv::SMatrix{N,N,T,N²}
    Δs::SVector{N,T}

    function DiagonalEntry{T,N}() where {T,N}
        N² = N^2
        D = @SMatrix zeros(T, N, N)
        Dinv = @SMatrix zeros(T, N, N)
        Δs = @SVector zeros(T, N)

        new{T,N,N²}(D, Dinv, Δs)
    end
end

mutable struct OffDiagonalEntry{T,N1,N2,N1N2} <: Entry{T}
    L::SMatrix{N2,N1,T,N1N2}
    U::SMatrix{N1,N2,T,N1N2}

    function OffDiagonalEntry{T,N1,N2}() where {T,N1,N2}
        L = @SMatrix zeros(T, N2, N1)
        U = @SMatrix zeros(T, N1, N2)

        new{T,N1,N2,N1 * N2}(L, U)
    end
end

mutable struct InequalityEntry{T,N} <: Entry{T}
    Δs::SVector{N,T}
    Δγ::SVector{N,T}

    function InequalityEntry{T,N}() where {T,N}
        Δs = @SVector zeros(T, N)
        Δγ = @SVector zeros(T, N)

        new{T,N}(Δs, Δγ)
    end
end

struct SparseLDU{T}
    diagonals::UnitDict{Base.OneTo{Int64},DiagonalEntry{T}}
    offdiagonals::Dict{Tuple{Int64,Int64},OffDiagonalEntry{T}}
    inequalities::UnitDict{UnitRange{Int64},InequalityEntry{T}}

    function SparseLDU(graph::Graph{N},bodies::Vector{Body{T}},eqcs::Vector{<:EqualityConstraint{T}},
        ineqcs::Vector{<:InequalityConstraint{T}},bdict::Dict,eqdict::Dict,ineqdict::Dict) where {T,N}

        diagonals = DiagonalEntry{T}[]
        for body in bodies
            push!(diagonals, DiagonalEntry{T,length(body)}())
        end
        for constraint in eqcs
            push!(diagonals, DiagonalEntry{T,length(constraint)}())
        end
        diagonals = UnitDict(diagonals)

        offdiagonals = Dict{Tuple{Int64,Int64},OffDiagonalEntry{T}}()
        for id in graph.dfslist
            haskey(bdict, id) ? node = bodies[bdict[id]] : node = eqcs[eqdict[id]]
            N1 = length(node)

            for cid in successors(graph, id)
                haskey(bdict, cid) ? cnode = bodies[bdict[cid]] : cnode = eqcs[eqdict[cid]]
                N2 = length(cnode)

                offdiagonals[(id, cid)] = OffDiagonalEntry{T,N2,N1}()
            end
        end

        inequalities = InequalityEntry{T}[]
        for ineq in ineqcs
            push!(inequalities, InequalityEntry{T,length(ineq)}())
        end
        if !isempty(inequalities)
            inequalities = UnitDict((ineqcs[1].id):(ineqcs[length(inequalities)].id), inequalities)
        else
            inequalities = UnitDict(0:0, inequalities)
        end

        new{T}(diagonals, offdiagonals, inequalities)
    end
end

@inline getentry(ldu::SparseLDU, id::Int64) = ldu.diagonals[id]
@inline getentry(ldu::SparseLDU, ids::Tuple{Int64,Int64}) = ldu.offdiagonals[ids]
@inline getineqentry(ldu::SparseLDU, id::Int64) = ldu.inequalities[id]


text:
using DiffEqBase, Catalyst, OrdinaryDiffEq, Sundials
using ReactionNetworkImporters
using TimerOutputs

# parameters
networkname = "testbcrbng"
tf          = 10000.
datadir     = joinpath(@__DIR__,"../data/bcr")
fname       = joinpath(datadir, "bcr.net")
build_jac   = true
sparse_jac  = false

const to = TimerOutput()
reset_timer!(to)

# get the reaction network
@timeit to "netgen" prnbng = loadrxnetwork(BNGNetwork(),string(networkname,"bng"), fname);
rn = prnbng.rn; u0 = prnbng.u₀; p = prnbng.p;
@timeit to "addodes" addodes!(rn; build_jac=build_jac, sparse_jac=sparse_jac, build_symfuncs=false, build_paramjac=false)
@timeit to "ODEProb" oprob = ODEProblem(rn, u0, (0.,tf), p)
u = copy(u0);
du = similar(u);
@timeit to "f1" rn.f(du,u,p,0.)
@timeit to "f2" rn.f(du,u,p,0.)
if build_jac
    J = zeros(length(u),length(u))
    @timeit to "J1" rn.jac(J,u,p,0.)
    @timeit to "J2" rn.jac(J,u,p,0.)
end
show(to)
println()

# note solvers run faster the second time

# Rodas4
#reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false,calck=false); end; show(to)
#reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false,calck=false); end; show(to)
#reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false); end;

# similarly CVODE_BDF with gmres
dense = true
calck = true
# @timeit to "CVODE_BDF-GMRES-1" begin sol = solve(oprob,CVODE_BDF(linear_solver=:GMRES),dense=dense, abstol=1e-8, reltol=1e-8); end;
# @timeit to "CVODE_BDF-GMRES-2" begin sol = solve(oprob,CVODE_BDF(linear_solver=:GMRES),dense=dense, abstol=1e-8, reltol=1e-8); end;

# CVODE_BDF with LU
@timeit to "CVODE_BDF-1" begin sol = solve(oprob,CVODE_BDF(),dense=dense, abstol=1e-8, reltol=1e-8); end; show(to)
@timeit to "CVODE_BDF-2" begin sol = solve(oprob,CVODE_BDF(),dense=dense, abstol=1e-8, reltol=1e-8); end; show(to)
@timeit to "KenCarp4-1" begin sol = solve(oprob,KenCarp4(autodiff=false),dense=dense, calck=calck, abstol=1e-8, reltol=1e-8); end; show(to)
@timeit to "KenCarp4-2" begin sol = solve(oprob,KenCarp4(autodiff=false),dense=dense, calck=calck, abstol=1e-8, reltol=1e-8); end;

show(to)


text:
include("exact_slater.jl")

function test_hydrogenic_slater_integrals!(data, tests,
                                           eq, term::Symbol, expected, rtol;
                                           print_res=false)
    eng = energy(eq, term)
    expectedf = convert(Float64, expected)
    δ = eng-expectedf
    test = :(isapprox($eng, $expected; rtol=$rtol))
    pass = @eval $test
    push!(data, [eq.orbital Dict(:onebody => "⟨h⟩", :direct => "⟨J⟩", :exchange => "⟨K⟩")[term] expected expectedf eng δ δ/abs(1e-10+abs(expectedf)) pass])
    push!(tests, test)
end

function test_hydrogenic_slater_integrals(fun::Function, fock::Fock, do_test::Bool=true)
    data = []
    tests = []

    for (i,eq) in enumerate(fock.equations.equations)
        cases = fun(i)
        for case in cases
            test_hydrogenic_slater_integrals!(data, tests, eq, case...)
        end
    end

    println("Orbital energies pre-optimization:")
    pass = Highlighter((data,i,j) -> j == 8 && data[i,j], bold = true, foreground = :green)
    fail = Highlighter((data,i,j) -> j == 8 && !data[i,j], bold = true, foreground = :red)
    pretty_table(vcat(data...), ["Orbital", "Term", "Expected", "Expected", "Actual", "Error", "Relative error", "Pass"],
                 highlighters=(pass,fail))
    do_test || return
    for test in tests
        @eval @test $test
    end
end


text: g(a,b) = 3a+2b


text:
# Based on https://gist.github.com/Flafla2/f0260a861be0ebdeef76

# Hash lookup table as defined by Ken Perlin
const PERMS1 = [151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7,
    225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247,
    120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
    88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134,
    139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220,
    105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80,
    73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86,
    164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38,
    147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189,
    28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101,
    155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232,
    178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12,
    191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181,
    199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236,
    205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180]

const PERMS = [PERMS1; PERMS1]

function fade{T<:AbstractFloat}(x::T)
    return 6x^5 - 15x^4 + 10x^3
end

function smoothstep{T<:AbstractFloat}(a0::T, a1::T, x::T)
    a0 != a1 || throw(ArgumentError("Arguments a0 and a1 cannot be equal"))
    x = clamp((x - a0) / (a1 - a0), 0.0, 1.0)
    return 3x^2 - 2x^3
end

function smootherstep{T<:AbstractFloat}(a0::T, a1::T, x::T)
    a0 != a1 || throw(ArgumentError("Arguments a0 and a1 cannot be equal"))
    x = clamp((x - a0) / (a1 - a0), 0.0, 1.0)
    return fade(x)
end

# Might refactor this down to a simple lerp
function interpolate{T<:AbstractFloat}(a0::T, a1::T, w::T, method::Symbol=:linear)
    0 ≤ w ≤ 1 || throw(ArgumentError("Expected 0 ≤ w ≤ 1, got $w"))
    if method ∉ [:linear, :smoothstep, :smootherstep]
        throw(ArgumentError("Unrecognized interpolation method $method"))
    end
    x = method == :linear ? w : eval(method)(a0, a1, w)
    return a0 + x * (a1 - a0)
end

function gradient{T<:AbstractFloat}(hash::Int, x::T, y::T, z::T)
    h = hash & 15
    u = h < 8 ? x : y
    v = h < 4 ? y : h == 12 || h == 14 ? x : z
    return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v)
end

function perlin{T<:AbstractFloat}(x::T, y::T, z::T)
    xi = trunc(Int, x) & 255 + 1
    yi = trunc(Int, y) & 255 + 1
    zi = trunc(Int, z) & 255 + 1

    xf = first(modf(x))
    yf = first(modf(y))
    zf = first(modf(z))

    u = fade(xf)
    v = fade(yf)
    w = fade(zf)

    aaa = PERMS[PERMS[PERMS[xi] + yi] + zi]
    aba = PERMS[PERMS[PERMS[xi] + yi + 1] + zi]
    aab = PERMS[PERMS[PERMS[xi] + yi] + zi + 1]
    abb = PERMS[PERMS[PERMS[xi] + yi + 1] + zi + 1]
    baa = PERMS[PERMS[PERMS[xi + 1] + yi] + zi]
    bba = PERMS[PERMS[PERMS[xi + 1] + yi + 1] + zi]
    bab = PERMS[PERMS[PERMS[xi + 1] + yi] + zi + 1]
    bbb = PERMS[PERMS[PERMS[xi + 1] + yi + 1] + zi + 1]

    x1 = interpolate(gradient(aaa, xf, yf, zf), gradient(baa, xf - 1, yf, zf), u)
    x2 = interpolate(gradient(aba, xf, yf - 1, zf), gradient(bba, xf - 1, yf - 1, zf), u)
    y1 = interpolate(x1, x2, v)

    x1 = interpolate(gradient(aab, xf, yf, zf - 1), gradient(bab, xf - 1, yf, zf - 1), u)
    x2 = interpolate(gradient(abb, xf, yf - 1, zf - 1), gradient(bbb, xf - 1, yf - 1, zf - 1), u)
    y2 = interpolate(x1, x2, v)

    return (interpolate(y1, y2, w) + 1) / 2
end

function octaveperlin{T<:AbstractFloat}(x::T, y::T, z::T, octaves::Int, persistence::T)
    total = 0.0
    frequency = 1.0
    amplitude = 1.0
    maxval = 0.0
    for i = 1:octaves
        total += perlin(x * frequency, y * frequency, z * frequency) * amplitude
        maxval += amplitude
        amplitude *= persistence
        frequency *= 2
    end
    return total / maxval
end


text:
using Combinatorics

function is_anagram(w1, w2)
  return length(findin(w1, w2)) == length(w1) == length(w2)
end

valid_count = 0
for passphrase in eachline(STDIN)
  words = split(passphrase, " ")
  anagram_found = false

  for (w1, w2) in combinations(words, 2)
    if is_anagram(w1, w2)
      anagram_found = true
      break
    end
  end

  if !anagram_found
    valid_count += 1
  end
end

println("Valid passphrases for part 2: $valid_count")


text:
using Statistics
"""
    CrossEntropyLoss(dims=1)
    (l::CrossEntropyLoss)(scores, answers; average=true)
Calculates negative log likelihood error on your predicted scores.
`answers` should be integers corresponding to correct class indices.
If an answer is 0, loss from that answer will not be included.
This is the masking feature when you are working with unequal length sequences.

if dims==1,
    First dimension is assumed to be predicted logits.
* size(scores) = D,[BATCH,T1,T2,...]
* size(answers)= [BATCH,T1,T2,...]

elseif dims==2,
    Last dimension is assumed to be predicted logits.
* size(scores) = [BATCH,T1,T2,...],D
* size(answers)= [BATCH,T1,T2,...]
"""
struct CrossEntropyLoss <: Loss
    dims::Integer
end
CrossEntropyLoss(;dims=1) = CrossEntropyLoss(dims)
@inline (l::CrossEntropyLoss)(y,answers::Array{<:Integer}; average=true) = nllmask(y, answers; dims=l.dims, average=average)

"""
    BCELoss(average=true)
    (l::BCELoss)(scores, answers)
    Computes binary cross entropy given scores(predicted values) and answer labels. answer values should be {0,1}, then it returns negative of
    mean|sum(answers * log(p) + (1-answers)*log(1-p)) where p is equal to 1/(1 + exp.(scores)). see `LogisticLoss`.
"""
struct BCELoss <: Loss end
@inline (l::BCELoss)(y,answers::Array{<:Integer})=bce(y,answers)

"""
    LogisticLoss(average=true)
    (l::LogisticLoss)(scores, answers)
    Computes logistic loss given scores(predicted values) and answer labels. answer values should be {-1,1}, then it returns mean|sum(log(1 +
    exp(-answers*scores))). See also `BCELoss`.
"""
struct LogisticLoss <: Loss end
@inline (l::LogisticLoss)(y,answers::Array{<:Integer}) = logistic(y,answers)

"""
    SigmoidCrossEntropyLoss()
    (l::SigmoidCrossEntropyLoss)(scores, labels)

Measures the error in discrete classification tasks in which each class is independent and not mutually exclusive.
`labels` should be same size with scores where existing classes pointed by one others zero.

if dims==1, First dimension is assumed to be predicted logits.

elseif dims==2, Last dimension is assumed to be predicted logits.
"""
struct SigmoidCrossEntropyLoss <: Loss
    dims::Int
end
SigmoidCrossEntropyLoss(;dims=1) = SigmoidCrossEntropyLoss(dims)
function (l::SigmoidCrossEntropyLoss)(x, z; average=true)
    y = sum(relu.(x) .- x .* z - log.(sigm.(abs.(x))), dims=:)
    average ? y ./ (prod(size(x))/prod(size(x)[l.dims])) : y
end

####
#### Utils
####
@inline function nllmask(y,a::AbstractArray{<:Integer}; dims=1, average=true)
    indices = findindices(y, a, dims=dims)
    lp = logp(y,dims=dims)[indices]
    average ? -mean(lp) : -sum(lp)
end

function findindices(y,a::AbstractArray{<:Integer}; dims=1)
    n       = length(a)
    nonmask = a .> 0
    indices = Vector{Int}(undef,sum(nonmask))
    if dims == 1                   # instances in first dimension
        y1 = size(y,1)
        y2 = div(length(y),y1)
        if n != y2; throw(DimensionMismatch()); end
        k = 1
        @inbounds for (j,v) in enumerate(nonmask)
            !v && continue
            indices[k] = (j-1)*y1 + a[j]
            k += 1
        end
    elseif dims == 2               # instances in last dimension
        y2 = size(y,ndims(y))
        y1 = div(length(y),y2)
        if n != y1; throw(DimensionMismatch()); end
        k = 1
        @inbounds for (j,v) in enumerate(nonmask)
            !v && continue
            indices[k] = (a[j]-1)*y1 + j
            k += 1
        end
    else
        error("findindices only supports dims = 1 or 2")
    end
    return indices
end


text:
"""
    rCIP.jl

This is a Julia port of a MATLAB implementation of batch and incremental
(Renyi's) representative Cross Information Potential (rCIP) Cluster Validity Index.

Authors:
MATLAB implementation: Leonardo Enzo Brito da Silva
Julia port: Sasha Petrenko <sap625@mst.edu>

REFERENCES
[1] L. E. Brito da Silva, N. M. Melton, and D. C. Wunsch II, "Incremental
Cluster Validity Indices for Hard Partitions: Extensions  and  Comparative
Study," ArXiv  e-prints, Feb 2019, arXiv:1902.06711v1 [cs.LG].
[2] E. Gokcay and J. C. Principe, "A new clustering evaluation function
using Renyi’s information potential," in Proc. Int. Conf. Acoust., Speech,
Signal Process. (ICASSP), vol. 6. Jun. 2000, pp. 3490–3493.
[3] E. Gokcay and J. C. Principe, "Information theoretic clustering," IEEE
Trans. Pattern Anal. Mach. Intell., vol. 24, no. 2, pp. 158–171, Feb. 2002.
[4] D. Araújo, A. D. Neto, and A. Martins, "Representative cross information
potential clustering," Pattern Recognit. Lett., vol. 34, no. 16,
pp. 2181–2191, Dec. 2013.
[5] D. Araújo, A. D. Neto, and A. Martins, "Information-theoretic clustering:
A representative and evolutionary approach," Expert Syst. Appl.,
vol. 40, no. 10, pp. 4190–4205, Aug. 2013.
[6] R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification, 2nd ed.
John Wiley & Sons, 2000.
"""

"""
    rCIP

The stateful information of the (Renyi's) representative Cross Information Potential (rCIP) Cluster Validity Index.

# References
1. L. E. Brito da Silva, N. M. Melton, and D. C. Wunsch II, "Incremental Cluster Validity Indices for Hard Partitions: Extensions  and  Comparative Study," ArXiv  e-prints, Feb 2019, arXiv:1902.06711v1 [cs.LG].
2. E. Gokcay and J. C. Principe, "A new clustering evaluation function using Renyi’s information potential," in Proc. Int. Conf. Acoust., Speech, Signal Process. (ICASSP), vol. 6. Jun. 2000, pp. 3490–3493.
3. E. Gokcay and J. C. Principe, "Information theoretic clustering," IEEE Trans. Pattern Anal. Mach. Intell., vol. 24, no. 2, pp. 158–171, Feb. 2002.
4. D. Araújo, A. D. Neto, and A. Martins, "Representative cross information potential clustering," Pattern Recognit. Lett., vol. 34, no. 16, pp. 2181–2191, Dec. 2013.
5. D. Araújo, A. D. Neto, and A. Martins, "Information-theoretic clustering: A representative and evolutionary approach," Expert Syst. Appl., vol. 40, no. 10, pp. 4190–4205, Aug. 2013.
6. R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification, 2nd ed. John Wiley & Sons, 2000.
"""
mutable struct rCIP <: CVI
    label_map::LabelMap
    dim::Integer
    n_samples::Integer
    n::IntegerVector            # dim
    v::RealMatrix               # dim x n_clusters
    sigma::Array{Float, 3}     # dim x dim x n_clusters
    constant::Float
    D::RealMatrix               # n_clusters x n_clusters
    delta_term::RealMatrix      # dim x dim
    n_clusters::Integer
    criterion_value::Float
end # rCIP <: CVI

"""
    rCIP()

Default constructor for the (Renyi's) representative Cross Information Potential (rCIP) Cluster Validity Index.
"""
function rCIP()
    rCIP(
        LabelMap(),                         # label_map
        0,                                  # dim
        0,                                  # n_samples
        Array{Integer, 1}(undef, 0),        # n
        Array{Float, 2}(undef, 0, 0),       # v
        Array{Float, 3}(undef, 0, 0, 0),    # sigma
        0.0,                                # constant
        Array{Float, 2}(undef, 0, 0),       # D
        Array{Float, 2}(undef, 0, 0),       # delta_term
        0,                                  # n_clusters
        0.0                                 # criterion_value
    )
end # rCIP()

"""
    setup!(cvi::rCIP, sample::Vector{T}) where {T<:RealFP}
"""
function setup!(cvi::rCIP, sample::Vector{T}) where {T<:RealFP}
    # Get the feature dimension
    cvi.dim = length(sample)
    # Initialize the 2-D and 3-D arrays with the correct feature dimension
    cvi.v = Array{T, 2}(undef, cvi.dim, 0)
    cvi.sigma = Array{T, 3}(undef, cvi.dim, cvi.dim, 0)
    # Calculate the delta term
    epsilon = 12
    delta = 10^(-epsilon/cvi.dim)
    cvi.delta_term = Matrix{Float64}(LinearAlgebra.I, cvi.dim, cvi.dim).*delta
end # setup!(cvi::rCIP, sample::Vector{T}) where {T<:RealFP}

function param_inc!(cvi::rCIP, sample::RealVector, label::Integer)
    # Get the internal label
    i_label = get_internal_label!(cvi.label_map, label)

    n_samples_new = cvi.n_samples + 1
    if isempty(cvi.v)
        setup!(cvi, sample)
    end
    cvi.constant = 1/sqrt((2*pi)^cvi.dim)

    if i_label > cvi.n_clusters
        n_new = 1
        v_new = sample
        sigma_new = cvi.delta_term
        if cvi.n_clusters == 0
            D_new = zeros(1, 1)
        else
            D_new = zeros(cvi.n_clusters + 1, cvi.n_clusters + 1)
            D_new[1:cvi.n_clusters, 1:cvi.n_clusters] = cvi.D
            d_column_new = zeros(cvi.n_clusters + 1)
            for jx = 1:cvi.n_clusters
                diff_m = v_new - cvi.v[:, jx]
                sigma_q = sigma_new + cvi.sigma[:,:,jx]
                d_column_new[jx] = cvi.constant * (1/sqrt(det(sigma_q)))*exp(-0.5*transpose(diff_m)*inv(sigma_q)*diff_m)
                # d_column_new[jx] = sum((v_new - cvi.v[:, jx]).^2)
            end
            D_new[:, i_label] = d_column_new
            D_new[i_label, :] = transpose(d_column_new)
        end
        # Update 1-D parameters with a push
        cvi.n_clusters += 1
        push!(cvi.n, n_new)
        # Update 2-D parameters with appending and reassignment
        cvi.v = [cvi.v v_new]
        cvi.D = D_new
        cvi.sigma = cat(cvi.sigma, sigma_new, dims=3)
    else
        n_new = cvi.n[i_label] + 1
        v_new = (1 - 1/n_new) .* cvi.v[:, i_label] + (1/n_new) .* sample
        diff_x_v = sample - cvi.v[:, i_label]
        if n_new > 1
            sigma_new = ((n_new - 2)/(n_new - 1))*(cvi.sigma[:,:,i_label] - cvi.delta_term) +
                (1/n_new)*(diff_x_v*transpose(diff_x_v)) + cvi.delta_term
        else
            sigma_new = cvi.delta_term
        end
        d_column_new = zeros(cvi.n_clusters)
        for jx = 1:cvi.n_clusters
            if jx == i_label
                continue
            end
            diff_m = v_new - cvi.v[:, jx]
            sigma_q = sigma_new + cvi.sigma[:,:,jx]
            d_column_new[jx] = cvi.constant*(1/sqrt(det(sigma_q)))*exp(-0.5*transpose(diff_m)*inv(sigma_q)*diff_m)
        end
        # Update parameters
        cvi.n[i_label] = n_new
        cvi.v[:, i_label] = v_new
        cvi.sigma[:,:,i_label] = sigma_new
        cvi.D[:, i_label] = d_column_new
        cvi.D[i_label, :] = transpose(d_column_new)
    end
    cvi.n_samples = n_samples_new
end # param_inc!(cvi::rCIP, sample::RealVector, label::Integer)

function param_batch!(cvi::rCIP, data::RealMatrix, labels::IntegerVector)
    cvi.dim, cvi.n_samples = size(data)
    cvi.constant = 1/sqrt((2*pi)^cvi.dim)
    # Calculate the delta term
    epsilon = 12
    delta = 10^(-epsilon/cvi.dim)
    cvi.delta_term = Matrix{Float64}(LinearAlgebra.I, cvi.dim, cvi.dim).*delta

    # Take the average across all samples, but cast to 1-D vector
    u = unique(labels)
    cvi.n_clusters = length(u)
    cvi.n = zeros(Integer, cvi.n_clusters)
    cvi.v = zeros(cvi.dim, cvi.n_clusters)
    cvi.sigma = zeros(cvi.dim, cvi.dim, cvi.n_clusters)
    cvi.D = zeros(cvi.n_clusters, cvi.n_clusters)

    for ix = 1:cvi.n_clusters
        subset = data[:, findall(x->x==u[ix], labels)]
        cvi.n[ix] = size(subset, 2)
        cvi.v[1:cvi.dim, ix] = mean(subset, dims=2)
        if cvi.n[ix] > 1
            cvi.sigma[:,:,ix] = (1/(cvi.n[ix] - 1)) *
                ((subset*transpose(subset)) - cvi.n[ix].*cvi.v[:,ix]*transpose(cvi.v[:,ix])) +
                cvi.delta_term
        else
            cvi.sigma[:,:,ix] = cvi.delta_term
        end
    end
    for ix = 1 : (cvi.n_clusters - 1)
        for jx = ix + 1 : cvi.n_clusters
            diff_m = cvi.v[:, ix] - cvi.v[:, jx]
            sigma_q = cvi.sigma[:,:,ix] + cvi.sigma[:,:,jx]
            cvi.D[jx, ix] = cvi.constant*(1/sqrt(det(sigma_q)))*exp(-0.5*transpose(diff_m)*inv(sigma_q)*diff_m)
        end
    end
    cvi.D = cvi.D + transpose(cvi.D)
end # param_batch!(cvi::rCIP, data::RealMatrix, labels::IntegerVector)

function evaluate!(cvi::rCIP)
    # Assume a symmetric dimension
    dim = size(cvi.D)[1]
    if dim > 1
        values = [cvi.D[i,j] for i = 1:dim, j=1:dim if j > i]
        cvi.criterion_value = sum(values)
    else
        cvi.criterion_value = 0.0
    end
end # evaluate(cvi::rCIP)


text:
"""
    trim(x; dims::Tuple, pad::Int)

Trim `missingval(x)` from `x` for axes in `dims`, returning a view of `x`.

# Arguments

- `x`: A `Raster` or `RasterStack`. For stacks, all layers must having
    missing values for a pixel for it to be trimmed.

# Keywords

- `dims`: By default `dims=(XDim, YDim)`, so that trimming keeps the area
    of `X` and `Y` that contains non-missing values along all other dimensions.
- `pad`: The trimmed size will be padded by `pad` on all sides, although
    padding will not be added beyond the original extent of the array.

`trim` does not accept `filename`/`suffix` arguments as it does not alter the underlying data.

# Example

Create trimmed layers of Australian habitat heterogeneity.

```jldoctest
using Rasters, Plots
layers = (:evenness, :range, :contrast, :correlation)
st = RasterStack(EarthEnv{HabitatHeterogeneity}, layers)
plot(st)

# Roughly cut out australia
ausbounds = X(Between(100, 160)), Y(Between(-10, -50))
aus = st[ausbounds...]
a = plot(aus)

# Trim missing values and plot
b = plot(trim(aus))

savefig(a, "build/trim_example_before.png")
savefig(b, "build/trim_example_after.png")
# output
```

### Before `trim`:

![before trim](trim_example_before.png)

### After `trim`:

![after trim](trim_example_after.png)

$EXPERIMENTAL
"""
function trim(x::RasterStackOrArray; dims::Tuple=(XDim, YDim), pad::Int=0)
    # Get the actual dimensions in their order in the array
    targetdims = commondims(x, dims)
    # Get the range of non-missing values for each dimension
    trackers = AxisTrackers(DD.dims(x), targetdims)
    _update!(trackers, x)
    ranges = _cropranges(trackers)
    padded_ranges = _pad(ranges, targetdims, pad)
    rangedims = map(rebuild, targetdims, padded_ranges)
    return view(x, rangedims...)
end
function trim(ser::AbstractRasterSeries; dims::Tuple=(XDim, YDim), pad::Int=0)
    x = first(ser)
    # Get the actual dimensions in their order in the array
    targetdims = commondims(x, dims)
    # Get the range of non-missing values for each dimension
    trackers = AxisTrackers(DD.dims(x), targetdims)
    map(x -> _update!(trackers, x), ser)
    # Get the ranges that contain all non-missing values
    ranges = _cropranges(trackers)
    padded_ranges = _pad(ranges, targetdims, pad)
    rangedims = map(rebuild, targetdims, padded_ranges)
    return map(x -> view(x, rangedims...), ser)
end

function _pad(ranges, dims, pad)
    map(ranges, size(dims)) do r, l
        max(first(r)-pad, 1):min(last(r)+pad, l)
    end
end

# Tracks the status of an index for some subset of dimensions of an Array
# This lets us track e.g. the X/Y indices that have only missing values
# accross all other dimensions.
# This is a hack to work with DiskArrays broadcast chunking without allocations.
struct AxisTrackers{N,Tr,D,TD} <: AbstractArray{Bool,N}
    tracking::Tr
    dims::D
    trackeddims::TD
end
function AxisTrackers(tracking::T, dims::D, trackeddims::TD) where {T,D,TD}
    AxisTrackers{length(dims),T,D,TD}(tracking, dims, trackeddims)
end
function AxisTrackers(dims::Tuple, trackeddims::Tuple)
    tracking = map(trackeddims) do td
        (_ -> false).(td)
    end
    return AxisTrackers(tracking, dims, trackeddims)
end

Base.axes(A::AxisTrackers) = map(d -> axes(d, 1), A.dims)
Base.size(A::AxisTrackers) = map(length, A.dims)
Base.getindex(A::AxisTrackers, I...) = map(getindex, A.tracking, _trackedinds(I)) |> any
function Base.setindex!(A::AxisTrackers, x, I::Int...)
    map(A.tracking, _trackedinds(A, I)) do axis, i
        axis[i] |= x
    end
end

function _trackedinds(A, I)
    # Wrap indices in dimensions so we can sort and filter them
    Id = map((d, i) -> DD.basetypeof(d)(i), A.dims, I)
    # Get just the tracked dimensions
    Itracked = dims(Id, A.trackeddims)
    # Get the indices for the tracked dimensions
    return map(val, Itracked)
end

function _cropranges(trackers::AxisTrackers)
    map(trackers.tracking) do a
        f = findfirst(a)
        l = findlast(a)
        f = f === nothing ? firstindex(a) : f
        l = l === nothing ? lastindex(a) : l
        f:l
    end
end

# Broadcast over the array and tracker to mark axis indices as being missing or not
_update!(tr::AxisTrackers, A::AbstractRaster) = tr .= A .!== missingval(A)
_update!(tr::AxisTrackers, st::AbstractRasterStack) = map(A -> tr .= A .!== missingval(A), st)


text:
using RData
using Test

include("basic.jl")
include("RDA.jl")
include("RDS.jl")


text:
    f = open(join([version_dir, "optim_benchmark.csv"], "/"), "w")
    write(f, join(["Problem", "Optimizer", "Converged", "Time", "Minimum", "Iterations", "f_calls", "g_calls", "f_hat", "f_error", "x_error"], ","))
    write(f, "\n")

    @showprogress 1 "Benchmarking..." for (name, problem) in Optim.UnconstrainedProblems.examples
        for (i,  algorithm) in enumerate(default_solvers)
            print_debuginfo && @show name, problem , algorithm
            try
            # Force compilation and obtain results
            results = optimize(problem.f, problem.g!, problem.h!,
                               problem.initial_x, algorithm, Optim.Options(g_tol = 1e-16))
            # Run each algorithm n times
            n = 10
            if algorithm == ParticleSwarm() && name == "Large Polynomial"
                n = 1
            end
            # Estimate run time in seconds
            run_time = minimum([@elapsed optimize(problem.f, problem.g!, problem.h!,
                                   problem.initial_x,
                                   algorithm, Optim.Options(g_tol = 1e-16)) for nn = 1:n])

            # Count iterations
            iterations = results.iterations

            # Print out results.
            write(f, join([problem.name,
                           default_names[i],
                           Optim.converged(results),
                           run_time,
                           Optim.minimum(results),
                           iterations,
                           Optim.f_calls(results),
                           Optim.g_calls(results),
                           problem.f(problem.solutions),
                           Optim.minimum(results)-problem.f(problem.solutions),
                           norm(Optim.minimizer(results)-problem.solutions, Inf)], ","))
            write(f, "\n")
            catch
                write(f, join([problem.name,
                               default_names[i],
                               "false",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf"], ","))
                write(f, "\n")
            end
        end
    end
    close(f)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ScenarioBasedInternal <: ForecastInternal
        label::String
        resolution::Dates.Period
        initial_time::Dates.DateTime
        scenario_count::Int64
        time_series_uuid::UUIDs.UUID
        horizon::Int
        internal::InfrastructureSystemsInternal
    end

A Discrete Scenario Based forecast for a particular data field in a Component.

# Arguments
- `label::String`: label of component parameter forecasted
- `resolution::Dates.Period`
- `initial_time::Dates.DateTime`: forecast availability time
- `scenario_count::Int64`: Number of scenarios
- `time_series_uuid::UUIDs.UUID`: reference to time series data; timestamp - scalingfactor
- `horizon::Int`: length of this forecast
- `internal::InfrastructureSystemsInternal`
"""
mutable struct ScenarioBasedInternal <: ForecastInternal
    "label of component parameter forecasted"
    label::String
    resolution::Dates.Period
    "forecast availability time"
    initial_time::Dates.DateTime
    "Number of scenarios"
    scenario_count::Int64
    "reference to time series data; timestamp - scalingfactor"
    time_series_uuid::UUIDs.UUID
    "length of this forecast"
    horizon::Int
    internal::InfrastructureSystemsInternal
end

function ScenarioBasedInternal(label, resolution, initial_time, scenario_count, time_series_uuid, horizon, )
    ScenarioBasedInternal(label, resolution, initial_time, scenario_count, time_series_uuid, horizon, InfrastructureSystemsInternal(), )
end

function ScenarioBasedInternal(; label, resolution, initial_time, scenario_count, time_series_uuid, horizon, )
    ScenarioBasedInternal(label, resolution, initial_time, scenario_count, time_series_uuid, horizon, )
end

"""Get ScenarioBasedInternal label."""
get_label(value::ScenarioBasedInternal) = value.label
"""Get ScenarioBasedInternal resolution."""
get_resolution(value::ScenarioBasedInternal) = value.resolution
"""Get ScenarioBasedInternal initial_time."""
get_initial_time(value::ScenarioBasedInternal) = value.initial_time
"""Get ScenarioBasedInternal scenario_count."""
get_scenario_count(value::ScenarioBasedInternal) = value.scenario_count
"""Get ScenarioBasedInternal time_series_uuid."""
get_time_series_uuid(value::ScenarioBasedInternal) = value.time_series_uuid
"""Get ScenarioBasedInternal horizon."""
get_horizon(value::ScenarioBasedInternal) = value.horizon
"""Get ScenarioBasedInternal internal."""
get_internal(value::ScenarioBasedInternal) = value.internal