From dbb901f7b4a7c129e7fc143e5fc68d64336d8760 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Tue, 17 May 2016 18:32:44 -0600 Subject: [PATCH 01/20] Informative names for simpleFunP --- bet/calculateP/simpleFunP.py | 82 ++++---- .../linearMap/linearMapUniformSampling.py | 196 +++++++----------- 2 files changed, 116 insertions(+), 162 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index f3292214..7b2fd431 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -20,7 +20,7 @@ class wrong_argument_type(Exception): """ -def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6): +def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rect_scale=0.2, M=50, num_d_emulate=1E6): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` where :math:`\rho_{\mathcal{D}}` is a uniform probability density on @@ -45,9 +45,9 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6): :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - play around with it and you can get reasonable results with a relatively small number here like 50. - :param bin_ratio: The ratio used to determine the width of the - uniform distributiion as ``bin_size = (data_max-data_min)*bin_ratio`` - :type bin_ratio: double or list() + :param rect_scale: The scale used to determine the support of the + uniform distribution as ``rect_size = (data_max-data_min)*rect_scale`` + :type rect_scale: double or list() :param int num_d_emulate: Number of samples used to emulate using an MC assumption :param data_set: Sample set that the probability measure is defined for. @@ -71,9 +71,10 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6): dim = data_set.shape[1] values = data_set else: - raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray") + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") - bin_size = (np.max(values, 0) - np.min(values, 0))*bin_ratio + rect_size = (np.max(values, 0) - np.min(values, 0))*rect_scale r''' @@ -100,7 +101,7 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6): ''' if comm.rank == 0: - d_distr_samples = 1.5*bin_size*(np.random.random((M, + d_distr_samples = 1.5*rect_size*(np.random.random((M, dim))-0.5)+Q_ref else: d_distr_samples = np.empty((M, dim)) @@ -120,7 +121,7 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6): ''' # Generate the samples from :math:`\rho_{\mathcal{D}}` num_d_emulate = int(num_d_emulate/comm.size)+1 - d_distr_emulate = bin_size*(np.random.random((num_d_emulate, + d_distr_emulate = rect_size*(np.random.random((num_d_emulate, dim))-0.5) + Q_ref # Bin these samples using nearest neighbor searches @@ -151,7 +152,7 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6): data_set._output_probability_set = s_set return s_set -def normal_normal(data_set, Q_ref, M, std, num_d_emulate=1E6): +def normal_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate=1E6): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability @@ -243,7 +244,8 @@ def normal_normal(data_set, Q_ref, M, std, num_d_emulate=1E6): if isinstance(data_set, samp.discretization): data_set._output_sample_set = s_set return s_set -def unif_normal(data_set, Q_ref, M, std, num_d_emulate=1E6): + +def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate=1E6): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability @@ -319,7 +321,7 @@ def unif_normal(data_set, Q_ref, M, std, num_d_emulate=1E6): data_set._output_probability_set = s_set return s_set -def uniform_hyperrectangle_user(data_set, domain, center_pts_per_edge=1): +def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domain, center_pts_per_edge=1): r""" Creates a simple function appoximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho{\mathcal{D}, M}` is a uniform probablity density over the @@ -332,9 +334,9 @@ def uniform_hyperrectangle_user(data_set, domain, center_pts_per_edge=1): :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` - :param domain: The domain overwhich :math:`\rho_\mathcal{D}` is + :param rect_domain: The domain overwhich :math:`\rho_\mathcal{D}` is uniform. - :type domain: :class:`numpy.ndarray` of shape (2, mdim) + :type rect_domain: :class:`numpy.ndarray` of shape (2, mdim) :param list() center_pts_per_edge: number of center points per edge and additional two points will be added to create the bounding layer @@ -356,17 +358,18 @@ def uniform_hyperrectangle_user(data_set, domain, center_pts_per_edge=1): dim = data_set.shape[1] values = data_set else: - raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray") + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") data = values - domain = util.fix_dimensions_data(domain, data.shape[1]) - domain_center = np.mean(domain, 0) - domain_lengths = np.max(domain, 0) - np.min(domain, 0) + rect_domain = util.fix_dimensions_data(rect_domain, data.shape[1]) + domain_center = np.mean(rect_domain, 0) + domain_lengths = np.max(rect_domain, 0) - np.min(rect_domain, 0) - return uniform_hyperrectangle_binsize(data_set, domain_center, domain_lengths, + return regular_partition_uniform_distribution_rectangle_size(data_set, domain_center, domain_lengths, center_pts_per_edge) -def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size, +def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_size, center_pts_per_edge=1): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` @@ -377,9 +380,9 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size, to represent it exactly with ``M = 3^mdim`` or rather ``len(d_distr_samples) == 3^mdim``. - :param bin_size: The size used to determine the width of the uniform + :param rect_size: The size used to determine the width of the uniform distribution - :type bin_size: double or list() + :type rect_size: double or list() :param int num_d_emulate: Number of samples used to emulate using an MC assumption :param data_set: Sample set that the probability measure is defined for. @@ -407,7 +410,8 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size, dim = data_set.shape[1] values = data_set else: - raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray") + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") data = values @@ -420,15 +424,15 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size, print 'Using 1 in each dimension.' if np.any(np.less(center_pts_per_edge, 0)): print 'Warning: center_pts_per_edge must be greater than 0' - if not isinstance(bin_size, collections.Iterable): - bin_size = bin_size*np.ones((dim,)) - if np.any(np.less(bin_size, 0)): + if not isinstance(rect_size, collections.Iterable): + rect_size = rect_size*np.ones((dim,)) + if np.any(np.less(rect_size, 0)): print 'Warning: center_pts_per_edge must be greater than 0' sur_domain = np.array([np.min(data, 0), np.max(data, 0)]).transpose() points, _, rect_domain = vHist.center_and_layer1_points_binsize\ - (center_pts_per_edge, Q_ref, bin_size, sur_domain) + (center_pts_per_edge, Q_ref, rect_size, sur_domain) edges = vHist.edges_regular(center_pts_per_edge, rect_domain, sur_domain) _, volumes, _ = vHist.histogramdd_volumes(edges, points) s_set = vHist.simple_fun_uniform(points, volumes, rect_domain) @@ -437,11 +441,11 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size, data_set._output_probability_set = s_set return s_set -def uniform_hyperrectangle(data_set, Q_ref, bin_ratio, center_pts_per_edge=1): +def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rect_scale, center_pts_per_edge=1): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density - centered at Q_ref with bin_ratio of the width + centered at Q_ref with rect_scale of the width of D. Since rho_D is a uniform distribution on a hyperrectanlge we should be able @@ -450,9 +454,9 @@ def uniform_hyperrectangle(data_set, Q_ref, bin_ratio, center_pts_per_edge=1): :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` - :param bin_ratio: The ratio used to determine the width of the - uniform distributiion as ``bin_size = (data_max-data_min)*bin_ratio`` - :type bin_ratio: double or list() + :param rect_scale: The scale used to determine the width of the + uniform distributiion as ``rect_size = (data_max-data_min)*rect_scale`` + :type rect_scale: double or list() :param int num_d_emulate: Number of samples used to emulate using an MC assumption :param Q_ref: :math:`Q(\lambda_{reference})` @@ -477,17 +481,18 @@ def uniform_hyperrectangle(data_set, Q_ref, bin_ratio, center_pts_per_edge=1): dim = data_set.shape[1] values = data_set else: - raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray") + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") data = values - if not isinstance(bin_ratio, collections.Iterable): - bin_ratio = bin_ratio*np.ones((dim, )) + if not isinstance(rect_scale, collections.Iterable): + rect_scale = rect_scale*np.ones((dim, )) - bin_size = (np.max(data, 0) - np.min(data, 0))*bin_ratio - return uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size, + rect_size = (np.max(data, 0) - np.min(data, 0))*rect_scale + return regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_size, center_pts_per_edge) -def uniform_data(data_set): +def uniform_partition_uniform_distribution_data_samples(data_set): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density over @@ -522,7 +527,8 @@ def uniform_data(data_set): s_set = samp.sample_set(dim = dim) s_set.set_values(values) else: - raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray") + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") s_set.set_probabilities(np.ones((num,), dtype=np.float)/num) diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index 1aca04a4..b91e5c23 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -3,9 +3,24 @@ # Copyright (C) 2014-2015 The BET Development Team """ -This example generates uniform samples on a 3D grid -and evaluates a linear map to a 2d space. Probabilities -in the paramter space are calculated using emulated points. +This example solves a stochastic inverse problem for a +linear 3-to-2 map. We refer to the map as the QoI map, +or just a QoI. We refer to the range of the QoI map as +the data space. +The 3-D input space is discretized with i.i.d. uniform +random samples. We refer to the input space as the +parameter space, and use parameter to refer to a particular +point (e.g., a particular random sample) in this space. +A reference parameter is used to define a reference QoI datum +and a uniform probability measure is defined on a small box +centered at this datum. +The measure on the data space is discretized either randomly +or deterministically, and this discretized measure is then +inverted by BET to determine a probability measure on the +parameter space whose support contains the measurable sets +of probable parameters. +We use emulation to estimate the measures of sets defined by +the random discretizations. 1D and 2D marginals are calculated, smoothed, and plotted. """ @@ -15,168 +30,101 @@ import bet.calculateP.simpleFunP as simpleFunP import bet.calculateP.calculateP as calculateP import bet.postProcess.plotP as plotP -import bet.sample as sample +import bet.sample as samp import bet.sampling.basicSampling as bsam +from myModel import my_model # Initialize 3-dimensional input parameter sample set object -input_samples = sample.sample_set(3) +input_samples = samp.sample_set(3) + # Set parameter domain input_samples.set_domain(np.repeat([[0.0, 1.0]], 3, axis=0)) -# Set reference parameter -ref_lam = [0.5, 0.5, 0.5] - +# Define the sampler that will be used to create the discretization +# object, which is the fundamental object used by BET to compute +# solutions to the stochastic inverse problem ''' Suggested changes for user: - -Try setting n0, n1, and n2 all to 10 and compare the results. - -Also, we can do uniform random sampling by setting - - random_sample = True - -If random_sample = True, consider defining - - n_samples = 1E3 - -Then also try n_samples = 1E4. What happens when n_samples = 1E2? -''' -random_sample = True - -if random_sample == False: - n0 = 30 # number of samples in lam0 direction - n1 = 30 # number of samples in lam1 direction - n2 = 30 # number of samples in lam2 direction - n_samples = n0*n1*n2 -else: - n_samples = 2E3 - -# Set the parameter samples -if random_sample == False: - sampler = bsam.sampler(None, n_samples) - - # Define a regular grid for the samples, eventually update to use basicSampling - vec0=list(np.linspace(lam_domain[0][0], lam_domain[0][1], n0)) - vec1 = list(np.linspace(lam_domain[1][0], lam_domain[1][1], n1)) - vec2 = list(np.linspace(lam_domain[2][0], lam_domain[2][1], n2)) - vecv0, vecv1, vecv2 = np.meshgrid(vec0, vec1, vec2, indexing='ij') - - input_samples.set_values( - np.vstack((vecv0.flat[:], - vecv1.flat[:], - vecv2.flat[:])).transpose() - ) -else: - # Use uniform i.i.d. random samples from the domain - sampler = bsam.sampler(None, n_samples) - input_samples = sampler.random_sample_set('random', input_samples) -# QoI map -Q_map = np.array([[0.506, 0.463],[0.253, 0.918], [0.085, 0.496]]) - -# reference QoI -Q_ref = np.array([0.422, 0.9385]) - -# calc data -data= np.dot(samples,Q_map) -np.savetxt('3to2_samples.txt.gz', samples) -np.savetxt('3to2_data.txt.gz', data) +Try num_samples = 1E3 and 1E4. What happens when num_samples = 1E2? +Try using 'lhs' instead of 'random' in the sampler. +''' +sampler = bsam.sampler(my_model, num_samples=1E4) +# Generate samples on the parameter space +input_samples = sampler.random_sample_set('random', input_samples) +# Estimate volumes of Voronoi cells associated with the parameter samples ''' Suggested changes for user: - -Try different ways of discretizing the probability measure on D defined as a uniform -probability measure on a rectangle (since D is 2-dimensional). - -unif_unif creates a uniform measure on a hyperbox with dimensions relative to the -size of the circumscribed hyperbox of the set D using the bin_ratio. A total of M samples -are drawn within a slightly larger scaled hyperbox to discretize this measure defining -M total generalized contour events in Lambda. The reason a slightly larger scaled hyperbox -is used to draw the samples to discretize D is because otherwise every generalized contour -event will have non-zero probability which obviously defeats the purpose of "localizing" -the probability within a subset of D. - -uniform_hyperrectangle uses the same measure defined in the same way as unif_unif, but the -difference is in the discretization which is on a regular grid defined by center_pts_per_edge. -If center_pts_per_edge = 1, then the contour event corresponding to the entire support of rho_D -is approximated as a single event. This is done by carefully placing a regular 3x3 grid (since D=2 in -this case) of points in D with the center point of the grid in the center of the support of -the measure and the other points placed outside of the rectangle defining the support to define -a total of 9 contour events with 8 of them having exactly zero probability. -''' -deterministic_discretize_D = False -if deterministic_discretize_D == True: - (d_distr_prob, d_distr_samples, d_Tree) = simpleFunP.uniform_hyperrectangle(data=data, - Q_ref=Q_ref, bin_ratio=0.2, center_pts_per_edge = 1) +Try different numbers of points to estimate the volume, and also +try comparing to the standard Monte Carlo assumption that all the +Voronoi cells have the same measure. +''' +MC_assumption = False +if MC_assumption is False: + input_samples.estimate_volume(n_mc_points=1E5) else: - (d_distr_prob, d_distr_samples, d_Tree) = simpleFunP.unif_unif(data=data, - Q_ref=Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E5) + input_samples.estimate_volume_mc() + +# Now create the discretization object using the input samples +my_discretization = sampler.create_random_discretization('random', input_samples, + savefile = '3to2_discretization.txt.gz') + +# Define the reference parameter +param_ref = np.array([0.5, 0.5, 0.5]) +# Compute the reference QoI +Q_ref = my_model(param_ref) ''' Suggested changes for user: - -If using a regular grid of sampling (if random_sample = False), we set - - lambda_emulate = samples - -Otherwise, play around with num_l_emulate. A value of 1E2 will probably -give poor results while results become fairly consistent with values -that are approximately 10x the number of samples. - -Note that you can always use - - lambda_emulate = samples - -and this simply will imply that a standard Monte Carlo assumption is -being used, which in a measure-theoretic context implies that each -Voronoi cell is assumed to have the same measure. This type of -approximation is more reasonable for large n_samples due to the slow -convergence rate of Monte Carlo (it converges like 1/sqrt(n_samples)). -''' -if random_sample == False: - lambda_emulate = samples -else: - lambda_emulate = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, num_l_emulate = 1E5) +Try different ways of discretizing the probability measure on D defined as a uniform +probability measure on a rectangle (since D is 2-dimensional) centered at Q_ref whose +size is determined by scaling the circumscribing box of D. +''' +simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.2, + center_pts_per_edge = 5) +#simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled( +# data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.2, +# M=50, num_d_emulate=1E5) # calculate probablities -(P, lambda_emulate, io_ptr, emulate_ptr) = calculateP.prob_emulated(samples=samples, - data=data, - rho_D_M=d_distr_prob, - d_distr_samples=d_distr_samples, - lambda_emulate=lambda_emulate, - d_Tree=d_Tree) +calculateP.prob(my_discretization) + # calculate 2d marginal probs ''' Suggested changes for user: - + At this point, the only thing that should change in the plotP.* inputs should be either the nbins values or sigma (which influences the kernel density estimation with smaller values implying a density estimate that looks more like a histogram and larger values smoothing out the values more). - + There are ways to determine "optimal" smoothing parameters (e.g., see CV, GCV, and other similar methods), but we have not incorporated these into the code as lower-dimensional marginal plots have limited value in understanding the structure of a high dimensional non-parametric probability measure. ''' -(bins, marginals2D) = plotP.calculate_2D_marginal_probs(P_samples = P, samples = lambda_emulate, lam_domain = lam_domain, nbins = [10, 10, 10]) +(bins, marginals2D) = plotP.calculate_2D_marginal_probs(input_samples, + nbins = 30) # smooth 2d marginals probs (optional) -marginals2D = plotP.smooth_marginals_2D(marginals2D,bins, sigma=0.1) +marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.1) # plot 2d marginals probs -plotP.plot_2D_marginal_probs(marginals2D, bins, lam_domain, filename = "linearMap", - plot_surface=False) +plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "linearMap", + file_extension = ".eps", plot_surface=False) # calculate 1d marginal probs -(bins, marginals1D) = plotP.calculate_1D_marginal_probs(P_samples = P, samples = lambda_emulate, lam_domain = lam_domain, nbins = [10, 10, 10]) +(bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, + nbins = 40) # smooth 1d marginal probs (optional) marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.1) # plot 2d marginal probs -plotP.plot_1D_marginal_probs(marginals1D, bins, lam_domain, filename = "linearMap") - +plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "linearMap", + file_extension = ".eps") From a2904610d5fed45938baeb6eb0ab05ea0b4ee9f3 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Thu, 19 May 2016 14:47:14 -0600 Subject: [PATCH 02/20] Added different discretizations of simpleFunP --- bet/calculateP/simpleFunP.py | 567 +++++++++++++++++++++-------------- 1 file changed, 348 insertions(+), 219 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 7b2fd431..ddc42900 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -18,25 +18,26 @@ class wrong_argument_type(Exception): Exception for when the argument for data_set is not one of the acceptible types. """ - -def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rect_scale=0.2, M=50, num_d_emulate=1E6): + +def uniform_partition_uniform_distribution_rectangle_size(data_set, Q_ref, + rect_size, M=50, + num_d_emulate=1E6): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` where :math:`\rho_{\mathcal{D}}` is a uniform probability density on a generalized rectangle centered at Q_ref. - The support of this density is defined by bin_ratio, which determines - the size of the generalized rectangle by scaling the circumscribing - generalized rectangle of :math:`\mathcal{D}`. - The simple function approximation is then defined by determining M + The support of this density is defined by rect_size, which determines + the size of the generalized rectangle. + The simple function approximation is then defined by determining M Voronoi cells (i.e., "bins") partitioning :math:`\mathcal{D}`. These bins are only implicitly defined by M samples in :math:`\mathcal{D}`. - Finally, the probabilities of each of these bins is computed by - sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor - searches to bin these samples in the M implicitly defined bins. + Finally, the probabilities of each of these bins is computed by + sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor + searches to bin these samples in the M implicitly defined bins. The result is the simple function approximation denoted by :math:`\rho_{\mathcal{D},M}`. - + Note that all computations in the measure-theoretic framework that follow from this are for the fixed simple function approximation :math:`\rho_{\mathcal{D},M}`. @@ -45,16 +46,16 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - play around with it and you can get reasonable results with a relatively small number here like 50. - :param rect_scale: The scale used to determine the support of the + :param rect_size: The scale used to determine the support of the uniform distribution as ``rect_size = (data_max-data_min)*rect_scale`` - :type rect_scale: double or list() + :type rect_size: double or list() :param int num_d_emulate: Number of samples used to emulate using an MC - assumption + assumption :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param Q_ref: :math:`Q(`\lambda_{reference})` :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) - + :rtype: :class:`~bet.sample.voronoi_sample_set` :returns: sample_set object defininng simple function approximation """ @@ -65,7 +66,7 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec elif isinstance(data_set, samp.discretization): num = data_set.check_nums() dim = data_set._output_sample_set._dim - values = data_set._output_sample_set._values + values = data_set._output_sample_set._values elif isinstance(data_set, np.ndarray): num = data_set.shape[0] dim = data_set.shape[1] @@ -74,39 +75,41 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " "bet.sample.discretization or np.ndarray") - rect_size = (np.max(values, 0) - np.min(values, 0))*rect_scale - + if not isinstance(rect_size, collections.Iterable): + rect_size = rect_size * np.ones((dim,)) + if np.any(np.less(rect_size, 0)): + print 'Warning: rect_size must be greater than 0' r''' - Create M samples defining M Voronoi cells (i.e., "bins") in D used to + Create M samples defining M Voronoi cells (i.e., "bins") in D used to define the simple function approximation :math:`\rho_{\mathcal{D},M}`. - + This does not have to be random, but here we assume this to be the case. We can choose these samples deterministically but that fails to scale with dimension efficiently. - + Note that these M samples are chosen for the sole purpose of determining the bins used to create the approximation to :math:`rho_{\mathcal{D}}`. - + We call these M samples "d_distr_samples" because they are samples on the data space and the distr implies these samples are chosen to create the approximation to the probability measure (distribution) on D. - + Note that we create these samples in a set containing the hyperrectangle in order to get output cells with zero probability. If all of the d_dstr_samples were taken from within the support of :math:`\rho_{\mathcal{D}}` then each of the M bins would have positive probability. This would in turn imply that the support of :math:`\rho_{\Lambda}` is all of :math:`\Lambda`. - ''' + ''' if comm.rank == 0: - d_distr_samples = 1.5*rect_size*(np.random.random((M, - dim))-0.5)+Q_ref + d_distr_samples = 1.5 * rect_size * (np.random.random((M, + dim)) - 0.5) + Q_ref else: d_distr_samples = np.empty((M, dim)) comm.Bcast([d_distr_samples, MPI.DOUBLE], root=0) - + # Initialize sample set object s_set = samp.voronoi_sample_set(dim) s_set.set_values(d_distr_samples) @@ -120,9 +123,9 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec :math:`\rho_{\mathcal{D}}`. ''' # Generate the samples from :math:`\rho_{\mathcal{D}}` - num_d_emulate = int(num_d_emulate/comm.size)+1 - d_distr_emulate = rect_size*(np.random.random((num_d_emulate, - dim))-0.5) + Q_ref + num_d_emulate = int(num_d_emulate / comm.size) + 1 + d_distr_emulate = rect_size * (np.random.random((num_d_emulate, + dim)) - 0.5) + Q_ref # Bin these samples using nearest neighbor searches (_, k) = s_set.query(d_distr_emulate) @@ -131,14 +134,13 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec for i in range(M): count_neighbors[i] = np.sum(np.equal(k, i)) - # Use the binning to define :math:`\rho_{\mathcal{D},M}` ccount_neighbors = np.copy(count_neighbors) comm.Allreduce([count_neighbors, MPI.INT], [ccount_neighbors, MPI.INT], - op=MPI.SUM) + op=MPI.SUM) count_neighbors = ccount_neighbors rho_D_M = count_neighbors.astype(np.float64) / \ - float(num_d_emulate*comm.size) + float(num_d_emulate * comm.size) s_set.set_probabilities(rho_D_M) ''' @@ -152,198 +154,104 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec data_set._output_probability_set = s_set return s_set -def normal_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate=1E6): +def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, + rect_scale=0.2, M=50, + num_d_emulate=1E6): r""" - Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` - where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability - density centered at Q_ref with standard deviation std using M bins sampled - from the given normal distribution. - - :param data_set: Sample set that the probability measure is defined for. - :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` - :param int M: Defines number M samples in D used to define - :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - - play around with it and you can get reasonable results with a - relatively small number here like 50. - :param int num_d_emulate: Number of samples used to emulate using an MC - assumption - :param Q_ref: :math:`Q(\lambda_{reference})` - :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) - :param std: The standard deviation of each QoI - :type std: :class:`~numpy.ndarray` of size (mdim,) - - :rtype: :class:`~bet.sample.voronoi_sample_set` - :returns: sample_set object defining simple function approximation - - """ - import scipy.stats as stats - r'''Create M smaples defining M bins in D used to define - :math:`\rho_{\mathcal{D},M}` rho_D is assumed to be a multi-variate normal - distribution with mean Q_ref and standard deviation std.''' - if not isinstance(Q_ref, collections.Iterable): - Q_ref = np.array([Q_ref]) - if not isinstance(std, collections.Iterable): - std = np.array([std]) - - covariance = std**2 - - d_distr_samples = np.zeros((M, len(Q_ref))) - print "d_distr_samples.shape", d_distr_samples.shape - print "Q_ref.shape", Q_ref.shape - print "std.shape", std.shape - - if comm.rank == 0: - for i in range(len(Q_ref)): - d_distr_samples[:, i] = np.random.normal(Q_ref[i], std[i], M) - comm.Bcast([d_distr_samples, MPI.DOUBLE], root=0) - - # Initialize sample set object - s_set = samp.voronoi_sample_set(len(Q_ref)) - s_set.set_values(d_distr_samples) - s_set.set_kdtree() - - - r'''Now compute probabilities for :math:`\rho_{\mathcal{D},M}` by sampling - from rho_D First generate samples of rho_D - I sometimes call this - emulation''' - num_d_emulate = int(num_d_emulate/comm.size)+1 - d_distr_emulate = np.zeros((num_d_emulate, len(Q_ref))) - for i in range(len(Q_ref)): - d_distr_emulate[:, i] = np.random.normal(Q_ref[i], std[i], - num_d_emulate) - - # Now bin samples of rho_D in the M bins of D to compute rho_{D, M} - if len(d_distr_samples.shape) == 1: - d_distr_samples = np.expand_dims(d_distr_samples, axis=1) - - (_, k) = s_set.query(d_distr_emulate) - count_neighbors = np.zeros((M,), dtype=np.int) - volumes = np.zeros((M,)) - for i in range(M): - Itemp = np.equal(k, i) - count_neighbors[i] = np.sum(Itemp) - volumes[i] = np.sum(1.0/stats.multivariate_normal.pdf\ - (d_distr_emulate[Itemp, :], Q_ref, covariance)) - # Now define probability of the d_distr_samples - # This together with d_distr_samples defines :math:`\rho_{\mathcal{D},M}` - ccount_neighbors = np.copy(count_neighbors) - comm.Allreduce([count_neighbors, MPI.INT], [ccount_neighbors, MPI.INT], - op=MPI.SUM) - count_neighbors = ccount_neighbors - cvolumes = np.copy(volumes) - comm.Allreduce([volumes, MPI.DOUBLE], [cvolumes, MPI.DOUBLE], op=MPI.SUM) - volumes = cvolumes - rho_D_M = count_neighbors.astype(np.float64)*volumes - rho_D_M = rho_D_M/np.sum(rho_D_M) - s_set.set_probabilities(rho_D_M) - s_set.set_volumes(volumes) + Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` + where :math:`\rho_{\mathcal{D}}` is a uniform probability density on + a generalized rectangle centered at Q_ref. + The support of this density is defined by rect_scale, which determines + the size of the generalized rectangle by scaling the circumscribing + generalized rectangle of :math:`\mathcal{D}`. + The simple function approximation is then defined by determining M + Voronoi cells (i.e., "bins") partitioning :math:`\mathcal{D}`. These + bins are only implicitly defined by M samples in :math:`\mathcal{D}`. + Finally, the probabilities of each of these bins is computed by + sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor + searches to bin these samples in the M implicitly defined bins. + The result is the simple function approximation denoted by + :math:`\rho_{\mathcal{D},M}`. - # NOTE: The computation of q_distr_prob, q_distr_emulate, q_distr_samples - # above, while informed by the sampling of the map Q, do not require - # solving the model EVER! This can be done "offline" so to speak. - if isinstance(data_set, samp.discretization): - data_set._output_sample_set = s_set - return s_set - -def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate=1E6): - r""" - Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` - where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability - density centered at Q_ref with standard deviation std using M bins sampled - from a uniform distribution with a size 4 standard deviations in each - direction. + Note that all computations in the measure-theoretic framework that + follow from this are for the fixed simple function approximation + :math:`\rho_{\mathcal{D},M}`. - :param data_set: Sample set that the probability measure is defined for. - :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param int M: Defines number M samples in D used to define :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - play around with it and you can get reasonable results with a relatively small number here like 50. + :param rect_scale: The scale used to determine the support of the + uniform distribution as ``rect_size = (data_max-data_min)*rect_scale`` + :type rect_scale: double or list() :param int num_d_emulate: Number of samples used to emulate using an MC assumption - :param Q_ref: :math:`Q(\lambda_{reference})` + :param data_set: Sample set that the probability measure is defined for. + :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :param Q_ref: :math:`Q(`\lambda_{reference})` :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) - :param std: The standard deviation of each QoI - :type std: :class:`~numpy.ndarray` of size (mdim,) :rtype: :class:`~bet.sample.voronoi_sample_set` :returns: sample_set object defininng simple function approximation - """ - r'''Create M samples defining M bins in D used to define - :math:`\rho_{\mathcal{D},M}` rho_D is assumed to be a multi-variate normal - distribution with mean Q_ref and standard deviation std.''' - - bin_size = 4.0*std - d_distr_samples = np.zeros((M, len(Q_ref))) - if comm.rank == 0: - d_distr_samples = bin_size*(np.random.random((M, - len(Q_ref)))-0.5)+Q_ref - comm.Bcast([d_distr_samples, MPI.DOUBLE], root=0) - - # Initialize sample set object - s_set = samp.voronoi_sample_set(len(Q_ref)) - s_set.set_values(d_distr_samples) - s_set.set_kdtree() - - r'''Now compute probabilities for :math:`\rho_{\mathcal{D},M}` by sampling - from rho_D First generate samples of rho_D - I sometimes call this - emulation''' - num_d_emulate = int(num_d_emulate/comm.size)+1 - d_distr_emulate = np.zeros((num_d_emulate, len(Q_ref))) - for i in range(len(Q_ref)): - d_distr_emulate[:, i] = np.random.normal(Q_ref[i], std[i], - num_d_emulate) + if isinstance(data_set, samp.sample_set_base): + num = data_set.check_num() + dim = data_set._dim + values = data_set._values + elif isinstance(data_set, samp.discretization): + num = data_set.check_nums() + dim = data_set._output_sample_set._dim + values = data_set._output_sample_set._values + elif isinstance(data_set, np.ndarray): + num = data_set.shape[0] + dim = data_set.shape[1] + values = data_set + else: + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") - # Now bin samples of rho_D in the M bins of D to compute rho_{D, M} - if len(d_distr_samples.shape) == 1: - d_distr_samples = np.expand_dims(d_distr_samples, axis=1) + rect_size = (np.max(values, 0) - np.min(values, 0))*rect_scale - (_, k) = s_set.query(d_distr_emulate) - count_neighbors = np.zeros((M,), dtype=np.int) - #volumes = np.zeros((M,)) - for i in range(M): - Itemp = np.equal(k, i) - count_neighbors[i] = np.sum(Itemp) - - r'''Now define probability of the d_distr_samples This together with - d_distr_samples defines :math:`\rho_{\mathcal{D},M}`''' - ccount_neighbors = np.copy(count_neighbors) - comm.Allreduce([count_neighbors, MPI.INT], [ccount_neighbors, MPI.INT], - op=MPI.SUM) - count_neighbors = ccount_neighbors - rho_D_M = count_neighbors.astype(np.float64)/float(comm.size*num_d_emulate) - s_set.set_probabilities(rho_D_M) - # NOTE: The computation of q_distr_prob, q_distr_emulate, q_distr_samples - # above, while informed by the sampling of the map Q, do not require - # solving the model EVER! This can be done "offline" so to speak. - if isinstance(data_set, samp.discretization): - data_set._output_probability_set = s_set - return s_set + return uniform_partition_uniform_distribution_rectangle_size(data_set, Q_ref, + rect_size, M, + num_d_emulate) -def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domain, center_pts_per_edge=1): +def uniform_partition_uniform_distribution_rectangle_domain(data_set, rect_domain, + M=50, num_d_emulate=1E6): r""" - Creates a simple function appoximation of :math:`\rho_{\mathcal{D},M}` - where :math:`\rho{\mathcal{D}, M}` is a uniform probablity density over the - hyperrectangular domain specified by domain. + Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` + where :math:`\rho_{\mathcal{D}}` is a uniform probability density on + a generalized rectangle defined by rect_domain. + The simple function approximation is then defined by determining M + Voronoi cells (i.e., "bins") partitioning :math:`\mathcal{D}`. These + bins are only implicitly defined by M samples in :math:`\mathcal{D}`. + Finally, the probabilities of each of these bins is computed by + sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor + searches to bin these samples in the M implicitly defined bins. + The result is the simple function approximation denoted by + :math:`\rho_{\mathcal{D},M}`. - Since :math:`\rho_\mathcal{D}` is a uniform distribution on a - hyperrectangle we should we able to represent it exactly with - :math:`M=3^{m}` where m is the dimension of the data space or rather - ``len(d_distr_samples) == 3**mdim``. + Note that all computations in the measure-theoretic framework that + follow from this are for the fixed simple function approximation + :math:`\rho_{\mathcal{D},M}`. + :param int M: Defines number M samples in D used to define + :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - + play around with it and you can get reasonable results with a + relatively small number here like 50. + :param rect_domain: The support of the density + :type rect_domain: double or list() + :param int num_d_emulate: Number of samples used to emulate using an MC + assumption :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` - :param rect_domain: The domain overwhich :math:`\rho_\mathcal{D}` is - uniform. - :type rect_domain: :class:`numpy.ndarray` of shape (2, mdim) - :param list() center_pts_per_edge: number of center points per edge and - additional two points will be added to create the bounding layer + :param Q_ref: :math:`Q(`\lambda_{reference})` + :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) :rtype: :class:`~bet.sample.voronoi_sample_set` :returns: sample_set object defininng simple function approximation - """ + # make sure the shape of the data and the domain are correct if isinstance(data_set, samp.sample_set_base): num = data_set.check_num() @@ -352,7 +260,7 @@ def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domai elif isinstance(data_set, samp.discretization): num = data_set.check_nums() dim = data_set._output_sample_set._dim - values = data_set._output_sample_set._values + values = data_set._output_sample_set._values elif isinstance(data_set, np.ndarray): num = data_set.shape[0] dim = data_set.shape[1] @@ -361,34 +269,36 @@ def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domai raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " "bet.sample.discretization or np.ndarray") - data = values + data = values rect_domain = util.fix_dimensions_data(rect_domain, data.shape[1]) domain_center = np.mean(rect_domain, 0) domain_lengths = np.max(rect_domain, 0) - np.min(rect_domain, 0) - - return regular_partition_uniform_distribution_rectangle_size(data_set, domain_center, domain_lengths, - center_pts_per_edge) + + return uniform_partition_uniform_distribution_rectangle_size(data_set, domain_center, + domain_lengths, M, + num_d_emulate) + def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_size, - center_pts_per_edge=1): + center_pts_per_edge=1): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density - centered at Q_ref with bin_size of the width of D. + centered at Q_ref with rect_size of the width of a hyperrectangle. Since rho_D is a uniform distribution on a hyperrectanlge we should be able to represent it exactly with ``M = 3^mdim`` or rather ``len(d_distr_samples) == 3^mdim``. :param rect_size: The size used to determine the width of the uniform - distribution + distribution :type rect_size: double or list() - :param int num_d_emulate: Number of samples used to emulate using an MC - assumption + :param int num_d_emulate: Number of samples used to emulate using an MC + assumption :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` - :param Q_ref: :math:`Q(\lambda_{reference})` - :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) + :param Q_ref: :math:`Q(\lambda_{reference})` + :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) :param list() center_pts_per_edge: number of center points per edge and additional two points will be added to create the bounding layer @@ -404,7 +314,7 @@ def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_ elif isinstance(data_set, samp.discretization): num = data_set.check_nums() dim = data_set._output_sample_set._dim - values = data_set._output_sample_set._values + values = data_set._output_sample_set._values elif isinstance(data_set, np.ndarray): num = data_set.shape[0] dim = data_set.shape[1] @@ -412,35 +322,84 @@ def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_ else: raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " "bet.sample.discretization or np.ndarray") - + data = values if not isinstance(center_pts_per_edge, collections.Iterable): center_pts_per_edge = np.ones((dim,)) * center_pts_per_edge else: - if not len(center_pts_per_edge) == dim: + if not len(center_pts_per_edge) == dim: center_pts_per_edge = np.ones((dim,)) print 'Warning: center_pts_per_edge dimension mismatch.' print 'Using 1 in each dimension.' if np.any(np.less(center_pts_per_edge, 0)): print 'Warning: center_pts_per_edge must be greater than 0' if not isinstance(rect_size, collections.Iterable): - rect_size = rect_size*np.ones((dim,)) + rect_size = rect_size * np.ones((dim,)) if np.any(np.less(rect_size, 0)): - print 'Warning: center_pts_per_edge must be greater than 0' + print 'Warning: rect_size must be greater than 0' sur_domain = np.array([np.min(data, 0), np.max(data, 0)]).transpose() - points, _, rect_domain = vHist.center_and_layer1_points_binsize\ - (center_pts_per_edge, Q_ref, rect_size, sur_domain) - edges = vHist.edges_regular(center_pts_per_edge, rect_domain, sur_domain) + points, _, rect_domain = vHist.center_and_layer1_points_binsize \ + (center_pts_per_edge, Q_ref, rect_size, sur_domain) + edges = vHist.edges_regular(center_pts_per_edge, rect_domain, sur_domain) _, volumes, _ = vHist.histogramdd_volumes(edges, points) - s_set = vHist.simple_fun_uniform(points, volumes, rect_domain) + s_set = vHist.simple_fun_uniform(points, volumes, rect_domain) if isinstance(data_set, samp.discretization): data_set._output_probability_set = s_set return s_set + +def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domain, center_pts_per_edge=1): + r""" + Creates a simple function appoximation of :math:`\rho_{\mathcal{D},M}` + where :math:`\rho{\mathcal{D}, M}` is a uniform probablity density over the + hyperrectangular domain specified by domain. + + Since :math:`\rho_\mathcal{D}` is a uniform distribution on a + hyperrectangle we should we able to represent it exactly with + :math:`M=3^{m}` where m is the dimension of the data space or rather + ``len(d_distr_samples) == 3**mdim``. + + :param data_set: Sample set that the probability measure is defined for. + :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :param rect_domain: The domain overwhich :math:`\rho_\mathcal{D}` is + uniform. + :type rect_domain: :class:`numpy.ndarray` of shape (2, mdim) + :param list() center_pts_per_edge: number of center points per edge and + additional two points will be added to create the bounding layer + + :rtype: :class:`~bet.sample.voronoi_sample_set` + :returns: sample_set object defininng simple function approximation + + """ + # make sure the shape of the data and the domain are correct + if isinstance(data_set, samp.sample_set_base): + num = data_set.check_num() + dim = data_set._dim + values = data_set._values + elif isinstance(data_set, samp.discretization): + num = data_set.check_nums() + dim = data_set._output_sample_set._dim + values = data_set._output_sample_set._values + elif isinstance(data_set, np.ndarray): + num = data_set.shape[0] + dim = data_set.shape[1] + values = data_set + else: + raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, " + "bet.sample.discretization or np.ndarray") + + data = values + rect_domain = util.fix_dimensions_data(rect_domain, data.shape[1]) + domain_center = np.mean(rect_domain, 0) + domain_lengths = np.max(rect_domain, 0) - np.min(rect_domain, 0) + + return regular_partition_uniform_distribution_rectangle_size(data_set, domain_center, + domain_lengths, center_pts_per_edge) + def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rect_scale, center_pts_per_edge=1): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` @@ -535,3 +494,173 @@ def uniform_partition_uniform_distribution_data_samples(data_set): if isinstance(data_set, samp.discretization): data_set._output_sample_set = s_set return s_set + + +def normal_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate=1E6): + r""" + Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` + where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability + density centered at Q_ref with standard deviation std using M bins sampled + from the given normal distribution. + + :param data_set: Sample set that the probability measure is defined for. + :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :param int M: Defines number M samples in D used to define + :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - + play around with it and you can get reasonable results with a + relatively small number here like 50. + :param int num_d_emulate: Number of samples used to emulate using an MC + assumption + :param Q_ref: :math:`Q(\lambda_{reference})` + :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) + :param std: The standard deviation of each QoI + :type std: :class:`~numpy.ndarray` of size (mdim,) + + :rtype: :class:`~bet.sample.voronoi_sample_set` + :returns: sample_set object defining simple function approximation + + """ + import scipy.stats as stats + r'''Create M smaples defining M bins in D used to define + :math:`\rho_{\mathcal{D},M}` rho_D is assumed to be a multi-variate normal + distribution with mean Q_ref and standard deviation std.''' + if not isinstance(Q_ref, collections.Iterable): + Q_ref = np.array([Q_ref]) + if not isinstance(std, collections.Iterable): + std = np.array([std]) + + covariance = std ** 2 + + d_distr_samples = np.zeros((M, len(Q_ref))) + print "d_distr_samples.shape", d_distr_samples.shape + print "Q_ref.shape", Q_ref.shape + print "std.shape", std.shape + + if comm.rank == 0: + for i in range(len(Q_ref)): + d_distr_samples[:, i] = np.random.normal(Q_ref[i], std[i], M) + comm.Bcast([d_distr_samples, MPI.DOUBLE], root=0) + + # Initialize sample set object + s_set = samp.voronoi_sample_set(len(Q_ref)) + s_set.set_values(d_distr_samples) + s_set.set_kdtree() + + r'''Now compute probabilities for :math:`\rho_{\mathcal{D},M}` by sampling + from rho_D First generate samples of rho_D - I sometimes call this + emulation''' + num_d_emulate = int(num_d_emulate / comm.size) + 1 + d_distr_emulate = np.zeros((num_d_emulate, len(Q_ref))) + for i in range(len(Q_ref)): + d_distr_emulate[:, i] = np.random.normal(Q_ref[i], std[i], + num_d_emulate) + + # Now bin samples of rho_D in the M bins of D to compute rho_{D, M} + if len(d_distr_samples.shape) == 1: + d_distr_samples = np.expand_dims(d_distr_samples, axis=1) + + (_, k) = s_set.query(d_distr_emulate) + count_neighbors = np.zeros((M,), dtype=np.int) + volumes = np.zeros((M,)) + for i in range(M): + Itemp = np.equal(k, i) + count_neighbors[i] = np.sum(Itemp) + volumes[i] = np.sum(1.0 / stats.multivariate_normal.pdf \ + (d_distr_emulate[Itemp, :], Q_ref, covariance)) + # Now define probability of the d_distr_samples + # This together with d_distr_samples defines :math:`\rho_{\mathcal{D},M}` + ccount_neighbors = np.copy(count_neighbors) + comm.Allreduce([count_neighbors, MPI.INT], [ccount_neighbors, MPI.INT], + op=MPI.SUM) + count_neighbors = ccount_neighbors + cvolumes = np.copy(volumes) + comm.Allreduce([volumes, MPI.DOUBLE], [cvolumes, MPI.DOUBLE], op=MPI.SUM) + volumes = cvolumes + rho_D_M = count_neighbors.astype(np.float64) * volumes + rho_D_M = rho_D_M / np.sum(rho_D_M) + s_set.set_probabilities(rho_D_M) + s_set.set_volumes(volumes) + + # NOTE: The computation of q_distr_prob, q_distr_emulate, q_distr_samples + # above, while informed by the sampling of the map Q, do not require + # solving the model EVER! This can be done "offline" so to speak. + if isinstance(data_set, samp.discretization): + data_set._output_sample_set = s_set + return s_set + + +def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate=1E6): + r""" + Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` + where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability + density centered at Q_ref with standard deviation std using M bins sampled + from a uniform distribution with a size 4 standard deviations in each + direction. + + :param data_set: Sample set that the probability measure is defined for. + :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :param int M: Defines number M samples in D used to define + :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - + play around with it and you can get reasonable results with a + relatively small number here like 50. + :param int num_d_emulate: Number of samples used to emulate using an MC + assumption + :param Q_ref: :math:`Q(\lambda_{reference})` + :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) + :param std: The standard deviation of each QoI + :type std: :class:`~numpy.ndarray` of size (mdim,) + + :rtype: :class:`~bet.sample.voronoi_sample_set` + :returns: sample_set object defininng simple function approximation + + """ + r'''Create M samples defining M bins in D used to define + :math:`\rho_{\mathcal{D},M}` rho_D is assumed to be a multi-variate normal + distribution with mean Q_ref and standard deviation std.''' + + bin_size = 4.0 * std + d_distr_samples = np.zeros((M, len(Q_ref))) + if comm.rank == 0: + d_distr_samples = bin_size * (np.random.random((M, + len(Q_ref))) - 0.5) + Q_ref + comm.Bcast([d_distr_samples, MPI.DOUBLE], root=0) + + # Initialize sample set object + s_set = samp.voronoi_sample_set(len(Q_ref)) + s_set.set_values(d_distr_samples) + s_set.set_kdtree() + + r'''Now compute probabilities for :math:`\rho_{\mathcal{D},M}` by sampling + from rho_D First generate samples of rho_D - I sometimes call this + emulation''' + num_d_emulate = int(num_d_emulate / comm.size) + 1 + d_distr_emulate = np.zeros((num_d_emulate, len(Q_ref))) + for i in range(len(Q_ref)): + d_distr_emulate[:, i] = np.random.normal(Q_ref[i], std[i], + num_d_emulate) + + # Now bin samples of rho_D in the M bins of D to compute rho_{D, M} + if len(d_distr_samples.shape) == 1: + d_distr_samples = np.expand_dims(d_distr_samples, axis=1) + + (_, k) = s_set.query(d_distr_emulate) + count_neighbors = np.zeros((M,), dtype=np.int) + # volumes = np.zeros((M,)) + for i in range(M): + Itemp = np.equal(k, i) + count_neighbors[i] = np.sum(Itemp) + + r'''Now define probability of the d_distr_samples This together with + d_distr_samples defines :math:`\rho_{\mathcal{D},M}`''' + ccount_neighbors = np.copy(count_neighbors) + comm.Allreduce([count_neighbors, MPI.INT], [ccount_neighbors, MPI.INT], + op=MPI.SUM) + count_neighbors = ccount_neighbors + rho_D_M = count_neighbors.astype(np.float64) / float(comm.size * num_d_emulate) + s_set.set_probabilities(rho_D_M) + # NOTE: The computation of q_distr_prob, q_distr_emulate, q_distr_samples + # above, while informed by the sampling of the map Q, do not require + # solving the model EVER! This can be done "offline" so to speak. + if isinstance(data_set, samp.discretization): + data_set._output_probability_set = s_set + return s_set From d458e3e413a40f33a744090742293e60865cba2e Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Thu, 19 May 2016 18:50:47 -0600 Subject: [PATCH 03/20] Added regular sampling, updates to linear example --- bet/postProcess/plotDomains.py | 25 ++--- bet/sampling/basicSampling.py | 103 +++++++++++++++++- .../linearMap/linearMapUniformSampling.py | 83 ++++++++++---- 3 files changed, 173 insertions(+), 38 deletions(-) diff --git a/bet/postProcess/plotDomains.py b/bet/postProcess/plotDomains.py index 772b1c8f..64eb0d99 100644 --- a/bet/postProcess/plotDomains.py +++ b/bet/postProcess/plotDomains.py @@ -264,7 +264,7 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, the other markers. :param sample_obj: Object containing the samples to plot - :type sample_obj: :class:`~bet.sample.sample_set` + :type sample_obj: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization` :param list sample_nos: sample numbers to plot :param rho_D: probability density on D :type rho_D: callable function that takes a :class:`np.array` and returns a @@ -280,7 +280,9 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, pairwise or tripletwise data sample scatter plots in 2 or 3 dimensions - """ + """ + if isinstance(sample_obj, sample.discretization): + sample_obj = sample_obj._output_sample_set # If there is density function given determine the pointwise probability # values of each sample based on the value in the data space. Otherwise, # color the samples in numerical order. @@ -522,7 +524,6 @@ def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None, def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True, interactive=False): r""" - Creates two-dimensional projections of scatter plots of samples. :param sample_obj: Object containing the samples to plot @@ -555,10 +556,9 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True xlabel = r'$\lambda_{'+str(showdim+1)+r'}$' ylabel = r'$\lambda_{'+str(i+1)+r'}$' - filenames = [img_folder+'domain_l'+str(showdim+1)+'_l'+\ - str(i+1)+'.eps', img_folder+'l'+str(showdim+1)+\ - '_l'+str(i+1)+'_domain_L_cs.eps'] - filename = filenames[0] + filename = [img_folder+'domain_l'+str(showdim+1)+'_l'+\ + str(i+1)+'.eps'] + plt.scatter(sample_obj.get_values()[:, 0], sample_obj.get_values()[:, 1]) if save: plt.autoscale(tight=True) @@ -576,10 +576,9 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True xlabel = r'$\lambda_{'+str(x+1)+r'}$' ylabel = r'$\lambda_{'+str(y+1)+r'}$' - filenames = [img_folder+'domain_l'+str(x+1)+'_l'+\ - str(y+1)+'.eps', img_folder+'l'+str(x+1)+\ - '_l'+str(y+1)+'_domain_L_cs.eps'] - filename = filenames[0] + filename = [img_folder+'domain_l'+str(x+1)+'_l'+\ + str(y+1)+'.eps'] + plt.scatter(sample_obj.get_values()[:, x], sample_obj.get_values()[:, y]) if save: plt.autoscale(tight=True) @@ -641,7 +640,7 @@ def scatter2D_multi(sample_obj, color=None, p_ref=None, img_folder='figs/', sample_obj_temp = sample.sample_set(2) sample_obj_temp.set_values(sample_obj.get_values()[:, [showdim, i]]) - if p_ref: + if p_ref is not None: scatter_2D(sample_obj_temp, sample_nos=None, color=color, p_ref=p_ref[[showdim, i]], save=True, interactive=False, xlabel=xlabel, ylabel=ylabel, @@ -664,7 +663,7 @@ def scatter2D_multi(sample_obj, color=None, p_ref=None, img_folder='figs/', sample_obj_temp = sample.sample_set(2) sample_obj_temp.set_values(sample_obj.get_values()[:, [x, y]]) - if p_ref: + if p_ref is not None: scatter_2D(sample_obj_temp, sample_nos=None, color=color, p_ref=p_ref[[x, y]], save=True, interactive=False, xlabel=xlabel, ylabel=ylabel, filename=myfilename) diff --git a/bet/sampling/basicSampling.py b/bet/sampling/basicSampling.py index e0a3082b..ec84b857 100644 --- a/bet/sampling/basicSampling.py +++ b/bet/sampling/basicSampling.py @@ -15,6 +15,8 @@ from pyDOE import lhs from bet.Comm import comm import bet.sample as sample +import bet.util as util +import collections def loadmat(save_file, disc_name=None, model=None): """ @@ -218,7 +220,106 @@ def random_sample_set_dimension(self, sample_type, input_dim, return self.random_sample_set(sample_type, input_sample_set, num_samples, criterion, parallel) - def compute_QoI_and_create_discretization(self, input_sample_set, savefile=None, parallel=False): + def regular_sample_set(self, input_sample_set, num_samples_per_dim=1): + """ + Sampling algorithm for generating a regular grid of samples taken + on the domain present with input_sample_set (a default unit hypercube + is used if no domain has been specified) + + :param input_sample_set: samples to evaluate the model at + :type input_sample_set: :class:`~bet.sample.sample_set` with + num_smaples + :param num_samples_per_dim: number of samples per dimension + :type num_samples_per_dim: :class: `~numpy.ndarray` of dimension + (input_sample_set._dim,) + + :rtype: :class:`~bet.sample.discretization` + :returns: :class:`~bet.sample.discretization` object which contains + input and output of ``num_samples`` + """ + + # Create N samples + dim = input_sample_set.get_dim() + + if not isinstance(num_samples_per_dim, collections.Iterable): + num_samples_per_dim = num_samples_per_dim * np.ones((dim,)) + if np.any(np.less(num_samples_per_dim, 0)): + print 'Warning: num_smaples_per_dim must be greater than 0' + + self.num_samples = np.product(num_samples_per_dim) + + if input_sample_set.get_domain() is None: + # create the domain + input_domain = np.array([[0., 1.]] * dim) + input_sample_set.set_domain(input_domain) + else: + input_domain = input_sample_set.get_domain() + # update the bounds based on the number of samples + input_sample_set.update_bounds(self.num_samples) + input_values = np.copy(input_sample_set._width) + + vec_samples_dimension = np.empty((dim), dtype=object) + for i in np.arange(0, dim): + vec_samples_dimension[i] = list(np.linspace( + input_domain[i,0], input_domain[i,1], + num_samples_per_dim[i]+2))[1:num_samples_per_dim[i]+1] + + arrays_samples_dimension = np.meshgrid( + *[vec_samples_dimension[i] for i in np.arange(0, dim)], indexing='ij') + + for i in np.arange(0, dim): + input_values[:,i:i+1] = np.vstack(arrays_samples_dimension[i].flat[:]) + + input_sample_set.set_values(input_values) + + return input_sample_set + + def regular_sample_set_domain(self, input_domain, num_samples_per_dim=1): + """ + Sampling algorithm for generating a regular grid of samples taken + on the domain present with input_sample_set (a default unit hypercube + is used if no domain has been specified) + + :param input_domain: min and max bounds for the input values, + ``min = input_domain[:, 0]`` and ``max = input_domain[:, 1]`` + :type input_domain: :class:`numpy.ndarray` of shape (ndim, 2) + :param num_samples_per_dim: number of samples per dimension + :type num_samples_per_dim: :class: `~numpy.ndarray` of dimension + (input_sample_set._dim,) + + :rtype: :class:`~bet.sample.discretization` + :returns: :class:`~bet.sample.discretization` object which contains + input and output of ``num_samples`` + + """ + # Create N samples + input_sample_set = sample.sample_set(input_domain.shape[0]) + input_sample_set.set_domain(input_domain) + + return self.regular_sample_set(input_sample_set, num_samples_per_dim) + + def regular_sample_set_dimension(self, input_dim, num_samples_per_dimension): + """ + Sampling algorithm for generating a regular grid of samples taken + on a unit hypercube of dimension input_dim + + :param int input_dim: the dimension of the input space + :param num_samples_per_dim: number of samples per dimension + :type num_samples_per_dim: :class: `~numpy.ndarray` of dimension + (input_sample_set._dim,) + + :rtype: :class:`~bet.sample.discretization` + :returns: :class:`~bet.sample.discretization` object which contains + input and output of ``num_samples`` + + """ + # Create N samples + input_sample_set = sample.sample_set(input_dim) + + return self.regular_sample_set(input_sample_set, num_samples_per_dimension) + + def compute_QoI_and_create_discretization(self, input_sample_set, savefile=None, + parallel=False): """ Samples the model at ``input_sample_set`` and saves the results. diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index b91e5c23..54bd539b 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -30,6 +30,7 @@ import bet.calculateP.simpleFunP as simpleFunP import bet.calculateP.calculateP as calculateP import bet.postProcess.plotP as plotP +import bet.postProcess.plotDomains as plotD import bet.sample as samp import bet.sampling.basicSampling as bsam from myModel import my_model @@ -43,39 +44,67 @@ # Define the sampler that will be used to create the discretization # object, which is the fundamental object used by BET to compute # solutions to the stochastic inverse problem +sampler = bsam.sampler(my_model) + ''' Suggested changes for user: -Try num_samples = 1E3 and 1E4. What happens when num_samples = 1E2? -Try using 'lhs' instead of 'random' in the sampler. -''' -sampler = bsam.sampler(my_model, num_samples=1E4) +Try with and without random sampling. + +If using random sampling, try num_samples = 1E3 and 1E4. +What happens when num_samples = 1E2? +Try using 'lhs' instead of 'random' in the random_sample_set. +If using regular sampling, try different numbers of samples +per dimension. +''' # Generate samples on the parameter space -input_samples = sampler.random_sample_set('random', input_samples) -# Estimate volumes of Voronoi cells associated with the parameter samples +randomSampling = False +if randomSampling is True: + sampler.random_sample_set('random', input_samples, num_samples=1E3) +else: + sampler.regular_sample_set(input_samples, num_samples_per_dim=[15, 15, 10]) + ''' Suggested changes for user: -Try different numbers of points to estimate the volume, and also -try comparing to the standard Monte Carlo assumption that all the -Voronoi cells have the same measure. +A standard Monte Carlo (MC) assumption is that every Voronoi cell +has the same volume. If a regular grid of samples was used, then +the standard MC assumption is true. + +See what happens if the MC assumption is not assumed to be true, and +if different numbers of points are used to estimate the volumes of +the Voronoi cells. ''' MC_assumption = False +# Estimate volumes of Voronoi cells associated with the parameter samples if MC_assumption is False: input_samples.estimate_volume(n_mc_points=1E5) else: input_samples.estimate_volume_mc() -# Now create the discretization object using the input samples -my_discretization = sampler.create_random_discretization('random', input_samples, +# Create the discretization object using the input samples +my_discretization = sampler.compute_QoI_and_create_discretization(input_samples, savefile = '3to2_discretization.txt.gz') +''' +Suggested changes for user: + +Try different reference parameters. +''' # Define the reference parameter param_ref = np.array([0.5, 0.5, 0.5]) +#param_ref = np.array([0.75, 0.75, 0.5]) +#param_ref = np.array([0.75, 0.75, 0.75]) +#param_ref = np.array([0.5, 0.5, 0.75]) + # Compute the reference QoI Q_ref = my_model(param_ref) +# Create some plots of input and output discretizations +plotD.scatter2D_multi(input_samples, p_ref = param_ref, showdim = 'all', filename = 'linearMapParameterSamples') +plotD.show_data(my_discretization, Q_ref = Q_ref) + ''' Suggested changes for user: @@ -83,17 +112,22 @@ probability measure on a rectangle (since D is 2-dimensional) centered at Q_ref whose size is determined by scaling the circumscribing box of D. ''' -simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( - data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.2, - center_pts_per_edge = 5) -#simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled( -# data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.2, -# M=50, num_d_emulate=1E5) +randomDataDiscretization = False +if randomDataDiscretization is False: + simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.25, + center_pts_per_edge = 3) +else: + simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.25, + M=50, num_d_emulate=1E5) # calculate probablities calculateP.prob(my_discretization) -# calculate 2d marginal probs +######################################## +# Post-process the results +######################################## ''' Suggested changes for user: @@ -105,13 +139,14 @@ There are ways to determine "optimal" smoothing parameters (e.g., see CV, GCV, and other similar methods), but we have not incorporated these into the code -as lower-dimensional marginal plots have limited value in understanding the -structure of a high dimensional non-parametric probability measure. +as lower-dimensional marginal plots generally have limited value in understanding +the structure of a high dimensional non-parametric probability measure. ''' +# calculate 2d marginal probs (bins, marginals2D) = plotP.calculate_2D_marginal_probs(input_samples, - nbins = 30) + nbins = [10, 10, 10]) # smooth 2d marginals probs (optional) -marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.1) +marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.2) # plot 2d marginals probs plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "linearMap", @@ -119,9 +154,9 @@ # calculate 1d marginal probs (bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, - nbins = 40) + nbins = [10, 10, 10]) # smooth 1d marginal probs (optional) -marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.1) +marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.2) # plot 2d marginal probs plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "linearMap", file_extension = ".eps") From 4d571e12e819734b0791569feabd08ac429283a4 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Fri, 20 May 2016 08:09:02 -0600 Subject: [PATCH 04/20] Final linear example --- bet/postProcess/plotDomains.py | 238 +++++++++--------- .../linearMap/linearMapUniformSampling.py | 2 +- 2 files changed, 127 insertions(+), 113 deletions(-) diff --git a/bet/postProcess/plotDomains.py b/bet/postProcess/plotDomains.py index 64eb0d99..2aeec9b5 100644 --- a/bet/postProcess/plotDomains.py +++ b/bet/postProcess/plotDomains.py @@ -9,8 +9,8 @@ import matplotlib.tri as tri import numpy as np import matplotlib.pyplot as plt -#plt.rc('text', usetex=True) -#plt.rc('font', family='serif') +# plt.rc('text', usetex=True) +# plt.rc('font', family='serif') from matplotlib.lines import Line2D from itertools import combinations from mpl_toolkits.mplot3d import Axes3D @@ -28,19 +28,22 @@ colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k') + class dim_not_matching(Exception): """ Exception for when the dimension is inconsistent. """ + class bad_object(Exception): """ Exception for when the wrong type of object is used. """ + def scatter_2D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True, interactive=False, xlabel='x', ylabel='y', - filename='scatter2d'): + filename='scatter2d'): r""" Creates a two-dimensional scatter plot of the samples within the sample object colored by ``color`` (usually an array of pointwise probability density values). @@ -92,18 +95,19 @@ def scatter_2D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True, cbar.set_label(r'$\rho_\mathcal{D}(q)$') # if there is a reference value plot it with a notiable marker if p_ref is not None: - plt.scatter(p_ref[0], p_ref[1], c='m', s=2*markersize) + plt.scatter(p_ref[0], p_ref[1], c='m', s=2 * markersize) if save: plt.autoscale(tight=True) plt.xlabel(xlabel) plt.ylabel(ylabel) plt.savefig(filename, bbox_inches='tight', transparent=True, - pad_inches=0) + pad_inches=0) if interactive: plt.show() else: plt.close() + def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True, interactive=False, xlabel='x', ylabel='y', zlabel='z', filename="scatter3d"): @@ -113,7 +117,7 @@ def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True, A reference sample (``p_ref``) can be chosen by the user. This reference sample will be plotted as a mauve circle twice the size of the other markers. - + :param sample_obj: Object containing the samples to plot :type sample_obj: :class:`~bet.sample.sample_set` :param list sample_nos: indicies of the samples to plot @@ -159,24 +163,25 @@ def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True, # samples are colored by the pointwise probability density on the data # space cbar = fig.colorbar(p) - cbar.set_label(r'$\rho_\mathcal{D}(q)$') + cbar.set_label(r'$\rho_\mathcal{D}(q)$') # if there is a reference value plot it with a notiable marker if p_ref is not None: - ax.scatter(p_ref[0], p_ref[1], p_ref[2], c='m', s=2*markersize) + ax.scatter(p_ref[0], p_ref[1], p_ref[2], c='m', s=2 * markersize) ax.autoscale(tight=True) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_zlabel(zlabel) if save: plt.savefig(filename, bbox_inches='tight', transparent=True, - pad_inches=0) + pad_inches=0) if interactive: plt.show() else: plt.close() - + + def show_param(sample_disc, rho_D=None, p_ref=None, sample_nos=None, - save=True, interactive=False, lnums=None, showdim=None): + save=True, interactive=False, lnums=None, showdim=None): r""" Create scatter plots of samples within the sample object colored by ``color`` (usually an array of pointwise probability density values). @@ -220,42 +225,43 @@ def show_param(sample_disc, rho_D=None, p_ref=None, sample_nos=None, # (e.g. i, where \lambda_i is a coordinate in the parameter space), then # set them to be the the counting numbers. if lnums is None: - lnums = 1+np.array(range(sample_obj.get_values().shape[1])) + lnums = 1 + np.array(range(sample_obj.get_values().shape[1])) # Create the labels based on the user selected parameter coordinates - xlabel = r'$\lambda_{'+str(lnums[0])+'}$' - ylabel = r'$\lambda_{'+str(lnums[1])+'}$' + xlabel = r'$\lambda_{' + str(lnums[0]) + '}$' + ylabel = r'$\lambda_{' + str(lnums[1]) + '}$' savename = 'param_samples_cs.eps' # Plot 2 or 3 dimensional scatter plots of the samples colored by rD. if sample_obj.get_dim() == 2: scatter_2D(sample_obj, sample_nos, rD, p_ref, save, interactive, xlabel, ylabel, savename) elif sample_obj.get_dim() == 3: - zlabel = r'$\lambda_{'+str(lnums[2])+'}$' + zlabel = r'$\lambda_{' + str(lnums[2]) + '}$' scatter_3D(sample_obj, sample_nos, rD, p_ref, save, interactive, xlabel, ylabel, zlabel, savename) elif sample_obj.get_dim() > 2 and showdim == 2: temp_obj = sample.sample_set(2) for x, y in combinations(lnums, 2): - xlabel = r'$\lambda_{'+str(x)+'}$' - ylabel = r'$\lambda_{'+str(y)+'}$' - savename = 'param_samples_l'+str(x)+'l'+str(y)+'_cs.eps' - temp_obj.set_values(sample_obj.get_values()[:, [x-1, y-1]]) + xlabel = r'$\lambda_{' + str(x) + '}$' + ylabel = r'$\lambda_{' + str(y) + '}$' + savename = 'param_samples_l' + str(x) + 'l' + str(y) + '_cs.eps' + temp_obj.set_values(sample_obj.get_values()[:, [x - 1, y - 1]]) scatter_2D(temp_obj, sample_nos, rD, p_ref, save, interactive, xlabel, ylabel, savename) elif sample_obj.get_dim() > 3 and showdim == 3: temp_obj = sample.sample_set(3) for x, y, z in combinations(lnums, 3): - xlabel = r'$\lambda_{'+str(x)+'}$' - ylabel = r'$\lambda_{'+str(y)+'}$' - zlabel = r'$\lambda_{'+str(z)+'}$' - savename = 'param_samples_l'+str(x)+'l'+str(y)+'l'+str(z)+\ + xlabel = r'$\lambda_{' + str(x) + '}$' + ylabel = r'$\lambda_{' + str(y) + '}$' + zlabel = r'$\lambda_{' + str(z) + '}$' + savename = 'param_samples_l' + str(x) + 'l' + str(y) + 'l' + str(z) + \ '_cs.eps' - temp_obj.set_values(sample_obj.get_values()[:, [x-1, y-1, z-1]]) + temp_obj.set_values(sample_obj.get_values()[:, [x - 1, y - 1, z - 1]]) scatter_3D(temp_obj, sample_nos, rD, p_ref, save, interactive, xlabel, ylabel, zlabel, savename) + def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, - save=True, interactive=False, Q_nums=None, showdim=None): + save=True, interactive=False, Q_nums=None, showdim=None): r""" Create scatter plots of data within the sample_obj colored by ``color`` (usually an array of pointwise probability density values). @@ -264,7 +270,7 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, the other markers. :param sample_obj: Object containing the samples to plot - :type sample_obj: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization` + :type sample_obj: :class:`~bet.sample.sample_set` :param list sample_nos: sample numbers to plot :param rho_D: probability density on D :type rho_D: callable function that takes a :class:`np.array` and returns a @@ -283,6 +289,7 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, """ if isinstance(sample_obj, sample.discretization): sample_obj = sample_obj._output_sample_set + # If there is density function given determine the pointwise probability # values of each sample based on the value in the data space. Otherwise, # color the samples in numerical order. @@ -296,8 +303,8 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, if Q_nums is None: Q_nums = range(sample_obj.get_dim()) # Create the labels based on the user selected data coordinates - xlabel = r'$q_{'+str(Q_nums[0]+1)+'}$' - ylabel = r'$q_{'+str(Q_nums[1]+1)+'}$' + xlabel = r'$q_{' + str(Q_nums[0] + 1) + '}$' + ylabel = r'$q_{' + str(Q_nums[1] + 1) + '}$' savename = 'data_samples_cs.eps' # Plot 2 or 3 dimensional scatter plots of the data colored by rD. if sample_obj.get_dim() == 2: @@ -305,18 +312,18 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, if isinstance(Q_ref, np.ndarray): q_ref = Q_ref[Q_nums[:2]] scatter_2D(sample_obj, sample_nos, rD, q_ref, save, interactive, xlabel, - ylabel, savename) + ylabel, savename) elif sample_obj.get_dim() == 3: - zlabel = r'$q_{'+str(Q_nums[2]+1)+'}$' + zlabel = r'$q_{' + str(Q_nums[2] + 1) + '}$' if isinstance(Q_ref, np.ndarray): q_ref = Q_ref[Q_nums[:3]] scatter_3D(sample_obj, sample_nos, rD, q_ref, save, interactive, xlabel, - ylabel, zlabel, savename) + ylabel, zlabel, savename) elif sample_obj.get_dim() > 2 and showdim == 2: for x, y in combinations(Q_nums, 2): - xlabel = r'$q_{'+str(x+1)+'}$' - ylabel = r'$q_{'+str(y+1)+'}$' - savename = 'data_samples_q'+str(x+1)+'q'+str(y+1)+'_cs.eps' + xlabel = r'$q_{' + str(x + 1) + '}$' + ylabel = r'$q_{' + str(y + 1) + '}$' + savename = 'data_samples_q' + str(x + 1) + 'q' + str(y + 1) + '_cs.eps' q_ref = None if isinstance(Q_ref, np.ndarray): q_ref = Q_ref[[x, y]] @@ -325,32 +332,33 @@ def show_data(sample_obj, rho_D=None, Q_ref=None, sample_nos=None, sample_obj_temp.set_values(sample_obj.get_values()[:, [x, y]]) scatter_2D(sample_obj_temp, sample_nos, rD, q_ref, save, - interactive, xlabel, ylabel, savename) + interactive, xlabel, ylabel, savename) elif sample_obj.get_dim() > 3 and showdim == 3: for x, y, z in combinations(Q_nums, 3): - xlabel = r'$q_{'+str(x+1)+'}$' - ylabel = r'$q_{'+str(y+1)+'}$' - zlabel = r'$q_{'+str(z+1)+'}$' + xlabel = r'$q_{' + str(x + 1) + '}$' + ylabel = r'$q_{' + str(y + 1) + '}$' + zlabel = r'$q_{' + str(z + 1) + '}$' q_ref = None if isinstance(Q_ref, np.ndarray): q_ref = Q_ref[[x, y, z]] - savename = 'data_samples_q'+str(x+1)+'q'+str(y+1)+'q'\ - +str(z+1)+'_cs.eps' + savename = 'data_samples_q' + str(x + 1) + 'q' + str(y + 1) + 'q' \ + + str(z + 1) + '_cs.eps' sample_obj_temp = sample.sample_set(3) sample_obj_temp.set_values(sample_obj.get_values()[:, [x, y, z]]) scatter_3D(sample_obj_temp, sample_nos, rD, q_ref, save, - interactive, xlabel, ylabel, zlabel, savename) + interactive, xlabel, ylabel, zlabel, savename) + def show_data_domain_multi(sample_disc, Q_ref=None, Q_nums=None, - img_folder='figs/', ref_markers=None, - ref_colors=None, showdim=None): + img_folder='figs/', ref_markers=None, + ref_colors=None, showdim=None): r""" Plots 2-D projections of the data domain D using a triangulation based on the first two coordinates (parameters) of the generating samples where :math:`Q={q_1, q_i}` for ``i=Q_nums``, with a marker for various - :math:`Q_{ref}`. + :math:`Q_{ref}`. :param sample_disc: Object containing the samples to plot :type sample_disc: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` @@ -403,35 +411,35 @@ def show_data_domain_multi(sample_disc, Q_ref=None, Q_nums=None, # Create plots of the showdim^th QoI (q_{showdim}) with all other QoI (q_i) if isinstance(showdim, int): for i in Q_nums: - xlabel = r'$q_{'+str(showdim+1)+r'}$' - ylabel = r'$q_{'+str(i+1)+r'}$' + xlabel = r'$q_{' + str(showdim + 1) + r'}$' + ylabel = r'$q_{' + str(i + 1) + r'}$' - filenames = [img_folder+'domain_q'+str(showdim+1)+'_q'+\ - str(i+1)+'.eps', img_folder+'q'+str(showdim+1)+\ - '_q'+str(i+1)+'_domain_Q_cs.eps'] + filenames = [img_folder + 'domain_q' + str(showdim + 1) + '_q' + \ + str(i + 1) + '.eps', img_folder + 'q' + str(showdim + 1) + \ + '_q' + str(i + 1) + '_domain_Q_cs.eps'] data_obj_temp = sample.sample_set(2) data_obj_temp.set_values(data_obj.get_values()[:, [showdim, i]]) sample_disc_temp = sample.discretization(sample_obj, data_obj_temp) if Q_ref is not None: - show_data_domain_2D(sample_disc_temp, Q_ref[:,[showdim, i]], - ref_markers, ref_colors, xlabel=xlabel, - ylabel=ylabel, triangles=triangles, save=True, - interactive=False, filenames=filenames) + show_data_domain_2D(sample_disc_temp, Q_ref[:, [showdim, i]], + ref_markers, ref_colors, xlabel=xlabel, + ylabel=ylabel, triangles=triangles, save=True, + interactive=False, filenames=filenames) else: show_data_domain_2D(sample_disc_temp, None, - ref_markers, ref_colors, xlabel=xlabel, ylabel=ylabel, - triangles=triangles, save=True, interactive=False, - filenames=filenames) - # Create plots of all combinations of QoI in 2D + ref_markers, ref_colors, xlabel=xlabel, ylabel=ylabel, + triangles=triangles, save=True, interactive=False, + filenames=filenames) + # Create plots of all combinations of QoI in 2D elif showdim == 'all' or showdim == 'ALL': for x, y in combinations(Q_nums, 2): - xlabel = r'$q_{'+str(x+1)+r'}$' - ylabel = r'$q_{'+str(y+1)+r'}$' + xlabel = r'$q_{' + str(x + 1) + r'}$' + ylabel = r'$q_{' + str(y + 1) + r'}$' - filenames = [img_folder+'domain_q'+str(x+1)+'_q'+str(y+1)+'.eps', - img_folder+'q'+str(x+1)+'_q'+str(y+1)+'_domain_Q_cs.eps'] + filenames = [img_folder + 'domain_q' + str(x + 1) + '_q' + str(y + 1) + '.eps', + img_folder + 'q' + str(x + 1) + '_q' + str(y + 1) + '_domain_Q_cs.eps'] data_obj_temp = sample.sample_set(2) data_obj_temp.set_values(data_obj.get_values()[:, [x, y]]) @@ -439,18 +447,19 @@ def show_data_domain_multi(sample_disc, Q_ref=None, Q_nums=None, if Q_ref is not None: show_data_domain_2D(sample_disc_temp, Q_ref[:, [x, y]], - ref_markers, ref_colors, xlabel=xlabel, ylabel=ylabel, - triangles=triangles, save=True, interactive=False, - filenames=filenames) + ref_markers, ref_colors, xlabel=xlabel, ylabel=ylabel, + triangles=triangles, save=True, interactive=False, + filenames=filenames) else: show_data_domain_2D(sample_disc_temp, None, - ref_markers, ref_colors, xlabel=xlabel, ylabel=ylabel, - triangles=triangles, save=True, interactive=False, - filenames=filenames) + ref_markers, ref_colors, xlabel=xlabel, ylabel=ylabel, + triangles=triangles, save=True, interactive=False, + filenames=filenames) + def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None, - ref_colors=None, xlabel=r'$q_1$', ylabel=r'$q_2$', - triangles=None, save=True, interactive=False, filenames=None): + ref_colors=None, xlabel=r'$q_1$', ylabel=r'$q_2$', + triangles=None, save=True, interactive=False, filenames=None): r""" Plots 2-D a single data domain D using a triangulation based on the first two coordinates (parameters) of the generating samples where :math:`Q={q_1, @@ -494,7 +503,7 @@ def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None, # Set default file names if filenames == None: filenames = ['domain_q1_q2_cs.eps', 'q1_q2_domain_Q_cs.eps'] - + # Make sure the shape of Q_ref is correct if Q_ref is not None: Q_ref = util.fix_dimensions_data(Q_ref, 2) @@ -502,30 +511,32 @@ def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None, # Create figure plt.tricontourf(data_obj.get_values()[:, 0], data_obj.get_values()[:, 1], np.zeros((data_obj.get_values().shape[0],)), - triangles=triangles, colors='grey') + triangles=triangles, colors='grey') plt.autoscale(tight=True) plt.xlabel(xlabel) plt.ylabel(ylabel) plt.savefig(filenames[0], bbox_inches='tight', transparent=True, - pad_inches=.2) + pad_inches=.2) # Add truth markers if Q_ref is not None: for i in xrange(Q_ref.shape[0]): plt.scatter(Q_ref[i, 0], Q_ref[i, 1], s=60, c=ref_colors[i], - marker=ref_markers[i]) + marker=ref_markers[i]) if save: plt.savefig(filenames[1], bbox_inches='tight', transparent=True, - pad_inches=.2) + pad_inches=.2) if interactive: plt.show() else: plt.close() + def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True, - interactive=False): + interactive=False): r""" + Creates two-dimensional projections of scatter plots of samples. - + :param sample_obj: Object containing the samples to plot :type sample_obj: :class:`~bet.sample.sample_set` :param bool save: flag whether or not to save the figure @@ -549,23 +560,24 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True # Create list of all the parameter coordinates L_nums = range(sample_obj.get_dim()) - # Create plots of the showdim^th parameter (\lambda_{showdim}) with all the - # other parameters + # Create plots of the showdim^th parameter (\lambda_{showdim}) with all the + # other parameters if isinstance(showdim, int): for i in L_nums: - xlabel = r'$\lambda_{'+str(showdim+1)+r'}$' - ylabel = r'$\lambda_{'+str(i+1)+r'}$' - - filename = [img_folder+'domain_l'+str(showdim+1)+'_l'+\ - str(i+1)+'.eps'] + xlabel = r'$\lambda_{' + str(showdim + 1) + r'}$' + ylabel = r'$\lambda_{' + str(i + 1) + r'}$' + filenames = [img_folder + 'domain_l' + str(showdim + 1) + '_l' + \ + str(i + 1) + '.eps', img_folder + 'l' + str(showdim + 1) + \ + '_l' + str(i + 1) + '_domain_L_cs.eps'] + filename = filenames[0] plt.scatter(sample_obj.get_values()[:, 0], sample_obj.get_values()[:, 1]) if save: plt.autoscale(tight=True) plt.xlabel(xlabel) plt.ylabel(ylabel) plt.savefig(filename, bbox_inches='tight', transparent=True, - pad_inches=0) + pad_inches=0) if interactive: plt.show() else: @@ -573,27 +585,29 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True # Create plots of all of the possible pairwise combinations of parameters elif showdim == 'all' or showdim == 'ALL': for x, y in combinations(L_nums, 2): - xlabel = r'$\lambda_{'+str(x+1)+r'}$' - ylabel = r'$\lambda_{'+str(y+1)+r'}$' - - filename = [img_folder+'domain_l'+str(x+1)+'_l'+\ - str(y+1)+'.eps'] + xlabel = r'$\lambda_{' + str(x + 1) + r'}$' + ylabel = r'$\lambda_{' + str(y + 1) + r'}$' + filenames = [img_folder + 'domain_l' + str(x + 1) + '_l' + \ + str(y + 1) + '.eps', img_folder + 'l' + str(x + 1) + \ + '_l' + str(y + 1) + '_domain_L_cs.eps'] + filename = filenames[0] plt.scatter(sample_obj.get_values()[:, x], sample_obj.get_values()[:, y]) if save: plt.autoscale(tight=True) plt.xlabel(xlabel) plt.ylabel(ylabel) plt.savefig(filename, bbox_inches='tight', transparent=True, - pad_inches=0) + pad_inches=0) if interactive: plt.show() else: plt.close() + def scatter2D_multi(sample_obj, color=None, p_ref=None, img_folder='figs/', filename="scatter2Dm", label_char=r'$\lambda', - showdim=None): + showdim=None): r""" Creates two-dimensional projections of scatter plots of samples colored by ``color`` (usually an array of pointwise probability density values). A @@ -624,50 +638,50 @@ def scatter2D_multi(sample_obj, color=None, p_ref=None, img_folder='figs/', # Create a folder for these figures if it doesn't already exist if not os.path.isdir(img_folder): os.mkdir(img_folder) - # Create list of all the parameter coordinates + # Create list of all the parameter coordinates p_nums = range(sample_obj.get_dim()) - # Create plots of the showdim^th parameter (\lambda_{showdim}) with all the - # other parameters + # Create plots of the showdim^th parameter (\lambda_{showdim}) with all the + # other parameters if isinstance(showdim, int): for i in p_nums: - xlabel = label_char+r'_{'+str(showdim+1)+r'}$' - ylabel = label_char+r'_{'+str(i+1)+r'}$' + xlabel = label_char + r'_{' + str(showdim + 1) + r'}$' + ylabel = label_char + r'_{' + str(i + 1) + r'}$' - postfix = '_d'+str(showdim+1)+'_d'+str(i+1)+'.eps' - myfilename = os.path.join(img_folder, filename+postfix) + postfix = '_d' + str(showdim + 1) + '_d' + str(i + 1) + '.eps' + myfilename = os.path.join(img_folder, filename + postfix) sample_obj_temp = sample.sample_set(2) sample_obj_temp.set_values(sample_obj.get_values()[:, [showdim, i]]) if p_ref is not None: scatter_2D(sample_obj_temp, sample_nos=None, - color=color, p_ref=p_ref[[showdim, i]], save=True, - interactive=False, xlabel=xlabel, ylabel=ylabel, - filename=myfilename) + color=color, p_ref=p_ref[[showdim, i]], save=True, + interactive=False, xlabel=xlabel, ylabel=ylabel, + filename=myfilename) else: scatter_2D(sample_obj_temp, sample_nos=None, - color=color, p_ref=None, save=True, - interactive=False, xlabel=xlabel, ylabel=ylabel, - filename=myfilename) + color=color, p_ref=None, save=True, + interactive=False, xlabel=xlabel, ylabel=ylabel, + filename=myfilename) # Create plots of all of the possible pairwise combinations of parameters elif showdim == 'all' or showdim == 'ALL': for x, y in combinations(p_nums, 2): - xlabel = label_char+r'_{'+str(x+1)+r'}$' - ylabel = label_char+r'_{'+str(y+1)+r'}$' + xlabel = label_char + r'_{' + str(x + 1) + r'}$' + ylabel = label_char + r'_{' + str(y + 1) + r'}$' - postfix = '_d'+str(x+1)+'_d'+str(y+1)+'.eps' - myfilename = os.path.join(img_folder, filename+postfix) + postfix = '_d' + str(x + 1) + '_d' + str(y + 1) + '.eps' + myfilename = os.path.join(img_folder, filename + postfix) sample_obj_temp = sample.sample_set(2) sample_obj_temp.set_values(sample_obj.get_values()[:, [x, y]]) if p_ref is not None: scatter_2D(sample_obj_temp, sample_nos=None, color=color, - p_ref=p_ref[[x, y]], save=True, interactive=False, - xlabel=xlabel, ylabel=ylabel, filename=myfilename) + p_ref=p_ref[[x, y]], save=True, interactive=False, + xlabel=xlabel, ylabel=ylabel, filename=myfilename) else: scatter_2D(sample_obj_temp, sample_nos=None, color=color, - p_ref=None, save=True, interactive=False, - xlabel=xlabel, ylabel=ylabel, filename=myfilename) + p_ref=None, save=True, interactive=False, + xlabel=xlabel, ylabel=ylabel, filename=myfilename) diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index 54bd539b..f75622cc 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -76,7 +76,7 @@ if different numbers of points are used to estimate the volumes of the Voronoi cells. ''' -MC_assumption = False +MC_assumption = True # Estimate volumes of Voronoi cells associated with the parameter samples if MC_assumption is False: input_samples.estimate_volume(n_mc_points=1E5) From 1f48875b148c1db495047cf944a07a4dfb7efffa Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Sun, 22 May 2016 15:04:55 -0600 Subject: [PATCH 05/20] updating tests to use new simpleFunP functions --- test/test_calculateP/test_calculateP.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/test/test_calculateP/test_calculateP.py b/test/test_calculateP/test_calculateP.py index 376df369..00ae0d95 100644 --- a/test/test_calculateP/test_calculateP.py +++ b/test/test_calculateP/test_calculateP.py @@ -21,6 +21,7 @@ from bet.Comm import comm data_path = os.path.dirname(bet.__file__) + "/../test/test_calculateP/datafiles" + class TestEmulateIIDLebesgue(unittest.TestCase): """ Test :meth:`bet.calculateP.calculateP.emulate_iid_lebesgue`. @@ -144,7 +145,8 @@ def setUp(self): self.inputs.set_values(np.loadtxt(data_path + "/3to2_samples.txt.gz")) self.outputs.set_values(np.loadtxt(data_path + "/3to2_data.txt.gz")) Q_ref = np.array([0.422, 0.9385]) - self.output_prob = simpleFunP.uniform_hyperrectangle(self.outputs, Q_ref = Q_ref, bin_ratio=0.2, center_pts_per_edge=1) + self.output_prob = simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + self.outputs, Q_ref = Q_ref, rect_scale=0.2, center_pts_per_edge=1) self.inputs.set_domain(np.array([[0.0, 1.0], [0.0, 1.0], @@ -208,7 +210,8 @@ def setUp(self): self.inputs.set_values(np.loadtxt(data_path + "/3to2_samples.txt.gz")) self.outputs.set_values(np.loadtxt(data_path + "/3to2_data.txt.gz")[:,0]) Q_ref = np.array([0.422]) - self.output_prob = simpleFunP.uniform_hyperrectangle(self.outputs, Q_ref = Q_ref, bin_ratio=0.2, center_pts_per_edge=1) + self.output_prob = simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + self.outputs, Q_ref = Q_ref, rect_scale=0.2, center_pts_per_edge=1) self.inputs.set_domain(np.array([[0.0, 1.0], [0.0, 1.0], @@ -281,7 +284,8 @@ def setUp(self): self.outputs.set_values(np.dot(self.inputs._values, rnd.rand(10, 4))) Q_ref = np.mean(self.outputs._values, axis=0) self.inputs_emulated = calcP.emulate_iid_lebesgue(self.inputs.get_domain(), num_l_emulate=1001, globalize=True) - self.output_prob = simpleFunP.uniform_hyperrectangle(self.outputs, Q_ref = Q_ref, bin_ratio=0.2, center_pts_per_edge=1) + self.output_prob = simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + self.outputs, Q_ref = Q_ref, rect_scale=0.2, center_pts_per_edge=1) self.disc = samp.discretization(input_sample_set=self.inputs, output_sample_set=self.outputs, output_probability_set=self.output_prob, @@ -353,7 +357,8 @@ def setUp(self): self.inputs_emulated = calcP.emulate_iid_lebesgue(self.lam_domain, self.num_l_emulate, globalize = True) - self.output_prob = simpleFunP.uniform_hyperrectangle(self.outputs, Q_ref = Q_ref, bin_ratio=0.2, center_pts_per_edge=1) + self.output_prob = simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + self.outputs, Q_ref = Q_ref, rect_scale=0.2, center_pts_per_edge=1) self.disc = samp.discretization(input_sample_set=self.inputs, output_sample_set=self.outputs, output_probability_set=self.output_prob, From 0c134fda2fb4378ac759b00a59301da31effb6dc Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Sun, 22 May 2016 16:11:03 -0600 Subject: [PATCH 06/20] updating existing tests for new function names --- test/test_calculateP/test_simpleFunP.py | 393 +++++++++++++----------- 1 file changed, 216 insertions(+), 177 deletions(-) diff --git a/test/test_calculateP/test_simpleFunP.py b/test/test_calculateP/test_simpleFunP.py index 77a63346..c2c29268 100644 --- a/test/test_calculateP/test_simpleFunP.py +++ b/test/test_calculateP/test_simpleFunP.py @@ -135,15 +135,16 @@ def createData(self): self.data_domain = np.array([[0.0, 10.0], [0.0, 10.0], [0.0, 10.0]]) self.mdim = 3 -class unif_unif(prob_uniform): +class uniform_partition_uniform_distribution_rectangle_scaled(prob_uniform): """ - Set up :meth:`bet.calculateP.simpleFunP.unif_unif` on data domain. + Set up :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled` on data domain. """ def setUp(self): """ Set up problem. """ - self.data_prob = sFun.unif_unif(self.data, self.Q_ref, M=67, bin_ratio=0.1, num_d_emulate=1E3) + self.data_prob = sFun.uniform_partition_uniform_distribution_rectangle_scaled( + self.data, self.Q_ref, rect_scale=0.1, M=67, num_d_emulate=1E3) self.d_distr_samples = self.data_prob.get_values() self.rho_D_M = self.data_prob.get_probabilities() @@ -183,55 +184,59 @@ def test_domain(self): print np.sum(self.rho_D_M[np.logical_not(inside)] == 0.0) assert np.sum(self.rho_D_M[np.logical_not(inside)] == 0.0)<100 -class test_unif_unif_01D(data_01D, unif_unif): +class test_uniform_partition_uniform_distribution_rectangle_scaled_01D(data_01D, + uniform_partition_uniform_distribution_rectangle_scaled): """ - Tests :meth:`bet.calculateP.simpleFunP.unif_unif` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_unif_unif_01D, self).createData() - super(test_unif_unif_01D, self).setUp() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_01D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_01D, self).setUp() -class test_unif_unif_1D(data_1D, unif_unif): +class test_uniform_partition_uniform_distribution_rectangle_scaled_1D(data_1D, + uniform_partition_uniform_distribution_rectangle_scaled): """ - Tests :meth:`bet.calculateP.simpleFunP.unif_unif` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_unif_unif_1D, self).createData() - super(test_unif_unif_1D, self).setUp() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_1D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_1D, self).setUp() -class test_unif_unif_2D(data_2D, unif_unif): +class test_uniform_partition_uniform_distribution_rectangle_scaled_2D(data_2D, + uniform_partition_uniform_distribution_rectangle_scaled): """ - Tests :meth:`bet.calculateP.simpleFunP.unif_unif` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_unif_unif_2D, self).createData() - super(test_unif_unif_2D, self).setUp() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_2D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_2D, self).setUp() -class test_unif_unif_3D(data_3D, unif_unif): +class test_uniform_partition_uniform_distribution_rectangle_scaled_3D(data_3D, + uniform_partition_uniform_distribution_rectangle_scaled): """ - Tests :meth:`bet.calculateP.simpleFunP.unif_unif` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_unif_unif_3D, self).createData() - super(test_unif_unif_3D, self).setUp() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_3D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_scaled_3D, self).setUp() -class normal_normal(prob): +class normal_partition_normal_distribution(prob): """ - Set up :meth:`bet.calculateP.simpleFunP.normal_normal` on data domain. + Set up :meth:`bet.calculateP.simpleFunP.normal_partition_normal_distribution` on data domain. """ def setUp(self): """ @@ -241,7 +246,7 @@ def setUp(self): std = 1.0 else: std = np.ones(self.Q_ref.shape) - self.data_prob = sFun.normal_normal(None, self.Q_ref, M=67, std=std, num_d_emulate=1E3) + self.data_prob = sFun.normal_partition_normal_distribution(None, self.Q_ref, std=std, M=67, num_d_emulate=1E3) self.d_distr_samples = self.data_prob.get_values() self.rho_D_M = self.data_prob.get_probabilities() @@ -252,51 +257,51 @@ def test_M(self): """ assert len(self.rho_D_M) == 67 -class test_normal_normal_01D(data_01D, normal_normal): +class test_normal_partition_normal_distribution_01D(data_01D, normal_partition_normal_distribution): """ - Tests :meth:`bet.calculateP.simpleFunP.normal_normal` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.normal_partition_normal_distribution` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_normal_normal_01D, self).createData() - super(test_normal_normal_01D, self).setUp() + super(test_normal_partition_normal_distribution_01D, self).createData() + super(test_normal_partition_normal_distribution_01D, self).setUp() -class test_normal_normal_1D(data_1D, normal_normal): +class test_normal_partition_normal_distribution_1D(data_1D, normal_partition_normal_distribution): """ - Tests :meth:`bet.calculateP.simpleFunP.normal_normal` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.normal_partition_normal_distribution` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_normal_normal_1D, self).createData() - super(test_normal_normal_1D, self).setUp() + super(test_normal_partition_normal_distribution_1D, self).createData() + super(test_normal_partition_normal_distribution_1D, self).setUp() -class test_normal_normal_2D(data_2D, normal_normal): +class test_normal_partition_normal_distribution_2D(data_2D, normal_partition_normal_distribution): """ - Tests :meth:`bet.calculateP.simpleFunP.normal_normal` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.normal_partition_normal_distribution` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_normal_normal_2D, self).createData() - super(test_normal_normal_2D, self).setUp() + super(test_normal_partition_normal_distribution_2D, self).createData() + super(test_normal_partition_normal_distribution_2D, self).setUp() -class test_normal_normal_3D(data_3D, normal_normal): +class test_normal_partition_normal_distribution_3D(data_3D, normal_partition_normal_distribution): """ - Tests :meth:`bet.calculateP.simpleFunP.normal_normal` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.normal_partition_normal_distribution` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_normal_normal_3D, self).createData() - super(test_normal_normal_3D, self).setUp() + super(test_normal_partition_normal_distribution_3D, self).createData() + super(test_normal_partition_normal_distribution_3D, self).setUp() class uniform_hyperrectangle_base(prob_uniform): @@ -338,9 +343,9 @@ def setUp(self): """ self.center_pts_per_edge = 2*np.ones((self.mdim,), dtype=np.int) -class uniform_hyperrectangle_user_int(uniform_hyperrectangle_int): +class regular_partition_uniform_distribution_rectangle_domain_int(uniform_hyperrectangle_int): """ - Set up :met:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user` with an + Set up :met:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` with an int type of value fo r``center_pts_per_edge`` """ @@ -348,7 +353,7 @@ def setUp(self): """ Set up problem. """ - super(uniform_hyperrectangle_user_int, self).setUp() + super(regular_partition_uniform_distribution_rectangle_domain_int, self).setUp() if type(self.Q_ref) != np.array: Q_ref = np.array([self.Q_ref]) else: @@ -364,13 +369,14 @@ def setUp(self): self.rect_domain[:, 0] = Q_ref - .5*r_width self.rect_domain[:, 1] = Q_ref + .5*r_width - self.data_prob = sFun.uniform_hyperrectangle_user(self.data, self.rect_domain.transpose(), self.center_pts_per_edge) + self.data_prob = sFun.regular_partition_uniform_distribution_rectangle_domain( + self.data, self.rect_domain.transpose(), self.center_pts_per_edge) self.rho_D_M = self.data_prob._probabilities self.d_distr_samples = self.data_prob._values -class uniform_hyperrectangle_user_list(uniform_hyperrectangle_list): +class regular_partition_uniform_distribution_rectangle_domain_list(uniform_hyperrectangle_list): """ - Set up :met:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user` with an + Set up :met:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` with an int type of value fo r``center_pts_per_edge`` """ @@ -378,7 +384,7 @@ def setUp(self): """ Set up problem. """ - super(uniform_hyperrectangle_user_list, self).setUp() + super(regular_partition_uniform_distribution_rectangle_domain_list, self).setUp() if type(self.Q_ref) != np.array: Q_ref = np.array([self.Q_ref]) else: @@ -394,108 +400,117 @@ def setUp(self): self.rect_domain[:, 0] = Q_ref - .5*r_width self.rect_domain[:, 1] = Q_ref + .5*r_width - self.data_prob = sFun.uniform_hyperrectangle_user(self.data, self.rect_domain.transpose(), self.center_pts_per_edge) + self.data_prob = sFun.regular_partition_uniform_distribution_rectangle_domain( + self.data, self.rect_domain.transpose(), self.center_pts_per_edge) self.rho_D_M = self.data_prob._probabilities self.d_distr_samples = self.data_prob._values -class test_uniform_hyperrectangle_user_int_01D(data_01D, uniform_hyperrectangle_user_int): +class test_regular_partition_uniform_distribution_rectangle_domain_int_01D(data_01D, + regular_partition_uniform_distribution_rectangle_domain_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_int` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_int_01D, self).createData() - super(test_uniform_hyperrectangle_user_int_01D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_01D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_01D, self).setUp() -class test_uniform_hyperrectangle_user_int_1D(data_1D, uniform_hyperrectangle_user_int): +class test_regular_partition_uniform_distribution_rectangle_domain_int_1D(data_1D, + regular_partition_uniform_distribution_rectangle_domain_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_int` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_int_1D, self).createData() - super(test_uniform_hyperrectangle_user_int_1D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_1D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_1D, self).setUp() -class test_uniform_hyperrectangle_user_int_2D(data_2D, uniform_hyperrectangle_user_int): +class test_regular_partition_uniform_distribution_rectangle_domain_int_2D(data_2D, + regular_partition_uniform_distribution_rectangle_domain_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_int` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_int_2D, self).createData() - super(test_uniform_hyperrectangle_user_int_2D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_2D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_2D, self).setUp() -class test_uniform_hyperrectangle_user_int_3D(data_3D, uniform_hyperrectangle_user_int): +class test_regular_partition_uniform_distribution_rectangle_domain_int_3D(data_3D, + regular_partition_uniform_distribution_rectangle_domain_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_int` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_int_3D, self).createData() - super(test_uniform_hyperrectangle_user_int_3D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_3D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_int_3D, self).setUp() -class test_uniform_hyperrectangle_user_list_01D(data_01D, uniform_hyperrectangle_user_list): +class test_regular_partition_uniform_distribution_rectangle_domain_list_01D(data_01D, + regular_partition_uniform_distribution_rectangle_domain_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_list` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_list_01D, self).createData() - super(test_uniform_hyperrectangle_user_list_01D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_01D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_01D, self).setUp() -class test_uniform_hyperrectangle_user_list_1D(data_1D, uniform_hyperrectangle_user_list): +class test_regular_partition_uniform_distribution_rectangle_domain_list_1D(data_1D, + regular_partition_uniform_distribution_rectangle_domain_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_list` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_list_1D, self).createData() - super(test_uniform_hyperrectangle_user_list_1D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_1D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_1D, self).setUp() -class test_uniform_hyperrectangle_user_list_2D(data_2D, uniform_hyperrectangle_user_list): +class test_regular_partition_uniform_distribution_rectangle_domain_list_2D(data_2D, + regular_partition_uniform_distribution_rectangle_domain_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_list` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_list_2D, self).createData() - super(test_uniform_hyperrectangle_user_list_2D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_2D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_2D, self).setUp() -class test_uniform_hyperrectangle_user_list_3D(data_3D, uniform_hyperrectangle_user_list): +class test_regular_partition_uniform_distribution_rectangle_domain_list_3D(data_3D, + regular_partition_uniform_distribution_rectangle_domain_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_user_list` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_domain` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_user_list_3D, self).createData() - super(test_uniform_hyperrectangle_user_list_3D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_3D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_domain_list_3D, self).setUp() -class uniform_hyperrectangle_size_int(uniform_hyperrectangle_int): +class regular_partition_uniform_distribution_rectangle_size_int(uniform_hyperrectangle_int): """ - Set up :met:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size` with an + Set up :met:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` with an int type of value fo r``center_pts_per_edge`` """ @@ -503,7 +518,7 @@ def setUp(self): """ Set up problem. """ - super(uniform_hyperrectangle_size_int, self).setUp() + super(regular_partition_uniform_distribution_rectangle_size_int, self).setUp() if type(self.Q_ref) != np.array: Q_ref = np.array([self.Q_ref]) else: @@ -520,13 +535,14 @@ def setUp(self): self.rect_domain[:, 0] = Q_ref - .5*r_width self.rect_domain[:, 1] = Q_ref + .5*r_width - self.data_prob = sFun.uniform_hyperrectangle_binsize(self.data, self.Q_ref,binsize, self.center_pts_per_edge) + self.data_prob = sFun.regular_partition_uniform_distribution_rectangle_size( + self.data, self.Q_ref, binsize, self.center_pts_per_edge) self.rho_D_M = self.data_prob._probabilities self.d_distr_samples = self.data_prob._values -class uniform_hyperrectangle_size_list(uniform_hyperrectangle_list): +class regular_partition_uniform_distribution_rectangle_size_list(uniform_hyperrectangle_list): """ - Set up :met:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size` with an + Set up :met:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` with an int type of value fo r``center_pts_per_edge`` """ @@ -534,7 +550,7 @@ def setUp(self): """ Set up problem. """ - super(uniform_hyperrectangle_size_list, self).setUp() + super(regular_partition_uniform_distribution_rectangle_size_list, self).setUp() if type(self.Q_ref) != np.array: Q_ref = np.array([self.Q_ref]) else: @@ -551,107 +567,116 @@ def setUp(self): self.rect_domain[:, 0] = Q_ref - .5*r_width self.rect_domain[:, 1] = Q_ref + .5*r_width - self.data_prob = sFun.uniform_hyperrectangle_binsize(self.data, self.Q_ref,binsize, self.center_pts_per_edge) + self.data_prob = sFun.regular_partition_uniform_distribution_rectangle_size( + self.data, self.Q_ref, binsize, self.center_pts_per_edge) self.rho_D_M = self.data_prob._probabilities self.d_distr_samples = self.data_prob._values -class test_uniform_hyperrectangle_size_int_01D(data_01D, uniform_hyperrectangle_size_int): +class test_regular_partition_uniform_distribution_rectangle_size_int_01D(data_01D, + regular_partition_uniform_distribution_rectangle_size_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_int` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_int_01D, self).createData() - super(test_uniform_hyperrectangle_size_int_01D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_int_01D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_int_01D, self).setUp() -class test_uniform_hyperrectangle_size_int_1D(data_1D, uniform_hyperrectangle_size_int): +class test_regular_partition_uniform_distribution_rectangle_size_int_1D(data_1D, + regular_partition_uniform_distribution_rectangle_size_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_int` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_int_1D, self).createData() - super(test_uniform_hyperrectangle_size_int_1D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_int_1D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_int_1D, self).setUp() -class test_uniform_hyperrectangle_size_int_2D(data_2D, uniform_hyperrectangle_size_int): +class test_regular_partition_uniform_distribution_rectangle_size_int_2D(data_2D, + regular_partition_uniform_distribution_rectangle_size_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_int` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_int_2D, self).createData() - super(test_uniform_hyperrectangle_size_int_2D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_int_2D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_int_2D, self).setUp() -class test_uniform_hyperrectangle_size_int_3D(data_3D, uniform_hyperrectangle_size_int): +class test_regular_partition_uniform_distribution_rectangle_size_int_3D(data_3D, + regular_partition_uniform_distribution_rectangle_size_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_int` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_int_3D, self).createData() - super(test_uniform_hyperrectangle_size_int_3D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_int_3D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_int_3D, self).setUp() -class test_uniform_hyperrectangle_size_list_01D(data_01D, uniform_hyperrectangle_size_list): +class test_regular_partition_uniform_distribution_rectangle_size_list_01D(data_01D, + regular_partition_uniform_distribution_rectangle_size_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_list` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_list_01D, self).createData() - super(test_uniform_hyperrectangle_size_list_01D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_list_01D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_list_01D, self).setUp() -class test_uniform_hyperrectangle_size_list_1D(data_1D, uniform_hyperrectangle_size_list): +class test_regular_partition_uniform_distribution_rectangle_size_list_1D(data_1D, + regular_partition_uniform_distribution_rectangle_size_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_list` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_list_1D, self).createData() - super(test_uniform_hyperrectangle_size_list_1D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_list_1D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_list_1D, self).setUp() -class test_uniform_hyperrectangle_size_list_2D(data_2D, uniform_hyperrectangle_size_list): +class test_regular_partition_uniform_distribution_rectangle_size_list_2D(data_2D, + regular_partition_uniform_distribution_rectangle_size_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_list` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_list_2D, self).createData() - super(test_uniform_hyperrectangle_size_list_2D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_list_2D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_list_2D, self).setUp() -class test_uniform_hyperrectangle_size_list_3D(data_3D, uniform_hyperrectangle_size_list): +class test_regular_partition_uniform_distribution_rectangle_size_list_3D(data_3D, + regular_partition_uniform_distribution_rectangle_size_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_size_list` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_size` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_size_list_3D, self).createData() - super(test_uniform_hyperrectangle_size_list_3D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_size_list_3D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_size_list_3D, self).setUp() -class uniform_hyperrectangle_ratio_int(uniform_hyperrectangle_int): +class regular_partition_uniform_distribution_rectangle_scaled_int(uniform_hyperrectangle_int): """ - Set up :met:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio` with an + Set up :met:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` with an int type of value fo r``center_pts_per_edge`` """ @@ -659,7 +684,7 @@ def setUp(self): """ Set up problem. """ - super(uniform_hyperrectangle_ratio_int, self).setUp() + super(regular_partition_uniform_distribution_rectangle_scaled_int, self).setUp() if type(self.Q_ref) != np.array: Q_ref = np.array([self.Q_ref]) else: @@ -676,13 +701,14 @@ def setUp(self): self.rect_domain[:, 0] = Q_ref - .5*r_width self.rect_domain[:, 1] = Q_ref + .5*r_width - self.data_prob = sFun.uniform_hyperrectangle(self.data, self.Q_ref, binratio, self.center_pts_per_edge) + self.data_prob = sFun.regular_partition_uniform_distribution_rectangle_scaled( + self.data, self.Q_ref, binratio, self.center_pts_per_edge) self.rho_D_M = self.data_prob._probabilities self.d_distr_samples = self.data_prob._values -class uniform_hyperrectangle_ratio_list(uniform_hyperrectangle_list): +class regular_partition_uniform_distribution_rectangle_scaled_list(uniform_hyperrectangle_list): """ - Set up :met:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio` with an + Set up :met:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` with an int type of value fo r``center_pts_per_edge`` """ @@ -690,7 +716,7 @@ def setUp(self): """ Set up problem. """ - super(uniform_hyperrectangle_ratio_list, self).setUp() + super(regular_partition_uniform_distribution_rectangle_scaled_list, self).setUp() if type(self.Q_ref) != np.array: Q_ref = np.array([self.Q_ref]) else: @@ -707,113 +733,122 @@ def setUp(self): self.rect_domain[:, 0] = Q_ref - .5*r_width self.rect_domain[:, 1] = Q_ref + .5*r_width - self.data_prob = sFun.uniform_hyperrectangle(self.data, self.Q_ref, binratio, self.center_pts_per_edge) + self.data_prob = sFun.regular_partition_uniform_distribution_rectangle_scaled( + self.data, self.Q_ref, binratio, self.center_pts_per_edge) self.rho_D_M = self.data_prob._probabilities self.d_distr_samples = self.data_prob._values -class test_uniform_hyperrectangle_ratio_int_01D(data_01D, uniform_hyperrectangle_ratio_int): +class test_regular_partition_uniform_distribution_rectangle_scaled_int_01D(data_01D, + regular_partition_uniform_distribution_rectangle_scaled_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_int` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_int_01D, self).createData() - super(test_uniform_hyperrectangle_ratio_int_01D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_01D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_01D, self).setUp() -class test_uniform_hyperrectangle_ratio_int_1D(data_1D, uniform_hyperrectangle_ratio_int): +class test_regular_partition_uniform_distribution_rectangle_scaled_int_1D(data_1D, + regular_partition_uniform_distribution_rectangle_scaled_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_int` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_int_1D, self).createData() - super(test_uniform_hyperrectangle_ratio_int_1D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_1D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_1D, self).setUp() -class test_uniform_hyperrectangle_ratio_int_2D(data_2D, uniform_hyperrectangle_ratio_int): +class test_regular_partition_uniform_distribution_rectangle_scaled_int_2D(data_2D, + regular_partition_uniform_distribution_rectangle_scaled_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_int` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_int_2D, self).createData() - super(test_uniform_hyperrectangle_ratio_int_2D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_2D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_2D, self).setUp() -class test_uniform_hyperrectangle_ratio_int_3D(data_3D, uniform_hyperrectangle_ratio_int): +class test_regular_partition_uniform_distribution_rectangle_scaled_int_3D(data_3D, + regular_partition_uniform_distribution_rectangle_scaled_int): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_int` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_int_3D, self).createData() - super(test_uniform_hyperrectangle_ratio_int_3D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_3D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_int_3D, self).setUp() -class test_uniform_hyperrectangle_ratio_list_01D(data_01D, uniform_hyperrectangle_ratio_list): +class test_regular_partition_uniform_distribution_rectangle_scaled_list_01D(data_01D, + regular_partition_uniform_distribution_rectangle_scaled_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_list` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled_list` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_list_01D, self).createData() - super(test_uniform_hyperrectangle_ratio_list_01D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_01D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_01D, self).setUp() -class test_uniform_hyperrectangle_ratio_list_1D(data_1D, uniform_hyperrectangle_ratio_list): +class test_regular_partition_uniform_distribution_rectangle_scaled_list_1D(data_1D, + regular_partition_uniform_distribution_rectangle_scaled_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_list` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_list_1D, self).createData() - super(test_uniform_hyperrectangle_ratio_list_1D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_1D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_1D, self).setUp() -class test_uniform_hyperrectangle_ratio_list_2D(data_2D, uniform_hyperrectangle_ratio_list): +class test_regular_partition_uniform_distribution_rectangle_scaled_list_2D(data_2D, + regular_partition_uniform_distribution_rectangle_scaled_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_list` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_list_2D, self).createData() - super(test_uniform_hyperrectangle_ratio_list_2D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_2D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_2D, self).setUp() -class test_uniform_hyperrectangle_ratio_list_3D(data_3D, uniform_hyperrectangle_ratio_list): +class test_regular_partition_uniform_distribution_rectangle_scaled_list_3D(data_3D, + regular_partition_uniform_distribution_rectangle_scaled_list): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_hyperrectangle_ratio_list` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.regular_partition_uniform_distribution_rectangle_scaled` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_hyperrectangle_ratio_list_3D, self).createData() - super(test_uniform_hyperrectangle_ratio_list_3D, self).setUp() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_3D, self).createData() + super(test_regular_partition_uniform_distribution_rectangle_scaled_list_3D, self).setUp() -class uniform_data(prob_uniform): +class uniform_partition_uniform_distribution_data_samples(prob_uniform): """ - Set up :meth:`bet.calculateP.simpleFunP.uniform_data` on data domain. + Set up :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_data_samples` on data domain. """ def setUp(self): """ Set up problem. """ - self.data_prob = sFun.uniform_data(self.data) + self.data_prob = sFun.uniform_partition_uniform_distribution_data_samples(self.data) self.d_distr_samples = self.data_prob.get_values() self.rho_D_M = self.data_prob.get_probabilities() self.data = self.data._values @@ -831,50 +866,54 @@ def test_M(self): """ assert len(self.rho_D_M) == self.data.shape[0] -class test_uniform_data_01D(data_01D, uniform_data): +class test_uniform_partition_uniform_distribution_data_samples_01D(data_01D, + uniform_partition_uniform_distribution_data_samples): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_data` on 01D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_data_samples` on 01D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_data_01D, self).createData() - super(test_uniform_data_01D, self).setUp() + super(test_uniform_partition_uniform_distribution_data_samples_01D, self).createData() + super(test_uniform_partition_uniform_distribution_data_samples_01D, self).setUp() -class test_uniform_data_1D(data_1D, uniform_data): +class test_uniform_partition_uniform_distribution_data_samples_1D(data_1D, + uniform_partition_uniform_distribution_data_samples): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_data` on 1D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_data_samples` on 1D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_data_1D, self).createData() - super(test_uniform_data_1D, self).setUp() + super(test_uniform_partition_uniform_distribution_data_samples_1D, self).createData() + super(test_uniform_partition_uniform_distribution_data_samples_1D, self).setUp() -class test_uniform_data_2D(data_2D, uniform_data): +class test_uniform_partition_uniform_distribution_data_samples_2D(data_2D, + uniform_partition_uniform_distribution_data_samples): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_data` on 2D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_data_samples` on 2D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_data_2D, self).createData() - super(test_uniform_data_2D, self).setUp() + super(test_uniform_partition_uniform_distribution_data_samples_2D, self).createData() + super(test_uniform_partition_uniform_distribution_data_samples_2D, self).setUp() -class test_uniform_data_3D(data_3D, uniform_data): +class test_uniform_partition_uniform_distribution_data_samples_3D(data_3D, + uniform_partition_uniform_distribution_data_samples): """ - Tests :meth:`bet.calculateP.simpleFunP.uniform_data` on 3D data domain. + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_data_samples` on 3D data domain. """ def setUp(self): """ Set up problem. """ - super(test_uniform_data_3D, self).createData() - super(test_uniform_data_3D, self).setUp() + super(test_uniform_partition_uniform_distribution_data_samples_3D, self).createData() + super(test_uniform_partition_uniform_distribution_data_samples_3D, self).setUp() From b5c6e2002c58781a82b9ff1b09e122372cce42c1 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Sun, 22 May 2016 17:38:42 -0600 Subject: [PATCH 07/20] Updates to examples --- .../linearMap/linearMapUniformSampling.py | 3 +- examples/nonlinearMap/myModel.py | 70 +++++ .../nonlinearMapUniformSampling.py | 243 +++++++----------- test/test_calculateP/test_simpleFunP.py | 218 ++++++++++++++++ 4 files changed, 378 insertions(+), 156 deletions(-) create mode 100644 examples/nonlinearMap/myModel.py diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index f75622cc..b6ff6d4f 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -8,7 +8,8 @@ or just a QoI. We refer to the range of the QoI map as the data space. The 3-D input space is discretized with i.i.d. uniform -random samples. We refer to the input space as the +random samples or a regular grid of samples. +We refer to the input space as the parameter space, and use parameter to refer to a particular point (e.g., a particular random sample) in this space. A reference parameter is used to define a reference QoI datum diff --git a/examples/nonlinearMap/myModel.py b/examples/nonlinearMap/myModel.py new file mode 100644 index 00000000..e7fda009 --- /dev/null +++ b/examples/nonlinearMap/myModel.py @@ -0,0 +1,70 @@ +# Copyright (C) 2016 The BET Development Team + +# -*- coding: utf-8 -*- +import numpy as np +import math as m + +''' +Suggested changes for user: + +Try setting QoI_num = 2. + +Play around with the x1, y1, and/or, x2, y2 values to try and +"optimize" the QoI to give the highest probability region +on the reference parameter above. + +Hint: Try using QoI_num = 1 and systematically varying the +x1 and y1 values to find QoI with contour structures (as inferred +through the 2D marginal plots) that are nearly orthogonal. + +Some interesting pairs of QoI to compare are: +(x1,y1)=(0.5,0.5) and (x2,y2)=(0.25,0.25) +(x1,y1)=(0.5,0.5) and (x2,y2)=(0.15,0.15) +(x1,y1)=(0.5,0.5) and (x2,y2)=(0.25,0.15) +''' +# Choose the number of QoI +QoI_num = 2 + +# Specify the spatial points to take measurements of solution defining the QoI +if QoI_num == 1: + x1 = 0.5 + y1 = 0.5 + x = np.array([x1]) + y = np.array([y1]) +else: + x1 = 0.5 + y1 = 0.15 + x2 = 0.15 + y2 = 0.25 + x = np.array([x1, x2]) + y = np.array([y1, y2]) + +class QoI_component(object): + def __init__(self, x, y): + self.x = x + self.y = y + def eval(self, parameter_samples): + if parameter_samples.shape == (2,): + lam1 = parameter_samples[0] + lam2 = parameter_samples[1] + else: + lam1 = parameter_samples[:,0] + lam2 = parameter_samples[:,1] + z = np.sin(m.pi * self.x * lam1) * np.sin(m.pi * self.y * lam2) + return z + +# Specify the QoI maps +if QoI_num == 1: + def QoI_map(parameter_samples): + Q1 = QoI_component(x[0], y[0]) + return np.array([Q1.eval(parameter_samples)]).transpose() +else: + def QoI_map(parameter_samples): + Q1 = QoI_component(x[0], y[0]) + Q2 = QoI_component(x[1], y[1]) + return np.array([Q1.eval(parameter_samples), Q2.eval(parameter_samples)]).transpose() + +# Define a model that is the QoI map +def my_model(parameter_samples): + QoI_samples = QoI_map(parameter_samples) + return QoI_samples \ No newline at end of file diff --git a/examples/nonlinearMap/nonlinearMapUniformSampling.py b/examples/nonlinearMap/nonlinearMapUniformSampling.py index 00ac9909..a487f2b2 100644 --- a/examples/nonlinearMap/nonlinearMapUniformSampling.py +++ b/examples/nonlinearMap/nonlinearMapUniformSampling.py @@ -19,216 +19,149 @@ and :math:`\Omega=[0,1]\times[0,1]`. Probabilities -in the paramter space are calculated using emulated points. +in the parameter space are calculated using emulated points. 1D and 2D marginals are calculated, smoothed, and plotted. """ + import numpy as np -import math as m import bet.calculateP as calculateP import bet.postProcess as postProcess import bet.calculateP.simpleFunP as simpleFunP import bet.calculateP.calculateP as calculateP import bet.postProcess.plotP as plotP +import bet.postProcess.plotDomains as plotD +import bet.sample as samp +import bet.sampling.basicSampling as bsam +from myModel import my_model + + +# Initialize 3-dimensional input parameter sample set object +input_samples = samp.sample_set(2) -# parameter domain -lam_domain= np.array([[3.0, 6.0], - [1.0, 5.0]]) +# Set parameter domain +input_samples.set_domain(np.array([[3.0, 6.0], + [1.0, 5.0]])) -# reference parameters -ref_lam = [5.5, 4.5] +# Define the sampler that will be used to create the discretization +# object, which is the fundamental object used by BET to compute +# solutions to the stochastic inverse problem +sampler = bsam.sampler(my_model) ''' Suggested changes for user: - -Try setting n0 and n1 both to 10 and compare the results. - -Also, we can do uniform random sampling by setting - random_sample = True - -If random_sample = True, consider defining - - n_samples = 1E3 - -Then also try n_samples = 1E4. What happens when n_samples = 1E2? -''' -random_sample = False +Try with and without random sampling. -if random_sample == False: - n0 = 50 # number of samples in lam0 direction - n1 = 50 # number of samples in lam1 direction - n_samples = n0*n1 -else: - n_samples = 1E3 - -#set up samples -if random_sample == False: - vec0 = list(np.linspace(lam_domain[0][0], lam_domain[0][1], n0)) - vec1 = list(np.linspace(lam_domain[1][0], lam_domain[1][1], n1)) - vecv0, vecv1 = np.meshgrid(vec0, vec1, indexing='ij') - samples = np.vstack((vecv0.flat[:], vecv1.flat[:])).transpose() -else: - samples = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, - num_l_emulate = n_samples) +If using random sampling, try num_samples = 1E3 and 1E4. +What happens when num_samples = 1E2? +Try using 'lhs' instead of 'random' in the random_sample_set. -# QoI function -def QoI(x,y,lam1,lam2): - z = np.sin(m.pi*x*lam1)*np.sin(m.pi*y*lam2) - return z #np.vstack(z.flat[:]).transpose() +If using regular sampling, try different numbers of samples +per dimension. +''' +# Generate samples on the parameter space +randomSampling = False +if randomSampling is True: + sampler.random_sample_set('random', input_samples, num_samples=1E4) +else: + sampler.regular_sample_set(input_samples, num_samples_per_dim=[50, 50]) ''' Suggested changes for user: -Try setting QoI_num = 2. - -Play around with the x1, y1, and/or, x2, y2 values to try and -"optimize" the QoI to give the highest probability region -on the reference parameter above. +A standard Monte Carlo (MC) assumption is that every Voronoi cell +has the same volume. If a regular grid of samples was used, then +the standard MC assumption is true. -Hint: Try using QoI_num = 1 and systematically varying the -x1 and y1 values to find QoI with contour structures (as inferred -through the 2D marginal plots) that are nearly orthogonal. - -Some interesting pairs of QoI to compare are: -(x1,y1)=(0.5,0.5) and (x2,y2)=(0.25,0.25) -(x1,y1)=(0.5,0.5) and (x2,y2)=(0.15,0.15) -(x1,y1)=(0.5,0.5) and (x2,y2)=(0.25,0.15) +See what happens if the MC assumption is not assumed to be true, and +if different numbers of points are used to estimate the volumes of +the Voronoi cells. ''' -# Choose the QoI and define Q_ref -QoI_num = 2 - -if QoI_num == 1: - x1 = 0.5 - y1 = 0.5 - x = np.array([x1]) - y = np.array([y1]) - Q_ref = np.array([QoI(x[0],y[0],ref_lam[0],ref_lam[1])]) -else: - x1 = 0.5 - y1 = 0.15 - x2 = 0.15 - y2 = 0.25 - x = np.array([x1,x2]) - y = np.array([y1,y2]) - Q_ref = np.array([QoI(x[0],y[0],ref_lam[0],ref_lam[1]), - QoI(x[1],y[1],ref_lam[0],ref_lam[1])]) - -if QoI_num == 1: - def QoI_map(x,y,lam1,lam2): - Q1 = QoI(x[0],y[0],lam1,lam2) - z = np.array([Q1]).transpose() - return z +MC_assumption = True +# Estimate volumes of Voronoi cells associated with the parameter samples +if MC_assumption is False: + input_samples.estimate_volume(n_mc_points=1E5) else: - def QoI_map(x,y,lam1,lam2): - Q1 = QoI(x[0],y[0],lam1,lam2) - Q2 = QoI(x[1],y[1],lam1,lam2) - z = np.array([Q1,Q2]).transpose() - return z - -# calc data -data = QoI_map(x,y,samples[:,0],samples[:,1]) + input_samples.estimate_volume_mc() + +# Create the discretization object using the input samples +my_discretization = sampler.compute_QoI_and_create_discretization(input_samples, + savefile = 'NonlinearExample.txt.gz') ''' Suggested changes for user: - -Try different ways of discretizing the probability measure on D defined -as a uniform probability measure on a rectangle (if QoI_num = 2) or on -an interval (if QoI_num = 1). - -unif_unif creates a uniform measure on a hyperbox with dimensions -relative to the size of the circumscribed hyperbox of the set D using -the bin_ratio. A total of M samples are drawn within a slightly larger -scaled hyperbox to discretize this measure defining M total generalized -contour events in Lambda. The reason a slightly larger scaled hyperbox -is used to draw the samples to discretize D is because otherwise every -generalized contour event will have non-zero probability which obviously -defeats the purpose of "localizing" the probability within a subset of D. - -uniform_hyperrectangle uses the same measure defined in the same way as -unif_unif, but the difference is in the discretization which is on a -regular grid defined by center_pts_per_edge. If center_pts_per_edge = 1, -then the contour event corresponding to the entire support of rho_D is -approximated as a single event. This is done by carefully placing a -regular 3x3 grid (for the D=2 case) of points in D with the center -point of the grid in the center of the support of the measure and the -other points placed outside of the rectangle defining the support to -define a total of 9 contour events with 8 of them with zero probability. + +Try different reference parameters. ''' -deterministic_discretize_D = True +# Define the reference parameter +param_ref = np.array([5.5, 4.5]) +#param_ref = np.array([4.5, 3.0]) +#param_ref = np.array([3.5, 1.5]) -if deterministic_discretize_D == True: - (d_distr_prob, d_distr_samples, d_Tree) = simpleFunP.uniform_hyperrectangle(data=data, - Q_ref=Q_ref, bin_ratio=0.2, center_pts_per_edge = 1) -else: - (d_distr_prob, d_distr_samples, d_Tree) = simpleFunP.unif_unif(data=data, - Q_ref=Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E5) +# Compute the reference QoI +Q_ref = my_model(param_ref) +# Create some plots of input and output discretizations +plotD.scatter_2D(input_samples, p_ref = param_ref, filename = 'nonlinearMapParameterSamples.eps') +if Q_ref.size == 2: + plotD.show_data(my_discretization, Q_ref = Q_ref) -# create emulated points ''' Suggested changes for user: -If using a regular grid of sampling (if random_sample = False), we set - - lambda_emulate = samples - -Otherwise, play around with num_l_emulate. A value of 1E2 will probably -give poor results while results become fairly consistent with values -that are approximately 10x the number of samples. - -Note that you can always use - - lambda_emulate = samples - -and this simply will imply that a standard Monte Carlo assumption is -being used, which in a measure-theoretic context implies that each -Voronoi cell is assumed to have the same measure. This type of -approximation is more reasonable for large n_samples due to the slow -convergence rate of Monte Carlo (it converges like 1/sqrt(n_samples)). +Try different ways of discretizing the probability measure on D defined +as a uniform probability measure on a rectangle or interval depending +on choice of QoI_num in myModel.py. ''' -if random_sample == False: - lambda_emulate = samples +randomDataDiscretization = False +if randomDataDiscretization is False: + simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.25, + center_pts_per_edge = 3) else: - lambda_emulate = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, num_l_emulate = 1E5) - + simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref, rect_scale=0.25, + M=50, num_d_emulate=1E5) # calculate probablities -(P, lambda_emulate, io_ptr, emulate_ptr) = calculateP.prob_emulated(samples=samples, - data=data, rho_D_M = d_distr_prob, d_distr_samples = d_distr_samples, - lambda_emulate=lambda_emulate, d_Tree=d_Tree) -# calculate 2d marginal probs +calculateP.prob(my_discretization) + +######################################## +# Post-process the results +######################################## ''' Suggested changes for user: - + At this point, the only thing that should change in the plotP.* inputs should be either the nbins values or sigma (which influences the kernel density estimation with smaller values implying a density estimate that looks more like a histogram and larger values smoothing out the values more). - + There are ways to determine "optimal" smoothing parameters (e.g., see CV, GCV, and other similar methods), but we have not incorporated these into the code -as lower-dimensional marginal plots have limited value in understanding the -structure of a high dimensional non-parametric probability measure. +as lower-dimensional marginal plots generally have limited value in understanding +the structure of a high dimensional non-parametric probability measure. ''' -(bins, marginals2D) = plotP.calculate_2D_marginal_probs(P_samples = P, samples = lambda_emulate, lam_domain = lam_domain, nbins = [20, 20]) +# calculate 2d marginal probs +(bins, marginals2D) = plotP.calculate_2D_marginal_probs(input_samples, + nbins = [20, 20]) # smooth 2d marginals probs (optional) -marginals2D = plotP.smooth_marginals_2D(marginals2D,bins, sigma=0.5) +marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.5) # plot 2d marginals probs -plotP.plot_2D_marginal_probs(marginals2D, bins, lam_domain, filename = "nonlinearMap", - plot_surface=False) +plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "nomlinearMap", + file_extension = ".eps", plot_surface=False) - # calculate 1d marginal probs -(bins, marginals1D) = plotP.calculate_1D_marginal_probs(P_samples = P, samples = lambda_emulate, lam_domain = lam_domain, nbins = [20, 20]) +(bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, + nbins = [20, 20]) # smooth 1d marginal probs (optional) marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.5) -# plot 1d marginal probs -plotP.plot_1D_marginal_probs(marginals1D, bins, lam_domain, filename = "nonlinearMap") - - +# plot 2d marginal probs +plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "nonlinearMap", + file_extension = ".eps") diff --git a/test/test_calculateP/test_simpleFunP.py b/test/test_calculateP/test_simpleFunP.py index c2c29268..e49645c9 100644 --- a/test/test_calculateP/test_simpleFunP.py +++ b/test/test_calculateP/test_simpleFunP.py @@ -917,3 +917,221 @@ def setUp(self): super(test_uniform_partition_uniform_distribution_data_samples_3D, self).setUp() +class uniform_partition_uniform_distribution_rectangle_size(prob_uniform): + """ + Set up :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_size` on data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + self.data_prob = sFun.uniform_partition_uniform_distribution_rectangle_size( + self.data, self.Q_ref, rect_size=1.0, M=67, num_d_emulate=1E3) + self.d_distr_samples = self.data_prob.get_values() + self.rho_D_M = self.data_prob.get_probabilities() + + if type(self.Q_ref) != np.array: + self.Q_ref = np.array([self.Q_ref]) + if len(self.data_domain.shape) == 1: + self.data_domain = np.expand_dims(self.data_domain, axis=0) + + self.rect_domain = np.zeros((self.data_domain.shape[0], 2)) + + binsize = 1.0 + r_width = binsize * np.ones(self.data_domain[:, 1].shape) + + self.rect_domain[:, 0] = self.Q_ref - .5 * r_width + self.rect_domain[:, 1] = self.Q_ref + .5 * r_width + + def test_M(self): + """ + Test that the right number of d_distr_samples are used to create + rho_D_M. + """ + assert len(self.rho_D_M) == 67 + + def test_domain(self): + """ + Test that the probabilities within the prescribed domain are non-zero + and that the probabilities outside of the prescribed domain are zero. + """ + # d_distr_samples are (mdim, M) + # rect_domain is (mdim, 2) + inside = np.logical_and(np.all(np.greater_equal(self.d_distr_samples, + self.rect_domain[:, 0]), axis=1), + np.all(np.less_equal(self.d_distr_samples, + self.rect_domain[:, 1]), axis=1)) + msg = "Due to the inherent randomness of this method, this may fail." + print msg + print np.sum(self.rho_D_M[inside] >= 0.0) + assert np.sum(self.rho_D_M[inside] >= 0.0) < 100 + print np.sum(self.rho_D_M[np.logical_not(inside)] == 0.0) + assert np.sum(self.rho_D_M[np.logical_not(inside)] == 0.0) < 100 + + +class test_uniform_partition_uniform_distribution_rectangle_size_01D(data_01D, + uniform_partition_uniform_distribution_rectangle_size): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_size` on 01D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_size_01D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_size_01D, self).setUp() + + +class test_uniform_partition_uniform_distribution_rectangle_size_1D(data_1D, + uniform_partition_uniform_distribution_rectangle_size): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_size` on 1D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_size_1D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_size_1D, self).setUp() + + +class test_uniform_partition_uniform_distribution_rectangle_size_2D(data_2D, + uniform_partition_uniform_distribution_rectangle_size): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_size` on 2D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_size_2D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_size_2D, self).setUp() + + +class test_uniform_partition_uniform_distribution_rectangle_size_3D(data_3D, + uniform_partition_uniform_distribution_rectangle_size): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_size` on 3D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_size_3D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_size_3D, self).setUp() + + +class uniform_partition_uniform_distribution_rectangle_domain(prob_uniform): + """ + Set up :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_domain` on data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + if type(self.Q_ref) != np.array: + Q_ref = np.array([self.Q_ref]) + else: + Q_ref = self.Q_ref + if len(self.data_domain.shape) == 1: + data_domain = np.expand_dims(self.data_domain, axis=0) + else: + data_domain = self.data_domain + + self.rect_domain = np.zeros((data_domain.shape[0], 2)) + r_width = 0.1 * data_domain[:, 1] + + self.rect_domain[:, 0] = Q_ref - .5 * r_width + self.rect_domain[:, 1] = Q_ref + .5 * r_width + + self.data_prob = sFun.uniform_partition_uniform_distribution_rectangle_domain( + self.data, self.rect_domain.transpose(), M=67, num_d_emulate=1E3) + self.d_distr_samples = self.data_prob.get_values() + self.rho_D_M = self.data_prob.get_probabilities() + + def test_M(self): + """ + Test that the right number of d_distr_samples are used to create + rho_D_M. + """ + assert len(self.rho_D_M) == 67 + + def test_domain(self): + """ + Test that the probabilities within the prescribed domain are non-zero + and that the probabilities outside of the prescribed domain are zero. + """ + # d_distr_samples are (mdim, M) + # rect_domain is (mdim, 2) + inside = np.logical_and(np.all(np.greater_equal(self.d_distr_samples, + self.rect_domain[:, 0]), axis=1), + np.all(np.less_equal(self.d_distr_samples, + self.rect_domain[:, 1]), axis=1)) + msg = "Due to the inherent randomness of this method, this may fail." + print msg + print np.sum(self.rho_D_M[inside] >= 0.0) + assert np.sum(self.rho_D_M[inside] >= 0.0) < 100 + print np.sum(self.rho_D_M[np.logical_not(inside)] == 0.0) + assert np.sum(self.rho_D_M[np.logical_not(inside)] == 0.0) < 100 + + +class test_uniform_partition_uniform_distribution_rectangle_domain_01D(data_01D, + uniform_partition_uniform_distribution_rectangle_domain): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_domain` on 01D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_domain_01D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_domain_01D, self).setUp() + + +class test_uniform_partition_uniform_distribution_rectangle_domain_1D(data_1D, + uniform_partition_uniform_distribution_rectangle_domain): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_domain` on 1D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_domain_1D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_domain_1D, self).setUp() + + +class test_uniform_partition_uniform_distribution_rectangle_domain_2D(data_2D, + uniform_partition_uniform_distribution_rectangle_domain): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_domain` on 2D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_domain_2D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_domain_2D, self).setUp() + + +class test_uniform_partition_uniform_distribution_rectangle_domain_3D(data_3D, + uniform_partition_uniform_distribution_rectangle_domain): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_uniform_distribution_rectangle_domain` on 3D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_uniform_distribution_rectangle_domain_3D, self).createData() + super(test_uniform_partition_uniform_distribution_rectangle_domain_3D, self).setUp() From b59b4071d2855b6978fac53385e6d2d18358d444 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Sun, 22 May 2016 17:40:24 -0600 Subject: [PATCH 08/20] nonlinear Example default setup --- examples/nonlinearMap/myModel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/nonlinearMap/myModel.py b/examples/nonlinearMap/myModel.py index e7fda009..819257aa 100644 --- a/examples/nonlinearMap/myModel.py +++ b/examples/nonlinearMap/myModel.py @@ -23,7 +23,7 @@ (x1,y1)=(0.5,0.5) and (x2,y2)=(0.25,0.15) ''' # Choose the number of QoI -QoI_num = 2 +QoI_num = 1 # Specify the spatial points to take measurements of solution defining the QoI if QoI_num == 1: From 0388f453697cf00ff745cb10800b5684660ff997 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Sun, 22 May 2016 20:46:42 -0600 Subject: [PATCH 09/20] Final linear, nonlinear, and validation examples --- .../linearMap/linearMapUniformSampling.py | 2 +- .../nonlinearMapUniformSampling.py | 2 +- examples/validationExample/linearMap.py | 224 +++++++++--------- examples/validationExample/myModel.py | 11 + 4 files changed, 130 insertions(+), 109 deletions(-) create mode 100644 examples/validationExample/myModel.py diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index b6ff6d4f..77a766e2 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -1,6 +1,6 @@ #! /usr/bin/env python -# Copyright (C) 2014-2015 The BET Development Team +# Copyright (C) 2014-2016 The BET Development Team """ This example solves a stochastic inverse problem for a diff --git a/examples/nonlinearMap/nonlinearMapUniformSampling.py b/examples/nonlinearMap/nonlinearMapUniformSampling.py index a487f2b2..70d4857a 100644 --- a/examples/nonlinearMap/nonlinearMapUniformSampling.py +++ b/examples/nonlinearMap/nonlinearMapUniformSampling.py @@ -1,6 +1,6 @@ #! /usr/bin/env python -# Copyright (C) 2014-2015 The BET Development Team +# Copyright (C) 2014-2016 The BET Development Team r""" This example generates samples on a 2D grid and evaluates diff --git a/examples/validationExample/linearMap.py b/examples/validationExample/linearMap.py index 6010cfbe..b9175b64 100644 --- a/examples/validationExample/linearMap.py +++ b/examples/validationExample/linearMap.py @@ -1,6 +1,6 @@ #! /usr/bin/env python -# Copyright (C) 2014-2015 Lindley Graham and Steven Mattis +# Copyright (C) 2014-2016 The BET Development Team """ This 2D linear example verifies that geometrically distinct QoI can @@ -8,170 +8,180 @@ used to define the output probability measure. """ +from bet.Comm import comm, MPI import numpy as np import bet.calculateP as calculateP import bet.postProcess as postProcess import bet.calculateP.simpleFunP as simpleFunP import bet.calculateP.calculateP as calculateP import bet.postProcess.plotP as plotP +import bet.postProcess.plotDomains as plotD +import bet.sample as samp +import bet.sampling.basicSampling as bsam import scipy.spatial as spatial +from myModel import my_model -# parameter domain -lam_domain= np.array([[0.0, 1.0], - [0.0, 1.0]]) +# Initialize 3-dimensional input parameter sample set object +input_samples = samp.sample_set(2) + +# Set parameter domain +input_samples.set_domain(np.repeat([[0.0, 1.0]], 2, axis=0)) + +# Define the sampler that will be used to create the discretization +# object, which is the fundamental object used by BET to compute +# solutions to the stochastic inverse problem +sampler = bsam.sampler(my_model) ''' Suggested changes for user: - -Try setting n0 and n1 all to 10 and compare the results. - -Also, we can do uniform random sampling by setting - - random_sample = True - -If random_sample = True, consider defining - - n_samples = 2.5E3 - -Then also try n_samples = 1E4. What happens when n_samples = 1E2? -''' -random_sample = False -if random_sample == False: - n0 = 50 # number of samples in lam0 direction - n1 = 50 # number of samples in lam1 direction - n_samples = n0*n1 -else: - n_samples = 2.5E3 +Try with and without random sampling. +If using random sampling, try num_samples = 1E3 and 1E4. +What happens when num_samples = 1E2? +Try using 'lhs' instead of 'random' in the random_sample_set. -#set up samples -if random_sample == False: - vec0=list(np.linspace(lam_domain[0][0], lam_domain[0][1], n0)) - vec1 = list(np.linspace(lam_domain[1][0], lam_domain[1][1], n1)) - vecv0, vecv1 = np.meshgrid(vec0, vec1, indexing='ij') - samples=np.vstack((vecv0.flat[:], vecv1.flat[:])).transpose() +If using regular sampling, try different numbers of samples +per dimension. +''' +# Generate samples on the parameter space +randomSampling = False +if randomSampling is True: + sampler.random_sample_set('random', input_samples, num_samples=1E3) else: - samples = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, - num_l_emulate = n_samples) + sampler.regular_sample_set(input_samples, num_samples_per_dim=[30, 30]) + +''' +Suggested changes for user: -# QoI map -Q_map = np.array([[0.506, 0.463],[0.253, 0.918]]) +A standard Monte Carlo (MC) assumption is that every Voronoi cell +has the same volume. If a regular grid of samples was used, then +the standard MC assumption is true. +See what happens if the MC assumption is not assumed to be true, and +if different numbers of points are used to estimate the volumes of +the Voronoi cells. +''' +MC_assumption = True +# Estimate volumes of Voronoi cells associated with the parameter samples +if MC_assumption is False: + input_samples.estimate_volume(n_mc_points=1E5) +else: + input_samples.estimate_volume_mc() -# calc data -data = np.dot(samples,Q_map) +# Create the discretization object using the input samples +my_discretization = sampler.compute_QoI_and_create_discretization(input_samples, + savefile = 'Validation_discretization.txt.gz') ''' Compute the output distribution simple function approximation by propagating a different set of samples to implicitly define a Voronoi -discretization of the data space, and then propagating i.i.d. uniform -samples to bin into these cells. +discretization of the data space, corresponding to an implicitly defined +set of contour events defining a discretization of the input parameter +space. The probabilities of the Voronoi cells in the data space (and +thus the probabilities of the corresponding contour events in the +input parameter space) are determined by Monte Carlo sampling using +a set of i.i.d. uniform samples to bin into these cells. Suggested changes for user: -See the effect of using different values for d_distr_samples_num. -Choosing +See the effect of using different values for num_samples_discretize_D. +Choosing + + num_samples_discretize_D = 1 - d_distr_samples_num = 1 - produces exactly the right answer and is equivalent to assigning a -uniform probability to each data sample above (why?). +uniform probability to each data sample above (why?). -Try setting this to 2, 5, 10, 50, and 100. Can you explain what you +Try setting this to 2, 5, 10, 50, and 100. Can you explain what you are seeing? To see an exaggerated effect, try using random sampling -above with n_samples set to 1E2. +above with n_samples set to 1E2. ''' -d_distr_samples_num = 1 +num_samples_discretize_D = 5 +num_iid_samples = 1E5 -samples_discretize = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, - num_l_emulate = d_distr_samples_num) +SimpleFunction_set = samp.sample_set(2) +Monte_Carlo_set = samp.sample_set(2) -d_distr_samples = np.dot(samples_discretize, Q_map) +SimpleFunction_set.set_domain(np.repeat([[0.0, 1.0]], 2, axis=0)) +Monte_Carlo_set.set_domain(np.repeat([[0.0, 1.0]], 2, axis=0)) -d_Tree = spatial.KDTree(d_distr_samples) +SimpleFunction_discretization = sampler.create_random_discretization('random', + SimpleFunction_set, + num_samples=num_samples_discretize_D) -samples_distr_prob_num = d_distr_samples_num*1E3 +Monte_Carlo_discretization = sampler.create_random_discretization('random', + Monte_Carlo_set, + num_samples=num_iid_samples) -samples_distr_prob = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, - num_l_emulate = samples_distr_prob_num) +my_discretization._output_probability_set = SimpleFunction_discretization._output_sample_set +my_discretization._emulated_output_sample_set = Monte_Carlo_discretization._output_sample_set -data_prob = np.dot(samples_distr_prob, Q_map) +my_discretization.set_emulated_oo_ptr() -# Determine which data samples go to which d_distr_samples_num bins using the QoI -(_, oo_ptr) = d_Tree.query(data_prob) +my_discretization._output_probability_set.set_probabilities(np.zeros((num_samples_discretize_D,))) -# Calculate Probabilities of the d_distr_samples defined Voronoi cells -d_distr_prob = np.zeros((d_distr_samples_num,)) -for i in range(d_distr_samples_num): - Itemp = np.equal(oo_ptr, i) - Itemp_sum = float(np.sum(Itemp)) - d_distr_prob[i] = Itemp_sum / samples_distr_prob_num +for i in range(num_samples_discretize_D): + Itemp = np.equal(my_discretization.get_emulated_oo_ptr(), i) + Itemp_sum = float(np.sum(Itemp)) + Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM) + if Itemp_sum > 0: + my_discretization._output_probability_set._probabilities[i] = Itemp_sum / num_iid_samples -''' -Suggested changes for user: - -If using a regular grid of sampling (if random_sample = False), we set - - lambda_emulate = samples - -Otherwise, play around with num_l_emulate. A value of 1E2 will probably -give poor results while results become fairly consistent with values -that are approximately 10x the number of samples. - -Note that you can always use - - lambda_emulate = samples - -and this simply will imply that a standard Monte Carlo assumption is -being used, which in a measure-theoretic context implies that each -Voronoi cell is assumed to have the same measure. This type of -approximation is more reasonable for large n_samples due to the slow -convergence rate of Monte Carlo (it converges like 1/sqrt(n_samples)). -''' -if random_sample == False: - lambda_emulate = samples -else: - lambda_emulate = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, num_l_emulate = 1E5) +# Calculate probabilities +calculateP.prob(my_discretization) +# Show some plots of the different sample sets +plotD.scatter_2D(my_discretization._input_sample_set, filename = 'Parameter_Samples.eps') +plotD.scatter_2D(my_discretization._output_sample_set, filename = 'QoI_Samples.eps') +plotD.scatter_2D(my_discretization._output_probability_set, filename = 'Data_Space_Discretization.eps') -# calculate probablities -(P, lambda_emulate, io_ptr, emulate_ptr) = calculateP.prob_emulated(samples=samples, - data=data, - rho_D_M=d_distr_prob, - d_distr_samples=d_distr_samples, - lambda_emulate=lambda_emulate, - d_Tree=d_Tree) -# calculate 2d marginal probs +######################################## +# Post-process the results +######################################## ''' Suggested changes for user: - + At this point, the only thing that should change in the plotP.* inputs should be either the nbins values or sigma (which influences the kernel density estimation with smaller values implying a density estimate that looks more like a histogram and larger values smoothing out the values more). - + There are ways to determine "optimal" smoothing parameters (e.g., see CV, GCV, and other similar methods), but we have not incorporated these into the code -as lower-dimensional marginal plots have limited value in understanding the -structure of a high dimensional non-parametric probability measure. +as lower-dimensional marginal plots generally have limited value in understanding +the structure of a high dimensional non-parametric probability measure. ''' -(bins, marginals2D) = plotP.calculate_2D_marginal_probs(P_samples = P, samples = lambda_emulate, lam_domain = lam_domain, nbins = [10, 10]) +# calculate 2d marginal probs +(bins, marginals2D) = plotP.calculate_2D_marginal_probs(input_samples, + nbins = [30, 30]) + +# plot 2d marginals probs +plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "validation_raw", + file_extension = ".eps", plot_surface=False) + # smooth 2d marginals probs (optional) -#marginals2D = plotP.smooth_marginals_2D(marginals2D,bins, sigma=0.01) +marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.1) # plot 2d marginals probs -plotP.plot_2D_marginal_probs(marginals2D, bins, lam_domain, filename = "linearMapValidation", - plot_surface=True, interactive=True, lambda_label = None, fileExtension = ".png") +plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "validation_smooth", + file_extension = ".eps", plot_surface=False) # calculate 1d marginal probs -(bins, marginals1D) = plotP.calculate_1D_marginal_probs(P_samples = P, samples = lambda_emulate, lam_domain = lam_domain, nbins = [10, 10]) +(bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, + nbins = [30, 30]) + +# plot 2d marginal probs +plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "validation_raw", + file_extension = ".eps") + # smooth 1d marginal probs (optional) -#marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.01) -# plot 1d marginal probs -plotP.plot_1D_marginal_probs(marginals1D, bins, lam_domain, filename = "linearMapValidation", lambda_label = None, fileExtension = ".png") +marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.1) +# plot 2d marginal probs +plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "validation_smooth", + file_extension = ".eps") diff --git a/examples/validationExample/myModel.py b/examples/validationExample/myModel.py new file mode 100644 index 00000000..60c8b36f --- /dev/null +++ b/examples/validationExample/myModel.py @@ -0,0 +1,11 @@ +# Copyright (C) 2016 The BET Development Team + +# -*- coding: utf-8 -*- +import numpy as np + +# Define a model that is a linear QoI map +def my_model(parameter_samples): + Q_map = np.array([[0.506, 0.463], [0.253, 0.918]]) + QoI_samples = data= np.dot(parameter_samples,Q_map) + return QoI_samples + From 3168f991c93453e42a622b2180bd5880118f59a9 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Sun, 22 May 2016 21:20:34 -0600 Subject: [PATCH 10/20] Fixed a previous merge typo --- bet/calculateP/simpleFunP.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 6e3b6368..67064ed8 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -76,8 +76,6 @@ def uniform_partition_uniform_distribution_rectangle_size(data_set, Q_ref, msg += "bet.sample.discretization or np.ndarray" raise wrong_argument_type(msg) - bin_size = (np.max(values, 0) - np.min(values, 0))*bin_ratio - if not isinstance(rect_size, collections.Iterable): rect_size = rect_size * np.ones((dim,)) if np.any(np.less(rect_size, 0)): From 72925b2f7ba3a83c462e4720d1e47bfaa9ea4735 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Mon, 23 May 2016 11:24:29 -0600 Subject: [PATCH 11/20] Updating for user defined simpleFunP --- bet/calculateP/simpleFunP.py | 97 ++++++++++++++++++++++++- examples/validationExample/linearMap.py | 19 ++--- 2 files changed, 99 insertions(+), 17 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 67064ed8..5b5fdd9a 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -372,7 +372,8 @@ def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domai ``len(d_distr_samples) == 3**mdim``. :param data_set: Sample set that the probability measure is defined for. - :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :type data_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param rect_domain: The domain overwhich :math:`\rho_\mathcal{D}` is uniform. :type rect_domain: :class:`numpy.ndarray` of shape (2, mdim) @@ -505,7 +506,7 @@ def uniform_partition_uniform_distribution_data_samples(data_set): s_set.set_probabilities(np.ones((num,), dtype=np.float)/num) if isinstance(data_set, samp.discretization): - data_set._output_sample_set = s_set + data_set._output_probability_set = s_set return s_set @@ -599,7 +600,7 @@ def normal_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate= # above, while informed by the sampling of the map Q, do not require # solving the model EVER! This can be done "offline" so to speak. if isinstance(data_set, samp.discretization): - data_set._output_sample_set = s_set + data_set._output_probability_set = s_set return s_set @@ -678,3 +679,93 @@ def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate if isinstance(data_set, samp.discretization): data_set._output_probability_set = s_set return s_set + +def user_discretization_user_distribution(data_set, data_discretization_set, + data_distribution_set): + r""" + Creates a user defined simple function approximation of a user + defined distribution. The simple function discretization is + specified in the data_discretization_set, and the set of i.i.d. + samples from the distribution is specified in the + data_distribution_set. + + :param data_set: + :param data_discretization_set: + :param data_distribution_set: + :return: + """ + if isinstance(data_set, samp.sample_set_base): + s_set = data_set.copy() + dim = s_set._dim + elif isinstance(data_set, samp.discretization): + s_set = data_set._output_sample_set.copy() + dim = s_set._dim + elif isinstance(data_set, np.ndarray): + dim = data_set.shape[1] + values = data_set + s_set = samp.sample_set(dim=dim) + s_set.set_values(values) + else: + msg = "The first argument must be of type bet.sample.sample_set, " + msg += "bet.sample.discretization or np.ndarray" + raise wrong_argument_type(msg) + + if isinstance(data_discretization_set, samp.sample_set_base): + num_samples_discretize_D = data_discretization_set.check_num() + simpleFun_set = data_discretization_set.copy() + dim_simpleFun = simpleFun_set._dim + elif isinstance(data_discretization_set, samp.discretization): + num_samples_discretize_D = data_discretization_set.check_nums() + simpleFun_set = data_discretization_set._output_sample_set.copy() + dim_simpleFun = simpleFun_set._dim + elif isinstance(data_discretization_set, np.ndarray): + num_samples_discretize_D = data_discretization_set.shape[0] + dim_simpleFun = data_discretization_set.shape[1] + simpleFun_set = samp.sample_set(dim=dim_simpleFun) + simpleFun_set.set_values(values) + else: + msg = "The second argument must be of type bet.sample.sample_set, " + msg += "bet.sample.discretization or np.ndarray" + raise wrong_argument_type(msg) + + if isinstance(data_distribution_set, samp.sample_set_base): + MonteCarlo_set = data_distribution_set.copy() + dim_MonteCarlo = MonteCarlo_set._dim + num_iid_samples = data_distribution_set.check_num() + elif isinstance(data_distribution_set, samp.discretization): + MonteCarlo_set = data_distribution_set._output_sample_set.copy() + dim_MonteCarlo = MonteCarlo_set._dim + num_iid_samples = data_distribution_set.check_nums() + elif isinstance(data_distribution_set, np.ndarray): + num_iid_samples = data_distribution_set.shape[0] + dim_MonteCarlo = data_distribution_set.shape[1] + MonteCarlo_set = samp.sample_set(dim=dim_MonteCarlo) + MonteCarlo_set.set_values(values) + else: + msg = "The second argument must be of type bet.sample.sample_set, " + msg += "bet.sample.discretization or np.ndarray" + raise wrong_argument_type(msg) + + if np.not_equal(dim_MonteCarlo, dim) or np.not_equal(dim_simpleFun, dim): + msg = "The argument types have conflicting dimensions" + raise wrong_argument_type(msg) + + my_discretization = samp.discretization(input_sample_set=s_set, + output_sample_set=s_set, + emulated_output_sample_set=MonteCarlo_set, + output_probability_set=simpleFun_set) + + my_discretization.set_emulated_oo_ptr() + + my_discretization._output_probability_set.set_probabilities(np.zeros((num_samples_discretize_D,))) + + for i in range(num_samples_discretize_D): + Itemp = np.equal(my_discretization.get_emulated_oo_ptr(), i) + Itemp_sum = float(np.sum(Itemp)) + Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM) + if Itemp_sum > 0: + my_discretization._output_probability_set._probabilities[i] = Itemp_sum / num_iid_samples + + if isinstance(data_set, samp.discretization): + data_set._output_probability_set = my_discretization._output_probability_set + return s_set diff --git a/examples/validationExample/linearMap.py b/examples/validationExample/linearMap.py index b9175b64..f5f19acf 100644 --- a/examples/validationExample/linearMap.py +++ b/examples/validationExample/linearMap.py @@ -97,7 +97,7 @@ are seeing? To see an exaggerated effect, try using random sampling above with n_samples set to 1E2. ''' -num_samples_discretize_D = 5 +num_samples_discretize_D = 1 num_iid_samples = 1E5 SimpleFunction_set = samp.sample_set(2) @@ -114,19 +114,10 @@ Monte_Carlo_set, num_samples=num_iid_samples) -my_discretization._output_probability_set = SimpleFunction_discretization._output_sample_set -my_discretization._emulated_output_sample_set = Monte_Carlo_discretization._output_sample_set - -my_discretization.set_emulated_oo_ptr() - -my_discretization._output_probability_set.set_probabilities(np.zeros((num_samples_discretize_D,))) - -for i in range(num_samples_discretize_D): - Itemp = np.equal(my_discretization.get_emulated_oo_ptr(), i) - Itemp_sum = float(np.sum(Itemp)) - Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM) - if Itemp_sum > 0: - my_discretization._output_probability_set._probabilities[i] = Itemp_sum / num_iid_samples +# Compute the simple function approximation to the distribution on the data space +simpleFunP.user_discretization_user_distribution(my_discretization, + SimpleFunction_discretization, + Monte_Carlo_discretization) # Calculate probabilities calculateP.prob(my_discretization) From 5ae6dc4cdca9726dd32a153e59f41e40eae24dd3 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Mon, 23 May 2016 11:30:45 -0600 Subject: [PATCH 12/20] Updating to some commenting --- bet/calculateP/simpleFunP.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 5b5fdd9a..60cd59cd 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -45,8 +45,8 @@ def uniform_partition_uniform_distribution_rectangle_size(data_set, Q_ref, :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - play around with it and you can get reasonable results with a relatively small number here like 50. - :param rect_size: The scale used to determine the support of the - uniform distribution as ``rect_size = (data_max-data_min)*rect_scale`` + :param rect_size: Determines the size of the support of the + uniform distribution on a generalized rectangle :type rect_size: double or list() :param int num_d_emulate: Number of samples used to emulate using an MC assumption @@ -248,7 +248,8 @@ def uniform_partition_uniform_distribution_rectangle_domain(data_set, rect_domai :param int num_d_emulate: Number of samples used to emulate using an MC assumption :param data_set: Sample set that the probability measure is defined for. - :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :type data_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param Q_ref: :math:`Q(`\lambda_{reference})` :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) @@ -279,9 +280,8 @@ def uniform_partition_uniform_distribution_rectangle_domain(data_set, rect_domai domain_lengths = np.max(rect_domain, 0) - np.min(rect_domain, 0) - return uniform_partition_uniform_distribution_rectangle_size(data_set, domain_center, - domain_lengths, M, - num_d_emulate) + return uniform_partition_uniform_distribution_rectangle_size(data_set, + domain_center, domain_lengths, M, num_d_emulate) def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_size, @@ -408,9 +408,10 @@ def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domai domain_lengths = np.max(rect_domain, 0) - np.min(rect_domain, 0) return regular_partition_uniform_distribution_rectangle_size(data_set, domain_center, - domain_lengths, center_pts_per_edge) + domain_lengths, center_pts_per_edge) -def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rect_scale, center_pts_per_edge=1): +def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, + rect_scale, center_pts_per_edge=1): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density @@ -422,7 +423,8 @@ def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec ``len(d_distr_samples) == 3^mdim``. :param data_set: Sample set that the probability measure is defined for. - :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :type data_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param rect_scale: The scale used to determine the width of the uniform distributiion as ``rect_size = (data_max-data_min)*rect_scale`` :type rect_scale: double or list() @@ -460,8 +462,8 @@ def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, rec rect_scale = rect_scale*np.ones((dim, )) rect_size = (np.max(data, 0) - np.min(data, 0))*rect_scale - return regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_size, - center_pts_per_edge) + return regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, + rect_size, center_pts_per_edge) def uniform_partition_uniform_distribution_data_samples(data_set): r""" @@ -613,7 +615,8 @@ def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate direction. :param data_set: Sample set that the probability measure is defined for. - :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :type data_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param int M: Defines number M samples in D used to define :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" - play around with it and you can get reasonable results with a @@ -757,14 +760,16 @@ def user_discretization_user_distribution(data_set, data_discretization_set, my_discretization.set_emulated_oo_ptr() - my_discretization._output_probability_set.set_probabilities(np.zeros((num_samples_discretize_D,))) + my_discretization._output_probability_set.set_probabilities( + np.zeros((num_samples_discretize_D,))) for i in range(num_samples_discretize_D): Itemp = np.equal(my_discretization.get_emulated_oo_ptr(), i) Itemp_sum = float(np.sum(Itemp)) Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM) if Itemp_sum > 0: - my_discretization._output_probability_set._probabilities[i] = Itemp_sum / num_iid_samples + my_discretization._output_probability_set._probabilities[i] = \ + Itemp_sum / num_iid_samples if isinstance(data_set, samp.discretization): data_set._output_probability_set = my_discretization._output_probability_set From 502a0f42c473d376e84a2fe949bd97879e95df4d Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Mon, 23 May 2016 11:55:30 -0600 Subject: [PATCH 13/20] Small changes to formatting --- bet/calculateP/simpleFunP.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 60cd59cd..789f42a9 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -360,7 +360,8 @@ def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_ return s_set -def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domain, center_pts_per_edge=1): +def regular_partition_uniform_distribution_rectangle_domain(data_set, + rect_domain, center_pts_per_edge=1): r""" Creates a simple function appoximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho{\mathcal{D}, M}` is a uniform probablity density over the @@ -407,11 +408,11 @@ def regular_partition_uniform_distribution_rectangle_domain(data_set, rect_domai domain_center = np.mean(rect_domain, 0) domain_lengths = np.max(rect_domain, 0) - np.min(rect_domain, 0) - return regular_partition_uniform_distribution_rectangle_size(data_set, domain_center, - domain_lengths, center_pts_per_edge) + return regular_partition_uniform_distribution_rectangle_size(data_set, + domain_center, domain_lengths, center_pts_per_edge) def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, - rect_scale, center_pts_per_edge=1): + rect_scale, center_pts_per_edge=1): r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density From f555e59a27e32e761d56d1b711c0eeb55cc535a0 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Mon, 23 May 2016 12:13:58 -0600 Subject: [PATCH 14/20] Mainly updating to commenting --- bet/calculateP/simpleFunP.py | 70 +++++++++++++++++++++--------------- 1 file changed, 41 insertions(+), 29 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 789f42a9..b0ff707b 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -25,15 +25,15 @@ def uniform_partition_uniform_distribution_rectangle_size(data_set, Q_ref, r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` where :math:`\rho_{\mathcal{D}}` is a uniform probability density on - a generalized rectangle centered at Q_ref. - The support of this density is defined by rect_size, which determines + a generalized rectangle centered at ``Q_ref``. + The support of this density is defined by ``rect_size``, which determines the size of the generalized rectangle. - The simple function approximation is then defined by determining M + The simple function approximation is then defined by determining ``M`` Voronoi cells (i.e., "bins") partitioning :math:`\mathcal{D}`. These - bins are only implicitly defined by M samples in :math:`\mathcal{D}`. + bins are only implicitly defined by ``M`` samples in :math:`\mathcal{D}`. Finally, the probabilities of each of these bins is computed by sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor - searches to bin these samples in the M implicitly defined bins. + searches to bin these samples in the ``M`` implicitly defined bins. The result is the simple function approximation denoted by :math:`\rho_{\mathcal{D},M}`. @@ -161,16 +161,16 @@ def uniform_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` where :math:`\rho_{\mathcal{D}}` is a uniform probability density on - a generalized rectangle centered at Q_ref. - The support of this density is defined by rect_scale, which determines + a generalized rectangle centered at ``Q_ref``. + The support of this density is defined by ``rect_scale``, which determines the size of the generalized rectangle by scaling the circumscribing generalized rectangle of :math:`\mathcal{D}`. - The simple function approximation is then defined by determining M + The simple function approximation is then defined by determining ``M `` Voronoi cells (i.e., "bins") partitioning :math:`\mathcal{D}`. These - bins are only implicitly defined by M samples in :math:`\mathcal{D}`. + bins are only implicitly defined by ``M`` samples in :math:`\mathcal{D}`. Finally, the probabilities of each of these bins is computed by sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor - searches to bin these samples in the M implicitly defined bins. + searches to bin these samples in the ``M`` implicitly defined bins. The result is the simple function approximation denoted by :math:`\rho_{\mathcal{D},M}`. @@ -225,13 +225,13 @@ def uniform_partition_uniform_distribution_rectangle_domain(data_set, rect_domai r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D}}` where :math:`\rho_{\mathcal{D}}` is a uniform probability density on - a generalized rectangle defined by rect_domain. - The simple function approximation is then defined by determining M + a generalized rectangle defined by ``rect_domain``. + The simple function approximation is then defined by determining ``M`` Voronoi cells (i.e., "bins") partitioning :math:`\mathcal{D}`. These - bins are only implicitly defined by M samples in :math:`\mathcal{D}`. + bins are only implicitly defined by ``M ``samples in :math:`\mathcal{D}`. Finally, the probabilities of each of these bins is computed by sampling from :math:`\rho{\mathcal{D}}` and using nearest neighbor - searches to bin these samples in the M implicitly defined bins. + searches to bin these samples in the ``M`` implicitly defined bins. The result is the simple function approximation denoted by :math:`\rho_{\mathcal{D},M}`. @@ -289,7 +289,7 @@ def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_ r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density - centered at Q_ref with rect_size of the width of a hyperrectangle. + centered at ``Q_ref`` with ``rect_size`` of the width of a hyperrectangle. Since rho_D is a uniform distribution on a hyperrectanlge we should be able to represent it exactly with ``M = 3^mdim`` or rather @@ -365,7 +365,7 @@ def regular_partition_uniform_distribution_rectangle_domain(data_set, r""" Creates a simple function appoximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho{\mathcal{D}, M}` is a uniform probablity density over the - hyperrectangular domain specified by domain. + hyperrectangular domain specified by ``rect_domain``. Since :math:`\rho_\mathcal{D}` is a uniform distribution on a hyperrectangle we should we able to represent it exactly with @@ -416,7 +416,7 @@ def regular_partition_uniform_distribution_rectangle_scaled(data_set, Q_ref, r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a uniform probability density - centered at Q_ref with rect_scale of the width + centered at ``Q_ref`` with ``rect_scale`` of the width of D. Since rho_D is a uniform distribution on a hyperrectanlge we should be able @@ -517,8 +517,8 @@ def normal_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate= r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability - density centered at Q_ref with standard deviation std using M bins sampled - from the given normal distribution. + density centered at ``Q_ref`` with standard deviation ``std`` using + ``M`` bins sampled from the given normal distribution. :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or @@ -611,9 +611,9 @@ def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate r""" Creates a simple function approximation of :math:`\rho_{\mathcal{D},M}` where :math:`\rho_{\mathcal{D},M}` is a multivariate normal probability - density centered at Q_ref with standard deviation std using M bins sampled - from a uniform distribution with a size 4 standard deviations in each - direction. + density centered at ``Q_ref`` with standard deviation ``std`` using + ``M`` bins sampled from a uniform distribution with a size 4 standard + deviations in each direction. :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or @@ -689,14 +689,26 @@ def user_discretization_user_distribution(data_set, data_discretization_set, r""" Creates a user defined simple function approximation of a user defined distribution. The simple function discretization is - specified in the data_discretization_set, and the set of i.i.d. + specified in the ``data_discretization_set``, and the set of i.i.d. samples from the distribution is specified in the - data_distribution_set. + ``data_distribution_set``. - :param data_set: - :param data_discretization_set: - :param data_distribution_set: - :return: + :param data_set: Sample set that the probability measure is defined for. + :type data_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :param data_discretization_set: Sample set defining the discretization + of the data space into Voronoi cells for which a simple function + is defined upon. + :type data_discretization_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + :param data_distribution_set: Sample set containing the i.i.d. samples + from the distribution on the data space that are binned within the + Voronoi cells implicitly defined by the data_discretization_set. + :type data_distribution_set: :class:`~bet.sample.discretization` or + :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` + + :rtype: :class:`~bet.sample.voronoi_sample_set` + :returns: sample_set object defininng simple function approximation """ if isinstance(data_set, samp.sample_set_base): s_set = data_set.copy() @@ -774,4 +786,4 @@ def user_discretization_user_distribution(data_set, data_discretization_set, if isinstance(data_set, samp.discretization): data_set._output_probability_set = my_discretization._output_probability_set - return s_set + return my_discretization._output_probability_set From 3b01efb8b8133b5739d8f068ca2936f3d5206fdd Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Mon, 23 May 2016 14:59:25 -0600 Subject: [PATCH 15/20] FEniCS example done --- examples/FEniCS/BET_script.py | 162 ++++ examples/FEniCS/Lshaped.xml | 1275 +++++++++++++++++++++++++++ examples/FEniCS/meshDS.py | 83 ++ examples/FEniCS/myModel.py | 182 ++++ examples/FEniCS/poissonRandField.py | 16 + examples/FEniCS/projectKL.py | 168 ++++ 6 files changed, 1886 insertions(+) create mode 100644 examples/FEniCS/BET_script.py create mode 100644 examples/FEniCS/Lshaped.xml create mode 100644 examples/FEniCS/meshDS.py create mode 100644 examples/FEniCS/myModel.py create mode 100644 examples/FEniCS/poissonRandField.py create mode 100644 examples/FEniCS/projectKL.py diff --git a/examples/FEniCS/BET_script.py b/examples/FEniCS/BET_script.py new file mode 100644 index 00000000..86464998 --- /dev/null +++ b/examples/FEniCS/BET_script.py @@ -0,0 +1,162 @@ +#! /usr/bin/env python + +# Copyright (C) 2014-2016 The BET Development Team + +r""" +An installation of FEniCS using the same python as used for +installing BET is required to run this example. + +This example generates samples for a KL expansion associated with +a covariance defined by ``cov`` in myModel.py that on an L-shaped +mesh defining the permeability field for a Poisson equation. + +The quantities of interest (QoI) are defined as two spatial +averages of the solution to the PDE. + +The user defines the dimension of the parameter space (corresponding +to the number of KL terms) and the number of samples in this space. +""" + +import numpy as np +import bet.calculateP as calculateP +import bet.postProcess as postProcess +import bet.calculateP.simpleFunP as simpleFunP +import bet.calculateP.calculateP as calculateP +import bet.postProcess.plotP as plotP +import bet.postProcess.plotDomains as plotD +import bet.sample as samp +import bet.sampling.basicSampling as bsam +from myModel import my_model + +# Initialize input parameter sample set object +num_KL_terms = 2 +input_samples = samp.sample_set(2) + +# Set parameter domain +KL_term_min = -3.0 +KL_term_max = 3.0 +input_samples.set_domain(np.repeat([[KL_term_min, KL_term_max]], + num_KL_terms, + axis=0)) + +# Define the sampler that will be used to create the discretization +# object, which is the fundamental object used by BET to compute +# solutions to the stochastic inverse problem +sampler = bsam.sampler(my_model) + +''' +Suggested changes for user: + +Try with and without random sampling. + +If using random sampling, try num_samples = 1E3 and 1E4. +What happens when num_samples = 1E2? +Try using 'lhs' instead of 'random' in the random_sample_set. + +If using regular sampling, try different numbers of samples +per dimension. +''' +# Generate samples on the parameter space +randomSampling = False +if randomSampling is True: + sampler.random_sample_set('random', input_samples, num_samples=1E4) +else: + sampler.regular_sample_set(input_samples, num_samples_per_dim=[50, 50]) + +''' +Suggested changes for user: + +A standard Monte Carlo (MC) assumption is that every Voronoi cell +has the same volume. If a regular grid of samples was used, then +the standard MC assumption is true. + +See what happens if the MC assumption is not assumed to be true, and +if different numbers of points are used to estimate the volumes of +the Voronoi cells. +''' +MC_assumption = True +# Estimate volumes of Voronoi cells associated with the parameter samples +if MC_assumption is False: + input_samples.estimate_volume(n_mc_points=1E5) +else: + input_samples.estimate_volume_mc() + +# Create the discretization object using the input samples +my_discretization = sampler.compute_QoI_and_create_discretization(input_samples, + savefile='FEniCS_Example.txt.gz') + +''' +Suggested changes for user: + +Try different reference parameters. +''' +# Define the reference parameter +#param_ref = np.zeros((1,num_KL_terms)) +param_ref = np.ones((1,num_KL_terms)) + +# Compute the reference QoI +Q_ref = my_model(param_ref) + +# Create some plots of input and output discretizations +plotD.scatter_2D(input_samples, p_ref=param_ref[0,:], filename='FEniCS_ParameterSamples.eps') +if Q_ref.size == 2: + plotD.show_data(my_discretization, Q_ref=Q_ref[0,:]) + +''' +Suggested changes for user: + +Try different ways of discretizing the probability measure on D defined +as a uniform probability measure on a rectangle or interval depending +on choice of QoI_num in myModel.py. +''' +randomDataDiscretization = False +if randomDataDiscretization is False: + simpleFunP.regular_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref[0,:], rect_scale=0.1, + center_pts_per_edge=3) +else: + simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled( + data_set=my_discretization, Q_ref=Q_ref[0,:], rect_scale=0.1, + M=50, num_d_emulate=1E5) + +# calculate probablities +calculateP.prob(my_discretization) + +######################################## +# Post-process the results +######################################## +''' +Suggested changes for user: + +At this point, the only thing that should change in the plotP.* inputs +should be either the nbins values or sigma (which influences the kernel +density estimation with smaller values implying a density estimate that +looks more like a histogram and larger values smoothing out the values +more). + +There are ways to determine "optimal" smoothing parameters (e.g., see CV, GCV, +and other similar methods), but we have not incorporated these into the code +as lower-dimensional marginal plots generally have limited value in understanding +the structure of a high dimensional non-parametric probability measure. +''' +# calculate 2d marginal probs +(bins, marginals2D) = plotP.calculate_2D_marginal_probs(input_samples, + nbins=20) +# smooth 2d marginals probs (optional) +marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.5) + +# plot 2d marginals probs +plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename="FEniCS", + file_extension=".eps", plot_surface=False) + +# calculate 1d marginal probs +(bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, + nbins=20) +# smooth 1d marginal probs (optional) +marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.5) +# plot 2d marginal probs +plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename="FEniCS", + file_extension=".eps") + + + diff --git a/examples/FEniCS/Lshaped.xml b/examples/FEniCS/Lshaped.xml new file mode 100644 index 00000000..4cfe598b --- /dev/null +++ b/examples/FEniCS/Lshaped.xml @@ -0,0 +1,1275 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/FEniCS/meshDS.py b/examples/FEniCS/meshDS.py new file mode 100644 index 00000000..cadfccf1 --- /dev/null +++ b/examples/FEniCS/meshDS.py @@ -0,0 +1,83 @@ +#!/usr/bin/en python + +from dolfin import * +from numpy import * + + +class meshDS(object): + + """Docstring for meshDS. """ + + def __init__(self,mesh): + """TODO: to be defined1. + + :mesh: reads a fenics mesh object + + """ + self._mesh = mesh + self.node_elem = {} # empty dictionary of node to elements connectivity + self.edges_elem = {} # empty dictionary of edges to elements connectivity + + # initialize the mesh and read in the values + self._mesh.init() + self._dim = self._mesh.topology().dim() + self.num_nodes = self._mesh.num_vertices() + self.num_elements = self._mesh.num_cells() + self.num_edges = self._mesh.num_edges() + + def getNodes(self): + """TODO: Docstring for getNodes. + :returns: num of nodes in the mesh + + """ + return self.num_nodes + def getElements(self): + """TODO: Docstring for getElements. + :returns: number of elements in the mesh + + """ + return self.num_elements + + def getEdges(self): + """TODO: Docstring for getElements. + :returns: number of elements in the mesh + + """ + return self.num_edges + def getElemToNodes(self): + """TODO: Docstring for getElemToNodes. + :returns: Elements - Nodes Connectivity array of array + + """ + return self._mesh.cells() + def getNodesToElem(self): + """TODO: Docstring for getNodesToElem. + :returns: returns Nodes to Element connectivity as a dictionary + where nodes_elem[i] is an array of all the elements attached to node i + + """ + for nodes in entities(self._mesh,0): + self.node_elem[nodes.index()] = nodes.entities(self._dim) + return self.node_elem + def getElemVCArray(self): + + """TODO: Docstring for getElemVCArray. + :returns: array of element volume and and an array of element centroid object + Thus elem_centroid_array[i][0] means the x co-ordinate of the centroid for element number i + Thus elem_centroid_array[i][1] means the y co-ordinate of the centroid for element number i + """ + + elem_vol_array = empty((self.num_elements),dtype=float) + elem_centroid_array = empty((self.num_elements),dtype=object) + + cell_indx = 0 + for node_list in self._mesh.cells(): + # First get the cell object corresponding to the cell_indx + cell_obj = Cell(self._mesh,cell_indx) + # Find the cell volume and cell centroid + elem_vol_array[cell_indx] = cell_obj.volume() + elem_centroid_array[cell_indx] = cell_obj.midpoint() + # update cell index + cell_indx = cell_indx + 1 + return elem_vol_array,elem_centroid_array + diff --git a/examples/FEniCS/myModel.py b/examples/FEniCS/myModel.py new file mode 100644 index 00000000..d72f7d60 --- /dev/null +++ b/examples/FEniCS/myModel.py @@ -0,0 +1,182 @@ +# Copyright (C) 2016 The BET Development Team + +# -*- coding: utf-8 -*- +import numpy as np +from dolfin import * +from meshDS import meshDS +from projectKL import projectKL +from poissonRandField import solvePoissonRandomField +import scipy.io as sio + +def my_model(parameter_samples): + # number of parameter samples + numSamples = parameter_samples.shape[0] + + # number of KL expansion terms. + numKL = parameter_samples.shape[1] + + # the samples are the coefficients of the KL expansion typically denoted by xi_k + xi_k = parameter_samples + + ''' + ++++++++++++++++ Steps in Computing the Numerical KL Expansion ++++++++++ + We proceed by loading the mesh and defining the function space for which + the eigenfunctions are defined upon. + + Then, we define the covariance kernel which requires correlation lengths + and a standard deviation. + + We then compute the truncated KL expansion. + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ''' + + # Step 1: Set up the Mesh and Function Space + + mesh = Mesh("Lshaped.xml") + #mesh = RectangleMesh(0,0,10,10,20,20) + + # Plot the mesh for visual check + #plot(mesh,interactive=True) + + # initialize the mesh to generate connectivity + mesh.init() + + # Random field is projected on the space of Hat functions in the mesh + V = FunctionSpace(mesh, "CG", 1) + + # Step 2: Project covariance in the mesh and get the eigenfunctions + + # Initialize the projectKL object with the mesh + Lmesh = projectKL(mesh) + + # Create the covariance expression to project on the mesh. + etaX = 10.0 + etaY = 10.0 + C = 1 + + # Pick your favorite covariance. Popular choices are Gaussian (of course), + # Exponential, triangular (has finite support which is nice). Check out + # Ghanem and Spanos' book for more classical options. + + # A Gaussian Covariance + ''' + cov = Expression("C*exp(-((x[0]-x[1]))*((x[0]-x[1]))/ex - \ + ((x[2]-x[3]))*((x[2]-x[3]))/ey)", + ex=etaX,ey=etaY, C=C) + ''' + # An Exponential Covariance + cov = Expression("C*exp(-fabs(x[0]-x[1])/ex - fabs(x[2]-x[3])/ey)",ex=etaX,ey=etaY, C=C) + + # Solve the discrete covariance relation on the mesh + Lmesh.projectCovToMesh(numKL,cov) + + # Get the eigenfunctions and eigenvalues + eigen_func = Lmesh.eigen_funcs + eigen_val = Lmesh.eigen_vals + + #print 'eigen_vals' + #print eigen_val + #print eigen_val.sum() + + ''' + ++++++++++++++++ Steps in Solving Poisson with the KL fields ++++++++++++ + First set up the necessary variables and boundary conditions for the + problem. + + Then create the QoI maps defined by average values over some part of the + physical domain. + + Loop through the sample fields and create the permeability defined by the + exponential of the KL field (i.e., the KL expansion represents the log of + the permeability). + + For each sample field, call the PoissonRandomField solver. + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ''' + # permeability + perm_k = Function(V) + + # Create Boundary Conditions -- Dirichlet on left and bottom boundary. + # Left Dirichlet Bc + def left_boundary(x,on_boundary): + """TODO: Docstring for left_boundary. + + :x: TODO + :on_boundary: TODO + :returns: TODO + + """ + tol = 1e-14 + return on_boundary and abs(x[0]) < tol + Gamma_0 = DirichletBC(V,Constant(0.0),left_boundary) + + def bottom_boundary(x,on_boundary): + """TODO: Docstring for left_boundary. + + :x: TODO + :on_boundary: TODO + :returns: TODO + + """ + tol = 1e-14 + return on_boundary and abs(x[1]) < tol + Gamma_1 = DirichletBC(V,Constant(0.0),bottom_boundary) + bcs = [Gamma_0,Gamma_1] + + # Setup adjoint boundary conditions. + Gamma_adj_0 = DirichletBC(V,Constant(0.0),left_boundary) + Gamma_adj_1 = DirichletBC(V,Constant(0.0),bottom_boundary) + bcs_adj = [Gamma_adj_0, Gamma_adj_1] + + # Setup the QoI class + class CharFunc(Expression): + def __init__(self, region): + self.a = region[0] + self.b = region[1] + self.c = region[2] + self.d = region[3] + def eval(self, v, x): + v[0] = 0 + if (x[0] >= self.a) & (x[0] <= self.b) & (x[1] >= self.c) & (x[1] <= self.d): + v[0] = 1 + return v + + # Define the QoI maps + Chi_1 = CharFunc([0.75, 1.25, 7.75, 8.25]) + Chi_2 = CharFunc([7.75, 8.25, 0.75, 1.25]) + + QoI_samples = np.zeros([numSamples,2]) + + QoI_deriv_1 = np.zeros([numSamples,numKL]) + QoI_deriv_2 = np.zeros([numSamples,numKL]) + + # For each sample solve the PDE + f = Constant(-1.0) # forcing of Poisson + + for i in range(0, numSamples): + + print "Sample point number: %g" % i + + # create a temp array to store logPerm as sum of KL expansions + # logPerm is log permeability + logPerm = np.zeros((mesh.num_vertices()), dtype=float) + for kl in range(0, numKL): + logPerm += xi_k[i, kl] * \ + sqrt(eigen_val[kl]) * eigen_func[kl].vector().array() + + # permiability is the exponential of log permeability logPerm + perm_k_array = 0.1 + np.exp(logPerm) + # print "Mean value of the random field is: %g" % perm_k_array.mean() + + ## use dof_to_vertex map to map values to the function space + perm_k.vector()[:] = perm_k_array + + # solve Poisson with this random field using FEM + u = solvePoissonRandomField(perm_k, mesh, 1, f, bcs) + + # Compute QoI + QoI_samples[i, 0] = assemble(u * Chi_1 * dx) + QoI_samples[i, 1] = assemble(u * Chi_2 * dx) + + + return QoI_samples \ No newline at end of file diff --git a/examples/FEniCS/poissonRandField.py b/examples/FEniCS/poissonRandField.py new file mode 100644 index 00000000..626474a6 --- /dev/null +++ b/examples/FEniCS/poissonRandField.py @@ -0,0 +1,16 @@ +from dolfin import* + +def solvePoissonRandomField(rand_field,mesh,poly_order,f,bcs): + """ + Solves the poisson equation with a random field : + (\grad \dot (rand_field \grad(u)) = -f) + """ + # create the function space + V = FunctionSpace(mesh, "CG", poly_order) + u = TrialFunction(V) + v = TestFunction(V) + L = f*v*dx + a = inner(rand_field*nabla_grad(u),nabla_grad(v))*dx + u = Function(V) + solve(a == L,u,bcs) + return u diff --git a/examples/FEniCS/projectKL.py b/examples/FEniCS/projectKL.py new file mode 100644 index 00000000..3d57c88c --- /dev/null +++ b/examples/FEniCS/projectKL.py @@ -0,0 +1,168 @@ +from dolfin import * +import numpy as np +import petsc4py +from petsc4py import PETSc +from slepc4py import SLEPc +from meshDS import* +# initialize petsc +petsc4py.init() + +class projectKL(object): + + """Docstring for projectKL. """ + + def __init__(self,mesh): + """TODO: to be defined1. """ + # create meshDS obect + self._mesh = mesh + self.domain = meshDS(mesh) + self.c_volume_array,self.c_centroid_array = self.domain.getElemVCArray() + self.node_to_elem = self.domain.getNodesToElem() + self.flag = False + def getCovMat(self, cov_expr): + """TODO: Docstring for getCovMat. + + :cov_expr: Expression (dolfin) as a function of + :returns: covariance PETSC matrix cov_mat + + """ + # store the expression + self.expr = cov_expr + # create a PETSC matrix for cov_mat + cov_mat = PETSc.Mat().create() + cov_mat.setType('aij') + cov_mat.setSizes(self.domain.getNodes(),self.domain.getNodes()) + cov_mat.setUp() + + cov_ij = np.empty((1),dtype=float) # scalar valued function is evaluated in this variable + xycor = np.empty((4),dtype=float) # the points to evalute the expression + + print '---------------------------' + print '---------------------------' + print ' Building Covariance Matrix' + print '---------------------------' + print '---------------------------' + # Loop through global nodes and build the matrix for i < j because of symmetric nature. + for node_i in range(0,self.domain.getNodes()): + # global node node_i + for node_j in range(node_i,self.domain.getNodes()): + # global node node_j + temp_cov_ij = 0 + for elem_i in self.node_to_elem[node_i]: + # elem_i : element attached to node_i + # x1 : x co-ordinate of the centroid of element elem_i + x1 = self.c_centroid_array[elem_i].x() + # y1 : x co-ordinate of the centroid of element elem_i + y1 = self.c_centroid_array[elem_i].y() + for elem_j in self.node_to_elem[node_j]: + # elem_j : element attached to node_j + # x2 : x co-ordinate for the centroid of element elem_j + x2 = self.c_centroid_array[elem_j].x() + # y2 : y co-ordinate for the centroid of element elem_j + y2 = self.c_centroid_array[elem_j].y() + xycor[0] = x1 + xycor[1] = x2 + xycor[2] = y1 + xycor[3] = y2 + # evaluate the expression + cov_expr.eval(cov_ij,xycor) + if cov_ij[0] > 0: + temp_cov_ij += (1.0/3)*(1.0/3)*cov_ij[0]*self.c_volume_array[elem_i]* \ + self.c_volume_array[elem_j] + cov_mat.setValue(node_i,node_j,temp_cov_ij) + cov_mat.setValue(node_j,node_i,temp_cov_ij) + cov_mat.assemblyBegin() + cov_mat.assemblyEnd() + print '---------------------------' + print '---------------------------' + print ' Finished Covariance Matrix' + print '---------------------------' + print '---------------------------' + + return cov_mat + + def _getBMat(self): + """TODO: Docstring for getBmat. We are solving for CX = BX where C is the covariance matrix + and B is just a mass matrix. Here we assemble B. This is a private function. DONT call this + unless debuging. + + :returns: PETScMatrix B + """ + + # B matrix is just a mass matrix, can be easily assembled through fenics + # however, the ordering in fenics is not the mesh ordering. so we build a temp matrix + # then use the vertex to dof map to get the right ordering interms of our mesh nodes + V = FunctionSpace(self._mesh, "CG", 1) + # Define basis and bilinear form + u = TrialFunction(V) + v = TestFunction(V) + a = u*v*dx + B_temp = assemble(a) + + B = PETSc.Mat().create() + B.setType('aij') + B.setSizes(self.domain.getNodes(),self.domain.getNodes()) + B.setUp() + + B_ij = B_temp.array() + + v_to_d_map = vertex_to_dof_map(V) + + print '---------------------------' + print '---------------------------' + print ' Building Mass Matrix ' + print '---------------------------' + print '---------------------------' + for node_i in range(0, self.domain.getNodes()): + for node_j in range(node_i, self.domain.getNodes()): + B_ij_nodes = B_ij[v_to_d_map[node_i],v_to_d_map[node_j]] + if B_ij_nodes > 0: + B.setValue(node_i,node_j,B_ij_nodes) + B.setValue(node_j,node_i,B_ij_nodes) + + B.assemblyBegin() + B.assemblyEnd() + print '---------------------------' + print '---------------------------' + print ' Finished Mass Matrix ' + print '---------------------------' + print '---------------------------' + return B + + def projectCovToMesh(self,num_kl,cov_expr): + """TODO: Docstring for projectCovToMesh. Solves CX = BX where C is the covariance matrix + :num_kl : number of kl exapansion terms needed + :returns: TODO + + """ + # turn the flag to true + self.flag = True + # get C,B matrices + C = PETScMatrix(self.getCovMat(cov_expr)) + B = PETScMatrix(self._getBMat()) + # Solve the generalized eigenvalue problem + eigensolver = SLEPcEigenSolver(C,B) + eigensolver.solve(num_kl) + # Get the number of eigen values that converged. + nconv = eigensolver.get_number_converged() + + # Get N eigenpairs where N is the number of KL expansion and check if N < nconv otherwise you had + # really bad matrix + + # create numpy array of vectors and eigenvalues + self.eigen_funcs = np.empty((num_kl),dtype=object) + self.eigen_vals = np.empty((num_kl),dtype=float) + + # store the eigenvalues and eigen functions + V = FunctionSpace(self._mesh, "CG", 1) + for eigen_pairs in range(0,num_kl): + lambda_r, lambda_c, x_real, x_complex = eigensolver.get_eigenpair(eigen_pairs) + self.eigen_funcs[eigen_pairs] = Function(V) + # use dof_to_vertex map to map values to the function space + self.eigen_funcs[eigen_pairs].vector()[:] = x_real[dof_to_vertex_map(V)]#*np.sqrt(lambda_r) + # divide by norm to make the unit norm again + self.eigen_funcs[eigen_pairs].vector()[:] = self.eigen_funcs[eigen_pairs].vector()[:] / \ + norm(self.eigen_funcs[eigen_pairs]) + self.eigen_vals[eigen_pairs] = lambda_r + + From 9f2c77d7d6ab63c02360fb395da6a93572ed85d4 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Mon, 23 May 2016 20:33:07 -0600 Subject: [PATCH 16/20] Updated contaminant example and plot param_ref --- bet/postProcess/postTools.py | 4 +- examples/FEniCS/BET_script.py | 5 +- examples/contaminantTransport/contaminant.py | 103 ++++++++++-------- .../linearMap/linearMapUniformSampling.py | 4 +- .../nonlinearMapUniformSampling.py | 4 +- 5 files changed, 66 insertions(+), 54 deletions(-) diff --git a/bet/postProcess/postTools.py b/bet/postProcess/postTools.py index 07a71379..b8510054 100644 --- a/bet/postProcess/postTools.py +++ b/bet/postProcess/postTools.py @@ -72,7 +72,7 @@ def sort_by_rho(sample_set): sample_set_out._input_sample_set.set_values(samples) sample_set_out._input_sample_set.set_probabilities(P_samples) sample_set_out._input_sample_set.set_volumes(lam_vol) - sample_set_out._output_sample_set.set_data(data) + sample_set_out._output_sample_set.set_values(data) else: sample_set_out = sample.sample_set(sample_set.get_dim()) sample_set_out.set_values(samples) @@ -161,7 +161,7 @@ def sample_prob(percentile, sample_set, sort=True, descending=False): sample_set_out._input_sample_set.set_values(samples) sample_set_out._input_sample_set.set_probabilities(P_samples) sample_set_out._input_sample_set.set_volumes(lam_vol) - sample_set_out._output_sample_set.set_data(data) + sample_set_out._output_sample_set.set_values(data) else: sample_set_out = sample.sample_set(sample_set.get_dim()) sample_set_out.set_values(samples) diff --git a/examples/FEniCS/BET_script.py b/examples/FEniCS/BET_script.py index 86464998..4df77161 100644 --- a/examples/FEniCS/BET_script.py +++ b/examples/FEniCS/BET_script.py @@ -147,7 +147,8 @@ # plot 2d marginals probs plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename="FEniCS", - file_extension=".eps", plot_surface=False) + lam_ref=param_ref[0,:], file_extension=".eps", + plot_surface=False) # calculate 1d marginal probs (bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, @@ -156,7 +157,7 @@ marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.5) # plot 2d marginal probs plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename="FEniCS", - file_extension=".eps") + lam_ref=param_ref[0,:], file_extension=".eps") diff --git a/examples/contaminantTransport/contaminant.py b/examples/contaminantTransport/contaminant.py index f59f0d26..37329520 100644 --- a/examples/contaminantTransport/contaminant.py +++ b/examples/contaminantTransport/contaminant.py @@ -26,92 +26,103 @@ import bet.postProcess.plotP as plotP import bet.postProcess.plotDomains as plotD import bet.postProcess.postTools as postTools - +import bet.sample as samp # Labels and descriptions of the uncertain parameters labels = ['Source $y$ coordinate [L]', 'Source $x$ coordinate [L]', 'Dispersivity x [L]', 'Flow Angle [degrees]', 'Contaminant flux [M/T]'] # Load data from files -lam_domain = np.loadtxt("files/lam_domain.txt.gz") #parameter domain -ref_lam = np.loadtxt("files/lam_ref.txt.gz") #reference parameter set -Q_ref = np.loadtxt("files/Q_ref.txt.gz") #reference QoI set -samples = np.loadtxt("files/samples.txt.gz") # uniform samples in parameter domain -dataf = np.loadtxt("files/data.txt.gz") # data from model - -QoI_indices=[0,1,2,3] # Indices for output data with which you want to invert -bin_ratio = 0.25 #ratio of length of data region to invert - -data = dataf[:,QoI_indices] -Q_ref=Q_ref[QoI_indices] - -dmax = data.max(axis=0) -dmin = data.min(axis=0) -dscale = bin_ratio*(dmax-dmin) -Qmax = Q_ref + 0.5*dscale -Qmin = Q_ref -0.5*dscale -def rho_D(x): - return np.all(np.logical_and(np.greater(x,Qmin), np.less(x,Qmax)),axis=1) +# First obtain info on the parameter domain +parameter_domain = np.loadtxt("files/lam_domain.txt.gz") #parameter domain +parameter_dim = parameter_domain.shape[0] +# Create input sample set +input_samples = samp.sample_set(parameter_dim) +input_samples.set_domain(parameter_domain) +input_samples.set_values(np.loadtxt("files/samples.txt.gz")) +input_samples.estimate_volume_mc() # Use standard MC estimate of volumes +# Choose which QoI to use and create output sample set +QoI_indices_observe = np.array([0,1,2,3]) +output_samples = samp.sample_set(QoI_indices_observe.size) +output_samples.set_values(np.loadtxt("files/data.txt.gz")[:,QoI_indices_observe]) + +# Create discretization object +my_discretization = samp.discretization(input_sample_set=input_samples, + output_sample_set=output_samples) + +# Load the reference parameter and QoI values +param_ref = np.loadtxt("files/lam_ref.txt.gz") #reference parameter set +Q_ref = np.loadtxt("files/Q_ref.txt.gz")[QoI_indices_observe] #reference QoI set # Plot the data domain -plotD.show_data(data, Q_ref = Q_ref, rho_D=rho_D, showdim=2) +plotD.show_data(my_discretization, Q_ref = Q_ref, showdim=2) # Whether or not to use deterministic description of simple function approximation of # ouput probability deterministic_discretize_D = True if deterministic_discretize_D == True: - (d_distr_prob, d_distr_samples, d_Tree) = simpleFunP.uniform_hyperrectangle(data=data, - Q_ref=Q_ref, - bin_ratio=bin_ratio, - center_pts_per_edge = 1) + simpleFunP.regular_partition_uniform_distribution_rectangle_scaled(data_set=my_discretization, + Q_ref=Q_ref, + rect_scale=0.25, + center_pts_per_edge = 1) else: - (d_distr_prob, d_distr_samples, d_Tree) = simpleFunP.unif_unif(data=data, - Q_ref=Q_ref, - M=50, - bin_ratio=bin_ratio, - num_d_emulate=1E5) + simpleFunP.uniform_partition_uniform_distribution_rectangle_scaled(data_set=my_discretization, + Q_ref=Q_ref, + rect_scale=0.25, + M=50, + num_d_emulate=1E5) # calculate probablities making Monte Carlo assumption -(P, lam_vol, io_ptr) = calculateP.prob(samples=samples, - data=data, - rho_D_M=d_distr_prob, - d_distr_samples=d_distr_samples) +calculateP.prob(my_discretization) # calculate 2D marginal probabilities -(bins, marginals2D) = plotP.calculate_2D_marginal_probs(P_samples = P, samples = samples, lam_domain = lam_domain, nbins = 10) +(bins, marginals2D) = plotP.calculate_2D_marginal_probs(my_discretization, nbins = 10) # smooth 2D marginal probabilites for plotting (optional) -marginals2D = plotP.smooth_marginals_2D(marginals2D,bins, sigma=1.0) +marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=1.0) # plot 2D marginal probabilities -plotP.plot_2D_marginal_probs(marginals2D, bins, lam_domain, filename = "contaminant_map", +plotP.plot_2D_marginal_probs(marginals2D, bins, my_discretization, filename = "contaminant_map", plot_surface=False, - lam_ref = ref_lam, + lam_ref = param_ref, lambda_label=labels, interactive=False) # calculate 1d marginal probs -(bins, marginals1D) = plotP.calculate_1D_marginal_probs(P_samples = P, samples = samples, lam_domain = lam_domain, nbins = 20) +(bins, marginals1D) = plotP.calculate_1D_marginal_probs(my_discretization, nbins = 20) # smooth 1d marginal probs (optional) marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=1.0) # plot 1d marginal probs -plotP.plot_1D_marginal_probs(marginals1D, bins, lam_domain, filename = "contaminant_map", interactive=False, lam_ref=ref_lam, lambda_label=labels) +plotP.plot_1D_marginal_probs(marginals1D, bins, my_discretization, + filename = "contaminant_map", + interactive=False, + lam_ref=param_ref, + lambda_label=labels) percentile = 1.0 # Sort samples by highest probability density and sample highest percentile percent samples -(num_samples, P_high, samples_high, lam_vol_high, data_high)= postTools.sample_highest_prob(top_percentile=percentile, P_samples=P, samples=samples, lam_vol=lam_vol,data = data,sort=True) +(num_samples, my_discretization_highP, indices)= postTools.sample_highest_prob( + percentile, my_discretization, sort=True) # print the number of samples that make up the highest percentile percent samples and # ratio of the volume of the parameter domain they take up -print (num_samples, np.sum(lam_vol_high)) +print (num_samples, np.sum(my_discretization_highP._input_sample_set.get_volumes())) + +# Choose unused QoI as prediction QoI and propagate measure onto predicted QoI data space +QoI_indices_predict = np.array([7]) +output_samples_predict = samp.sample_set(QoI_indices_predict.size) +output_samples_predict.set_values(np.loadtxt("files/data.txt.gz")[:,QoI_indices_predict]) +output_samples_predict.set_probabilities(input_samples.get_probabilities()) -# Propogate the probability measure through a different QoI map -(_, P_pred, _, _ , data_pred)= postTools.sample_highest_prob(top_percentile=percentile, P_samples=P, samples=samples, lam_vol=lam_vol,data = dataf[:,7],sort=True) +# Determine range of predictions and store as domain for plotting purposes +output_samples_predict.set_domain(output_samples_predict.get_bounding_box()) # Plot 1D pdf of predicted QoI # calculate 1d marginal probs -(bins_pred, marginals1D_pred) = plotP.calculate_1D_marginal_probs(P_samples = P_pred, samples = data_pred, lam_domain = np.array([[np.min(data_pred),np.max(data_pred)]]), nbins = 20) +(bins_pred, marginals1D_pred) = plotP.calculate_1D_marginal_probs(output_samples_predict, + nbins = 20) # plot 1d pdf -plotP.plot_1D_marginal_probs(marginals1D_pred, bins_pred, lam_domain= np.array([[np.min(data_pred),np.max(data_pred)]]), filename = "contaminant_prediction", interactive=False) +plotP.plot_1D_marginal_probs(marginals1D_pred, bins_pred, output_samples_predict, + filename = "contaminant_prediction", interactive=False) diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index 77a766e2..260ff77b 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -151,7 +151,7 @@ # plot 2d marginals probs plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "linearMap", - file_extension = ".eps", plot_surface=False) + lam_ref=param_ref, file_extension = ".eps", plot_surface=False) # calculate 1d marginal probs (bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, @@ -160,7 +160,7 @@ marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.2) # plot 2d marginal probs plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "linearMap", - file_extension = ".eps") + lam_ref=param_ref, file_extension = ".eps") diff --git a/examples/nonlinearMap/nonlinearMapUniformSampling.py b/examples/nonlinearMap/nonlinearMapUniformSampling.py index 70d4857a..97cb3f37 100644 --- a/examples/nonlinearMap/nonlinearMapUniformSampling.py +++ b/examples/nonlinearMap/nonlinearMapUniformSampling.py @@ -152,7 +152,7 @@ # plot 2d marginals probs plotP.plot_2D_marginal_probs(marginals2D, bins, input_samples, filename = "nomlinearMap", - file_extension = ".eps", plot_surface=False) + lam_ref = param_ref, file_extension = ".eps", plot_surface=False) # calculate 1d marginal probs (bins, marginals1D) = plotP.calculate_1D_marginal_probs(input_samples, @@ -161,7 +161,7 @@ marginals1D = plotP.smooth_marginals_1D(marginals1D, bins, sigma=0.5) # plot 2d marginal probs plotP.plot_1D_marginal_probs(marginals1D, bins, input_samples, filename = "nonlinearMap", - file_extension = ".eps") + lam_ref = param_ref, file_extension = ".eps") From bd35758181c1fe542c6cda97028d45dd7643d501 Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Tue, 24 May 2016 14:53:29 -0600 Subject: [PATCH 17/20] Tests for all basicSampling done --- bet/sampling/basicSampling.py | 50 +++--- test/test_sampling/test_basicSampling.py | 197 +++++++++++++++++++++++ 2 files changed, 225 insertions(+), 22 deletions(-) diff --git a/bet/sampling/basicSampling.py b/bet/sampling/basicSampling.py index ac1ebdeb..691edc7e 100644 --- a/bet/sampling/basicSampling.py +++ b/bet/sampling/basicSampling.py @@ -169,9 +169,9 @@ def random_sample_set_domain(self, sample_type, input_domain, :param string criterion: latin hypercube criterion see `PyDOE `_ - :rtype: :class:`~bet.sample.discretization` - :returns: :class:`~bet.sample.discretization` object which contains - input and output of ``num_samples`` + :rtype: :class:`~bet.sample.sample_set` + :returns: :class:`~bet.sample.sample_Set` object which contains + input ``num_samples`` """ # Create N samples @@ -203,9 +203,9 @@ def random_sample_set_dimension(self, sample_type, input_dim, :param string criterion: latin hypercube criterion see `PyDOE `_ - :rtype: :class:`~bet.sample.discretization` - :returns: :class:`~bet.sample.discretization` object which contains - input and output of ``num_samples`` + :rtype: :class:`~bet.sample.sample_set` + :returns: :class:`~bet.sample.sample_Set` object which contains + input ``num_samples`` """ # Create N samples @@ -227,9 +227,9 @@ def regular_sample_set(self, input_sample_set, num_samples_per_dim=1): :type num_samples_per_dim: :class: `~numpy.ndarray` of dimension (input_sample_set._dim,) - :rtype: :class:`~bet.sample.discretization` - :returns: :class:`~bet.sample.discretization` object which contains - input and output of ``num_samples`` + :rtype: :class:`~bet.sample.sample_set` + :returns: :class:`~bet.sample.sample_Set` object which contains + input ``num_samples`` """ # Create N samples @@ -237,7 +237,7 @@ def regular_sample_set(self, input_sample_set, num_samples_per_dim=1): if not isinstance(num_samples_per_dim, collections.Iterable): num_samples_per_dim = num_samples_per_dim * np.ones((dim,)) - if np.any(np.less(num_samples_per_dim, 0)): + if np.any(np.less_equal(num_samples_per_dim, 0)): print 'Warning: num_smaples_per_dim must be greater than 0' self.num_samples = np.product(num_samples_per_dim) @@ -258,11 +258,17 @@ def regular_sample_set(self, input_sample_set, num_samples_per_dim=1): input_domain[i,0], input_domain[i,1], num_samples_per_dim[i]+2))[1:num_samples_per_dim[i]+1] - arrays_samples_dimension = np.meshgrid( - *[vec_samples_dimension[i] for i in np.arange(0, dim)], indexing='ij') + if np.equal(dim, 1): + arrays_samples_dimension = np.array([vec_samples_dimension]) + else: + arrays_samples_dimension = np.meshgrid( + *[vec_samples_dimension[i] for i in np.arange(0, dim)], indexing='ij') - for i in np.arange(0, dim): - input_values[:,i:i+1] = np.vstack(arrays_samples_dimension[i].flat[:]) + if np.equal(dim, 1): + input_values = arrays_samples_dimension.transpose() + else: + for i in np.arange(0, dim): + input_values[:,i:i+1] = np.vstack(arrays_samples_dimension[i].flat[:]) input_sample_set.set_values(input_values) @@ -281,9 +287,9 @@ def regular_sample_set_domain(self, input_domain, num_samples_per_dim=1): :type num_samples_per_dim: :class: `~numpy.ndarray` of dimension (input_sample_set._dim,) - :rtype: :class:`~bet.sample.discretization` - :returns: :class:`~bet.sample.discretization` object which contains - input and output of ``num_samples`` + :rtype: :class:`~bet.sample.sample_set` + :returns: :class:`~bet.sample.sample_Set` object which contains + input ``num_samples`` """ # Create N samples @@ -292,7 +298,7 @@ def regular_sample_set_domain(self, input_domain, num_samples_per_dim=1): return self.regular_sample_set(input_sample_set, num_samples_per_dim) - def regular_sample_set_dimension(self, input_dim, num_samples_per_dimension): + def regular_sample_set_dimension(self, input_dim, num_samples_per_dim=1): """ Sampling algorithm for generating a regular grid of samples taken on a unit hypercube of dimension input_dim @@ -302,15 +308,15 @@ def regular_sample_set_dimension(self, input_dim, num_samples_per_dimension): :type num_samples_per_dim: :class: `~numpy.ndarray` of dimension (input_sample_set._dim,) - :rtype: :class:`~bet.sample.discretization` - :returns: :class:`~bet.sample.discretization` object which contains - input and output of ``num_samples`` + :rtype: :class:`~bet.sample.sample_set` + :returns: :class:`~bet.sample.sample_Set` object which contains + input ``num_samples`` """ # Create N samples input_sample_set = sample.sample_set(input_dim) - return self.regular_sample_set(input_sample_set, num_samples_per_dimension) + return self.regular_sample_set(input_sample_set, num_samples_per_dim) def compute_QoI_and_create_discretization(self, input_sample_set, savefile=None, parallel=False): diff --git a/test/test_sampling/test_basicSampling.py b/test/test_sampling/test_basicSampling.py index 5da856d3..abb90623 100644 --- a/test/test_sampling/test_basicSampling.py +++ b/test/test_sampling/test_basicSampling.py @@ -15,6 +15,7 @@ import bet.sample from bet.sample import sample_set from bet.sample import discretization as disc +import collections local_path = os.path.join(".") @@ -358,6 +359,158 @@ def verify_random_sample_set(sampler, sample_type, input_sample_set, nptest.assert_array_equal(test_sample_set._values, my_sample_set._values) +def verify_regular_sample_set(sampler, input_sample_set, + num_samples_per_dim): + + test_sample_set = input_sample_set + dim = input_sample_set.get_dim() + # recreate the samples + if num_samples_per_dim is None: + num_samples_per_dim = 5 + + if not isinstance(num_samples_per_dim, collections.Iterable): + num_samples_per_dim = num_samples_per_dim * np.ones((dim,), dtype='int') + + sampler.num_samples = np.product(num_samples_per_dim) + + test_domain = test_sample_set.get_domain() + if test_domain is None: + test_domain = np.repeat([[0, 1]], test_sample_set.get_dim(), axis=0) + + test_values = np.zeros((sampler.num_samples, test_sample_set.get_dim())) + + vec_samples_dimension = np.empty((dim), dtype=object) + for i in np.arange(0, dim): + vec_samples_dimension[i] = list(np.linspace( + test_domain[i, 0], test_domain[i, 1], + num_samples_per_dim[i] + 2))[1:num_samples_per_dim[i] + 1] + + if np.equal(dim, 1): + arrays_samples_dimension = np.array([vec_samples_dimension]) + else: + arrays_samples_dimension = np.meshgrid( + *[vec_samples_dimension[i] for i in np.arange(0, dim)], indexing='ij') + + if np.equal(dim, 1): + test_values = arrays_samples_dimension.transpose() + else: + for i in np.arange(0, dim): + test_values[:, i:i + 1] = np.vstack(arrays_samples_dimension[i].flat[:]) + + test_sample_set.set_values(test_values) + + # create the sample set from sampler + my_sample_set = sampler.regular_sample_set(input_sample_set, + num_samples_per_dim=num_samples_per_dim) + + # compare the samples + nptest.assert_array_equal(test_sample_set._values, + my_sample_set._values) + +def verify_regular_sample_set_domain(sampler, input_domain, + num_samples_per_dim): + + input_sample_set = sample_set(input_domain.shape[0]) + input_sample_set.set_domain(input_domain) + + test_sample_set = input_sample_set + dim = input_sample_set.get_dim() + # recreate the samples + if num_samples_per_dim is None: + num_samples_per_dim = 5 + + if not isinstance(num_samples_per_dim, collections.Iterable): + num_samples_per_dim = num_samples_per_dim * np.ones((dim,), dtype='int') + + sampler.num_samples = np.product(num_samples_per_dim) + + test_domain = test_sample_set.get_domain() + if test_domain is None: + test_domain = np.repeat([[0, 1]], test_sample_set.get_dim(), axis=0) + + test_values = np.zeros((sampler.num_samples, test_sample_set.get_dim())) + + vec_samples_dimension = np.empty((dim), dtype=object) + for i in np.arange(0, dim): + vec_samples_dimension[i] = list(np.linspace( + test_domain[i, 0], test_domain[i, 1], + num_samples_per_dim[i] + 2))[1:num_samples_per_dim[i] + 1] + + if np.equal(dim, 1): + arrays_samples_dimension = np.array([vec_samples_dimension]) + else: + arrays_samples_dimension = np.meshgrid( + *[vec_samples_dimension[i] for i in np.arange(0, dim)], indexing='ij') + + if np.equal(dim, 1): + test_values = arrays_samples_dimension.transpose() + else: + for i in np.arange(0, dim): + test_values[:, i:i + 1] = np.vstack(arrays_samples_dimension[i].flat[:]) + + test_sample_set.set_values(test_values) + + # create the sample set from sampler + my_sample_set = sampler.regular_sample_set_domain(input_domain, + num_samples_per_dim=num_samples_per_dim) + + # compare the samples + nptest.assert_array_equal(test_sample_set._values, + my_sample_set._values) + +def verify_regular_sample_set_dimension(sampler, input_dim, + num_samples_per_dim): + + input_domain = np.repeat([[0, 1]], input_dim, axis=0) + input_sample_set = sample_set(input_dim) + input_sample_set.set_domain(input_domain) + + test_sample_set = input_sample_set + dim = input_dim + # recreate the samples + if num_samples_per_dim is None: + num_samples_per_dim = 5 + + if not isinstance(num_samples_per_dim, collections.Iterable): + num_samples_per_dim = num_samples_per_dim * np.ones((dim,), dtype='int') + + sampler.num_samples = np.product(num_samples_per_dim) + + test_domain = test_sample_set.get_domain() + if test_domain is None: + test_domain = np.repeat([[0, 1]], test_sample_set.get_dim(), axis=0) + + test_values = np.zeros((sampler.num_samples, test_sample_set.get_dim())) + + vec_samples_dimension = np.empty((dim), dtype=object) + for i in np.arange(0, dim): + vec_samples_dimension[i] = list(np.linspace( + test_domain[i, 0], test_domain[i, 1], + num_samples_per_dim[i] + 2))[1:num_samples_per_dim[i] + 1] + + if np.equal(dim, 1): + arrays_samples_dimension = np.array([vec_samples_dimension]) + else: + arrays_samples_dimension = np.meshgrid( + *[vec_samples_dimension[i] for i in np.arange(0, dim)], indexing='ij') + + if np.equal(dim, 1): + test_values = arrays_samples_dimension.transpose() + else: + for i in np.arange(0, dim): + test_values[:, i:i + 1] = np.vstack(arrays_samples_dimension[i].flat[:]) + + test_sample_set.set_values(test_values) + + # create the sample set from sampler + my_sample_set = sampler.regular_sample_set_dimension(input_dim, + num_samples_per_dim=num_samples_per_dim) + + # compare the samples + nptest.assert_array_equal(test_sample_set._values, + my_sample_set._values) + + class Test_basic_sampler(unittest.TestCase): """ Test :class:`bet.sampling.basicSampling.sampler`. @@ -508,6 +661,50 @@ def test_random_sample_set_dim(self): verify_random_sample_set_dimension(sampler, sample_type, input_dim, num_samples) + def test_regular_sample_set(self): + """ + Test :meth:`bet.sampling.basicSampling.sampler.regular_sample_set` + for six different sample sets + """ + input_sample_set_list = [self.input_sample_set1, + self.input_sample_set2, + self.input_sample_set4, + self.input_sample_set5] + + test_list = zip(self.samplers, input_sample_set_list) + + for sampler, input_sample_set in test_list: + for num_samples_per_dim in [None, 10]: + verify_regular_sample_set(sampler, input_sample_set, num_samples_per_dim) + + def test_regular_sample_set_domain(self): + """ + Test :meth:`bet.sampling.basicSampling.sampler.regular_sample_set_domain` + for six different sample sets + """ + input_domain_list= [self.input_domain1, + self.input_domain3] + + test_list = zip(self.samplers, input_domain_list) + + for sampler, input_domain in test_list: + for num_samples_per_dim in [None, 10]: + verify_regular_sample_set_domain(sampler, input_domain, num_samples_per_dim) + + def test_regular_sample_set_dimension(self): + """ + Test :meth:`bet.sampling.basicSampling.sampler.regular_sample_set_dimension` + for six different sample sets + """ + input_dimension_list = [self.input_dim1, + self.input_dim2] + + test_list = zip(self.samplers, input_dimension_list) + + for sampler, input_dim in test_list: + for num_samples_per_dim in [None, 10]: + verify_regular_sample_set_dimension(sampler, input_dim, num_samples_per_dim) + def test_create_random_discretization(self): """ Test :meth:`bet.sampling.basicSampling.sampler.create_random_discretization` From 92b181b0c1ce608b5b0921c9ef1fdfed550fa52c Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Tue, 24 May 2016 16:14:00 -0600 Subject: [PATCH 18/20] New tests for new simpleFunP functions --- bet/calculateP/simpleFunP.py | 30 ++--- examples/validationExample/linearMap.py | 14 +-- test/test_calculateP/test_simpleFunP.py | 148 ++++++++++++++++++++++++ 3 files changed, 172 insertions(+), 20 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index b0ff707b..d0f84fae 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -636,6 +636,10 @@ def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate r'''Create M samples defining M bins in D used to define :math:`\rho_{\mathcal{D},M}` rho_D is assumed to be a multi-variate normal distribution with mean Q_ref and standard deviation std.''' + if not isinstance(Q_ref, collections.Iterable): + Q_ref = np.array([Q_ref]) + if not isinstance(std, collections.Iterable): + std = np.array([std]) bin_size = 4.0 * std d_distr_samples = np.zeros((M, len(Q_ref))) @@ -684,22 +688,22 @@ def uniform_partition_normal_distribution(data_set, Q_ref, std, M, num_d_emulate data_set._output_probability_set = s_set return s_set -def user_discretization_user_distribution(data_set, data_discretization_set, +def user_partition_user_distribution(data_set, data_partition_set, data_distribution_set): r""" Creates a user defined simple function approximation of a user defined distribution. The simple function discretization is - specified in the ``data_discretization_set``, and the set of i.i.d. + specified in the ``data_partition_set``, and the set of i.i.d. samples from the distribution is specified in the ``data_distribution_set``. :param data_set: Sample set that the probability measure is defined for. :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` - :param data_discretization_set: Sample set defining the discretization + :param data_partition_set: Sample set defining the discretization of the data space into Voronoi cells for which a simple function is defined upon. - :type data_discretization_set: :class:`~bet.sample.discretization` or + :type data_partition_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray` :param data_distribution_set: Sample set containing the i.i.d. samples from the distribution on the data space that are binned within the @@ -726,17 +730,17 @@ def user_discretization_user_distribution(data_set, data_discretization_set, msg += "bet.sample.discretization or np.ndarray" raise wrong_argument_type(msg) - if isinstance(data_discretization_set, samp.sample_set_base): - num_samples_discretize_D = data_discretization_set.check_num() - simpleFun_set = data_discretization_set.copy() + if isinstance(data_partition_set, samp.sample_set_base): + num_samples_discretize_D = data_partition_set.check_num() + simpleFun_set = data_partition_set.copy() dim_simpleFun = simpleFun_set._dim - elif isinstance(data_discretization_set, samp.discretization): - num_samples_discretize_D = data_discretization_set.check_nums() - simpleFun_set = data_discretization_set._output_sample_set.copy() + elif isinstance(data_partition_set, samp.discretization): + num_samples_discretize_D = data_partition_set.check_nums() + simpleFun_set = data_partition_set._output_sample_set.copy() dim_simpleFun = simpleFun_set._dim - elif isinstance(data_discretization_set, np.ndarray): - num_samples_discretize_D = data_discretization_set.shape[0] - dim_simpleFun = data_discretization_set.shape[1] + elif isinstance(data_partition_set, np.ndarray): + num_samples_discretize_D = data_partition_set.shape[0] + dim_simpleFun = data_partition_set.shape[1] simpleFun_set = samp.sample_set(dim=dim_simpleFun) simpleFun_set.set_values(values) else: diff --git a/examples/validationExample/linearMap.py b/examples/validationExample/linearMap.py index f5f19acf..389b9c92 100644 --- a/examples/validationExample/linearMap.py +++ b/examples/validationExample/linearMap.py @@ -100,14 +100,14 @@ num_samples_discretize_D = 1 num_iid_samples = 1E5 -SimpleFunction_set = samp.sample_set(2) +Partition_set = samp.sample_set(2) Monte_Carlo_set = samp.sample_set(2) -SimpleFunction_set.set_domain(np.repeat([[0.0, 1.0]], 2, axis=0)) +Partition_set.set_domain(np.repeat([[0.0, 1.0]], 2, axis=0)) Monte_Carlo_set.set_domain(np.repeat([[0.0, 1.0]], 2, axis=0)) -SimpleFunction_discretization = sampler.create_random_discretization('random', - SimpleFunction_set, +Partition_discretization = sampler.create_random_discretization('random', + Partition_set, num_samples=num_samples_discretize_D) Monte_Carlo_discretization = sampler.create_random_discretization('random', @@ -115,9 +115,9 @@ num_samples=num_iid_samples) # Compute the simple function approximation to the distribution on the data space -simpleFunP.user_discretization_user_distribution(my_discretization, - SimpleFunction_discretization, - Monte_Carlo_discretization) +simpleFunP.user_partition_user_distribution(my_discretization, + Partition_discretization, + Monte_Carlo_discretization) # Calculate probabilities calculateP.prob(my_discretization) diff --git a/test/test_calculateP/test_simpleFunP.py b/test/test_calculateP/test_simpleFunP.py index e49645c9..86e1fa7c 100644 --- a/test/test_calculateP/test_simpleFunP.py +++ b/test/test_calculateP/test_simpleFunP.py @@ -304,6 +304,83 @@ def setUp(self): super(test_normal_partition_normal_distribution_3D, self).setUp() +class uniform_partition_normal_distribution(prob): + """ + Set up :meth:`bet.calculateP.simpleFunP.uniform_partition_normal_distribution` on data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + if type(self.Q_ref) != np.array and type(self.Q_ref) != np.ndarray: + std = 1.0 + else: + std = np.ones(self.Q_ref.shape) + self.data_prob = sFun.uniform_partition_normal_distribution(None, self.Q_ref, std=std, M=67, num_d_emulate=1E3) + self.d_distr_samples = self.data_prob.get_values() + self.rho_D_M = self.data_prob.get_probabilities() + + def test_M(self): + """ + Test that the right number of d_distr_samples are used to create + rho_D_M. + """ + assert len(self.rho_D_M) == 67 + + +class test_uniform_partition_normal_distribution_01D(data_01D, uniform_partition_normal_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_normal_distribution` on 01D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_normal_distribution_01D, self).createData() + super(test_uniform_partition_normal_distribution_01D, self).setUp() + + +class test_uniform_partition_normal_distribution_1D(data_1D, uniform_partition_normal_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_normal_distribution` on 1D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_normal_distribution_1D, self).createData() + super(test_uniform_partition_normal_distribution_1D, self).setUp() + + +class test_uniform_partition_normal_distribution_2D(data_2D, uniform_partition_normal_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_normal_distribution` on 2D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_normal_distribution_2D, self).createData() + super(test_uniform_partition_normal_distribution_2D, self).setUp() + + +class test_uniform_partition_normal_distribution_3D(data_3D, uniform_partition_normal_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.uniform_partition_normal_distribution` on 3D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_uniform_partition_normal_distribution_3D, self).createData() + super(test_uniform_partition_normal_distribution_3D, self).setUp() + + class uniform_hyperrectangle_base(prob_uniform): """ Provides set up and a test to check the number of ``d_distr_samples`` for @@ -1135,3 +1212,74 @@ def setUp(self): """ super(test_uniform_partition_uniform_distribution_rectangle_domain_3D, self).createData() super(test_uniform_partition_uniform_distribution_rectangle_domain_3D, self).setUp() + + +class user_partition_user_distribution(prob): + """ + Set up :meth:`bet.calculateP.simpleFunP.user_partition_user_distribution` on data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + self.data_prob = sFun.user_partition_user_distribution(self.data, + self.data, + self.data) + self.rho_D_M = self.data_prob.get_probabilities() + self.d_distr_samples = self.data_prob.get_values() + +class test_user_partition_user_distribution_01D(data_01D, + user_partition_user_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.user_partition_user_distribution` on 01D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_user_partition_user_distribution_01D, self).createData() + super(test_user_partition_user_distribution_01D, self).setUp() + + +class test_user_partition_user_distribution_1D(data_1D, + user_partition_user_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.user_partition_user_distribution` on 1D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_user_partition_user_distribution_1D, self).createData() + super(test_user_partition_user_distribution_1D, self).setUp() + + +class test_user_partition_user_distribution_2D(data_2D, + user_partition_user_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.user_partition_user_distribution` on 2D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_user_partition_user_distribution_2D, self).createData() + super(test_user_partition_user_distribution_2D, self).setUp() + + +class test_user_partition_user_distribution_3D(data_3D, + user_partition_user_distribution): + """ + Tests :meth:`bet.calculateP.simpleFunP.user_partition_user_distribution` on 3D data domain. + """ + + def setUp(self): + """ + Set up problem. + """ + super(test_user_partition_user_distribution_3D, self).createData() + super(test_user_partition_user_distribution_3D, self).setUp() From ebed516baecd00e2f19aebd6959fe1ae61a557ac Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Tue, 24 May 2016 16:20:05 -0600 Subject: [PATCH 19/20] Raising exceptions for incorrect argument types --- bet/calculateP/simpleFunP.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index d0f84fae..5ccb8c28 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -78,8 +78,9 @@ def uniform_partition_uniform_distribution_rectangle_size(data_set, Q_ref, if not isinstance(rect_size, collections.Iterable): rect_size = rect_size * np.ones((dim,)) - if np.any(np.less(rect_size, 0)): - print 'Warning: rect_size must be greater than 0' + if np.any(np.less_equal(rect_size, 0)): + msg = 'rect_size must be greater than 0' + raise wrong_argument_type(msg) r''' Create M samples defining M Voronoi cells (i.e., "bins") in D used to @@ -340,12 +341,14 @@ def regular_partition_uniform_distribution_rectangle_size(data_set, Q_ref, rect_ msg = 'center_pts_per_edge dimension mismatch.' msg += 'Using 1 in each dimension.' logging.warning(msg) - if np.any(np.less(center_pts_per_edge, 0)): - print 'Warning: center_pts_per_edge must be greater than 0' + if np.any(np.less_equal(center_pts_per_edge, 0)): + msg = 'center_pts_per_edge must be greater than 0' + raise wrong_argument_type(msg) if not isinstance(rect_size, collections.Iterable): rect_size = rect_size * np.ones((dim,)) - if np.any(np.less(rect_size, 0)): - print 'Warning: rect_size must be greater than 0' + if np.any(np.less_equal(rect_size, 0)): + msg = 'rect_size must be greater than 0' + raise wrong_argument_type(msg) sur_domain = np.array([np.min(data, 0), np.max(data, 0)]).transpose() From 79a2547fcfc11b837714433723f509495c9b937a Mon Sep 17 00:00:00 2001 From: Troy Butler Date: Tue, 24 May 2016 18:08:03 -0600 Subject: [PATCH 20/20] Fixed user simpleFunP to work in parallel --- bet/calculateP/simpleFunP.py | 68 +++++++++---------- .../linearMap/linearMapUniformSampling.py | 1 + 2 files changed, 35 insertions(+), 34 deletions(-) diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py index 5ccb8c28..4ba7c2f6 100644 --- a/bet/calculateP/simpleFunP.py +++ b/bet/calculateP/simpleFunP.py @@ -734,36 +734,34 @@ def user_partition_user_distribution(data_set, data_partition_set, raise wrong_argument_type(msg) if isinstance(data_partition_set, samp.sample_set_base): - num_samples_discretize_D = data_partition_set.check_num() - simpleFun_set = data_partition_set.copy() - dim_simpleFun = simpleFun_set._dim + M = data_partition_set.check_num() + d_distr_samples = data_partition_set._values + dim_simpleFun = d_distr_samples.shape[1] elif isinstance(data_partition_set, samp.discretization): - num_samples_discretize_D = data_partition_set.check_nums() - simpleFun_set = data_partition_set._output_sample_set.copy() - dim_simpleFun = simpleFun_set._dim + M = data_partition_set.check_nums() + d_distr_samples = data_partition_set._output_sample_set._values + dim_simpleFun = d_distr_samples.shape[1] elif isinstance(data_partition_set, np.ndarray): - num_samples_discretize_D = data_partition_set.shape[0] + M = data_partition_set.shape[0] dim_simpleFun = data_partition_set.shape[1] - simpleFun_set = samp.sample_set(dim=dim_simpleFun) - simpleFun_set.set_values(values) + d_distr_samples = data_partition_set else: msg = "The second argument must be of type bet.sample.sample_set, " msg += "bet.sample.discretization or np.ndarray" raise wrong_argument_type(msg) if isinstance(data_distribution_set, samp.sample_set_base): - MonteCarlo_set = data_distribution_set.copy() - dim_MonteCarlo = MonteCarlo_set._dim - num_iid_samples = data_distribution_set.check_num() + d_distr_emulate = data_distribution_set._values + dim_MonteCarlo = d_distr_emulate.shape[1] + num_d_emulate = data_distribution_set.check_num() elif isinstance(data_distribution_set, samp.discretization): - MonteCarlo_set = data_distribution_set._output_sample_set.copy() - dim_MonteCarlo = MonteCarlo_set._dim - num_iid_samples = data_distribution_set.check_nums() + d_distr_emulate = data_distribution_set._output_sample_set._values + dim_MonteCarlo = d_distr_emulate.shape[1] + num_d_emulate = data_distribution_set.check_nums() elif isinstance(data_distribution_set, np.ndarray): - num_iid_samples = data_distribution_set.shape[0] + num_d_emulate = data_distribution_set.shape[0] dim_MonteCarlo = data_distribution_set.shape[1] - MonteCarlo_set = samp.sample_set(dim=dim_MonteCarlo) - MonteCarlo_set.set_values(values) + d_distr_emulate = data_distribution_set else: msg = "The second argument must be of type bet.sample.sample_set, " msg += "bet.sample.discretization or np.ndarray" @@ -773,24 +771,26 @@ def user_partition_user_distribution(data_set, data_partition_set, msg = "The argument types have conflicting dimensions" raise wrong_argument_type(msg) - my_discretization = samp.discretization(input_sample_set=s_set, - output_sample_set=s_set, - emulated_output_sample_set=MonteCarlo_set, - output_probability_set=simpleFun_set) + # Initialize sample set object + s_set = samp.sample_set(dim) + s_set.set_values(d_distr_samples) + s_set.set_kdtree() - my_discretization.set_emulated_oo_ptr() + (_, k) = s_set.query(d_distr_emulate) - my_discretization._output_probability_set.set_probabilities( - np.zeros((num_samples_discretize_D,))) + count_neighbors = np.zeros((M,), dtype=np.int) + for i in range(M): + count_neighbors[i] = np.sum(np.equal(k, i)) - for i in range(num_samples_discretize_D): - Itemp = np.equal(my_discretization.get_emulated_oo_ptr(), i) - Itemp_sum = float(np.sum(Itemp)) - Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM) - if Itemp_sum > 0: - my_discretization._output_probability_set._probabilities[i] = \ - Itemp_sum / num_iid_samples + # Use the binning to define :math:`\rho_{\mathcal{D},M}` + ccount_neighbors = np.copy(count_neighbors) + comm.Allreduce([count_neighbors, MPI.INT], [ccount_neighbors, MPI.INT], + op=MPI.SUM) + count_neighbors = ccount_neighbors + rho_D_M = count_neighbors.astype(np.float64) / \ + float(num_d_emulate * comm.size) + s_set.set_probabilities(rho_D_M) if isinstance(data_set, samp.discretization): - data_set._output_probability_set = my_discretization._output_probability_set - return my_discretization._output_probability_set + data_set._output_probability_set = s_set + return s_set diff --git a/examples/linearMap/linearMapUniformSampling.py b/examples/linearMap/linearMapUniformSampling.py index 260ff77b..b578d18b 100644 --- a/examples/linearMap/linearMapUniformSampling.py +++ b/examples/linearMap/linearMapUniformSampling.py @@ -146,6 +146,7 @@ # calculate 2d marginal probs (bins, marginals2D) = plotP.calculate_2D_marginal_probs(input_samples, nbins = [10, 10, 10]) + # smooth 2d marginals probs (optional) marginals2D = plotP.smooth_marginals_2D(marginals2D, bins, sigma=0.2)