Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • wenxuan.jia/pygwinc
  • sean-leavey/pygwinc
  • sebastian.steinlechner/pygwinc
  • nicholas.demos/pygwinc
  • chris.whittle/pygwinc
  • raymond.robie/pygwinc
  • mateusz.bawaj/pygwinc
  • anchal.gupta/pygwinc
  • 40m/pygwinc
  • evan.hall/pygwinc
  • kevin.kuns/pygwinc
  • geoffrey-lovelace/pygwinc
  • brittany.kamai/pygwinc
  • daniel-brown/pygwinc
  • lee-mcculler/pygwinc
  • jameson.rollins/pygwinc
  • gwinc/pygwinc
17 results
Show changes
Commits on Source (15)
Hide whitespace changes
Inline Side-by-side
Showing
with 554 additions and 499 deletions
.gitattributes
View file @ b3adea0a
gwinc/test/matlab.pkl filter=lfs diff=lfs merge=lfs -text
gwinc/test/A+.pkl filter=lfs diff=lfs merge=lfs -text
gwinc/test/aLIGO.pkl filter=lfs diff=lfs merge=lfs -text
*.h5 filter=lfs diff=lfs merge=lfs -text
.gitlab-ci.yml
View file @ b3adea0a
......@@ -10,15 +10,14 @@ test:
- echo $CI_COMMIT_SHA | cut -b1-8 > gitID.txt
script:
- apt-get update -qq
- apt-get install -y -qq python3-yaml python3-scipy python3-matplotlib python3-ipython lalsimulation-python3
- apt-get install -y -qq python3-yaml python3-scipy python3-matplotlib python3-ipython lalsimulation-python3 python3-pypdf2
- git clone https://gitlab-ci-token:ci_token@git.ligo.org/gwinc/inspiral_range.git
- export PYTHONPATH=inspiral_range
- export MPLBACKEND=agg
- python3 -m gwinc aLIGO -s aLIGO.png
- python3 -m gwinc Aplus -s Aplus.png
- python3 -m gwinc Voyager -s Voyager.png
- python3 -m gwinc.test aLIGO -t 10e-6 -k Seismic -k "Substrate Thermo-Elastic" -p -s aLIGO_test.png
- python3 -m gwinc.test -t 100e-6 -k Seismic -k "Substrate Thermo-Elastic" Aplus -p -s Aplus_test.png
- for ifo in aLIGO Aplus Voyager CE1 CE2; do
- python3 -m gwinc $ifo -s $ifo.png
- done
- python3 -m gwinc.test -r gwinc_test_report.pdf
after_script:
- rm gitID.txt
cache:
......@@ -31,8 +30,9 @@ test:
- aLIGO.png
- Aplus.png
- Voyager.png
- aLIGO_test.png
- Aplus_test.png
- CE1.png
- CE2.png
- gwinc_test_report.pdf
pages:
stage: deploy
......@@ -40,11 +40,10 @@ pages:
- test
script:
- mkdir public
- mv aLIGO.png public/
- mv Aplus.png public/
- mv Voyager.png public/
- mv aLIGO_test.png public/ || true
- mv Aplus_test.png public/ || true
- for ifo in aLIGO Aplus Voyager CE1 CE2; do
- mv $ifo.png public/
- done
- mv gwinc_test_report.pdf public/ || true
- apt-get install -y -qq python3-pip python3-dev make
- pip3 install sphinx sphinx-rtd-theme
- cd docs
......
README.md
View file @ b3adea0a
......@@ -16,10 +16,10 @@ description is loaded, the noise budget can be calculated and plotted:
>>> import gwinc
>>> import numpy as np
>>> freq = np.logspace(1, 3, 1000)
>>> Budget, ifo, freq_, plot_style = gwinc.load_ifo('aLIGO')
>>> ifo = gwinc.precompIFO(freq, ifo)
>>> traces = Budget(freq, ifo=ifo).calc_trace()
>>> fig = gwinc.plot_noise(freq, traces, **plot_style)
>>> Budget = gwinc.load_budget('aLIGO')
>>> ifo = gwinc.precompIFO(freq, Budget.ifo)
>>> traces = Budget(freq, ifo=ifo).run()
>>> fig = gwinc.plot_noise(freq, traces)
>>> fig.show()
```
......
gwinc/__init__.py
View file @ b3adea0a
from .ifo import available_ifos, load_ifo
from .struct import load_struct
from __future__ import division
import os
import logging
import importlib
import numpy as np
from .ifo import available_ifos
from .struct import Struct
from .precomp import precompIFO
from .plot import plot_noise
from .io import load_hdf5, save_hdf5
def _load_module(name_or_path):
"""Load module from name or path.
Return loaded module and module path.
"""
if os.path.exists(name_or_path):
path = name_or_path.rstrip('/')
modname = os.path.splitext(os.path.basename(path))[0]
if os.path.isdir(path):
path = os.path.join(path, '__init__.py')
spec = importlib.util.spec_from_file_location(modname, path)
mod = importlib.util.module_from_spec(spec)
spec.loader.exec_module(mod)
else:
mod = importlib.import_module(name_or_path)
try:
path = mod.__path__[0]
except AttributeError:
path = mod.__file__
return mod, path
def load_budget(name_or_path):
"""Load GWINC IFO Budget by name or from path.
Named IFOs should correspond to one of the IFOs available in the
gwinc package (see gwinc.available_ifos()). If a path is provided
it should either be a budget package (directory) or module (ending
in .py), or an IFO struct definition (see gwinc.Struct).
If a budget package path is provided, and the package includes an
'ifo.yaml' file, that file will be loaded into a Struct and
assigned as an attribute to the returned Budget class.
If a bare struct is provided the base aLIGO budget definition will
be assumed.
Returns a Budget class with 'ifo', 'freq', and 'plot_style', ifo
structure, frequency array, and plot style dictionary, with the
last three being None if they are not defined in the budget.
"""
ifo = None
if os.path.exists(name_or_path):
path = name_or_path.rstrip('/')
bname, ext = os.path.splitext(os.path.basename(path))
if ext in Struct.STRUCT_EXT:
logging.info("loading struct {}...".format(path))
ifo = Struct.from_file(path)
bname = 'aLIGO'
modname = 'gwinc.ifo.aLIGO'
logging.info("loading budget {}...".format(modname))
else:
modname = path
logging.info("loading module path {}...".format(modname))
else:
if name_or_path not in available_ifos():
raise RuntimeError("Unknonw IFO '{}' (available IFOs: {}).".format(
name_or_path,
available_ifos(),
))
bname = name_or_path
modname = 'gwinc.ifo.'+name_or_path
logging.info("loading module {}...".format(modname))
mod, modpath = _load_module(modname)
Budget = getattr(mod, bname)
ifopath = os.path.join(modpath, 'ifo.yaml')
if not ifo and ifopath:
ifo = Struct.from_file(ifopath)
Budget.ifo = ifo
return Budget
def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
"""Calculate strain noise budget for a specified interferometer model.
Argument `freq` is the frequency array for which the noises will
be calculated, and `ifoin` is the IFO model (see the `load_ifo()`
function).
If `source` structure provided, so evaluates the sensitivity of
the detector to several potential gravitational wave
sources.
If `plot` is True a plot of the budget will be created.
Returns tuple of (score, noises, ifo)
"""
# assume generic aLIGO configuration
# FIXME: how do we allow adding arbitrary addtional noise sources
# from just ifo description, without having to specify full budget
Budget = load_budget('aLIGO')
ifo = precompIFO(freq, ifo, PRfixed)
traces = Budget(freq, ifo=ifo).run()
plot_style = getattr(Budget, 'plot_style', {})
# construct matgwinc-compatible noises structure
noises = {}
for name, (data, style) in traces.items():
noises[style.get('label', name)] = data
noises['Freq'] = freq
pbs = ifo.gwinc.pbs
parm = ifo.gwinc.parm
finesse = ifo.gwinc.finesse
prfactor = ifo.gwinc.prfactor
if ifo.Laser.Power * prfactor != pbs:
pass
#warning(['Thermal lensing limits input power to ' num2str(pbs/prfactor, 3) ' W']);
# report astrophysical scores if so desired
score = None
if source:
logging.warning("Source score calculation currently not supported. See `inspiral-range` package for similar functionality:")
logging.warning("https://git.ligo.org/gwinc/inspiral-range")
# score = int731(freq, noises['Total'], ifo, source)
# score.Omega = intStoch(freq, noises['Total'], 0, ifo, source)
# --------------------------------------------------------
# output graphics
if plot:
logging.info('Laser Power: %7.2f Watt' % ifo.Laser.Power)
logging.info('SRM Detuning: %7.2f degree' % (ifo.Optics.SRM.Tunephase*180/pi))
logging.info('SRM transmission: %9.4f' % ifo.Optics.SRM.Transmittance)
logging.info('ITM transmission: %9.4f' % ifo.Optics.ITM.Transmittance)
logging.info('PRM transmission: %9.4f' % ifo.Optics.PRM.Transmittance)
logging.info('Finesse: %7.2f' % finesse)
logging.info('Power Recycling Gain: %7.2f' % prfactor)
logging.info('Arm Power: %7.2f kW' % (parm/1000))
logging.info('Power on BS: %7.2f W' % pbs)
# coating and substrate thermal load on the ITM
PowAbsITM = (pbs/2) * \
np.hstack([(finesse*2/np.pi) * ifo.Optics.ITM.CoatingAbsorption,
(2 * ifo.Materials.MassThickness) * ifo.Optics.ITM.SubstrateAbsorption])
logging.info('Thermal load on ITM: %8.3f W' % sum(PowAbsITM))
logging.info('Thermal load on BS: %8.3f W' %
(ifo.Materials.MassThickness*ifo.Optics.SubstrateAbsorption*pbs))
if (ifo.Laser.Power*prfactor != pbs):
logging.info('Lensing limited input power: %7.2f W' % (pbs/prfactor))
if score:
logging.info('BNS Inspiral Range: ' + str(score.effr0ns) + ' Mpc/ z = ' + str(score.zHorizonNS))
logging.info('BBH Inspiral Range: ' + str(score.effr0bh) + ' Mpc/ z = ' + str(score.zHorizonBH))
logging.info('Stochastic Omega: %4.1g Universes' % score.Omega)
plot_noise(ifo, traces, **plot_style)
return score, noises, ifo
gwinc/__main__.py
View file @ b3adea0a
......@@ -9,9 +9,8 @@ import logging
logging.basicConfig(format='%(message)s',
level=os.getenv('LOG_LEVEL', logging.INFO))
from .ifo import available_ifos, load_ifo
from . import available_ifos, load_budget, plot_noise
from .precomp import precompIFO
from . import plot_noise
from . import io
##################################################
......@@ -97,10 +96,12 @@ def main():
plot_style = attrs
else:
Budget, ifo, freq, plot_style = load_ifo(args.IFO)
Budget = load_budget(args.IFO)
ifo = Budget.ifo
# FIXME: this should be done only if specified, to allow for
# using any FREQ specified in budget module by default
# using any FREQ specified in the Budget
freq = np.logspace(np.log10(args.flo), np.log10(args.fhi), args.npoints)
plot_style = getattr(Budget, 'plot_style', {})
traces = None
if args.yaml:
......@@ -117,7 +118,7 @@ def main():
if not ifo:
parser.exit(2, "no IFO structure available.")
fmt = '{:30} {:>20} {:>20}'
Budget, ifoo, freq, plot_style = load_ifo(args.diff)
Budget = load_ifo(args.diff)
diffs = ifo.diff(ifoo)
if diffs:
print(fmt.format('', args.IFO, args.diff))
......@@ -178,10 +179,7 @@ def main():
logging.info("precomputing ifo...")
ifo = precompIFO(freq, ifo)
logging.info("calculating budget...")
budget = Budget(freq=freq, ifo=ifo)
budget.load()
budget.update()
traces = budget.calc_trace()
traces = Budget(freq=freq, ifo=ifo).run()
# logging.info('recycling factor: {: >0.3f}'.format(ifo.gwinc.prfactor))
# logging.info('BS power: {: >0.3f} W'.format(ifo.gwinc.pbs))
......
gwinc/gwinc.py deleted 100644 → 0
View file @ 160b4a1d
from __future__ import division
import numpy as np
from numpy import pi, sqrt
from collections import OrderedDict
import logging
from . import load_ifo
from .precomp import precompIFO
from . import noise
from .plot import plot_noise
def gwinc(freq, ifoin, source=None, plot=False, PRfixed=True):
"""Calculate strain noise budget for a specified interferometer model.
Argument `freq` is the frequency array for which the noises will
be calculated, and `ifoin` is the IFO model (see the `load_ifo()`
function).
If `source` structure provided, so evaluates the sensitivity of
the detector to several potential gravitational wave
sources.
If `plot` is True a plot of the budget will be created.
Returns tuple of (score, noises, ifo)
"""
# assume generic aLIGO configuration
# FIXME: how do we allow adding arbitrary addtional noise sources
# from just ifo description, without having to specify full budget
Budget, ifo_, freq_, plot_style = load_ifo('aLIGO')
# add some precomputed info to the ifo struct
# this implicitly deepcopies and the return value is the copy
ifo = precompIFO(freq, ifoin, PRfixed)
traces = Budget(freq, ifo=ifo).calc_trace()
# construct matgwinc-compatible noises structure
noises = {}
for name, (data, style) in traces.items():
noises[style.get('label', name)] = data
noises['Freq'] = freq
pbs = ifo.gwinc.pbs
parm = ifo.gwinc.parm
finesse = ifo.gwinc.finesse
prfactor = ifo.gwinc.prfactor
if ifo.Laser.Power * prfactor != pbs:
pass
#warning(['Thermal lensing limits input power to ' num2str(pbs/prfactor, 3) ' W']);
#TODO decide if all this below this should remain, since it is already inside of __main__
# report astrophysical scores if so desired
score = None
if source:
score = int73(freq, noises['Total'], ifo, source)
score.Omega = intStoch(freq, noises['Total'], 0, ifo, source)
# --------------------------------------------------------
# output graphics
if plot:
# Report input parameters
if ifo.Optics.Type == 'DualCarrier_new': #include the case for Dual carrier
finesseA = 2*pi/ifo.Optics.ITM.TransmittanceD1
finesseB = 2*pi/ifo.Optics.ITM.TransmittanceD2
pbsA = ifo.Laser.PBSD1
pbsB = ifo.Laser.PBSD2
logging.info('Finesse for carrier A: %7.2f' % finesseA)
logging.info('Finesse for carrier B: %7.2f' % finesseB)
logging.info('Power Recycling Factor: %7.2f' % ifo.PRCgain)
logging.info('Arm power for carrier A:%7.2f kW' % (finesseA*2/pi*pbsA/2/1000))
logging.info('Arm power for carrier B:%7.2f kW' % (finesseB*2/pi*pbsB/2/1000))
logging.info('Power on beam splitter for carrier A: %7.2f W' % pbsA)
logging.info('Power on beam splitter for carrier B: %7.2f W' % pbsB)
logging.info('Laser Power for Carrier A: %7.2f Watt' % ifo.LP1)
logging.info('Laser Power for Carrier B: %7.2f Watt' % ifo.LP2)
logging.info('SRM Detuning for Carrier A: %7.2f degree' % (ifo.Optics.SRM.TunephaseD1*180/pi))
logging.info('SRM Detuning for Carrier B: %7.2f degree' % (ifo.Optics.SRM.TunephaseD2*180/pi))
logging.info('SRM transmission for Carrier A:%9.4f' % ifo.Optics.SRM.TransmittanceD1)
logging.info('SRM transmission for Carrier B:%9.4f' % ifo.Optics.SRM.TransmittanceD2)
logging.info('ITM transmission for Carrier A:%9.4f' % ifo.Optics.ITM.TransmittanceD1)
logging.info('ITM transmission for Carrier B:%9.4f' % ifo.Optics.ITM.TransmittanceD2)
logging.info('PRM transmission for both: %9.4f' % ifo.Optics.PRM.Transmittance)
else:
logging.info('Laser Power: %7.2f Watt' % ifo.Laser.Power)
logging.info('SRM Detuning: %7.2f degree' % (ifo.Optics.SRM.Tunephase*180/pi))
logging.info('SRM transmission: %9.4f' % ifo.Optics.SRM.Transmittance)
logging.info('ITM transmission: %9.4f' % ifo.Optics.ITM.Transmittance)
logging.info('PRM transmission: %9.4f' % ifo.Optics.PRM.Transmittance)
logging.info('Finesse: %7.2f' % finesse)
logging.info('Power Recycling Gain: %7.2f' % prfactor)
logging.info('Arm Power: %7.2f kW' % (parm/1000))
logging.info('Power on BS: %7.2f W' % pbs)
# coating and substrate thermal load on the ITM
PowAbsITM = (pbs/2) * \
np.hstack([(finesse*2/pi) * ifo.Optics.ITM.CoatingAbsorption,
(2 * ifo.Materials.MassThickness) * ifo.Optics.ITM.SubstrateAbsorption])
# Stefan's Mysterious TCS Section ~~~~ `~~ ` ` 324@#%@#$ !
M = np.array([[ifo.TCS.s_cc, ifo.TCS.s_cs], [ifo.TCS.s_cs, ifo.TCS.s_ss]])
S_uncorr = PowAbsITM.T*M*PowAbsITM
TCSeff = 1-sqrt(ifo.TCS.SRCloss/S_uncorr)
logging.info('Thermal load on ITM: %8.3f W' % sum(PowAbsITM))
logging.info('Thermal load on BS: %8.3f W' %
(ifo.Materials.MassThickness*ifo.Optics.SubstrateAbsorption*pbs))
if (ifo.Laser.Power*prfactor != pbs):
logging.info('Lensing limited input power: %7.2f W' % (pbs/prfactor))
if source:
logging.info('BNS Inspiral Range: ' + str(score.effr0ns) + ' Mpc/ z = ' + str(score.zHorizonNS))
logging.info('BBH Inspiral Range: ' + str(score.effr0bh) + ' Mpc/ z = ' + str(score.zHorizonBH))
logging.info('Stochastic Omega: %4.1g Universes' % score.Omega)
plot_noise(ifo, traces, **plot_style)
return score, noises, ifo
gwinc/gwinc_matlab.py
View file @ b3adea0a
......@@ -120,7 +120,7 @@ def ifo_matlab_transform(ifo):
"""
# add constants
CONSTS = {k:v for k, v in const.__dict__ if not k.startswith('__')}
ifo.Constants = Struct.from_dict(CONSTS)
ifo.Constants = Struct(CONSTS)
# copy tempurature into Constants
ifo.Constants.Temp = ifo.Infrastructure.Temp
......
gwinc/ifo/__init__.py
View file @ b3adea0a
import os
import logging
from ..struct import load_struct, STRUCT_EXT
from ..util import load_module
PLOT_STYLE = dict(
......@@ -10,6 +6,20 @@ PLOT_STYLE = dict(
)
def lpath(file0, file1):
"""Return path of file1 when expressed relative to file0.
For instance, if file0 is "/path/to/file0" and file1 is
"../for/file1" then what is returned is "/path/for/file1".
This is useful for resolving paths within packages with e.g.:
rpath = lpath(__file__, '../resource.yaml')
"""
return os.path.abspath(os.path.join(os.path.dirname(file0), file1))
def available_ifos():
"""List available included pre-defined IFOs
......@@ -20,57 +30,3 @@ def available_ifos():
if os.path.isdir(os.path.join(root, f)) and f[0] != '_':
ifos.append(f)
return sorted(ifos)
def load_ifo(name_or_path):
"""Load GWINC IFO Budget by name or from file.
Named IFOs should correspond to one of the IFOs available in the
gwinc package (see gwinc.available_ifos()). If a path is provided
it should either be a budget package (directory) or module (ending
in .py), or an IFO struct (see gwinc.load_struct()). In the
latter case the base aLIGO budget definition will be used.
Returns primary Budget class, ifo structure, frequency array, and
plot style dictionary, with the last three being None if they are
not defined in the budget.
"""
ifo = None
if os.path.exists(name_or_path):
path = name_or_path.rstrip('/')
bname, ext = os.path.splitext(os.path.basename(path))
if ext in STRUCT_EXT:
logging.info("loading struct {}...".format(path))
ifo = load_struct(path)
bname = 'aLIGO'
modname = 'gwinc.ifo.aLIGO'
logging.info("loading budget {}...".format(modname))
else:
modname = path
logging.info("loading module path {}...".format(modname))
else:
if name_or_path not in available_ifos():
raise RuntimeError("Unknonw IFO '{}' (available IFOs: {}).".format(
name_or_path,
available_ifos(),
))
bname = name_or_path
modname = 'gwinc.ifo.'+name_or_path
logging.info("loading module {}...".format(modname))
mod, modpath = load_module(modname)
Budget = getattr(mod, bname)
ifo = getattr(mod, 'IFO', ifo)
ifopath = os.path.join(modpath, 'ifo.yaml')
if not ifo and ifopath:
ifo = load_struct(ifopath)
freq = getattr(mod, 'FREQ', None)
plot_style = getattr(mod, 'PLOT_STYLE', PLOT_STYLE)
return Budget, ifo, freq, plot_style
gwinc/nb.py
View file @ b3adea0a
......@@ -32,7 +32,7 @@ class BudgetItem:
return None
def update(self, **kwargs):
"""Overload method for updating data needed to calculate final PSD.
"""Overload method for updating data.
By default any keyword arguments provided are written directly
as attribute variables (as with __init__).
......@@ -42,7 +42,7 @@ class BudgetItem:
setattr(self, key, val)
def calc(self):
"""Overload method for calculation of final PSD.
"""Overload method for final PSD calculation.
Should return an array of power-referenced values evaluated at
all evaluation frequencies (self.freq).
......@@ -152,6 +152,17 @@ class Noise(BudgetItem):
data = None
return data, self.style
def run(self, **kwargs):
"""Convenience method to load, update, and return calc traces.
Equivalent of load(), update(), calc_traces() run in sequence.
Keyword arguments are passed to update().
"""
self.load()
self.update(**kwargs)
return self.calc_trace()
class Budget(Noise):
"""GWINC Budget class
......
gwinc/noise/substratethermal.py
View file @ b3adea0a
......@@ -21,13 +21,12 @@ def carrierdensity_adiabatic(f, ifo):
diffHole = ifo.Materials.Substrate.HoleDiffusion
cdDens = ifo.Materials.Substrate.CarrierDensity
r0 = ifo.Optics.ITM.BeamRadius/np.sqrt(2)
L = ifo.Infrastructure.Length
gPhase = ifo.gwinc.finesse*2/pi
psdElec = 4*H*gammaElec**2*cdDens*diffElec/(pi*r0**4*Omega**2) # units are meters
psdHole = 4*H*gammaHole**2*cdDens*diffHole/(pi*r0**4*Omega**2) # units are meters
psdMeters = 2 * (psdElec + psdHole) # electrons and holes for two ITMs
n = psdMeters / (gPhase*L)**2
n = psdMeters / (gPhase)**2 * ifo.gwinc.dhdl_sqr
return n
......@@ -38,7 +37,6 @@ def carrierdensity_exact(f, ifo):
"""
w = ifo.Optics.ITM.BeamRadius
L = ifo.Infrastructure.Length
H = ifo.Materials.MassThickness
kBT = const.kB * ifo.Materials.Substrate.Temp
hbar = const.hbar
......@@ -71,7 +69,7 @@ def carrierdensity_exact(f, ifo):
psdMeters = 2 * (psdElec + psdHole)
n = psdMeters / (gPhase*L)**2
n = psdMeters / (gPhase)**2 * ifo.gwinc.dhdl_sqr
return n
......@@ -91,12 +89,11 @@ def thermorefractiveITM_adiabatic(f, ifo):
Temp = ifo.Materials.Substrate.Temp
kBT = const.kB * Temp
r0 = ifo.Optics.ITM.BeamRadius/np.sqrt(2)
L = ifo.Infrastructure.Length
gPhase = ifo.gwinc.finesse*2/pi
psd = 4*H*beta**2*kappa*kBT*Temp/(pi*r0**4*Omega**2*(rho*C)**2) # units are meters
psdMeters = 2*psd # two ITMs
n = psdMeters / (gPhase*L)**2
n = psdMeters / (gPhase)**2 * ifo.gwinc.dhdl_sqr
return n
......@@ -108,7 +105,6 @@ def thermorefractiveITM_exact(f, ifo):
"""
w = ifo.Optics.ITM.BeamRadius
L = ifo.Infrastructure.Length
H = ifo.Materials.MassThickness
kBT = const.kB * ifo.Materials.Substrate.Temp
Temp = ifo.Materials.Substrate.Temp
......@@ -137,7 +133,7 @@ def thermorefractiveITM_exact(f, ifo):
psdMeters = 2*psd # two itms
n = psdMeters / (gPhase*L)**2
n = psdMeters / (gPhase)**2 * ifo.gwinc.dhdl_sqr
return n
......@@ -156,7 +152,6 @@ def subbrownian(f, ifo):
c2 = ifo.Materials.Substrate.c2
n = ifo.Materials.Substrate.MechanicalLossExponent
alphas = ifo.Materials.Substrate.Alphas
L = ifo.Infrastructure.Length
kBT = const.kB * ifo.Materials.Substrate.Temp
# Bulk substrate contribution
......
gwinc/struct.py
View file @ b3adea0a
......@@ -57,7 +57,10 @@ class Struct(object):
"""
# FIXME: This would be a way to allow setting nested struct
STRUCT_EXT = ['.yaml', '.yml', '.mat', '.m']
"""accepted extension types for struct files"""
# FIXME: There should be a way to allow setting nested struct
# attributes, e.g.:
#
# >>> s = Struct()
......@@ -75,11 +78,15 @@ class Struct(object):
##########
def __init__(self, **kwargs):
"""Arguments can pre-fill the structure
def __init__(self, *args, **kwargs):
"""Initialize Struct object
Initializes similar to dict(), taking a single dict or mapping
argument, or keyword arguments to initially populate the
Struct.
"""
self.__dict__.update(kwargs)
self.update(dict(*args, **kwargs))
def __getitem__(self, key):
"""Get a (possibly nested) value from the struct.
......@@ -96,7 +103,7 @@ class Struct(object):
else:
return self.__dict__[key]
def get(self, key, default):
def get(self, key, default=None):
"""Get a (possibly nested) value from the struct, or default.
"""
......@@ -115,6 +122,34 @@ class Struct(object):
def setdefault(self, key, default):
return self.__dict__.setdefault(key, default)
def update(self, other):
"""Update Struct from other Struct or dict.
"""
if isinstance(other, Struct):
d = other.__dict__
else:
d = dict(other)
for k, v in d.items():
if k in self:
if isinstance(self[k], Struct) \
and isinstance(v, (dict, Struct)):
self[k].update(v)
continue
try:
delattr(self, k)
except AttributeError:
delattr(self.__class__, k)
if isinstance(v, dict):
self[k] = Struct(v)
elif isinstance(v, (list, tuple)):
try:
self[k] = list(map(Struct, v))
except TypeError:
self[k] = v
else:
self[k] = v
def items(self):
return self.__dict__.items()
......@@ -215,6 +250,9 @@ class Struct(object):
diffs.append((k, v, ov))
return diffs
def __eq__(self, other):
"""True if structs have all equal values"""
return not bool(self.diff(other))
def to_txt(self, path=None, fmt='0.6e', delimiter=': ', end=''):
"""Return text represenation of Struct, one element per line.
......@@ -251,30 +289,13 @@ class Struct(object):
return txt.getvalue()
@classmethod
def from_dict(cls, d):
"""Create Struct from nested dict.
"""
c = cls()
for k,v in d.items():
if type(v) == dict:
c.__dict__[k] = Struct.from_dict(v)
else:
try:
c.__dict__[k] = list(map(Struct.from_dict, v))
except (AttributeError, TypeError):
c.__dict__[k] = v
return c
@classmethod
def from_yaml(cls, y):
"""Create Struct from YAML string.
"""
d = yaml.load(y)
return cls.from_dict(d)
d = yaml.load(y, Loader=yaml_loader)
return cls(d)
@classmethod
......@@ -310,45 +331,39 @@ class Struct(object):
def from_file(cls, path):
"""Load Struct from .yaml or MATLAB .mat file.
File type will be determined by extension.
Accepted file types are .yaml, .mat, or .m.
For .m files, the file is expected to include either an object
or function that corresponds to the basename of the file. The
MATLAB engine will be invoked to execute the .m code and
extract the resultant IFO data.
If `path` is a tuple, all elements will be joined ala
os.path.join, with the first element resolved to it's absolute
dirname. This is useful for loading package-relative files
with e.g.:
Struct.from_file((__file__, 'myifo.yaml'))
"""
(root, ext) = os.path.splitext(path)
if type(path) == tuple:
path = os.path.join(os.path.abspath(os.path.dirname(path[0])), *path[1:])
base, ext = os.path.splitext(path)
if ext == '.m':
from ..gwinc_matlab import Matlab
matlab = Matlab()
matlab.addpath(os.path.dirname(path))
func_name = os.path.basename(base)
matlab.eval("ifo = {};".format(func_name), nargout=0)
ifo = matlab.extract('ifo')
return Struct.from_matstruct(ifo)
with open(path, 'r') as f:
if ext in ['.yaml', '.yml']:
d = yaml.load(f, Loader=yaml_loader)
return cls.from_dict(d)
return cls.from_yaml(f)
elif ext == '.mat':
s = loadmat(f, squeeze_me=True, struct_as_record=False)
return cls.from_matstruct(s)
else:
raise IOError("Unknown file type: {}".format(ext))
def load_struct(path):
"""Load struct from YAML or MATLAB file.
Files may be either .yaml, .mat or .m. For .m files, the file is
expected to include either an object or function that corresponds
to the basename of the file. The MATLAB engine will be invoked to
execute the .m code and extract the resultant IFO data.
"""
root, ext = os.path.splitext(path)
if ext == '.m':
from ..gwinc_matlab import Matlab
matlab = Matlab()
matlab.addpath(os.path.dirname(path))
func_name = os.path.basename(root)
matlab.eval("ifo = {};".format(func_name), nargout=0)
ifo = matlab.extract('ifo')
return Struct.from_matstruct(ifo)
else:
return Struct.from_file(path)
# accepted extension types for struct files
STRUCT_EXT = ['.yaml', '.yml', '.mat', '.m']
gwinc/test/Aplus.pkl deleted 100644 → 0
View file @ 160b4a1d
File deleted
gwinc/test/__main__.py
View file @ b3adea0a
import os
import sys
import glob
import signal
import pickle
import subprocess
import hashlib
import logging
import tempfile
import argparse
import numpy as np
import matplotlib.pyplot as plt
import argparse
from collections import OrderedDict
from collections.abc import Mapping
from PyPDF2 import PdfFileReader, PdfFileWriter
import logging
logging.basicConfig(format='%(message)s',
level=os.getenv('LOG_LEVEL', logging.INFO))
from .. import load_ifo
from ..gwinc import gwinc
from ..gwinc_matlab import gwinc_matlab
from .. import available_ifos, load_budget
from .. import precompIFO
from .. import load_hdf5, save_hdf5
try:
import inspiral_range
except ImportError:
inspiral_range = None
FLO = 5
FHI = 6000
NPOINTS = 3000
TOLERANCE = 1e-6
CACHE_PATH = os.path.join(os.path.dirname(__file__), 'cache')
def path_hash(path):
"""Calculate SHA1 hash of path, either directory or file"""
if not path or not os.path.exists(path):
return
path = os.path.expanduser(path)
if os.path.isdir(path):
d = path
f = '.'
else:
d = os.path.dirname(path)
f = os.path.basename(path)
CWD = os.getcwd()
os.chdir(d)
cmd = 'find {} -type f ! -wholename "*/.*" -print0 | sort -z | xargs -0 sha1sum | sha1sum'.format(f)
sha1sum_out = subprocess.check_output(cmd, shell=True)
sha1sum = sha1sum_out.split()[0]
os.chdir(CWD)
return sha1sum.decode()
##################################################
def main():
parser = argparse.ArgumentParser()
parser.add_argument('--plot', '-p', action='store_true', help='plot differences')
parser.add_argument('--save', '-s', help='save plot to file')
parser.add_argument('--recalc', '-r', action='store_true', help='recalculate all traces')
parser.add_argument('--tolerance', '-t', help='fractional tolerance', type=float, default=1e-6)
parser.add_argument('--skip', '-k', action='append', help='traces to skip comparing')
parser.add_argument('IFO', help='IFO name or description file')
args = parser.parse_args()
def walk_traces(traces, root=()):
"""recursively walk through traces
logging.info("loading IFO '{}'...".format(args.IFO))
Budget, ifo, freq, plot_style = load_ifo(args.IFO)
yields (key_tuple, value), where key_tuple is a tuple of nested
dict keys of the traces dict locating the given value.
freq = np.logspace(np.log10(FLO), np.log10(FHI), NPOINTS)
##############################
# matgwinc processing
mdata_pkl = os.path.join(os.path.dirname(__file__), '{}.pkl'.format(args.IFO))
ifo_hash = hashlib.sha1(ifo.to_txt().encode()).hexdigest()
gwinc_hash = path_hash(os.getenv('GWINCPATH'))
if not gwinc_hash:
logging.warning("GWINCPATH not specified or does not exist; skipping check for changes to matgwinc code.")
mrecalc = args.recalc
"""
for key, val in traces.items():
r = root+(key,)
if isinstance(val, Mapping):
for k, v in walk_traces(val, root=r):
yield k, v
continue
else:
yield r, val
if os.path.exists(mdata_pkl) and not mrecalc:
mrecalc = False
logging.info("loading matgwinc data {}...".format(mdata_pkl))
with open(mdata_pkl, 'rb') as f:
if sys.version_info.major > 2:
mdata = pickle.load(f, encoding='latin1')
else:
mdata = pickle.load(f)
if mdata['ifo_hash'] != ifo_hash:
logging.info("ifo hash has changed: {}".format(ifo_hash))
mrecalc = True
if gwinc_hash and mdata['gwinc_hash'] != gwinc_hash:
logging.info("matgwinc hash has changed: {}".format(gwinc_hash))
mrecalc = True
def zip_noises(tracesA, tracesB, skip):
"""zip matching noises from traces"""
for keys, (noiseA, _) in walk_traces(tracesA):
# keys is a tuple of nested keys for noise
name = keys[-1]
# skip keys we'll add at the end
if name in ['Total', 'int73']:
continue
if skip and name in skip:
logging.warning("SKIPPING TEST: '{}'".format(name))
continue
if mrecalc:
logging.info("calculating matgwinc noises...")
try:
mscore, mnoises, mifo = gwinc_matlab(freq, ifo)
except (ImportError, OSError):
sys.exit("MATLAB engine not available.")
mdata = dict(score=mscore, noises=mnoises, ifo=mifo, ifo_hash=ifo_hash, gwinc_hash=gwinc_hash)
with open(mdata_pkl, 'wb') as f:
pickle.dump(mdata, f)
t = tracesB
for key in keys:
t = t[key]
noiseB, style = t
except (KeyError, TypeError):
logging.warning("MISSING B TRACE: '{}'".format(name))
continue
mnoises = mdata['noises']
yield name, noiseA, noiseB
##############################
# pygwinc processing
yield 'Total', tracesA['Total'][0], tracesB['Total'][0]
logging.info("calculating pygwinc noises...")
score, noises, ifo = gwinc(freq, ifo)
if 'int73' in tracesA:
yield 'int73', tracesA['int73'][0], tracesB['int73'][0]
##############################
# calc inspiral ranges
if inspiral_range:
logging.info("calculating inspiral ranges...")
def compare_traces(tracesA, tracesB, tolerance=TOLERANCE, skip=None):
"""Compare two gwinc trace dictionaries
range_func = inspiral_range.range
H = inspiral_range.waveform.CBCWaveform(freq)
Noises listed in `skip` will not be compared.
mfom = range_func(freq, mnoises['Total'], H=H)
_, mnoises['int73'] = inspiral_range.int73(freq, mnoises['Total'])
logging.info("matgwinc range: {:.2f} Mpc".format(mfom))
Returns a dictionary of noises that differ fractionally (on a
point-by-point basis) by more that `tolerance` between `tracesA`
and `tracesB`.
fom = range_func(freq, noises['Total'], H=H)
_, noises['int73'] = inspiral_range.int73(freq, noises['Total'])
logging.info("pygwinc range: {:.2f} Mpc".format(fom))
"""
#name_width = max([len(n[0][-1]) for n in walk_traces(tracesA)])
name_width = 15
diffs = OrderedDict()
for name, noiseA, noiseB in zip_noises(tracesA, tracesB, skip):
logging.debug("comparing {}...".format(name))
fom_title = """inspiral {func} {m1}/{m2} Msol:
matgwinc: {mfom:.2f} Mpc
pygwinc: {fom:.2f} Mpc""".format(
func=range_func.__name__,
m1=H.params['m1'],
m2=H.params['m2'],
mfom=mfom,
fom=fom,
)
ampA = np.sqrt(noiseA)
ampB = np.sqrt(noiseB)
diff = ampB - ampA
frac = abs(diff / ampA)
else:
fom_title = ''
if max(frac) < tolerance:
continue
##############################
# find differences
logging.warning("EXCESSIVE DIFFERENCE: {:{w}} {:6.1f} ppm".format(
name, max(frac)*1e6, w=name_width))
diffs[name] = (noiseA, noiseB, frac)
return diffs
def plot_diffs(freq, diffs, tolerance,
name, labelA, labelB, fom_title='',
save=None):
spec = (len(diffs)+1, 2)
sharex = None
for i, nname in enumerate(diffs):
noiseA, noiseB, frac = diffs[nname]
axl = plt.subplot2grid(spec, (i, 0), sharex=None)
axl.loglog(freq, np.sqrt(noiseA), label=labelA)
axl.loglog(freq, np.sqrt(noiseB), label=labelB)
axl.grid()
axl.legend(loc='upper right')
axl.set_ylabel(nname)
if i == 0:
axl.set_title("noise value")
if i == 0:
sharex = axl
axr = plt.subplot2grid(spec, (i, 1), sharex=sharex)
axr.loglog(freq, frac)
axr.grid()
axr.axhline(y=max(frac), color='r', linestyle='--')
axr.text(max(freq)+4000, max(frac), '{:.1f} ppm'.format(max(frac)*1e6),
horizontalalignment='left', verticalalignment='center',
color='red')
if i == 0:
axr.set_title("fractional difference")
plt.suptitle('''{} {}/{} noise comparison
(noises that differ by more than {} ppm)
{}'''.format(name, labelA, labelB, tolerance*1e6, fom_title))
axl.set_xlabel("frequency [Hz]")
axr.set_xlabel("frequency [Hz]")
plt.subplots_adjust(top=0.8, right=0.85, wspace=0.3)
if save:
pwidth = 10
pheight = (len(diffs) * 5) + 2
plt.gcf().set_size_inches(pwidth, pheight)
plt.savefig(save)
else:
plt.show()
skip = args.skip
fractional_tolerance = args.tolerance
##################################################
diffs = {}
for name, noise in noises.items():
if name in ['Freq']:
continue
if skip and name in skip:
logging.warning("SKIPPING TEST: '{}'".format(name))
continue
def main():
parser = argparse.ArgumentParser()
parser.add_argument('--tolerance', '-t', type=float, default=TOLERANCE,
help='fractional tolerance [{}]'.format(TOLERANCE))
parser.add_argument('--skip', '-k', metavar='NOISE', action='append',
help='traces to skip in comparison (multiple may be specified)')
parser.add_argument('--cache', '-c', metavar='PATH', default=CACHE_PATH,
help='specify alternate IFO traces cache path')
rgroup = parser.add_mutually_exclusive_group()
rgroup.add_argument('--plot', '-p', action='store_true',
help='plot differences')
rgroup.add_argument('--report', '-r', metavar='REPORT.pdf',
help='create PDF report of test results (only created if differences found)')
rgroup.add_argument('--gen-cache', action='store_true',
help='update/create IFO traces cache directory')
parser.add_argument('ifo', metavar='IFO', nargs='*',
help='specific ifos to test (default all)')
args = parser.parse_args()
if args.gen_cache:
try:
mnoise = mnoises[name]
except KeyError:
logging.warning("skipping missing MATGWINC trace: {}".format(name))
os.makedirs(args.cache)
except FileExistsError:
pass
freq = np.logspace(np.log10(5), np.log10(6000), 3000)
for name in available_ifos():
Budget = load_budget(name)
ifo = precompIFO(freq, Budget.ifo)
traces = Budget(freq, ifo=ifo).run()
path = os.path.join(args.cache, name+'.h5')
save_hdf5(path, freq, traces)
return
if args.report:
base, ext = os.path.splitext(args.report)
if ext != '.pdf':
parser.error("Test reports only support PDF format.")
outdir = tempfile.TemporaryDirectory()
# find all cached IFOs
logging.info("loading cache {}...".format(args.cache))
cached_ifos = {}
for f in sorted(os.listdir(args.cache)):
name, ext = os.path.splitext(f)
if ext != '.h5':
continue
# logging.info("compare: {}".format(name))
mn = mnoise
pn = noise
# _, mn = inspiral_range.int73(freq, mnoise)
# _, pn = inspiral_range.int73(freq, noise)
diff = np.sqrt(mn) - np.sqrt(pn)
frac = abs(diff / np.sqrt(pn))
cached_ifos[name] = os.path.join(args.cache, f)
if max(frac) < fractional_tolerance:
continue
# select
if args.ifo:
ifos = {name:cached_ifos[name] for name in args.ifo}
else:
ifos = cached_ifos
labelA = 'cache'
labelB = 'head'
fail = False
# compare
for name, path in ifos.items():
logging.info("{} tests...".format(name))
freq, tracesA, attrs = load_hdf5(path)
Budget = load_budget(name)
ifo = precompIFO(freq, Budget.ifo)
tracesB = Budget(freq, ifo=ifo).run()
if inspiral_range:
totalA = tracesA['Total'][0]
totalB = tracesB['Total'][0]
range_func = inspiral_range.range
H = inspiral_range.waveform.CBCWaveform(freq)
fomA = range_func(freq, totalA, H=H)
tracesA['int73'] = inspiral_range.int73(freq, totalA)[1], None
fomB = range_func(freq, totalB, H=H)
tracesB['int73'] = inspiral_range.int73(freq, totalB)[1], None
fom_summary = """
inspiral {func} {m1}/{m2} Msol:
{labelA}: {fomA:.2f} Mpc
{labelB}: {fomB:.2f} Mpc
""".format(
func=range_func.__name__,
m1=H.params['m1'],
m2=H.params['m2'],
labelA=labelA,
fomA=fomA,
labelB=labelB,
fomB=fomB,
)
else:
fom_summary = ''
logging.warning("EXCESSIVE DIFFERENCE: {:{w}} {:6.1f} ppm".format(
name, max(frac)*1e6, w=max([len(n) for n in noises])))
# logging.warning(" max: {:e}, min: {:e}".format(max(frac), min(frac)))
diffs[name] = (mn, pn, frac)
##############################
# plot
if args.plot:
spec = (len(diffs)+1, 2)
sharex = None
for i, name in enumerate(diffs):
mn, pn, frac = diffs[name]
axl = plt.subplot2grid(spec, (i, 0), sharex=sharex)
axl.loglog(freq, np.sqrt(pn), label='pygwinc')
axl.loglog(freq, np.sqrt(mn), label='matgwinc')
axl.grid()
axl.legend(loc='upper right')
axl.set_ylabel(name)
if i == 0:
sharex = axl
axr = plt.subplot2grid(spec, (i, 1), sharex=sharex)
axr.loglog(freq, frac)
axr.grid()
axr.axhline(y=max(frac), color='r', linestyle='--')
axr.text(max(freq)+4000, max(frac), '{:.1f} ppm'.format(max(frac)*1e6),
horizontalalignment='left', verticalalignment='center',
color='red')
diffs = compare_traces(tracesA, tracesB, args.tolerance, args.skip)
if diffs:
axl.set_xlabel("frequency [Hz]")
axr.set_xlabel("frequency [Hz]")
plt.suptitle("""{} mat/py gwinc noise comparison
noises that differ by more than {} ppm [(mat-py)/py]
{}""".format(args.IFO, fractional_tolerance*1e6, fom_title))
if args.save:
plt.gcf().set_size_inches(11, (len(diffs)+1)*4)
plt.savefig(args.save)
else:
plt.show()
logging.warning("{} tests FAIL".format(name))
fail |= True
if args.plot or args.report:
if args.report:
save = os.path.join(outdir.name, name+'.pdf')
else:
save = None
plot_diffs(
freq, diffs, args.tolerance,
name, labelA, labelB, fom_summary,
save=save,
)
else:
logging.warning("All tests passed, so no plot was generated")
##############################
if len(diffs) > 0:
return 1
return 0
logging.info("{} tests pass.".format(name))
if not fail:
logging.info("all tests pass.")
return 0
if args.report:
logging.info("generating report {}...".format(args.report))
pdf_writer = PdfFileWriter()
for path in glob.glob(os.path.join(outdir.name, '*.pdf')):
pdf_reader = PdfFileReader(path)
for page in range(pdf_reader.getNumPages()):
pdf_writer.addPage(pdf_reader.getPage(page))
with open(args.report, 'wb') as f:
pdf_writer.write(f)
outdir.cleanup()
logging.info("TESTS FAILED.")
return 1
##################################################
......
gwinc/test/aLIGO.pkl → gwinc/test/cache/Aplus.h5 LFS
View file @ b3adea0a
No preview for this file type
gwinc/test/cache/CE1.h5 0 → 100644 LFS
View file @ b3adea0a
File added
gwinc/test/cache/CE2.h5 0 → 100644 LFS
View file @ b3adea0a
File added
gwinc/test/cache/Voyager.h5 0 → 100644 LFS
View file @ b3adea0a
File added
gwinc/test/cache/aLIGO.h5 0 → 100644 LFS
View file @ b3adea0a
File added
gwinc/util.py deleted 100644 → 0
View file @ 160b4a1d
import os
import importlib
def lpath(file0, file1):
"""Return path of file1 when expressed relative to file0.
For instance, if file0 is "/path/to/file0" and file1 is
"../for/file1" then what is returned is "/path/for/file1".
This is useful for resolving paths within packages with e.g.:
rpath = lpath(__file__, '../resource.yaml')
"""
return os.path.abspath(os.path.join(os.path.dirname(file0), file1))
def load_module(name_or_path):
"""Load module from name or path.
Return loaded module and module path.
"""
if os.path.exists(name_or_path):
path = name_or_path.rstrip('/')
modname = os.path.splitext(os.path.basename(path))[0]
if os.path.isdir(path):
path = os.path.join(path, '__init__.py')
spec = importlib.util.spec_from_file_location(modname, path)
mod = importlib.util.module_from_spec(spec)
spec.loader.exec_module(mod)
else:
mod = importlib.import_module(name_or_path)
try:
path = mod.__path__[0]
except AttributeError:
path = mod.__file__
return mod, path
setup.py 100644 → 100755
View file @ b3adea0a
......@@ -30,8 +30,10 @@ setup_args = dict(
install_requires = [
'h5py',
'ipython',
'matplotlib',
'numpy',
'PyPDF2',
'pyyaml',
'scipy',
],
......