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
|
#include "Cleanup.hh"
#include "algorithms/eliminate_metric.hh"
#include "properties/Metric.hh"
#include "properties/InverseMetric.hh"
//#define DEBUG 1
using namespace cadabra;
eliminate_metric::eliminate_metric(const Kernel& k, Ex& e, Ex& pref, bool redundant)
: eliminate_converter(k, e, pref, redundant)
{
}
bool eliminate_metric::is_conversion_object(iterator fit) const
{
const Metric *vb=kernel.properties.get<Metric>(fit);
const InverseMetric *ivb=kernel.properties.get<InverseMetric>(fit);
if(vb || ivb) return true;
else return false;
}
eliminate_converter::eliminate_converter(const Kernel& k, Ex& e, Ex& pref, bool redun)
: Algorithm(k, e), preferred(pref), redundant(redun)
{
}
bool eliminate_converter::can_apply(iterator it)
{
if(*it->name=="\\prod") return true;
return false;
}
bool eliminate_converter::handle_one_index(index_iterator ind1, index_iterator ind2, iterator fit, sibling_iterator objs)
{
bool replaced=false;
#ifdef DEBUG
std::cerr << "handling " << ind1 << ", " << ind2 << ", " << fit << std::endl;
#endif
// For a conversion to be possible, we need one upper and one lower index
// (for position!=free) or two indices (for position=free).
#ifdef DEBUG
std::cerr << "dummies: " << std::endl;
for(const auto& d: ind_dummy)
std::cerr << d.second << std::endl;
#endif
auto locs= ind_dummy.equal_range(Ex(ind1));
int num1=std::distance(locs.first, locs.second);
Ex other(ind1);
other.begin()->flip_parent_rel();
locs = ind_dummy.equal_range(other);
int num2=std::distance(locs.first, locs.second);
#ifdef DEBUG
std::cerr << "num1 = " << num1 << ", num2 = " << num2 << std::endl;
#endif
if(num1==1 && num2==1) {
while(locs.first!=locs.second) {
if(locs.first->second!=(iterator)ind1) {
// Does this index sit on an object in the "preferred form" list?
// (if there is no preferred form, always eliminate)
if(separated_by_derivative(locs.first->second, ind2, fit)==false) {
if(objs==preferred.end()) { // no
#ifdef DEBUG
std::cerr << "Eliminate regardless" << std::endl;
#endif
tr.move_ontop(locs.first->second, iterator(ind2))->fl.parent_rel=ind2->fl.parent_rel;;
fit=tr.erase(fit);
replaced=true;
}
else { // yes
#ifdef DEBUG
std::cerr << "Index on preferred form tensor" << std::endl;
#endif
iterator par=tr.parent(locs.first->second);
sibling_iterator prefit=tr.begin(objs);
while(prefit!=tr.end(objs)) {
if(subtree_equal(&kernel.properties, prefit, par, -1, false)) {
tr.move_ontop(locs.first->second, iterator(ind2))->fl.parent_rel=ind2->fl.parent_rel;;
fit=tr.erase(fit);
replaced=true;
break;
}
++prefit;
}
}
}
}
if(replaced)
return true;
++locs.first;
}
}
return false;
}
Algorithm::result_t eliminate_converter::apply(iterator& it)
{
result_t res=result_t::l_no_action;
// Put arguments in canonical form.
sibling_iterator objs=preferred.begin();
if(objs!=preferred.end())
if(*objs->name!="\\comma")
objs=tr.wrap(objs, str_node("\\comma"));
ind_free.clear();
ind_dummy.clear();
classify_indices(it, ind_free, ind_dummy);
// Run over all factors, find metrics, figure out whether they can
// be used to turn the indices on the other tensors to preferred type.
sibling_iterator fit=tr.begin(it);
while(fit!=tr.end(it)) {
if(is_conversion_object(fit)) {
index_iterator ind1=index_iterator::begin(kernel.properties, fit), ind2=ind1;
++ind2;
if (!redundant || (ind_free.find(Ex(ind1)) == ind_free.end() &&
ind_free.find(Ex(ind2)) == ind_free.end())) {
// 1st index to 2nd index conversion?
if(handle_one_index(ind1, ind2, fit, objs)) {
res=result_t::l_applied;
break;
}
// 2nd index to 1st index conversion?
if(handle_one_index(ind2, ind1, fit, objs)) {
res=result_t::l_applied;
break;
}
}
}
++fit;
}
if(res==result_t::l_applied)
cleanup_dispatch(kernel, tr, it);
return res;
}
|
