def _fit_RSA_UV(self, X, Y, X_base,
                    scan_onsets=None, coords=None, inten=None):
        """ The major utility of fitting Bayesian RSA.
            Note that there is a naming change of variable. X in fit()
            is changed to Y here, and design in fit() is changed to X here.
            This is because we follow the tradition that X expresses the
            variable defined (controlled) by the experimenter, i.e., the
            time course of experimental conditions convolved by an HRF,
            and Y expresses data.
            However, in wrapper function fit(), we follow the naming
            routine of scikit-learn.
        """
        GP_inten = self.GP_inten
        GP_space = self.GP_space
        rank = self.rank
        n_V = np.size(Y, axis=1)
        n_T = np.size(Y, axis=0)
        n_C = np.size(X, axis=1)
        l_idx, rank = self._chol_idx(n_C, rank)
        n_l = np.size(l_idx[0])  # the number of parameters for L

        t_start = time.time()

        D, F, run_TRs, n_run = self._prepare_DF(
            n_T, scan_onsets=scan_onsets)
        XTY, XTDY, XTFY, YTY_diag, YTDY_diag, YTFY_diag, XTX, \
            XTDX, XTFX = self._prepare_data_XY(X, Y, D, F)

        X0TX0, X0TDX0, X0TFX0, XTX0, XTDX0, XTFX0, \
            X0TY, X0TDY, X0TFY, X0, X_base, n_X0, idx_DC = \
            self._prepare_data_XYX0(
                X, Y, X_base, None, D, F, run_TRs, no_DC=False)
        # Prepare the data for fitting. These pre-calculated matrices
        # will be re-used a lot in evaluating likelihood function and
        # gradient.
        # DC component will be added to the nuisance regressors.
        # In later steps, we do not need to add DC components again

        dist2, inten_diff2, space_smooth_range, inten_smooth_range,\
            n_smooth = self._calc_dist2_GP(
                coords=coords, inten=inten,
                GP_space=GP_space, GP_inten=GP_inten)
        # Calculating the distance between voxel locations and betweeen
        # voxel intensities. These are used if a Gaussian Process prior
        # is requested to regularize log(SNR^2)

        idx_param_sing, idx_param_fitU, idx_param_fitV = \
            self._build_index_param(n_l, n_V, n_smooth)
        # Indexes to find each parameter in a combined parameter vector.

        current_GP = np.zeros(n_smooth)

        # We will perform the fitting in 2~3 steps:
        # (1) A preliminary fitting assuming all voxels share
        # exactly the same temporal covariance matrix for their noise.
        # SNR is assumed to be 1 for all voxels in this fitting.
        # Therefore, there are only n_l+2 free parameters.
        # (2) (optional) A fitting which allows each voxel to have their
        # own pseudo-SNR and AR(1) coefficients. But no Gaussian Process
        # prior is imposed on log(SNR). This step is neglected if GP
        # prior is not requested. This step allows the SNR parameters to
        # move closer to their correct values before GP is introduced.
        # This step alternately fits the shared covariance and voxel-
        # specific variance. It fits for init_iter steps and the
        # tolerance is also increased by a factor of 5 to speed up
        # fitting.
        # (3) Final fitting. If GP prior is requested, it will be
        # introduced in this step. Otherwise, just fit as the previous
        # step, but using un-altered tolerance setting, and n_iter
        # as the number of iteration.

        # Step 1 fitting, with a simplified model
        current_vec_U_chlsk_l, current_a1, current_logSigma2 = \
            self._initial_fit_singpara(
                XTX, XTDX, XTFX, YTY_diag, YTDY_diag, YTFY_diag,
                XTY, XTDY, XTFY, X0TX0, X0TDX0, X0TFX0,
                XTX0, XTDX0, XTFX0, X0TY, X0TDY, X0TFY,
                X, Y, X0, idx_param_sing,
                l_idx, n_C, n_T, n_V, n_l, n_run, n_X0, rank)

        current_logSNR2 = -current_logSigma2
        norm_factor = np.mean(current_logSNR2)
        current_logSNR2 = current_logSNR2 - norm_factor
        X_res = None

        # Step 2 fitting, which only happens if
        # GP prior is requested
        if GP_space:
            current_vec_U_chlsk_l, current_a1, current_logSNR2, X_res\
                = self._fit_diagV_noGP(
                    XTY, XTDY, XTFY, YTY_diag, YTDY_diag, YTFY_diag,
                    XTX, XTDX, XTFX, X, Y, X_base, X_res, D, F, run_TRs,
                    current_vec_U_chlsk_l,
                    current_a1, current_logSNR2,
                    idx_param_fitU, idx_param_fitV,
                    l_idx, n_C, n_T, n_V, n_l, n_run, n_X0, rank)

            current_GP[0] = np.log(np.min(
                dist2[np.tril_indices_from(dist2, k=-1)]))
            # We start fitting the model with GP prior with a small
            # length scale: the size of voxels.
            # Alternatively, initialize with a large distance.
            # Further testing of initial parameters need to be done.

#             current_GP[0] = np.log(np.max(dist2)/4.0)

            logger.debug('current GP[0]:{}'.format(current_GP[0]))
            if GP_inten:
                current_GP[1] = np.log(np.maximum(
                    np.percentile(inten_diff2[np.tril_indices_from(
                        inten_diff2, k=-1)], 2), 0.5))
                logger.debug(
                    'current GP[1]:{}'.format(current_GP[1]))
                # We start the length scale for intensity with
                # a small value. A heuristic is 2 percentile of
                # all the square differences. But it should not be
                # smaller than 0.5. This limit is set in case
                # many voxels have close to equal intensities,
                # which might render 2 percentile to 0.

        # Step 3 fitting. GP prior is imposed if requested.
        # In this step, unless auto_nuisance is set to False, X_res
        # will be re-estimated from the residuals after each step
        # of fitting. And X0 will be concatenation of X_base and X_res
        logger.debug('indexing:{}'.format(idx_param_fitV))
        logger.debug('initial GP parameters:{}'.format(current_GP))
        current_vec_U_chlsk_l, current_a1, current_logSNR2,\
            current_GP, X_res = self._fit_diagV_GP(
                XTY, XTDY, XTFY, YTY_diag, YTDY_diag, YTFY_diag,
                XTX, XTDX, XTFX, X, Y, X_base, X_res, D, F, run_TRs,
                current_vec_U_chlsk_l,
                current_a1, current_logSNR2, current_GP, n_smooth,
                idx_param_fitU, idx_param_fitV,
                l_idx, n_C, n_T, n_V, n_l, n_run, n_X0, rank,
                GP_space, GP_inten, dist2, inten_diff2,
                space_smooth_range, inten_smooth_range)

        estU_chlsk_l_AR1_UV = np.zeros([n_C, rank])
        estU_chlsk_l_AR1_UV[l_idx] = current_vec_U_chlsk_l

        est_cov_AR1_UV = np.dot(estU_chlsk_l_AR1_UV, estU_chlsk_l_AR1_UV.T)

        est_rho1_AR1_UV = 2 / np.pi * np.arctan(current_a1)
        est_SNR_AR1_UV = np.exp(current_logSNR2 / 2.0)

        # Calculating est_sigma_AR1_UV, est_sigma_AR1_UV,
        # est_beta_AR1_UV and est_beta0_AR1_UV
        X0TX0, X0TDX0, X0TFX0, XTX0, XTDX0, XTFX0, \
            X0TY, X0TDY, X0TFY, X0, X_base, n_X0, _ \
            = self._prepare_data_XYX0(
                X, Y, X_base, X_res, D, F, run_TRs, no_DC=True)

        X0TAX0, XTAX0, X0TAY, X0TAX0_i, \
            XTAcorrX, XTAcorrY, YTAcorrY, LTXTAcorrY, XTAcorrXL, LTXTAcorrXL\
            = self._precompute_ar1_quad_forms(XTY, XTDY, XTFY,
                                              YTY_diag, YTDY_diag, YTFY_diag,
                                              XTX, XTDX, XTFX, X0TX0, X0TDX0,
                                              X0TFX0, XTX0, XTDX0, XTFX0, X0TY,
                                              X0TDY, X0TFY,
                                              estU_chlsk_l_AR1_UV,
                                              est_rho1_AR1_UV, n_V, n_X0)
        LL, LAMBDA_i, LAMBDA, YTAcorrXL_LAMBDA, sigma2 \
            = self._calc_LL(est_rho1_AR1_UV, LTXTAcorrXL, LTXTAcorrY, YTAcorrY,
                            X0TAX0, est_SNR_AR1_UV**2,
                            n_V, n_T, n_run, rank, n_X0)
        est_sigma_AR1_UV = sigma2**0.5
        est_beta_AR1_UV = est_SNR_AR1_UV**2 \
            * np.dot(estU_chlsk_l_AR1_UV, YTAcorrXL_LAMBDA.T)
        est_beta_AR1_UV_latent =  \
            est_SNR_AR1_UV**2 * YTAcorrXL_LAMBDA.T
        # the latent term means that X*L multiplied by this term
        # is the same as X*beta. This will be used for decoding
        # and cross-validating, in case L is low-rank
        est_beta0_AR1_UV = np.einsum(
            'ijk,ki->ji', X0TAX0_i,
            (X0TAY - np.einsum('ikj,ki->ji', XTAX0, est_beta_AR1_UV)))

        # Now we want to collapse all beta0 corresponding to DC components
        # of different runs to a single map, and preserve only one DC component
        # across runs. This is because they should express the same component
        # and the new data to transform do not necessarily have the same
        # numbers of runs as the training data.

        if idx_DC.size > 1:
            collapsed_DC = np.sum(X0[:, idx_DC], axis=1)
            X0 = np.insert(np.delete(X0, idx_DC, axis=1), 0,
                           collapsed_DC, axis=1)
            collapsed_beta0 = np.mean(est_beta0_AR1_UV[idx_DC, :], axis=0)
            est_beta0_AR1_UV = np.insert(
                np.delete(est_beta0_AR1_UV, idx_DC, axis=0),
                0, collapsed_beta0, axis=0)
        t_finish = time.time()
        logger.info(
            'total time of fitting: {} seconds'.format(t_finish - t_start))
        logger.debug('final GP parameters:{}'.format(current_GP))
        if GP_space:
            est_space_smooth_r = np.exp(current_GP[0] / 2.0)
            if GP_inten:
                est_intensity_kernel_r = np.exp(current_GP[1] / 2.0)
                K_major = np.exp(- (dist2 / est_space_smooth_r**2 +
                                    inten_diff2 / est_intensity_kernel_r**2)
                                 / 2.0)
            else:
                est_intensity_kernel_r = None
                K_major = np.exp(- dist2 / est_space_smooth_r**2 / 2.0)
            K = K_major + np.diag(np.ones(n_V) * self.eta)
            invK_tilde_log_SNR = np.linalg.solve(K, current_logSNR2) / 2
            log_SNR_invK_tilde_log_SNR = np.dot(current_logSNR2,
                                                invK_tilde_log_SNR) / 2
            tau2, _ = self.tau2_prior(log_SNR_invK_tilde_log_SNR, n_V,
                                      self.tau_range)
            est_std_log_SNR = tau2 ** 0.5
        else:
            est_space_smooth_r = None
            est_intensity_kernel_r = None
            est_std_log_SNR = None
        return est_cov_AR1_UV, estU_chlsk_l_AR1_UV, est_SNR_AR1_UV, \
            est_beta_AR1_UV, est_beta0_AR1_UV, est_beta_AR1_UV_latent,\
            est_sigma_AR1_UV, est_rho1_AR1_UV, est_space_smooth_r, \
            est_std_log_SNR, est_intensity_kernel_r, X0
The major utility of fitting Bayesian RSA.
            Note that there is a naming change of variable. X in fit()
            is changed to Y here, and design in fit() is changed to X here.
            This is because we follow the tradition that X expresses the
            variable defined (controlled) by the experimenter, i.e., the
            time course of experimental conditions convolved by an HRF,
            and Y expresses data.
            However, in wrapper function fit(), we follow the naming
            routine of scikit-learn.