1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
|
#include "algorithms/expand_delta.hh"
#include "properties/KroneckerDelta.hh"
#include "properties/TableauBase.hh"
#include "Cleanup.hh"
using namespace cadabra;
#define ASYM 1
expand_delta::expand_delta(const Kernel& k, Ex& tr)
: Algorithm(k, tr)
{
}
bool expand_delta::can_apply(iterator st)
{
const KroneckerDelta *kr=kernel.properties.get<KroneckerDelta>(st);
if(kr)
if(tr.number_of_children(st)>2)
return true;
return false;
}
Algorithm::result_t expand_delta::apply(iterator& st)
{
// std::cerr << "applying " << Ex(st) << std::endl;
Ex rep("\\sum");
iterator sum=rep.begin();
combin::combinations<Ex> ci;
#ifdef ASYM
// Determine implicit anti-symmetrisation by looking for anti-symmetric
// tensors in the product (if any). Store these sets in `implicit_sets`.
std::vector<std::vector<Ex>> implicit_sets;
if(tr.is_head(st)==false && *(tr.parent(st)->name)=="\\prod") {
sibling_iterator oth=tr.begin(tr.parent(st));
while(oth!=tr.end(tr.parent(st))) {
if(oth!=st) {
const TableauBase *tb=kernel.properties.get<TableauBase>(oth);
if(tb) {
// std::cerr << "object " << oth << " has symmetry" << std::endl;
for(unsigned int i=0; i<tb->size(kernel.properties, tr, oth); ++i) {
// std::cerr << "object " << oth << " tab " << i << std::endl;
TableauBase::tab_t tmptab(tb->get_tab(kernel.properties, tr, oth, i));
for(unsigned int asset=0; asset<tmptab.row_size(0); ++asset) {
if(tmptab.column_size(asset)>1) {
// std::cerr << " asym size " << tmptab.column_size(asset) << std::endl;
std::vector<Ex> iset;
for(unsigned int asel=0; asel<tmptab.column_size(asset); ++asel) {
int indexnum=tmptab(asel,asset);
auto oth2=begin_index(oth);
oth2+=indexnum;
iset.push_back(*oth2);
}
if(iset.size()>1)
implicit_sets.push_back(iset);
}
}
}
}
}
++oth;
}
}
#endif
// std::cerr << implicit_sets.size() << std::endl;
// for(const auto& implicit_set: implicit_sets) {
// for(const auto& implicit_index: implicit_set) {
// std::cerr << *implicit_index.begin()->name << " ";
// }
// std::cerr << std::endl;
// }
// Determine the overlap factors. For each index set
// which is to be kept implicitly anti-symmetrised, we need
// to determine how many of those indices actually sit on
// the Kronecker delta. E.g. \delta^a_m^b_n A_{m n p} will
// have an implicit set [m,n,p], but there are only two
// indices overlapping, so the reduction in number of terms
// is 2, not 6.
multiplier_t divfactor1=1, divfactor2=1;
#ifdef ASYM
Ex_comparator comp(kernel.properties);
for(const auto& implicit_set: implicit_sets) {
unsigned int overlap1=0, overlap2=0;
for(const auto& implicit_index: implicit_set) {
sibling_iterator ind=tr.begin(st);
for(unsigned int pairnum=0; pairnum<tr.number_of_children(st); pairnum+=2) {
// std::cerr << "does " << ind << " match " << implicit_index.begin() << std::endl;
comp.clear();
auto m1=comp.equal_subtree(ind, implicit_index.begin(), Ex_comparator::useprops_t::never, true);
// std::cerr << static_cast<int>(m1) << std::endl;
if(m1<=Ex_comparator::match_t::subtree_match) ++overlap1;
++ind;
comp.clear();
auto m2=comp.equal_subtree(ind, implicit_index.begin(), Ex_comparator::useprops_t::never, true);
if(m2<=Ex_comparator::match_t::subtree_match) ++overlap2;
++ind;
}
}
divfactor1*=combin::fact(overlap1);
divfactor2*=combin::fact(overlap2);
}
#endif
// std::cerr << divfactor1 << ", " << divfactor2 << std::endl;
bool permute_second_set=(divfactor2>divfactor1);
// Construct the original permutation.
sibling_iterator ind=tr.begin(st);
for(unsigned int pairnum=0; pairnum<tr.number_of_children(st); pairnum+=2) {
if(permute_second_set)
++ind;
if(std::find(ci.original.begin(), ci.original.end(), *ind)!=ci.original.end()) {
zero(st->multiplier);
cleanup_dispatch(kernel, tr, st);
return result_t::l_applied;
}
ci.original.push_back(*ind);
++ind;
if(!permute_second_set)
++ind;
ci.sublengths.push_back(1);
}
#ifdef ASYM
// Now construct the output asym ranges.
for(const auto& implicit_set: implicit_sets) {
combin::range_t asymrange1;
for(const auto& implicit_index: implicit_set) {
for(unsigned int i1=0; i1<ci.original.size(); ++i1) {
comp.clear();
auto m = comp.equal_subtree(ci.original[i1].begin(), implicit_index.begin(), Ex_comparator::useprops_t::never, true);
if(m<=Ex_comparator::match_t::subtree_match) {
// std::cerr << "add to asymrange: " << i1 << std::endl;
asymrange1.push_back(i1);
break;
}
}
}
ci.input_asym.push_back(asymrange1);
}
#endif
// Now generate all the permutations of the indices.
ci.permute();
// std::cerr << "number of permutations: " << ci.size() << std::endl;
// For each permutation, add a product of deltas with two indices.
for(unsigned int permnum=0; permnum<ci.size(); ++permnum) {
// FIXME: replace with progress indicator
// if(permnum%10000==0)
// std::cerr << "permutation " << permnum << std::endl;
iterator prod=tr.append_child(sum, str_node("\\prod"));
int sgn=combin::ordersign(ci[permnum].begin(), ci[permnum].end(),
ci.original.begin(), ci.original.end());
multiplier_t mult=multiplier_t(ci.multiplier(permnum))/multiplier_t(combin::fact((unsigned long)ci.original.size()))*sgn;
// std::cerr << "multipliers: " << mult << " = " << ci.multiplier(permnum)
// << " / " << combin::fact(ci.original.size()) << " * " << sgn << std::endl;
multiply(prod->multiplier, mult);
multiply(prod->multiplier, *st->multiplier);
sibling_iterator ind=tr.begin(st);
for(unsigned int pairnum=0; pairnum<ci.original.size(); ++pairnum) {
iterator delta=tr.append_child(prod, str_node(st->name));//"\\delta"));
if(permute_second_set)
tr.append_child(delta, (iterator)(ind));
tr.append_child(delta, ci[permnum][pairnum].begin());
++ind;
if(!permute_second_set)
tr.append_child(delta, (iterator)(ind));
++ind;
}
}
// std::cerr << "flatten" << std::endl;
if(rep.number_of_children(sum)==1) { // flatten \sum node if there was only one term
rep.flatten(sum);
rep.erase(sum);
}
st=tr.replace(st, rep.begin());
// std::cerr << "cleanup" << std::endl;
cleanup_dispatch(kernel, tr, st);
// std::cerr << "done" << std::endl;
// After the return, we need to indicate that no new dummies were
// introduced, otherwise we go into rename_replacement_dummies,
// which can be extremely time-consuming.
// std::cerr << "####### returning new value " << std::endl;
return result_t::l_applied_no_new_dummies;
}
|
