Bug(?) in reposition_dummy_indices: test case and patch
Hello! This simple program #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main(int argc, char** argv) { varidx mu1(symbol("mu1"), 4); varidx mu2(symbol("mu2"), 4); varidx mu3(symbol("mu3"), 4); varidx mu4(symbol("mu4"), 4); varidx mu5(symbol("mu5"), 4); varidx mu6(symbol("mu6"), 4); exvector ev2; ev2.push_back(mu3.toggle_variance()); ev2.push_back(mu6); ev2.push_back(mu5.toggle_variance()); ev2.push_back(mu6.toggle_variance()); ev2.push_back(mu5); ev2.push_back(mu3); ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2); test_cycl = test_cycl.simplify_indexed(); if (test_cycl.get_free_indices().size() != 0) { std::cerr << e << std::endl; throw std::logic_error("Inconsistent indices in expression"); } return 0; } fails (both with GiNaC 1.4 CVS and 1.3.x). I believe that the reason is bug in indexed.cpp:reposition_dummy_indices() and propose attached patch to fix it. Note that patch *seems* to be correct, but IMHO it is somewhat ugly and probably inefficient. Could anyone suggest a better solution? Best regards, Alexei. -- All science is either physics or stamp collecting.
Dear Alexei, On Tue, 29 Aug 2006, Sheplyakov Alexei wrote:
This simple program
[...]
fails (both with GiNaC 1.4 CVS and 1.3.x). I believe that the reason is bug in indexed.cpp:reposition_dummy_indices() and propose attached patch to fix it. Note that patch *seems* to be correct, but IMHO it is somewhat ugly and probably inefficient. Could anyone suggest a better solution?
What if we simply don't do this kind of toggling of the variance in the case that both indices of the dummy pair belong to the same indexed object and the indexed object has a symmetry in which the dummy pair is involved? Seems simpler to me. Best wishes, Chris
On Wed, Aug 30, 2006 at 12:08:26PM +0200, Chris Dams wrote:
This simple program
[...]
fails (both with GiNaC 1.4 CVS and 1.3.x). I believe that the reason is bug in indexed.cpp:reposition_dummy_indices() and propose attached patch to fix it. Note that patch *seems* to be correct, but IMHO it is somewhat ugly and probably inefficient. Could anyone suggest a better solution?
What if we simply don't do this kind of toggling of the variance in the case that both indices of the dummy pair belong to the same indexed object and the indexed object has a symmetry in which the dummy pair is involved? Then
T.mu.nu~mu~nu - T~nu~mu.nu.mu where T has cyclic symmetry, would not evaluate to zero. Which is (IMHO) very bad thing...
Seems simpler to me.
There exist a lot of simple *incorrect* solutions :) Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, On Wed, 30 Aug 2006, Sheplyakov Alexei wrote:
On Wed, Aug 30, 2006 at 12:08:26PM +0200, Chris Dams wrote:
What if we simply don't do this kind of toggling of the variance in the case that both indices of the dummy pair belong to the same indexed object and the indexed object has a symmetry in which the dummy pair is involved?
Then
T.mu.nu~mu~nu - T~nu~mu.nu.mu
where T has cyclic symmetry, would not evaluate to zero. Which is (IMHO) very bad thing...
Yes. That is true. However, your patch does not simplify T~mu.mu~nu.nu - T~mu~nu.nu.mu to zero if T has cyclic symmetry. Neither does the code that is in CVS at the moment. So, currently we cannot guarantee much about the simplification of this kind of expression. What about modifying the canonicalization code in symmetry.cpp, then? We could at first decide the canonical form based only on the values and dimensions of the indices. Not on their variance. If that does not give a prefered ordering we would order them the same way as is done now. Besides that, we should reject a toggling of indices in reposition_dummy indices if it results in an expression with reordered indices. Does this sound any better? Best wishes, Chris
Dear all, On Wed, 30 Aug 2006, Chris Dams wrote:
What about modifying the canonicalization code in symmetry.cpp, then? We could at first decide the canonical form based only on the values and dimensions of the indices. Not on their variance. If that does not give a prefered ordering we would order them the same way as is done now. Besides that, we should reject a toggling of indices in reposition_dummy indices if it results in an expression with reordered indices. Does this sound any better?
I now patched the function reposition_dummy_indices to fix this. The idea is that the problem Alexei pointed out can only occur if there are varidxes involved in a cyclic symmetry. In that case the function reposition_dummy_indices now starts by trying out all variations obtained by raising/lowering varidxes that occur in dummy pairs of the tensor object that is given to reposition_dummy_indices. Of these posibilities the one that is the smallest in the sense of ex_is_less is chosen. The dummy indices that have taken part in this trying out of all posibilities are removed from the vector variant_dummy_indices. After that there shouldn't be a problem anymore using the rest of the function reposition_dummy_indices as it always has been. Best wishes, Chris
Hello, On Fri, Sep 08, 2006 at 06:02:13PM +0200, Chris Dams wrote:
What about modifying the canonicalization code in symmetry.cpp, then? We could at first decide the canonical form based only on the values and dimensions of the indices. Not on their variance. If that does not give a prefered ordering we would order them the same way as is done now. Besides that, we should reject a toggling of indices in reposition_dummy indices if it results in an expression with reordered indices. Does this sound any better?
I now patched the function reposition_dummy_indices to fix this. The idea is that the problem Alexei pointed out can only occur if there are varidxes involved in a cyclic symmetry. [snipped]
I don't think your analysis is completely correct. The problem is that e = e.subs(something) calls e.eval(), and eval can change the original expression so that iterators become invalid. Actually, after the evaluation expression can be *anything*, even not necessarily indexed. For example, the expression in question could evaluate to zero (if someone implements derived class such that for some combination of indices the thing evaluates to zero). What your patch does is work-around for one of such (non-trivial eval) cases, but I think the code is still terribly broken. Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, On Fri, 6 Oct 2006, Sheplyakov Alexei wrote:
I don't think your analysis is completely correct. The problem is that
e = e.subs(something)
calls e.eval(), and eval can change the original expression so that iterators become invalid. Actually, after the evaluation expression can be *anything*, even not necessarily indexed. For example, the expression in question could evaluate to zero (if someone implements derived class such that for some combination of indices the thing evaluates to zero).
It is true that an evaluation function can be written such that this goes wrong. However, I was assuming that evaluation functions are doing things to expressions that make sense. Assuming that the metric is invertible, it is not possible that A~mu~nu~sigma~tau is nonzero and that A.mu~nu~sigma~tau is zero. Hence, an evaluation function that turns one of them into zero should do this to the other as well. Do you know of any evaluation that does make sense where we would encounter this kind of problem? I thought about this, but came to the conclusion that this probably cannot happen. On the other hand, as a matter of principle you are right that it is best to assume as little about the rest of the code as possbile. For this reason I would accept a patch that would do the substitutions as they are currently but does not cause eval to be called. Best wishes, Chris
Hello! On Tue, Oct 10, 2006 at 01:16:21PM +0200, Chris Dams wrote:
I don't think your analysis is completely correct. The problem is that
e = e.subs(something)
calls e.eval(), and eval can change the original expression so that iterators become invalid. Actually, after the evaluation expression can be *anything*, even not necessarily indexed. For example, the expression in question could evaluate to zero (if someone implements derived class such that for some combination of indices the thing evaluates to zero).
It is true that an evaluation function can be written such that this goes wrong. However, I was assuming that evaluation functions are doing things to expressions that make sense.
First of all, I think it is perfectly sensible for eval() to return zero (e.g., what about traceless tensors?). Secondly, your patch does not catch this: [I admit that example is a little bit weird] #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main(int argc, char** argv) { symbol d("d"); varidx mu(symbol("mu"), d), nu(symbol("nu"), d); ex e = indexed(symbol("T"), sy_symm(3, 0, 1, 2), mu, nu, mu.toggle_variance(), nu.toggle_variance()); cout << e << " == > "; e = e.simplify_indexed(); cout << e << endl; return 0; } On my system, this gives: T~mu.nu~nu.mu ==> T~mu.nu.nu.mu Note that we got invalid expression after the evaluation. [snipped]
Do you know of any evaluation that does make sense where we would encounter this kind of problem? I thought about this, but came to the conclusion that this probably cannot happen.
I was trying to implement two classes (derived from indexed) which represent Riemann and Ricci tensors. I did R~k.i.k.j => r.i.j substitution in R::eval(). But my code did not work. After a while I found a way to reproduce the problem with classes which are part of GiNaC.
On the other hand, as a matter of principle you are right that it is best to assume as little about the rest of the code as possbile. For this reason I would accept a patch that would do the substitutions as they are currently but does not cause eval to be called.
The patch I've posted does that. Is there something wrong/stupid/etc with it? Best regards, Alexei. -- All science is either physics or stamp collecting.
participants (2)
-
Chris Dams
-
varg＠theor.jinr.ru