Skip to content

Commit e0fdf55

Browse files
kpeeterskpeeters
committed
Merge branch 'cbehan-patch5'
2 parents 039991d + 7b414a0 commit e0fdf55

File tree

2 files changed

+109
-30
lines changed

2 files changed

+109
-30
lines changed

‎core/algorithms/combine.cc‎

Lines changed: 92 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ Algorithm::result_t combine::apply(iterator& it)
2020
{
2121
sibling_iterator sib=tr.begin(it);
2222
index_map_t ind_free, ind_dummy;
23+
std::vector<Ex::iterator> dummies;
2324
while(sib!=tr.end(it)) { // iterate over all factors in the product
2425
sibling_iterator ch=tr.begin(sib);
2526
while(ch!=tr.end(sib)) { // iterate over all indices of this factor
@@ -33,72 +34,135 @@ Algorithm::result_t combine::apply(iterator& it)
3334
}
3435
if(ind_dummy.size()==0) return result_t::l_no_action;
3536

36-
index_map_t::iterator dums1=ind_dummy.begin(), dums2;
37-
while(dums1!=ind_dummy.end()) {
37+
while(ind_dummy.begin()!=ind_dummy.end()) {
38+
bool found=false;
39+
index_map_t::iterator start=ind_dummy.begin();
40+
while(!found && start!=ind_dummy.end()) {
41+
iterator parent=tr.parent(start->second);
42+
sibling_iterator ch=tr.begin(parent), last_part;
43+
while(ch!=tr.end(parent)) {
44+
auto fnd=ind_dummy.find((Ex::iterator)ch);
45+
if(fnd!=ind_dummy.end()) last_part=ch;
46+
++ch;
47+
}
48+
if(last_part==start->second) {
49+
// We are on a rightmost contracted index
50+
found=true;
51+
}
52+
else ++start;
53+
}
54+
bool paired=true;
55+
while(paired && start!=ind_dummy.end()) {
56+
iterator parent=tr.parent(start->second);
57+
sibling_iterator ch=tr.begin(parent), last_part;
58+
while(ch!=tr.end(parent) && ind_dummy.size()>0) {
59+
auto fnd2=ind_dummy.equal_range((Ex::iterator)ch);
60+
auto fnd1=fnd2.first;
61+
if(fnd1->second!=ch) ++fnd1;
62+
if(fnd1->second==ch) {
63+
dummies.insert(dummies.end(), ch);
64+
ind_dummy.erase(fnd1);
65+
last_part=ch;
66+
}
67+
++ch;
68+
}
69+
auto fnd=ind_dummy.find((Ex::iterator)last_part);
70+
if(fnd==ind_dummy.end()) {
71+
last_part->flip_parent_rel();
72+
fnd=ind_dummy.find((Ex::iterator)last_part);
73+
last_part->flip_parent_rel();
74+
}
75+
if(fnd==ind_dummy.end()) {
76+
// Contraction ends because we are on a vector
77+
// It could also be a trace if we removed the paired index more than one iteration ago
78+
paired=false;
79+
}
80+
else {
81+
start=fnd;
82+
index_map_t::iterator check=ind_dummy.end();
83+
iterator parent=tr.parent(start->second);
84+
sibling_iterator ch=tr.begin(parent), first_part;
85+
while(check==ind_dummy.end()) {
86+
first_part=ch;
87+
check=ind_dummy.find((Ex::iterator)ch);
88+
++ch;
89+
}
90+
if(first_part!=start->second) {
91+
throw NotYetImplemented("Evaluation requires transposing a matrix.");
92+
return result_t::l_no_action;
93+
}
94+
}
95+
}
96+
}
97+
98+
std::vector<Ex::iterator>::iterator dums1=dummies.begin(), dums2;
99+
dums2=dums1;
100+
++dums2;
101+
while(dums1!=dummies.end() && dums2!=dummies.end()) {
38102
// txtout << "analysing " << std::endl;
39103
// txtout << *(dums1->second->name) << std::endl;
40-
dums2=dums1;
41-
++dums2;
42104

43105
bool isbrack1=false, isbrack2=false;
44106
bool ismatorvec1=false, ismatorvec2=false;
45-
bool diffparents=tr.parent(dums1->second)!=tr.parent(dums2->second);
46-
const Matrix *mat1=kernel.properties.get<Matrix>(tr.parent(dums1->second));
107+
// These are both to recognize traces
108+
bool diffparents=tr.parent(*dums1)!=tr.parent(*dums2);
109+
bool consecutive=*(*dums1)->name==*(*dums2)->name;
110+
const Matrix *mat1=kernel.properties.get<Matrix>(tr.parent(*dums1));
47111
if(mat1)
48112
ismatorvec1=true;
49-
else if(*(tr.parent(dums1->second)->name)=="\\indexbracket") {
113+
else if(*(tr.parent(*dums1)->name)=="\\indexbracket") {
50114
ismatorvec1=true;
51115
isbrack1=true;
52116
}
53-
else if(tr.number_of_children(tr.parent(dums1->second))==1)
117+
else if(tr.number_of_children(tr.parent(*dums1))==1)
54118
ismatorvec1=true;
55-
const Matrix *mat2=kernel.properties.get<Matrix>(tr.parent(dums2->second));
119+
const Matrix *mat2=kernel.properties.get<Matrix>(tr.parent(*dums2));
56120
if(mat2)
57121
ismatorvec2=true;
58-
else if(*(tr.parent(dums2->second)->name)=="\\indexbracket") {
122+
else if(*(tr.parent(*dums2)->name)=="\\indexbracket") {
59123
ismatorvec2=true;
60124
isbrack2=true;
61125
}
62-
else if(tr.number_of_children(tr.parent(dums2->second))==1)
126+
else if(tr.number_of_children(tr.parent(*dums2))==1)
63127
ismatorvec2=true;
64128

65-
if(ismatorvec1 && ismatorvec2 && diffparents) {
129+
if(ismatorvec1 && ismatorvec2 && diffparents && consecutive) {
66130
// txtout << "gluing " << *(dums2->second->name) << std::endl;
67131
// create new indexbracket with product node
68-
iterator outerbrack=tr.insert(tr.parent(dums1->second), str_node("\\indexbracket"));
132+
iterator outerbrack=tr.insert(tr.parent(*dums1), str_node("\\indexbracket"));
69133
iterator brackprod=tr.append_child(outerbrack, str_node("\\prod"));
70-
iterator parn1=tr.parent(dums1->second);
71-
iterator parn2=tr.parent(dums2->second);
134+
iterator parn1=tr.parent(*dums1);
135+
iterator parn2=tr.parent(*dums2);
72136
// remove the dummy index from these two objects, and move
73137
// other (dummy or not) indices to the outer indexbracket.
74-
sibling_iterator ind1=tr.begin(tr.parent(dums1->second));
75-
sibling_iterator stop1=tr.end(tr.parent(dums1->second));
138+
sibling_iterator ind1=tr.begin(tr.parent(*dums1));
139+
sibling_iterator stop1=tr.end(tr.parent(*dums1));
76140
if(isbrack1)
77141
++ind1;
78142
while(ind1!=stop1) {
79-
if(ind1!=dums1->second) {
143+
if(ind1!=*dums1) {
80144
sibling_iterator nxt=ind1;
81145
++nxt;
82146
tr.reparent(outerbrack, ind1, nxt);
83147
}
84148
++ind1;
85149
// ind1=tr.erase(ind1);
86150
}
87-
tr.erase(dums1->second);
88-
sibling_iterator ind2=tr.begin(tr.parent(dums2->second));
89-
sibling_iterator stop2=tr.end(tr.parent(dums2->second));
151+
tr.erase(*dums1);
152+
sibling_iterator ind2=tr.begin(tr.parent(*dums2));
153+
sibling_iterator stop2=tr.end(tr.parent(*dums2));
90154
if(isbrack2)
91155
++ind2;
92156
while(ind2!=stop2) {
93-
if(ind2!=dums2->second) {
157+
if(ind2!=*dums2) {
94158
sibling_iterator nxt=ind2;
95159
++nxt;
96160
tr.reparent(outerbrack, ind2, nxt);
97161
}
98162
++ind2;
99163
// ind2=tr.erase(ind2);
100164
}
101-
tr.erase(dums2->second);
165+
tr.erase(*dums2);
102166

103167
// put both objects inside the indexbracket.
104168
if(isbrack1) {
@@ -130,8 +194,12 @@ Algorithm::result_t combine::apply(iterator& it)
130194
tr.reparent(brackprod,parn2,nxt);
131195
}
132196
}
197+
if(consecutive) {
198+
++dums1;
199+
++dums2;
200+
}
133201
++dums1;
134-
++dums1;
202+
++dums2;
135203
}
136204

137205
//std::cerr << it << std::endl;

‎tests/implicit.cdb‎

Lines changed: 17 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ def test02():
2121
A::ImplicitIndex(A_{a b});
2222
B::ImplicitIndex(B_{a b});
2323
ex:= Tr( A B A C B ) - Tr( C );
24-
untrace(_);
24+
untrace(_)
2525
tst:= C Tr( A B A B ) - C Tr(1) - @(ex);
2626
assert(tst==0)
2727
print("Test 02 passed")
@@ -34,7 +34,7 @@ def test03():
3434
\sigma^{\mu}::ImplicitIndex(\sigma^{\mu a}_{b});
3535
Tr{#}::Trace(indices=spinor);
3636
ex:= Tr( \sigma^{\mu} q \sigma^{\nu} );
37-
untrace(_);
37+
untrace(_)
3838
tst:= q Tr( \sigma^{\mu} \sigma^{\nu} ) - @(ex);
3939
assert(tst==0)
4040
print("Test 03 passed")
@@ -54,7 +54,7 @@ def test05():
5454
__cdbkernel__=create_scope()
5555
{a,b,c,d}::Indices(vector);
5656
ex:=(A)_{a b} (B)_{b c} (C)_{c d};
57-
combine(_);
57+
combine(_)
5858
tst:= (A B C)_{a d} - @(ex);
5959
assert(tst==0)
6060
print("Test 05 passed")
@@ -65,8 +65,8 @@ def test06():
6565
__cdbkernel__=create_scope()
6666
{a,b}::Indices(vector);
6767
ex:=(A)_{a b} (B)_{b a};
68-
combine(_);
69-
tst:= (A B)_{b b} - @(ex);
68+
combine(_)
69+
tst:= (B A)_{b b} - @(ex);
7070
assert(tst==0)
7171
print("Test 06 passed")
7272

@@ -79,10 +79,21 @@ def test07():
7979
\Gamma^{m}::ImplicitIndex(\Gamma^{m a}_{b});
8080
\Gamma^{m}::Matrix;
8181
ex:=(\Gamma^{m} \Gamma^{n})^{a}_{c};
82-
expand(_);
82+
expand(_)
8383
tst:= (\Gamma^{m})^{a}_{b} (\Gamma^{n})^{b}_{c} - @(ex);
8484
assert(tst==0)
8585
print("Test 07 passed")
8686

8787
test07()
8888

89+
def test08():
90+
__cdbkernel__=create_scope()
91+
{a,b}::Indices(vector);
92+
ex:=(u)_{a} (M)^{a b} (v)_{b};
93+
combine(_)
94+
tst:= \indexbracket(u M v) - @(ex);
95+
assert(tst==0)
96+
print("Test 08 passed")
97+
98+
test08()
99+

0 commit comments

Comments
 (0)