smolyak.cc 5.93 KB
Newer Older
1 2 3
// Copyright 2005, Ondra Kamenik

#include "smolyak.hh"
4
#include "symmetry.hh"
5

6 7
#include <iostream>
#include <iomanip>
8 9 10 11 12 13

/* This constructs a beginning of |isum| summand in |smolq|. We must be
   careful here, since |isum| can be past-the-end, so no reference to
   vectors in |smolq| by |isum| must be done in this case. */

smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum)
14
  : smolq(q), isummand(isum),
15
    jseq(q.dimen(), 0),
16 17
    sig{q.dimen()},
    p{q.dimen()}
18 19
{
  if (isummand < q.numSummands())
20
    setPointAndWeight();
21 22 23 24 25
}

bool
smolpit::operator==(const smolpit &spit) const
{
26
  return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq;
27 28 29 30 31 32 33 34 35 36
}

/* We first try to increase index within the current summand. If we are
   at maximum, we go to a subsequent summand. Note that in this case all
   indices in |jseq| will be zero, so no change is needed. */

smolpit &
smolpit::operator++()
{
  // todo: throw if |smolq==NULL| or |jseq==NULL| or |sig==NULL|
37 38 39 40
  const IntSequence &levpts = smolq.levpoints[isummand];
  int i = smolq.dimen()-1;
  jseq[i]++;
  while (i >= 0 && jseq[i] == levpts[i])
41
    {
42
      jseq[i] = 0;
43 44
      i--;
      if (i >= 0)
45
        jseq[i]++;
46
    }
47
  sig.signalAfter(std::max(i, 0));
48 49 50 51

  if (i < 0)
    isummand++;

52
  if (isummand < smolq.numSummands())
53 54 55 56 57 58 59 60 61 62 63 64
    setPointAndWeight();

  return *this;
}

/* Here we set the point coordinates according to |jseq| and
   |isummand|. Also the weight is set here. */

void
smolpit::setPointAndWeight()
{
  // todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or
65 66 67 68
  // |p==NULL| or |isummand>=smolq.numSummands()|
  int l = smolq.level;
  int d = smolq.dimen();
  int sumk = (smolq.levels[isummand]).sum();
69 70
  int m1exp = l + d - sumk - 1;
  w = (2*(m1exp/2) == m1exp) ? 1.0 : -1.0;
71
  w *= smolq.psc.noverk(d-1, sumk-l);
72 73
  for (int i = 0; i < d; i++)
    {
74 75 76
      int ki = (smolq.levels[isummand])[i];
      p[i] = (smolq.uquad).point(ki, jseq[i]);
      w *= (smolq.uquad).weight(ki, jseq[i]);
77 78 79 80 81 82 83
    }
}

/* Debug print. */
void
smolpit::print() const
{
84 85 86 87 88 89 90 91 92 93 94 95 96
  auto ff = std::cout.flags();
  std::cout << "isum=" << std::left << std::setw(3) << isummand << std::right << ": [";
  for (int i = 0; i < smolq.dimen(); i++)
    std::cout << std::setw(2) << (smolq.levels[isummand])[i] << ' ';
  std::cout << "] j=[";
  for (int i = 0; i < smolq.dimen(); i++)
    std::cout << std::setw(2) << jseq[i] << ' ';
  std::cout << std::showpos << std::fixed << std::setprecision(3)
            << "] " << std::setw(4) << w << "*(";
  for (int i = 0; i < smolq.dimen()-1; i++)
    std::cout << std::setw(4) << p[i] << ' ';
  std::cout << std::setw(4) << p[smolq.dimen()-1] << ')' << std::endl;
  std::cout.flags(ff);
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 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 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
}

/* Here is the constructor of |SmolyakQuadrature|. We have to setup
   |levels|, |levpoints| and |cumevals|. We have to go through all
   $d$-dimensional sequences $k$, such that $l\leq \vert k\vert\leq
   l+d-1$ and all $k_i$ are positive integers. This is equivalent to
   going through all $k$ such that $l-d\leq\vert k\vert\leq l-1$ and all
   $k_i$ are non-negative integers. This is equivalent to going through
   $d+1$ dimensional sequences $(k,x)$ such that $\vert(k,x)\vert =l-1$
   and $x=0,\ldots,d-1$. The resulting sequence of positive integers is
   obtained by adding $1$ to all $k_i$. */

SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
  : QuadratureImpl<smolpit>(d), level(l), uquad(uq), psc(d-1, d-1)
{
  // 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)
    {
      if ((*si)[d] <= d-1)
        {
          IntSequence lev((const IntSequence &)*si, 0, d);
          lev.add(1);
          levels.push_back(lev);
          IntSequence levpts(d);
          for (int i = 0; i < d; i++)
            levpts[i] = uquad.numPoints(lev[i]);
          levpoints.push_back(levpts);
          cum += levpts.mult();
          cumevals.push_back(cum);
        }
    }
}

/* Here we return a number of evalutions of the quadrature for the
   given level. If the given level is the current one, we simply return
   the maximum cumulative number of evaluations. Otherwise we call costly
   |calcNumEvaluations| method. */

int
SmolyakQuadrature::numEvals(int l) const
{
  if (l != level)
    return calcNumEvaluations(l);
  else
    return cumevals[numSummands()-1];
}

/* This divides all the evaluations to |tn| approximately equal groups,
   and returns the beginning of the specified group |ti|. The granularity
   of divisions are summands as listed by |levels|. */

smolpit
SmolyakQuadrature::begin(int ti, int tn, int l) const
{
  // todo: raise is |level!=l|
  if (ti == tn)
    return smolpit(*this, numSummands());

  int totevals = cumevals[numSummands()-1];
  int evals = (totevals*ti)/tn;
  unsigned int isum = 0;
  while (isum+1 < numSummands() && cumevals[isum+1] < evals)
    isum++;
  return smolpit(*this, isum);
}

/* This is the same in a structure as |@<|SmolyakQuadrature| constructor@>|.
   We have to go through all summands and calculate
   a number of evaluations in each summand. */

int
SmolyakQuadrature::calcNumEvaluations(int lev) const
{
  int cum = 0;
  SymmetrySet ss(lev-1, dim+1);
  for (symiterator si(ss); !si.isEnd(); ++si)
    {
      if ((*si)[dim] <= dim-1)
        {
          IntSequence lev((const IntSequence &)*si, 0, dim);
          lev.add(1);
          IntSequence levpts(dim);
          for (int i = 0; i < dim; i++)
            levpts[i] = uquad.numPoints(lev[i]);
          cum += levpts.mult();
        }
    }
  return cum;
}

/* This returns a maximum level such that the number of evaluations is
   less than the given number. */

void
SmolyakQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
{
  int last_evals;
  evals = 1;
  lev = 1;
  do
    {
      lev++;
      last_evals = evals;
      evals = calcNumEvaluations(lev);
    }
  while (lev < uquad.numLevels() && evals <= max_evals);
  lev--;
  evals = last_evals;
}