Programming with Redberry

Basics

Redberry interface is written in Groovy and is intended to be used within the Groovy environment. Groovy is a general-purpose programming language and one can use all features and programming language constructs that are available in Groovy: looping, branching, functions, lambda-expressions, lists, classes etc. Besides, Redberry provides a specialized domain-specific programming constructs which are useful for computer algebra purposes.

The documentation on general programming constructs can be found on Groovy website. The most useful things are:

One can find a comprehensive documentation on other Groovy features on Groovy website.

Examples

Let us consider some applications of Redberry that involve programming.

Basic constructs

Consider some basic programming constructs:

• Looping

Loop over expression (or list/set/array etc.):

def t = 'a + b + c'.t
for(def a in t)
    println a
//or equivalently
for(int i=0; i < t.size(); ++i)
    println t[i]
//one more way
t.each { a->
    println a
}
//using while
def i =0
while(i < t.size()){
    println t[i++]
}
For further details see Groovy documentation:

• Logical branching

Typical if-else statement:

def t = 'a + b + c'.t
if( t.size() > 3){
    //do something
} else if (t.size() == 3){
    //do something else
} else{
    //do something  else else
}
It is important to note, that in order to compare tensors, one should use .equals() method instead of == operator (this issue will be fixed in Groovy 3.0 and one can use == for comparison everywhere):
def expr1 = 'a_i + b_i + c_i'.t
def expr2 = 'a_j + b_j + c_j'.t

if ( expr1.equals(expr2) ){
    //somth
}

For further details see Groovy documentation:

• Functions

Inside a Groovy script one can define function using closure:

def pow3 = { x -> x**3 }
println pow3(3) //gives 27

def max = { x, y -> x > y ? x : y }
println max(2, 3) //gives 3
The last statement inside a closure is automatically considered as return statement.

In Redberry a function that transform expression to another expression is a Transformation. In order to convert a closure to a transformation one can do:

def tr = { expr ->
    //invert indices
    (expr.indices % expr.indices.inverted) >> expr
} as Transformation

println tr >> 't_ab'.t // gives t^ab
//use & to join transformation 
println((EliminateMetrics & tr) >> 'g^ab*t_ac'.t) //gives t_b^c

For further details see Groovy documentation:

• Iterables

Lists and expressions are iterable, which means that one can use for example .each { } in order to iterate over expression elements. Additionally, there are some other useful methods that allow to iterate over a list (or expression) and select elements:

def t = 'a_i + b_i+ c_i'.t
t.eachWithIndex { e, i -> println "$i: $e" } // gives 0: a_i 1: b_i 2: c_i

println t.collect { 2 * it } // gives [2*a_i, 2*b_i, 2*c_i]

println t.find { (it % 'a_k'.t).first != null } // gives a_i

println t.findAll { it.indices.lower.size() < 2 } // gives [a_i, b_i, c_i]
To convert a list to Sum or Product one can use .sum() or .multiply() methods:
 // gives a_i + b_i + c_i
println t.findAll( { it.indices.lower.size() < 2 } ).sum()
For further details see Groovy documentation:

Toy example

Let us consider function that selects all elements with size less than n from Sum and write them to list:

def trunc = { expr, n ->
    //if expr is not Sum --- return empty List
    if (expr.class != Sum)
        return []
    //resulting list
    def r = []
    //loop over summands
    for (def s in expr) {
        //check summand size
        if (s.class == Product && s.size() >= n)
            continue;
        r << s
    }
    return r
}

println trunc('a + a*b + a*b*c + a*b*c*d'.t, 3)
   > [a, a*b]
println trunc('a + a*b + a*b*c + a*b*c*d'.t, 4)
   > [a, a*b, a*b*c]

The for-loop at lines 8-13 can be rewritten equivalently using each closure:

expr.each { s ->
    //check summand size
    if (s.class != Product || s.size() < n)
        r << s
}
Convert resulting list to a new Sum:
def selected = trunc('a + a*b + a*b*c + a*b*c*d'.t, 4)
def newSum = selected.sum()
println newSum
   > a + a*b + a*b*c
Implement Redberry Transformation that removes all elements with size greater or equals than 4 from Sum:
//transformation that removes elements from sum with size >=4
def truncTr4 = { expr ->
    expr.class == Sum ? trunc(expr, 4).sum() : expr
} as Transformation
//apply transformation using >>
println truncTr4 >> 'a + a*b + a*b*c + a*b*c*d'.t
   > a + a*b + a*b*c
Here we omitted the return keyword at line 24. Let us generalise Transformation truncTr4 to work with arbitrary n:
//transformation that removes elements from sum with size >=n
def truncTr = { n ->
    { expr ->
        expr.class == Sum ? trunc(expr, n).sum() : expr
    } as Transformation
}
//apply transformation using >>
println truncTr(4) >> 'a + a*b + a*b*c + a*b*c*d'.t
   > a + a*b + a*b*c
Transformation truncTr will apply only to the top algebraic level. In order to implement Transformation that will change all sums in expression, one can use Tree traversal:
//transformation that applies to each part of expression
def truncAll = { n ->
    { expr ->
        expr.transformParentAfterChild { e ->
            e.class == Sum ? truncTr(n) >> e : e
        }
    } as Transformation
}
//this will do nothing
println truncTr(4) >> 'x*(a + a*b + a*b*c + a*b*c*d)'.t
   > x*(a + a*b + a*b*c + a*b*c*d)
//this will be applied
println truncAll(4) >> 'x*(a + a*b + a*b*c + a*b*c*d)'.t
   > x*(a + a*b + a*b*c)

Advanced example: Fourier transform of Lagrangian

Let us write a function that performs Fourier transform of Langrangian. So, for example:

\[ \int d^4 x\left( -\frac{1}{4}\left(\partial_\mu A_\nu (x) - \partial_\nu A_\mu (x)\right)^2 \right)= \int d^4 p\left( \frac{1}{2}A_\mu(p) \left(g^{\mu\nu} p^2 - p^\mu p^\nu \right) A_\nu(-p) \right) \]

Let's implement Redberry Transformation that transforms l.h.s. integrand to r.h.s integrand.

Obviously this transformation will separately transform each term in Lagrangian. So, let's first implement Fourier transform of a single Product. For simplicity we assume that we have only one field. The algorithm sequentially goes through Product multipliers; for each field (e.g. $\partial_{a}\partial_{b} A_c(x_a)$) it generates a new momentum (e.g. $p2_a$) and replaces field argument with momentum and partial derivatives with products of momentum (e.g. $\partial_{a}\partial_{b} A_c(x_a) \to i \times p2_a \times i \times p2_b \times A_c(p_a)$). The last generated momentum should be replaced with negated sum of all other momentums (e.g. $p3_i \to -p0_i - p1_i -p2_i$). Here's the implementation:

//transform product of tensors
def fourierProduct = { product ->
    //list of generated momentums
    def momentums = []
    //the result
    def result = product.builder
    //counter of momentums
    def i = 0
    //let's transform each term in product
    for (def term in product) {
        //transform those terms that are functions
        if (term.class == TensorField) {
            //generate next momentum
            def momentum = "p${i++}".t
            momentums << momentum
            //replace function argument with momentum
            // (e.g. f~(2)_{a bc}[x_a] -> f~(2)_{a bc}[p_a])
            term = "${term[0]} = $momentum${term[0].indices}".t >> term
            //in case of derivative we need
            // to replace partials with momentums
            if (term.isDerivative()) {
                //indices of differentiating variables
                def dIndices = term.partitionOfIndices[1]
                //extract just parent field from derivative
                // (e.g. f~(2)_{a bc}[p_a] -> f_a[p_a])
                term = term.parentField
                //multiply by momentums
                // (e.g. f~(2)_{a bc}[p_a] -> I*p_b*I*p_c*f_a[p_a])
                dIndices.each { indices ->
                    term *= "I * $momentum $indices".t
                }
            }
        }
        //put transformed term to new product
        result << term
    }
    //result
    def r = result.build()
    //we must replace the last momentum with -(sum of other momentums)
    def rhs = '0'.t
    //sum generated momentums except last one
    momentums.eachWithIndex { momentum, c ->
        if (c != momentums.size() - 1)
            rhs -= "${momentum}_a".t
    }
    //replace last momentum with sum of others and return
    "${momentums[momentums.size() - 1]}_a = $rhs".t >> r
} as Transformation

The final transformation for a sum of terms is:

//transform sum of products
def fourierSum = { expr ->
    //expand all brackets and unfold powers of scalar tensors
    expr = (ExpandAndEliminate & PowerUnfold) >> expr
    //apply fourierProduct to each summand:
    expr = expr.collect { s -> fourierProduct >> s }.sum()
    //simplify and return:
    ExpandAndEliminate >> expr
} as Transformation

Consider, for example, Lagrangian of electromagnetic field :

def emTensor = 'F_ab := A~(1)_ab[x_a] - A~(1)_ba[x_a]'.t
def lagrangian = '-(1/4)*F_ab*F^ab - 1/(2*x)*(A~(1)_a^a[x_a])**2'.t
def fourier = (fourierSum & Collect['A_a[p_a]'.t, Factor]) >> lagrangian
println fourier
   > A_{a}[p0_{a}]*A^{c}[-p0_{a}]*
           ((1/2)*x**(-1)*(-1+x)*p0^{a}*p0_{c}-(1/2)*p0_{b}*p0^{b}*d^{a}_{c})
One can use this result in order to inverse the quadratic form in the Lagrangian (using Reduce function) and obtain propagator:
def q2 = 'K^a_c := (1/2)*x**(-1)*(-1+x)*p0^a*p0_c-(1/2)*p0_b*p0^b*d^a_c'.t
def options = [ExternalSolver : [
                       Solver: 'Mathematica',
                       Path  : '/usr/local/bin']]
def emProp = Reduce(['2*K^a_i*P^i_b = d^a_b'.t], ['P_ab'], options)
println emProp
   > [[P_{ab} =
           (p0^{c}*p0_{c})**(-2)*(1-x)*p0_{a}*p0_{b}-(p0^{c}*p0_{c})**(-1)*g_{ab}]]
which is a well known result:

\[ \frac{1}{p^2} \left( -g_{ab} + (1-x)\frac{p_a p_b}{p^2}\right) \]

Another example: let's find a propagator for a Fierz-Pauli Lagrangian in $D$ dimensions:

\[ L = -\frac{1}{2} \partial_l h_{mn} \partial^l h^{mn} + \partial_m h_{nl} \partial^n h^{ml} - \partial^m h_{ml} \partial^l h^n_n + \frac{1}{2} \partial_l h^m_m \partial^l h^n_n -\frac{1}{2} m^2 \left (h_{mn} h^{mn} - h^a_a h^b_b \right) \]

def FierzPauli = '''
             -(1/2)*h~(1)_mnl[x_a]*h~(1)^mnl[x_a]
             + h~(1)_nlm[x_a]*h~(1)^mln[x_a]
             - h~(1)_ml^m[x_a]*h~(1)_n^nl[x_a]
             + (1/2)*h~(1)^m_ml[x_a]*h~(1)^n_n^l[x_a]
             -(1/2)*m**2*(h_mn[x_a]*h^mn[x_a] - h^a_a[x_a]*h^b_b[x_a])
             '''.t
fourier = (fourierSum & Collect['h_ab[p_a]'.t, Factor]) >> FierzPauli
println fourier
   > h_{ml}[-p0_{a}]*h_{n}^{a}[p0_{a}]*
          (-p0^{l}*p0^{m}*d^{n}_{a}+p0^{m}*p0^{n}*d_{a}^{l}
          -(1/2)*(m**2+p0_{b}*p0^{b})*d_{a}^{l}*g^{nm}
          +(1/2)*(m**2+p0_{b}*p0^{b})*d_{a}^{n}*g^{lm})
Given this quadratic form, one can find a propagator in the following way:
//quadratic form
q2 = '''(-p0^{l}*p0^{m}*d^{n}_{a}+p0^{m}*p0^{n}*d_{a}^{l}
         -(1/2)*(m**2+p0_{b}*p0^{b})*d_{a}^{l}*g^{nm}
         +(1/2)*(m**2+p0_{b}*p0^{b})*d_{a}^{n}*g^{lm})'''.t
//making symmetric with respect to field indices
def p1 = '^ml'.si
p1.symmetries.setSymmetric()
def p2 = '^n_a'.si
p2.symmetries.setSymmetric()
q2 = (Symmetrize[p1] & Symmetrize[p2]) >> q2
//make a substitution
"K^mln_a := $q2".t

//the propagator symmetries
addSymmetry 'P^ab_mn', [[0, 1]].p
addSymmetry 'P^ab_mn', [[0, 2], [1, 3]].p
options = [Transformations: 'd_n^n = D'.t,
           ExternalSolver : [
                   Solver: 'Mathematica',
                   Path  : '/usr/local/bin']]
//equation
def eq = '(K_abcd + K_cdab)*P^abmn = (d_c^m*d_d^n+d_c^n*d_d^m)/2'.t
def grProp = Reduce([eq], ['P_abmn'], options)
println grProp
   > [[P_{abmn} = -(1/2)*(m**2+p0^{e}*p0_{e})**(-1)*g_{an}*g_{bm}
          -(1/2)*(m**2+p0^{e}*p0_{e})**(-1)*g_{am}*g_{bn}
          -(1/2)*m**(-2)*(m**2+p0^{e}*p0_{e})**(-1)*p0_{a}*p0_{m}*g_{bn}
          -(1/2)*m**(-2)*(m**2+p0^{e}*p0_{e})**(-1)*p0_{b}*p0_{n}*g_{am}
          -(1/2)*m**(-2)*(m**2+p0^{e}*p0_{e})**(-1)*p0_{b}*p0_{m}*g_{an}
          -(1/2)*m**(-2)*(m**2+p0^{e}*p0_{e})**(-1)*p0_{a}*p0_{n}*g_{bm}
          +(m**2+p0^{e}*p0_{e})**(-1)*(-1+D)**(-1)*g_{ab}*g_{mn}
          -m**(-4)*(m**2+p0^{e}*p0_{e})**(-1)*(-1+D)**(-1)*(-2+D)
          *p0_{a}*p0_{b}*p0_{m}*p0_{n}+m**(-2)*(m**2+p0^{e}*p0_{e})**(-1)
          *(-1+D)**(-1)*p0_{m}*p0_{n}*g_{ab}+m**(-2)*(m**2+p0^{e}*p0_{e})**(-1)
          *(-1+D)**(-1)*p0_{a}*p0_{b}*g_{mn}]]
Which is also a well known result that can be equivalently written as \[ P_{absl} = -\frac{1}{m^2 + p^2} \left( \frac{1}{2}(P_{as} P_{bl} + P_{al} P_{bs}) - \frac{1}{D-1} P_{ab} P_{sl} \right), \] where $P_{ab} = g_{ab} + p_a p_b / m^2$.