...
 
Commits (5)
......@@ -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(ifo, 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 .ifo import load_ifo
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,10 +11,8 @@ logging.basicConfig(format='%(message)s',
from .ifo import available_ifos, load_ifo
from .precomp import precompIFO
from .gwinc import gwinc as pygwinc
from . import gwinc_matlab
from . import plot_noise
from . import util
from . import io
##################################################
......@@ -31,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,
......@@ -52,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,
......@@ -63,17 +63,13 @@ parser.add_argument('--npoints', '-n', default=NPOINTS,
help="number of frequency points [{}]".format(NPOINTS))
parser.add_argument('--title', '-t',
help="plot title")
parser.add_argument('--matlab', '-m', action='store_true',
help="use MATLAB gwinc engine to calculate noises")
parser.add_argument('--fom',
help="calculate inspiral range for resultant spectrum ('func[:param=val,param=val]')")
parser.add_argument('--no-displacement', '-nd', action='store_false', dest='displacement',
help="suppress displacement sensitivity axis")
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',
......@@ -82,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():
......@@ -94,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))
......@@ -123,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
......@@ -140,11 +152,6 @@ def main():
logging.warning("no display, plotting disabled.")
args.plot = False
if args.matlab:
gwinc = gwinc_matlab.gwinc_matlab
else:
gwinc = pygwinc
if args.fom:
import inspiral_range
try:
......@@ -166,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))
......@@ -191,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.:
......@@ -222,22 +236,21 @@ You may interact with plot using "plt." methods, e.g.:
>>> plt.savefig("foo.pdf")
''')
ipshell.enable_pylab()
ipshell.run_code("plot_noise(ifo, 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(
ifo,
noises,
freq,
traces,
ax=ax,
displacement=args.displacement,
**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',
)
def calc(self):
return noise.coatingthermal.coatbrownian(self.freq, self.ifo)
class CoatingThermoOptic(nb.Noise):
"""Coating Thermo-Optic
"""
style = dict(
label='Coating Thermo-Optic',
color='#02ccfe',
linestyle='--',
)
def calc(self):
return noise.coatingthermal.thermooptic(self.freq, self.ifo)
class ITMThermoRefractive(nb.Noise):
"""ITM Thermo-Refractive
"""
style = dict(
label='ITM Thermo-Refractive',
color='#448ee4',
linestyle='--',
)
def calc(self):
return noise.substratethermal.thermorefractiveITM(self.freq, self.ifo)
class ITMCarrierDensity(nb.Noise):
"""ITM Carrier Density
"""
style = dict(
label='ITM Carrier Density',
color='#929591',
linestyle='--',
)
def calc(self):
return noise.substratethermal.carrierdensity(self.freq, self.ifo)
class SubstrateBrownian(nb.Noise):
"""Substrate Brownian
"""
style = dict(
label='Substrate Brownian',
color='#fb7d07',
linestyle='--',
)
def calc(self):
return noise.substratethermal.subbrownian(self.freq, self.ifo)
class SubstrateThermoElastic(nb.Noise):
"""Substrate Thermo-Elastic
"""
style = dict(
label='Substrate Thermo-Elastic',
color='#f5bf03',
linestyle='--',
)
def calc(self):
return noise.substratethermal.subtherm(self.freq, self.ifo)
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)
import datetime
import h5py
SCHEMA = 'GWINC noise budget'
SCHEMA_VERSION = 1
def _write_trace_recursive(grp, traces):
for name, trace in traces.items():
if isinstance(trace, dict):
tgrp = grp.create_group(name)
_write_trace_recursive(tgrp, trace)
else:
data, style = trace
dset = grp.create_dataset(name, data=data)
for key, val in style.items():
dset.attrs[key] = val
def save_hdf5(path, freq, traces, **kwargs):
"""Save GWINC traces dict to HDF5 file.
See HDF5_SCHEMA.
"""
with h5py.File(path, 'w') as f:
f.attrs['SCHEMA'] = SCHEMA
f.attrs['SCHEMA_VERSION'] = SCHEMA_VERSION
# FIXME: add budget code hash or something
f.attrs['date'] = datetime.datetime.now().isoformat()
for key, val in kwargs.items():
f.attrs[key] = val
f.create_dataset('Freq', data=freq)
tgrp = f.create_group('traces')
_write_trace_recursive(tgrp, traces)
def _read_trace_recursive(element):
trace = {}
for name, item in element.items():
if isinstance(item, h5py.Group):
trace[name] = _read_trace_recursive(item)
else:
trace[name] = item.value, dict(item.attrs.items())
return trace
def load_hdf5(path):
"""Load traces from HDF5 file.
Returns a recursive traces dictionary. See HDF5_SCHEMA.
"""
with h5py.File(path, 'r') as f:
# FIXME: check SCHEMA name/version
freq = f['Freq'].value
traces = _read_trace_recursive(f['/traces'])
attrs = dict(f.attrs.items())
return freq, traces, attrs
This diff is collapsed.
from numpy import sqrt
PRIORITY_MAP = {
'Quantum Vacuum': 0,
'Seismic': 1,
'Newtonian Gravity': 2,
'Suspension Thermal': 3,
'Coating Brownian': 4,
'Coating Thermo-Optic': 5,
'ITM Thermo-Refractive': 6,
'ITM Carrier Density': 7,
'Substrate Brownian': 8,
'Substrate Thermo-Elastic': 9,
'Excess Gas': 10,
}
STYLE_MAP = {
'Quantum Vacuum': dict(
# color = 'xkcd:vibrant purple',
color = '#ad03de',
),
'Seismic': dict(
# color = 'xkcd:brown',
color = '#653700',
),
'Newtonian Gravity': dict(
# color = 'xkcd:green',
color = '#15b01a',
),
'Atmospheric Infrasound': dict(
# color = 'xkcd:neon green',
color = '#0cff0c',
ls = '--',
),
'Suspension Thermal': dict(
# color = 'xkcd:deep sky blue',
color = '#0d75f8',
),
'Coating Brownian': dict(
# color = 'xkcd:fire engine red',
color = '#fe0002',
),
'Coating Thermo-Optic': dict(
# color = 'xkcd:bright sky blue',
color = '#02ccfe',
ls = '--',
),
'ITM Thermo-Refractive': dict(
# color = 'xkcd:dark sky blue',
color = '#448ee4',
ls = '--',
),
'ITM Carrier Density': dict(
# color = 'xkcd:grey',
color = '#929591',
ls = '--',
),
'Substrate Brownian': dict(
# color = 'xkcd:pumpkin orange',
color = '#fb7d07',
ls = '--',
),
'Substrate Thermo-Elastic': dict(
# color = 'xkcd:golden',
color = '#f5bf03',
ls = '--',
),
'Excess Gas': dict(
# color = 'xkcd:baby shit brown',
color = '#ad900d',
ls = '--',
),
}
def plot_noise(
ifo,
noises,
freq,
traces,
ax=None,
displacement=True,
**kwargs
):
"""Plot a GWINC noise budget from calculated noises
If an axis handle is provided it will be used for the plot.
`displacement` may be set to False to supress the right had
displacement axis.
Returns the figure handle used.
Returns the figure handle.
"""
f = noises['Freq']
if ax is None:
import matplotlib.pyplot as plt
fig = plt.figure()
......@@ -97,31 +21,27 @@ def plot_noise(
else:
fig = ax.figure
def plot_dict(noises):
#use sorted to force a consistent ordering
#The key lambda in sorted gets the (name, noise) pair and so nn[0] returns the name
#the return tuple causes sorting by priority, with a fallback to lexical sort on the name
for name, noise in sorted(noises.items(), key = lambda nn : (PRIORITY_MAP.get(nn[0], 100), nn[0])):
if name in ['Freq', 'Total']:
continue
if isinstance(noise, dict):
plot_dict(noise)
else:
noise = noises[name]
stylekw = dict(
color = (0, 0, 0),
label = name,
lw = 3,
)
try:
style = STYLE_MAP[name]
stylekw.update(style)
except KeyError:
pass
ax.loglog(f, sqrt(noise), **stylekw)
plot_dict(noises)
ax.loglog(f, sqrt(noises['Total']), color='#000000', alpha=0.6, label='Total', lw=4)
for name, trace in traces.items():
if isinstance(trace, dict):
data, style = trace['Total']
else:
data, style = trace
# assuming all data is PSD
data = sqrt(data)
if name == 'Total':
style = dict(
color='#000000',
alpha=0.6,
lw=4,
)
ylim = [min(data)/10, max(data)]
if 'label' not in style:
style['label'] = name
if 'linewidth' in style:
style['lw'] = style['linewidth']
elif 'lw' not in style:
style['lw'] = 3
ax.loglog(freq, data, **style)
ax.grid(
True,
......@@ -136,16 +56,17 @@ def plot_noise(
fontsize='small',
)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel(u"Strain [1/\u221AHz]")
ax.set_xlim([min(f), max(f)])
ax.set_ylim([3e-25, 1e-21])
ax.autoscale(enable=True, axis='y', tight=True)
if 'ylim' in kwargs:
ax.set_ylim(kwargs['ylim'])
else:
ax.set_ylim(ylim)
ax.set_xlim(freq[0], freq[-1])
if displacement:
ax_d = ax.twinx()
ax_d.set_yscale('log')
y1, y2 = ax.get_ylim()
ax_d.set_ylim(y1*ifo.Infrastructure.Length, y2*ifo.Infrastructure.Length)
ax_d.set_ylabel(u"Displacement [m/\u221AHz]")
ax.set_xlabel('Frequency [Hz]')
if 'ylabel' in kwargs:
ax.set_ylabel(kwargs['ylabel'])
if 'title' in kwargs:
ax.set_title(kwargs['title'])
return fig
......@@ -324,3 +324,31 @@ class Struct(object):
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']
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
......@@ -12,7 +12,8 @@ import logging
logging.basicConfig(format='%(message)s',
level=os.getenv('LOG_LEVEL', logging.INFO))
from .. import load_ifo, gwinc
from .. import load_ifo
from ..gwinc import gwinc
from ..gwinc_matlab import gwinc_matlab
try:
import inspiral_range
......@@ -57,7 +58,7 @@ def main():
args = parser.parse_args()
logging.info("loading IFO '{}'...".format(args.IFO))
ifo = load_ifo(args.IFO)
Budget, ifo, freq, plot_style = load_ifo(args.IFO)
freq = np.logspace(np.log10(FLO), np.log10(FHI), NPOINTS)
......
import datetime
import h5py
import os
import importlib
from . import struct
def lpath(file0, file1):
"""Return path of file1 when expressed relative to file0.
def load_hdf5(path):
"""Load IFO and noises from HDF5 file.
For instance, if file0 is "/path/to/file0" and file1 is
"../for/file1" then what is returned is "/path/for/file1".
Returns (name, ifo, noises) tuple load from file.
This is useful for resolving paths within packages with e.g.:
rpath = lpath(__file__, '../resource.yaml')
"""
with h5py.File(path, 'r') as f:
name = f.attrs['name']
ifo = struct.Struct.from_yaml(f.attrs['IFO'])
noises = {}
for name, data in f['/traces'].items():
noises[name] = data.value
return name, ifo, noises
return os.path.abspath(os.path.join(os.path.dirname(file0), file1))
def load_module(name_or_path):
"""Load module from name or path.
def save_hdf5(name, ifo, noises, path):
"""Save IFO and noises to HDF5 file.
Return loaded module and module path.
"""
with h5py.File(path, 'w') as f:
f.attrs['schema'] = 'GWINC noise budget'
# FIXME: add GWINC version
f.attrs['date'] = datetime.datetime.now().isoformat()
f.attrs['name'] = name
f.attrs['IFO'] = ifo.to_yaml()
tgrp = f.create_group('traces')
for noise, data in noises.items():
tgrp.create_dataset(noise, data=data)
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