File: Adjform.hh

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 (229 lines) | stat: -rw-r--r-- 8,323 bytes parent folder | download | duplicates (2)
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
#pragma once

#include <vector>
#include <map>
#include <iosfwd>
#include <cstdint>
#include <string>
#include "Compare.hh"
#include "Kernel.hh"
#include "Storage.hh"

namespace cadabra {

	class IndexMap;

	/// Representation of the index structure of a tensor monomial,
	/// using a storage format which resembles an adjacency matrix. The
	/// structure is stored as a vector of integers.  Negative integers
	/// denote free indices. Positive indices denote a contraction with
	/// an index at the indicated position.
	///
	/// Example:
	///
	///    A_{m n} B_{p n}   ->   -1 3 -2 1
	///       0 1     2 3
	///
	/// The tensor names themselves ('A' and 'B' above) are not stored.
	
	class Adjform
		{
		public:
			/// The maximal number of index slots is set by `value_type`: for a
			/// short, the maximal number of slots is 127.
			using value_type = short;
			using size_type = value_type;
			using difference_type = value_type;
			using array_type = std::vector<value_type>;
			using const_reference = array_type::const_reference;
			using const_iterator = array_type::const_iterator;

			Adjform();
			template <typename IndexIterator>
			Adjform(IndexIterator beg, IndexIterator end, IndexMap& index_map, const Kernel& kernel);
			template <typename ValueIterator>
			Adjform(ValueIterator beg, ValueIterator end, bool push_as_coordinates);

			const_iterator begin() const;
			const_iterator end() const;

			// Return the position of the given index starting from the given
			// offset (but returning the position from the actual start). Returns
			// this->size() if not found.
			size_type index_of(value_type index, size_type offset = 0) const;

			bool operator < (const Adjform& other) const;
			bool operator == (const Adjform& other) const;
			bool operator != (const Adjform& other) const;

			const_reference operator [] (size_type idx) const;

			size_type size() const;
			size_type max_size() const;
			bool empty() const;

			// Return true if the value at 'pos' is < 0
			bool is_free_index(size_type pos) const;
			// Return true if the value at pos is >= 0
			bool is_dummy_index(size_type pos) const;
	
			size_type n_free_indices() const;
			size_type n_dummy_indices() const;

			// Contract pairs of this type of index, e.g. with value=-1 '-1 -2 -3 -1' -> '3 -2 -3 0'.
			// Returns true if the index was found and replaced.
			bool resolve_dummy(value_type value);

			// Push value(s) to the end and contract it if a eqivalent is found, e.g. with
			// value=-1 '-1 -2 -3' -> '3 -2 -3 0'
			void push_index(value_type value);
			void push_indices(const Adjform& other);

			// PUsh value to the end, but do not attempt to contract it, e.g. with
			// value=-1 '-1 -2 -3' -> '-1 -2 -3 -1'
			void push_coordinate(value_type value);
			void push_coordinates(const Adjform& other);

			// Push an iterator by calling IndexMap::is_coordinate to see if push_index or 
			// push_coordinate should be used
			void push(Ex::iterator iterator, IndexMap& index_map, const Kernel& kernel);

			// Swap the values at positions a and b
			void swap(size_type a, size_type b);

			// Cycle the values by moving the last element to the front n times,
			// e.g. with n=1 '3 -2 -3 0' -> '1 0 -2 -3'
			void rotate(size_type n);

			// Sort indices so that all free indices are at the beginning and all contracted
			// pairs at the end and next to each other, e.g. '3 -2 -3 0' -> '-3 -2 3 2'
			void sort();

			// Produce a unique integer to represent the current permutation
			uint64_t to_lehmer_code() const;
			// The total number of possible permutations of values
			uint64_t max_lehmer_code() const;
			// String representation where indices are named 'a-z' in the order they appear
			std::string to_string() const;

		private:
			// Storage of an index configuration; has a maximal size
			// equal to size_type (not size_t).
			array_type data;
		};

	/// To ensure consistency when creating adjforms out of two
	/// different Ex objects an IndexMap object is required which
	/// keeps track of which numeric index represents which index
	/// name.
	
	class IndexMap
		{
		public:
			IndexMap(const Kernel& kernel);
			~IndexMap();
			// Return a negative integer for each unique index
			Adjform::value_type get_free_index(Ex::iterator index);
			// Return true if index has the Coordinate or Symbol property, or is an integer
			static bool is_coordinate(const Kernel& kernel, Ex::iterator index);
		private:
			std::unique_ptr<Ex_comparator> comp;
			std::unique_ptr<Ex> data;
		};

	/// Representation of a sum of tensor monomials, each having the
	/// same tensor names, but with different index positions. As with
	/// AdjForm, the names of the tensors are not stored, only the
	/// index structure and the coefficient of each term.
	
	class ProjectedAdjform
		{
		public:
			using integer_type = int32_t;
			using map_t = std::map<Adjform, integer_type>;
			using iterator = map_t::iterator;
			using const_iterator = map_t::const_iterator;

			ProjectedAdjform();
			ProjectedAdjform(const Adjform& adjform, const integer_type& value = 1);

			// Add all contributions from 'other' into 'this'
			void combine(const ProjectedAdjform& other); 
			void combine(const ProjectedAdjform& other, integer_type factor);
			ProjectedAdjform& operator += (const ProjectedAdjform& other);
			friend ProjectedAdjform operator + (ProjectedAdjform lhs, const ProjectedAdjform& rhs);

			// Multiply all terms by a scalar factor
			void multiply(const integer_type& k);
			ProjectedAdjform& operator *= (const integer_type& k);
			friend ProjectedAdjform operator * (ProjectedAdjform lhs, const integer_type& rhs);

			iterator begin();
			const_iterator begin() const;
			iterator end();
			const_iterator end() const;

			void clear(); // Remove all entries
			size_t size() const; // Number of entries
			size_t max_size() const; // Returns the number of terms there would be if fully symmetrized
			size_t n_indices() const; // Returns the number of indices each adjform has
			bool empty() const; // True if there are no entries

			// Get the value of the term, or zero if it doesn't exist
			const integer_type& get(const Adjform& adjform) const;
			// Sets the given term to value, creating/removing the term if required
			void set(const Adjform& adjform, const integer_type& value = 1);
			// Adds value to the given term, creating/removing the term if required
			void add(const Adjform& adjform, const integer_type& value = 1);

			// Symmetrise or anti-symmetrise in the given indices
			// e.g. if the only term is abcd then
			//        apply_young_symmetry({0, 1, 2}, false) -> abcd + acbd + bacd + bcad + cabd + cbad
			//        apply_young_symmetry({2, 3, 4}, true) -> abcd - abdc - acbd + acdb - adcb + adbc
			void apply_young_symmetry(const std::vector<size_t>& indices, bool antisymmetric);

			// Symmetrize in indices starting at the indices in 'positions' with each group
			// 'n_indices' long.
			// e.g. if the only term is abcdefg then apply_ident_symmetry({0, 2, 4}, 2) ->
			//      abcdefg + abefcdg + cdabefg + cdefabg + efabcdg + efcdabg
			void apply_ident_symmetry(const std::vector<size_t>& positions, size_t n_indices);
			void apply_ident_symmetry(const std::vector<size_t>& positions, size_t n_indices, const std::vector<std::vector<int>>& commutation_matrix);

			// Symmetrize cyclically so abc -> abc + bca + cab
			void apply_cyclic_symmetry();

		private:
			// Unsafe (but faster) versions of the public functions
			void set_(const Adjform& adjform, const integer_type& value = 1);
			void add_(const Adjform& adjform, const integer_type& value = 1);

			map_t data;
			static integer_type zero;
		};

	template <typename IndexIterator>
	Adjform::Adjform(IndexIterator beg, IndexIterator end, IndexMap& index_map, const Kernel& kernel)
		{
		while (beg != end) {
			push(beg, index_map, kernel);
			++beg;
			}
		}

	template <typename ValueIterator>
	Adjform::Adjform(ValueIterator beg, ValueIterator end, bool push_as_coordinates)
		{
		while (beg != end) {
			auto val = *beg;
			if (push_as_coordinates)
				push_coordinate(val);
			else
				push_index(val);
			++beg;
			}
		}
	}


std::ostream& operator << (std::ostream& os, const cadabra::Adjform& adjform);
std::ostream& operator << (std::ostream& os, const cadabra::ProjectedAdjform& adjex);