Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
MCMC-Laplace.h
Go to the documentation of this file.
1 
304  *
305  *
306  *
307  * #include <deal.II/base/timer.h>
308  * #include <deal.II/grid/tria.h>
309  * #include <deal.II/dofs/dof_handler.h>
310  * #include <deal.II/grid/grid_generator.h>
311  * #include <deal.II/grid/tria_accessor.h>
312  * #include <deal.II/grid/tria_iterator.h>
313  * #include <deal.II/dofs/dof_accessor.h>
314  * #include <deal.II/fe/fe_q.h>
315  * #include <deal.II/dofs/dof_tools.h>
316  * #include <deal.II/fe/fe_values.h>
317  * #include <deal.II/base/quadrature_lib.h>
318  * #include <deal.II/base/function.h>
319  * #include <deal.II/numerics/vector_tools.h>
320  * #include <deal.II/numerics/matrix_tools.h>
321  * #include <deal.II/lac/vector.h>
322  * #include <deal.II/lac/full_matrix.h>
323  * #include <deal.II/lac/sparse_matrix.h>
324  * #include <deal.II/lac/dynamic_sparsity_pattern.h>
325  * #include <deal.II/lac/solver_cg.h>
326  * #include <deal.II/lac/precondition.h>
327  * #include <deal.II/lac/sparse_direct.h>
328  *
329  * #include <deal.II/numerics/data_out.h>
330  *
331  * #include <fstream>
332  * #include <iostream>
333  * #include <random>
334  *
335  * #include <deal.II/base/logstream.h>
336  *
337  * using namespace dealii;
338  *
339  *
340  * @endcode
341  *
342  * The following is a namespace in which we define the solver of the PDE.
343  * The main class implements an abstract `Interface` class declared at
344  * the top, which provides for an `evaluate()` function that, given
345  * a coefficient vector, solves the PDE discussed in the Readme file
346  * and then evaluates the solution at the 169 mentioned points.
347  *
348 
349  *
350  * The solver follows the basic layout of @ref step_4 "step-4", though it precomputes
351  * a number of things in the `setup_system()` function, such as the
352  * evaluation of the matrix that corresponds to the point evaluations,
353  * as well as the local contributions to matrix and right hand side.
354  *
355 
356  *
357  * Rather than commenting on everything in detail, in the following
358  * we will only document those things that are not already clear from
359  * @ref step_4 "step-4" and a small number of other tutorial programs.
360  *
361  * @code
362  * namespace ForwardSimulator
363  * {
364  * class Interface
365  * {
366  * public:
367  * virtual Vector<double> evaluate(const Vector<double> &coefficients) = 0;
368  * };
369  *
370  *
371  *
372  * template <int dim>
373  * class PoissonSolver : public Interface
374  * {
375  * public:
376  * PoissonSolver(const unsigned int global_refinements,
377  * const unsigned int fe_degree,
378  * const std::string &dataset_name);
379  * virtual Vector<double>
380  * evaluate(const Vector<double> &coefficients) override;
381  *
382  * private:
383  * void make_grid(const unsigned int global_refinements);
384  * void setup_system();
385  * void assemble_system(const Vector<double> &coefficients);
386  * void solve();
387  * void output_results(const Vector<double> &coefficients) const;
388  *
390  * FE_Q<dim> fe;
391  * DoFHandler<dim> dof_handler;
392  *
394  * Vector<double> cell_rhs;
395  * std::map<types::global_dof_index,double> boundary_values;
396  *
397  * SparsityPattern sparsity_pattern;
398  * SparseMatrix<double> system_matrix;
399  *
400  * Vector<double> solution;
401  * Vector<double> system_rhs;
402  *
403  * std::vector<Point<dim>> measurement_points;
404  *
405  * SparsityPattern measurement_sparsity;
406  * SparseMatrix<double> measurement_matrix;
407  *
408  * TimerOutput timer;
409  * unsigned int nth_evaluation;
410  *
411  * const std::string &dataset_name;
412  * };
413  *
414  *
415  *
416  * template <int dim>
417  * PoissonSolver<dim>::PoissonSolver(const unsigned int global_refinements,
418  * const unsigned int fe_degree,
419  * const std::string &dataset_name)
420  * : fe(fe_degree)
421  * , dof_handler(triangulation)
422  * , timer(std::cout, TimerOutput::summary, TimerOutput::cpu_times)
423  * , nth_evaluation(0)
424  * , dataset_name(dataset_name)
425  * {
426  * make_grid(global_refinements);
427  * setup_system();
428  * }
429  *
430  *
431  *
432  * template <int dim>
433  * void PoissonSolver<dim>::make_grid(const unsigned int global_refinements)
434  * {
435  * Assert(global_refinements >= 3,
436  * ExcMessage("This program makes the assumption that the mesh for the "
437  * "solution of the PDE is at least as fine as the one used "
438  * "in the definition of the coefficient."));
440  * triangulation.refine_global(global_refinements);
441  *
442  * std::cout << " Number of active cells: " << triangulation.n_active_cells()
443  * << std::endl;
444  * }
445  *
446  *
447  *
448  * template <int dim>
449  * void PoissonSolver<dim>::setup_system()
450  * {
451  * @endcode
452  *
453  * First define the finite element space:
454  *
455  * @code
456  * dof_handler.distribute_dofs(fe);
457  *
458  * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
459  * << std::endl;
460  *
461  * @endcode
462  *
463  * Then set up the main data structures that will hold the discrete problem:
464  *
465  * @code
466  * {
467  * DynamicSparsityPattern dsp(dof_handler.n_dofs());
468  * DoFTools::make_sparsity_pattern(dof_handler, dsp);
469  * sparsity_pattern.copy_from(dsp);
470  *
471  * system_matrix.reinit(sparsity_pattern);
472  *
473  * solution.reinit(dof_handler.n_dofs());
474  * system_rhs.reinit(dof_handler.n_dofs());
475  * }
476  *
477  * @endcode
478  *
479  * And then define the tools to do point evaluation. We choose
480  * a set of 13x13 points evenly distributed across the domain:
481  *
482  * @code
483  * {
484  * const unsigned int n_points_per_direction = 13;
485  * const double dx = 1. / (n_points_per_direction + 1);
486  *
487  * for (unsigned int x = 1; x <= n_points_per_direction; ++x)
488  * for (unsigned int y = 1; y <= n_points_per_direction; ++y)
489  * measurement_points.emplace_back(x * dx, y * dx);
490  *
491  * @endcode
492  *
493  * First build a full matrix of the evaluation process. We do this
494  * even though the matrix is really sparse -- but we don't know
495  * which entries are nonzero. Later, the `copy_from()` function
496  * calls build a sparsity pattern and a sparse matrix from
497  * the dense matrix.
498  *
499  * @code
500  * Vector<double> weights(dof_handler.n_dofs());
501  * FullMatrix<double> full_measurement_matrix(n_points_per_direction *
502  * n_points_per_direction,
503  * dof_handler.n_dofs());
504  *
505  * for (unsigned int index = 0; index < measurement_points.size(); ++index)
506  * {
507  * VectorTools::create_point_source_vector(dof_handler,
508  * measurement_points[index],
509  * weights);
510  * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
511  * full_measurement_matrix(index, i) = weights(i);
512  * }
513  *
514  * measurement_sparsity.copy_from(full_measurement_matrix);
515  * measurement_matrix.reinit(measurement_sparsity);
516  * measurement_matrix.copy_from(full_measurement_matrix);
517  * }
518  *
519  * @endcode
520  *
521  * Next build the mapping from cell to the index in the 64-element
522  * coefficient vector:
523  *
524  * @code
525  * for (const auto &cell : triangulation.active_cell_iterators())
526  * {
527  * const unsigned int i = std::floor(cell->center()[0] * 8);
528  * const unsigned int j = std::floor(cell->center()[1] * 8);
529  *
530  * const unsigned int index = i + 8 * j;
531  *
532  * cell->set_user_index(index);
533  * }
534  *
535  * @endcode
536  *
537  * Finally prebuild the building blocks of the linear system as
538  * discussed in the Readme file :
539  *
540  * @code
541  * {
542  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
543  *
544  * cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
545  * cell_rhs.reinit(dofs_per_cell);
546  *
547  * const QGauss<dim> quadrature_formula(fe.degree+1);
548  * const unsigned int n_q_points = quadrature_formula.size();
549  *
550  * FEValues<dim> fe_values(fe,
551  * quadrature_formula,
552  * update_values | update_gradients |
553  * update_JxW_values);
554  *
555  * fe_values.reinit(dof_handler.begin_active());
556  *
557  * for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
558  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
559  * {
560  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
561  * cell_matrix(i, j) +=
562  * (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
563  * fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
564  * fe_values.JxW(q_index)); // dx
565  *
566  * cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
567  * 10.0 * // f(x_q)
568  * fe_values.JxW(q_index)); // dx
569  * }
570  *
571  * VectorTools::interpolate_boundary_values(dof_handler,
572  * 0,
573  * ZeroFunction<dim>(),
574  * boundary_values);
575  * }
576  * }
577  *
578  *
579  *
580  * @endcode
581  *
582  * Given that we have pre-built the matrix and right hand side contributions
583  * for a (representative) cell, the function that assembles the matrix is
584  * pretty short and straightforward:
585  *
586  * @code
587  * template <int dim>
588  * void PoissonSolver<dim>::assemble_system(const Vector<double> &coefficients)
589  * {
590  * Assert(coefficients.size() == 64, ExcInternalError());
591  *
592  * system_matrix = 0;
593  * system_rhs = 0;
594  *
595  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
596  *
597  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
598  *
599  * for (const auto &cell : dof_handler.active_cell_iterators())
600  * {
601  * const double coefficient = coefficients(cell->user_index());
602  *
603  * cell->get_dof_indices(local_dof_indices);
604  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
605  * {
606  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
607  * system_matrix.add(local_dof_indices[i],
608  * local_dof_indices[j],
609  * coefficient * cell_matrix(i, j));
610  *
611  * system_rhs(local_dof_indices[i]) += cell_rhs(i);
612  * }
613  * }
614  *
615  * MatrixTools::apply_boundary_values(boundary_values,
616  * system_matrix,
617  * solution,
618  * system_rhs);
619  * }
620  *
621  *
622  * @endcode
623  *
624  * The same is true for the function that solves the linear system:
625  *
626  * @code
627  * template <int dim>
628  * void PoissonSolver<dim>::solve()
629  * {
630  * SparseDirectUMFPACK solver;
631  * solver.factorize(system_matrix);
632  * solver.vmult(solution, system_rhs);
633  * }
634  *
635  *
636  *
637  * @endcode
638  *
639  * The following function outputs graphical data for the most recently
640  * used coefficient and corresponding solution of the PDE. Collecting
641  * the coefficient values requires translating from the 64-element
642  * coefficient vector and the cells that correspond to each of these
643  * elements. The rest remains pretty obvious, with the exception
644  * of including the number of the current sample into the file name.
645  *
646  * @code
647  * template <int dim>
648  * void
649  * PoissonSolver<dim>::output_results(const Vector<double> &coefficients) const
650  * {
651  * Vector<float> coefficient_values(triangulation.n_active_cells());
652  * for (const auto &cell : triangulation.active_cell_iterators())
653  * coefficient_values[cell->active_cell_index()] =
654  * coefficients(cell->user_index());
655  *
656  * DataOut<dim> data_out;
657  *
658  * data_out.attach_dof_handler(dof_handler);
659  * data_out.add_data_vector(solution, "solution");
660  * data_out.add_data_vector(coefficient_values, "coefficient");
661  *
662  * data_out.build_patches();
663  *
664  * std::ofstream output("solution-" +
665  * Utilities::int_to_string(nth_evaluation, 10) + ".vtu");
666  * data_out.write_vtu(output);
667  * }
668  *
669  *
670  *
671  * @endcode
672  *
673  * The following is the main function of this class: Given a coefficient
674  * vector, it assembles the linear system, solves it, and then evaluates
675  * the solution at the measurement points by applying the measurement
676  * matrix to the solution vector. That vector of "measured" values
677  * is then returned.
678  *
679 
680  *
681  * The function will also output the solution in a graphical format
682  * if you un-comment the corresponding statement in the third
683  * code block. However, you may end up with a very large amount
684  * of data: This code is producing, at the minimum, 10,000 samples
685  * and creating output for each one of them is surely more data
686  * than you ever want to see!
687  *
688 
689  *
690  * At the end of the function, we output some timing information
691  * every 10,000 samples.
692  *
693  * @code
694  * template <int dim>
695  * Vector<double>
696  * PoissonSolver<dim>::evaluate(const Vector<double> &coefficients)
697  * {
698  * {
699  * TimerOutput::Scope section(timer, "Building linear systems");
700  * assemble_system(coefficients);
701  * }
702  *
703  * {
704  * TimerOutput::Scope section(timer, "Solving linear systems");
705  * solve();
706  * }
707  *
708  * Vector<double> measurements(measurement_matrix.m());
709  * {
710  * TimerOutput::Scope section(timer, "Postprocessing");
711  *
712  * measurement_matrix.vmult(measurements, solution);
713  * Assert(measurements.size() == measurement_points.size(),
714  * ExcInternalError());
715  *
716  * /* output_results(coefficients); */
717  * }
718  *
719  * ++nth_evaluation;
720  * if (nth_evaluation % 10000 == 0)
721  * timer.print_summary();
722  *
723  * return std::move(measurements);
724  * }
725  * } // namespace ForwardSimulator
726  *
727  *
728  * @endcode
729  *
730  * The following namespaces define the statistical properties of the Bayesian
731  * inverse problem. The first is about the definition of the measurement
732  * statistics (the "likelihood"), which we here assume to be a normal
733  * distribution @f$N(\mu,\sigma I)@f$ with mean value @f$\mu@f$ given by the
734  * actual measurement vector (passed as an argument to the constructor
735  * of the `Gaussian` class and standard deviation @f$\sigma@f$.
736  *
737 
738  *
739  * For reasons of numerical accuracy, it is useful to not return the
740  * actual likelihood, but its logarithm. This is because these
741  * values can be very small, occasionally on the order of @f$e^{-100}@f$,
742  * for which it becomes very difficult to compute accurate
743  * values.
744  *
745  * @code
746  * namespace LogLikelihood
747  * {
748  * class Interface
749  * {
750  * public:
751  * virtual double log_likelihood(const Vector<double> &x) const = 0;
752  * };
753  *
754  *
755  * class Gaussian : public Interface
756  * {
757  * public:
758  * Gaussian(const Vector<double> &mu, const double sigma);
759  *
760  * virtual double log_likelihood(const Vector<double> &x) const override;
761  *
762  * private:
763  * const Vector<double> mu;
764  * const double sigma;
765  * };
766  *
767  * Gaussian::Gaussian(const Vector<double> &mu, const double sigma)
768  * : mu(mu)
769  * , sigma(sigma)
770  * {}
771  *
772  *
773  * double Gaussian::log_likelihood(const Vector<double> &x) const
774  * {
775  * Vector<double> x_minus_mu = x;
776  * x_minus_mu -= mu;
777  *
778  * return -x_minus_mu.norm_sqr() / (2 * sigma * sigma);
779  * }
780  * } // namespace LogLikelihood
781  *
782  *
783  * @endcode
784  *
785  * Next up is the "prior" imposed on the coefficients. We assume
786  * that the logarithms of the entries of the coefficient vector
787  * are all distributed as a Gaussian with given mean and standard
788  * deviation. If the logarithms of the coefficients are normally
789  * distributed, then this implies in particular that the coefficients
790  * can only be positive, which is a useful property to ensure the
791  * well-posedness of the forward problem.
792  *
793 
794  *
795  * For the same reasons as for the likelihood above, the interface
796  * for the prior asks for returning the *logarithm* of the prior,
797  * instead of the prior probability itself.
798  *
799  * @code
800  * namespace LogPrior
801  * {
802  * class Interface
803  * {
804  * public:
805  * virtual double log_prior(const Vector<double> &x) const = 0;
806  * };
807  *
808  *
809  * class LogGaussian : public Interface
810  * {
811  * public:
812  * LogGaussian(const double mu, const double sigma);
813  *
814  * virtual double log_prior(const Vector<double> &x) const override;
815  *
816  * private:
817  * const double mu;
818  * const double sigma;
819  * };
820  *
821  * LogGaussian::LogGaussian(const double mu, const double sigma)
822  * : mu(mu)
823  * , sigma(sigma)
824  * {}
825  *
826  *
827  * double LogGaussian::log_prior(const Vector<double> &x) const
828  * {
829  * double log_of_product = 0;
830  *
831  * for (const auto &el : x)
832  * log_of_product +=
833  * -(std::log(el) - mu) * (std::log(el) - mu) / (2 * sigma * sigma);
834  *
835  * return log_of_product;
836  * }
837  * } // namespace LogPrior
838  *
839  *
840  *
841  * @endcode
842  *
843  * The Metropolis-Hastings algorithm requires a method to create a new sample
844  * given a previous sample. We do this by perturbing the current (coefficient)
845  * sample randomly using a Gaussian distribution centered at the current
846  * sample. To ensure that the samples' individual entries all remain
847  * positive, we use a Gaussian distribution in logarithm space -- in other
848  * words, instead of *adding* a small perturbation with mean value zero,
849  * we *multiply* the entries of the current sample by a factor that
850  * is the exponential of a random number with mean zero. (Because the
851  * exponential of zero is one, this means that the most likely factors
852  * to multiply the existing sample entries by are close to one. And
853  * because the exponential of a number is always positive, we never
854  * get negative samples this way.)
855  *
856 
857  *
858  * But the Metropolis-Hastings sampler doesn't just need a perturbed
859  * sample @f$y@f$ location given the current sample location @f$x@f$. It also
860  * needs to know the ratio of the probability of reaching @f$y@f$ from
861  * @f$x@f$, divided by the probability of reaching @f$x@f$ from @f$y@f$. If we
862  * were to use a symmetric proposal distribution (e.g., a Gaussian
863  * distribution centered at @f$x@f$ with a width independent of @f$x@f$), then
864  * these two probabilities would be the same, and the ratio one. But
865  * that's not the case for the Gaussian in log space. It's not
866  * terribly difficult to verify that in that case, for a single
867  * component the ratio of these probabilities is @f$y_i/x_i@f$, and
868  * consequently for all components of the vector together, the
869  * probability is the product of these ratios.
870  *
871  * @code
872  * namespace ProposalGenerator
873  * {
874  * class Interface
875  * {
876  * public:
877  * virtual
878  * std::pair<Vector<double>,double>
879  * perturb(const Vector<double> &current_sample) const = 0;
880  * };
881  *
882  *
883  * class LogGaussian : public Interface
884  * {
885  * public:
886  * LogGaussian(const unsigned int random_seed, const double log_sigma);
887  *
888  * virtual
889  * std::pair<Vector<double>,double>
890  * perturb(const Vector<double> &current_sample) const;
891  *
892  * private:
893  * const double log_sigma;
894  * mutable std::mt19937 random_number_generator;
895  * };
896  *
897  *
898  *
899  * LogGaussian::LogGaussian(const unsigned int random_seed,
900  * const double log_sigma)
901  * : log_sigma(log_sigma)
902  * {
903  * random_number_generator.seed(random_seed);
904  * }
905  *
906  *
907  * std::pair<Vector<double>,double>
908  * LogGaussian::perturb(const Vector<double> &current_sample) const
909  * {
910  * Vector<double> new_sample = current_sample;
911  * double product_of_ratios = 1;
912  * for (auto &x : new_sample)
913  * {
914  * const double rnd = std::normal_distribution<>(0, log_sigma)(random_number_generator);
915  * const double exp_rnd = std::exp(rnd);
916  * x *= exp_rnd;
917  * product_of_ratios *= exp_rnd;
918  * }
919  *
920  * return {new_sample, product_of_ratios};
921  * }
922  *
923  * } // namespace ProposalGenerator
924  *
925  *
926  * @endcode
927  *
928  * The last main class is the Metropolis-Hastings sampler itself.
929  * If you understand the algorithm behind this method, then
930  * the following implementation should not be too difficult
931  * to read. The only thing of relevance is that descriptions
932  * of the algorithm typically ask whether the *ratio* of two
933  * probabilities (the "posterior" probabilities of the current
934  * and the previous samples, where the "posterior" is the product of the
935  * likelihood and the prior probability) is larger or smaller than a
936  * randomly drawn number. But because our interfaces return the
937  * *logarithms* of these probabilities, we now need to take
938  * the ratio of appropriate exponentials -- which is made numerically
939  * more stable by considering the exponential of the difference of
940  * the log probabilities. The only other slight complication is that
941  * we need to multiply this ratio by the ratio of proposal probabilities
942  * since we use a non-symmetric proposal distribution.
943  *
944 
945  *
946  * Finally, we note that the output is generated with 7 digits of
947  * accuracy. (The C++ default is 6 digits.) We do this because,
948  * as shown in the paper, we can determine the mean value of the
949  * probability distribution we are sampling here to at least six
950  * digits of accuracy, and do not want to be limited by the precision
951  * of the output.
952  *
953  * @code
954  * namespace Sampler
955  * {
956  * class MetropolisHastings
957  * {
958  * public:
959  * MetropolisHastings(ForwardSimulator::Interface & simulator,
960  * const LogLikelihood::Interface & likelihood,
961  * const LogPrior::Interface & prior,
962  * const ProposalGenerator::Interface &proposal_generator,
963  * const unsigned int random_seed,
964  * const std::string & dataset_name);
965  *
966  * void sample(const Vector<double> &starting_guess,
967  * const unsigned int n_samples);
968  *
969  * private:
970  * ForwardSimulator::Interface & simulator;
971  * const LogLikelihood::Interface & likelihood;
972  * const LogPrior::Interface & prior;
973  * const ProposalGenerator::Interface &proposal_generator;
974  *
975  * std::mt19937 random_number_generator;
976  *
977  * unsigned int sample_number;
978  * unsigned int accepted_sample_number;
979  *
980  * std::ofstream output_file;
981  *
982  * void write_sample(const Vector<double> &current_sample,
983  * const double current_log_likelihood);
984  * };
985  *
986  *
987  * MetropolisHastings::MetropolisHastings(
988  * ForwardSimulator::Interface & simulator,
989  * const LogLikelihood::Interface & likelihood,
990  * const LogPrior::Interface & prior,
991  * const ProposalGenerator::Interface &proposal_generator,
992  * const unsigned int random_seed,
993  * const std::string & dataset_name)
994  * : simulator(simulator)
995  * , likelihood(likelihood)
996  * , prior(prior)
997  * , proposal_generator(proposal_generator)
998  * , sample_number(0)
999  * , accepted_sample_number(0)
1000  * {
1001  * output_file.open("samples-" + dataset_name + ".txt");
1002  * output_file.precision(7);
1003  *
1004  * random_number_generator.seed(random_seed);
1005  * }
1006  *
1007  *
1008  * void MetropolisHastings::sample(const Vector<double> &starting_guess,
1009  * const unsigned int n_samples)
1010  * {
1011  * std::uniform_real_distribution<> uniform_distribution(0, 1);
1012  *
1013  * Vector<double> current_sample = starting_guess;
1014  * double current_log_posterior =
1015  * (likelihood.log_likelihood(simulator.evaluate(current_sample)) +
1016  * prior.log_prior(current_sample));
1017  *
1018  * ++sample_number;
1019  * ++accepted_sample_number;
1020  * write_sample(current_sample, current_log_posterior);
1021  *
1022  * for (unsigned int k = 1; k < n_samples; ++k, ++sample_number)
1023  * {
1024  * std::pair<Vector<double>,double>
1025  * perturbation = proposal_generator.perturb(current_sample);
1026  * const Vector<double> trial_sample = std::move (perturbation.first);
1027  * const double perturbation_probability_ratio = perturbation.second;
1028  *
1029  * const double trial_log_posterior =
1030  * (likelihood.log_likelihood(simulator.evaluate(trial_sample)) +
1031  * prior.log_prior(trial_sample));
1032  *
1033  * if (std::exp(trial_log_posterior - current_log_posterior) * perturbation_probability_ratio
1034  * >=
1035  * uniform_distribution(random_number_generator))
1036  * {
1037  * current_sample = trial_sample;
1038  * current_log_posterior = trial_log_posterior;
1039  *
1040  * ++accepted_sample_number;
1041  * }
1042  *
1043  * write_sample(current_sample, current_log_posterior);
1044  * }
1045  * }
1046  *
1047  *
1048  *
1049  * void MetropolisHastings::write_sample(const Vector<double> &current_sample,
1050  * const double current_log_posterior)
1051  * {
1052  * output_file << current_log_posterior << '\t';
1053  * output_file << accepted_sample_number << '\t';
1054  * for (const auto &x : current_sample)
1055  * output_file << x << ' ';
1056  * output_file << '\n';
1057  *
1058  * output_file.flush();
1059  * }
1060  * } // namespace Sampler
1061  *
1062  *
1063  * @endcode
1064  *
1065  * The final function is `main()`, which simply puts all of these pieces
1066  * together into one. The "exact solution", i.e., the "measurement values"
1067  * we use for this program are tabulated to make it easier for other
1068  * people to use in their own implementations of this benchmark. These
1069  * values created using the same main class above, but using 8 mesh
1070  * refinements and using a Q3 element -- i.e., using a much more accurate
1071  * method than the one we use in the forward simulator for generating
1072  * samples below (which uses 5 global mesh refinement steps and a Q1
1073  * element). If you wanted to regenerate this set of numbers, then
1074  * the following code snippet would do that:
1075  * <div class=CodeFragmentInTutorialComment>
1076  * @code
1077  * /* Set the exact coefficient: */
1078  * Vector<double> exact_coefficients(64);
1079  * for (auto &el : exact_coefficients)
1080  * el = 1.;
1081  * exact_coefficients(9) = exact_coefficients(10) = exact_coefficients(17) =
1082  * exact_coefficients(18) = 0.1;
1083  * exact_coefficients(45) = exact_coefficients(46) = exact_coefficients(53) =
1084  * exact_coefficients(54) = 10.;
1085  *
1086 
1087  * /* Compute the "correct" solution vector: */
1088  * const Vector<double> exact_solution =
1089  * ForwardSimulator::PoissonSolver<2>(/* global_refinements = */ 8,
1090  * /* fe_degree = */ 3,
1091  * /* prefix = */ "exact")
1092  * .evaluate(exact_coefficients);
1093  * @endcode
1094  * </div>
1095  *
1096  * @code
1097  * int main()
1098  * {
1099  * const bool testing = true;
1100  *
1101  * @endcode
1102  *
1103  * Run with one thread, so as to not step on other processes
1104  * doing the same at the same time. It turns out that the problem
1105  * is also so small that running with more than one thread
1106  * *increases* the runtime.
1107  *
1108  * @code
1109  * MultithreadInfo::set_thread_limit(1);
1110  *
1111  * const unsigned int random_seed = (testing ? 1U : std::random_device()());
1112  * const std::string dataset_name = Utilities::to_string(random_seed, 10);
1113  *
1114  * const Vector<double> exact_solution(
1115  * { 0.06076511762259369, 0.09601910120848481,
1116  * 0.1238852517838584, 0.1495184117375201,
1117  * 0.1841596127549784, 0.2174525028261122,
1118  * 0.2250996160898698, 0.2197954769002993,
1119  * 0.2074695698370926, 0.1889996477663016,
1120  * 0.1632722532153726, 0.1276782480038186,
1121  * 0.07711845915789312, 0.09601910120848552,
1122  * 0.2000589533367983, 0.3385592591951766,
1123  * 0.3934300024647806, 0.4040223892461541,
1124  * 0.4122329537843092, 0.4100480091545554,
1125  * 0.3949151637189968, 0.3697873264791232,
1126  * 0.33401826235924, 0.2850397806663382,
1127  * 0.2184260032478671, 0.1271121156350957,
1128  * 0.1238852517838611, 0.3385592591951819,
1129  * 0.7119285162766475, 0.8175712861756428,
1130  * 0.6836254116578105, 0.5779452419831157,
1131  * 0.5555615956136897, 0.5285181561736719,
1132  * 0.491439702849224, 0.4409367494853282,
1133  * 0.3730060082060772, 0.2821694983395214,
1134  * 0.1610176733857739, 0.1495184117375257,
1135  * 0.3934300024647929, 0.8175712861756562,
1136  * 0.9439154625527653, 0.8015904115095128,
1137  * 0.6859683749254024, 0.6561235366960599,
1138  * 0.6213197201867315, 0.5753611315000049,
1139  * 0.5140091754526823, 0.4325325506354165,
1140  * 0.3248315148915482, 0.1834600412730086,
1141  * 0.1841596127549917, 0.4040223892461832,
1142  * 0.6836254116578439, 0.8015904115095396,
1143  * 0.7870119561144977, 0.7373108331395808,
1144  * 0.7116558878070463, 0.6745179049094283,
1145  * 0.6235300574156917, 0.5559332704045935,
1146  * 0.4670304994474178, 0.3499809143811,
1147  * 0.19688263746294, 0.2174525028261253,
1148  * 0.4122329537843404, 0.5779452419831566,
1149  * 0.6859683749254372, 0.7373108331396063,
1150  * 0.7458811983178246, 0.7278968022406559,
1151  * 0.6904793535357751, 0.6369176452710288,
1152  * 0.5677443693743215, 0.4784738764865867,
1153  * 0.3602190632823262, 0.2031792054737325,
1154  * 0.2250996160898818, 0.4100480091545787,
1155  * 0.5555615956137137, 0.6561235366960938,
1156  * 0.7116558878070715, 0.727896802240657,
1157  * 0.7121928678670187, 0.6712187391428729,
1158  * 0.6139157775591492, 0.5478251665295381,
1159  * 0.4677122687599031, 0.3587654911000848,
1160  * 0.2050734291675918, 0.2197954769003094,
1161  * 0.3949151637190157, 0.5285181561736911,
1162  * 0.6213197201867471, 0.6745179049094407,
1163  * 0.690479353535786, 0.6712187391428787,
1164  * 0.6178408289359514, 0.5453605027237883,
1165  * 0.489575966490909, 0.4341716881061278,
1166  * 0.3534389974779456, 0.2083227496961347,
1167  * 0.207469569837099, 0.3697873264791366,
1168  * 0.4914397028492412, 0.5753611315000203,
1169  * 0.6235300574157017, 0.6369176452710497,
1170  * 0.6139157775591579, 0.5453605027237935,
1171  * 0.4336604929612851, 0.4109641743019312,
1172  * 0.3881864790111245, 0.3642640090182592,
1173  * 0.2179599909280145, 0.1889996477663011,
1174  * 0.3340182623592461, 0.4409367494853381,
1175  * 0.5140091754526943, 0.5559332704045969,
1176  * 0.5677443693743304, 0.5478251665295453,
1177  * 0.4895759664908982, 0.4109641743019171,
1178  * 0.395727260284338, 0.3778949322004734,
1179  * 0.3596268271857124, 0.2191250268948948,
1180  * 0.1632722532153683, 0.2850397806663325,
1181  * 0.373006008206081, 0.4325325506354207,
1182  * 0.4670304994474315, 0.4784738764866023,
1183  * 0.4677122687599041, 0.4341716881061055,
1184  * 0.388186479011099, 0.3778949322004602,
1185  * 0.3633362567187364, 0.3464457261905399,
1186  * 0.2096362321365655, 0.1276782480038148,
1187  * 0.2184260032478634, 0.2821694983395252,
1188  * 0.3248315148915535, 0.3499809143811097,
1189  * 0.3602190632823333, 0.3587654911000799,
1190  * 0.3534389974779268, 0.3642640090182283,
1191  * 0.35962682718569, 0.3464457261905295,
1192  * 0.3260728953424643, 0.180670595355394,
1193  * 0.07711845915789244, 0.1271121156350963,
1194  * 0.1610176733857757, 0.1834600412730144,
1195  * 0.1968826374629443, 0.2031792054737354,
1196  * 0.2050734291675885, 0.2083227496961245,
1197  * 0.2179599909279998, 0.2191250268948822,
1198  * 0.2096362321365551, 0.1806705953553887,
1199  * 0.1067965550010013 });
1200  *
1201  * @endcode
1202  *
1203  * Now run the forward simulator for samples:
1204  *
1205  * @code
1206  * ForwardSimulator::PoissonSolver<2> laplace_problem(
1207  * /* global_refinements = */ 5,
1208  * /* fe_degree = */ 1,
1209  * dataset_name);
1210  * LogLikelihood::Gaussian log_likelihood(exact_solution, 0.05);
1211  * LogPrior::LogGaussian log_prior(0, 2);
1212  * ProposalGenerator::LogGaussian proposal_generator(
1213  * random_seed, 0.09); /* so that the acceptance ratio is ~0.24 */
1214  * Sampler::MetropolisHastings sampler(laplace_problem,
1215  * log_likelihood,
1216  * log_prior,
1217  * proposal_generator,
1218  * random_seed,
1219  * dataset_name);
1220  *
1221  * Vector<double> starting_coefficients(64);
1222  * for (auto &el : starting_coefficients)
1223  * el = 1.;
1224  * sampler.sample(starting_coefficients,
1225  * (testing ? 250 * 40 /* takes 40 seconds */
1226  * :
1227  * 100000000 /* takes 6 days */
1228  * ));
1229  * }
1230  * @endcode
1231 
1232 
1233 */
DoFRenumbering::random
void random(DoFHandlerType &dof_handler)
Definition: dof_renumbering.cc:2102
FE_Q
Definition: fe_q.h:554
dealii
Definition: namespace_dealii.h:25
std::log
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5390
TimerOutput::cpu_times
@ cpu_times
Definition: timer.h:645
Triangulation< dim >
DataOutBase::dx
@ dx
Definition: data_out_base.h:1562
SparseMatrix< double >
TimerOutput::summary
@ summary
Definition: timer.h:605
DoFTools::always
@ always
Definition: dof_tools.h:236
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)
Definition: advection.h:80
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
TimerOutput
Definition: timer.h:546
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
VectorTools::mean
@ mean
Definition: vector_tools_common.h:79
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
SparsityPattern
Definition: sparsity_pattern.h:865
value
static const bool value
Definition: dof_tools_constraints.cc:433
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
FullMatrix< double >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
Vector< double >