|
| 1 | +# To run: mpiexec -np 16 python pkrun.py --type ELG... |
| 2 | + |
| 3 | +import os |
| 4 | +import argparse |
| 5 | +import logging |
| 6 | +import numpy as np |
| 7 | +from astropy.table import Table, vstack |
| 8 | +from matplotlib import pyplot as plt |
| 9 | + |
| 10 | +# To install pypower: python -m pip install git+https://github.com/adematti/pypower |
| 11 | +# To install mockfactory: python -m pip install git+https://github.com/adematti/mockfactory |
| 12 | +from pypower import CatalogFFTPower, utils, setup_logging |
| 13 | +from mockfactory import Catalog |
| 14 | +from LSS.tabulated_cosmo import TabulatedDESI |
| 15 | +cosmo = TabulatedDESI() |
| 16 | +distance = cosmo.comoving_radial_distance |
| 17 | + |
| 18 | +os.environ['NUMEXPR_MAX_THREADS'] = '8' |
| 19 | +parser = argparse.ArgumentParser() |
| 20 | +parser.add_argument("--type", help="tracer type to be selected") |
| 21 | +parser.add_argument("--basedir", help="where to find catalogs",default='/global/cfs/cdirs/desi/survey/catalogs') |
| 22 | +parser.add_argument("--version", help="catalog version",default='test') |
| 23 | +parser.add_argument("--verspec",help="version for redshifts",default='everest') |
| 24 | +parser.add_argument("--survey",help="e.g., SV3 or main",default='SV3') |
| 25 | +parser.add_argument("--nran",help="number of random files to combine together (1-18 available)",default=10) |
| 26 | +parser.add_argument("--weight_type",help="types of weights to use; use angular_bitwise for PIP; default just uses WEIGHT column",default='default') |
| 27 | + |
| 28 | +#only relevant for reconstruction |
| 29 | +parser.add_argument("--rectype",help="IFT or MG supported so far",default='IFT') |
| 30 | +parser.add_argument("--convention",help="recsym or reciso supported so far",default='reciso') |
| 31 | + |
| 32 | +setup_logging() |
| 33 | +args = parser.parse_args() |
| 34 | + |
| 35 | +ttype = args.type |
| 36 | +basedir = args.basedir |
| 37 | +version = args.version |
| 38 | +specrel = args.verspec |
| 39 | +survey = args.survey |
| 40 | +nran = int(args.nran) |
| 41 | +weight_type = args.weight_type |
| 42 | + |
| 43 | +edges = {'min':0., 'step':0.005} |
| 44 | + |
| 45 | +dirpk = os.environ['CSCRATCH']+'/'+survey+'pk/' |
| 46 | + |
| 47 | +utils.mkdir(dirpk) |
| 48 | +print('made '+dirpk) |
| 49 | + |
| 50 | +lssdir = basedir+'/'+survey+'/LSS/'+specrel+'/LSScats/' |
| 51 | +dirname = lssdir + version |
| 52 | +#dirout = svdir+'LSScats/'+version+'/' |
| 53 | + |
| 54 | +zmask = [''] |
| 55 | +minn = 0 |
| 56 | + |
| 57 | +subt = None |
| 58 | + |
| 59 | +if ttype[:3] == 'LRG': |
| 60 | + zl = [0.4,0.6,0.8,1.1] |
| 61 | + |
| 62 | + |
| 63 | +if ttype[:3] == 'ELG':# or type == 'ELG_HIP': |
| 64 | + #minn = 5 |
| 65 | + zl = [0.8,1.1,1.5] |
| 66 | + #zmask = ['','_zmask'] |
| 67 | + |
| 68 | + #zmin = 0.8 |
| 69 | + #zmax = 1.6 |
| 70 | + |
| 71 | + |
| 72 | +if ttype == 'QSO': |
| 73 | + zl = [0.8,1.1,1.5,2.1] |
| 74 | + #zmin = 1. |
| 75 | + #zmax = 2.1 |
| 76 | + |
| 77 | +if ttype == 'QSOh': |
| 78 | + zl = [2.1,3.5] |
| 79 | + ttype = 'QSO' |
| 80 | + #zmin = 1. |
| 81 | + #zmax = 2.1 |
| 82 | + |
| 83 | +if ttype[:3] == 'BGS': |
| 84 | + #minn = 2 |
| 85 | + zl = [0.1,0.3,0.5] |
| 86 | + #zmin = 0.1 |
| 87 | + #zmax = 0.5 |
| 88 | + |
| 89 | +if ttype[:3] == 'BGS' and ttype[-1] == 'l': |
| 90 | + #minn = 2 |
| 91 | + zl = [0.1,0.3] |
| 92 | + ttype = ttype[:-1] |
| 93 | + #zmin = 0.1 |
| 94 | + #zmax = 0.5 |
| 95 | + |
| 96 | +if ttype[:3] == 'BGS' and ttype[-1] == 'h': |
| 97 | + #minn = 2 |
| 98 | + zl = [0.3,0.5] |
| 99 | + ttype = ttype[:-1] |
| 100 | + #zmin = 0.1 |
| 101 | + #zmax = 0.5 |
| 102 | + |
| 103 | + |
| 104 | +wa = '' |
| 105 | +if survey in ['main', 'DA02']: |
| 106 | + wa = 'zdone' |
| 107 | + |
| 108 | +def compute_power_spectrum(edges, tracer='LRG', region='_N', nrandoms=4, zlim=(0., np.inf), weight_type=None, ells=(0, 2, 4), boxsize=5000., nmesh=1024, dtype='f4'): |
| 109 | + if ttype == 'ELGrec' or ttype == 'LRGrec': |
| 110 | + data_fn = os.path.join(dirname, tracer+wa+ region+'_clustering_'+args.rectype+args.convention+'.dat.fits') |
| 111 | + data = Catalog.load_fits(data_fn) |
| 112 | + |
| 113 | + randoms_fn = os.path.join(dirname, tracer+wa+ region+'_clustering_'+args.rectype+args.convention+'.ran.fits') |
| 114 | + randoms = Catalog.load_fits(randoms_fn) |
| 115 | + else: |
| 116 | + data_fn = os.path.join(dirname, '{}{}_clustering.dat.fits'.format(tracer+wa, region)) |
| 117 | + data = Catalog.load_fits(data_fn) |
| 118 | + |
| 119 | + randoms_fn = [os.path.join(dirname, '{}{}_{:d}_clustering.ran.fits'.format(tracer+wa, region, iran)) for iran in range(nrandoms)] |
| 120 | + randoms = Catalog.concatenate(*(Catalog.load_fits(fn) for fn in randoms_fn), keep_order=False) |
| 121 | + |
| 122 | + def get_positions_weights(catalog, name='data'): |
| 123 | + mask = (catalog['Z'] >= zlim[0]) & (catalog['Z'] < zlim[1]) |
| 124 | + positions = [catalog['RA'][mask], catalog['DEC'][mask], distance(catalog['Z'][mask])] |
| 125 | + nmask = catalog.mpicomm.allreduce(mask.sum()) |
| 126 | + if catalog.mpicomm.rank == 0: |
| 127 | + catalog.log_info('Using {} rows for {}'.format(nmask, name)) |
| 128 | + #if weight_type is None: |
| 129 | + # weights = None |
| 130 | + #else: |
| 131 | + weights = np.ones_like(positions[0]) |
| 132 | + if name == 'data': |
| 133 | + if 'photometric' in weight_type: |
| 134 | + rfweight = RFWeight(tracer=tracer) |
| 135 | + weights *= rfweight(positions[0], positions[1]) |
| 136 | + if 'zfail' in weight_type: |
| 137 | + weights *= catalog['WEIGHT_ZFAIL'][mask] |
| 138 | + if 'default' in weight_type: |
| 139 | + weights *= catalog['WEIGHT'][mask] |
| 140 | + if 'completeness' in weight_type: |
| 141 | + weights *= catalog['WEIGHT'][mask]/catalog['WEIGHT_ZFAIL'][mask] |
| 142 | + elif 'bitwise' in weight_type: |
| 143 | + weights = list(catalog['BITWEIGHTS'][mask].T) + [weights] |
| 144 | + return positions, weights |
| 145 | + |
| 146 | + data_positions, data_weights = get_positions_weights(data, name='data') |
| 147 | + randoms_positions, randoms_weights = get_positions_weights(randoms, name='randoms') |
| 148 | + |
| 149 | + result = CatalogFFTPower(data_positions1=data_positions, data_weights1=data_weights, |
| 150 | + randoms_positions1=randoms_positions, randoms_weights1=randoms_weights, |
| 151 | + edges=edges, ells=ells, boxsize=boxsize, nmesh=nmesh, resampler='tsc', interlacing=2, |
| 152 | + position_type='rdd', dtype=dtype) |
| 153 | + return result |
| 154 | + |
| 155 | +ranwt1=False |
| 156 | + |
| 157 | +regl = ['_N','_S'] |
| 158 | + |
| 159 | +tcorr = ttype |
| 160 | +tw = ttype |
| 161 | +if survey == 'main': |
| 162 | + regl = ['_DN','_DS','_N','_S'] |
| 163 | + if ttype == 'LRGrec': |
| 164 | + regl = ['_DN','_N'] |
| 165 | + tcorr = 'LRG' |
| 166 | + tw = 'LRG'+args.rectype+args.convention |
| 167 | + if ttype == 'ELGrec': |
| 168 | + regl = ['_DN','_N'] |
| 169 | + tcorr = 'ELG' |
| 170 | + tw = 'ELG'+args.rectype+args.convention |
| 171 | + |
| 172 | + |
| 173 | +nzr = len(zl) |
| 174 | +if len(zl) == 2: |
| 175 | + nzr = len(zl)-1 |
| 176 | +for i in range(0,nzr): |
| 177 | + if i == len(zl)-1: |
| 178 | + zmin=zl[0] |
| 179 | + zmax=zl[-1] |
| 180 | + else: |
| 181 | + zmin = zl[i] |
| 182 | + zmax = zl[i+1] |
| 183 | + print(zmin,zmax) |
| 184 | + for reg in regl: |
| 185 | + print(reg) |
| 186 | + result = compute_power_spectrum(edges, tracer=tcorr, region=reg, nrandoms=args.nran, zlim=(zmin,zmax), weight_type=weight_type) |
| 187 | + poles = result.poles |
| 188 | + fo = open(dirpk+'pk024'+tw+survey+reg+'_'+str(zmin)+str(zmax)+version+'_'+weight_type+'lin.dat','w') |
| 189 | + fo.write('#norm = {}\n'.format(poles.wnorm)) |
| 190 | + fo.write('#shotnoise = {}\n'.format(poles.shotnoise)) |
| 191 | + for k,p in zip(poles.k, poles.power.T.real): |
| 192 | + fo.write(' '.join([str(k)] + [str(p_) for p_ in p]) + '\n') |
| 193 | + fo.close() |
0 commit comments