kNNpy.Auxiliary.TPCF.3DTPCF_Tracer-Field

  1####################################################################################################
  2
  3#-------------------------------  Importing the required libraries  --------------------------------
  4
  5import numpy as np
  6import time
  7import MAS_library as MASL
  8import sys
  9import os
 10
 11#Necessary for relative imports (see https://stackoverflow.com/questions/34478398/import-local-function-from-a-module-housed-in-another-directory-with-relative-im)
 12module_path = os.path.abspath(os.path.join('../../'))
 13'''
 14@private
 15'''
 16if module_path not in sys.path:
 17    sys.path.append(module_path)
 18
 19from kNNpy.HelperFunctions import create_smoothed_field_dict_3D
 20
 21####################################################################################################
 22
 23# Function that returns the the 2-point Cross-Correlation Function between a set of discrete tracers and a continuous field using the stacking method
 24
 25def CrossCorr2pt(boxsize, bins, QueryPos, TracerPos, delta, thickness, R, kmin=None, kmax=None, Verbose=False):
 26    '''
 27    Calculates the Two-point Cross-correlation function between a set of tracers and a field. The interpolation can only be done using the 
 28    CIC-mass assignment scheme.
 29
 30    Parameters
 31    ----------
 32    boxsize: float
 33        The length (in Mpc/h) of the cubic box containing the tracers and the field
 34
 35    bins : float array of shape (m,)
 36        Set of m radial distances at which the 2PCF will be computed
 37        
 38    QueryPos : float array of shape (n_query,3) where n_query is the number of query points
 39            3D positions of the random query points inside the box, given in Cartesian coordinates (x, y, z) within the range         [0, boxsize]
 40
 41    TracerPos : float array of shape (n_pos, 3) where n_pos is the number of discrete tracers
 42        3D positions of the n tracers (e.g., galaxies) inside the box, given in Cartesian coordinates (x, y, z) within the range         [0, boxsize]
 43
 44    delta : float array of shape (ngrid, ngrid, ngrid)
 45        Smoothed overdensity field defined on a uniform grid with ngrid³ points
 46
 47    thickness : float, optional
 48        the thickness of the shell used for smoothing. Only use this parameter when 'Shell' filter is used. The smoothing is done using a shell with inner radius R-thickness/2 and outer radius R+thickness/2.
 49
 50    R : float, optional
 51        radial scale (in Mpc/h) at which the field is to be smoothed. Only use this parameter for real space smoothing.
 52
 53    kmin : float, optional
 54        the minimum value of the wavenumber. Only use this parameter when 'Top-Hat-k' filter is used.
 55
 56    kmax : float, optional
 57        the maximum value of the wavenumber. Only use this parameter when 'Top-Hat-k' filter is used.
 58    
 59    Verbose : bool, optional
 60        if set to True, the time taken to complete each step of the calculation will be printed, by default False.
 61    
 62    Returns
 63    -------
 64    xi : float array of shape (m,)
 65        The 2-point cross-correlation function (2PCF) between the tracer positions and the field, computed at each of the m radial bins.
 66    
 67    Raises
 68    ------
 69    ValueError
 70        if the given query points are not on a three-dimensional grid.
 71    ValueError
 72        if x,y, or z coordinate of any of the query points is not in ``(0, boxsize)``.
 73    ValueError
 74        if x,y, or z coordinate of any of the tracer points is not in ``(0, boxsize)``..
 75    ValueError
 76        if the given tracer points are not on a three-dimensional grid.
 77    ValueError
 78        if the given field is not a cubical three-dimensional array.
 79    '''
 80    #-----------------------------------------------------------------------------------------------
 81
 82    if Verbose: total_start_time = time.perf_counter()
 83
 84    #-----------------------------------------------------------------------------------------------
 85        
 86    # Check all inputs are consistent with the function requirement
 87
 88    if Verbose: print('Checking inputs ...')
 89
 90    if QueryPos.shape[1]!=3: 
 91        raise ValueError('Incorrect spatial dimension for query points: array containing the query point positions must be of shape (n_query, 3), ' \
 92        'where n_query is the number of query points.')
 93
 94    if np.any(QueryPos[:, 0] <= 0 or QueryPos[:, 0] >= boxsize):
 95        raise ValueError('Invalid query point position(s): please ensure 0 < x < boxsize.')
 96
 97    if np.any(QueryPos[:, 1] <= 0 or QueryPos[:, 1] >= boxsize):
 98        raise ValueError('Invalid query point position(s): please ensure 0 < y < boxsize.')
 99
100    if np.any(QueryPos[:, 2] <= 0 or QueryPos[:, 2] >= boxsize):
101        raise ValueError('Invalid query point position(s): please ensure 0 < z < boxsize.')
102
103    if np.any(TracerPos[:, 0] <= 0 or TracerPos[:, 0] >= boxsize):
104        raise ValueError('Invalid tracer point position(s): please ensure 0 < x < boxsize.')
105
106    if np.any(TracerPos[:, 1]<= 0 or TracerPos[:, 1]>= boxsize):
107        raise ValueError('Invalid tracer point position(s): please ensure 0 < y < boxsize.')
108
109    if np.any(TracerPos[:, 2]<= 0 or TracerPos[:, 2]>= boxsize):
110        raise ValueError('Invalid tracer point position(s): please ensure 0 < z < boxsize.')
111
112    if TracerPos.shape[1]!=3: 
113        raise ValueError('Incorrect spatial dimension for tracers: array containing the tracer positions must be of shape (n_tracer, 3), ' \
114        'where n_tracer is the number of tracers.')
115
116    if Verbose: print('\tdone.')
117
118    #-----------------------------------------------------------------------------------------------
119    
120    # Calculating the number of grid points along each axis of the field delta
121    shape = np.shape(delta)
122    if len(shape) != 3 or shape[0] != shape[1] or shape[1] != shape[2]:
123        raise ValueError("Error: Input array is not cubical (n, n, n).")
124    
125    smoothed_delta_dict= create_smoothed_field_dict_3D(field=delta, Filter='Shell', grid=QueryPos, Boxsize=boxsize, bins=bins, kmin=kmin, kmax=kmax, thickness=thickness, Verbose=Verbose)
126    
127    delta_interp = np.zeros((len(bins), len(TracerPos)))  # Shape (number_of_bins, number_of_tracers)
128
129    # Perform interpolation for each smoothed field
130    for i, R in enumerate(bins):
131        density_interpolated = np.zeros(TracerPos.shape[0], dtype=np.float32)
132        MASL.CIC_interp(smoothed_delta_dict[str(R)][i], boxsize, TracerPos, density_interpolated)
133        delta_interp[i] = density_interpolated
134
135    # Computing the 2-point Cross-Correlation Function by averaging over the interpolated field at the tracer positions
136    xi = np.zeros_like(bins)
137    for i, R in enumerate(bins):
138        xi[i] = np.mean(delta_interp[i])
139    
140    return xi
141
142#---------------------------------------------------------------------------------------------------