11
22#include " explicit_indices.hh"
33
4+ #include " algorithms/substitute.hh"
45#include " properties/ImplicitIndex.hh"
56
67using namespace cadabra ;
@@ -17,13 +18,117 @@ bool explicit_indices::can_apply(iterator st)
1718 // All this because when we generate free indices, we need to
1819 // ensure that all terms and all sides of an equals node use the
1920 // same index names.
21+
22+ if (*st->name ==" \\ equals" return false ; // switch
23+ if (*st->name ==" \\ sum" return true ;
24+ if (is_termlike (st)) {
25+ if (tr.is_head (st)) return true ;
26+ if (*tr.parent (st)->name ==" \\ sum" return false ;
27+ return true ;
28+ }
2029
2130 return false ;
2231 }
2332
24- Algorithm::result_t explicit_indices::apply (iterator& prod )
33+ Algorithm::result_t explicit_indices::apply (iterator& it )
2534 {
35+ result_t res=result_t ::l_no_action;
36+
37+ // Ensure that we are always working on a sum, even
38+ // if there is only one term.
39+ if (is_termlike (it))
40+ force_node_wrap (it, " \\ sum"
41+
42+ // Classify all free and dummy indices already present. Any new
43+ // indices cannot be taken from these.
44+ index_map_t ind_free_sum, ind_dummy_sum;
45+ classify_indices (it, ind_free_sum, ind_dummy_sum);
46+ for (auto & k: ind_free_sum)
47+ std::cerr << k.first << std::endl;
48+ for (auto & k: ind_dummy_sum)
49+ std::cerr << k.first << std::endl;
50+
51+ sibling_iterator term=tr.begin (it);
52+ while (term!=tr.end (it)) {
53+ iterator tmp=term;
54+ prod_wrap_single_term (tmp);
55+ term=tmp;
56+
57+ // For each index set, keep track of the last used index in
58+ // building the explicit index line.
59+ std::map<const Indices *, Ex::iterator> index_lines;
60+
61+ sibling_iterator factor=tr.begin (term);
62+ while (factor!=tr.end (term)) {
63+ int tmp;
64+ auto ii = kernel.properties .get_with_pattern <ImplicitIndex>(factor, tmp);
65+ if (ii.first ) {
66+ // Determine indices on this factor. Use a copy because we
67+ // need this object later.
68+ Ex orig (factor);
69+ index_map_t factor_ind_free, factor_ind_dummy;
70+ classify_indices (orig.begin (), factor_ind_free, factor_ind_dummy);
71+
72+ // Substitute explcit replacement and rename the indices
73+ // already present on the implicit factor.
74+ Ex replacement (" \\ arrow"
75+ replacement.append_child (replacement.begin (), ii.second ->obj .begin ());
76+ replacement.append_child (replacement.begin (), ii.first ->explicit_form .begin ());
77+ substitute subs (kernel, tr, replacement);
78+ iterator factor_tmp=factor;
79+ if (subs.can_apply (factor_tmp))
80+ subs.apply (factor_tmp);
81+ else
82+ throw InternalError (" Internal inconsistency encountered, aborting."
83+ factor=factor_tmp;
84+
85+ // Determine indices on the replacement.
86+ index_map_t repl_ind_free, repl_ind_dummy;
87+ classify_indices (factor, repl_ind_free, repl_ind_dummy);
88+
89+ // Which indices have been added?
90+ index_map_t ind_same;
91+ IndexClassifier ic (kernel);
92+ ic.determine_intersection (factor_ind_free, repl_ind_free, ind_same, true );
93+
94+ // Go through indices in order of appearance, determine if they have
95+ // been added, and rename. Note: indices do not appear in order in the
96+ // index maps above.
97+ auto iit=index_iterator::begin (kernel.properties , factor);
98+ while (iit!=index_iterator::end (kernel.properties , factor)) {
99+ auto search=repl_ind_free.begin ();
100+ bool found=false ;
101+ while (search!=repl_ind_free.end ()) {
102+ if (search->second == (iterator)iit) {
103+ found=true ;
104+ break ;
105+ }
106+ ++search;
107+ }
108+ if (found) {
109+ // This index was added.
110+ // Get a new free index.
111+
112+ auto ip = kernel.properties .get <Indices>(search->second );
113+ if (!ip)
114+ throw InternalError (" Do not have Indices property for all implicit indices."
115+
116+ std::cerr << " getting dummy index"
117+ auto di = ic.get_dummy (ip, &ind_free_sum, &ind_dummy_sum);
118+ std::cerr << di << std::endl;
119+ tr.replace (search->second , di.begin ());
120+ }
121+
122+ ++iit;
123+ }
124+ }
125+ ++factor;
126+ }
26127
128+ tmp=term;
129+ prod_unwrap_single_term (tmp);
130+ ++term;
131+ }
27132
28- return result_t ::l_applied ;
133+ return res ;
29134 }
0 commit comments