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
|
#pragma once
#include "Compare.hh"
/// \ingroup algorithms
///
/// Evaluate a tensorial expression to components, performing all
/// dummy index sums.
/**
See tests/components.cdb for basic samples.
Components nodes have the following structure:
\components_{a b ... n}({t,t,...r}=value, {t,t,...r}=value, ...)
\verbatim
\components
_{m} // index names/types
_{n}
\comma // last child
\equals
\comma // list of index values
r
t
[value] // value of tensor for above values of indices
\equals
...
\equals
...
ex:= A_{m n} ( A_{m n} + B_{m n} );
crds:= m -> { t, r, \phi, \theta };
vals:= { A_{t t} -> r, A_{r r} -> r**2, A_{t \phi} -> t, B_{t t} -> 1 };
evaluate(ex, crds, vals);
tst:= r * ( r + 1 ) + r**2 * r**2 + t*t - @(ex);
\endverbatim
The algorithm performs
the following steps. First, all free and dummy indices are
classified using members of Algorithm (maybe later?). For all indices, the 'crds'
argument to the algorithm is scanned for a
\verbatim
m -> { a,b,c }
\endverbatim
type pattern.
Do a dept-first scan until the first sum or product node.
In a sum node, all terms must have the same index structure, so
we can now make a map from the subtree pattern to all possible
index values. For each such index value set, lookup if there is a
replacement pattern, if not replace by zero. This map can be
reused later for equal sub-expressions (here it is important that
a normal tensor expression can be used as a pattern immediately,
though I would search both on explicit iterator (for things that
we have just seen) and on pattern (for things that occur again
later in the tree)).
In a product node, find all dummies and free indices. Create a
set of all index value sets, e.g.
\verbatim
{m, n, p} ->
{ {t, t, t}, {r,t,t}, {t,r,t}, ... }
\endverbatim
For each such index value set, lookup replacement rules.
KEY: quick way to lookup, for a given set of index values on a
(found) pattern/pattern pointed to by iterator, whether and what is the
value of the pattern.
*/
#include "Algorithm.hh"
#include "properties/Indices.hh"
namespace cadabra {
class evaluate : public Algorithm {
public:
evaluate(const Kernel&, Ex&, const Ex& component_values, bool rhs=false, bool simplify=true);
virtual bool can_apply(iterator) override;
virtual result_t apply(iterator&) override;
/// Merge the information in two 'components' nodes at the given
/// iterators, moving all out of the second one into the first
/// one.
void merge_components(iterator it1, iterator it2);
private:
const Ex& components;
bool only_rhs, call_sympy;
bool is_component(iterator it) const;
bool is_scalar_function(iterator it) const;
iterator handle_components(iterator it);
iterator handle_sum(iterator it);
iterator handle_prod(iterator it);
iterator handle_derivative(iterator it);
iterator handle_epsilon(iterator it);
/// Replace a single factor with a 'components' ...
/// The full_ind_free argument can contain a list of indices in the order
/// in which values should be stored in index value sets.
iterator handle_factor(sibling_iterator sib, const index_map_t& full_ind_free);
/// Expand a tensor factor into a components node with all components
/// written out explicitly. Used when there is no sparse rule matching
/// this factor.
iterator dense_factor(iterator sib, const index_map_t& ind_free, const index_map_t& ind_dummy);
/// Merge entries in a single 'components' node when they are for the
/// same index value(s).
void merge_component_children(iterator it);
/// Cleanup all components in a 'components' node; that is, call the
/// cleanup_dispatch function on them.
void cleanup_components(iterator it1);
/// Simplify all components of a 'components' node by running sympy's simplify
/// on them. Returns a replacement iterator to the top.
void simplify_components(iterator);
/// Wrap a non-component scalar node in a 'components' node.
iterator wrap_scalar_in_components_node(iterator sib);
/// Inverse of the above.
void unwrap_scalar_in_components_node(iterator sib);
/// Determine all the Coordinate dependencies of the object at 'it'. For the
/// time being this can only be a 'components' node.
std::set<Ex, tree_exact_less_obj> dependencies(iterator it);
};
}
|
