Merge pull request #195 from eecsu/v2_master

v2 master -- updates to basicSampling, simpleFunP, some tests, and some examples

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)

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):
     """
@@ -167,9 +169,9 @@ def random_sample_set_domain(self, sample_type, input_domain,
         :param string criterion: latin hypercube criterion see
             `PyDOE <http://pythonhosted.org/pyDOE/randomized.html>`_
 
-        :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
@@ -201,9 +203,9 @@ def random_sample_set_dimension(self, sample_type, input_dim,
         :param string criterion: latin hypercube criterion see
             `PyDOE <http://pythonhosted.org/pyDOE/randomized.html>`_
 
-        :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
@@ -212,8 +214,112 @@ def random_sample_set_dimension(self, sample_type, input_dim,
         return self.random_sample_set(sample_type, input_sample_set,
                                        num_samples, criterion)
 
+    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.sample_set`
+        :returns: :class:`~bet.sample.sample_Set` object which contains
+            input ``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_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)
+
+        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]
+
+        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):
+            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)
+
+        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.sample_set`
+        :returns: :class:`~bet.sample.sample_Set` object which contains
+            input ``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_dim=1):
+        """
+        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.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_dim)
+
     def compute_QoI_and_create_discretization(self, input_sample_set,
-            savefile=None, parallel=False): 
+            savefile=None, parallel=False):
         """
         Samples the model at ``input_sample_set`` and saves the results.
 

examples/FEniCS/BET_script.py

@@ -0,0 +1,163 @@
+#! /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",
+                             lam_ref=param_ref[0,:], 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",
+                             lam_ref=param_ref[0,:], file_extension=".eps")
+
+
+