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;
}
|