SimiLie
Loading...
Searching...
No Matches
young_tableau.hpp
1// SPDX-FileCopyrightText: 2024 Baptiste Legouix
2// SPDX-License-Identifier: MIT
3
4#pragma once
5
6#if defined BUILD_YOUNG_TABLEAU
7
8#include <ddc/ddc.hpp>
9
10#include <similie/csr/csr.hpp>
11#include <similie/csr/csr_dynamic.hpp>
12#include <similie/misc/macros.hpp>
13#include <similie/misc/permutation_parity.hpp>
14#include <similie/misc/specialization.hpp>
15#include <similie/misc/stride.hpp>
16#include <similie/tensor/dummy_index.hpp>
17#include <similie/tensor/prime.hpp>
18#include <similie/tensor/tensor_impl.hpp>
19
20namespace sil {
21
22namespace young_tableau {
23
24namespace detail {
25
26template <class T, T... Args>
27constexpr T sum(std::integer_sequence<T, Args...> = {})
28{
29 return (Args + ...);
30}
31
32} // namespace detail
33
34template <class... Row>
35class YoungTableauSeq
36{
37public:
38 using shape = std::index_sequence<Row::size()...>;
39
40 static constexpr std::size_t rank = detail::sum(shape());
41};
42
43namespace detail {
44
45// Extract Row from a YoungTableauSeq at index I
46template <std::size_t I, std::size_t Id, class TableauSeq>
47struct ExtractRow;
48
49template <std::size_t I, std::size_t Id, class HeadRow, class... TailRow>
50struct ExtractRow<I, Id, YoungTableauSeq<HeadRow, TailRow...>>
51{
52 using type = std::conditional_t<
53 I == Id,
54 HeadRow,
55 typename ExtractRow<I, Id + 1, YoungTableauSeq<TailRow...>>::type>;
56};
57
58template <std::size_t I, std::size_t Id>
59struct ExtractRow<I, Id, YoungTableauSeq<>>
60{
61 using type = std::index_sequence<>;
62};
63
64template <std::size_t I, class TableauSeq>
65using extract_row_t = ExtractRow<I, 0, TableauSeq>::type;
66
67// Override the row of YoungTableauSeq at index I with RowToSet
68template <
69 std::size_t I,
70 std::size_t Id,
71 class RowToSet,
72 class HeadTableauSeq,
73 class InterestTableauSeq,
74 class TailTableauSeq>
75struct OverrideRow;
76
77template <
78 std::size_t I,
79 std::size_t Id,
80 class RowToSet,
81 class... HeadRow,
82 class InterestRow,
83 class HeadOfTailRow,
84 class... TailOfTailRow>
85struct OverrideRow<
86 I,
87 Id,
88 RowToSet,
89 YoungTableauSeq<HeadRow...>,
90 YoungTableauSeq<InterestRow>,
91 YoungTableauSeq<HeadOfTailRow, TailOfTailRow...>>
92{
93 using type = std::conditional_t<
94 I == Id,
95 YoungTableauSeq<HeadRow..., RowToSet, HeadOfTailRow, TailOfTailRow...>,
96 typename OverrideRow<
97 I,
98 Id + 1,
99 RowToSet,
100 YoungTableauSeq<HeadRow..., InterestRow>,
101 YoungTableauSeq<HeadOfTailRow>,
102 YoungTableauSeq<TailOfTailRow...>>::type>;
103};
104
105template <std::size_t I, std::size_t Id, class RowToSet, class... HeadRow, class InterestRow>
106struct OverrideRow<
107 I,
108 Id,
109 RowToSet,
110 YoungTableauSeq<HeadRow...>,
111 YoungTableauSeq<InterestRow>,
112 YoungTableauSeq<>>
113{
114 using type = YoungTableauSeq<HeadRow..., RowToSet>;
115};
116
117template <std::size_t I, class RowToSet, class TableauSeq>
118struct OverrideRowHelper;
119
120template <std::size_t I, class RowToSet, class InterestRow, class... TailRow>
121struct OverrideRowHelper<I, RowToSet, YoungTableauSeq<InterestRow, TailRow...>>
122{
123 using type = OverrideRow<
124 I,
125 0,
126 RowToSet,
127 YoungTableauSeq<>,
128 YoungTableauSeq<InterestRow>,
129 YoungTableauSeq<TailRow...>>::type;
130};
131
132template <std::size_t I, class RowToSet, class TableauSeq>
133using override_row_t = OverrideRowHelper<I, RowToSet, TableauSeq>::type;
134
135// Add cell with value Value at the end of a row
136template <class Row, std::size_t Value>
137struct AddCellToRow;
138
139template <std::size_t... Elem, std::size_t Value>
140struct AddCellToRow<std::index_sequence<Elem...>, Value>
141{
142 using type = std::index_sequence<Elem..., Value>;
143};
144
145template <class Row, std::size_t Value>
146using add_cell_to_row_t = AddCellToRow<Row, Value>::type;
147
148// Add cell with value Value at the end of the row at index I
149template <class Tableau, std::size_t I, std::size_t Value>
150struct AddCellToTableau;
151
152template <class... Row, std::size_t I, std::size_t Value>
153struct AddCellToTableau<YoungTableauSeq<Row...>, I, Value>
154{
155 using type = detail::override_row_t<
156 I,
157 add_cell_to_row_t<extract_row_t<I, YoungTableauSeq<Row...>>, Value>,
158 YoungTableauSeq<Row...>>;
159};
160
161template <class Tableau, std::size_t I, std::size_t Value>
162using add_cell_to_tableau_t = AddCellToTableau<Tableau, I, Value>::type;
163
164// Compute dual of YoungTableauSeq
165template <class TableauDual, class Tableau, std::size_t I, std::size_t J>
166struct Dual;
167
168template <
169 class TableauDual,
170 std::size_t HeadElemOfHeadRow,
171 std::size_t... TailElemOfHeadRow,
172 class... TailRow,
173 std::size_t I,
174 std::size_t J>
175struct Dual<
176 TableauDual,
177 YoungTableauSeq<std::index_sequence<HeadElemOfHeadRow, TailElemOfHeadRow...>, TailRow...>,
178 I,
179 J>
180{
181 using type = std::conditional_t<
182 sizeof...(TailRow) == 0 && sizeof...(TailElemOfHeadRow) == 0,
183 add_cell_to_tableau_t<TableauDual, J, HeadElemOfHeadRow>,
184 std::conditional_t<
185 sizeof...(TailElemOfHeadRow) == 0,
186 typename Dual<
187 add_cell_to_tableau_t<TableauDual, J, HeadElemOfHeadRow>,
188 YoungTableauSeq<TailRow...>,
189 I + 1,
190 0>::type,
191 typename Dual<
192 add_cell_to_tableau_t<TableauDual, J, HeadElemOfHeadRow>,
193 YoungTableauSeq<std::index_sequence<TailElemOfHeadRow...>, TailRow...>,
194 I,
195 J + 1>::type>>;
196};
197
198template <class TableauDual, class... TailRow, std::size_t I, std::size_t J>
199struct Dual<TableauDual, YoungTableauSeq<std::index_sequence<>, TailRow...>, I, J>
200{
201 using type = YoungTableauSeq<>;
202};
203
204template <class TableauDual, std::size_t I, std::size_t J>
205struct Dual<TableauDual, YoungTableauSeq<>, I, J>
206{
207 using type = TableauDual;
208};
209
210template <class Tableau>
211struct DualHelper;
212
213template <std::size_t... ElemOfHeadRow, class... TailRow>
214struct DualHelper<YoungTableauSeq<std::index_sequence<ElemOfHeadRow...>, TailRow...>>
215{
216 using type
217 = Dual<YoungTableauSeq<std::conditional_t<
218 ElemOfHeadRow == -1,
219 std::index_sequence<>,
220 std::index_sequence<>>...>,
221 YoungTableauSeq<std::index_sequence<ElemOfHeadRow...>, TailRow...>,
222 0,
223 0>::type;
224};
225
226template <>
227struct DualHelper<YoungTableauSeq<>>
228{
229 using type = YoungTableauSeq<>;
230};
231
232template <class Tableau>
233using dual_t = DualHelper<Tableau>::type;
234
235// Compute hooks
236template <class TableauHooks, class Tableau, class Shape, std::size_t I, std::size_t J>
237struct Hooks;
238
239template <
240 class TableauHooks,
241 std::size_t HeadElemOfHeadRow,
242 std::size_t... TailElemOfHeadRow,
243 class... TailRow,
244 std::size_t HeadRowSize,
245 std::size_t... TailRowSize,
246 std::size_t I,
247 std::size_t J>
248struct Hooks<
249 TableauHooks,
250 YoungTableauSeq<std::index_sequence<HeadElemOfHeadRow, TailElemOfHeadRow...>, TailRow...>,
251 std::index_sequence<HeadRowSize, TailRowSize...>,
252 I,
253 J>
254{
255 using type = std::conditional_t<
256 sizeof...(TailRow) == 0 && sizeof...(TailElemOfHeadRow) == 0,
257 add_cell_to_tableau_t<TableauHooks, I, HeadRowSize - J>,
258 std::conditional_t<
259 sizeof...(TailElemOfHeadRow) == 0,
260 typename Hooks<
261 add_cell_to_tableau_t<TableauHooks, I, HeadRowSize - J>,
262 YoungTableauSeq<TailRow...>,
263 std::index_sequence<TailRowSize...>,
264 I + 1,
265 0>::type,
266 typename Hooks<
267 add_cell_to_tableau_t<TableauHooks, I, HeadRowSize - J>,
268 YoungTableauSeq<std::index_sequence<TailElemOfHeadRow...>, TailRow...>,
269 std::index_sequence<HeadRowSize, TailRowSize...>,
270 I,
271 J + 1>::type>>;
272};
273
274template <class TableauHooks, class... TailRow, class Shape, std::size_t I, std::size_t J>
275struct Hooks<TableauHooks, YoungTableauSeq<std::index_sequence<>, TailRow...>, Shape, I, J>
276{
277 using type = YoungTableauSeq<>;
278};
279
280template <class TableauHooks, class Shape, std::size_t I, std::size_t J>
281struct Hooks<TableauHooks, YoungTableauSeq<>, Shape, I, J>
282{
283 using type = TableauHooks;
284};
285
286template <class Tableau>
287struct HooksHelper;
288
289template <class... Row>
290struct HooksHelper<YoungTableauSeq<Row...>>
291{
292 using type
293 = Hooks<YoungTableauSeq<std::conditional_t<true, std::index_sequence<>, Row>...>,
294 YoungTableauSeq<Row...>,
295 typename YoungTableauSeq<Row...>::shape,
296 0,
297 0>::type;
298};
299
300template <>
301struct HooksHelper<YoungTableauSeq<>>
302{
303 using type = YoungTableauSeq<>;
304};
305
306template <class Tableau>
307using hooks_t = HooksHelper<Tableau>::type;
308
309// Sum the partial contributions of hooks to get hook lengths (= hooks+hooks_of_dual^T-1)
310template <class TableauHookLengths, class Tableau1, class Tableau2, std::size_t I>
311struct HookLengths;
312
313template <
314 class TableauHookLengths,
315 std::size_t HeadElemOfHeadRow1,
316 std::size_t... TailElemOfHeadRow1,
317 class... TailRow1,
318 std::size_t HeadElemOfHeadRow2,
319 std::size_t... TailElemOfHeadRow2,
320 class... TailRow2,
321 std::size_t I>
322struct HookLengths<
323 TableauHookLengths,
324 YoungTableauSeq<
325 std::index_sequence<HeadElemOfHeadRow1, TailElemOfHeadRow1...>,
326 TailRow1...>,
327 YoungTableauSeq<
328 std::index_sequence<HeadElemOfHeadRow2, TailElemOfHeadRow2...>,
329 TailRow2...>,
330 I>
331{
332 using type = std::conditional_t<
333 sizeof...(TailRow1) == 0 && sizeof...(TailElemOfHeadRow1) == 0,
334 add_cell_to_tableau_t<
335 TableauHookLengths,
336 I,
337 HeadElemOfHeadRow1 + HeadElemOfHeadRow2 - 1>,
338 std::conditional_t<
339 sizeof...(TailElemOfHeadRow1) == 0,
340 typename HookLengths<
341 add_cell_to_tableau_t<
342 TableauHookLengths,
343 I,
344 HeadElemOfHeadRow1 + HeadElemOfHeadRow2 - 1>,
345 YoungTableauSeq<TailRow1...>,
346 YoungTableauSeq<TailRow2...>,
347 I + 1>::type,
348 typename HookLengths<
349 add_cell_to_tableau_t<
350 TableauHookLengths,
351 I,
352 HeadElemOfHeadRow1 + HeadElemOfHeadRow2 - 1>,
353 YoungTableauSeq<
354 std::index_sequence<TailElemOfHeadRow1...>,
355 TailRow1...>,
356 YoungTableauSeq<
357 std::index_sequence<TailElemOfHeadRow2...>,
358 TailRow2...>,
359 I>::type>>;
360};
361
362template <class TableauHookLengths, class... TailRow1, class... TailRow2, std::size_t I>
363struct HookLengths<
364 TableauHookLengths,
365 YoungTableauSeq<std::index_sequence<>, TailRow1...>,
366 YoungTableauSeq<std::index_sequence<>, TailRow2...>,
367 I>
368{
369 using type = YoungTableauSeq<>;
370};
371
372template <class TableauHookLengths, std::size_t I>
373struct HookLengths<TableauHookLengths, YoungTableauSeq<>, YoungTableauSeq<>, I>
374{
375 using type = TableauHookLengths;
376};
377
378template <class Tableau1, class Tableau2>
379struct HookLengthsHelper;
380
381template <class... Row1, class Tableau2>
382struct HookLengthsHelper<YoungTableauSeq<Row1...>, Tableau2>
383{
384 using type = HookLengths<
385 YoungTableauSeq<std::conditional_t<true, std::index_sequence<>, Row1>...>,
386 YoungTableauSeq<Row1...>,
387 Tableau2,
388 0>::type;
389};
390
391template <>
392struct HookLengthsHelper<YoungTableauSeq<>, YoungTableauSeq<>>
393{
394 using type = YoungTableauSeq<>;
395};
396
397template <class Tableau1, class Tableau2>
398using hook_lengths_t = HookLengthsHelper<Tableau1, Tableau2>::type;
399
400// Compute irreductible representations dimension
401template <std::size_t Dimension, class TableauHookLengths, std::size_t I, std::size_t J>
402struct IrrepDim;
403
404template <
405 std::size_t Dimension,
406 std::size_t HeadElemOfHeadRow,
407 std::size_t... TailElemOfHeadRow,
408 class... TailRow,
409 std::size_t I,
410 std::size_t J>
411struct IrrepDim<
412 Dimension,
413 YoungTableauSeq<std::index_sequence<HeadElemOfHeadRow, TailElemOfHeadRow...>, TailRow...>,
414 I,
415 J>
416{
417 static consteval std::size_t run(double prod)
418 {
419 prod *= Dimension + J - I;
420 prod /= HeadElemOfHeadRow;
421 if constexpr (sizeof...(TailRow) == 0 && sizeof...(TailElemOfHeadRow) == 0) {
422 return prod;
423 } else if constexpr (sizeof...(TailElemOfHeadRow) == 0) {
424 return IrrepDim<Dimension, YoungTableauSeq<TailRow...>, I + 1, 0>::run(prod);
425 } else {
426 return IrrepDim<
427 Dimension,
428 YoungTableauSeq<std::index_sequence<TailElemOfHeadRow...>, TailRow...>,
429 I,
430 J + 1>::run(prod);
431 }
432 }
433};
434
435} // namespace detail
436
440template <std::size_t Dimension, misc::Specialization<YoungTableauSeq> TableauSeq>
441class YoungTableau
442{
443public:
444 using tableau_seq = TableauSeq;
445 using shape = typename TableauSeq::shape;
446
447private:
448 static constexpr std::size_t s_d = Dimension;
449 static constexpr std::size_t s_r = TableauSeq::rank;
450
451public:
452 using dual = YoungTableau<s_d, detail::dual_t<tableau_seq>>;
453 using hook_lengths = detail::hook_lengths_t<
454 detail::hooks_t<tableau_seq>,
455 detail::dual_t<detail::hooks_t<detail::dual_t<tableau_seq>>>>;
456
457private:
458 static constexpr std::size_t s_irrep_dim = detail::IrrepDim<s_d, hook_lengths, 0, 0>::run(1);
459
460 static constexpr std::array<char, 64> generate_tag_array();
461
462 static constexpr std::array<char, 64> s_tag_array = generate_tag_array();
463
464 static constexpr std::size_t s_tag_size
465 = std::find(s_tag_array.begin(), s_tag_array.end(), '\0') - s_tag_array.begin();
466
467 static constexpr std::string_view s_tag = {s_tag_array.data(), s_tag_size};
468
469 static consteval auto load_irrep();
470
471 static constexpr auto s_irrep = load_irrep();
472
473public:
474 YoungTableau();
475
476 static constexpr std::size_t dimension()
477 {
478 return s_d;
479 }
480
481 static constexpr std::size_t rank()
482 {
483 return s_r;
484 }
485
486 static constexpr std::size_t irrep_dim()
487 {
488 return s_irrep_dim;
489 }
490
491 static std::string tag()
492 {
493 return std::string(s_tag);
494 }
495
496private:
497 static constexpr std::size_t n_nonzeros_in_irrep()
498 {
499 return std::get<2>(std::get<0>(s_irrep)).size();
500 }
501
502public:
503 template <tensor::TensorNatIndex... Id>
504 using projector_domain = ddc::DiscreteDomain<tensor::prime<Id>..., Id...>;
505
506 template <tensor::TensorNatIndex... Id>
507 static auto projector();
508
509 template <class BasisId, class... Id>
510 static constexpr csr::Csr<n_nonzeros_in_irrep(), BasisId, Id...> u(
511 ddc::DiscreteDomain<Id...> restricted_domain)
512 {
513 if constexpr (n_nonzeros_in_irrep() != 0) {
514 ddc::DiscreteDomain<BasisId, Id...>
515 domain(ddc::DiscreteDomain<BasisId>(
516 ddc::DiscreteElement<BasisId>(0),
517 ddc::DiscreteVector<BasisId>(n_nonzeros_in_irrep())),
518 restricted_domain);
519
520 return csr::Csr<n_nonzeros_in_irrep(), BasisId, Id...>(
521 domain,
522 std::get<0>(std::get<0>(s_irrep)),
523 std::get<1>(std::get<0>(s_irrep)),
524 std::get<2>(std::get<0>(s_irrep)));
525 } else {
526 ddc::DiscreteDomain<BasisId, Id...>
527 domain(ddc::DiscreteDomain<BasisId>(
528 ddc::DiscreteElement<BasisId>(0),
529 ddc::DiscreteVector<BasisId>(0)),
530 restricted_domain);
531
532 return csr::Csr<n_nonzeros_in_irrep(), BasisId, Id...>(
533 domain,
534 std::array<std::size_t, BasisId::mem_size() + 1> {},
535 std::get<1>(std::get<0>(s_irrep)),
536 std::get<2>(std::get<0>(s_irrep)));
537 }
538 }
539
540 template <class BasisId, class... Id>
541 static constexpr csr::Csr<n_nonzeros_in_irrep(), BasisId, Id...> v(
542 ddc::DiscreteDomain<Id...> restricted_domain)
543 {
544 if constexpr (n_nonzeros_in_irrep() != 0) {
545 ddc::DiscreteDomain<BasisId, Id...>
546 domain(ddc::DiscreteDomain<BasisId>(
547 ddc::DiscreteElement<BasisId>(0),
548 ddc::DiscreteVector<BasisId>(n_nonzeros_in_irrep())),
549 restricted_domain);
550
551 return csr::Csr<n_nonzeros_in_irrep(), BasisId, Id...>(
552 domain,
553 std::get<0>(std::get<1>(s_irrep)),
554 std::get<1>(std::get<1>(s_irrep)),
555 std::get<2>(std::get<1>(s_irrep)));
556 } else {
557 ddc::DiscreteDomain<BasisId, Id...>
558 domain(ddc::DiscreteDomain<BasisId>(
559 ddc::DiscreteElement<BasisId>(0),
560 ddc::DiscreteVector<BasisId>(0)),
561 restricted_domain);
562
563 return csr::Csr<n_nonzeros_in_irrep(), BasisId, Id...>(
564 domain,
565 std::array<std::size_t, BasisId::mem_size() + 1> {},
566 std::get<1>(std::get<1>(s_irrep)),
567 std::get<2>(std::get<1>(s_irrep)));
568 }
569 }
570};
571
572namespace detail {
573
574template <std::size_t Dimension, class... TailRow, std::size_t I, std::size_t J>
575struct IrrepDim<Dimension, YoungTableauSeq<std::index_sequence<>, TailRow...>, I, J>
576{
577 static consteval std::size_t run(double prod)
578 {
579 return prod;
580 }
581};
582
583template <std::size_t Dimension, std::size_t I, std::size_t J>
584struct IrrepDim<Dimension, YoungTableauSeq<>, I, J>
585{
586 static consteval std::size_t run(double prod)
587 {
588 return prod;
589 }
590};
591
592// Gram-Schmidt approach to produce from vec an orthogonal vector to the vector space generated by basis
593template <
594 class ElementType,
595 class... Id,
596 class LayoutStridedPolicy,
597 class MemorySpace,
598 class BasisId>
599tensor::Tensor<ElementType, ddc::DiscreteDomain<Id...>, LayoutStridedPolicy, MemorySpace>
600orthogonalize(
601 tensor::Tensor<ElementType, ddc::DiscreteDomain<Id...>, LayoutStridedPolicy, MemorySpace>
602 tensor,
603 csr::CsrDynamic<BasisId, Id...> basis,
604 std::size_t max_basis_id)
605{
606 ddc::Chunk eigentensor_alloc(
607 ddc::DiscreteDomain<BasisId, Id...>(basis.domain()),
608 ddc::HostAllocator<double>());
609 tensor::Tensor eigentensor(eigentensor_alloc);
610 for (ddc::DiscreteElement<BasisId> elem : ddc::DiscreteDomain<BasisId>(
611 ddc::DiscreteElement<BasisId>(0),
612 ddc::DiscreteVector<BasisId>(max_basis_id))) {
613 csr::csr2dense(eigentensor, basis.get(elem));
614
615 ddc::Chunk scalar_prod_alloc(ddc::DiscreteDomain<> {}, ddc::HostAllocator<double>());
616 tensor::Tensor scalar_prod(scalar_prod_alloc);
617 tensor::tensor_prod(scalar_prod, tensor, eigentensor[ddc::DiscreteElement<BasisId>(0)]);
618 ddc::Chunk norm_squared_alloc(ddc::DiscreteDomain<> {}, ddc::HostAllocator<double>());
619 tensor::Tensor norm_squared(norm_squared_alloc);
620 tensor::tensor_prod(norm_squared, eigentensor, eigentensor);
621
622 eigentensor *= -1. * scalar_prod(ddc::DiscreteElement<>())
623 / norm_squared(ddc::DiscreteElement<>());
624 tensor += eigentensor[ddc::DiscreteElement<BasisId>(0)];
625 }
626 return tensor;
627}
628
629// Function to find the number at the given index `n` with increasing Hamming weight
630std::vector<bool> index_hamming_weight_code(std::size_t index, std::size_t length)
631{
632 std::size_t count = 0;
633 std::vector<bool> bits(length);
634 for (std::size_t hamming_weight = 0; hamming_weight < length; ++hamming_weight) {
635 std::fill(bits.begin(), bits.begin() + hamming_weight, true);
636 std::fill(bits.begin() + hamming_weight, bits.end(), false);
637 do {
638 if (count == index) {
639 return bits;
640 }
641 count++;
642 } while (std::prev_permutation(bits.begin(), bits.end()));
643 }
644 assert(false && "hamming weight code not found for index");
645 return bits;
646}
647
648// Dummy tag used by OrthonormalBasisSubspaceEigenvalueOne (coalescent dimension of the CsrDynamic storage)
649struct BasisId : tensor::TensorNaturalIndex<>
650{
651};
652
653template <class Ids>
654struct OrthonormalBasisSubspaceEigenvalueOne;
655
656template <class... Id>
657struct OrthonormalBasisSubspaceEigenvalueOne<tensor::TensorFullIndex<Id...>>
658{
659 template <class YoungTableau>
660 static std::pair<csr::CsrDynamic<BasisId, Id...>, csr::CsrDynamic<BasisId, Id...>> run(
661 YoungTableau tableau)
662 {
663 auto [proj_alloc, proj] = tableau.template projector<Id...>();
664
665 tensor::TensorAccessor<Id...> candidate_accessor;
666 ddc::DiscreteDomain<Id...> candidate_dom = candidate_accessor.domain();
667 ddc::Chunk candidate_alloc(candidate_dom, ddc::HostAllocator<double>());
668 tensor::Tensor candidate(candidate_alloc);
669 ddc::Chunk prod_alloc(
670 tensor::natural_tensor_prod_domain(proj.domain(), candidate.domain()),
671 ddc::HostAllocator<double>());
672 tensor::Tensor prod(prod_alloc);
673
674 ddc::DiscreteDomain<BasisId, Id...> basis_dom(
675 ddc::DiscreteDomain<BasisId>(
676 ddc::DiscreteElement<BasisId>(0),
677 ddc::DiscreteVector<BasisId>(tableau.irrep_dim())),
678 candidate_dom);
679 csr::CsrDynamic<BasisId, Id...> u(basis_dom);
680 csr::CsrDynamic<BasisId, Id...> v(basis_dom);
681 std::size_t n_irreps = 0;
682 std::size_t index = 0;
683 while (n_irreps < tableau.irrep_dim()) {
684 index++;
685 std::vector<bool> hamming_weight_code
686 = index_hamming_weight_code(index, (Id::size() * ...));
687
688 SIMILIE_DEBUG_LOG("similie_init_young_tableau_candidate");
689 ddc::parallel_for_each(
690 "similie_init_young_tableau_candidate",
691 Kokkos::DefaultHostExecutionSpace(),
692 candidate.domain(),
693 [&](ddc::DiscreteElement<Id...> elem) {
694 candidate(elem) = hamming_weight_code[(
695 (misc::detail::stride<Id, Id...>() * elem.template uid<Id>())
696 + ...)];
697 });
698
699 tensor::tensor_prod(prod, proj, candidate);
700 SIMILIE_DEBUG_LOG("similie_compute_young_tableau_candidate_copy");
701 Kokkos::deep_copy(
702 candidate.allocation_kokkos_view(),
703 prod.allocation_kokkos_view()); // We rely on Kokkos::deep_copy in place of ddc::parallel_deepcopy to avoid type verification of the type dimensions
704
705 orthogonalize(candidate, v, n_irreps);
706
707 if (ddc ::host_transform_reduce(
708 candidate.domain(),
709 false,
710 ddc::reducer::lor<bool>(),
711 [&](ddc::DiscreteElement<Id...> elem) { return candidate(elem) > 1e-6; })) {
712 ddc::Chunk
713 norm_squared_alloc(ddc::DiscreteDomain<> {}, ddc::HostAllocator<double>());
714 tensor::Tensor norm_squared(norm_squared_alloc);
715 tensor::tensor_prod(norm_squared, candidate, candidate);
716 SIMILIE_DEBUG_LOG("similie_normalize_young_tableau_candidate");
717 ddc::parallel_for_each(
718 "similie_normalize_young_tableau_candidate",
719 Kokkos::DefaultHostExecutionSpace(),
720 candidate.domain(),
721 [&](ddc::DiscreteElement<Id...> elem) {
722 candidate(elem) /= Kokkos::sqrt(norm_squared(ddc::DiscreteElement<>()));
723 });
724 // Not sure if u = v is correct in any case (ie. complex tensors ?)
725 u.push_back(candidate);
726 v.push_back(candidate);
727 n_irreps++;
728 std::cout << n_irreps << "/" << tableau.irrep_dim()
729 << " eigentensors found associated to the eigenvalue 1 for the Young "
730 "projector labelized "
731 << tableau.tag() << std::endl;
732 }
733 }
734 return std::pair<csr::CsrDynamic<BasisId, Id...>, csr::CsrDynamic<BasisId, Id...>>(u, v);
735 }
736};
737
738} // namespace detail
739
740template <std::size_t Dimension, misc::Specialization<YoungTableauSeq> TableauSeq>
741YoungTableau<Dimension, TableauSeq>::YoungTableau()
742{
743 // Check if the irrep is available in the dictionnary
744 {
745 std::ifstream file(IRREPS_DICT_PATH, std::ios::out | std::ios::binary);
746 std::string line;
747 while (!file.eof()) {
748 getline(file, line);
749 if (line == s_tag) {
750 file.close();
751 if (n_nonzeros_in_irrep() == 0) {
752 std::cout << "\033[1;31mIrrep " << s_tag << " in dimension " << s_d
753 << " required and found in dictionnary " << IRREPS_DICT_PATH
754 << " but the executable has been compiled without it. Please "
755 "recompile.\033[0m"
756 << std::endl;
757 }
758 return;
759 }
760 }
761 }
762
763 // If the current irrep is not found in the dictionnary, compute and dump it
764 std::cout << "\033[1;31mIrrep " << s_tag << " corresponding to the Young Tableau:\033[0m\n"
765 << *this << "\n\033[1;31min dimension " << s_d
766 << " required but not found in dictionnary " << IRREPS_DICT_PATH
767 << ". It will be computed, and you will have to recompile once it is done.\033[0m"
768 << std::endl;
769
770 auto [u, v]
771 = detail::OrthonormalBasisSubspaceEigenvalueOne<tensor::dummy_index_t<s_d, s_r>>::run(
772 *this);
773
774 std::ofstream file(IRREPS_DICT_PATH, std::ios::app | std::ios::binary);
775 if (!file) {
776 std::cerr << "Error opening file: " << IRREPS_DICT_PATH << std::endl;
777 return;
778 }
779 file << s_tag << "\n";
780 u.write(file);
781 v.write(file);
782 file << "\n";
783 file.close();
784 if (!file.good()) {
785 std::cerr << "Error occurred while writing to file " << IRREPS_DICT_PATH
786 << " while adding irrep " << s_tag << std::endl;
787 } else {
788 std::cout << "\033[1;32mIrrep " << s_tag << " added to the dictionnary " << IRREPS_DICT_PATH
789 << ".\033[0m \033[1;31mPlease recompile.\033[0m" << std::endl;
790 }
791}
792
793namespace detail {
794
795// Build index for symmetrizer (such that sym*proj is properly defined)
796template <class OId, class... Id>
797using symmetrizer_index_t = std::conditional_t<
798 (ddc::type_seq_rank_v<OId, ddc::detail::TypeSeq<Id...>> < (sizeof...(Id) / 2)),
799 tensor::prime<OId>,
800 ddc::type_seq_element_t<
801 static_cast<std::size_t>(
802 std::
803 max(static_cast<std::ptrdiff_t>(0),
804 static_cast<std::ptrdiff_t>(
805 ddc::type_seq_rank_v<OId, ddc::detail::TypeSeq<Id...>>
806 - (sizeof...(Id) / 2)))),
807 ddc::detail::TypeSeq<Id...>>>;
808
809// Functor to fill identity or transpose projectors
810template <std::size_t Dimension, std::size_t Rank, class... NaturalId>
811class TrFunctor
812{
813private:
814 using TensorType = tensor::Tensor<
815 double,
816 ddc::DiscreteDomain<NaturalId...>,
817 Kokkos::layout_right,
818 Kokkos::DefaultHostExecutionSpace::memory_space>;
819
820 TensorType t;
821 std::array<std::size_t, Rank> idx_to_permute;
822
823public:
824 TrFunctor(TensorType t_, std::array<std::size_t, Rank> idx_to_permute_)
825 : t(t_)
826 , idx_to_permute(idx_to_permute_)
827 {
828 }
829
830 void operator()(ddc::DiscreteElement<NaturalId...> elem) const
831 {
832 std::array<std::size_t, sizeof...(NaturalId)> elem_array = ddc::detail::array(elem);
833 bool has_to_be_one = true;
834 for (std::size_t i = 0; i < Rank; ++i) {
835 has_to_be_one
836 = has_to_be_one && (elem_array[i] == elem_array[idx_to_permute[i] + Rank]);
837 }
838 if (has_to_be_one) {
839 t(elem) = 1;
840 }
841 }
842};
843
844template <std::size_t Dimension, bool AntiSym, class... Id>
845static tensor::Tensor<
846 double,
847 ddc::DiscreteDomain<Id...>,
848 Kokkos::layout_right,
849 Kokkos::DefaultHostExecutionSpace::memory_space>
850fill_symmetrizer(
851 tensor::Tensor<
852 double,
853 ddc::DiscreteDomain<Id...>,
854 Kokkos::layout_right,
855 Kokkos::DefaultHostExecutionSpace::memory_space> sym,
856 std::array<std::size_t, sizeof...(Id) / 2> idx_to_permute)
857{
858 ddc::Chunk tr_alloc(sym.domain(), ddc::HostAllocator<double>());
859 tensor::Tensor tr(tr_alloc);
860 ddc::parallel_fill(tr, 0);
861 ddc::host_for_each(
862 tr.domain(),
863 TrFunctor<Dimension, sizeof...(Id) / 2, Id...>(tr, idx_to_permute));
864
865 if constexpr (AntiSym) {
866 tr *= static_cast<double>(misc::permutation_parity(idx_to_permute));
867 }
868 sym += tr;
869
870 return sym;
871}
872
873/*
874 Compute all permutations on a subset of indices in idx_to_permute, keep the
875 rest of the indices where they are.
876 */
877template <std::size_t Nt, std::size_t Ns>
878static std::vector<std::array<std::size_t, Nt>> permutations_subset(
879 std::array<std::size_t, Nt> t,
880 std::array<std::size_t, Ns> subset_values)
881{
882 std::array<std::size_t, Ns> subset_indices;
883 std::array<std::size_t, Ns> elements_to_permute;
884 int j = 0;
885 for (std::size_t i = 0; i < t.size(); ++i) {
886 if (std::find(subset_values.begin(), subset_values.end(), t[i] + 1)
887 != std::end(subset_values)) {
888 subset_indices[j] = i;
889 elements_to_permute[j++] = t[i];
890 }
891 }
892
893 std::vector<std::array<std::size_t, Nt>> result;
894 do {
895 std::array<std::size_t, Nt> tmp = t;
896 for (std::size_t i = 0; i < Ns; ++i) {
897 tmp[subset_indices[i]] = elements_to_permute[i];
898 }
899 result.push_back(tmp);
900 } while (std::next_permutation(elements_to_permute.begin(), elements_to_permute.end()));
901
902 return result;
903}
904
905// Compute projector
906template <class PartialTableauSeq, std::size_t Dimension, bool AntiSym = false>
907struct Projector;
908
909template <std::size_t... ElemOfHeadRow, class... TailRow, std::size_t Dimension, bool AntiSym>
910struct Projector<
911 YoungTableauSeq<std::index_sequence<ElemOfHeadRow...>, TailRow...>,
912 Dimension,
913 AntiSym>
914{
915 template <class... Id>
916 static tensor::Tensor<
917 double,
918 ddc::DiscreteDomain<Id...>,
919 Kokkos::layout_right,
920 Kokkos::DefaultHostExecutionSpace::memory_space>
921 run(tensor::Tensor<
922 double,
923 ddc::DiscreteDomain<Id...>,
924 Kokkos::layout_right,
925 Kokkos::DefaultHostExecutionSpace::memory_space> proj)
926 {
927 if constexpr (sizeof...(ElemOfHeadRow) >= 2) {
928 // Allocate & build a symmetric projector for the row
929 tensor::TensorAccessor<symmetrizer_index_t<Id, Id...>...> sym_accessor;
930 ddc::DiscreteDomain<symmetrizer_index_t<Id, Id...>...> sym_dom = sym_accessor.domain();
931
932 ddc::Chunk sym_alloc(sym_dom, ddc::HostAllocator<double>());
933 tensor::Tensor sym(sym_alloc);
934 ddc::parallel_fill(sym, 0);
935 std::array<std::size_t, sizeof...(Id) / 2> idx_to_permute;
936 for (std::size_t i = 0; i < sizeof...(Id) / 2; ++i) {
937 idx_to_permute[i] = i;
938 }
939 std::array<std::size_t, sizeof...(ElemOfHeadRow)> row_values {ElemOfHeadRow...};
940 auto idx_permutations = detail::permutations_subset(
941 idx_to_permute,
942 row_values); // TODO check https://indico.cern.ch/event/814040/contributions/3452485/attachments/1860434/3057354/psr_alcock.pdf page 6, it may be incomplete
943 for (std::size_t i = 0; i < idx_permutations.size(); ++i) {
944 fill_symmetrizer<
945 Dimension,
946 AntiSym,
947 symmetrizer_index_t<Id, Id...>...>(sym, idx_permutations[i]);
948 }
949
950 // Extract the symmetric part (for the row) of the projector (requires an intermediate prod tensor)
951 ddc::Chunk prod_alloc(
952 natural_tensor_prod_domain(sym.domain(), proj.domain()),
953 ddc::HostAllocator<double>());
954 tensor::Tensor prod(prod_alloc);
955 tensor::tensor_prod(prod, sym, proj);
956 Kokkos::deep_copy(
957 proj.allocation_kokkos_view(),
958 prod.allocation_kokkos_view()); // We rely on Kokkos::deep_copy in place of ddc::parallel_deepcopy to avoid type verification of the type dimensions
959 }
960 if constexpr (sizeof...(TailRow) == 0) {
961 return proj;
962 } else {
963 return Projector<YoungTableauSeq<TailRow...>, Dimension, AntiSym>::run(proj);
964 }
965 }
966};
967
968} // namespace detail
969
970template <std::size_t Dimension, misc::Specialization<YoungTableauSeq> TableauSeq>
971template <tensor::TensorNatIndex... Id>
972auto YoungTableau<Dimension, TableauSeq>::projector()
973{
974 static_assert(sizeof...(Id) == s_r);
975 tensor::TensorAccessor<tensor::prime<Id>..., Id...> proj_accessor;
976 ddc::DiscreteDomain<tensor::prime<Id>..., Id...> proj_dom = proj_accessor.domain();
977
978 // Allocate a projector and fill it as an identity tensor
979 ddc::Chunk proj_alloc(proj_dom, ddc::HostAllocator<double>());
980 tensor::Tensor proj(proj_alloc);
981 ddc::parallel_fill(proj, 0);
982 std::array<std::size_t, s_r> idx_to_permute;
983 for (std::size_t i = 0; i < s_r; ++i) {
984 idx_to_permute[i] = i;
985 }
986 ddc::host_for_each(
987 proj.domain(),
988 detail::TrFunctor<s_d, s_r, tensor::prime<Id>..., Id...>(proj, idx_to_permute));
989
990 // Build the projector
991 detail::Projector<tableau_seq, s_d>::run(proj);
992 detail::Projector<typename dual::tableau_seq, s_d, true>::run(proj);
993 return std::make_tuple(std::move(proj_alloc), proj);
994}
995
996// Load binary files and build u and v static constexpr Csr at compile-time
997namespace detail {
998
999template <std::size_t Line>
1000consteval std::string_view load_irrep_line_for_tag(std::string_view const tag)
1001{
1002 constexpr static unsigned char raw[] = {
1003#embed IRREPS_DICT_PATH
1004 };
1005
1006 constexpr static auto raw_chars = []() consteval {
1007 std::array<char, sizeof(raw)> chars {};
1008 for (std::size_t i = 0; i < sizeof(raw); ++i) {
1009 chars[i] = static_cast<char>(raw[i]);
1010 }
1011 return chars;
1012 }();
1013
1014 constexpr static std::string_view str(raw_chars.data(), raw_chars.size());
1015
1016 size_t tagPos = str.find(tag);
1017
1018 if (tagPos != std::string::npos) {
1019 std::size_t endOfTagLine = str.find('\n', tagPos);
1020 if (endOfTagLine != std::string::npos) {
1021 std::size_t block_size_pos = endOfTagLine + 1;
1022 for (std::size_t i = 0; i < Line; ++i) {
1023 if (block_size_pos + sizeof(std::size_t) > str.size()) {
1024 return "";
1025 }
1026 std::array<char, sizeof(std::size_t)> size_chars {};
1027 for (std::size_t j = 0; j < size_chars.size(); ++j) {
1028 size_chars[j] = str[block_size_pos + j];
1029 }
1030 std::size_t const block_size = std::bit_cast<std::size_t>(size_chars);
1031 block_size_pos += sizeof(std::size_t) + block_size;
1032 }
1033 if (block_size_pos + sizeof(std::size_t) > str.size()) {
1034 return "";
1035 }
1036 std::array<char, sizeof(std::size_t)> size_chars {};
1037 for (std::size_t j = 0; j < size_chars.size(); ++j) {
1038 size_chars[j] = str[block_size_pos + j];
1039 }
1040 std::size_t const block_size = std::bit_cast<std::size_t>(size_chars);
1041 std::size_t const block_start = block_size_pos + sizeof(std::size_t);
1042 if (block_start + block_size > str.size()) {
1043 return "";
1044 }
1045
1046 return std::string_view(str.data() + block_start, block_size);
1047 }
1048 }
1049
1050 return "";
1051}
1052
1053template <class Ids, std::size_t Offset>
1054struct LoadIrrepIdxForTag;
1055
1056template <std::size_t... I, std::size_t Offset>
1057struct LoadIrrepIdxForTag<std::index_sequence<I...>, Offset>
1058{
1059 static consteval std::array<std::string_view, sizeof...(I)> run(std::string_view const tag)
1060 {
1061 return std::array<std::string_view, sizeof...(I)> {
1062 load_irrep_line_for_tag<I + Offset>(tag)...};
1063 }
1064};
1065
1066template <class T, std::size_t N, std::size_t I = 0>
1067consteval std::array<T, N> bit_cast_array(
1068 std::array<T, N> vec,
1069 std::string_view const str,
1070 std::string_view const tag)
1071{
1072 if constexpr (I == N) {
1073 return vec;
1074 } else {
1075 std::array<char, sizeof(T) / sizeof(char)> chars = {};
1076 for (std::size_t j = 0; j < sizeof(T) / sizeof(char); ++j) {
1077 chars[j] = str[sizeof(T) / sizeof(char) * I + j];
1078 }
1079
1080 vec[I] = std::bit_cast<T>(
1081 chars); // We rely on std::bit_cast because std::reinterprest_cast is not constexpr
1082 return bit_cast_array<T, N, I + 1>(vec, str, tag);
1083 }
1084}
1085
1086template <class T, std::size_t N, std::size_t I = 0>
1087consteval std::array<T, N> bit_cast_array(std::string_view const str, std::string_view const tag)
1088{
1089 std::array<T, N> vec {};
1090 return bit_cast_array<T, N>(vec, str, tag);
1091}
1092
1093template <class T, std::size_t N, class Ids>
1094struct BitCastArrayOfArrays;
1095
1096template <class T, std::size_t N, std::size_t... I>
1097struct BitCastArrayOfArrays<T, N, std::index_sequence<I...>>
1098{
1099 static consteval std::array<std::array<T, N>, sizeof...(I)> run(
1100 std::array<std::string_view, sizeof...(I)> const str,
1101 std::string_view const tag)
1102 {
1103 return std::array<std::array<T, N>, sizeof...(I)> {bit_cast_array<T, N>(str[I], tag)...};
1104 }
1105};
1106
1107} // namespace detail
1108
1109template <std::size_t Dimension, misc::Specialization<YoungTableauSeq> TableauSeq>
1110consteval auto YoungTableau<Dimension, TableauSeq>::load_irrep()
1111{
1112 static constexpr std::string_view str_u_coalesc_idx(detail::load_irrep_line_for_tag<0>(s_tag));
1113 static constexpr std::array<std::string_view, s_r> str_u_idx(
1114 detail::LoadIrrepIdxForTag<std::make_index_sequence<s_r>, 1>::run(s_tag));
1115 static constexpr std::string_view str_u_values(detail::load_irrep_line_for_tag<s_r + 1>(s_tag));
1116 static constexpr std::string_view str_v_coalesc_idx(
1117 detail::load_irrep_line_for_tag<s_r + 2>(s_tag));
1118 static constexpr std::array<std::string_view, s_r> str_v_idx(
1119 detail::LoadIrrepIdxForTag<std::make_index_sequence<s_r>, s_r + 3>::run(s_tag));
1120 static constexpr std::string_view str_v_values(
1121 detail::load_irrep_line_for_tag<2 * s_r + 3>(s_tag));
1122
1123 if constexpr (str_u_values.size() != 0) {
1124 static constexpr std::array u_coalesc_idx = detail::bit_cast_array<
1125 std::size_t,
1126 str_u_coalesc_idx.size() / sizeof(std::size_t)>(str_u_coalesc_idx, s_tag);
1127 static constexpr std::array u_idx = detail::BitCastArrayOfArrays<
1128 std::size_t,
1129 str_u_idx[0].size() / sizeof(std::size_t),
1130 std::make_index_sequence<s_r>>::run(str_u_idx, s_tag);
1131 static constexpr std::array u_values = detail::
1132 bit_cast_array<double, str_u_values.size() / sizeof(double)>(str_u_values, s_tag);
1133 static constexpr std::array v_coalesc_idx = detail::bit_cast_array<
1134 std::size_t,
1135 str_v_coalesc_idx.size() / sizeof(std::size_t)>(str_v_coalesc_idx, s_tag);
1136 static constexpr std::array v_idx = detail::BitCastArrayOfArrays<
1137 std::size_t,
1138 str_v_idx[0].size() / sizeof(std::size_t),
1139 std::make_index_sequence<s_r>>::run(str_v_idx, s_tag);
1140 static constexpr std::array v_values = detail::
1141 bit_cast_array<double, str_v_values.size() / sizeof(double)>(str_v_values, s_tag);
1142 return std::make_pair(
1143 std::make_tuple(u_coalesc_idx, u_idx, u_values),
1144 std::make_tuple(v_coalesc_idx, v_idx, v_values));
1145 } else {
1146 return std::make_pair(
1147 std::make_tuple(
1148 std::array<std::size_t, 1> {0},
1149 std::array<std::array<std::size_t, 0>, s_r> {},
1150 std::array<double, 0> {}),
1151 std::make_tuple(
1152 std::array<std::size_t, 1> {0},
1153 std::array<std::array<std::size_t, 0>, s_r> {},
1154 std::array<double, 0> {}));
1155 }
1156}
1157
1158namespace detail {
1159
1160// Produce tag as a string (std::string features are limited at compile-time that's why we manipulate char arrays)
1161template <class Row>
1162struct YoungTableauRowToArray;
1163
1164template <std::size_t... RowElement>
1165struct YoungTableauRowToArray<std::index_sequence<RowElement...>>
1166{
1167 static constexpr auto run()
1168 {
1169 static constexpr std::array row = {RowElement...};
1170 return row;
1171 }
1172};
1173
1174template <class TableauSeq>
1175struct YoungTableauToArray;
1176
1177template <class... Row>
1178struct YoungTableauToArray<YoungTableauSeq<Row...>>
1179{
1180 static constexpr auto run()
1181 {
1182 static constexpr std::tuple tableau = {YoungTableauRowToArray<Row>::run()...};
1183 return tableau;
1184 }
1185};
1186
1187template <bool RowDelimiter>
1188struct RowToString
1189{
1190 template <std::size_t N>
1191 static constexpr std::array<char, 2 * N - !RowDelimiter> run(
1192 const std::array<std::size_t, N>& array)
1193 {
1194 std::array<char, 2 * N - !RowDelimiter> buf = {};
1195
1196 std::size_t current_pos = 0;
1197 for (std::size_t i = 0; i < N; ++i) {
1198 char temp[20];
1199 auto [ptr, ec] = std::to_chars(temp, temp + sizeof(temp), array[i]);
1200 if (ec != std::errc()) {
1201 throw std::runtime_error("Failed to convert number to string");
1202 }
1203
1204 for (char* p = temp; p < ptr; ++p) {
1205 buf[current_pos++] = *p;
1206 }
1207
1208 if (i < N - 1) {
1209 buf[current_pos++] = '_';
1210 }
1211 }
1212
1213 if constexpr (RowDelimiter) {
1214 buf[2 * N - 1] = 'l';
1215 }
1216
1217 return buf;
1218 }
1219};
1220
1221template <std::size_t... sizes>
1222constexpr auto concatenate(const std::array<char, sizes>&... arrays)
1223{
1224 std::array<char, (sizes + ...)> result;
1225 std::size_t index {};
1226
1227 ((std::copy_n(arrays.begin(), sizes, result.begin() + index), index += sizes), ...);
1228
1229 return result;
1230}
1231
1232template <class RowIdx>
1233struct ArrayToString;
1234
1235template <std::size_t... RowId>
1236struct ArrayToString<std::index_sequence<RowId...>>
1237{
1238 template <class Tuple>
1239 static constexpr auto run(Tuple const tableau)
1240 {
1241 return concatenate(
1242 RowToString < RowId != sizeof...(RowId) - 1 > ::run(std::get<RowId>(tableau))...);
1243 }
1244};
1245
1246template <std::size_t size>
1247constexpr auto add_dimension(const std::array<char, size>& array, std::size_t d)
1248{
1249 std::array<char, size + 2> result;
1250 char temp[1];
1251 std::to_chars(temp, temp + 1, d);
1252 std::copy_n(array.begin(), size, result.begin());
1253 result[size] = 'd';
1254 result[size + 1] = temp[0];
1255
1256 return result;
1257}
1258
1259} // namespace detail
1260
1261template <std::size_t Dimension, misc::Specialization<YoungTableauSeq> TableauSeq>
1262constexpr std::array<char, 64> YoungTableau<Dimension, TableauSeq>::generate_tag_array()
1263{
1264 static constexpr std::tuple tableau = detail::YoungTableauToArray<tableau_seq>::run();
1265 constexpr auto row_str_wo_dimension
1266 = detail::ArrayToString<std::make_index_sequence<tableau_seq::shape::size()>>::run(
1267 tableau);
1268 constexpr auto row_str = detail::add_dimension(row_str_wo_dimension, s_d);
1269 std::array<char, 64> tag {};
1270 static_assert(row_str.size() < tag.size());
1271 std::copy_n(row_str.begin(), row_str.size(), tag.begin());
1272 return tag;
1273}
1274
1275namespace detail {
1276
1277// Print Young tableau in a string
1278template <class TableauSeq>
1279struct PrintYoungTableauSeq;
1280
1281template <std::size_t HeadRowHeadElement, std::size_t... HeadRowTailElement, class... TailRow>
1282struct PrintYoungTableauSeq<
1283 YoungTableauSeq<std::index_sequence<HeadRowHeadElement, HeadRowTailElement...>, TailRow...>>
1284{
1285 static std::string run(std::string str)
1286 {
1287 str += std::to_string(HeadRowHeadElement) + " ";
1288 if constexpr (sizeof...(TailRow) == 0 && sizeof...(HeadRowTailElement) == 0) {
1289 } else if constexpr (sizeof...(HeadRowTailElement) == 0) {
1290 str += "\n";
1291 str = PrintYoungTableauSeq<YoungTableauSeq<TailRow...>>::run(str);
1292 } else {
1293 str = PrintYoungTableauSeq<
1294 YoungTableauSeq<std::index_sequence<HeadRowTailElement...>, TailRow...>>::
1295 run(str);
1296 }
1297 return str;
1298 }
1299};
1300
1301template <>
1302struct PrintYoungTableauSeq<YoungTableauSeq<>>
1303{
1304 static std::string run(std::string str)
1305 {
1306 return str;
1307 }
1308};
1309
1310} // namespace detail
1311
1312template <std::size_t Dimension, misc::Specialization<YoungTableauSeq> TableauSeq>
1313std::ostream& operator<<(std::ostream& os, YoungTableau<Dimension, TableauSeq> const& tableau)
1314{
1315 std::string str = "";
1316 os << detail::PrintYoungTableauSeq<TableauSeq>::run(str);
1317 return os;
1318}
1319
1320} // namespace young_tableau
1321
1322} // namespace sil
1323
1324#endif
sil::csr::operator<<
std::ostream & operator<<(std::ostream &os, Csr< N, TensorIndex... > const &csr)
Definition csr.hpp:191
sil::tensor::natural_tensor_prod_domain
natural_tensor_prod_domain_t< Dom1, Dom2 > natural_tensor_prod_domain(Dom1 dom1, Dom2 dom2)
Definition tensor_impl.hpp:812
sil
The top-level namespace of SimiLie.
Definition csr.hpp:15