From aaab1adbfa5bc36146298800af2bdd7d976f4656 Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 10:57:00 -0400
Subject: [PATCH 1/9] fixed type checking in postProcess and basicSampling, see
 issue #200

---
 bet/postProcess/plotDomains.py | 26 +++++++++++++-------------
 bet/postProcess/plotP.py       | 16 ++++++++--------
 bet/postProcess/postTools.py   | 14 +++++++-------
 bet/sampling/basicSampling.py  |  2 +-
 4 files changed, 29 insertions(+), 29 deletions(-)

diff --git a/bet/postProcess/plotDomains.py b/bet/postProcess/plotDomains.py
index fbc1e723..b0d8e155 100644
--- a/bet/postProcess/plotDomains.py
+++ b/bet/postProcess/plotDomains.py
@@ -52,7 +52,7 @@ def scatter_2D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
     the other markers.
 
     :param sample_obj: contains samples to create scatter plot
-    :type sample_obj: :class:`~bet.sample.sample_set`
+    :type sample_obj: :class:`~bet.sample.sample_set_base`
     :param list sample_nos: indicies of the samples to plot
     :param color: values to color the samples by
     :type color: :class:`numpy.ndarray`
@@ -65,7 +65,7 @@ def scatter_2D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
     :param string filename: filename to save the figure as
 
     """
-    if not isinstance(sample_obj, sample.sample_set):
+    if not isinstance(sample_obj, sample.sample_set_base):
         raise bad_object("Improper sample object")
     # check dimension of data to plot
     if sample_obj.get_dim() != 2:
@@ -119,7 +119,7 @@ def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
     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_base`
     :param list sample_nos: indicies of the samples to plot
     :param color: values to color the samples by
     :type color: :class:`numpy.ndarray`
@@ -133,7 +133,7 @@ def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
     :param string filename: filename to save the figure as
 
     """
-    if not isinstance(sample_obj, sample.sample_set):
+    if not isinstance(sample_obj, sample.sample_set_base):
         raise bad_object("Improper sample object")
     # check dimension of data to plot
     if sample_obj.get_dim() != 3:
@@ -191,7 +191,7 @@ def show_param(sample_disc, rho_D=None, p_ref=None, sample_nos=None,
 
     :param sample_disc: Object containing the samples to plot
     :type sample_disc: :class:`~bet.sample.discretization` 
-        or :class:`~bet.sample.sample_set`
+        or :class:`~bet.sample.sample_set_base`
     :param list sample_nos: sample numbers to plot
     :param rho_D: probability density function on D
     :type rho_D: callable function that takes a :class:`np.array` and returns a
@@ -217,7 +217,7 @@ def show_param(sample_disc, rho_D=None, p_ref=None, sample_nos=None,
     else:
         if isinstance(sample_disc, sample.discretization):
             sample_obj = sample_disc._input_sample_set
-        elif isinstance(sample_disc, sample.sample_set):
+        elif isinstance(sample_disc, sample.sample_set_base):
             sample_obj = sample_disc
         else:
             raise bad_object("Improper sample object")
@@ -271,7 +271,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_base`
     :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
@@ -363,7 +363,7 @@ def show_data_domain_multi(sample_disc, Q_ref=None, Q_nums=None,
 
     :param sample_disc: Object containing the samples to plot
     :type sample_disc: :class:`~bet.sample.discretization` or 
-        :class:`~bet.sample.sample_set`
+        :class:`~bet.sample.sample_set_base`
     :param Q_ref: reference data value
     :type Q_ref: :class:`numpy.ndarray` of shape (M, mdim)
     :param list Q_nums: dimensions of the QoI to plot
@@ -472,7 +472,7 @@ def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None,
 
     :param sample_disc: Object containing the samples to plot
     :type sample_disc: :class:`~bet.sample.discretization` 
-        or :class:`~bet.sample.sample_set`
+        or :class:`~bet.sample.sample_set_base`
     :param Q_ref: reference data value
     :type Q_ref: :class:`numpy.ndarray` of shape (M, 2)
     :param list ref_markers: list of marker types for :math:`Q_{ref}`
@@ -542,7 +542,7 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all',
     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`
+    :type sample_obj: :class:`~bet.sample.sample_set_base`
     :param bool save: flag whether or not to save the figure
     :param bool interactive: flag whether or not to show the figure
     :param string img_folder: folder to save the plots to
@@ -551,7 +551,7 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all',
     :type showdim: int or string
 
     """
-    if not isinstance(sample_obj, sample.sample_set):
+    if not isinstance(sample_obj, sample.sample_set_base):
         raise bad_object("Improper sample object")
 
     # If no specific coordinate number of choice is given set to be the first
@@ -622,7 +622,7 @@ def scatter2D_multi(sample_obj, color=None, p_ref=None, img_folder='figs/',
     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_base`
     :param color: values to color the ``samples`` by
     :type color: :class:`numpy.ndarray`
     :param string filename: filename to save the figure as
@@ -635,7 +635,7 @@ def scatter2D_multi(sample_obj, color=None, p_ref=None, img_folder='figs/',
     :type showdim: int or string
 
     """
-    if not isinstance(sample_obj, sample.sample_set):
+    if not isinstance(sample_obj, sample.sample_set_base):
         raise bad_object("Improper sample object")
     # If no specific coordinate number of choice is given set to be the first
     # coordinate direction.
diff --git a/bet/postProcess/plotP.py b/bet/postProcess/plotP.py
index 2fb33ebe..e3b804a7 100644
--- a/bet/postProcess/plotP.py
+++ b/bet/postProcess/plotP.py
@@ -32,7 +32,7 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
     that the probabilities to be plotted are from the input space.
 
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or 
+    :type sample_set: :class:`~bet.sample.sample_set_base` or 
         :class:`~bet.sample.discretization`
     :param nbins: Number of bins in each direction.
     :type nbins: :int or :class:`~numpy.ndarray` of shape (ndim,)
@@ -42,7 +42,7 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
     """
     if isinstance(sample_set, sample.discretization):
         sample_obj = sample_set._input_sample_set
-    elif isinstance(sample_set, sample.sample_set):
+    elif isinstance(sample_set, sample.sample_set_base):
         sample_obj = sample_set
     else:
         raise bad_object("Improper sample object")
@@ -78,7 +78,7 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
     that the probabilities to be plotted are from the input space.
 
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` 
+    :type sample_set: :class:`~bet.sample.sample_set_base` 
         or :class:`~bet.sample.discretization`
     :param nbins: Number of bins in each direction.
     :type nbins: :int or :class:`~numpy.ndarray` of shape (ndim,)
@@ -88,7 +88,7 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
     """
     if isinstance(sample_set, sample.discretization):
         sample_obj = sample_set._input_sample_set
-    elif isinstance(sample_set, sample.sample_set):
+    elif isinstance(sample_set, sample.sample_set_base):
         sample_obj = sample_set
     else:
         raise bad_object("Improper sample object")
@@ -139,7 +139,7 @@ def plot_1D_marginal_probs(marginals, bins, sample_set,
         calculating marginals
     :type bins: :class:`~numpy.ndarray` of shape (nbins+1,)
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` 
+    :type sample_set: :class:`~bet.sample.sample_set_base` 
         or :class:`~bet.sample.discretization`
     :param filename: Prefix for output files.
     :type filename: str
@@ -153,7 +153,7 @@ def plot_1D_marginal_probs(marginals, bins, sample_set,
     """
     if isinstance(sample_set, sample.discretization):
         sample_obj = sample_set._input_sample_set
-    elif isinstance(sample_set, sample.sample_set):
+    elif isinstance(sample_set, sample.sample_set_base):
         sample_obj = sample_set
     else:
         raise bad_object("Improper sample object")
@@ -202,7 +202,7 @@ def plot_2D_marginal_probs(marginals, bins, sample_set,
     :param bins: Endpoints of bins used in calculating marginals
     :type bins: :class:`~numpy.ndarray` of shape (nbins+1,2)
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` 
+    :type sample_set: :class:`~bet.sample.sample_set_base` 
         or :class:`~bet.sample.discretization`
     :param filename: Prefix for output files.
     :type filename: str
@@ -216,7 +216,7 @@ def plot_2D_marginal_probs(marginals, bins, sample_set,
     """
     if isinstance(sample_set, sample.discretization):
         sample_obj = sample_set._input_sample_set
-    elif isinstance(sample_set, sample.sample_set):
+    elif isinstance(sample_set, sample.sample_set_base):
         sample_obj = sample_set
     else:
         raise bad_object("Improper sample object")
diff --git a/bet/postProcess/postTools.py b/bet/postProcess/postTools.py
index b8510054..33be8d52 100644
--- a/bet/postProcess/postTools.py
+++ b/bet/postProcess/postTools.py
@@ -28,7 +28,7 @@ def sort_by_rho(sample_set):
     are also sorted.
 
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or
+    :type sample_set: :class:`~bet.sample.sample_set_base` or
         :class:`~bet.sample.discretization` 
     :param indices: sorting indices
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
@@ -45,7 +45,7 @@ def sort_by_rho(sample_set):
         P_samples = sample_set._input_sample_set.get_probabilities()
         lam_vol = sample_set._input_sample_set.get_volumes()
         data = sample_set._output_sample_set.get_values()
-    elif isinstance(sample_set, sample.sample_set):
+    elif isinstance(sample_set, sample.sample_set_base):
         samples = sample_set.get_values()
         P_samples = sample_set.get_probabilities()
         lam_vol = sample_set.get_volumes()
@@ -95,7 +95,7 @@ def sample_prob(percentile, sample_set, sort=True, descending=False):
     :param percentile: ratio of highest probability samples to select
     :type percentile: float
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or 
+    :type sample_set: :class:`~bet.sample.sample_set_base` or 
         :class:`~bet.sample.discretization`
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param indices: sorting indices
@@ -114,7 +114,7 @@ def sample_prob(percentile, sample_set, sort=True, descending=False):
         P_samples = sample_set._input_sample_set.get_probabilities()
         lam_vol = sample_set._input_sample_set.get_volumes()
         data = sample_set._output_sample_set.get_values()
-    elif isinstance(sample_set, sample.sample_set):
+    elif isinstance(sample_set, sample.sample_set_base):
         samples = sample_set.get_values()
         P_samples = sample_set.get_probabilities()
         lam_vol = sample_set.get_volumes()
@@ -129,7 +129,7 @@ def sample_prob(percentile, sample_set, sort=True, descending=False):
             P_samples = sample_set._input_sample_set.get_probabilities()
             lam_vol = sample_set._input_sample_set.get_volumes()
             data = sample_set._output_sample_set.get_values()
-        elif isinstance(sample_set, sample.sample_set):
+        elif isinstance(sample_set, sample.sample_set_base):
             samples = sample_set.get_values()
             P_samples = sample_set.get_probabilities()
             lam_vol = sample_set.get_volumes()
@@ -182,7 +182,7 @@ def sample_highest_prob(top_percentile, sample_set, sort=True):
     :param top_percentile: ratio of highest probability samples to select
     :type top_percentile: float
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` 
+    :type sample_set: :class:`~bet.sample.sample_set_base` 
         or :class:`~bet.sample.discretization`
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param indices: sorting indices
@@ -208,7 +208,7 @@ def sample_lowest_prob(bottom_percentile, sample_set, sort=True):
     :param top_percentile: ratio of highest probability samples to select
     :type top_percentile: float
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` 
+    :type sample_set: :class:`~bet.sample.sample_set_base` 
         or :class:`~bet.sample.discretization`
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param indices: sorting indices of unsorted ``P_samples``
diff --git a/bet/sampling/basicSampling.py b/bet/sampling/basicSampling.py
index 691edc7e..74b7942c 100644
--- a/bet/sampling/basicSampling.py
+++ b/bet/sampling/basicSampling.py
@@ -415,7 +415,7 @@ def create_random_discretization(self, sample_type, input_obj,
         if num_samples is None:
             num_samples = self.num_samples
 
-        if isinstance(input_obj, sample.sample_set):
+        if isinstance(input_obj, sample.sample_set_base):
             input_sample_set = self.random_sample_set(sample_type, input_obj,
                     num_samples, criterion)
         elif isinstance(input_obj, np.ndarray):

From 689ca0ff95c31ffbb4ac2f0b4f2fe2c213338e14 Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 15:03:43 -0400
Subject: [PATCH 2/9] added todo and updated some api documentation

---
 bet/sample.py        | 7 +++++++
 doc/bet.sampling.rst | 8 ++++++++
 2 files changed, 15 insertions(+)

diff --git a/bet/sample.py b/bet/sample.py
index 4f336193..97ef66fd 100644
--- a/bet/sample.py
+++ b/bet/sample.py
@@ -610,6 +610,13 @@ def estimate_volume(self, n_mc_points=int(1E4)):
         Calculate the volume faction of cells approximately using Monte
         Carlo integration. 
 
+        .. todo::
+
+           This currently presumes a uniform Lesbegue measure on the
+           ``domain``. Currently the way this is written
+           ``emulated_input_sample_set`` is NOT used to calculate the volume.
+           This should at least be an option. 
+
         :param int n_mc_points: If estimate is True, number of MC points to use
         """
         num = self.check_num()
diff --git a/doc/bet.sampling.rst b/doc/bet.sampling.rst
index 7a53b1b9..b443ee02 100644
--- a/doc/bet.sampling.rst
+++ b/doc/bet.sampling.rst
@@ -4,6 +4,14 @@ bet.sampling package
 Submodules
 ----------
 
+bet.sampling.LpGeneralizedSamples module
+----------------------------------------
+
+.. automodule:: bet.sampling.LpGeneralizedSamples
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 bet.sampling.adaptiveSampling module
 ------------------------------------
 

From 122bbcb03f92317c8a2be86f149a05ef73990c6c Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 15:04:01 -0400
Subject: [PATCH 3/9] updating examples

---
 examples/fromFile_ADCIRCMap/Q_1D_serial.py    | 55 ++++++++---------
 examples/fromFile_ADCIRCMap/Q_2D_parallel.py  | 61 ++++++++-----------
 examples/fromFile_ADCIRCMap/Q_2D_serial.py    | 56 ++++++++---------
 examples/fromFile_ADCIRCMap/Q_3D_serial.py    | 32 +++++-----
 examples/fromFile_ADCIRCMap/fromFile2D.py     |  6 +-
 examples/fromFile_ADCIRCMap/fromFile3D.py     |  6 +-
 examples/fromFile_ADCIRCMap/plotDomains2D.py  | 24 +++++---
 examples/fromFile_ADCIRCMap/plotDomains3D.py  | 22 ++++---
 .../fromFile_ADCIRCMap/sandbox_test_2D.py     | 17 +++---
 .../fromFile_ADCIRCMap/sandbox_test_3D.py     | 16 +++--
 10 files changed, 141 insertions(+), 154 deletions(-)

diff --git a/examples/fromFile_ADCIRCMap/Q_1D_serial.py b/examples/fromFile_ADCIRCMap/Q_1D_serial.py
index 6288d928..8857f309 100644
--- a/examples/fromFile_ADCIRCMap/Q_1D_serial.py
+++ b/examples/fromFile_ADCIRCMap/Q_1D_serial.py
@@ -4,16 +4,23 @@
 import bet.calculateP.simpleFunP as sfun
 import numpy as np
 import scipy.io as sio
+import bet.sample as sample
 
 # Import "Truth"
-mdat = sio.loadmat('Q_2D')
+mdat = sio.loadmat('../matfiles/Q_2D')
 Q = mdat['Q']
 Q_ref = mdat['Q_true']
 
 # Import Data
-samples = mdat['points'].transpose()
+points = mdat['points']
 lam_domain = np.array([[0.07, .15], [0.1, 0.2]])
 
+# Create input, output, and discretization from data read from file
+input_sample_set = sample.sample_set(points.shape[0])
+input_sample_set.set_values(points.transpose())
+input_sample_set.set_domain(lam_domain)
+
+
 print "Finished loading data"
 
 def postprocess(station_nums, ref_num):
@@ -24,55 +31,41 @@ def postprocess(station_nums, ref_num):
     filename += '_ref_'+str(ref_num+1)
 
     data = Q[:, station_nums]
+    output_sample_set = sample.sample_set(data.shape[1])
+    output_sample_set.set_values(data)
     q_ref = Q_ref[ref_num, station_nums]
 
     # Create Simple function approximation
     # Save points used to parition D for simple function approximation and the
     # approximation itself (this can be used to make close comparisions...)
-    (rho_D_M, d_distr_samples, d_Tree) = sfun.uniform_hyperrectangle(data,
-            q_ref, bin_ratio=0.15,
+    output_probability_set = sfun.regular_partition_uniform_distribution_rectangle_scaled(\
+            output_sample_set, q_ref, rect_scale=0.15,
             center_pts_per_edge=np.ones((data.shape[1],)))
 
     num_l_emulate = 1e6
-    lambda_emulate = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
-    print "Finished emulating lambda samples"
+    set_emulated = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
+    my_disc = sample.discretization(input_sample_set, output_sample_set,
+            output_probability_set, emulated_input_sample_set=set_emulated)
 
-    mdict = dict()
-    mdict['rho_D_M'] = rho_D_M
-    mdict['d_distr_samples'] = d_distr_samples 
-    mdict['num_l_emulate'] = num_l_emulate
-    mdict['lambda_emulate'] = lambda_emulate
+    print "Finished emulating lambda samples"
 
     # Calculate P on lambda emulate
-    (P0, lem0, io_ptr0, emulate_ptr0) = calcP.prob_emulated(samples, data,
-            rho_D_M, d_distr_samples, lambda_emulate, d_Tree)
     print "Calculating prob_emulated"
-    mdict['P0'] = P0
-    mdict['lem0'] = lem0
-    mdict['io_ptr0'] = io_ptr0
-    mdict['emulate_ptr0'] = emulate_ptr0
+    calcP.prob_emulated(my_disc)
+    sample.save_discretization(my_disc, filename, "prob_emulated_solution")
 
     # Calclate P on the actual samples with assumption that voronoi cells have
     # equal size
-    (P1, lam_vol1, io_ptr1) = calcP.prob(samples, data,
-            rho_D_M, d_distr_samples, d_Tree)
+    input_sample_set.estimate_volume_mc()
     print "Calculating prob"
-    mdict['P1'] = P1
-    mdict['lam_vol1'] = lam_vol1
-    mdict['lem1'] = samples
-    mdict['io_ptr1'] = io_ptr1
+    calcP.prob(my_disc)
+    sample.save_discretization(my_disc, filename, "prob_solution")
 
     # Calculate P on the actual samples estimating voronoi cell volume with MC
     # integration
-    (P3, lam_vol3, lambda_emulate3, io_ptr3, emulate_ptr3) = calcP.prob_mc(samples,
-            data, rho_D_M, d_distr_samples, lambda_emulate, d_Tree)
+    calcP.prob_mc(my_disc)
     print "Calculating prob_mc"
-    mdict['P3'] = P3
-    mdict['lam_vol3'] = lam_vol3
-    mdict['io_ptr3'] = io_ptr3
-    mdict['emulate_ptr3'] = emulate_ptr3
-    # Export P
-    sio.savemat(filename, mdict, do_compression=True)
+    sample.save_discretization(my_disc, filename, "prob_mc_solution")
 
 # Post-process and save P and emulated points
 ref_nums = [6, 11, 15] # 7, 12, 16
diff --git a/examples/fromFile_ADCIRCMap/Q_2D_parallel.py b/examples/fromFile_ADCIRCMap/Q_2D_parallel.py
index e05c8fac..bdad804f 100644
--- a/examples/fromFile_ADCIRCMap/Q_2D_parallel.py
+++ b/examples/fromFile_ADCIRCMap/Q_2D_parallel.py
@@ -6,16 +6,21 @@
 import scipy.io as sio
 import bet.util as util
 from bet.Comm import comm
+import bet.sample as sample
 
 # Import "Truth"
-mdat = sio.loadmat('Q_2D')
+mdat = sio.loadmat('../matfiles/Q_2D')
 Q = mdat['Q']
 Q_ref = mdat['Q_true']
 
 # Import Data
-samples = mdat['points'].transpose()
+points = mdat['points']
 lam_domain = np.array([[0.07, .15], [0.1, 0.2]])
 
+# Create input, output, and discretization from data read from file
+input_sample_set = sample.sample_set(points.shape[0])
+input_sample_set.set_values(points.transpose())
+input_sample_set.set_domain(lam_domain)
 print "Finished loading data"
 
 def postprocess(station_nums, ref_num):
@@ -26,58 +31,44 @@ def postprocess(station_nums, ref_num):
     filename += '_ref_'+str(ref_num+1)
 
     data = Q[:, station_nums]
+    output_sample_set = sample.sample_set(data.shape[1])
+    output_sample_set.set_values(data)
     q_ref = Q_ref[ref_num, station_nums]
 
     # Create Simple function approximation
     # Save points used to parition D for simple function approximation and the
     # approximation itself (this can be used to make close comparisions...)
-    (rho_D_M, d_distr_samples, d_Tree) = sfun.uniform_hyperrectangle(data,
-            q_ref, bin_ratio=0.15,
+    output_probability_set = sfun.regular_partition_uniform_distribution_rectangle_scaled(\
+            output_sample_set, q_ref, rect_scale=0.15,
             center_pts_per_edge=np.ones((data.shape[1],)))
 
     num_l_emulate = 1e6
-    lambda_emulate = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
-    
-    if comm.rank == 0:
-        print "Finished emulating lambda samples"
-        mdict = dict()
-        mdict['rho_D_M'] = rho_D_M
-        mdict['d_distr_samples'] = d_distr_samples 
-        mdict['num_l_emulate'] = num_l_emulate
+    set_emulated = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
+    my_disc = sample.discretization(input_sample_set, output_sample_set,
+            output_probability_set, emulated_input_sample_set=set_emulated)
+
+    print "Finished emulating lambda samples"
 
     # Calculate P on lambda emulate
-    (P0, lem0, io_ptr0, emulate_ptr0) = calcP.prob_emulated(samples, data,
-            rho_D_M, d_distr_samples, lambda_emulate, d_Tree)
+    print "Calculating prob_emulated"
+    calcP.prob_emulated(my_disc)
     if comm.rank == 0:
-        print "Calculating prob_emulated"
-        mdict['P0'] = P0
-        mdict['lem0'] = lem0
-        mdict['io_ptr0'] = io_ptr0
-        mdict['emulate_ptr0'] = emulate_ptr0
+        sample.save_discretization(my_disc, filename, "prob_emulated_solution")
 
     # Calclate P on the actual samples with assumption that voronoi cells have
     # equal size
-    (P1, lam_vol1, io_ptr1) = calcP.prob(samples, data,
-            rho_D_M, d_distr_samples, d_Tree)
+    input_sample_set.estimate_volume_mc()
+    print "Calculating prob"
+    calcP.prob(my_disc)
     if comm.rank == 0:
-        print "Calculating prob"
-        mdict['P1'] = P1
-        mdict['lam_vol1'] = lam_vol1
-        mdict['lem1'] = samples
-        mdict['io_ptr1'] = io_ptr1
+        sample.save_discretization(my_disc, filename, "prob_solution")
 
     # Calculate P on the actual samples estimating voronoi cell volume with MC
     # integration
-    (P3, lam_vol3, lambda_emulate3, io_ptr3, emulate_ptr3) = calcP.prob_mc(samples,
-            data, rho_D_M, d_distr_samples, lambda_emulate, d_Tree)
+    calcP.prob_mc(my_disc)
+    print "Calculating prob_mc"
     if comm.rank == 0:
-        print "Calculating prob_mc"
-        mdict['P3'] = P3
-        mdict['lam_vol3'] = lam_vol3
-        mdict['io_ptr3'] = io_ptr3
-        mdict['emulate_ptr3'] = emulate_ptr3
-        # Export P
-        sio.savemat(filename, mdict, do_compression=True)
+        sample.save_discretization(my_disc, filename, "prob_mc_solution")
 
 # Post-process and save P and emulated points
 ref_nums = [6, 11, 15] # 7, 12, 16
diff --git a/examples/fromFile_ADCIRCMap/Q_2D_serial.py b/examples/fromFile_ADCIRCMap/Q_2D_serial.py
index d3d50c3c..86862a11 100644
--- a/examples/fromFile_ADCIRCMap/Q_2D_serial.py
+++ b/examples/fromFile_ADCIRCMap/Q_2D_serial.py
@@ -4,16 +4,21 @@
 import bet.calculateP.simpleFunP as sfun
 import numpy as np
 import scipy.io as sio
+import bet.sample as sample
 
 # Import "Truth"
-mdat = sio.loadmat('Q_2D')
+mdat = sio.loadmat('../matfiles/Q_2D')
 Q = mdat['Q']
 Q_ref = mdat['Q_true']
 
 # Import Data
-samples = mdat['points'].transpose()
+points = mdat['points']
 lam_domain = np.array([[0.07, .15], [0.1, 0.2]])
 
+# Create input, output, and discretization from data read from file
+input_sample_set = sample.sample_set(points.shape[0])
+input_sample_set.set_values(points.transpose())
+input_sample_set.set_domain(lam_domain)
 print "Finished loading data"
 
 def postprocess(station_nums, ref_num):
@@ -24,55 +29,44 @@ def postprocess(station_nums, ref_num):
     filename += '_ref_'+str(ref_num+1)
 
     data = Q[:, station_nums]
+    output_sample_set = sample.sample_set(data.shape[1])
+    output_sample_set.set_values(data)
     q_ref = Q_ref[ref_num, station_nums]
 
     # Create Simple function approximation
     # Save points used to parition D for simple function approximation and the
     # approximation itself (this can be used to make close comparisions...)
-    (rho_D_M, d_distr_samples, d_Tree) = sfun.uniform_hyperrectangle(data,
-            q_ref, bin_ratio=0.15,
+    output_probability_set = sfun.regular_partition_uniform_distribution_rectangle_scaled(\
+            output_sample_set, q_ref, rect_scale=0.15,
             center_pts_per_edge=np.ones((data.shape[1],)))
 
     num_l_emulate = 1e6
-    lambda_emulate = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
-    print "Finished emulating lambda samples"
+    set_emulated = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
+    my_disc = sample.discretization(input_sample_set, output_sample_set,
+            output_probability_set, emulated_input_sample_set=set_emulated)
 
-    mdict = dict()
-    mdict['rho_D_M'] = rho_D_M
-    mdict['d_distr_samples'] = d_distr_samples 
-    mdict['num_l_emulate'] = num_l_emulate
-    mdict['lambda_emulate'] = lambda_emulate
+    print "Finished emulating lambda samples"
 
     # Calculate P on lambda emulate
-    (P0, lem0, io_ptr0, emulate_ptr0) = calcP.prob_emulated(samples, data,
-            rho_D_M, d_distr_samples, lambda_emulate, d_Tree)
     print "Calculating prob_emulated"
-    mdict['P0'] = P0
-    mdict['lem0'] = lem0
-    mdict['io_ptr0'] = io_ptr0
-    mdict['emulate_ptr0'] = emulate_ptr0
+    calcP.prob_emulated(my_disc)
+    if comm.rank == 0:
+        sample.save_discretization(my_disc, filename, "prob_emulated_solution")
 
     # Calclate P on the actual samples with assumption that voronoi cells have
     # equal size
-    (P1, lam_vol1, io_ptr1) = calcP.prob(samples, data,
-            rho_D_M, d_distr_samples, d_Tree)
+    input_sample_set.estimate_volume_mc()
     print "Calculating prob"
-    mdict['P1'] = P1
-    mdict['lam_vol1'] = lam_vol1
-    mdict['lem1'] = samples
-    mdict['io_ptr1'] = io_ptr1
+    calcP.prob(my_disc)
+    if comm.rank == 0:
+        sample.save_discretization(my_disc, filename, "prob_solution")
 
     # Calculate P on the actual samples estimating voronoi cell volume with MC
     # integration
-    (P3, lam_vol3, lambda_emulate3, io_ptr3, emulate_ptr3) = calcP.prob_mc(samples,
-            data, rho_D_M, d_distr_samples, lambda_emulate, d_Tree)
+    calcP.prob_mc(my_disc)
     print "Calculating prob_mc"
-    mdict['P3'] = P3
-    mdict['lam_vol3'] = lam_vol3
-    mdict['io_ptr3'] = io_ptr3
-    mdict['emulate_ptr3'] = emulate_ptr3
-    # Export P 
-    sio.savemat(filename, mdict, do_compression=True)
+    if comm.rank == 0:
+        sample.save_discretization(my_disc, filename, "prob_mc_solution")
 
 # Post-process and save P and emulated points
 ref_nums = [6, 11, 15] # 7, 12, 16
diff --git a/examples/fromFile_ADCIRCMap/Q_3D_serial.py b/examples/fromFile_ADCIRCMap/Q_3D_serial.py
index b9656cf3..d4dfd8c5 100644
--- a/examples/fromFile_ADCIRCMap/Q_3D_serial.py
+++ b/examples/fromFile_ADCIRCMap/Q_3D_serial.py
@@ -4,9 +4,10 @@
 import bet.calculateP.simpleFunP as sfun
 import numpy as np
 import scipy.io as sio
+import bet.sample as sample
 
 # Import "Truth"
-mdat = sio.loadmat('Q_3D')
+mdat = sio.loadmat('../matfiles/Q_3D')
 Q = mdat['Q']
 Q_ref = mdat['Q_true']
 
@@ -14,6 +15,11 @@
 samples = mdat['points'].transpose()
 lam_domain = np.array([[-900, 1200], [0.07, .15], [0.1, 0.2]])
 
+# Create input, output, and discretization from data read from file
+points = mdat['points']
+input_sample_set = sample.sample_set(points.shape[0])
+input_sample_set.set_values(points.transpose())
+input_sample_set.set_domain(lam_domain)
 print "Finished loading data"
 
 def postprocess(station_nums, ref_num):
@@ -24,30 +30,26 @@ def postprocess(station_nums, ref_num):
     filename += '_ref_'+str(ref_num+1)
 
     data = Q[:, station_nums]
+    output_sample_set = sample.sample_set(data.shape[1])
+    output_sample_set.set_values(data)
     q_ref = Q_ref[ref_num, station_nums]
 
     # Create Simple function approximation
     # Save points used to parition D for simple function approximation and the
     # approximation itself (this can be used to make close comparisions...)
-    (rho_D_M, d_distr_samples, d_Tree) = sfun.uniform_hyperrectangle(data,
-            q_ref, bin_ratio=0.15,
+    output_probability_set = sfun.regular_partition_uniform_distribution_rectangle_scaled(\
+            output_sample_set, q_ref, rect_scale=0.15,
             center_pts_per_edge=np.ones((data.shape[1],)))
-    mdict = dict()
-    mdict['rho_D_M'] = rho_D_M
-    mdict['d_distr_samples'] = d_distr_samples
+
+    my_disc = sample.discretization(input_sample_set, output_sample_set,
+            output_probability_set)
 
     # Calclate P on the actual samples with assumption that voronoi cells have
     # equal size
-    (P1, lam_vol1, io_ptr1) = calcP.prob(samples, data, rho_D_M,
-            d_distr_samples, d_Tree)
+    input_sample_set.estimate_volume_mc()
     print "Calculating prob"
-    mdict['P1'] = P1
-    mdict['lam_vol1'] = lam_vol1
-    mdict['lem1'] = samples
-    mdict['io_ptr1'] = io_ptr1
-
-    # Export P and compare to MATLAB solution visually
-    sio.savemat(filename, mdict, do_compression=True)
+    calcP.prob(my_disc)
+    sample.save_discretization(my_disc, filename, "prob_solution")
 
 # Post-process and save P and emulated points
 ref_num = 14
diff --git a/examples/fromFile_ADCIRCMap/fromFile2D.py b/examples/fromFile_ADCIRCMap/fromFile2D.py
index c3b6e628..fdfa9797 100644
--- a/examples/fromFile_ADCIRCMap/fromFile2D.py
+++ b/examples/fromFile_ADCIRCMap/fromFile2D.py
@@ -12,8 +12,6 @@
 
 # Set minima and maxima
 lam_domain = np.array([[.07, .15], [.1, .2]])
-param_min = lam_domain[:, 0]
-param_max = lam_domain[:, 1]
 
 # Select only the stations I care about this will lead to better sampling
 station_nums = [0, 5] # 1, 6
@@ -22,7 +20,7 @@
 transition_set = asam.transition_set(.5, .5**5, 1.0)
 
 # Read in Q_ref and Q to create the appropriate rho_D 
-mdat = sio.loadmat('Q_2D')
+mdat = sio.loadmat('../matfiles/Q_2D')
 Q = mdat['Q']
 Q = Q[:, station_nums]
 Q_ref = mdat['Q_true']
@@ -60,7 +58,7 @@ def rho_D(outputs):
 
 # Get samples
 inital_sample_type = "lhs"
-(samples, data, all_step_ratios) = sampler.generalized_chains(param_min, param_max,
+(my_disc, data, all_step_ratios) = sampler.generalized_chains(lam_domain,
         transition_set, kernel_rD, sample_save_file, inital_sample_type)
 
 # Read in points_ref and plot results
diff --git a/examples/fromFile_ADCIRCMap/fromFile3D.py b/examples/fromFile_ADCIRCMap/fromFile3D.py
index c662b246..e5ae925f 100644
--- a/examples/fromFile_ADCIRCMap/fromFile3D.py
+++ b/examples/fromFile_ADCIRCMap/fromFile3D.py
@@ -18,8 +18,6 @@
 ymax = 1500
 wall_height = -2.5
 
-param_min = param_domain[:, 0]
-param_max = param_domain[:, 1]
 
 # Select only the stations I care about this will lead to better
 # sampling
@@ -29,7 +27,7 @@
 transition_set = asam.transition_set(.5, .5**5, 0.5)
 
 # Read in Q_ref and Q to create the appropriate rho_D 
-mdat = sio.loadmat('Q_3D')
+mdat = sio.loadmat('../matfiles/Q_3D')
 Q = mdat['Q']
 Q = Q[:, station_nums]
 Q_ref = mdat['Q_true']
@@ -67,7 +65,7 @@ def rho_D(outputs):
 
 # Get samples
 inital_sample_type = "lhs"
-(samples, data, all_step_ratios) = sampler.generalized_chains(param_min, param_max,
+(my_disc, all_step_ratios) = sampler.generalized_chains(param_domain,
         transition_set, kernel_rD, sample_save_file, inital_sample_type)
 
 # Read in points_ref and plot results
diff --git a/examples/fromFile_ADCIRCMap/plotDomains2D.py b/examples/fromFile_ADCIRCMap/plotDomains2D.py
index 5dbf13b4..71508668 100644
--- a/examples/fromFile_ADCIRCMap/plotDomains2D.py
+++ b/examples/fromFile_ADCIRCMap/plotDomains2D.py
@@ -6,17 +6,16 @@
 import numpy as np
 import bet.postProcess.plotDomains as pDom
 import scipy.io as sio
+import bet.sample as sample
 
 # Set minima and maxima
 lam_domain = np.array([[.07, .15], [.1, .2]])
-param_min = lam_domain[:, 0]
-param_max = lam_domain[:, 1]
 
 # Select only the stations I care about this will lead to better sampling
 station_nums = [0, 5] # 1, 6
 
 # Read in Q_ref and Q to create the appropriate rho_D 
-mdat = sio.loadmat('Q_2D')
+mdat = sio.loadmat('../matfiles/Q_2D.mat')
 Q = mdat['Q']
 Q = Q[:, station_nums]
 Q_ref = mdat['Q_true']
@@ -39,15 +38,24 @@ def rho_D(outputs):
 p_ref = mdat['points_true']
 p_ref = p_ref[5:7, 15]
 
+# Create input, output, and discretization from data read from file
+points = mdat['points']
+input_sample_set = sample.sample_set(points.shape[0])
+input_sample_set.set_values(points.transpose())
+input_sample_set.set_domain(lam_domain)
+output_sample_set = sample.sample_set(Q.shape[1])
+output_sample_set.set_values(Q)
+my_disc = sample.discretization(input_sample_set, output_sample_set)
+
 # Show the samples in the parameter space
-pDom.show_param(samples=points.transpose(), data=Q, rho_D=rho_D, p_ref=p_ref)
+pDom.show_param(my_disc, rho_D=rho_D, p_ref=p_ref)
 # Show the corresponding samples in the data space
-pDom.show_data(data=Q, rho_D=rho_D, Q_ref=Q_ref)
+pDom.show_data(output_sample_set, rho_D=rho_D, Q_ref=Q_ref)
 # Show the data domain that corresponds with the convex hull of samples in the
 # parameter space
-pDom.show_data_domain_2D(samples=points.transpose(), data=Q, Q_ref=Q_ref)
+pDom.show_data_domain_2D(my_disc, Q_ref=Q_ref)
 
 # Show multiple data domains that correspond with the convex hull of samples in
 # the parameter space
-pDom.show_data_domain_multi(samples=points.transpose(), data=mdat['Q'],
-        Q_ref=mdat['Q_true'][15], Q_nums=[1,2,5], showdim='all')
+pDom.show_data_domain_multi(my_disc, Q_ref=mdat['Q_true'][15], Q_nums=[1,2,5],
+        showdim='all')
diff --git a/examples/fromFile_ADCIRCMap/plotDomains3D.py b/examples/fromFile_ADCIRCMap/plotDomains3D.py
index b9bdec86..94d8771c 100644
--- a/examples/fromFile_ADCIRCMap/plotDomains3D.py
+++ b/examples/fromFile_ADCIRCMap/plotDomains3D.py
@@ -7,6 +7,7 @@
 import bet.postProcess.plotDomains as pDom
 import scipy.io as sio
 from scipy.interpolate import griddata
+import bet.sample as sample
 
 # Set minima and maxima
 param_domain = np.array([[-900, 1500], [.07, .15], [.1, .2]])
@@ -15,15 +16,13 @@
 xmax = 1580
 ymax = 1500
 
-param_min = param_domain[:, 0]
-param_max = param_domain[:, 1]
 
 # Select only the stations I care about this will lead to better
 # sampling
 station_nums = [0, 4, 1] # 1, 5, 2
 
 # Read in Q_ref and Q to create the appropriate rho_D 
-mdat = sio.loadmat('Q_3D')
+mdat = sio.loadmat('../matfiles/Q_3D')
 Q = mdat['Q']
 Q = Q[:, station_nums]
 Q_ref = mdat['Q_true']
@@ -48,12 +47,21 @@ def rho_D(outputs):
 p_ref = mdat['points_true']
 p_ref = p_ref[:, 14]
 
+
+# Create input, output, and discretization from data read from file
+input_sample_set = sample.sample_set(points.shape[0])
+input_sample_set.set_values(points.transpose())
+input_sample_set.set_domain(param_domain)
+output_sample_set = sample.sample_set(Q.shape[1])
+output_sample_set.set_values(Q)
+my_disc = sample.discretization(input_sample_set, output_sample_set)
+
 # Show the samples in the parameter space
-pDom.show_param(samples=points.transpose(), data=Q, rho_D=rho_D, p_ref=p_ref)
+pDom.show_param(my_disc, rho_D=rho_D, p_ref=p_ref)
 # Show the corresponding samples in the data space
-pDom.show_data(data=Q, rho_D=rho_D, Q_ref=Q_ref)
+pDom.show_data(output_sample_set, rho_D=rho_D, Q_ref=Q_ref)
 
 # Show multiple data domains that correspond with the convex hull of samples in
 # the parameter space
-pDom.show_data_domain_multi(samples=points.transpose(), data=mdat['Q'],
-        Q_ref=mdat['Q_true'][15], Q_nums=[1,2,5], showdim='all')       
+pDom.show_data_domain_multi(my_disc, Q_ref=mdat['Q_true'][15], Q_nums=[1,2,5],
+        showdim='all')       
diff --git a/examples/fromFile_ADCIRCMap/sandbox_test_2D.py b/examples/fromFile_ADCIRCMap/sandbox_test_2D.py
index 1c097e0a..9d65d91f 100644
--- a/examples/fromFile_ADCIRCMap/sandbox_test_2D.py
+++ b/examples/fromFile_ADCIRCMap/sandbox_test_2D.py
@@ -21,9 +21,6 @@
 ymax = 1500
 wall_height = -2.5
 
-param_min = lam_domain[:, 0]
-param_max = lam_domain[:, 1]
-
 
 # Select only the stations I care about this will lead to better sampling
 station_nums = [0, 5] # 1, 6
@@ -33,7 +30,7 @@
 transition_set = asam.transition_set(.5, .5**5, 1.0)
 
 # Read in Q_ref and Q to create the appropriate rho_D 
-mdat = sio.loadmat('Q_2D')
+mdat = sio.loadmat('../matfiles/Q_2D')
 Q = mdat['Q']
 Q = Q[:, station_nums]
 Q_ref = mdat['Q_true']
@@ -75,24 +72,24 @@ def rho_D(outputs):
 
 # Get samples
 # Run with varying kernels
-gen_results = sampler.run_gen(kern_list, rho_D, maximum, param_min,
-        param_max, transition_set, sample_save_file)
-#run_reseed_results = sampler.run_gen(kern_list, rho_D, maximum, param_min,
-#        param_max, t_kernel, sample_save_file, reseed=3)
+gen_results = sampler.run_gen(kern_list, rho_D, maximum, lam_domain,
+        transition_set, sample_save_file) 
+#run_reseed_results = sampler.run_gen(kern_list, rho_D, maximum, lam_domain,
+#       t_kernel, sample_save_file, reseed=3)
 
 # Run with varying transition sets bounds
 init_ratio = [0.1, 0.25, 0.5]
 min_ratio = [2e-3, 2e-5, 2e-8]
 max_ratio = [.5, .75, 1.0]
 tk_results = sampler.run_tk(init_ratio, min_ratio, max_ratio, rho_D,
-        maximum, param_min, param_max, kernel_rD, sample_save_file)
+        maximum, lam_domain, kernel_rD, sample_save_file)
 
 # Run with varying increase/decrease ratios and tolerances for a rhoD_kernel
 increase = [1.0, 2.0, 4.0]
 decrease = [0.5, 0.5e2, 0.5e3]
 tolerance = [1e-4, 1e-6, 1e-8]
 incdec_results = sampler.run_inc_dec(increase, decrease, tolerance, rho_D,
-        maximum, param_min, param_max, transition_set, sample_save_file)
+        maximum, lam_domain, transition_set, sample_save_file)
 
 # Compare the quality of several sets of samples
 print "Compare yield of sample sets with various kernels"
diff --git a/examples/fromFile_ADCIRCMap/sandbox_test_3D.py b/examples/fromFile_ADCIRCMap/sandbox_test_3D.py
index 61a7de2f..75fd0fcb 100644
--- a/examples/fromFile_ADCIRCMap/sandbox_test_3D.py
+++ b/examples/fromFile_ADCIRCMap/sandbox_test_3D.py
@@ -21,8 +21,6 @@
 ymax = 1500
 wall_height = -2.5
 
-param_min = param_domain[:, 0]
-param_max = param_domain[:, 1]
 
 # Select only the stations I care about this will lead to better sampling
 station_nums = [0, 4, 1] # 1, 5, 2
@@ -31,7 +29,7 @@
 transition_set = asam.transition_set(.5, .5**5, 0.5)
 
 # Read in Q_ref and Q to create the appropriate rho_D 
-mdat = sio.loadmat('Q_3D')
+mdat = sio.loadmat('../matfiles/Q_3D')
 Q = mdat['Q']
 Q = Q[:, station_nums]
 Q_ref = mdat['Q_true']
@@ -73,24 +71,24 @@ def rho_D(outputs):
 
 # Get samples
 # Run with varying kernels
-gen_results = sampler.run_gen(heur_list, rho_D, maximum, param_min,
-        param_max, transition_set, sample_save_file)
-#run_reseed_results = sampler.run_gen(heur_list, rho_D, maximum, param_min,
-#        param_max, t_kernel, sample_save_file, reseed=3)
+gen_results = sampler.run_gen(heur_list, rho_D, maximum, param_domain,
+        transition_set, sample_save_file) 
+#run_reseed_results = sampler.run_gen(heur_list, rho_D, maximum, param_domain,
+#        t_kernel, sample_save_file, reseed=3)
 
 # Run with varying transition sets bounds
 init_ratio = [0.1, 0.25, 0.5]
 min_ratio = [2e-3, 2e-5, 2e-8]
 max_ratio = [.5, .75, 1.0]
 tk_results = sampler.run_tk(init_ratio, min_ratio, max_ratio, rho_D,
-        maximum, param_min, param_max, kernel_rD, sample_save_file)
+        maximum, param_domain, kernel_rD, sample_save_file)
 
 # Run with varying increase/decrease ratios and tolerances for a rhoD_kernel
 increase = [1.0, 2.0, 4.0]
 decrease = [0.5, 0.5e2, 0.5e3]
 tolerance = [1e-4, 1e-6, 1e-8]
 incdec_results = sampler.run_inc_dec(increase, decrease, tolerance, rho_D,
-        maximum, param_min, param_max, transition_set, sample_save_file)
+        maximum, param_domain, transition_set, sample_save_file)
 
 # Compare the quality of several sets of samples
 result_list = [gen_results, tk_results, incdec_results]

From 05794feba8075f177ecf38717b69ecbb97b6b700 Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 15:50:49 -0400
Subject: [PATCH 4/9] updated to so that emulated samples need never be
 globalized even in pserial

---
 bet/calculateP/calculateP.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/bet/calculateP/calculateP.py b/bet/calculateP/calculateP.py
index c3a47420..61b93214 100644
--- a/bet/calculateP/calculateP.py
+++ b/bet/calculateP/calculateP.py
@@ -159,7 +159,10 @@ def prob_mc(discretization):
     cvol = np.copy(vol)
     comm.Allreduce([vol, MPI.DOUBLE], [cvol, MPI.DOUBLE], op=MPI.SUM)
     vol = cvol
-    vol = vol/float(discretization._emulated_input_sample_set._values.shape[0])
+    num_l_emulate = discretization._emulated_input_sample_set.\
+            _values_local.shape[0]
+    num_l_emulate = comm.allreduce(num_l_emulate, op=MPI.SUM)
+    vol = vol/float(num_l_emulate)
     discretization._input_sample_set._volumes = vol
     discretization._input_sample_set.global_to_local()
 

From a6ce4ee109f8ff33bde9bd47db881cf861e32cd5 Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 15:51:17 -0400
Subject: [PATCH 5/9] lowered num_l_emulate for serial examples

---
 examples/fromFile_ADCIRCMap/Q_1D_serial.py | 2 +-
 examples/fromFile_ADCIRCMap/Q_2D_serial.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/fromFile_ADCIRCMap/Q_1D_serial.py b/examples/fromFile_ADCIRCMap/Q_1D_serial.py
index 8857f309..6447f465 100644
--- a/examples/fromFile_ADCIRCMap/Q_1D_serial.py
+++ b/examples/fromFile_ADCIRCMap/Q_1D_serial.py
@@ -42,7 +42,7 @@ def postprocess(station_nums, ref_num):
             output_sample_set, q_ref, rect_scale=0.15,
             center_pts_per_edge=np.ones((data.shape[1],)))
 
-    num_l_emulate = 1e6
+    num_l_emulate = 1e4
     set_emulated = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
     my_disc = sample.discretization(input_sample_set, output_sample_set,
             output_probability_set, emulated_input_sample_set=set_emulated)
diff --git a/examples/fromFile_ADCIRCMap/Q_2D_serial.py b/examples/fromFile_ADCIRCMap/Q_2D_serial.py
index 86862a11..2bb1a987 100644
--- a/examples/fromFile_ADCIRCMap/Q_2D_serial.py
+++ b/examples/fromFile_ADCIRCMap/Q_2D_serial.py
@@ -40,7 +40,7 @@ def postprocess(station_nums, ref_num):
             output_sample_set, q_ref, rect_scale=0.15,
             center_pts_per_edge=np.ones((data.shape[1],)))
 
-    num_l_emulate = 1e6
+    num_l_emulate = 1e4
     set_emulated = calcP.emulate_iid_lebesgue(lam_domain, num_l_emulate)
     my_disc = sample.discretization(input_sample_set, output_sample_set,
             output_probability_set, emulated_input_sample_set=set_emulated)

From 331d40782f6207718600b4fe0843ae9e63dc63ed Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 16:43:29 -0400
Subject: [PATCH 6/9] addressing issue #199

---
 bet/calculateP/calculateP.py        |  6 ++--
 bet/postProcess/plotP.py            | 50 ++++++++++++++++++++++++-----
 test/test_postProcess/test_plotP.py | 38 +++++++++++++---------
 3 files changed, 66 insertions(+), 28 deletions(-)

diff --git a/bet/calculateP/calculateP.py b/bet/calculateP/calculateP.py
index 61b93214..63bea40d 100644
--- a/bet/calculateP/calculateP.py
+++ b/bet/calculateP/calculateP.py
@@ -50,10 +50,6 @@ def prob_emulated(discretization, globalize=True):
     ``num_l_emulate`` iid samples :math:`(\lambda_{emulate})`.
     This is added to the emulated input sample set object.
 
-    .. todo::
-        
-        @smattis the way this is written globalize does nothing
-
     :param discretization: An object containing the discretization information.
     :type class:`bet.sample.discretization`
     :param bool globalize: Makes local variables global.
@@ -86,6 +82,8 @@ def prob_emulated(discretization, globalize=True):
                         _probabilities[i]/Itemp_sum
     
     discretization._emulated_input_sample_set._probabilities_local = P
+    if globalize:
+        discretization._emulated_input_sample_set.local_to_global()
     pass
 
 def prob(discretization): 
diff --git a/bet/postProcess/plotP.py b/bet/postProcess/plotP.py
index e3b804a7..d1dd39ff 100644
--- a/bet/postProcess/plotP.py
+++ b/bet/postProcess/plotP.py
@@ -23,13 +23,23 @@ class bad_object(Exception):
     Exception for when the wrong type of object is used.
     """
 
+class missing_attribute(Exception):
+    """
+    Exception for missing attribute.
+    """
+
 def calculate_1D_marginal_probs(sample_set, nbins=20):
         
     r"""
     This calculates every single marginal of the probability measure
     described by the probabilities within the sample_set object.
     If the sample_set object is a discretization object, we assume
-    that the probabilities to be plotted are from the input space.
+    that the probabilities to be plotted are from the input space on the
+    emulated samples
+    (``discretization._emulated_input_sample_set._probabilties_local``).
+
+    This assumes that the user has already run
+    :meth:`~bet.calculateP.calculateP.prob_emulated`.
 
     :param sample_set: Object containing samples and probabilities
     :type sample_set: :class:`~bet.sample.sample_set_base` or 
@@ -41,12 +51,21 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
 
     """
     if isinstance(sample_set, sample.discretization):
-        sample_obj = sample_set._input_sample_set
+        sample_obj = sample_set._emulated_input_sample_set
+        if sample_obj is None:
+            raise missing_attribute("Missing emulated_input_sample_set")
     elif isinstance(sample_set, sample.sample_set_base):
         sample_obj = sample_set
     else:
         raise bad_object("Improper sample object")
 
+    # Check for local probabilities
+    if sample_obj.probabilities_local is None:
+        if sample_obj.probabilities is None:
+            raise missing_attribute("Missing probabilities")
+        else:
+            sample_obj.global_to_local()
+
     # Make list of bins if only an integer is given
     if isinstance(nbins, int):
         nbins = nbins*np.ones(sample_obj.get_dim(), dtype=np.int)
@@ -61,8 +80,8 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
     # Calculate marginals
     marginals = {}
     for i in range(sample_obj.get_dim()):
-        [marg, _] = np.histogram(sample_obj.get_values()[:, i], bins=bins[i],
-                                 weights=sample_obj.get_probabilities())
+        [marg, _] = np.histogram(sample_obj.get_values_local()[:, i], bins=bins[i],
+                                 weights=sample_obj.get_probabilities_local())
         marg_temp = np.copy(marg)
         comm.Allreduce([marg, MPI.DOUBLE], [marg_temp, MPI.DOUBLE], op=MPI.SUM)
         marginals[i] = marg_temp
@@ -75,7 +94,13 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
     This calculates every pair of marginals (or joint in 2d case) of
     input probability measure defined on a rectangular grid.
     If the sample_set object is a discretization object, we assume
-    that the probabilities to be plotted are from the input space.
+    that the probabilities to be plotted are from the input space on the
+    emulated samples
+    (``discretization._emulated_input_sample_set._probabilties_local``).
+
+    This assumes that the user has already run
+    :meth:`~bet.calculateP.calculateP.prob_emulated`.
+
 
     :param sample_set: Object containing samples and probabilities
     :type sample_set: :class:`~bet.sample.sample_set_base` 
@@ -87,12 +112,21 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
 
     """
     if isinstance(sample_set, sample.discretization):
-        sample_obj = sample_set._input_sample_set
+        sample_obj = sample_set._emulated_input_sample_set
+        if sample_obj is None:
+            raise missing_attribute("Missing emulated_input_sample_set")
     elif isinstance(sample_set, sample.sample_set_base):
         sample_obj = sample_set
     else:
         raise bad_object("Improper sample object")
 
+    # Check for local probabilities
+    if sample_obj.probabilities_local is None:
+        if sample_obj.probabilities is None:
+            raise missing_attribute("Missing probabilities")
+        else:
+            sample_obj.global_to_local()
+
     if sample_obj.get_dim() < 2:
         raise dim_not_matching("Incompatible dimensions of sample set"
                                " for plotting")
@@ -112,9 +146,9 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
     marginals = {}
     for i in range(sample_obj.get_dim()):
         for j in range(i+1, sample_obj.get_dim()):
-            (marg, _) = np.histogramdd(sample_obj.get_values()[:, [i, j]],
+            (marg, _) = np.histogramdd(sample_obj.get_values_local()[:, [i, j]],
                     bins=[bins[i], bins[j]],
-                    weights=sample_obj.get_probabilities())
+                    weights=sample_obj.get_probabilities_local())
             marg = np.ascontiguousarray(marg)
             marg_temp = np.copy(marg)
             comm.Allreduce([marg, MPI.DOUBLE], [marg_temp, MPI.DOUBLE],
diff --git a/test/test_postProcess/test_plotP.py b/test/test_postProcess/test_plotP.py
index aa138f7c..d29c31df 100644
--- a/test/test_postProcess/test_plotP.py
+++ b/test/test_postProcess/test_plotP.py
@@ -29,21 +29,22 @@ def setUp(self):
         """
         Set up problem.
         """
-        input_samples = sample.sample_set(1)
-        input_samples.set_domain(np.array([[0.0, 1.0]]))
+        emulated_input_samples = sample.sample_set(1)
+        emulated_input_samples.set_domain(np.array([[0.0, 1.0]]))
 
         num_samples=1000
 
-        input_samples.set_values(np.linspace(input_samples.get_domain()[0][0],
-                                             input_samples.get_domain()[0][1],
+        emulated_input_samples.set_values_local(np.linspace(emulated_input_samples.get_domain()[0][0],
+                                             emulated_input_samples.get_domain()[0][1],
                                              num_samples+1))
 
-        input_samples.set_probabilities(1.0/float(comm.size)*(1.0/float(input_samples.get_values().shape[0]))
-                                        *np.ones((input_samples.get_values().shape[0],)))
+        emulated_input_samples.set_probabilities_local(1.0/float(comm.size)*(1.0/float(\
+                emulated_input_samples.get_values_local().shape[0]))\
+                *np.ones((emulated_input_samples.get_values_local().shape[0],)))
 
-        input_samples.check_num()
+        emulated_input_samples.check_num()
 
-        self.samples = input_samples
+        self.samples = emulated_input_samples
 
     def test_1_bin(self):
         """
@@ -67,23 +68,28 @@ def test_10_bins(self):
 
 class Test_calc_marg_2D(unittest.TestCase):
     """
-    Test :meth:`bet.postProcess.plotP.calculate_1D_marginal_probs` and  :meth:`bet.postProcess.plotP.calculate_2D_marginal_probs` for a 2D
+    Test :meth:`bet.postProcess.plotP.calculate_1D_marginal_probs` and
+    :meth:`bet.postProcess.plotP.calculate_2D_marginal_probs` for a 2D
     parameter space.
     """
     def setUp(self):
         """
         Set up problem.
         """
-        input_samples = sample.sample_set(2)
-        input_samples.set_domain(np.array([[0.0,1.0],[0.0,1.0]]))
+        emulated_input_samples = sample.sample_set(2)
+        emulated_input_samples.set_domain(np.array([[0.0,1.0],[0.0,1.0]]))
 
-        input_samples.set_values(util.meshgrid_ndim((np.linspace(input_samples.get_domain()[0][0], input_samples.get_domain()[0][1], 10),
-                                                     np.linspace(input_samples.get_domain()[1][0], input_samples.get_domain()[1][1], 10))))
+        emulated_input_samples.set_values_local(util.meshgrid_ndim((np.linspace(emulated_input_samples.get_domain()[0][0],
+            emulated_input_samples.get_domain()[0][1], 10),
+            np.linspace(emulated_input_samples.get_domain()[1][0],
+                emulated_input_samples.get_domain()[1][1], 10))))
 
-        input_samples.set_probabilities(1.0/float(comm.size)*(1.0/float(input_samples.get_values().shape[0]))*np.ones((input_samples.get_values().shape[0],)))
-        input_samples.check_num()
+        emulated_input_samples.set_probabilities_local(1.0/float(comm.size)*\
+                (1.0/float(emulated_input_samples.get_values_local().shape[0]))*\
+                np.ones((emulated_input_samples.get_values_local().shape[0],)))
+        emulated_input_samples.check_num()
 
-        self.samples = input_samples
+        self.samples = emulated_input_samples
 
     def test_1_bin_1D(self):
         """ 

From d22322463641250fd54f6f61b679135c3e4a16fa Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 16:50:17 -0400
Subject: [PATCH 7/9] fixed VisibleDeprecationWarning by recasting
 num_l_emulate to int

---
 bet/calculateP/calculateP.py               | 4 ++--
 examples/fromFile_ADCIRCMap/Q_2D_serial.py | 9 +++------
 2 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/bet/calculateP/calculateP.py b/bet/calculateP/calculateP.py
index 63bea40d..dd96a1ad 100644
--- a/bet/calculateP/calculateP.py
+++ b/bet/calculateP/calculateP.py
@@ -29,8 +29,8 @@ def emulate_iid_lebesgue(domain, num_l_emulate, globalize=False):
     :returns: a set of samples for emulation
 
     """
-    num_l_emulate = (num_l_emulate/comm.size) + \
-            (comm.rank < num_l_emulate%comm.size)
+    num_l_emulate = int((num_l_emulate/comm.size) + \
+            (comm.rank < num_l_emulate%comm.size))
     lam_width = domain[:, 1] - domain[:, 0]
     lambda_emulate = lam_width*np.random.random((num_l_emulate,
         domain.shape[0]))+domain[:, 0]
diff --git a/examples/fromFile_ADCIRCMap/Q_2D_serial.py b/examples/fromFile_ADCIRCMap/Q_2D_serial.py
index 2bb1a987..85e4eab0 100644
--- a/examples/fromFile_ADCIRCMap/Q_2D_serial.py
+++ b/examples/fromFile_ADCIRCMap/Q_2D_serial.py
@@ -50,23 +50,20 @@ def postprocess(station_nums, ref_num):
     # Calculate P on lambda emulate
     print "Calculating prob_emulated"
     calcP.prob_emulated(my_disc)
-    if comm.rank == 0:
-        sample.save_discretization(my_disc, filename, "prob_emulated_solution")
+    sample.save_discretization(my_disc, filename, "prob_emulated_solution")
 
     # Calclate P on the actual samples with assumption that voronoi cells have
     # equal size
     input_sample_set.estimate_volume_mc()
     print "Calculating prob"
     calcP.prob(my_disc)
-    if comm.rank == 0:
-        sample.save_discretization(my_disc, filename, "prob_solution")
+    sample.save_discretization(my_disc, filename, "prob_solution")
 
     # Calculate P on the actual samples estimating voronoi cell volume with MC
     # integration
     calcP.prob_mc(my_disc)
     print "Calculating prob_mc"
-    if comm.rank == 0:
-        sample.save_discretization(my_disc, filename, "prob_mc_solution")
+    sample.save_discretization(my_disc, filename, "prob_mc_solution")
 
 # Post-process and save P and emulated points
 ref_nums = [6, 11, 15] # 7, 12, 16

From 73ecb28a5d38a4443dd321334eb84c20b2aa9083 Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 17:07:25 -0400
Subject: [PATCH 8/9] fixes name error, and example

---
 bet/postProcess/plotP.py                  | 4 ++--
 examples/fromFile_ADCIRCMap/fromFile2D.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/bet/postProcess/plotP.py b/bet/postProcess/plotP.py
index d1dd39ff..3883e88f 100644
--- a/bet/postProcess/plotP.py
+++ b/bet/postProcess/plotP.py
@@ -60,7 +60,7 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
         raise bad_object("Improper sample object")
 
     # Check for local probabilities
-    if sample_obj.probabilities_local is None:
+    if sample_obj._probabilities_local is None:
         if sample_obj.probabilities is None:
             raise missing_attribute("Missing probabilities")
         else:
@@ -121,7 +121,7 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
         raise bad_object("Improper sample object")
 
     # Check for local probabilities
-    if sample_obj.probabilities_local is None:
+    if sample_obj._probabilities_local is None:
         if sample_obj.probabilities is None:
             raise missing_attribute("Missing probabilities")
         else:
diff --git a/examples/fromFile_ADCIRCMap/fromFile2D.py b/examples/fromFile_ADCIRCMap/fromFile2D.py
index fdfa9797..180a2e18 100644
--- a/examples/fromFile_ADCIRCMap/fromFile2D.py
+++ b/examples/fromFile_ADCIRCMap/fromFile2D.py
@@ -58,7 +58,7 @@ def rho_D(outputs):
 
 # Get samples
 inital_sample_type = "lhs"
-(my_disc, data, all_step_ratios) = sampler.generalized_chains(lam_domain,
+(my_disc,  all_step_ratios) = sampler.generalized_chains(lam_domain,
         transition_set, kernel_rD, sample_save_file, inital_sample_type)
 
 # Read in points_ref and plot results

From 8a31ef6ed00a00e12f779d244715714ae455dc9d Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Wed, 25 May 2016 17:28:35 -0400
Subject: [PATCH 9/9] all examples in fromFile_ADCIRCMap have been tested
 except for Q_2D_parallel, and sandbox_test_[23]D.py

---
 examples/fromFile_ADCIRCMap/plotDomains2D.py | 2 +-
 examples/fromFile_ADCIRCMap/plotDomains3D.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/fromFile_ADCIRCMap/plotDomains2D.py b/examples/fromFile_ADCIRCMap/plotDomains2D.py
index 71508668..e4a00e63 100644
--- a/examples/fromFile_ADCIRCMap/plotDomains2D.py
+++ b/examples/fromFile_ADCIRCMap/plotDomains2D.py
@@ -57,5 +57,5 @@ def rho_D(outputs):
 
 # Show multiple data domains that correspond with the convex hull of samples in
 # the parameter space
-pDom.show_data_domain_multi(my_disc, Q_ref=mdat['Q_true'][15], Q_nums=[1,2,5],
+pDom.show_data_domain_multi(my_disc, Q_ref=mdat['Q_true'][15], 
         showdim='all')
diff --git a/examples/fromFile_ADCIRCMap/plotDomains3D.py b/examples/fromFile_ADCIRCMap/plotDomains3D.py
index 94d8771c..2b15b5ac 100644
--- a/examples/fromFile_ADCIRCMap/plotDomains3D.py
+++ b/examples/fromFile_ADCIRCMap/plotDomains3D.py
@@ -63,5 +63,5 @@ def rho_D(outputs):
 
 # Show multiple data domains that correspond with the convex hull of samples in
 # the parameter space
-pDom.show_data_domain_multi(my_disc, Q_ref=mdat['Q_true'][15], Q_nums=[1,2,5],
+pDom.show_data_domain_multi(my_disc, Q_ref=mdat['Q_true'][15],
         showdim='all')