diff --git a/bet/postProcess/plot_voronoi_probs.m b/bet/postProcess/plot_voronoi_probs.m
deleted file mode 100644
index cc557d62..00000000
--- a/bet/postProcess/plot_voronoi_probs.m
+++ /dev/null
@@ -1,187 +0,0 @@
-function plot_voronoi_probs(P_samples,samples,lam_dim,lam_domain,post_process, nbins)
-%
-% This makes plots of every pair of marginals (or joint in 2d case) of
-% input probability measure defined by P_samples
-%
-% post_process - is an input that only applies to the 2d case if you want
-% to see marginals on a regular grid instead of w.r.t. the Voronoi cells.
-%
-lam_domain = lam_domain';
-if nargin < 5
-    post_process = 0;
-end
-if nargin == 5
-    if isempty(post_process)
-        post_process = 0;
-    end
-end
-if nargin < 6
-    nbins = 20;  % number of bins in each direction
-end
-if nargin >= 5
-    if isempty(nbins)
-        nbins = 20;  % number of bins in each direction
-    end
-end
-
-if lam_dim == 2 % Plot Voronoi tesselations, otherwise plot 2d projections/marginals of the joint inverse measure
-    [num_samples,~] = size(samples);
-    fig_num = 5; % Starting figure number
-    
-    fprintf('Plot Voronoi tesselations \n')
-    % Add fake samples outside of lam_domain to close Voronoi tesselations at
-    % infinity
-    midpt = mean(lam_domain);
-    for i=1:lam_dim
-        new_pt = midpt;
-        new_pt(i) = lam_domain(1,i)-10;
-        samples = [samples; new_pt];
-        new_pt = midpt;
-        new_pt(i) = lam_domain(2,i)+10;
-        samples = [samples; new_pt];
-    end
-    
-    indices = nchoosek(1:lam_dim,2);
-    
-    if ~post_process
-        % Get convex hulls for Voronoi cells in each plane of input domain
-        vv = [];
-        cc = [];
-        for i=1:nchoosek(lam_dim,2)
-            [Vor(i).vertices,Vor(i).cells]=voronoin(samples(:,indices(i,:)));
-            %     fprintf('Length of cc = %6i \n', length(cc))
-            if length(Vor(i).cells) ~= num_samples+2*lam_dim
-                fprintf('plot_voronoi: Houston, we have a problem \n')
-                return
-            end
-        end
-        for i=1:nchoosek(lam_dim,2)
-            figure(fig_num);
-            clf
-            fig_num = fig_num+1;            
-            hold on
-            for l=1:length(Vor(i).cells)
-                if all (Vor(i).cells{l}-1 ~= 0)
-                    patch(Vor(i).vertices(Vor(i).cells{l},1),Vor(i).vertices(Vor(i).cells{l},2),P_samples(l));
-                end
-            end
-            axis([lam_domain(1,indices(i,1)) lam_domain(2,indices(i,1)) lam_domain(1,indices(i,2)) lam_domain(2,indices(i,2))])
-            %title('Voronoi Tesselation colored by probability', 'FontSize', 15)
-            s = ['$\lambda_' int2str(indices(i,1)) '$'];
-            xlabel(s, 'interpreter', 'latex', 'FontSize', 15)
-            s = ['$\lambda_' int2str(indices(i,2)) '$'];
-            ylabel(s, 'interpreter', 'latex', 'FontSize', 15)
-            colorbar
-            hold off
-        end
-    else
-        [num_samples,~] = size(samples);
-        % Init the data struct holding info used for marginals
-        Lambda_Info = struct('bin_vals',zeros(lam_dim,nbins+1),'bin_ptr',zeros(num_samples,lam_dim));
-        for k=1:lam_dim
-            % Compute bins in each direction
-            Lambda_Info.bin_vals(k,:) = linspace(lam_domain(1,k),lam_domain(2,k),nbins+1);
-            % Compute bin pointer for each sample
-            [~,Lambda_Info.bin_ptr(:,k)] = histc(samples(:,k),Lambda_Info.bin_vals(k,:));
-        end
-        indices = nchoosek(1:lam_dim,2);
-        Lambda_pairs = struct('pair',struct());
-        for i=1:nchoosek(lam_dim,2) % Compute marginals for each 2d pair of coordinates
-            Lambda_pairs.pair(i).terms = indices(i,:);
-            Lambda_pairs.pair(i).marginal = zeros(nbins+1,nbins+1);
-        end
-        
-        for i = 1:nchoosek(lam_dim,2) % Compute marginals
-            for j=1:num_samples
-                flag = 0;
-                k = 1;
-                while k<nbins+1 && flag == 0
-                    if Lambda_Info.bin_ptr(j,indices(i,1)) == k
-                        kk = 1;
-                        while kk<nbins+1 && flag == 0
-                            if Lambda_Info.bin_ptr(j,indices(i,2)) == kk
-                                Lambda_pairs.pair(i).marginal(kk+1,k+1) = Lambda_pairs.pair(i).marginal(kk+1,k+1) + P_samples(j);
-                                flag = 1;
-                            end
-                            kk = kk+1;
-                        end
-                    end
-                    k = k+1;
-                end
-            end
-        end
-        
-        fig_num = 10;
-        for i = 1:nchoosek(lam_dim,2)
-            figure(fig_num);
-            zz = [Lambda_pairs.pair(i).marginal(2:end, 2:end) zeros(nbins, 1); zeros(1,nbins+1)];
-            surf(Lambda_Info.bin_vals(indices(i,1),:),Lambda_Info.bin_vals(indices(i,2),:),zz)
-            shading flat
-            axis([lam_domain(1,indices(i,1)) lam_domain(2,indices(i,1)) lam_domain(1,indices(i,2)) lam_domain(2,indices(i,2))])
-            s = ['Marginal in $\lambda_' int2str(indices(i,1)) '$ x $\lambda_' int2str(indices(i,2)) '$'];
-            %title(s, 'interpreter', 'latex', 'FontSize', 15)
-            s = ['$\lambda_' int2str(indices(i,1)) '$'];
-            xlabel(s, 'interpreter', 'latex', 'FontSize', 15)
-            s = ['$\lambda_' int2str(indices(i,2)) '$'];
-            ylabel(s, 'interpreter', 'latex', 'FontSize', 15)
-            colorbar
-            fig_num = fig_num+1;
-        end
-    end
-    
-else % Higher dimension case
-    [num_samples,~] = size(samples);
-    % Init the data struct holding info used for marginals
-    Lambda_Info = struct('bin_vals',zeros(lam_dim,nbins+1),'bin_ptr',zeros(num_samples,lam_dim));
-    for k=1:lam_dim
-        % Compute bins in each direction
-        Lambda_Info.bin_vals(k,:) = linspace(lam_domain(1,k),lam_domain(2,k),nbins+1);
-        % Compute bin pointer for each sample
-        [~,Lambda_Info.bin_ptr(:,k)] = histc(samples(:,k),Lambda_Info.bin_vals(k,:));
-    end
-    indices = nchoosek(1:lam_dim,2);
-    Lambda_pairs = struct('pair',struct());
-    for i=1:nchoosek(lam_dim,2) % Compute marginals for each 2d pair of coordinates
-        Lambda_pairs.pair(i).terms = indices(i,:);
-        Lambda_pairs.pair(i).marginal = zeros(nbins+1,nbins+1);
-        Lambda_pairs.pair(i).box_area = ( Lambda_Info.bin_vals(indices(i,1),2) - Lambda_Info.bin_vals(indices(i,1),1) ) * ( Lambda_Info.bin_vals(indices(i,2),2) - Lambda_Info.bin_vals(indices(i,2),1) );
-    end
-    
-    for i = 1:nchoosek(lam_dim,2) % Compute marginals
-        for j=1:num_samples
-            flag = 0;
-            k = 1;
-            while k<nbins+1 && flag == 0
-                if Lambda_Info.bin_ptr(j,indices(i,1)) == k
-                    kk = 1;
-                    while kk<nbins+1 && flag == 0
-                        if Lambda_Info.bin_ptr(j,indices(i,2)) == kk
-                            Lambda_pairs.pair(i).marginal(kk+1,k+1) = Lambda_pairs.pair(i).marginal(kk+1,k+1) + P_samples(j);
-                            flag = 1;
-                        end
-                        kk = kk+1;
-                    end
-                end
-                k = k+1;
-            end
-        end
-    end
-    
-    fig_num = 5;
-    for i = 1:nchoosek(lam_dim,2)
-        figure(fig_num);
-        zz = [Lambda_pairs.pair(i).marginal(2:end, 2:end) zeros(nbins, 1); zeros(1,nbins+1)];
-        surf(Lambda_Info.bin_vals(indices(i,1),:),Lambda_Info.bin_vals(indices(i,2),:),zz)
-        shading flat
-        axis([lam_domain(1,indices(i,1)) lam_domain(2,indices(i,1)) lam_domain(1,indices(i,2)) lam_domain(2,indices(i,2))])
-        s = ['Marginal in $\lambda_' int2str(indices(i,1)) '$ x $\lambda_' int2str(indices(i,2)) '$'];
-        %title(s, 'interpreter', 'latex', 'FontSize', 15)
-        s = ['$\lambda_' int2str(indices(i,1)) '$'];
-        xlabel(s, 'interpreter', 'latex', 'FontSize', 15)
-        s = ['$\lambda_' int2str(indices(i,2)) '$'];
-        ylabel(s, 'interpreter', 'latex', 'FontSize', 15)
-        colorbar
-        fig_num = fig_num+1;
-    end
-    
-end
diff --git a/bet/sampling/adaptiveSampling.py b/bet/sampling/adaptiveSampling.py
index fe2bfeb5..90e4a554 100644
--- a/bet/sampling/adaptiveSampling.py
+++ b/bet/sampling/adaptiveSampling.py
@@ -277,7 +277,7 @@ def generalized_chains(self, param_min, param_max, t_set, kern,
         # hypercube/space-filling curve to fully explore parameter space - not
         # necessarily random). Call these Samples_old.
         (samples_old, data_old) = super(sampler, self).random_samples(
-                initial_sample_type, param_min, param_max, psavefile,
+                initial_sample_type, param_min, param_max, savefile,
                 self.num_chains, criterion)
         self.num_samples = self.chain_length * self.num_chains
         comm.Barrier()
@@ -326,8 +326,10 @@ def generalized_chains(self, param_min, param_max, t_set, kern,
             mdat['step_ratios'] = all_step_ratios
             mdat['samples'] = samples
             mdat['data'] = data
-            super(sampler, self).save(mdat, psavefile)
-
+            if size > 1:
+                super(sampler, self).save(mdat, psavefile)
+            else:
+                super(sampler, self).save(mdat, savefile)
             MYsamples_old = samples_new
 
         # collect everything
diff --git a/setup.py b/setup.py
index 340beff2..41003a6e 100644
--- a/setup.py
+++ b/setup.py
@@ -5,15 +5,18 @@
 '''
 The python script for building the BET package and subpackages.
 '''
-from distutils.core import setup
+try:
+  from setuptools import setup
+except ImportError:
+    from distutils.core import setup
 
 setup(name='bet',
       version='0.2.0',
       description='Butler, Estep, Tavener method',
       author='Steven Mattis',
       author_email='steve.a.mattis@gmail.com',
-      liscense='GNU LGPL',
+      license='GNU LGPL',
       url='https://github.com/UT-CHG/BET',
       packages=['bet', 'bet.sampling', 'bet.calculateP', 'bet.postProcess'],
-      install_requires=['matplotlib', 'mpl_toolkits', 'pyDOE', 'scipy',
+      install_requires=['matplotlib', 'pyDOE', 'scipy',
           'numpy'])