Re: Bug(?) in reposition_dummy_indices: test case
Dear Alexei,
First of all, I think it is perfectly sensible for eval() to return zero (e.g., what about traceless tensors?).
Yes. However, if we have a traceless tensor both T~mu.mu and T.mu~mu should evaluated themselves into 0. It does not make sense to only have one of them evaluate into zero. The point I was trying to make is that tensors with some index lowered/raised always (?) have the same properties so that an evaluate step that aplies to the version with an up-index should also aply to the one with a down-index.
Secondly, your patch does not catch this: [I admit that example is a little bit weird]
Your example is perfectly fine. I should have deleted all indices that are dummy indices of the indexed object e that is given as a parameter to reposition_dummy_indices in any case, not just when a cyclic symmetry is involved. Because you also mentioned the possibility that users would want to implement things like the Riemann curvature tensor, it makes sense that the evaluation of an indexed object could do any weird kind of symmetrization it would want, so I also removed the condition that a cyclic symmetry should be involved when trying all combinations of repositioning local dummy indices.
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?
Oh well, the thing seemed rather complicated to me and didn't do all simplifications possible, as noted in http://www.ginac.de/pipermail/ginac-devel/2006-August/001055.html. I now more or less applied it, but with some modifications. I hope that I haven't broken anything ;-). (1) Now that the raising and lowering of indices is done on an exvector, we do not need to raise/lower both indices at the same time. Another reason why this isn't necessary anymore is that indices that occur in dummy pairs within the indexed object e should already have been removed from the vector variant_dummy_indices. I think this means that the for_each that you had can be omitted. (2) It is not correct to return the indexed object by doing indexed ei_ = ex_to<indexed>(e); ei_.seq = seq; e = ei_; This should be e = ex_to<indexed>(e).thiscontainer(seq); because we may be dealing with something that inherits from indexed rather than an indexed object. I also simplified the clifford exam by removing a huge macro. A patch is in CVS. Best wishes, Chris
Hi, Chris! On Fri, Oct 20, 2006 at 10:51:21AM +0200, Chris Dams wrote:
Yes. However, if we have a traceless tensor both T~mu.mu and T.mu~mu should evaluated themselves into 0. It does not make sense to only have one of them evaluate into zero. The point I was trying to make is that tensors with some index lowered/raised always (?) have the same properties so that an evaluate step that aplies to the version with an up-index should also aply to the one with a down-index.
I think it *does* make sense to replace only one of combinations and rely on simplify_indexed to get indices canonicalized. Otherwise one need to duplicate quite a large portion of code from indexed.cpp and symmetry.cpp (to decide what combinations are equivalent). [snipped]
Oh well, the thing seemed rather complicated to me and didn't do all simplifications possible, as noted in http://www.ginac.de/pipermail/ginac-devel/2006-August/001055.html. I now more or less applied it, but with some modifications. I hope that I haven't broken anything ;-). (1) Now that the raising and lowering of indices is done on an exvector, we do not need to raise/lower both indices at the same time. Another reason why this isn't necessary anymore is that indices that occur in dummy pairs within the indexed object e should already have been removed from the vector variant_dummy_indices. I think this means that the for_each that you had can be omitted.
I agree.
(2) It is not correct to return the indexed object by doing
indexed ei_ = ex_to<indexed>(e); ei_.seq = seq; e = ei_;
This should be e = ex_to<indexed>(e).thiscontainer(seq); because we may be dealing with something that inherits from indexed rather than an indexed object.
I didn't like that chunk myself exactly for this reason, but I didn't know the correct way.
I also simplified the clifford exam by removing a huge macro.
A patch is in CVS. Works fine for me. I've re-diffed it for 1.3 branch, could you please apply it?
Best regards, Alexei. [PATCH] reposition_dummy_indices: fix bugs w.r.t. raising/lowering dummy indices. Now reposition_dummy_indices works fine if the argument has non-trivial symmetry properties and/or evaluation method (which triggers re-evaluation in course of exchanging of dummy indices). Thus, simplify_indexed gives correct results for such expressions. Backported from the main branch. --- check/exam_paranoia.cpp | 28 ++++++++++++++++ ginac/indexed.cpp | 82 ++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 97 insertions(+), 13 deletions(-) diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp index 7c82e48..0ab6491 100644 --- a/check/exam_paranoia.cpp +++ b/check/exam_paranoia.cpp @@ -437,6 +437,33 @@ static unsigned exam_paranoia16() return result; } +// Bug in reposition_dummy_indices() could result in correct expression +// turned into one with inconsistent indices. Fixed on Aug 29, 2006 +static unsigned exam_paranoia17() +{ + 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); + // notice: all indices are contracted ... + + ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2); + test_cycl = test_cycl.simplify_indexed(); + // ... so there should be zero free indices in the end. + return test_cycl.get_free_indices().size(); +} + + unsigned exam_paranoia() { unsigned result = 0; @@ -460,6 +487,7 @@ unsigned exam_paranoia() result += exam_paranoia14(); cout << '.' << flush; result += exam_paranoia15(); cout << '.' << flush; result += exam_paranoia16(); cout << '.' << flush; + result += exam_paranoia17(); cout << '.' << flush; if (!result) { cout << " passed " << endl; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 05c3b90..2de726b 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -632,10 +632,60 @@ bool reposition_dummy_indices(ex & e, ex { bool something_changed = false; + // Find dummy symbols that occur twice in the same indexed object. + exvector local_var_dummies; + local_var_dummies.reserve(e.nops()/2); + for (size_t i=1; i<e.nops(); ++i) { + if (!is_a<varidx>(e.op(i))) + continue; + for (size_t j=i+1; j<e.nops(); ++j) { + if (is_dummy_pair(e.op(i), e.op(j))) { + local_var_dummies.push_back(e.op(i)); + for (exvector::iterator k = variant_dummy_indices.begin(); + k!=variant_dummy_indices.end(); ++k) { + if (e.op(i).op(0) == k->op(0)) { + variant_dummy_indices.erase(k); + break; + } + } + break; + } + } + } + + // In the case where a dummy symbol occurs twice in the same indexed object + // we try all posibilities of raising/lowering and keep the least one in + // the sense of ex_is_less. + ex optimal_e = e; + size_t numpossibs = 1 << local_var_dummies.size(); + for (size_t i=0; i<numpossibs; ++i) { + ex try_e = e; + for (size_t j=0; j<local_var_dummies.size(); ++j) { + exmap m; + if (1<<j & i) { + ex curr_idx = local_var_dummies[j]; + ex curr_toggle = ex_to<varidx>(curr_idx).toggle_variance(); + m[curr_idx] = curr_toggle; + m[curr_toggle] = curr_idx; + } + try_e = e.subs(m, subs_options::no_pattern); + } + if(ex_is_less()(try_e, optimal_e)) + { optimal_e = try_e; + something_changed = true; + } + } + e = optimal_e; + + if (!is_a<indexed>(e)) + return true; + + exvector seq = ex_to<indexed>(e).seq; + // If a dummy index is encountered for the first time in the // product, pull it up, otherwise, pull it down - exvector::const_iterator it2, it2start, it2end; - for (it2start = ex_to<indexed>(e).seq.begin(), it2end = ex_to<indexed>(e).seq.end(), it2 = it2start + 1; it2 != it2end; ++it2) { + for (exvector::iterator it2 = seq.begin()+1, it2end = seq.end(); + it2 != it2end; ++it2) { if (!is_exactly_a<varidx>(*it2)) continue; @@ -643,14 +693,20 @@ bool reposition_dummy_indices(ex & e, ex for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to<varidx>(*it2).is_covariant()) { - e = e.subs(lst( - *it2 == ex_to<varidx>(*it2).toggle_variance(), - ex_to<varidx>(*it2).toggle_variance() == *it2 - ), subs_options::no_pattern); + /* + * N.B. we don't want to use + * + * e = e.subs(lst( + * *it2 == ex_to<varidx>(*it2).toggle_variance(), + * ex_to<varidx>(*it2).toggle_variance() == *it2 + * ), subs_options::no_pattern); + * + * since this can trigger non-trivial repositioning of indices, + * e.g. due to non-trivial symmetry properties of e, thus + * invalidating iterators + */ + *it2 = ex_to<varidx>(*it2).toggle_variance(); something_changed = true; - it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start); - it2start = ex_to<indexed>(e).seq.begin(); - it2end = ex_to<indexed>(e).seq.end(); } moved_indices.push_back(*vit); variant_dummy_indices.erase(vit); @@ -661,11 +717,8 @@ bool reposition_dummy_indices(ex & e, ex for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to<varidx>(*it2).is_contravariant()) { - e = e.subs(*it2 == ex_to<varidx>(*it2).toggle_variance(), subs_options::no_pattern); + *it2 = ex_to<varidx>(*it2).toggle_variance(); something_changed = true; - it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start); - it2start = ex_to<indexed>(e).seq.begin(); - it2end = ex_to<indexed>(e).seq.end(); } goto next_index; } @@ -674,6 +727,9 @@ bool reposition_dummy_indices(ex & e, ex next_index: ; } + if (something_changed) + e = ex_to<indexed>(e).thiscontainer(seq); + return something_changed; } -- 1.4.2.3 -- All science is either physics or stamp collecting.
participants (2)
-
Chris Dams
-
varg＠theor.jinr.ru