{'fraction_non_alphanumeric': 0.037398426260077, 'fraction_numerical': 0.009436009548534364, 'mean_word_length': 4.706462413429344, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian model calibration for block copolymer self-assembly: Likelihood-free inference and expected information gain computation via measure transport

We consider the Bayesian calibration of models describing the phenomenon of block copolymer (BCP) self-assembly using image data produced by microscopy or X-ray scattering techniques. To account for the random long-range disorder in BCP equilibrium structures, we introduce auxiliary variables to represent this aleatory uncertainty. These variables, however, result in an integrated likelihood for high-dimensional image data that is generally intractable to evaluate. We tackle this challenging Bayesian inference problem using a likelihood-free approach based on measure transport together with the construction of summary statistics for the image data. We also show that expected information gains (EIGs) from the observed data about the model parameters can be computed with no significant additional cost. Lastly, we present a numerical case study based on the Ohta--Kawasaki model for diblock copolymer thin film self-assembly and top-down microscopy characterization. For calibration, we introduce several domain-specific energy- and Fourier-based summary statistics, and quantify their informativeness using EIG. We demonstrate the power of the proposed approach to study the effect of data corruptions and experimental designs on the calibration results.


The self-assembly of block copolymers: Application and modeling
Block copolymer (BCP) melts are large collections of one or more types of BCPs, or polymers consisting of several blocks of distinct monomers joined together by covalent bonds. Originally mixed in high temperature equilibrium, the distinct monomer blocks in a BCP melt spontaneously segregate upon thermal or solvent annealing as a result of their thermodynamic incompatibility. Periodically ordered nanoscale structures emerge as the melt achieves a new equilibrium configuration. This physical phenomenon is known as A B C Figure 1: (left) The experimental phase diagram for the self-assembly of Polyisoprene-Polystyrene diblock copolymers. The parameter f P I denotes the volume ratio of the Polyisoprene block, χ denotes the Flory-Huggins interaction parameter, and N denotes the degree of polymerization. Reprinted (adapted) with permission from [13]. Copyright ©1995 American Chemical Society. (right) Top-down scanning electron microscope images of lamella phase PS-b-PMMA copolymer films with (A) a chemically neutral surface; (B) a chemically nanopatterned surface with period smaller than the natural period of the copolymer, thus leading to defects; (C) a chemically nanopatterned surface with period close to the natural period of the copolymer. Reprinted (adapted) with permission from [4]. Copyright ©2003 Springer Nature.
the self-assembly of BCPs [1,2]. Examples of the equilibrium structures for the self-assembly of diblock copolymers (Di-BCPs), i.e., BCPs with two distinct blocks, are lamellae, spheres, cylinders, etc.; see Fig. 1 (left) along with an experimental phase diagram. The massive attention given to BCP self-assembly stems from the flexibility in their achievable equilibrium structures when subject to external guidance. Specifically, carefully engineered techniques, such as substrate-based techniques [3][4][5] or techniques using thermal and electric fields [6][7][8], can influence the self-assembly process to form equilibrium structures with desired features. The application of these techniques is known as directed self-assembly (DSA) of BCPs. The DSA of BCPs has quickly attracted attention for its wide-ranging commercial applications, the most notable being patterning for lithography at the sub-10 nm scale [9][10][11][12].
Computer simulation of BCP self-assembly based on mathematical models has played an important role in the study of self-assembly. These simulations have proven to be valuable in determining the effects of material properties, film thickness, and polymer-substrate interactions on the self-assembly [14]. The popular models and simulation techniques for self-assembly are predominantly derived from statistical polymer theory, such as Monte Carlo simulations [15,16], theoretically-informed coarse-grained (TICG) simulation [17], self-consistent field theory (SCFT) [18,19], and density functional theory (DFT) [20,21]-listed in decreasing order by granularity 2 . While model-based studies of the self-assembly mainly serve as aides to interpreting empirical findings, there is an increased interest in using models to make potentially consequential predictions. Most noticeably, recent efforts have been made on the inverse design 3 of DSA, where model simulations and iterative optimization schemes suggest designs of external guidance and BCP material properties that are best suited for achieving a target equilibrium pattern [22][23][24][25][26]. Model-based inverse design is often considered the critical pathway towards transitioning DSA technology from the laboratory into commercial manufacturing of integrated circuit devices [27].

Calibration of block copolymer self-assembly models: Uncertainties and the Bayesian approach
With the increasing usage and growing impact of model-based predictions for BCP self-assembly, it is increasingly important to ensure the reliability of models. Model parameter estimation, or model calibration, is necessary in order to produce reliable model-based predictions [28] and is the focus of this work. In particular, we formulate and propose a suitable methodology for calibrating models of BCP self-assembly. Models are calibrated according to observational data obtained by experimental characterizations [29,Vol. 2], via microscopy or X-ray scattering techniques (see, for example, Fig. 2), of BCP melts, and various intrinsic uncertainties are also explicitly accounted for.
One source of intrinsic uncertainty stems from corruptions introduced in the characterization process. This is often modeled as stochastic sensor noise. We also address an aleatoric uncertainty, related to the well-known phenomenon that unguided or weakly guided self-assembly of BCPs results in a range of possible equilibrium morphologies with long-range disorder. A characteristic example for such long-range disorder is the 'finger-print' patterns for the lamellae phase of Di-BCPs thin films under unguided self-assembly. For DSA, the aleatoric uncertainty results in defects in the equilibrium patterns. Examples of a 'finger-print' pattern and DSA pattern defects are shown in Fig. 1 (right). Pattern defects have significant implications on product yield and quality in the nanolithography application of BCP self-assembly [11,30].
Statistical polymer theory [31] accounts for the impact of aleatoric uncertainty by describing the variety of achievable equilibrium morphologies for a particular polymer melt system via a statistical ensemble. The statistical ensemble assigns a probability of occurrence for each achievable equilibrium morphology [32]. In this work, we consider all BCP self-assembly models as statistical models, through which the computer-generated solutions are samples from a model-specific statistical ensemble. We assume the stochastic component of the models is accounted for by auxiliary random variables which are possibly infinite-dimensional and/or parameter-dependent.
To perform model calibration that accounts for these uncertainties, we advocate for a Bayesian probabilistic approach that uses probability distributions to describe the uncertainty in model parameters. Under such an approach, the output of the model calibration is a posterior probability distribution of the model parameters conditioned on the observed data. In addition to performing parameter estimation, we are also interested in studying the effects of instrument-specific parameters that describe the strength of various kinds of corruptions introduced in the characterization process. These effects can be quantified by the amount of information content, i.e., the Shannon entropy [33], preserved by the characterization process. Doing so informs decision-making for appropriate specifications of characterization instruments, so that any experimental data are as effective as possible for model calibration. For the numerical case study presented in Sections 7 and 8, we focus on top-down microscopy characterization of Di-BCP thin films, and study the effect of (i) various degrees of image blurring and/or noise intensities, and (ii) various levels of magnification on the calibration of the Ohta-Kawasaki (OK) model [21,34,35], a DFT model, for Di-BCP self-assembly.

Likelihood-free inference and expected information gain computation via measure transport
To solve the proposed Bayesian model calibration problem, we must first tackle key challenges that may seem insurmountable, especially if addressed via conventional likelihood-based methods for Bayesian inference. Specifically, analytical evaluations of integrated likelihood functions [36] are deemed intractable due to (i) the high-dimensional auxiliary random variables, and (ii) the potentially complex forms of sensor noise present in the characterization data. To address these challenges, we adopt the measure transport approach to likelihood-free inference (LFI) that was considered in [37,38]. This approach addresses the challenges of intractable likelihoods, given that it only requires samples from the joint distribution of the parameters and model-generated data. Moreover, the generation of these samples can be performed offline in a scalable manner.
Performing LFI with high-dimensional image data can be challenging, however, as the regions of the parameter space with large likelihood values are typically highly concentrated. Therefore, we advocate constructing low-dimensional summary statistics of the images to capture information in the data that is relevant to the model parameters of interest. This approach is commonly used for the analysis of the image data for BCP self-assembly [39][40][41]. Furthermore, using summary statistics is natural for the model calibration of strongly guided Di-BCP thin film self-assembly using X-ray scattering image data [24,42], as these scattering images are inherently low-dimensional. In our numerical case study in Section 7, we propose several energy-based and Fourier-based summary statistics of top-down microscope images of Di-BCP thin film self-assembly, and their effectiveness for model calibration is quantitatively studied.
The ability to produce posterior densities for Bayesian calibration of BCP models is a prerequisite for the more computationally difficult task of computing the expected information gain (EIG) from summary statistics of image data under various forms of corruptions introduced by the characterization process. To do so, we present a novel likelihood-free estimator for EIG based on measure transport estimators for the posterior distribution. While several EIG estimators based on variational approximations have been recently proposed [43,44], the estimators we use in this work possess nice theoretical properties, such as having the potential of being consistent in the infinite sample setting. Furthermore, we show that the estimation of the EIG can be computed using the proposed LFI approach with no significant additional cost.

Related work
Parameter estimation of models for self-assembly based on experimental characterization of the equilibrium patterns has been previously studied. Khaira et. al [45] determined material and geometric parameters in the TICG model for chemical substrate guided self-assembly of lamellae phase Di-BCP thin films. The parameter estimation is based on data obtained from resonant critical-dimension small-angle X-ray scattering (res-CDSAXS) and grazing-incidence small-angle X-ray scattering (GISAXS) characterization data. Hannon et. al [42] accomplished similar calibration tasks for the SCFT model with res-CDSAXS data.
The parameter estimation in previous works generally follow a deterministic approach. Since deterministic approaches do not quantify the impact of various intrinsic uncertainties on calibration confidence, the reliability of calibrated models may be poor. In particular, by calibrating only to characterization data of Di-BCP equilibrium states strongly influenced by DSA external guidance, the aleatoric uncertainty can be reasonably omitted, since it has a negligible effect on the data. However, in the scenario where the selfassembly is unguided or weakly guided, the aleatoric uncertainty cannot be neglected, as model predictions and experimental characterizations can differ in nontrivial ways due to the long-range disorder or abundance of defected states.
We present related works in LFI and estimation of EIG in Sections 4.3 and 6.1, respectively.

Paper layout
The paper is organized as follows. In Section 2, the relevant components of probabilistic data models for BCP self-assembly are defined. These components include the forward operator of computer simulations and the statistical model for experimental characterization. In Section 3, we introduce the Bayesian framework for the calibration of BCP models and provide a short review of relevant concepts from information theory. In Section 4, we outline the challenges in the Bayesian calibration of BCP models and discuss how these challenges can be addressed with the LFI approach and constructing summary statistics of image data. In Section 5, we present a measure transport approach to LFI for estimating the posterior distribution. In Section 6, we introduce a methodology for estimating EIG via triangular transport maps. In Section 7, we describe the data models for the numerical case study: Bayesian calibration of the OK model for the self-assembly of Di-BCP thin films with top-down microscopy image data. Several energy-and Fourier-based summary statistics of the images are proposed. The prior distribution is defined via a careful analysis of an admissible parameter regime for the OK model. Lastly, in Section 8 we present numerical results for the case study obtained through the proposed LFI and EIG estimation methods.

Contributions
The contributions of this work are summarized, in order of appearance, as follows.
1. We introduce a Bayesian approach for calibration of a BCP self-assembly model that accounts for both measurement and aleatoric uncertainties. 2. We overview a quantitative approach for understanding the informativeness of experimental characterizations about parameters of BCP self-assembly models. Reprinted (adapted) with permission from [46]. Copyright ©2007 American Chemical Society. (right) A schematic drawing of res-CDSAXS experimental setup and a reconstructed scattering image of PS-b-PMMA film strongly guided by a chemically nanopatterned substrate [40].
3. We highlight the computational challenges of integrated likelihood evaluations for calibrating a BCP self-assembly model and demonstrate that a LFI approach is suitable for overcoming such challenges. 4. We overview one methodology for LFI, based on measure transport, that is well-suited for inference problems with generative models, given that it only requires a limited number of offline realizations of parameters and model evaluations. 5. We introduce the notion of extracting summary statistics from BCP characterization image data and propose an efficient approach to measure the utility of selected summary statistics in the LFI setting via a noval EIG estimator that is, in principle, consistent. 6. We present a numerical case study on the Bayesian calibration and EIG computation of the OK model for Di-BCP thin film self-assembly with top-down microscopy characterization.

A probabilistic data model of block copolymer self-assembly
In this section, we conceptualize and formulate generic probabilistic data models of BCP self-assembly that consist of (i) statistical models of experimental characterization, introduced in Section 2.1, and (ii) statistical models of BCP self-assembly, introduced in Section 2.2. The former account for measurement uncertainty in the characterization data, while the latter account for the aleatoric uncertainty. A particular example of the data model that is used for a numerical case study is provided in Section 7.

Statistical models of experimental characterization
We consider the process that generates an image, y ∈ R y , that characterizes a latent BCP melt u ∈ V u (D) in a physical domain D ⊂ R 3 . The material state u may be described by discrete spatial positions of monomers or coarse-grained monomer segments. On the meso-or continuum scale, the material state can also be modeled as the local fractions or normalized densities fields of distinct monomers. We denote the set of all possible Di-BCP material states as V u (D).
Let the characterization process refer to the map from u → y. We model this process as the composition of a state-to-observable map J and a statistical sensor operator S, respectively representing signal formation and image formation in the characterization process. That is, where n represents the stochastic noise introduced during the characterization process.

6
The state-to-observable map can be further decomposed into a restriction operator R and an imaging operator I, i.e., J := I • R. The restriction operator 4 R : V u (D) → V u R(D)) extracts the observed portion of the material state in the subdomain R(D) ⊂ D. For microscopy characterization, R(D) may represent regions of top or cross-section surfaces in the field of view of the microscope. For X-ray scattering characterization, it may represent the subdomain around the incident X-ray beam that affects its trajectory.
The operator I : V u R(D) → V o models the imaging mechanism, which depends on the experimental setup and certain properties of the latent BCP material, such as the refractive indices of the distinct monomers that are relevant for X-ray scattering imaging [48]. Most microscopy imaging methods are realspace methods, where the imaging operator is often constructed using identity maps. One may also include effects related to image resolution, contrast, brightness, depth of field, etc., in the microscopy imaging operator [49]. For X-ray scattering characterization, one may additionally incorporate transformations to the reciprocal space (e.g. via Fourier transform) [39,[50][51][52] in the imaging operator for reproducing scattering patterns from given material states.
The statistical sensor operator S : V o × R n → R y models the discretization and noise corruptions introduced by image-forming sensors. The sensor noise n ∈ R n is a realization of the random variable N , often attributed to infinitesimal fluctuations in sensor voltage. A simple additive model with independent noise is often considered: where the components of the stochastic noise are independent and identically distributed (i.i.d.), and D : V o → R y represents the local averaging of signals introduced by image-forming sensors. Let us note that realistic sensor noise models also account for more complex spatial and signal intensity dependence [53][54][55].
2.2. The aleatoric uncertainty and statistical models of block copolymer self-assembly While the self-assembly of block copolymers leads to the formation of ordered periodic structures, such as those depicted in the experimental phase diagram in Fig. 1 (left), they are typically not realizable in experiments without strong external guidance, even though they are often considered to be the most stable configurations. Instead, a rich variety of states that exhibit long-range disorder, often referred to as defected states or metasable states are observed; see Fig. 1 (right).
To account for the impact of aleatoric uncertainty on the experimentally observed equilibrium states of BCP self-assembly, we assume that the various states attainable under the same experimental scenario S E to be realizations of a V u -valued random variable, U (S E ). We refer to U (S E ) as the state random variable. The state random variable gives rise to a probability distribution of the achievable states, whose structure depends on the material, geometrical, thermal, and various other controllable macroscopic factors specified by S E . When the self-assembly is strongly guided, the distribution of the state random variable is close to a Dirac mass, meaning that only a handful of BCP equilibrium states are attainable. When the self-assembly is unguided or weakly guided, the distribution can be quite complex, yet it typically concentrates on states in the same phase or mixture of phases that have similar structures, such as similar radius and periodicity length of spheres in the sphere phase.
A statistical model of BCP self-assembly can be generalized as a mathematical model derived from statistical polymer theory for specifying the probability distribution, i.e., the statistical ensemble, of U (S E ) for a range of experimental scenarios. It typically contains a model of polymer chains, models of the inter-polymer and polymer-environment interactions, and possibly a model of thermal fluctuations and/or dynamics of polymers during the self-assembly. A comprehensive review of these models can be found in [29,Vol. 1], [19,31,56], and references therein.
Computer simulations of BCP self-assembly can be treated as stochastic simulators for generating independent samples from model-specific statistical ensembles of equilibrium states. Let us denote the forward operator of a statistical model of BCP self-assembly as F : R x × V z → V u . This model assumes the equality in distribution where x ∈ R x are parameters 5 that represents the material, geometrical, thermal, and/or other macroscopic factors of the experimental scenario S E that are considered in the model. The V z -valued random variable Z represents the stochastic component of the forward operator. The origin of the stochastic component Z may be purely numerical. For example, Z can partly represent random initial configurations of the polymeric system and the sampling mechanism for accepting perturbed configurations in Monte Carlo simulations of equilibrium BCP melts [58]. It can also be high-dimensional or even infinite-dimensional. For example, Z can be specified by the stochastic forcing term that resides in infinite-dimensional function spaces for dynamical DFT models [20]. We refer to Z as an auxiliary random variable, or nuisance variable, as it is either non-physical or represents uncontrollable microscopic factors in experiments that are difficult to observe. For the numerical case study in Section 7 and 8, we use the OK model, in which the stochastic component is specified by random initial states for free energy minimization.

A generic data model of block copolymer self-assembly
We call the combination of the model of BCP self-assembly and model of experimental characterizations of equilibrium states as the data model. For notational convenience, we first define the parameter-to-observable map O := J • F, where O(X, Z) is the observable random variable. For example, for X-ray scattering characterization, O(x, Z) may be the (noise-free) model-predicted latent scattering pattern where variability comes from the model-predicted random BCP equilibrium structures. We define the image data as a R yvalued random variable arising from both the measurement uncertainty and the aleatoric uncertainty. The data model G : R x × V z × R n → R y summarizes the relationship between model parameters and data as

Bayesian inference and quantifying uncertainty
In this section, we present the Bayesian inference framework for calibrating models for BCP self-assembly. In Section 3.1 we introduce the posterior distributions that arise from the likelihood and the, so called, conditional likelihood. To summarize the uncertainty in calibrated parameters, we introduce some informationtheoretic quantities in Section 3.2. Lastly, in Section 3.3, we present the expected information gain as a natural quantity to consider in a Bayesian inference context for measuring the information contained in the data about the unknown model parameters.

The posterior and joint posterior
For a Bayesian formulation to model calibration, we introduce a random variable for the parameters X ∈ R x and seek to characterize its distribution given an observation of Y . Initially, the parameter X is distributed according to the prior measure ν X ∈ P(R x ) with density π X . The prior measure represents our epistemic uncertainty in the model parameters before collecting observations. Similarly, the (possibly infinite-dimensional) auxiliary random variable is distributed according to the conditional prior measure ν Z|X that may depend on the model parameter. In our numerical case study, a prior distribution for the OK model parameters is constructed in Section 7.5 and the distribution for the auxiliary variable is described in Appendix A.
The data model in Section 2.3 defines the conditional distribution for the observations. Let Y = G(X, Z, N ) be a data random variable and assume the sensor noise random variable N ∼ π N is independent of (X, Z). This noise yields a conditional density for the data given the inputs x and z that we denote by π Y |X,Z (·|x, z). For example, if the sensor operator in the data model (4) is given by the additive form in (2), then this conditional density is given by π N (· − D(O(x, z))).
For the BCP model calibration problem, however, we are not interested in inferring the auxiliary variable Z. Thus, we seek to describe the conditional distribution of the data given only the model parameters. The uncertainties in the sensor noise and the auxiliary variable define a marginal conditional distribution for Y given X, whose density we denote by π Y |X (·|x). For these conditional distributions of the data, we define both the conditional likelihood and the likelihood for an observation y * as the functions which arise when fixing the data to y * in the conditional densities. That is, π y * |X,Z (x, z) := π Y |X,Z (y * |x, z) is the conditional likelihood for both the model parameters and the auxiliary variable, while π y * |X (x) := π Y |X (y * |x) is the likelihood for only the model parameters. The two likelihood functions represent, respectively, the likelihood of observing an image y * for (i) a particular model-predicted equilibrium state generated at (x, z), and (ii) a variety of random model-predicted equilibrium states generated at x. We will discuss how both are defined in depth in Section 4.
Given a prior and likelihood function for an observation y * , Bayes' theorem provides a rule for producing the conditional distributions of the parameters. The joint posterior density for the parameters and the auxiliary variable is given by 6 where C y * = π Y (y * ) = R x V z π y * |X,Z (x, z)dν Z,X (z|x)π X (x)dx > 0 is the normalization constant for the posterior density and is often referred to as the model evidence or the marginal likelihood. To calibrate the model parameters, our interest is in the marginal posterior density of the parameters alone given by Many computational tasks in Bayesian inference rely on evaluating these densities to draw samples from the posterior distributions. While evaluations of (5) depend on evaluations of the conditional likelihood π y * |X,Z , evaluations of (6) inherit difficulties from evaluating the likelihood π y * |X that stem from the influence of the nuisance variable Z on the observations. We discuss these computational difficulties in Section 4. For a thorough treatment on these challenges in Bayesian inference, in general, see [59,60].

Entropy and information gain
To understand and quantify the effectiveness of Bayesian model calibration, we turn to information theory that formally quantifies the information content in random variables. One information-theoretic measure for the uncertainty of a parameter X is the Shannon entropy H : P(R x ) → R given by The entropy is larger for a parameter that is less certain. For example, if a parameter is known, i.e., its distribution is a Dirac mass, then its entropy is zero. After conditioning on the data y * , we can measure the entropy of the parameters by computing the conditional entropy H(X|Y = y * ). If the entropy decreases after conditioning, i.e., H(X|Y = y * ) < H(X), we say that observing the data y * has decreased our uncertainty or increased our confidence in the values of model parameters. The expectation of H(X|Y = y) over the distribution of the data random variable defines the conditional entropy of the parameters given the data H(X|Y ) and it always satisfies H(X|Y ) ≤ H(X). If the conditional entropy is less than the prior's entropy, we say that the data is, on average, informative of the parameters.
To quantify the information gained from performing Bayesian inference, we measure the change in the entropy of the parameters X before and after conditioning on data. This change quantifies the value of acquiring certain noisy, corrupted, and/or partial observations of the system to learn about the unknown model parameters. In the context of Bayesian inference, this change compares the prior density π X and the posterior density π X|y * for the parameters. An information-theoretic quantity for measuring information gain based on the change in entropy is the Kullback-Leibler (KL) divergence D KL . The KL divergence from the prior to the posterior is given by The KL divergence is always non-negative and is zero if and only if the posterior is equal to the prior almost everywhere. Furthermore, a large value for the KL divergence implies that the posterior has lower entropy relative to the prior and its density is typically more concentrated than the prior density. A large departure of the posterior from the prior indicates that we gained more information from the observations about the parameters.

Expected information gain
Before observing any particular data, we can measure the information gained about the parameter in a Bayesian inference problem by considering the expected information gain (EIG) I : P(R x ) × P(R y ) → R + between the random variables X and Y . We say that I(X; Y ) represents the the EIG for learning about X from Y . The EIG can be defined by the integration of the KL divergence in (8) with respect the data y. That is, Within the field of Bayesian optimal experimental design (BOED), the EIG is commonly used to assess the utility of an experimental measurement for reducing uncertainty in the model parameters [61]. In this context, the EIG is typically a function of the experimental design variables, and is maximized to find the optimal design. An example of BOED is the optimal design of a discrete observation operator, e.g., choosing sensor locations that maximize the EIG on certain quantities of interest [62]. In the context of microscopy characterization of BCP melts, for example, we can measure I(X; Y ) for different restriction operators R, introduced in Section 2.1, to determine the optimal level of magnification for inferring a particular subset of the model parameters. In our numerical case study in Section 7 and 8, we will demonstrate results on the computation of EIG that enables such a type of BOED study. We note that the expression for the EIG I(X; Y ) in (9) is equivalent to the mutual information (MI) between random variables X and Y . The MI is always a non-negative random variable. While a zero value implies that the two random variables are independent, a larger value indicates that the data is more informative on average for inferring the parameter. Using the form of the KL divergence in (8) and Bayes' rule in (6), the EIG can also be written in the two equivalent ways While the first expression uses the ratio of the posterior and prior densities, the second expression uses the ratio of the conditional likelihood and the marginal density for the observations. We note that as compared to standard Bayesian inference problems that characterize the posterior for a single realization of the data, computing the EIG requires evaluating the posterior densities for many possible realizations of the data. In Section 6, we propose an efficient estimator for the EIG in (10a) that is based on a variational approximation for all of the posterior densities.

Addressing the challenges of auxiliary variables and high-dimensional image data via likelihoodfree inference
Performing Bayesian calibration and computing EIG requires characterizing posteriors, i.e., producing representative samples from the the posterior distribution or providing an exact or approximate expression of the density π X|y * . In this section, we outline distinctive challenges for the Bayesian calibration of BCP self-assembly models. In Section 4.1 we describe intractable integrated likelihood functions that arise when accounting for the aleatoric uncertainty in the data model. Sections 4.2 and 4.3 discuss computational methods for dealing with integrated likelihoods, including LFI methods. Lastly, in Section 4.4 we address inference with high-dimensional image data by presenting an approach that constructs low-dimensional summary statistics of the data to capture relevant information about the model parameters of interest.

Evaluating integrated likelihoods
To calibrate the model parameters x given image data y * , it is necessary to prescribe the likelihood function π y * |X in (6). This likelihood is derived from the stochastic procedure for generating the image data using the data model for BCP self-assembly in (4). The procedure depends on, not only the sensor noise N , but also the auxiliary variable Z entering the parameter-to-observable map. In particular, the likelihood of an observed image y * ∈ R y given model parameters x accounts for all realizations of the auxiliary variable and stochastic noise that could possibly generate the particular image. Mathematically, this is achieved by summing the contributions of the conditional likelihoods π y * |X,Z (x, z) weighted by the probability of generating each auxiliary variable z. Hence, the likelihood function for the parameters is defined as the following integrated likelihood [36] that marginalizes out the auxiliary variables Using this likelihood leads to the posterior density for the parameters x in (6) that is a low-dimensional marginal of the joint posterior distribution for (X, Z) 7 . The integrated likelihood in (11) has a closed-form expression only in very special situations. One setting is when both the conditional distribution for the auxiliary variables ν Z|X (·|x) and the conditional density for the observations π Y |X,Z (·|x, z) are Gaussian, and the parameter-to-observable map is affine. In general, however, the integral we consider over the infinite-dimensional auxiliary variables z is unavailable in closed form 8 . In these cases, we say that the integrated likelihood is unavailable and consider it intractable to evaluate.
Despite their intractability, integrated likelihoods are frequently encountered in many fields. Other interesting inference problems that involve integrated likelihoods include time, phase, and distance marginalization in gravitational wave inference [64], boundary condition marginalization for linear parabolic PDEs [65], and hyper-parameter discovery for hierarchical models in image deblurring and X-ray computational tomography [66]. Therefore, developing efficient computational inference methods that work with integrated likelihoods is an important ongoing area of research.

Inference methods for intractable likelihoods
In settings with intractable likelihoods, it is not possible to compute posteriors using likelihood-based sampling algorithms, such Markov chain Monte Carlo (MCMC), or other variational methods which require evaluating the likelihood at multiple values of the parameters. To satisfy this constraint, three classes of inference methods have been developed that do not require explicit likelihood evaluations. These classes include surrogate likelihood functions, pseudo-marginal Monte Carlo methods [67,68], which assume tractable conditional likelihood function evaluations, and likelihood-free inference, often referred to as simulation-based inference [69].
The profile likelihood method [70] is one popular choice of surrogate that replaces integration in (11) with a maximization over z. Specifically, the profile likelihood π y * |X is defined as the conditional likelihood evaluated at the maximum likelihood estimate (MLE) for z as a function of x, i.e., z MLE (x) = arg max z∈V z π y * |X,Z (x, z). Thus, if a well-defined unique z MLE exists, the profile likelihood can be tractably evaluated as π y * |X (x) := π y * |X,Z (x, z MLE (x)). This approach, however, replaces z by a point estimate and ignores the impact of its uncertainty on the likelihood. As a result, the profile likelihood is inconsistent and often requires adding factors to correctly characterize the posterior uncertainty [71].
If conditional likelihood evaluations are available, pseudo-marginal methods can perform inference using Monte Carlo estimates of the integrated likelihood function in (11). These approaches generate N samples of the auxiliary variable If the variance Var(π y * |X,Z (x, ·)) is finite, the standard deviation of this estimator decreases with more samples at the rate O(1/ √ N ). Thus, we can reduce the sample variance of π y * |X by increasing the number of samples until sufficient accuracy of the integration has been reached. While unbiased Monte Carlo estimators of the integrated likelihood have been shown to yield consistent posterior sampling [67], these methods may necessitate a large number of samples to achieve a reasonable computational efficiency in practice [72]. On the other hand, the evaluation of the conditional likelihood function at each sample of the auxiliary variable requires evaluating the forward operator, which can be computationally expensive. Furthermore, estimates of the integrated likelihood are not shared for different parameter values and need to be recomputed at each iteration of some sampling methods, such as MCMC.
In addition to the challenge of marginalization, the conditional likelihood π y * |X,Z can also become analytically intractable to evaluate for many non-trivial sensor noise models, S(·, N ), which rules out the use of pseudo-marginal methods and the profile likelihood. The analytical expression for the conditional likelihood stems from a change of variables between the noise random variable and the data random variable. For example, the sensor operator with multiplicative noise, i.e., S(·, N ) = D(·)N , with noise N ∼ π N has the conditional likelihood In general, evaluating the conditional likelihood requires computing the inverse and Jacobian inverse of the sensor noise model, which is often unfeasible when it has a complicated form 9 . When the data arises from a "black-box" simulator, i.e., when the data model G relies on legacy scientific code whose mathematical expressions may be unavailable, it is not possible to derive the conditional likelihood analytically. In polymer modeling and characterization, complex mathematical expressions and scientific codes underlay both the forward operators F and realistic sensor operators S. Therefore, we consider the evaluation of conditional likelihoods to be generally unfeasible.
On the other hand, likelihood-free inference (LFI) methods enable statistical inference in generative models. They rely only on sampling from the joint distribution of model parameters and model-generated data. In the context of the integrated likelihood function in (11), this is accomplished by 1. sampling the model parameters, : R n → R y is presumed to be a generally non-monotonic differentiable transformation, and π Y |X,Z (y|x, z) is a weighted sum of densities over a partition of the input space A k ⊂ R n into regions for which Sx,z is monotonic. Unless this reduces to a simpler closed-form expression with explicit computations, the evaluation of π Y |X,Z (y|x, z) can be quite complex for arbitrary S, due to the difficulties of evaluating generalized inverses S −1 x,z , computing Jacobian determinants (generally O(n 3 ) computational complexity, where n 1 for our sensor noise), and determining the partitions {A k }.
4. sampling the sensor noise n i ∼ π N , and 5. evaluating the statistical sensor model, The resulting samples (x i , y i ), with z i and n i discarded, are samples from the marginal distribution of (X, Y ). This collection of marginal samples avoids the step of explicitly marginalizing the conditional likelihood π Y |X,Z to evaluate the integrated likelihood in (11). Given a set of i.i.d. joint samples (x i , y i ) that is large enough to be representative of their joint distribution, the goal of LFI methods is to approximate the posterior density π X|y * or to sample from it, without evaluating the likelihood function. Thus, LFI is agnostic to the form of the likelihood and applicable for inference with various sources of intractability, including the high-dimensional integration with respect to the auxiliary variable. This makes LFI appealing for the BCP model calibration problem.

Likelihood-free inference: Related work and practical desiderata
The family of LFI methods for posterior sampling includes both classical techniques such as approximate Bayesian computation (ABC) and synthetic likelihoods, as well as modern machine learning methods. ABC methods [73][74][75] define a non-parametric approximation to the likelihood function, and thereby the posterior, using an accept/reject sampling scheme. For each pair of joint samples (x i , y i ), this scheme accepts parameters whose simulated data y i is close to the observed data y * with respect to a discrepancy function d and tolerance level , i.e., ABC accepts x i if d(y i , y * ) < . While ABC can be shown to be a consistent approximation of the posterior [76] in the limit as → 0, small values for lead to rejection of many candidate parameters, and often an impractical number of simulations is required to accurately estimate posterior statistics. Furthermore, the quality of the non-parametric approximation in ABC deteriorates as the data dimension increases [77].
An alternative to ABC that is biased, but typically more efficient, is a parametric approximation to the likelihood function. As an example, the Bayesian synthetic likelihood (BSL) method uses samples y i ∼ ν Y |X (·|x) to create a multivariate Gaussian approximation to the likelihood at each parameter value x given by π Y |X = N (y(x), C(x)) where y(x) and C(x) denote the estimated mean and covariance matrix. The approximate likelihood evaluations π y * |x (x) = π Y |x (y * ) can then be efficiently used within MCMC algorithms [78,79]. While BSL can be more efficient than ABC for high-dimensional data, this approach requires sampling multiple data, i.e., evaluating the data model, for each parameter in order to determine y(x) and C(x). As a result, both ABC and synthetic likelihood approaches are less suitable for applications that involve computationally expensive forward operators. Lastly, the inference results of both methods are applicable only to the particular set of observed data y * and thus require repeating the computation for any new set of data.
To amortize the cost of performing inference across many sets of observed data, recent machine learning methods propose to learn models that are functions of both the parameters and the data [69]. These methods use the collection of joint samples (x i , y i ) to learn variational approximations for either the posterior density π X|Y [37,80], the likelihood function π Y |X [81], or the likelihood ratio π X|Y /π X [82]. The resulting models are then used in an inference step to define a tractable approximation to the posterior density or to generate approximate posterior samples for any observed data y * . This enables the inference results to be used for downstream tasks that require access to the posteriors for all data realizations, such as EIG computation and BOED. To avoid constraining the structural forms for these models, many techniques use flexible neural network architectures that can capture complex posterior distributions and can learn the structure of highdimensional and domain-specific data (e.g., using architectures designed for images). Many of these models, however, often require a large upfront cost to learn them reliably and adapting these methods to small-sample settings remains challenging. To select a LFI method for the calibration of BCP self-assembly models, we have the following desiderata: 1. Sample efficiency: Method provides reasonable posterior approximations in settings where it is computationally expensive to generate simulated data. 2. Consistency: Method does not impose arbitrary approximations on the form of the likelihood or the posterior distribution in order to capture non-Gaussianity without biases from the inference procedure.
3. Amortized inference: Method provides posterior densities π X|y * or samples for almost all values of the data without needing to repeat inference for each realization of y * .
To produce consistent posterior approximations, without discarding or rejecting samples, we consider ABC and BSL to be less appealing for our application. To also benefit from amortization similarly to the machine learning methods, we propose to use a measure transport approach for LFI in this work, which relies on a careful parameterization for settings with limited sample sizes. The details of the proposed approach are discussed in Section 5.

Informative summary statistics of high-dimensional image data
When performing LFI with high-dimensional data, it is common to construct and work with lowdimensional summary statistics of the data [83], Q = Q(Y ) ∈ R q . Working with this type of summary statistics has two inherent advantages: (i) it can improve the quality and tractability of LFI, and (ii) it can provide insight into what information contained in the data is relevant for inference.
In principle, arbitrarily high-dimensional data could be used to perform inference. This data, however, can lead to concentrated regions of the parameter space with large likelihood values, which hampers the effectiveness of many LFI methods. Here we provide two numerical examples to visualize this phenomenon.
For both examples, we consider the calibration of the OK model for Di-BCP thin film self-assembly using top-down microscopy image data. We consider the characterization model with additive Gaussian noise N ∼ N (0, C N ) as described in (2). That is where the state-to-observable map is taken as the identity resulting in O ≡ F. The OK model, its associated forward operator F, and sensor signal discretization D are defined in Section 7. Here we use the squared noise-covariance-weighted l 2 data misfit 10 M f defined as to compare the synthetic data y * that is generated from the fixed latent inputs (x * , z * ) via (12), with model-predicted patterns at other inputs (x, z). We note that the data misfit is proportional to the negative log-likelihood generated by the sensor noise in (12). In the first numerical example, the data misfit is evaluated at x * and 10 3 i.i.d. samples of the auxiliary variable, z i ∼ ν Z|X (·|x * ). In Fig. 3, a histogram of the squared l 2 data misfit evaluations is shown along with the synthetic data y * , the most likely observable corresponding to the lowest l 2 misfit, and the least likely observable corresponding to the highest l 2 misfit. In the second numerical example, the data misfit is evaluated at z * and 10 3 i.i.d. samples of the model parameter x i ∼ N (x * , δ 2 Id x ), where δ is chosen to be small 11 . Similar plots for the second example are shown in Fig. 4.
In both figures, we observe that, while the true parameters reproduce observables identical to the image up to the additive noise, most realizations of the nearby parameters or auxiliary variables produce observables that are drastically different from the image data according to the l 2 data misfit. As a result, the misfit values for most of the samples are very large, relative to misfit values for the latent input that generated the synthetic data. This leads to a likelihood concentration phenomenon that results in LFI methods requiring many joint samples to recover the true parameters.
Furthermore, similar features in characterization images y * of BCP self-assembly indicate that certain characteristics of their latent BCP equilibrium structures are similar. This should, in general, imply certain closeness in their underlying model parameters. However, the likelihood function π y * |X (x) measures the . Likelihood-free inference via triangular transport maps In this section, we present a measure transport approach for characterizing the family of conditional ensities for the parameters given the observations ⇡ X|Q using only samples from the joint distribution ⌫ X,Q . n 5.1 and 5.2 we introduce the notion of measure transport and triangular transport maps, respectively. n 5.3 we show how to estimate the maps from samples. Then, in 5.4 we demonstrate how these maps can e used for likelihood-free inference. A numerical implementation of this approach will be available in an pcoming release of the TransportMaps Python package 10 , a software library implementing algorithms for ransport map-based variational inference.

.1. Measure transport
Let W 2 R w be a target random variable with measure ⌫ W 2 P(R w ) we would like to characterize, nd let R 2 R w be a reference random variable with measure ⌫ R 2 P(R w ) that we can easily sample from, or example, a standard Gaussian. Our goal is to find a measurable function S : R w ! R w that defines a eterministic coupling between these two variables such that the law of S(W ) is equal to the law of R. For uch a map S, we say that S pushes forward ⌫ W to ⌫ R and we denote this by S ] ⌫ W = ⌫ R . Equivalently, we ay that S pulls back ⌫ R to ⌫ W and we denote this by S ] ⌫ R = ⌫ W . For measures ⌫ W and ⌫ R with densities W and ⇡ R , respectively, similar relationships apply, S ] ⇡ W = ⇡ R and S ] ⇡ R = ⇡ W . If S is a diffeomorphism i.e., a differentiable bijection with a differentiable inverse), we can evaluate the pullback density S ] ⌫ ⇢ using he change-of-variables formula as herefore, having the transport map S that couples W and R and the reference density ⇡ R enables us to valuate the target density ⇡ W . Furthermore, we can use this map to generate i.i.d. samples from the target istribution by sampling from the reference measure r i ⇠ ⌫ R and applying the inverse map to generate amples S 1 (r i ) ⇠ ⌫ W . Figure 5 provides an example of using S to characterize the pullback density and enerate samples from a 2-dimensional "banana" target distribution. The existence of a map S that pushes forward ⇡ W to ⇡ R is guaranteed under very mild conditions such s ⇡ W and ⇡ R being absolutely continuous with respect to the Lebesgue measure, i.e., the distributions do ot contain atoms. We refer the reader to [? ] for how to relax these conditions.

.2. The Knothe-Rosenblatt rearrangement
In general, many maps S can satisfy the measure transport conditions above [? ]. For instance, the ptimal transport literature looks for maps S that couple two distributions and also minimize an integrated ransportation distance, such as the Wasserstein metric [? ]. In this work, we consider the Knothe-Rosenblatt KR) rearrangement for S, the unique monotone and lower triangular map that pushes forward ⇡ W to ⇡ R [? ? ]. A monotone and lower triangular map is a multivariate function of the form 10 transportmaps.mit.edu

14
The synthetic data y * The most likely observable The least likely observable

Likelihood-free inference via triangular transport maps
In this section, we present a measure transport approach for characterizing the family of conditional sities for the parameters given the observations ⇡ X|Q using only samples from the joint distribution ⌫ X,Q . .1 and 5.2 we introduce the notion of measure transport and triangular transport maps, respectively. .3 we show how to estimate the maps from samples. Then, in 5.4 we demonstrate how these maps can used for likelihood-free inference. A numerical implementation of this approach will be available in an oming release of the TransportMaps Python package 10 , a software library implementing algorithms for sport map-based variational inference.

Measure transport
Let W 2 R w be a target random variable with measure ⌫ W 2 P(R w ) we would like to characterize, let R 2 R w be a reference random variable with measure ⌫ R 2 P(R w ) that we can easily sample from, example, a standard Gaussian. Our goal is to find a measurable function S : R w ! R w that defines a erministic coupling between these two variables such that the law of S(W ) is equal to the law of R. For h a map S, we say that S pushes forward ⌫ W to ⌫ R and we denote this by S ] ⌫ W = ⌫ R . Equivalently, we that S pulls back ⌫ R to ⌫ W and we denote this by S ] ⌫ R = ⌫ W . For measures ⌫ W and ⌫ R with densities and ⇡ R , respectively, similar relationships apply, S ] ⇡ W = ⇡ R and S ] ⇡ R = ⇡ W . If S is a diffeomorphism ., a differentiable bijection with a differentiable inverse), we can evaluate the pullback density S ] ⌫ ⇢ using change-of-variables formula as refore, having the transport map S that couples W and R and the reference density ⇡ R enables us to luate the target density ⇡ W . Furthermore, we can use this map to generate i.i.d. samples from the target ribution by sampling from the reference measure r i ⇠ ⌫ R and applying the inverse map to generate ples S 1 (r i ) ⇠ ⌫ W . Figure 5 provides an example of using S to characterize the pullback density and erate samples from a 2-dimensional "banana" target distribution. The existence of a map S that pushes forward ⇡ W to ⇡ R is guaranteed under very mild conditions such W and ⇡ R being absolutely continuous with respect to the Lebesgue measure, i.e., the distributions do contain atoms. We refer the reader to [? ] for how to relax these conditions.
The synthetic data y * The most likely observable The least likely observable probability of the model parameterized by x matching a specific image. As a result, the observables that are visually similar to the image data, but are mismatched in entry-wise values, will not be weighted heavily by the likelihood function, as demonstrated through the two numerical examples above. This leads to unnecessarily ineffective parameter inference. Although using entire images for inference ensures that all of the available information about the latent state is preserved, some information may be irrelevant and/or redundant for the purpose of model calibration and prediction. For example, the specific orientation of phase interfaces in a top-down microscope image of unguided lamella phase BCP thin films (i.e. the 'figer-print' patterns) are irrelevant (when ignoring boundary effects) if our interest is to predict the natural periodicity length of the BCP material. In other words, compressing irrelevant information in the data should not affect the posterior density for the parameters.
The notion of data reduction is very common for image data. While the l 2 distance can be easily computed, it does not capture the visual similarity between images. For instance, the distance can be arbitrarily large when the two images are identical, up to a rotation or reflection. Therefore, it is often recommended to measure the discrepancy between images in a more natural way. Different metrics that mitigate the issues for high-dimensional images have been considered by many authors [84][85][86] and some image processing techniques have been applied to image data for BCP self-assembly in the past [26,39,41]. These metrics and the corresponding likelihood functions are designed to be relatively insensitive to certain transformations that do not change the fundamental characteristics of the image related to the model parameters of interest.
In this work we advocate for constructing low-dimensional summary statistics of BCP characterization images that quantitatively describe image features that are related to characteristics of their latent BCP equilibrium structures. Ideally, the summary statistics Q are chosen to be sufficient statistics. This means that they do not introduce any bias in the inference procedure, i.e., they satisfy π X|Y (x|y) = π X|Q (x|Q(y)) for all x, y. Furthermore, we can also look for minimal sufficient statistics, i.e., the summary statistics that most reduce the dimension of the data without discarding any relevant information. Unfortunately, sufficient statistics are not readily available outside of exponential family models, and consequently many methods have been proposed to select summary (but not necessarily sufficient) statistics for LFI in practice [87]. These summary statistics may discard some information during the inference, as a trade-off for suffering less from the concentration phenomenon of the likelihood for high-dimensional images [88].
It is important to note that working with summary statistics of the data remains feasible with LFI methods. While the conditional distributions for the summary statistics π Q|X (·|x) may not be analytically available or feasible to evaluate, LFI methods only require joint samples of the parameters x i and summary statistics q i . These samples can be acquired by an additional step of computing the summary statistics q i = Q(y i ) for each data sample y i .

Likelihood-free inference via triangular transport maps
In this section, we present a measure transport approach for characterizing the family of conditional densities for the parameters given the observations π X|Q using only samples from the joint distribution ν X,Q . In Section 5.1 and 5.2 we introduce the notion of measure transport and triangular transport maps, respectively. In Section 5.3 we show how to estimate the maps from samples. Then, in Section 5.4, we demonstrate how these maps can be used for likelihood-free inference. A numerical implementation of this approach is available online in a MATLAB package for transport map estimation 12 .

Measure transport
Let W ∈ R w be a target random variable with measure ν W ∈ P(R w ) we would like to characterize, and let R ∈ R w be a reference random variable with measure ν R ∈ P(R w ) that we can easily sample from, for example, a standard Gaussian. Our goal is to find a measurable function S : R w → R w that defines a deterministic coupling between these two variables such that the law of S(W ) is equal to the law of R. For such a map S, we say that S pushes forward ν W to ν R and we denote this by S ν W = ν R . Equivalently, we say that S pulls back ν R to ν W and we denote this by S ν R = ν W . For measures ν W and ν R with densities π W and π R , respectively, similar relationships apply, S π W = π R and S π R = π W . If S is a diffeomorphism (i.e., a differentiable bijection with a differentiable inverse), we can evaluate the pullback density S ν R using the change-of-variables formula as Therefore, having the transport map S that couples W and R and the reference density π R enables us to evaluate the target density π W . Furthermore, we can use this map to generate i.i.d. samples from the target Figure 5: The target density π W is represented as the pullback of a standard Gaussian reference density π R through the transport map S (left). Samples from ν R are generated by mapping reference samples through the inverse of S (right).
Similarly we can represent the reference density as the pushforward of the target density through S and sample from ν W by mapping target samples through S.
distribution by sampling from the reference measure r i ∼ ν R and applying the inverse map to generate samples S −1 (r i ) ∼ ν W . In Fig. 5, we provide an example of using S to characterize the pullback density and generate samples from a 2-dimensional "banana" target distribution. The existence of a map S that pushes forward π W to π R is guaranteed under very mild conditions such as ν W and ν R being absolutely continuous with respect to the Lebesgue measure, i.e., the distributions do not contain atoms and their densities exist. We refer the reader to [89] for approaches that relax these conditions.

The Knothe-Rosenblatt rearrangement
In general, many maps S can satisfy the measure transport conditions above [90]. For instance, the optimal transport literature looks for maps S that couple two distributions and also minimize an integrated transportation distance, such as the Wasserstein metric [91]. In this work, we consider the Knothe-Rosenblatt (KR) rearrangement for S, i.e., the unique monotone and lower triangular map that pushes forward π W to π R [92][93][94]. A monotone and lower triangular map is a multivariate function of the form where ξ → S k (w 1 , . . . , w k−1 , ξ) : R → R, i.e., the restriction of the map component S k to its last variable w k , by fixing its first k − 1 variables, is monotone increasing for all (w 1 , . . . , w k−1 ) ∈ R k−1 . We denote the space of triangular and monotone maps on R w by S (R w ).
Working with the KR rearrangement, rather than other (non-triangular) transport maps, has several computational benefits [38]. First, monotone triangular maps S ∈ S (R w ) can be easily inverted to sample from ν W by solving the nonlinear system S(w i ) = r i , r i ∼ ν R via a sequence of 1-dimensional nonlinear solves. Furthermore, each nonlinear solve has a unique solution due to the injectivity of the map. Second, the Jacobian determinant of a triangular map det ∇S(w) is tractable to compute as the product of the partial derivatives of each component with respect to the last variables, i.e., det ∇S(w) = w k=1 ∂ k S k (w 1:k ). Third, these maps are easier to parameterize than maps derived via optimal transport and can be computed quickly using gradient-based optimization techniques [95,96].
Most importantly in our context, these maps provide direct access to the conditional densities of the target distribution. For a reference distribution with independent components, such as the standard Gaussian, the k map component S k pushes forward the kth marginal conditional density π W k |W1,...,W k−1 in the factorization of the target density π W = π W1 π W2|W1 . . . π Ww|W1,...,Ww−1 , to the kth marginal of the reference density π R [97, App. A]. In Section 5.4, we use this property to extract the posterior density, i.e., the conditional density of the parameters given the data, from their joint density.

Estimating the Knothe-Rosenblatt rearrangement
The exact KR rearrangement S KR that couples two densities is, in general, analytically unavailable. However, we consider approximations of the KR map in this work. We follow the variational approach in [38,96,97] of minimizing the KL divergence, 13 from the pullback density S π R to the target density π W , which has the form The KR rearrangement is given by the global minimizer of the KL divergence over the space of lower triangular and monotone maps on R w . That is, In all of our numerical experiments, we let the reference distribution be standard Gaussian. Substituting the form of the reference density π R = N (0, Id w ) in (16), the optimization problem can be written as where each term in the sum only depends on one map component. Thus, optimizing for S can be decoupled into w optimization problems for each map component 14 , which can each be solved independently and in parallel to determine each S k . More information on the parametrization and computation of monotone and triangular maps is provided in Appendix C.

Conditional density estimation via triangular transport maps
To estimate the conditional densities π X|Q and any posterior π X|q * in particular, we begin by considering the KR rearrangement S KR that pushes forward the target joint density π X,Q to a reference density with independent components π R1,R2 = π R1 π R2 , that are both defined on R x+q . Given that S KR is lower triangular, we can partition the map into two blocks as where S Q : R q → R q and S X (q, ·) : R x → R x are monotone functions for all q ∈ R q .From the measure transport condition, then S Q pulls back the marginal π R1 to the marginal for the summary statistics π Q , and S X (q, ·) pulls back the marginal π R2 to the conditional π X|Q (·|q) for each q ∈ R q .
From the decoupling of the optimization problem over the map components, we only need to solve x independent optimization problems for the components in S X to approximate the conditional density π X|Q . Moreover, the conditional density can be estimated given N i.i.d. samples from the joint distribution, We use these samples to define a Monte Carlo approximation of the expectation in (17) and to find the kth component of S X by solving under the monotonicity constraint ∂ k (S X ) k (q, x 1:k ) > 0 for all (q, x 1:k ) ∈ R q+k . Then, collecting the components for k = 1, . . . , x in S X produces an estimator of the conditional density π X|Q (·|q) for any q ∈ R q as An example of the estimated conditional densities π X|Q (·|q) extracted from the transport map S X (q, ·) is shown in Fig. 6. We note that a single map S X that is a function of the summary statistics q can characterize conditional densities with different properties across the data space, such as unimodal and bimodal behavior. In addition to having a functional form for the conditional density π X|Q (·|q), we can also use the estimated map to generate i.i.d. samples x i from the conditional distribution ν X|Q (·|q) by sampling r i from the reference measure and solving for the parameter x i that is mapped to r i by the function S X (q, ·), i.e., To summarize the approach outlined in this section, we use a sample representation (x i , q i ) of the joint distribution ν X,Q to compute an approximate map S X that couples a reference density to the conditional densities π X|Q ; this map enables both evaluation of approximate posterior densities π X|q * for any and all summary statistics q * , via (20), and simulation of i.i.d. samples from their respective distributions ν X|Q (·|q * ), via (21).

Estimating expected information gain via triangular transport maps
A core advantage of estimating the conditional densities for all values of q is that it enables the computation of expectations over the data random variables. One of these quantities is the EIG defined in Section 3.3. An example of using EIG to compare the informativeness of two summary statistics for parameter inference is shown in Fig. 7. In Section 6.1 we describe existing estimators for EIG and in Section 6.2 we present a lower bound for the EIG based on the transport maps computed in Section 5.  Fig. 6 given two summary statistics Q 1 = −X + N (left) and Q 2 = X 2 + N (right) for N ∼ N (0, 0.01). The posteriors for Q 1 are concentrated as compared to the bi-modal posteriors for Q 2 , indicating that Q 1 is more informative for predicting X. Integrating the KL divergence between the prior for X and the posteriors over the different values of Q i results in EIGs of I(X; Q 1 ) = 1.64 ± 0.01 and I(X; Q 2 ) = 0.86 ± 0.01. This provides a quantitative approach to measure and rank the informativeness of different summary statistics.

Related work
Given conditional densities π Q|X that are tractable to evaluate, several nested Monte Carlo (NMC) estimators have been proposed to estimate the EIG using samples from the joint distribution of parameters and observables [98][99][100]. The EIG in (10a) and (10b) contains an expectation over ν X,Q and an unknown marginal density for the summary statistics that is typically unavailable in closed form. To estimate the EIG, NMC uses and an independent set of N 2 i.i.d samples {x j } N2 j=1 ∼ ν X to approximate the integral for the marginal density of the summary statistics at q i as π Q (q i ) ≈ 1 N2 N2 j=1 π Q|X (q i |x j ). Let us note that estimates for the marginal density are often numerically zero from arithmetic underflow when evaluating the likelihood function in high dimensions.
To address this, Beck et. al [101] proposed a Laplace approximation-based importance sampling scheme to estimate the marginal likelihood. Although NMC is asymptotically unbiased for estimating EIG as (N 1 , N 2 ) → ∞, NMC is often considered too computationally costly. Even with an optimal allocation of samples to the inner and outer estimators, the mean squared error of a NMC estimator for EIG converges at a rate of O(N −2/3 ), where N is the total number of samples. However, in order to achieve this rate, one needs to have N 2 = c √ N 1 , and finding what c is involves some tuning [102].
To alleviate this cost, Wu et. al [103] considered EIG estimation based on Laplace approximations that exploit likelihood-informed subspaces and data dimension reduction. The authors numerically demonstrate that various biased estimators are often sufficient to perform BOED with unimodal posteriors that arise from affine or mildly nonlinear forward operators, unimodal priors, and additive Gaussian observation errors with small variance. The downside, however, is that accurate approximations or bounds on EIG are no longer guaranteed, and (gradient) information of the forward operator is required.
In a likelihood-free setting, estimating EIG requires us to characterize two densities, the prior and posterior, or the likelihood and marginal density of the data, which can be computationally expensive. One of the most well-known methods proposed in [104] is to estimate mutual information (MI) using knearest neighbor entropy estimators. Another method is to estimate MI using kernel density estimators (KDE) [105]. While the former requires the number of samples to scale exponentially with the MI, the latter KDE often performs poorly in high dimensions. To circumvent these challenges, Suzuki et. al [106] proposed to approximate EIG by directly estimating a ratios of densities. Thomas et. al [82] proposed an approach using logistic regression to classify prior and posterior samples in order to estimate the ratio of the prior and posterior densities. Although this method allows us to estimate the density ratio directly, it requires solving an optimization problem for each observation. Alternatively, several approaches have been proposed in the machine learning literature for constructing variational bounds for MI as proposed in [44,107]. In [44], the approximate posterior is obtained by optimizing over a parameterized family of distributions, e.g., the mean-field approximation. The choice can introduce bias in the EIG estimation.

Likelihood-free estimators for expectation information gain
We now describe an approach to estimate the EIG that uses the variational approximation π X|Q in (20) for the posterior density π X|Q . A simple calculation [44, Appx. A] shows that a lower bound for the EIG in (10a) is: Moreover, the gap between the EIG and the lower bound is equal to the expectation of the KL divergence from the variational approximation to the true posterior, i.e., E ν Q [D KL (π X|Q (·|q)|| π X|Q (·|q))]. Thus, the lower bound for the EIG provided by the variational approximation approaches the true EIG as π X|Q (·|q) better approximates the true posterior density. Given N i.i.d. samples from the joint density of parameters and summary statistics, an unbiased estimate for the lower bound in (22) is where (23) is the entropy of the prior density for the parameters. When the prior density is specified in a parametric family, this quantity is often available in closed form. As an example, for a standard Gaussian prior density π X = N (0, Id x ) where X ∈ R x , the entropy is given by E ν X [− log π X ] = 1 2 log(2πex). We note that in the general likelihood-free setting when the prior density π X is also analytically unavailable, we can also find a variational approximation π X for π X and define an estimator for the EIG as where However, as shown in [44], the estimator in (24) does not yield a bound on EIG. This is problematic if the goal is to maximize the EIG with respect to different design variables in the context of optimal experimental design.
An important property of MI, and hence EIG, is its invariance under marginal transformations of the random variables [108]. If we have two bijective functions f and g, then I(X; Q) = I(f (X); g(Q)). The property allows us to compute the EIG by finding variational approximations on transformed random variables so that the resulting EIG estimator has better statistical properties. In our numerical experiments, we transform the model parameters X defined on a bounded subset of R 3 + to the unbounded domain R 3 using the map described in Section 7.5. The transformed random variables have a standard Gaussian prior distribution allowing us to use the EIG estimator in (23). Moreover, the distributions for the transformed parameters satisfy the absolute continuity assumptions required by our approach to compute the transport maps in Section 5. A similar transformation can be defined for the summary statistics. For instance, the function g(Q) = log(Q) transforms variables fully supported on the positive numbers R + to be fully supported on R.

Case study: Model definitions
In the next two sections, we perform Bayesian calibration of the OK model for the unguided selfassembly of Di-BCP thin films with top-down microscopy image data. In this section, we focus on defining the data model and the prior distribution for the model parameters used by the LFI and EIG computation methodologies described above. In the next section, we present numerical results of model calibration and EIG computation.
We first introduce the OK model for unguided Di-BCP thin film self-assembly and the model parameters in Section 7.1. The associated forward operator and auxiliary variable are defined in Section 7.2. We then define a model for top-down microscopy characterization that contains restriction operators that mimic microscope magnification, a blurring convolution, and additive sensor noise in Section 7.3. Several lowdimensional summary statistics that quantitatively describe relevant features in the image data are proposed in Section 4.4, which are categorized into energy-based statistics and Fourier-based statistics. Lastly, we construct a model-specific prior for the parameters based on assumptions of experimental conditions in Section 7.5.

The Ohta-Kawasaki model
Consider a Di-BCP melt, composed of two types of monomers that are labeled A and B respectively, in a physical domain D ⊂ R 3 . Let u A and u B : D → [0, 1], denote the normalized segment densities for the two monomers. These segment densities satisfy the incompressibility constraint u A + u B = 1. Then, our state variable u := u A − u B ∈ [−1, 1], known as the order parameter, represents the difference between the normalized local fractions of the two phases. Under the thin film approximation, we may only consider the variation of the order parameter parallel to the substrate; thus D ⊂ R 2 is assumed in our case study 15 .
The free energy of a Di-BCP melt, E OK , given by the OK model [21,34] is assumed to be the sum of three contributing energies: the double well energy D OK , the interfacial energy I OK , and the non-local energy N OK , scaled by the model parameters (κ, , σ) ∈ R 3 + , respectively. That is, The double well energy D OK quantifies the free energy associated with the mixing of monomers and favors pure monomer phases represented by values ±1. The interfacial energy I OK is a H 1 (D) seminorm that characterizes the short-range interaction between monomers and controls the interfacial length between different monomer phases. Lastly, the nonlocal energy N OK characterizes the long-range interaction between monomers and controls the expansion of phase boundaries. The competition between the three energies during free-energy minimization forms a length scale selection mechanism that leads to various structures of equilibrium Di-BCP melts depicted by the experimental phase diagram in Fig. 1 (left). It has been shown numerically that most periodic structures observed in experiments are minimizers of the OK free energy [110,111]. The form of these contributing energies is given by: where m ∈ [−1, 1] is the spatial mean M of the order parameter, i.e., and W(u) = (1 − u 2 ) 2 /4 is the double well potential that is plotted in Fig. 8. The action of the inverse Laplacian is defined by the solution of a Poisson equation with homogeneous Neumann boundary conditions 16 : where n is the unit vector normal to ∂D. Since −1 and 1 refer to normalized pure densities of u A and u B , respectively, when the order parameter is outside of (−1, 1), the model does not represent a physically meaningful quantity. In other words, this double well potential W is non-degenerate. A degenerate, or obstacle, double well potential, such as described in [113], produces an infinite wall at these boundaries. The parameters in the OK model are connected to the commonly used material parameters of diblock copolymers. In particular, an appeal to the higher-fidelity SCFT model in [35] reveals the relationships: where N is the degree of polymerization of the polymer, χ is the Flory-Huggins parameter describing the strength of repulsion between the two monomers, l is the statistical segment length of the polymer 17 , and f ∈ (0, 1) represents the volume fraction of monomer A. In this study, we take the common assumption of κ = 1 and calibrate the collection of model parameters x = (m, , σ). We refer to the three parameters as the spatial average parameter, the interfacial parameter, and the nonlocal parameter, respectively.

The forward operator of the Ohta-Kawasaki model
The forward operator F is defined as the map from parameters x and a random initial state z to a local minimizer of the OK energy, E OK x (u) := E OK (u; 1, m, , σ), for the given parameters. That is, where the state space is defined as As a result of the non-convexity of the double well potential W, the OK energy is non-convex and does not have a unique minimizer. Thus, the local minimizers in (29) consist of an ensemble of numerous achievable equilibrium patterns for Di-BCP thin film self-assembly. We note that the states satisfying the first order condition for the minimization problem in (29) can be found by solving the following mixed problem [114].
where (·, ·) L 2 denotes L 2 (D) or (L 2 (D)) 2 inner product. The evaluation of the forward operator F(x, z) is accomplished by minimizing the energy E OK x via a variational procedure, starting from an initial state z. Conventionally, an H −1 gradient flow approach is used, which involves mimicking a time-dependent process by solving the nonlocal Cahn-Hilliard equation [23,110,111], a fourth-order time-evolving partial differential equation (PDE), whose steady-state solutions satisfy (31). According to the dynamical DFT [115], the gradient flow approach may faithfully represent the kinetic pathway of the self-assembly if an appropriate mobility operator [116] and/or a stochastic noise that satisfies the fluctuation-dissipation theorem [20] are/is incorporated. Instead, the OK model is more commonly treated as an equilibrium model, i.e., the thermal fluctuations and dynamics of Di-BCPs are not modeled.
Several authors have considered a diverse set of strategies to evaluate the forward operator via more efficient optimization schemes 18 . Although there are few studies on the direct optimization of the OK energy, quasi-Newton methods, such as Anderson acceleration [117] and semi-implicit Seidel schemes [118], are widely used in context of the SCFT model, which is also formulated as a free energy minimization problem. Similarly, a modified Newton method for the direct minimization of the OK energy was recently proposed [114] and shown to be typically three orders of magnitude faster than the gradient flow approach. The latter Newton method is adopted as the optimization procedure in our numerical case study.
Given the presence of many local minimizers for OK energy, an approach is needed to model the random initial state. A commonly adopted initial state is the realization of a uniform white noise random field [23,110,111]. A more rigorous approach that uses transformed Gaussian random fields is proposed in [114]. These random fields play the role of the auxiliary variable Z in the forward operator. When coupled with the optimization procedure, the random fields produce a distribution of equilibrium patterns. As described in Section 2.2, a physically consistent model should satisfy the equality in distribution in (3), meaning that the probability associated with each model solution matches the probability of realizing the corresponding equilibrium pattern in an experiment. The validation of these model assumptions underlying the forward operator is an important topic that we leave for future work.
We provide a description of our model for the auxiliary variable and finite element discretization of the state space in Appendix A and Appendix B, respectively. In conjunction with the variational minimization procedure for the OK model and its numerical implementation detailed in [114], this completes the definition of the forward operator in our case study. See Fig. 9 for examples of Di-BCP patterns generated by the forward operator.

A model of top-down microscopy characterization
In this subsection, we define the state-to-observable map J for top-down microscopy characterization of Di-BCP thin films. We present the restriction operator R, the imaging operator I, and the statistical sensor operator S used in our case study. Recall that J = I • R, as defined in Section 2.1.
We assume a set of experimental scenarios, in which microscope images are taken in fields of view of various sizes, mimicking the effect of varying levels of magnification. In terms of the characterization model, we first define a reference field of view D = [0, L] 2 , where L denotes its side length, corresponds to a relatively large region of the specimen top surface. Then, we define the restrictions to its centered subdomains R j (D), where j = 2, 4, . . . , 10, 20, . . . 90 denotes the percentage of the area of D retained by the restriction. The restriction operator is thus defined by R j : H 1 (D) → H 1 (R j (D)). The effect of these restriction operators, i.e., varying levels of magnification, on model calibration will be quantitatively studied in the next section. The numerical implementation of the restriction operators is discussed Appendix B and a numerical example of the synthetic data generated by these restriction operators is provided in Fig. 10.
reference field of view 70% 40% 10% For the microscopy imaging operator I, we only consider the blurring effect, which is often modeled as a point-spread function (PSF) in the context of electron or optical microscopy. The PSF is a linear integral operator which degrades the fine-scale features in the observed specimen to resemble lens-aberrations. Numerous PSFs for microscopy imaging have been considered in the literature [119][120][121], such as Airy disks, Debye diffraction integral, and Gaussian convolution integral. For our numerical case study, we employ the convolution with a Gaussian kernel 19 as a representative prototype: where I : H 1 (R(D)) → H 1 (R(D)) and σ blur is a blurring strength parameter. 19 We apply zero-padding to extend u beyond the boundary of the reference field of view ∂D.
Since the sensor noise in microscope images are generally 'fine-grained', we assume there is i.i.d. additive noise for each pixel. Other authors have considered Poisson, Gaussian or mixed Poisson-Gaussian noise random variable in a combination of additive and multiplicative form [54,120,122]; we employ additive Gaussian noise as a prototype for simplicity, although more complex and accurate models could have been employed within the introduced inference framework without additional difficulty.
In our study, we assume a 256 × 256 grid array of sensors. The sensor signal is mathematically defined as local averages in grid cells, denoted as D jk , centered at the sensor positions: 1, 2, . . . , 256). (33) We thus arrive at the following characterization model: where J = I • R as defined above.
In the next section, we examine, jointly and separately, the effect of various levels of blurring, σ blur = (None, 5L × 10 −3 , L × 10 −2 ), and/or the effect of various sensor noise variances σ 2 N , based on a signal-tonoise (SNR) ratio of (∞, 5.0, 3.0, 2.5), on the informativeness of the data. SNR is related to the sensor noise variance σ 2 N by, naively, σ −1 N = SNR.

Summary statistics for top-down microscope images of diblock copolymer thin films
In this section, several summary statistics for top-down microscopy characterization images of Di-BCP thin films, i.e. Q = Q(Y ) ∈ R q , are proposed. First, in Section 7.4.1 we introduce summary statistics defined by the composition of energy functionals and an H 1 -interpolation operator. Then, in Section 7.4.2, we present summary statistics based on the magnitude spectrum of these images obtained by discrete Fourier transform. For both cases, we discuss connections between these summary statistics and features in the image data, and provide reasoning for why we expect these summary statistics to be informative of the model parameters.

Energy-based summary statistics
We introduce five summary statistics, Q i = Q i (Y ), i = 0, 1, . . . , 4, where Q i is defined by the composition of an energy functional and an H 1 interpolation operator. The H 1 interpolation operator, denoted as Π, interpolate the image y ∈ R y into the function space H 1 (D u ) using a linear finite element basis, where D u is a unit square domain 20 . The details of the interpolation operator for our numerical case study are given in Appendix B. The energy functionals act on the H 1 -interpolated image data. Their names, notations, and mathematical forms are listed in Tab. 1.
Three summary statistics, Q 1 , Q 2 , Q 3 , are inspired by the OK free energy. They are defined using the functionals in the OK energy, i.e., D OK , I OK , and N OK , evaluated at the interpolated image data random variable Π(Y ). The motivation behind choosing these summary statistics are our perceptions that (i) the magnitude of the three energies are sensitive to the model parameters and relatively insensitive to the auxiliary variable, and (ii) they characterize informative geometric features of Di-BCP patterns, such the total area of phase interface, degree of phase separation, etc. Since the microscopy characterization visually preserves the geometry of phase interfaces, we suspect that the negative impact of sensor blurring and noise on the informativeness of these summary statistics for inferring the model parameters is relatively minor. Therefore, we hypothesize that the parameters are well-informed by these energy-based summary statistics.
We also define two additional energy-based summary statistics, given by the empirical spatial average, Q 0 = Q 0 (Y ) := (M • Π)(Y ), and a total-variation-like functional, Q 4 = Q 4 (Y ) := (T V • Π)(Y ). We expect the empirical spatial average to be informative of the material parameter f that represents the segment Notation Name number ratio of monomer A. For the OK model, f is directly related to the model parameter m through (28). The total-variation functional T V formally captures the average ratio between the interfacial length and the inclusion domain size, such as the radius of spots for spot patterns, of Di-BCP patterns. It is also relatively insensitive to the sensor noise with edge preserving property, which makes it popular in image denoising problems [123]. We visualize integrands of the functionals defining the energy-based summary statistics for synthetic data with lamella and spot latent patterns in Fig. 11.

Fourier-based summary statistics
Di-BCP thin film equilibrium patterns are well characterized by their intrinsic length scales. For example, the interfacial thickness and the periodicity length of Di-BCP patterns are often extracted from microscopy images for quantitative analysis [41,124]. In particular, the azimuthally-averaged magnitude spectrum is widely used to process and analyze the microscopy characterization images of BCP equilibrium patterns [125,126]. In the absence of external guidance or boundary geometry effects, there is no orientational preference in the latent patterns, which suggests that the orientational averaging of the magnitude spectrum preserves most of the information for inferring the model parameters. We thus hypothesize that the azimuthallyaveraged magnitude spectrum of microscopy image data and its summary statistics may be effective for model parameter inference.
We note that the azimuthally-averaged magnitude spectrum is similar to scattering intensity profiles in small angle X-ray scattering. The scattering intensity profiles are modeled as noisy Fourier domain representations of the latent states [39,127]. In this sense, the use of azimuthally-averaged magnitude spectrum and its related summary statistics for parameter inference and EIG computation can be seen as a rough prototype for the same task performed with small angle X-ray scattering data 21 .
We now define the procedure for obtaining the azimuthally-averaged magnitude spectrum from an image y * ∈ R N ×M . The discrete Fourier transform of the image outputs a complex Fourier coefficient matrix y * ∈ C N ×M . We first translate the entries to center the zero frequency and examine the magnitude spectrum: For azimuthal averaging of the magnitude spectrum, a change-of-variables to radial coordinates centered at the zero frequency is performed. In the radial coordinate, the zero-centered magnitude spectrum has N r unique radial frequencies, (r 1 , . . . , r Nr ) ∈ (0, r max ] Nr . We only study the magnitude spectrum for frequencies no larger than some right-closed intervals, as higher radial frequencies are not well-resolved and typically uninformative. Finally, we obtain the azimuthally-averaged magnitude spectrum, P a : R N ×M → R Nr , by averaging over orientations θ. Specifically, the average magnitude P a (y * )(r k ) is computed for N θ (r k ) number of entries in P(y * ) with the identical radial frequency r k : We hypothesize that not much, if any, information that is relevant toward inferring the length scale dependent parameters is lost after the azimuthal averaging. The number of unique radial frequencies N r is typically large. To further reduce the dimension of the azimuthally-averaged magnitude spectrum data, we also apply averaging over each unit interval of the radial frequency to obtain r max -dimensional averaged spectrumP a (y * ): P(y * )(r k , θ i ) , j = 1, 2, . . . , r max .
As shown in Fig. 12, the zero-centered magnitude spectrum of synthetic microscope images have 'rings' of frequencies with dominant magnitude values. They represent the periodicity length of the latent Di-BCP pattern. Furthermore, the overall shape of the magnitude spectrum distribution is sensitive to the change in model parameters and relatively insensitive to the auxiliary variable. Therefore, we further extract summary statistics such as dominant radial frequency, or some quantitative descriptions of the shape of the azimuthally-averaged magnitude spectrum, such as moments of various orders. Moreover, using Q a =P a (Y ) as summary statistics of microscopy data allows us to further investigate and compare the informativeness of the averaged magnitude spectrum in different ranges of the radial frequency for learning the model parameters.

Prior distribution
We now define the prior distribution ν X , with density π X , for the parameters used in our case study. The parameters of interest, again, are x = (m, , σ), where m ∈ (−1, 1), ∈ R + , and σ ∈ R + . We label the individual parameter random variables as X = (M, E, Σ) ∼ ν X . ν X is constructed based on hypotheses on the experimental setting related to the parameters of the OK model. These assumptions put restrictions on the feasible or admissible portions in parameter space. We Synthetic image y * ln(P(y * )) P a (y * )P a (y * ) (1) Radial frequencies Radial frequencies Interval frequencies (4) Radial frequencies Interval frequencies Radial frequencies 2 0 20 40 Radial frequencies Interval frequencies

Radial frequencies
Interval frequencies 2 Figure 12: Two synthetic top-down microscopy characterization images, y * , of Di-BCP thin films, their log-scale magnitude spectrum, ln(P(y * )), their azimuthally-averaged magnitude spectrum, Pa(y * ), with rmax = 50, and their interval-averaged magnitude spectrum,Pa(y * ). The synthetic images are generated by the OK model at the same model parameter.
denote this admissible set of parameters by A R 3 + . Then, under the principal of maximum entropy [128], we consider a uniform prior distribution 22 over this admissible set, whose density is given by: Our uniform prior is defined hierarchically, that is, the parameters are defined conditionally on the others given the variable ordering (M, E, Σ). This hierarchical distribution implies a factorization of the prior density: π X (m, , σ) = π M (m)π E|M ( |m)π Σ|M,E (σ|m, ).
The marginal density π M , for the spatial average parameter M , known as the hyper-parameter prior or hyperprior, is uniform and is independent of the other two parameter random variables (E, Σ); the conditional densities π E|M (·|m) and π Σ|M,E (·|m, ) are also uniform. We implicitly define the factors in (39) by describing the constraints that form the boundary ∂A of the admissible set A. First, we consider constraints on the spatial average parameter M . In our experimental setting, we assume that the latent material states of observed images exhibit phase separation with a high probability. Note that the homogeneously mixed state are trivial critical points of the OK energy. In this case, the length scale related parameters E and Σ do not play a role in the forward operator, and thus not identifiable from data. To define constraint boundaries for the hyper-prior π M based on this assumption, we consider a result from the analysis of the OK energy functional: Observation 1 (Constant state: Global minimum or metastable). The constant state u = m with no phase seperation is the global minimum of E OK x (25) whenever |m| > 1/ √ 3. Otherwise, it is a local minimum, i.e. a metasable state. [110,111,114].
It is often assumed that global minimum states of E OK x are relatively more stable configurations, i.e., more probable. Therefore, we exclude portion of parameter space in which the constant states are global minima, and define the admissible set for the hyper-prior for M as (0, 1/ √ 3), and the marginal prior density for M as where we consider only positive numbers in our study, without loss of generality, due to symmetry of the forward operator about m = 0.
To define admissibility constraints for the conditional density of E and Σ given the parameter M , we employ two assumptions on the equilibrium patterns encountered in our experimental scenario. These assumptions put restrictions on the size of the field of view of the microscope relative to different length scales of latent Di-BCP patterns. We note that these assumptions are based on a rescaled computational domain D c = [0, L] 2 , L = 40 , with the material and model parameters rescaled according to the relations in (28). See a discussion on the computational domain in Appendix B.
Assumption 1 (Resolving phase interface). We assume the phase interface of the latent Di-BCP patterns can be sufficiently resolved if given little or no image corruptions, which correspond to > min = 0.01.
Assumption 2 (Observing relatively large surfaces). We assume the periodicity of Di-BCP patterns can be sufficiently captured in the image data [110,111], if given little or no image corruptions: and we thus assume σ > σ min = 0.01.
To complete the admissibility constraints, we emulate the work of [110,111], in which the authors characterize a phase diagram for the OK model that describes the morphology of equilibrium states for given model parameters ( , σ). By analyzing the (linear) stability of global minimizers of the OK energy under periodic boundary conditions through Fourier expansions, they derive "order/disorder transition (ODT) curves" that partition the parameter space of ( , σ) into: (i) a region in which the OK energy functional possess the constant states as the global minimizers, (ii) a region in which the constant states are not the global minimizers, and (iii) a transition region where the analysis is not applicable. It is hypothesized in [111], based on numerical evidence, that the actual region in parameter space for which the constant states are not the global minimizers is actually larger, therefore expanding the ODT curve further into the transition region.
We follow these authors in producing a numerically discovered empirical boundary for the ODT curve based on the OK model. 23 The result of the numerical procedure is two parameterized curved boundary surfaces σ max ( ; m) and σ min ( ; m) which define the upper and lower boundary for the parameter σ conditioned on the values for m and . A visualization of the curved line "slices" of the surfaces, for three values of m are provided in Fig. 13. Intersecting these surfaces with the fixed bounds min = 0.01, σ min = 0.01 from Assumptions 1, 2, and M ∈ (0, 1/ √ 3) from Observation 1, produces the boundary ∂A.

Case study: Numerical results
In this section, we present the numerical results for our case study. These results demonstrate the insights one can extract from performing calibration tasks via the measure transport approach. In Section 8.1, we provide representative posterior samples to visually demonstrate the reduction of the model parameter uncertainty after performing the model calibration with the energy-based summary statistics. Samples of the observables, i.e., the latent Di-BCP patterns, at these parameter values are also generated and plotted to qualitatively validate the inference results. In Section 8.2, we present results on estimating EIG of model parameters from various combinations of the energy-based summary statistics with varying levels of noise . The prior distribution for the parameters (E, Σ) conditioned on the hyper-parameter M , ν E,Σ|M (·, ·|m), is uniformly distributed over the admissible set. We also plot 10,000 independently drawn samples from each conditional distribution. and blurring corruptions. The similar EIG estimation results from the Fourier-based summary statistics are provided in Section 8.3. In Section 8.4, we study the impact of various levels of magnification for microscopy imaging on the calibration results. Finally, in Section 8.5, we study the EIG estimation for the interfacial and nonlocal parameters (E, Σ) when the spatial average parameter M = m is known.
To perform our numerical experiments, we first collect n = 10 5 samples from the joint distribution of parameters and the model-generated data {(x i , y i )} n i=1 , as described in Section 4.2. Then, the samples of model-generated data are used to compute samples of the proposed summary statistics {(x i , q i = Q(y i ))} n i=1 . We use 60% of the samples as a training set for estimating the posterior probability densities using the triangular transport map approach introduced in Section 5 and use the remaining samples to estimate the EIG using the likelihood-free estimator introduced in Section 6.
A very small subset (∼ 200) of the samples are numerically problematic due to (i) failure to resolve phase interface geometries in the observable and (ii) failure to evaluate the forward operator within a reasonable residual tolerance. In both cases, the numerical anomalies primarily occur near the boundary of the admissible set ∂A and are caused by relatively coarse spatial discretization of the state space. In our numerical case study, we choose to discard these samples. To detect these numerically problematic samples, an empirical approach is adopted based on the value of the KL divergence betweenP a (y i ) and a reference set of values. We observe that the numerically problematic samples invariably possess a large value of such a KL-divergence. This ansatz is used an automated criterion to discard a majority of these numerically problematic samples. Additionally, to deal with outlier samples, we compute the 0.05 and 0.995 quantiles of the log posterior evaluations π X|Q (x i |q i ) and discard samples that are outside of these quantiles for the EIG estimation.

Posterior and posterior predictives
To visualize posterior predictive samples, we first generate representative posterior quadrature points. In particular, we take Gauss-Hermite quadrature points r GH which are used for accurately integrating functions with respect to the standard normal distribution, i.e, the reference distribution ν R2 , and transform them through the inverse of the posterior transport map, i.e. S X (q * , ·) −1 (r GH ), as described in (20). These points are depicted within the admissible set for the prior in Tab. 2. While these points, paired with appropriate weights, can be used for both visualizing and accurately integrating functions with respect to the posterior distribution, we only present numerical results on the visualization of the model predictions at posterior samples in this study. We note that the points include the inverse transformation of the center quadrature point, x 0 = S X (q * , ·) −1 (0), as a representative sample that is typically in the region with highest posterior probability 24 . In addition, we include the transformation of quadrature points in the tails of each reference marginal N (0, 1), in order to study the 'tails', or regions of low probability, of the posterior. These points are labeled as x ±δm := S X (q * , ·) −1 (0 ± δm), x ±δ := S X (q * , ·) −1 (0 ± δ ), and x ±δσ := S X (q * , ·) −1 (0 ± δσ) where δm = ∆e 1 , δ = ∆e 2 and δσ = ∆e 3 for some scalar ∆ > 0. In our studies we set ∆ = 2, which corresponds to the 0.977 quantile of a univariate standard normal random variable.
In Tab Next, we present posterior predictive samples in Tab. 3 in order to examine the downstream effects of Bayesian calibration on model predictions. That is, we generate samples from the posterior predictive distribution The random variable U pred = F(X|q * , Z) ∼ ν U pred |Q (·|q * ) represents the Di-BCP patterns predicted by the forward operator F(·, Z) given samples of the calibrated parameters. These posterior predictive samples help us to visually understand the gain in prediction confidence from model calibration. Generating these samples is accomplished by simulating from the forward operator x → F(x, z i ) where z i ∼ ν Z|M (·|m) and x are the quadrature points in Tab. 2. The Di-BCP pattern predictions of models calibrated according to the summary statistics Q = (Q 2 , Q 4 ) (top row for each u * ) are much more visually dissimilar than the predictions of models calibrated according to the summary statistics Q = Q 1:4 (bottom row for each u * ), regardless of the phases of latent Di-BCP patterns.
Overall, we see that the energy-based summary statistics Q 1:4 are highly informative of the OK model parameters, despite having a much lower dimension than the image data. We note that sampling from all three posterior distributions in Tab. 2 and 3 is easily accomplished after finding the map S X . Furthermore, evaluating the posterior densities can also be easily accomplished using the map, as discussed in Section 5.

Expected information gain: Energy-based summary statistics
In this subsection, we present results on the estimation of EIGs for inferring the model parameters from the energy-based summary statistics introduced in Section 4.4.
In Fig. 14, we plot the estimation of EIGs for inferring each marginal parameter random variable in X = (M, E, Σ) from the energy-based summary statistics when observing raw images in the reference field of view, i.e., assuming no noise or blurring corruption in the characterization process. The EIG calculations imply that the empirical spatial average summary statistics M strongly informs the spatial average parameter M , as expected (left). Furthermore, other energy-based summary statistics also jointly inform the spatial Table 2: (left column) The latent material states u * , described by order parameters, of three synthetic microscopy image y * data of Di-BCP thin films. (middle column) The Gauss-Hermite quadrature points of the reference density pulled back to the model parameter posterior, illustrating posterior concentration conditioned on q * = (Q 2 (y * ), Q 4 (y * )). (right column) The same as the middle column, but conditioned on q * = (Q 1 (y * ), . . . , Q 4 (y * )).
average parameter. This corroborates what we learn from the phase diagram in Fig. 1: the spatial average parameter, or the volume ratio parameter f , can be learned from the morphological phase of the Di-BCP material. The interfacial energy Q 2 and/or the total variation energy Q 4 strongly inform the interfacial parameter E (center ), and, finally, the double well energy Q 1 and nonlocal energy Q 3 jointly strongly inform the nonlocal parameter Σ (right). The EIG calculations also reveal that even though, marginally, Table 3: Posterior predictive samples, F (x, z i ), where each image on the right column is generated by independent draws from z i ∼ ν Z|M (·|m), The model parameter samples x 0 , x ±δm , x ±δ , and x ±δσ are transformed quadrature points shown in Table 2.
the double well energy Q 1 does not strongly inform Σ, Q 1 and Q 3 are jointly very informative of Σ. In fact, the estimated EIG from these two summary statistics, I(Σ; (Q 1 , Q 3 )), is larger than the sum of the estimated EIGs from each summary statistics marginally, I(Σ; Q 1 ) + I(Σ; Q 3 ) 25 . Similarly, we observe that Q 1 and Q 3 jointly inform the interfacial parameter E 'almost' as much as Q 2 or Q 4 each informs E, even though Q 1 and Q 3 marginally each provide very little information for learning E. Next, to investigate the effects of additive noise and blurring corruptions on the informativeness of the energy-based summary statistics, we plot in Fig. 15 the estimated EIGs for inferring each model parameter from its most informative summary statistics according to Fig. 14, with different levels of additive noise and blurring corruptions specified in Section 7.3. First, we observe that while the empirical spatial average summary statistics M is highly informative of the spatial average parameter in the absence of the noise and blurring corruptions, its informativeness deteriorates by half with the presence of additive noise (left), despite the noise being zero mean. Thus, we may have difficulty directly pinpointing the spatial average parameter from the empirical spatial average statistics M of noisy characterization images. On the other hand, the informativeness of (Q 1 , Q 3 ) for inferring Σ is not significantly degraded by the various degrees of noise and blurring corruptions (right), while the informativeness of (Q 2 , Q 4 ) for inferring E is more significantly degraded by the blurring corruption than by the additive noise corruption (center ). The variables of these grids conceptually represent the specifications of the characterization instrument one might be able to choose or control. Ultimately, the EIG computation can help one determine the design of characterization experiments, trading-off the cost of experiments and model calibration accuracy. While some of these quantitative results can sometimes be theoretically deduced without numerical simulations, often these relationships are complex and intractable. Thus it is highly useful to be able to empirically extract these relationships with numerical simulations instead. Using the framework for EIG computation introduced in Section 6, each EIG value in the entries of the 3 × 4 grids is estimated by a map estimate S X , given the parameter and summary statistics joint samples (x i , q i ) ∼ ν X,Q .
Lastly, we note that the EIGs for learning model parameters are examined marginally, although the same computations could have been conducted to estimate the EIG for the three parameters jointly. However, since we were mostly interested in the impact of various summary statistics on learning each individual parameter, we leave EIG studies involving the joint parameter X = (M, E, Σ) for future study.

Expected information gain: Fourier-based summary statistics
Next, we computed the EIG from the Fourier-based summary statistics for learning each of the three parameters. In particular, we consider Q a =P a (Y ) ∈ R 40 , as defined in (37), where each entry corresponds to the averaged magnitude spectrum value in a unit interval of the radial frequency. Let Q a,i denote the ith component of Q a . In Fig. 16 and 17, we show the estimated EIGs for inferring marginally each of the three parameters from each component of Q a , with various degrees of additive noise and blurring corruptions. With the presence of the additive noise corruption, the informativeness of the averaged magnitude spectrum decreases, while the blurring corruption primarily reduces the information contained in the averaged magnitude spectrum at higher radial frequencies. Both the blurring and noise corruptions tend to more strongly impact the EIG of the spectrum at frequencies that are most informative for each parameter, with the exception of the EIG for Σ subjects to blurring corruptions (Fig. 16, right). Moreover, the spectrum at higher frequencies generally tend to inform M and E more while the spectrum at lower frequencies (∼ 10) inform Σ more, so that the spectrum at both higher and lower frequencies play a non-negligible role in informing the model parameters. For the remainder of this subsection, we address the difference in EIG approximation quality when using the parametric transport framework laid out in Section 6, compared to a traditional Gaussian approximation based EIG approximation. In Fig. 18, the EIG computed with our non-Gaussian parametric approach is compared to the EIG computed via a Gaussian approximation of the posterior densities. As expected from expression (22), the EIG estimator with a variational approximation to the posterior density is negatively biased and this bias increases prediction for an approximate density that has a larger discrepancy with respect to the true posterior (see Section 6.2). For some components of the averaged magnitude spectrum, this results in an estimated EIG based on the Gaussian approximation (red ) that is near zero, while the actual information gain in those components (blue) is quite larger as seen from the non-Gaussian approximation. This kind of bias can hamper experimental design decision-making that relies on EIG estimations. We further demonstrate how poor the approximation can be by visualizing the Gaussian approximation of the joint density that underrepresents the mutual information between two random variables. Specifically, we consider the joint density of Σ and Q a,17 . In Fig. 19, we plot the joint samples from this distribution on the left along with a Gaussian and non-Gaussian approximation to the joint density. While we considered the EIG from the individual components of the averaged magnitude spectrum for inferring each parameter, interesting future work would be to determine the subset of spectrum data that jointly inform the parameters most. Techniques for doing so are found within the field of active learning or optimal experimental design-these fields can benefit from the framework laid out in the paper, since they often rely on EIG estimates in order to evaluate optimization objectives.

Expected information gain: Levels of magnification
In this subsection, we consider the effects of varying levels of magnification, specified by the restriction operator R defined in Section 7.3, on the informativeness of energy-based summary statistics, M , (Q 1 , Q 3 ) or (Q 2 , Q 4 ), about X = M , E, or Σ, respectively. In particular, we plot the estimated EIG as a function of magnification levels for a fixed level of blurring corruption and varying levels of noise corruptions in Fig. 20. Similarly, we plot the same quantities in Fig. 21 for a fixed level of noise corruption and varying levels of blurring corruptions.  For learning the spatial average parameter M (left in Fig. 20 and 21), the estimated EIG from the empirical mean summary statistics M decreases as the magnification level increases. Moreover, the estimated EIG is insensitive to changes in noise and blurring corruptions, caused by the spatial average-preserving property of the Gaussian blur and the zero-mean property of the noise. We expect this behaviour to change for more complex characterization models that does not satisfy these properties. These numerical results suggest that, under the assumptions of the data model and parameter prior distribution, the action of zooming in will not lead to a decrease in uncertainty in M on average. These result may be interpreted by the loss of phase information related to the material parameter f (see Fig. 1) when we only observe a small portion of the top surface of a Di-BCP thin film. For example, one might only observe several spots in a pattern that contains a mixture of spots and strips.
For learning the interfacial parameter E (middle in Fig. 20 and 21), the estimated EIG from the energybased summary statistics (Q 2 , Q 4 ) decreases as the level of magnification increases in the presence of noise and blurring corruptions. This also suggests that the action of zooming in will not lead to a decrease in uncertainty in the E parameter on average. We see that in an idealistic noise-free setting, however, the estimated EIG increases slightly when the level of magnification increases. We hypothesize that the increase 38 in image pixel density via magnification leads to better resolved phase interfaces only in the presence of the blurring corruption, and thus leads to a small increase in the estimated EIG for learning E.
For learning the nonlocal parameter Σ (right in Fig. 20 and 21), we observe an interesting response of the estimated EIG from the energy-based summary statistics (Q 1 , Q 3 ) with respect to changes in the levels of magnification and image corruptions. In the high SNR setting, the estimated EIG decreases as the magnification level increases. On the contrary, we see an increase in the estimated EIG in the low SNR setting. Similarly, in the low blurring setting, the estimated EIGs decreases as the magnification level increases. We also observe an increase in the estimated EIG in the high blurring setting. These results suggest that whether the action of zooming in is beneficial for inferring Σ depends on the level of image corruptions.
Overall, while some of the figures behave as one might expect, the unexpected behavior of the EIG for certain summary statistics and model parameters would otherwise be difficult to determine quantitatively without these numerical experiments. These empirical findings could aid and inspire further theoretical investigations about the effects of various image corruptions and experimental designs on model calibration and prediction.

Expected information gain: Known volume ratio of monomer blocks
In many experimental settings, the volume ratio of the two monomer blocks in Di-BCPs can be determined with high certainty. In particular, we have shown in Section 8.2 that the spatial average summary statistics is exceptional in learning the OK parameter m. Therefore, in this subsection, we consider the inference problem for two length scale related parameters in the OK model (E, Σ) given image data of Di-BCP patterns with a fixed, or known, spatial average parameter. Figures 22,23, and 24 plot the estimated EIGs for learning the parameters (E, Σ) jointly from (Q 1 , Q 3 ), (Q 2 , Q 4 ), and (Q 1 , . . . , Q 4 ), respectively, conditioned on three different spatial averages m = (0, 0.3, 1/ √ 3). In general, we observe a decrease in the estimated EIGs with with increasing level of noise and blurring corruptions. We also note that the EIGs from (Q 2 , Q 4 ) are more sensitive to noise than those from (Q 1 , Q 3 ) or conditioning on all observables, particularly for the spatial averages m = 0 and m = 0.3.

Conclusion
In this work, we explored the Bayesian model calibration problem for block copolymer self-assembly using imaging data. To set up the problem, we introduced a generic probabilistic data model that includes both the measurement uncertainty in the characterization process and the aleatoric uncertainty in the self-assembly procedure. To calibrate the model parameters we characterize their posterior distribution, where the likelihood function accounts for the aleatoric uncertainty by marginalizing over high-dimensional auxiliary variables. Moreover, we quantify the expected information gain to measure the informativeness of data for calibrating the parameters.  We argue and demonstrate that, as a result of the generally intractable integrated likelihood evaluations and the high dimensional data, many conventional Bayesian inference methodologies are unsuitable for performing the proposed calibration task. To overcome these challenges, we propose to (i) employ a likelihood-free inference method that only requires joint samples of data and model parameters, and (ii) construct informative low-dimensional summary statistics for the data. In particular, we use an inference methodology based on triangular transport maps, which allows us to (i) consistently characterize the joint distribution of model parameters and summary statistics with limited samples, (ii) perform amortized inference for any realization of the data with a single collection of joint samples, and (iii) consistently compute expected information gain with no significant additional cost, which enables us to study the effects of different summary statistics and corruptions in the data on the calibration results.
Lastly, we demonstrated the power of the proposed inference approach and the computation of expected information gain on a numerical case study of calibrating the Ohta-Kawasaki model parameters using top-down microscopy characterization images. We introduced several energy-and Fourier-based summary statistics to reduce the dimension of the image data. A hierarchical uniform prior distribution for the model parameters is constructed by analyzing the forward operator. First, we observe the concentration of the posterior distribution in the parameter space and visualize posterior predictive samples, which show that the proposed energy-based summary statistics are effective for calibrating the parameters with low uncertainty. We then estimate the expected information gains from different summary statistics for learning the model parameters, under various experimental scenarios with different levels of (i) noise and blurring corruptions, and (ii) magnifications of the microscope. These numerical results allow us to investigate the effects of experimental settings on the calibration results for this representative model, and enable further optimal experimental design studies of interest.
where K 1 is the order 1 modified Bessel function of the second kind and σ 2 = Γ(1) Γ(2)4πκ [131]. In this work, κ = 200 is chosen such that the spatial correlation is below the length scales of the features of all equilibrium states in our experimental scenario modeled by the OK model (25). A is isotropic and stationary since the function c depends on Euclidean distance |s 1 − s 2 | alone. For this choice of covariance, we have almost surely continuous realizations [57,Sec. 6.5]. Realizations A i ∼ ν A can be numerically approximated given a sufficiently fine spatial discretization. After discretization, these fields become Gauss-Markov random fields [132], which are efficiently implemented using the hIPPylib python package [133], where finite element approximations of these PDEs are solved in the FEniCS library [134,135].. A description of finite element discretizations for similar random fields can be found in [47]. Finally, we define the transformation Z = Z(x, A) such that the random field Z is a pointwise uniform a.s. continuous field that satisfies certain physical constraints. Specifically, given the point-wise variance σ 2 , the pointwise transformation ensures that the pointwise values of the auxiliary random field realizations z i = Z(X, A i ) for A i ∼ ν A are within the physical order parameter bounds (−1, 1) and have the pointwise mean m. Furthermore, its pointwise marginal densities are π Z(s)|M (·|m) = U(−1+2m, 1) ∀s ∈ R 2 . Thus, they are physically reasonable random fields to employ as starting points Z for the variational procedure that defines F(x, Z).
To exactly match the specified spatial average m for the order parameter (27), an additional transformation is employed, in which we choose the closest point satisfying the correct spatial average within the chosen finite element function space. Since this projection does not account for the pointwise constraint z i (s) ∈ (−1, 1), transformed random field samples may actually violate the constraint in regions of the domain-therefore, our samples only approximately satisfy this physical constraint. However, this has a negligible effect on forward solutions, since the double well potential D OK heavily penalizes deviation from non-physical regions in state space, so that any variational procedure will correct for this within a few iterations.
Appendix B. The discretization of the state space, the restriction operators, and image interpolation The state space V u for th OK model isH 1 m D , as in (30), with D corresponds to the reference field of view, as presented in Section 7.3, and m corresponds to the spatial average model parameter. However, when evaluating the forward operator, we solve the equivalent mixed problem as in (31), thus the state space is transformed into the product space (H 1 (D)) 2 . Furthermore, we transform the physical domain D into a computational domain D c , with a side length of 40, thus equivalently rescale the model parameters according to (28). The purpose of the domain scaling is to improve the numerical conditioning of the forward operator by balancing the relative magnitude of the three terms in the OK free energy (25).
The computational domain D c is discretized with uniform triangular mesh, denoted by P h , such that D c = K∈P hK . Each element is the image of a scaling map F K of the reference triangular elementK. We consider a linear finite element basis: let P 1 (K) denote the space of complete linear polynomials, then the finite-dimensional spaces S h,1 := S h,1 (P h ) by: In our case study, we consider P h with same number of cells as the number of sensors, which makes the evaluation of the sensor signal discretization operator D (33) relatively straight forward. The restriction operators R j , as introduced in Section 7.3, are applied by projecting the space S h,1 to spaces defined with denser uniform triangular meshes P hj . The denser meshes have the same cells arrangement as P h in the subdomain domains R j (D c ). We then apply the imaging operator I and the statistical sensor operator S to the order parameter restricted to R j (D c ).
The interpolation of images into H 1 (D u ) is implemented through the similar finite element basis as define in (B.1). The DoFs of the finite element basis match the the dimension of the image and the spatial arrangement of the pixels is preserved.
The finite element discretization, the restriction operators, and image interpolation are implemented through the FEniCS library.