Repository: scikit-learn/scikit-learn
File: sklearn/gaussian_process/tests/test_gpc.py (Python)

Code before:
            else:
                result.append(
                    [gpc.log_marginal_likelihood([theta]) for theta in kernel.theta]
                )
        if length_scales.ndim == 1:
            return np.stack(result)
        elif length_scales.ndim == 2:
            return np.stack(result).reshape(
                length_scales.shape[0], length_scales.shape[1]
            )

    lml_gradient = []
    length_scale_param_name = next(
        name for name in kernel.get_params() if name.endswith("length_scale")
    )
    for i, length_scale in enumerate(length_scales):
        kernel.set_params(**{length_scale_param_name: length_scale})
        lml_gradient.append(
            gpc.log_marginal_likelihood(kernel.theta, eval_gradient=True)[1][0]
        )

    lml_gradient_approx = derivative(
        evaluate_grad_at_length_scales,
        length_scales,
        maxiter=20,
    ).df

    assert_almost_equal(np.stack(lml_gradient), lml_gradient_approx, 3)


def test_random_starts(global_random_seed):
    # Test that an increasing number of random-starts of GP fitting only
    # increases the log marginal likelihood of the chosen theta.
    n_samples, n_features = 25, 2
    rng = np.random.RandomState(global_random_seed)
    X = rng.randn(n_samples, n_features) * 2 - 1
    y = (np.sin(X).sum(axis=1) + np.sin(3 * X).sum(axis=1)) > 0

    kernel = C(1.0, (1e-2, 1e2)) * RBF(
        length_scale=[1e-3] * n_features, length_scale_bounds=[(1e-4, 1e2)] * n_features
    )
    last_lml = -np.inf
    for n_restarts_optimizer in range(5):
        gp = GaussianProcessClassifier(
            kernel=kernel,
            n_restarts_optimizer=n_restarts_optimizer,
            random_state=global_random_seed,
        ).fit(X, y)
        lml = gp.log_marginal_likelihood(gp.kernel_.theta)
        assert lml > last_lml - np.finfo(np.float32).eps
        last_lml = lml

Review comment (suggestion):
```suggestion
    derivative_results = derivative(
         evaluate_grad_at_length_scales,
         length_scales,
         maxiter=20,
    )
    # Check that the numerical estimation converged within tolerances.
    assert np.all(derivative_results.status == 0)
    lml_gradient_approx = derivative_results.df
```

Code after:
    # Test that an increasing number of random-starts of GP fitting only
    # increases the log marginal likelihood of the chosen theta.
    n_samples, n_features = 25, 2
    rng = np.random.RandomState(global_random_seed)
    X = rng.randn(n_samples, n_features) * 2 - 1
    y = (np.sin(X).sum(axis=1) + np.sin(3 * X).sum(axis=1)) > 0

    kernel = C(1.0, (1e-2, 1e2)) * RBF(
        length_scale=[1e-3] * n_features, length_scale_bounds=[(1e-4, 1e2)] * n_features
    )
    last_lml = -np.inf
    for n_restarts_optimizer in range(5):
        gp = GaussianProcessClassifier(
            kernel=kernel,
            n_restarts_optimizer=n_restarts_optimizer,
            random_state=global_random_seed,
        ).fit(X, y)
        lml = gp.log_marginal_likelihood(gp.kernel_.theta)
        assert lml > last_lml - np.finfo(np.float32).eps
        last_lml = lml


@pytest.mark.parametrize("kernel", non_fixed_kernels)
def test_custom_optimizer(kernel, global_random_seed):
    # Test that GPC can use externally defined optimizers.
    # Define a dummy optimizer that simply tests 10 random hyperparameters
    def optimizer(obj_func, initial_theta, bounds):
        rng = np.random.RandomState(global_random_seed)
        theta_opt, func_min = (
            initial_theta,
            obj_func(initial_theta, eval_gradient=False),
        )
        for _ in range(10):
            theta = np.atleast_1d(
                rng.uniform(np.maximum(-2, bounds[:, 0]), np.minimum(1, bounds[:, 1]))
            )
            f = obj_func(theta, eval_gradient=False)
            if f < func_min:
                theta_opt, func_min = theta, f
        return theta_opt, func_min

    gpc = GaussianProcessClassifier(kernel=kernel, optimizer=optimizer)
    gpc.fit(X, y_mc)
    # Checks that optimizer improved marginal likelihood
    assert gpc.log_marginal_likelihood(
        gpc.kernel_.theta
    ) >= gpc.log_marginal_likelihood(kernel.theta)


@pytest.mark.parametrize("kernel", kernels)
def test_multi_class(kernel):