def observe(self, obs_mesh, obs_V, output_type='r'):
        """
        Observes self on obs_mesh with observation variance obs_V.
        Output_type controls the information returned:

        'r' : returns information needed by Realization objects.
        'o' : returns information needed by function observe.
        's' : returns information needed by the Gaussian process
              submodel.
        """

        # print 'C.observe called'

        # Number of spatial dimensions.
        ndim = obs_mesh.shape[1]

        if self.ndim is not None:
            if not ndim == self.ndim:
                raise ValueError(
                    "Dimension of observation mesh is not equal to dimension of base mesh.")
        else:
            self.ndim = ndim

        # print ndim

        # =====================================
        # = If self hasn't been observed yet: =
        # =====================================
        if not self.observed:

            # If self has not been observed, get the Cholesky factor of self(obs_mesh, obs_mesh)
            # and the side information and store it.

            # Rank so far is 0.
            m_old = 0

            # Number of observation points so far is 0.
            N_old = 0

            if output_type != 's':
                obs_dict = self.cholesky(
                    obs_mesh,
                    apply_pivot=False,
                    nugget=obs_V,
                    regularize=False,
                    rank_limit=self.rank_limit)
            else:
                C_eval = self.__call__(obs_mesh, obs_mesh, regularize=False)
                U = C_eval.copy('F')
                for i in xrange(U.shape[0]):
                    U[i, i] += obs_V[i]
                info = dpotrf_wrap(U)
                if info > 0:
                    raise LinAlgError(
                        "Matrix does not appear to be positive definite by row %i. Could not observe with assume_full_rank=True." %
                        info)
                obs_dict = {
                    'U': U,
                    'pivots': arange(
                        U.shape[0]),
                    'U_new': U,
                    'C_eval': C_eval}
            obs_dict_new = obs_dict

            # Rank of self(obs_mesh, obs_mesh)
            m_new = obs_dict['U'].shape[0]

            # Upper-triangular Cholesky factor of self(obs_mesh, obs_mesh)
            self.full_Uo = obs_dict['U']
            # print (self.full_Uo[:,argsort(obs_dict['pivots'])].T*self.full_Uo[:,argsort(obs_dict['pivots'])] - self(obs_mesh,obs_mesh)).max()

            # Upper-triangular square Cholesky factor of self(obs_mesh_*,
            # obs_mesh_*). See documentation.
            self.Uo = obs_dict['U'][:, :m_new]

            # Pivots.
            piv_new = obs_dict['pivots']
            self.full_piv = piv_new
            self.obs_piv = piv_new[:m_new]

            # Remember full observation mesh.
            self.full_obs_mesh = obs_mesh

            # relevant slice is the positive-definite indices, which get into
            # obs_mesh_*. See documentation.
            relevant_slice = self.obs_piv

            # obs_mesh_new is obs_mesh_* from documentation.
            obs_mesh_new = obs_mesh[relevant_slice, :]

            self.obs_mesh = obs_mesh_new
            self.obs_V = obs_V[piv_new]
            self.obs_len = m_new

        # =======================================
        # = If self has been observed already:  =
        # =======================================
        else:

            # If self has been observed, get the Cholesky factor of the _full_ observation mesh (new
            # and old observations) using continue_cholesky, along with side
            # information, and store it.

            # Extract information from self's existing attributes related to
            # the observation mesh..
            obs_old, piv_old = self.Uo, self.obs_piv

            # Rank of self's evaluation on the observation mesh so far.
            m_old = len(self.obs_piv)

            # Number of observations so far.
            N_old = self.full_obs_mesh.shape[0]

            # Number of new observations.
            N_new = obs_mesh.shape[0]

            # Call to self.continue_cholesky.
            obs_dict_new = self.continue_cholesky(x=obs_mesh,
                                                  x_old=self.full_obs_mesh,
                                                  chol_dict_old={
                                                      'U': self.full_Uo,
                                                      'pivots': self.full_piv},
                                                  apply_pivot=False,
                                                  observed=False,
                                                  regularize=False,
                                                  nugget=obs_V,
                                                  assume_full_rank=output_type == 's',
                                                  rank_limit=self.rank_limit)

            if output_type == 's':
                C_eval = obs_dict_new['C_eval']

            # Full Cholesky factor of self(obs_mesh, obs_mesh), where obs_mesh
            # is the combined observation mesh.
            self.full_Uo = obs_dict_new['U']

            # Rank of self(obs_mesh, obs_mesh)
            m_new = self.full_Uo.shape[0]

            # Square upper-triangular Cholesky factor of self(obs_mesh_*,
            # obs_mesh_*). See documentation.
            self.Uo = self.full_Uo[:, :m_new]

            # Pivots.
            piv_new = obs_dict_new['pivots']
            self.obs_piv = piv_new[:m_new]
            self.full_piv = piv_new

            # Concatenate old and new observation meshes.
            self.full_obs_mesh = vstack((self.full_obs_mesh, obs_mesh))
            relevant_slice = piv_new[m_old:m_new] - N_old
            obs_mesh_new = obs_mesh[relevant_slice, :]

            # Remember obs_mesh_* and corresponding observation variances.
            self.obs_mesh = vstack((self.obs_mesh, obs_mesh[relevant_slice, :]))
            self.obs_V = hstack((self.obs_V, obs_V[relevant_slice]))

            # Length of obs_mesh_*.
            self.obs_len = m_new

        self.observed = True
        # Output expected by Realization
        if output_type == 'r':
            return relevant_slice, obs_mesh_new, self.full_Uo[
                m_old:, argsort(piv_new)[N_old:]], self.full_Uo[:m_old, argsort(piv_new)[N_old:]]

        # Ouptut expected by observe
        if output_type == 'o':
            return relevant_slice, obs_mesh_new

        # Output expected by the GP submodel
        if output_type == 's':
            return obs_dict_new['U_new'], obs_dict_new[
                'C_eval'], self.full_Uo[:m_old, argsort(piv_new)[N_old:]]
Observes self on obs_mesh with observation variance obs_V.
        Output_type controls the information returned:

        'r' : returns information needed by Realization objects.
        'o' : returns information needed by function observe.
        's' : returns information needed by the Gaussian process
              submodel.