ttNetwork.cpp 59.1 KB
Newer Older
Sebastian Wolf's avatar
Sebastian Wolf committed
1
// Xerus - A General Purpose Tensor Library
2
// Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf. 
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// 
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
// 
// Xerus is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org 
// or contact us at contact@libXerus.org.

20
/**
21 22 23
* @file
* @brief Implementation of the TTNetwork class (and thus TTTensor and TTOperator).
*/
24

Sebastian Wolf's avatar
Sebastian Wolf committed
25
#include <algorithm>
26
#include <memory>
Sebastian Wolf's avatar
Sebastian Wolf committed
27

28
#include <xerus/ttNetwork.h>
Sebastian Wolf's avatar
Sebastian Wolf committed
29

Sebastian Wolf's avatar
Sebastian Wolf committed
30
#include <xerus/misc/check.h>
31
#include <xerus/misc/math.h>
Sebastian Wolf's avatar
Sebastian Wolf committed
32
#include <xerus/misc/performanceAnalysis.h>
33
#include <xerus/misc/internal.h>
Sebastian Wolf's avatar
Sebastian Wolf committed
34

35
#include <xerus/basic.h>
36
#include <xerus/misc/basicArraySupport.h>
37
#include <xerus/index.h>
Sebastian Wolf's avatar
Sebastian Wolf committed
38
#include <xerus/tensor.h>
Sebastian Wolf's avatar
Sebastian Wolf committed
39
#include <xerus/ttStack.h>
40
#include <xerus/indexedTensorList.h>
41
#include <xerus/indexedTensorMoveable.h>
42

43
namespace xerus {
Sebastian Wolf's avatar
Sebastian Wolf committed
44
	/*- - - - - - - - - - - - - - - - - - - - - - - - - - Constructors - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
Sebastian Wolf's avatar
Sebastian Wolf committed
45
	template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
46
	TTNetwork<isOperator>::TTNetwork() : TensorNetwork(), canonicalized(false) {}
Sebastian Wolf's avatar
Sebastian Wolf committed
47 48
	
	
49 50 51 52 53
	template<bool isOperator>
	TTNetwork<isOperator>::TTNetwork(const Tensor& _tensor, const double _eps, const size_t _maxRank) :
		TTNetwork(_tensor, _eps, std::vector<size_t>(_tensor.degree() == 0 ? 0 : _tensor.degree()/N-1, _maxRank)) {}
	
	
54
	template<bool isOperator>
55 56 57
	TTNetwork<isOperator>::TTNetwork(const size_t _degree) : TTNetwork(std::vector<size_t>(_degree, 1)) { }
	
	template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
58
	TTNetwork<isOperator>::TTNetwork(Tensor::DimensionTuple _dimensions) : TensorNetwork(ZeroNode::None), canonicalized(true), corePosition(0) {
59 60 61
		dimensions = std::move(_dimensions);
		
		REQUIRE(dimensions.size()%N==0, "Illegal degree for TTOperator.");
Sebastian Wolf's avatar
Sebastian Wolf committed
62
		REQUIRE(!misc::contains(dimensions, size_t(0)), "Zero is no valid dimension.");
63
		const size_t numComponents = dimensions.size()/N;
64 65
		
		if (numComponents == 0) {
66
			nodes.emplace_back(std::make_unique<Tensor>());
67 68 69
			return;
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
70
		// ExternalLinks
71 72 73
		externalLinks.reserve(dimensions.size());
		for (size_t i = 0; i < numComponents; ++i) {
			externalLinks.emplace_back(i+1, 1, dimensions[i], false);
74
		}
Sebastian Wolf's avatar
Sebastian Wolf committed
75
		
76
		if (isOperator) {
77 78
			for (size_t i = 0; i < numComponents; ++i) {
				externalLinks.emplace_back(i+1, 2, dimensions[numComponents+i], false);
79 80 81
			}
		}
		
82
		std::vector<TensorNetwork::Link> neighbors;
83
		
84
		neighbors.emplace_back(1, 0, 1, false);
85
		
86
		nodes.emplace_back(std::make_unique<Tensor>(Tensor::ones({1})), std::move(neighbors));
Sebastian Wolf's avatar
Sebastian Wolf committed
87
		
Sebastian Wolf's avatar
Sebastian Wolf committed
88
		for (size_t i = 0; i < numComponents; ++i) {
89
			neighbors.clear();
90 91 92
			neighbors.emplace_back(i, i==0 ? 0 : N+1, 1, false);
			neighbors.emplace_back(-1, i, dimensions[i], true);
			if(isOperator) { neighbors.emplace_back(-1, i+numComponents, dimensions[numComponents+i], true); }
93
			neighbors.emplace_back(i+2, 0, 1, false);
94
			
95
			if(!isOperator) {
96
				nodes.emplace_back( std::make_unique<Tensor>(Tensor::dirac({1, dimensions[i], 1}, 0)), std::move(neighbors) );
97
			} else {
98
				nodes.emplace_back( std::make_unique<Tensor>(Tensor::dirac({1, dimensions[i], dimensions[numComponents+i], 1}, 0)), std::move(neighbors) );
99
			}
100
		}
Sebastian Wolf's avatar
Sebastian Wolf committed
101
		
102 103
		neighbors.clear();
		neighbors.emplace_back(numComponents, N+1, 1, false);
104
		nodes.emplace_back( std::make_unique<Tensor>(Tensor::ones({1})), std::move(neighbors));
105 106 107
		
		// Make a Zero Tensor (at core)
		(*nodes[1].tensorObject)[0] = 0;
108 109
	}
	
110
	
111
	template<bool isOperator>
112
	TTNetwork<isOperator>::TTNetwork(const Tensor& _tensor, const double _eps, const TensorNetwork::RankTuple& _maxRanks): TTNetwork(_tensor.degree()) {
113
		REQUIRE(_tensor.degree()%N==0, "Number of indicis must be even for TTOperator");
114
		REQUIRE(_eps >= 0 && _eps < 1, "_eps must be positive and smaller than one. " << _eps << " was given.");
Sebastian Wolf's avatar
Sebastian Wolf committed
115
		REQUIRE(_maxRanks.size() == num_ranks(), "We need " << num_ranks() <<" ranks but " << _maxRanks.size() << " where given");
Sebastian Wolf's avatar
Sebastian Wolf committed
116
		REQUIRE(!misc::contains(_maxRanks, size_t(0)), "Maximal ranks must be strictly positive. Here: " << _maxRanks);
117
		
118
		const size_t numComponents = degree()/N;
119
		
120
		if (_tensor.degree() == 0) {
121 122
			*nodes[0].tensorObject = _tensor;
			return;
123 124
		}
		
125
		dimensions = _tensor.dimensions;
126
		
127 128 129 130 131 132
		Tensor remains;
		if(isOperator) {
			std::vector<size_t> shuffle(_tensor.degree());
			for(size_t i = 0; i < numComponents; ++i) {
				shuffle[i] = 2*i;
				shuffle[numComponents + i] = 2*i+1; 
133
			}
Benjamin Huber's avatar
Benjamin Huber committed
134
			
135
			xerus::reshuffle(remains, _tensor, shuffle);
136
		} else {
137
			remains = _tensor;
138 139
		}
		
140 141 142 143 144 145 146 147
		// Add ghost dimensions used in the nodes
		std::vector<size_t> extDimensions;
		extDimensions.reserve(remains.degree()+2);
		extDimensions.emplace_back(1);
		extDimensions.insert(extDimensions.end(), remains.dimensions.begin(), remains.dimensions.end());
		extDimensions.emplace_back(1);
		remains.reinterpret_dimensions(extDimensions);
		
Sebastian Wolf's avatar
Sebastian Wolf committed
148
		
149 150 151 152
		Tensor singularValues, newNode;
		for(size_t position = numComponents-1; position > 0; --position) {
			calculate_svd(remains, singularValues, newNode, remains, 1+position*N, _maxRanks[position-1], _eps);
			
Sebastian Wolf's avatar
Sebastian Wolf committed
153 154
			set_component(position, std::move(newNode)); 
			newNode.reset();
155 156
			xerus::contract(remains, remains, false, singularValues, false, 1);
		}
157
		
158 159
		set_component(0, remains);
		assume_core_position(0);
160
	}
161 162
	

163 164 165 166
// 	template<bool isOperator>
// 	TTNetwork<isOperator>::TTNetwork(const TensorNetwork &_network, double _eps) : TTNetwork(Tensor(_network)) {
// 		LOG(warning, "Cast of arbitrary tensor network to TT not yet supported. Casting to Tensor first"); // TODO
// 	}
167

Sebastian Wolf's avatar
Sebastian Wolf committed
168
	
169 170
	template<bool isOperator>
	TTNetwork<isOperator> TTNetwork<isOperator>::ones(const std::vector<size_t>& _dimensions) {
171
		REQUIRE(_dimensions.size()%N == 0, "Illegal number of dimensions for ttOperator");
Sebastian Wolf's avatar
Sebastian Wolf committed
172
		REQUIRE(!misc::contains(_dimensions, size_t(0)), "Trying to construct a TTTensor with dimension 0 is not possible.");
173
		
174 175 176 177
		if(_dimensions.empty()) {
			return TTNetwork(Tensor::ones({}));
		}
		
178 179
		TTNetwork result(_dimensions.size());
		const size_t numNodes = _dimensions.size()/N;
180
		
Sebastian Wolf's avatar
Sebastian Wolf committed
181
		std::vector<size_t> dimensions(isOperator ? 4 : 3, 1);
182
		for(size_t i = 0; i < numNodes; ++i) {
Sebastian Wolf's avatar
Sebastian Wolf committed
183
			dimensions[1] = _dimensions[i];
184
			if (isOperator) {
Sebastian Wolf's avatar
Sebastian Wolf committed
185
				dimensions[2] = _dimensions[i+numNodes];
186 187 188
			}
			result.set_component(i, Tensor::ones(dimensions));
		}
Sebastian Wolf's avatar
Sebastian Wolf committed
189
		result.canonicalize_left();
190 191
		return result;
	}
192
	
Sebastian Wolf's avatar
Sebastian Wolf committed
193
	
Sebastian Wolf's avatar
Sebastian Wolf committed
194
	template<> template<>
195
	TTNetwork<true> TTNetwork<true>::identity(const std::vector<size_t>& _dimensions) {
196
		REQUIRE(_dimensions.size()%N==0, "Illegal number of dimensions for ttOperator");
Sebastian Wolf's avatar
Sebastian Wolf committed
197
		REQUIRE(!misc::contains(_dimensions, size_t(0)), "Trying to construct a TTTensor with dimension 0 is not possible.");
198
		
199 200 201 202
		if(_dimensions.empty()) {
			return TTNetwork(Tensor::ones({}));
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
203 204
		const size_t numComponents = _dimensions.size()/N;
		
Sebastian Wolf's avatar
Sebastian Wolf committed
205 206 207 208 209 210
		TTNetwork result(_dimensions.size());
		
		std::vector<size_t> constructionVector(4, 1);
		for (size_t i = 0; i < numComponents; ++i) {
			constructionVector[1] = _dimensions[i];
			constructionVector[2] = _dimensions[i+numComponents];
Sebastian Wolf's avatar
Sebastian Wolf committed
211
			result.set_component(i, Tensor(constructionVector, [](const std::vector<size_t> &_idx){
212 213 214
				if (_idx[1] == _idx[2]) {
					return 1.0;
				}
215
				return 0.0;
Sebastian Wolf's avatar
Sebastian Wolf committed
216
			}));
217 218
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
219
		result.canonicalize_left();
220 221 222
		return result;
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
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
	template<bool isOperator>
	TTNetwork<isOperator> TTNetwork<isOperator>::kronecker(const std::vector<size_t>& _dimensions) {
		REQUIRE(_dimensions.size()%N == 0, "Illegal number of dimensions for ttOperator");
		REQUIRE(!misc::contains(_dimensions, size_t(0)), "Trying to construct a TTNetwork with dimension 0 is not possible.");
		
		if(_dimensions.empty()) {
			return TTNetwork(Tensor::kronecker({}));
		}
		
		TTNetwork result(_dimensions.size());
		const size_t numNodes = _dimensions.size()/N;
		
		const auto minN = misc::min(_dimensions);
		
		// All nodes are simply kronecker tensors themself
		std::vector<size_t> dimensions;
		for(size_t i = 0; i < numNodes; ++i) {
			dimensions.reserve(4);
			if(i > 0) { dimensions.push_back(minN); }
			dimensions.push_back(_dimensions[i]);
			if (isOperator) { dimensions.push_back(_dimensions[i+numNodes]); }
			if(i+1 < numNodes) { dimensions.push_back(minN); }
			auto newComp = Tensor::kronecker(dimensions);
			if(i == 0) { dimensions.insert(dimensions.begin(), 1); }
			if(i+1 == numNodes) { dimensions.insert(dimensions.end(), 1); }
			if(i == 0 || i+1 == numNodes) { newComp.reinterpret_dimensions(std::move(dimensions)); }
			result.set_component(i, std::move(newComp));
			dimensions.clear();
		}
Sebastian Wolf's avatar
Sebastian Wolf committed
252
		result.canonicalize_left();
Sebastian Wolf's avatar
Sebastian Wolf committed
253 254 255
		return result;
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278
	template<bool isOperator>
	TTNetwork<isOperator> TTNetwork<isOperator>::dirac(std::vector<size_t> _dimensions, const std::vector<size_t>& _position) {
		REQUIRE(_dimensions.size()%N==0, "Illegal number of dimensions for ttOperator");
		REQUIRE(!misc::contains(_dimensions, size_t(0)), "Trying to construct a TTTensor with dimension 0 is not possible.");
		REQUIRE(_dimensions.size() == _position.size(), "Inconsitend number of entries in _dimensions and _position.");
		
		const size_t numComponents = _dimensions.size()/N;
		
		if(numComponents <= 1) {
			return TTNetwork(Tensor::dirac(_dimensions, _position));
		}
		
		TTNetwork result(_dimensions);
		
		for (size_t i = 0; i < numComponents; ++i) {
			if(isOperator) {
				result.set_component(i, Tensor::dirac({1, result.dimensions[i], result.dimensions[numComponents+i], 1}, _position[i]*result.dimensions[numComponents+i] + _position[numComponents+i]));
			} else {
				result.set_component(i, Tensor::dirac({1, result.dimensions[i], 1}, _position[i]));
			}
		}
		return result;
	}
Sebastian Wolf's avatar
Sebastian Wolf committed
279
	
Sebastian Wolf's avatar
Sebastian Wolf committed
280 281 282 283
	template<bool isOperator>
	TTNetwork<isOperator> TTNetwork<isOperator>::dirac(std::vector<size_t> _dimensions, const size_t _position) {
		return dirac(_dimensions, Tensor::position_to_multiIndex(_position, _dimensions));
	}
Sebastian Wolf's avatar
Sebastian Wolf committed
284
	
285 286
	
	/*- - - - - - - - - - - - - - - - - - - - - - - - - - Internal helper functions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
287
	
288
	#ifndef XERUS_DISABLE_RUNTIME_CHECKS
289
		template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
290
		void TTNetwork<isOperator>::require_correct_format() const {
291 292
			require_valid_network(); // Network must at least be valid.
			
293
			const size_t numComponents = degree()/N;
294
			const size_t numNodes = degree() == 0 ? 1 : degree()/N + 2;
295
			REQUIRE(nodes.size() == numNodes, "Wrong number of nodes: " << nodes.size() << " expected " << numNodes << ".");
Sebastian Wolf's avatar
Sebastian Wolf committed
296
			REQUIRE(!canonicalized || (degree() == 0 && corePosition == 0) || corePosition < numComponents, "Invalid corePosition: " << corePosition << " there are only " << numComponents << " components.");
297
			
298
			// Per external link
Sebastian Wolf's avatar
Sebastian Wolf committed
299
			for (size_t n = 0; n < externalLinks.size(); ++n) {
300
				const TensorNetwork::Link &l = externalLinks[n];
301
				REQUIRE(l.other == (n%numComponents)+1, "The " << n << "-th external link must point the the " << (n%numComponents) << "-th component (i.e. the " << (n%numComponents)+1 << "-th node) but does point to the " << l.other << "-th node.");
302 303
			}
			
304
			// Virtual nodes
305
			if(degree() > 0) {
306 307 308 309 310 311
				REQUIRE(nodes.front().degree() == 1, "The left virtual node must have degree 1, but has size " << nodes.front().degree());
				REQUIRE(nodes.front().neighbors[0].dimension == 1, "The left virtual node's single dimension must be 1, but is " << nodes.front().neighbors[0].dimension);
				REQUIRE(nodes.front().neighbors[0].other == 1, "The left virtual node's single link must be to node 1, but is towards node " << nodes.front().neighbors[0].other);
				REQUIRE(nodes.front().neighbors[0].indexPosition == 0, "The left virtual node's single link must link at indexPosition 0, but link at " << nodes.front().neighbors[0].indexPosition);
				REQUIRE(misc::hard_equal((*nodes.front().tensorObject)[0], 1.0), "The left virtual node's single entry must be 1.0, but it is " << (*nodes.front().tensorObject)[0]);
				REQUIRE(!nodes.front().tensorObject->has_factor(), "The left virtual node must no carry a non-trivial factor.");
312
				
313 314 315 316 317 318
				REQUIRE(nodes.back().degree() == 1, "The right virtual node must have degree 1, but has size " << nodes.back().degree());
				REQUIRE(nodes.back().neighbors[0].dimension == 1, "The right virtual node's single dimension must be 1, but is " << nodes.back().neighbors[0].dimension);
				REQUIRE(nodes.back().neighbors[0].other == numNodes-2, "The right virtual node's single link must be to node " << numNodes-2 << ", but is towards node " << nodes.back().neighbors[0].other);
				REQUIRE(nodes.back().neighbors[0].indexPosition == N+1, "The right virtual node's single link must link at indexPosition " << N+1 << ", but link at " << nodes.back().neighbors[0].indexPosition);
				REQUIRE(misc::hard_equal((*nodes.back().tensorObject)[0], 1.0), "The right virtual node's single entry must be 1.0, but it is " << (*nodes.back().tensorObject)[0]);
				REQUIRE(!nodes.back().tensorObject->has_factor(), "The right virtual node must no carry a non-trivial factor.");
319 320
			}
			
321
			// Per component
Sebastian Wolf's avatar
Sebastian Wolf committed
322
			for (size_t n = 0; n < numComponents; ++n) {
323 324
				const TensorNode& node = nodes[n+1];
				
Sebastian Wolf's avatar
Sebastian Wolf committed
325
				REQUIRE(!canonicalized || n == corePosition || !node.tensorObject->has_factor(), "In canonicalized TTNetworks only the core may carry a non-trivial factor. Violated by component " << n << " factor: " << node.tensorObject->factor << " corepos: " << corePosition);
326 327 328 329 330 331 332 333 334 335 336 337 338 339
				
				REQUIRE(node.degree() == N+2, "Every TT-Component must have degree " << N+2 << ", but component " << n << " has degree " << node.degree());
				REQUIRE(!node.neighbors[0].external, "The first link of each TT-Component must not be external. Violated by component " << n);
				REQUIRE(node.neighbors[0].other == n, "The first link of each TT-Component must link to the previous node. Violated by component " << n << ", which instead links to node " << node.neighbors[0].other << " (expected " << n << ").");
				REQUIRE(node.neighbors[0].indexPosition == (n==0?0:N+1), "The first link of each TT-Component must link to the last last index of the previous node. Violated by component " << n << ", which instead links to index " << node.neighbors[0].indexPosition << " (expected " << (n==0?0:N+1) << ").");
				
				REQUIRE(node.neighbors[1].external, "The second link of each TT-Component must be external. Violated by component " << n << ".");
				REQUIRE(node.neighbors[1].indexPosition == n, "The second link of each TT-Component must link to the external dimension equal to the component position. Violated by component " << n << " which links at " << node.neighbors[1].indexPosition);
				REQUIRE(!isOperator || node.neighbors[2].external, "The third link of each TTO-Component must be external. Violated by component " << n << ".");
				REQUIRE(!isOperator || node.neighbors[2].indexPosition == numComponents+n, "The third link of each TTO-Component must link to the external dimension equal to the component position + numComponents. Violated by component " << n << " which links at " << node.neighbors[2].indexPosition << " (expected " << numComponents+n << ").");
				
				REQUIRE(!node.neighbors.back().external, "The last link of each TT-Component must not be external. Violated by component " << n);
				REQUIRE(node.neighbors.back().other == n+2, "The last link of each TT-Component must link to the next node. Violated by component " << n << ", which instead links to node " << node.neighbors.back().other << " (expected " << n+2 << ").");
				REQUIRE(node.neighbors.back().indexPosition == 0, "The last link of each TT-Component must link to the first index of the next node. Violated by component " << n << ", which instead links to index " << node.neighbors.back().indexPosition << " (expected 0).");
340 341 342 343
			}
		}
	#else
		template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
344
		void TTNetwork<isOperator>::require_correct_format() const { }
345
	#endif
346
	
Sebastian Wolf's avatar
Sebastian Wolf committed
347
	
348 349
	template<bool isOperator>
	bool TTNetwork<isOperator>::exceeds_maximal_ranks() const {
Sebastian Wolf's avatar
Sebastian Wolf committed
350 351
		const size_t numComponents = dimensions.size()/N;
		for (size_t i = 0; i < numComponents; ++i) {
352
			const Tensor& comp = get_component(i);
Sebastian Wolf's avatar
Sebastian Wolf committed
353
			const size_t extDim = isOperator ? comp.dimensions[1]*comp.dimensions[2] : comp.dimensions[1];
354
			if (comp.dimensions.front() > extDim * comp.dimensions.back() || comp.dimensions.back() > extDim * comp.dimensions.front()) {
355 356 357 358 359 360
				return true;
			}
		}
		return false;
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
361 362 363 364 365 366
	
	template<bool isOperator>
	size_t TTNetwork<isOperator>::num_ranks() const {
		return degree() == 0 ? 0 : degree()/N-1;
	}
	
367
	/*- - - - - - - - - - - - - - - - - - - - - - - - - - Miscellaneous - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
368 369
	
	template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
370
	std::vector<size_t> TTNetwork<isOperator>::reduce_to_maximal_ranks(std::vector<size_t> _ranks, const std::vector<size_t>& _dimensions) {
371
		const size_t numComponents = _dimensions.size()/N;
Benjamin Huber's avatar
Benjamin Huber committed
372
		REQUIRE(_dimensions.size()%N == 0, "invalid number of dimensions for TTOperator");
373
		REQUIRE(numComponents == _ranks.size()+1, "Invalid number of ranks ("<<_ranks.size()<<") or dimensions ("<<_dimensions.size()<<") given.");
374
		
375 376 377 378 379 380 381 382 383 384
		// Left to right sweep
		size_t currMax = 1;
		for (size_t i = 0; i+1 < numComponents; ++i) {
			currMax *= _dimensions[i];
			if (isOperator) { currMax *= _dimensions[numComponents+i]; }
			
			if (currMax < _ranks[i]) { 
				_ranks[i] = currMax;
			} else {
				currMax = _ranks[i];
385
			}
386 387 388 389 390 391 392 393 394 395 396 397
		}
	
		// Right to left sweep
		currMax = 1;
		for (size_t i = 1; i < numComponents; ++i) {
			currMax *= _dimensions[numComponents-i];
			if (isOperator) { currMax *= _dimensions[2*numComponents-i]; }
			
			if (currMax < _ranks[numComponents-i-1]) {
				_ranks[numComponents-i-1] = currMax;
			} else {
				currMax = _ranks[numComponents-i-1];
398 399 400
			}
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
401
		return _ranks;
402
	}
403
	
404
	
Benjamin Huber's avatar
Benjamin Huber committed
405 406
	template<bool isOperator>
	size_t TTNetwork<isOperator>::degrees_of_freedom(const std::vector<size_t> &_dimensions, const std::vector<size_t> &_ranks) {
407
		if (_dimensions.empty()) { return 1; }
Benjamin Huber's avatar
Benjamin Huber committed
408 409 410 411 412 413 414
		const size_t numComponents = _dimensions.size()/N;
		REQUIRE(_dimensions.size()%N == 0, "invalid number of dimensions for TTOperator");
		REQUIRE(numComponents == _ranks.size()+1, "Invalid number of ranks ("<<_ranks.size()<<") or dimensions ("<<_dimensions.size()<<") given.");
		size_t result = 0;
		for (size_t i=0; i<numComponents; ++i) {
			size_t component = i==0? 1 : _ranks[i-1];
			component *= _dimensions[i];
415 416
			if (isOperator) { component *= _dimensions[i+numComponents]; }
			if (i<_ranks.size()) { component *= _ranks[i]; }
Benjamin Huber's avatar
Benjamin Huber committed
417 418
			result += component;
		}
419 420
		for (const auto r : _ranks) {
			result -= misc::sqr(r);
Benjamin Huber's avatar
Benjamin Huber committed
421 422 423 424 425 426 427 428 429 430
		}
		return result;
	}
	
	template<bool isOperator>
	size_t TTNetwork<isOperator>::degrees_of_freedom() {
		return degrees_of_freedom(dimensions, ranks());
	}
	
	
431
	template<bool isOperator>
432 433 434
	void TTNetwork<isOperator>::fix_mode(const size_t _mode, const size_t _slatePosition) {
		REQUIRE(!isOperator, "fix_mode(), does not work for TTOperators, if applicable cast to TensorNetwork first");
		TensorNetwork::fix_mode(_mode, _slatePosition);
435 436
	}
	
437
	template<bool isOperator>
438 439
	void TTNetwork<isOperator>::resize_mode(const size_t _mode, const size_t _newDim, const size_t _cutPos) {
		TensorNetwork::resize_mode(_mode, _newDim, _cutPos);
Sebastian Wolf's avatar
Sebastian Wolf committed
440
		if(canonicalized && _newDim != corePosition) {
441 442
			const size_t oldCorePosition = corePosition;
			const size_t numComponents = degree()/N;
443
			move_core(_mode%numComponents);
444 445 446 447
			move_core(oldCorePosition);
		}
	}
	
448
	
449 450 451 452 453 454 455
	template<bool isOperator>
	void TTNetwork<isOperator>::use_dense_representations() {
		for (size_t i = 0; i < degree(); ++i) { 
			component(i).use_dense_representation(); 
		}
	}
	
456
	template<bool isOperator>
457
	Tensor& TTNetwork<isOperator>::component(const size_t _idx) {
458
		REQUIRE(_idx == 0 || _idx < degree()/N, "Illegal index " << _idx <<" in TTNetwork::component, as there are onyl " << degree()/N << " components.");
459
		return *nodes[degree() == 0 ? 0 : _idx+1].tensorObject;
460
	}
461
	
462
	
463
	template<bool isOperator>
464
	const Tensor& TTNetwork<isOperator>::get_component(const size_t _idx) const {
465 466
		REQUIRE(_idx == 0 || _idx < degree()/N, "Illegal index " << _idx <<" in TTNetwork::get_component.");
		return *nodes[degree() == 0 ? 0 : _idx+1].tensorObject;
467 468
	}
	
469
	
470
	template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
471
	void TTNetwork<isOperator>::set_component(const size_t _idx, Tensor _T) {
472 473 474 475 476 477
		if(degree() == 0) {
			REQUIRE(_idx == 0, "Illegal index " << _idx <<" in TTNetwork::set_component");
			REQUIRE(_T.degree() == 0, "Component of degree zero TTNetwork must have degree zero. Given: " << _T.degree());
			*nodes[0].tensorObject = std::move(_T);
		} else {
			REQUIRE(_idx < degree()/N, "Illegal index " << _idx <<" in TTNetwork::set_component");
Sebastian Wolf's avatar
Sebastian Wolf committed
478
            REQUIRE(_T.degree() == N+2, "Component " << _idx << " must have degree " << N+2 << ". Given: " << _T.degree());
479 480 481 482 483 484 485 486 487
			
			TensorNode& currNode = nodes[_idx+1];
			*currNode.tensorObject = std::move(_T);
			for (size_t i = 0; i < N+2; ++i) {
				currNode.neighbors[i].dimension = currNode.tensorObject->dimensions[i];
				if (currNode.neighbors[i].external) {
					externalLinks[currNode.neighbors[i].indexPosition].dimension = currNode.tensorObject->dimensions[i];
					dimensions[currNode.neighbors[i].indexPosition] = currNode.tensorObject->dimensions[i];
				}
488 489
			}
		}
Sebastian Wolf's avatar
Sebastian Wolf committed
490
		
Sebastian Wolf's avatar
Sebastian Wolf committed
491
		canonicalized = canonicalized && (corePosition == _idx);
492 493
	}
	
494
	
495 496 497 498 499 500 501 502
	template<bool isOperator>
	std::vector<std::vector<std::tuple<size_t, size_t, value_t>>> get_grouped_entries(const Tensor& _component) {
		REQUIRE(_component.is_sparse(), "Not usefull (and not implemented) for dense Tensors.");
		
		const size_t externalDim = isOperator ? _component.dimensions[1] * _component.dimensions[2] : _component.dimensions[1];
		
		std::vector<std::vector<std::tuple<size_t, size_t, value_t>>> groups(externalDim);
		
503
		for(const auto& entry : _component.get_unsanitized_sparse_data()) {
504 505 506 507 508 509 510 511
			const size_t r2 = entry.first%_component.dimensions.back();
			const size_t n = (entry.first/_component.dimensions.back())%externalDim;
			const size_t r1 = (entry.first/_component.dimensions.back())/externalDim;
			groups[n].emplace_back(r1, r2, _component.factor*entry.second);
		}
		
		return groups;
	}
512
	
513
	
514 515
	template<bool isOperator>
	std::pair<TensorNetwork, TensorNetwork> TTNetwork<isOperator>::chop(const size_t _position) const {
Sebastian Wolf's avatar
Sebastian Wolf committed
516
		require_correct_format();
517 518 519 520 521
		
		const size_t numComponents = degree()/N;
		REQUIRE(_position < numComponents, "Can't split a " << numComponents << " component TTNetwork at position " << _position);
		
		// Create the resulting TNs
Sebastian Wolf's avatar
Sebastian Wolf committed
522 523
		TensorNetwork left(ZeroNode::None);
		TensorNetwork right(ZeroNode::None);
524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
		
		left.nodes.push_back(nodes[0]);
		for (size_t i = 0; i < _position; ++i) {
			left.dimensions.push_back(dimensions[i]);
			left.externalLinks.push_back(externalLinks[i]);
			left.nodes.push_back(nodes[i+1]);
		}
		if(isOperator) {
			for(size_t i = 0; i < _position; ++i) {
				left.dimensions.push_back(dimensions[i+numComponents]);
				left.externalLinks.push_back(externalLinks[i+numComponents]);
			}
		}
		left.dimensions.push_back(left.nodes.back().neighbors.back().dimension);
		left.externalLinks.emplace_back(_position, _position==0?0:N+1, left.nodes.back().neighbors.back().dimension , false);
		left.nodes.back().neighbors.back().external = true;
		left.nodes.back().neighbors.back().indexPosition = isOperator ? 2*_position-1 : _position;
Sebastian Wolf's avatar
Sebastian Wolf committed
541
		
542 543
		right.dimensions.push_back(nodes[_position+2].neighbors.front().dimension);
		right.externalLinks.emplace_back(_position+2, 0, nodes[_position+2].neighbors.front().dimension , false); // NOTE other will be corrected to 0 in the following steps
544 545

		for(size_t i = _position+1; i < numComponents; ++i) {
546 547 548 549 550 551 552 553 554 555
			right.dimensions.push_back(dimensions[i]);
			right.externalLinks.push_back(externalLinks[i]);
			right.nodes.push_back(nodes[i+1]);
		}
		if(isOperator) {
			for(size_t i = _position+1; i < numComponents+1; ++i) {
				right.dimensions.push_back(dimensions[i+numComponents]);
				right.externalLinks.push_back(externalLinks[i+numComponents]);
			}
		}
556 557 558
		// The last node
		right.nodes.push_back(nodes.back());
		
559 560 561 562
		right.nodes.front().neighbors.front().external = true;
		right.nodes.front().neighbors.front().indexPosition = _position; // NOTE indexPosition will be corrected to 0 in the following steps
		
		// Account for the fact that the first _position+2 nodes do not exist
563
		for(TensorNetwork::Link& link : right.externalLinks) {
564 565 566 567
			link.other -= _position+2;
		}
		
		for(TensorNode& node : right.nodes) {
568
			for(TensorNetwork::Link& link : node.neighbors) {
569
				if(link.external) {
570
					link.indexPosition -= _position;
571 572 573 574 575 576 577 578 579
				} else {
					link.other -= _position+2;
				}
			}
		}
		
		return std::pair<TensorNetwork, TensorNetwork>(std::move(left), std::move(right));
	}
	
580 581 582 583
	
	template<bool isOperator>
	void TTNetwork<isOperator>::move_core(const size_t _position, const bool _keepRank) {
		const size_t numComponents = degree()/N;
584
		REQUIRE(_position < numComponents || (_position == 0 && degree() == 0), "Illegal core-position " << _position << " chosen for TTNetwork with " << numComponents << " components");
585 586
		require_correct_format();
		
Sebastian Wolf's avatar
Sebastian Wolf committed
587
		if(canonicalized) {
588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603
			// Move right?
			for (size_t n = corePosition; n < _position; ++n) {
				transfer_core(n+1, n+2, !_keepRank);
			}
			
			// Move left?
			for (size_t n = corePosition; n > _position; --n) {
				transfer_core(n+1, n, !_keepRank);
			}
		} else {
			// Move right?
			for (size_t n = 0; n < _position; ++n) {
				transfer_core(n+1, n+2, !_keepRank);
			}
			
			// Move left?
604 605
			for (size_t n = numComponents; n > _position+1; --n) {
				transfer_core(n, n-1, !_keepRank);
606 607 608 609 610 611 612 613 614 615
			}
		}
		
		while (exceeds_maximal_ranks()) {
			// Move left from given CorePosition
			for (size_t n = _position; n > 0; --n) {
				transfer_core(n+1, n, !_keepRank);
			}
			
			// Move to the most right
616
			for (size_t n = 0; n+1 < numComponents; ++n) {
617 618 619 620
				transfer_core(n+1, n+2, !_keepRank);
			}
			
			// Move back left to given CorePosition
621 622
			for (size_t n = numComponents; n > _position+1; --n) {
				transfer_core(n, n-1, !_keepRank);
623 624 625
			}
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
626
		canonicalized = true;
627 628 629 630 631
		corePosition = _position;
	}
	
	
	template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
632
	void TTNetwork<isOperator>::canonicalize_left() {
633 634 635 636 637
		move_core(0);
	}
	
	
	template<bool isOperator>
Sebastian Wolf's avatar
Sebastian Wolf committed
638
	void TTNetwork<isOperator>::canonicalize_right() {
639
		move_core(degree() == 0 ? 0 : degree()/N-1);
640 641 642
	}
	
	
643
	template<bool isOperator>
644
	void TTNetwork<isOperator>::round(const std::vector<size_t>& _maxRanks, const double _eps) {
645
		require_correct_format();
646 647
		const size_t numComponents = degree()/N;
		REQUIRE(_eps < 1, "_eps must be smaller than one. " << _eps << " was given.");
648
		REQUIRE(_maxRanks.size()+1 == numComponents || (_maxRanks.empty() && numComponents == 0), "There must be exactly degree/N-1 maxRanks. Here " << _maxRanks.size() << " instead of " << numComponents-1 << " are given.");
Sebastian Wolf's avatar
Sebastian Wolf committed
649
		REQUIRE(!misc::contains(_maxRanks, size_t(0)), "Trying to round a TTTensor to rank 0 is not possible.");
650
		
Sebastian Wolf's avatar
Sebastian Wolf committed
651
		const bool initialCanonicalization = canonicalized;
652 653
		const size_t initialCorePosition = corePosition;
		
Sebastian Wolf's avatar
Sebastian Wolf committed
654
		canonicalize_right();
655 656
		
		for(size_t i = 0; i+1 < numComponents; ++i) {
Ben Huber's avatar
Ben Huber committed
657
			round_edge(numComponents-i, numComponents-i-1, _maxRanks[numComponents-i-2], _eps, 0.0);
658 659
		}
		
660
		assume_core_position(0);
661 662 663
		
		if(initialCanonicalization) {
			move_core(initialCorePosition);
664
		}
665
	}
666
	
Sebastian Wolf's avatar
Sebastian Wolf committed
667
	
668
	template<bool isOperator>
669
	void TTNetwork<isOperator>::round(const size_t _maxRank) {
Sebastian Wolf's avatar
Sebastian Wolf committed
670
		round(std::vector<size_t>(num_ranks(), _maxRank), EPSILON);
671
	}
672
	
Sebastian Wolf's avatar
Sebastian Wolf committed
673
	
674
	template<bool isOperator>
675 676 677
	void TTNetwork<isOperator>::round(const int _maxRank) {
		REQUIRE( _maxRank > 0, "MaxRank must be positive");
		round(size_t(_maxRank));
678 679
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
680
	
681
	template<bool isOperator>
682
	void TTNetwork<isOperator>::round(const value_t _eps) {
Sebastian Wolf's avatar
Sebastian Wolf committed
683
		round(std::vector<size_t>(num_ranks(), std::numeric_limits<size_t>::max()), _eps);
684
	}
685

686
	
687
	template<bool isOperator>
Ben Huber's avatar
Ben Huber committed
688
	void TTNetwork<isOperator>::soft_threshold(const std::vector<double> &_taus, const bool  /*_preventZero*/) {
689
		const size_t numComponents = degree()/N;
690
		REQUIRE(_taus.size()+1 == numComponents || (_taus.empty() && numComponents == 0), "There must be exactly degree/N-1 taus. Here " << _taus.size() << " instead of " << numComponents-1 << " are given.");
Sebastian Wolf's avatar
Sebastian Wolf committed
691
		require_correct_format();
692
		
Sebastian Wolf's avatar
Sebastian Wolf committed
693
		const bool initialCanonicalization = canonicalized;
694 695
		const size_t initialCorePosition = corePosition;
		
Sebastian Wolf's avatar
Sebastian Wolf committed
696
		canonicalize_right();
697
		
698
		for(size_t i = 0; i+1 < numComponents; ++i) {
Ben Huber's avatar
Ben Huber committed
699
			round_edge(numComponents-i, numComponents-i-1, std::numeric_limits<size_t>::max(), 0.0, _taus[i]);
700 701
		}
		
702
		assume_core_position(0);
703
		
704 705 706 707
		if(initialCanonicalization) {
			move_core(initialCorePosition);
		}
	}
708
	
709
	
710
	template<bool isOperator>
711
	void TTNetwork<isOperator>::soft_threshold(const double _tau, const bool _preventZero) {
Sebastian Wolf's avatar
Sebastian Wolf committed
712
		soft_threshold(std::vector<double>(num_ranks(), _tau), _preventZero);
713
	}
714
	
715

716 717 718
	template<bool isOperator>
	std::vector<size_t> TTNetwork<isOperator>::ranks() const {
		std::vector<size_t> res;
Sebastian Wolf's avatar
Sebastian Wolf committed
719
		res.reserve(num_ranks());
720
		for (size_t n = 1; n+2 < nodes.size(); ++n) {
721 722 723 724 725
			res.push_back(nodes[n].neighbors.back().dimension);
		}
		return res;
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
726
	
727
	template<bool isOperator>
728
	size_t TTNetwork<isOperator>::rank(const size_t _i) const {
729
		REQUIRE(_i+1 < degree()/N, "Requested illegal rank " << _i);
730 731 732
		return nodes[_i+1].neighbors.back().dimension;
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
733
	
734 735
	template<bool isOperator>
	void TTNetwork<isOperator>::assume_core_position(const size_t _pos) {
736
		REQUIRE(_pos < degree()/N || (degree() == 0 && _pos == 0), "Invalid core position.");
737
		corePosition = _pos;
Sebastian Wolf's avatar
Sebastian Wolf committed
738
		canonicalized = true;
739 740
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
741
	
742 743 744 745 746
	template<bool isOperator>
	TensorNetwork* TTNetwork<isOperator>::get_copy() const {
		return new TTNetwork(*this);
	}
	
747 748 749 750 751 752
	template<bool isOperator>
	void TTNetwork<isOperator>::contract_unconnected_subnetworks() {
		if(degree() == 0) {
			std::set<size_t> all;
			for(size_t i = 0; i < nodes.size(); ++i) { all.emplace_hint(all.end(), i); }
			contract(all);
Sebastian Wolf's avatar
Sebastian Wolf committed
753
			canonicalized = false;
754 755 756 757
		} else {
			REQUIRE(nodes.size() > 2, "Invalid TTNetwork");
			const size_t numComponents = nodes.size()-2;
			
758 759 760 761
			for(size_t i = 0; i+1 < numComponents; ++i) {
				if(nodes[i+1].degree() == 2) {
						// If we are the core, everything is fine, we contract ourself to the next node, then get removed and the corePositions stays. If the next Node is the core, we have to change the corePosition to ours, because we will be removed. In all other cases cannonicalization is destroyed.
						if(corePosition == i+1) { corePosition = i; }
Sebastian Wolf's avatar
Sebastian Wolf committed
762
						else if(corePosition != i) { canonicalized = false; }
763
						contract(i+1, i+2);
764 765 766 767 768
				}
			}
			
			// Extra treatment for last component to avoid contraction to the pseudo-node.
			if(nodes[numComponents].degree() == 2) {
769
				if(corePosition == numComponents-1) { corePosition = numComponents-2; }
Sebastian Wolf's avatar
Sebastian Wolf committed
770
				else if(corePosition != numComponents-2) { canonicalized = false; }
771 772 773 774
				contract(numComponents-1, numComponents);
			}
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
775
		INTERNAL_CHECK(corePosition < degree() || !canonicalized, "Woot");
776
		
777 778 779
		sanitize();
	}