Overwrite checkgrad method to check whole block instead of looping through

        Shows diagnostics using matshow instead

        :param verbose: If True, print a "full" checking of each parameter
        :type verbose: bool
        :param step: The size of the step around which to linearise the objective
        :type step: float (default 1e-6)
        :param tolerance: the tolerance allowed (see note)
        :type tolerance: float (default 1e-3)

        Note:-
           The gradient is considered correct if the ratio of the analytical
           and numerical gradients is within <tolerance> of unity.

def checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3, block_indices=None, plot=False):
        """
        Overwrite checkgrad method to check whole block instead of looping through

        Shows diagnostics using matshow instead

        :param verbose: If True, print a "full" checking of each parameter
        :type verbose: bool
        :param step: The size of the step around which to linearise the objective
        :type step: float (default 1e-6)
        :param tolerance: the tolerance allowed (see note)
        :type tolerance: float (default 1e-3)

        Note:-
           The gradient is considered correct if the ratio of the analytical
           and numerical gradients is within <tolerance> of unity.
        """
        try:
            import numdifftools as nd
        except:
            raise ImportError("Don't have numdifftools package installed, it is not a GPy dependency as of yet, it is only used for hessian tests")

        if target_param:
            raise NotImplementedError('Only basic functionality is provided with this gradchecker')

        #Repeat for each parameter, not the nicest but shouldn't be many cases where there are many
        #variables
        current_index = 0
        for name, shape in zip(self.names, self.shapes):
            current_size = numpy.prod(shape)
            x = self.optimizer_array.copy()
            #x = self._get_params_transformed().copy()
            x = x[current_index:current_index + current_size].reshape(shape)

            # Check gradients
            analytic_hess = self._ddf(x)
            if analytic_hess.shape[1] == 1:
                analytic_hess = numpy.diagflat(analytic_hess)

            #From the docs:
            #x0 : vector location
            #at which to differentiate fun
            #If x0 is an N x M array, then fun is assumed to be a function
            #of N*M variables., thus we must have it flat, not (N,1), but just (N,)
            #numeric_hess_partial = nd.Hessian(self._f, vectorized=False)
            numeric_hess_partial = nd.Jacobian(self._df, vectorized=False)
            #numeric_hess_partial = nd.Derivative(self._df, vectorized=True)
            numeric_hess = numeric_hess_partial(x)

            check_passed = self.checkgrad_block(analytic_hess, numeric_hess, verbose=verbose, step=step, tolerance=tolerance, block_indices=block_indices, plot=plot)
            current_index += current_size
        return check_passed