491 const unsigned int n_theta = 64;
492 for (
unsigned int test=0; test<10; ++test)
494 std::cout <<
"Generating output for test " << test << std::endl;
497 std::ifstream
test_input (
"../testing/input." + std::to_string(test) +
".txt");
501 for (
unsigned int i=0; i<
n_theta; ++i)
507 std::ofstream
test_output_z (
"output." + std::to_string(test) +
".z.txt");
516 std::ofstream
test_output_prior (
"output." + std::to_string(test) +
".logprior.txt");
530---------------------------------------
543---------------------------------------
564<a name=
"ann-Matlab/exact_values.m"></a>
567%% -----------------------------------------------------------------------------
574%% -----------------------------------------------------------------------------
576%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
578%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
580z_hat = [0.06076511762259369;
752<a name=
"ann-Matlab/forward_solver.m"></a>
755%% -----------------------------------------------------------------------------
762%% -----------------------------------------------------------------------------
764%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
765%%%%%%%%%%%%%%%%%%%%%% forward solver function %%%%%%%%%%%%%%%%%%%%%%%%%%%%
766%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
775%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
778%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
805%solve linear
equation for coefficients,
U
813<a name=
"ann-Matlab/get_statistics.m"></a>
816%% -----------------------------------------------------------------------------
823%% -----------------------------------------------------------------------------
825%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
827%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
831%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
836%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
864%compute mean of autocovariance matrix
865autocovar = autocovar(1:64,1:64,L:2*L-1)/N;
869<a name="ann-Matlab/log_probability.m"></a>
870<h1>Annotated version of Matlab/log_probability.m</h1>
872%% -----------------------------------------------------------------------------
874%% SPDX-License-Identifier: LGPL-2.1-or-later
875%% Copyright (C) 2022 by Wolfgang Bangerth
877%% This file is part of the deal.II code gallery.
879%% -----------------------------------------------------------------------------
881%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
882%%%%%%%%%%%%%%%%% compute log probability, log pi %%%%%%%%%%%%%%%%%%%%%%%%%
883%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
885%theta = current 8x8 parameter matrix
886%z = current vector of measurements
887%z_hat = vector of "exact" measurements
888%sig = standard deviation parameter in likelihood
889%sig_pr = standard deviation parameter in prior
890%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
892%log_pi = logarithm of posterior probability
893%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
895function log_pi = log_probability(theta,z,z_hat,sig,sig_pr)
897%compute log likelihood
898log_L = -sum((z-z_hat).^2)/(2*sig^2);
901log_pi_pr = -sum(log(theta).^2,'all
')/(2*sig_pr^2);
903%compute log posterior
904log_pi = log_L+log_pi_pr;
908<a name="ann-Matlab/main.m"></a>
909<h1>Annotated version of Matlab/main.m</h1>
911%% -----------------------------------------------------------------------------
913%% SPDX-License-Identifier: LGPL-2.1-or-later
914%% Copyright (C) 2022 by Wolfgang Bangerth
916%% This file is part of the deal.II code gallery.
918%% -----------------------------------------------------------------------------
920%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
921%%%%%%%%% Run MCMC sampler to estimate posterior distribution %%%%%%%%%%%%%
922%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
924%define number of chains, chain length, and lag time
931%open Matlab parallel pool
935load precomputations.mat
937%define lag time and data matrix
938data = zeros(8,8,L,N); %data matrix of samples at lag times
939theta_means = zeros(8,8,N); %overall mean of theta
945 %set initial theta, theta mean, and z values of chain
954 %define proposal, theta_tilde
955 xi = normrnd(0,sig_prop,[8 8]);
956 theta_tilde = theta.*exp(xi);
958 %compute new z values
959 z_tilde = forward_solver_(theta_tilde);
961 %compute posterior log probability of theta_tilde
962 log_pi_tilde = log_probability_(theta_tilde,z_tilde);
963 log_pi = log_probability_(theta,z);
965 %compute acceptance probability; accept proposal appropriately
966 accept = exp(log_pi_tilde-log_pi) ...
967 *prod(theta_tilde./theta,'all
');
969 theta = theta_tilde; %accept new theta values
970 z = z_tilde; %record associated measurements
973 %update mean of theta
974 theta_mean = theta_mean + theta;
979 data(:,:,m,n) = theta;
984 theta_means(:,:,n) = theta_mean/N_L;
990%compute statistics on data set
991[theta_mean,covars,autocovar] = get_statistics(data,theta_means);
993%save data to Matlab workspace, labeled by N and N_L
998<a name="ann-Matlab/plot_solution.m"></a>
999<h1>Annotated version of Matlab/plot_solution.m</h1>
1001%% -----------------------------------------------------------------------------
1003%% SPDX-License-Identifier: LGPL-2.1-or-later
1004%% Copyright (C) 2022 by Wolfgang Bangerth
1006%% This file is part of the deal.II code gallery.
1008%% -----------------------------------------------------------------------------
1010%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1011%%%%%%%%%%%%%%% plots solution, u, to Poisson equation %%%%%%%%%%%%%%%%%%%%
1012%%%%%%%%%%%%%%%% associated to theta and a 32x32 mesh %%%%%%%%%%%%%%%%%%%%%
1013%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1015%input matrix theta for plotting (e.g., theta = theta_hat)
1018%construct mass matrix, M_plot, for plotting
1021Mp = zeros(n,n,33^2);
1026 Mp(i,j,k) = phi(xsp(i)-h*c(1),xsp(j)-h*c(2));
1030Mp = reshape(Mp,[n^2 33^2]);
1032%run forward solver on mean of theta
1033A = zeros(33^2,33^2);
1035 for j=0:31 %build A by summing over contribution from each cell
1037 %find local coefficient in 8x8 grid
1038 theta_loc = theta(floor(i/4)+1,floor(j/4)+1);
1040 %update A by including contribution from cell (i,j)
1041 dof = [lbl(i,j),lbl(i,j+1),lbl(i+1,j+1),lbl(i+1,j)];
1042 A(dof,dof) = A(dof,dof) + theta_loc*A_loc;
1046%enforce boundary condition
1047A(boundaries,:) = Id(boundaries,:);
1052%solve linear equation for coefficients, U
1055%close all current plots
1060zs = reshape(Mp*U,[n n]);
1065<a name="ann-Matlab/precomputations.m"></a>
1066<h1>Annotated version of Matlab/precomputations.m</h1>
1068%% -----------------------------------------------------------------------------
1070%% SPDX-License-Identifier: LGPL-2.1-or-later
1071%% Copyright (C) 2022 by Wolfgang Bangerth
1073%% This file is part of the deal.II code gallery.
1075%% -----------------------------------------------------------------------------
1077%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1078%%%%%%% do all precomputations necessary for MCMC simulations %%%%%%%%%%%%%
1079%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1084%define characteristic function of unit square
1085S = @(x,y) heaviside(x).*heaviside(y) ...
1086 .*(1-heaviside(x-h)).*(1-heaviside(y-h));
1088%define tent function on the domain [-h,h]x[-h,h]
1089phi = @(x,y) ((x+h).*(y+h).*S(x+h,y+h) + (h-x).*(h-y).*S(x,y) ...
1090 + (x+h).*(h-y).*S(x+h,y) + (h-x).*(y+h).*S(x,y+h))/h^2;
1092%define function that converts from (i,j) to dof, and its inverse
1093lbl = @(i,j) 33*j+i+1;
1094inv_lbl = @(k) [k-1-33*floor((k-1)/33),floor((k-1)/33)];
1096%construct measurement matrix, M
1097xs = 1/14:1/14:13/14; %measurement points
1098M = zeros(13,13,33^2);
1103 M(i,j,k) = phi(xs(i)-h*c(1),xs(j)-h*c(2));
1107M = reshape(M,[13^2 33^2]);
1109%construct exact coefficient matrix, theta_hat
1110theta_hat = ones(8,8);
1111theta_hat(2:3,2:3) = 0.1;
1112theta_hat(6:7,6:7) = 10;
1114%construct local overlap matrix, A_loc, and identity matrix Id
1115A_loc = [2/3 -1/6 -1/3 -1/6;
1118 -1/6 -1/3 -1/6 2/3];
1121%locate boundary labels
1122boundaries = [lbl(0:1:32,0),lbl(0:1:32,32),lbl(0,1:1:31),lbl(32,1:1:31)];
1124%define RHS of FEM linear system, AU = b
1125b = ones(33^2,1)*10*h^2;
1126b(boundaries) = zeros(128,1); %enforce boundary conditions on b
1128%load exact z_hat values
1131%set global parameters and functions for simulation
1132sig = 0.05; %likelihood standard deviation
1133sig_pr = 2; %prior (log) standard deviation
1134sig_prop = 0.0725; %proposal (log) standard deviation
1135theta0 = ones(8,8); %initial theta values
1136forward_solver_ = @(theta) ...
1137 forward_solver(theta,lbl,A_loc,Id,boundaries,b,M);
1138log_probability_ = @(theta,z) log_probability(theta,z,z_hat,sig,sig_pr);
1140%find initial z values
1141z0 = forward_solver_(theta0);
1143save precomputations.mat
1147<a name="ann-Python/forward_solver.py"></a>
1148<h1>Annotated version of Python/forward_solver.py</h1>
1150## -----------------------------------------------------------------------------
1152## SPDX-License-Identifier: LGPL-2.1-or-later
1153## Copyright (C) 2022 by Wolfgang Bangerth
1155## This file is part of the deal.II code gallery.
1157## -----------------------------------------------------------------------------
1161from scipy.sparse.linalg import spsolve
1164###########################################################################
1165############ list of "exact" measurement values, z_hat ####################
1166###########################################################################
1169 [0.06076511762259369, 0.09601910120848481,
1170 0.1238852517838584, 0.1495184117375201,
1171 0.1841596127549784, 0.2174525028261122,
1172 0.2250996160898698, 0.2197954769002993,
1173 0.2074695698370926, 0.1889996477663016,
1174 0.1632722532153726, 0.1276782480038186,
1175 0.07711845915789312, 0.09601910120848552,
1176 0.2000589533367983, 0.3385592591951766,
1177 0.3934300024647806, 0.4040223892461541,
1178 0.4122329537843092, 0.4100480091545554,
1179 0.3949151637189968, 0.3697873264791232,
1180 0.33401826235924, 0.2850397806663382,
1181 0.2184260032478671, 0.1271121156350957,
1182 0.1238852517838611, 0.3385592591951819,
1183 0.7119285162766475, 0.8175712861756428,
1184 0.6836254116578105, 0.5779452419831157,
1185 0.5555615956136897, 0.5285181561736719,
1186 0.491439702849224, 0.4409367494853282,
1187 0.3730060082060772, 0.2821694983395214,
1188 0.1610176733857739, 0.1495184117375257,
1189 0.3934300024647929, 0.8175712861756562,
1190 0.9439154625527653, 0.8015904115095128,
1191 0.6859683749254024, 0.6561235366960599,
1192 0.6213197201867315, 0.5753611315000049,
1193 0.5140091754526823, 0.4325325506354165,
1194 0.3248315148915482, 0.1834600412730086,
1195 0.1841596127549917, 0.4040223892461832,
1196 0.6836254116578439, 0.8015904115095396,
1197 0.7870119561144977, 0.7373108331395808,
1198 0.7116558878070463, 0.6745179049094283,
1199 0.6235300574156917, 0.5559332704045935,
1200 0.4670304994474178, 0.3499809143811,
1201 0.19688263746294, 0.2174525028261253,
1202 0.4122329537843404, 0.5779452419831566,
1203 0.6859683749254372, 0.7373108331396063,
1204 0.7458811983178246, 0.7278968022406559,
1205 0.6904793535357751, 0.6369176452710288,
1206 0.5677443693743215, 0.4784738764865867,
1207 0.3602190632823262, 0.2031792054737325,
1208 0.2250996160898818, 0.4100480091545787,
1209 0.5555615956137137, 0.6561235366960938,
1210 0.7116558878070715, 0.727896802240657,
1211 0.7121928678670187, 0.6712187391428729,
1212 0.6139157775591492, 0.5478251665295381,
1213 0.4677122687599031, 0.3587654911000848,
1214 0.2050734291675918, 0.2197954769003094,
1215 0.3949151637190157, 0.5285181561736911,
1216 0.6213197201867471, 0.6745179049094407,
1217 0.690479353535786, 0.6712187391428787,
1218 0.6178408289359514, 0.5453605027237883,
1219 0.489575966490909, 0.4341716881061278,
1220 0.3534389974779456, 0.2083227496961347,
1221 0.207469569837099, 0.3697873264791366,
1222 0.4914397028492412, 0.5753611315000203,
1223 0.6235300574157017, 0.6369176452710497,
1224 0.6139157775591579, 0.5453605027237935,
1225 0.4336604929612851, 0.4109641743019312,
1226 0.3881864790111245, 0.3642640090182592,
1227 0.2179599909280145, 0.1889996477663011,
1228 0.3340182623592461, 0.4409367494853381,
1229 0.5140091754526943, 0.5559332704045969,
1230 0.5677443693743304, 0.5478251665295453,
1231 0.4895759664908982, 0.4109641743019171,
1232 0.395727260284338, 0.3778949322004734,
1233 0.3596268271857124, 0.2191250268948948,
1234 0.1632722532153683, 0.2850397806663325,
1235 0.373006008206081, 0.4325325506354207,
1236 0.4670304994474315, 0.4784738764866023,
1237 0.4677122687599041, 0.4341716881061055,
1238 0.388186479011099, 0.3778949322004602,
1239 0.3633362567187364, 0.3464457261905399,
1240 0.2096362321365655, 0.1276782480038148,
1241 0.2184260032478634, 0.2821694983395252,
1242 0.3248315148915535, 0.3499809143811097,
1243 0.3602190632823333, 0.3587654911000799,
1244 0.3534389974779268, 0.3642640090182283,
1245 0.35962682718569, 0.3464457261905295,
1246 0.3260728953424643, 0.180670595355394,
1247 0.07711845915789244, 0.1271121156350963,
1248 0.1610176733857757, 0.1834600412730144,
1249 0.1968826374629443, 0.2031792054737354,
1250 0.2050734291675885, 0.2083227496961245,
1251 0.2179599909279998, 0.2191250268948822,
1252 0.2096362321365551, 0.1806705953553887,
1253 0.1067965550010013])
1256###########################################################################
1257####### do all precomputations necessary for MCMC simulations #############
1258###########################################################################
1260# Define the mesh width
1263# Define characteristic function of unit square
1271 return heaviside(x)*heaviside(y) * (1-heaviside(x-h))*(1-heaviside(y-h));
1273# Define tent function on the domain [0,2h]x[0,2h]
1275 return ((x+h)*(y+h)*S(x+h,y+h) + (h-x)*(h-y)*S(x,y)
1276 + (x+h)*(h-y)*S(x+h,y) + (h-x)*(y+h)*S(x,y+h))/h**2
1290M =
np.zeros((13,13,33**2));
1299# Construct local overlap matrix, A_loc, and identity matrix Id
1300A_loc =
np.array([[2./3, -1./6, -1./3, -1./6],
1301 [-1./6, 2./3, -1./6, -1./3],
1302 [-1./3, -1./6, 2./3, -1./6],
1303 [-1./6, -1./3, -1./6, 2./3]])
1304Id =
np.eye(33**2,33**2)
1306# Locate boundary labels
1312# Define RHS of FEM linear system, AU = b
1313b =
np.ones(33**2)*10*h**2
1320###########################################################################
1321###################### forward solver function ############################
1322###########################################################################
1358###########################################################################
1360###########################################################################
1377###########################################################################
1378#############
A function
to test
against known output #####################
1379###########################################################################
1384 print (
"Verifying against data set", i)
1392 # solver
twice),
but we don't care about efficiency here since
1393 # we are only computing with ten samples
1394 this_z = forward_solver(theta)
1395 this_log_likelihood = log_likelihood(theta)
1396 this_log_prior = log_prior(theta)
1398 # Then also read the reference output generated by the C++ program:
1399 f_output_z = open ("../testing/output.{}.z.txt".format(i), 'r')
1400 f_output_likelihood = open ("../testing/output.{}.loglikelihood.txt".format(i), 'r')
1401 f_output_prior = open ("../testing/output.{}.logprior.txt".format(i), 'r')
1403 reference_z = np.fromfile(f_output_z, count=13**2, sep=" ")
1404 reference_log_likelihood = float(f_output_likelihood.read())
1405 reference_log_prior = float(f_output_prior.read())
1407 print (" || z-z_ref || : ",
1408 np.linalg.norm(this_z - reference_z))
1409 print (" log likelihood : ",
1410 "Python value=", this_log_likelihood,
1411 "(C++ reference value=", reference_log_likelihood,
1412 ", error=", abs(this_log_likelihood - reference_log_likelihood),
1414 print (" log prior : ",
1415 "Python value=", this_log_prior,
1416 "(C++ reference value=", reference_log_prior,
1417 ", error=", abs(this_log_prior - reference_log_prior),
1421def time_forward_solver() :
1425 for i in range(n_runs) :
1426 # Create a random vector (with entries between 0 and 1), scale
1427 # it by a factor of 4, subtract 2, then take the exponential
1428 # of each entry to get random entries between e^{-2} and
1430 theta = np.exp(np.random.rand(64) * 4 - 2)
1431 z = forward_solver(theta)
1433 print ("Time per forward evaluation:", (end-begin)/n_runs)
1435verify_against_stored_tests()
1436time_forward_solver()
1440<a name="ann-mcmc-laplace.cc"></a>
1441<h1>Annotated version of mcmc-laplace.cc</h1>
1447 * /* -----------------------------------------------------------------------------
1449 * * SPDX-License-Identifier: LGPL-2.1-or-later
1450 * * Copyright (C) 2019 by Wolfgang Bangerth
1452 * * This file is part of the deal.II code gallery.
1454 * * -----------------------------------------------------------------------------
1456 * * Author: Wolfgang Bangerth, Colorado State University, 2019.
1460 * #include <deal.II/base/timer.h>
1461 * #include <deal.II/base/multithread_info.h>
1462 * #include <deal.II/grid/tria.h>
1463 * #include <deal.II/dofs/dof_handler.h>
1464 * #include <deal.II/grid/grid_generator.h>
1465 * #include <deal.II/grid/tria_accessor.h>
1466 * #include <deal.II/grid/tria_iterator.h>
1467 * #include <deal.II/dofs/dof_accessor.h>
1468 * #include <deal.II/fe/fe_q.h>
1469 * #include <deal.II/dofs/dof_tools.h>
1470 * #include <deal.II/fe/fe_values.h>
1471 * #include <deal.II/base/quadrature_lib.h>
1472 * #include <deal.II/base/function.h>
1473 * #include <deal.II/numerics/vector_tools.h>
1474 * #include <deal.II/numerics/matrix_tools.h>
1475 * #include <deal.II/lac/vector.h>
1476 * #include <deal.II/lac/full_matrix.h>
1477 * #include <deal.II/lac/sparse_matrix.h>
1478 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
1479 * #include <deal.II/lac/solver_cg.h>
1480 * #include <deal.II/lac/precondition.h>
1481 * #include <deal.II/lac/sparse_ilu.h>
1483 * #include <deal.II/numerics/data_out.h>
1485 * #include <fstream>
1486 * #include <iostream>
1489 * #include <deal.II/base/logstream.h>
1491 * using namespace dealii;
1496 * The following is a namespace in which we define the solver of the PDE.
1497 * The main class implements an abstract `Interface` class declared at
1498 * the top, which provides for an `evaluate()` function that, given
1499 * a coefficient vector, solves the PDE discussed in the Readme file
1500 * and then evaluates the solution at the 169 mentioned points.
1504 * The solver follows the basic layout of @ref step_4 "step-4", though it precomputes
1505 * a number of things in the `setup_system()` function, such as the
1506 * evaluation of the matrix that corresponds to the point evaluations,
1507 * as well as the local contributions to matrix and right hand side.
1511 * Rather than commenting on everything in detail, in the following
1512 * we will only document those things that are not already clear from
1513 * @ref step_4 "step-4" and a small number of other tutorial programs.
1516 * namespace ForwardSimulator
1521 * virtual Vector<double> evaluate(const Vector<double> &coefficients) = 0;
1523 * virtual ~Interface() = default;
1528 * template <int dim>
1529 * class PoissonSolver : public Interface
1532 * PoissonSolver(const unsigned int global_refinements,
1533 * const unsigned int fe_degree,
1534 * const std::string &dataset_name);
1535 * virtual Vector<double>
1536 * evaluate(const Vector<double> &coefficients) override;
1539 * void make_grid(const unsigned int global_refinements);
1540 * void setup_system();
1541 * void assemble_system(const Vector<double> &coefficients);
1543 * void output_results(const Vector<double> &coefficients) const;
1545 * Triangulation<dim> triangulation;
1547 * DoFHandler<dim> dof_handler;
1549 * FullMatrix<double> cell_matrix;
1550 * Vector<double> cell_rhs;
1551 * std::map<types::global_dof_index,double> boundary_values;
1553 * SparsityPattern sparsity_pattern;
1554 * SparseMatrix<double> system_matrix;
1556 * Vector<double> solution;
1557 * Vector<double> system_rhs;
1559 * std::vector<Point<dim>> measurement_points;
1561 * SparsityPattern measurement_sparsity;
1562 * SparseMatrix<double> measurement_matrix;
1564 * TimerOutput timer;
1565 * unsigned int nth_evaluation;
1567 * const std::string &dataset_name;
1572 * template <int dim>
1573 * PoissonSolver<dim>::PoissonSolver(const unsigned int global_refinements,
1574 * const unsigned int fe_degree,
1575 * const std::string &dataset_name)
1577 * , dof_handler(triangulation)
1578 * , timer(std::cout, TimerOutput::summary, TimerOutput::cpu_times)
1579 * , nth_evaluation(0)
1580 * , dataset_name(dataset_name)
1582 * make_grid(global_refinements);
1588 * template <int dim>
1589 * void PoissonSolver<dim>::make_grid(const unsigned int global_refinements)
1591 * Assert(global_refinements >= 3,
1592 * ExcMessage("This program makes the assumption that the mesh for the "
1593 * "solution of the PDE is at least as fine as the one used "
1594 * "in the definition of the coefficient."));
1595 * GridGenerator::hyper_cube(triangulation, 0, 1);
1596 * triangulation.refine_global(global_refinements);
1598 * std::cout << " Number of active cells: " << triangulation.n_active_cells()
1604 * template <int dim>
1605 * void PoissonSolver<dim>::setup_system()
1609 * First define the finite element space:
1612 * dof_handler.distribute_dofs(fe);
1614 * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1619 * Then set up the main data structures that will hold the discrete problem:
1623 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
1624 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
1625 * sparsity_pattern.copy_from(dsp);
1627 * system_matrix.reinit(sparsity_pattern);
1629 * solution.reinit(dof_handler.n_dofs());
1630 * system_rhs.reinit(dof_handler.n_dofs());
1635 * And then define the tools to do point evaluation. We choose
1636 * a set of 13x13 points evenly distributed across the domain:
1640 * const unsigned int n_points_per_direction = 13;
1641 * const double dx = 1. / (n_points_per_direction + 1);
1643 * for (unsigned int x = 1; x <= n_points_per_direction; ++x)
1644 * for (unsigned int y = 1; y <= n_points_per_direction; ++y)
1645 * measurement_points.emplace_back(x * dx, y * dx);
1649 * First build a full matrix of the evaluation process. We do this
1650 * even though the matrix is really sparse -- but we don't
know
1656 *
Vector<double> weights(dof_handler.n_dofs());
1659 *
dof_handler.n_dofs());
1666 *
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1681 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1684 *
const unsigned int j = std::floor(cell->center()[1] * 8);
1686 *
const unsigned int index = i + 8 *
j;
1688 *
cell->set_user_index(index);
1698 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1700 *
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1711 *
fe_values.reinit(dof_handler.begin_active());
1713 *
for (
unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1714 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1716 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1718 *
(fe_values.shape_grad(i, q_index) *
1719 *
fe_values.shape_grad(
j, q_index) *
1720 *
fe_values.JxW(q_index));
1722 *
cell_rhs(i) += (fe_values.shape_value(i, q_index) *
1724 *
fe_values.JxW(q_index));
1743 *
template <
int dim>
1746 *
Assert(coefficients.size() == 64, ExcInternalError());
1748 *
system_matrix = 0;
1751 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1753 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1755 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1759 *
cell->get_dof_indices(local_dof_indices);
1760 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1762 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1763 *
system_matrix.add(local_dof_indices[i],
1764 *
local_dof_indices[
j],
1783 *
template <
int dim>
1787 *
ilu.initialize(system_matrix);
1805 *
template <
int dim>
1810 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1812 *
coefficients(cell->user_index());
1816 *
data_out.attach_dof_handler(dof_handler);
1817 *
data_out.add_data_vector(solution,
"solution");
1820 *
data_out.build_patches();
1822 *
std::ofstream output(
"solution-" +
1824 *
data_out.write_vtu(output);
1852 *
template <
int dim>
1872 *
ExcInternalError());
1879 *
timer.print_summary();
1909 * virtual double log_likelihood(const Vector<double> &x) const = 0;
1911 * virtual ~Interface() = default;
1915 *
class Gaussian :
public Interface
1924 *
const double sigma;
1967 *
virtual ~Interface() =
default;
1980 *
const double sigma;
1993 *
for (
const auto &el :
x)
2009 * positive, we use a Gaussian distribution in logarithm space -- in other
2010 * words, instead of *adding* a small perturbation with mean value zero,
2011 * we *multiply* the entries of the current sample by a factor that
2012 * is the exponential of a random number with mean zero. (Because the
2013 * exponential of zero is one, this means that the most likely factors
2014 * to multiply the existing sample entries by are close to one. And
2015 * because the exponential of a number is always positive, we never
2016 * get negative samples this way.)
2027 *
that's not the case for the Gaussian in log space. It's
not
2040 *
std::pair<Vector<double>,
double>
2043 *
virtual ~Interface() =
default;
2053 *
std::pair<Vector<double>,
double>
2071 *
std::pair<Vector<double>,
double>
2124 *
const LogLikelihood::Interface &
likelihood,
2125 *
const LogPrior::Interface &
prior,
2134 *
ForwardSimulator::Interface &
simulator;
2135 *
const LogLikelihood::Interface &
likelihood;
2136 *
const LogPrior::Interface &
prior;
2152 *
ForwardSimulator::Interface &
simulator,
2153 *
const LogLikelihood::Interface &
likelihood,
2154 *
const LogPrior::Interface &
prior,
2188 *
std::pair<Vector<double>,
double>
2279 *
{ 0.06076511762259369, 0.09601910120848481,
2280 *
0.1238852517838584, 0.1495184117375201,
2281 *
0.1841596127549784, 0.2174525028261122,
2282 *
0.2250996160898698, 0.2197954769002993,
2283 *
0.2074695698370926, 0.1889996477663016,
2284 *
0.1632722532153726, 0.1276782480038186,
2285 *
0.07711845915789312, 0.09601910120848552,
2286 *
0.2000589533367983, 0.3385592591951766,
2287 *
0.3934300024647806, 0.4040223892461541,
2288 *
0.4122329537843092, 0.4100480091545554,
2289 *
0.3949151637189968, 0.3697873264791232,
2290 *
0.33401826235924, 0.2850397806663382,
2291 *
0.2184260032478671, 0.1271121156350957,
2292 *
0.1238852517838611, 0.3385592591951819,
2293 *
0.7119285162766475, 0.8175712861756428,
2294 *
0.6836254116578105, 0.5779452419831157,
2295 *
0.5555615956136897, 0.5285181561736719,
2296 *
0.491439702849224, 0.4409367494853282,
2297 *
0.3730060082060772, 0.2821694983395214,
2298 *
0.1610176733857739, 0.1495184117375257,
2299 *
0.3934300024647929, 0.8175712861756562,
2300 *
0.9439154625527653, 0.8015904115095128,
2301 *
0.6859683749254024, 0.6561235366960599,
2302 *
0.6213197201867315, 0.5753611315000049,
2303 *
0.5140091754526823, 0.4325325506354165,
2304 *
0.3248315148915482, 0.1834600412730086,
2305 *
0.1841596127549917, 0.4040223892461832,
2306 *
0.6836254116578439, 0.8015904115095396,
2307 *
0.7870119561144977, 0.7373108331395808,
2308 *
0.7116558878070463, 0.6745179049094283,
2309 *
0.6235300574156917, 0.5559332704045935,
2310 *
0.4670304994474178, 0.3499809143811,
2311 *
0.19688263746294, 0.2174525028261253,
2312 *
0.4122329537843404, 0.5779452419831566,
2313 *
0.6859683749254372, 0.7373108331396063,
2314 *
0.7458811983178246, 0.7278968022406559,
2315 *
0.6904793535357751, 0.6369176452710288,
2316 *
0.5677443693743215, 0.4784738764865867,
2317 *
0.3602190632823262, 0.2031792054737325,
2318 *
0.2250996160898818, 0.4100480091545787,
2319 *
0.5555615956137137, 0.6561235366960938,
2320 *
0.7116558878070715, 0.727896802240657,
2321 *
0.7121928678670187, 0.6712187391428729,
2322 *
0.6139157775591492, 0.5478251665295381,
2323 *
0.4677122687599031, 0.3587654911000848,
2324 *
0.2050734291675918, 0.2197954769003094,
2325 *
0.3949151637190157, 0.5285181561736911,
2326 *
0.6213197201867471, 0.6745179049094407,
2327 *
0.690479353535786, 0.6712187391428787,
2328 *
0.6178408289359514, 0.5453605027237883,
2329 *
0.489575966490909, 0.4341716881061278,
2330 *
0.3534389974779456, 0.2083227496961347,
2331 *
0.207469569837099, 0.3697873264791366,
2332 *
0.4914397028492412, 0.5753611315000203,
2333 *
0.6235300574157017, 0.6369176452710497,
2334 *
0.6139157775591579, 0.5453605027237935,
2335 *
0.4336604929612851, 0.4109641743019312,
2336 *
0.3881864790111245, 0.3642640090182592,
2337 *
0.2179599909280145, 0.1889996477663011,
2338 *
0.3340182623592461, 0.4409367494853381,
2339 *
0.5140091754526943, 0.5559332704045969,
2340 *
0.5677443693743304, 0.5478251665295453,
2341 *
0.4895759664908982, 0.4109641743019171,
2342 *
0.395727260284338, 0.3778949322004734,
2343 *
0.3596268271857124, 0.2191250268948948,
2344 *
0.1632722532153683, 0.2850397806663325,
2345 *
0.373006008206081, 0.4325325506354207,
2346 *
0.4670304994474315, 0.4784738764866023,
2347 *
0.4677122687599041, 0.4341716881061055,
2348 *
0.388186479011099, 0.3778949322004602,
2349 *
0.3633362567187364, 0.3464457261905399,
2350 *
0.2096362321365655, 0.1276782480038148,
2351 *
0.2184260032478634, 0.2821694983395252,
2352 *
0.3248315148915535, 0.3499809143811097,
2353 *
0.3602190632823333, 0.3587654911000799,
2354 *
0.3534389974779268, 0.3642640090182283,
2355 *
0.35962682718569, 0.3464457261905295,
2356 *
0.3260728953424643, 0.180670595355394,
2357 *
0.07711845915789244, 0.1271121156350963,
2358 *
0.1610176733857757, 0.1834600412730144,
2359 *
0.1968826374629443, 0.2031792054737354,
2360 *
0.2050734291675885, 0.2083227496961245,
2361 *
0.2179599909279998, 0.2191250268948822,
2362 *
0.2096362321365551, 0.1806705953553887,
2363 *
0.1067965550010013 });
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
std::vector< index_type > data
Expression floor(const Expression &x)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int one
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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation