utils.py 28.8 KB
Newer Older
Moritz Huebner's avatar
Moritz Huebner committed
1
from __future__ import division
MoritzThomasHuebner's avatar
MoritzThomasHuebner committed
2

Moritz Huebner's avatar
Moritz Huebner committed
3
4
5
6
import logging
import os
from math import fmod
import argparse
7
import traceback
MoritzThomasHuebner's avatar
MoritzThomasHuebner committed
8
import inspect
9
import subprocess
Colm Talbot's avatar
Colm Talbot committed
10
import json
11
import multiprocessing
Moritz Huebner's avatar
Moritz Huebner committed
12

MoritzThomasHuebner's avatar
MoritzThomasHuebner committed
13
import numpy as np
Colm Talbot's avatar
Colm Talbot committed
14
from scipy.interpolate import interp2d
15
from scipy.special import logsumexp
Colm Talbot's avatar
Colm Talbot committed
16
import pandas as pd
MoritzThomasHuebner's avatar
MoritzThomasHuebner committed
17

Colm Talbot's avatar
Colm Talbot committed
18
logger = logging.getLogger('bilby')
Gregory Ashton's avatar
Gregory Ashton committed
19

Moritz Huebner's avatar
Moritz Huebner committed
20

21
22
23
24
25
# Constants: values taken from astropy v3.0.4
speed_of_light = 299792458.0  # m/s
parsec = 3.0856775814671916e+16  # m
solar_mass = 1.9884754153381438e+30  # Kg
radius_of_earth = 6378100.0  # m
Moritz Huebner's avatar
Moritz Huebner committed
26

Moritz Huebner's avatar
Moritz Huebner committed
27
28
_TOL = 14

Moritz Huebner's avatar
Moritz Huebner committed
29

30
def infer_parameters_from_function(func):
31
32
33
34
35
36
37
38
39
40
    """ Infers the arguments of a function
        (except the first arg which is assumed to be the dep. variable).

        Throws out *args and **kwargs type arguments

        Can deal with type hinting!

        Returns
        ---------
        list: A list of strings with the parameters
41
    """
42
43
44
45
46
47
48
49
50
51
52
53
54
    return _infer_args_from_function_except_for_first_arg(func=func)


def infer_args_from_method(method):
    """ Infers all arguments of a method except for 'self'

    Throws out *args and **kwargs type arguments.

    Can deal with type hinting!

    Returns
    ---------
    list: A list of strings with the parameters
55
    """
56
57
58
59
60
61
62
63
    return _infer_args_from_function_except_for_first_arg(func=method)


def _infer_args_from_function_except_for_first_arg(func):
    try:
        parameters = inspect.getfullargspec(func).args
    except AttributeError:
        parameters = inspect.getargspec(func).args
64
65
66
67
    parameters.pop(0)
    return parameters


Moritz Huebner's avatar
Moritz Huebner committed
68
def get_sampling_frequency(time_array):
Moritz Huebner's avatar
Moritz Huebner committed
69
70
    """
    Calculate sampling frequency from a time series
71

Moritz Huebner's avatar
Moritz Huebner committed
72
73
74
75
76
    Attributes:
    -------
    time_array: array_like
        Time array to get sampling_frequency from

77
78
    Returns
    -------
Moritz Huebner's avatar
Moritz Huebner committed
79
    Sampling frequency of the time series: float
80
81
82

    Raises
    -------
83
84
    ValueError: If the time series is not evenly sampled.

Moritz Huebner's avatar
Moritz Huebner committed
85
86
    """
    tol = 1e-10
Moritz Huebner's avatar
Moritz Huebner committed
87
    if np.ptp(np.diff(time_array)) > tol:
Moritz Huebner's avatar
Moritz Huebner committed
88
89
        raise ValueError("Your time series was not evenly sampled")
    else:
Moritz Huebner's avatar
Moritz Huebner committed
90
        return np.round(1. / (time_array[1] - time_array[0]), decimals=_TOL)
Moritz Huebner's avatar
Moritz Huebner committed
91
92


93
94
95
96
def get_sampling_frequency_and_duration_from_time_array(time_array):
    """
    Calculate sampling frequency and duration from a time array

Moritz Huebner's avatar
Moritz Huebner committed
97
98
99
100
101
    Attributes:
    -------
    time_array: array_like
        Time array to get sampling_frequency/duration from: array_like

102
103
    Returns
    -------
Moritz Huebner's avatar
Moritz Huebner committed
104
    sampling_frequency, duration: float, float
105
106
107
108
109
110
111
112

    Raises
    -------
    ValueError: If the time_array is not evenly sampled.

    """

    sampling_frequency = get_sampling_frequency(time_array)
Moritz Huebner's avatar
Moritz Huebner committed
113
    duration = len(time_array) / sampling_frequency
114
115
116
    return sampling_frequency, duration


117
118
119
120
def get_sampling_frequency_and_duration_from_frequency_array(frequency_array):
    """
    Calculate sampling frequency and duration from a frequency array

Moritz Huebner's avatar
Moritz Huebner committed
121
122
123
124
125
    Attributes:
    -------
    frequency_array: array_like
        Frequency array to get sampling_frequency/duration from: array_like

126
127
    Returns
    -------
Moritz Huebner's avatar
Moritz Huebner committed
128
    sampling_frequency, duration: float, float
129
130
131
132
133
134
135
136
137
138
139
140
141

    Raises
    -------
    ValueError: If the frequency_array is not evenly sampled.

    """

    tol = 1e-10
    if np.ptp(np.diff(frequency_array)) > tol:
        raise ValueError("Your frequency series was not evenly sampled")

    number_of_frequencies = len(frequency_array)
    delta_freq = frequency_array[1] - frequency_array[0]
Moritz Huebner's avatar
Moritz Huebner committed
142
143
144
    duration = np.round(1 / delta_freq, decimals=_TOL)

    sampling_frequency = np.round(2 * (number_of_frequencies - 1) / duration, decimals=14)
145
146
147
    return sampling_frequency, duration


Moritz Huebner's avatar
Moritz Huebner committed
148
def create_time_series(sampling_frequency, duration, starting_time=0.):
149
150
151
152
153
154
155
156
157
158
159
160
161
    """

    Parameters
    ----------
    sampling_frequency: float
    duration: float
    starting_time: float, optional

    Returns
    -------
    float: An equidistant time series given the parameters

    """
Moritz Huebner's avatar
Moritz Huebner committed
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
    _check_legal_sampling_frequency_and_duration(sampling_frequency, duration)
    number_of_samples = int(duration * sampling_frequency)
    return np.linspace(start=starting_time,
                       stop=duration + starting_time - 1 / sampling_frequency,
                       num=number_of_samples)


def create_frequency_series(sampling_frequency, duration):
    """ Create a frequency series with the correct length and spacing.

    Parameters
    -------
    sampling_frequency: float
    duration: float

    Returns
    -------
    array_like: frequency series

    """
    _check_legal_sampling_frequency_and_duration(sampling_frequency, duration)
    number_of_samples = int(np.round(duration * sampling_frequency))
    number_of_frequencies = int(np.round(number_of_samples / 2) + 1)

    return np.linspace(start=0,
                       stop=sampling_frequency / 2,
                       num=number_of_frequencies)


def _check_legal_sampling_frequency_and_duration(sampling_frequency, duration):
    """ By convention, sampling_frequency and duration have to multiply to an integer

    This will check if the product of both parameters multiplies reasonably close
    to an integer.

    Parameters
    -------
    sampling_frequency: float
    duration: float

    """
    num = sampling_frequency * duration
204
    if np.abs(num - np.round(num)) > 10**(-_TOL):
Moritz Huebner's avatar
Moritz Huebner committed
205
206
207
208
209
210
211
212
        raise IllegalDurationAndSamplingFrequencyException(
            '\nYour sampling frequency and duration must multiply to a number'
            'up to (tol = {}) decimals close to an integer number. '
            '\nBut sampling_frequency={} and  duration={} multiply to {}'.format(
                _TOL, sampling_frequency, duration,
                sampling_frequency * duration
            )
        )
Moritz Huebner's avatar
Moritz Huebner committed
213
214
215


def ra_dec_to_theta_phi(ra, dec, gmst):
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
    """ Convert from RA and DEC to polar coordinates on celestial sphere

    Parameters
    -------
    ra: float
        right ascension in radians
    dec: float
        declination in radians
    gmst: float
        Greenwich mean sidereal time of arrival of the signal in radians

    Returns
    -------
    float: zenith angle in radians
    float: azimuthal angle in radians

Moritz Huebner's avatar
Moritz Huebner committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    """
    phi = ra - gmst
    theta = np.pi / 2 - dec
    return theta, phi


def gps_time_to_gmst(gps_time):
    """
    Convert gps time to Greenwich mean sidereal time in radians

    This method assumes a constant rotation rate of earth since 00:00:00, 1 Jan. 2000
    A correction has been applied to give the exact correct value for 00:00:00, 1 Jan. 2018
    Error accumulates at a rate of ~0.0001 radians/decade.

246
247
248
249
250
251
252
253
254
    Parameters
    -------
    gps_time: float
        gps time

    Returns
    -------
    float: Greenwich mean sidereal time in radians

Moritz Huebner's avatar
Moritz Huebner committed
255
256
257
258
259
260
261
262
263
264
265
    """
    omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
    gps_2000 = 630720013.
    gmst_2000 = (6 + 39. / 60 + 51.251406103947375 / 3600) * np.pi / 12
    correction_2018 = -0.00017782487379358614
    sidereal_time = omega_earth * (gps_time - gps_2000) + gmst_2000 + correction_2018
    gmst = fmod(sidereal_time, 2 * np.pi)
    return gmst


def create_white_noise(sampling_frequency, duration):
266
    """ Create white_noise which is then coloured by a given PSD
Moritz Huebner's avatar
Moritz Huebner committed
267

268
269
270
271
272
    Parameters
    -------
    sampling_frequency: float
    duration: float
        duration of the data
Moritz Huebner's avatar
Moritz Huebner committed
273

274
275
    Returns
    -------
276
277
    array_like: white noise
    array_like: frequency array
Moritz Huebner's avatar
Moritz Huebner committed
278
279
280
281
282
    """

    number_of_samples = duration * sampling_frequency
    number_of_samples = int(np.round(number_of_samples))

Gregory Ashton's avatar
Gregory Ashton committed
283
    delta_freq = 1. / duration
Moritz Huebner's avatar
Moritz Huebner committed
284
285
286

    frequencies = create_frequency_series(sampling_frequency, duration)

Gregory Ashton's avatar
Gregory Ashton committed
287
    norm1 = 0.5 * (1. / delta_freq)**0.5
Moritz Huebner's avatar
Moritz Huebner committed
288
289
    re1 = np.random.normal(0, norm1, len(frequencies))
    im1 = np.random.normal(0, norm1, len(frequencies))
Gregory Ashton's avatar
Gregory Ashton committed
290
    htilde1 = re1 + 1j * im1
Moritz Huebner's avatar
Moritz Huebner committed
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308

    # convolve data with instrument transfer function
    otilde1 = htilde1 * 1.
    # set DC and Nyquist = 0
    otilde1[0] = 0
    # no Nyquist frequency when N=odd
    if np.mod(number_of_samples, 2) == 0:
        otilde1[-1] = 0

    # normalise for positive frequencies and units of strain/rHz
    white_noise = otilde1
    # python: transpose for use with infft
    white_noise = np.transpose(white_noise)
    frequencies = np.transpose(frequencies)

    return white_noise, frequencies


309
310
def nfft(time_domain_strain, sampling_frequency):
    """ Perform an FFT while keeping track of the frequency bins. Assumes input
Gregory Ashton's avatar
Gregory Ashton committed
311
        time series is real (positive frequencies only).
Moritz Huebner's avatar
Moritz Huebner committed
312

313
    Parameters
314
315
316
317
318
    ----------
    time_domain_strain: array_like
        Time series of strain data.
    sampling_frequency: float
        Sampling frequency of the data.
Moritz Huebner's avatar
Moritz Huebner committed
319

320
321
    Returns
    -------
Moritz Huebner's avatar
Moritz Huebner committed
322
    frequency_domain_strain, frequency_array: (array_like, array_like)
323
        Single-sided FFT of time domain strain normalised to units of
Gregory Ashton's avatar
Gregory Ashton committed
324
        strain / Hz, and the associated frequency_array.
325

Moritz Huebner's avatar
Moritz Huebner committed
326
    """
327
    frequency_domain_strain = np.fft.rfft(time_domain_strain)
Moritz Huebner's avatar
Moritz Huebner committed
328
    frequency_domain_strain /= sampling_frequency
Moritz Huebner's avatar
Moritz Huebner committed
329

Moritz Huebner's avatar
Moritz Huebner committed
330
331
    frequency_array = np.linspace(
        0, sampling_frequency / 2, len(frequency_domain_strain))
Moritz Huebner's avatar
Moritz Huebner committed
332

Moritz Huebner's avatar
Moritz Huebner committed
333
    return frequency_domain_strain, frequency_array
Moritz Huebner's avatar
Moritz Huebner committed
334
335


336
def infft(frequency_domain_strain, sampling_frequency):
Gregory Ashton's avatar
Gregory Ashton committed
337
    """ Inverse FFT for use in conjunction with nfft.
338
339

    Parameters
340
341
342
    ----------
    frequency_domain_strain: array_like
        Single-sided, normalised FFT of the time-domain strain data (in units
Gregory Ashton's avatar
Gregory Ashton committed
343
        of strain / Hz).
Moritz Huebner's avatar
Moritz Huebner committed
344
    sampling_frequency: int, float
345
        Sampling frequency of the data.
346
347
348

    Returns
    -------
Moritz Huebner's avatar
Moritz Huebner committed
349
    time_domain_strain: array_like
350
        An array of the time domain strain
Moritz Huebner's avatar
Moritz Huebner committed
351
    """
352
353
354
    time_domain_strain_norm = np.fft.irfft(frequency_domain_strain)
    time_domain_strain = time_domain_strain_norm * sampling_frequency
    return time_domain_strain
Moritz Huebner's avatar
Moritz Huebner committed
355
356


357
def setup_logger(outdir=None, label=None, log_level='INFO', print_version=False):
Moritz Huebner's avatar
Moritz Huebner committed
358
359
360
361
362
363
    """ Setup logging output: call at the start of the script to use

    Parameters
    ----------
    outdir, label: str
        If supplied, write the logging output to outdir/label.log
364
365
366
    log_level: str, optional
        ['debug', 'info', 'warning']
        Either a string from the list above, or an integer as specified
Moritz Huebner's avatar
Moritz Huebner committed
367
        in https://docs.python.org/2/library/logging.html#logging-levels
Gregory Ashton's avatar
Gregory Ashton committed
368
369
    print_version: bool
        If true, print version information
Moritz Huebner's avatar
Moritz Huebner committed
370
371
372
373
    """

    if type(log_level) is str:
        try:
374
            level = getattr(logging, log_level.upper())
Moritz Huebner's avatar
Moritz Huebner committed
375
376
377
        except AttributeError:
            raise ValueError('log_level {} not understood'.format(log_level))
    else:
378
        level = int(log_level)
Moritz Huebner's avatar
Moritz Huebner committed
379

Colm Talbot's avatar
Colm Talbot committed
380
    logger = logging.getLogger('bilby')
Gregory Ashton's avatar
Gregory Ashton committed
381
    logger.propagate = False
382
    logger.setLevel(level)
Moritz Huebner's avatar
Moritz Huebner committed
383

Gregory Ashton's avatar
Gregory Ashton committed
384
385
386
    if any([type(h) == logging.StreamHandler for h in logger.handlers]) is False:
        stream_handler = logging.StreamHandler()
        stream_handler.setFormatter(logging.Formatter(
Gregory Ashton's avatar
Gregory Ashton committed
387
            '%(asctime)s %(name)s %(levelname)-8s: %(message)s', datefmt='%H:%M'))
388
        stream_handler.setLevel(level)
Gregory Ashton's avatar
Gregory Ashton committed
389
390
391
392
393
394
395
396
397
398
399
400
401
        logger.addHandler(stream_handler)

    if any([type(h) == logging.FileHandler for h in logger.handlers]) is False:
        if label:
            if outdir:
                check_directory_exists_and_if_not_mkdir(outdir)
            else:
                outdir = '.'
            log_file = '{}/{}.log'.format(outdir, label)
            file_handler = logging.FileHandler(log_file)
            file_handler.setFormatter(logging.Formatter(
                '%(asctime)s %(levelname)-8s: %(message)s', datefmt='%H:%M'))

402
            file_handler.setLevel(level)
Gregory Ashton's avatar
Gregory Ashton committed
403
            logger.addHandler(file_handler)
Moritz Huebner's avatar
Moritz Huebner committed
404

Gregory Ashton's avatar
Gregory Ashton committed
405
    for handler in logger.handlers:
406
        handler.setLevel(level)
Gregory Ashton's avatar
Gregory Ashton committed
407

Gregory Ashton's avatar
Gregory Ashton committed
408
    if print_version:
Matthew David Pitkin's avatar
Matthew David Pitkin committed
409
        version = get_version_information()
Colm Talbot's avatar
Colm Talbot committed
410
        logger.info('Running bilby version: {}'.format(version))
Moritz Huebner's avatar
Moritz Huebner committed
411
412


Matthew David Pitkin's avatar
Matthew David Pitkin committed
413
414
415
416
417
418
def get_version_information():
    version_file = os.path.join(
        os.path.dirname(os.path.dirname(__file__)), '.version')
    try:
        with open(version_file, 'r') as f:
            return f.readline().rstrip()
419
420
    except EnvironmentError:
        print("No version information file '.version' found")
Matthew David Pitkin's avatar
Matthew David Pitkin committed
421
422


Moritz Huebner's avatar
Moritz Huebner committed
423
def get_progress_bar(module='tqdm'):
424
425
426
    """
    TODO: Write proper docstring
    """
Moritz Huebner's avatar
Moritz Huebner committed
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
    if module in ['tqdm']:
        try:
            from tqdm import tqdm
        except ImportError:
            def tqdm(x, *args, **kwargs):
                return x
        return tqdm
    elif module in ['tqdm_notebook']:
        try:
            from tqdm import tqdm_notebook as tqdm
        except ImportError:
            def tqdm(x, *args, **kwargs):
                return x
        return tqdm


def spherical_to_cartesian(radius, theta, phi):
444
445
446
447
448
449
450
451
452
453
    """ Convert from spherical coordinates to cartesian.

    Parameters
    -------
    radius: float
        radial coordinate
    theta: float
        axial coordinate
    phi: float
        azimuthal coordinate
Moritz Huebner's avatar
Moritz Huebner committed
454

455
456
457
    Returns
    -------
    list: cartesian vector
Moritz Huebner's avatar
Moritz Huebner committed
458
459
460
461
462
463
    """
    cartesian = [radius * np.sin(theta) * np.cos(phi), radius * np.sin(theta) * np.sin(phi), radius * np.cos(theta)]
    return cartesian


def check_directory_exists_and_if_not_mkdir(directory):
464
465
466
467
468
469
470
471
    """ Checks if the given directory exists and creates it if it does not exist

    Parameters
    ----------
    directory: str
        Name of the directory

    """
Moritz Huebner's avatar
Moritz Huebner committed
472
473
    if not os.path.exists(directory):
        os.makedirs(directory)
Gregory Ashton's avatar
Gregory Ashton committed
474
        logger.debug('Making directory {}'.format(directory))
Moritz Huebner's avatar
Moritz Huebner committed
475
    else:
Gregory Ashton's avatar
Gregory Ashton committed
476
        logger.debug('Directory {} exists'.format(directory))
Moritz Huebner's avatar
Moritz Huebner committed
477
478
479


def set_up_command_line_arguments():
480
    """ Sets up command line arguments that can be used to modify how scripts are run.
481
482
483

    Returns
    -------
484
485
486
487
488
489
490
491
492
    command_line_args, command_line_parser: tuple
        The command_line_args is a Namespace of the command line arguments while
        the command_line_parser can be given to a new `argparse.ArgumentParser`
        as a parent object from which to inherit.

    Notes
    -----
        The command line arguments are passed initially at runtime, but this parser
        does not have a `--help` option (i.e., the command line options are
Colm Talbot's avatar
Colm Talbot committed
493
        available for any script which includes `import bilby`, but no help command
494
495
496
497
498
499
        is available. This is done to avoid conflicts with child argparse routines
        (see the example below).

    Example
    -------
    In the following example we demonstrate how to setup a custom command line for a
Colm Talbot's avatar
Colm Talbot committed
500
    project which uses bilby.
501

Colm Talbot's avatar
Colm Talbot committed
502
503
        # Here we import bilby, which initialses and parses the default command-line args
        >>> import bilby
504
        # The command line arguments can then be accessed via
Colm Talbot's avatar
Colm Talbot committed
505
        >>> bilby.core.utils.command_line_args
506
507
508
        Namespace(clean=False, log_level=20, quite=False)
        # Next, we import argparse and define a new argparse object
        >>> import argparse
Colm Talbot's avatar
Colm Talbot committed
509
        >>> parser = argparse.ArgumentParser(parents=[bilby.core.utils.command_line_parser])
510
511
512
513
        >>> parser.add_argument('--argument', type=int, default=1)
        >>> args = parser.parse_args()
        Namespace(clean=False, log_level=20, quite=False, argument=1)

Colm Talbot's avatar
Colm Talbot committed
514
    Placing these lines into a script, you'll be able to pass in the usual bilby default
515
516
    arguments, in addition to `--argument`. To see a list of all options, call the script
    with `--help`.
517
518

    """
519
520
521
522
523
524
525
526
    try:
        parser = argparse.ArgumentParser(
            description="Command line interface for bilby scripts",
            add_help=False, allow_abbrev=False)
    except TypeError:
        parser = argparse.ArgumentParser(
            description="Command line interface for bilby scripts",
            add_help=False)
Moritz Huebner's avatar
Moritz Huebner committed
527
528
529
    parser.add_argument("-v", "--verbose", action="store_true",
                        help=("Increase output verbosity [logging.DEBUG]." +
                              " Overridden by script level settings"))
Colm Talbot's avatar
Colm Talbot committed
530
    parser.add_argument("-q", "--quiet", action="store_true",
Moritz Huebner's avatar
Moritz Huebner committed
531
532
533
534
535
536
                        help=("Decrease output verbosity [logging.WARNING]." +
                              " Overridden by script level settings"))
    parser.add_argument("-c", "--clean", action="store_true",
                        help="Force clean data, never use cached data")
    parser.add_argument("-u", "--use-cached", action="store_true",
                        help="Force cached data and do not check its validity")
Gregory Ashton's avatar
Gregory Ashton committed
537
538
    parser.add_argument("--sampler-help", nargs='?', default=False,
                        const='None', help="Print help for given sampler")
539
    parser.add_argument("--bilby-test-mode", action="store_true",
Moritz Huebner's avatar
Moritz Huebner committed
540
541
                        help=("Used for testing only: don't run full PE, but"
                              " just check nothing breaks"))
Gregory Ashton's avatar
Gregory Ashton committed
542
543
544
    parser.add_argument("--bilby-zero-likelihood-mode", action="store_true",
                        help=("Used for testing only: don't run full PE, but"
                              " just check nothing breaks"))
545
    args, unknown_args = parser.parse_known_args()
Colm Talbot's avatar
Colm Talbot committed
546
    if args.quiet:
Moritz Huebner's avatar
Moritz Huebner committed
547
548
549
550
551
552
        args.log_level = logging.WARNING
    elif args.verbose:
        args.log_level = logging.DEBUG
    else:
        args.log_level = logging.INFO

553
    return args, parser
Moritz Huebner's avatar
Moritz Huebner committed
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
def derivatives(vals, func, releps=1e-3, abseps=None, mineps=1e-9, reltol=1e-3,
                epsscale=0.5, nonfixedidx=None):
    """
    Calculate the partial derivatives of a function at a set of values. The
    derivatives are calculated using the central difference, using an iterative
    method to check that the values converge as step size decreases.

    Parameters
    ----------
    vals: array_like
        A set of values, that are passed to a function, at which to calculate
        the gradient of that function
    func:
        A function that takes in an array of values.
    releps: float, array_like, 1e-3
        The initial relative step size for calculating the derivative.
    abseps: float, array_like, None
        The initial absolute step size for calculating the derivative.
        This overrides `releps` if set.
        `releps` is set then that is used.
    mineps: float, 1e-9
        The minimum relative step size at which to stop iterations if no
        convergence is achieved.
    epsscale: float, 0.5
        The factor by which releps if scaled in each iteration.
    nonfixedidx: array_like, None
        An array of indices in `vals` that are _not_ fixed values and therefore
        can have derivatives taken. If `None` then derivatives of all values
        are calculated.
Gregory Ashton's avatar
Gregory Ashton committed
585

586
587
588
589
590
591
592
593
    Returns
    -------
    grads: array_like
        An array of gradients for each non-fixed value.
    """

    if nonfixedidx is None:
        nonfixedidx = range(len(vals))
Gregory Ashton's avatar
Gregory Ashton committed
594

595
596
    if len(nonfixedidx) > len(vals):
        raise ValueError("To many non-fixed values")
Gregory Ashton's avatar
Gregory Ashton committed
597

598
599
600
601
602
603
604
605
606
607
608
    if max(nonfixedidx) >= len(vals) or min(nonfixedidx) < 0:
        raise ValueError("Non-fixed indexes contain non-existant indices")

    grads = np.zeros(len(nonfixedidx))

    # maximum number of times the gradient can change sign
    flipflopmax = 10.

    # set steps
    if abseps is None:
        if isinstance(releps, float):
Gregory Ashton's avatar
Gregory Ashton committed
609
610
611
            eps = np.abs(vals) * releps
            eps[eps == 0.] = releps  # if any values are zero set eps to releps
            teps = releps * np.ones(len(vals))
612
613
614
615
616
617
618
619
620
621
        elif isinstance(releps, (list, np.ndarray)):
            if len(releps) != len(vals):
                raise ValueError("Problem with input relative step sizes")
            eps = np.multiply(np.abs(vals), releps)
            eps[eps == 0.] = np.array(releps)[eps == 0.]
            teps = releps
        else:
            raise RuntimeError("Relative step sizes are not a recognised type!")
    else:
        if isinstance(abseps, float):
Gregory Ashton's avatar
Gregory Ashton committed
622
            eps = abseps * np.ones(len(vals))
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
        elif isinstance(abseps, (list, np.ndarray)):
            if len(abseps) != len(vals):
                raise ValueError("Problem with input absolute step sizes")
            eps = np.array(abseps)
        else:
            raise RuntimeError("Absolute step sizes are not a recognised type!")
        teps = eps

    # for each value in vals calculate the gradient
    count = 0
    for i in nonfixedidx:
        # initial parameter diffs
        leps = eps[i]
        cureps = teps[i]

        flipflop = 0

        # get central finite difference
        fvals = np.copy(vals)
        bvals = np.copy(vals)

        # central difference
Gregory Ashton's avatar
Gregory Ashton committed
645
646
647
        fvals[i] += 0.5 * leps  # change forwards distance to half eps
        bvals[i] -= 0.5 * leps  # change backwards distance to half eps
        cdiff = (func(fvals) - func(bvals)) / leps
648
649

        while 1:
Gregory Ashton's avatar
Gregory Ashton committed
650
651
            fvals[i] -= 0.5 * leps  # remove old step
            bvals[i] += 0.5 * leps
652
653
654
655
656
657
658
659
660
661
662

            # change the difference by a factor of two
            cureps *= epsscale
            if cureps < mineps or flipflop > flipflopmax:
                # if no convergence set flat derivative (TODO: check if there is a better thing to do instead)
                logger.warning("Derivative calculation did not converge: setting flat derivative.")
                grads[count] = 0.
                break
            leps *= epsscale

            # central difference
Gregory Ashton's avatar
Gregory Ashton committed
663
664
665
            fvals[i] += 0.5 * leps  # change forwards distance to half eps
            bvals[i] -= 0.5 * leps  # change backwards distance to half eps
            cdiffnew = (func(fvals) - func(bvals)) / leps
666
667
668
669
670
671

            if cdiffnew == cdiff:
                grads[count] = cdiff
                break

            # check whether previous diff and current diff are the same within reltol
Gregory Ashton's avatar
Gregory Ashton committed
672
            rat = (cdiff / cdiffnew)
673
674
            if np.isfinite(rat) and rat > 0.:
                # gradient has not changed sign
Gregory Ashton's avatar
Gregory Ashton committed
675
                if np.abs(1. - rat) < reltol:
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
                    grads[count] = cdiffnew
                    break
                else:
                    cdiff = cdiffnew
                    continue
            else:
                cdiff = cdiffnew
                flipflop += 1
                continue

        count += 1

    return grads


691
692
693
694
695
696
697
698
def logtrapzexp(lnf, dx):
    """
    Perform trapezium rule integration for the logarithm of a function on a regular grid.

    Parameters
    ----------
    lnf: array_like
        A :class:`numpy.ndarray` of values that are the natural logarithm of a function
Moritz Huebner's avatar
Moritz Huebner committed
699
    dx: Union[array_like, float]
700
701
702
703
704
705
706
707
708
709
        A :class:`numpy.ndarray` of steps sizes between values in the function, or a
        single step size value.

    Returns
    -------
    The natural logarithm of the area under the function.
    """
    return np.log(dx / 2.) + logsumexp([logsumexp(lnf[:-1]), logsumexp(lnf[1:])])


710
711
class SamplesSummary(object):
    """ Object to store a set of samples and calculate summary statistics
Sylvia Biscoveanu's avatar
Sylvia Biscoveanu committed
712
713
714

    Parameters
    ----------
715
716
717
718
719
720
721
    samples: array_like
        Array of samples
    average: str {'median', 'mean'}
        Use either a median average or mean average when calculating relative
        uncertainties
    level: float
        The default confidence interval level, defaults t0 0.9
Sylvia Biscoveanu's avatar
Sylvia Biscoveanu committed
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
    def __init__(self, samples, average='median', confidence_level=.9):
        self.samples = samples
        self.average = average
        self.confidence_level = confidence_level

    @property
    def samples(self):
        return self._samples

    @samples.setter
    def samples(self, samples):
        self._samples = samples

    @property
    def confidence_level(self):
        return self._confidence_level

    @confidence_level.setter
    def confidence_level(self, confidence_level):
        if 0 < confidence_level and confidence_level < 1:
            self._confidence_level = confidence_level
        else:
            raise ValueError("Confidence level must be between 0 and 1")

    @property
    def average(self):
        if self._average == 'mean':
            return self.mean
        elif self._average == 'median':
            return self.median

    @average.setter
    def average(self, average):
        allowed_averages = ['mean', 'median']
        if average in allowed_averages:
            self._average = average
        else:
            raise ValueError("Average {} not in allowed averages".format(average))

    @property
    def median(self):
        return np.median(self.samples, axis=0)

    @property
    def mean(self):
        return np.mean(self.samples, axis=0)

    @property
    def _lower_level(self):
        """ The credible interval lower quantile value """
        return (1 - self.confidence_level) / 2.

    @property
    def _upper_level(self):
        """ The credible interval upper quantile value """
        return (1 + self.confidence_level) / 2.

    @property
    def lower_absolute_credible_interval(self):
        """ Absolute lower value of the credible interval """
        return np.quantile(self.samples, self._lower_level, axis=0)

    @property
    def upper_absolute_credible_interval(self):
        """ Absolute upper value of the credible interval """
        return np.quantile(self.samples, self._upper_level, axis=0)

    @property
    def lower_relative_credible_interval(self):
        """ Relative (to average) lower value of the credible interval """
        return self.lower_absolute_credible_interval - self.average

    @property
    def upper_relative_credible_interval(self):
        """ Relative (to average) upper value of the credible interval """
        return self.upper_absolute_credible_interval - self.average
Sylvia Biscoveanu's avatar
Sylvia Biscoveanu committed
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
def run_commandline(cl, log_level=20, raise_error=True, return_output=True):
    """Run a string cmd as a subprocess, check for errors and return output.

    Parameters
    ----------
    cl: str
        Command to run
    log_level: int
        See https://docs.python.org/2/library/logging.html#logging-levels,
        default is '20' (INFO)

    """

    logger.log(log_level, 'Now executing: ' + cl)
    if return_output:
        try:
            out = subprocess.check_output(
                cl, stderr=subprocess.STDOUT, shell=True,
                universal_newlines=True)
        except subprocess.CalledProcessError as e:
            logger.log(log_level, 'Execution failed: {}'.format(e.output))
            if raise_error:
                raise
            else:
                out = 0
        os.system('\n')
        return(out)
    else:
        process = subprocess.Popen(cl, shell=True)
        process.communicate()


834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
class Counter(object):
    """
    General class to count number of times a function is Called, returns total
    number of function calls
    Parameters
    ----------
    initalval : int, 0
    number to start counting from
    """
    def __init__(self, initval=0):
        self.val = multiprocessing.RawValue('i', initval)
        self.lock = multiprocessing.Lock()

    def increment(self):
        with self.lock:
            self.val.value += 1

    @property
    def value(self):
        return self.val.value


Colm Talbot's avatar
Colm Talbot committed
856
857
class UnsortedInterp2d(interp2d):

Moritz Huebner's avatar
Moritz Huebner committed
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
    def __call__(self, x, y, dx=0, dy=0, assume_sorted=False):
        """  Wrapper to scipy.interpolate.interp2d which preserves the input ordering.

        See https://stackoverflow.com/questions/44941271/scipy-interp2d-returned-function-sorts-input-argument-automatically-and-undesira
        for the implementation details.


        Parameters
        ----------
        x: See superclass
        y: See superclass
        dx: See superclass
        dy: See superclass
        assume_sorted: bool, optional
            This is just a place holder to prevent a warning.
            Overwriting this will not do anything

        Returns
        ----------
        array_like: See superclass
Colm Talbot's avatar
Colm Talbot committed
878

Moritz Huebner's avatar
Moritz Huebner committed
879
        """
Colm Talbot's avatar
Colm Talbot committed
880
        unsorted_idxs = np.argsort(np.argsort(x))
Moritz Huebner's avatar
Moritz Huebner committed
881
        return interp2d.__call__(self, x, y, dx=dx, dy=dy, assume_sorted=False)[unsorted_idxs]
Colm Talbot's avatar
Colm Talbot committed
882
883


884
885
886
887
#  Instantiate the default argument parser at runtime
command_line_args, command_line_parser = set_up_command_line_arguments()
#  Instantiate the default logging
setup_logger(print_version=True, log_level=command_line_args.log_level)
Moritz Huebner's avatar
Moritz Huebner committed
888
889

if 'DISPLAY' in os.environ:
Gregory Ashton's avatar
Gregory Ashton committed
890
    logger.debug("DISPLAY={} environment found".format(os.environ['DISPLAY']))
Moritz Huebner's avatar
Moritz Huebner committed
891
892
    pass
else:
Gregory Ashton's avatar
Gregory Ashton committed
893
    logger.debug('No $DISPLAY environment variable found, so importing \
894
                   matplotlib.pyplot with non-interactive "Agg" backend.')
Moritz Huebner's avatar
Moritz Huebner committed
895
    import matplotlib
Colm Talbot's avatar
Colm Talbot committed
896
    import matplotlib.pyplot as plt
897

898
899
900
    non_gui_backends = matplotlib.rcsetup.non_interactive_bk
    for backend in non_gui_backends:
        try:
Gregory Ashton's avatar
Gregory Ashton committed
901
            logger.debug("Trying backend {}".format(backend))
902
            matplotlib.use(backend, warn=False)
Colm Talbot's avatar
Colm Talbot committed
903
            plt.switch_backend(backend)
904
            break
905
        except Exception:
906
            print(traceback.format_exc())
Moritz Huebner's avatar
Moritz Huebner committed
907
908


Colm Talbot's avatar
Colm Talbot committed
909
910
911
912
913
914
class BilbyJsonEncoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.ndarray):
            return {'__array__': True, 'content': obj.tolist()}
        if isinstance(obj, complex):
            return {'__complex__': True, 'real': obj.real, 'imag': obj.imag}
Moritz Huebner's avatar
Moritz Huebner committed
915
        if isinstance(obj, pd.DataFrame):
Colm Talbot's avatar
Colm Talbot committed
916
917
918
919
920
921
922
923
924
925
926
927
928
929
            return {'__dataframe__': True, 'content': obj.to_dict(orient='list')}
        return json.JSONEncoder.default(self, obj)


def decode_bilby_json(dct):
    if dct.get("__array__", False):
        return np.asarray(dct["content"])
    if dct.get("__complex__", False):
        return complex(dct["real"], dct["imag"])
    if dct.get("__dataframe__", False):
        return pd.DataFrame(dct['content'])
    return dct


Moritz Huebner's avatar
Moritz Huebner committed
930
931
class IllegalDurationAndSamplingFrequencyException(Exception):
    pass