function [A, B, ell, d] = forward(x,xb, p,pb, theta, k0, c, sigma2_n, N, y, simPar, forceSingleKernel)
    % This function returns the GPR loss and the derivative of the loss w.r.t. each parameter. 

    sigma2_mu = 3000;
    sigma2_v = 3000;

    if simPar.channelEstimator.GPR.mixedKernel == false || forceSingleKernel
        [K, grads] = get_K(x,x, p,p, theta, k0, c);         % Compute the kernel K(x, p, theta). 
        Ky = K + sigma2_n*eye(simPar.N*simPar.T);                       % Compute the noisy kernel Ky. 
        %         Ky_inv_y = Ky\y;                                
        % R = chol(Ky);  
                [Ky_inv, logdet] = complex_cholesky_ops(Ky);
        Ky_inv_y = Ky_inv*y;

%         Ky_inv_y = Ky\y;  B = -real(log(det(Ky))); 

A = real(-y'*Ky_inv_y);   B = -real(logdet); 
%         A = real(-y'*Ky_inv_y);   B = -real(log(det(Ky)));
%         ell = A+B; 
        ell = A+B;                                          % The GPR target function ell = log p(y|x, theta). (29)
    
    
        d_ell_K = conj(Ky_inv_y)*(Ky_inv_y.') - conj(inv(Ky));    % The derivative of ell w.r.t. K.  (aaH − K−1 y )(31)
        
        gvs     = [grads.v];
        gmus    = [grads.mu];
        % grhos   = [grads.rho];
        
        d_ell_mu    = 2*real(gmus*d_ell_K(:));              % Compute the derivative of ell w.r.t. mu. (32)
        d_ell_v     = 2*real(gvs*d_ell_K(:));
        % d_ell_rho   = 2*real(grhos*d_ell_K(:));
        
    
        % Normalization on mu. 
        mu = theta.mu;
        v = theta.v;
        
        
        % sigma2_mu = 3e9; 
        
        % Derivatives. 
        d = struct("d_mu", d_ell_mu - 2*mu/sigma2_mu,"d_v", d_ell_v - 2*v/sigma2_v); 
    
        % Loss function. 
        ell = ell - norm(mu)^2/sigma2_mu- norm(v)^2/sigma2_v; %?

    else
        N_kernels = simPar.channelEstimator.GPR.mixKernelNum;
        Ks = cell(N_kernels, 1);
        grads = cell(N_kernels, 1);
        w = theta.kernelWeights;

        Ky = sigma2_n * eye(N); 
        for idx = 1:N_kernels
            [Ks{idx}, grads{idx}] = get_K(x,x, p,p, theta.kernelParams(idx), k0, c); 
            Ky = Ky + w(idx)*Ks{idx};     % mixed kernel. 
        end

        %         Ky_inv_y = Ky\y;                                
        % R = chol(Ky);  
                [Ky_inv, logdet] = complex_cholesky_ops(Ky);
        Ky_inv_y = Ky_inv*y;

%         Ky_inv_y = Ky\y;  B = -real(log(det(Ky))); 

A = real(-y'*Ky_inv_y);   B = -real(logdet); 
%         A = real(-y'*Ky_inv_y);   B = -real(log(det(Ky)));
        ell = A+B; 

        d_ell_K = conj(Ky_inv_y)*(Ky_inv_y.') - conj(inv(Ky)); 

        
        d_mu = zeros(3, N_kernels);
        d_w = zeros(1, N_kernels);

        for idx = 1:N_kernels
            gmus    = [grads{idx}.mu];
%             d_mu(:, idx)    = 2*w(idx)*real(gmus*d_ell_K(:)) - 2*theta.kernelParams(idx).mu/sigma2_mu; 

            K_component = Ks{idx}; 
            d_w(idx) = 2*real((K_component(:)).'*d_ell_K(:)); %和trace效果相同

            ell = ell - norm(theta.kernelParams(idx).mu)^2/sigma2_mu; 
        end

        d = struct("d_mu", d_mu, "d_w", d_w); 
    end

end

