Merge remote-tracking branch 'origin/master' into add_bbh_wfg

.gitlab-ci.yml

@@ -38,8 +38,13 @@ python-3:
     - pip install -r optional_requirements.txt
     - pip install 'coverage>=4.5'
     - pip install coverage-badge
+    - pip install flake8
   script:
     - python setup.py install
+
+    # Run pyflakes
+    - flake8 .
+
     # Run tests and collect coverage data
     - coverage --version
     - coverage erase

cli_tupak/plot_multiple_posteriors.py

@@ -4,7 +4,7 @@ import argparse
 def setup_command_line_args():
     parser = argparse.ArgumentParser(
         description="Plot corner plots from results files")
-    parser.add_argument("-r", "--results",  nargs='+',
+    parser.add_argument("-r", "--results", nargs='+',
                         help="List of results files to use.")
     parser.add_argument("-f", "--filename", default=None,
                         help="Output file name.")

examples/injection_examples/binary_neutron_star_example.py

+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation on a binary neutron star
+system taking into account tidal deformabilities.
+
+This example estimates the masses using a uniform prior in both component masses
+and also estimates the tidal deformabilities using a uniform prior in both
+tidal deformabilities
+"""
+
+from __future__ import division, print_function
+
+import numpy as np
+
+import tupak
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'bns_example'
+tupak.core.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility.  This is optional!
+np.random.seed(88170235)
+
+# We are going to inject a binary neutron star waveform.  We first establish a
+# dictionary of parameters that includes all of the different waveform
+# parameters, including masses of the two black holes (mass_1, mass_2),
+# spins of both black holes (a_1,a_2) , etc. 
+injection_parameters = dict(
+    mass_1=1.5, mass_2=1.3, a_1=0.0, a_2=0.0, luminosity_distance=50.,
+    iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
+    ra=1.375, dec=-1.2108, lambda_1=400, lambda_2=450)
+
+# Set the duration and sampling frequency of the data segment that we're going
+# to inject the signal into. For the
+# TaylorF2 waveform, we cut the signal close to the isco frequency
+duration = 8
+sampling_frequency = 2*1570.
+start_time = injection_parameters['geocent_time'] + 2 - duration
+
+# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
+waveform_arguments = dict(waveform_approximant='TaylorF2',
+                          reference_frequency=50., minimum_frequency=40.0)
+
+# Create the waveform_generator using a LAL Binary Neutron Star source function
+waveform_generator = tupak.gw.WaveformGenerator(
+    duration=duration, sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=tupak.gw.source.lal_binary_neutron_star,
+    waveform_arguments=waveform_arguments)
+
+# Set up interferometers.  In this case we'll use three interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
+# These default to their design sensitivity and start at 40 Hz.
+interferometers = tupak.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+for interferometer in interferometers:
+    interferometer.minimum_frequency = 40
+interferometers.set_strain_data_from_power_spectral_densities(
+    sampling_frequency=sampling_frequency, duration=duration,
+    start_time=start_time)
+interferometers.inject_signal(parameters=injection_parameters,
+                              waveform_generator=waveform_generator)
+
+priors = tupak.gw.prior.BNSPriorSet()
+for key in ['a_1', 'a_2', 'psi', 'geocent_time', 'ra', 'dec',
+            'iota', 'luminosity_distance', 'phase']:
+    priors[key] = injection_parameters[key]
+    
+# Initialise the likelihood by passing in the interferometer data (IFOs)
+# and the waveoform generator
+likelihood = tupak.gw.GravitationalWaveTransient(
+    interferometers=interferometers, waveform_generator=waveform_generator,
+    time_marginalization=False, phase_marginalization=False,
+    distance_marginalization=False, prior=priors)
+
+# Run sampler.  In this case we're going to use the `nestle` sampler
+result = tupak.run_sampler(
+    likelihood=likelihood, priors=priors, sampler='nestle', npoints=1000,
+    injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+result.plot_corner()
+

examples/injection_examples/using_gwin.py

+#!/bin/python
+"""
+An example of how to use tupak with `gwin` (https://github.com/gwastro/gwin) to
+perform CBC parameter estimation.
+
+To run this example, it is sufficient to use the pip-installable pycbc package,
+but the source installation of gwin. You can install this by cloning the
+repository (https://github.com/gwastro/gwin) and running
+
+$ python setup.py install
+
+A practical difference between gwin and tupak is that while fixed parameters
+are specified via the prior in tupak, in gwin, these are fixed at instantiation
+of the model. So, in the following, we only create priors for the parameters
+to be searched over.
+
+"""
+from __future__ import division, print_function
+import numpy as np
+import tupak
+
+import gwin
+from pycbc import psd as pypsd
+from pycbc.waveform.generator import (FDomainDetFrameGenerator,
+                                      FDomainCBCGenerator)
+
+label = 'using_gwin'
+
+# Search priors
+priors = dict()
+priors['distance'] = tupak.core.prior.Uniform(500, 2000, 'distance')
+priors['polarization'] = tupak.core.prior.Uniform(0, np.pi, 'iota')
+
+# Data variables
+seglen = 4
+sample_rate = 2048
+N = seglen * sample_rate / 2 + 1
+fmin = 30.
+
+# Injected signal variables
+injection_parameters = dict(mass1=38.6, mass2=29.3, spin1z=0, spin2z=0,
+                            tc=0, ra=3.1, dec=1.37, polarization=2.76,
+                            distance=1500)
+
+# These lines figure out which parameters are to be searched over
+variable_parameters = priors.keys()
+fixed_parameters = injection_parameters.copy()
+for key in priors:
+    fixed_parameters.pop(key)
+
+# These lines generate the `model` object - see https://gwin.readthedocs.io/en/latest/api/gwin.models.gaussian_noise.html
+generator = FDomainDetFrameGenerator(
+    FDomainCBCGenerator, 0.,
+    variable_args=variable_parameters, detectors=['H1', 'L1'],
+    delta_f=1. / seglen, f_lower=fmin,
+    approximant='IMRPhenomPv2', **fixed_parameters)
+signal = generator.generate(**injection_parameters)
+psd = pypsd.aLIGOZeroDetHighPower(int(N), 1. / seglen, 20.)
+psds = {'H1': psd, 'L1': psd}
+model = gwin.models.GaussianNoise(
+    variable_parameters, signal, generator, fmin, psds=psds)
+model.update(**injection_parameters)
+
+
+# This create a dummy class to convert the model into a tupak.likelihood object
+class GWINLikelihood(tupak.core.likelihood.Likelihood):
+    def __init__(self, model):
+        """ A likelihood to wrap around GWIN model objects
+
+        Parameters
+        ----------
+        model: gwin.model.GaussianNoise
+            A gwin model
+
+        """
+        self.model = model
+        self.parameters = {x: None for x in self.model.variable_params}
+
+    def log_likelihood(self):
+        self.model.update(**self.parameters)
+        return self.model.loglikelihood
+
+
+likelihood = GWINLikelihood(model)
+
+
+# Now run the inference
+result = tupak.run_sampler(
+    likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
+    label=label)
+result.plot_corner()
+

examples/other_examples/linear_regression_pymc3_custom_likelihood.py

@@ -99,7 +99,7 @@ class GaussianLikelihoodPyMC3(tupak.Likelihood):
         with sampler.pymc3_model:
             mdist = sampler.pymc3_priors['m']
             cdist = sampler.pymc3_priors['c']
-            
+
             mu = model(time, mdist, cdist)
 
             # set the likelihood distribution

setup.cfg

+[flake8]
+exclude = .git,docs,build,dist,test,examples,*__init__.py
+max-line-length = 160
+ignore = E129

setup.py

@@ -22,8 +22,8 @@ def write_version_file(version):
     try:
         git_log = subprocess.check_output(
             ['git', 'log', '-1', '--pretty=%h %ai']).decode('utf-8')
-        git_diff = (subprocess.check_output(['git', 'diff', '.'])
-                    + subprocess.check_output(
+        git_diff = (subprocess.check_output(['git', 'diff', '.']) +
+                    subprocess.check_output(
                         ['git', 'diff', '--cached', '.'])).decode('utf-8')
         if git_diff == '':
             git_status = '(CLEAN) ' + git_log
@@ -50,7 +50,7 @@ def get_long_description():
     return long_description
 
 
-version = '0.2.1'
+version = '0.2.2'
 version_file = write_version_file(version)
 long_description = get_long_description()
 

test/calibration_tests.py

@@ -11,6 +11,11 @@ class TestBaseClass(unittest.TestCase):
     def tearDown(self):
         del self.model
 
+    def test_repr(self):
+        expected = 'Recalibrate(prefix={})'.format('\'recalib_\'')
+        actual = repr(self.model)
+        self.assertEqual(expected, actual)
+
     def test_calibration_factor(self):
         frequency_array = np.linspace(20, 1024, 1000)
         cal_factor = self.model.get_calibration_factor(frequency_array)
@@ -20,14 +25,22 @@ class TestBaseClass(unittest.TestCase):
 class TestCubicSpline(unittest.TestCase):
 
     def setUp(self):
+        self.prefix = 'recalib_'
+        self.minimum_frequency = 20
+        self.maximum_frequency = 1024
+        self.n_points = 5
         self.model = calibration.CubicSpline(
-            prefix='recalib_', minimum_frequency=20, maximum_frequency=1024,
-            n_points=5)
+            prefix=self.prefix, minimum_frequency=self.minimum_frequency,
+            maximum_frequency=self.maximum_frequency, n_points=self.n_points)
         self.parameters = {'recalib_{}_{}'.format(param, ii): 0.0
                            for ii in range(5)
                            for param in ['amplitude', 'phase']}
 
     def tearDown(self):
+        del self.prefix
+        del self.minimum_frequency
+        del self.maximum_frequency
+        del self.n_points
         del self.model
         del self.parameters
 
@@ -37,6 +50,12 @@ class TestCubicSpline(unittest.TestCase):
                                                        **self.parameters)
         assert np.alltrue(cal_factor.real == np.ones_like(frequency_array))
 
+    def test_repr(self):
+        expected = 'CubicSpline(prefix=\'{}\', minimum_frequency={}, maximum_frequency={}, n_points={})'\
+            .format(self.prefix, self.minimum_frequency, self.maximum_frequency, self.n_points)
+        actual = repr(self.model)
+        self.assertEqual(expected, actual)
+
 
 class TestCubicSplineRequiresFourNodes(unittest.TestCase):
 

test/conversion_tests.py

@@ -8,14 +8,35 @@ import numpy as np
 class TestBasicConversions(unittest.TestCase):
 
     def setUp(self):
-        self.mass_1 = 20
-        self.mass_2 = 10
-        self.mass_ratio = 0.5
-        self.total_mass = 30
-        self.chirp_mass = 200**0.6 / 30**0.2
-        self.symmetric_mass_ratio = 2/9
+        self.mass_1 = 1.4
+        self.mass_2 = 1.3
+        self.mass_ratio = 13/14
+        self.total_mass = 2.7
+        self.chirp_mass = (1.4 * 1.3)**0.6 / 2.7**0.2
+        self.symmetric_mass_ratio = (1.4 * 1.3) / 2.7**2
         self.cos_angle = -1
         self.angle = np.pi
+        self.lambda_1 = 300
+        self.lambda_2 = 300 * (14 / 13)**5
+        self.lambda_tilde = 8 / 13 * (
+            (1 + 7 * self.symmetric_mass_ratio
+             - 31 * self.symmetric_mass_ratio**2)
+            * (self.lambda_1 + self.lambda_2)
+            + (1 - 4 * self.symmetric_mass_ratio)**0.5
+            * (1 + 9 * self.symmetric_mass_ratio
+               - 11 * self.symmetric_mass_ratio**2)
+            * (self.lambda_1 - self.lambda_2)
+        )
+        self.delta_lambda = 1 / 2 * (
+                (1 - 4 * self.symmetric_mass_ratio)**0.5
+                * (1 - 13272 / 1319 * self.symmetric_mass_ratio
+                   + 8944 / 1319 * self.symmetric_mass_ratio**2)
+                * (self.lambda_1 + self.lambda_2)
+                + (1 - 15910 / 1319 * self.symmetric_mass_ratio
+                   + 32850 / 1319 * self.symmetric_mass_ratio**2
+                   + 3380 / 1319 * self.symmetric_mass_ratio**3)
+                * (self.lambda_1 - self.lambda_2)
+        )
 
     def tearDown(self):
         del self.mass_1
@@ -27,7 +48,8 @@ class TestBasicConversions(unittest.TestCase):
 
     def test_total_mass_and_mass_ratio_to_component_masses(self):
         mass_1, mass_2 = tupak.gw.conversion.total_mass_and_mass_ratio_to_component_masses(self.mass_ratio, self.total_mass)
-        self.assertTupleEqual((mass_1, mass_2), (self.mass_1, self.mass_2))
+        self.assertTrue(all([abs(mass_1 - self.mass_1) < 1e-5,
+                             abs(mass_2 - self.mass_2) < 1e-5]))
 
     def test_symmetric_mass_ratio_to_mass_ratio(self):
         mass_ratio = tupak.gw.conversion.symmetric_mass_ratio_to_mass_ratio(self.symmetric_mass_ratio)
@@ -61,6 +83,20 @@ class TestBasicConversions(unittest.TestCase):
         mass_ratio = tupak.gw.conversion.mass_1_and_chirp_mass_to_mass_ratio(self.mass_1, self.chirp_mass)
         self.assertAlmostEqual(self.mass_ratio, mass_ratio)
 
+    def test_lambda_tilde_to_lambda_1_lambda_2(self):
+        lambda_1, lambda_2 =\
+            tupak.gw.conversion.lambda_tilde_to_lambda_1_lambda_2(
+                self.lambda_tilde, self.mass_1, self.mass_2)
+        self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
+                             abs(self.lambda_2 - lambda_2) < 1e-5]))
+
+    def test_lambda_tilde_delta_lambda_to_lambda_1_lambda_2(self):
+        lambda_1, lambda_2 =\
+            tupak.gw.conversion.lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
+                self.lambda_tilde, self.delta_lambda, self.mass_1, self.mass_2)
+        self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
+                             abs(self.lambda_2 - lambda_2) < 1e-5]))
+
 
 class TestConvertToLALBBHParams(unittest.TestCase):
 

test/detector_tests.py

@@ -8,6 +8,8 @@ from mock import patch
 import numpy as np
 import scipy.signal.windows
 import gwpy
+import os
+import logging
 
 
 class TestDetector(unittest.TestCase):
@@ -191,14 +193,6 @@ class TestDetector(unittest.TestCase):
             self.ifo.yarm_azimuth = 12
             self.assertEqual(self.ifo.detector_tensor, 0)
 
-    def test_amplitude_spectral_density_array(self):
-        self.ifo.power_spectral_density.power_spectral_density_interpolated = MagicMock(return_value=np.array([1, 4]))
-        self.assertTrue(np.array_equal(self.ifo.amplitude_spectral_density_array, np.array([1, 2])))
-
-    def test_power_spectral_density_array(self):
-        self.ifo.power_spectral_density.power_spectral_density_interpolated = MagicMock(return_value=np.array([1, 4]))
-        self.assertTrue(np.array_equal(self.ifo.power_spectral_density_array, np.array([1, 4])))
-
     def test_antenna_response_default(self):
         with mock.patch('tupak.gw.utils.get_polarization_tensor') as m:
             with mock.patch('numpy.einsum') as n:
@@ -310,6 +304,17 @@ class TestDetector(unittest.TestCase):
             self.assertTrue(np.array_equal(expected[0], actual[0]))  # array-like element has to be evaluated separately
             self.assertListEqual(expected[1], actual[1])
 
+    def test_repr(self):
+        expected = 'Interferometer(name=\'{}\', power_spectral_density={}, minimum_frequency={}, ' \
+                   'maximum_frequency={}, length={}, latitude={}, longitude={}, elevation={}, xarm_azimuth={}, ' \
+                   'yarm_azimuth={}, xarm_tilt={}, yarm_tilt={})' \
+            .format(self.name, self.power_spectral_density, float(self.minimum_frequency),
+                    float(self.maximum_frequency), float(self.length), float(self.latitude), float(self.longitude),
+                    float(self.elevation), float(self.xarm_azimuth), float(self.yarm_azimuth), float(self.xarm_tilt),
+                    float(self.yarm_tilt))
+        print(repr(self.ifo))
+        self.assertEqual(expected, repr(self.ifo))
+
 
 class TestInterferometerStrainData(unittest.TestCase):
 
@@ -542,10 +547,10 @@ class TestInterferometerStrainData(unittest.TestCase):
     def test_frequency_domain_strain_when_set(self):
         self.ifosd.sampling_frequency = 200
         self.ifosd.duration = 4
-        expected_strain = self.ifosd.frequency_array*self.ifosd.frequency_mask
+        expected_strain = self.ifosd.frequency_array * self.ifosd.frequency_mask
         self.ifosd._frequency_domain_strain = expected_strain
         self.assertTrue(np.array_equal(expected_strain,
-                                        self.ifosd.frequency_domain_strain))
+                                       self.ifosd.frequency_domain_strain))
 
     @patch('tupak.core.utils.nfft')
     def test_frequency_domain_strain_from_frequency_domain_strain(self, m):
@@ -790,5 +795,201 @@ class TestInterferometerList(unittest.TestCase):
         self.assertListEqual([self.ifo1.name, new_ifo.name, self.ifo2.name], names)
 
 
+class TestPowerSpectralDensityWithoutFiles(unittest.TestCase):
+
+    def setUp(self):
+        self.frequency_array = np.array([1., 2., 3.])
+        self.psd_array = np.array([16., 25., 36.])
+        self.asd_array = np.array([4., 5., 6.])
+
+    def tearDown(self):
+        del self.frequency_array
+        del self.psd_array
+        del self.asd_array
+
+    def test_init_with_asd_array(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, asd_array=self.asd_array)
+        self.assertTrue(np.array_equal(self.frequency_array, psd.frequency_array))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+
+    def test_init_with_psd_array(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, psd_array=self.psd_array)
+        self.assertTrue(np.array_equal(self.frequency_array, psd.frequency_array))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+
+    def test_setting_asd_array_after_init(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array)
+        psd.asd_array = self.asd_array
+        self.assertTrue(np.array_equal(self.frequency_array, psd.frequency_array))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+
+    def test_setting_psd_array_after_init(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array)
+        psd.psd_array = self.psd_array
+        self.assertTrue(np.array_equal(self.frequency_array, psd.frequency_array))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+
+    def test_power_spectral_density_interpolated_from_asd_array(self):
+        expected = np.array([25.])
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, asd_array = self.asd_array)
+        self.assertEqual(expected, psd.power_spectral_density_interpolated(2))
+
+    def test_power_spectral_density_interpolated_from_psd_array(self):
+        expected = np.array([25.])
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, psd_array = self.psd_array)
+        self.assertEqual(expected, psd.power_spectral_density_interpolated(2))
+
+    def test_from_amplitude_spectral_density_array(self):
+        actual = tupak.gw.detector.PowerSpectralDensity.from_amplitude_spectral_density_array(
+            frequency_array=self.frequency_array, asd_array=self.asd_array)
+        self.assertTrue(np.array_equal(self.psd_array, actual.psd_array))
+        self.assertTrue(np.array_equal(self.asd_array, actual.asd_array))
+
+    def test_from_power_spectral_density_array(self):
+        actual = tupak.gw.detector.PowerSpectralDensity.from_power_spectral_density_array(
+            frequency_array=self.frequency_array, psd_array=self.psd_array)
+        self.assertTrue(np.array_equal(self.psd_array, actual.psd_array))
+        self.assertTrue(np.array_equal(self.asd_array, actual.asd_array))
+
+    def test_repr(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, psd_array=self.psd_array)
+        expected = 'PowerSpectralDensity(frequency_array={}, psd_array={}, asd_array={})'.format(self.frequency_array,
+                                                                                                 self.psd_array,
+                                                                                                 self.asd_array)
+        self.assertEqual(expected, repr(psd))
+
+
+class TestPowerSpectralDensityWithFiles(unittest.TestCase):
+
+    def setUp(self):
+        self.dir = os.path.join(os.path.dirname(__file__), 'noise_curves')
+        os.mkdir(self.dir)
+        self.asd_file = os.path.join(os.path.dirname(__file__), 'noise_curves', 'asd_test_file.txt')
+        self.psd_file = os.path.join(os.path.dirname(__file__), 'noise_curves', 'psd_test_file.txt')
+        with open(self.asd_file, 'w') as f:
+            f.write('1.\t1.0e-21\n2.\t2.0e-21\n3.\t3.0e-21')
+        with open(self.psd_file, 'w') as f:
+            f.write('1.\t1.0e-42\n2.\t4.0e-42\n3.\t9.0e-42')
+        self.frequency_array = np.array([1.0, 2.0, 3.0])
+        self.asd_array = np.array([1.0e-21, 2.0e-21, 3.0e-21])
+        self.psd_array = np.array([1.0e-42, 4.0e-42, 9.0e-42])
+
+    def tearDown(self):
+        os.remove(self.asd_file)
+        os.remove(self.psd_file)
+        os.rmdir(self.dir)
+        del self.dir
+        del self.asd_array
+        del self.psd_array
+        del self.asd_file
+        del self.psd_file
+
+    def test_init_with_psd_file(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, psd_file=self.psd_file)
+        self.assertEqual(self.psd_file, psd.psd_file)
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+        self.assertTrue(np.allclose(self.asd_array, psd.asd_array, atol=1e-30))
+
+    def test_init_with_asd_file(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, asd_file=self.asd_file)
+        self.assertEqual(self.asd_file, psd.asd_file)
+        self.assertTrue(np.allclose(self.psd_array, psd.psd_array, atol=1e-60))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+
+    def test_setting_psd_array_after_init(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array)
+        psd.psd_file = self.psd_file
+        self.assertEqual(self.psd_file, psd.psd_file)
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+        self.assertTrue(np.allclose(self.asd_array, psd.asd_array, atol=1e-30))
+
+    def test_init_with_asd_array_after_init(self):
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array)
+        psd.asd_file = self.asd_file
+        self.assertEqual(self.asd_file, psd.asd_file)
+        self.assertTrue(np.allclose(self.psd_array, psd.psd_array, atol=1e-60))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+
+    def test_power_spectral_density_interpolated_from_asd_file(self):
+        expected = np.array([4.0e-42])
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, asd_file=self.asd_file)
+        self.assertTrue(np.allclose(expected, psd.power_spectral_density_interpolated(2), atol=1e-60))
+
+    def test_power_spectral_density_interpolated_from_psd_file(self):
+        expected = np.array([4.0e-42])
+        psd = tupak.gw.detector.PowerSpectralDensity(frequency_array=self.frequency_array, psd_file=self.psd_file)
+        self.assertAlmostEqual(expected, psd.power_spectral_density_interpolated(2))
+
+    def test_from_amplitude_spectral_density_file(self):
+        psd = tupak.gw.detector.PowerSpectralDensity.from_amplitude_spectral_density_file(asd_file=self.asd_file)
+        self.assertEqual(self.asd_file, psd.asd_file)
+        self.assertTrue(np.allclose(self.psd_array, psd.psd_array, atol=1e-60))
+        self.assertTrue(np.array_equal(self.asd_array, psd.asd_array))
+
+    def test_from_power_spectral_density_file(self):
+        psd = tupak.gw.detector.PowerSpectralDensity.from_power_spectral_density_file(psd_file=self.psd_file)
+        self.assertEqual(self.psd_file, psd.psd_file)
+        self.assertTrue(np.array_equal(self.psd_array, psd.psd_array))
+        self.assertTrue(np.allclose(self.asd_array, psd.asd_array, atol=1e-30))
+
+    def test_from_aligo(self):
+        psd = tupak.gw.detector.PowerSpectralDensity.from_aligo()
+        expected_filename = 'aLIGO_ZERO_DET_high_P_psd.txt'
+        expected = tupak.gw.detector.PowerSpectralDensity(psd_file=expected_filename)
+        actual_filename = psd.psd_file.split('/')[-1]
+        self.assertEqual(expected_filename, actual_filename)
+        self.assertTrue(np.allclose(expected.psd_array, psd.psd_array, atol=1e-60))
+        self.assertTrue(np.array_equal(expected.asd_array, psd.asd_array))
+
+    def test_check_file_psd_file_set_to_asd_file(self):
+        logger = logging.getLogger('tupak')
+        m = MagicMock()
+        logger.warning = m
+        psd = tupak.gw.detector.PowerSpectralDensity(psd_file=self.asd_file)
+        self.assertEqual(4, m.call_count)
+
+    def test_check_file_not_called_psd_file_set_to_psd_file(self):
+        logger = logging.getLogger('tupak')
+        m = MagicMock()
+        logger.warning = m