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
|
#include "IndexIterator.hh"
#include "Exceptions.hh"
#include "Cleanup.hh"
#include "properties/Indices.hh"
#include "properties/Integer.hh"
#include "algorithms/young_project_tensor.hh"
#include "algorithms/collect_terms.hh"
#include "algorithms/indexsort.hh"
using namespace cadabra;
young_project_tensor::young_project_tensor(const Kernel& k, Ex& tr, bool modmono)
: Algorithm(k, tr), modulo_monoterm(modmono)
{
}
bool young_project_tensor::can_apply(iterator it)
{
if(*it->name=="\\prod") return false;
tb=kernel.properties.get<TableauBase>(it);
if(tb) {
if(tb->size(kernel.properties, tr, it)>0)
return true;
}
return false;
}
Algorithm::result_t young_project_tensor::apply(iterator& it)
{
// std::cout << "at " << *it->name << std::endl;
assert(tb);
// txtout << typeid(*tb).name() << std::endl;
// FIXME: handle cases with multiple tabs.
TableauBase::tab_t tab=tb->get_tab(kernel.properties, tr, it, 0);
if(tab.number_of_rows()==0)
return result_t::l_no_action;
// txtout << tab << std::endl;
if(modulo_monoterm) {
if(tab.number_of_rows()==1) // Nothing happends with fully symmetric tensors modulo monoterm.
return result_t::l_no_action;
if(tab.row_size(0)==1 && tab.selfdual_column==0) // Ditto for fully anti-symmetric tensors & modmono.
return result_t::l_no_action;
}
// For non-trivial tableau shapes, apply the Young projector.
Ex rep;
rep.set_head(str_node("\\sum"));
if(tab.row_size(0)>0) {
sym.clear();
tab.projector(sym); //, modulo_monoterm);
// txtout << sym.size() << std::endl;
for(unsigned int i=0; i<sym.size(); ++i) {
Ex repfac(it);
for(unsigned int j=0; j<sym[i].size(); ++j) {
index_iterator src_fd=index_iterator::begin(kernel.properties, it);
index_iterator dst_fd=index_iterator::begin(kernel.properties, repfac.begin());
// txtout << sym[i][j] << " " << sym.original[j] << std::endl;
src_fd+=sym[i][j];
dst_fd+=sym.original[j];
// txtout << *src_fd->name << std::endl;
// txtout << *dst_fd->name << std::endl;
dst_fd->name=src_fd->name;
}
multiply(repfac.begin()->multiplier, sym.signature(i));
multiply(repfac.begin()->multiplier, tb->get_tab(kernel.properties, tr, it, 0).projector_normalisation());
iterator newtensor=rep.append_child(rep.begin(), repfac.begin());
if(modulo_monoterm) { // still necessary for column exchange
indexsort isort(kernel, rep);
[[maybe_unused]] auto res=isort.can_apply(newtensor); // to set tb
assert(res);
isort.apply(newtensor);
}
}
collect_terms cterms(kernel, rep);
iterator rephead=rep.begin();
cterms.apply(rephead);
}
else {
rep.append_child(rep.begin(), it);
}
// If there is a selfdual or anti-selfdual component, we need to add a term
// for each generated term, which contains an epsilon. In the generated term,
// find the positions of the indices which originally sat on the selfdual tensor,
// then replace these with dummies which are repeated on the epsilon tensor.
// We should do all of this _here_, not before the Young projector, because we
// want to avoid using indexsort.
if(tab.selfdual_column!=0) {
// Classify indices so we can insert dummies.
index_map_t one, two, three, four, added_dummies;
classify_indices_up(it, one, two);
classify_indices(it, three, four);
// Figure out the properties of the indices for which we want dummy partners.
index_iterator iit=index_iterator::begin(kernel.properties, it);
iit+=tab(0,abs(tab.selfdual_column)-1);
const Integer *itg=kernel.properties.get<Integer>(iit, true);
const Indices *ind=kernel.properties.get<Indices>(iit, true);
if(itg==0)
throw ConsistencyException("young_project_tensor: Need to know the range of the indices.");
if(ind==0)
throw ConsistencyException("young_project_tensor: Need to have a set of dummy indices.");
// Step through all generated terms and give each one an epsilon partner term.
sibling_iterator tt=rep.begin(rep.begin());
while(tt!=rep.end(rep.begin())) {
Ex repfac(tt);
iterator prodit=repfac.wrap(repfac.begin(), str_node("\\prod"));
iterator tensit=repfac.begin(prodit);
// FIXME: take care of Euclidean signature cases.
// repfac.insert(repfac.begin(prodit), str_node("I"));
iterator epsit =repfac.append_child(prodit, str_node("\\epsilon"));
// Normalise the epsilon term appropriately.
multiply(prodit->multiplier,
multiplier_t(1)/combin::factorial(to_long(*(itg->difference.begin()->multiplier)/2)));
if(tab.selfdual_column<0)
flip_sign(prodit->multiplier);
// Move the free indices to the epsilon factor.
for(unsigned int row=0; row<tab.column_size(abs(tab.selfdual_column)-1); ++row) {
// Find index of original in newly generated term (all indices are different,
// so this works).
iit=index_iterator::begin(kernel.properties, tensit);
index_iterator iit_orig=index_iterator::begin(kernel.properties, it);
iit_orig+=tab(row, abs(tab.selfdual_column)-1);
while(subtree_exact_equal(&kernel.properties, iit, iit_orig)==false)
++iit;
Ex dum=get_dummy(ind, &one, &two, &three, &four, &added_dummies);
repfac.append_child(epsit, iterator(iit)); // move index to eps
iterator repind=rep.replace_index(iterator(iit), dum.begin()); // replace index on tens
added_dummies.insert(index_map_t::value_type(dum, repind));
}
// Now insert the new dummies in the epsilon.
index_map_t::iterator adi=added_dummies.begin();
while(adi!=added_dummies.end()) {
iterator ni=repfac.append_child(epsit, adi->first.begin());
ni->fl.parent_rel=str_node::p_sub;
++adi;
}
added_dummies.clear();
++tt;
rep.insert_subtree(tt, repfac.begin());
}
}
// Final cleanup
it=tr.replace(it,rep.begin());
cleanup_dispatch(kernel, tr, it);
return result_t::l_applied;
}
|
