docs: def pca(onarray, offarray, n=10, exchs=None, pc=False, mode='mean'):
      """Apply Principal Component Analysis (PCA) method to estimate baselines at each time.
  
      Args:
          onarray (decode.array): Decode array of on-point observations.
          offarray (decode.array): Decode array of off-point observations.
          n (int): The number of pricipal components.
          pc (bool): When True, this function also returns eigen vectors and their coefficients.
          mode (None or str): The way of correcting offsets.
              'mean': Mean.
              'median': Median.
              None: No correction.
  
      Returns:
          filtered (decode.array): Baseline-subtracted array.
          When pc is True:
              Ps (list(np.ndarray)): Eigen vectors.
              Cs (list(np.ndarray)): Coefficients.
      """
      logger = getLogger('decode.models.pca')
      logger.info('n_components exchs mode')
      if exchs is None:
          exchs = [16, 44, 46]
      logger.info('{} {} {}'.format(n, exchs, mode))
  
      offid = np.unique(offarray.scanid)
      onid  = np.unique(onarray.scanid)
  
      onarray   = onarray.copy() # Xarray
      onarray[:, exchs] = 0
      onvalues  = onarray.values
      onscanid  = onarray.scanid.values
      offarray  = offarray.copy() # Xarray
      offarray[:, exchs] = 0
      offvalues = offarray.values
      offscanid = offarray.scanid.values
  
      Ps, Cs     = [], []
      Xatm       = dc.full_like(onarray, onarray)
      Xatmvalues = Xatm.values
      model = TruncatedSVD(n_components=n)
      for i in onid:
          leftid  = np.searchsorted(offid, i) - 1
          rightid = np.searchsorted(offid, i)
  
          Xon = onvalues[onscanid == i]
          if leftid == -1:
              Xoff   = offvalues[offscanid == offid[rightid]]
              Xoff_m = getattr(np, 'nan'+mode)(Xoff, axis=0) if mode is not None else 0
              Xon_m  = Xoff_m
              model.fit(Xoff - Xoff_m)
          elif rightid == len(offid):
              Xoff   = offvalues[offscanid == offid[leftid]]
              Xoff_m = getattr(np, 'nan'+mode)(Xoff, axis=0) if mode is not None else 0
              Xon_m  = Xoff_m
              model.fit(Xoff - Xoff_m)
          else:
              Xoff_l  = offvalues[offscanid == offid[leftid]]
              Xoff_lm = getattr(np, 'nan'+mode)(Xoff_l, axis=0) if mode is not None else 0
              Xoff_r  = offvalues[offscanid == offid[rightid]]
              Xoff_rm = getattr(np, 'nan'+mode)(Xoff_r, axis=0) if mode is not None else 0
              Xon_m   = getattr(np, 'nan'+mode)(np.vstack([Xoff_l, Xoff_r]), axis=0) if mode is not None else 0
              model.fit(np.vstack([Xoff_l - Xoff_lm, Xoff_r - Xoff_rm]))
          P = model.components_
          C = model.transform(Xon - Xon_m)
  
          Xatmvalues[onscanid == i] = C @ P + Xon_m
          # Xatms.append(dc.full_like(Xon, C @ P + Xon_m.values))
          Ps.append(P)
          Cs.append(C)
  
      if pc:
          return Xatm, Ps, Cs
      else:
          return Xatm