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

Code before:
    def evaluate_grad_at_length_scales(length_scales):
        length_scale_param_name = next(
            name for name in kernel.get_params() if name.endswith("length_scale")
        )
        result = []
        for i, length_scale in enumerate(length_scales.flatten()):
            kernel.set_params(**{length_scale_param_name: length_scale})
            if type(kernel) == Product or len(kernel.theta) == 1:
                result.append(gpc.log_marginal_likelihood(kernel.theta))

            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
    )

Review comment (suggestion):
```suggestion
    for length_scale in length_scales:
```

Code after:
    epsilon = 1e-9
    lml_gradient_approx = approx_fprime(
        kernel.theta.copy(),
        lambda theta: gpr.log_marginal_likelihood(theta, False),
        epsilon=epsilon,
    )
    assert_allclose(lml_gradient, lml_gradient_approx, rtol=1e-4, atol=epsilon * 100)


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


@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