Emri

nrutils/__init__.py

@@ -33,6 +33,7 @@
 '''
 
 # # Import all python files within this folder
+from __future__ import print_function
 from os.path import dirname, basename, isdir, realpath
 from commands import getstatusoutput as bash
 
@@ -45,7 +46,7 @@
 
 # Search recurssively within the config's sim_dir for files matching the config's metadata_id
 this_file = realpath(__file__)
-if verbose: print "The highest level init for nrutils is located at: %s" % this_file
+if verbose: print( "The highest level init for nrutils is located at: %s" % this_file)
 if this_file[-1] == 'c': this_file = this_file[:-1]
 cmd = 'find %s -maxdepth 2 -name "__init__.py"' % dirname(this_file)
 status, output = bash(cmd)
@@ -72,24 +73,24 @@
 
 
 # Import all modules from each package
-if verbose: print '\n>> Initiating nrutils ...'
+if verbose: print( '\n>> Initiating nrutils ...')
 
 # Let the people know
 if verbose:
-    print "\n>> Sub-Packages to be imported:"
+    print( "\n>> Sub-Packages to be imported:")
     for k in internal_packages:
-        print '   -> %s' % k
+        print( '   -> %s' % k)
 
 # Some other notes
-if verbose: print '>> Please note style conventions:\
+if verbose: print( '>> Please note style conventions:\
                   \n   * lower case function/method/variable names\
                   \n   * no underscore in names unless there are repeated letters, or counfounded syllables\
                   \n   * information is implicitely in time domain unless explicitely stated.\
-                  \n   * frequency domain information will start with "fd".\n'
+                  \n   * frequency domain information will start with "fd".\n')
 
-if verbose: print '%s:\n' % __name__
+if verbose: print( '%s:\n' % __name__)
 for p in internal_packages:
-    if verbose: print '  .%s: ' % p
+    if verbose: print( '  .%s: ' % p)
     exec r'import %s' % p
     # exec 'from %s import *' % p
 
@@ -100,13 +101,8 @@
 from nrutils.formula import *
 
 #
-if verbose: print ''
+if verbose: print( '')
 # Cleanup
 del cmd, bash, p, dir_list, status, output, this_file, basename, dirname, isdir, realpath
 
-# for d in modules:
-#     print '\n|&) Importing: %s' % d
-#     # __import__(d,locals(),globals())
-#     # cmd = 'from %s import *'
-#     cmd = 'from %s import *' % d
-#     exec cmd
+

nrutils/config/bradwr.ini

@@ -1,4 +1,4 @@
-
+    
 # Superficial label for the institute
 institute = gt
 

nrutils/config/emr.ini

@@ -0,0 +1,37 @@
+
+#
+institute = umass
+
+#
+metadata_id = .md
+
+# Directory where all smulation folders are stored. The assumtion is that each simulation folder contains a related metadata file with a filname which contains the metadata_id
+catalog_dir = /Volumes/lucidum/emr-waveforms-tmp/
+; catalog_dir = /Volumes/lucidum/emr-equatorial-waveforms-tmp/
+
+# NOTE standard parameter order for every simulation catalog
+# extraction_parameter, l, m, level
+# Also NOTE that only default extraction_parameter and level are to be specified in this file. The multipolar indeces l andm must be input at waveform loadeding via instances of the gwylm class.
+
+# NOTE that the ordering here is: extraction_parameter, l, m, level
+# NOTE that information on the string formatting here can be found at https://pyformat.info/
+data_file_name_format = rStrainByMass2_l{1:d}_m{2:d}_rInf.gz
+
+
+# NOTE that the ordering here is: default_extraction_parameter, default_level. NOTE that not all groups use level as a sepcifier and that not all levels exist for each single simulation.
+default_par_list = 0,0
+
+#
+is_extrapolated = 1
+
+# Toggle for whether waveform is scaled by extration radius or not
+is_rscaled = 1
+
+# Data type
+__disk_data_kind__ = strain
+
+# Tag for whether set is of the extreme massratio type and waveforms are scaled by smaller mass
+__is_extreme_mass_ratio__ = 1
+
+# Name of handler file within install directory: .../nrutils/handlers/[handler_name].py
+handler_name = emr

nrutils/config/hogshead.ini

@@ -6,8 +6,8 @@ institute = bam
 metadata_id = .bbh
 
 # Directory where all smulation folders are stored. The assumtion is that each simulation folder contains a related metadata file with a filname which contains the metadata_id
-#catalog_dir = ~/MOUNT/tmp/gluster/hogshead/NR_data/ReducedData-hogshead/dirac-nr-followup/
-catalog_dir = /Users/jthompson/hmnt/hogshead/NR_data/ReducedData-hogshead/,/Users/shrobanaghosh/gluster/NR_data/ReducedData-hogshead/,/gluster/hogshead/NR_data/ReducedData-hogshead/, ~/MOUNT/tmp/gluster/hogshead/NR_data/ReducedData-hogshead/,/Users/markodh/sshfs-atlas2/gluster/hogshead/NR_data/ReducedData-hogshead/, /Users/Eleanor/Documents/Jupyter/Precession/Precession_angles/minion_data/,/Users/eleanorhamilton/mnt/gw-data/NR_data/,/Users/alex/mount/minion/NR_data/ReducedData-hogshead/
+#catalog_dir = ~/MOUNT/tmp/gluster/gw-data/NR_data/ReducedData-hogshead/dirac-nr-followup/
+catalog_dir = /Users/jthompson/hmnt/hogshead/NR_data/ReducedData-hogshead/,/Users/shrobanaghosh/gluster/NR_data/ReducedData-hogshead/,/gluster/gw-data/NR_data/ReducedData-hogshead/, ~/MOUNT/tmp/gluster/gw-data/NR_data/ReducedData-hogshead/,/Users/markodh/sshfs-atlas2/gluster/gw-data/NR_data/ReducedData-hogshead/, /Users/Eleanor/Documents/Jupyter/Precession/Precession_angles/minion_data/,/Users/eleanorhamilton/mnt/gw-data/NR_data/,/Users/alex/mount/minion/NR_data/ReducedData-hogshead/
 
 # NOTE standard parameter order for every simulation catalog
 # extraction_parameter, l, m, level

nrutils/config/rit.ini

@@ -0,0 +1,36 @@
+
+#
+institute = rit
+
+#
+metadata_id = Metadata
+
+# Directory where all smulation folders are stored. The assumtion is that each simulation folder contains a related metadata file with a filname which contains the metadata_id
+catalog_dir = /Volumes/lucidum/rit-catalog/Data/
+
+# NOTE standard parameter order for every simulation catalog
+# extraction_parameter, l, m, level
+# Also NOTE that only default extraction_parameter and level are to be specified in this file. The multipolar indeces l andm must be input at waveform loadeding via instances of the gwylm class.
+
+# NOTE that the ordering here is: extraction_parameter, l, m, level
+# NOTE that information on the string formatting here can be found at https://pyformat.info/
+data_file_name_format = rPsi4_l{1:d}_m{2:d}_rInf.asc
+
+
+# NOTE that the ordering here is: default_extraction_parameter, default_level. NOTE that not all groups use level as a sepcifier and that not all levels exist for each single simulation.
+default_par_list = 0,0
+
+#
+is_extrapolated = 1
+
+# Toggle for whether waveform is scaled by extration radius or not
+is_rscaled = 1
+
+# Data type
+__disk_data_kind__ = psi4
+
+# Tag for whether set is of the extreme massratio type and waveforms are scaled by smaller mass
+__is_extreme_mass_ratio__ = 0
+
+# Name of handler file within install directory: .../nrutils/handlers/[handler_name].py
+handler_name = rit

nrutils/core/basics.py

@@ -149,7 +149,7 @@ def intrp_wfarr(wfarr,delta=None,domain=None,verbose = False):
 
 
 # Fucntion to pad wfarr with zeros. NOTE that this should only be applied to a time domain waveform that already begins and ends with zeros.
-def pad_wfarr(wfarr,new_length,where=None,verbose=None,extend=True):
+def pad_wfarr(wfarr,new_length,where=None,verbose=None,extend=True,__nowarn__=False):
 
     #
     from numpy import hstack,zeros,arange,pad,unwrap,angle,cos,sin
@@ -171,7 +171,7 @@ def pad_wfarr(wfarr,new_length,where=None,verbose=None,extend=True):
     # Warn the user if extend is false
     if not extend:
         msg = 'You have disabled the extend option. As a result the input padding length will be interpreted as the desired total length of the new waveform array. This course is discouraged in favor of e.g. using the fftlength option when taking fouorier transforms, OR simply inputting the desired pad amount.'
-        warning(msg)
+        warning(msg, verbose=(not __nowarn__) and verbose)
 
     #
     if where is None:
@@ -400,7 +400,7 @@ def recompose_wfarrs( wfarr_dict, theta, phi ):
                 error( 'invalid multipole eigenvalue found: %s'%[v] )
         # key values must be ndarray
         if not isinstance(wfarr_dict[k],ndarray):
-            print wfarr_dict[k].__class__
+            alert(wfarr_dict[k].__class__)
             error('key values must be ndarray')
 
     # Number of samples
@@ -518,6 +518,13 @@ def straighten_wfarr( wfarr, verbose=False ):
     # Import useful things
     from numpy import arange,sum,array,diff,isfinite,hstack,allclose,median
     thisfun = 'straighten_wfarr'
+
+    #
+    rw,cl = wfarr.shape 
+    if rw == 0:
+        error('input waveform array has zero rows')
+    if cl == 0:
+        error('input waveform array has zero columns')
 
     # check whether t is monotonically increasing
     isincreasing = allclose( wfarr[:,0], sorted(wfarr[:,0]), 1e-6 )
@@ -529,28 +536,28 @@ def straighten_wfarr( wfarr, verbose=False ):
         map_ = arange( len(wfarr[:,0]) )
         map_ = sorted( map_, key = lambda x: wfarr[x,0] )
         wfarr = wfarr[ map_, : ]
-        if allclose( wfarr[:,0], sorted(wfarr[:,0]), 1e-6 ) and verbose: warning(red('The waveform time series is now monotonic.'))
+        # if allclose( wfarr[:,0], sorted(wfarr[:,0]), 1e-6 ) and verbose: warning(red('The waveform time series is now monotonic.'))
 
     # Remove rows that contain non-finite data
     finite_mask = isfinite( sum( wfarr, 1 ) )
-    if sum(finite_mask)!=len(finite_mask):
-        if verbose: warning('Non-finite values found in waveform array. Corresponding rows will be removed.',thisfun)
+    # if sum(finite_mask)!=len(finite_mask):
+    #     if verbose: warning('Non-finite values found in waveform array. Corresponding rows will be removed.',thisfun)
     wfarr = wfarr[ finite_mask, : ]
 
     # Sort rows by the time series' values
     time = array( wfarr[:,0] )
     space = arange( wfarr.shape[0] )
     chart = sorted( space, key = lambda k: time[k] )
-    if (space != chart).all():
-        if verbose: warning('The waveform array was found to have nonmonotinicities in its time series. The array will now be straightened.',thisfun)
+    # if (space != chart).all():
+    #     if verbose: warning('The waveform array was found to have nonmonotinicities in its time series. The array will now be straightened.',thisfun)
     wfarr = wfarr[ chart, : ]
 
     # Remove rows with duplicate time series values
     time = array( wfarr[:,0] )
     dt = median( diff(time) )
     diff_mask = hstack( [ True, diff(time)/dt>1e-6 ] )
-    if sum(diff_mask)!=len(diff_mask):
-        if verbose: warning('Repeated time values were found in the array. Offending rows will be removed.',thisfun)
+    # if sum(diff_mask)!=len(diff_mask):
+    #     if verbose: warning('Repeated time values were found in the array. Offending rows will be removed.',thisfun)
     wfarr = wfarr[ diff_mask, : ]
 
     # The wfarr should now be straight
@@ -604,7 +611,7 @@ def alphamax(_gwylmo,dphi,plt=False,verbose=False,n=13):
         plot( dpsis, aspl, label=dpsi )
         plot( dpsis[argmax(aspl)], aspl[argmax(aspl)], 'or', mfc='none' )
         xlabel(r'$\psi$')
-    if verbose: print dpsi_opt,action([dphi,dpsi_opt])
+    if verbose: alert(dpsi_opt,action([dphi,dpsi_opt]))
     return [ dpsi_opt, action([dphi,dpsi_opt])    ]
 
 def betamax(_gwylmo,n=10,plt=False,opt=True,verbose=False):
@@ -651,10 +658,10 @@ def betamax(_gwylmo,n=10,plt=False,opt=True,verbose=False):
         xlabel(r'$\phi$')
         title(val_max)
 
-    if verbose:
-        print 'dphi_opt = ' + str(dphi_opt)
-        print 'dpsi_opt = ' + str(dpsi_opt)
-        print 'val_max = ' + str(val_max)
+    # if verbose:
+    #     print 'dphi_opt = ' + str(dphi_opt)
+    #     print 'dpsi_opt = ' + str(dpsi_opt)
+    #     print 'val_max = ' + str(val_max)
 
     return dphi_opt,dpsi_opt
 
@@ -672,7 +679,7 @@ def betamax2(_gwylmo,n=10,plt=False,opt=True,verbose=False):
         dphi_action = lambda _dphi: action( [_dphi,dpsi] )
         dphi = minimize( dphi_action, dphi, bounds=[(0,2*pi)] ).x[0]
         done = k>n
-        print '>> ',dphi,dpsi,action([dphi,dpsi])
+        alert(dphi,dpsi,action([dphi,dpsi]))
         k+=1
 
     return dphi,dpsi
@@ -848,6 +855,21 @@ def isgwylm( obj ):
     else:
         return False
 
+
+#
+def kind2texlabel(kind):
+    '''
+    Convert __disc_data_kind__ to latex for internal labeling of plots. See methods in gwylm class
+    '''
+    if kind=='psi4':
+        return '\psi'
+    elif kind=='strain':
+        return 'h'
+    elif kind=='news':
+        return 'f'
+    else:
+        error('unknown __disc_data_kind__ of %s'%red(kind))
+
 #00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%00%%#
 ''' Low level functions for rotating waveforms '''
 # https://arxiv.org/pdf/1304.3176.pdf
@@ -1152,7 +1174,9 @@ def calc_coprecessing_angles(multipole_dict,       # Dict of multipoles { ... l,
 
     #################################################
     # Reflect Y according to nrutils conventions    #
-    Y *= -1                         #
+    # NOTE that this step is needed due to the      #
+    # m_relative_sign_convention                    #
+    Y *= -1                         
     #################################################
     if IS_FD:
         alert('The domain values seem evenly split between positive and negative values. Thus, we will interpret the input as corresponding to ' +
@@ -1703,7 +1727,7 @@ def rotate_wfarrs_at_all_times( l,                          # the l of the new m
 
 
 # Careful function to find peak index for BBH waveform cases
-def find_amp_peak_index( t, amp, plot = False, return_jid=False ):
+def find_amp_peak_index( t, amp, phi, plot = False, return_jid=False ):
 
     '''
     Careful function to find peak index for BBH waveform cases
@@ -1783,13 +1807,36 @@ def find_amp_peak_index( t, amp, plot = False, return_jid=False ):
         pre[half_way:] *= 0
         dt = t[1]-t[0]
         jid = argmax( pre ) + int(100/dt)
+
+        #
+        if k_amp_max==0:
+
+            #warning('the data input may not peak near merger. we will use the second derivative of the phase to estimate the location of merger in the timeseries')
+
+            mask = amp>(1e-2*max(amp))
+            d2phi = abs(spline_diff(t[mask],phi[mask],n=2))
+
+            k_amp_max = argmax(d2phi)
+            ans = k_amp_max
+            jid=0
 
         # Return answer
         ans = k_amp_max
         if return_jid: ans = (k_amp_max,jid)
         return ans
     else:
-        warning('the data input here may be flat and zero')
-        ans = 0
+
+        # print sum(amp),sum(phi),sum(amp)<1e-6
+        if sum(phi)<1e-6:
+            ans = 0
+            if return_jid: ans = (k_amp_max,0)
+            warning('the data appears to be ZERO; please make sure you understand why')
+            return ans
+
+        warning('the data input here may be flat.')
+
+        k_amp_max = 0
+        jid=0
+        ans = k_amp_max
         if return_jid: ans = (0,0)
         return ans