Simplifying indexed expressions not working on multi-indices scalar_product
According to section 4.14.5 "Simplifying indexed expressions" of the tutorial, if I set: scalar_products sp; sp.add(A, B, 0); then: A~mu~nu*B.mu.nu should be replaced by zero when using .simplify_indexed(sp). This is not working. See my simple code below. Is this a bug?, any help will be appreciated Alejandro PS. I am using Ginac version 1.3.4 on debian this are the rpm I've got: ginac-tools_1.3.4-1_i386.deb libginac1.3c2a_1.3.4-1_i386.deb libginac-dev_1.3.4-1_i386.deb %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main (void) { symbol A("A"),B("B"),mu_sym("mu"),nu_sym("nu"); idx mu(mu_sym,3), nu(nu_sym,3); ex e, e1; scalar_products sp; sp.add(A,B,0); cout << "scalar_prod 1 index:" << endl; e= indexed(A,mu)*indexed(B,mu); cout <<e <<endl; cout << e.simplify_indexed(sp) << " ==> correct " << endl; cout << "scalar_prod 2 indices:" << endl; e1= indexed(A,mu,nu)*indexed(B,mu,nu); cout << e1 <<endl; cout << e1.simplify_indexed(sp) << " ==> wrong, it is not zero" << endl; return 0; }
Dear Alejandro, On Wed, 26 Jul 2006, Alejandro Limache wrote:
According to section 4.14.5 "Simplifying indexed expressions" of the tutorial, if I set: scalar_products sp; sp.add(A, B, 0); then: A~mu~nu*B.mu.nu should be replaced by zero when using .simplify_indexed(sp). This is not working. See my simple code below.
Ah yes, I see. Actually, I removed this possibility from the library some time ago. See: http://www.ginac.de/pipermail/ginac-devel/2005-December/000878.html Since I thought that nobody would be using this feature, I, in the end, did what I said in the "Note:" and now GiNaC does not try anymore to apply scalar_products in case that the indexed objects have more than one index. Of course, the tutorial should also be made consistent with this, so I will do that. On the other hand, you might also try to convince me that we really do need this kind of contraction... Thanks for pointing this out, anyway. Best wishes, Chris
Chris, Thank you very much for your prompt reply. On the issue if scalar_products defined for more than one indices are useful, I would say definitely yes. Scalar products also known as "inner products" appear in many areas of physics and they may involve higher order tensors (tensor with more than one index). For example, in the area of continuum mechanics, its very common to see inner products defined for 2nd order tensors (2 indices). In particular the "stress power P" is the inner product of the stress tensor T with the rate of deformation tensor D: P = T*D=T.i.j * D.i.j and measures the power produced by the inner (stress) forces. Three observations on inner products: 1) By definition an inner product must be a scalar so all tensor indices should get contracted. 2) One must be careful with the order of the indices. The inner product of a non zero tensor with itself must be positive, from this condition it follows that T.i.j * D.j.i is NOT an inner product 3) Covariant indices should contract with contravariant indices so T~a.mu*T~a~mu is not correct a should throw an exception or not get simplified, but T.a.mu*T~a~mu should get simplified. My first impression on GINAC has been very good. I am really hopeful to be able to make a unique code where I can mix symbolic calculations with numerical computations, that would be so great for my research. One of the things I like to do for example is to use symbolic math to compute quantities defined as C++ objects such as A = B+C+D, so they get expanded in indices A.i = B.i+C.i.D.i to be computed efficiently as they were C languaje. I hope GINAC can be helpful me on this kind task. Alejandro
From: Chris Dams <Chris.Dams@mi.infn.it> Reply-To: GiNaC discussion list <ginac-list@ginac.de> To: GiNaC discussion list <ginac-list@ginac.de> Subject: Re: [GiNaC-list] Simplifying indexed expressions not working onmulti-indices scalar_product Date: Thu, 27 Jul 2006 13:08:59 +0200 (CEST)
Dear Alejandro,
On Wed, 26 Jul 2006, Alejandro Limache wrote:
According to section 4.14.5 "Simplifying indexed expressions" of the tutorial, if I set: scalar_products sp; sp.add(A, B, 0); then: A~mu~nu*B.mu.nu should be replaced by zero when using .simplify_indexed(sp). This is not working. See my simple code below.
Ah yes, I see. Actually, I removed this possibility from the library some time ago. See: http://www.ginac.de/pipermail/ginac-devel/2005-December/000878.html Since I thought that nobody would be using this feature, I, in the end, did what I said in the "Note:" and now GiNaC does not try anymore to apply scalar_products in case that the indexed objects have more than one index. Of course, the tutorial should also be made consistent with this, so I will do that. On the other hand, you might also try to convince me that we really do need this kind of contraction...
Thanks for pointing this out, anyway.
Best wishes, Chris
_______________________________________________ GiNaC-list mailing list GiNaC-list@ginac.de https://www.cebix.net/mailman/listinfo/ginac-list
Dear Alejandro, On Thu, 27 Jul 2006, Alejandro Limache wrote:
For example, in the area of continuum mechanics, its very common to see inner products defined for 2nd order tensors (2 indices). In particular the "stress power P" is the inner product of the stress tensor T with the rate of deformation tensor D: P = T*D=T.i.j * D.i.j and measures the power produced by the inner (stress) forces.
Yes, such a thing could be useful. However....
1) By definition an inner product must be a scalar so all tensor indices should get contracted.
Yes, obviously.
2) One must be careful with the order of the indices. The inner product of a non zero tensor with itself must be positive, from this condition it follows that T.i.j * D.j.i is NOT an inner product
Yes, the old implementation that I deleted would do this. Looking at the old code again, I guess it would even do T.i.j.k.k*D.i.j -> P if the user were to define that the scalar product between T and D is P. *shrug* I can think of two problems with your solution (1) What if the indices are of different dimension? Should it be possible to specify a different dimension for every index postion? IMO, this would get a bit overly complicated. Maybe user-defined inner products should only be used if all dimensions of indices involved are the same. This, at least, seems to cover your case. (2) What if T happens to be symmetric and D is not? Then the product T.i.j*D.i.j might automatically evaluate into T.j.i*D.i.j and the contraction would no longer occur, depending on what the canonical ordering of i and j happens to be. Not good. Maybe looping over all posibilities for a symmetric tensor would need to be done. I'm not sure whether I will have time to implement this, though.
3) Covariant indices should contract with contravariant indices so T~a.mu*T~a~mu is not correct a should throw an exception or not get simplified, but T.a.mu*T~a~mu should get simplified.
Oh yes, of course I should have written T.a.mu*T.a~mu, where a was not a varidx but an ordinary idx. Never mind though, the problem was that if these indices a and mu were of incomparable dimension (in the sense of the function minimal_dim) an exception would be thrown. Best wishes, Chris
Chris,
(1) What if the indices are of different dimension? Should it be possible to specify a different dimension for every index postion? IMO, this would get a bit overly complicated. Maybe user-defined inner products should only be used if all dimensions of indices involved are the same. This, at least, seems to cover your case.
The cases of inner products of physical quantities that I can recall always involve indices of the same dimension, but maybe in some general optimization or minimization problems one might have different dimensions. These inner products usually measure some kind of energy or power. In the case of FEM (Finite Element Method) there are times where you can get products M.k.i.j*A.i.j of multidimensional arrays M,A where for example i and j have dimension of space (3) and k the dimension of the number of nodes or elements. However in this case there is no contraction of the k index.
(2) What if T happens to be symmetric and D is not? Then the product T.i.j*D.i.j might automatically evaluate into T.j.i*D.i.j and the contraction would no longer occur, depending on what the canonical ordering of i and j happens to be. Not good. Maybe looping over all posibilities for a symmetric tensor would need to be done.
Yes that is true and alas...T is most of the time symmetric in the example I gave. Maybe looping could be a good solution
I'm not sure whether I will have time to implement this, though.
Chris, I even could not say that all these things are extremely important. I my case most of the time I am interested in computing products such as T.i.j*D.i.j numerically. But I would like to be able to make use of symbolic math in order to: 1) to get some simplifications so as to have less terms to compute numerically 2) to expand tensor expressions (combinations of sums and products) in indicial form (as I think is done in some libraries that use expression templates like blitz++) in order to reduce the loss of performance of C++ with respect to fortran o C code. 3) to do some analytical differentiation of terms I think these types of things could make a code more general, powerful and efficient (faster). I hope Ginac could help me on that. Thank you, Alejandro
--- Chris Dams <Chris.Dams@mi.infn.it> wrote: [...]
On Thu, 27 Jul 2006, Alejandro Limache wrote: [...]
1) By definition an inner product must be a scalar so all tensor indices should get contracted.
Yes, obviously.
Hmmm..., seems to me that in general you would also like to have scalar products that may be dependent on other variables, including integers. --- Alejandro Limache <alejandrolimache@hotmail.com> wrote:
[--- Chris Dams <Chris.Dams@mi.infn.it> wrote:]
(1) What if the indices are of different dimension? [[...]] Maybe user-defined inner products should only be used if all dimensions of indices involved are the same. [[...]] [...] one might have different dimensions. [...] [...] there are times where you can get products M.k.i.j*A.i.j of multidimensional arrays M,A where for example i and j have dimension of space (3) and k the dimension of the number of nodes or elements. [...]
(M.k.i.j*A.i.j){k=1...N} here seems to be such an example of an *indexed scalar*. I also vaguely seem to recall seeing some treatments in general relativity where some indices of indexed quantities run over "space" dimensions only, whereas other indices run over all four space-time dimensions. My memory is hazy here, so I suppose I could be a bit mistaken on this. There is also some mumbo-jumbo concerning what indexed entities constitute a "tensor" -- something about the indices relating to coordinates and the indexed entities conforming to certain rules of transformation -- so that, I suppose, the tensors constitute a kind of equivalence class with such an indexed entity being a representative of that class. So it seems that one would like to be able to distinguish indices in various convenient ways; some indices would relate to the "tensor" properties and some would (perhaps) be extraneous to the concept of a tensor -- perhaps relevant to some other type of algebraic concept. And then there are some who like to work with "coordinate-free" representations of tensors. I'm a bit of a newbie to GiNaC; so I don't know yet how relevant all of this is to GiNaC. Richard Haney __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com
Hello people, I have the following question. Is there a way to SET and GET numeric values from indexed expressions? Consider the 2 following cases: 1) For matrices. In the tutorial Section 4.13 one sees that is possible to do: matrix A(2, 2), B(2, 2), C(2, 2); A = 1, 2, 3, 4; B = -1, 0, 2, 1; C = 8, 4, 2, 1; e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j); So: cout << e.simplify_indexed() << endl; // -> [[-13,-6],[1,2]].i.j Is there a way to get the numeric values of each of the components of e? I would like to be able to do something like double d; d = e(0,0); to get: d=-13 2) Can something similar be done for indexed expression not defined as matrices? For example given something like symbol A("A").... idx i(.., 3), j(.., 3), k(..,10); ex e = indexed(A,i,j,k) could I set for example: i=2, j=0, k=4 and set A(2,0,4)=3.5 ? could I then somehow get the value "3.5" back? Thank you for your help! Alejandro
Dear Alejandro, On Thu, 27 Jul 2006, Alejandro Limache wrote:
// -> [[-13,-6],[1,2]].i.j Is there a way to get the numeric values of each of the components of e? I would like to be able to do something like double d; d = e(0,0); to get: d=-13
ex_to<numeric>(e.subs(lst(i==0,j==0))).to_double() You should only do that if you are certain that the relevant matrix entry is numeric or all hell will break loose. Safer is to first do a test: if (is_a<numeric>(....))
2) Can something similar be done for indexed expression not defined as matrices? For example given something like symbol A("A").... idx i(.., 3), j(.., 3), k(..,10); ex e = indexed(A,i,j,k) could I set for example: i=2, j=0, k=4 and set A(2,0,4)=3.5 ? could I then somehow get the value "3.5" back?
Not directly. The closest that you can get to that is keeping an exmap m and say m[A(2,0,4)] = 3.5 and then do .subs(m) in an expression where you have unrolled the dummy indices. With the dimensions you give above this could get a costly though. If you don't have any patterns in the exmap you could use the option subs_options::no_pattern to speed it up a quite bit. Best wishes, Chris
Chris, Thank you for helping me to get a grasp on GINAC.
Do you mean the writing out of contractions as in A.i*B.i -> A.1*B.1 + ... + A.3*B.3. That can be done using the function ex expand_dummy_sum(const ex & e, bool subs_idx). Yes that is what I wanted to do. I did an example to do just that in the simple code below however instead of getting: A.i*B.i -> A.0*B.0 + A.1*B.1 I am getting: A.i*B.i -> A.0*B.i + A.1*B.i so why is it not working?. Does it has something to do with the flag: bool subs_idx I just put a "1" in there. What is the purpose of that flag?
Not directly. The closest that you can get to that is keeping an exmap m and say m[A(2,0,4)] = 3.5 and then do .subs(m) in an expression where you have unrolled the dummy indices. I tried to do this in the simple code below where I have the expanded expression: A.0*B.i + A.1*B.i I tried to replace A.1 by 7.5. It worked ok but I noticed that in order for this to work I needed to use an expression of the form: indexed(A,idx(1,)) in .subs Is there any alternative where I can just put some kind of string str "A.1" inside .subs and still work?
Alejandro %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #include <iostream> #include <ginac/ginac.h> #include <ginac/flags.h> using namespace std; using namespace GiNaC; int main (void) { symbol i_sym("i"); idx i(i_sym, 2); symbol A("A"),B("B"); ex exA = indexed(A,i); ex exB = indexed(B,i); ex exAB = exA*exB; cout << exAB << endl; ex exABexp = expand_dummy_sum(exAB,1); cout << exABexp << endl; // replacing numeric values for component A.1 ex exA1 = indexed(A,idx(1,2)); cout << exABexp.subs(exA1 == 7.5) << endl; return 0; }
Dear Alejandro, On Fri, 28 Jul 2006, Alejandro Limache wrote:
ex expand_dummy_sum(const ex & e, bool subs_idx).
Yes that is what I wanted to do. I did an example to do just that in the simple code below however instead of getting: A.i*B.i -> A.0*B.0 + A.1*B.1 I am getting: A.i*B.i -> A.0*B.i + A.1*B.i so why is it not working?.
Yes, this function is broken. The commentary in indexed.h states that it was intended to be used for all indices but actually it only works for varidxes. I will fix this. In a few days the CVS should have code that does work as expected. If you don't want to use CVS-code or are in a hurry, I suggest you write your own function to do this.
Not directly. The closest that you can get to that is keeping an exmap m and say m[A(2,0,4)] = 3.5 and then do .subs(m) in an expression where you have unrolled the dummy indices. I tried to do this in the simple code below where I have the expanded expression: A.0*B.i + A.1*B.i I tried to replace A.1 by 7.5. It worked ok but I noticed that in order for this to work I needed to use an expression of the form: indexed(A,idx(1,)) in .subs Is there any alternative where I can just put some kind of string str "A.1" inside .subs and still work?
That is not possible. You will have to use indexed(A,idx(1,2)). Best wishes, Chris
Dear Alejandro, On Thu, 27 Jul 2006, Alejandro Limache wrote:
The cases of inner products of physical quantities that I can recall always involve indices of the same dimension, but maybe in some general optimization or minimization problems one might have different dimensions.
It can quite easily get more hairy then that. In the standard model before symmetry breaking the left-handed quark field has four indices that all live in different spaces. In dimensional regularization there can be indices that live in the same space, yet have a different dimension. The latter case is the main reason that we really do care about the dimension.
1) to get some simplifications so as to have less terms to compute numerically
You can always use algebraic substitution: .subs( indexed(T, idx(wild(0),3), idx(wild(1),3)) * indexed(D, idx(wild(0),3), idx(wild(1),3)) == P, subs_options::subs_algebraic )
2) to expand tensor expressions (combinations of sums and products) in indicial form (as I think is done in some libraries that use expression templates like blitz++) in order to reduce the loss of performance of C++ with respect to fortran o C code.
Do you mean the writing out of contractions as in A.i*B.i -> A.1*B.1 + ... + A.3*B.3. That can be done using the function ex expand_dummy_sum(const ex & e, bool subs_idx).
3) to do some analytical differentiation of terms I think these types of things could make a code more general, powerful and efficient (faster).
Do you mean differentiation with respect to indexed objects? Yes, that certainly is very useful. It is not present in GiNaC at the moment. It is not so difficult, though, to write your own function to do this. Just write a function for it that uses a type switch. Below an example is given of some code that I once wrote for this purpose. It is probably more complicated then what you would need. Best wishes, Chris ex int_deriv(const ex&x,const ex&y) { if(is_a<indexed>(x)&&is_a<indexed>(y)&&x.nops()==y.nops()) { ex result=int_deriv(x.op(0),y.op(0)); if(result==0) return 0; for(int i=1;i<x.nops();++i) { if(is_a<varidx>(x.op(i))&&is_a<varidx>(y.op(i))) { ex yother=ex_to<varidx>(y.op(i)).toggle_variance(); result*=lorentz_g(x.op(i),yother); } else result*=delta_tensor(x.op(i),y.op(i)); } return result; } else if(is_a<function>(x)&&x.nops()==1) { ex deriv,repl; if(is_a<indexed>(y)) { deriv=y.op(0); repl=y.op(1); } else if(is_a<lst>(y)&&y.nops()==2) { deriv=y.op(1); repl=y.op(0); } else { deriv=y; repl=1; } if(is_a<function>(deriv)&&is_a<function>(x)) { unsigned ser=ex_to<function>(deriv).get_serial(); if(ser!=ex_to<function>(x).get_serial()) return 0; return x.subs(x==repl)*nullify(x.op(0)==deriv.op(0)); } } else if(is_a<mul>(x)||is_a<ncmul>(x)||is_a<lst>(x)|| x.match(inp(wild(0),wild(1)))) { ex result=0; for(int i=0;i<x.nops();++i) { ex factor=int_deriv(x.op(i),y); if(factor!=0) result+=x.subs(x.op(i)==factor).expand(); } return result; } else if(is_a<add>(x)) { exvector result; for(int i=0;i<x.nops();++i) { if(!(i%1000)&&x.nops()>1000) cerr << i << "/" << x.nops() << endl; result.push_back(int_deriv(x.op(i),y)); } return add(result); } return 0; }
participants (3)
-
Alejandro Limache
-
Chris Dams
-
Richard Haney