====== Programming with Redberry ======
Next topic: [[documentation:guide:notes_on_internal_architecture]]
----
=====Basics=====
Redberry interface is written in [[http://groovy.codehaus.org|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:
* Looping: [[http://groovy.codehaus.org/Looping|http://groovy.codehaus.org/Looping]]
* Branching: [[http://groovy.codehaus.org/Logical+Branching|http://groovy.codehaus.org/Logical+Branching]]
* Collections:
* Overview: [[http://groovy.codehaus.org/Collections|http://groovy.codehaus.org/Collections]]
* Detailed examples: [[http://groovy.codehaus.org/JN1015-Collections|http://groovy.codehaus.org/JN1015-Collections]]
* Functions: [[http://groovy.codehaus.org/Closures|http://groovy.codehaus.org/Closures]]
* Strings: [[http://groovy.codehaus.org/JN1525-Strings|http://groovy.codehaus.org/JN1525-Strings]]
One can find a comprehensive documentation on other Groovy features on [[http://groovy.codehaus.org|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:
* Looping: [[http://groovy.codehaus.org/Looping|http://groovy.codehaus.org/Looping]]
* Collections:
* Overview: [[http://groovy.codehaus.org/Collections|http://groovy.codehaus.org/Collections]]
* Detailed examples: [[http://groovy.codehaus.org/JN1015-Collections|http://groovy.codehaus.org/JN1015-Collections]]
=== • 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:
* Branching: [[http://groovy.codehaus.org/Logical+Branching|http://groovy.codehaus.org/Logical+Branching]]
=== • 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 [[documentation:ref: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:
* Functions: [[http://groovy.codehaus.org/Closures|http://groovy.codehaus.org/Closures]]
=== • 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 [[documentation:ref:Sum]] or [[documentation:ref: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:
* Collections: [[http://groovy.codehaus.org/JN1015-Collections|http://groovy.codehaus.org/JN1015-Collections]]
====Toy example====
Let us consider function that selects all elements with size less than ''n'' from [[documentation:ref: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'' [[http://groovy.codehaus.org/Closures|closure]]:
expr.each { s ->
//check summand size
if (s.class != Product || s.size() < n)
r << s
}
Convert resulting list to a new [[documentation:ref: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 [[documentation:ref:Transformation]] that removes all elements with size greater or equals than ''4'' from [[documentation:ref: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 [[documentation:ref: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
[[documentation:ref:Transformation]] ''truncTr'' will apply only to the top algebraic level. In order to implement [[documentation:ref: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 [[documentation:ref: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 [[documentation:ref:product]]. For simplicity we assume that we have only one field. The algorithm sequentially goes through [[documentation:ref: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 [[documentation:ref: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$.