code: def _bsa_dpmm(
      hrois,
      prior_h0,
      subjects,
      coords,
      sigma,
      prevalence_pval,
      prevalence_threshold,
      dof=10,
      alpha=0.5,
      n_iter=1000,
      burnin=100,
      algorithm="density",
  ):
      """Estimation of the population level model of activation density using
      dpmm and inference
  
      Parameters
      ----------
      hrois: list of nipy.labs.spatial_models.hroi.HierarchicalROI instances
             representing individual ROIs
             Let nr be the number of terminal regions across subjects
      prior_h0: array of shape (nr)
                mixture-based prior probability
                that the terminal regions are true positives
      subjects: array of shape (nr)
                subject index associated with the terminal regions
      coords: array of shape (nr, coord.shape[1])
              coordinates of the of the terminal regions
      sigma: float > 0,
             expected cluster scatter in the common space in units of coord
      prevalence_pval: float in the [0,1] interval, optional
                       p-value of the prevalence test
      prevalence_threshold: float in the rannge [0,nsubj]
                           null hypothesis on region prevalence
      dof: float > 0, optional,
           degrees of freedom of the prior
      alpha: float > 0, optional,
             creation parameter of the DPMM
      niter: int, optional,
             number of iterations of the DPMM
      burnin: int, optional,
              number of iterations of the DPMM
      algorithm: {'density', 'co_occurrence'}, optional,
                 algorithm used in the DPMM inference
  
      Returns
      -------
      landmarks: instance of sbf.LandmarkRegions
                 that describes the ROIs found in inter-subject inference
                 If no such thing can be defined landmarks is set to None
      hrois: List of nipy.labs.spatial_models.hroi.HierarchicalROI instances
             representing individual ROIs
      """
      from nipy.algorithms.graph.field import (
          field_from_coo_matrix_and_data,
      )
  
      domain = hrois[0].domain
      n_subjects = len(hrois)
  
      landmarks = None
      density = np.zeros(domain.size)
      if len(subjects) < 1:
          return landmarks, hrois
  
      null_density = 1.0 / domain.local_volume.sum()
  
      # prepare the DPMM
      dim = domain.em_dim
      prior_precision = 1.0 / (sigma**2) * np.ones((1, dim))
  
      # n_iter = number of iterations to estimate density
      if algorithm == "density":
          density, post_proba = _dpmm(
              coords,
              alpha,
              null_density,
              dof,
              prior_precision,
              prior_h0,
              subjects,
              domain.coord,
              n_iter=n_iter,
              burnin=burnin,
          )
          # associate labels with coords
          Fbeta = field_from_coo_matrix_and_data(
              domain.topology, density
          )
          _, label = Fbeta.custom_watershed(0, null_density)
          midx = np.array(
              [
                  np.argmin(np.sum((domain.coord - coord_) ** 2, 1))
                  for coord_ in coords
              ]
          )
          components = label[midx]
      elif algorithm == "co-occurrence":
          post_proba, density, co_clustering = _dpmm(
              coords,
              alpha,
              null_density,
              dof,
              prior_precision,
              prior_h0,
              subjects,
              n_iter=n_iter,
              burnin=burnin,
              co_clust=True,
          )
          contingency_graph = wgraph_from_coo_matrix(co_clustering)
          if contingency_graph.E > 0:
              contingency_graph.remove_edges(
                  contingency_graph.weights > 0.5
              )
  
          components = contingency_graph.cc()
          components[density < null_density] = (
              components.max()
              + 1
              + np.arange(np.sum(density < null_density))
          )
      else:
          raise ValueError("Unknown algorithm")
  
      # append some information to the hroi in each subject
      for subject in range(n_subjects):
          bfs = hrois[subject]
          if bfs is None:
              continue
          if bfs.k == 0:
              bfs.set_roi_feature("label", np.array([]))
              continue
  
          leaves_pos = [bfs.select_id(k) for k in bfs.get_leaves_id()]
          # save posterior proba
          post_proba_ = np.zeros(bfs.k)
          post_proba_[leaves_pos] = post_proba[subjects == subject]
          bfs.set_roi_feature("posterior_proba", post_proba_)
  
          # save prior proba
          prior_proba = np.zeros(bfs.k)
          prior_proba[leaves_pos] = 1 - prior_h0[subjects == subject]
          bfs.set_roi_feature("prior_proba", prior_proba)
  
          # assign labels to ROIs
          roi_label = -np.ones(bfs.k).astype(np.int_)
          roi_label[leaves_pos] = components[subjects == subject]
          # when parent regions has similarly labelled children,
          # include it also
          roi_label = bfs.make_forest().propagate_upward(roi_label)
          bfs.set_roi_feature("label", roi_label)
  
      # derive the group-level landmarks
      # with a threshold on the number of subjects
      # that are represented in each one
      landmarks, new_labels = build_landmarks(
          domain,
          coords,
          subjects,
          np.array(components),
          1 - prior_h0,
          prevalence_pval,
          prevalence_threshold,
          sigma,
      )
  
      # relabel the regions
      _update_hroi_labels(hrois, new_labels)
  
      return landmarks, hrois
  
file_path: nipy/nipy/labs/spatial_models/bayesian_structural_analysis.py
imports: import numpy
  import scipy.stats
  from nipy.algorithms.graph import wgraph_from_coo_matrix
  from ...algorithms.statistics.empirical_pvalue import NormalEmpiricalNull, gamma_gaussian_fit, three_classes_GMM_fit
  from .hroi import HROI_as_discrete_domain_blobs
  from .structural_bfls import build_landmarks
label_desc: Use image node_id 5 for calling a global function with example usage: _bsa_dpmm(hrois, prior_h0, subjects, coords, sigma, prevalence_pval, prevalence_threshold, dof, alpha, n_iter, burnin, algorithm) and returns: landmarks,  hrois, landmarks,  hrois