text: r"""
  Module containing the `~halotools.mock_observables.marked_tpcf` function used to
  calculate the marked two-point correlation function.
  """
  
  from __future__ import absolute_import, division, print_function, unicode_literals
  
  import numpy as np
  from astropy.utils.misc import NumpyRNGContext
  
  from .clustering_helpers import process_optional_input_sample2
  
  from ..mock_observables_helpers import (enforce_sample_has_correct_shape,
      get_separation_bins_array, get_period, get_num_threads)
  from ..pair_counters.mesh_helpers import _enforce_maximum_search_length
  
  from ..pair_counters import npairs_3d, marked_npairs_3d
  
  from ...custom_exceptions import HalotoolsError
  
  
  __all__ = ['marked_tpcf']
  __author__ = ['Duncan Campbell']
  
  np.seterr(divide='ignore', invalid='ignore')  # ignore divide by zero in e.g. DD/RR
  
  
  def marked_tpcf(sample1, rbins, sample2=None,
          marks1=None, marks2=None, period=None, do_auto=True, do_cross=True,
          num_threads=1, weight_func_id=1,
          normalize_by='random_marks', iterations=1, randomize_marks=None, seed=None):
      r"""
      Calculate the real space marked two-point correlation function, :math:`\mathcal{M}(r)`.
  
      Example calls to this function appear in the documentation below.
      See the :ref:`mock_obs_pos_formatting` documentation page for
      instructions on how to transform your coordinate position arrays into the
      format accepted by the ``sample1`` and ``sample2`` arguments.
  
      Parameters
      ----------
      sample1 : array_like
          Npts1 x 3 numpy array containing 3-D positions of points.
          See the :ref:`mock_obs_pos_formatting` documentation page, or the
          Examples section below, for instructions on how to transform
          your coordinate position arrays into the
          format accepted by the ``sample1`` and ``sample2`` arguments.
          Length units are comoving and assumed to be in Mpc/h, here and throughout Halotools.
  
      rbins : array_like
          array of boundaries defining the real space radial bins in which pairs are counted.
          Length units are comoving and assumed to be in Mpc/h, here and throughout Halotools.
  
      sample2 : array_like, optional
          Npts2 x 3 array containing 3-D positions of points.
          Passing ``sample2`` as an input permits the calculation of
          the cross-correlation function. Default is None, in which case only the
          auto-correlation function will be calculated.
  
      marks1 : array_like, optional
          len(sample1) x N_marks array of marks.  The supplied marks array must have the
          appropriate shape for the chosen ``weight_func_id`` (see Notes for requirements).  If this
          parameter is not specified, it is set to numpy.ones((len(sample1), N_marks)).
  
      marks2 : array_like, optional
          len(sample2) x N_marks array of marks.  The supplied marks array must have the
          appropriate shape for the chosen ``weight_func_id`` (see Notes for requirements).  If this
          parameter is not specified, it is set to numpy.ones((len(sample2), N_marks)).
  
      period : array_like, optional
          Length-3 sequence defining the periodic boundary conditions
          in each dimension. If you instead provide a single scalar, Lbox,
          period is assumed to be the same in all Cartesian directions.
          If set to None (the default option), PBCs are set to infinity.
          Length units are comoving and assumed to be in Mpc/h, here and throughout Halotools.
  
      do_auto : boolean, optional
          Boolean determines whether the auto-correlation function will
          be calculated and returned. Default is True.
  
      do_cross : boolean, optional
          Boolean determines whether the cross-correlation function will
          be calculated and returned. Only relevant when ``sample2`` is also provided.
          Default is True for the case where ``sample2`` is provided, otherwise False.
  
      num_threads : int, optional
          Number of threads to use in calculation, where parallelization is performed
          using the python ``multiprocessing`` module. Default is 1 for a purely serial
          calculation, in which case a multiprocessing Pool object will
          never be instantiated. A string 'max' may be used to indicate that
          the pair counters should use all available cores on the machine.
  
      weight_func_id : int, optional
          Integer ID indicating which marking function should be used.
          See notes for a list of available marking functions.
  
      normalize_by : string, optional
          A string indicating how to normailze the weighted pair counts in the marked
          correlation function calculation.  Options are: 'random_marks' or 'number_counts'.
          See Notes for more detail.
  
      iterations : int, optional
          integer indicating the number of times to calculate the random weights,
          taking the mean of the outcomes.  Only applicable if ``normalize_by`` is set
          to 'random_marks'.  See Notes for further explanation.
  
      randomize_marks : array_like, optional
          Boolean array of length N_marks indicating which elements should be randomized
          when calculating the random weighted pair counts.  Default is [True]*N_marks.
          This parameter is only applicable if ``normalize_by`` is set to 'random_marks'.
          See Notes for more detail.
  
      seed : int, optional
          Random number seed used to shuffle the marks
          and to randomly downsample data, if applicable.
          Default is None, in which case downsampling and shuffling will be stochastic.
  
      Returns
      -------
      marked_correlation_function(s) : numpy.array
          *len(rbins)-1* length array containing the marked correlation function
          :math:`\mathcal{M}(r)` computed in each of the bins defined by ``rbins``.
  
          .. math::
              \mathcal{M}(r) \equiv \mathrm{WW}(r) / \mathrm{XX}(r),
  
          where :math:`\mathrm{WW}(r)` is the weighted number of pairs with separations
          equal to :math:`r`, and :math:`\mathrm{XX}(r)` is dependent on the choice of the
          ``normalize_by`` parameter.  If ``normalize_by`` is 'random_marks'
          :math:`XX \equiv \mathcal{RR}`, the weighted pair counts where the marks have
          been randomized marks.  If ``normalize_by`` is 'number_counts'
          :math:`XX \equiv DD`, the unweighted pair counts.
          See Notes for more detail.
  
          If ``sample2`` is passed as input, three arrays of length *len(rbins)-1* are
          returned:
  
          .. math::
              \mathcal{M}_{11}(r), \ \mathcal{M}_{12}(r), \ \mathcal{M}_{22}(r),
  
          the autocorrelation of ``sample1``, the cross-correlation between ``sample1`` and
          ``sample2``, and the autocorrelation of ``sample2``.  If ``do_auto`` or
          ``do_cross`` is set to False, the appropriate result(s) is not returned.
  
      Notes
      -----
      Pairs are counted using
      `~halotools.mock_observables.pair_counters.marked_npairs_3d`.
  
      If the ``period`` argument is passed in, the ith coordinate of all points
      must be between 0 and period[i].
  
      ``normalize_by`` indicates how to calculate :math:`\mathrm{XX}`.  If ``normalize_by``
      is 'random_marks', then :math:`\mathrm{XX} \equiv \mathcal{RR}`, and
      :math:`\mathcal{RR}` is calculated by randomizing the marks among points according
      to the ``randomize_marks`` mask.  This marked correlation function is then:
  
      .. math::
          \mathcal{M}(r) \equiv \frac{\sum_{ij}f(m_i,m_j)}{\sum_{kl}f(m_k,m_l)}
  
      where the sum in the numerator is of pairs :math:`i,j` with separation :math:`r`,
      and marks :math:`m_i,m_j`.  :math:`f()` is the marking function, ``weight_func_id``.  The sum
      in the denominator is over an equal number of random pairs :math:`k,l`. The
      calculation of this sum can be done multiple times, by setting the ``iterations``
      parameter. The mean of the sum is then taken amongst iterations and used in the
      calculation.
  
      If ``normalize_by`` is 'number_counts', then :math:`\mathrm{XX} \equiv \mathrm{DD}`
      is calculated by counting total number of pairs using
      `~halotools.mock_observables.pair_counters.npairs_3d`.
      This is:
  
      .. math::
          \mathcal{M}(r) \equiv \frac{\sum_{ij}f(m_i,m_j)}{\sum_{ij} 1},
  
      There are multiple marking functions available.  In general, each requires a different
      number of marks per point, N_marks.  The marking function gets passed two vectors
      per pair, w1 and w2, of length N_marks and return a float.  The available marking
      functions, ``weight_func_id`` and the associated integer ID numbers are:
  
      #. multiplicaitive weights (N_marks = 1)
          .. math::
              f(w_1,w_2) = w_1[0] \times w_2[0]
  
      #. summed weights (N_marks = 1)
          .. math::
              f(w_1,w_2) = w_1[0] + w_2[0]
  
      #. equality weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_1[1]\times w_2[1] & : w_1[0] = w_2[0] \\
                      0.0 & : w_1[0] \neq w_2[0] \\
                  \end{array}
                  \right.
  
      #. inequality weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_1[1]\times w_2[1] & : w_1[0] \neq w_2[0] \\
                      0.0 & : w_1[0] = w_2[0] \\
                  \end{array}
                  \right.
  
      #. greater than weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_1[1]\times w_2[1] & : w_2[0] > w_1[0] \\
                      0.0 & : w_2[0] \leq w_1[0] \\
                  \end{array}
                  \right.
  
      #. less than weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_1[1]\times w_2[1] & : w_2[0] < w_1[0] \\
                      0.0 & : w_2[0] \geq w_1[0] \\
                  \end{array}
                  \right.
  
      #. greater than tolerance weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_2[1] & : w_2[0]>(w_1[0]+w_1[1]) \\
                      0.0 & : w_2[0] \leq (w_1[0]+w_1[1]) \\
                  \end{array}
                  \right.
  
      #. less than tolerance weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_2[1] & : w_2[0]<(w_1[0]+w_1[1]) \\
                      0.0 & : w_2[0] \geq (w_1[0]+w_1[1]) \\
                  \end{array}
                  \right.
  
      #. tolerance weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_2[1] & : |w_1[0]-w_2[0]|<w_1[1] \\
                      0.0 & : |w_1[0]-w_2[0]| \geq w_1[1] \\
                  \end{array}
                  \right.
  
      #. exclusion weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_2[1] & : |w_1[0]-w_2[0]|>w_1[1] \\
                      0.0 & : |w_1[0]-w_2[0]| \leq w_1[1] \\
                  \end{array}
                  \right.
  
      #. ratio weights (N_marks = 2)
          .. math::
              f(w_1,w_2) =
                  \left \{
                  \begin{array}{ll}
                      w_2[1] & : w2[0] > w1[0]*w1[1] \\
                      0.0 & : otherwise \\
                  \end{array}
                  \right.
  
      Examples
      --------
      For demonstration purposes we create a randomly distributed set of points within a
      periodic unit cube.
  
      >>> Npts = 1000
      >>> Lbox = 1.0
      >>> period = np.array([Lbox,Lbox,Lbox])
  
      >>> x = np.random.random(Npts)
      >>> y = np.random.random(Npts)
      >>> z = np.random.random(Npts)
  
      We transform our *x, y, z* points into the array shape used by the function by
      taking the transpose of the result of `numpy.vstack`. This boilerplate transformation
      is used throughout the `~halotools.mock_observables` sub-package:
  
      >>> coords = np.vstack((x,y,z)).T
  
      Assign random floats in the range [0,1] to the points to use as the marks:
  
      >>> marks = np.random.random(Npts)
  
      Use the multiplicative marking function:
  
      >>> rbins = np.logspace(-2,-1,10)
      >>> MCF = marked_tpcf(coords, rbins, marks1=marks, period=period, normalize_by='number_counts', weight_func_id=1)
  
      The result should be consistent with :math:`\langle {\rm mark}\rangle^2` at all *r*
      within the statistical errors.
      """
  
      # process parameters
      function_args = (sample1, rbins, sample2, marks1, marks2,
          period, do_auto, do_cross, num_threads,
          weight_func_id, normalize_by, iterations, randomize_marks, seed)
      sample1, rbins, sample2, marks1, marks2, period, do_auto, do_cross, num_threads,\
          weight_func_id, normalize_by, _sample1_is_sample2, PBCs,\
          randomize_marks = _marked_tpcf_process_args(*function_args)
  
      # calculate marked pairs
      W1W1, W1W2, W2W2 = marked_pair_counts(sample1, sample2, rbins, period,
          num_threads, do_auto, do_cross, marks1, marks2, weight_func_id, _sample1_is_sample2)
  
      if normalize_by == 'number_counts':
          R1R1, R1R2, R2R2 = pair_counts(sample1, sample2, rbins, period,
              num_threads, do_auto, do_cross, _sample1_is_sample2, None, None)
      # calculate randomized marked pairs
      elif normalize_by == 'random_marks':
          if iterations > 1:
              # create storage arrays of the right shape
              R1R1 = np.zeros((iterations, len(rbins)-1))
              R1R2 = np.zeros((iterations, len(rbins)-1))
              R2R2 = np.zeros((iterations, len(rbins)-1))
              for i in range(iterations):
                  # get arrays to randomize marks
                  with NumpyRNGContext(seed):
                      permutate1 = np.random.permutation(np.arange(0, len(sample1)))
                      permutate2 = np.random.permutation(np.arange(0, len(sample2)))
                  R1R1[i, :], R1R2[i, :], R2R2[i, :] = random_counts(
                      sample1, sample2, rbins, period, num_threads,
                      do_auto, do_cross, marks1, marks2, weight_func_id,
                      _sample1_is_sample2, permutate1, permutate2, randomize_marks)
  
              R1R1 = np.median(R1R1, axis=0)
              R1R2 = np.median(R1R2, axis=0)
              R2R2 = np.median(R2R2, axis=0)
          else:
              # get arrays to randomize marks
              with NumpyRNGContext(seed):
                  permutate1 = np.random.permutation(np.arange(0, len(sample1)))
                  permutate2 = np.random.permutation(np.arange(0, len(sample2)))
              R1R1, R1R2, R2R2 = random_counts(sample1, sample2, rbins, period,
                  num_threads, do_auto, do_cross, marks1, marks2, weight_func_id,
                  _sample1_is_sample2, permutate1, permutate2, randomize_marks)
  
      # return results
      if _sample1_is_sample2:
          M_11 = W1W1/R1R1
          return M_11
      else:
          if (do_auto is True) & (do_cross is True):
              M_11 = W1W1/R1R1
              M_12 = W1W2/R1R2
              M_22 = W2W2/R2R2
              return M_11, M_12, M_22
          elif (do_cross is True):
              M_12 = W1W2/R1R2
              return M_12
          elif (do_auto is True):
              M_11 = W1W1/R1R1
              M_22 = W2W2/R2R2
              return M_11, M_22
  
  
  def marked_pair_counts(sample1, sample2, rbins, period, num_threads,
          do_auto, do_cross, marks1, marks2, weight_func_id, _sample1_is_sample2):
      """
      Count weighted data pairs.
      """
  
      if do_auto is True:
          D1D1 = marked_npairs_3d(sample1, sample1, rbins,
              weights1=marks1, weights2=marks1,
              weight_func_id=weight_func_id, period=period, num_threads=num_threads)
          D1D1 = np.diff(D1D1)
      else:
          D1D1 = None
          D2D2 = None
  
      if _sample1_is_sample2:
          D1D2 = D1D1
          D2D2 = D1D1
      else:
          if do_cross is True:
              D1D2 = marked_npairs_3d(sample1, sample2, rbins,
                  weights1=marks1, weights2=marks2, weight_func_id=weight_func_id,
                  period=period, num_threads=num_threads)
              D1D2 = np.diff(D1D2)
          else:
              D1D2 = None
          if do_auto is True:
              D2D2 = marked_npairs_3d(sample2, sample2, rbins,
                  weights1=marks2, weights2=marks2, weight_func_id=weight_func_id,
                  period=period, num_threads=num_threads)
              D2D2 = np.diff(D2D2)
          else:
              D2D2 = None
  
      return D1D1, D1D2, D2D2
  
  
  def random_counts(sample1, sample2, rbins, period, num_threads,
          do_auto, do_cross, marks1, marks2, weight_func_id,
          _sample1_is_sample2, permutate1, permutate2, randomize_marks):
      """
      Count random weighted data pairs.
      """
  
      permuted_marks1 = marks1
      permuted_marks2 = marks2
      for i in range(marks1.shape[1]):
          if randomize_marks[i]:
              permuted_marks1[:, i] = marks1[permutate1, i]
      for i in range(marks2.shape[1]):
          if randomize_marks[i]:
              permuted_marks2[:, i] = marks2[permutate2, i]
  
      if do_auto is True:
          R1R1 = marked_npairs_3d(sample1, sample1, rbins,
              weights1=marks1, weights2=permuted_marks1,
              weight_func_id=weight_func_id, period=period, num_threads=num_threads)
          R1R1 = np.diff(R1R1)
      else:
          R1R1 = None
          R2R2 = None
  
      if _sample1_is_sample2:
          R1R2 = R1R1
          R2R2 = R1R1
      else:
          if do_cross is True:
              R1R2 = marked_npairs_3d(sample1, sample2, rbins,
                  weights1=permuted_marks1,
                  weights2=permuted_marks2,
                  weight_func_id=weight_func_id, period=period, num_threads=num_threads)
              R1R2 = np.diff(R1R2)
          else:
              R1R2 = None
          if do_auto is True:
              R2R2 = marked_npairs_3d(sample2, sample2, rbins,
                  weights1=marks2, weights2=permuted_marks2,
                  weight_func_id=weight_func_id, period=period, num_threads=num_threads)
              R2R2 = np.diff(R2R2)
          else:
              R2R2 = None
  
      return R1R1, R1R2, R2R2
  
  
  def pair_counts(sample1, sample2, rbins, period, num_threads, do_auto, do_cross,
          _sample1_is_sample2, approx_cell1_size, approx_cell2_size):
      """
      Count data-data pairs.
      """
      if do_auto is True:
          D1D1 = npairs_3d(sample1, sample1, rbins, period=period, num_threads=num_threads,
                        approx_cell1_size=approx_cell1_size,
                        approx_cell2_size=approx_cell1_size)
          D1D1 = np.diff(D1D1)
      else:
          D1D1 = None
          D2D2 = None
  
      if _sample1_is_sample2:
          D1D2 = D1D1
          D2D2 = D1D1
      else:
          if do_cross is True:
              D1D2 = npairs_3d(sample1, sample2, rbins, period=period,
                            num_threads=num_threads,
                            approx_cell1_size=approx_cell1_size,
                            approx_cell2_size=approx_cell2_size)
              D1D2 = np.diff(D1D2)
          else:
              D1D2 = None
          if do_auto is True:
              D2D2 = npairs_3d(sample2, sample2, rbins, period=period,
                            num_threads=num_threads,
                            approx_cell1_size=approx_cell2_size,
                            approx_cell2_size=approx_cell2_size)
              D2D2 = np.diff(D2D2)
          else:
              D2D2 = None
  
      return D1D1, D1D2, D2D2
  
  
  def _marked_tpcf_process_args(sample1, rbins, sample2, marks1, marks2,
          period, do_auto, do_cross, num_threads,
          wfunc, normalize_by, iterations, randomize_marks, seed):
      """
      Private method to do bounds-checking on the arguments passed to
      `~halotools.mock_observables.marked_tpcf`.
      """
  
      sample1 = enforce_sample_has_correct_shape(sample1)
  
      sample2, _sample1_is_sample2, do_cross = process_optional_input_sample2(
          sample1, sample2, do_cross)
  
      # process wfunc parameter
      try:
          int(wfunc) == wfunc
      except:
          msg = ("\n `wfunc` parameter must be an integer ID of the desired function.")
          raise ValueError(msg)
  
      # process normalize_by parameter
      if normalize_by not in ['random_marks', 'number_counts']:
          msg = ("\n `normalize_by` parameter not recognized.")
          raise ValueError(msg)
  
      # process marks
      if marks1 is not None:
          marks1 = np.atleast_1d(marks1).astype(float)
      else:
          marks1 = np.ones(len(sample1)).astype(float)
      if marks2 is not None:
          marks2 = np.atleast_1d(marks2).astype(float)
      else:
          marks2 = np.ones(len(sample2)).astype(float)
  
      if marks1.ndim == 1:
          npts1 = len(marks1)
          marks1 = marks1.reshape((npts1, 1))
      elif marks1.ndim == 2:
          pass
      else:
          ndim1 = marks1.ndim
          msg = ("\n You must either pass in a 1-D or 2-D array \n"
                 "for the input `marks1`. \n"
                 "The `pair_counters._wnpairs_process_weights` function received \n"
                 "a `marks1` array of dimension %i")
          raise HalotoolsError(msg % ndim1)
  
      if marks2.ndim == 1:
          npts2 = len(marks2)
          marks2 = marks2.reshape((npts2, 1))
      elif marks2.ndim == 2:
          pass
      else:
          ndim2 = marks2.ndim
          msg = ("\n You must either pass in a 1-D or 2-D array \n"
                 "for the input `marks2`. \n"
                 "The `pair_counters._wnpairs_process_weights` function received \n"
                 "a `marks2` array of dimension %i")
          raise HalotoolsError(msg % ndim2)
  
      # check for consistency between marks and samples
      if len(marks1) != len(sample1):
          msg = ("\n `marks1` must have same length as `sample1`.")
          raise HalotoolsError(msg)
      if len(marks2) != len(sample2):
          msg = ("\n `marks2` must have same length as `sample2`.")
          raise HalotoolsError(msg)
  
      if randomize_marks is not None:
          randomize_marks = np.atleast_1d(randomize_marks)
      else:
          randomize_marks = np.array([True]*marks1.shape[1])
  
      if randomize_marks.ndim == 1:
          if len(randomize_marks) != marks1.shape[1]:
              msg = ("\n `randomize_marks` must have same length \n"
                     " as the number of weights per point.")
              raise HalotoolsError(msg)
      else:
          msg = ("\n `randomize_marks` must be one dimensional.")
          raise HalotoolsError(msg)
  
      rbins = get_separation_bins_array(rbins)
      rmax = np.amax(rbins)
  
      period, PBCs = get_period(period)
  
      _enforce_maximum_search_length(rmax, period)
  
      try:
          assert do_auto == bool(do_auto)
          assert do_cross == bool(do_cross)
      except:
          msg = "`do_auto` and `do_cross` keywords must be boolean-valued."
          raise ValueError(msg)
  
      num_threads = get_num_threads(num_threads)
  
      return sample1, rbins, sample2, marks1, marks2, period, do_auto, do_cross,\
          num_threads, wfunc, normalize_by, _sample1_is_sample2, PBCs, randomize_marks
  