Merge branch 'port-to-python3-and-new-lalsuite' into 'master'

Port to python3 and new lalsuite

See merge request !19

examples/MCMC_examples/fully_coherent_search_using_MCMC.py

@@ -28,13 +28,13 @@ data.make_data()
 
 # The predicted twoF, given by lalapps_predictFstat can be accessed by
 twoF = data.predict_fstat()
-print 'Predicted twoF value: {}\n'.format(twoF)
+print('Predicted twoF value: {}\n'.format(twoF))
 
 DeltaF0 = 1e-7
 DeltaF1 = 1e-13
 VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
 VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
-print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
+print('\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1))
 
 theta_prior = {'F0': {'type': 'unif',
                       'lower': F0-DeltaF0/2.,

examples/MCMC_examples/semi_coherent_search_using_MCMC.py

@@ -28,13 +28,13 @@ data.make_data()
 
 # The predicted twoF, given by lalapps_predictFstat can be accessed by
 twoF = data.predict_fstat()
-print 'Predicted twoF value: {}\n'.format(twoF)
+print('Predicted twoF value: {}\n'.format(twoF))
 
 DeltaF0 = 1e-7
 DeltaF1 = 1e-13
 VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
 VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
-print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
+print('\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1))
 
 theta_prior = {'F0': {'type': 'unif',
                       'lower': F0-DeltaF0/2.,

examples/followup_examples/semi_coherent_directed_follow_up.py

@@ -28,7 +28,7 @@ data.make_data()
 
 # The predicted twoF, given by lalapps_predictFstat can be accessed by
 twoF = data.predict_fstat()
-print 'Predicted twoF value: {}\n'.format(twoF)
+print('Predicted twoF value: {}\n'.format(twoF))
 
 # Search
 VF0 = VF1 = 1e5

examples/glitch_examples/semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py

@@ -66,5 +66,5 @@ mcmc.plot_corner(label_offset=0.25, truths=[0, 0, 0, 0],
 
 mcmc.print_summary()
 
-print('Prior widths =', F0_width, F1_width)
-print("Actual run time = {}".format(dT))
+print(('Prior widths =', F0_width, F1_width))
+print(("Actual run time = {}".format(dT)))

examples/glitch_examples/semicoherent_glitch_robust_directed_grid_search_on_1_glitch.py

@@ -62,6 +62,6 @@ fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label),
             bbox_inches='tight')
 
 
-print('Prior widths =', F0_width, F1_width)
-print("Actual run time = {}".format(dT))
-print("Actual number of grid points = {}".format(search.data.shape[0]))
+print(('Prior widths =', F0_width, F1_width))
+print(("Actual run time = {}".format(dT)))
+print(("Actual number of grid points = {}".format(search.data.shape[0])))

examples/other_examples/twoF_cumulative.py

@@ -27,13 +27,13 @@ data.make_data()
 
 # The predicted twoF, given by lalapps_predictFstat can be accessed by
 twoF = data.predict_fstat()
-print 'Predicted twoF value: {}\n'.format(twoF)
+print('Predicted twoF value: {}\n'.format(twoF))
 
 DeltaF0 = 1e-7
 DeltaF1 = 1e-13
 VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
 VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
-print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
+print('\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1))
 
 theta_prior = {'F0': {'type': 'unif',
                       'lower': F0-DeltaF0/2.,

examples/other_examples/using_initialisation.py

@@ -28,13 +28,13 @@ data.make_data()
 
 # The predicted twoF, given by lalapps_predictFstat can be accessed by
 twoF = data.predict_fstat()
-print 'Predicted twoF value: {}\n'.format(twoF)
+print('Predicted twoF value: {}\n'.format(twoF))
 
 DeltaF0 = 1e-7
 DeltaF1 = 1e-13
 VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
 VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
-print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
+print('\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1))
 
 theta_prior = {'F0': {'type': 'unif',
                       'lower': F0-DeltaF0/2.,

pyfstat/__init__.py

-from __future__ import division as _division
+
 
 from .core import BaseSearchClass, ComputeFstat, SemiCoherentSearch, SemiCoherentGlitchSearch
 from .make_sfts import Writer, GlitchWriter, FrequencyModulatedArtifactWriter, FrequencyAmplitudeModulatedArtifactWriter

pyfstat/core.py

 """ The core tools used in pyfstat """
-from __future__ import division, absolute_import, print_function
+
 
 import os
 import logging
@@ -530,8 +530,17 @@ class ComputeFstat(BaseSearchClass):
                 self.injectSources))
             PPV = lalpulsar.CreatePulsarParamsVector(1)
             PP = PPV.data[0]
-            PP.Amp.h0 = self.injectSources['h0']
-            PP.Amp.cosi = self.injectSources['cosi']
+            h0 = self.injectSources['h0']
+            cosi = self.injectSources['cosi']
+            use_aPlus = ('aPlus' in dir(PP.Amp))
+            print("use_aPlus = {}".format(use_aPlus))
+            if use_aPlus:  # lalsuite interface changed in aff93c45
+                PP.Amp.aPlus = 0.5 * h0 * (1.0 + cosi**2)
+                PP.Amp.aCross = h0 * cosi
+            else:
+                PP.Amp.h0 = h0
+                PP.Amp.cosi = cosi
+
             PP.Amp.phi0 = self.injectSources['phi0']
             PP.Amp.psi = self.injectSources['psi']
             PP.Doppler.Alpha = self.injectSources['Alpha']
@@ -681,14 +690,14 @@ class ComputeFstat(BaseSearchClass):
                                argp=None):
         """ Returns twoF or ln(BSGL) fully-coherently at a single point """
         self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0])
-        self.PulsarDopplerParams.Alpha = Alpha
-        self.PulsarDopplerParams.Delta = Delta
+        self.PulsarDopplerParams.Alpha = float(Alpha)
+        self.PulsarDopplerParams.Delta = float(Delta)
         if self.binary:
-            self.PulsarDopplerParams.asini = asini
-            self.PulsarDopplerParams.period = period
-            self.PulsarDopplerParams.ecc = ecc
-            self.PulsarDopplerParams.tp = tp
-            self.PulsarDopplerParams.argp = argp
+            self.PulsarDopplerParams.asini = float(asini)
+            self.PulsarDopplerParams.period = float(period)
+            self.PulsarDopplerParams.ecc = float(ecc)
+            self.PulsarDopplerParams.tp = float(tp)
+            self.PulsarDopplerParams.argp = float(argp)
 
         lalpulsar.ComputeFstat(self.FstatResults,
                                self.FstatInput,
@@ -1026,14 +1035,14 @@ class SemiCoherentSearch(ComputeFstat):
         """ Returns twoF or ln(BSGL) semi-coherently at a single point """
 
         self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0])
-        self.PulsarDopplerParams.Alpha = Alpha
-        self.PulsarDopplerParams.Delta = Delta
+        self.PulsarDopplerParams.Alpha = float(Alpha)
+        self.PulsarDopplerParams.Delta = float(Delta)
         if self.binary:
-            self.PulsarDopplerParams.asini = asini
-            self.PulsarDopplerParams.period = period
-            self.PulsarDopplerParams.ecc = ecc
-            self.PulsarDopplerParams.tp = tp
-            self.PulsarDopplerParams.argp = argp
+            self.PulsarDopplerParams.asini = float(asini)
+            self.PulsarDopplerParams.period = float(period)
+            self.PulsarDopplerParams.ecc = float(ecc)
+            self.PulsarDopplerParams.tp = float(tp)
+            self.PulsarDopplerParams.argp = float(argp)
 
         lalpulsar.ComputeFstat(self.FstatResults,
                                self.FstatInput,

pyfstat/grid_based_searches.py

 """ Searches using grid-based methods """
-from __future__ import division, absolute_import, print_function
+
 
 import os
 import logging
@@ -12,7 +12,7 @@ import socket
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
-from scipy.misc import logsumexp
+from scipy.special import logsumexp
 
 import pyfstat.helper_functions as helper_functions
 from pyfstat.core import (BaseSearchClass, ComputeFstat,
@@ -342,7 +342,7 @@ class GridSearch(BaseSearchClass):
     def print_max_twoF(self):
         d = self.get_max_twoF()
         print('Max twoF values for {}:'.format(self.label))
-        for k, v in d.iteritems():
+        for k, v in d.items():
             print('  {}={}'.format(k, v))
 
     def set_out_file(self, extra_label=None):
@@ -1006,7 +1006,7 @@ class EarthTest(GridSearch):
         vals = [self.minStartTime, self.maxStartTime, self.F0, self.F1,
                 self.F2, self.Alpha, self.Delta]
         self.special_data = {'zero': [0, 0, 0]}
-        for key, (dR, dphi, dP) in self.special_data.iteritems():
+        for key, (dR, dphi, dP) in self.special_data.items():
             rescaleRadius = (1 + dR / lal.REARTH_SI)
             rescalePeriod = (1 + dP / lal.DAYSID_SI)
             lalpulsar.BarycenterModifyEarthRotation(

pyfstat/helper_functions.py

@@ -110,7 +110,7 @@ def get_ephemeris_files():
             logging.warning('No [earth/sun]_ephem found in '+config_file+'. '+please)
             earth_ephem = None
             sun_ephem = None
-    elif env_var in os.environ.keys():
+    elif env_var in list(os.environ.keys()):
         earth_ephem = os.path.join(os.environ[env_var],'earth00-40-DE421.dat.gz')
         sun_ephem = os.path.join(os.environ[env_var],'sun00-40-DE421.dat.gz')
         if not ( os.path.isfile(earth_ephem) and os.path.isfile(sun_ephem) ):

pyfstat/make_sfts.py

 """ pyfstat tools to generate sfts """
-from __future__ import division, absolute_import, print_function
+
 
 import numpy as np
 import logging
@@ -120,10 +120,10 @@ refTime = {:10.6f}""")
         template = (self.get_base_template(
             i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, F2, tref) + """
 transientWindowType = {:s}
-transientStartTime = {:10.3f}
-transientTauDays = {:1.3f}\n""")
+transientStartTime = {:10.0f}
+transientTau = {:10.0f}\n""")
         return template.format(i, Alpha, Delta, h0, cosi, psi, phi, F0, F1,
-                               F2, tref, window, tstart, duration_days)
+                               F2, tref, window, tstart, duration_days*86400)
 
     def get_single_config_line(self, i, Alpha, Delta, h0, cosi, psi, phi, F0,
                                F1, F2, tref, window, tstart, duration_days):
@@ -477,7 +477,7 @@ class FrequencyModulatedArtifactWriter(Writer):
 
         linePhi = 0
         lineFreq_old = 0
-        for i in tqdm(range(self.nsfts)):
+        for i in tqdm(list(range(self.nsfts))):
             mid_time = self.tstart + (i+.5)*self.Tsft
             lineFreq = self.get_frequency(mid_time)
 
@@ -517,7 +517,7 @@ class FrequencyModulatedArtifactWriter(Writer):
             logging.info('Using {} threads'.format(args.N))
             try:
                 with pathos.pools.ProcessPool(args.N) as p:
-                    list(tqdm(p.imap(self.make_ith_sft, range(self.nsfts)),
+                    list(tqdm(p.imap(self.make_ith_sft, list(range(self.nsfts))),
                               total=self.nsfts))
             except KeyboardInterrupt:
                 p.terminate()
@@ -525,7 +525,7 @@ class FrequencyModulatedArtifactWriter(Writer):
             logging.info(
                 "No multiprocessing requested or `pathos` not install, cont."
                 " without multiprocessing")
-            for i in tqdm(range(self.nsfts)):
+            for i in tqdm(list(range(self.nsfts))):
                 self.make_ith_sft(i)
 
         self.concatenate_sft_files()

pyfstat/mcmc_based_searches.py

 """ Searches using MCMC-based methods """
-from __future__ import division, absolute_import, print_function
+
 
 import sys
 import os
@@ -221,7 +221,7 @@ class MCMCSearch(core.BaseSearchClass):
 
         self.theta_keys = []
         fixed_theta_dict = {}
-        for key, val in self.theta_prior.iteritems():
+        for key, val in self.theta_prior.items():
             if type(val) is dict:
                 fixed_theta_dict[key] = 0
                 self.theta_keys.append(key)
@@ -953,7 +953,7 @@ class MCMCSearch(core.BaseSearchClass):
         See the pyfstat.core.plot_twoF_cumulative function for further details
         """
         d, maxtwoF = self.get_max_twoF()
-        for key, val in self.theta_prior.iteritems():
+        for key, val in self.theta_prior.items():
             if key not in d:
                 d[key] = val
 
@@ -1223,8 +1223,8 @@ class MCMCSearch(core.BaseSearchClass):
     def _generate_scattered_p0(self, p):
         """ Generate a set of p0s scattered about p """
         p0 = [[p + self.scatter_val * p * np.random.randn(self.ndim)
-               for i in xrange(self.nwalkers)]
-              for j in xrange(self.ntemps)]
+               for i in range(self.nwalkers)]
+              for j in range(self.ntemps)]
         return p0
 
     def _generate_initial_p0(self):
@@ -1349,7 +1349,7 @@ class MCMCSearch(core.BaseSearchClass):
                 setattr(self, key, new_d[key])
 
         mod_keys = []
-        for key in new_d.keys():
+        for key in list(new_d.keys()):
             if key in old_d:
                 if new_d[key] != old_d[key]:
                     mod_keys.append((key, old_d[key], new_d[key]))
@@ -1486,10 +1486,10 @@ class MCMCSearch(core.BaseSearchClass):
             if hasattr(self, 'theta0_index'):
                 f.write('theta0_index = {}\n'.format(self.theta0_idx))
             if method == 'med':
-                for key, val in median_std_d.iteritems():
+                for key, val in median_std_d.items():
                     f.write('{} = {:1.16e}\n'.format(key, val))
             if method == 'twoFmax':
-                for key, val in max_twoF_d.iteritems():
+                for key, val in max_twoF_d.items():
                     f.write('{} = {:1.16e}\n'.format(key, val))
 
     def generate_loudest(self):
@@ -1514,7 +1514,7 @@ class MCMCSearch(core.BaseSearchClass):
             f.write(r"\begin{tabular}{c l c} \hline" + '\n'
                     r"Parameter & & &  \\ \hhline{====}")
 
-            for key, prior in self.theta_prior.iteritems():
+            for key, prior in self.theta_prior.items():
                 if type(prior) is dict:
                     Type = prior['type']
                     if Type == "unif":
@@ -1546,10 +1546,10 @@ class MCMCSearch(core.BaseSearchClass):
         if hasattr(self, 'theta0_idx'):
             logging.info('theta0 index: {}'.format(self.theta0_idx))
         logging.info('Max twoF: {} with parameters:'.format(max_twoF))
-        for k in np.sort(max_twoFd.keys()):
+        for k in np.sort(list(max_twoFd.keys())):
             print('  {:10s} = {:1.9e}'.format(k, max_twoFd[k]))
         logging.info('Median +/- std for production values')
-        for k in np.sort(median_std_d.keys()):
+        for k in np.sort(list(median_std_d.keys())):
             if 'std' not in k:
                 logging.info('  {:10s} = {:1.9e} +/- {:1.9e}'.format(
                     k, median_std_d[k], median_std_d[k+'_std']))
@@ -1668,7 +1668,7 @@ class MCMCSearch(core.BaseSearchClass):
 
     def write_evidence_file_from_dict(self, EvidenceDict, evidence_file_name):
         with open(evidence_file_name, 'w+') as f:
-            for key, val in EvidenceDict.iteritems():
+            for key, val in EvidenceDict.items():
                 f.write('{} {} {}\n'.format(key, val[0], val[1]))
 
 
@@ -1801,7 +1801,7 @@ class MCMCGlitchSearch(MCMCSearch):
                                r'$\delta$'] + full_glitch_symbols)
         self.theta_keys = []
         fixed_theta_dict = {}
-        for key, val in self.theta_prior.iteritems():
+        for key, val in self.theta_prior.items():
             if type(val) is dict:
                 fixed_theta_dict[key] = 0
                 if key in glitch_keys:
@@ -1863,7 +1863,7 @@ class MCMCGlitchSearch(MCMCSearch):
 
         fig, ax = plt.subplots()
         d, maxtwoF = self.get_max_twoF()
-        for key, val in self.theta_prior.iteritems():
+        for key, val in self.theta_prior.items():
             if key not in d:
                 d[key] = val
 
@@ -2223,7 +2223,7 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
 
     def check_old_run_setup(self, old_setup, **kwargs):
         try:
-            truths = [val == old_setup[key] for key, val in kwargs.iteritems()]
+            truths = [val == old_setup[key] for key, val in kwargs.items()]
             if all(truths):
                 return True
             else:
@@ -2540,7 +2540,7 @@ class MCMCTransientSearch(MCMCSearch):
 
         self.theta_keys = []
         fixed_theta_dict = {}
-        for key, val in self.theta_prior.iteritems():
+        for key, val in self.theta_prior.items():
             if type(val) is dict:
                 fixed_theta_dict[key] = 0
                 self.theta_keys.append(key)

pyfstat/optimal_setup_functions.py

@@ -3,7 +3,7 @@
 Provides functions to aid in calculating the optimal setup for zoom follow up
 
 """
-from __future__ import division, absolute_import, print_function
+
 
 import logging
 import numpy as np

pyfstat/tcw_fstat_map_funcs.py

@@ -31,7 +31,7 @@ def _optional_import ( modulename, shorthand=None ):
         logging.debug('Successfully imported module %s%s.'
                       % (modulename, shorthandbit))
         success = True
-    except ImportError, e:
+    except ImportError as e:
         if e.message == 'No module named '+modulename:
             logging.debug('No module {:s} found.'.format(modulename))
             success = False
@@ -111,7 +111,7 @@ def init_transient_fstat_map_features ( wantCuda=False, cudaDeviceName=None ):
                       ' then checking all available devices...')
         try:
             context0 = pycuda.tools.make_default_context()
-        except pycuda._driver.LogicError, e:
+        except pycuda._driver.LogicError as e:
             if e.message == 'cuDeviceGet failed: invalid device ordinal':
                 devn = int(os.environ['CUDA_DEVICE'])
                 raise RuntimeError('Requested CUDA device number {} exceeds' \

setup.py

 #!/usr/bin/env python
 
-from distutils.core import setup
+from setuptools import setup, find_packages
+from os import path
+
+here = path.abspath(path.dirname(__file__))
+# Get the long description from the README file
+with open(path.join(here, 'README.md'), encoding='utf-8') as f:
+    long_description = f.read()
 
 setup(name='PyFstat',
       version='0.2',
       author='Gregory Ashton',
       author_email='gregory.ashton@ligo.org',
-      packages=['pyfstat'],
+      packages=find_packages(where="pyfstat"),
       include_package_data=True,
       package_data={'pyfstat': ['pyCUDAkernels/cudaTransientFstatExpWindow.cu',
                                 'pyCUDAkernels/cudaTransientFstatRectWindow.cu']},
-      )
+      install_requires=[
+          'matplotlib',
+          'scipy',
+          'ptemcee',
+          'corner',
+          'dill',
+          'tqdm',
+          'bashplotlib',
+          'peakutils',
+          'pathos',
+          'pycuda',
+      ],
+)

tests.py

@@ -93,7 +93,6 @@ class par(Test):
     label = 'TestPar'
 
     def test(self):
-        os.system('mkdir {}'.format(self.outdir))
         os.system(
             'echo "x=100\ny=10" > {}/{}.par'.format(self.outdir, self.label))
 
@@ -315,10 +314,10 @@ class SemiCoherentGlitchSearch(Test):
         Writer.tend = maxStartTime
         FSB = Writer.predict_fstat()
 
-        print FSA, FSB
+        print(FSA, FSB)
         predicted_FS = (FSA + FSB)
 
-        print(predicted_FS, FS)
+        print((predicted_FS, FS))
         self.assertTrue(np.abs((FS - predicted_FS))/predicted_FS < 0.3)
 
 
@@ -359,8 +358,8 @@ class MCMCSearch(Test):
         search.run(create_plots=False)
         _, FS = search.get_max_twoF()
 
-        print('Predicted twoF is {} while recovered is {}'.format(
-                predicted_FS, FS))
+        print(('Predicted twoF is {} while recovered is {}'.format(
+                predicted_FS, FS)))
         self.assertTrue(
             FS > predicted_FS or np.abs((FS-predicted_FS))/predicted_FS < 0.3)