diff --git a/dynare++/integ/cc/smolyak.cc b/dynare++/integ/cc/smolyak.cc index 6fec35e9c3a83fc6e80f6945ef8b07a523bfae2e..36e1c267cd8e1087230da11d0ef9496843b6346d 100644 --- a/dynare++/integ/cc/smolyak.cc +++ b/dynare++/integ/cc/smolyak.cc @@ -112,12 +112,11 @@ SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq) // todo: check |l>1|, |l>=d| // todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()| int cum = 0; - SymmetrySet ss(l-1, d+1); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(l-1, d+1)) { - if ((*si)[d] <= d-1) + if (si[d] <= d-1) { - IntSequence lev((const IntSequence &)*si, 0, d); + IntSequence lev((const IntSequence &) si, 0, d); lev.add(1); levels.push_back(lev); IntSequence levpts(d); @@ -171,12 +170,11 @@ int SmolyakQuadrature::calcNumEvaluations(int lev) const { int cum = 0; - SymmetrySet ss(lev-1, dim+1); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(lev-1, dim+1)) { - if ((*si)[dim] <= dim-1) + if (si[dim] <= dim-1) { - IntSequence lev((const IntSequence &)*si, 0, dim); + IntSequence lev((const IntSequence &) si, 0, dim); lev.add(1); IntSequence levpts(dim); for (int i = 0; i < dim; i++) diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc index 7f4be316ac72c53d65d724babbc327b6d65ec038..553f1eb690648d57e81d9d2e5e0ca97b486d0ee5 100644 --- a/dynare++/kord/korder.cc +++ b/dynare++/kord/korder.cc @@ -342,30 +342,27 @@ KOrder::switchToFolded() int maxdim = g<unfold>().getMaxDim(); for (int dim = 1; dim <= maxdim; dim++) - { - SymmetrySet ss(dim, 4); - for (symiterator si(ss); !si.isEnd(); ++si) - { - if ((*si)[2] == 0 && g<unfold>().check(*si)) - { - auto *ft = new FGSTensor(*(g<unfold>().get(*si))); - insertDerivative<fold>(ft); - if (dim > 1) - { - gss<unfold>().remove(*si); - gs<unfold>().remove(*si); - g<unfold>().remove(*si); - } - } - if (G<unfold>().check(*si)) - { - auto *ft = new FGSTensor(*(G<unfold>().get(*si))); - G<fold>().insert(ft); - if (dim > 1) - { - G<fold>().remove(*si); - } - } - } - } + for (auto &si : SymmetrySet(dim, 4)) + { + if (si[2] == 0 && g<unfold>().check(si)) + { + auto *ft = new FGSTensor(*(g<unfold>().get(si))); + insertDerivative<fold>(ft); + if (dim > 1) + { + gss<unfold>().remove(si); + gs<unfold>().remove(si); + g<unfold>().remove(si); + } + } + if (G<unfold>().check(si)) + { + auto *ft = new FGSTensor(*(G<unfold>().get(si))); + G<fold>().insert(ft); + if (dim > 1) + { + G<fold>().remove(si); + } + } + } } diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh index 7c504c00dfefc9781c1cd241e9288515d04cd601..784e0c8df3dff54bd477d0bbd63cbcf1e9fd77cb 100644 --- a/dynare++/kord/korder.hh +++ b/dynare++/kord/korder.hh @@ -871,12 +871,11 @@ KOrder::check(int dim) const } // check for $F_{y^iu^ju'^k}+D_{ijk}+E_{ijk}=0$ - SymmetrySet ss(dim, 3); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(dim, 3)) { - int i = (*si)[0]; - int j = (*si)[1]; - int k = (*si)[2]; + int i = si[0]; + int j = si[1]; + int k = si[2]; if (i+j > 0 && k > 0) { Symmetry sym{i, j, 0, k}; diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh index 65705afe329e5d3827bca8b96da330364de08fb3..109b1752fa54f89ef8225d5624a35df35723df94 100644 --- a/dynare++/kord/korder_stoch.hh +++ b/dynare++/kord/korder_stoch.hh @@ -516,18 +516,17 @@ KOrderStoch::performStep(int order) int maxd = g<t>().getMaxDim(); KORD_RAISE_IF(order-1 != maxd && (order != 1 || maxd != -1), "Wrong order for KOrderStoch::performStep"); - SymmetrySet ss(order, 4); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(order, 4)) { - if ((*si)[2] == 0) + if (si[2] == 0) { JournalRecordPair pa(journal); - pa << "Recovering symmetry " << *si << endrec; + pa << "Recovering symmetry " << si << endrec; - _Ttensor *G_sym = faaDiBrunoG<t>(*si); + _Ttensor *G_sym = faaDiBrunoG<t>(si); G<t>().insert(G_sym); - _Ttensor *g_sym = faaDiBrunoZ<t>(*si); + _Ttensor *g_sym = faaDiBrunoZ<t>(si); g_sym->mult(-1.0); matA.multInv(*g_sym); g<t>().insert(g_sym); diff --git a/dynare++/tl/cc/stack_container.cc b/dynare++/tl/cc/stack_container.cc index 311f94b058bf4596a114e4e1c445c4657e933e36..3bd2ef64f170ca74dca0e55af7b108485f266481 100644 --- a/dynare++/tl/cc/stack_container.cc +++ b/dynare++/tl/cc/stack_container.cc @@ -42,10 +42,10 @@ FoldedStackContainer::multAndAdd(int dim, const FGSContainer &c, FGSTensor &out) "Wrong symmetry length of container for FoldedStackContainer::multAndAdd"); sthread::detach_thread_group gr; - SymmetrySet ss(dim, c.num()); - for (symiterator si(ss); !si.isEnd(); ++si) - if (c.check(*si)) - gr.insert(std::make_unique<WorkerFoldMAADense>(*this, *si, c, out)); + + for (auto &si : SymmetrySet(dim, c.num())) + if (c.check(si)) + gr.insert(std::make_unique<WorkerFoldMAADense>(*this, si, c, out)); gr.run(); } @@ -395,10 +395,9 @@ UnfoldedStackContainer::multAndAdd(int dim, const UGSContainer &c, "Wrong symmetry length of container for UnfoldedStackContainer::multAndAdd"); sthread::detach_thread_group gr; - SymmetrySet ss(dim, c.num()); - for (symiterator si(ss); !si.isEnd(); ++si) - if (c.check(*si)) - gr.insert(std::make_unique<WorkerUnfoldMAADense>(*this, *si, c, out)); + for (auto &si : SymmetrySet(dim, c.num())) + if (c.check(si)) + gr.insert(std::make_unique<WorkerUnfoldMAADense>(*this, si, c, out)); gr.run(); } diff --git a/dynare++/tl/cc/symmetry.cc b/dynare++/tl/cc/symmetry.cc index 9b94df492c84ea305b815387c51edb275878cb28..23a7e878c608d6d645d4dcadc80d171a4c731096 100644 --- a/dynare++/tl/cc/symmetry.cc +++ b/dynare++/tl/cc/symmetry.cc @@ -51,53 +51,42 @@ Symmetry::isFull() const return count <= 1; } -/* Here we construct the beginning of the |symiterator|. The first - symmetry index is 0. If length is 2, the second index is the - dimension, otherwise we create the subordinal symmetry set and its - beginning as subordinal |symiterator|. */ +/* Construct a symiterator of given dimension, starting from the given + symmetry. */ -symiterator::symiterator(SymmetrySet &ss) - : s(ss), end_flag(false) +symiterator::symiterator(int dim_arg, Symmetry run_arg) + : dim{dim_arg}, run(std::move(run_arg)) { - s.sym()[0] = 0; - if (s.size() == 2) - s.sym()[1] = s.dimen(); - else - { - subs = std::make_unique<SymmetrySet>(s, s.dimen()); - subit = std::make_unique<symiterator>(*subs); - } } /* Here we move to the next symmetry. We do so only, if we are not at the end. If length is 2, we increase lower index and decrease upper index, otherwise we increase the subordinal symmetry. If we got to the end, we recreate the subordinal symmetry set and set the subordinal - iterator to the beginning. At the end we test, if we are not at the - end. This is recognized if the lowest index exceeded the dimension. */ + iterator to the beginning. */ symiterator & symiterator::operator++() { - if (!end_flag) + if (run[0] == dim) + run[0]++; // Jump to the past-the-end iterator + else if (run.size() == 2) { - if (s.size() == 2) - { - s.sym()[0]++; - s.sym()[1]--; - } - else + run[0]++; + run[1]--; + } + else + { + symiterator subit{dim-run[0], Symmetry(run, run.size()-1)}; + ++subit; + if (run[1] == dim-run[0]+1) { - ++(*subit); - if (subit->isEnd()) - { - s.sym()[0]++; - subs = std::make_unique<SymmetrySet>(s, s.dimen()-s.sym()[0]); - subit = std::make_unique<symiterator>(*subs); - } + run[0]++; + run[1] = 0; + /* subit is equal to the past-the-end iterator, so the range + 2..(size()-1) is already set to 0 */ + run[run.size()-1] = dim-run[0]; } - if (s.sym()[0] == s.dimen()+1) - end_flag = true; } return *this; } diff --git a/dynare++/tl/cc/symmetry.hh b/dynare++/tl/cc/symmetry.hh index 2a4976ce791a6c1817664f77788fca0149418c14..e5e97bc1387d21cdb797b9f29de3349666213a01 100644 --- a/dynare++/tl/cc/symmetry.hh +++ b/dynare++/tl/cc/symmetry.hh @@ -105,95 +105,77 @@ public: bool isFull() const; }; -/* The class |SymmetrySet| defines a set of symmetries of the given - length having given dimension. It does not store all the symmetries, - rather it provides a storage for one symmetry, which is changed as an - adjoint iterator moves. - - The iterator class is |symiterator|. It is implemented - recursively. The iterator object, when created, creates subordinal - iterator, which iterates over a symmetry set whose length is one less, - and dimension is the former dimension. When the subordinal iterator - goes to its end, the superordinal iterator increases left most index in - the symmetry, resets the subordinal symmetry set with different - dimension, and iterates through the subordinal symmetry set until its - end, and so on. That's why we provide also |SymmetrySet| constructor - for construction of a subordinal symmetry set. +/* This is an iterator that iterates over all symmetries of given length and + dimension (see the SymmetrySet class for details). - The typical usage of the abstractions for |SymmetrySet| and - |symiterator| is as follows: - - \kern0.3cm - \centerline{|for (symiterator si(SymmetrySet(6, 4)); !si.isEnd(); ++si) {body}|} - \kern0.3cm + The beginning iterator is (0, …, 0, dim). + Increasing it gives (0, … , 1, dim-1) + The just-before-end iterator is (dim, 0, …, 0) + The past-the-end iterator is (dim+1, 0, …, 0) - \noindent It goes through all symmetries of size 4 having dimension - 6. One can use |*si| as the symmetry in the body. */ + The constructor creates the iterator which starts from the given symmetry + symmetry (beginning). */ -class SymmetrySet +class symiterator { + const int dim; Symmetry run; - int dim; public: - SymmetrySet(int d, int length) - : run(length), dim(d) - { - } - SymmetrySet(SymmetrySet &s, int d) - : run(s.run, s.size()-1), dim(d) - { - } - int - dimen() const - { - return dim; - } + symiterator(int dim_arg, Symmetry run_arg); + ~symiterator() = default; + symiterator &operator++(); const Symmetry & - sym() const + operator*() const { return run; } - Symmetry & - sym() + bool + operator=(const symiterator &it) { - return run; + return dim == it.dim && run == it.run; } - int - size() const + bool + operator!=(const symiterator &it) { - return run.size(); + return !operator=(it); } }; -/* The logic of |symiterator| was described in |@<|SymmetrySet| class - declaration@>|. Here we only comment that: the class has a reference - to the |SymmetrySet| only to know dimension and for access of its - symmetry storage. Further we have pointers to subordinal |symiterator| - and its |SymmetrySet|. These are pointers, since the recursion ends at - length equal to 2, in which case these pointers are uninitialized. +/* The class |SymmetrySet| defines a set of symmetries of the given length + having given dimension (i.e. it represents all the lists of integers of + length "len" and of sum equal to "dim"). It does not store all the + symmetries, it is just a convenience class for iteration. - The constructor creates the iterator which initializes to the first - symmetry (beginning). */ + The typical usage of the abstractions for |SymmetrySet| and + |symiterator| is as follows: -class symiterator + for (auto &si : SymmetrySet(6, 4)) + + It goes through all symmetries of lenght 4 having dimension 6. One can use + "si" as the symmetry in the body. */ + +class SymmetrySet { - SymmetrySet &s; - std::unique_ptr<symiterator> subit; - std::unique_ptr<SymmetrySet> subs; - bool end_flag; public: - symiterator(SymmetrySet &ss); - ~symiterator() = default; - symiterator &operator++(); - bool - isEnd() const + const int len; + const int dim; + SymmetrySet(int dim_arg, int len_arg) + : len(len_arg), dim(dim_arg) { - return end_flag; } - const Symmetry & - operator*() const + symiterator + begin() const + { + Symmetry run(len); + run[len-1] = dim; + return { dim, run }; + } + symiterator + end() const { - return s.sym(); + Symmetry run(len); + run[0] = dim+1; + return { dim, run }; } }; diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc index 391f146d0a3b6a99e62104d9bca027d2bf0a4e65..a8a2c269d3fa118005f340a246fe671a3f3387fb 100644 --- a/dynare++/tl/testing/monoms.cc +++ b/dynare++/tl/testing/monoms.cc @@ -471,13 +471,12 @@ SparseDerivGenerator::SparseDerivGenerator( for (int dim = 1; dim <= maxdimen; dim++) { - SymmetrySet ss(dim, 4); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(dim, 4)) { - bigg->insert(bigg_m.deriv(*si)); - rcont->insert(r.deriv(*si)); - if ((*si)[2] == 0) - g->insert(g_m.deriv(*si)); + bigg->insert(bigg_m.deriv(si)); + rcont->insert(r.deriv(si)); + if (si[2] == 0) + g->insert(g_m.deriv(si)); } ts[dim-1] = f.deriv(dim); diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc index d574fcb01bb044f69ddc477c8933a8e9326ab2b6..3bb93092d4a51842eb042c0535ee268fa0299a52 100644 --- a/dynare++/tl/testing/tests.cc +++ b/dynare++/tl/testing/tests.cc @@ -362,21 +362,19 @@ TestRunnable::fold_zcont(int nf, int ny, int nu, int nup, int nbigg, for (int d = 2; d <= dim; d++) { - SymmetrySet ss(d, 4); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(d, 4)) { - printf("\tSymmetry: "); (*si).print(); - FGSTensor res(nf, TensorDimens(*si, nvs)); + printf("\tSymmetry: "); + si.print(); + FGSTensor res(nf, TensorDimens(si, nvs)); res.getData().zeros(); clock_t stime = clock(); - for (int l = 1; l <= (*si).dimen(); l++) - { - zc.multAndAdd(*(dg.ts[l-1]), res); - } + for (int l = 1; l <= si.dimen(); l++) + zc.multAndAdd(*(dg.ts[l-1]), res); stime = clock() - stime; printf("\t\ttime for symmetry: %8.4g\n", ((double) stime)/CLOCKS_PER_SEC); - const FGSTensor *mres = dg.rcont->get(*si); + const FGSTensor *mres = dg.rcont->get(si); res.add(-1.0, *mres); double normtmp = res.getData().getMax(); printf("\t\terror normMax: %10.6g\n", normtmp); @@ -419,22 +417,20 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg, for (int d = 2; d <= dim; d++) { - SymmetrySet ss(d, 4); - for (symiterator si(ss); !si.isEnd(); ++si) + for (auto &si : SymmetrySet(d, 4)) { - printf("\tSymmetry: "); (*si).print(); - UGSTensor res(nf, TensorDimens(*si, nvs)); + printf("\tSymmetry: "); + si.print(); + UGSTensor res(nf, TensorDimens(si, nvs)); res.getData().zeros(); clock_t stime = clock(); - for (int l = 1; l <= (*si).dimen(); l++) - { - zc.multAndAdd(*(dg.ts[l-1]), res); - } + for (int l = 1; l <= si.dimen(); l++) + zc.multAndAdd(*(dg.ts[l-1]), res); stime = clock() - stime; printf("\t\ttime for symmetry: %8.4g\n", ((double) stime)/CLOCKS_PER_SEC); FGSTensor fold_res(res); - const FGSTensor *mres = dg.rcont->get(*si); + const FGSTensor *mres = dg.rcont->get(si); fold_res.add(-1.0, *mres); double normtmp = fold_res.getData().getMax(); printf("\t\terror normMax: %10.6g\n", normtmp);