diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 01a9ea41..60ee438c 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -22,6 +22,7 @@ class GP(GPBase):
 
     """
     def __init__(self, X, Y, kernel, likelihood, inference_method=None, name='gp'):
+        super(GP, self).__init__(X, likelihood, kernel, normalize_X=normalize_X, name=name)
 
         if inference_method is None:
             if isinstance(likelihood, likelihoods.Gaussian):
diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index 9466a011..ce5cd311 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -45,6 +45,7 @@ class GPBase(Model):
         
         self.add_parameter(self.kern, gradient=self.dL_dtheta_K)
         self.add_parameter(self.likelihood, gradient=lambda:self.posterior.dL_dtheta_lik)
+        
 
     def posterior_samples_f(self,X,size=10,which_parts='all',full_cov=True):
         """
diff --git a/GPy/core/index_operations.py b/GPy/core/index_operations.py
index 2e1fc774..ea7dfe2d 100644
--- a/GPy/core/index_operations.py
+++ b/GPy/core/index_operations.py
@@ -9,9 +9,12 @@ from parameter import Param
 from collections import defaultdict
 
 class ParamDict(defaultdict):
-    def __init__(self, default=lambda: numpy.array([], dtype=int)):
-        defaultdict.__init__(self, default)
-    
+    def __init__(self):
+        """
+        Default will be self._default, if not set otherwise
+        """
+        defaultdict.__init__(self, self.default_factory)
+        
     def __getitem__(self, key):
         try:
             return defaultdict.__getitem__(self, key)
@@ -35,7 +38,14 @@ class ParamDict(defaultdict):
                 if numpy.all(a==key) and a._parent_index_==key._parent_index_:
                     return super(ParamDict, self).__setitem__(a, value)
         defaultdict.__setitem__(self, key, value)
-        
+
+class SetDict(ParamDict):
+    def default_factory(self):
+        return set()
+
+class IntArrayDict(ParamDict):
+    def default_factory(self):
+        return numpy.int_([])
 
 class ParameterIndexOperations(object):
     '''
@@ -52,11 +62,11 @@ class ParameterIndexOperations(object):
         #self._reverse = collections.defaultdict(list)
         
     def __getstate__(self):
-        return self._properties, self._reverse
+        return self._properties#, self._reverse
         
     def __setstate__(self, state):
         self._properties = state[0]
-        self._reverse = state[1]
+        # self._reverse = state[1]
 
     def iteritems(self):
         return self._properties.iteritems()
diff --git a/GPy/core/parameter.py b/GPy/core/parameter.py
index 78e53d5b..ae7061e8 100644
--- a/GPy/core/parameter.py
+++ b/GPy/core/parameter.py
@@ -6,8 +6,8 @@ Created on 4 Sep 2013
 import itertools
 import numpy
 from transformations import Logexp, NegativeLogexp, Logistic
-from parameterized import Nameable, Pickleable, Observable
-from GPy.core.parameterized import _adjust_name_for_printing
+from parameterized import Nameable, Pickleable, Observable, Parentable
+from parameterized import _adjust_name_for_printing
 
 ###### printing
 __constraints_name__ = "Constraint"
@@ -60,8 +60,95 @@ class ObservableArray(ListArray, Observable):
     def __setslice__(self, start, stop, val):
         return self.__setitem__(slice(start, stop), val)    
 
-    
-class Param(ObservableArray, Nameable, Pickleable):
+class Constrainable(Nameable):
+    def __init__(self, name):
+        super(Constrainable,self).__init__(name)
+    #===========================================================================
+    # Constrain operations -> done
+    #===========================================================================
+    def constrain(self, transform, warning=True, update=True):
+        """
+        :param transform: the :py:class:`GPy.core.transformations.Transformation`
+                          to constrain the this parameter to.
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain the parameter to the given
+        :py:class:`GPy.core.transformations.Transformation`.
+        """
+        if self.has_parent():
+            self._highest_parent_._add_constrain(self, transform, warning)
+            if update:
+                self._highest_parent_.parameters_changed()
+        else:
+            for p in self._parameters_:
+                self._add_constrain(p, transform, warning)
+            if update:
+                self.parameters_changed()
+    def constrain_positive(self, warning=True):
+        """
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain this parameter to the default positive constraint.
+        """
+        self.constrain(Logexp(), warning)
+
+    def constrain_negative(self, warning=True):
+        """
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain this parameter to the default negative constraint.
+        """
+        self.constrain(NegativeLogexp(), warning)
+
+    def constrain_bounded(self, lower, upper, warning=True):
+        """
+        :param lower, upper: the limits to bound this parameter to
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain this parameter to lie within the given range.
+        """
+        self.constrain(Logistic(lower, upper), warning)
+
+    def unconstrain(self, *transforms):
+        """
+        :param transforms: The transformations to unconstrain from.
+        
+        remove all :py:class:`GPy.core.transformations.Transformation` 
+        transformats of this parameter object.
+        """
+        if self.has_parent():
+            self._highest_parent_._remove_constrain(self, *transforms)
+        else:
+            for p in self._parameters_:
+                self._remove_constrain(p, *transforms)
+
+    def unconstrain_positive(self):
+        """
+        Remove positive constraint of this parameter. 
+        """
+        self.unconstrain(Logexp())
+
+    def unconstrain_negative(self):
+        """
+        Remove negative constraint of this parameter. 
+        """
+        self.unconstrain(NegativeLogexp())
+
+    def unconstrain_bounded(self, lower, upper):
+        """
+        :param lower, upper: the limits to unbound this parameter from
+        
+        Remove (lower, upper) bounded constrain from this parameter/
+        """
+        self.unconstrain(Logistic(lower, upper))
+
+class Float(numpy.float64, Constrainable):
+    def __init__(self, f, base):
+        super(Float,self).__init__(f)
+        self._base = base
+        
+        
+class Param(ObservableArray, Constrainable):
     """
     Parameter object for GPy models.
 
@@ -82,23 +169,29 @@ class Param(ObservableArray, Nameable, Pickleable):
     Fixing parameters will fix them to the value they are right now. If you change
     the fixed value, it will be fixed to the new value!
     
-    See :py:class:`GPy.core.parameterized.Parameterized` for more details.
+    See :py:class:`GPy.core.parameterized.Parameterized` for more details on constraining etc.
+
+    This ndarray can be stored in lists and checked if it is in.
+
+    >>> import numpy as np
+    >>> x = np.random.normal(size=(10,3))
+    >>> x in [[1], x, [3]]
+    True
+    
+    WARNING: This overrides the functionality of x==y!!!
+    Use numpy.equal(x,y) for element-wise equality testing.
     """
-    __array_priority__ = -numpy.inf # Never give back Param
+    __array_priority__ = 0 # Never give back Param
     _fixes_ = None
     def __new__(cls, name, input_array, *args, **kwargs):
         obj = numpy.atleast_1d(super(Param, cls).__new__(cls, input_array=input_array))
-        obj._direct_parent_ = None
-        #obj.name = name
-        obj._parent_index_ = None
-        obj._highest_parent_ = None
         obj._current_slice_ = (slice(obj.shape[0]),)
         obj._realshape_ = obj.shape
         obj._realsize_ = obj.size
         obj._realndim_ = obj.ndim
         obj._updated_ = False
-        from index_operations import ParamDict
-        obj._tied_to_me_ = ParamDict(set)
+        from index_operations import SetDict
+        obj._tied_to_me_ = SetDict()
         obj._tied_to_ = []
         obj._original_ = True
         return obj
@@ -110,11 +203,10 @@ class Param(ObservableArray, Nameable, Pickleable):
         # see InfoArray.__array_finalize__ for comments
         if obj is None: return
         super(Param, self).__array_finalize__(obj)
-        self.name = getattr(obj, 'name', None)
-        self._current_slice_ = getattr(obj, '_current_slice_', None)
         self._direct_parent_ = getattr(obj, '_direct_parent_', None)
         self._parent_index_ = getattr(obj, '_parent_index_', None)
         self._highest_parent_ = getattr(obj, '_highest_parent_', None)
+        self._current_slice_ = getattr(obj, '_current_slice_', None)
         self._tied_to_me_ = getattr(obj, '_tied_to_me_', None)
         self._tied_to_ = getattr(obj, '_tied_to_', None)
         self._realshape_ = getattr(obj, '_realshape_', None)
@@ -122,13 +214,14 @@ class Param(ObservableArray, Nameable, Pickleable):
         self._realndim_ = getattr(obj, '_realndim_', None)
         self._updated_ = getattr(obj, '_updated_', None)
         self._original_ = getattr(obj, '_original_', None)
+        self._name = getattr(obj, 'name', None)
         
     def __array_wrap__(self, out_arr, context=None):
         return out_arr.view(numpy.ndarray)
     #===========================================================================
     # Pickling operations
     #===========================================================================
-    def __reduce__(self):
+    def __reduce_ex__(self):
         func, args, state = super(Param, self).__reduce__()
         return func, args, (state, 
                             (self.name,
@@ -207,83 +300,7 @@ class Param(ObservableArray, Nameable, Pickleable):
         self._highest_parent_._unfix(self)
     unfix = unconstrain_fixed
     #===========================================================================
-    # Constrain operations -> done
-    #===========================================================================
-    def constrain(self, transform, warning=True, update=True):
-        """
-        :param transform: the :py:class:`GPy.core.transformations.Transformation`
-                          to constrain the this parameter to.
-        :param warning: print a warning if re-constraining parameters.
-        
-        Constrain the parameter to the given
-        :py:class:`GPy.core.transformations.Transformation`.
-        """
-        if self._original_: # this happens when indexing created a copy of the array
-            self.__setitem__(slice(None), transform.initialize(self), update=False)
-        else:
-            self._direct_parent_._get_original(self).__setitem__(self._current_slice_, transform.initialize(self), update=False)
-        self._highest_parent_._add_constrain(self, transform, warning)
-        for t in self._tied_to_me_.iterkeys():
-            if transform not in self._highest_parent_._constraints_for_collect(t, t._raveled_index()):
-                t._direct_parent_._get_original(t)[t._current_slice_].constrain(transform, warning, update)
-        if update:
-            self._highest_parent_.parameters_changed()        
-
-    def constrain_positive(self, warning=True):
-        """
-        :param warning: print a warning if re-constraining parameters.
-        
-        Constrain this parameter to the default positive constraint.
-        """
-        self.constrain(Logexp(), warning)
-
-    def constrain_negative(self, warning=True):
-        """
-        :param warning: print a warning if re-constraining parameters.
-        
-        Constrain this parameter to the default negative constraint.
-        """
-        self.constrain(NegativeLogexp(), warning)
-
-    def constrain_bounded(self, lower, upper, warning=True):
-        """
-        :param lower, upper: the limits to bound this parameter to
-        :param warning: print a warning if re-constraining parameters.
-        
-        Constrain this parameter to lie within the given range.
-        """
-        self.constrain(Logistic(lower, upper), warning)
-
-    def unconstrain(self, *transforms):
-        """
-        :param transforms: The transformations to unconstrain from.
-        
-        remove all :py:class:`GPy.core.transformations.Transformation` 
-        transformats of this parameter object.
-        """
-        self._highest_parent_._remove_constrain(self, *transforms)
-
-    def unconstrain_positive(self):
-        """
-        Remove positive constraint of this parameter. 
-        """
-        self.unconstrain(Logexp())
-
-    def unconstrain_negative(self):
-        """
-        Remove negative constraint of this parameter. 
-        """
-        self.unconstrain(NegativeLogexp())
-
-    def unconstrain_bounded(self, lower, upper):
-        """
-        :param lower, upper: the limits to unbound this parameter from
-        
-        Remove (lower, upper) bounded constrain from this parameter/
-        """
-        self.unconstrain(Logistic(lower, upper))
-    #===========================================================================
-    # Tying operations -> done
+    # Tying operations -> bugged, TODO
     #===========================================================================
     def tie_to(self, param):
         """
@@ -667,6 +684,8 @@ class ParamConcatenation(object):
         return "\n".join(map(repr,self.params))
     
 if __name__ == '__main__':
+    
+
     from GPy.core.parameterized import Parameterized
     from GPy.core.parameter import Param
 
diff --git a/GPy/core/parameterized.py b/GPy/core/parameterized.py
index 36ae547e..67bc18f5 100644
--- a/GPy/core/parameterized.py
+++ b/GPy/core/parameterized.py
@@ -11,17 +11,21 @@ from re import compile, _pattern_type
 import re
 
 class Parentable(object):
-    _direct_parent_ = None
-    _parent_index_ = None
-    
+    def __init__(self, direct_parent=None, highest_parent=None, parent_index=None):
+        super(Parentable,self).__init__()        
+        self._direct_parent_ = direct_parent
+        self._parent_index_ = parent_index
+        self._highest_parent_ = highest_parent
+        
     def has_parent(self):
         return self._direct_parent_ is not None
     
 class Nameable(Parentable):
     _name = None
-    def __init__(self, name):
+    def __init__(self, name, direct_parent=None, highest_parent=None, parent_index=None):
+        super(Nameable,self).__init__(direct_parent, highest_parent, parent_index)
         self._name = name or self.__class__.__name__
-        self.name = name
+
     @property
     def name(self):
         return self._name
@@ -68,7 +72,7 @@ def _adjust_name_for_printing(name):
     if name is not None:
         return name.replace(" ", "_").replace(".", "_").replace("-","").replace("+","").replace("!","").replace("*","").replace("/","")
     return ''
-from parameter import ParamConcatenation, Param
+from parameter import ParamConcatenation, Param, Constrainable
 from index_operations import ParameterIndexOperations,\
     index_empty
 
@@ -83,7 +87,7 @@ FIXED = False
 UNFIXED = True
 #===============================================================================
 
-class Parameterized(Nameable, Pickleable, Observable):
+class Parameterized(Constrainable, Pickleable, Observable):
     """
     Parameterized class
     
@@ -92,10 +96,9 @@ class Parameterized(Nameable, Pickleable, Observable):
     Printing parameters:
     
         - print m:           prints a nice summary over all parameters
-        - print m.name:      prints details for all the parameters 
-                             which start with name
-        - print m['.*name']: prints details for all the parameters 
-                             which contain "name"
+        - print m.name:      prints details for parameter with name 'name'
+        - print m[regexp]: prints details for all the parameters 
+                             which match (!) regexp
         - print m['']:       prints details for all parameters
     
         Fields:
@@ -108,11 +111,10 @@ class Parameterized(Nameable, Pickleable, Observable):
             Tied_to:    which paramter it is tied to.
     
     Getting and setting parameters:
-        
-        Two ways to get parameters:
-            
-            - m.name regular expression matches all parameters beginning with name
-            - m['name'] regular expression matches all parameters with name
+
+        Set all values in parameter to one:
+
+            m.name.to.parameter = 1
     
     Handling of constraining, fixing and tieing parameters:
         
@@ -120,15 +122,15 @@ class Parameterized(Nameable, Pickleable, Observable):
         
             - m.name[:,1].constrain_positive()
             - m.name[0].tie_to(m.name[1])
-            
+
         Fixing parameters will fix them to the value they are right now. If you change
         the parameters value, the parameter will be fixed to the new value!
-        
+
         If you want to operate on all parameters use m[''] to wildcard select all paramters 
         and concatenate them. Printing m[''] will result in printing of all parameters in detail.
     """
     def __init__(self, name=None):
-        super(Parameterized, self).__init__(name)
+        super(Parameterized, self).__init__(name=name)
         self._in_init_ = True
         self._constraints_ = None#ParameterIndexOperations()
         if not hasattr(self, "_parameters_"):
@@ -189,7 +191,7 @@ class Parameterized(Nameable, Pickleable, Observable):
         
         
         Add all parameters to this parameter class, you can insert parameters 
-        at any given index using the :py:func:`list.insert` syntax 
+        at any given index using the :func:`list.insert` syntax 
         """
         if parameter in self._parameters_ and index is not None:
             # make sure fixes and constraints are indexed right
@@ -279,10 +281,8 @@ class Parameterized(Nameable, Pickleable, Observable):
             return
         i = 0
         sizes = [0]
-        #self.size = sum(p.size for p in self._parameters_)
         self._param_slices_ = []
         for p in self._parameters_:
-            #if p._parent_ is None:
             p._direct_parent_ = self
             p._parent_index_ = i
             i += 1
@@ -291,10 +291,8 @@ class Parameterized(Nameable, Pickleable, Observable):
             not_unique = []
             sizes.append(p.size+sizes[-1])
             self._param_slices_.append(slice(sizes[-2], sizes[-1]))
-#             if p._fixes_ is not None:
-#                 self._fixes_[p._raveled_index_for(p)] = p._fixes_
-#                 p._fixes_ = None
             pname = _adjust_name_for_printing(p.name)  
+            # and makes sure to not delete programmatically added parameters
             if pname in self.__dict__:
                 if isinstance(self.__dict__[pname], (Parameterized, Param)):
                     if not p is self.__dict__[pname]:
@@ -303,16 +301,6 @@ class Parameterized(Nameable, Pickleable, Observable):
             elif not (pname in not_unique):
                 self.__dict__[pname] = p
                 self._added_names_.add(pname)
-#         for p in self._parameters_:
-#             if hasattr(p, '_constraints_') and p._constraints_ is not None:
-#                 for t, ind in p._constraints_.iteritems():
-#                     self.constraints.add(t, ind+self._offset_for(p))
-#                 p._constraints_.clear()
-#         if np.all(self._fixes_): # ==UNFIXED
-#             self._fixes_= None
-#         else:
-#             self.constraints.add(__fixed__, np.nonzero(~self._fixes_)[0])
-#         self.parameters_changed()
     #===========================================================================
     # Pickling operations
     #===========================================================================
@@ -355,20 +343,18 @@ class Parameterized(Nameable, Pickleable, Observable):
         return [
                 self._fixes_,
                 self._constraints_,
-                self._priors_,
                 self._parameters_,
                 self._name,
-                self.gradient_mapping,
+                #self.gradient_mapping,
                 self._added_names_,
                 ]
         
     def setstate(self, state):
         self._added_names_ = state.pop()
-        self.gradient_mapping = state.pop(),
+        #self.gradient_mapping = state.pop()
         self._name = state.pop()
         self._parameters_ = state.pop()
         self._connect_parameters()
-        self._priors = state.pop()
         self._constraints_ = state.pop()
         self._fixes_ = state.pop()
         self.parameters_changed()
@@ -379,20 +365,10 @@ class Parameterized(Nameable, Pickleable, Observable):
         if self.has_parent():
             return g
         x = self._get_params()
-        #g = g.copy()
-        #for constraint, index in self.constraints.iteritems():
-        #    if constraint != __fixed__:
-        #        g[index] = g[index] * constraint.gradfactor(x[index])
         [numpy.put(g, i, g[i]*c.gradfactor(x[i])) for c,i in self.constraints.iteritems() if c != __fixed__]
-        #[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
         for p in self.flattened_parameters:
             for t,i in p._tied_to_me_.iteritems():
                 g[self._offset_for(p) + numpy.array(list(i))] += g[self._raveled_index_for(t)]
-        #[g[self._offset_for(t) + numpy.array(list(i))].__iadd__(v) for i, v in [[i, g[self._raveled_index_for(p)].sum()] for p in self.flattened_parameters for t,i in p._tied_to_me_.iteritems()]]
-#         if len(self.tied_indices) or len(self.fixed_indices):
-#             to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
-#             return np.delete(g, to_remove)
-#         else:
         if self._has_fixes(): return g[self._fixes_]
         return g
     #===========================================================================
@@ -405,18 +381,20 @@ class Parameterized(Nameable, Pickleable, Observable):
         return n
     def _get_params(self):
         # don't overwrite this anymore!
-        return numpy.hstack([x._get_params() for x in self._parameters_])#numpy.fromiter(itertools.chain(*itertools.imap(lambda x: x._get_params(), self._parameters_)), dtype=numpy.float64, count=sum(self._parameter_sizes_))    
+        return numpy.hstack([x._get_params() for x in self._parameters_])
     def _set_params(self, params, update=True):
         # don't overwrite this anymore!
         [p._set_params(params[s], update=update) for p,s in itertools.izip(self._parameters_,self._param_slices_)]
         self.parameters_changed()
     def _get_params_transformed(self):
+        # transformed parameters (apply transformation rules)
         p = self._get_params()
         [numpy.put(p, ind, c.finv(p[ind])) for c,ind in self.constraints.iteritems() if c != __fixed__]
         if self._has_fixes():
             return p[self._fixes_]
         return p
     def _set_params_transformed(self, p):
+        # inverse apply transformations for parameters and set the resulting parameters
         p = p.copy()
         if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp
         [numpy.put(p, ind, c.f(p[ind])) for c,ind in self.constraints.iteritems() if c != __fixed__]
@@ -448,9 +426,18 @@ class Parameterized(Nameable, Pickleable, Observable):
         return 0
     
     def _raveled_index_for(self, param):
+        """
+        get the raveled index for a parameter
+        that is an int array, containing the indexes for the flattened
+        parameter inside this parameterized logic.
+        """
         return param._raveled_index() + self._offset_for(param)
     
     def _raveled_index(self):
+        """
+        get the raveled index for this object,
+        this is not in the global view of things!
+        """
         return numpy.r_[:self.size]
     #===========================================================================
     # Handle ties:
@@ -519,6 +506,8 @@ class Parameterized(Nameable, Pickleable, Observable):
         reconstrained = self._remove_constrain(param, index=rav_i) # remove constraints before
         # if removing constraints before adding new is not wanted, just delete the above line!        
         self.constraints.add(transform, rav_i)
+        param = self._get_original(param)
+        param._set_params(transform.initialize(param._get_params()))
         if warning and any(reconstrained):
             # if you want to print the whole params object, which was reconstrained use:
             # m = str(param[self._backtranslate_index(param, reconstrained)])
@@ -610,6 +599,12 @@ class Parameterized(Nameable, Pickleable, Observable):
     #===========================================================================
     # Printing:
     #===========================================================================
+    def _short(self):
+        # short string to print
+        if self.has_parent():
+            return self._direct_parent_.hirarchy_name() + _adjust_name_for_printing(self.name)
+        else:
+            return _adjust_name_for_printing(self.name)
     def _parameter_names(self, add_name=False):
         if add_name:
             return [_adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x._parameter_names(add_name=True)]
@@ -639,9 +634,11 @@ class Parameterized(Nameable, Pickleable, Observable):
     def _ties_str(self):
         return [','.join(x._ties_str) for x in self.flattened_parameters]
     def __str__(self, header=True):
+        
+        name  = _adjust_name_for_printing(self.name) + "."
         constrs = self._constraints_str; ts = self._ties_str
         desc = self._description_str; names = self.parameter_names
-        nl = max([len(str(x)) for x in names + [_adjust_name_for_printing(self.name)]])
+        nl = max([len(str(x)) for x in names + [name]])
         sl = max([len(str(x)) for x in desc + ["Value"]])
         cl = max([len(str(x)) if x else 0 for x in constrs  + ["Constraint"]])
         tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]])
@@ -652,7 +649,7 @@ class Parameterized(Nameable, Pickleable, Observable):
         #to_print = [format_spec.format(p=p, const=c, t=t) if isinstance(p, Param) else p.__str__(header=False) for p, c, t in itertools.izip(self._parameters_, constrs, ts)]
         sep = '-'*(nl+sl+cl+tl+8*2+3)
         if header:
-            header = "  {{0:<{0}s}}  |  {{1:^{1}s}}  |  {{2:^{2}s}}  |  {{3:^{3}s}}".format(nl, sl, cl, tl).format(_adjust_name_for_printing(self.name), "Value", "Constraint", "Tied to")
+            header = "  {{0:<{0}s}}  |  {{1:^{1}s}}  |  {{2:^{2}s}}  |  {{3:^{3}s}}".format(nl, sl, cl, tl).format(name, "Value", "Constraint", "Tied to")
             #header += '\n' + sep
             to_print.insert(0, header)
         return '\n'.format(sep).join(to_print)
diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py
index a7be7334..2d15ff0c 100644
--- a/GPy/core/sparse_gp.py
+++ b/GPy/core/sparse_gp.py
@@ -29,8 +29,8 @@ class SparseGP(GPBase):
 
     """
 
-    def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False):
-        GPBase.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, name="sparse GP")
+    def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False, name='sparse gp'):
+        GPBase.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, name=name)
 
         self.Z = Z
         self.num_inducing = Z.shape[0]
diff --git a/GPy/core/variational.py b/GPy/core/variational.py
index 74287dcf..510d36da 100644
--- a/GPy/core/variational.py
+++ b/GPy/core/variational.py
@@ -3,8 +3,10 @@ Created on 6 Nov 2013
 
 @author: maxz
 '''
+import numpy as np
 from parameterized import Parameterized
 from parameter import Param
+from ..util.misc import param_to_array
 
 class Normal(Parameterized):
     '''
@@ -12,8 +14,53 @@ class Normal(Parameterized):
     
     holds the means and variances for a factorizing multivariate normal distribution
     '''
-    def __init__(self, name, means, variances):
+    def __init__(self, means, variances, name='latent space'):
         Parameterized.__init__(self, name=name)
         self.means = Param("mean", means)
         self.variances = Param('variance', variances)
-        self.add_parameters(self.means, self.variances)
\ No newline at end of file
+        self.add_parameters(self.means, self.variances)
+
+    def plot(self, fignum=None, ax=None, colors=None):
+        """
+        Plot latent space X in 1D:
+
+            - if fig is given, create input_dim subplots in fig and plot in these
+            - if ax is given plot input_dim 1D latent space plots of X into each `axis`
+            - if neither fig nor ax is given create a figure with fignum and plot in there
+
+        colors:
+            colors of different latent space dimensions input_dim
+
+        """
+        import pylab
+        if ax is None:
+            fig = pylab.figure(num=fignum, figsize=(8, min(12, (2 * self.means.shape[1]))))
+        if colors is None:
+            colors = pylab.gca()._get_lines.color_cycle
+            pylab.clf()
+        else:
+            colors = iter(colors)
+        plots = []
+        means, variances = param_to_array(self.means, self.variances)
+        x = np.arange(means.shape[0])
+        for i in range(means.shape[1]):
+            if ax is None:
+                a = fig.add_subplot(means.shape[1], 1, i + 1)
+            elif isinstance(ax, (tuple, list)):
+                a = ax[i]
+            else:
+                raise ValueError("Need one ax per latent dimnesion input_dim")
+            a.plot(means, c='k', alpha=.3)
+            plots.extend(a.plot(x, means.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
+            a.fill_between(x,
+                            means.T[i] - 2 * np.sqrt(variances.T[i]),
+                            means.T[i] + 2 * np.sqrt(variances.T[i]),
+                            facecolor=plots[-1].get_color(),
+                            alpha=.3)
+            a.legend(borderaxespad=0.)
+            a.set_xlim(x.min(), x.max())
+            if i < means.shape[1] - 1:
+                a.set_xticklabels('')
+        pylab.draw()
+        fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
+        return fig
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 4c286a9e..badef423 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -285,6 +285,7 @@ def bgplvm_simulation(optimize='scg',
     k = kern.linear(input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2)) # + kern.bias(input_dim)
     m = BayesianGPLVM(Y, input_dim, init="PCA", num_inducing=num_inducing, kernel=k)
 
+    import ipdb; ipdb.set_trace()
     # m.constrain('variance|noise', LogexpClipped())
     m['gaussian'] = Y.var() / 100.
 
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 2d3a78c6..782cfa1b 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -59,7 +59,7 @@ class kern(Parameterized):
         Get the current state of the class,
         here just all the indices, rest can get recomputed
         """
-        return Parameterized.getstate(self) + [self._parameters_,
+        return Parameterized.getstate(self) + [#self._parameters_,
                 #self.num_params,
                 self.input_dim,
                 self.input_slices,
@@ -71,7 +71,7 @@ class kern(Parameterized):
         self.input_slices = state.pop()
         self.input_dim = state.pop()
         #self.num_params = state.pop()
-        self._parameters_ = state.pop()
+        #self._parameters_ = state.pop()
         Parameterized.setstate(self, state)
 
 
diff --git a/GPy/kern/parts/linear.py b/GPy/kern/parts/linear.py
index c8b2f618..c3bbf837 100644
--- a/GPy/kern/parts/linear.py
+++ b/GPy/kern/parts/linear.py
@@ -2,12 +2,12 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from kernpart import Kernpart
 import numpy as np
-from ...util.linalg import tdot
-from ...util.misc import fast_array_equal
 from scipy import weave
-from GPy.core.parameter import Param
+from kernpart import Kernpart
+from ...util.linalg import tdot
+from ...util.misc import fast_array_equal, param_to_array
+from ...core.parameter import Param
 
 class Linear(Kernpart):
     """
@@ -235,8 +235,8 @@ class Linear(Kernpart):
         weave_options = {'headers'           : ['<omp.h>'],
                          'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
                          'extra_link_args'   : ['-lgomp']}
-
-        N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
+        
+        N,num_inducing,input_dim,mu = mu.shape[0],Z.shape[0],mu.shape[1],param_to_array(mu)
         weave.inline(code, support_code=support_code, libraries=['gomp'],
                      arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
@@ -271,6 +271,7 @@ class Linear(Kernpart):
                          'extra_link_args'   : ['-lgomp']}
 
         N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
+        mu, AZA, target, dL_dpsi2 = param_to_array(mu, AZA, target, dL_dpsi2)
         weave.inline(code, support_code=support_code, libraries=['gomp'],
                      arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py
index 8bd95dc5..c6e827b0 100644
--- a/GPy/kern/parts/rbf.py
+++ b/GPy/kern/parts/rbf.py
@@ -2,12 +2,12 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from kernpart import Kernpart
 import numpy as np
 from scipy import weave
+from kernpart import Kernpart
 from ...util.linalg import tdot
-from ...util.misc import fast_array_equal
-from GPy.core.parameter import Param
+from ...util.misc import fast_array_equal, param_to_array
+from ...core.parameter import Param
 
 class RBF(Kernpart):
     """
@@ -63,10 +63,11 @@ class RBF(Kernpart):
         #self._X, self._X2, self._params_save = np.empty(shape=(3, 1))
 
         # a set of optional args to pass to weave
-        self.weave_options = {'headers'           : ['<omp.h>'],
-                              'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
-                              'extra_link_args'   : ['-lgomp']}
-    
+        # self.weave_options = {'headers'           : ['<omp.h>'],
+        #                       'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
+        #                       'extra_link_args'   : ['-lgomp']}
+        self.weave_options = {}
+        
     def on_input_change(self, X):
         #self._K_computations(X, None)
         pass
@@ -133,7 +134,8 @@ class RBF(Kernpart):
                 }
                 """
                 num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim
-                weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
+                X = param_to_array(X)
+                weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
             else:
                 code = """
                 int q,i,j;
@@ -150,6 +152,7 @@ class RBF(Kernpart):
                 """
                 num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim
                 # [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)]
+                X, X2 = param_to_array(X, X2)
                 weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
         else:
             target[1] += (self.variance / self.lengthscale) * np.sum(self._K_dvar * self._K_dist2 * dL_dK)
diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py
index 9f19393b..5223397f 100644
--- a/GPy/models/bayesian_gplvm.py
+++ b/GPy/models/bayesian_gplvm.py
@@ -27,7 +27,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
 
     """
     def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
-                 Z=None, kernel=None, **kwargs):
+                 Z=None, kernel=None, name='bayesian gplvm', **kwargs):
         if type(likelihood_or_Y) is np.ndarray:
             likelihood = Gaussian(likelihood_or_Y)
         else:
@@ -47,9 +47,8 @@ class BayesianGPLVM(SparseGP, GPLVM):
         if kernel is None:
             kernel = kern.rbf(input_dim) # + kern.white(input_dim)
 
-        SparseGP.__init__(self, X, likelihood, kernel, Z=Z, X_variance=X_variance, **kwargs)
-        
-        self.q = Normal('latent space', self.X, self.X_variance)
+        SparseGP.__init__(self, X=X, likelihood=likelihood, kernel=kernel, Z=Z, X_variance=X_variance, name=name, **kwargs)
+        self.q = Normal(self.X, self.X_variance)
         self.add_parameter(self.q, gradient=self._dbound_dmuS, index=0)
         self.ensure_default_constraints()
 
@@ -252,50 +251,6 @@ class BayesianGPLVM(SparseGP, GPLVM):
                 controller.deactivate()
         return controller.view
 
-    def plot_X_1d(self, fignum=None, ax=None, colors=None):
-        """
-        Plot latent space X in 1D:
-
-            - if fig is given, create input_dim subplots in fig and plot in these
-            - if ax is given plot input_dim 1D latent space plots of X into each `axis`
-            - if neither fig nor ax is given create a figure with fignum and plot in there
-
-        colors:
-            colors of different latent space dimensions input_dim
-
-        """
-        import pylab
-        if ax is None:
-            fig = pylab.figure(num=fignum, figsize=(8, min(12, (2 * self.X.shape[1]))))
-        if colors is None:
-            colors = pylab.gca()._get_lines.color_cycle
-            pylab.clf()
-        else:
-            colors = iter(colors)
-        plots = []
-        x = np.arange(self.X.shape[0])
-        for i in range(self.X.shape[1]):
-            if ax is None:
-                a = fig.add_subplot(self.X.shape[1], 1, i + 1)
-            elif isinstance(ax, (tuple, list)):
-                a = ax[i]
-            else:
-                raise ValueError("Need one ax per latent dimnesion input_dim")
-            a.plot(self.X, c='k', alpha=.3)
-            plots.extend(a.plot(x, self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
-            a.fill_between(x,
-                            self.X.T[i] - 2 * np.sqrt(self.X_variance.T[i]),
-                            self.X.T[i] + 2 * np.sqrt(self.X_variance.T[i]),
-                            facecolor=plots[-1].get_color(),
-                            alpha=.3)
-            a.legend(borderaxespad=0.)
-            a.set_xlim(x.min(), x.max())
-            if i < self.X.shape[1] - 1:
-                a.set_xticklabels('')
-        pylab.draw()
-        fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
-        return fig
-
 def latent_cost_and_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
     """
     objective function for fitting the latent variables for test points
diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py
index 28367be9..6b5ac07f 100644
--- a/GPy/models/gplvm.py
+++ b/GPy/models/gplvm.py
@@ -28,15 +28,15 @@ class GPLVM(GP):
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, Y, input_dim, init='PCA', X=None, kernel=None, normalize_Y=False):
+    def __init__(self, Y, input_dim, init='PCA', X=None, kernel=None, normalize_Y=False, name="gplvm"):
         if X is None:
             X = self.initialise_latent(init, input_dim, Y)
         if kernel is None:
             kernel = kern.rbf(input_dim, ARD=input_dim > 1) + kern.bias(input_dim, np.exp(-2))
         likelihood = Gaussian(Y, normalize=normalize_Y, variance=np.exp(-2.))
-        GP.__init__(self, X, likelihood, kernel, normalize_X=False)
+        GP.__init__(self, X, likelihood, kernel, normalize_X=False, name=name)
         self.X = Param('q_mean', self.X)
-        self.add_parameter(self.X, self.dK_dX, 0)
+        self.add_parameter(self.X, gradient=self.dK_dX, index=0)
         #self.set_prior('.*X', Gaussian_prior(0, 1))
         self.ensure_default_constraints()
 
diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py
index 8e1471d2..2b9ef805 100644
--- a/GPy/models/mrd.py
+++ b/GPy/models/mrd.py
@@ -323,9 +323,6 @@ class MRD(Model):
         else:
             return pylab.gcf()
 
-    def plot_X_1d(self, *a, **kw):
-        return self.gref.plot_X_1d(*a, **kw)
-
     def plot_X(self, fignum=None, ax=None):
         fig = self._handle_plotting(fignum, ax, lambda i, g, ax: ax.imshow(g.X))
         return fig
diff --git a/GPy/util/misc.py b/GPy/util/misc.py
index 1cb4c182..8cf5806f 100644
--- a/GPy/util/misc.py
+++ b/GPy/util/misc.py
@@ -32,6 +32,18 @@ def chain_3(d3f_dg3, dg_dx, d2f_dg2, d2g_dx2, df_dg, d3g_dx3):
     """
     return d3f_dg3*(dg_dx**3) + 3*d2f_dg2*dg_dx*d2g_dx2 + df_dg*d3g_dx3
 
+### make a parameter to its corresponding array:
+def param_to_array(*param):
+    """
+    Convert an arbitrary number of parameters to :class:ndarray class objects. This is for
+    converting parameter objects to numpy arrays, when using scipy.weave.inline routine.
+    In scipy.weave.blitz there is no automatic array detection (even when the array inherits
+    from :class:ndarray)"""
+    assert len(param) > 0, "At least one parameter needed"
+    if len(param) == 1:
+        return param[0].view(np.ndarray)
+    return map(lambda x: x.view(np.ndarray), param)
+
 def opt_wrapper(m, **kwargs):
     """
     This function just wraps the optimization procedure of a GPy
@@ -159,6 +171,7 @@ def fast_array_equal(A, B):
     elif ((A == None) and (B != None)) or ((A != None) and (B == None)):
         return False
     elif A.shape == B.shape:
+        A, B = param_to_array(A, B)
         if A.ndim == 2:
             N, D = [int(i) for i in A.shape]
             value = weave.inline(code2, support_code=support_code,
@@ -174,7 +187,6 @@ def fast_array_equal(A, B):
 
     return value
 
-
 if __name__ == '__main__':
     import pylab as plt
     X = np.linspace(1,10, 100)[:, None]
diff --git a/GPy/util/plot_latent.py b/GPy/util/plot_latent.py
index 62442650..cecb811c 100644
--- a/GPy/util/plot_latent.py
+++ b/GPy/util/plot_latent.py
@@ -38,9 +38,11 @@ def plot_latent(model, labels=None, which_indices=None,
 
     input_1, input_2 = most_significant_input_dimensions(model, which_indices)
 
+    X = np.asarray(model.X)
+    
     # first, plot the output variance as a function of the latent space
-    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution)
-    Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1]))
+    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(X[:, [input_1, input_2]], resolution=resolution)
+    Xtest_full = np.zeros((Xtest.shape[0], X.shape[1]))
 
     def plot_function(x):
         Xtest_full[:, [input_1, input_2]] = x
@@ -48,7 +50,7 @@ def plot_latent(model, labels=None, which_indices=None,
         var = var[:, :1]
         return np.log(var)
     view = ImshowController(ax, plot_function,
-                            tuple(model.X.min(0)[:, [input_1, input_2]]) + tuple(model.X.max(0)[:, [input_1, input_2]]),
+                            tuple(X[:, [input_1, input_2]].min(0)) + tuple(X[:, [input_1, input_2]].max(0)),
                             resolution, aspect=aspect, interpolation='bilinear',
                             cmap=pb.cm.binary)
     
@@ -74,11 +76,11 @@ def plot_latent(model, labels=None, which_indices=None,
 
         index = np.nonzero(labels == ul)[0]
         if model.input_dim == 1:
-            x = model.X[index, input_1]
+            x = X[index, input_1]
             y = np.zeros(index.size)
         else:
-            x = model.X[index, input_1]
-            y = model.X[index, input_2]
+            x = X[index, input_1]
+            y = X[index, input_2]
         ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label)
 
     ax.set_xlabel('latent dimension %i' % input_1)
@@ -117,16 +119,17 @@ def plot_magnification(model, labels=None, which_indices=None,
         labels = np.ones(model.num_data)
 
     input_1, input_2 = most_significant_input_dimensions(model, which_indices)
-
+    X = np.asarray(model.X)
+    
     # first, plot the output variance as a function of the latent space
-    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution)
-    Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1]))
+    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(X[:, [input_1, input_2]], resolution=resolution)
+    Xtest_full = np.zeros((Xtest.shape[0], X.shape[1]))
     def plot_function(x):
         Xtest_full[:, [input_1, input_2]] = x
         mf=model.magnification(Xtest_full)
         return mf
     view = ImshowController(ax, plot_function,
-                            tuple(model.X.min(0)[:, [input_1, input_2]]) + tuple(model.X.max(0)[:, [input_1, input_2]]),
+                            tuple(X.min(0)[:, [input_1, input_2]]) + tuple(X.max(0)[:, [input_1, input_2]]),
                             resolution, aspect=aspect, interpolation='bilinear',
                             cmap=pb.cm.gray)
 
@@ -149,11 +152,11 @@ def plot_magnification(model, labels=None, which_indices=None,
 
         index = np.nonzero(labels == ul)[0]
         if model.input_dim == 1:
-            x = model.X[index, input_1]
+            x = X[index, input_1]
             y = np.zeros(index.size)
         else:
-            x = model.X[index, input_1]
-            y = model.X[index, input_2]
+            x = X[index, input_1]
+            y = X[index, input_2]
         ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label)
 
     ax.set_xlabel('latent dimension %i' % input_1)
diff --git a/GPy/util/visualize.py b/GPy/util/visualize.py
index 683c6c67..b13c100c 100644
--- a/GPy/util/visualize.py
+++ b/GPy/util/visualize.py
@@ -92,7 +92,7 @@ class lvm(matplotlib_show):
         :param latent_axes: the axes where the latent visualization should be plotted.
         """
         if vals == None:
-            vals = model.X[0]
+            vals = np.asarray(model.X[0])
 
         matplotlib_show.__init__(self, vals, axes=latent_axes)
 
@@ -171,21 +171,21 @@ class lvm_subplots(lvm):
     latent_axes is a np array of dimension np.ceil(input_dim/2),
     one for each pair of the latent dimensions.
     """
-    def __init__(self, vals, Model, data_visualize, latent_axes=None, sense_axes=None):
-        self.nplots = int(np.ceil(Model.input_dim/2.))+1
+    def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None):
+        self.nplots = int(np.ceil(model.input_dim/2.))+1
         assert len(latent_axes)==self.nplots
         if vals==None:
-            vals = Model.X[0, :]
+            vals = np.asarray(model.X[0, :])
         self.latent_values = vals 
 
         for i, axis in enumerate(latent_axes):
             if i == self.nplots-1:
-                if self.nplots*2!=Model.input_dim:
+                if self.nplots*2!=model.input_dim:
                     latent_index = [i*2, i*2]
-                lvm.__init__(self, self.latent_vals, Model, data_visualize, axis, sense_axes, latent_index=latent_index)
+                lvm.__init__(self, self.latent_vals, model, data_visualize, axis, sense_axes, latent_index=latent_index)
             else:
                 latent_index = [i*2, i*2+1]
-                lvm.__init__(self, self.latent_vals, Model, data_visualize, axis, latent_index=latent_index)
+                lvm.__init__(self, self.latent_vals, model, data_visualize, axis, latent_index=latent_index)