'Score each estimator on each task.
The number of tasks in X should match the number of tasks/estimators
given at fit time, i.e. we need
``X.shape[-1] == len(self.estimators_)``.
Parameters
X : array, shape (n_samples, nd_features, n_tasks)
The input samples. For each data slice, the corresponding estimator
scores the prediction, e.g.:
``[estimators[ii].score(X[..., ii], y) for ii in range(n_estimators)]``.
The feature dimension can be multidimensional e.g.
X.shape = (n_samples, n_features_1, n_features_2, n_tasks)
y : array, shape (n_samples,) | (n_samples, n_targets)
The target values.
Returns
score : array, shape (n_samples, n_estimators)
Score for each estimator/task.'

def score(self, X, y):

    from sklearn.metrics.scorer import check_scoring
    self._check_Xy(X)
    if (X.shape[(-1)] != len(self.estimators_)):
        raise ValueError('The    number    of    estimators    does    not    match    X.shape[-1]')
    scoring = check_scoring(self.base_estimator, self.scoring)
    y = _fix_auc(scoring, y)
    (parallel, p_func, n_jobs) = parallel_func(_sl_score, self.n_jobs)
    n_jobs = min(n_jobs, X.shape[(-1)])
    X_splits = np.array_split(X, n_jobs, axis=(-1))
    est_splits = np.array_split(self.estimators_, n_jobs)
    score = parallel((p_func(est, scoring, x, y) for (est, x) in zip(est_splits, X_splits)))
    score = np.concatenate(score, axis=0)
    return score
