dynare/dynare++/integ/src/quadrature-points.cc

// Copyright (C) 2008-2011, Ondra Kamenik

#include "parser/cc/matrix_parser.h"
#include "utils/cc/memory_file.h"
#include "utils/cc/exception.h"
#include "sylv/cc/GeneralMatrix.h"
#include "sylv/cc/Vector.h"
#include "sylv/cc/SymSchurDecomp.h"
#include "sylv/cc/SylvException.h"
#include "integ/cc/quadrature.hh"
#include "integ/cc/smolyak.hh"
#include "integ/cc/product.hh"

#include <getopt.h>
#include <cstdio>

#include <cmath>

struct QuadParams
{
  const char *outname;
  const char *vcovname;
  int max_level;
  double discard_weight;
  QuadParams(int argc, char **argv);
  void check_consistency() const;
private:
  enum {opt_max_level, opt_discard_weight, opt_vcov};
};

QuadParams::QuadParams(int argc, char **argv)
  : outname(NULL), vcovname(NULL), max_level(3), discard_weight(0.0)
{
  if (argc == 1)
    {
      // print the help and exit
      exit(1);
    }

  outname = argv[argc-1];
  argc--;

  struct option const opts [] = {
    {"max-level", required_argument, NULL, opt_max_level},
    {"discard-weight", required_argument, NULL, opt_discard_weight},
    {"vcov", required_argument, NULL, opt_vcov},
    {NULL, 0, NULL, 0}
  };

  int ret;
  int index;
  while (-1 != (ret = getopt_long(argc, argv, "", opts, &index)))
    {
      switch (ret)
        {
        case opt_max_level:
          if (1 != sscanf(optarg, "%d", &max_level))
            fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
          break;
        case opt_discard_weight:
          if (1 != sscanf(optarg, "%lf", &discard_weight))
            fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
          break;
        case opt_vcov:
          vcovname = optarg;
          break;
        }
    }

  check_consistency();
}

void
QuadParams::check_consistency() const
{
  if (outname == NULL)
    {
      fprintf(stderr, "Error: output name not set\n");
      exit(1);
    }

  if (vcovname == NULL)
    {
      fprintf(stderr, "Error: vcov file name not set\n");
      exit(1);
    }
}

/** Utility class for ordering pointers to vectors according their
 * ordering. */
struct OrderVec
{
  bool
  operator()(const Vector *a, const Vector *b) const
  {
    return *a < *b;
  }
};

int
main(int argc, char **argv)
{
  QuadParams params(argc, argv);

  // open output file for writing
  FILE *fout;
  if (NULL == (fout = fopen(params.outname, "w")))
    {
      fprintf(stderr, "Could not open %s for writing\n", params.outname);
      exit(1);
    }

  try
    {

      // open memory file for vcov
      ogu::MemoryFile vcov_mf(params.vcovname);

      // parse the vcov matrix
      ogp::MatrixParser mp;
      mp.parse(vcov_mf.length(), vcov_mf.base());
      if (mp.nrows() != mp.ncols())
        throw ogu::Exception(__FILE__, __LINE__,
                             "VCOV matrix not square");
      // and put to the GeneralMatrix
      GeneralMatrix vcov(mp.nrows(), mp.ncols());
      vcov.zeros();
      for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
        vcov.get(it.row(), it.col()) = *it;

      // calculate the factor A of vcov, so that A*A^T=VCOV
      GeneralMatrix A(vcov.numRows(), vcov.numRows());
      SymSchurDecomp ssd(vcov);
      ssd.getFactor(A);

      // construct Gauss-Hermite quadrature
      GaussHermite ghq;
      // construct Smolyak quadrature
      int level = params.max_level;
      SmolyakQuadrature sq(vcov.numRows(), level, ghq);

      printf("Dimension:                %d\n", vcov.numRows());
      printf("Maximum level:            %d\n", level);
      printf("Total number of nodes:    %d\n", sq.numEvals(level));

      // put the points to the vector
      std::vector<Vector *> points;
      for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
        points.push_back(new Vector((const Vector &) qit.point()));
      // sort and uniq
      OrderVec ordvec;
      std::sort(points.begin(), points.end(), ordvec);
      std::vector<Vector *>::iterator new_end = std::unique(points.begin(), points.end());
      for (std::vector<Vector *>::iterator it = new_end; it != points.end(); ++it)
        delete *it;
      points.erase(new_end, points.end());

      printf("Duplicit nodes removed:   %lu\n", (unsigned long) (sq.numEvals(level)-points.size()));

      // calculate weights and mass
      double mass = 0.0;
      std::vector<double> weights;
      for (int i = 0; i < (int) points.size(); i++)
        {
          weights.push_back(std::exp(-points[i]->dot(*(points[i]))));
          mass += weights.back();
        }

      // calculate discarded mass
      double discard_mass = 0.0;
      for (int i = 0; i < (int) weights.size(); i++)
        if (weights[i]/mass < params.discard_weight)
          discard_mass += weights[i];

      printf("Total mass discarded:     %f\n", discard_mass/mass);

      // dump the results
      int npoints = 0;
      double upscale_weight = 1/(mass-discard_mass);
      Vector x(vcov.numRows());
      for (int i = 0; i < (int) weights.size(); i++)
        if (weights[i]/mass >= params.discard_weight)
          {
            // print the upscaled weight
            fprintf(fout, "%20.16g", upscale_weight*weights[i]);
            // multiply point with the factor A and sqrt(2)
            A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
            // print the coordinates
            for (int j = 0; j < x.length(); j++)
              fprintf(fout, " %20.16g", x[j]);
            fprintf(fout, "\n");
            npoints++;
          }

      printf("Final number of points:   %d\n", npoints);

      fclose(fout);

    }
  catch (const SylvException &e)
    {
      e.printMessage();
      return 1;
    }
  catch (const ogu::Exception &e)
    {
      e.print();
      return 1;
    }

  return 0;
}