fixed formatting conflict from merge

bet/postProcess/plotP.py

@@ -2,42 +2,11 @@
 This module provides methods for plotting probabilities. 
 """
 
-# NOTE: We need nbins+1 because of plotting considerations.
-""" 
-TODO: Make sure plots are transparent when saving them. Copy or emulate code
-for saving the plots from plotDomains.py
-"""
-
 from bet.Comm import *
 import matplotlib.pyplot as plt
 import numpy as np
 import copy, math
 
-def plot_voronoi_probs(P_samples, samples, lam_domain, nbins=20,
-        plot_surface=False):
-    """
-    This makes plots of the joint probabilies of input probability measure
-    defined by P_samples for 2d cases. post_process - is an input that only
-    applies to the 2d case  w.r.t. the Voronoi cells.
-
-    :param P_samples: Probabilities.
-    :type P_samples: :class:'~numpy.ndarray' of shape (num_samples,)
-    :param samples: The samples in parameter space for which the model was run.
-    :type samples: :class:'~numpy.ndarray' of shape (num_samples, ndim)
-    :param lam_domain: The domain for each parameter for the model.
-    :type lam_domain: :class:'~numpy.ndarray' of shape (ndim, 2)
-    :param nbins: Number of bins in each direction.
-    :type nbins: :int
-
-    """
-    lam_dim = lam_domain.shape[0]
-
-    if lam_dim == 2: # Plot Voronoi tesselations, otherwise plot 2d 
-        #projections/marginals of the joint inverse measure
-        num_samples = samples.shape[0]
-        #Add fake samples outside of lam_domain to close Voronoi 
-        #tesselations at infinity
-    # TODO: finish this method
 
 def calculate_1D_marginal_probs(P_samples, samples, lam_domain, nbins=20):
 
@@ -66,35 +35,18 @@ def calculate_1D_marginal_probs(P_samples, samples, lam_domain, nbins=20):
     if isinstance(nbins, int):
         nbins = nbins*np.ones(num_dim, dtype=np.int)
 
-    """TODO: determine if np.histogram or np.searchsorted is faster and change
-    to appropiate method. 
-    """
-
     # Create bins
     bins = []
     for i in range(num_dim):
-        bins.append(np.linspace(lam_domain[i][0]-np.finfo(float).eps, lam_domain[i][1]+np.finfo(float).eps,
-            nbins[i]+1))
-    bin_ptr = np.zeros((num_samples, num_dim), dtype=np.int)
-    # Bin samples
-    for j in range(num_dim):
-        bin_ptr[:, j] = np.searchsorted(bins[j], samples[:, j])
-    bin_ptr -= 1
-    bins = []
-    for i in range(num_dim):
-        bins.append(np.linspace(lam_domain[i][0], lam_domain[i][1],
-                                nbins[i]+1))
-
-    # Calculate marginal probabilities 
+        bins.append(np.linspace(lam_domain[i][0], lam_domain[i][1], nbins[i]+1))
+
+    # Calculate marginals
     marginals = {}
     for i in range(num_dim):
-        marg = np.zeros((nbins[i]+1,))
-        # This may be sped up with logical indices
-        for k in range(num_samples):
-            marg[bin_ptr[k][i]] += P_samples[k]
+        [marg, _] = np.histogram(samples[:,i], bins=bins[i], weights=P_samples)
         marg_temp = np.copy(marg)
         comm.Allreduce([marg, MPI.DOUBLE], [marg_temp, MPI.DOUBLE], op=MPI.SUM)
-        marginals[i] = marg_temp[:-1]
+        marginals[i] = marg_temp
 
     return (bins, marginals)
 
@@ -116,47 +68,35 @@ def calculate_2D_marginal_probs(P_samples, samples, lam_domain, nbins=20):
     :returns: (bins, marginals)
 
     """
-
+    if len(samples.shape) == 1:
+        samples = np.expand_dims(samples, axis=1)
     num_samples = samples.shape[0]
     num_dim = samples.shape[1]
 
     # Make list of bins if only an integer is given
     if isinstance(nbins, int):
         nbins = nbins*np.ones(num_dim, dtype=np.int)
-    """TODO: determine if np.histogram or np.searchsorted is faster and change
-    to appropiate method. 
-    """
 
     # Create bins
-    bins = []
-    for i in range(num_dim):
-        bins.append(np.linspace(lam_domain[i][0]-np.finfo(float).eps, lam_domain[i][1]+np.finfo(float).eps, nbins[i]+1))
-    bin_ptr = np.zeros((num_samples, num_dim), dtype=np.int)
-    # Bin samples
-    for j in range(num_dim):
-        bin_ptr[:, j] = np.searchsorted(bins[j], samples[:, j])
-    bin_ptr -= 1
-
     bins = []
     for i in range(num_dim):
         bins.append(np.linspace(lam_domain[i][0], lam_domain[i][1], nbins[i]+1))
-         
-    # Calculate marginal probabilities 
+
+    # Calculate marginals
     marginals = {}
     for i in range(num_dim):
         for j in range(i+1, num_dim):
-            marg = np.zeros((nbins[i]+1, nbins[j]+1))
-            # TODO: This may be sped up with logical indices
-            for k in range(num_samples):
-                marg[bin_ptr[k][i]][bin_ptr[k][j]] += P_samples[k]
+            (marg, _) = np.histogramdd(samples[:,[i,j]], bins = [bins[i], bins[j]], weights = P_samples)
+            marg=np.ascontiguousarray(marg)
             marg_temp = np.copy(marg)
             comm.Allreduce([marg, MPI.DOUBLE], [marg_temp, MPI.DOUBLE],
                     op=MPI.SUM) 
-            marginals[(i, j)] = marg_temp[:-1, :-1]
+            marginals[(i, j)] = marg_temp
+
     return (bins, marginals)
 
 def plot_1D_marginal_probs(marginals, bins, lam_domain,
-        filename="file", lam_ref=None, interactive=True,
+        filename="file", lam_ref=None, interactive=False,
         lambda_label=None):
 
     """
@@ -197,14 +137,15 @@ def plot_1D_marginal_probs(marginals, bins, lam_domain,
                 label1 = lambda_label[i]
             ax.set_xlabel(label1) 
             ax.set_ylabel(r'$\rho$')
-            fig.savefig(filename + "_1D_" + str(i) + ".eps")
+            fig.savefig(filename + "_1D_" + str(i) + ".eps", transparent=True)
             if interactive:
                 plt.show()
             else:
                 plt.close()
+            plt.clf()
 
 def plot_2D_marginal_probs(marginals, bins, lam_domain,
-        filename="file", lam_ref=None, plot_surface=False, interactive=True,
+        filename="file", lam_ref=None, plot_surface=False, interactive=False,
         lambda_label=None):
 
     """
@@ -260,7 +201,7 @@ def plot_2D_marginal_probs(marginals, bins, lam_domain,
             fig.colorbar(quadmesh, ax=ax, label=label_cbar)
             plt.axis([lam_domain[i][0], lam_domain[i][1], lam_domain[j][0],
                 lam_domain[j][1]]) 
-            fig.savefig(filename + "_2D_" + str(i) + "_" + str(j) + ".eps")
+            fig.savefig(filename + "_2D_" + str(i) + "_" + str(j) + ".eps", transparent=True)
             if interactive:
                 plt.show()
             else:
@@ -283,11 +224,12 @@ def plot_2D_marginal_probs(marginals, bins, lam_domain,
                 ax.set_zlabel(r'$P$')
                 plt.backgroundcolor = 'w'
                 fig.colorbar(surf, shrink=0.5, aspect=5, label=r'$P$')
-                fig.savefig(filename + "_surf_"+str(i)+"_"+str(j)+".eps")
+                fig.savefig(filename + "_surf_"+str(i)+"_"+str(j)+".eps", transparent=True)
                 if interactive:
                     plt.show()
                 else:
                     plt.close()
+                plt.clf()
 
 def smooth_marginals_1D(marginals, bins, sigma=10.0):
     """
@@ -351,8 +293,7 @@ def smooth_marginals_2D(marginals, bins, sigma=10.0):
     marginals_smooth = {}
     pairs = copy.deepcopy(marginals.keys())
     pairs.sort()
-    #TODO: k is an unused variable
-    for k, (i, j) in enumerate(pairs):   
+    for (i, j) in pairs:   
         nx = len(bins[i])-1
         ny = len(bins[j])-1
         dx = bins[i][1] - bins[i][0]

bet/postProcess/postTools.py

@@ -1,7 +1,5 @@
 """
-TODO: rewrite to match what the module acutally does
-
-This module provides methods for plotting probabilities. 
+This module provides methods for postprocessing probabilities and data. 
 """
 from bet.Comm import rank
 import numpy as np
@@ -70,8 +68,6 @@ def sample_highest_prob(top_percentile, P_samples, samples, lam_vol=None,
     """
     if len(samples.shape) == 1:
         samples = np.expand_dims(samples, axis=1)
-    # TODO: should we trust the user or should we do a quick check of part of
-    # the array to see if it's already sorted?
     if sort:
         (P_samples, samples, lam_vol, data) = sort_by_rho(P_samples, samples,
                 lam_vol, data)

examples/linearMap/linearMapUniformSampling.py

@@ -69,6 +69,8 @@
 
 # calc data
 data= np.dot(samples,Q_map)
+np.savetxt('3to2_samples.txt.gz', samples)
+np.savetxt('3to2_data.txt.gz', data)
 
 '''
 Suggested changes for user:
@@ -122,7 +124,7 @@
     Monte Carlo (it converges like 1/sqrt(n_samples)).
 '''
 #lambda_emulate = samples
-lambda_emulate = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, num_l_emulate = 1E5)
+lambda_emulate = calculateP.emulate_iid_lebesgue(lam_domain=lam_domain, num_l_emulate = 1E4)
 
 
 # calculate probablities

test/test_postProcess/test_plotP.py

@@ -30,7 +30,7 @@ def setUp(self):
         self.lam_domain=np.array([[0.0,1.0]])
         num_samples=1000
         self.samples = np.linspace(self.lam_domain[0][0], self.lam_domain[0][1], num_samples+1)
-        self.P_samples = (1.0/float(self.samples.shape[0]))*np.ones((self.samples.shape[0],))
+        self.P_samples = 1.0/float(size)*(1.0/float(self.samples.shape[0]))*np.ones((self.samples.shape[0],))
 
     def test_1_bin(self):
         """
@@ -65,7 +65,7 @@ def setUp(self):
         """
         self.lam_domain=np.array([[0.0,1.0],[0.0,1.0]])
         self.samples=util.meshgrid_ndim((np.linspace(self.lam_domain[0][0], self.lam_domain[0][1], 10),np.linspace(self.lam_domain[1][0], self.lam_domain[1][1], 10)))
-        self.P_samples = (1.0/float(self.samples.shape[0]))*np.ones((self.samples.shape[0],))
+        self.P_samples = 1.0/float(size)*(1.0/float(self.samples.shape[0]))*np.ones((self.samples.shape[0],))
 
     def test_1_bin_1D(self):
         """