File: unwrap.cc

package info (click to toggle)
cadabra2 2.4.3.2-2
  • links: PTS, VCS
  • area: main
  • in suites: sid
  • size: 78,732 kB
  • sloc: ansic: 133,450; cpp: 92,064; python: 1,530; javascript: 203; sh: 184; xml: 182; objc: 53; makefile: 51
file content (340 lines) | stat: -rw-r--r-- 9,546 bytes parent folder | download | duplicates (3)
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
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340

#include "Compare.hh"
#include "Cleanup.hh"
#include "algorithms/unwrap.hh"
#include "properties/Derivative.hh"
#include "properties/Accent.hh"
#include "properties/DiracBar.hh"
#include "properties/Spinor.hh"
#include "properties/DifferentialForm.hh"
#include "properties/GammaMatrix.hh"
#include "properties/DependsBase.hh"
//#include "algorithms/prodcollectnum.hh"

// #define DEBUG 1

using namespace cadabra;

unwrap::unwrap(const Kernel& k, Ex& tr, Ex& w)
	: Algorithm(k, tr)
	{
	if(w.begin()!=w.end()) {
		if(*w.begin()->name!="\\comma")
			wrappers.push_back(w);
		else {
			auto sib=w.begin(w.begin());
			while(sib!=w.end(w.begin())) {
				wrappers.push_back(Ex(sib));
				++sib;
				}
			}
		}
	}

bool unwrap::can_apply(iterator it)
	{
	const Derivative *der=kernel.properties.get<Derivative>(it);
	const Accent     *acc=kernel.properties.get<Accent>(it);
	if(der || acc) {
		Ex_comparator comp(kernel.properties);
		if(wrappers.size()>0) {
			bool found=false;
			for(auto&w : wrappers) {
#ifdef DEBUG
				std::cerr << "comparing " << w << " to " << Ex(it) << std::endl;
#endif
				if(comp.equal_subtree(w.begin(), it)==Ex_comparator::match_t::subtree_match) {
#ifdef DEBUG
					std::cerr << "yes" << std::endl;
#endif
					found=true;
					break;
					}
				}
			if(!found) return false;
			}
		return true;
		}

	if(*it->name=="\\wedge") return true;

	return false;
	}

Algorithm::result_t unwrap::apply_on_wedge(iterator& it)
	{
	result_t res=result_t::l_no_action;

	auto prodwrap = tr.wrap(it, str_node("\\prod"));
	
	sibling_iterator sib=tr.begin(it);
	while(sib!=tr.end(it)) {
		if(*sib->name=="\\prod") {
			sibling_iterator fac=tr.begin(sib);
			while(fac!=tr.end(sib)) {
				const DifferentialFormBase *diff = kernel.properties.get<DifferentialFormBase>(fac);
				sibling_iterator nxt=fac;
				++nxt;
				if(diff==0 || diff->degree(kernel.properties, fac).begin()->is_zero() ) {
					// Figure out the sign. First move to front of existing
					// simple product. Then move through all other factors
					// in the wedge product.
					Ex_comparator comp(kernel.properties);
					int sign=comp.can_move_to_front(tr, sib, fac);
					if(sign!=0) {
						sibling_iterator prevwedfac=sib;
						if(prevwedfac!=tr.begin(it)) {
							do {
								--prevwedfac;
								auto stc=comp.equal_subtree(prevwedfac, fac);
								sign*=comp.can_swap(prevwedfac, fac, stc);
								} while(prevwedfac!=tr.begin(it));
							}
						if(sign!=0) {
							tr.move_before(it, fac);
							multiply(prodwrap->multiplier, sign);
							}
						}
					}
				fac=nxt;
				}
			}
		sibling_iterator nxt=sib;
		++nxt;
		iterator tmp=sib;
		cleanup_dispatch(kernel, tr, tmp);
		sib=nxt;
		}
	cleanup_dispatch(kernel, tr, prodwrap);	
	
	return res;
	}

// \del{a*f*A*b*C*d*e} -> a*f*\del{A*b*C}*d*e
//
// locate first dependent factor
// locate first following independent factor
// ...
//
// Should also work for brackets, like
//
//    \poisson( A )( B ) -> 0
//
// if either A or B has vanishing poisson bracket, and
//
//    \poisson( N1 D1 )( N2 D2 ) -> N1 N2 \poisson( D1 )( D2 ).
//
// So all "derivatives"

Algorithm::result_t unwrap::apply(iterator& it)
	{
#ifdef DEBUG
	std::cerr << "Applying unwrap at " << Ex(it) << std::endl;
#endif

	if(*it->name=="\\wedge")
		return apply_on_wedge(it);
	
	result_t res = result_t::l_no_action;

	bool is_accent=kernel.properties.get<Accent>(it);
	bool is_diracbar=kernel.properties.get<DiracBar>(it);

	// Wrap the 'derivative' in a product node so we can take
	// child nodes out and stuff them inside the product.

	iterator old_it=it;
	it=tr.wrap(it, str_node("\\prod"));

	bool all_arguments_moved_out=true;
	sibling_iterator acton=tr.begin(old_it);
	while(acton!=tr.end(old_it)) {
		// Only look at child nodes which are not indices.
		if(acton->is_index()==false) {
			sibling_iterator derarg=acton;
			++acton; // don't use this anymore this loop

			if(*derarg->name=="\\sum") {
				all_arguments_moved_out=false;
				continue; // FIXME: Don't know how to handle this yet.
				}

#ifdef DEBUG
			std::cerr << "doing " << *derarg->name << std::endl;
#endif

			// If the argument of the derivative is not a product, make
			// into one, so we can handle everything using the same code.
			if(*derarg->name!="\\prod")
				derarg=tr.wrap(derarg, str_node("\\prod"));

			// Iterate over all arguments of the product sitting inside
			// the derivative (but see the comment above).
			sibling_iterator factor=tr.begin(derarg);
			while(factor!=tr.end(derarg)) {
				// std::cerr << "checking " << Ex(factor) << std::endl;

				sibling_iterator nxt=factor;
				++nxt;
				bool move_out=true;

				// An object pattern like 'A??' should always be assumed to have
				// dependence, because we don't yet know what it will match.
				if(move_out) {
					if(factor->is_name_wildcard() || factor->is_object_wildcard())
						move_out=false;
					}

				if(move_out) {
					if(is_diracbar)
						if(kernel.properties.get<Spinor>(factor) || kernel.properties.get<GammaMatrix>(factor))
							move_out=false;
					}

				// Then figure out whether there is implicit dependence on the operator.
				// or on the coordinate.
				if(move_out) {
					const DependsBase *dep=kernel.properties.get<DependsBase>(factor);
					if(dep!=0) {
#ifdef DEBUG
						std::cerr << *factor->name << " acted on by " << *old_it->name << "; depends" << std::endl;
#endif
						auto derivative=kernel.properties.get<Derivative>(old_it);
						Ex deps=dep->dependencies(kernel, factor /* it */);
						sibling_iterator depobjs=deps.begin(deps.begin());
						while(depobjs!=deps.end(deps.begin())) {
#ifdef DEBUG
							std::cerr << "?" << *old_it->name << " == " << *depobjs->name << std::endl;
#endif
							// FIXME: need to compare more than the name
							if(old_it->name == depobjs->name) {
								move_out=false;
								break;
								}
							else if(derivative && derivative->with_respect_to.size()>0 && derivative->with_respect_to.begin()->name==depobjs->name) {
								move_out=false;
								break;
								}
							else {
								// compare all indices
#ifdef DEBUG
								std::cerr << "comparing indices" << std::endl;
#endif
								sibling_iterator indit=tr.begin(old_it);
								while(indit!=tr.end(old_it)) {
									if(indit->is_index()) {
#ifdef DEBUG
										std::cerr << "compare " << *indit->name << " to " << *depobjs->name << std::endl;
#endif
										if(subtree_compare(&kernel.properties, indit, depobjs, 0, false, 0, false)==0) {
#ifdef DEBUG
											std::cerr << "not moving out" << std::endl;
#endif
											move_out=false;
											break;
											}
										}
									++indit;
									}
								if(!move_out) break;
								}
							++depobjs;
							}
						}
					}

				// Finally, there may also be explicit dependence.
				if(move_out) {
					// FIXME: This certainly does not handle Y(a,b) correctly
					sibling_iterator chldit=tr.begin(factor);
					while(chldit!=tr.end(factor)) {
						if(chldit->is_index()==false) {
							sibling_iterator indit=tr.begin(old_it);
							while(indit!=tr.end(old_it)) {
								if(subtree_exact_equal(&kernel.properties, chldit, indit, 0)) {
									move_out=false;
									break;
									}
								++indit;
								}
							if(!move_out) break;
							}
						++chldit;
						}
					}

				// If no dependence found, move this child out of the derivative.
				if(move_out) {
					// FIXME: Does does not look at the commutativity
					// property of the index/indices wrt. which the
					// derivative is taken.
					int sign=1;
					if(factor!=tr.begin(derarg)) {
						Ex_comparator compare(kernel.properties);
						sign=compare.can_move_to_front(tr, derarg, factor); //, Ex_comparator::match_t::no_match_less);
						}
					if(sign!=0) {
						// If the sign *is* zero, it means that we are trying to move a factor
						// through another one, but do not know how to do that. In that case,
						// you cannot move the factor out.
						res=result_t::l_applied;
						tr.move_before(old_it, factor);
						multiply(it->multiplier, sign);
						}
					}

				factor=nxt;
				}

			// std::cerr << "after step " << Ex(it) << std::endl;

			// All factors in this argument have been handled now, let's see what's left.
			unsigned int derarg_num_chldr=tr.number_of_children(derarg);
			if(derarg_num_chldr==0) {
				// Empty accents should simply be ignored, but empty derivatives vanish.
				if(!is_accent) {
					zero(it->multiplier);
					break; // we can stop now, the entire expression is zero.
					}
				}
			else {
				all_arguments_moved_out=false;
				if(derarg_num_chldr==1) {
					derarg=tr.flatten_and_erase(derarg);
					}
				}
			}
		else ++acton;
		}


	// All non-index arguments have now been handled.
	if(all_arguments_moved_out && is_accent) {
		zero(it->multiplier);
		}
	else if(*it->multiplier!=0) {
		if(tr.number_of_children(it)==1) { // nothing was moved out
			tr.flatten(it);
			it=tr.erase(it);
			}
		else {
			// Moving factors around has potentially led to a top-level product
			// which contains children with non-unit multiplier.
			cleanup_dispatch(kernel, tr, it);

			// If the derivative acts on another derivative, we need
			// to un-nest the argument of the outer (and this situation
			// can only happen if there is only one non-index child node)
			iterator itarg=tr.begin(it);
			while(itarg->is_index())
				++itarg;

			cleanup_dispatch(kernel, tr, itarg);
			}
		}
	cleanup_dispatch(kernel, tr, it);
	
//	std::cerr << "unwrap done " << it << std::endl;

	return res;
	}