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
|
#include "algorithms/epsilon_to_delta.hh"
#include "algorithms/reduce_delta.hh"
#include "properties/EpsilonTensor.hh"
#include "properties/Metric.hh"
using namespace cadabra;
epsilon_to_delta::epsilon_to_delta(const Kernel& k, Ex& tr, bool r)
: Algorithm(k, tr), reduce(r)
{
}
bool epsilon_to_delta::can_apply(iterator st)
{
if(*st->name!="\\prod")
return false;
epsilons.clear();
// Find the two epsilon tensors in the product.
std::multimap<std::string,Ex::iterator> emap;
sibling_iterator it=tr.begin(st);
while(it!=tr.end(st)) {
const EpsilonTensor *eps=kernel.properties.get<EpsilonTensor>(it);
if(eps) emap.insert(std::pair<std::string,Ex::iterator>(get_index_set_name(begin_index(it)),it));
++it;
}
signature=1;
std::multimap<std::string,Ex::iterator>::iterator eit=emap.begin();
while(epsilons.size()<2 && eit!=emap.end()) {
if(emap.count(eit->first)>1) {
epsilons.push_back(eit->second);
++eit;
epsilons.push_back(eit->second);
const EpsilonTensor *eps=kernel.properties.get<EpsilonTensor>(eit->second);
if(eps->metric.begin()!=eps->metric.end()) {
const Metric *met=kernel.properties.get<Metric>(eps->metric.begin());
if(met)
signature=met->signature;
}
if(eps->krdelta.begin()!=eps->krdelta.end())
repdelta=eps->krdelta;
}
++eit;
}
if(epsilons.size()<2)
return false;
if(repdelta.begin()==repdelta.end())
return false;
else repdelta.erase_children(repdelta.begin());
return true;
}
Algorithm::result_t epsilon_to_delta::apply(iterator& st)
{
Ex rep(repdelta);
iterator delta=rep.begin();
// std::cerr << Ex(st) << std::endl;
sibling_iterator eps1=tr.begin(epsilons[0]);
sibling_iterator eps2=tr.begin(epsilons[1]);
while(eps1!=tr.end(epsilons[0])) {
rep.append_child(delta, *eps1);
rep.append_child(delta, *eps2);
++eps1;
++eps2;
}
// std::cerr << "result: " << rep.begin() << std::endl;
multiply(st->multiplier, *epsilons[0]->multiplier);
multiply(st->multiplier, *epsilons[1]->multiplier);
tr.erase(epsilons[0]);
// std::cerr << tr.number_of_children(epsilons[1]) << std::endl;
// std::cerr << " -> " << *st->multiplier << " * " << combin::fact(multiplier_t(tr.number_of_children(epsilons[1]))) << std::endl;
multiply(st->multiplier, combin::fact(multiplier_t(tr.number_of_children(epsilons[1]))));
// std::cerr << "A:" << *st->multiplier << std::endl;
multiply(st->multiplier, signature);
// std::cerr << "B:" << *st->multiplier << std::endl;
iterator gend=tr.replace(epsilons[1], rep.begin());
// std::cerr << "B2:" << *st->multiplier << std::endl;
if(reduce) {
// std::cerr << "reducing" << std::endl;
reduce_delta rg(kernel, tr);
if(rg.can_apply(gend))
rg.apply(gend);
if(*gend->multiplier==0) {
zero(st->multiplier);
return result_t::l_applied;
}
}
if(*gend->multiplier!=1) {
// std::cerr << "B3:" << *st->multiplier << std::endl;
// std::cerr << "B3:" << *gend->multiplier << std::endl;
multiply(tr.parent(gend)->multiplier, *gend->multiplier);
one(gend->multiplier);
// std::cerr << "C:" << *st->multiplier << std::endl;
}
if(tr.number_of_children(st)==1) {
multiply(tr.begin(st)->multiplier, *st->multiplier);
tr.flatten(st);
st=tr.erase(st);
}
// std::cerr << Ex(st) << std::endl;
return result_t::l_applied;
}
|
