korder_stoch.hh 17.9 KB
Newer Older
1 2 3 4 5 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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
// Copyright 2005, Ondra Kamenik

// Higher order at stochastic steady

/* This file defines a number of classes of which |KOrderStoch| is the
   main purpose. Basically, |KOrderStoch| calculates first and higher
   order Taylor expansion of a policy rule at $\sigma>0$ with explicit
   forward $g^{**}$. More formally, we have to solve a policy rule $g$
   from
   $$E_t\left[f(g^{**}(g^*(y^*_t,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y^*,u_t)\right]$$
   As an introduction in {\tt approximation.hweb} argues, $g^{**}$ at
   tine $t+1$ must be given from outside. Let the explicit
   $E_t(g^{**}(y^*,u_{t+1},\sigma)$ be equal to $h(y^*,\sigma)$. Then we
   have to solve
   $$f(h(g^*(y^*,u,\sigma),\sigma),g(y,u,\sigma),y,u),$$
   which is much easier than fully implicit system for $\sigma=0$.

   Besides the class |KOrderStoch|, we declare here also classes for the
   new containers corresponding to
   $f(h(g^*(y^*,u,\sigma),\sigma),g(y,u,\sigma),y,u)$. Further, we
   declare |IntegDerivs| and |StochForwardDerivs| classes which basically
   calculate $h$ as an extrapolation based on an approximation to $g$ at
   lower $\sigma$. */

#include "korder.hh"
#include "faa_di_bruno.hh"
#include "journal.hh"

/* This class is a container, which has a specialized constructor
   integrating the policy rule at given $\sigma$. */

template <int t>
class IntegDerivs : public ctraits<t>::Tgss
{
public:
  IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const __Tm &mom,
              double at_sigma);
};

/* This constructor integrates a rule (namely its $g^{**}$ part) with
   respect to $u=\tilde\sigma\eta$, and stores to the object the
   derivatives of this integral $h$ at $(y^*,u,\sigma)=(\tilde
   y^*,0,\tilde\sigma)$. The original container of $g^{**}$, the moments of
   the stochastic shocks |mom| and the $\tilde\sigma$ are input.

   The code follows the following derivation
   \def\lims{\vbox{\baselineskip=0pt\lineskip=1pt
   \setbox0=\hbox{$\scriptstyle n+k=p$}\hbox to\wd0{\hss$\scriptstyle m=0$\hss}\box0}}
   $$
   \eqalign{h(y,\sigma)&=E_t\left[g(y,u',\sigma)\right]=\cr
   &=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+j+k=d}\pmatrix{d\cr i,j,k}\left[g_{y^iu^j\sigma^k}\right]
   (y^*-\tilde y^*)^i\sigma^j\Sigma^j(\sigma-\tilde\sigma)^k\cr
   &=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+m+n+k=d}\pmatrix{d\cr i,m+n,k}
   \left[g_{y^iu^{m+n}\sigma^k}\right]
   \hat y^{*i}\Sigma^{m+n}\pmatrix{m+n\cr m,n}{\tilde\sigma}^m\hat\sigma^{k+n}\cr
   &=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+m+n+k=d}\pmatrix{d\cr i,m,n,k}
   \left[g_{y^iu^{m+n}\sigma^k}\right]
   \Sigma^{m+n}{\tilde\sigma}^m\hat y^{*i}\hat\sigma^{k+n}\cr
   &=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+p=d}\sum_{\lims}\pmatrix{d\cr i,m,n,k}
   \left[g_{y^iu^{m+n}\sigma^k}\right]
   \Sigma^{m+n}{\tilde\sigma}^m\hat y^{*i}\hat\sigma^{k+n}\cr
   &=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+p=d}\pmatrix{d\cr i,p}
   \left[\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
   \left[g_{y^iu^{m+n}\sigma^k}\right]
   \Sigma^{m+n}{\tilde\sigma}^m\right]\hat y^{*i}\hat\sigma^{k+n},
   }
   $$
   where $\pmatrix{a\cr b_1,\ldots, b_n}$ is a generalized combination
   number, $p=k+n$, $\hat\sigma=\sigma-\tilde\sigma$, $\hat
   y^*=y^*-\tilde y$, and we dropped writing the multidimensional indexes
   in Einstein summation.

   This implies that
   $$h_{y^i\sigma^p}=\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
   \left[g_{y^iu^{m+n}\sigma^k}\right]
   \Sigma^{m+n}{\tilde\sigma}^m$$
   and this is exactly what the code does. */

template <int t>
IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const __Tm &mom,
                            double at_sigma)
  : ctraits<t>::Tgss(4)
{
  int maxd = g.getMaxDim();
  for (int d = 1; d <= maxd; d++)
    {
      for (int i = 0; i <= d; i++)
        {
          int p = d-i;
90
          Symmetry sym{i, 0, 0, p};
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
          _Ttensor *ten = new _Ttensor(r, TensorDimens(sym, nvs));

          // calculate derivative $h_{y^i\sigma^p}$
          /* This code calculates
             $$h_{y^i\sigma^p}=\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
             \left[g_{y^iu^{m+n}\sigma^k}\right]
             \Sigma^{m+n}{\tilde\sigma}^m$$
             and stores it in |ten|. */
          ten->zeros();
          for (int n = 0; n <= p; n++)
            {
              int k = p-n;
              int povern = Tensor::noverk(p, n);
              int mfac = 1;
              for (int m = 0; i+m+n+k <= maxd; m++, mfac *= m)
                {
                  double mult = (pow(at_sigma, m)*povern)/mfac;
108
                  Symmetry sym_mn{i, m+n, 0, k};
109 110 111 112 113 114
                  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);
115
                      gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry{m+n})));
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 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
                    }
                }
            }

          this->insert(ten);
        }
    }
}

/* This class calculates an extrapolation of expectation of forward
   derivatives. It is a container, all calculations are done in a
   constructor.

   The class calculates derivatives of $E[g(y*,u,\sigma)]$ at $(\bar
   y^*,\bar\sigma)$. The derivatives are extrapolated based on
   derivatives at $(\tilde y^*,\tilde\sigma)$. */

template <int t>
class StochForwardDerivs : public ctraits<t>::Tgss
{
public:
  StochForwardDerivs(const PartitionY &ypart, int nu,
                     const _Tgss &g, const __Tm &m,
                     const Vector &ydelta, double sdelta,
                     double at_sigma);
};

/* This is the constructor which performs the integration and the
   extrapolation. Its parameters are: |g| is the container of derivatives
   at $(\tilde y,\tilde\sigma)$; |m| are the moments of stochastic
   shocks; |ydelta| is a difference of the steady states $\bar y-\tilde
   y$; |sdelta| is a difference between new sigma and old sigma
   $\bar\sigma-\tilde\sigma$, and |at_sigma| is $\tilde\sigma$. There is
   no need of inputing the $\tilde y$.

   We do the operation in four steps:
   \orderedlist
   \li Integrate $g^{**}$, the derivatives are at $(\tilde y,\tilde\sigma)$
   \li Form the (full symmetric) polynomial from the derivatives stacking
   $\left[\matrix{y^*\cr\sigma}\right]$
   \li Centralize this polynomial about $(\bar y,\bar\sigma)$
   \li Recover general symmetry tensors from the (full symmetric) polynomial
   \endorderedlist */

template <int t>
StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
                                          const _Tgss &g, const __Tm &m,
                                          const Vector &ydelta, double sdelta,
                                          double at_sigma)
  : ctraits<t>::Tgss(4)
{
  int maxd = g.getMaxDim();
  int r = ypart.nyss();

  // make |g_int| be integral of $g^{**}$ at $(\tilde y,\tilde\sigma)$
  /* This simply constructs |IntegDerivs| class. Note that the |nvs| of
     the tensors has zero dimensions for shocks, this is because we need to
     make easily stacks of the form $\left[\matrix{y^*\cr\sigma}\right]$ in
     the next step. */
  IntSequence nvs(4);
  nvs[0] = ypart.nys(); nvs[1] = 0; nvs[2] = 0; nvs[3] = 1;
  IntegDerivs<t> g_int(r, nvs, g, m, at_sigma);

  // make |g_int_sym| be full symmetric polynomial from |g_int|
  /* Here we just form a polynomial whose unique variable corresponds to
     $\left[\matrix{y^*\cr\sigma}\right]$ stack. */
  _Tpol g_int_sym(r, ypart.nys()+1);
  for (int d = 1; d <= maxd; d++)
    {
185
      auto *ten = new _Ttensym(r, ypart.nys()+1, d);
186 187 188 189
      ten->zeros();
      for (int i = 0; i <= d; i++)
        {
          int k = d-i;
190 191
          if (g_int.check(Symmetry{i, 0, 0, k}))
            ten->addSubTensor(*(g_int.get(Symmetry{i, 0, 0, k})));
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
        }
      g_int_sym.insert(ten);
    }

  // make |g_int_cent| the centralized polynomial about $(\bar y,\bar\sigma)$
  /* Here we centralize the polynomial to $(\bar y,\bar\sigma)$ knowing
     that the polynomial was centralized about $(\tilde
     y,\tilde\sigma)$. This is done by derivating and evaluating the
     derivated polynomial at $(\bar y-\tilde
     y,\bar\sigma-\tilde\sigma)$. The stack of this vector is |delta| in
     the code. */
  Vector delta(ypart.nys()+1);
  Vector dy(delta, 0, ypart.nys());
  ConstVector dy_in(ydelta, ypart.nstat, ypart.nys());
  dy = dy_in;
  delta[ypart.nys()] = sdelta;
  _Tpol g_int_cent(r, ypart.nys()+1);
  for (int d = 1; d <= maxd; d++)
    {
      g_int_sym.derivative(d-1);
      _Ttensym *der = g_int_sym.evalPartially(d, delta);
      g_int_cent.insert(der);
    }

  // pull out general symmetry tensors from |g_int_cent|
  /* Here we only recover the general symmetry derivatives from the full
     symmetric polynomial. Note that the derivative get the true |nvs|. */
  IntSequence ss(4);
  ss[0] = ypart.nys(); ss[1] = 0; ss[2] = 0; ss[3] = 1;
  IntSequence pp(4);
  pp[0] = 0; pp[1] = 1; pp[2] = 2; pp[3] = 3;
  IntSequence true_nvs(nvs);
  true_nvs[1] = nu; true_nvs[2] = nu;
  for (int d = 1; d <= maxd; d++)
    {
227
      if (g_int_cent.check(Symmetry{d}))
228 229 230
        {
          for (int i = 0; i <= d; i++)
            {
231
              Symmetry sym{i, 0, 0, d-i};
232
              IntSequence coor(sym, pp);
233
              _Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry{d})), ss, coor,
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249
                                           TensorDimens(sym, true_nvs));
              this->insert(ten);
            }
        }
    }
}

/* This container corresponds to $h(g^*(y,u,\sigma),\sigma)$. Note that
   in our application, the $\sigma$ as a second argument to $h$ will be
   its fourth variable in symmetry, so we have to do four member stack
   having the second and third stack dummy. */

template <class _Ttype>
class GXContainer : public GContainer<_Ttype>
{
public:
250 251 252
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename StackContainer<_Ttype>::_Ctype;
  using itype = typename StackContainer<_Ttype>::itype;
253 254 255 256
  GXContainer(const _Ctype *gs, int ngs, int nu)
    : GContainer<_Ttype>(gs, ngs, nu)
  {
  }
257
  itype getType(int i, const Symmetry &s) const override;
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
};

/* This routine corresponds to this stack:
   $$\left[\matrix{g^*(y,u,\sigma)\cr dummy\cr dummy\cr\sigma}\right]$$ */

template <class _Ttype>
typename GXContainer<_Ttype>::itype
GXContainer<_Ttype>::getType(int i, const Symmetry &s) const
{
  if (i == 0)
    if (s[2] > 0)
      return _Stype::zero;
    else
      return _Stype::matrix;
  if (i == 1)
    return _Stype::zero;
  if (i == 2)
    return _Stype::zero;
  if (i == 3)
277
    if (s == Symmetry{0, 0, 0, 1})
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
      return _Stype::unit;
    else
      return _Stype::zero;

  KORD_RAISE("Wrong stack index in GXContainer::getType");
}

/* This container corresponds to $f(H(y,u,\sigma),g(y,u,sigma),y,u)$,
   where the $H$ has the size (number of rows) as $g^{**}$. Since it is
   very simmilar to |ZContainer|, we inherit form it and override only
   |getType| method. */

template <class _Ttype>
class ZXContainer : public ZContainer<_Ttype>
{
public:
294 295 296
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename StackContainer<_Ttype>::_Ctype;
  using itype = typename StackContainer<_Ttype>::itype;
297 298 299 300
  ZXContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng, int ny, int nu)
    : ZContainer<_Ttype>(gss, ngss, g, ng, ny, nu)
  {
  }
301
  itype getType(int i, const Symmetry &s) const override;
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
};

/* This |getType| method corresponds to this stack:
   $$\left[\matrix{H(y,u,\sigma)\cr g(y,u,\sigma)\cr y\cr u}\right]$$ */

template <class _Ttype>
typename ZXContainer<_Ttype>::itype
ZXContainer<_Ttype>::getType(int i, const Symmetry &s) const
{
  if (i == 0)
    if (s[2] > 0)
      return _Stype::zero;
    else
      return _Stype::matrix;
  if (i == 1)
    if (s[2] > 0)
      return _Stype::zero;
    else
      return _Stype::matrix;
  if (i == 2)
322
    if (s == Symmetry{1, 0, 0, 0})
323 324 325 326
      return _Stype::unit;
    else
      return _Stype::zero;
  if (i == 3)
327
    if (s == Symmetry{0, 1, 0, 0})
328 329 330 331 332 333 334 335 336 337
      return _Stype::unit;
    else
      return _Stype::zero;

  KORD_RAISE("Wrong stack index in ZXContainer::getType");
}

class UnfoldedGXContainer : public GXContainer<UGSTensor>, public UnfoldedStackContainer
{
public:
338
  using _Ctype = TensorContainer<UGSTensor>;
339 340 341 342 343 344 345 346 347
  UnfoldedGXContainer(const _Ctype *gs, int ngs, int nu)
    : GXContainer<UGSTensor>(gs, ngs, nu)
  {
  }
};

class FoldedGXContainer : public GXContainer<FGSTensor>, public FoldedStackContainer
{
public:
348
  using _Ctype = TensorContainer<FGSTensor>;
349 350 351 352 353 354 355 356 357
  FoldedGXContainer(const _Ctype *gs, int ngs, int nu)
    : GXContainer<FGSTensor>(gs, ngs, nu)
  {
  }
};

class UnfoldedZXContainer : public ZXContainer<UGSTensor>, public UnfoldedStackContainer
{
public:
358
  using _Ctype = TensorContainer<UGSTensor>;
359 360 361 362 363 364 365 366 367
  UnfoldedZXContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng, int ny, int nu)
    : ZXContainer<UGSTensor>(gss, ngss, g, ng, ny, nu)
  {
  }
};

class FoldedZXContainer : public ZXContainer<FGSTensor>, public FoldedStackContainer
{
public:
368
  using _Ctype = TensorContainer<FGSTensor>;
369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494
  FoldedZXContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng, int ny, int nu)
    : ZXContainer<FGSTensor>(gss, ngss, g, ng, ny, nu)
  {
  }
};

/* This matrix corresponds to
   $$\left[f_{y}\right]+ \left[0
   \left[f_{y^{**}_+}\right]\cdot\left[h^{**}_{y^*}\right] 0\right]$$
   This is very the same as |MatrixA|, the only difference that the
   |MatrixA| is constructed from whole $h_{y^*}$, not only from
   $h^{**}_{y^*}$, hence the new abstraction. */

class MatrixAA : public PLUMatrix
{
public:
  MatrixAA(const FSSparseTensor &f, const IntSequence &ss,
           const TwoDMatrix &gyss, const PartitionY &ypart);
};

/* This class calculates derivatives of $g$ given implicitly by
   $f(h(g^*(y,u,\sigma),\sigma),g(y,u,\sigma),y,u)$, where $h(y,\sigma)$
   is given from outside.

   Structurally, the class is very similar to |KOrder|, but calculations
   are much easier. The two constructors construct an object from sparse
   derivatives of $f$, and derivatives of $h$. The caller must ensure
   that the both derivatives are done at the same point.

   The calculation for order $k$ (including $k=1$) is done by a call
   |performStep(k)|. The derivatives can be retrived by |getFoldDers()|
   or |getUnfoldDers()|. */

class KOrderStoch
{
protected:
  IntSequence nvs;
  PartitionY ypart;
  Journal &journal;
  UGSContainer _ug;
  FGSContainer _fg;
  UGSContainer _ugs;
  FGSContainer _fgs;
  UGSContainer _uG;
  FGSContainer _fG;
  const UGSContainer *_uh;
  const FGSContainer *_fh;
  UnfoldedZXContainer _uZstack;
  FoldedZXContainer _fZstack;
  UnfoldedGXContainer _uGstack;
  FoldedGXContainer _fGstack;
  const TensorContainer<FSSparseTensor> &f;
  MatrixAA matA;
public:
  KOrderStoch(const PartitionY &ypart, int nu, const TensorContainer<FSSparseTensor> &fcont,
              const FGSContainer &hh, Journal &jr);
  KOrderStoch(const PartitionY &ypart, int nu, const TensorContainer<FSSparseTensor> &fcont,
              const UGSContainer &hh, Journal &jr);
  template <int t>
  void performStep(int order);
  const FGSContainer &
  getFoldDers() const
  {
    return _fg;
  }
  const UGSContainer &
  getUnfoldDers() const
  {
    return _ug;
  }
protected:
  template <int t>
  _Ttensor *faaDiBrunoZ(const Symmetry &sym) const;
  template <int t>
  _Ttensor *faaDiBrunoG(const Symmetry &sym) const;

  // convenience access methods
  template<int t>
  _Tg&g();
  template<int t>
  const _Tg&g() const;
  template<int t>
  _Tgs&gs();
  template<int t>
  const _Tgs&gs() const;
  template<int t>
  const _Tgss&h() const;
  template<int t>
  _TG&G();
  template<int t>
  const _TG&G() const;
  template<int t>
  _TZXstack&Zstack();
  template<int t>
  const _TZXstack&Zstack() const;
  template<int t>
  _TGXstack&Gstack();
  template<int t>
  const _TGXstack&Gstack() const;
};

/* This calculates a derivative of $f(G(y,u,\sigma),g(y,u,\sigma),y,u)$
   of a given symmetry. */

template <int t>
_Ttensor *
KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
{
  JournalRecordPair pa(journal);
  pa << "Faa Di Bruno ZX container for " << sym << endrec;
  _Ttensor *res = new _Ttensor(ypart.ny(), TensorDimens(sym, nvs));
  FaaDiBruno bruno(journal);
  bruno.calculate(Zstack<t>(), f, *res);
  return res;
}

/* This calculates a derivative of
   $G(y,u,\sigma)=h(g^*(y,u,\sigma),\sigma)$ of a given symmetry. */

template <int t>
_Ttensor *
KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
{
  JournalRecordPair pa(journal);
  pa << "Faa Di Bruno GX container for " << sym << endrec;
  TensorDimens tdims(sym, nvs);
495
  auto *res = new _Ttensor(ypart.nyss(), tdims);
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
  FaaDiBruno bruno(journal);
  bruno.calculate(Gstack<t>(), h<t>(), *res);
  return res;
}

/* This retrieves all $g$ derivatives of a given dimension from implicit
   $f(h(g^*(y,u,\sigma),\sigma),g(y,u,\sigma),y,u)$. It suppose that all
   derivatives of smaller dimensions have been retrieved.

   So, we go through all symmetries $s$, calculate $G_s$ conditional on
   $g_s=0$, insert the derivative to the $G$ container, then calculate
   $F_s$ conditional on $g_s=0$. This is a righthand side. The left hand
   side is $matA\cdot g_s$. The $g_s$ is retrieved as
   $$g_s=-matA^{-1}\cdot RHS.$$ Finally we have to update $G_s$ by
   calling |Gstack<t>().multAndAdd(1, h<t>(), *G_sym)|. */

template <int t>
void
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)
    {
      if ((*si)[2] == 0)
        {
          JournalRecordPair pa(journal);
          pa << "Recovering symmetry " << *si << endrec;

          _Ttensor *G_sym = faaDiBrunoG<t>(*si);
          G<t>().insert(G_sym);

          _Ttensor *g_sym = faaDiBrunoZ<t>(*si);
          g_sym->mult(-1.0);
          matA.multInv(*g_sym);
          g<t>().insert(g_sym);
          gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *g_sym));

          Gstack<t>().multAndAdd(1, h<t>(), *G_sym);
        }
    }
}