33
44#include " algorithms/substitute.hh"
55#include " properties/ImplicitIndex.hh"
6+ #include " properties/PartialDerivative.hh"
67
78using namespace cadabra ;
89
@@ -41,7 +42,8 @@ Algorithm::result_t explicit_indices::apply(iterator& it)
4142
4243 // Classify all free and dummy indices already present. Any new
4344 // indices cannot be taken from these.
44- index_map_t ind_free_sum, ind_dummy_sum;
45+ ind_free_sum.clear ();
46+ ind_dummy_sum.clear ();
4547 classify_indices (it, ind_free_sum, ind_dummy_sum);
4648
4749 sibling_iterator term=tr.begin (it);
@@ -52,81 +54,25 @@ Algorithm::result_t explicit_indices::apply(iterator& it)
5254
5355 // For each index set, keep track of the last used index in
5456 // building the explicit index line.
55- index_map_t added_this_term;
56- std::map< const Indices *, Ex::iterator> index_lines;
57+ added_this_term. clear () ;
58+ index_lines. clear () ;
5759
5860 sibling_iterator factor=tr.begin (term);
5961 while (factor!=tr.end (term)) {
60- int tmp;
61- auto ii = kernel.properties .get_with_pattern <ImplicitIndex>(factor, tmp);
62- if (ii.first ) {
63- // Determine indices on this factor. Use a copy because we
64- // need this object later.
65- Ex orig (factor);
66- index_map_t factor_ind_free, factor_ind_dummy;
67- classify_indices (orig.begin (), factor_ind_free, factor_ind_dummy);
68-
69- // Substitute explcit replacement and rename the indices
70- // already present on the implicit factor.
71- Ex replacement (" \\ arrow" );
72- replacement.append_child (replacement.begin (), ii.second ->obj .begin ());
73- replacement.append_child (replacement.begin (), ii.first ->explicit_form .begin ());
74- substitute subs (kernel, tr, replacement);
75- iterator factor_tmp=factor;
76- if (subs.can_apply (factor_tmp))
77- subs.apply (factor_tmp);
78- else
79- throw InternalError (" Internal inconsistency encountered, aborting." );
80- factor=factor_tmp;
81-
82- // Determine indices on the replacement.
83- index_map_t repl_ind_free, repl_ind_dummy;
84- classify_indices (factor, repl_ind_free, repl_ind_dummy);
85-
86- // Which indices have been added?
87- index_map_t ind_same;
88- IndexClassifier ic (kernel);
89- ic.determine_intersection (factor_ind_free, repl_ind_free, ind_same, true );
90-
91- // Go through indices in order of appearance, determine if they have
92- // been added, and rename. Note: indices do not appear in order in the
93- // index maps above.
94- auto iit=index_iterator::begin (kernel.properties , factor);
95- while (iit!=index_iterator::end (kernel.properties , factor)) {
96- auto search=repl_ind_free.begin ();
97- bool found=false ;
98- while (search!=repl_ind_free.end ()) {
99- if (search->second == (iterator)iit) {
100- found=true ;
101- break ;
102- }
103- ++search;
104- }
105- ++iit; // Update now, we may be replacing this index.
106- if (found) {
107- // This index was added.
108- const Indices *ip = kernel.properties .get <Indices>(search->second );
109- if (!ip)
110- throw InternalError (" Do not have Indices property for all implicit indices." );
111-
112- // Determine if we have an 'active' index line for
113- // this index type.
114- auto line = index_lines.find (ip);
115- if (line==index_lines.end ()) {
116- // No active line. Get a new free index.
117- auto di = ic.get_dummy (ip, &ind_free_sum, &ind_dummy_sum, &added_this_term);
118- auto loc = tr.replace_index (search->second , di.begin (), true );
119- added_this_term.insert (index_map_t::value_type (di, loc));
120- index_lines[ip]=loc;
121- }
122- else {
123- // Use the active line index, then unset the active line.
124- auto loc = tr.replace_index (search->second , line->second , true );
125- index_lines.erase (line);
126- }
62+ const PartialDerivative *pd = kernel.properties .get <PartialDerivative>(factor);
63+ if (pd) {
64+ sibling_iterator args=tr.begin (factor);
65+ while (args!=tr.end (factor)) {
66+ if (args->fl .parent_rel ==str_node::p_none) {
67+ handle_factor (args);
68+ break ;
12769 }
70+ ++args;
12871 }
12972 }
73+ else {
74+ handle_factor (factor);
75+ }
13076 ++factor;
13177 }
13278
@@ -137,3 +83,83 @@ Algorithm::result_t explicit_indices::apply(iterator& it)
13783
13884 return res;
13985 }
86+
87+ void explicit_indices::handle_factor (sibling_iterator& factor)
88+ {
89+ int tmp;
90+ auto ii = kernel.properties .get_with_pattern <ImplicitIndex>(factor, tmp);
91+ if (ii.first ) {
92+ // Determine indices on this factor. Use a copy because we
93+ // need this object later.
94+ Ex orig (factor);
95+ index_map_t factor_ind_free, factor_ind_dummy;
96+ classify_indices (orig.begin (), factor_ind_free, factor_ind_dummy);
97+
98+ // Substitute explcit replacement and rename the indices
99+ // already present on the implicit factor.
100+ Ex replacement (" \\ arrow" );
101+ replacement.append_child (replacement.begin (), ii.second ->obj .begin ());
102+ replacement.append_child (replacement.begin (), ii.first ->explicit_form .begin ());
103+ substitute subs (kernel, tr, replacement);
104+ iterator factor_tmp=factor;
105+ if (subs.can_apply (factor_tmp))
106+ subs.apply (factor_tmp);
107+ else
108+ throw InternalError (" Internal inconsistency encountered, aborting." );
109+ factor=factor_tmp;
110+
111+ // Determine indices on the replacement.
112+ index_map_t repl_ind_free, repl_ind_dummy;
113+ classify_indices (factor, repl_ind_free, repl_ind_dummy);
114+
115+ // Which indices have been added?
116+ index_map_t ind_same;
117+ IndexClassifier ic (kernel);
118+ ic.determine_intersection (factor_ind_free, repl_ind_free, ind_same, true );
119+
120+ // Keep track of indices already added in this factor, to avoid making
121+ // an index trace on a single factor.
122+ std::map<const Indices *, Ex::iterator> index_lines_factor;
123+
124+ // Go through indices in order of appearance, determine if they have
125+ // been added, and rename. Note: indices do not appear in order in the
126+ // index maps above.
127+ auto iit=index_iterator::begin (kernel.properties , factor);
128+ while (iit!=index_iterator::end (kernel.properties , factor)) {
129+ auto search=repl_ind_free.begin ();
130+ bool found=false ;
131+ while (search!=repl_ind_free.end ()) {
132+ if (search->second == (iterator)iit) {
133+ found=true ;
134+ break ;
135+ }
136+ ++search;
137+ }
138+ ++iit; // Update now, we may be replacing this index.
139+ if (found) {
140+ // This index was added.
141+ const Indices *ip = kernel.properties .get <Indices>(search->second );
142+ if (!ip)
143+ throw InternalError (" Do not have Indices property for all implicit indices." );
144+
145+ // Determine if we have an 'active' index line for
146+ // this index type.
147+ auto line = index_lines.find (ip);
148+ auto line_factor = index_lines_factor.find (ip);
149+ if (line==index_lines.end () || line_factor!=index_lines_factor.end ()) {
150+ // No active line. Get a new free index.
151+ auto di = ic.get_dummy (ip, &ind_free_sum, &ind_dummy_sum, &added_this_term);
152+ auto loc = tr.replace_index (search->second , di.begin (), true );
153+ added_this_term.insert (index_map_t::value_type (di, loc));
154+ index_lines[ip]=loc;
155+ index_lines_factor[ip]=loc;
156+ }
157+ else {
158+ // Use the active line index, then unset the active line.
159+ auto loc = tr.replace_index (search->second , line->second , true );
160+ index_lines.erase (line);
161+ }
162+ }
163+ }
164+ }
165+ }
0 commit comments