Commit 39d430dd authored by Jameson Graef Rollins's avatar Jameson Graef Rollins

new nb noise budget module

This patch provides a new nb sub-module that defines classes for managing
and calculating noise budgets.  It provides the following overridable
classes:

nb.Calibration  A noise calibration
nb.Noise        A noise source
nb.Budget       A budget of noises

The Budget class includes a calc_trace() method that will return a traces
dictionary that includes data and trace plot styling for every noise term
in the budget recursively.

The existing included interferometers are updated to define their budgets
using this new interface, and the plot_noises function is updated to
accept the new traces dictionary.

An HDF5_SCHEMA describes how trace dictionaries are encoded into HDF5 files.
The new io module includes functions for writing traces to HDF5 files, and
for reading traces stored in this format.

The command line interface is updated to handle this new structure.
parent 1ae0f5cb
Pipeline #74044 passed with stages
in 1 minute and 12 seconds
......@@ -15,10 +15,10 @@ test:
- export PYTHONPATH=inspiral_range
- export MPLBACKEND=agg
- python3 -m gwinc aLIGO -s aLIGO.png
- python3 -m gwinc A+ -s A+.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" A+ -p -s A+_test.png
- python3 -m gwinc.test -t 100e-6 -k Seismic -k "Substrate Thermo-Elastic" Aplus -p -s Aplus_test.png
after_script:
- rm gitID.txt
cache:
......@@ -29,10 +29,10 @@ test:
expire_in: 4w
paths:
- aLIGO.png
- A+.png
- Aplus.png
- Voyager.png
- aLIGO_test.png
- A+_test.png
- Aplus_test.png
pages:
stage: deploy
......@@ -41,10 +41,10 @@ pages:
script:
- mkdir public
- mv aLIGO.png public/
- mv A+.png public/
- mv Aplus.png public/
- mv Voyager.png public/
- mv aLIGO_test.png public/ || true
- mv A+_test.png public/ || true
- mv Aplus_test.png public/ || true
- apt-get install -y -qq python3-pip python3-dev make
- pip3 install sphinx sphinx-rtd-theme
- cd docs
......
# HDF5 Schema for GWINC noise trace storage
This file describes a schemata for HDF5 storage of noise trace data
and plot styling GWINC noise budgets.
HDF5 is a heirarchical, structured data storage format [0]. Content
is organized into a heirarchical folder-like structure, with two
types of named objects:
* groups: holder of other objects (like folders)
* datasets: holder of data arrays (like files)
Objects can also have attributes as (almost) arbitrary key:value
pairs.
Bindings are available for most platforms including Python [1] and
Matlab.
[0] https://en.wikipedia.org/wiki/Hierarchical_Data_Format
[1] http://www.h5py.org/
## version history
v1
- first versioned schema release
## schema
The following describes the noise budget schema. Specific strings are
enclosed in single quotes (''), and variables are described in
brackets (<>). Group objects are indicated by a closing '/'
separator, data set are indicated by a closing ':' followed by a
specification of their length and type (e.g. "(N),float"), and
attributes are specified in the .attrs[] dictionary format. Optional
elements are enclosed in parentheses.
A single trace is a length N array (with optional plot style specified
in attributes:
```
/<trace>: (N),float
(/<trace>.attrs['label'] = <label>)
(/<trace>.attrs['color] = <color>)
...
```
A budget item, i.e. a collection of noises is structured as follows:
```
/<budget>/
/'Total': (N),float
/<trace_0>: (N),float
(/<trace_1>: (N),float)
```
<!-- ``` -->
<!-- /'noises'/ -->
<!-- /*<noise_0>: (N),float -->
<!-- ... -->
<!-- /\*'references'/ -->
<!-- /*<ref_0>: (N),float -->
<!-- ... -->
<!-- ``` -->
## Top-level Budget
The following two root attributes expected: a string describing the schema,
and an int schema version number:
```
/.attrs['SCHEMA'] = 'GWINC Noise Budget'
/.attrs['SCHEMA_VERSION'] = 1
```
Top-level attributes are generally used for specifying overall plot styling, but the
following root attributes are typically defined:
```
/.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
/.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>
```
The budget frequency array is defined in a 'Freq' dataset:
```
/'Freq': (N), float
```
The budget traces are defined a traces group. The overall structure
looks something like this:
```
/.attrs['SCHEMA'] = 'GWINC Noise Budget'
/.attrs['SCHEMA_VERSION'] = 1
/.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
/.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>
/'Freq': (N), float
/traces/
/'Total': (N),float
/<noise_0>: (N),float
/<noise_1>: (N),float
/<noise_2>/
/'Total': (N),float
/<noise_3>: (N),float
/<noise_4>: (N),float
...
```
......@@ -16,14 +16,11 @@ description is loaded, the noise budget can be calculated and plotted:
>>> import gwinc
>>> import numpy as np
>>> freq = np.logspace(1, 3, 1000)
>>> ifo = gwinc.load_ifo('aLIGO')
>>> Budget, ifo, freq_, plot_style = gwinc.load_ifo('aLIGO')
>>> ifo = gwinc.precompIFO(freq, ifo)
>>> noises = gwinc.noise_calc(freq, ifo)
>>> gwinc.plot_noise(noises)
```
Or the `gwinc` convenience function can be used to handle it all:
```
>>> score, data, ifo = gwinc.gwinc(freq, ifo, plot=True)
>>> traces = Budget(freq, ifo=ifo).calc_trace()
>>> fig = gwinc.plot_noise(freq, traces, **plot_style)
>>> fig.show()
```
......@@ -32,39 +29,7 @@ Or the `gwinc` convenience function can be used to handle it all:
You can make gwinc plots directly from the command line by executing
the package directly:
```shell
$ python3 -m gwinc -h
usage: gwinc [-h] [--flo FLO] [--fhi FHI] [--npoints NPOINTS] [--title TITLE]
[--matlab] [--fom FOM] [--dump | --save SAVE | --interactive]
[IFO]
Plot GWINC noise budget for specified IFO.
If the inspiral_range package is installed, various figures of merit
can be calculated for the resultant spectrum with the --fom argument,
e.g.:
gwinc --fom horizon ...
gwinc --fom range:m1=20,m2=20 ...
See documentation for inspiral_range package for details.
positional arguments:
IFO IFO name or description file path (.yaml or .mat)
optional arguments:
-h, --help show this help message and exit
--flo FLO, -fl FLO lower frequency bound in Hz [5]
--fhi FHI, --fh FHI upper frequency bound in Hz [6000]
--npoints NPOINTS, -n NPOINTS
number of frequency points [3000]
--title TITLE, -t TITLE
plot title
--matlab, -m use MATLAB gwinc engine to calculate noises
--fom FOM calculate inspiral range for resultant spectrum
('func[:param=val,param=val]')
--dump, -d print IFO parameters to stdout and exit
--save SAVE, -s SAVE save figure to file
--interactive, -i open interactive shell when plotting
$ python3 -m gwinc aLIGO
```
......@@ -80,6 +45,51 @@ various detectors:
* [Voyager.yaml](https://git.ligo.org/gwinc/pygwinc/blob/master/gwinc/ifo/Voyager.yaml)
## noise budgets
GWINC provides an `nb` package for defining arbitrary noise budgets:
```python
import numpy as np
from gwinc import nb
from gwinc import noise
class ExcessGas(nb.Noise):
"""Excess gas"""
style = dict(
label='Excess Gas',
color='#ad900d',
linestyle='--',
)
def calc(self):
return noise.residualgas.gas(self.freq, self.ifo)
class MeasuredNoise(nb.Noise):
"""My measured noise"""
style = dict(
label='Measured Noise,
color='#838209,
linestyle='-',
)
def load(self):
psd, freq = np.loadtxt('/path/to/measured/psd.txt')
self.data = self.interpolate(f, psd)
def calc(self):
return self.data
class MyBudget(nb.Budget):
noises = [
ExcessGas,
MeasuredNoise,
]
```
## comparison with MATLAB gwinc
`pygwinc` includes the ability use MATLAB gwinc directly via the
......
from .struct import load_struct
from .ifo import available_ifos, load_ifo
from .struct import load_struct
from .precomp import precompIFO
from .gwinc import noise_calc
from .gwinc import gwinc
from .plot import plot_noise
from .io import load_hdf5, save_hdf5
......@@ -11,9 +11,8 @@ logging.basicConfig(format='%(message)s',
from .ifo import available_ifos, load_ifo
from .precomp import precompIFO
from . import gwinc
from . import plot_noise
from . import util
from . import io
##################################################
......@@ -30,10 +29,11 @@ By default a GWINC noise budget of the specified IFO will calculated,
and plotted with an interactive plotter. If the --save option is
specified the plot will be saved directly to a file (without display)
(various formats are supported, indicated by file extension). If the
requested extension is "hdf5" then the noise traces and IFO parameters
will be saved to an HDF5 file (without plotting). The input file
(IFO) can be an HDF5 file saved from a previous call, in which all
noise traces and IFO parameters will be loaded from that file.
requested extension is 'hdf5' or 'h5' then the noise traces and IFO
parameters will be saved to an HDF5 file (without plotting). The
input file (IFO) can be an HDF5 file saved from a previous call, in
which case all noise traces and IFO parameters will be loaded from
that file.
If the inspiral_range package is installed, various figures of merit
can be calculated for the resultant spectrum with the --fom argument,
......@@ -51,9 +51,10 @@ FLO = 5
FHI = 6000
NPOINTS = 3000
parser = argparse.ArgumentParser(prog='gwinc',
description=description,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser = argparse.ArgumentParser(
prog='gwinc',
description=description,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('--flo', '-fl', default=FLO, type=float,
help="lower frequency bound in Hz [{}]".format(FLO))
parser.add_argument('--fhi', '--fh', default=FHI, type=float,
......@@ -68,7 +69,7 @@ group = parser.add_mutually_exclusive_group()
group.add_argument('--interactive', '-i', action='store_true',
help="interactive plot with interactive shell")
group.add_argument('--save', '-s',
help="save noise (hdf5) or figure (pdf/png/svg) to file")
help="save budget traces (.hdf5/.h5) or plot (.pdf/.png/.svg) to file")
group.add_argument('--yaml', '-y', action='store_true',
help="print IFO as yaml to stdout and exit")
group.add_argument('--text', '-x', action='store_true',
......@@ -77,8 +78,8 @@ group.add_argument('--diff', '-d', metavar='IFO',
help="show differences table between another IFO description")
group.add_argument('--no-plot', '-np', action='store_false', dest='plot',
help="supress plotting")
parser.add_argument('IFO', default=IFO,
help="IFO name or description file path (.yaml, .mat, .m), or full HDF5 data file")
parser.add_argument('IFO',
help="IFO name, description file path (.yaml, .mat, .m), budget module (.py), or HDF5 data file (.hdf5, .h5)")
def main():
......@@ -89,26 +90,34 @@ def main():
##########
# initial arg processing
if os.path.splitext(args.IFO)[1] in ['.hdf5', '.h5']:
title, ifo, noises = util.load_hdf5(args.IFO)
if os.path.splitext(os.path.basename(args.IFO))[1] in ['.hdf5', '.h5']:
Budget = None
freq, traces, attrs = io.load_hdf5(args.IFO)
ifo = getattr(attrs, 'IFO', None)
plot_style = attrs
else:
ifo = load_ifo(args.IFO)
noises = None
if args.title:
title = args.title
else:
title = '{} GWINC Noise Budget'.format(args.IFO)
Budget, ifo, freq, plot_style = load_ifo(args.IFO)
# FIXME: this should be done only if specified, to allow for
# using any FREQ specified in budget module by default
freq = np.logspace(np.log10(args.flo), np.log10(args.fhi), args.npoints)
traces = None
if args.yaml:
if not ifo:
parser.exit(2, "no IFO structure available.")
print(ifo.to_yaml(), end='')
return
if args.text:
if not ifo:
parser.exit(2, "no IFO structure available.")
print(ifo.to_txt(), end='')
return
if args.diff:
if not ifo:
parser.exit(2, "no IFO structure available.")
fmt = '{:30} {:>20} {:>20}'
ifoo = load_ifo(args.diff)
Budget, ifoo, freq, plot_style = load_ifo(args.diff)
diffs = ifo.diff(ifoo)
if diffs:
print(fmt.format('', args.IFO, args.diff))
......@@ -118,6 +127,14 @@ def main():
ov = repr(p[2])
print(fmt.format(k, v, ov))
return
if args.title:
plot_style['title'] = args.title
elif Budget:
plot_style['title'] = "GWINC Noise Budget: {}".format(Budget.name)
else:
plot_style['title'] = "GWINC Noise Budget: {}".format(args.IFO)
if args.plot:
if args.save:
# FIXME: this silliness seems to be the only way to have
......@@ -156,18 +173,20 @@ def main():
##########
# main calculations
if noises:
freq = noises['Freq']
else:
logging.info("calculating noises...")
freq = np.logspace(np.log10(args.flo), np.log10(args.fhi), args.npoints)
score, noises, ifo = gwinc(freq, ifo)
logging.info('recycling factor: {: >0.3f}'.format(ifo.gwinc.prfactor))
logging.info('BS power: {: >0.3f} W'.format(ifo.gwinc.pbs))
logging.info('arm finesse: {: >0.3f}'.format(ifo.gwinc.finesse))
logging.info('arm power: {: >0.3f} kW'.format(ifo.gwinc.parm/1000))
if not traces:
if ifo:
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()
# logging.info('recycling factor: {: >0.3f}'.format(ifo.gwinc.prfactor))
# logging.info('BS power: {: >0.3f} W'.format(ifo.gwinc.pbs))
# logging.info('arm finesse: {: >0.3f}'.format(ifo.gwinc.finesse))
# logging.info('arm power: {: >0.3f} kW'.format(ifo.gwinc.parm/1000))
if args.fom:
logging.info("calculating inspiral {}...".format(range_func))
......@@ -181,30 +200,35 @@ def main():
m2=H.params['m2'],
fom=fom,
)
title += '\n{}'.format(fom_title)
plot_style['title'] += '\n{}'.format(fom_title)
##########
# output
# save noise traces to HDF5 file
if args.save and os.path.splitext(args.save)[1] in ['.hdf5', '.h5']:
util.save_hdf5(
title,
ifo,
noises,
args.save,
logging.info("saving budget traces {}...".format(args.save))
if ifo:
plot_style['IFO'] = ifo.to_yaml()
io.save_hdf5(
path=args.save,
freq=freq,
traces=traces,
**plot_style
)
# interactive shell plotting
elif args.interactive:
ipshell = InteractiveShellEmbed(
user_ns={'freq': freq,
'noises': noises,
'ifo': ifo,
'plot_noise': plot_noise,
user_ns={
'freq': freq,
'traces': traces,
'ifo': ifo,
'plot_style': plot_style,
'plot_noise': plot_noise,
},
banner1='''
PYGWINC interactive plotter
GWINC interactive plotter
You may interact with plot using "plt." methods, e.g.:
......@@ -212,20 +236,21 @@ You may interact with plot using "plt." methods, e.g.:
>>> plt.savefig("foo.pdf")
''')
ipshell.enable_pylab()
ipshell.run_code("plot_noise(noises)")
ipshell.run_code("plt.title('{}')".format(title))
ipshell.run_code("plot_noise(freq, traces, **plot_style)")
ipshell.run_code("plt.title('{}')".format(plot_style['title']))
ipshell()
# standard plotting
elif args.plot:
logging.info("plotting...")
logging.info("plotting noises...")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plot_noise(
noises,
freq,
traces,
ax=ax,
**plot_style
)
ax.set_title(title)
fig.tight_layout()
if args.save:
fig.savefig(
......
from gwinc.ifo.noises import *
class Aplus(nb.Budget):
name = 'A+'
noises = [
QuantumVacuum,
Seismic,
Newtonian,
SuspensionThermal,
CoatingBrownian,
CoatingThermoOptic,
SubstrateBrownian,
SubstrateThermoElastic,
ExcessGas,
]
from gwinc.ifo.noises import *
class CE(nb.Budget):
name = 'Cosmic Explorer'
noises = [
QuantumVacuum,
Seismic,
Newtonian,
SuspensionThermal,
CoatingBrownian,
CoatingThermoOptic,
ITMThermoRefractive,
ITMCarrierDensity,
SubstrateBrownian,
SubstrateThermoElastic,
ExcessGas,
]
from gwinc.ifo.noises import *
class Voyager(nb.Budget):
name = 'Voyager'
noises = [
QuantumVacuum,
Seismic,
Newtonian,
SuspensionThermal,
CoatingBrownian,
CoatingThermoOptic,
ITMThermoRefractive,
ITMCarrierDensity,
SubstrateBrownian,
SubstrateThermoElastic,
ExcessGas,
]
import os
import fnmatch
import logging
from ..struct import Struct
from ..gwinc_matlab import Matlab
from ..struct import load_struct, STRUCT_EXT
from ..util import load_module
PLOT_STYLE = dict(
ylabel=u"Strain [1/\u221AHz]",
)
def available_ifos():
"""List available included IFO files"""
"""List available included pre-defined IFOs
"""
ifos = []
for f in os.listdir(os.path.dirname(__file__)):
if fnmatch.fnmatch(f, '*.yaml'):
ifos.append(f.split('.')[0])
return ifos
root = os.path.dirname(__file__)
for f in os.listdir(root):
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 IFO by name or from file.
"""Load GWINC IFO Budget by name or from file.
Named IFOs should correspond to the basename of .yaml IFO
definition files included with pygwinc (see available_ifos()
above).
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.
When specifying by path 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.
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
else:
path = os.path.join(os.path.dirname(__file__),
name_or_path+'.yaml')
path = name_or_path.rstrip('/')
bname, ext = os.path.splitext(os.path.basename(path))
(root, ext) = os.path.splitext(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))
if ext == '.m':
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:
modname = path
logging.info("loading module path {}...".format(modname))
else:
return Struct.from_file(path)
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
from gwinc.ifo.noises import *
class aLIGO(nb.Budget):
name = 'Advanced LIGO'
noises = [
QuantumVacuum,
Seismic,
Newtonian,
SuspensionThermal,
CoatingBrownian,
CoatingThermoOptic,
SubstrateBrownian,
SubstrateThermoElastic,
ExcessGas,
]
from .. import nb
from .. import noise
class QuantumVacuum(nb.Noise):
"""Quantum Vacuum
"""
style = dict(
label='Quantum',
color='#ad03de',
)
def calc(self):
return noise.quantum.shotrad(self.freq, self.ifo)
class Seismic(nb.Noise):
"""Seismic
"""
style = dict(
label='Seismic',
color='#653700',
)
def calc(self):
return noise.seismic.seismic(self.freq, self.ifo)
class Newtonian(nb.Noise):
"""Newtonian Gravity
"""
style = dict(
label='Newtonian Gravity',
color='#15b01a',
)
def calc(self):
return noise.newtonian.gravg(self.freq, self.ifo)
class SuspensionThermal(nb.Noise):
"""Suspension Thermal
"""
style = dict(
label='Suspension Thermal',
color='#0d75f8',
)
def calc(self):
return noise.suspensionthermal.susptherm(self.freq, self.ifo)
class CoatingBrownian(nb.Noise):
"""Coating Brownian
"""
style = dict(
label='Coating Brownian',
color='#fe0002',
)