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>
329 * #include <deal.II/numerics/data_out.h>
332 * #include <iostream>
335 * #include <deal.II/base/logstream.h>
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.
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.
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.
362 *
namespace ForwardSimulator
373 *
class PoissonSolver :
public Interface
376 * PoissonSolver(
const unsigned int global_refinements,
377 *
const unsigned int fe_degree,
378 *
const std::string &dataset_name);
383 *
void make_grid(
const unsigned int global_refinements);
384 *
void setup_system();
395 * std::map<types::global_dof_index,double> boundary_values;
403 * std::vector<Point<dim>> measurement_points;
409 *
unsigned int nth_evaluation;
411 *
const std::string &dataset_name;
417 * PoissonSolver<dim>::PoissonSolver(
const unsigned int global_refinements,
418 *
const unsigned int fe_degree,
419 *
const std::string &dataset_name)
423 * , nth_evaluation(0)
424 * , dataset_name(dataset_name)
426 * make_grid(global_refinements);
433 *
void PoissonSolver<dim>::make_grid(
const unsigned int global_refinements)
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."));
442 * std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
449 *
void PoissonSolver<dim>::setup_system()
453 * First define the finite element space:
456 * dof_handler.distribute_dofs(fe);
458 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
463 * Then
set up the main data structures that will hold the discrete problem:
469 * sparsity_pattern.copy_from(dsp);
471 * system_matrix.reinit(sparsity_pattern);
473 * solution.reinit(dof_handler.n_dofs());
474 * system_rhs.reinit(dof_handler.n_dofs());
479 * And then define the tools to
do point evaluation. We choose
480 * a
set of 13x13 points evenly distributed across the domain:
484 *
const unsigned int n_points_per_direction = 13;
485 *
const double dx = 1. / (n_points_per_direction + 1);
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);
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
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());
505 * for (unsigned int index = 0; index < measurement_points.size(); ++index)
507 * VectorTools::create_point_source_vector(dof_handler,
508 * measurement_points[index],
510 * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
511 * full_measurement_matrix(index, i) = weights(i);
514 * measurement_sparsity.copy_from(full_measurement_matrix);
515 * measurement_matrix.reinit(measurement_sparsity);
516 * measurement_matrix.copy_from(full_measurement_matrix);
521 * Next build the mapping from cell to the index in the 64-element
522 * coefficient vector:
525 * for (const auto &cell : triangulation.active_cell_iterators())
527 * const unsigned int i = std::floor(cell->center()[0] * 8);
528 * const unsigned int j = std::floor(cell->center()[1] * 8);
530 * const unsigned int index = i + 8 * j;
532 * cell->set_user_index(index);
537 * Finally prebuild the building blocks of the linear system as
538 * discussed in the Readme file :
542 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
544 * cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
545 * cell_rhs.reinit(dofs_per_cell);
547 * const QGauss<dim> quadrature_formula(fe.degree+1);
548 * const unsigned int n_q_points = quadrature_formula.size();
550 * FEValues<dim> fe_values(fe,
551 * quadrature_formula,
552 * update_values | update_gradients |
553 * update_JxW_values);
555 * fe_values.reinit(dof_handler.begin_active());
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)
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
566 * cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
568 * fe_values.JxW(q_index)); // dx
571 * VectorTools::interpolate_boundary_values(dof_handler,
573 * ZeroFunction<dim>(),
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:
588 * void PoissonSolver<dim>::assemble_system(const Vector<double> &coefficients)
590 * Assert(coefficients.size() == 64, ExcInternalError());
595 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
597 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
599 * for (const auto &cell : dof_handler.active_cell_iterators())
601 * const double coefficient = coefficients(cell->user_index());
603 * cell->get_dof_indices(local_dof_indices);
604 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
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));
611 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
615 * MatrixTools::apply_boundary_values(boundary_values,
624 * The same is true for the function that solves the linear system:
628 * void PoissonSolver<dim>::solve()
630 * SparseDirectUMFPACK solver;
631 * solver.factorize(system_matrix);
632 * solver.vmult(solution, system_rhs);
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.
649 * PoissonSolver<dim>::output_results(const Vector<double> &coefficients) const
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());
656 * DataOut<dim> data_out;
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");
662 * data_out.build_patches();
664 * std::ofstream output("solution-" +
665 * Utilities::int_to_string(nth_evaluation, 10) + ".vtu");
666 * data_out.write_vtu(output);
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
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!
690 * At the end of the function, we output some timing information
691 * every 10,000 samples.
696 * PoissonSolver<dim>::evaluate(const Vector<double> &coefficients)
699 * TimerOutput::Scope section(timer, "Building linear systems");
700 * assemble_system(coefficients);
704 * TimerOutput::Scope section(timer, "Solving linear systems");
708 * Vector<double> measurements(measurement_matrix.m());
710 * TimerOutput::Scope section(timer, "Postprocessing");
712 * measurement_matrix.vmult(measurements, solution);
713 * Assert(measurements.size() == measurement_points.size(),
714 * ExcInternalError());
716 * /* output_results(coefficients); */
720 * if (nth_evaluation % 10000 == 0)
721 * timer.print_summary();
723 * return std::move(measurements);
725 * } // namespace ForwardSimulator
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$.
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
746 * namespace LogLikelihood
751 * virtual double log_likelihood(const Vector<double> &x) const = 0;
755 * class Gaussian : public Interface
758 * Gaussian(const Vector<double> &mu, const double sigma);
760 * virtual double log_likelihood(const Vector<double> &x) const override;
763 * const Vector<double> mu;
764 * const double sigma;
767 * Gaussian::Gaussian(const Vector<double> &mu, const double sigma)
773 * double Gaussian::log_likelihood(const Vector<double> &x) const
775 * Vector<double> x_minus_mu = x;
778 * return -x_minus_mu.norm_sqr() / (2 * sigma * sigma);
780 * } // namespace LogLikelihood
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.
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.
805 * virtual double log_prior(const Vector<double> &x) const = 0;
809 * class LogGaussian : public Interface
812 * LogGaussian(const double mu, const double sigma);
814 * virtual double log_prior(const Vector<double> &x) const override;
818 * const double sigma;
821 * LogGaussian::LogGaussian(const double mu, const double sigma)
827 * double LogGaussian::log_prior(const Vector<double> &x) const
829 * double log_of_product = 0;
831 * for (const auto &el : x)
833 * -(std::log(el) - mu) * (std::log(el) - mu) / (2 * sigma * sigma);
835 * return log_of_product;
837 * } // namespace LogPrior
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.)
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.
872 * namespace ProposalGenerator
878 * std::pair<Vector<double>,double>
879 * perturb(const Vector<double> ¤t_sample) const = 0;
883 * class LogGaussian : public Interface
886 * LogGaussian(const unsigned int random_seed, const double log_sigma);
889 * std::pair<Vector<double>,double>
890 * perturb(const Vector<double> ¤t_sample) const;
893 * const double log_sigma;
894 * mutable std::mt19937 random_number_generator;
899 * LogGaussian::LogGaussian(const unsigned int random_seed,
900 * const double log_sigma)
901 * : log_sigma(log_sigma)
903 * random_number_generator.seed(random_seed);
907 * std::pair<Vector<double>,double>
908 * LogGaussian::perturb(const Vector<double> ¤t_sample) const
910 * Vector<double> new_sample = current_sample;
911 * double product_of_ratios = 1;
912 * for (auto &x : new_sample)
914 * const double rnd = std::normal_distribution<>(0, log_sigma)(random_number_generator);
915 * const double exp_rnd = std::exp(rnd);
917 * product_of_ratios *= exp_rnd;
920 * return {new_sample, product_of_ratios};
923 * } // namespace ProposalGenerator
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.
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
956 * class MetropolisHastings
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);
966 * void sample(const Vector<double> &starting_guess,
967 * const unsigned int n_samples);
970 * ForwardSimulator::Interface & simulator;
971 * const LogLikelihood::Interface & likelihood;
972 * const LogPrior::Interface & prior;
973 * const ProposalGenerator::Interface &proposal_generator;
975 * std::mt19937 random_number_generator;
977 * unsigned int sample_number;
978 * unsigned int accepted_sample_number;
980 * std::ofstream output_file;
982 * void write_sample(const Vector<double> ¤t_sample,
983 * const double current_log_likelihood);
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)
997 * , proposal_generator(proposal_generator)
999 * , accepted_sample_number(0)
1001 * output_file.open("samples-" + dataset_name + ".txt");
1002 * output_file.precision(7);
1004 * random_number_generator.seed(random_seed);
1008 * void MetropolisHastings::sample(const Vector<double> &starting_guess,
1009 * const unsigned int n_samples)
1011 * std::uniform_real_distribution<> uniform_distribution(0, 1);
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));
1019 * ++accepted_sample_number;
1020 * write_sample(current_sample, current_log_posterior);
1022 * for (unsigned int k = 1; k < n_samples; ++k, ++sample_number)
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;
1029 * const double trial_log_posterior =
1030 * (likelihood.log_likelihood(simulator.evaluate(trial_sample)) +
1031 * prior.log_prior(trial_sample));
1033 * if (std::exp(trial_log_posterior - current_log_posterior) * perturbation_probability_ratio
1035 * uniform_distribution(random_number_generator))
1037 * current_sample = trial_sample;
1038 * current_log_posterior = trial_log_posterior;
1040 * ++accepted_sample_number;
1043 * write_sample(current_sample, current_log_posterior);
1049 * void MetropolisHastings::write_sample(const Vector<double> ¤t_sample,
1050 * const double current_log_posterior)
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
';
1058 * output_file.flush();
1060 * } // namespace Sampler
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>
1077 * /* Set the exact coefficient: */
1078 * Vector<double> exact_coefficients(64);
1079 * for (auto &el : exact_coefficients)
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.;
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);
1099 * const bool testing = true;
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.
1109 * MultithreadInfo::set_thread_limit(1);
1111 * const unsigned int random_seed = (testing ? 1U : std::random_device()());
1112 * const std::string dataset_name = Utilities::to_string(random_seed, 10);
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 });
1203 * Now run the forward simulator for samples:
1206 * ForwardSimulator::PoissonSolver<2> laplace_problem(
1207 * /* global_refinements = */ 5,
1208 * /* fe_degree = */ 1,
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,
1217 * proposal_generator,
1221 * Vector<double> starting_coefficients(64);
1222 * for (auto &el : starting_coefficients)
1224 * sampler.sample(starting_coefficients,
1225 * (testing ? 250 * 40 /* takes 40 seconds */
1227 * 100000000 /* takes 6 days */