tests.cc 14.9 KB
Newer Older
1 2 3
/* $Id: tests.cpp 431 2005-08-16 15:41:01Z kamenik $ */
/* Copyright 2005, Ondra Kamenik */

4
#include "GeneralMatrix.hh"
5
#include <dynlapack.h>
6
#include "SylvException.hh"
7

8 9
#include "rfs_tensor.hh"
#include "normal_moments.hh"
10

11 12 13 14 15
#include "vector_function.hh"
#include "quadrature.hh"
#include "smolyak.hh"
#include "product.hh"
#include "quasi_mcarlo.hh"
16

17 18
#include <iomanip>
#include <chrono>
19
#include <cmath>
20 21 22 23 24
#include <iostream>
#include <utility>
#include <array>
#include <memory>
#include <cstdlib>
25 26 27 28 29 30 31 32 33 34 35 36 37

// evaluates unfolded (Dx)^k power, where x is a vector, D is a
// Cholesky factor (lower triangular)
class MomentFunction : public VectorFunction
{
  GeneralMatrix D;
  int k;
public:
  MomentFunction(const GeneralMatrix &inD, int kk)
    : VectorFunction(inD.numRows(), UFSTensor::calcMaxOffset(inD.numRows(), kk)),
      D(inD), k(kk)
  {
  }
38 39 40
  MomentFunction(const MomentFunction &func) = default;

  std::unique_ptr<VectorFunction>
41
  clone() const override
42
  {
43
    return std::make_unique<MomentFunction>(*this);
44
  }
45
  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
46 47 48 49 50 51 52
};

void
MomentFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
  if (point.length() != indim() || out.length() != outdim())
    {
53 54
      std::cerr << "Wrong length of vectors in MomentFunction::eval" << std::endl;
      std::exit(EXIT_FAILURE);
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
    }
  Vector y(point);
  y.zeros();
  D.multaVec(y, point);
  URSingleTensor ypow(y, k);
  out.zeros();
  out.add(1.0, ypow.getData());
}

class TensorPower : public VectorFunction
{
  int k;
public:
  TensorPower(int nvar, int kk)
    : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk)
  {
  }
72 73 74
  TensorPower(const TensorPower &func) = default;

  std::unique_ptr<VectorFunction>
75
  clone() const override
76
  {
77
    return std::make_unique<TensorPower>(*this);
78
  }
79
  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
80 81 82 83 84 85 86
};

void
TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
  if (point.length() != indim() || out.length() != outdim())
    {
87 88
      std::cerr << "Wrong length of vectors in TensorPower::eval" << std::endl;
      std::exit(EXIT_FAILURE);
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
    }
  URSingleTensor ypow(point, k);
  out.zeros();
  out.add(1.0, ypow.getData());
}

// evaluates (1+1/d)^d*(x_1*...*x_d)^(1/d), its integral over <0,1>^d
// is 1.0, and its variation grows exponetially
// d = dim
class Function1 : public VectorFunction
{
  int dim;
public:
  Function1(int d)
    : VectorFunction(d, 1), dim(d)
  {
  }
  Function1(const Function1 &f)
    : VectorFunction(f.indim(), f.outdim()), dim(f.dim)
  {
  }
110
  std::unique_ptr<VectorFunction>
111
  clone() const override
112
  {
113
    return std::make_unique<Function1>(*this);
114
  }
115
  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
116 117 118 119 120 121 122
};

void
Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
  if (point.length() != dim || out.length() != 1)
    {
123 124
      std::cerr << "Wrong length of vectors in Function1::eval" << std::endl;
      std::exit(EXIT_FAILURE);
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
    }
  double r = 1;
  for (int i = 0; i < dim; i++)
    r *= point[i];
  r = pow(r, 1.0/dim);
  r *= pow(1.0 + 1.0/dim, (double) dim);
  out[0] = r;
}

// evaluates Function1 but with transformation x_i=0.5(y_i+1)
// this makes the new function integrate over <-1,1>^d to 1.0
class Function1Trans : public Function1
{
public:
  Function1Trans(int d)
    : Function1(d)
  {
  }
143 144 145
  Function1Trans(const Function1Trans &func) = default;

  std::unique_ptr<VectorFunction>
146
  clone() const override
147
  {
148
    return std::make_unique<Function1Trans>(*this);
149
  }
150
  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
};

void
Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
  Vector p(point.length());
  for (int i = 0; i < p.length(); i++)
    p[i] = 0.5*(point[i]+1);
  Function1::eval(p, sig, out);
  out.mult(pow(0.5, indim()));
}

// WallTimer class. Constructor saves the wall time, destructor
// cancels the current time from the saved, and prints the message
// with time information
class WallTimer
{
168
  std::string mes;
169
  std::chrono::time_point<std::chrono::high_resolution_clock> start;
170 171
  bool new_line;
public:
172
  WallTimer(std::string m, bool nl = true)
173
    : mes{m}, start{std::chrono::high_resolution_clock::now()}, new_line{nl}
174 175 176 177
  {
  }
  ~WallTimer()
  {
178 179 180
    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    std::cout << mes << std::setw(8) << std::setprecision(4) << duration.count();
181
    if (new_line)
182
      std::cout << std::endl;
183 184 185 186 187 188 189 190 191
  }
};

/****************************************************/
/*     declaration of TestRunnable class            */
/****************************************************/
class TestRunnable
{
public:
192
  const std::string name;
193 194
  int dim; // dimension of the solved problem
  int nvar; // number of variable of the solved problem
195
  TestRunnable(std::string name_arg, int d, int nv)
196
    : name{move(name_arg)}, dim(d), nvar(nv)
197 198
  {
  }
199
  virtual ~TestRunnable() = default;
200 201 202 203 204 205 206 207 208 209 210 211 212 213
  bool test() const;
  virtual bool run() const = 0;
protected:
  static bool smolyak_normal_moments(const GeneralMatrix &m, int imom, int level);
  static bool product_normal_moments(const GeneralMatrix &m, int imom, int level);
  static bool qmc_normal_moments(const GeneralMatrix &m, int imom, int level);
  static bool smolyak_product_cube(const VectorFunction &func, const Vector &res,
                                   double tol, int level);
  static bool qmc_cube(const VectorFunction &func, double res, double tol, int level);
};

bool
TestRunnable::test() const
{
214
  std::cout << "Running test <" << name << ">" << std::endl;
215 216 217 218 219 220 221
  bool passed;
  {
    WallTimer tim("Wall clock time ", false);
    passed = run();
  }
  if (passed)
    {
222
      std::cout << "............................ passed" << std::endl << std::endl;
223 224 225 226
      return passed;
    }
  else
    {
227
      std::cout << "............................ FAILED" << std::endl << std::endl;
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
      return passed;
    }
}

/****************************************************/
/*     definition of TestRunnable static methods    */
/****************************************************/
bool
TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level)
{
  // first make m*m' and then Cholesky factor
  GeneralMatrix mtr(m, "transpose");
  GeneralMatrix msq(m, mtr);

  // make vector function
  int dim = m.numRows();
  TensorPower tp(dim, imom);
  GaussConverterFunction func(tp, msq);

  // smolyak quadrature
  Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
  {
    WallTimer tim("\tSmolyak quadrature time:         ");
    GaussHermite gs;
    SmolyakQuadrature quad(dim, level, gs);
253
    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, smol_out);
254
    std::cout << "\tNumber of Smolyak evaluations:    " << quad.numEvals(level) << std::endl;
255 256 257 258
  }

  // check against theoretical moments
  UNormalMoments moments(imom, msq);
259
  smol_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
260
  std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
  return smol_out.getMax() < 1.e-7;
}

bool
TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level)
{
  // first make m*m' and then Cholesky factor
  GeneralMatrix mtr(m, "transpose");
  GeneralMatrix msq(m, mtr);

  // make vector function
  int dim = m.numRows();
  TensorPower tp(dim, imom);
  GaussConverterFunction func(tp, msq);

  // product quadrature
  Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
  {
    WallTimer tim("\tProduct quadrature time:         ");
    GaussHermite gs;
    ProductQuadrature quad(dim, gs);
282
    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, prod_out);
283
    std::cout << "\tNumber of product evaluations:    " << quad.numEvals(level) << std::endl;
284 285 286 287
  }

  // check against theoretical moments
  UNormalMoments moments(imom, msq);
288
  prod_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
289
  std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
290 291 292 293 294 295 296 297 298
  return prod_out.getMax() < 1.e-7;
}

bool
TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res,
                                   double tol, int level)
{
  if (res.length() != func.outdim())
    {
299 300
      std::cerr << "Incompatible dimensions of check value and function." << std::endl;
      std::exit(EXIT_FAILURE);
301 302 303 304 305 306 307 308 309
    }

  GaussLegendre glq;
  Vector out(func.outdim());
  double smol_error;
  double prod_error;
  {
    WallTimer tim("\tSmolyak quadrature time:         ");
    SmolyakQuadrature quad(func.indim(), level, glq);
310
    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
311 312
    out.add(-1.0, res);
    smol_error = out.getMax();
313 314
    std::cout << "\tNumber of Smolyak evaluations:    " << quad.numEvals(level) << std::endl;
    std::cout << "\tError:                            " << std::setw(16) << std::setprecision(12) << smol_error << std::endl;
315 316 317 318
  }
  {
    WallTimer tim("\tProduct quadrature time:         ");
    ProductQuadrature quad(func.indim(), glq);
319
    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
320 321
    out.add(-1.0, res);
    prod_error = out.getMax();
322 323
    std::cout << "\tNumber of product evaluations:    " << quad.numEvals(level) << std::endl;
    std::cout << "\tError:                            " << std::setw(16) << std::setprecision(12) << prod_error << std::endl;
324 325 326 327 328 329 330 331 332 333 334 335 336 337 338
  }

  return smol_error < tol && prod_error < tol;
}

bool
TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int level)
{
  Vector r(1);
  double error1;
  {
    WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time:  ");
    WarnockPerScheme wps;
    QMCarloCubeQuadrature qmc(func.indim(), level, wps);
    //		qmc.savePoints("warnock.txt", level);
339
    qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
340
    error1 = std::max(res - r[0], r[0] - res);
341
    std::cout << "\tQuasi-Monte Carlo (Warnock scrambling) error: " << std::setw(16) << std::setprecision(12) << error1 << std::endl;
342 343 344 345 346 347 348
  }
  double error2;
  {
    WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time:  ");
    ReversePerScheme rps;
    QMCarloCubeQuadrature qmc(func.indim(), level, rps);
    //		qmc.savePoints("reverse.txt", level);
349
    qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
350
    error2 = std::max(res - r[0], r[0] - res);
351
    std::cout << "\tQuasi-Monte Carlo (reverse scrambling) error: " << std::setw(16) << std::setprecision(12) << error2 << std::endl;
352 353 354 355 356 357 358
  }
  double error3;
  {
    WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time:       ");
    IdentityPerScheme ips;
    QMCarloCubeQuadrature qmc(func.indim(), level, ips);
    //		qmc.savePoints("identity.txt", level);
359
    qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
360
    error3 = std::max(res - r[0], r[0] - res);
361
    std::cout << "\tQuasi-Monte Carlo (no scrambling) error:      " << std::setw(16) << std::setprecision(12) << error3 << std::endl;
362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
  }

  return error1 < tol && error2 < tol && error3 < tol;
}

/****************************************************/
/*     definition of TestRunnable subclasses        */
/****************************************************/
class SmolyakNormalMom1 : public TestRunnable
{
public:
  SmolyakNormalMom1()
    : TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2)
  {
  }

  bool
379
  run() const override
380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
  {
    GeneralMatrix m(2, 2);
    m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
    return smolyak_normal_moments(m, 4, 4);
  }
};

class SmolyakNormalMom2 : public TestRunnable
{
public:
  SmolyakNormalMom2()
    : TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3)
  {
  }

  bool
396
  run() const override
397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414
  {
    GeneralMatrix m(3, 3);
    m.zeros();
    m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
    m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
    return smolyak_normal_moments(m, 8, 8);
  }
};

class ProductNormalMom1 : public TestRunnable
{
public:
  ProductNormalMom1()
    : TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2)
  {
  }

  bool
415
  run() const override
416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
  {
    GeneralMatrix m(2, 2);
    m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
    return product_normal_moments(m, 4, 4);
  }
};

class ProductNormalMom2 : public TestRunnable
{
public:
  ProductNormalMom2()
    : TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3)
  {
  }

  bool
432
  run() const override
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
  {
    GeneralMatrix m(3, 3);
    m.zeros();
    m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
    m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
    return product_normal_moments(m, 8, 8);
  }
};

// note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
class F1GaussLegendre : public TestRunnable
{
public:
  F1GaussLegendre()
    : TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1)
  {
  }

  bool
452
  run() const override
453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
  {
    Function1Trans f1(6);
    Vector res(1); res[0] = 1.0;
    return smolyak_product_cube(f1, res, 1e-2, 13);
  }
};

class F1QuasiMCarlo : public TestRunnable
{
public:
  F1QuasiMCarlo()
    : TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1)
  {
  }

  bool
469
  run() const override
470 471 472 473 474 475 476 477 478
  {
    Function1 f1(6);
    return qmc_cube(f1, 1.0, 1.e-4, 1000000);
  }
};

int
main()
{
479
  std::vector<std::unique_ptr<TestRunnable>> all_tests;
480
  // fill in vector of all tests
481 482 483 484
  all_tests.push_back(std::make_unique<SmolyakNormalMom1>());
  all_tests.push_back(std::make_unique<SmolyakNormalMom2>());
  all_tests.push_back(std::make_unique<ProductNormalMom1>());
  all_tests.push_back(std::make_unique<ProductNormalMom2>());
485 486 487
  all_tests.push_back(std::make_unique<F1GaussLegendre>());
  all_tests.push_back(std::make_unique<F1QuasiMCarlo>());

488 489 490
  // find maximum dimension and maximum nvar
  int dmax = 0;
  int nvmax = 0;
491
  for (const auto &test : all_tests)
492
    {
493 494 495 496
      if (dmax < test->dim)
        dmax = test->dim;
      if (nvmax < test->nvar)
        nvmax = test->nvar;
497 498 499 500 501
    }
  tls.init(dmax, nvmax); // initialize library

  // launch the tests
  int success = 0;
502
  for (const auto &test : all_tests)
503 504 505
    {
      try
        {
506
          if (test->test())
507 508 509 510
            success++;
        }
      catch (const TLException &e)
        {
511
          std::cout << "Caught TL exception in <" << test->name << ">:" << std::endl;
512 513 514 515
          e.print();
        }
      catch (SylvException &e)
        {
516
          std::cout << "Caught Sylv exception in <" << test->name << ">:" << std::endl;
517 518 519 520
          e.printMessage();
        }
    }

521 522
  int nfailed = all_tests.size() - success;
  std::cout << "There were " << nfailed << " tests that failed out of "
523
            << all_tests.size() << " tests run." << std::endl;
524

525 526 527 528
  if (nfailed)
    return EXIT_FAILURE;
  else
    return EXIT_SUCCESS;
529
}