math_algos.h 111 KB
Newer Older
1
2
3
4
/**
 * future math library, developed from scratch to eventually replace tlibs(2)
 * container-agnostic math algorithms
 * @author Tobias Weber <tweber@ill.fr>
5
 * @date dec-2017 -- 2019
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
 * @license GPLv3, see 'LICENSE' file
 * @desc The present version was forked on 8-Nov-2018 from the privately developed "magtools" project (https://github.com/t-weber/magtools).
 */

#ifndef __TOBIS_MATH_ALGOS_H__
#define __TOBIS_MATH_ALGOS_H__

#include <cstddef>
#include <cmath>
#include <cassert>
#include <complex>
#include <tuple>
#include <vector>
#include <limits>
#include <algorithm>
#include <iterator>
#include <numeric>
#include <iostream>
#include <iomanip>
#include <boost/algorithm/string.hpp>


// separator tokens
#define COLSEP ';'
#define ROWSEP '|'


namespace m {

Tobias WEBER's avatar
Tobias WEBER committed
35
36

// constants
37
38
39
40
41
42
43
template<typename T> constexpr T pi = T(M_PI);
template<typename T> T golden = T(0.5) + std::sqrt(T(5))/T(2);


// ----------------------------------------------------------------------------
// concepts
// ----------------------------------------------------------------------------
Tobias WEBER's avatar
Tobias WEBER committed
44
45
46
47
48
49
50
/**
 * requirements for having a value_type
 */
template<class T>
concept bool has_value_type = requires { typename T::value_type; };


51
52
53
54
/**
 * requirements for a scalar type
 */
template<class T>
55
concept bool is_scalar =
56
57
58
59
60
61
62
63
64
65
66
	std::is_floating_point_v<T> || std::is_integral_v<T> /*|| std::is_arithmetic_v<T>*/;


/**
 * requirements for a basic vector container like std::vector
 */
template<class T>
concept bool is_basic_vec = requires(const T& a)
{
	a.size();					// must have a size() member function
	a.operator[](1);			// must have an operator[]
Tobias WEBER's avatar
Tobias WEBER committed
67
} && has_value_type<T>;
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

/**
 * requirements of a vector type with a dynamic size
 */
template<class T>
concept bool is_dyn_vec = requires(const T& a)
{
	T(3);						// constructor
};

/**
 * requirements for a vector container
 */
template<class T>
concept bool is_vec = requires(const T& a)
{
	a+a;						// operator+
	a-a;						// operator-
	a[0]*a;						// operator*
	a*a[0];
	a/a[0];						// operator/
} && is_basic_vec<T>;


/**
 * requirements for a basic matrix container
 */
template<class T>
concept bool is_basic_mat = requires(const T& a)
{
	a.size1();					// must have a size1() member function
	a.size2();					// must have a size2() member function
	a.operator()(1,1);			// must have an operator()
Tobias WEBER's avatar
Tobias WEBER committed
101
} && has_value_type<T>;
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

/**
 * requirements of a matrix type with a dynamic size
 */
template<class T>
concept bool is_dyn_mat = requires(const T& a)
{
	T(3,3);						// constructor
};

/**
 * requirements for a matrix container
 */
template<class T>
concept bool is_mat = requires(const T& a)
{
	a+a;						// operator+
	a-a;						// operator-
	a(0,0)*a;					// operator*
	a*a(0,0);
	a/a(0,0);					// operator/
} && is_basic_mat<T>;




/**
 * requirements for a complex number
 */
template<class T>
concept bool is_complex = requires(const T& a)
{
	std::conj(a);
	a.real();			// must have a real() member function
	a.imag();			// must have an imag() member function

	a+a;
	a-a;
	a*a;
	a/a;
Tobias WEBER's avatar
Tobias WEBER committed
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
} && has_value_type<T>;
// ----------------------------------------------------------------------------





// ----------------------------------------------------------------------------
// type traits
// ----------------------------------------------------------------------------
/**
 * if it exists, get T's underlying value_type, otherwise use the type T directly
 */
template<class T, bool b = has_value_type<T>> struct _underlying_value_type {};
template<class T> struct _underlying_value_type<T, false> { using ty = T; };
template<class T> struct _underlying_value_type<T, true> { using ty = typename T::value_type; };
template<class T> using underlying_value_type = typename _underlying_value_type<T>::ty;
159
160
161
162
// ----------------------------------------------------------------------------



Tobias WEBER's avatar
Tobias WEBER committed
163
164
165
166
167
168
169
170
171
172
173
174

// ----------------------------------------------------------------------------
// forward declarations
// ----------------------------------------------------------------------------
template<class t_mat, class t_vec = std::vector<typename t_mat::value_type>>
std::tuple<bool, t_mat, t_mat> qr(const t_mat& mat)
requires is_mat<t_mat> && is_vec<t_vec>;
// ----------------------------------------------------------------------------




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
// ----------------------------------------------------------------------------
// adapters
// ----------------------------------------------------------------------------
template<typename size_t, size_t N, typename T, template<size_t, size_t, class...> class t_mat_base>
class qvec_adapter : public t_mat_base<1, N, T>
{
public:
	// types
	using base_type = t_mat_base<1, N, T>;
	using size_type = size_t;
	using value_type = T;

	// constructors
	using base_type::base_type;
	qvec_adapter(const base_type& vec) : base_type{vec} {}

	constexpr size_t size() const { return N; }

	T& operator[](size_t i) { return base_type::operator()(i,0); }
	const T operator[](size_t i) const { return base_type::operator()(i,0); }
};

template<typename size_t, size_t ROWS, size_t COLS, typename T, template<size_t, size_t, class...> class t_mat_base>
class qmat_adapter : public t_mat_base<COLS, ROWS, T>
{
public:
	// types
	using base_type = t_mat_base<COLS, ROWS, T>;
	using size_type = size_t;
	using value_type = T;

	// constructors
	using base_type::base_type;
	qmat_adapter(const base_type& mat) : base_type{mat} {}

	size_t size1() const { return ROWS; }
	size_t size2() const { return COLS; }
};


template<typename size_t, size_t N, typename T, class t_vec_base>
class qvecN_adapter : public t_vec_base
{
public:
	// types
	using base_type = t_vec_base;
	using size_type = size_t;
	using value_type = T;

	// constructors
	using base_type::base_type;
	qvecN_adapter(const base_type& vec) : base_type{vec} {}

	constexpr size_t size() const { return N; }

	T& operator[](size_t i) { return static_cast<base_type&>(*this)[i]; }
	const T operator[](size_t i) const { return static_cast<const base_type&>(*this)[i]; }
};

template<typename size_t, size_t ROWS, size_t COLS, typename T, class t_mat_base>
class qmatNN_adapter : public t_mat_base
{
public:
	// types
	using base_type = t_mat_base;
	using size_type = size_t;
	using value_type = T;

	// constructors
	using base_type::base_type;
	qmatNN_adapter(const base_type& mat) : base_type{mat} {}

Tobias WEBER's avatar
Tobias WEBER committed
247
248
249
250
251
252
253
254
255
	// convert from a different matrix type
	template<class t_matOther> qmatNN_adapter(const t_matOther& matOther)
		requires is_basic_mat<t_matOther>
	{
		const std::size_t minRows = std::min(static_cast<std::size_t>(size1()), static_cast<std::size_t>(matOther.size1()));
		const std::size_t minCols = std::min(static_cast<std::size_t>(size2()), static_cast<std::size_t>(matOther.size2()));

		for(std::size_t i=0; i<minRows; ++i)
			for(std::size_t j=0; j<minCols; ++j)
Tobias WEBER's avatar
Tobias WEBER committed
256
				(*this)(i,j) = static_cast<value_type>(matOther(i,j));
Tobias WEBER's avatar
Tobias WEBER committed
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
	size_t size1() const { return ROWS; }
	size_t size2() const { return COLS; }
};
// ----------------------------------------------------------------------------



// ----------------------------------------------------------------------------
// matrix
// ----------------------------------------------------------------------------

template<class T=double, template<class...> class t_cont = std::vector>
requires is_basic_vec<t_cont<T>> && is_dyn_vec<t_cont<T>>
class mat
{
public:
	using value_type = T;
	using container_type = t_cont<T>;

	mat() = default;
	mat(std::size_t ROWS, std::size_t COLS) : m_data(ROWS*COLS), m_rowsize(ROWS), m_colsize(COLS) {}
	~mat() = default;

	std::size_t size1() const { return m_rowsize; }
	std::size_t size2() const { return m_colsize; }
	const T& operator()(std::size_t row, std::size_t col) const { return m_data[row*m_colsize + col]; }
	T& operator()(std::size_t row, std::size_t col) { return m_data[row*m_colsize + col]; }

private:
	container_type m_data;
	std::size_t m_rowsize, m_colsize;
};
// ----------------------------------------------------------------------------



// ----------------------------------------------------------------------------
// n-dim algos
// ----------------------------------------------------------------------------

Tobias WEBER's avatar
Tobias WEBER committed
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
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
/**
 * are two scalars equal within an epsilon range?
 */
template<class T>
bool equals(T t1, T t2, T eps = std::numeric_limits<T>::epsilon())
requires is_scalar<T>
{
	return std::abs(t1 - t2) <= eps;
}

/**
 * are two complex numbers equal within an epsilon range?
 */
template<class T>
bool equals(const T& t1, const T& t2,
	typename T::value_type eps = std::numeric_limits<typename T::value_type>::epsilon())
requires is_complex<T>
{
	return (std::abs(t1.real() - t2.real()) <= eps) &&
		(std::abs(t1.imag() - t2.imag()) <= eps);
}

/**
 * are two vectors equal within an epsilon range?
 */
template<class t_vec>
bool equals(const t_vec& vec1, const t_vec& vec2,
	typename t_vec::value_type eps = std::numeric_limits<typename t_vec::value_type>::epsilon())
requires is_basic_vec<t_vec>
{
	using T = typename t_vec::value_type;

	// size has to be equal
	if(vec1.size() != vec2.size())
		return false;

	// check each element
	for(std::size_t i=0; i<vec1.size(); ++i)
	{
		if constexpr(is_complex<decltype(eps)>)
		{
			if(!equals<T>(vec1[i], vec2[i], eps.real()))
				return false;
		}
		else
		{
			if(!equals<T>(vec1[i], vec2[i], eps))
				return false;
		}
	}

	return true;
}

/**
 * are two matrices equal within an epsilon range?
 */
template<class t_mat>
bool equals(const t_mat& mat1, const t_mat& mat2,
	typename t_mat::value_type eps = std::numeric_limits<typename t_mat::value_type>::epsilon())
requires is_mat<t_mat>
{
	using T = typename t_mat::value_type;

	if(mat1.size1() != mat2.size1() || mat1.size2() != mat2.size2())
		return false;

	for(std::size_t i=0; i<mat1.size1(); ++i)
	{
		for(std::size_t j=0; j<mat1.size2(); ++j)
		{
			if(!equals<T>(mat1(i,j), mat2(i,j), eps))
				return false;
		}
	}

	return true;
}


/**
 * set submatrix to unit
 */
template<class t_mat>
void unit(t_mat& mat, std::size_t rows_begin, std::size_t cols_begin, std::size_t rows_end, std::size_t cols_end)
Tobias WEBER's avatar
Tobias WEBER committed
385
requires is_basic_mat<t_mat>
386
387
388
389
390
391
392
393
394
395
396
397
{
	for(std::size_t i=rows_begin; i<rows_end; ++i)
		for(std::size_t j=cols_begin; j<cols_end; ++j)
			mat(i,j) = (i==j ? 1 : 0);
}


/**
 * unit matrix
 */
template<class t_mat>
t_mat unit(std::size_t N1, std::size_t N2)
Tobias WEBER's avatar
Tobias WEBER committed
398
requires is_basic_mat<t_mat>
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
{
	t_mat mat;
	if constexpr(is_dyn_mat<t_mat>)
		mat = t_mat(N1, N2);

	unit<t_mat>(mat, 0,0, mat.size1(),mat.size2());
	return mat;
}


/**
 * unit matrix
 */
template<class t_mat>
t_mat unit(std::size_t N=0)
Tobias WEBER's avatar
Tobias WEBER committed
414
requires is_basic_mat<t_mat>
415
416
417
418
419
420
421
422
423
424
{
	return unit<t_mat>(N,N);
}


/**
 * zero matrix
 */
template<class t_mat>
t_mat zero(std::size_t N1, std::size_t N2)
Tobias WEBER's avatar
Tobias WEBER committed
425
requires is_basic_mat<t_mat>
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
{
	t_mat mat;
	if constexpr(is_dyn_mat<t_mat>)
		mat = t_mat(N1, N2);

	for(std::size_t i=0; i<mat.size1(); ++i)
		for(std::size_t j=0; j<mat.size2(); ++j)
			mat(i,j) = 0;

	return mat;
}

/**
 * zero matrix
 */
template<class t_mat>
t_mat zero(std::size_t N=0)
Tobias WEBER's avatar
Tobias WEBER committed
443
requires is_basic_mat<t_mat>
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
{
	return zero<t_mat>(N, N);
}


/**
 * zero vector
 */
template<class t_vec>
t_vec zero(std::size_t N=0)
requires is_basic_vec<t_vec>
{
	t_vec vec;
	if constexpr(is_dyn_vec<t_vec>)
		vec = t_vec(N);

	for(std::size_t i=0; i<vec.size(); ++i)
		vec[i] = 0;

	return vec;
}


Tobias WEBER's avatar
Tobias WEBER committed
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
/**
 * diagonal matrix
 */
template<class t_mat, class t_vec=std::vector<typename t_mat::value_type>>
t_mat diag(const t_vec& vals)
requires is_basic_mat<t_mat> && is_basic_vec<t_vec>
{
	const std::size_t N = vals.size();
	t_mat mat = zero<t_mat>(N);

	// static matrix does not necessarily have the required size!
	if constexpr(!m::is_dyn_mat<t_mat>)
		assert(mat.size1() == mat.size1() && mat.size1() == N);

	for(std::size_t i=0; i<std::min(mat.size1(), N); ++i)
		mat(i,i) = vals[i];

	return mat;
}


488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
/**
 * tests for zero vector
 */
template<class t_vec>
bool equals_0(const t_vec& vec,
	typename t_vec::value_type eps = std::numeric_limits<typename t_vec::value_type>::epsilon())
requires is_basic_vec<t_vec>
{
	return equals<t_vec>(vec, zero<t_vec>(vec.size()), eps);
}

/**
 * tests for zero matrix
 */
template<class t_mat>
bool equals_0(const t_mat& mat,
	typename t_mat::value_type eps = std::numeric_limits<typename t_mat::value_type>::epsilon())
requires is_mat<t_mat>
{
	return equals<t_mat>(mat, zero<t_mat>(mat.size1(), mat.size2()), eps);
}


/**
 * transpose matrix
 * WARNING: not possible for static non-square matrix!
 */
template<class t_mat>
t_mat trans(const t_mat& mat)
requires is_mat<t_mat>
{
	t_mat mat2;
	if constexpr(is_dyn_mat<t_mat>)
		mat2 = t_mat(mat.size2(), mat.size1());

	for(std::size_t i=0; i<mat.size1(); ++i)
		for(std::size_t j=0; j<mat.size2(); ++j)
			mat2(j,i) = mat(i,j);

	return mat2;
}


/**
 * create vector from initializer_list
 */
template<class t_vec>
t_vec create(const std::initializer_list<typename t_vec::value_type>& lst)
requires is_basic_vec<t_vec>
{
	t_vec vec;
	if constexpr(is_dyn_vec<t_vec>)
		vec = t_vec(lst.size());

	auto iterLst = lst.begin();
	for(std::size_t i=0; i<vec.size(); ++i)
	{
		if(iterLst != lst.end())
		{
			vec[i] = *iterLst;
			std::advance(iterLst, 1);
		}
		else	// vector larger than given list?
		{
			vec[i] = 0;
		}
	}

	return vec;
}


/**
 * create matrix from nested initializer_lists in columns/rows order
 */
template<class t_mat,
	template<class...> class t_cont_outer = std::initializer_list,
	template<class...> class t_cont = std::initializer_list>
t_mat create(const t_cont_outer<t_cont<typename t_mat::value_type>>& lst)
requires is_mat<t_mat>
{
	const std::size_t iCols = lst.size();
	const std::size_t iRows = lst.begin()->size();

	t_mat mat = unit<t_mat>(iRows, iCols);

	auto iterCol = lst.begin();
	for(std::size_t iCol=0; iCol<iCols; ++iCol)
	{
		auto iterRow = iterCol->begin();
		for(std::size_t iRow=0; iRow<iRows; ++iRow)
		{
			mat(iRow, iCol) = *iterRow;
			std::advance(iterRow, 1);
		}

		std::advance(iterCol, 1);
	}

	return mat;
}


/**
 * create matrix from column (or row) vectors
 */
template<class t_mat, class t_vec, template<class...> class t_cont_outer = std::initializer_list>
t_mat create(const t_cont_outer<t_vec>& lst, bool bRow = false)
requires is_mat<t_mat> && is_basic_vec<t_vec>
{
	const std::size_t iCols = lst.size();
	const std::size_t iRows = lst.begin()->size();

	t_mat mat = unit<t_mat>(iRows, iCols);

	auto iterCol = lst.begin();
	for(std::size_t iCol=0; iCol<iCols; ++iCol)
	{
		for(std::size_t iRow=0; iRow<iRows; ++iRow)
			mat(iRow, iCol) = (*iterCol)[iRow];
		std::advance(iterCol, 1);
	}

	if(bRow) mat = trans<t_mat>(mat);
	return mat;
}


/**
 * create matrix from initializer_list in column/row order
 */
template<class t_mat>
t_mat create(const std::initializer_list<typename t_mat::value_type>& lst)
requires is_mat<t_mat>
{
	const std::size_t N = std::sqrt(lst.size());

	t_mat mat = unit<t_mat>(N, N);

	auto iter = lst.begin();
	for(std::size_t iRow=0; iRow<N; ++iRow)
	{
		for(std::size_t iCol=0; iCol<N; ++iCol)
		{
			mat(iRow, iCol) = *iter;
			std::advance(iter, 1);
		}
	}

	return mat;
}


Tobias WEBER's avatar
Tobias WEBER committed
641
642
643
/**
 * convert between vector types
 */
Tobias WEBER's avatar
Tobias WEBER committed
644
template<class t_vecTo, class t_vecFrom>
Tobias WEBER's avatar
Tobias WEBER committed
645
646
647
648
649
650
651
652
653
654
t_vecTo convert(const t_vecFrom& vec)
requires is_basic_vec<t_vecFrom> && is_basic_vec<t_vecTo>
{
	using t_ty = typename t_vecTo::value_type;

	t_vecTo vecRet;
	if constexpr(is_dyn_vec<t_vecTo>)
		vecRet = t_vecTo(vec.size());

	for(std::size_t i=0; i<vec.size(); ++i)
Tobias WEBER's avatar
Tobias WEBER committed
655
		vecRet[i] = static_cast<t_ty>(vec[i]);
Tobias WEBER's avatar
Tobias WEBER committed
656
657
658
659
660

	return vecRet;
}


661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
/**
 * get a column vector from a matrix
 */
template<class t_mat, class t_vec>
t_vec col(const t_mat& mat, std::size_t col)
requires is_mat<t_mat> && is_basic_vec<t_vec>
{
	t_vec vec;
	if constexpr(is_dyn_vec<t_vec>)
		vec = t_vec(mat.size1());

	for(std::size_t i=0; i<mat.size1(); ++i)
		vec[i] = mat(i, col);

	return vec;
}

/**
 * get a row vector from a matrix
 */
template<class t_mat, class t_vec>
t_vec row(const t_mat& mat, std::size_t row)
requires is_mat<t_mat> && is_basic_vec<t_vec>
{
	t_vec vec;
	if constexpr(is_dyn_vec<t_vec>)
		vec = t_vec(mat.size2());

	for(std::size_t i=0; i<mat.size2(); ++i)
		vec[i] = mat(row, i);

	return vec;
}


/**
 * inner product <vec1|vec2>
 */
template<class t_vec>
typename t_vec::value_type inner(const t_vec& vec1, const t_vec& vec2)
requires is_basic_vec<t_vec>
{
	typename t_vec::value_type val(0);

	for(std::size_t i=0; i<vec1.size(); ++i)
	{
		if constexpr(is_complex<typename t_vec::value_type>)
			val += std::conj(vec1[i]) * vec2[i];
		else
			val += vec1[i] * vec2[i];
	}

	return val;
}


/**
 * inner product between two vectors of different type
 */
template<class t_vec1, class t_vec2>
typename t_vec1::value_type inner(const t_vec1& vec1, const t_vec2& vec2)
requires is_basic_vec<t_vec1> && is_basic_vec<t_vec2>
{
	if(vec1.size()==0 || vec2.size()==0)
		return typename t_vec1::value_type{};

	// first element
	auto val = vec1[0]*vec2[0];

	// remaining elements
	for(std::size_t i=1; i<std::min(vec1.size(), vec2.size()); ++i)
	{
		if constexpr(is_complex<typename t_vec1::value_type>)
		{
			auto prod = std::conj(vec1[i]) * vec2[i];
			val = val + prod;
		}
		else
		{
			auto prod = vec1[i]*vec2[i];
			val = val + prod;
		}
	}

	return val;
}


Tobias WEBER's avatar
Tobias WEBER committed
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
/**
 * matrix-matrix product
 */
template<class t_mat>
t_mat prod(const t_mat& mat1, const t_mat& mat2, bool assert_sizes=true)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
{
	// if not asserting sizes, the inner size will use the minimum of the two matrix sizes
	if(assert_sizes)
	{
		if constexpr(m::is_dyn_mat<t_mat>)
			assert((mat1.size2() == mat2.size1()));
		else
			static_assert(mat1.size2() == mat2.size1());
	}


	t_mat matRet(mat1.size1(), mat2.size2());
	const std::size_t innersize = std::min(mat1.size2(), mat2.size1());

	for(std::size_t row=0; row<matRet.size1(); ++row)
	{
		for(std::size_t col=0; col<matRet.size2(); ++col)
		{
			matRet(row, col) = 0;
			for(std::size_t i=0; i<innersize; ++i)
				matRet(row, col) += mat1(row, i) * mat2(i, col);
		}
	}

	return matRet;
}



784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
/**
 * 2-norm
 */
template<class t_vec>
typename t_vec::value_type norm(const t_vec& vec)
requires is_basic_vec<t_vec>
{
	return std::sqrt(inner<t_vec>(vec, vec));
}


/**
 * outer product
 */
template<class t_mat, class t_vec>
t_mat outer(const t_vec& vec1, const t_vec& vec2)
requires is_basic_vec<t_vec> && is_mat<t_mat>
{
	const std::size_t N1 = vec1.size();
	const std::size_t N2 = vec2.size();

	t_mat mat;
	if constexpr(is_dyn_mat<t_mat>)
		mat = t_mat(N1, N2);

	for(std::size_t n1=0; n1<N1; ++n1)
		for(std::size_t n2=0; n2<N2; ++n2)
		{
			if constexpr(is_complex<typename t_vec::value_type>)
				mat(n1, n2) = std::conj(vec1[n1]) * vec2[n2];
			else
				mat(n1, n2) = vec1[n1]*vec2[n2];
		}

	return mat;
}
Tobias WEBER's avatar
Tobias WEBER committed
820
}
821
822


Tobias WEBER's avatar
Tobias WEBER committed
823
824
825
namespace m_ops {
// ----------------------------------------------------------------------------
// vector operators
826
827
828
// ----------------------------------------------------------------------------

/**
Tobias WEBER's avatar
Tobias WEBER committed
829
 * unary +
830
 */
Tobias WEBER's avatar
Tobias WEBER committed
831
832
833
template<class t_vec>
const t_vec& operator+(const t_vec& vec1)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
834
{
Tobias WEBER's avatar
Tobias WEBER committed
835
	return vec1;
836
837
838
839
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
840
 * unary -
841
 */
Tobias WEBER's avatar
Tobias WEBER committed
842
843
844
template<class t_vec>
t_vec operator-(const t_vec& vec1)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
845
{
Tobias WEBER's avatar
Tobias WEBER committed
846
	t_vec vec(vec1.size());
847

Tobias WEBER's avatar
Tobias WEBER committed
848
849
	for(std::size_t i=0; i<vec1.size(); ++i)
		vec[i] = -vec1[i];
850

Tobias WEBER's avatar
Tobias WEBER committed
851
	return vec;
852
853
854
855
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
856
 * binary +
857
 */
Tobias WEBER's avatar
Tobias WEBER committed
858
859
860
template<class t_vec>
t_vec operator+(const t_vec& vec1, const t_vec& vec2)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
861
{
Tobias WEBER's avatar
Tobias WEBER committed
862
863
864
865
	if constexpr(m::is_dyn_vec<t_vec>)
		assert((vec1.size() == vec2.size()));
	else
		static_assert(vec1.size() == vec2.size());
866

Tobias WEBER's avatar
Tobias WEBER committed
867
	t_vec vec(vec1.size());
868

Tobias WEBER's avatar
Tobias WEBER committed
869
870
871
872
	for(std::size_t i=0; i<vec1.size(); ++i)
		vec[i] = vec1[i] + vec2[i];

	return vec;
873
874
875
876
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
877
 * binary -
878
 */
Tobias WEBER's avatar
Tobias WEBER committed
879
880
881
template<class t_vec>
t_vec operator-(const t_vec& vec1, const t_vec& vec2)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
882
{
Tobias WEBER's avatar
Tobias WEBER committed
883
	return vec1 + (-vec2);
884
885
886
887
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
888
 * vector * scalar
889
 */
Tobias WEBER's avatar
Tobias WEBER committed
890
891
892
template<class t_vec>
t_vec operator*(const t_vec& vec1, typename t_vec::value_type d)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
893
{
Tobias WEBER's avatar
Tobias WEBER committed
894
895
896
897
	t_vec vec(vec1.size());

	for(std::size_t i=0; i<vec1.size(); ++i)
		vec[i] = vec1[i] * d;
898

Tobias WEBER's avatar
Tobias WEBER committed
899
900
	return vec;
}
901
902
903


/**
Tobias WEBER's avatar
Tobias WEBER committed
904
 * scalar * vector
905
 */
Tobias WEBER's avatar
Tobias WEBER committed
906
907
908
909
template<class t_vec>
t_vec operator*(typename t_vec::value_type d, const t_vec& vec)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
	//&& !m::is_basic_mat<typename t_vec::value_type>	// hack!
910
{
Tobias WEBER's avatar
Tobias WEBER committed
911
	return vec * d;
912
913
914
}

/**
Tobias WEBER's avatar
Tobias WEBER committed
915
 * vector / scalar
916
917
 */
template<class t_vec>
Tobias WEBER's avatar
Tobias WEBER committed
918
919
t_vec operator/(const t_vec& vec, typename t_vec::value_type d)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
920
{
Tobias WEBER's avatar
Tobias WEBER committed
921
922
	using T = typename t_vec::value_type;
	return vec * (T(1)/d);
923
924
925
926
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
927
 * vector += vector
928
929
 */
template<class t_vec>
Tobias WEBER's avatar
Tobias WEBER committed
930
931
t_vec& operator+=(t_vec& vec1, const t_vec& vec2)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
932
{
Tobias WEBER's avatar
Tobias WEBER committed
933
934
	vec1 = vec1 + vec2;
	return vec1;
935
936
937
}

/**
Tobias WEBER's avatar
Tobias WEBER committed
938
 * vector -= vector
939
940
 */
template<class t_vec>
Tobias WEBER's avatar
Tobias WEBER committed
941
942
t_vec& operator-=(t_vec& vec1, const t_vec& vec2)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
943
{
Tobias WEBER's avatar
Tobias WEBER committed
944
945
	vec1 = vec1 - vec2;
	return vec1;
946
947
948
949
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
950
 * vector *= scalar
951
 */
Tobias WEBER's avatar
Tobias WEBER committed
952
953
954
template<class t_vec>
t_vec& operator*=(t_vec& vec1, typename t_vec::value_type d)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
955
{
Tobias WEBER's avatar
Tobias WEBER committed
956
957
	vec1 = vec1 * d;
	return vec1;
958
959
960
}

/**
Tobias WEBER's avatar
Tobias WEBER committed
961
 * vector /= scalar
962
 */
Tobias WEBER's avatar
Tobias WEBER committed
963
964
965
template<class t_vec>
t_vec& operator/=(t_vec& vec1, typename t_vec::value_type d)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
966
{
Tobias WEBER's avatar
Tobias WEBER committed
967
968
	vec1 = vec1 / d;
	return vec1;
969
970
971
}


Tobias WEBER's avatar
Tobias WEBER committed
972

973
/**
Tobias WEBER's avatar
Tobias WEBER committed
974
 * operator <<
975
 */
Tobias WEBER's avatar
Tobias WEBER committed
976
977
978
template<class t_vec>
std::ostream& operator<<(std::ostream& ostr, const t_vec& vec)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
979
980
981
{
	const std::size_t N = vec.size();

Tobias WEBER's avatar
Tobias WEBER committed
982
983
984
985
986
987
	for(std::size_t i=0; i<N; ++i)
	{
		ostr << vec[i];
		if(i < N-1)
			ostr << COLSEP << " ";
	}
988

Tobias WEBER's avatar
Tobias WEBER committed
989
	return ostr;
990
991
992
993
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
994
 * operator >>
995
 */
Tobias WEBER's avatar
Tobias WEBER committed
996
997
998
template<class t_vec>
std::istream& operator>>(std::istream& istr, t_vec& vec)
requires m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
999
{
Tobias WEBER's avatar
Tobias WEBER committed
1000
	vec.clear();
1001

Tobias WEBER's avatar
Tobias WEBER committed
1002
1003
	std::string str;
	std::getline(istr, str);
1004

Tobias WEBER's avatar
Tobias WEBER committed
1005
1006
1007
1008
	std::vector<std::string> vecstr;
	boost::split(vecstr, str, [](auto c)->bool { return c==COLSEP; }, boost::token_compress_on);

	for(auto& tok : vecstr)
1009
	{
Tobias WEBER's avatar
Tobias WEBER committed
1010
1011
1012
1013
1014
1015
		boost::trim(tok);
		std::istringstream istr(tok);

		typename t_vec::value_type c{};
		istr >> c;
		vec.emplace_back(std::move(c));
1016
1017
	}

Tobias WEBER's avatar
Tobias WEBER committed
1018
	return istr;
1019
}
Tobias WEBER's avatar
Tobias WEBER committed
1020
// ----------------------------------------------------------------------------
1021
1022
1023



Tobias WEBER's avatar
Tobias WEBER committed
1024
1025
1026
// ----------------------------------------------------------------------------
// matrix operators
// ----------------------------------------------------------------------------
1027
1028

/**
Tobias WEBER's avatar
Tobias WEBER committed
1029
 * unary +
1030
 */
Tobias WEBER's avatar
Tobias WEBER committed
1031
1032
1033
template<class t_mat>
const t_mat& operator+(const t_mat& mat1)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1034
{
Tobias WEBER's avatar
Tobias WEBER committed
1035
	return mat1;
1036
1037
1038
1039
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1040
 * unary -
1041
 */
Tobias WEBER's avatar
Tobias WEBER committed
1042
1043
1044
template<class t_mat>
t_mat operator-(const t_mat& mat1)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1045
{
Tobias WEBER's avatar
Tobias WEBER committed
1046
	t_mat mat(mat1.size1(), mat1.size2());
1047

Tobias WEBER's avatar
Tobias WEBER committed
1048
1049
1050
1051
1052
	for(std::size_t i=0; i<mat1.size1(); ++i)
		for(std::size_t j=0; j<mat1.size2(); ++j)
			mat(i,j) = -mat1(i,j);

	return mat;
1053
1054
1055
1056
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1057
 * binary +
1058
 */
Tobias WEBER's avatar
Tobias WEBER committed
1059
1060
1061
template<class t_mat>
t_mat operator+(const t_mat& mat1, const t_mat& mat2)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1062
{
Tobias WEBER's avatar
Tobias WEBER committed
1063
1064
1065
1066
	if constexpr(m::is_dyn_mat<t_mat>)
		assert((mat1.size1() == mat2.size1() && mat1.size2() == mat2.size2()));
	else
		static_assert(mat1.size1() == mat2.size1() && mat1.size2() == mat2.size2());
1067

Tobias WEBER's avatar
Tobias WEBER committed
1068
	t_mat mat(mat1.size1(), mat1.size2());
1069

Tobias WEBER's avatar
Tobias WEBER committed
1070
1071
1072
	for(std::size_t i=0; i<mat1.size1(); ++i)
		for(std::size_t j=0; j<mat1.size2(); ++j)
			mat(i,j) = mat1(i,j) + mat2(i,j);
1073

Tobias WEBER's avatar
Tobias WEBER committed
1074
	return mat;
1075
1076
1077
1078
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1079
 * binary -
1080
 */
Tobias WEBER's avatar
Tobias WEBER committed
1081
1082
1083
template<class t_mat>
t_mat operator-(const t_mat& mat1, const t_mat& mat2)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1084
{
Tobias WEBER's avatar
Tobias WEBER committed
1085
	return mat1 + (-mat2);
1086
1087
1088
1089
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1090
 * matrix * scalar
1091
 */
Tobias WEBER's avatar
Tobias WEBER committed
1092
1093
1094
template<class t_mat>
t_mat operator*(const t_mat& mat1, typename t_mat::value_type d)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1095
{
Tobias WEBER's avatar
Tobias WEBER committed
1096
	t_mat mat(mat1.size1(), mat1.size2());
1097

Tobias WEBER's avatar
Tobias WEBER committed
1098
1099
1100
	for(std::size_t i=0; i<mat1.size1(); ++i)
		for(std::size_t j=0; j<mat1.size2(); ++j)
			mat(i,j) = mat1(i,j) * d;
1101

Tobias WEBER's avatar
Tobias WEBER committed
1102
	return mat;
1103
1104
1105
}

/**
Tobias WEBER's avatar
Tobias WEBER committed
1106
 * scalar * matrix
1107
 */
Tobias WEBER's avatar
Tobias WEBER committed
1108
1109
1110
template<class t_mat>
t_mat operator*(typename t_mat::value_type d, const t_mat& mat)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1111
{
Tobias WEBER's avatar
Tobias WEBER committed
1112
	return mat * d;
1113
1114
1115
}

/**
Tobias WEBER's avatar
Tobias WEBER committed
1116
 * matrix / scalar
1117
1118
 */
template<class t_mat>
Tobias WEBER's avatar
Tobias WEBER committed
1119
1120
t_mat operator/(const t_mat& mat, typename t_mat::value_type d)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1121
1122
{
	using T = typename t_mat::value_type;
Tobias WEBER's avatar
Tobias WEBER committed
1123
	return mat * (T(1)/d);
1124
1125
1126
1127
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1128
 * matrix-matrix product
1129
1130
 */
template<class t_mat>
Tobias WEBER's avatar
Tobias WEBER committed
1131
1132
t_mat operator*(const t_mat& mat1, const t_mat& mat2)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1133
{
Tobias WEBER's avatar
Tobias WEBER committed
1134
	return m::prod<t_mat>(mat1, mat2);
1135
1136
1137
1138
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1139
 * matrix *= scalar
1140
1141
 */
template<class t_mat>
Tobias WEBER's avatar
Tobias WEBER committed
1142
1143
t_mat& operator*=(t_mat& mat1, typename t_mat::value_type d)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1144
{
Tobias WEBER's avatar
Tobias WEBER committed
1145
1146
1147
	mat1 = mat1 * d;
	return mat1;
}
1148

Tobias WEBER's avatar
Tobias WEBER committed
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
/**
 * matrix /= scalar
 */
template<class t_mat>
t_mat& operator/=(t_mat& mat1, typename t_mat::value_type d)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
{
	mat1 = mat1 / d;
	return mat1;
}
1159
1160


Tobias WEBER's avatar
Tobias WEBER committed
1161
1162
1163
1164
1165
1166
1167
1168
1169
/**
 * operator <<
 */
template<class t_mat>
std::ostream& operator<<(std::ostream& ostr, const t_mat& mat)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
{
	const std::size_t ROWS = mat.size1();
	const std::size_t COLS = mat.size2();
1170

Tobias WEBER's avatar
Tobias WEBER committed
1171
	for(std::size_t row=0; row<ROWS; ++row)
1172
	{
Tobias WEBER's avatar
Tobias WEBER committed
1173
		for(std::size_t col=0; col<COLS; ++col)
1174
		{
Tobias WEBER's avatar
Tobias WEBER committed
1175
1176
1177
			ostr << mat(row, col);
			if(col < COLS-1)
				ostr << COLSEP << " ";
1178
		}
Tobias WEBER's avatar
Tobias WEBER committed
1179
1180
1181

		if(row < ROWS-1)
			ostr << ROWSEP << " ";
1182
1183
	}

Tobias WEBER's avatar
Tobias WEBER committed
1184
	return ostr;
1185
1186
1187
1188
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1189
 * prints matrix in nicely formatted form
1190
 */
Tobias WEBER's avatar
Tobias WEBER committed
1191
1192
1193
template<class t_mat>
std::ostream& niceprint(std::ostream& ostr, const t_mat& mat)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
1194
{
Tobias WEBER's avatar
Tobias WEBER committed
1195
1196
	const std::size_t ROWS = mat.size1();
	const std::size_t COLS = mat.size2();
1197

Tobias WEBER's avatar
Tobias WEBER committed
1198
	for(std::size_t i=0; i<ROWS; ++i)
1199
	{
Tobias WEBER's avatar
Tobias WEBER committed
1200
1201
1202
1203
1204
1205
1206
		ostr << "(";
		for(std::size_t j=0; j<COLS; ++j)
			ostr << std::setw(ostr.precision()*1.5) << std::right << mat(i,j);
		ostr << ")";

		if(i < ROWS-1)
			ostr << "\n";
1207
1208
	}

Tobias WEBER's avatar
Tobias WEBER committed
1209
	return ostr;
1210
}
Tobias WEBER's avatar
Tobias WEBER committed
1211
1212
1213
// ----------------------------------------------------------------------------


1214

Tobias WEBER's avatar
Tobias WEBER committed
1215
1216
1217
// ----------------------------------------------------------------------------
// mixed operators
// ----------------------------------------------------------------------------
1218
1219

/**
Tobias WEBER's avatar
Tobias WEBER committed
1220
 * matrix-vector product
1221
 */
Tobias WEBER's avatar
Tobias WEBER committed
1222
1223
1224
1225
template<class t_mat, class t_vec>
t_vec operator*(const t_mat& mat, const t_vec& vec)
requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
	&& m::is_basic_vec<t_vec> && m::is_dyn_vec<t_vec>
1226
{
Tobias WEBER's avatar
Tobias WEBER committed
1227
1228
1229
1230
	if constexpr(m::is_dyn_mat<t_mat>)
		assert((mat.size2() == vec.size()));
	else
		static_assert(mat.size2() == vec.size());
1231
1232


Tobias WEBER's avatar
Tobias WEBER committed
1233
1234
1235
1236
1237
1238
	t_vec vecRet(mat.size1());

	for(std::size_t row=0; row<mat.size1(); ++row)
	{
		vecRet[row] = typename t_vec::value_type{/*0*/};
		for(std::size_t col=0; col<mat.size2(); ++col)
1239
		{
Tobias WEBER's avatar
Tobias WEBER committed
1240
1241
			auto elem = mat(row, col) * vec[col];
			vecRet[row] = vecRet[row] + elem;
1242
1243
1244
		}
	}

Tobias WEBER's avatar
Tobias WEBER committed
1245
1246
1247
	return vecRet;
}
// ----------------------------------------------------------------------------
1248
1249
1250
}


Tobias WEBER's avatar
Tobias WEBER committed
1251
1252
1253
1254
1255

namespace m{
// ----------------------------------------------------------------------------
// with metric

1256
/**
Tobias WEBER's avatar
Tobias WEBER committed
1257
1258
 * covariant metric tensor
 * g_{i,j} = e_i * e_j
1259
 */
Tobias WEBER's avatar
Tobias WEBER committed
1260
1261
1262
template<class t_mat, class t_vec, template<class...> class t_cont=std::initializer_list>
t_mat metric(const t_cont<t_vec>& basis_co)
requires is_basic_mat<t_mat> && is_basic_vec<t_vec>
1263
{
Tobias WEBER's avatar
Tobias WEBER committed
1264
	const std::size_t N = basis_co.size();
1265

Tobias WEBER's avatar
Tobias WEBER committed
1266
1267
1268
1269
1270
1271
	t_mat g_co;
	if constexpr(is_dyn_mat<t_mat>)
		g_co = t_mat(N, N);

	auto iter_i = basis_co.begin();
	for(std::size_t i=0; i<N; ++i)
1272
	{
Tobias WEBER's avatar
Tobias WEBER committed
1273
1274
1275
1276
1277
1278
1279
		auto iter_j = basis_co.begin();
		for(std::size_t j=0; j<N; ++j)
		{
			g_co(i,j) = inner<t_vec>(*iter_i, *iter_j);
			std::advance(iter_j, 1);
		}
		std::advance(iter_i, 1);
1280
1281
	}

Tobias WEBER's avatar
Tobias WEBER committed
1282
	return g_co;
1283
1284
1285
1286
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1287
 * lower index using metric
1288
 */
Tobias WEBER's avatar
Tobias WEBER committed
1289
1290
1291
template<class t_mat, class t_vec>
t_vec lower_index(const t_mat& metric_co, const t_vec& vec_contra)
requires is_basic_mat<t_mat> && is_basic_vec<t_vec>
1292
{
Tobias WEBER's avatar
Tobias WEBER committed
1293
1294
	const std::size_t N = vec_contra.size();
	t_vec vec_co = zero<t_vec>(N);
1295

Tobias WEBER's avatar
Tobias WEBER committed
1296
1297
1298
	for(std::size_t i=0; i<N; ++i)
		for(std::size_t j=0; j<N; ++j)
			vec_co[i] += metric_co(i,j) * vec_contra[j];
1299

Tobias WEBER's avatar
Tobias WEBER committed
1300
	return vec_co;
1301
1302
1303
1304
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1305
 * raise index using metric
1306
 */
Tobias WEBER's avatar
Tobias WEBER committed
1307
1308
1309
template<class t_mat, class t_vec>
t_vec raise_index(const t_mat& metric_contra, const t_vec& vec_co)
requires is_basic_mat<t_mat> && is_basic_vec<t_vec>
1310
{
Tobias WEBER's avatar
Tobias WEBER committed
1311
1312
	const std::size_t N = vec_co.size();
	t_vec vec_contra = zero<t_vec>(N);
1313

Tobias WEBER's avatar
Tobias WEBER committed
1314
1315
1316
	for(std::size_t i=0; i<N; ++i)
		for(std::size_t j=0; j<N; ++j)
			vec_contra[i] += metric_contra(i,j) * vec_co[j];
1317

Tobias WEBER's avatar
Tobias WEBER committed
1318
1319
	return vec_contra;
}
1320
1321


Tobias WEBER's avatar
Tobias WEBER committed
1322
1323
1324
1325
1326
1327
1328
1329
1330
/**
 * inner product using metric
 */
template<class t_mat, class t_vec>
typename t_vec::value_type inner(const t_mat& metric_co, const t_vec& vec1_contra, const t_vec& vec2_contra)
requires is_basic_mat<t_mat> && is_basic_vec<t_vec>
{
	t_vec vec2_co = lower_index<t_mat, t_vec>(metric_co, vec2_contra);
	return inner<t_vec>(vec1_contra, vec2_co);
1331
1332
1333
1334
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1335
 * 2-norm using metric
1336
 */
Tobias WEBER's avatar
Tobias WEBER committed
1337
1338
1339
template<class t_mat, class t_vec>
typename t_vec::value_type norm(const t_mat& metric_co, const t_vec& vec_contra)
requires is_basic_vec<t_vec>
1340
{
Tobias WEBER's avatar
Tobias WEBER committed
1341
1342
1343
	return std::sqrt(inner<t_mat, t_vec>(metric_co, vec_contra, vec_contra));
}
// ----------------------------------------------------------------------------
1344
1345
1346



Tobias WEBER's avatar
Tobias WEBER committed
1347
1348
1349
1350
1351
1352
1353
1354
1355
/**
 * matrix to project onto vector: P = |v><v|
 * from: |x'> = <v|x> * |v> = |v><v|x> = |v><v| * |x>
 */
template<class t_mat, class t_vec>
t_mat projector(const t_vec& vec, bool bIsNormalised = true)
requires is_vec<t_vec> && is_mat<t_mat>
{
	if(bIsNormalised)
1356
	{
Tobias WEBER's avatar
Tobias WEBER committed
1357
1358
1359
1360
1361
1362
1363
		return outer<t_mat, t_vec>(vec, vec);
	}
	else
	{
		const auto len = norm<t_vec>(vec);
		t_vec _vec = vec / len;
		return outer<t_mat, t_vec>(_vec, _vec);
1364
1365
1366
1367
1368
	}
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1369
 * project vector vec onto another vector vecProj
1370
 */
Tobias WEBER's avatar
Tobias WEBER committed
1371
1372
1373
template<class t_vec>
t_vec project(const t_vec& vec, const t_vec& vecProj, bool bIsNormalised = true)
requires is_vec<t_vec>
1374
{
Tobias WEBER's avatar
Tobias WEBER committed
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
	if(bIsNormalised)
	{
		return inner<t_vec>(vec, vecProj) * vecProj;
	}
	else
	{
		const auto len = norm<t_vec>(vecProj);
		const t_vec _vecProj = vecProj / len;
		return inner<t_vec>(vec, _vecProj) * _vecProj;
	}
1385
1386
1387
1388
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1389
1390
 * project vector vec onto another vector vecProj
 * don't multiply with direction vector
1391
1392
 */
template<class t_vec>
Tobias WEBER's avatar
Tobias WEBER committed
1393
1394
typename t_vec::value_type
project_scalar(const t_vec& vec, const t_vec& vecProj, bool bIsNormalised = true)
1395
1396
requires is_vec<t_vec>
{
Tobias WEBER's avatar
Tobias WEBER committed
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
	if(bIsNormalised)
	{
		return inner<t_vec>(vec, vecProj);
	}
	else
	{
		const auto len = norm<t_vec>(vecProj);
		const t_vec _vecProj = vecProj / len;
		return inner<t_vec>(vec, _vecProj);
	}
}
1408
1409


Tobias WEBER's avatar
Tobias WEBER committed
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
/**
 * project vector vec onto the line lineOrigin + lam*lineDir
 * shift line to go through origin, calculate projection and shift back
 * returns [closest point, distance]
 */
template<class t_vec>
std::tuple<t_vec, typename t_vec::value_type> project_line(const t_vec& vec,
	const t_vec& lineOrigin, const t_vec& lineDir, bool bIsNormalised = true)
requires is_vec<t_vec>
{
	const t_vec ptShifted = vec - lineOrigin;
	const t_vec ptProj = project<t_vec>(ptShifted, lineDir, bIsNormalised);
	const t_vec ptNearest = lineOrigin + ptProj;
1423

Tobias WEBER's avatar
Tobias WEBER committed
1424
1425
1426
	const typename t_vec::value_type dist = norm<t_vec>(vec - ptNearest);
	return std::make_tuple(ptNearest, dist);
}
1427
1428


Tobias WEBER's avatar
Tobias WEBER committed
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
/**
 * matrix to project onto orthogonal complement (plane perpendicular to vector): P = 1-|v><v|
 * from completeness relation: 1 = sum_i |v_i><v_i| = |x><x| + |y><y| + |z><z|
 */
template<class t_mat, class t_vec>
t_mat ortho_projector(const t_vec& vec, bool bIsNormalised = true)
requires is_vec<t_vec> && is_mat<t_mat>
{
	const std::size_t iSize = vec.size();
	return unit<t_mat>(iSize) -
		projector<t_mat, t_vec>(vec, bIsNormalised);
}
1441
1442


Tobias WEBER's avatar
Tobias WEBER committed
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
/**
 * matrix to mirror on plane perpendicular to vector: P = 1 - 2*|v><v|
 * subtracts twice its projection onto the plane normal from the vector
 */
template<class t_mat, class t_vec>
t_mat ortho_mirror_op(const t_vec& vec, bool bIsNormalised = true)
requires is_vec<t_vec> && is_mat<t_mat>
{
	using T = typename t_vec::value_type;
	const std::size_t iSize = vec.size();
1453

Tobias WEBER's avatar
Tobias WEBER committed
1454
1455
	return unit<t_mat>(iSize) -
		T(2)*projector<t_mat, t_vec>(vec, bIsNormalised);
1456
1457
1458
1459
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
1460
 * matrix to mirror [a, b, c, ...] into, e.g.,  [a, b', 0, 0]
1461
 */
Tobias WEBER's avatar
Tobias WEBER committed
1462
1463
1464
template<class t_mat, class t_vec>
t_mat ortho_mirror_zero_op(const t_vec& vec, std::size_t row)
requires is_vec<t_vec> && is_mat<t_mat>
1465
1466
{
	using T = typename t_vec::value_type;
Tobias WEBER's avatar
Tobias WEBER committed
1467
	const std::size_t N = vec.size();
1468

Tobias WEBER's avatar
Tobias WEBER committed
1469
1470
1471
	t_vec vecSub = zero<t_vec>(N);
	for(std::size_t i=0; i<row; ++i)
		vecSub[i] = vec[i];
1472

Tobias WEBER's avatar
Tobias WEBER committed
1473
1474
1475
1476
1477
	// norm of rest vector
	T n = T(0);
	for(std::size_t i=row; i<N; ++i)
		n += vec[i]*vec[i];
	vecSub[row] = std::sqrt(n);