/
mtl.py
990 lines (840 loc) · 37.6 KB
/
mtl.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
"""
desitarget.mtl
==============
Merged target lists.
"""
import os
import numpy as np
import healpy as hp
import numpy.lib.recfunctions as rfn
import sys
from astropy.table import Table
from astropy.io import ascii
import fitsio
from time import time
from datetime import datetime
from glob import glob
from . import __version__ as dt_version
from desitarget.targetmask import obsmask, obsconditions
from desitarget.targets import calc_priority, calc_numobs_more
from desitarget.targets import main_cmx_or_sv, switch_main_cmx_or_sv
from desitarget.targets import set_obsconditions
from desitarget.geomask import match, match_to
from desitarget.internal import sharedmem
from desitarget import io
# ADM set up the DESI default logger.
from desiutil.log import get_logger
log = get_logger()
# ADM the data model for MTL. Note that the _TARGET columns will have
# ADM to be changed on the fly for SV1_, SV2_, etc. files.
# ADM OBSCONDITIONS is formatted as just 'i4' for backward compatibility.
mtldatamodel = np.array([], dtype=[
('RA', '>f8'), ('DEC', '>f8'), ('PARALLAX', '>f4'),
('PMRA', '>f4'), ('PMDEC', '>f4'), ('REF_EPOCH', '>f4'),
('DESI_TARGET', '>i8'), ('BGS_TARGET', '>i8'), ('MWS_TARGET', '>i8'),
('SCND_TARGET', '>i8'), ('TARGETID', '>i8'),
('SUBPRIORITY', '>f8'), ('OBSCONDITIONS', 'i4'),
('PRIORITY_INIT', '>i8'), ('NUMOBS_INIT', '>i8'), ('PRIORITY', '>i8'),
('NUMOBS', '>i8'), ('NUMOBS_MORE', '>i8'), ('Z', '>f8'), ('ZWARN', '>i8'),
('TIMESTAMP', 'S19'), ('VERSION', 'S14'), ('TARGET_STATE', 'S16'),
('ZTILEID', '>i4')
])
zcatdatamodel = np.array([], dtype=[
('RA', '>f8'), ('DEC', '>f8'), ('TARGETID', '>i8'),
('NUMOBS', '>i4'), ('Z', '>f8'), ('ZWARN', '>i8'), ('ZTILEID', '>i4')
])
mtltilefiledm = np.array([], dtype=[
('TILEID', '>i4'), ('TIMESTAMP', 'S19'),
('VERSION', 'S14'), ('PROGRAM', 'S6')
])
# ADM when using basic or csv ascii writes, specifying the formats of
# ADM float32 columns can make things easier on the eye.
mtlformatdict = {"PARALLAX": '%16.8f', 'PMRA': '%16.8f', 'PMDEC': '%16.8f'}
def get_utc_date():
"""Convenience function to grab the UTC date.
Returns
-------
:class:`str`
The UTC data, appropriate to make a TIMESTAMP.
Notes
-----
- This is spun off into its own function to have a consistent way to
record time across the entire desitarget package.
"""
return datetime.utcnow().isoformat(timespec='seconds')
def get_mtl_dir(mtldir=None):
"""Convenience function to grab the $MTL_DIR environment variable.
Parameters
----------
mtldir : :class:`str`, optional, defaults to $MTL_DIR
If `mtldir` is passed, it is returned from this function. If it's
not passed, the $MTL_DIR environment variable is returned.
Returns
-------
:class:`str`
If `mtldir` is passed, it is returned from this function. If it's
not passed, the directory stored in the $MTL_DIR environment
variable is returned.
"""
if mtldir is None:
mtldir = os.environ.get('MTL_DIR')
# ADM check that the $MTL_DIR environment variable is set.
if mtldir is None:
msg = "Pass mtldir or set $MTL_DIR environment variable!"
log.critical(msg)
raise ValueError(msg)
return mtldir
def get_zcat_dir(zcatdir=None):
"""Convenience function to grab the $ZCAT_DIR environment variable.
Parameters
----------
zcatdir : :class:`str`, optional, defaults to $ZCAT_DIR
If `zcatdir` is passed, it is returned from this function. If it's
not passed, the $ZCAT_DIR environment variable is returned.
Returns
-------
:class:`str`
If `zcatdir` is passed, it is returned from this function. If it's
not passed, the directory stored in the $ZCAT_DIR environment
variable is returned.
"""
if zcatdir is None:
zcatdir = os.environ.get('ZCAT_DIR')
# ADM check that the $ZCAT_DIR environment variable is set.
if zcatdir is None:
msg = "Pass zcatdir or set $ZCAT_DIR environment variable!"
log.critical(msg)
raise ValueError(msg)
return zcatdir
def get_tile_full_path(tilefn=None):
"""Convenience function to grab the $TILE_FN environment variable.
Parameters
----------
tilefn : :class:`str`, optional, defaults to $TILE_FN
If `tilefn` is passed, it is returned from this function. If it's
not passed, the $TILE_FN environment variable is returned.
Returns
-------
:class:`str`
If `tilefn` is passed, it is returned from this function. If it's
not passed, the filename stored in the $TILE_FN environment
variable is returned.
"""
if tilefn is None:
tilefn = os.environ.get('TILE_FN')
# ADM check that the $TILE_FN environment variable is set.
if tilefn is None:
msg = "Pass tilefn or set $TILE_FN environment variable!"
log.critical(msg)
raise ValueError(msg)
return tilefn
def get_mtl_tile_file_name():
"""Convenience function to grab the name of the MTL tile file.
Returns
-------
:class:`str`
The name of the MTL tile file.
"""
fn = "mtl-done-tiles.ecsv"
return fn
def _get_mtl_nside():
"""Grab the HEALPixel nside to be used for MTL ledger files.
Returns
-------
:class:`int`
The HEALPixel nside number for MTL file creation and retrieval.
"""
# from desitarget.geomask import pixarea2nside
# ADM the nside (16) appropriate to a 7 sq. deg. field.
# nside = pixarea2nside(7)
nside = 32
return nside
def get_mtl_ledger_format():
"""Grab the file format for MTL ledger files.
Returns
-------
:class:`str`
The file format for MTL ledgers. Should be "ecsv" or "fits".
"""
# ff = "fits"
ff = "ecsv"
return ff
def make_mtl(targets, obscon, zcat=None, scnd=None,
trim=False, trimcols=False, trimtozcat=False):
"""Adds fiberassign and zcat columns to a targets table.
Parameters
----------
targets : :class:`~numpy.array` or `~astropy.table.Table`
A numpy rec array or astropy Table with at least the columns
``TARGETID``, ``DESI_TARGET``, ``NUMOBS_INIT``, ``PRIORITY_INIT``.
or the corresponding columns for SV or commissioning.
obscon : :class:`str`
A combination of strings that are in the desitarget bitmask yaml
file (specifically in `desitarget.targetmask.obsconditions`), e.g.
"DARK|GRAY". Governs the behavior of how priorities are set based
on "obsconditions" in the desitarget bitmask yaml file.
zcat : :class:`~astropy.table.Table`, optional
Redshift catalog table with columns ``TARGETID``, ``NUMOBS``,
``Z``, ``ZWARN``, ``ZTILEID``.
scnd : :class:`~numpy.array`, `~astropy.table.Table`, optional
A set of secondary targets associated with the `targets`. As with
the `target` must include at least ``TARGETID``, ``NUMOBS_INIT``,
``PRIORITY_INIT`` or the corresponding SV columns.
The secondary targets will be padded to have the same columns
as the targets, and concatenated with them.
trim : :class:`bool`, optional
If ``True`` (default), don't include targets that don't need
any more observations. If ``False``, include every input target.
trimcols : :class:`bool`, optional, defaults to ``False``
Only pass through columns in `targets` that are actually needed
for fiberassign (see `desitarget.mtl.mtldatamodel`).
trimtozcat : :class:`bool`, optional, defaults to ``False``
Only return targets that have been UPDATED (i.e. the targets with
a match in `zcat`). Returns all targets if `zcat` is ``None``.
Returns
-------
:class:`~astropy.table.Table`
MTL Table with targets columns plus:
* NUMOBS_MORE - number of additional observations requested
* PRIORITY - target priority (larger number = higher priority)
* TARGET_STATE - the observing state that corresponds to PRIORITY
* OBSCONDITIONS - replaces old GRAYLAYER
* TIMESTAMP - time that (this) make_mtl() function was run
* VERSION - version of desitarget used to run make_mtl()
"""
start = time()
# ADM set up the default logger.
from desiutil.log import get_logger
log = get_logger()
# ADM if trimcols was passed, reduce input target columns to minimal.
if trimcols:
mtldm = switch_main_cmx_or_sv(mtldatamodel, targets)
cullcols = list(set(targets.dtype.names) - set(mtldm.dtype.names))
if isinstance(targets, Table):
targets.remove_columns(cullcols)
else:
targets = rfn.drop_fields(targets, cullcols)
# ADM determine whether the input targets are main survey, cmx or SV.
colnames, masks, survey = main_cmx_or_sv(targets, scnd=True)
# ADM set the first column to be the "desitarget" column
desi_target, desi_mask = colnames[0], masks[0]
scnd_target = colnames[-1]
# ADM if secondaries were passed, concatenate them with the targets.
if scnd is not None:
nrows = len(scnd)
log.info('Pad {} primary targets with {} secondaries...t={:.1f}s'.format(
len(targets), nrows, time()-start))
padit = np.zeros(nrows, dtype=targets.dtype)
sharedcols = set(targets.dtype.names).intersection(set(scnd.dtype.names))
for col in sharedcols:
padit[col] = scnd[col]
targets = np.concatenate([targets, padit])
# APC Propagate a flag on which targets came from scnd
is_scnd = np.repeat(False, len(targets))
is_scnd[-nrows:] = True
log.info('Done with padding...t={:.1f}s'.format(time()-start))
# Trim targets from zcat that aren't in original targets table.
if zcat is not None:
ok = np.in1d(zcat['TARGETID'], targets['TARGETID'])
num_extra = np.count_nonzero(~ok)
if num_extra > 0:
log.warning("Ignoring {} zcat entries that aren't "
"in the input target list".format(num_extra))
zcat = zcat[ok]
n = len(targets)
# ADM if a redshift catalog was passed, order it to match the input targets
# ADM catalog on 'TARGETID'.
if zcat is not None:
# ADM find where zcat matches target array.
zmatcher = match_to(targets["TARGETID"], zcat["TARGETID"])
ztargets = zcat
if ztargets.masked:
unobs = ztargets['NUMOBS'].mask
ztargets['NUMOBS'][unobs] = 0
unobsz = ztargets['Z'].mask
ztargets['Z'][unobsz] = -1
unobszw = ztargets['ZWARN'].mask
ztargets['ZWARN'][unobszw] = -1
else:
ztargets = Table()
ztargets['TARGETID'] = targets['TARGETID']
ztargets['NUMOBS'] = np.zeros(n, dtype=np.int32)
ztargets['Z'] = -1 * np.ones(n, dtype=np.float32)
ztargets['ZWARN'] = -1 * np.ones(n, dtype=np.int32)
ztargets['ZTILEID'] = -1 * np.ones(n, dtype=np.int32)
# ADM if zcat wasn't passed, there is a one-to-one correspondence
# ADM between the targets and the zcat.
zmatcher = np.arange(n)
# ADM extract just the targets that match the input zcat.
targets_zmatcher = targets[zmatcher]
# ADM update the number of observations for the targets.
ztargets['NUMOBS_MORE'] = calc_numobs_more(targets_zmatcher, ztargets, obscon)
# ADM assign priorities. Only things in the zcat can have changed
# ADM priorities. Anything else is assigned PRIORITY_INIT, below.
priority, target_state = calc_priority(
targets_zmatcher, ztargets, obscon, state=True)
# If priority went to 0==DONOTOBSERVE or 1==OBS or 2==DONE, then
# NUMOBS_MORE should also be 0.
# ## mtl['NUMOBS_MORE'] = ztargets['NUMOBS_MORE']
ii = (priority <= 2)
log.info('{:d} of {:d} targets have priority zero, setting N_obs=0.'.format(
np.sum(ii), n))
ztargets['NUMOBS_MORE'][ii] = 0
# - Set the OBSCONDITIONS mask for each target bit.
obsconmask = set_obsconditions(targets)
# APC obsconmask will now be incorrect for secondary-only targets. Fix this
# APC using the mask on secondary targets.
if scnd is not None:
obsconmask[is_scnd] = set_obsconditions(targets[is_scnd], scnd=True)
# ADM set up the output mtl table.
mtl = Table(targets)
mtl.meta['EXTNAME'] = 'MTL'
# ADM add a placeholder for the secondary bit-mask, if it isn't there.
if scnd_target not in mtl.dtype.names:
mtl[scnd_target] = np.zeros(len(mtl),
dtype=mtldatamodel["SCND_TARGET"].dtype)
# ADM initialize columns to avoid zero-length/missing/format errors.
zcols = ["NUMOBS_MORE", "NUMOBS", "Z", "ZWARN"]
for col in zcols + ["TARGET_STATE", "TIMESTAMP", "VERSION"]:
mtl[col] = np.empty(len(mtl), dtype=mtldatamodel[col].dtype)
# ADM any target that wasn't matched to the ZCAT should retain its
# ADM original (INIT) value of PRIORITY and NUMOBS.
mtl['NUMOBS_MORE'] = mtl['NUMOBS_INIT']
mtl['PRIORITY'] = mtl['PRIORITY_INIT']
mtl['TARGET_STATE'] = "UNOBS"
# ADM add the time and version of the desitarget code that was run.
mtl["TIMESTAMP"] = get_utc_date()
mtl["VERSION"] = dt_version
# ADM now populate the new mtl columns with the updated information.
mtl['OBSCONDITIONS'] = obsconmask
mtl['PRIORITY'][zmatcher] = priority
mtl['TARGET_STATE'][zmatcher] = target_state
for col in zcols:
mtl[col][zmatcher] = ztargets[col]
# ADM also add the ZTILEID column, if passed, otherwise we're likely
# ADM to be working with non-ledger-based mocks and can let it slide.
if "ZTILEID" in ztargets:
mtl["ZTILEID"][zmatcher] = ztargets["ZTILEID"]
else:
mtl["ZTILEID"] = -1
# Filter out any targets marked as done.
if trim:
notdone = mtl['NUMOBS_MORE'] > 0
log.info('{:d} of {:d} targets are done, trimming these'.format(
len(mtl) - np.sum(notdone), len(mtl))
)
mtl = mtl[notdone]
# Filtering can reset the fill_value, which is just wrong wrong wrong
# See https://github.com/astropy/astropy/issues/4707
# and https://github.com/astropy/astropy/issues/4708
mtl['NUMOBS_MORE'].fill_value = -1
# ADM assert the data model is complete.
# ADM turning this off for now, useful for testing.
# mtltypes = [mtl[i].dtype.type for i in mtl.dtype.names]
# mtldmtypes = [mtldm[i].dtype.type for i in mtl.dtype.names]
# assert set(mtl.dtype.names) == set(mtldm.dtype.names)
# assert mtltypes == mtldmtypes
log.info('Done...t={:.1f}s'.format(time()-start))
if trimtozcat:
return mtl[zmatcher]
return mtl
def make_ledger_in_hp(targets, outdirname, nside, pixlist,
obscon="DARK", indirname=None, verbose=True):
"""
Make an initial MTL ledger file for targets in a set of HEALPixels.
Parameters
----------
targets : :class:`~numpy.array`
Targets made by, e.g. `desitarget.cuts.select_targets()`.
outdirname : :class:`str`
Output directory to which to write the MTLs (the file names are
constructed on the fly).
nside : :class:`int`
(NESTED) HEALPixel nside that corresponds to `pixlist`.
pixlist : :class:`list` or `int`
HEALPixels at `nside` at which to write the MTLs.
obscon : :class:`str`, optional, defaults to "DARK"
A string matching ONE obscondition in the desitarget bitmask yaml
file (i.e. in `desitarget.targetmask.obsconditions`), e.g. "GRAY"
Governs how priorities are set when merging targets. Also governs
the sub-directory to which the ledger is written.
indirname : :class:`str`
A directory associated with the targets. Written to the headers
of the output MTL files.
verbose : :class:`bool`, optional, defaults to ``True``
If ``True`` then log target and file information.
Returns
-------
Nothing, but writes the `targets` out to `outdirname` split across
each HEALPixel in `pixlist`.
"""
t0 = time()
# ADM in case an integer was passed.
pixlist = np.atleast_1d(pixlist)
# ADM execute MTL.
mtl = make_mtl(targets, obscon, trimcols=True)
# ADM the HEALPixel within which each target in the MTL lies.
theta, phi = np.radians(90-mtl["DEC"]), np.radians(mtl["RA"])
mtlpix = hp.ang2pix(nside, theta, phi, nest=True)
# ADM write the MTLs.
_, _, survey = main_cmx_or_sv(mtl)
for pix in pixlist:
inpix = mtlpix == pix
ecsv = get_mtl_ledger_format() == "ecsv"
nt, fn = io.write_mtl(
outdirname, mtl[inpix].as_array(), indir=indirname, ecsv=ecsv,
survey=survey, obscon=obscon, nsidefile=nside, hpxlist=pix)
if verbose:
log.info('{} targets written to {}...t={:.1f}s'.format(
nt, fn, time()-t0))
return
def make_ledger(hpdirname, outdirname, pixlist=None, obscon="DARK", numproc=1):
"""
Make initial MTL ledger files for HEALPixels, in parallel.
Parameters
----------
hpdirname : :class:`str`
Full path to either a directory containing targets that
have been partitioned by HEALPixel (i.e. as made by
`select_targets` with the `bundle_files` option). Or the
name of a single file of targets.
outdirname : :class:`str`
Output directory to which to write the MTL (the file name is
constructed on the fly).
pixlist : :class:`list` or `int`, defaults to ``None``
(Nested) HEALPixels at which to write the MTLs at the default
`nside` (which is `_get_mtl_nside()`). Defaults to ``None``,
which runs all of the pixels at `_get_mtl_nside()`.
obscon : :class:`str`, optional, defaults to "DARK"
A string matching ONE obscondition in the desitarget bitmask yaml
file (i.e. in `desitarget.targetmask.obsconditions`), e.g. "GRAY"
Governs how priorities are set based on "obsconditions". Also
governs the sub-directory to which the ledger is written.
numproc : :class:`int`, optional, defaults to 1 for serial
Number of processes to parallelize across.
Returns
-------
Nothing, but writes the full HEALPixel-split ledger to `outdirname`.
Notes
-----
- For _get_mtl_nside()=32, takes about 25 minutes with `numproc=12`.
`numproc>12` can run into memory issues.
- For _get_mtl_nside()=16, takes about 50 minutes with `numproc=8`.
`numproc>8` can run into memory issues.
"""
# ADM grab information regarding how the targets were constructed.
hdr, dt = io.read_targets_header(hpdirname, dtype=True)
# ADM check the obscon for which the targets were made is
# ADM consistent with the requested obscon.
oc = hdr["OBSCON"]
if obscon not in oc:
msg = "File is type {} but requested behavior is {}".format(oc, obscon)
log.critical(msg)
raise ValueError(msg)
# ADM the MTL datamodel must reflect the target flavor (SV, etc.).
mtldm = switch_main_cmx_or_sv(mtldatamodel, np.array([], dt))
# ADM speed-up by only reading the necessary columns.
cols = list(set(mtldm.dtype.names).intersection(dt.names))
# ADM optimal nside for reading in the targeting files.
nside = hdr["FILENSID"]
npixels = hp.nside2npix(nside)
pixels = np.arange(npixels)
# ADM the nside at which to write the MTLs.
mtlnside = _get_mtl_nside()
# ADM default to running all pixels.
if pixlist is None:
mtlnpixels = hp.nside2npix(mtlnside)
pixlist = np.arange(mtlnpixels)
# ADM check that the nside for writing MTLs is not at a lower
# ADM resolution than the nside at which the files are stored.
msg = "Ledger nside ({}) must be higher than file nside ({})!!!".format(
mtlnside, nside)
assert mtlnside >= nside, msg
from desitarget.geomask import nside2nside
# ADM the common function that is actually parallelized across.
def _make_ledger_in_hp(pixnum):
"""make initial ledger in a single HEALPixel"""
# ADM construct a list of all pixels in pixnum at the MTL nside.
setpix = set(nside2nside(nside, mtlnside, pixnum))
pix = [p for p in pixlist if p in setpix]
if len(pix) == 0:
return
# ADM read in the needed columns from the targets.
targs = io.read_targets_in_hp(hpdirname, nside, pixnum, columns=cols)
if len(targs) == 0:
return
# ADM write MTLs for the targs split over HEALPixels in pixlist.
return make_ledger_in_hp(
targs, outdirname, mtlnside, pix,
obscon=obscon, indirname=hpdirname, verbose=False)
# ADM this is just to count pixels in _update_status.
npix = np.ones((), dtype='i8')
t0 = time()
def _update_status(result):
"""wrap key reduction operation on the main parallel process"""
if npix % 2 == 0 and npix > 0:
rate = (time() - t0) / npix
log.info('{}/{} HEALPixels; {:.1f} secs/pixel...t = {:.1f} mins'.
format(npix, npixels, rate, (time()-t0)/60.))
npix[...] += 1
return result
# ADM Parallel process across HEALPixels.
if numproc > 1:
pool = sharedmem.MapReduce(np=numproc)
with pool:
pool.map(_make_ledger_in_hp, pixels, reduce=_update_status)
else:
for pixel in pixels:
_update_status(_make_ledger_in_hp(pixel))
log.info("Done writing ledger to {}...t = {:.1f} mins".format(
outdirname, (time()-t0)/60.))
return
def update_ledger(hpdirname, zcat, targets=None, obscon="DARK",
numobs_from_ledger=False):
"""
Update relevant HEALPixel-split ledger files for some targets.
Parameters
----------
hpdirname : :class:`str`
Full path to a directory containing an MTL ledger that has been
partitioned by HEALPixel (i.e. as made by `make_ledger`).
zcat : :class:`~astropy.table.Table`, optional
Redshift catalog table with columns ``TARGETID``, ``NUMOBS``,
``Z``, ``ZWARN``, ``ZTILEID``.
targets : :class:`~numpy.array` or `~astropy.table.Table`, optional, defaults to ``None``
A numpy rec array or astropy Table with at least the columns
``RA``, ``DEC``, ``TARGETID``, ``DESI_TARGET``, ``NUMOBS_INIT``,
and ``PRIORITY_INIT``. If ``None``, then assume the `zcat`
includes ``RA`` and ``DEC`` and look up `targets` in the ledger.
obscon : :class:`str`, optional, defaults to "DARK"
A string matching ONE obscondition in the desitarget bitmask yaml
file (i.e. in `desitarget.targetmask.obsconditions`), e.g. "GRAY"
Governs how priorities are set using "obsconditions". Basically a
check on whether the files in `hpdirname` are as expected.
numobs_from_ledger : :class:`bool`, optional, defaults to ``True``
If ``True`` then inherit the number of observations so far from
the ledger rather than expecting it to have a reasonable value
in the `zcat.`
Returns
-------
Nothing, but relevant ledger files are updated.
"""
# ADM find the general format for the ledger files in `hpdirname`.
# ADM also returning the obsconditions.
fileform, oc = io.find_mtl_file_format_from_header(hpdirname, returnoc=True)
# ADM check the obscondition is as expected.
if obscon != oc:
msg = "File is type {} but requested behavior is {}".format(oc, obscon)
log.critical(msg)
raise ValueError(msg)
# ADM if targets wasn't sent, that means the zcat includes
# ADM coordinates and we can read relevant targets from the ledger.
if targets is None:
nside = _get_mtl_nside()
theta, phi = np.radians(90-zcat["DEC"]), np.radians(zcat["RA"])
pixnum = hp.ang2pix(nside, theta, phi, nest=True)
# ADM we'll read in too many targets, here, but that's OK as
# ADM make_mtl(trimtozcat=True) only returns the updated targets.
targets, fndict = io.read_mtl_in_hp(hpdirname, nside, pixnum,
unique=True, returnfn=True)
# ADM if requested, use the previous values in the ledger to set
# ADM NUMOBS in the zcat.
if numobs_from_ledger:
# ADM match the zcat to the targets.
tii, zii = match(targets["TARGETID"], zcat["TARGETID"])
# ADM update NUMOBS in the zcat for matches.
zcat["NUMOBS"][zii] = targets["NUMOBS"][tii] + 1
# ADM run MTL, only returning the targets that are updated.
mtl = make_mtl(targets, oc, zcat=zcat, trimtozcat=True, trimcols=True)
# ADM this is redundant if targets wasn't sent, but it's quick.
nside = _get_mtl_nside()
theta, phi = np.radians(90-mtl["DEC"]), np.radians(mtl["RA"])
pixnum = hp.ang2pix(nside, theta, phi, nest=True)
# ADM loop through the pixels and update the ledger, depending
# ADM on whether we're working with .fits or .ecsv files.
ender = get_mtl_ledger_format()
for pix in set(pixnum):
# ADM grab the targets in the pixel.
ii = pixnum == pix
mtlpix = mtl[ii]
# ADM sorting on TARGETID is important for
# ADM io.read_mtl_ledger(unique=True)
mtlpix = mtlpix[np.argsort(mtlpix["TARGETID"])]
# ADM the correct filename for this pixel number.
fn = fileform.format(pix)
# ADM if we're working with .ecsv, simply append to the ledger.
if ender == 'ecsv':
f = open(fn, "a")
ascii.write(mtlpix, f, format='no_header', formats=mtlformatdict)
f.close()
# ADM otherwise, for FITS, we'll have to read in the whole file.
else:
ledger, hd = fitsio.read(fn, extname="MTL", header=True)
done = np.concatenate([ledger, mtlpix.as_array()])
fitsio.write(fn+'.tmp', done, extname='MTL', header=hd, clobber=True)
os.rename(fn+'.tmp', fn)
return
def inflate_ledger(mtl, hpdirname, columns=None, header=False, strictcols=False,
quick=False):
"""Add a fuller set of target columns to an MTL.
Parameters
----------
mtl : :class:`~numpy.array` or `~astropy.table.Table`
A Merged Target List in array or Table form. Must contain the
columns "RA", "DEC" and "TARGETID"
hpdirname : :class:`str`
Full path to a directory containing targets that have been
partitioned by HEALPixel (i.e. as made by `select_targets`
with the `bundle_files` option).
columns : :class:`list` or :class:`str`, optional
Only return these target columns. If ``None`` or not passed
then return all of the target columns.
header : :class:`bool`, optional, defaults to ``False``
If ``True`` then also return the header of the last file read
from the `hpdirname` directory.
strictcols : :class:`bool`, optional, defaults to ``False``
If ``True`` then strictly return only the columns in `columns`,
otherwise, inflate the ledger with the new columns.
quick : :class:`bool`, optional, defaults to ``False``
If ``True``, call the "quick" version of reading targets
which is much faster but makes fewer error checks.
Returns
-------
:class:`~numpy.array`
The original MTL with the fuller set of columns.
Notes
-----
- Will run more quickly if the targets in `mtl` are clustered.
- TARGETID is always returned, even if it isn't in `columns`.
"""
# ADM if a table was passed convert it to a numpy array.
if isinstance(mtl, Table):
mtl = mtl.as_array()
# ADM if a string was passed for the columns, convert it to a list.
if isinstance(columns, str):
columns = [columns]
# ADM we have to have TARGETID, even if it wasn't a passed column.
if columns is not None:
origcols = columns.copy()
if "TARGETID" not in columns:
columns.append("TARGETID")
# ADM look up the optimal nside for reading targets.
nside, _ = io.check_hp_target_dir(hpdirname)
# ADM which pixels do we need to read.
theta, phi = np.radians(90-mtl["DEC"]), np.radians(mtl["RA"])
pixnums = hp.ang2pix(nside, theta, phi, nest=True)
pixlist = list(set(pixnums))
# ADM read in targets in the required pixels.
targs = io.read_targets_in_hp(hpdirname, nside, pixlist, columns=columns,
header=header, quick=quick)
if header:
targs, hdr = targs
# ADM match the mtl back to the targets on TARGETID.
ii = match_to(targs["TARGETID"], mtl["TARGETID"])
# ADM extract just the targets that match the mtl.
targs = targs[ii]
# ADM create an array to contain the fuller set of target columns.
# ADM start with the data model for the target columns.
dt = targs.dtype.descr
# ADM add the unique columns from the mtl.
xtracols = [nam for nam in mtl.dtype.names if nam not in targs.dtype.names]
for col in xtracols:
dt.append((col, mtl[col].dtype.str))
# ADM remove columns from the data model that weren't requested.
if columns is not None and strictcols:
dt = [(name, form) for name, form in dt if name in origcols]
# ADM create the output array.
done = np.empty(len(mtl), dtype=dt)
# ADM populate the output array with the fuller target columns.
for col in targs.dtype.names:
if col in done.dtype.names:
done[col] = targs[col]
# ADM populate the output array with the unique MTL columns.
for col in xtracols:
if col in done.dtype.names:
done[col] = mtl[col]
if header:
return done, hdr
return done
def tiles_to_be_processed(zcatdir, mtltilefn, tilefn, obscon):
"""Execute full MTL loop, including reading files, updating ledgers.
Parameters
----------
zcatdir : :class:`str`
Full path to the directory that hosts the redshift catalogs.
mtltilefn : :class:`str`
Full path to the file of tiles that have been processed by MTL.
tilefn : :class:`str`, optional, defaults to ``None``
Full path to the name of the tile file. This file is used to link
TILEIDs to observing conditions.
obscon : :class:`str`
A string matching ONE obscondition in the desitarget bitmask yaml
file (i.e. in `desitarget.targetmask.obsconditions`), e.g. "DARK"
Governs how priorities are set when merging targets.
Returns
-------
:class:`~numpy.array`
An array of tiles that have not yet been processed and written to
the mtl tile file.
"""
# ADM read in the tile file.
tilelookup = Table.read(tilefn)
# ADM determine the tiles to be processed.
alltiles = []
for fn in glob(os.path.join(zcatdir, "*")):
# ADM sometimes there are weird tile names like "temp".
try:
tileid = int(os.path.basename(fn))
except ValueError:
pass
alltiles.append(tileid)
# ADM read in the tile file, guarding against it not having being
# ADM created yet.
donetiles = None
if os.path.isfile(mtltilefn):
donetiles = io.read_mtl_tile_file(mtltilefn)
# ADM extract the updated tiles.
if donetiles is None:
tilesallcon = np.array(alltiles)
else:
tilesallcon = np.array(list(set(alltiles) - set(donetiles["TILEID"])))
obstiles = tilelookup[tilelookup["PROGRAM"] == obscon]["TILEID"]
tileids = list(set(obstiles).intersection(set(tilesallcon)))
# ADM initialize the output array and add the tiles.
donetiles = np.zeros(len(tileids), dtype=mtltilefiledm.dtype)
donetiles["TILEID"] = tileids
# ADM look up the time.
donetiles["TIMESTAMP"] = get_utc_date()
# ADM add the version of desitarget.
donetiles["VERSION"] = dt_version
# ADM add the program/obscon.
donetiles["PROGRAM"] = obscon
return donetiles
def make_zcat_rr_backstop(zcatdir, tiles):
"""Make the simplest possible zcat using SV1-era redrock outputs.
Parameters
----------
zcatdir : :class:`str`
Full path to the directory that hosts the redshift catalogs.
tiles : :class:`~numpy.array`
List of TILEIDs in `zcatdir` from which to construct the `zcat`.
Returns
-------
:class:`~astropy.table.Table`
A zcat in the official format (`zcatdatamodel`) compiled from
the `tiles` in `zcatdir`.
Notes
-----
- How the `zcat` is constructed could certainly change once we have
the final schema in place.
"""
# ADM for each tile, read in the spectroscopic and targeting info.
allzs = []
allfms = []
for tile in tiles:
zbestfns = os.path.join(zcatdir, "{}".format(tile), "*", "zbest*")
for zbestfn in glob(zbestfns):
allzs.append(fitsio.read(zbestfn, "ZBEST"))
fm = fitsio.read(zbestfn, "FIBERMAP")
allfms.append(fm)
# ADM check the correct TILEID was written in the fibermap.
if set(fm["TILEID"]) != set([tile]):
msg = "Directory and fibermap don't match for tile".format(tile)
log.critical(msg)
raise ValueError(msg)
zs = np.concatenate(allzs)
fms = np.concatenate(allfms)
# ADM currently, the spectroscopic files aren't coadds, so aren't
# ADM unique. We therefore need to look up (any) coordinates for
# ADM each z in the fibermap.
zid = match_to(fms["TARGETID"], zs["TARGETID"])
# ADM write out the zcat as a file with the correct data model.
zcat = Table(np.zeros(len(zs), dtype=zcatdatamodel.dtype))
zcat["RA"] = fms[zid]["TARGET_RA"]
zcat["DEC"] = fms[zid]["TARGET_DEC"]
zcat["ZTILEID"] = fms[zid]["TILEID"]
zcat["NUMOBS"] = zs["NUMTILE"]
for col in set(zcat.dtype.names) - set(['RA', 'DEC', 'NUMOBS', 'ZTILEID']):
zcat[col] = zs[col]
return zcat
def loop_ledger(obscon, survey='main', zcatdir=None, mtldir=None,
tilefn=None, numobs_from_ledger=True):
"""Execute full MTL loop, including reading files, updating ledgers.
Parameters
----------
obscon : :class:`str`
A string matching ONE obscondition in the desitarget bitmask yaml
file (i.e. in `desitarget.targetmask.obsconditions`), e.g. "DARK"
Governs how priorities are set when merging targets.
survey : :class:`str`, optional, defaults to "main"
Used to look up the correct filename for `obscon`. Options are
``'main'`` and ``'svX``' (where X is 1, 2, 3 etc.) for the main
survey and different iterations of SV, respectively.
zcatdir : :class:`str`, optional, defaults to ``None``
Full path to the directory that hosts the redshift catalogs.
If this is ``None``, look up the redshift catalog directory from
the $ZCAT_DIR environment variable.
mtldir : :class:`str`, optional, defaults to ``None``
Full path to the directory that hosts the MTL ledgers and the MTL
tile file. If this ``None``, look up the MTL directory from the
$MTL_DIR environment variable.
tilefn : :class:`str`, optional, defaults to ``None``
Full path to the name of the tile file. This file is used to link
TILEIDs to observing conditions. If passed as ``None``, looked up
from the $TILE_FN environment variable.
numobs_from_ledger : :class:`bool`, optional, defaults to ``True``
If ``True`` then inherit the number of observations so far from
the ledger rather than expecting it to have a reasonable value
in the `zcat.`
Returns
-------
:class:`str`
The directory containing the ledger that was updated.
:class:`str`
The name of the MTL tile file that was updated.
:class:`str`
The name of the tile file that was used to link TILEIDs to
observing conditions.
:class:`~numpy.array`
Information for the tiles that were processed.
Notes
-----
- Assumes all of the relevant ledgers have already been made by,
e.g., `make_ledger()`.
"""
# ADM first grab all of the relevant files.
# ADM grab the MTL directory (in case we're relying on $MTL_DIR).
mtldir = get_mtl_dir(mtldir)
# ADM construct the full path to the mtl tile file.
mtltilefn = os.path.join(mtldir, get_mtl_tile_file_name())
# ADM grab the TILE filename (in case we're relying on $TILE_FN).
tilefn = get_tile_full_path(tilefn)
# ADM construct the relevant sub-directory for this survey and
# ADM set of observing conditions..
form = get_mtl_ledger_format()
hpdirname = io.find_target_files(mtldir, flavor="mtl",
survey=survey, obscon=obscon, ender=form)
# ADM grab the zcat directory (in case we're relying on $ZCAT_DIR).
zcatdir = get_zcat_dir(zcatdir)
# ADM grab an array of tiles that are yet to be processed.
tiles = tiles_to_be_processed(zcatdir, mtltilefn, tilefn, obscon)
# ADM stop if there are no tiles to process.
if len(tiles) == 0:
return hpdirname, mtltilefn, tilefn, tiles
# ADM create the zcat: This will likely change, but for now let's
# ADM just use redrock in the SV1-era format.
zcat = make_zcat_rr_backstop(zcatdir, tiles["TILEID"])
# ADM insist that for an MTL loop with real observations, the zcat
# ADM must conform to the data model. In particular, it must include
# ADM ZTILEID, which may not be needed for non-ledger simulations.
if zcat.dtype.descr != zcatdatamodel.dtype.descr:
msg = "zcat data model must be {} not {}!".format(
zcatdatamodel.dtype.descr, zcat.dtype.descr)
log.critical(msg)
raise ValueError(msg)
# ADM update the appropriate ledger.
update_ledger(hpdirname, zcat, obscon=obscon,
numobs_from_ledger=numobs_from_ledger)
# ADM write the processed tiles to the MTL tile file.
io.write_mtl_tile_file(mtltilefn, tiles)
return hpdirname, mtltilefn, tilefn, tiles