Skip to content

lyapower.utils_fitter

Logging, error models and array rebinning helpers.

lyapower.utils_fitter

Utilities for the lyapower fitter: logging, error models and rebinning.

Provides a small logging helper (:class:Logger), matplotlib/LaTeX number formatting, the pluggable power-spectrum error estimators used by the fits (uncorrelated / constant / computed / computed_epsilon), and n-dimensional array rebinning helpers (including on-disk rebinning of a Nyx/HDF5 simulation field).

Author: Corentin Ravoux

Logger

Bases: object

Thin wrapper over the standard :mod:logging module.

Source code in lyapower/utils_fitter.py 
 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
class Logger(object):
    """Thin wrapper over the standard :mod:`logging` module."""

    def __init__(self, name="Python_Report", log_level="info"):
        """Store the logger name and level.

        Args:
            name (str, optional): Logger / report-file name.
            log_level (str, optional): ``"info"``, ``"debug"`` or ``"warning"``.
        """
        self.name = name
        self.log_level = log_level

    def setup_logging(self):
        """
        Taken from https://nbodykit.readthedocs.io/
        Turn on logging, with the specified level.
        Parameters
        ----------
        log_level : 'info', 'debug', 'warning'
                the logging level to set; logging below this level is ignored
        """

        # This gives:
        #
        # [ 000000.43 ]   0: 06-28 14:49  measurestats	INFO	 Nproc = [2, 1, 1]
        # [ 000000.43 ]   0: 06-28 14:49  measurestats	INFO	 Rmax = 120

        levels = {
            "info": logging.INFO,
            "debug": logging.DEBUG,
            "warning": logging.WARNING,
        }

        logger = logging.getLogger()
        t0 = time.time()

        class Formatter(logging.Formatter):
            def format(self, record):
                s1 = "[ %09.2f ]: " % (time.time() - t0)
                return s1 + logging.Formatter.format(self, record)

        fmt = Formatter(
            fmt="%(asctime)s %(name)-15s %(levelname)-8s %(message)s",
            datefmt="%m-%d %H:%M ",
        )

        global _logging_handler
        if _logging_handler is None:
            _logging_handler = logging.StreamHandler()
            logger.addHandler(_logging_handler)

        _logging_handler.setFormatter(fmt)
        logger.setLevel(levels[self.log_level])

    def setup_report_logging(self):
        """Configure logging to write to the report file ``self.name``."""
        levels = {
            "info": logging.INFO,
            "debug": logging.DEBUG,
            "warning": logging.WARNING,
        }
        logging.basicConfig(
            filename=self.name,
            filemode="w",
            level=levels[self.log_level],
            format="%(asctime)s :: %(levelname)s :: %(message)s",
        )

    @staticmethod
    def add(line, level="info"):
        """Log a line at the given level.

        Args:
            line (str): Message to log.
            level (str, optional): ``"info"``, ``"warning"`` or ``"debug"``.
        """
        if level == "info":
            logging.info(line)
        if level == "warning":
            logging.warning(line)
        if level == "debug":
            logging.debug(line)

    @staticmethod
    def close():
        """Shut down the logging system."""
        logging.shutdown()

setup_logging 

setup_logging()

Taken from https://nbodykit.readthedocs.io/ Turn on logging, with the specified level. Parameters

log_level : 'info', 'debug', 'warning' the logging level to set; logging below this level is ignored

Source code in lyapower/utils_fitter.py 
 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
def setup_logging(self):
    """
    Taken from https://nbodykit.readthedocs.io/
    Turn on logging, with the specified level.
    Parameters
    ----------
    log_level : 'info', 'debug', 'warning'
            the logging level to set; logging below this level is ignored
    """

    # This gives:
    #
    # [ 000000.43 ]   0: 06-28 14:49  measurestats	INFO	 Nproc = [2, 1, 1]
    # [ 000000.43 ]   0: 06-28 14:49  measurestats	INFO	 Rmax = 120

    levels = {
        "info": logging.INFO,
        "debug": logging.DEBUG,
        "warning": logging.WARNING,
    }

    logger = logging.getLogger()
    t0 = time.time()

    class Formatter(logging.Formatter):
        def format(self, record):
            s1 = "[ %09.2f ]: " % (time.time() - t0)
            return s1 + logging.Formatter.format(self, record)

    fmt = Formatter(
        fmt="%(asctime)s %(name)-15s %(levelname)-8s %(message)s",
        datefmt="%m-%d %H:%M ",
    )

    global _logging_handler
    if _logging_handler is None:
        _logging_handler = logging.StreamHandler()
        logger.addHandler(_logging_handler)

    _logging_handler.setFormatter(fmt)
    logger.setLevel(levels[self.log_level])

setup_report_logging 

setup_report_logging()

Configure logging to write to the report file self.name.

Source code in lyapower/utils_fitter.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
def setup_report_logging(self):
    """Configure logging to write to the report file ``self.name``."""
    levels = {
        "info": logging.INFO,
        "debug": logging.DEBUG,
        "warning": logging.WARNING,
    }
    logging.basicConfig(
        filename=self.name,
        filemode="w",
        level=levels[self.log_level],
        format="%(asctime)s :: %(levelname)s :: %(message)s",
    )

add staticmethod 

add(line, level='info')

Log a line at the given level.

Parameters:

Name Type Description Default
line str

Message to log.

 required
level str

"info", "warning" or "debug".

 'info'
Source code in lyapower/utils_fitter.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
@staticmethod
def add(line, level="info"):
    """Log a line at the given level.

    Args:
        line (str): Message to log.
        level (str, optional): ``"info"``, ``"warning"`` or ``"debug"``.
    """
    if level == "info":
        logging.info(line)
    if level == "warning":
        logging.warning(line)
    if level == "debug":
        logging.debug(line)

close staticmethod 

close()

Shut down the logging system.

Source code in lyapower/utils_fitter.py 
156
157
158
159
@staticmethod
def close():
    """Shut down the logging system."""
    logging.shutdown()

create_log 

create_log(name='Python_Report', log_level='info')

Create and initialise a stream :class:Logger.

Parameters:

Name Type Description Default
name str

Logger name.

 'Python_Report'
log_level str

"info", "debug" or "warning".

 'info'

Returns:

Name Type Description
Logger

The configured logger.
Source code in lyapower/utils_fitter.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
def create_log(name="Python_Report", log_level="info"):
    """Create and initialise a stream :class:`Logger`.

    Args:
        name (str, optional): Logger name.
        log_level (str, optional): ``"info"``, ``"debug"`` or ``"warning"``.

    Returns:
        Logger: The configured logger.
    """
    log = Logger(name=name, log_level=log_level)
    log.setup_logging()
    return log

create_report_log 

create_report_log(name='Python_Report', log_level='info')

Create a :class:Logger that writes to a report file.

Parameters:

Name Type Description Default
name str

Output report file name.

 'Python_Report'
log_level str

"info", "debug" or "warning".

 'info'

Returns:

Name Type Description
Logger

The configured file logger.
Source code in lyapower/utils_fitter.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
def create_report_log(name="Python_Report", log_level="info"):
    """Create a :class:`Logger` that writes to a report file.

    Args:
        name (str, optional): Output report file name.
        log_level (str, optional): ``"info"``, ``"debug"`` or ``"warning"``.

    Returns:
        Logger: The configured file logger.
    """
    log = Logger(name=name, log_level=log_level)
    log.setup_report_logging()
    return log

latex_float 

latex_float(float_input, decimals_input='{0:.2g}')

example use: import matplotlib.pyplot as plt plt.figure(),plt.clf() plt.plot(np.array([1,2.]),'ko-',label="$P_0="+latex_float(7.63e-5)+'$'), plt.legend()

Source code in lyapower/utils_fitter.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def latex_float(float_input, decimals_input="{0:.2g}"):
    """
    example use:
    import matplotlib.pyplot as plt
    plt.figure(),plt.clf()
    plt.plot(np.array([1,2.]),'ko-',label="$P_0="+latex_float(7.63e-5)+'$'),
    plt.legend()
    """
    float_str = decimals_input.format(float_input)
    if "e" in float_str:
        base, exponent = float_str.split("e")
        return r"{0} \times 10^{{{1}}}".format(base, int(exponent))
    else:
        return float_str

return_key 

return_key(dictionary, string, default_value)

Return dictionary[string] if present, else default_value.

Parameters:

Name Type Description Default
dictionary dict

Source mapping.

 required
string

Key to look up.

 required
default_value

Value returned when the key is absent.

 required

Returns:

Type Description

The value at string or default_value.
Source code in lyapower/utils_fitter.py 
178
179
180
181
182
183
184
185
186
187
188
189
def return_key(dictionary, string, default_value):
    """Return ``dictionary[string]`` if present, else ``default_value``.

    Args:
        dictionary (dict): Source mapping.
        string: Key to look up.
        default_value: Value returned when the key is absent.

    Returns:
        The value at ``string`` or ``default_value``.
    """
    return dictionary[string] if string in dictionary.keys() else default_value

error_estimator 

error_estimator(power, model='uncorrelated', **kwargs)

Estimate power-spectrum errors with the chosen model.

Parameters:

Name Type Description Default
power ndarray

Power-spectrum values.

 required
model str

"uncorrelated", "constant", "computed" or "computed_epsilon".

 'uncorrelated'
**kwargs

Model inputs (bin_count, epsilon).

 {}

Returns:

Type Description

numpy.ndarray: The estimated errors.

Raises:

Type Description
KeyError

If the model name is unknown.
Source code in lyapower/utils_fitter.py 
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
def error_estimator(power, model="uncorrelated", **kwargs):
    """Estimate power-spectrum errors with the chosen model.

    Args:
        power (numpy.ndarray): Power-spectrum values.
        model (str, optional): ``"uncorrelated"``, ``"constant"``,
            ``"computed"`` or ``"computed_epsilon"``.
        **kwargs: Model inputs (``bin_count``, ``epsilon``).

    Returns:
        numpy.ndarray: The estimated errors.

    Raises:
        KeyError: If the model name is unknown.
    """
    if model == "uncorrelated":
        return error_estimator_uncorrelated(power, **kwargs)
    elif model == "constant":
        return error_estimator_constant(power, **kwargs)
    elif model == "computed":
        return error_estimator_computed(power, **kwargs)
    elif model == "computed_epsilon":
        return error_estimator_computed_epsilon(power, **kwargs)
    else:
        raise KeyError("model of error estimator not available")

error_estimator_uncorrelated 

error_estimator_uncorrelated(power, **kwargs)

Error assuming uncorrelated bins: power * (1/sqrt(N) + epsilon).

Parameters:

Name Type Description Default
power ndarray

Power-spectrum values.

 required
**kwargs

bin_count (counts per bin) and epsilon (floor).

 {}

Returns:

Type Description

numpy.ndarray: The estimated errors.

Raises:

Type Description
KeyError

If bin_count or epsilon is missing.
Source code in lyapower/utils_fitter.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
def error_estimator_uncorrelated(power, **kwargs):
    """Error assuming uncorrelated bins: ``power * (1/sqrt(N) + epsilon)``.

    Args:
        power (numpy.ndarray): Power-spectrum values.
        **kwargs: ``bin_count`` (counts per bin) and ``epsilon`` (floor).

    Returns:
        numpy.ndarray: The estimated errors.

    Raises:
        KeyError: If ``bin_count`` or ``epsilon`` is missing.
    """
    epsilon = return_key(kwargs, "epsilon", None)
    bin_count = return_key(kwargs, "bin_count", None)
    if (bin_count is None) | (epsilon is None):
        raise KeyError("Need bin_count and epsilon")
    return power * ((1 / np.sqrt(bin_count)) + epsilon)

error_estimator_constant 

error_estimator_constant(power, **kwargs)

Error as a constant fraction of the power: power * epsilon.

Parameters:

Name Type Description Default
power ndarray

Power-spectrum values.

 required
**kwargs

epsilon (fractional error).

 {}

Returns:

Type Description

numpy.ndarray: The estimated errors.

Raises:

Type Description
KeyError

If epsilon is missing.
Source code in lyapower/utils_fitter.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def error_estimator_constant(power, **kwargs):
    """Error as a constant fraction of the power: ``power * epsilon``.

    Args:
        power (numpy.ndarray): Power-spectrum values.
        **kwargs: ``epsilon`` (fractional error).

    Returns:
        numpy.ndarray: The estimated errors.

    Raises:
        KeyError: If ``epsilon`` is missing.
    """
    epsilon = return_key(kwargs, "epsilon", None)
    if epsilon is None:
        raise KeyError("Need bin_count and epsilon")
    return power * epsilon

error_estimator_computed 

error_estimator_computed(power, **kwargs)

Error taken directly from a precomputed per-bin error array.

Parameters:

Name Type Description Default
power ndarray

Power-spectrum values (unused).

 required
**kwargs

bin_count (here the precomputed errors).

 {}

Returns:

Type Description

numpy.ndarray: The provided bin_count errors.

Raises:

Type Description
KeyError

If bin_count is missing.
Source code in lyapower/utils_fitter.py 
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
def error_estimator_computed(power, **kwargs):
    """Error taken directly from a precomputed per-bin error array.

    Args:
        power (numpy.ndarray): Power-spectrum values (unused).
        **kwargs: ``bin_count`` (here the precomputed errors).

    Returns:
        numpy.ndarray: The provided ``bin_count`` errors.

    Raises:
        KeyError: If ``bin_count`` is missing.
    """
    bin_count = return_key(kwargs, "bin_count", None)
    if bin_count is None:
        raise KeyError("Need bin_count")
    return bin_count

error_estimator_computed_epsilon 

error_estimator_computed_epsilon(power, **kwargs)

Precomputed error plus a fractional floor: bin_count + epsilon*power.

Parameters:

Name Type Description Default
power ndarray

Power-spectrum values.

 required
**kwargs

bin_count (precomputed errors) and epsilon (floor).

 {}

Returns:

Type Description

numpy.ndarray: The estimated errors.

Raises:

Type Description
KeyError

If bin_count is missing.
Source code in lyapower/utils_fitter.py 
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
def error_estimator_computed_epsilon(power, **kwargs):
    """Precomputed error plus a fractional floor: ``bin_count + epsilon*power``.

    Args:
        power (numpy.ndarray): Power-spectrum values.
        **kwargs: ``bin_count`` (precomputed errors) and ``epsilon`` (floor).

    Returns:
        numpy.ndarray: The estimated errors.

    Raises:
        KeyError: If ``bin_count`` is missing.
    """
    bin_count = return_key(kwargs, "bin_count", None)
    epsilon = return_key(kwargs, "epsilon", None)
    if bin_count is None:
        raise KeyError("Need bin_count")
    return bin_count + epsilon * power

bin_ndarray 

bin_ndarray(ndarray, new_shape, operation='mean')

From : https://stackoverflow.com/questions/8090229/resize-with-averaging-or-rebin-a-numpy-2d-array/29042041 Bins an ndarray in all axes based on the target shape, by summing or averaging. Number of output dimensions must match number of input dimensions. Example

m = np.arange(0,100,1).reshape((10,10)) n = bin_ndarray(m, new_shape=(5,5), operation='sum') print(n) [[ 22 30 38 46 54][102 110 118 126 134] [182 190 198 206 214][262 270 278 286 294] [342 350 358 366 374]]

Source code in lyapower/utils_fitter.py 
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
def bin_ndarray(ndarray, new_shape, operation="mean"):
    """
    From : https://stackoverflow.com/questions/8090229/resize-with-averaging-or-rebin-a-numpy-2d-array/29042041
    Bins an ndarray in all axes based on the target shape, by summing or
    averaging.
    Number of output dimensions must match number of input dimensions.
    Example
    -------
    >>> m = np.arange(0,100,1).reshape((10,10))
    >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum')
    >>> print(n)
    [[ 22  30  38  46  54]
    [102 110 118 126 134]
    [182 190 198 206 214]
    [262 270 278 286 294]
    [342 350 358 366 374]]
    """
    if not operation.lower() in ["sum", "mean", "average", "avg", "gauss"]:
        raise ValueError("Operation {} not supported.".format(operation))
    if ndarray.ndim != len(new_shape):
        raise ValueError("Shape mismatch: {} -> {}".format(ndarray.shape, new_shape))
    compression_pairs = [(d, c // d) for d, c in zip(new_shape, ndarray.shape)]
    flattened = [l for p in compression_pairs for l in p]
    ndarray = ndarray.reshape(flattened)
    for i in range(len(new_shape)):
        if operation.lower() == "sum":
            ndarray = ndarray.sum(-1 * (i + 1))
        elif operation.lower() in ["mean", "average", "avg"]:
            ndarray = ndarray.mean(-1 * (i + 1))
        elif operation.lower() in ["gauss"]:
            if i != 0:
                raise KeyError("gaussian mean is not available for dim higher than 1")
            from scipy import signal

            newndarray = np.zeros(new_shape)
            gaussian_weights = signal.gaussian(
                int(ndarray.shape[1]), int(ndarray.shape[1]) / 4
            )
            for j in range(len(ndarray)):
                newndarray[j] = np.average(ndarray[j], axis=0, weights=gaussian_weights)
            ndarray = newndarray
    return ndarray

rebin_slice 

rebin_slice(sim_name, new_shape_slice, index_rescaling, index, operation='mean', first_field='derived_fields', second_field='tau_red', transform=None)

Read and rebin one x-slice of a Nyx/HDF5 simulation field.

Parameters:

Name Type Description Default
sim_name str

HDF5 simulation file path.

 required
new_shape_slice tuple[int]

Target shape of the rebinned slice.

 required
index_rescaling int

Number of input x-planes per output plane.

 required
index int

Output slice index.

 required
operation str

Reduction ("mean"/"sum" ...).

 'mean'
first_field str

HDF5 group holding the field.

 'derived_fields'
second_field str

HDF5 dataset name.

 'tau_red'
transform callable

Applied to the slice before rebinning.

 None

Returns:

Type Description

numpy.ndarray: The rebinned slice.
Source code in lyapower/utils_fitter.py 
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
def rebin_slice(
    sim_name,
    new_shape_slice,
    index_rescaling,
    index,
    operation="mean",
    first_field="derived_fields",
    second_field="tau_red",
    transform=None,
):
    """Read and rebin one x-slice of a Nyx/HDF5 simulation field.

    Args:
        sim_name (str): HDF5 simulation file path.
        new_shape_slice (tuple[int]): Target shape of the rebinned slice.
        index_rescaling (int): Number of input x-planes per output plane.
        index (int): Output slice index.
        operation (str, optional): Reduction (``"mean"``/``"sum"`` ...).
        first_field (str, optional): HDF5 group holding the field.
        second_field (str, optional): HDF5 dataset name.
        transform (callable, optional): Applied to the slice before rebinning.

    Returns:
        numpy.ndarray: The rebinned slice.
    """
    print("Treating ", index)
    sim = h5py.File(sim_name)
    full_slice = sim[first_field][second_field][
        index * index_rescaling : (index + 1) * index_rescaling, :, :
    ]
    if transform is not None:
        full_slice = transform(full_slice)
    return bin_ndarray(full_slice, new_shape_slice, operation=operation)

rebin_simulation 

rebin_simulation(sim_name, index_rescaling, operation='mean', number_worker=1, first_field='derived_fields', second_field='tau_red', transform_name=None)

Rebin a whole Nyx/HDF5 simulation field slice by slice.

Processes the field one x-slice at a time (optionally in a multiprocessing pool) to keep memory bounded, applying an optional transform first.

Parameters:

Name Type Description Default
sim_name str

HDF5 simulation file path.

 required
index_rescaling int

Downsampling factor along each axis.

 required
operation str

Reduction ("mean"/"sum" ...).

 'mean'
number_worker int

Number of processes (1 = serial).

 1
first_field str

HDF5 group holding the field.

 'derived_fields'
second_field str

HDF5 dataset name.

 'tau_red'
transform_name str

"exp" applies exp(-x) (tau -> transmitted flux) before rebinning.

 None

Returns:

Type Description

numpy.ndarray: The rebinned 3D field.
Source code in lyapower/utils_fitter.py 
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
def rebin_simulation(
    sim_name,
    index_rescaling,
    operation="mean",
    number_worker=1,
    first_field="derived_fields",
    second_field="tau_red",
    transform_name=None,
):
    """Rebin a whole Nyx/HDF5 simulation field slice by slice.

    Processes the field one x-slice at a time (optionally in a multiprocessing
    pool) to keep memory bounded, applying an optional transform first.

    Args:
        sim_name (str): HDF5 simulation file path.
        index_rescaling (int): Downsampling factor along each axis.
        operation (str, optional): Reduction (``"mean"``/``"sum"`` ...).
        number_worker (int, optional): Number of processes (1 = serial).
        first_field (str, optional): HDF5 group holding the field.
        second_field (str, optional): HDF5 dataset name.
        transform_name (str, optional): ``"exp"`` applies ``exp(-x)`` (tau ->
            transmitted flux) before rebinning.

    Returns:
        numpy.ndarray: The rebinned 3D field.
    """
    sim = h5py.File(sim_name)
    shape_sim = sim["domain"].attrs["shape"]
    shape_rebinned_field = (shape_sim / index_rescaling).astype(int)
    rebinned_field = np.zeros(shape_rebinned_field)
    new_shape_slice = (1, shape_rebinned_field[1], shape_rebinned_field[2])
    if transform_name == "exp":
        transform = lambda x: np.exp(-x)
    else:
        transform = None
    func = partial(
        rebin_slice,
        sim_name,
        new_shape_slice,
        index_rescaling,
        operation=operation,
        first_field=first_field,
        second_field=second_field,
        transform=transform,
    )
    if number_worker == 1:
        for i in range(shape_rebinned_field[0]):
            print("Treating ", i)
            rebinned_field[i, :, :] = func(i)
    else:

        with mp.Pool(number_worker) as p:
            results = p.map(func, np.arange(shape_rebinned_field[0]))
        rebinned_field = np.array(results)
    return rebinned_field