Differences

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

Link to this comparison view

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>​
 +<div class="​text-right"​ style="​font-size:​ 15px; ">
 +</​html>​
 +Next topic: [[documentation:​guide:​notes_on_internal_architecture]]
 +<​html>​
 +<span class="​glyphicon glyphicon-arrow-right"></​span>​
 +</​div>​
 +</​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.):
 +<sxh groovy; gutter: false>
 +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++]
 +}
 +</​sxh>​
 +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:
 +<sxh groovy; gutter: false>
 +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
 +}
 +</​sxh>​
 +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):​
 +<sxh groovy; gutter: false>
 +def expr1 = 'a_i + b_i + c_i'.t
 +def expr2 = 'a_j + b_j + c_j'.t
 +
 +if ( expr1.equals(expr2) ){
 +    //somth
 +}
 +</​sxh>​
 +
 +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:
 +<sxh groovy; gutter: false>
 +def pow3 = { x -> x**3 }
 +println pow3(3) //gives 27
 +
 +def max = { x, y -> x > y ? x : y }
 +println max(2, 3) //gives 3
 +</​sxh>​
 +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:
 +<sxh groovy; gutter: false>
 +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
 +</​sxh>​
 +
 +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:
 +<sxh groovy; gutter: false>
 +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]
 +</​sxh>​
 + To convert a list to [[documentation:​ref:​Sum]] or [[documentation:​ref:​Product]] one can use ''​.sum()''​ or ''​.multiply()''​ methods:
 +<sxh groovy; gutter: false>
 + // gives a_i + b_i + c_i
 +println t.findAll( { it.indices.lower.size() < 2 } ).sum()
 +</​sxh>​
 +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:
 +<sxh groovy; gutter: true>
 +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)
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > [a, a*b]
 +</​sxh>​
 +<sxh groovy; gutter: true; first-line: 18 >
 +println trunc('​a + a*b + a*b*c + a*b*c*d'​.t,​ 4)
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > [a, a*b, a*b*c]
 +</​sxh>​
 +
 +The for-loop at lines 8-13 can be rewritten equivalently using ''​each''​ [[http://​groovy.codehaus.org/​Closures|closure]]:​
 +<sxh groovy; gutter: false>
 +expr.each { s ->
 +    //check summand size
 +    if (s.class != Product || s.size() < n)
 +        r << s
 +}
 +</​sxh>​
 +Convert resulting list to a new [[documentation:​ref:​Sum]]:​
 +<sxh groovy; gutter: true; first-line: 19>
 +def selected = trunc('​a + a*b + a*b*c + a*b*c*d'​.t,​ 4)
 +def newSum = selected.sum()
 +println newSum
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > a + a*b + a*b*c
 +</​sxh>​
 +Implement Redberry [[documentation:​ref:​Transformation]] that removes all elements with size greater or equals than ''​4''​ from [[documentation:​ref:​Sum]]:​
 +<sxh groovy; gutter: true; first-line: 22>
 +//​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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > a + a*b + a*b*c
 +</​sxh>​
 +Here we omitted the ''​return''​ keyword at line 24. Let us generalise [[documentation:​ref:​Transformation]] ''​truncTr4''​ to work with arbitrary ''​n'':​
 +<sxh groovy; gutter: true; first-line: 28>
 +//​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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > a + a*b + a*b*c
 +</​sxh>​
 +[[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]]:​
 +<sxh groovy; gutter: true; first-line: 36>
 +//​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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > x*(a + a*b + a*b*c + a*b*c*d)
 +</​sxh>​
 +<sxh groovy; gutter: true; first-line: 46>
 +//this will be applied
 +println truncAll(4) >> 'x*(a + a*b + a*b*c + a*b*c*d)'​.t
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > x*(a + a*b + a*b*c)
 +</​sxh>​
 +
 +====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:​
 +<sxh groovy; gutter: true>
 +//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
 +</​sxh>​
 +
 +The final transformation for a sum of terms is:
 +<sxh groovy; gutter: true; first-line: 49>
 +//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
 +</​sxh>​
 +
 +Consider, for example, Lagrangian of electromagnetic field :
 +
 +<sxh groovy; gutter: true; first-line: 58>
 +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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > 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})
 +</​sxh>​
 +One can use this result in order to inverse the quadratic form in the Lagrangian (using [[documentation:​ref:​Reduce]] function) and obtain propagator:
 +<sxh groovy; gutter: true; first-line: 62>
 +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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > [[P_{ab} =
 +           ​(p0^{c}*p0_{c})**(-2)*(1-x)*p0_{a}*p0_{b}-(p0^{c}*p0_{c})**(-1)*g_{ab}]]
 +</​sxh>​
 +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)
 +\]
 +
 +<sxh groovy; gutter: true; first-line: 68>
 +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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > 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})
 +</​sxh>​
 +Given this quadratic form, one can find a propagator in the following way:
 +<sxh groovy; gutter: true; first-line: 77>
 +//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
 +</​sxh>​
 +<sxh plain; gutter: false>
 +   > [[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}]]
 +</​sxh>​
 +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$.