diff --git a/.gitmodules b/.gitmodules
index 3d19f46050..8ade2a5c15 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -53,3 +53,6 @@
 [submodule "external/base/samplers/looplessFluxSampler"]
 	path = external/base/samplers/looplessFluxSampler
 	url = https://github.com/rmtfleming/looplessFluxSampler
+[submodule "external/base/samplers/constrained-logconcave-sampler"]
+	path = external/base/samplers/constrained-logconcave-sampler
+	url = https://github.com/Bounciness/constrained-logconcave-sampler
diff --git a/external/base/samplers/constrained-logconcave-sampler b/external/base/samplers/constrained-logconcave-sampler
new file mode 160000
index 0000000000..6ba5af35a7
--- /dev/null
+++ b/external/base/samplers/constrained-logconcave-sampler
@@ -0,0 +1 @@
+Subproject commit 6ba5af35a72ce59750c438556f4bd58e9a34e1b2
diff --git a/external/base/utilities/histogram_distance/build.m b/external/base/utilities/histogram_distance/build.m
new file mode 100755
index 0000000000..7771e24264
--- /dev/null
+++ b/external/base/utilities/histogram_distance/build.m
@@ -0,0 +1,11 @@
+% Build the C/C++ files provided in the package
+%
+% @author: B. Schauerte
+% @date:   2009
+% @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+ 
+cpp_files=dir('*.cpp');
+for i=1:length(cpp_files)
+  fprintf('Building %d of %d: %s\n',i,length(cpp_files),cpp_files(i).name);
+  mex(cpp_files(i).name);
+end
\ No newline at end of file
diff --git a/external/base/utilities/histogram_distance/chi_square_statistics.m b/external/base/utilities/histogram_distance/chi_square_statistics.m
new file mode 100755
index 0000000000..ff44d93326
--- /dev/null
+++ b/external/base/utilities/histogram_distance/chi_square_statistics.m
@@ -0,0 +1,55 @@
+function d=chi_square_statistics(XI,XJ)
+  % Implementation of the Chi^2 distance to use with pdist
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+  for i=1:m
+    for j=1:p
+      m=(XI(1,j) + XJ(i,j)) / 2;
+			if m ~= 0 % if m == 0, then xi and xj are both 0 ... this way we avoid the problem with (xj - m)^2 / m = (0 - 0)^2 / 0 = 0 / 0 = ?
+				d(i,1) = d(i,1) + ((XI(1,j) - m)^2 / m); % XJ is the model! makes it possible to determine each "likelihood" that XI was drawn from each of the models in XJ
+			end
+    end
+  end
\ No newline at end of file
diff --git a/external/base/utilities/histogram_distance/chi_square_statistics_fast.cpp b/external/base/utilities/histogram_distance/chi_square_statistics_fast.cpp
new file mode 100755
index 0000000000..f89b6521f6
--- /dev/null
+++ b/external/base/utilities/histogram_distance/chi_square_statistics_fast.cpp
@@ -0,0 +1,101 @@
+/**
+ * Copyright 2009 B. Schauerte. All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without 
+ * modification, are permitted provided that the following conditions are 
+ * met:
+ * 
+ *    1. Redistributions of source code must retain the above copyright 
+ *       notice, this list of conditions and the following disclaimer.
+ * 
+ *    2. Redistributions in binary form must reproduce the above copyright 
+ *       notice, this list of conditions and the following disclaimer in 
+ *       the documentation and/or other materials provided with the 
+ *       distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+ * DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+ * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+ * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *  
+ * The views and conclusions contained in the software and documentation
+ * are those of the authors and should not be interpreted as representing 
+ * official policies, either expressed or implied, of B. Schauerte.
+ */
+
+/**
+ * chi_square_statistics_fast
+ *
+ * Fast C/C++ calculation of the chi-square statistics (compatible with pdist).
+ * (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+ *       Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+ *
+ * @author: B. Schauerte
+ * @date:   2009
+ * @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+ */
+#include "mex.h"
+
+#define SQR(x) ((x)*(x))
+
+void
+mexFunction (int nlhs, mxArray* plhs[], 
+	int nrhs, const mxArray* prhs[])
+{
+  mwSize i = 0, j = 0; /* variables for for-loops */
+  
+  /* Check number of input parameters */
+	if (nrhs != 2) 
+  {
+  	mexErrMsgTxt("Two inputs required.");
+  } 
+  else 
+		if (nlhs > 1) 
+	  {
+  		mexErrMsgTxt("Wrong number of output arguments.");
+	  }    
+
+  /* Check type of input parameters */
+	if (!mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1])) 
+		mexErrMsgTxt("Input should be double.\n");
+	
+	/* Input data */
+	const mxArray* XI = prhs[0];
+	const mxArray* XJ = prhs[1];
+	const double* XI_data = mxGetPr(XI);
+	const double* XJ_data = mxGetPr(XJ);
+	/* some helper variables */
+	const mwSize m = mxGetM(XJ); /* number of samples of p */
+  const mwSize p = mxGetN(XI); /* dimension of samples */
+	if (p != mxGetN(XJ)) 
+		mexErrMsgTxt("Dimension mismatch (1).\n");
+	if (1 != mxGetM(XI)) 
+		mexErrMsgTxt("Dimension mismatch. XI has to be an (1,n) vector.\n");
+	/* Output data */
+	mxArray* OUT = mxCreateNumericMatrix (m, 1, mxDOUBLE_CLASS, mxREAL);
+	plhs[0] = OUT;
+	double* out_data = mxGetPr(OUT);
+	
+	for (i = 0; i < m; i++) /* initialize output array */
+		out_data[i] = 0;
+
+	for (j = 0; j < p; j++)
+	{
+		const double xi = XI_data[j];
+		for (i = 0; i < m; i++)
+		{
+			const double mean = (xi + *XJ_data++) / 2.0;
+			if (mean != 0)
+				out_data[i] += SQR(xi - mean) / mean;
+		}
+	}
+	
+	return;
+}
diff --git a/external/base/utilities/histogram_distance/hist_dist_example.m b/external/base/utilities/histogram_distance/hist_dist_example.m
new file mode 100755
index 0000000000..2fcf41d8c8
--- /dev/null
+++ b/external/base/utilities/histogram_distance/hist_dist_example.m
@@ -0,0 +1,143 @@
+% An example of how to use the histogram distance functions for image
+% matching.
+%
+% Please note that this is a demo to show case the usage of the histogram
+% functions. But, in general, matching images solely based on their color
+% histograms ist - imho - not the best idea, unless you have a really large
+% image database.
+%
+% Some of the histogram distance functions have been used for outlier
+% reduction when learning color term/name models from web images, see:
+%
+% [1] B. Schauerte, G. A. Fink, "Web-based Learning of Naturalized Color 
+%     Models for Human-Machine Interaction". In Proceedings of the 12th 
+%     International Conference on Digital Image Computing: Techniques and 
+%     Applications (DICTA), IEEE, Sydney, Australia, December 1-3, 2010. 
+% [2] B. Schauerte, R. Stiefelhagen, "Learning Robust Color Name Models 
+%     from Web Images". In Proceedings of the 21st International Conference
+%     on Pattern Recognition (ICPR), Tsukuba, Japan, November 11-15, 2012
+%
+% If you use and like this code, you are kindly requested to cite some of 
+% the work above. 
+%
+% Anyway, I hope it saves you some work. Have fun with it ;)
+%
+% @author: B. Schauerte
+% @date:   2012,2013
+% @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+
+%%
+% Build the .cpp files, if necessary
+if ~exist('chi_square_statistics_fast','file') && exist('./build.m')
+  build;
+end
+
+% Download some random sample images from the Google-512 dataset. For 
+% information about the dataset see: 
+%
+% [1] B. Schauerte, G. A. Fink, "Web-based Learning of Naturalized Color 
+%     Models for Human-Machine Interaction". In Proceedings of the 12th 
+%     International Conference on Digital Image Computing: Techniques and 
+%     Applications (DICTA), IEEE, Sydney, Australia, December 1-3, 2010. 
+% [2] B. Schauerte, R. Stiefelhagen, "Learning Robust Color Name Models 
+%     from Web Images". In Proceedings of the 21st International Conference
+%     on Pattern Recognition (ICPR), Tsukuba, Japan, November 11-15, 2012
+
+%colornames={'red','green','blue','yellow', ...
+%  'pink','purple','brown','orange', ...
+%  'black','grey','white'};
+colornames={'red','green','blue','yellow'};
+fendings={'jpeg','png','gif'};
+tmp_foldername='google-512-samples';
+n_samples = 100;
+
+% download the images in a temporary folder
+if ~exist(tmp_foldername,'dir'), mkdir(tmp_foldername); end % create temporary directory
+filenames=cell(n_samples,1);
+for i=1:n_samples
+  colorname=colornames{randi(numel(colornames))};
+  %colorname=colornames{1};
+  for j=1:numel(fendings)
+    url=sprintf('https://cvhci.anthropomatik.kit.edu/~bschauer/datasets/google-512/images-resized-128/%s+color/%d.%s',colorname,i,fendings{j});
+    filename=sprintf('%s_%d.%s',colorname,i,fendings{j});
+    [~,status] = urlwrite(url,fullfile(tmp_foldername,filename));
+    if status
+      filenames{i} = filename;
+      break; 
+    end
+  end
+end
+
+%%
+% We simply use all files that have already been downloaded
+filenames=dir(fullfile(tmp_foldername,'*_*.*'));
+filenames={filenames.name};
+n_samples=numel(filenames);
+
+%%
+% calculate color image histograms
+n_bins=4;
+edges=(0:(n_bins-1))/n_bins;
+histograms=zeros(n_samples,n_bins*n_bins*n_bins);
+for i=1:n_samples
+  I=imread(fullfile(tmp_foldername,filenames{i}));
+  IR=imresize(I,[64 64]);
+  IR=im2double(IR);
+  
+  [~,r_bins] = histc(reshape(IR(:,:,1),1,[]),edges); r_bins = r_bins + 1;
+  [~,g_bins] = histc(reshape(IR(:,:,1),1,[]),edges); g_bins = g_bins + 1;
+  [~,b_bins] = histc(reshape(IR(:,:,1),1,[]),edges); b_bins = b_bins + 1;
+  
+  histogram=zeros(n_bins,n_bins,n_bins);
+  for j=1:numel(r_bins)
+    histogram(r_bins(j),g_bins(j),b_bins(j)) = histogram(r_bins(j),g_bins(j),b_bins(j)) + 1;
+  end
+  histograms(i,:) = reshape(histogram,1,[]) / sum(histogram(:)); % normalize, better for all probabilistic methods
+end
+
+%%
+% match histograms and show best matching pairs
+dist_func=@chi_square_statistics_fast;
+% 1. You can use pdist to calculate the distances, iff the distance measure
+%    is symmetric
+%D=squareform(pdist(histograms,dist_func)); % use pdist to calculate the distance for all image pairs 
+% 2. Use the following loop to calculate the distances, iff the measure is
+%    not symmetric
+% D=zeros(size(histograms,1),size(histograms,1));
+% for i=1:size(histograms,1)
+%   for j=1:size(histograms,1)
+%     D(i,j) = dist_func(histograms(i,:),histograms(j,:));
+%   end
+% end
+% 2. ... alternatively, use pdist2
+D=pdist2(histograms,histograms,dist_func);
+
+D(D == 0) = NaN;
+n_show_samples=5; % number of samples for the illustration
+figure('name','Random images (left) with their best (middle) and worst (right) match');
+c = 1;
+rand_indices=randperm(numel(filenames));
+for i=1:n_show_samples
+  % image we want to match
+  I=imread(fullfile(tmp_foldername,filenames{rand_indices(i)}));
+  if numel(size(I)) > 3, I=I(:,:,1:3); end
+  subplot(n_show_samples,3,c); imshow(I); c = c + 1;
+  
+  % best match
+  %[d,j]=min(D(rand_indices(i),:)); % if distances are not symmetric, then
+  % it might be useful to try the other order, see below, depending on the
+  % definition of the metric
+  [d,j]=min(D(:,rand_indices(i)));
+  I=imread(fullfile(tmp_foldername,filenames{j}));
+  if numel(size(I)) > 3, I=I(:,:,1:3); end
+  subplot(n_show_samples,3,c); imshow(I); title(sprintf('Dist: %.3f',d*100)); c = c + 1;
+  
+  % worst match
+  %[d,j]=max(D(rand_indices(i),:)); % if distances are not symmetric, then
+  % it might be useful to try the other order, see below, depending on the
+  % definition of the metric
+  [d,j]=max(D(:,rand_indices(i)));
+  I=imread(fullfile(tmp_foldername,filenames{j}));
+  if numel(size(I)) > 3, I=I(:,:,1:3); end
+  subplot(n_show_samples,3,c); imshow(I); title(sprintf('Dist: %.3f',d*100)); c = c + 1;
+end
\ No newline at end of file
diff --git a/external/base/utilities/histogram_distance/histogram_intersection.m b/external/base/utilities/histogram_distance/histogram_intersection.m
new file mode 100755
index 0000000000..37bb6b7491
--- /dev/null
+++ b/external/base/utilities/histogram_distance/histogram_intersection.m
@@ -0,0 +1,51 @@
+function d=histogram_intersection(XI,XJ)
+  % Implementation of the histogram intersection distance to use with pdist
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+  sxi=sum(XI);
+  for i=1:m
+    d(i,1) = 1 - (sum(min(XI, XJ(i,:))) / sxi);
+  end
\ No newline at end of file
diff --git a/external/base/utilities/histogram_distance/jeffrey_divergence.m b/external/base/utilities/histogram_distance/jeffrey_divergence.m
new file mode 100755
index 0000000000..11652ca2c6
--- /dev/null
+++ b/external/base/utilities/histogram_distance/jeffrey_divergence.m
@@ -0,0 +1,55 @@
+function d=jeffrey_divergence(XI,XJ)
+  % Implementation of the Jeffrey Divergence
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+  for i=1:m
+    for j=1:p
+      m=(XJ(i,j) + XI(1,j)) / 2;
+			if m ~= 0  % if m == 0, then xi == xj == 0
+				d(i,1) = d(i,1) + (XI(1,j) * log(XI(1,j) / m)) + (XJ(i,j) * log(XJ(i,j) / m));
+			end
+    end
+  end
diff --git a/external/base/utilities/histogram_distance/jensen_shannon_divergence.m b/external/base/utilities/histogram_distance/jensen_shannon_divergence.m
new file mode 100755
index 0000000000..4b604a4ed7
--- /dev/null
+++ b/external/base/utilities/histogram_distance/jensen_shannon_divergence.m
@@ -0,0 +1,41 @@
+function d=jensen_shannon_divergence(XI,XJ)
+  % Implementation of the Jensen-Shannon Divergence to use with pdist
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+  
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  d=jeffrey_divergence(XI,XJ);
+  d=d / 2;
diff --git a/external/base/utilities/histogram_distance/kolmogorov_smirnov_distance.m b/external/base/utilities/histogram_distance/kolmogorov_smirnov_distance.m
new file mode 100755
index 0000000000..6e4b0a74e8
--- /dev/null
+++ b/external/base/utilities/histogram_distance/kolmogorov_smirnov_distance.m
@@ -0,0 +1,55 @@
+function d=kolmogorov_smirnov_distance(XI,XJ)
+  % Implementation of the Kolmogorov-Smirnov distance to use with pdist
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+  
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+  cxi=cumsum(XI,2); % cumulative histograms
+  cxj=cumsum(XJ,2); % cumulative histograms
+  
+  for i=1:m
+    for j=1:p
+      d(i,1) = max(d(i,1), abs(cxi(1,j) - cxj(i,j)));
+    end
+  end
diff --git a/external/base/utilities/histogram_distance/kullback_leibler_divergence.m b/external/base/utilities/histogram_distance/kullback_leibler_divergence.m
new file mode 100755
index 0000000000..80e8f30add
--- /dev/null
+++ b/external/base/utilities/histogram_distance/kullback_leibler_divergence.m
@@ -0,0 +1,55 @@
+function d=kullback_leibler_divergence(XI,XJ)
+  % Implementation of the Kullback-Leibler Divergence to use with pdist
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+  
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+  for i=1:m
+    for j=1:p
+      %d(i,1) = d(i,1) + (XJ(i,j) * log(XJ(i,j) / XI(1,j))); % XI is the model!
+			if XI(1,j) ~= 0
+				d(i,1) = d(i,1) + (XI(1,j) * log(XI(1,j) / XJ(i,j))); % XJ is the model! makes it possible to determine each "likelihood" that XI was drawn from each of the models in XJ
+			end
+    end
+  end
diff --git a/external/base/utilities/histogram_distance/license.txt b/external/base/utilities/histogram_distance/license.txt
new file mode 100644
index 0000000000..5842fc9539
--- /dev/null
+++ b/external/base/utilities/histogram_distance/license.txt
@@ -0,0 +1,24 @@
+Copyright (c) 2012, Boris Schauerte
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in
+      the documentation and/or other materials provided with the distribution
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
diff --git a/external/base/utilities/histogram_distance/match_distance.m b/external/base/utilities/histogram_distance/match_distance.m
new file mode 100755
index 0000000000..b30e627996
--- /dev/null
+++ b/external/base/utilities/histogram_distance/match_distance.m
@@ -0,0 +1,55 @@
+function d=match_distance(XI,XJ)
+  % Implementation of the match distance to use with pdist
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+  
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+  cxi=cumsum(XI,2); % cumulative histograms
+  cxj=cumsum(XJ,2); % cumulative histograms
+  
+  for i=1:m
+    for j=1:p
+      d(i,1) = d(i,1) + abs(cxi(1,j) - cxj(i,j));
+    end
+  end
diff --git a/external/base/utilities/histogram_distance/quadratic_form_distance.m b/external/base/utilities/histogram_distance/quadratic_form_distance.m
new file mode 100755
index 0000000000..6c8189afee
--- /dev/null
+++ b/external/base/utilities/histogram_distance/quadratic_form_distance.m
@@ -0,0 +1,56 @@
+function d=quadratic_form_distance(XI,XJ,A)
+  % Implementation of the Quadratic form distance
+  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
+  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
+  %
+  % @author: B. Schauerte
+  % @date:   2009
+  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/
+
+  % Copyright 2009 B. Schauerte. All rights reserved.
+  % 
+  % Redistribution and use in source and binary forms, with or without 
+  % modification, are permitted provided that the following conditions are 
+  % met:
+  % 
+  %    1. Redistributions of source code must retain the above copyright 
+  %       notice, this list of conditions and the following disclaimer.
+  % 
+  %    2. Redistributions in binary form must reproduce the above copyright 
+  %       notice, this list of conditions and the following disclaimer in 
+  %       the documentation and/or other materials provided with the 
+  %       distribution.
+  % 
+  % THIS SOFTWARE IS PROVIDED BY B. SCHAUERTE ''AS IS'' AND ANY EXPRESS OR 
+  % IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+  % WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+  % DISCLAIMED. IN NO EVENT SHALL B. SCHAUERTE OR CONTRIBUTORS BE LIABLE 
+  % FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+  % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+  % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
+  % BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+  % WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+  % OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+  % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+  % 
+  % The views and conclusions contained in the software and documentation
+  % are those of the authors and should not be interpreted as representing 
+  % official policies, either expressed or implied, of B. Schauerte.
+  
+	if nargin < 3
+		error('Not enough arguments, i.e. similarity matrix is missing')
+	end
+
+  m=size(XJ,1); % number of samples of p
+  p=size(XI,2); % dimension of samples
+  
+  assert(p == size(XJ,2)); % equal dimensions
+  assert(size(XI,1) == 1); % pdist requires XI to be a single sample
+  
+  d=zeros(m,1); % initialize output array
+  
+	diff=zeros(1,p);
+  for i=1:m
+		diff=(XI - XJ(i,:));
+		d(i) = sqrt( diff * A * diff' );
+  end
\ No newline at end of file
diff --git a/initCobraToolbox.m b/initCobraToolbox.m
index dafd08530c..5dd1e26e8d 100644
--- a/initCobraToolbox.m
+++ b/initCobraToolbox.m
@@ -388,6 +388,9 @@ function initCobraToolbox(updateToolbox)
         [solverOK,solverInstalled] = changeCobraSolver(supportedSolversNames{i},...
                                      SOLVERS.(supportedSolversNames{i}).type{1},...
                                      0, 2);
+        if strcmp(SOLVERS.(supportedSolversNames{i}),'gurobi')
+            disp(SOLVERS.(supportedSolversNames{i}));
+        end
         if solverOK
             SOLVERS.(supportedSolversNames{i}).working = true;
         end
diff --git a/src/analysis/topology/getCorrespondingRows.m b/src/analysis/topology/getCorrespondingRows.m
index a1c8adfc17..2d3900d535 100644
--- a/src/analysis/topology/getCorrespondingRows.m
+++ b/src/analysis/topology/getCorrespondingRows.m
@@ -11,7 +11,7 @@
 %    S:                     `m x n` stoichiometric matrix
 %    rowBool:               `m x 1` boolean vector
 %    colBool:               `n x 1` boolean vector
-%    mode:                  'exclusive' or 'inclusive'
+%    mode:                  'exclusive' , 'inclusive' or 'partial'
 %
 % OUTPUT:
 %    restrictedRowBool:     `m x 1` boolean vector
diff --git a/src/base/io/utilities/convertOldCouplingFormat.m b/src/base/io/utilities/convertOldCouplingFormat.m
new file mode 100644
index 0000000000..2f7978840b
--- /dev/null
+++ b/src/base/io/utilities/convertOldCouplingFormat.m
@@ -0,0 +1,135 @@
+function model = convertOldCouplingFormat(model)
+%Converts an old style model implementation of coupling constraints to a
+%new style
+%
+% INPUT
+% model     model with model.A but without model.d   
+%
+% OUTPUT
+%    model:     A COBRA model structure with the following fields
+%
+%                * `.S` - The stoichiometric matrix
+%                * `.c` - Objective coeff vector
+%                * `.lb` - Lower bound vector
+%                * `.ub` - Upper bound vector              
+%                  * `.b`: accumulation/depletion vector (default 0 for each metabolite).
+%                  * `.C`: the Constraint matrix;
+%                  * `.d`: the right hand side vector for C;
+%                  * `.dsense`: the constraint sense vector;
+
+if isfield(model,'A') && isfield(model,'S')
+%     if printLevel >=1
+%         warning('The inserted Model contains an old style coupling matrix (A). The MAtrix will be converted into a Coupling Matrix (C) and fields will be adapted.')
+%     end
+    [nMets,nRxns] = size(model.S);
+    if size(model.S,1) == size(model.A,1)
+        bool=strncmp('slack_',model.mets,length('slack_'));
+        if any(bool)
+            %this is a draft harvey or harvetta
+            nMetNames = length(model.metNames);
+            if nMets~=nMetNames
+                model.metNames=cell(nnz(~bool),1);
+            end
+            
+            nMetFormulas = length(model.metFormulas);
+            if nMets~=nMetFormulas
+                model.metFormulas=cell(nnz(~bool),1);
+            end
+            
+            nMetCharges = length(model.metCharge);
+            if nMets~=nMetCharges
+                model.metCharges=ones(nnz(~bool),1)*NaN;
+            end
+            model=rmfield(model,'metCharge');
+            %         %metCharges is an empty character array but needs to be a numeric vector
+            %         model.metCharges=ones(size(model.mets,1),1)*NaN;
+            
+            mets=model.mets(~bool);
+            ctrs=model.mets(bool);
+            
+            S=model.A(~bool,:);
+            csense=columnVector(model.csense(~bool));
+            b=columnVector(model.b(~bool));
+            C=model.A(bool,:);
+            dsense=columnVector(model.csense(bool));
+            d=columnVector(model.b(bool));
+                
+            %replace some fields
+            model.rxnNotes=cell(nRxns,1);
+            model.rxnECNumbers=cell(nRxns,1);
+            model.rxnReferences=cell(nRxns,1);
+            model=rmfield(model,'rules');
+            model.rxnNames=cell(nRxns,1);
+            
+            model.rxnECNumbers = columnVector(model.rxnECNumbers);
+            model.grRulesNotes = columnVector(model.grRulesNotes);
+            model.rxnECNumbers = columnVector(model.rxnECNumbers);
+            model.rxnReferences = columnVector(model.rxnReferences);
+            
+            %grRule{x} needs to be a character array, but some grRules{x} are cells
+            bool=cellfun(@(y) ischar(y) , model.grRules);
+            for i=1:length(bool)
+                if ~bool(i)
+                    %replace cell with contents of cell
+                    tmp=model.grRules{i};
+                    model.grRules{i}=tmp{1};
+                end
+            end
+            bool=cellfun(@(y) ischar(y) , model.grRules);
+            if ~any(bool)
+                error('convertOldCouplingFormat: model.grRules format still wrong')
+            end
+            
+            %add any missing model.genes that were present in model.grRules
+            %uses code from generateRules.m
+            [~,genes] = preparseGPR(model.grRules);  % preparse all model.grRules
+            allGenes =  unique([genes{~cellfun(@isempty,genes)}]); %Get the unique gene list
+            if (~isfield(model, 'genes'))
+                newGenes = allGenes;
+            else
+                %         C = setdiff(A,B) for vectors A and B, returns the values in A that
+                %         are not in B with no repetitions. C will be sorted.
+                newGenes = setdiff(allGenes,model.genes);
+            end
+            if ~isempty(newGenes)
+                model = addGenes(model,newGenes);
+            end
+                       
+        else
+            error('not clear what type of old model this is')
+        end
+    else
+        % get the Constraint data
+        C = model.A(nMets+1:end,:);
+        ctrs = columnVector(model.mets(nMets+1:end));
+        if isempty(ctrs)
+            ctrs=cell(size(model.A,1)-size(model.S,1),1);
+        end
+        dsense =  columnVector(model.csense(nMets+1:end));
+        d = columnVector(model.b(nMets+1:end));
+        
+        % now, we assume, that those are the only modified fields, if not,
+        % something is seriously broken.
+        mets = columnVector(model.mets(1:nMets));
+        b = columnVector(model.b(1:nMets));
+        csense = columnVector(model.csense(1:nMets));
+    end
+
+
+    % set the constraint data
+    model.S = S;
+    model.mets = mets;
+    model.csense = csense;
+    model.b = b;
+    
+    model.C = C;
+    model.ctrs = ctrs;
+    model.dsense = dsense;
+    model.d = d;
+    
+    %remove model.A
+    model = rmfield(model,'A');
+end
+
+model = orderfields(model);
+
diff --git a/src/base/io/utilities/convertOldStyleModel.m b/src/base/io/utilities/convertOldStyleModel.m
index dbe44103ce..3661daeb37 100644
--- a/src/base/io/utilities/convertOldStyleModel.m
+++ b/src/base/io/utilities/convertOldStyleModel.m
@@ -104,31 +104,8 @@
 % get the defined field properties.    
 definedFields = getDefinedFieldProperties();
     
-    
-if isfield(model,'A') && isfield(model,'S')
-    if printLevel >=1
-        warning('The inserted Model contains an old style coupling matrix (A). The MAtrix will be converted into a Coupling Matrix (C) and fields will be adapted.')
-    end
-    nMets = size(model.S,1);
-    % get the Constraint data
-    C = model.A(nMets+1:end,:);     
-    ctrs = columnVector(model.mets(nMets+1:end));    
-    dsense =  columnVector(model.csense(nMets+1:end));
-    d = columnVector(model.b(nMets+1:end));
-    % set the constraint data
-    model.C = C;
-    model.ctrs = ctrs;
-    model.dsense = dsense;
-    model.d = d;
-    % now, we assume, that those are the only modified fields, if not,
-    % something is seriously broken.    
-    model.mets = columnVector(model.mets(1:nMets));
-    model.b = columnVector(model.b(1:nMets));
-    model.csense = columnVector(model.csense(1:nMets));
-    model = rmfield(model,'A');
-    
-end
-    
+%convert from old coupling constraints if necessary
+model = convertOldCouplingFormat(model);
 
 % convert old fields to current fields. 
 for i = 1:numel(oldFields)
@@ -330,4 +307,4 @@
     model.genes = cell(0,1);
 end
 
-    
\ No newline at end of file
+model = orderfields(model);
\ No newline at end of file
diff --git a/src/base/preparseGPR.m b/src/base/preparseGPR.m
index 80df8cca2d..e05ce300a9 100644
--- a/src/base/preparseGPR.m
+++ b/src/base/preparseGPR.m
@@ -14,6 +14,18 @@
 %
 % .. Author: -  Laurent Heirendt - December 2017
 
+%using regexprep all cells must be char row vectors.
+glt=length(grRules);
+for g=1:glt
+    if ~ischar(grRules{g})
+        if iscell(grRules{g})
+            tmp=grRules{g};
+            grRules{g}=tmp{1};
+        else
+            error(['grRules ' int2str(g)  ' is not a character array'])
+        end
+    end
+end
     preParsedGrRules = regexprep(grRules, '[\]\}]',')'); %replace other brackets by parenthesis.
     preParsedGrRules = regexprep(preParsedGrRules, '[\[\{]','('); %replace other brackets by parenthesis.
     preParsedGrRules = regexprep(preParsedGrRules,'([\(])\s*','$1'); %replace all spaces after opening parenthesis
diff --git a/src/base/solvers/buildLPproblemFromModel.m b/src/base/solvers/buildLPproblemFromModel.m
index 0741ff2e21..31c8c080c7 100644
--- a/src/base/solvers/buildLPproblemFromModel.m
+++ b/src/base/solvers/buildLPproblemFromModel.m
@@ -1,6 +1,13 @@
 function LPproblem = buildLPproblemFromModel(model, checked)
 % Builds an COBRA Toolbox LP problem structure from a COBRA Toolbox model structure.
 %
+% 
+%.. math::
+%
+%    max/min  ~& c^T x \\
+%    s.t.     ~& [S, E; C, D] x <=> b ~~~~~~~~~~~:y \\
+%             ~& lb \leq x \leq ub~~~~:w
+%
 % USAGE:
 %
 %    LPproblem = buildLPproblemFromModel(model)
@@ -45,6 +52,8 @@
     checked = true;
 end
 
+%backward compatibility with old formulation of coupling constraints
+
 %Build some fields, if they don't exist
 
 optionalFields = {'C','d','dsense','E','evarlb','evarub','evarc','D'};
diff --git a/src/base/utilities/dispMatrix.m b/src/base/utilities/dispMatrix.m
new file mode 100644
index 0000000000..e4dc641d01
--- /dev/null
+++ b/src/base/utilities/dispMatrix.m
@@ -0,0 +1,34 @@
+function dispMatrix(A,mode)
+% display a matrix A in a tight format
+%
+% INPUT
+% A     m x n matrix
+% mode  'Z' integers (default)
+%       'N' natural number
+%
+
+% from: https://stackoverflow.com/questions/7919004/tightening-the-display-of-matrices-in-matlab
+
+if ~exist('mode','var')
+    mode='Z';
+end
+
+switch mode
+    case 'disp'
+        disp(A)
+    case 'N'
+        fprintf([repmat('%d ',1,size(A,2)) '\n'],A');
+    case 'Z'
+        fprintf([repmat(sprintf('%% %dd',max(floor(log10(abs(A(:)))))+2+any(A(:)<0)),1,size(A,2)) '\n'],A');
+    case 'nonzeroZ'
+        A = removeZeroRowsCols(A);
+        fprintf([repmat(sprintf('%% %dd',max(floor(log10(abs(A(:)))))+2+any(A(:)<0)),1,size(A,2)) '\n'],A');
+    case 'nonzeroN'
+        A = removeZeroRowsCols(A);
+        fprintf([repmat('%d ',1,size(A,2)) '\n'],A');
+        
+        
+end
+
+end
+
diff --git a/src/base/utilities/removeZeroRowsCols.m b/src/base/utilities/removeZeroRowsCols.m
new file mode 100644
index 0000000000..41b9d5ec31
--- /dev/null
+++ b/src/base/utilities/removeZeroRowsCols.m
@@ -0,0 +1,10 @@
+function A = removeZeroRowsCols(A)
+%removes all zero rows and columns from a matrix
+
+% Remove zero rows
+A( all(~A,2), : ) = [];
+% Remove zero columns
+A( :, all(~A,1) ) = [];
+
+end
+
diff --git a/src/reconstruction/modelGeneration/modelVerification/checkDatabaseIDs.m b/src/reconstruction/modelGeneration/modelVerification/checkDatabaseIDs.m
new file mode 100644
index 0000000000..a0a926c09d
--- /dev/null
+++ b/src/reconstruction/modelGeneration/modelVerification/checkDatabaseIDs.m
@@ -0,0 +1,61 @@
+function results = checkDatabaseIDs(model,results)
+% Checks the model for validity of database identifiers
+%
+% USAGE:
+%
+%    results = checkDatabaseIDs(model,results)
+%
+% INPUT:
+%    model:       a structure that represents the COBRA model.
+%    results:     the results structure for this test
+%
+% OUTPUT:
+%
+%    results:     a struct with problematic database ids added.
+%
+% .. Authors:
+%       - Thomas Pfau, May 2017
+
+dbMappings = getDefinedFieldProperties('Database',true);
+
+for i= 1:size(dbMappings,1)
+    if isfield(model,dbMappings{i,3})
+        fits = cellfun(@(x) isempty(x) || checkID(x,dbMappings{i,5}),model.(dbMappings{i,3}));
+        %Only add the something if we really have wrong IDs.
+        if any(~fits)
+            if ~isfield(results,'checkDatabaseIDs')
+                results.checkDatabaseIDs = struct();
+            end
+            if ~isfield(results.checkDatabaseIDs,'invalidIDs')
+                results.checkDatabaseIDs.invalidIDs = struct();
+            end
+            results.checkDatabaseIDs.invalidIDs.(dbMappings{i,3}) = cell(size(model.(dbMappings{i,3})));
+            results.checkDatabaseIDs.invalidIDs.(dbMappings{i,3})(:) = {'valid'};
+            results.checkDatabaseIDs.invalidIDs.(dbMappings{i,3})(~fits) = model.(dbMappings{i,3})(~fits);
+        end
+    end
+end
+end
+
+
+function accepted = checkID(id,pattern)
+% Checks the the given id(s), i.e. strings split by ; versus the pattern
+%
+% USAGE:
+%
+%    accepted = checkID(id,pattern)
+%
+% INPUT:
+%    id:          A String representing ids (potentially separated by ;)
+%    pattern:     The pattern to check the id(s) against.
+%
+% OUTPUT:
+%
+%    accepted:     Whether all ids are ok.
+%
+% .. Authors:
+%       - Thomas Pfau, May 2017
+ids = strsplit(id,';');
+matches = regexp(ids,pattern);
+accepted = all(~cellfun(@isempty,matches));
+end
diff --git a/src/reconstruction/modelGeneration/modelVerification/checkPresentFields.m b/src/reconstruction/modelGeneration/modelVerification/checkPresentFields.m
new file mode 100644
index 0000000000..934c225011
--- /dev/null
+++ b/src/reconstruction/modelGeneration/modelVerification/checkPresentFields.m
@@ -0,0 +1,99 @@
+function results = checkPresentFields(fieldProperties,model, results)
+% Check the model fields for consistency with the given fieldProperties and
+% update the results struct.
+% The desired properties of each field are described here:
+% https://github.com/opencobra/cobratoolbox/blob/master/docs/source/notes/COBRAModelFields.md
+%
+% USAGE:
+%
+%    results = checkPresentFields(fieldProperties,model, results)
+%
+% INPUT:
+%    fieldProperties:  field properties as obtained by
+%                      getDefinedFieldProperties
+%    model:            a structure that represents the COBRA model.
+%    results:          the results structure for this test
+%
+% OUTPUT:
+%
+%    results:          The updated results struct.
+%
+% .. Authors:
+%       - Thomas Pfau, May 2017
+
+presentFields = find(ismember(fieldProperties(:,1),fieldnames(model)));
+
+%Check all Field Sizes
+for i = 1:numel(presentFields)
+    testedField = fieldProperties{presentFields(i),1};
+    [x_size,y_size] = size(model.(testedField));
+    xFieldMatch = fieldProperties{presentFields(i),2};
+    yFieldMatch = fieldProperties{presentFields(i),3};
+    
+    checkX = ~isnan(xFieldMatch);
+    checkY = ~isnan(yFieldMatch);
+    if checkX
+        if ischar(xFieldMatch)
+            if ~isfield(model,xFieldMatch)
+                x_pres = 0;
+                if ~isfield(results.Errors,'missingFields')
+                    results.Errors.missingFields = {};
+                end
+                results.Errors.missingFields(end+1) = {xFieldMatch};
+            else
+                x_pres = numel(model.(xFieldMatch));
+            end
+            errorMessage = sprintf('%s: Size of %s does not match elements in %s', xFieldMatch,testedField,xFieldMatch);
+        elseif isnumeric(xFieldMatch)
+            errorMessage = sprintf('X Size of %s was %i. Expected %i',testedField, x_size,x_pres);
+            x_pres = xFieldMatch;
+        end
+        if x_pres ~= x_size
+            if ~isfield(results.Errors,'inconsistentFields')
+                results.Errors.inconsistentFields = struct();
+            end
+            results.Errors.inconsistentFields.(testedField) = errorMessage;
+        end
+    end
+    if checkY
+        if ischar(yFieldMatch)
+            if ~isfield(model,yFieldMatch)
+                y_pres = 0;
+                if ~isfield(results.Errors,'missingFields')
+                    results.Errors.missingFields = {};
+                end
+                results.Errors.missingFields(end+1) = {yFieldMatch};
+            else
+                y_pres = numel(model.(yFieldMatch));
+            end
+            errorMessage = sprintf('%s: Size of %s does not match elements in %s', yFieldMatch,testedField,yFieldMatch);
+        elseif isnumeric(yFieldMatch)
+            y_pres = yFieldMatch;
+            errorMessage = sprintf('Y Size of %s was %i. Expected %i',testedField, y_size,y_pres);
+        end
+        if y_pres ~= y_size
+            if ~isfield(results.Errors,'inconsistentFields')
+                results.Errors.inconsistentFields = struct();
+            end
+            results.Errors.inconsistentFields.(testedField) = errorMessage;
+        end
+    end
+    %Test the field content properties
+    %x is necessary here, since it is used for the eval below!
+    x = model.(testedField);
+    try
+        propertiesMatched = eval(fieldProperties{presentFields(i),4});
+    catch
+        propertiesMatched = false;
+    end
+    if ~propertiesMatched
+        if ~isfield(results.Errors,'propertiesNotMatched')
+            results.Errors.propertiesNotMatched = struct();
+        end
+        %results.Errors.propertiesNotMatched.(testedField) = 'Field does not match the required properties.;
+        results.Errors.propertiesNotMatched.(testedField) = ['Field does not match the required properties: model.' testedField ' is ' class(x) ' but must satisfy: ' fieldProperties{presentFields(i),4}];
+    end
+    
+end
+%disp('done')
+end
diff --git a/src/reconstruction/modelGeneration/modelVerification/verifyModel.m b/src/reconstruction/modelGeneration/modelVerification/verifyModel.m
index 20c15816ea..af811a0ef9 100644
--- a/src/reconstruction/modelGeneration/modelVerification/verifyModel.m
+++ b/src/reconstruction/modelGeneration/modelVerification/verifyModel.m
@@ -25,7 +25,7 @@
 %                   * 'checkDatabaseIDs', check whether the database identifiers in specified fields (please have a look at the documentation), match to the expected patterns for those databases.
 %                   * 'silentCheck', do not print any information. Only applies to the model structure check. (default is to print info)
 %                   * 'restrictToFields' restricts the check to the listed fields. This will lead to requiredFields being reduced to those fields present in the restricted fields. If an empty cell array is provided no restriction is applied. (default: {})
-%                   * 'FBAOnly' checks only fields relevant for FBA (default: false) 
+%                   * 'FBAOnly' checks only fields relevant for FBA (default: false)
 %
 % OUTPUT:
 %
@@ -33,7 +33,7 @@
 %                 additional field `Errors` indicating the problems with the
 %                 model structure detected by the `verifyModel` function.
 %                 Results of additional options are returned in fields with
-%                 the respective names. 
+%                 the respective names.
 %
 % EXAMPLES:
 %    1) Do a simple Field check with the default required fields.
@@ -119,21 +119,21 @@
 if massBalance || chargeBalance
     doChargeBalance = false;
     if chargeBalance
-         if checkFields(results,'metCharges',model)
-             doChargeBalance = true;
-         else
-             results.chargeBalance = struct();
-             results.chargeBalance.missingFields = 'The metCharges field is missing. Cannot determine the charge Balance';
-         end
+        if checkFields(results,'metCharges',model)
+            doChargeBalance = true;
+        else
+            results.chargeBalance = struct();
+            results.chargeBalance.missingFields = 'The metCharges field is missing. Cannot determine the charge Balance';
+        end
     end
     doMassBalance = false;
     if massBalance
-         if checkFields(results,'metFormulas',model)
-             doMassBalance = true;
-         else
-             results.massBalance = struct();
-             results.massBalance.missingFields = 'The metFormulas field is missing. Cannot determine the mass Balance';
-         end
+        if checkFields(results,'metFormulas',model)
+            doMassBalance = true;
+        else
+            results.massBalance = struct();
+            results.massBalance.missingFields = 'The metFormulas field is missing. Cannot determine the mass Balance';
+        end
     end
     
     if doMassBalance || doChargeBalance
@@ -149,17 +149,17 @@
         [massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, Elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(model,0);
         %put the fields, if simpleCheck is active and there are imbalanced,
         %or if simplecheck is not active.
-        if doMassBalance && (~simpleCheck || any(imBalancedRxnBool))            
+        if doMassBalance && (~simpleCheck || any(imBalancedRxnBool))
             results.massBalance = struct();
             results.massBalance.massImbalance = massImbalance;
             results.massBalance.imBalancedMass = imBalancedMass;
             results.massBalance.imBalancedRxnBool = imBalancedRxnBool;
             results.massBalance.Elements = Elements;
             results.massBalance.missingFormulaeBool = missingFormulaeBool;
-            results.massBalance.balancedMetBool = balancedMetBool;            
+            results.massBalance.balancedMetBool = balancedMetBool;
         end
-        %Add the fields, if its either not a simple Check or if 
-        if doChargeBalance && (~simpleCheck || any(imBalancedCharge ~= 0))            
+        %Add the fields, if its either not a simple Check or if
+        if doChargeBalance && (~simpleCheck || any(imBalancedCharge ~= 0))
             results.chargeBalance = struct();
             results.chargeBalance.imBalanceCharge = imBalancedCharge;
         end
@@ -177,7 +177,7 @@
         [mins,maxs] = fluxVariability(model);
         %if this is not a simple check, or we have reactions which can't
         %carry flux.
-        if ~simpleCheck || any( ~(abs(mins) > 1e-12)| (abs(maxs) > 1e-12)) 
+        if ~simpleCheck || any( ~(abs(mins) > 1e-12)| (abs(maxs) > 1e-12))
             results.fluxConsistency = struct();
             if isfield(model,'rxns')
                 results.fluxConsistency.consistentReactions = model.rxns((abs(mins) > 1e-12) | (abs(maxs) > 1e-12));
@@ -185,7 +185,7 @@
             results.fluxConsistency.consistentReactionBool = ( (abs(mins) > 1e-12)| (abs(maxs) > 1e-12));
         end
     end
-
+    
 end
 
 if deadEndMetabolites
@@ -219,18 +219,18 @@
 if ~isempty(results) && ~silentCheck
     if isfield(results, 'Errors')
         problems = fieldnames(results.Errors);
-        disp('The following problems have been encountered in the model structure')
+        disp('verifyModel: The following problems have been encountered in the model structure')
         for i = 1:numel(problems)
             fprintf('%s:\n',problems{i})
-            problem_data = results.Errors.(problems{i});            
-            if isstruct(problem_data)             
+            problem_data = results.Errors.(problems{i});
+            if isstruct(problem_data)
                 problematic_fields = fieldnames(problem_data);
                 for field = 1: numel(problematic_fields)
-                    fprintf('%s: %s\n',problematic_fields{field}, results.Errors.(problems{i}).(problematic_fields{field}));                
+                    fprintf('%s: %s\n',problematic_fields{field}, results.Errors.(problems{i}).(problematic_fields{field}));
                 end
             else
                 for field = 1:numel(problem_data)
-                    fprintf('%s\n',problem_data{field});                
+                    fprintf('%s\n',problem_data{field});
                 end
             end
         end
@@ -248,67 +248,6 @@
 end
 
 
-function results = checkDatabaseIDs(model,results)
-% Checks the model for validity of database identifiers
-%
-% USAGE:
-%
-%    results = checkDatabaseIDs(model,results)
-%
-% INPUT:
-%    model:       a structure that represents the COBRA model.
-%    results:     the results structure for this test
-%
-% OUTPUT:
-%
-%    results:     a struct with problematic database ids added.
-%
-% .. Authors:
-%       - Thomas Pfau, May 2017
-
-dbMappings = getDefinedFieldProperties('Database',true);
-
-for i= 1:size(dbMappings,1)
-    if isfield(model,dbMappings{i,3})
-        fits = cellfun(@(x) isempty(x) || checkID(x,dbMappings{i,5}),model.(dbMappings{i,3}));
-        %Only add the something if we really have wrong IDs.
-        if any(~fits)
-            if ~isfield(results,'checkDatabaseIDs')
-                results.checkDatabaseIDs = struct();
-            end
-            if ~isfield(results.checkDatabaseIDs,'invalidIDs')
-                results.checkDatabaseIDs.invalidIDs = struct();
-            end
-            results.checkDatabaseIDs.invalidIDs.(dbMappings{i,3}) = cell(size(model.(dbMappings{i,3})));
-            results.checkDatabaseIDs.invalidIDs.(dbMappings{i,3})(:) = {'valid'};
-            results.checkDatabaseIDs.invalidIDs.(dbMappings{i,3})(~fits) = model.(dbMappings{i,3})(~fits);
-        end
-    end
-end
-end
-
-
-function accepted = checkID(id,pattern)
-% Checks the the given id(s), i.e. strings split by ; versus the pattern
-%
-% USAGE:
-%
-%    accepted = checkID(id,pattern)
-%
-% INPUT:
-%    id:          A String representing ids (potentially separated by ;)
-%    pattern:     The pattern to check the id(s) against.
-%
-% OUTPUT:
-%
-%    accepted:     Whether all ids are ok. 
-%
-% .. Authors:
-%       - Thomas Pfau, May 2017
-    ids = strsplit(id,';');
-    matches = regexp(ids,pattern);
-    accepted = all(~cellfun(@isempty,matches));    
-end
 
 function valid = checkFields(results,FieldNames,model)
 % Checks the given fields of the model in the results struct for
@@ -320,7 +259,7 @@
 %
 % INPUT:
 %    results:     the results structure for this test
-%    FieldNames:  The names of the fields to check. 
+%    FieldNames:  The names of the fields to check.
 %    model:       a structure that represents the COBRA model.
 %
 % OUTPUT:
@@ -341,98 +280,3 @@
 end
 
 
-function results = checkPresentFields(fieldProperties,model, results)
-% Check the model fields for consistency with the given fieldProperties and
-% update the results struct.
-%
-% USAGE:
-%
-%    results = checkPresentFields(fieldProperties,model, results)
-%
-% INPUT:
-%    fieldProperties:  field properties as obtained by
-%                      getDefinedFieldProperties
-%    model:            a structure that represents the COBRA model.
-%    results:          the results structure for this test
-%
-% OUTPUT:
-%
-%    results:          The updated results struct. 
-%
-% .. Authors:
-%       - Thomas Pfau, May 2017
-
-presentFields = find(ismember(fieldProperties(:,1),fieldnames(model)));
-
-%Check all Field Sizes
-for i = 1:numel(presentFields)
-    testedField = fieldProperties{presentFields(i),1};
-    [x_size,y_size] = size(model.(testedField));
-    xFieldMatch = fieldProperties{presentFields(i),2};
-    yFieldMatch = fieldProperties{presentFields(i),3};
-
-    checkX = ~isnan(xFieldMatch);
-    checkY = ~isnan(yFieldMatch);
-    if checkX
-        if ischar(xFieldMatch)
-            if ~isfield(model,xFieldMatch)
-                x_pres = 0;
-                if ~isfield(results.Errors,'missingFields')
-                    results.Errors.missingFields = {};
-                end
-                results.Errors.missingFields(end+1) = {xFieldMatch};
-            else
-                x_pres = numel(model.(xFieldMatch));
-            end
-            errorMessage = sprintf('%s: Size of %s does not match elements in %s', xFieldMatch,testedField,xFieldMatch);
-        elseif isnumeric(xFieldMatch)
-            errorMessage = sprintf('X Size of %s was %i. Expected %i',testedField, x_size,x_pres);
-            x_pres = xFieldMatch;
-        end
-        if x_pres ~= x_size
-            if ~isfield(results.Errors,'inconsistentFields')
-                results.Errors.inconsistentFields = struct();
-            end
-            results.Errors.inconsistentFields.(testedField) = errorMessage;
-        end
-    end
-    if checkY
-        if ischar(yFieldMatch)
-            if ~isfield(model,yFieldMatch)
-                y_pres = 0;
-                if ~isfield(results.Errors,'missingFields')
-                    results.Errors.missingFields = {};
-                end
-                results.Errors.missingFields(end+1) = {yFieldMatch};
-            else
-                y_pres = numel(model.(yFieldMatch));
-            end
-            errorMessage = sprintf('%s: Size of %s does not match elements in %s', yFieldMatch,testedField,yFieldMatch);
-        elseif isnumeric(yFieldMatch)
-            y_pres = yFieldMatch;
-            errorMessage = sprintf('Y Size of %s was %i. Expected %i',testedField, y_size,y_pres);
-        end
-        if y_pres ~= y_size
-            if ~isfield(results.Errors,'inconsistentFields')
-                results.Errors.inconsistentFields = struct();
-            end
-            results.Errors.inconsistentFields.(testedField) = errorMessage;
-        end
-    end
-    %Test the field content properties
-    %x is necessary here, since it is used for the eval below!
-    x = model.(testedField);
-    try
-        propertiesMatched = eval(fieldProperties{presentFields(i),4});
-    catch
-        propertiesMatched = false;
-    end
-    if ~propertiesMatched
-        if ~isfield(results.Errors,'inconsistentFields')
-            results.Errors.propertiesNotMatched = struct();
-        end
-        results.Errors.propertiesNotMatched.(testedField) = 'Field does not match the required properties';
-    end
-
-end
-end
diff --git a/src/reconstruction/refinement/generateRules.m b/src/reconstruction/refinement/generateRules.m
index 148a6132a0..cb10422de5 100644
--- a/src/reconstruction/refinement/generateRules.m
+++ b/src/reconstruction/refinement/generateRules.m
@@ -27,11 +27,13 @@
     if (~isfield(model, 'genes'))
         newGenes = allGenes;
     else
+%         C = setdiff(A,B) for vectors A and B, returns the values in A that 
+%         are not in B with no repetitions. C will be sorted.
         newGenes = setdiff(allGenes,model.genes);
     end
     if ~isempty(newGenes)
         if printLevel
-            warning('Found the following genes not present in the original model:\n%s\nAdding them to the model.',strjoin(newGenes,'\n'));
+            warning('Found the following genes in grRules that were not present in model.genes:\n%s\nAdding them to the model.',strjoin(newGenes,'\n'));
         end
         model = addGenes(model,newGenes);
     end