   for (auto i = 0; i < trees.tree_limit; i++) {
     TreeEnsemble tree;
     trees.get_tree(tree, i);
-    std::vector<gpu_treeshap::PathElement> tmp_path;
-    tmp_path.reserve(tree.max_depth);
-    tmp_path.emplace_back(0, -1, 0, -inf, inf, false, 1.0, 0.0f);
-    RecurseTree(0, tree, &tmp_path, &paths, &path_idx);
+    for (auto j = 0ull; j < tree.num_outputs; j++) {
+      std::vector<gpu_treeshap::PathElement> tmp_path;
+      tmp_path.reserve(tree.max_depth);
+      tmp_path.emplace_back(0, -1, 0, -inf, inf, false, 1.0, 0.0f);
+      RecurseTree(0, tree, &tmp_path, &paths, &path_idx, j, tree.num_outputs);
+    }
   }
   return paths;
 }

  DeviceExplanationDataset::DenseDatasetWrapper X =
      device_data.GetDeviceAccessor();

  thrust::device_vector<float> phis((X.NumCols() + 1) * X.NumRows());
  gpu_treeshap::GPUTreeShap(X, paths.begin(), paths.end(), 1, phis.data().get(),
                            phis.size());
  // Add the base offset term to bias
  auto base_offset = trees.base_offset[0];
  auto counting = thrust::make_counting_iterator(size_t(0));
  auto d_phis = phis.data().get();
  thrust::for_each(counting, counting + X.NumRows(),
                   [=] __device__(size_t row_idx) {
                     auto phi_idx = gpu_treeshap::IndexPhi(
                         row_idx, 1, 0, X.NumCols(), X.NumCols());
                     d_phis[phi_idx] += base_offset;
                   });
  thrust::copy(phis.begin(), phis.end(), out_contribs);
}

inline void dense_tree_path_dependent_interactions_gpu(

    int num_rows;
    int num_cols;

  public:
    DenseDatasetWrapper() = default;
    DenseDatasetWrapper(const float *data, const bool *missing, int num_rows,
                        int num_cols)
        : data(data), missing(missing), num_rows(num_rows), num_cols(num_cols) {
    }
    __device__ float GetElement(size_t row_idx, size_t col_idx) const {
      auto idx = row_idx * num_cols + col_idx;
      if (missing[idx]) {
