# Differences

This shows you the differences between two versions of the page.

 documentation:guide:programming_with_redberry [2015/11/21 12:33] documentation:guide:programming_with_redberry [2015/11/21 12:33] (current) Line 1: Line 1: + ====== Programming with Redberry ====== + <​html>​ +
+ ​ + Next topic: [[documentation:​guide:​notes_on_internal_architecture]] + <​html>​ + ​ + ​ + ​ + + ---- + =====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$.