• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

igmhub / picca / 25863215785

14 May 2026 01:38PM UTC coverage: 62.644% (-0.06%) from 62.708%
25863215785

Pull #1138

gihtub

iprafols
added safety "/" in case in_dir does not finish with it
Pull Request #1138: Fix Issue 1137

1388 of 2827 branches covered (49.1%)

Branch coverage included in aggregate %.

26 of 58 new or added lines in 14 files covered. (44.83%)

411 existing lines in 4 files now uncovered.

6544 of 9835 relevant lines covered (66.54%)

2.0 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

51.62
/py/picca/data.py
1
"""This module defines data structure to deal with line of sight data.
2

3
This module provides with three classes (QSO, Forest, Delta)
4
to manage the line-of-sight data.
5
See the respective docstrings for more details
6
"""
7
import numpy as np
3✔
8
from itertools import repeat
3✔
9
import warnings
3✔
10

11
from . import constants
3✔
12
from .utils import userprint
3✔
13

14
class QSO(object):
3✔
15
    """Class to represent quasar objects.
16

17
    Attributes:
18
        ra: float
19
            Right-ascension of the quasar (in radians).
20
        dec: float
21
            Declination of the quasar (in radians).
22
        z_qso: float
23
            Redshift of the quasar.
24
        plate: integer
25
            Plate number of the observation.
26
        fiberid: integer
27
            Fiberid of the observation.
28
        mjd: integer
29
            Modified Julian Date of the observation.
30
        thingid: integer
31
            Thingid of the observation.
32
        x_cart: float
33
            The x coordinate when representing ra, dec in a cartesian
34
            coordinate system.
35
        y_cart: float
36
            The y coordinate when representing ra, dec in a cartesian
37
            coordinate system.
38
        z_cart: float
39
            The z coordinate when representing ra, dec in a cartesian
40
            coordinate system.
41
        cos_dec: float
42
            Cosine of the declination angle.
43
        weights: float
44
            Weight assigned to object
45
        r_comov: float or None
46
            Comoving distance to the object
47
        dist_m: float or None
48
            Angular diameter distance to object
49
        log_lambda: float or None
50
            Wavelength associated with the quasar redshift
51

52
    Note that plate-fiberid-mjd is a unique identifier
53
    for the quasar.
54

55
    Methods:
56
        __init__: Initialize class instance.
57
        get_angle_between: Computes the angular separation between two quasars.
58
    """
59

60
    def __init__(self, los_id, ra, dec, z_qso, plate, mjd, fiberid):
3✔
61
        """Initializes class instance.
62

63
        Args:
64
            thingid: integer
65
                Thingid of the observation.
66
            ra: float
67
                Right-ascension of the quasar (in radians).
68
            dec: float
69
                Declination of the quasar (in radians).
70
            z_qso: float
71
                Redshift of the quasar.
72
            plate: integer
73
                Plate number of the observation.
74
            mjd: integer
75
                Modified Julian Date of the observation.
76
            fiberid: integer
77
                Fiberid of the observation.
78
        """
79
        self.ra = ra
3✔
80
        self.dec = dec
3✔
81

82
        self.plate = plate
3✔
83
        self.mjd = mjd
3✔
84
        self.fiberid = fiberid
3✔
85

86
        ## cartesian coordinates
87
        self.x_cart = np.cos(ra) * np.cos(dec)
3✔
88
        self.y_cart = np.sin(ra) * np.cos(dec)
3✔
89
        self.z_cart = np.sin(dec)
3✔
90
        self.cos_dec = np.cos(dec)
3✔
91

92
        self.z_qso = z_qso
3✔
93
        self.los_id = los_id
3✔
94
        #this is for legacy purposes only
95
        self.thingid = los_id
3✔
96
        warnings.warn("currently a thingid entry is created in QSO.__init__, this feature will be removed", DeprecationWarning)
3✔
97

98
        # variables computed in function io.read_objects
99
        self.weight = None
3✔
100
        self.r_comov = None
3✔
101
        self.dist_m = None
3✔
102

103
        # variables computed in modules bin.picca_xcf_angl and bin.picca_xcf1d
104
        self.log_lambda = None
3✔
105

106
    def get_angle_between(self, data):
3✔
107
        """Computes the angular separation between two quasars.
108

109
        Args:
110
            data: QSO or list of QSO
111
                Objects with which the angular separation will
112
                be computed.
113

114
        Returns
115
            A float or an array (depending on input data) with the angular
116
            separation between this quasar and the object(s) in data.
117
        """
118
        # case 1: data is list-like
119
        try:
3✔
120
            x_cart = np.array([d.x_cart for d in data])
3✔
121
            y_cart = np.array([d.y_cart for d in data])
3✔
122
            z_cart = np.array([d.z_cart for d in data])
3✔
123
            ra = np.array([d.ra for d in data])
3✔
124
            dec = np.array([d.dec for d in data])
3✔
125

126
            cos = x_cart * self.x_cart + y_cart * self.y_cart + z_cart * self.z_cart
3✔
127
            w = cos >= 1.
3✔
128
            if w.sum() != 0:
3!
129
                userprint('WARNING: {} pairs have cos>=1.'.format(w.sum()))
×
130
                cos[w] = 1.
×
131
            w = cos <= -1.
3✔
132
            if w.sum() != 0:
3!
133
                userprint('WARNING: {} pairs have cos<=-1.'.format(w.sum()))
×
134
                cos[w] = -1.
×
135
            angl = np.arccos(cos)
3✔
136

137
            w = ((np.absolute(ra - self.ra) < constants.SMALL_ANGLE_CUT_OFF) &
3✔
138
                 (np.absolute(dec - self.dec) < constants.SMALL_ANGLE_CUT_OFF))
139
            if w.sum() != 0:
3!
140
                angl[w] = np.sqrt((dec[w] - self.dec)**2 +
×
141
                                  (self.cos_dec * (ra[w] - self.ra))**2)
142
        # case 2: data is a QSO
143
        except TypeError:
3✔
144
            x_cart = data.x_cart
3✔
145
            y_cart = data.y_cart
3✔
146
            z_cart = data.z_cart
3✔
147
            ra = data.ra
3✔
148
            dec = data.dec
3✔
149

150
            cos = x_cart * self.x_cart + y_cart * self.y_cart + z_cart * self.z_cart
3✔
151
            if cos >= 1.:
3!
152
                userprint('WARNING: 1 pair has cosinus>=1.')
×
153
                cos = 1.
×
154
            elif cos <= -1.:
3!
155
                userprint('WARNING: 1 pair has cosinus<=-1.')
×
156
                cos = -1.
×
157
            angl = np.arccos(cos)
3✔
158
            if ((np.absolute(ra - self.ra) < constants.SMALL_ANGLE_CUT_OFF) &
3!
159
                    (np.absolute(dec - self.dec) < constants.SMALL_ANGLE_CUT_OFF)):
160
                angl = np.sqrt((dec - self.dec)**2 + (self.cos_dec *
×
161
                                                      (ra - self.ra))**2)
162
        return angl
3✔
163

164

165
class Forest(QSO):
3✔
166
    """Class to represent a Lyman alpha (or other absorption) forest
167

168
    This is not a proper forest class anymore, it's just the minimal set of variables to not make the corr_function, stacking etc fail
169
    """
170
    log_lambda_min = None
3✔
171
    log_lambda_max = None
3✔
172
    delta_log_lambda = None
3✔
173

174
    @classmethod
3✔
175
    def get_var_lss(cls, log_lambda):
3✔
176
        """Interpolates the pixel variance due to the Large Scale Strucure on
177
        the wavelength array.
178

179
        Empty function to be loaded at run-time.
180

181
        Args:
182
            log_lambda: array of float
183
                Array containing the logarithm of the wavelengths (in Angs)
184

185
        Returns:
186
            An array with the correction
187

188
        Raises:
189
            NotImplementedError: Function was not specified
190
        """
191
        raise NotImplementedError("Function should be specified at run-time")
192

193
    @classmethod
3✔
194
    def get_eta(cls, log_lambda):
3✔
195
        """Interpolates the correction factor to the contribution of the
196
        pipeline estimate of the instrumental noise to the variance on the
197
        wavelength array.
198

199
        See equation 4 of du Mas des Bourboux et al. 2020 for details.
200

201
        Empty function to be loaded at run-time.
202

203
        Args:
204
            log_lambda: array of float
205
                Array containing the logarithm of the wavelengths (in Angs)
206

207
        Returns:
208
            An array with the correction
209

210
        Raises:
211
            NotImplementedError: Function was not specified
212
        """
213
        raise NotImplementedError("Function should be specified at run-time")
214

215
    @classmethod
3✔
216
    def get_fudge(cls, log_lambda):
3✔
217
        """Interpolates the fudge contribution to the variance on the
218
        wavelength array.
219

220
        See function epsilon in equation 4 of du Mas des Bourboux et al.
221
        2020 for details.
222

223
        Args:
224
            log_lambda: array of float
225
                Array containing the logarithm of the wavelengths (in Angs)
226

227
        Returns:
228
            An array with the correction
229

230
        Raises:
231
            NotImplementedError: Function was not specified
232
        """
233
        raise NotImplementedError("Function should be specified at run-time")
234

235
   
236

237

238
class Delta(QSO):
3✔
239
    """Class to represent the mean transimission fluctuation field (delta)
240

241
    This class stores the information for the deltas for a given line of sight
242

243
    Attributes:
244
        ## Inherits from QSO ##
245
        log_lambda : array of floats
246
            Array containing the logarithm of the wavelengths (in Angs)
247
        weights : array of floats
248
            Weights associated to pixel. Overloaded from parent class
249
        cont: array of floats
250
            Quasar continuum
251
        delta: array of floats
252
            Mean transmission fluctuation (delta field)
253
        order: 0 or 1
254
            Order of the log10(lambda) polynomial for the continuum fit
255
        ivar: array of floats
256
            Inverse variance associated to each flux
257
        exposures_diff: array of floats
258
            Difference between exposures
259
        mean_snr: float
260
            Mean signal-to-noise ratio in the forest
261
        mean_reso: float
262
            Mean resolution of the forest in units of velocity (FWHM)
263
        mean_z: float
264
            Mean redshift of the forest
265
        mean_reso_pix: float
266
            Mean resolution of the forest in units of pixels (FWHM)
267
        mean_resolution_matrix: array of floats or None
268
            Mean (over wavelength) resolution matrix for that forest
269
        resolution_matrix: 2d array of floats or None
270
            Wavelength dependent resolution matrix for that forest
271
        delta_log_lambda: float
272
            Variation of the logarithm of the wavelength between two pixels
273
        z: array of floats or None
274
            Redshift of the abosrption
275
        r_comov: array of floats or None
276
            Comoving distance to the object. Overloaded from parent class
277
        dist_m: array of floats or None
278
            Angular diameter distance to object. Overloaded from parent
279
            class
280
        neighbours: list of Delta or QSO or None
281
            Neighbouring deltas/quasars
282
        fname: string or None
283
            String identifying Delta as part of a group
284

285
    Methods:
286
        __init__: Initializes class instances.
287
        from_fitsio: Initialize instance from a fits file.
288
        from_ascii: Initialize instance from an ascii file.
289
        from_image: Initialize instance from an ascii file.
290
        project: Project the delta field.
291

292
    """
293

294
    def __init__(self, los_id, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
3✔
295
                 weights, cont, delta, order, ivar, exposures_diff, mean_snr,
296
                 mean_reso, mean_z, resolution_matrix=None,
297
                 mean_resolution_matrix=None, mean_reso_pix=None):
298
        """Initializes class instances.
299

300
        Args:
301
            los_id: integer
302
                Thingid or Targetid of the observation.
303
            ra: float
304
                Right-ascension of the quasar (in radians).
305
            dec: float
306
                Declination of the quasar (in radians).
307
            z_qso: float
308
                Redshift of the quasar.
309
            plate: integer
310
                Plate number of the observation.
311
            mjd: integer
312
                Modified Julian Date of the observation.
313
            fiberid: integer
314
                Fiberid of the observation.
315
            log_lambda: array of floats
316
                Logarithm of the wavelengths (in Angs)
317
            weights: array of floats
318
                Pixel weights
319
            cont: array of floats
320
                Quasar continuum
321
            delta: array of floats
322
                Mean transmission fluctuation (delta field)
323
            order: 0 or 1
324
                Order of the log10(lambda) polynomial for the continuum fit
325
            ivar: array of floats
326
                Inverse variance associated to each flux
327
            exposures_diff: array of floats
328
                Difference between exposures
329
            mean_snr: float
330
                Mean signal-to-noise ratio in the forest
331
            mean_reso: float
332
                Mean resolution of the forest
333
            mean_z: float
334
                Mean redshift of the forest
335
            mean_reso_pix: float
336
                Mean resolution of the forest in units of pixels (FWHM)
337
            mean_resolution_matrix: array of floats or None
338
                Mean (over wavelength) resolution matrix for that forest
339
            resolution_matrix: 2d array of floats or None
340
                Wavelength dependent resolution matrix for that forest
341
            delta_log_lambda: float
342
                Variation of the logarithm of the wavelength between two pixels
343
        """
344
        QSO.__init__(self, los_id, ra, dec, z_qso, plate, mjd, fiberid)
3✔
345
        self.log_lambda = log_lambda
3✔
346
        self.weights = weights
3✔
347
        self.cont = cont
3✔
348
        self.delta = delta
3✔
349
        self.order = order
3✔
350
        self.ivar = ivar
3✔
351
        self.exposures_diff = exposures_diff
3✔
352
        self.mean_snr = mean_snr
3✔
353
        self.mean_reso = mean_reso
3✔
354
        self.mean_z = mean_z
3✔
355
        self.resolution_matrix = resolution_matrix
3✔
356
        self.mean_resolution_matrix = mean_resolution_matrix
3✔
357
        self.mean_reso_pix = mean_reso_pix
3✔
358

359
        # variables computed in function io.read_deltas
360
        self.z = None
3✔
361
        self.r_comov = None
3✔
362
        self.dist_m = None
3✔
363

364
        # variables computed in function cf.fill_neighs or xcf.fill_neighs
365
        self.neighbours = None
3✔
366

367
        # variables used in function cf.compute_wick_terms and
368
        # main from bin.picca_wick
369
        self.fname = None
3✔
370

371
    @classmethod
3✔
372
    def from_fitsio(cls, hdu, pk1d_type=False, order=None):
3✔
373
        """Initialize instance from a fits file.
374

375
        Args:
376
            hdu: fitsio.hdu.table.TableHDU
377
                A Header Data Unit opened with fitsio
378
            pk1d_type: bool - default: False
379
                Specifies if the fits file is formatted for the 1D Power
380
                Spectrum analysis
381
            order: int - default: None
382
                Order of the polynomial used for the continuum fitting
383
        Returns:
384
            a Delta instance
385
        """
386
        header = hdu.read_header()
3✔
387

388
        # new runs of picca_deltas should have a blinding keyword
389
        if "BLINDING" in header:
3!
UNCOV
390
            blinding = header["BLINDING"]
×
391
        # older runs are not from DESI main survey and should not be blinded
392
        else:
393
            blinding = "none"
3✔
394

395
        if blinding != "none":
3!
396
            delta_name = "DELTA_BLIND"
×
397
        else:
398
            delta_name = "DELTA"
3✔
399

400
        delta = hdu[delta_name][:].astype(float)
3✔
401

402
        if 'LOGLAM' in hdu.get_colnames():
3!
403
            log_lambda = hdu['LOGLAM'][:].astype(float)
3✔
UNCOV
404
        elif 'LAMBDA' in hdu.get_colnames():
×
UNCOV
405
            log_lambda = np.log10(hdu['LAMBDA'][:].astype(float))
×
406
        else:
UNCOV
407
            raise KeyError("Did not find LOGLAM or LAMBDA in delta file")
×
408

409
        if pk1d_type:
3✔
410
            ivar = hdu['IVAR'][:].astype(float)
3✔
411
            try:
3✔
412
                exposures_diff = hdu['DIFF'][:].astype(float)
3✔
413
            except (KeyError, ValueError):
×
UNCOV
414
                userprint('WARNING: no DIFF in hdu while pk1d_type=True, filling with zeros.')
×
UNCOV
415
                exposures_diff = np.zeros(delta.shape)
×
416
            mean_snr = header['MEANSNR']
3✔
417
            mean_reso = header['MEANRESO']
3✔
418
            try:
3✔
419
                mean_reso_pix = header['MEANRESO_PIX']
3✔
420
            except (KeyError, ValueError):
3✔
421
                mean_reso_pix = None
3✔
422

423
            mean_z = header['MEANZ']
3✔
424
            try:
3✔
425
                #transposing here gives back the actual reso matrix which has been stored transposed
426
                resolution_matrix = hdu['RESOMAT'][:].T.astype(float)
3✔
UNCOV
427
                if resolution_matrix is not None:
×
UNCOV
428
                    mean_resolution_matrix = np.mean(resolution_matrix, axis=1)
×
429
                else:
UNCOV
430
                    mean_resolution_matrix = None
×
431
            except (KeyError, ValueError):
3✔
432
                resolution_matrix = None
3✔
433
                mean_resolution_matrix = None
3✔
434
            weights = None
3✔
435
            cont = None
3✔
436
        else:
437
            ivar = None
3✔
438
            exposures_diff = None
3✔
439
            mean_snr = None
3✔
440
            mean_reso = None
3✔
441
            mean_z = None
3✔
442
            resolution_matrix = None
3✔
443
            mean_resolution_matrix = None
3✔
444
            mean_reso_pix = None
3✔
445
            weights = hdu['WEIGHT'][:].astype(float)
3✔
446
            cont = hdu['CONT'][:].astype(float)
3✔
447

448
        if 'THING_ID' in header:
3!
449
            los_id = header['THING_ID']
3✔
450
            plate = header['PLATE']
3✔
451
            mjd = header['MJD']
3✔
452
            fiberid = header['FIBERID']
3✔
UNCOV
453
        elif 'LOS_ID' in header:
×
UNCOV
454
            los_id = header['LOS_ID']
×
UNCOV
455
            plate=los_id
×
UNCOV
456
            mjd=los_id
×
UNCOV
457
            fiberid=los_id
×
458
        else:
459
            raise Exception("Could not find THING_ID or LOS_ID")
×
460

461
        ra = header['RA']
3✔
462
        dec = header['DEC']
3✔
463
        z_qso = header['Z']
3✔
464

465
        return cls(los_id, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
3✔
466
                   weights, cont, delta, order, ivar, exposures_diff, mean_snr,
467
                   mean_reso, mean_z, resolution_matrix,
468
                   mean_resolution_matrix, mean_reso_pix)
469

470
    @classmethod
3✔
471
    def from_ascii(cls, line):
3✔
472
        """Initialize instance from an ascii file.
473

474
        Args:
475
            line: string
476
                A line of the ascii file containing information from a line
477
                of sight
478

479
        Returns:
480
            a Delta instance
481
        """
482

UNCOV
483
        cols = line.split()
×
UNCOV
484
        plate = int(cols[0])
×
UNCOV
485
        mjd = int(cols[1])
×
UNCOV
486
        fiberid = int(cols[2])
×
UNCOV
487
        ra = float(cols[3])
×
UNCOV
488
        dec = float(cols[4])
×
489
        z_qso = float(cols[5])
×
490
        mean_z = float(cols[6])
×
491
        mean_snr = float(cols[7])
×
492
        mean_reso = float(cols[8])
×
493
        delta_log_lambda = float(cols[9])
×
494

495
        num_pixels = int(cols[10])
×
496
        delta = np.array(cols[11:11 + num_pixels]).astype(float)
×
497
        log_lambda = np.array(cols[11 + num_pixels:11 +
×
498
                                   2 * num_pixels]).astype(float)
499
        ivar = np.array(cols[11 + 2 * num_pixels:11 +
×
500
                             3 * num_pixels]).astype(float)
501
        exposures_diff = np.array(cols[11 + 3 * num_pixels:11 +
×
502
                                       4 * num_pixels]).astype(float)
503

UNCOV
504
        thingid = 0
×
505
        order = 0
×
UNCOV
506
        weights = None
×
507
        cont = None
×
508

UNCOV
509
        return cls(thingid, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
×
510
                   weights, cont, delta, order, ivar, exposures_diff, mean_snr,
511
                   mean_reso, mean_z, delta_log_lambda)
512

513
    @classmethod
3✔
514
    def from_image(cls, hdul, pk1d_type=False, z_min_qso=0, z_max_qso=10, order=None):
3✔
515
        """Initialize instance from an ascii file.
516

517
        Args:
518
            hdu: fitsio.hdu.table.TableHDU
519
                A Header Data Unit opened with fitsio
520
            pk1d_type: bool - default: False
521
                Specifies if the fits file is formatted for the 1D Power
522
                Spectrum analysis
523
            z_min_qso: float - default: 0
524
                Specifies the minimum redshift for QSOs
525
            z_max_qso: float - default: 10
526
                Specifies the maximum redshift for QSOs
527
            order: int - default: None
528
                Order of the polynomial used for the continuum fitting
529
        Returns:
530
            a Delta instance
531
        """
UNCOV
532
        if pk1d_type:
×
533
            raise ValueError("ImageHDU format not implemented for Pk1D forests.")
×
534

535
        header = hdul["METADATA"].read_header()
×
536
        N_forests = hdul["METADATA"].get_nrows()
×
UNCOV
537
        Nones = np.full(N_forests, None)
×
538

539
        # new runs of picca_deltas should have a blinding keyword
540
        if "BLINDING" in header:
×
541
            blinding = header["BLINDING"]
×
542
        else:
543
            blinding = "none"
×
544

545
        if blinding != "none":
×
UNCOV
546
            delta_name = "DELTA_BLIND"
×
547
        else:
548
            delta_name = "DELTA"
×
549

UNCOV
550
        delta = hdul[delta_name].read().astype(float)
×
551

UNCOV
552
        if "LOGLAM" in hdul:
×
553
            log_lambda = hdul["LOGLAM"][:].astype(float)
×
554
        elif "LAMBDA" in hdul:
×
UNCOV
555
            log_lambda = np.log10(hdul["LAMBDA"][:].astype(float))
×
556
        else:
UNCOV
557
            raise KeyError("Did not find LOGLAM or LAMBDA in delta file")
×
558

UNCOV
559
        ivar = Nones
×
560
        exposures_diff = Nones
×
561
        mean_snr = Nones
×
562
        mean_reso = Nones
×
563
        mean_z = Nones
×
UNCOV
564
        resolution_matrix = Nones
×
565
        mean_resolution_matrix = Nones
×
UNCOV
566
        mean_reso_pix = Nones
×
567
        weights = hdul["WEIGHT"].read().astype(float)
×
568
        w = weights > 0
×
569
        cont = hdul["CONT"].read().astype(float)
×
570

571
        if "THING_ID" in hdul["METADATA"].get_colnames():
×
572
            los_id = hdul["METADATA"]["THING_ID"][:]
×
573
            plate = hdul["METADATA"]["PLATE"][:]
×
574
            mjd = hdul["METADATA"]["MJD"][:]
×
575
            fiberid=hdul["METADATA"]["FIBERID"][:]
×
576
        elif "LOS_ID" in hdul["METADATA"].get_colnames():
×
577
            los_id = hdul["METADATA"]["LOS_ID"][:]
×
UNCOV
578
            plate=los_id
×
579
            mjd=los_id
×
580
            fiberid=los_id
×
581
        else:
582
            raise Exception("Could not find THING_ID or LOS_ID")
×
583

584
        ra = hdul["METADATA"]["RA"][:]
×
585
        dec = hdul["METADATA"]["DEC"][:]
×
586
        z_qso = hdul["METADATA"]["Z"][:]
×
587
        
588
        deltas = []
×
UNCOV
589
        for (los_id_i, ra_i, dec_i, z_qso_i, plate_i, mjd_i, fiberid_i, log_lambda,
×
590
            weights_i, cont_i, delta_i, ivar_i, exposures_diff_i, mean_snr_i,
591
            mean_reso_i, mean_z_i, resolution_matrix_i,
592
            mean_resolution_matrix_i, mean_reso_pix_i, w_i
593
        ) in zip(los_id, ra, dec, z_qso, plate, mjd, fiberid, repeat(log_lambda),
594
                   weights, cont, delta, ivar, exposures_diff, mean_snr,
595
                   mean_reso, mean_z, resolution_matrix,
596
                   mean_resolution_matrix, mean_reso_pix, w):
597
            if z_qso_i >= z_min_qso and z_qso_i <= z_max_qso:        
×
598
                deltas.append(cls(
×
599
                    los_id_i, ra_i, dec_i, z_qso_i, plate_i, mjd_i, fiberid_i, log_lambda[w_i],
600
                    weights_i[w_i] if weights_i is not None else None, 
601
                    cont_i[w_i], 
602
                    delta_i[w_i],
603
                    order, 
604
                    ivar_i[w_i] if ivar_i is not None else None,
605
                    exposures_diff_i[w_i] if exposures_diff_i is not None else None, 
606
                    mean_snr_i, mean_reso_i, mean_z_i,
607
                    resolution_matrix_i if resolution_matrix_i is not None else None,
608
                    mean_resolution_matrix_i if mean_resolution_matrix_i is not None else None,
609
                    mean_reso_pix_i,
610
                ))
611

UNCOV
612
        return deltas
×
613

614
    def project(self):
3✔
615
        """Project the delta field.
616

617
        The projection gets rid of the distortion caused by the continuum
618
        fitiing. See equations 5 and 6 of du Mas des Bourboux et al. 2020
619
        """
620
        # 2nd term in equation 6
621
        sum_weights = np.sum(self.weights)
3✔
622
        if sum_weights > 0.0:
3!
623
            mean_delta = np.average(self.delta, weights=self.weights)
3✔
624
        else:
625
            # should probably write a warning
UNCOV
626
            return
×
627

628
        # 3rd term in equation 6
629
        res = 0
3✔
630
        if (self.order == 1) and self.delta.shape[0] > 1:
3!
631
            mean_log_lambda = np.average(self.log_lambda, weights=self.weights)
3✔
632
            meanless_log_lambda = self.log_lambda - mean_log_lambda
3✔
633
            mean_delta_log_lambda = (
3✔
634
                np.sum(self.weights * self.delta * meanless_log_lambda) /
635
                np.sum(self.weights * meanless_log_lambda**2))
636
            res = mean_delta_log_lambda * meanless_log_lambda
3✔
UNCOV
637
        elif self.order == 1:
×
UNCOV
638
            res = self.delta
×
639

640
        self.delta -= mean_delta + res
3✔
641

642
    def rebin(self, factor, dwave=0.8):
3✔
643
        """Rebin deltas by an integer factor
644

645
        Args:
646
            factor: int
647
                Factor to rebin deltas (new_bin_size = factor * old_bin_size)
648
            dwave: float
649
                Delta lambda of original deltas
650
        """
UNCOV
651
        wave = 10**np.array(self.log_lambda)
×
652

653
        start = wave.min() - dwave / 2
×
UNCOV
654
        num_bins = np.ceil(((wave[-1] - wave[0]) / dwave + 1) / factor)
×
655

UNCOV
656
        edges = np.arange(num_bins) * dwave * factor + start
×
657

UNCOV
658
        new_indx = np.searchsorted(edges, wave)
×
659

UNCOV
660
        binned_delta = np.bincount(new_indx, weights=self.delta*self.weights,
×
661
                                   minlength=edges.size+1)[1:-1]
UNCOV
662
        binned_weight = np.bincount(new_indx, weights=self.weights, minlength=edges.size+1)[1:-1]
×
663

UNCOV
664
        mask = binned_weight != 0
×
UNCOV
665
        binned_delta[mask] /= binned_weight[mask]
×
666

UNCOV
667
        new_wave = (edges[1:] + edges[:-1]) / 2
×
668

669
        self.log_lambda = np.log10(new_wave[mask])
×
UNCOV
670
        self.delta = binned_delta[mask]
×
671
        self.weights = binned_weight[mask]
×
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc