Dynare++ tensor library: modernize the Symmetry class

We now use a initializer list constructor for creating symmetries of the form
$y^n$, $y^n u^m$, $y^nu^m\sigma^k$.

The constructor taking a single integer is used to initialize a symmetry of a
given length.

Similar changes are made to IntSequence.

This behavior is similar to std::vector.

dynare++/integ/cc/product.cc

@@ -10,7 +10,7 @@
 
 prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
   : prodq(q), level(l), npoints(q.uquad.numPoints(l)),
-    jseq{q.dimen(), 0},
+    jseq(q.dimen(), 0),
     end_flag(false),
     sig{q.dimen()},
     p{q.dimen()}

dynare++/integ/cc/smolyak.cc

@@ -12,7 +12,7 @@
 
 smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum)
   : smolq(q), isummand(isum),
-    jseq{q.dimen(), 0},
+    jseq(q.dimen(), 0),
     sig{q.dimen()},
     p{q.dimen()}
 {

dynare++/integ/testing/tests.cc

@@ -256,7 +256,7 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level
 
   // check against theoretical moments
   UNormalMoments moments(imom, msq);
-  smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
+  smol_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
   std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
   return smol_out.getMax() < 1.e-7;
 }
@@ -285,7 +285,7 @@ TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level
 
   // check against theoretical moments
   UNormalMoments moments(imom, msq);
-  prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
+  prod_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
   std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
   return prod_out.getMax() < 1.e-7;
 }

dynare++/kord/approximation.cc

@@ -80,7 +80,7 @@ Approximation::approxAtSteady()
 {
   model.calcDerivativesAtSteady();
   FirstOrder fo(model.nstat(), model.npred(), model.nboth(), model.nforw(),
-                model.nexog(), *(model.getModelDerivatives().get(Symmetry(1))),
+                model.nexog(), *(model.getModelDerivatives().get(Symmetry{1})),
                 journal, qz_criterium);
   KORD_RAISE_IF_X(!fo.isStable(),
                   "The model is not Blanchard-Kahn stable",
@@ -293,21 +293,21 @@ Approximation::calcStochShift(Vector &out, double at_sigma) const
     {
       if (KOrder::is_even(d))
         {
-          Symmetry sym(0, d, 0, 0);
+          Symmetry sym{0, d, 0, 0};
 
           // calculate $F_{u'^d}$ via |ZAuxContainer|
           FGSTensor *ten = new FGSTensor(ypart.ny(), TensorDimens(sym, nvs));
           ten->zeros();
           for (int l = 1; l <= d; l++)
             {
-              const FSSparseTensor *f = model.getModelDerivatives().get(Symmetry(l));
+              const FSSparseTensor *f = model.getModelDerivatives().get(Symmetry{l});
               zaux.multAndAdd(*f, *ten);
             }
 
           // multiply with shocks and add to result
-          FGSTensor *tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry(0, 0, 0, 0), nvs));
+          FGSTensor *tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry{0, 0, 0, 0}, nvs));
           tmp->zeros();
-          ten->contractAndAdd(1, *tmp, *(mom.get(Symmetry(d))));
+          ten->contractAndAdd(1, *tmp, *(mom.get(Symmetry{d})));
 
           out.add(pow(at_sigma, d)/dfac, tmp->getData());
           delete ten;
@@ -360,8 +360,8 @@ Approximation::check(double at_sigma) const
 TwoDMatrix *
 Approximation::calcYCov() const
 {
-  const TwoDMatrix &gy = *(rule_ders->get(Symmetry(1, 0, 0, 0)));
-  const TwoDMatrix &gu = *(rule_ders->get(Symmetry(0, 1, 0, 0)));
+  const TwoDMatrix &gy = *(rule_ders->get(Symmetry{1, 0, 0, 0}));
+  const TwoDMatrix &gu = *(rule_ders->get(Symmetry{0, 1, 0, 0}));
   TwoDMatrix G(model.numeq(), model.numeq());
   G.zeros();
   G.place(gy, 0, model.nstat());

dynare++/kord/decision_rule.hh

@@ -217,11 +217,11 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
           int j = d-i;
           int kfact = 1;
           _Ttensor tmp(ypart.ny(),
-                       TensorDimens(Symmetry(i, j), tns));
+                       TensorDimens(Symmetry{i, j}, tns));
           tmp.zeros();
           for (int k = 0; k+d <= g.getMaxDim(); k++, kfact *= k)
             {
-              Symmetry sym(i, j, 0, k);
+              Symmetry sym{i, j, 0, k};
               if (g.check(sym))
                 {
                   double mult = pow(sigma, k)/dfact/kfact;
@@ -539,7 +539,7 @@ DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
   fillTensors(g, sigma);
   _Tparent yspol(ypart.nstat, ypart.nys(), *this);
   bigf = new _Tparent((const _Tparent &) yspol);
-  _Ttensym *frst = bigf->get(Symmetry(1));
+  _Ttensym *frst = bigf->get(Symmetry{1});
   for (int i = 0; i < ypart.nys(); i++)
     frst->get(i, i) = frst->get(i, i) - 1;
   bigfder = new _Tparent(*bigf, 0);
@@ -572,9 +572,9 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
       int kfact = 1;
       for (int k = 0; d+k <= g.getMaxDim(); k++, kfact *= k)
         {
-          if (g.check(Symmetry(d, 0, 0, k)))
+          if (g.check(Symmetry{d, 0, 0, k}))
             {
-              const _Ttensor *ten = g.get(Symmetry(d, 0, 0, k));
+              const _Ttensor *ten = g.get(Symmetry{d, 0, 0, k});
               double mult = pow(sigma, k)/dfact/kfact;
               g_yd->add(mult, *ten);
             }

dynare++/kord/first_order.hh

@@ -72,10 +72,10 @@ public:
   {
     IntSequence nvs(4);
     nvs[0] = fo.ypart.nys(); nvs[1] = fo.nu; nvs[2] = fo.nu; nvs[3] = 1;
-    _Ttensor *ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(1, 0, 0, 0), nvs));
+    _Ttensor *ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
     ten->zeros(); ten->add(1.0, fo.gy);
     this->insert(ten);
-    ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(0, 1, 0, 0), nvs));
+    ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
     ten->zeros(); ten->add(1.0, fo.gu);
     this->insert(ten);
   }

dynare++/kord/global_check.cc

@@ -95,9 +95,9 @@ ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
   hss = new FTensorPolynomial(dr_ss, ytmp_star);
   ConstVector ysteady_ss(dr.c->getSteady(), model->nstat()+model->npred(),
                          model->nboth()+model->nforw());
-  if (hss->check(Symmetry(0)))
+  if (hss->check(Symmetry{0}))
     {
-      hss->get(Symmetry(0))->getData().add(1.0, ysteady_ss);
+      hss->get(Symmetry{0})->getData().add(1.0, ysteady_ss);
     }
   else
     {

dynare++/kord/korder.cc

@@ -239,9 +239,9 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     _um(maxk, v), _fm(_um), f(fcont),
-    matA(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),
-    matS(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),
-    matB(*(f.get(Symmetry(1))), _uZstack.getStackSizes()),
+    matA(*(f.get(Symmetry{1})), _uZstack.getStackSizes(), gy, ypart),
+    matS(*(f.get(Symmetry{1})), _uZstack.getStackSizes(), gy, ypart),
+    matB(*(f.get(Symmetry{1})), _uZstack.getStackSizes()),
     journal(jr)
 {
   KORD_RAISE_IF(gy.ncols() != ypart.nys(),
@@ -265,20 +265,20 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
   // put $g_y$ and $g_u$ to the container
   /* Note that $g_\sigma$ is zero by the nature and we do not insert it to
      the container. We insert a new physical copies. */
-  UGSTensor *tgy = new UGSTensor(ny, TensorDimens(Symmetry(1, 0, 0, 0), nvs));
+  UGSTensor *tgy = new UGSTensor(ny, TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
   tgy->getData() = gy.getData();
   insertDerivative<unfold>(tgy);
-  UGSTensor *tgu = new UGSTensor(ny, TensorDimens(Symmetry(0, 1, 0, 0), nvs));
+  UGSTensor *tgu = new UGSTensor(ny, TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
   tgu->getData() = gu.getData();
   insertDerivative<unfold>(tgu);
 
   // put $G_y$, $G_u$ and $G_{u'}$ to the container
   /* Also note that since $g_\sigma$ is zero, so $G_\sigma$. */
-  UGSTensor *tGy = faaDiBrunoG<unfold>(Symmetry(1, 0, 0, 0));
+  UGSTensor *tGy = faaDiBrunoG<unfold>(Symmetry{1, 0, 0, 0});
   G<unfold>().insert(tGy);
-  UGSTensor *tGu = faaDiBrunoG<unfold>(Symmetry(0, 1, 0, 0));
+  UGSTensor *tGu = faaDiBrunoG<unfold>(Symmetry{0, 1, 0, 0});
   G<unfold>().insert(tGu);
-  UGSTensor *tGup = faaDiBrunoG<unfold>(Symmetry(0, 0, 1, 0));
+  UGSTensor *tGup = faaDiBrunoG<unfold>(Symmetry{0, 0, 1, 0});
   G<unfold>().insert(tGup);
 }
 
@@ -306,7 +306,7 @@ KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
     {
       KORD_RAISE_IF(!der.isFinite(),
                     "RHS of Sylverster is not finite");
-      TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry(1, 0, 0, 0))));
+      TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry{1, 0, 0, 0})));
       GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(),
                             ypart.nstat+ypart.npred,
                             matA.getData(), matB.getData(),

dynare++/kord/korder.hh

@@ -461,7 +461,7 @@ template<int t>
 void
 KOrder::recover_y(int i)
 {
-  Symmetry sym(i, 0, 0, 0);
+  Symmetry sym{i, 0, 0, 0};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -475,7 +475,7 @@ KOrder::recover_y(int i)
 
   insertDerivative<t>(g_yi);
 
-  _Ttensor *gss_y = gss<t>().get(Symmetry(1, 0, 0, 0));
+  _Ttensor *gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
   gs<t>().multAndAdd(*gss_y, *G_yi);
   _Ttensor *gss_yi = gss<t>().get(sym);
   gs<t>().multAndAdd(*gss_yi, *G_yi);
@@ -496,7 +496,7 @@ template <int t>
 void
 KOrder::recover_yu(int i, int j)
 {
-  Symmetry sym(i, j, 0, 0);
+  Symmetry sym{i, j, 0, 0};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -508,7 +508,7 @@ KOrder::recover_yu(int i, int j)
   matA.multInv(*g_yiuj);
   insertDerivative<t>(g_yiuj);
 
-  gs<t>().multAndAdd(*(gss<t>().get(Symmetry(1, 0, 0, 0))), *G_yiuj);
+  gs<t>().multAndAdd(*(gss<t>().get(Symmetry{1, 0, 0, 0})), *G_yiuj);
 }
 
 /* Here we solve
@@ -535,7 +535,7 @@ template <int t>
 void
 KOrder::recover_ys(int i, int j)
 {
-  Symmetry sym(i, 0, 0, j);
+  Symmetry sym{i, 0, 0, j};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -595,7 +595,7 @@ template <int t>
 void
 KOrder::recover_yus(int i, int j, int k)
 {
-  Symmetry sym(i, j, 0, k);
+  Symmetry sym{i, j, 0, k};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -662,7 +662,7 @@ template <int t>
 void
 KOrder::recover_s(int i)
 {
-  Symmetry sym(0, 0, 0, i);
+  Symmetry sym{0, 0, 0, i};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -710,7 +710,7 @@ KOrder::fillG(int i, int j, int k)
     {
       if (is_even(k-m))
         {
-          _Ttensor *G_yiujupms = faaDiBrunoG<t>(Symmetry(i, j, m, k-m));
+          _Ttensor *G_yiujupms = faaDiBrunoG<t>(Symmetry{i, j, m, k-m});
           G<t>().insert(G_yiujupms);
         }
     }
@@ -727,12 +727,12 @@ template <int t>
 _Ttensor *
 KOrder::calcD_ijk(int i, int j, int k) const
 {
-  _Ttensor *res =       new _Ttensor(ny, TensorDimens(Symmetry(i, j, 0, 0), nvs));
+  _Ttensor *res =       new _Ttensor(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
   res->zeros();
   if (is_even(k))
     {
-      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry(i, j, k, 0));
-      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(k))));
+      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry{i, j, k, 0});
+      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{k})));
       delete tmp;
     }
   return res;
@@ -749,13 +749,13 @@ template <int t>
 _Ttensor *
 KOrder::calcE_ijk(int i, int j, int k) const
 {
-  _Ttensor *res = new _Ttensor(ny, TensorDimens(Symmetry(i, j, 0, 0), nvs));
+  _Ttensor *res = new _Ttensor(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
   res->zeros();
   for (int n = 2; n <= k-1; n += 2)
     {
-      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry(i, j, n, k-n));
+      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry{i, j, n, k-n});
       tmp->mult((double) (Tensor::noverk(k, n)));
-      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(n))));
+      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{n})));
       delete tmp;
     }
   return res;
@@ -861,7 +861,7 @@ KOrder::check(int dim) const
   // check for $F_{y^iu^j}=0
   for (int i = 0; i <= dim; i++)
     {
-      Symmetry sym(dim-i, i, 0,  0);
+      Symmetry sym{dim-i, i, 0,  0};
       _Ttensor *r = faaDiBrunoZ<t>(sym);
       double err = r->getData().getMax();
       JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
@@ -879,7 +879,7 @@ KOrder::check(int dim) const
       int k = (*si)[2];
       if (i+j > 0 && k > 0)
         {
-          Symmetry sym(i, j, 0, k);
+          Symmetry sym{i, j, 0, k};
           _Ttensor *r = faaDiBrunoZ<t>(sym);
           _Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
           r->add(1.0, *D_ijk);
@@ -894,7 +894,7 @@ KOrder::check(int dim) const
     }
 
   // check for $F_{\sigma^i}+D_i+E_i=0
-  _Ttensor *r = faaDiBrunoZ<t>(Symmetry(0, 0, 0, dim));
+  _Ttensor *r = faaDiBrunoZ<t>(Symmetry{0, 0, 0, dim});
   _Ttensor *D_k = calcD_k<t>(dim);
   r->add(1.0, *D_k);
   delete D_k;
@@ -902,7 +902,7 @@ KOrder::check(int dim) const
   r->add(1.0, *E_k);
   delete E_k;
   double err = r->getData().getMax();
-  Symmetry sym(0, 0, 0, dim);
+  Symmetry sym{0, 0, 0, dim};
   JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
   if (err > maxerror)
     maxerror = err;

dynare++/kord/korder_stoch.cc

@@ -35,7 +35,7 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     f(fcont),
-    matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1, 0, 0, 0))),
+    matA(*(fcont.get(Symmetry{1})), _uZstack.getStackSizes(), *(hh.get(Symmetry{1, 0, 0, 0})),
          ypart)
 {
   nvs[0] = ypart.nys();
@@ -56,7 +56,7 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     f(fcont),
-    matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1, 0, 0, 0))),
+    matA(*(fcont.get(Symmetry{1})), _uZstack.getStackSizes(), *(hh.get(Symmetry{1, 0, 0, 0})),
          ypart)
 {
   nvs[0] = ypart.nys();

dynare++/kord/korder_stoch.hh

@@ -87,7 +87,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
       for (int i = 0; i <= d; i++)
         {
           int p = d-i;
-          Symmetry sym(i, 0, 0, p);
+          Symmetry sym{i, 0, 0, p};
           _Ttensor *ten = new _Ttensor(r, TensorDimens(sym, nvs));
 
           // calculate derivative $h_{y^i\sigma^p}$
@@ -105,14 +105,14 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
               for (int m = 0; i+m+n+k <= maxd; m++, mfac *= m)
                 {
                   double mult = (pow(at_sigma, m)*povern)/mfac;
-                  Symmetry sym_mn(i, m+n, 0, k);
+                  Symmetry sym_mn{i, m+n, 0, k};
                   if (m+n == 0 && g.check(sym_mn))
                     ten->add(mult, *(g.get(sym_mn)));
                   if (m+n > 0 && KOrder::is_even(m+n) && g.check(sym_mn))
                     {
                       _Ttensor gtmp(*(g.get(sym_mn)));
                       gtmp.mult(mult);
-                      gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry(m+n))));
+                      gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry{m+n})));
                     }
                 }
             }
@@ -187,8 +187,8 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
       for (int i = 0; i <= d; i++)
         {
           int k = d-i;
-          if (g_int.check(Symmetry(i, 0, 0, k)))
-            ten->addSubTensor(*(g_int.get(Symmetry(i, 0, 0, k))));
+          if (g_int.check(Symmetry{i, 0, 0, k}))
+            ten->addSubTensor(*(g_int.get(Symmetry{i, 0, 0, k})));
         }
       g_int_sym.insert(ten);
     }
@@ -224,13 +224,13 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
   true_nvs[1] = nu; true_nvs[2] = nu;
   for (int d = 1; d <= maxd; d++)
     {
-      if (g_int_cent.check(Symmetry(d)))
+      if (g_int_cent.check(Symmetry{d}))
         {
           for (int i = 0; i <= d; i++)
             {
-              Symmetry sym(i, 0, 0, d-i);
+              Symmetry sym{i, 0, 0, d-i};
               IntSequence coor(sym, pp);
-              _Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry(d))), ss, coor,
+              _Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry{d})), ss, coor,
                                            TensorDimens(sym, true_nvs));
               this->insert(ten);
             }
@@ -274,7 +274,7 @@ GXContainer<_Ttype>::getType(int i, const Symmetry &s) const
   if (i == 2)
     return _Stype::zero;
   if (i == 3)
-    if (s == Symmetry(0, 0, 0, 1))
+    if (s == Symmetry{0, 0, 0, 1})
       return _Stype::unit;
     else
       return _Stype::zero;
@@ -319,12 +319,12 @@ ZXContainer<_Ttype>::getType(int i, const Symmetry &s) const
     else
       return _Stype::matrix;
   if (i == 2)
-    if (s == Symmetry(1, 0, 0, 0))
+    if (s == Symmetry{1, 0, 0, 0})
       return _Stype::unit;
     else
       return _Stype::zero;
   if (i == 3)
-    if (s == Symmetry(0, 1, 0, 0))
+    if (s == Symmetry{0, 1, 0, 0})
       return _Stype::unit;
     else
       return _Stype::zero;

dynare++/kord/tests.cc

@@ -287,7 +287,7 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
   for (int d = 1; d <= maxdim; d++)
     {
       printf("\ttensor fill for dim=%d is:   %3.2f %%\n",
-             d, c.get(Symmetry(d))->getFillFactor()*100.0);
+             d, c.get(Symmetry{d})->getFillFactor()*100.0);
     }
   Journal jr("out.txt");
   KOrder kord(nstat, npred, nboth, nforw, c, gy, gu, v, jr);

dynare++/src/dynare3.cc

@@ -353,7 +353,7 @@ DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms &a,
 void
 DynareDerEvalLoader::load(int i, int iord, const int *vars, double res)
 {
-  FSSparseTensor *t = md.get(Symmetry(iord));
+  FSSparseTensor *t = md.get(Symmetry{iord});
   IntSequence s(iord, 0);
   for (int j = 0; j < iord; j++)