diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 24e4f03c..bd316eba 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -161,7 +161,7 @@ def BGPLVM_oil(optimize=True, N=200, Q=10, num_inducing=15, max_iters=150, plot=
     # optimize
     if optimize:
         m.constrain_fixed('noise')
-        m.optimize('scg', messages=1, max_iters=100, gtol=.05)
+        m.optimize('scg', messages=1, max_iters=200, gtol=.05)
         m.constrain_positive('noise')
         m.optimize('scg', messages=1, max_iters=max_iters, gtol=.05)
 
@@ -277,7 +277,6 @@ def bgplvm_simulation(optimize='scg',
     from GPy import kern
     reload(mrd); reload(kern)
 
-
     Y = Ylist[0]
 
     k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
@@ -286,7 +285,6 @@ def bgplvm_simulation(optimize='scg',
     # m.constrain('variance|noise', logexp_clipped())
     m['noise'] = Y.var() / 100.
 
-
     if optimize:
         print "Optimizing model:"
         m.optimize(optimize, max_iters=max_iters,
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 5cd90749..a09a94f3 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -66,12 +66,26 @@ class kern(Parameterized):
         Parameterized.setstate(self, state)
 
 
-    def plot_ARD(self, fignum=None, ax=None, title=None):
-        """If an ARD kernel is present, it bar-plots the ARD parameters"""
+    def plot_ARD(self, fignum=None, ax=None, title='', legend=False):
+        """If an ARD kernel is present, it bar-plots the ARD parameters,
+        :param fignum: figure number of the plot
+        :param ax: matplotlib axis to plot on
+        :param title: 
+            title of the plot, 
+            pass '' to not print a title
+            pass None for a generic title
+        """
         if ax is None:
             fig = pb.figure(fignum)
             ax = fig.add_subplot(111)
+        from GPy.util import Tango
+        from matplotlib.textpath import TextPath
+        Tango.reset()
+        xticklabels = []
+        bars = []
+        x0 = 0
         for p in self.parts:
+            c = Tango.nextMedium()
             if hasattr(p, 'ARD') and p.ARD:
                 if title is None:
                     ax.set_title('ARD parameters, %s kernel' % p.name)
@@ -82,10 +96,32 @@ class kern(Parameterized):
                 else:
                     ard_params = 1. / p.lengthscale
 
-                x = np.arange(len(ard_params))
-                ax.bar(x - 0.4, ard_params)
-                ax.set_xticks(x)
-                ax.set_xticklabels([r"${}$".format(i) for i in x])
+                x = np.arange(x0, x0 + len(ard_params))
+                bars.append(ax.bar(x, ard_params, align='center', color=c, edgecolor='k', linewidth=1.2, label=p.name))
+                xticklabels.extend([r"$\mathrm{{{name}}}\ {x}$".format(name=p.name, x=i) for i in np.arange(len(ard_params))])
+                x0 += len(ard_params)
+        x = np.arange(x0)
+        for bar in bars:
+            for patch, num in zip(bar.patches, np.arange(len(bar.patches))):
+                height = patch.get_height()
+                xi = patch.get_x() + patch.get_width() / 2.
+                va = 'top'
+                c = 'w'
+                t = TextPath((0, 0), "${xi}$".format(xi=xi), rotation=0, usetex=True, ha='center')
+                if patch.get_extents().height <= t.get_extents().height + 2:
+                    va = 'bottom'
+                    c = 'k'
+                ax.text(xi, height, "${xi}$".format(xi=int(num)), color=c, rotation=0, ha='center', va=va)
+        # for xi, t in zip(x, xticklabels):
+        #    ax.text(xi, maxi / 2, t, rotation=90, ha='center', va='center')
+        # ax.set_xticklabels(xticklabels, rotation=17)
+        ax.set_xticks([])
+        ax.set_xlim(-.5, x0 - .5)
+        if title is '':
+            ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
+                      ncol=max(2, len(bars)), mode="expand", borderaxespad=0.)
+        else:
+            ax.legend()
         return ax
 
     def _transform_gradients(self, g):
diff --git a/GPy/util/misc.py b/GPy/util/misc.py
index f931ae6a..0914adb8 100644
--- a/GPy/util/misc.py
+++ b/GPy/util/misc.py
@@ -59,7 +59,7 @@ def kmm_init(X, m = 10):
     return X[inducing]
 
 def fast_array_equal(A, B):
-    code="""
+    code2="""
     int i, j;
     return_val = 1;
 
@@ -74,6 +74,23 @@ def fast_array_equal(A, B):
     }
     """
 
+    code3="""
+    int i, j, z;
+    return_val = 1;
+
+    #pragma omp parallel for private(i, j, z)
+    for(i=0;i<N;i++){
+       for(j=0;j<D;j++){
+         for(z=0;z<Q;z++){
+            if(A(i, j, z) != B(i, j, z)){
+               return_val = 0;
+               break;
+            }
+          }
+       }
+    }
+    """
+
     support_code = """
     #include <omp.h>
     #include <math.h>
@@ -93,9 +110,14 @@ def fast_array_equal(A, B):
     elif A.shape == B.shape:
         if len(A.shape) == 2:
             N, D = A.shape
-            value = weave.inline(code, support_code=support_code, libraries=['gomp'],
+            value = weave.inline(code2, support_code=support_code, libraries=['gomp'],
                                  arg_names=['A', 'B', 'N', 'D'],
                                  type_converters=weave.converters.blitz,**weave_options)
+        elif len(A.shape) == 3:
+            N, D, Q = A.shape
+            value = weave.inline(code3, support_code=support_code, libraries=['gomp'],
+                                 arg_names=['A', 'B', 'N', 'D', 'Q'],
+                                 type_converters=weave.converters.blitz,**weave_options)
         else:
             value = np.array_equal(A,B)