docstr: Bayesian parcel-based analysis.
  
  Given a sequence of independent images registered to a common
  space (for instance, a set of contrast images from a
  first-level fMRI analysis), perform a second-level analysis
  assuming constant effects throughout parcels defined from a
  given label image in reference space. Specifically, a model of
  the following form is assumed:
  
  Y = X * beta + variability,
  
  where Y denotes the input image sequence, X is a design
  matrix, and beta are parcel-wise parameter vectors. The
  algorithm computes the Bayesian posterior probability of beta
  in each parcel using an expectation propagation scheme.
  
  Parameters
  ----------
  con_imgs: sequence of nipy-like images
    Images input to the group analysis.
  parcel_img: nipy-like image
    Label image where each label codes for a parcel.
  parcel_info: sequence of arrays, optional
    A sequence of two arrays with same length equal to the
    number of distinct parcels consistently with the
    `parcel_img` argument. The first array gives parcel names
    and the second, parcel values, i.e., corresponding
    intensities in the associated parcel image. By default,
    parcel values are taken as
    `np.unique(parcel_img.get_fdata())` and parcel names are
    these values converted to strings.
  msk_img: nipy-like image, optional
    Binary mask to restrict analysis. By default, analysis is
    carried out on all parcels with nonzero value.
  vcon_imgs: sequence of nipy-like images, optional
    First-level variance estimates corresponding to
    `con_imgs`. This is useful if the input images are
    "noisy". By default, first-level variances are assumed to be
    zero.
  design_matrix: array, optional
    If None, a one-sample analysis model is used. Otherwise, an
    array with shape (n, p) where `n` matches the number of
    input scans, and `p` is the number of regressors.
  cvect: array, optional
    Contrast vector of interest. The method makes an inference
    on the contrast defined as the dot product cvect'*beta,
    where beta are the unknown parcel-wise effects. If None,
    `cvect` is assumed to be np.array((1,)). However, the
    `cvect` argument is mandatory if `design_matrix` is
    provided.
  fwhm: float, optional
    A parameter that represents the localization uncertainty in
    reference space in terms of the full width at half maximum
    of an isotropic Gaussian kernel.
  smooth_method: str, optional
    One of 'default' and 'spm'. Setting `smooth_method=spm`
    results in simply smoothing the input images using a
    Gaussian kernel, while the default method involves more
    complex smoothing in order to propagate spatial uncertainty
    into the inference process.
  res_path: str, optional
    An existing path to write output images. If None, no output
    is written.
  write_smoothed_images: bool, optional
    Specify whether smoothed images computed throughout the
    inference process are to be written on disk in `res_path`.
docstr_tok: ["Bayesian","parcel-based","analysis",".","Given","a","sequence","of","independent","images","registered","to","a","common","space","(","for","instance",",","a","set","of","contrast","images","from","a","first-level","fMRI","analysis",")",",","perform","a","second-level","analysis","assuming","constant","effects","throughout","parcels","defined","from","a","given","label","image","in","reference","space",".","Specifically",",","a","model","of","the","following","form","is","assumed",":","Y","=","X","*","beta","+","variability",",","where","Y","denotes","the","input","image","sequence",",","X","is","a","design","matrix",",","and","beta","are","parcel-wise","parameter","vectors",".","The","algorithm","computes","the","Bayesian","posterior","probability","of","beta","in","each","parcel","using","an","expectation","propagation","scheme",".","Parameters","--","--","--","--","--","con_imgs",":","sequence","of","nipy-like","images","Images","input","to","the","group","analysis",".","parcel_img",":","nipy-like","image","Label","image","where","each","label","codes","for","a","parcel",".","parcel_info",":","sequence","of","arrays",",","optional","A","sequence","of","two","arrays","with","same","length","equal","to","the","number","of","distinct","parcels","consistently","with","the","`","parcel_img","`","argument",".","The","first","array","gives","parcel","names","and","the","second",",","parcel","values",",","i.e.",",","corresponding","intensities","in","the","associated","parcel","image",".","By","default",",","parcel","values","are","taken","as","`","np.unique","(","parcel_img.get_fdata","(",")",")","`","and","parcel","names","are","these","values","converted","to","strings",".","msk_img",":","nipy-like","image",",","optional","Binary","mask","to","restrict","analysis",".","By","default",",","analysis","is","carried","out","on","all","parcels","with","nonzero","value",".","vcon_imgs",":","sequence","of","nipy-like","images",",","optional","First-level","variance","estimates","corresponding","to","`","con_imgs","`",".","This","is","useful","if","the","input","images","are","``","noisy","''",".","By","default",",","first-level","variances","are","assumed","to","be","zero",".","design_matrix",":","array",",","optional","If","None",",","a","one-sample","analysis","model","is","used",".","Otherwise",",","an","array","with","shape","(","n",",","p",")","where","`","n","`","matches","the","number","of","input","scans",",","and","`","p","`","is","the","number","of","regressors",".","cvect",":","array",",","optional","Contrast","vector","of","interest",".","The","method","makes","an","inference","on","the","contrast","defined","as","the","dot","product","cvect","'","*","beta",",","where","beta","are","the","unknown","parcel-wise","effects",".","If","None",",","`","cvect","`","is","assumed","to","be","np.array","(","(","1",",",")",")",".","However",",","the","`","cvect","`","argument","is","mandatory","if","`","design_matrix","`","is","provided",".","fwhm",":","float",",","optional","A","parameter","that","represents","the","localization","uncertainty","in","reference","space","in","terms","of","the","full","width","at","half","maximum","of","an","isotropic","Gaussian","kernel",".","smooth_method",":","str",",","optional","One","of","'default","'","and","'spm","'",".","Setting","`","smooth_method=spm","`","results","in","simply","smoothing","the","input","images","using","a","Gaussian","kernel",",","while","the","default","method","involves","more","complex","smoothing","in","order","to","propagate","spatial","uncertainty","into","the","inference","process",".","res_path",":","str",",","optional","An","existing","path","to","write","output","images",".","If","None",",","no","output","is","written",".","write_smoothed_images",":","bool",",","optional","Specify","whether","smoothed","images","computed","throughout","the","inference","process","are","to","be","written","on","disk","in","`","res_path","`","."]
code: def __init__(
      self,
      con_imgs,
      parcel_img,
      parcel_info=None,
      msk_img=None,
      vcon_imgs=None,
      design_matrix=None,
      cvect=None,
      fwhm=8,
      smooth_method="default",
      res_path=None,
      write_smoothed_images=False,
  ):
      """
      Bayesian parcel-based analysis.
  
      Given a sequence of independent images registered to a common
      space (for instance, a set of contrast images from a
      first-level fMRI analysis), perform a second-level analysis
      assuming constant effects throughout parcels defined from a
      given label image in reference space. Specifically, a model of
      the following form is assumed:
  
      Y = X * beta + variability,
  
      where Y denotes the input image sequence, X is a design
      matrix, and beta are parcel-wise parameter vectors. The
      algorithm computes the Bayesian posterior probability of beta
      in each parcel using an expectation propagation scheme.
  
      Parameters
      ----------
      con_imgs: sequence of nipy-like images
        Images input to the group analysis.
      parcel_img: nipy-like image
        Label image where each label codes for a parcel.
      parcel_info: sequence of arrays, optional
        A sequence of two arrays with same length equal to the
        number of distinct parcels consistently with the
        `parcel_img` argument. The first array gives parcel names
        and the second, parcel values, i.e., corresponding
        intensities in the associated parcel image. By default,
        parcel values are taken as
        `np.unique(parcel_img.get_fdata())` and parcel names are
        these values converted to strings.
      msk_img: nipy-like image, optional
        Binary mask to restrict analysis. By default, analysis is
        carried out on all parcels with nonzero value.
      vcon_imgs: sequence of nipy-like images, optional
        First-level variance estimates corresponding to
        `con_imgs`. This is useful if the input images are
        "noisy". By default, first-level variances are assumed to be
        zero.
      design_matrix: array, optional
        If None, a one-sample analysis model is used. Otherwise, an
        array with shape (n, p) where `n` matches the number of
        input scans, and `p` is the number of regressors.
      cvect: array, optional
        Contrast vector of interest. The method makes an inference
        on the contrast defined as the dot product cvect'*beta,
        where beta are the unknown parcel-wise effects. If None,
        `cvect` is assumed to be np.array((1,)). However, the
        `cvect` argument is mandatory if `design_matrix` is
        provided.
      fwhm: float, optional
        A parameter that represents the localization uncertainty in
        reference space in terms of the full width at half maximum
        of an isotropic Gaussian kernel.
      smooth_method: str, optional
        One of 'default' and 'spm'. Setting `smooth_method=spm`
        results in simply smoothing the input images using a
        Gaussian kernel, while the default method involves more
        complex smoothing in order to propagate spatial uncertainty
        into the inference process.
      res_path: str, optional
        An existing path to write output images. If None, no output
        is written.
      write_smoothed_images: bool, optional
        Specify whether smoothed images computed throughout the
        inference process are to be written on disk in `res_path`.
      """
      self.smooth_method = smooth_method
      self.con_imgs = con_imgs
      self.vcon_imgs = vcon_imgs
      self.n_subjects = len(con_imgs)
      if self.vcon_imgs is not None:
          if not self.n_subjects == len(vcon_imgs):
              raise ValueError(
                  "List of contrasts and variances"
                  " do not have the same length"
              )
      if msk_img is None:
          self.msk = None
      else:
          self.msk = msk_img.get_fdata().astype(bool).squeeze()
      self.res_path = res_path
  
      # design matrix
      if design_matrix is None:
          self.design_matrix = np.ones(self.n_subjects)
          self.cvect = np.ones((1,))
          if cvect is not None:
              raise ValueError("No contrast vector expected")
      else:
          self.design_matrix = np.asarray(design_matrix)
          if cvect is None:
              raise ValueError(
                  "`cvect` cannot be None with"
                  " provided design matrix"
              )
          self.cvect = np.asarray(cvect)
          if not self.design_matrix.shape[0] == self.n_subjects:
              raise ValueError(
                  "Design matrix shape is inconsistent"
                  " with number of input images"
              )
          if not len(self.cvect) == self.design_matrix.shape[1]:
              raise ValueError(
                  "Design matrix shape is inconsistent"
                  " with provided `cvect`"
              )
  
      # load the parcellation and resample it at the appropriate
      # resolution
      self.reference = parcel_img.reference
      self.parcel_full_res = (
          parcel_img.get_fdata().astype("uintp").squeeze()
      )
      self.affine_full_res = xyz_affine(parcel_img)
      parcel_img = make_xyz_image(
          self.parcel_full_res, self.affine_full_res, self.reference
      )
      self.affine = xyz_affine(self.con_imgs[0])
      parcel_img_rsp = resample(
          parcel_img,
          reference=(self.con_imgs[0].shape, self.affine),
          interp_order=0,
      )
      self.parcel = parcel_img_rsp.get_fdata().astype("uintp").squeeze()
      if self.msk is None:
          self.msk = self.parcel > 0
  
      # get parcel labels and values
      if parcel_info is None:
          self._parcel_values = np.unique(self.parcel)
          self._parcel_labels = self._parcel_values.astype(str)
      else:
          self._parcel_labels = np.asarray(parcel_info[0]).astype(str)
          self._parcel_values = np.asarray(parcel_info[1])
  
      # determine smoothing kernel size, which involves converting
      # the input full-width-at-half-maximum parameter given in mm
      # to standard deviation in voxel units.
      orient = io_orientation(self.affine)[:, 0].astype(int)
      # `orient` is an array, so this slicing leads to advanced indexing.
      voxsize = np.abs(self.affine[orient, list(range(3))])
      self.sigma = np.maximum(fwhm2sigma(fwhm) / voxsize, SIGMA_MIN)
  
      # run approximate belief propagation
      self._smooth_images(write_smoothed_images)
      self._voxel_level_inference()
      self._parcel_level_inference()
  
imports: import warnings
  from os.path import join
  import numpy
  import scipy.ndimage
  import scipy.stats
  from nibabel import io_orientation
  from ...None import save_image
  from ...core.image.image_spaces import make_xyz_image, xyz_affine
  from ..kernel_smooth import fwhm2sigma
  from ..registration import resample
  from ..statistics.bayesian_mixed_effects import two_level_glm
  from ..statistics.histogram import histogram
label_desc: Use image node_id 1 to create a new ParcelAnalysis object  with example: obj = ParcelAnalysis(con_imgs, parcel_img, parcel_info, msk_img, vcon_imgs, design_matrix, cvect, fwhm, smooth_method, res_path, write_smoothed_images)