code: def _optim_hparcel(
      feature,
      domain,
      graphs,
      nb_parcel,
      lamb=1.0,
      dmax=10.0,
      niter=5,
      initial_mask=None,
      chunksize=1.0e5,
      verbose=0,
  ):
      """Core function of the heirrachical parcellation procedure.
  
      Parameters
      ----------
      feature: list of subject-related feature arrays
      Pa : parcellation instance that is updated
      graphs: graph that represents the topology of the parcellation
      anat_coord: array of shape (nvox,3) space defining set of coordinates
      nb_parcel: int
                 the number of desired parcels
      lamb=1.0: parameter to weight position
                and feature impact on the algorithm
      dmax = 10: locality parameter (in the space of anat_coord)
                to limit surch volume (CPU save)
      chunksize = int, optional
      niter = 5: number of iterations in the algorithm
      verbose=0: verbosity level
  
      Returns
      -------
      U: list of arrays of length nsubj
         subject-dependent parcellations
      Proto_anat: array of shape (nvox) labelling of the common space
                  (template parcellation)
      """
      nb_subj = len(feature)
  
      # a1. perform a rough clustering of the data to make prototype
      indiv_coord = np.array(
          [
              domain.coord[initial_mask[:, s] > -1]
              for s in range(nb_subj)
          ]
      )
      reduced_anat, reduced_feature = _reduce_and_concatenate(
          indiv_coord, feature, chunksize
      )
  
      _, labs, _ = kmeans(
          reduced_feature, nb_parcel, Labels=None, maxiter=10
      )
      proto_anat = [
          np.mean(reduced_anat[labs == k], 0) for k in range(nb_parcel)
      ]
      proto_anat = np.array(proto_anat)
      proto = [
          np.mean(reduced_feature[labs == k], 0)
          for k in range(nb_parcel)
      ]
      proto = np.array(proto)
  
      # a2. topological model of the parcellation
      # group-level part
      spatial_proto = Field(nb_parcel)
      spatial_proto.set_field(proto_anat)
      spatial_proto.voronoi_diagram(proto_anat, domain.coord)
      spatial_proto.set_gaussian(proto_anat)
      spatial_proto.normalize()
  
      for git in range(niter):
          LP = []
          LPA = []
          U = []
          Energy = 0
          for s in range(nb_subj):
              # b.subject-specific instances of the model
              # b.0 subject-specific information
              Fs = feature[s]
              lac = indiv_coord[s]
              target = proto_anat.copy()
              lseeds = np.zeros(nb_parcel, np.int_)
              aux = np.argsort(rand(nb_parcel))
              toto = np.zeros(lac.shape[0])
              for j in range(nb_parcel):
                  # b.1 speed-up :only take a small ball
                  i = aux[j]
                  dx = lac - target[i]
                  iz = np.nonzero(np.sum(dx**2, 1) < dmax**2)
                  iz = np.reshape(iz, np.size(iz))
                  if np.size(iz) == 0:
                      iz = np.array([np.argmin(np.sum(dx**2, 1))])
  
                  # b.2: anatomical constraints
                  lanat = np.reshape(
                      lac[iz], (np.size(iz), domain.coord.shape[1])
                  )
                  pot = np.zeros(np.size(iz))
                  JM, rmin = _exclusion_map(
                      i, spatial_proto, target, lanat
                  )
                  pot[JM < 0] = np.inf
                  pot[JM >= 0] = -JM[JM >= 0]
  
                  # b.3: add feature discrepancy
                  df = Fs[iz] - proto[i]
                  df = np.reshape(df, (np.size(iz), proto.shape[1]))
                  pot += lamb * np.sum(df**2, 1)
  
                  # b.4: solution
                  if np.sum(np.isinf(pot)) == np.size(pot):
                      pot = np.sum(dx[iz] ** 2, 1)
  
                  sol = iz[np.argmin(pot)]
                  target[i] = lac[sol]
                  lseeds[i] = sol
                  toto[sol] = 1
  
              if verbose > 1:
                  jm = _field_gradient_jac(spatial_proto, target)
                  print(jm.min(), jm.max(), np.sum(toto > 0))
  
              # c.subject-specific parcellation
              g = graphs[s]
              f = Field(g.V, g.edges, g.weights, Fs)
              U.append(f.constrained_voronoi(lseeds))
  
              Energy += np.sum((Fs - proto[U[-1]]) ** 2) / np.sum(
                  initial_mask[:, s] > -1
              )
              # recompute the prototypes
              # (average in subject s)
              lproto = [
                  np.mean(Fs[U[-1] == k], 0) for k in range(nb_parcel)
              ]
              lproto = np.array(lproto)
              lproto_anat = np.array(
                  [
                      np.mean(lac[U[-1] == k], 0)
                      for k in range(nb_parcel)
                  ]
              )
  
              LP.append(lproto)
              LPA.append(lproto_anat)
  
          # recompute the prototypes across subjects
          proto_mem = proto.copy()
          proto = np.mean(np.array(LP), 0)
          proto_anat = np.mean(np.array(LPA), 0)
          displ = np.sqrt(np.sum((proto_mem - proto) ** 2, 1).max())
          if verbose:
              print("energy", Energy, "displacement", displ)
  
          # recompute the topological model
          spatial_proto.set_field(proto_anat)
          spatial_proto.voronoi_diagram(proto_anat, domain.coord)
          spatial_proto.set_gaussian(proto_anat)
          spatial_proto.normalize()
  
          if displ < 1.0e-4 * dmax:
              break
      return U, proto_anat
  
file_path: nipy/nipy/labs/spatial_models/hierarchical_parcellation.py
imports: from warnings import warn
  import numpy
  from numpy.random import rand
  from nipy.algorithms.clustering.utils import kmeans, voronoi
  from nipy.algorithms.graph.field import Field
  from nipy.algorithms.graph.graph import wgraph_from_coo_matrix
  from .parcellation import MultiSubjectParcellation
label_desc: Use image node_id 8 for calling a global function with example usage: _optim_hparcel(feature, domain, graphs, nb_parcel, lamb, dmax, niter, initial_mask, chunksize, verbose) and returns: U,  proto_anat