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 (40)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1195 additions and 653 deletions
HDF5_SCHEMATA → HDF5_SCHEMATA.md
View file @ 169b1afa
# 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. Python functions for writing
budget data to and reading budget data from this format are included
......@@ -19,53 +18,22 @@ 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
* [0] https://en.wikipedia.org/wiki/Hierarchical_Data_Format
* [1] http://www.h5py.org/
## 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 -->
<!-- ... -->
<!-- ``` -->
brackets (<>). Group objects are indicated by a closing '/', data
sets 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.
## Top-level Budget
## top-level attributes
The following two root attributes expected: a string describing the schema,
and an int schema version number:
......@@ -74,34 +42,81 @@ and an int schema version number:
/.attrs['SCHEMA_VERSION'] = 1
```
Top-level attributes can be used for storing any relevant metadata,
such as "ifo" parameters in YAML format or overall plot styling.
Typically at least the following root attributes are defined:
The following root root attributes are defined:
```
/.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
/.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>
/.attrs['ifo'] = <IFO Struct object, YAML serialized>
```
The remaining root level attributes usually pertain to the plot style.
The budget frequency array is defined in a top level 'Freq' dataset:
```
/'Freq': (N),float
```
## version history
### v1
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>
...
```
The budget frequency array is defined in a 'Freq' dataset:
A budget item, i.e. a collection of noises is structured as follows:
```
/'Freq': (N), float
/<budget>/
/<budget>/Total': (N),float
/<budget>/<trace_0>: (N),float
/<budget>/...
```
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
...
/traces/'Total': (N),float
/traces/<noise_0>: (N),float
/traces/<noise_1>: (N),float
/traces/<noise_2>/
/traces/<noise_2>/'Total': (N),float
/traces/<noise_2>/<noise_3>: (N),float
/traces/<noise_2>/...
/traces/...
```
### v2
Each trace is given the following structure:
```
/<trace>/
/<trace>.attrs['style'] = <YAML trace style dict>
/<trace>/'PSD': (N),float
/<trace>/budget/
/<trace>/budget/<subtrace_0>/
/<trace>/budget/<subtrace_1>/
/<trace>/budget/...
```
The overall structure is:
```
/<budget>/
/<budget>/.attrs['plot_style'] = <YAML plot style dict>
/<budget>/.attrs['style'] = <YAML "Total" trace style dict>
/<budget>/.attrs[*] = <arbitrary data>
/<budget>/'Freq': (N),float
/<budget>/'PSD': (N),float
/<budget>/budget/
/<budget>/budget/<trace_0>/...
/<budget>/budget/<trace_1>/...
/<budget>/budget/...
```
README.md
View file @ 169b1afa
......@@ -39,8 +39,8 @@ data.
The [`inspiral_range`](https://git.ligo.org/gwinc/inspiral-range)
package can be used to calculate various common "inspiral range"
figures of merit for gravitational wave detector budgets. See
[figures of merit](#figures-of-merit) section below.
figures of merit for gravitational wave detector budgets. See the
[inspiral range](#inspiral-range) section below.
## usage
......@@ -48,12 +48,12 @@ figures of merit for gravitational wave detector budgets. See
### command line interface
`pygwinc` provides a command line interface that can be used to
calculate and plot noise budgets for generic noise budgets or the
various canonical IFOs described above, save/plot hdf5 trace data, and
dump budget IFO parameters:
calculate and plot the various canonical IFO noise budgets described
above, as well as custom noise budgets (see below):
```shell
$ python3 -m gwinc aLIGO
```
Budget plots save/plot hdf5 trace data, and dump budget IFO parameters:
You can play with IFO parameters and see the effects on the budget by
dumping the pre-defined parameters to a [YAML-formatted parameter
......@@ -74,7 +74,7 @@ command line:
$ python3 -m gwinc aLIGO --ifo Optics.SRM.Tunephase=3.14
```
Stand-alone YAML files will always assume the nominal ['aLIGO' budget
Stand-alone YAML files assume the nominal ['aLIGO' budget
description](gwinc/ifo/aLIGO).
[Custom budgets](#budget-interface) may also be processed by providing
......@@ -84,18 +84,19 @@ $ python3 -m gwinc path/to/mybudget
```
The command line interface also includes an "interactive" mode which
provides an [IPython](https://ipython.org/) shell for interacting with a processed budget:
provides an [IPython](https://ipython.org/) shell for interacting with
a processed budget:
```shell
$ python3 -m gwinc -i Aplus
GWINC interactive shell
The 'ifo' Struct and 'traces' data are available for inspection.
The 'ifo' Struct and 'trace' data are available for inspection.
Use the 'whos' command to view the workspace.
You may interact with the plot using the 'plt' functions, e.g.:
In [.]: plt.title("foo")
In [.]: plt.savefig("foo.pdf")
In [.]: plt.title("My Special Budget")
In [.]: plt.savefig("mybudget.pdf")
In [1]:
```
......@@ -106,7 +107,7 @@ $ python3 -m gwinc -h
```
### python library
### library interface
For custom plotting, parameter optimization, etc. all functionality can be
accessed directly through the `gwinc` library interface:
......@@ -114,31 +115,32 @@ accessed directly through the `gwinc` library interface:
>>> import gwinc
>>> import numpy as np
>>> freq = np.logspace(1, 3, 1000)
>>> Budget = gwinc.load_budget('aLIGO')
>>> traces = Budget(freq).run()
>>> fig = gwinc.plot_noise(freq, traces)
>>> budget = gwinc.load_budget('aLIGO', freq)
>>> trace = budget.run()
>>> fig = gwinc.plot_budget(trace)
>>> fig.show()
```
The `load_budget()` function takes most of the same inputs as the
command line interface (e.g. IFO names, budget module paths, YAML
parameter files), and returns the un-instantiated `Budget` class
defined in the specified budget module (see [budget
interface](#budget-interface) below).
parameter files), and returns the instantiated `Budget` object defined
in the specified budget module (see [budget
interface](#budget-interface) below). The budget `ifo` `gwinc.Struct`
is available in the `budget.ifo` attribute.
The budget `run()` method is a convenience method that calculates all
budget noises and the noise total and returns a (possibly) nested
dictionary of a noise data, in the form of a `(data, style)` tuple
where 'data' is the PSD data and 'style' is a plot style dictionary
for the trace. The dictionary will be nested if the budget includes
any sub-budgets.
The budget `run()` method calculates all budget noises and the noise
total and returns a `BudgetTrace` object with `freq`, `psd`, and `asd`
properties. The budget sub-traces are available through a dictionary
(`trace['QuantumVacuum']`) interface and via attributes
(`trace.QuantumVacumm`).
## noise functions
`pygwinc` noise functions are available in the `gwinc.noise` package.
This package includes multiple sub-modules for the different types of
noises, e.g. `suspensionthermal`, `coatingthermal`, `quantum`, etc.)
The `pygwinc` analytical noise functions are available in the
`gwinc.noise` package. This package includes multiple sub-modules for
the different types of noises, e.g. `suspensionthermal`,
`coatingthermal`, `quantum`, etc.)
The various noise functions need many different parameters to
calculate their noise outputs. Many parameters are expected to be in
......@@ -161,7 +163,7 @@ def coating_brownian(f, materials, wavelength, wBeam, dOpt):
```
### `gwinc.Struct` objects
## `gwinc.Struct` objects
`pygwinc` provides a `Struct` class that can hold parameters in
attributes and additionally acts like a dictionary, for passing to the
......@@ -232,8 +234,12 @@ are:
* `update(**kwargs)`: update data/attributes
* `calc()`: return final data array
See the built-in documentation for more info (e.g. `pydoc3
gwinc.nb.BudgetItem`)
Generally these methods are not called directly. Instead, the `Noise`
and `Budget` classes include a `run` method that smartly executes them
in sequence and returns a `BudgetTrace` object for the budget.
See the built-in `BudgetItem` documentation for more info
(e.g. `pydoc3 gwinc.nb.BudgetItem`)
### budget module definition
......@@ -313,11 +319,10 @@ The `style` attributes of the various `Noise` classes define plot
style for the noise.
This budget can be loaded with the `gwinc.load_budget()` function, and
processed with the `Budget.run()` method:
processed as usual with the `Budget.run()` method:
```python
Budget = load_budget('/path/to/MyBudget')
budget = Budget(freq)
traces = budget.run()
budget = load_budget('/path/to/MyBudget', freq)
trace = budget.run()
```
Other than the necessary `freq` initialization argument that defines
......@@ -333,24 +338,73 @@ suspension Struct is extracted from the `self.ifo` Struct at
If a budget module defined as a package includes an `ifo.yaml`
[parameter file](#parameter-files) in the package directory, the
`load_budget()` function will automatically load the YAML data into a
`gwinc.Struct` and include it as an `Budget.ifo` attribute in the
returned `Budget` class. This would provide the `self.ifo` needed in
the `SuspensionThermal` Noise class above and is therefore a
convenient way to provide parameter structures in budget packages.
Otherwise it would need to be created/loaded in some other way and
passed to the budget at instantiation, e.g.:
`load_budget()` function will automatically load the YAML data into an
`ifo` `gwinc.Struct` and assign it to the `budget.ifo` attribute.
Alternate ifos can be specified at run time:
```python
Budget = load_budget('/path/to/MyBudget')
budget = load_budget('/path/to/MyBudget', freq)
ifo = Struct.from_file('/path/to/MyBudget.ifo')
budget = Budget(freq, ifo=ifo)
traces = budget.run()
trace = budget.run(ifo=ifo)
```
The IFOs included in `gwinc.ifo` provide examples of the use of the
budget interface (e.g. [gwinc.ifo.aLIGO](gwinc/ifo/aLIGO)).
### the "precomp" decorator
The `BudgetItem` supports "precomp" functions that can be used to
calculate common derived values needed in multiple `BudgetItems`.
They are specified using the `nb.precomp` decorator applied to the
`BudgetItem.calc()` method. These functions are executed during the
`update()` method call, supplied with the budget `freq` and `ifo`
attributes as input arguments, and are expected to update the `ifo`
struct. The execution state of the precomp functions is cached at the
Budget level, to prevent needlessly re-calculating them multiple
times. For example:
```python
from gwinc import nb
def precomp_foo(freq, ifo):
pc = Struct()
...
return pc
def precomp_bar(freq, ifo):
pc = Struct()
...
return pc
class Noise0(nb.Noise):
@nb.precomp(foo=precomp_foo)
def calc(self, foo):
...
class Noise1(nb.Noise):
@nb.precomp(foo=precomp_foo)
@nb.precomp(bar=precomp_bar)
def calc(self, foo, bar):
...
class MyBudget(nb.Budget):
noises = [
Noise0,
Noise1,
]
```
When `MyBudget.run()` is called, all the `precomp` functions will be
executed, e.g.:
```python
precomp_foo(self.freq, self.ifo)
precomp_bar(self.freq, self.ifo)
```
But the `precomp_foo` function will only be calculated once even
though it's specified as needed by both `Noise0` and `Noise1`.
### extracting single noise terms
There are various way to extract single noise terms from the Budget
......@@ -358,10 +412,9 @@ interface. The most straightforward way is to run the full budget,
and extract the noise data the output traces dictionary:
```python
Budget = load_budget('/path/to/MyBudget')
budget = Budget(freq)
traces = budget.calc_traces()
data, plot_style = traces['QuantumVacuum']
budget = load_budget('/path/to/MyBudget', freq)
trace = budget.run()
quantum_trace = trace['QuantumVacuum']
```
You can also calculate the final calibrated output noise for just a
......@@ -379,7 +432,7 @@ data = budget['QuantumVacuum'].calc()
```
# figures of merit
# inspiral range
The [`inspiral_range`](https://git.ligo.org/gwinc/inspiral-range)
package can be used to calculate various common "inspiral range"
......@@ -392,18 +445,15 @@ import inspiral_range
import numpy as np
freq = np.logspace(1, 3, 1000)
Budget = gwinc.load_budget('Aplus')
traces = Budget(freq).run()
budget = gwinc.load_budget('Aplus', freq)
trace = budget.run()
range = inspiral_range.range(
freq, traces['Total'][0],
freq, trace.psd,
m1=30, m2=30,
)
```
Note you need to extract the zeroth element of the `traces['Total']`
tuple, which is the actual PSD data.
See the [`inspiral_range`](https://git.ligo.org/gwinc/inspiral-range)
package for more details.
......
gwinc/__init__.py
View file @ 169b1afa
......@@ -7,6 +7,7 @@ import numpy as np
from .ifo import IFOS
from .struct import Struct
from .plot import plot_budget
from .plot import plot_noise
......@@ -37,24 +38,24 @@ def _load_module(name_or_path):
return mod, path
def load_budget(name_or_path):
"""Load GWINC IFO Budget by name or from path.
def load_budget(name_or_path, freq=None):
"""Load GWINC Budget
Named IFOs should correspond to one of the IFOs available in the
gwinc package (see gwinc.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).
Accepts either the name of a built-in canonical budget (see
gwinc.IFOS), the path to a budget package (directory) or module
(ending in .py), or the path to an IFO struct definition file (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 the budget is a package directory which includes an 'ifo.yaml'
file the ifo Struct will be loaded from that file and assigned to
the budget.ifo attribute. If a struct definition file is provided
the base aLIGO budget definition and ifo will be assumed.
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.
Returns a Budget object instantiated with the provided frequency
array (if specified), and with any included ifo assigned as an
object attribute. If a frequency array is not provided and the
budget class does not define it's own, the frequency array must be
provided at budget update() or run() time.
"""
ifo = None
......@@ -68,11 +69,9 @@ def load_budget(name_or_path):
ifo = Struct.from_file(path)
bname = 'aLIGO'
modname = 'gwinc.ifo.aLIGO'
logger.info("loading budget {}...".format(modname))
else:
modname = path
logger.info("loading module path {}...".format(modname))
else:
if name_or_path not in IFOS:
......@@ -82,17 +81,14 @@ def load_budget(name_or_path):
))
bname = name_or_path
modname = 'gwinc.ifo.'+name_or_path
logger.info("loading module {}...".format(modname))
logger.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 os.path.exists(ifopath):
ifo = Struct.from_file(ifopath)
Budget.ifo = ifo
return Budget
return Budget(freq=freq, ifo=ifo)
def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
......@@ -114,15 +110,16 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
# 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')
traces = Budget(freq, ifo=ifo).run()
plot_style = getattr(Budget, 'plot_style', {})
budget = load_budget('aLIGO', freq)
traces = budget.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
for name, trace in traces.items():
noises[name] = trace.psd
noises['Total'] = traces.psd
noises['Freq'] = traces.freq
pbs = ifo.gwinc.pbs
parm = ifo.gwinc.parm
......@@ -169,6 +166,6 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
logger.info('BBH Inspiral Range: ' + str(score.effr0bh) + ' Mpc/ z = ' + str(score.zHorizonBH))
logger.info('Stochastic Omega: %4.1g Universes' % score.Omega)
plot_noise(ifo, traces, **plot_style)
plot_budget(traces, **plot_style)
return score, noises, ifo
gwinc/__main__.py
View file @ 169b1afa
......@@ -5,10 +5,10 @@ import logging
import argparse
import numpy as np
from . import IFOS, load_budget, plot_noise, logger
from . import IFOS, load_budget, plot_budget, logger
logger.setLevel(os.getenv('LOG_LEVEL', 'WARNING').upper())
formatter = logging.Formatter('%(name)s: %(message)s')
formatter = logging.Formatter('%(message)s')
handler = logging.StreamHandler()
handler.setFormatter(formatter)
logger.addHandler(handler)
......@@ -61,7 +61,7 @@ parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument(
'--freq', '-f', metavar='FLO:[NPOINTS:]FHI',
help="frequency array specification in Hz [{}]".format(FREQ))
help="logarithmic frequency array specification in Hz [{}]".format(FREQ))
parser.add_argument(
'--ifo', '-o', metavar='PARAM=VAL',
#nargs='+', action='extend',
......@@ -71,7 +71,7 @@ parser.add_argument(
'--title', '-t',
help="plot title")
parser.add_argument(
'--fom', metavar='FUNC[:PARAM=VAL,...]',
'--fom', metavar='FUNC[:PARAM=VAL,...]', default=FOM,
help="use inspiral_range.FUNC to calculate range figure-of-merit on resultant spectrum")
group = parser.add_mutually_exclusive_group()
group.add_argument(
......@@ -118,29 +118,28 @@ def main():
if os.path.splitext(os.path.basename(args.IFO))[1] in DATA_SAVE_FORMATS:
if args.freq:
parser.exit(2, "Can not specify frequency array when loading traces from file.\n")
parser.exit(2, "Frequency specification not allowed when loading traces from file.\n")
if args.ifo:
parser.exit(2, "Can not override ifo parameters when loading traces from file.\n")
parser.exit(2, "IFO parameter specification not allowed when loading traces from file.\n")
from .io import load_hdf5
Budget = None
freq, traces, attrs = load_hdf5(args.IFO)
ifo = attrs.get('ifo')
# FIXME: deprecate 'IFO'
ifo = attrs.get('IFO', ifo)
plot_style = attrs
budget = None
trace = load_hdf5(args.IFO)
freq = trace.freq
ifo = trace.ifo
plot_style = trace.plot_style
else:
Budget = load_budget(args.IFO)
ifo = Budget.ifo
budget = load_budget(args.IFO)
ifo = budget.ifo
if args.freq:
try:
freq = freq_from_spec(args.freq)
except IndexError:
parser.error("improper frequency specification '{}'".format(args.freq))
parser.exit(2, "Improper frequency specification: {}\n".format(args.freq))
else:
freq = getattr(Budget, 'freq', freq_from_spec(FREQ))
plot_style = getattr(Budget, 'plot_style', {})
traces = None
freq = getattr(budget, 'freq', freq_from_spec(FREQ))
plot_style = getattr(budget, 'plot_style', {})
trace = None
if args.ifo:
for paramval in args.ifo:
......@@ -149,19 +148,19 @@ def main():
if args.yaml:
if not ifo:
parser.exit(2, "No 'ifo' structure available.\n")
parser.exit(2, "IFO structure not provided.\n")
print(ifo.to_yaml(), end='')
return
if args.text:
if not ifo:
parser.exit(2, "No 'ifo' structure available.\n")
parser.exit(2, "IFO structure not provided.\n")
print(ifo.to_txt(), end='')
return
if args.diff:
if not ifo:
parser.exit(2, "No 'ifo' structure available.\n")
Budget = load_budget(args.diff)
diffs = ifo.diff(Budget.ifo)
parser.exit(2, "IFO structure not provided.\n")
dbudget = load_budget(args.diff)
diffs = ifo.diff(dbudget.ifo)
if diffs:
w = max([len(d[0]) for d in diffs])
fmt = '{{:{}}} {{:>20}} {{:>20}}'.format(w)
......@@ -204,22 +203,21 @@ def main():
import inspiral_range
logger_ir = logging.getLogger('inspiral_range')
logger_ir.setLevel(logger.getEffectiveLevel())
logger_ir.removeHandler(logger_ir.handlers[0])
logger_ir.addHandler(logger.handlers[0])
if not args.fom:
args.fom = FOM
handler = logging.StreamHandler()
handler.setFormatter(logging.Formatter('%(name)s: %(message)s'))
logger_ir.addHandler(handler)
except ModuleNotFoundError:
if args.fom:
raise
logger.warning("inspiral_range package not available, figure of merit will not be calculated...")
logger.warning("WARNING: inspiral_range package not available, figure of merit will not be calculated.")
args.fom = None
if args.fom:
try:
range_func, fargs = args.fom.split(':')
range_func, range_func_args = args.fom.split(':')
except ValueError:
range_func = args.fom
fargs = ''
range_func_args = ''
range_params = {}
for param in fargs.split(','):
for param in range_func_args.split(','):
if not param:
continue
p, v = param.split('=')
......@@ -234,37 +232,33 @@ def main():
##########
# main calculations
if not traces:
if not trace:
logger.info("calculating budget...")
traces = Budget(freq=freq, ifo=ifo).run()
# logger.info('recycling factor: {: >0.3f}'.format(ifo.gwinc.prfactor))
# logger.info('BS power: {: >0.3f} W'.format(ifo.gwinc.pbs))
# logger.info('arm finesse: {: >0.3f}'.format(ifo.gwinc.finesse))
# logger.info('arm power: {: >0.3f} kW'.format(ifo.gwinc.parm/1000))
trace = budget.run(freq=freq)
if args.title:
plot_style['title'] = args.title
elif 'title' in plot_style:
pass
elif Budget:
plot_style['title'] = "GWINC Noise Budget: {}".format(Budget.name)
elif budget:
plot_style['title'] = "GWINC Noise Budget: {}".format(budget.name)
else:
plot_style['title'] = "GWINC Noise Budget: {}".format(args.IFO)
if args.fom:
logger.info("calculating inspiral {}...".format(range_func))
H = inspiral_range.CBCWaveform(freq, **range_params)
logger.debug("params: {}".format(H.params))
fom = eval('inspiral_range.{}'.format(range_func))(freq, traces['Total'][0], H=H)
logger.info("{}({}) = {:.2f} Mpc".format(range_func, fargs, fom))
fom_title = 'inspiral {func} {m1}/{m2} $\mathrm{{M}}_\odot$: {fom:.0f} Mpc'.format(
logger.debug("waveform params: {}".format(H.params))
fom = eval('inspiral_range.{}'.format(range_func))(freq, trace.psd, H=H)
logger.info("{}({}) = {:.2f} Mpc".format(range_func, range_func_args, fom))
subtitle = 'inspiral {func} {m1}/{m2} $\mathrm{{M}}_\odot$: {fom:.0f} Mpc'.format(
func=range_func,
m1=H.params['m1'],
m2=H.params['m2'],
fom=fom,
)
plot_style['title'] += '\n{}'.format(fom_title)
else:
subtitle = None
##########
# interactive
......@@ -273,14 +267,14 @@ def main():
if args.interactive:
banner = """GWINC interactive shell
The 'ifo' Struct and 'traces' data are available for inspection.
Use the 'whos' command to view the workspace.
The 'ifo' Struct and 'budget' trace data objects are available for
inspection. Use the 'whos' command to view the workspace.
"""
if not args.plot:
banner += """
You may plot the budget using the 'plot_noise()' function:
You may plot the budget using the 'plot_budget()' function:
In [.]: plot_noise(freq, traces, **plot_style)
In [.]: plot_budget(budget, **plot_style)
"""
banner += """
You may interact with the plot using the 'plt' functions, e.g.:
......@@ -292,34 +286,31 @@ In [.]: plt.savefig("foo.pdf")
ipshell = InteractiveShellEmbed(
banner1=banner,
user_ns={
'freq': freq,
'traces': traces,
'budget': trace,
'ifo': ifo,
'plot_style': plot_style,
'plot_noise': plot_noise,
'plot_budget': plot_budget,
},
)
ipshell.enable_pylab()
if args.plot:
ipshell.ex("fig = plot_noise(freq, traces, **plot_style)")
ipshell.ex("plt.title('{}')".format(plot_style['title']))
ipshell.ex("fig = plot_budget(budget, **plot_style)")
ipshell.ex("plt.title(plot_style['title'])")
ipshell()
##########
# output
# save noise traces to HDF5 file
# save noise trace to HDF5 file
if out_data_files:
from .io import save_hdf5
attrs = dict(plot_style)
for path in out_data_files:
logger.info("saving budget traces: {}".format(path))
logger.info("saving budget trace: {}".format(path))
save_hdf5(
path=path,
freq=freq,
traces=traces,
trace=trace,
ifo=ifo,
**attrs,
plot_style=plot_style,
)
# standard plotting
......@@ -327,10 +318,10 @@ In [.]: plt.savefig("foo.pdf")
logger.debug("plotting noises...")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plot_noise(
freq,
traces,
plot_budget(
trace,
ax=ax,
subtitle=subtitle,
**plot_style
)
fig.tight_layout()
......
gwinc/ifo/CE2/__init__.py
View file @ 169b1afa
......@@ -54,7 +54,6 @@ class Substrate(nb.Budget):
noises = [
ITMThermoRefractive,
ITMCarrierDensity,
SubstrateBrownian,
SubstrateThermoElastic,
]
......
gwinc/ifo/CE2/ifo.yaml
View file @ 169b1afa
......@@ -227,15 +227,6 @@ Materials:
RefractiveIndex: 3.5 # 3.38 * (1 + 4e-5 * T) (ioffe)
dndT: 1e-4 # ~123K & 1900 nm : http://arxiv.org/abs/physics/0606168
Temp: 123 # mirror temperature [K]
## parameters for semiconductor optics
isSemiConductor: True # we are doing semiconductor optics
CarrierDensity: 1e19 # 1/m^3; carrier density for phosphorous-doped silicon
ElectronDiffusion: 9.7e-3 # m^2/s; electron diffusion coefficient for silicon at 120 K
HoleDiffusion: 3.5e-3 # m^2/s; hole diffusion coefficient for silicon at 120 K
ElectronEffMass: 9.747e-31 # kg; effective mass of each electron 1.07*m_e
HoleEffMass: 8.016e-31 # kg; effective mass of each hole 0.88*m_e
ElectronIndexGamma: -8.8e-28 # m**3; dependence of index of refraction on electron carrier density
HoleIndexGamma: -1.02e-27 # m**3; dependence of index of refraction on hole carrier density
MassRadius: 0.4 # m; 80 cm mCZ silicon
MassThickness: 0.286
......@@ -257,7 +248,6 @@ Optics:
# zz = loadmat('CryogenicLIGO/Sensitivity/GWINC/optRuns/' + qopt_mat)
ITM:
SubstrateAbsorption: 1e-3 # 1/m; 10 ppm/cm for MCZ Si
BeamRadius: 0.141 # m; 1/e^2 power radius w1
CoatingAbsorption: 0.5e-6 # absorption of ITM
Transmittance: 0.014 # zz['x'][0][3]
#CoatingThicknessLown: 0.308
......@@ -278,7 +268,6 @@ Optics:
- 0.3867382
- 0.08814237
ETM:
BeamRadius: 0.141 # m; 1/e^2 power radius w2
Transmittance: 5e-6 # Transmittance of ETM
#CoatingThicknessLown: 0.27
#CoatingThicknessCap: 0.5
......
gwinc/ifo/Voyager/__init__.py
View file @ 169b1afa
......@@ -14,7 +14,6 @@ class Voyager(nb.Budget):
CoatingBrownian,
CoatingThermoOptic,
ITMThermoRefractive,
ITMCarrierDensity,
SubstrateBrownian,
SubstrateThermoElastic,
ExcessGas,
......
gwinc/ifo/Voyager/ifo.yaml
View file @ 169b1afa
......@@ -252,15 +252,6 @@ Materials:
RefractiveIndex: 3.5 # 3.38 * (1 + 4e-5 * T) (ioffe)
dndT: 1e-4 # ~123K & 1900 nm : http://arxiv.org/abs/physics/0606168
Temp: 123 # mirror temperature [K]
## parameters for semiconductor optics
isSemiConductor: True # we are doing semiconductor optics
CarrierDensity: 1e19 # 1/m^3; carrier density for phosphorous-doped silicon
ElectronDiffusion: 9.7e-3 # m^2/s; electron diffusion coefficient for silicon at 120 K
HoleDiffusion: 3.5e-3 # m^2/s; hole diffusion coefficient for silicon at 120 K
ElectronEffMass: 9.747e-31 # kg; effective mass of each electron 1.07*m_e
HoleEffMass: 8.016e-31 # kg; effective mass of each hole 0.88*m_e
ElectronIndexGamma: -8.8e-28 # m**3; dependence of index of refraction on electron carrier density
HoleIndexGamma: -1.02e-27 # m**3; dependence of index of refraction on hole carrier density
MassRadius: 0.225 # m; 45 cm mCZ silicon
MassThickness: 0.55
......@@ -282,7 +273,6 @@ Optics:
# zz = loadmat('CryogenicLIGO/Sensitivity/GWINC/optRuns/' + qopt_mat)
ITM:
SubstrateAbsorption: 1e-3 # 1/m; 10 ppm/cm for MCZ Si
BeamRadius: 0.0585 # m; 1/e^2 power radius w1
CoatingAbsorption: 1e-6 # absorption of ITM
Transmittance: 2e-3 # zz['x'][0][3]
#CoatingThicknessLown: 0.308
......@@ -303,7 +293,6 @@ Optics:
- 0.386738199718736
- 0.088142365969070
ETM:
BeamRadius: 0.0835 # m; 1/e^2 power radius w2
Transmittance: 5e-6 # Transmittance of ETM
#CoatingThicknessLown: 0.27
#CoatingThicknessCap: 0.5
......
gwinc/ifo/__main__.py
View file @ 169b1afa
......@@ -2,8 +2,6 @@ import argparse
import numpy as np
import matplotlib.pyplot as plt
import inspiral_range
from . import IFOS, PLOT_STYLE
from .. import load_budget
......@@ -31,20 +29,27 @@ def main():
budgets = {}
range_pad = 0
for ifo in IFOS:
budget = load_budget(ifo)(freq)
budget = load_budget(ifo, freq)
name = budget.name
budgets[name] = budget
range_pad = max(len(name), range_pad)
for name, budget in budgets.items():
data = budget.calc()
BNS_range = inspiral_range.range(freq, data)
label = '{name:<{pad}} {bns:>6.0f} Mpc'.format(
trace = budget.run()
try:
import inspiral_range
label_range = ' {:>6.0f} Mpc'.format(
inspiral_range.range(freq, trace.psd),
)
except ModuleNotFoundError:
label_range = ''
label = '{name:<{pad}}{range}'.format(
name=name,
pad=range_pad,
bns=BNS_range,
range=label_range,
)
ax.loglog(freq, np.sqrt(data), label=label)
ax.loglog(freq, trace.asd, label=label, lw=2)
ax.grid(
True,
......
gwinc/ifo/aLIGO/ifo.yaml
View file @ 169b1afa
......@@ -188,7 +188,7 @@ Suspension:
## Optic Material -------------------------------------------------------
Materials:
MassRadius: 0.17 # m;
MassRadius: 0.17 # m;
MassThickness: 0.200 # m; Peter F 8/11/2005
## Dielectric coating material parameters----------------------------------
......@@ -251,13 +251,11 @@ Optics:
dc: 1.5707963 # pi/2 # demod/detection/homodyne phase
ITM:
# BeamRadius: 0.055 # m, 1/e^2 power radius, now in precompIFO
Transmittance: 0.014
CoatingThicknessLown: 0.308
CoatingThicknessCap: 0.5
CoatingAbsorption: 0.5e-6
ETM:
# BeamRadius: 0.062 # m, 1/e^2 power radius, now in precompIFO
Transmittance: 5e-6
CoatingThicknessLown: 0.27
CoatingThicknessCap: 0.5
......
gwinc/ifo/noises.py
View file @ 169b1afa
import copy
import numpy as np
from numpy import pi, sin, exp, sqrt
......@@ -10,14 +11,33 @@ from .. import suspension
##################################################
def coating_thickness(ifo, optic):
optic = ifo.Optics.get(optic)
def coating_thickness(materials, optic):
if 'CoatLayerOpticalThickness' in optic:
return optic['CoatLayerOpticalThickness']
T = optic.Transmittance
dL = optic.CoatingThicknessLown
dCap = optic.CoatingThicknessCap
return noise.coatingthermal.getCoatDopt(ifo.Materials, T, dL, dCap=dCap)
return noise.coatingthermal.getCoatDopt(materials, T, dL, dCap=dCap)
def mirror_struct(ifo, tm):
"""Create a "mirror" Struct for a LIGO core optic
This is a copy of the ifo.Materials Struct, containing Substrate
and Coating sub-Structs, as well as some basic geometrical
properties of the optic.
"""
# NOTE: we deepcopy this struct since we'll be modifying it (and
# it would otherwise be stored by reference)
mirror = copy.deepcopy(ifo.Materials)
optic = ifo.Optics.get(tm)
mirror.Coating.dOpt = coating_thickness(ifo.Materials, optic)
mirror.update(optic)
mirror.MassVolume = pi * mirror.MassRadius**2 * mirror.MassThickness
mirror.MirrorMass = mirror.MassVolume * mirror.Substrate.MassDensity
return mirror
def arm_cavity(ifo):
......@@ -51,7 +71,6 @@ def ifo_power(ifo, PRfixed=True):
"""Compute power on beamsplitter, finesse, and power recycling factor.
"""
c = const.c
pin = ifo.Laser.Power
t1 = sqrt(ifo.Optics.ITM.Transmittance)
r1 = sqrt(1 - ifo.Optics.ITM.Transmittance)
......@@ -65,7 +84,7 @@ def ifo_power(ifo, PRfixed=True):
# Finesse, effective number of bounces in cavity, power recycling factor
finesse = 2*pi / (t1**2 + 2*loss) # arm cavity finesse
neff = 2 * finesse / pi
neff = 2 * finesse / pi
# Arm cavity reflectivity with finite loss
garm = t1 / (1 - r1*r2*sqrt(1-2*loss)) # amplitude gain wrt input field
......@@ -79,12 +98,12 @@ def ifo_power(ifo, PRfixed=True):
r5 = sqrt(1 - Tpr)
prfactor = t5**2 / (1 + r5 * rarm * sqrt(1-bsloss))**2
pbs = pin * prfactor # BS power from input power
pbs = pin * prfactor # BS power from input power
parm = pbs * garm**2 / 2 # arm power from BS power
thickness = ifo.Optics.ITM.get('Thickness', ifo.Materials.MassThickness)
asub = 1.3 * 2 * thickness * ifo.Optics.SubstrateAbsorption
pbsl = 2 *pcrit / (asub+acoat*neff) # bs power limited by thermal lensing
pbsl = 2 * pcrit / (asub+acoat*neff) # bs power limited by thermal lensing
if pbs > pbsl:
logger.warning('P_BS exceeds BS Thermal limit!')
......@@ -119,34 +138,46 @@ def dhdl(f, armlen):
##################################################
def precomp_mirror(f, ifo):
ifo.Materials.MirrorVolume = \
pi * ifo.Materials.MassRadius**2 \
* ifo.Materials.MassThickness
ifo.Materials.MirrorMass = \
ifo.Materials.MirrorVolume \
* ifo.Materials.Substrate.MassDensity
def precomp_suspension(f, ifo):
pc = Struct()
pc.VHCoupling = Struct()
if 'VHCoupling' in ifo.Suspension:
pc.VHCoupling.theta = ifo.Suspension.VHCoupling.theta
else:
pc.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension)
pc.hForce = hForce
pc.vForce = vForce
pc.hTable = hTable
pc.vTable = vTable
return pc
def precomp_gwinc(f, ifo):
ifo.gwinc = Struct()
pbs, parm, finesse, prfactor, Tpr = ifo_power(ifo)
ifo.gwinc.pbs = pbs
ifo.gwinc.parm = parm
ifo.gwinc.finesse = finesse
ifo.gwinc.prfactor = prfactor
def precomp_mirror(f, ifo):
pc = Struct()
pc.ITM = mirror_struct(ifo, 'ITM')
pc.ETM = mirror_struct(ifo, 'ETM')
return pc
def precomp_suspension(f, ifo):
if 'VHCoupling' not in ifo.Suspension:
ifo.Suspension.VHCoupling = Struct()
ifo.Suspension.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension)
ifo.Suspension.hForce = hForce
ifo.Suspension.vForce = vForce
ifo.Suspension.hTable = hTable
ifo.Suspension.vTable = vTable
def precomp_power(f, ifo):
pc = Struct()
pbs, parm, finesse, prfactor, Tpr = ifo_power(ifo)
pc.pbs = pbs
pc.parm = parm
pc.finesse = finesse
pc.gPhase = finesse * 2/np.pi
pc.prfactor = prfactor
return pc
def precomp_cavity(f, ifo):
pc = Struct()
w0, wBeam_ITM, wBeam_ETM = arm_cavity(ifo)
pc.w0 = w0
pc.wBeam_ITM = wBeam_ITM
pc.wBeam_ETM = wBeam_ETM
return pc
##################################################
......@@ -166,10 +197,10 @@ class QuantumVacuum(nb.Noise):
color='#ad03de',
)
def calc(self):
precomp_mirror(self.freq, self.ifo)
precomp_gwinc(self.freq, self.ifo)
return noise.quantum.shotrad(self.freq, self.ifo)
@nb.precomp(mirror=precomp_mirror)
@nb.precomp(power=precomp_power)
def calc(self, mirror, power):
return noise.quantum.shotrad(self.freq, self.ifo, mirror.ETM.MirrorMass, power)
class StandardQuantumLimit(nb.Noise):
......@@ -179,11 +210,12 @@ class StandardQuantumLimit(nb.Noise):
style = dict(
label="Standard Quantum Limit",
color="#000000",
linestyle=":", # Dotted.
linestyle=":",
)
def calc(self):
return 8 * const.hbar / (self.ifo.Materials.MirrorMass * (2 * np.pi * self.freq) ** 2)
@nb.precomp(mirror=precomp_mirror)
def calc(self, mirror):
return 8 * const.hbar / (mirror.ETM.MirrorMass * (2 * np.pi * self.freq) ** 2)
class Seismic(nb.Noise):
......@@ -195,8 +227,8 @@ class Seismic(nb.Noise):
color='#855700',
)
def calc(self):
precomp_suspension(self.freq, self.ifo)
@nb.precomp(sustf=precomp_suspension)
def calc(self, sustf):
if 'PlatformMotion' in self.ifo.Seismic:
if self.ifo.Seismic.PlatformMotion == 'BSC':
nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
......@@ -207,7 +239,7 @@ class Seismic(nb.Noise):
else:
nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
n, nh, nv = noise.seismic.seismic_suspension_fitered(
self.ifo.Suspension, nt)
self.ifo.Suspension, sustf, nt)
return n * 4
......@@ -277,9 +309,10 @@ class SuspensionThermal(nb.Noise):
color='#0d75f8',
)
def calc(self):
precomp_suspension(self.freq, self.ifo)
n = noise.suspensionthermal.suspension_thermal(self.freq, self.ifo.Suspension)
@nb.precomp(sustf=precomp_suspension)
def calc(self, sustf):
n = noise.suspensionthermal.suspension_thermal(
self.freq, self.ifo.Suspension, sustf)
return n * 4
......@@ -292,16 +325,17 @@ class CoatingBrownian(nb.Noise):
color='#fe0002',
)
def calc(self):
@nb.precomp(mirror=precomp_mirror)
@nb.precomp(cavity=precomp_cavity)
@nb.precomp(power=precomp_power)
def calc(self, mirror, cavity, power):
wavelength = self.ifo.Laser.Wavelength
materials = self.ifo.Materials
w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
dOpt_ITM = coating_thickness(self.ifo, 'ITM')
dOpt_ETM = coating_thickness(self.ifo, 'ETM')
nITM = noise.coatingthermal.coating_brownian(
self.freq, materials, wavelength, wBeam_ITM, dOpt_ITM)
self.freq, mirror.ITM, wavelength, cavity.wBeam_ITM, power.parm,
)
nETM = noise.coatingthermal.coating_brownian(
self.freq, materials, wavelength, wBeam_ETM, dOpt_ETM)
self.freq, mirror.ETM, wavelength, cavity.wBeam_ETM, power.parm,
)
return (nITM + nETM) * 2
......@@ -315,16 +349,17 @@ class CoatingThermoOptic(nb.Noise):
linestyle='--',
)
def calc(self):
@nb.precomp(mirror=precomp_mirror)
@nb.precomp(cavity=precomp_cavity)
def calc(self, mirror, cavity):
wavelength = self.ifo.Laser.Wavelength
materials = self.ifo.Materials
w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
dOpt_ITM = coating_thickness(self.ifo, 'ITM')
dOpt_ETM = coating_thickness(self.ifo, 'ETM')
nITM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic(
self.freq, materials, wavelength, wBeam_ITM, dOpt_ITM[:])
self.freq, mirror.ITM, wavelength, cavity.wBeam_ITM,
)
nETM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic(
self.freq, materials, wavelength, wBeam_ETM, dOpt_ETM[:])
self.freq, mirror.ETM, wavelength, cavity.wBeam_ETM,
)
return (nITM + nETM) * 2
......@@ -338,32 +373,12 @@ class ITMThermoRefractive(nb.Noise):
linestyle='--',
)
def calc(self):
finesse = ifo_power(self.ifo)[2]
gPhase = finesse * 2/np.pi
w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
@nb.precomp(cavity=precomp_cavity)
@nb.precomp(power=precomp_power)
def calc(self, cavity, power):
n = noise.substratethermal.substrate_thermorefractive(
self.freq, self.ifo.Materials, wBeam_ITM)
return n * 2 / gPhase**2
class ITMCarrierDensity(nb.Noise):
"""ITM Carrier Density
"""
style = dict(
label='ITM Carrier Density',
color='#929591',
linestyle='--',
)
def calc(self):
finesse = ifo_power(self.ifo)[2]
gPhase = finesse * 2/np.pi
w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
n = noise.substratethermal.substrate_carrierdensity(
self.freq, self.ifo.Materials, wBeam_ITM)
return n * 2 / gPhase**2
self.freq, self.ifo.Materials, cavity.wBeam_ITM)
return n * 2 / power.gPhase**2
class SubstrateBrownian(nb.Noise):
......@@ -376,12 +391,12 @@ class SubstrateBrownian(nb.Noise):
linestyle='--',
)
def calc(self):
w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
@nb.precomp(cavity=precomp_cavity)
def calc(self, cavity):
nITM = noise.substratethermal.substrate_brownian(
self.freq, self.ifo.Materials, wBeam_ITM)
self.freq, self.ifo.Materials, cavity.wBeam_ITM)
nETM = noise.substratethermal.substrate_brownian(
self.freq, self.ifo.Materials, wBeam_ETM)
self.freq, self.ifo.Materials, cavity.wBeam_ETM)
return (nITM + nETM) * 2
......@@ -395,12 +410,12 @@ class SubstrateThermoElastic(nb.Noise):
linestyle='--',
)
def calc(self):
w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
@nb.precomp(cavity=precomp_cavity)
def calc(self, cavity):
nITM = noise.substratethermal.substrate_thermoelastic(
self.freq, self.ifo.Materials, wBeam_ITM)
self.freq, self.ifo.Materials, cavity.wBeam_ITM)
nETM = noise.substratethermal.substrate_thermoelastic(
self.freq, self.ifo.Materials, wBeam_ETM)
self.freq, self.ifo.Materials, cavity.wBeam_ETM)
return (nITM + nETM) * 2
......
gwinc/io.py
View file @ 169b1afa
......@@ -4,25 +4,23 @@ import datetime
from . import logger
from . import Struct
from .trace import BudgetTrace
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 _write_trace_recursive(f, trace):
f.attrs['style'] = yaml.dump(trace.style)
f.create_dataset('PSD', data=trace.psd)
if trace.budget:
bgrp = f.create_group('budget')
for t in trace:
tgrp = bgrp.create_group(t.name)
_write_trace_recursive(tgrp, t)
def save_hdf5(path, freq, traces, **kwargs):
def save_hdf5(path, trace, ifo=None, plot_style=None, **kwargs):
"""Save GWINC budget data to an HDF5 file.
The `freq` argument should be the frequency array, and `traces`
......@@ -36,46 +34,107 @@ def save_hdf5(path, freq, traces, **kwargs):
"""
with h5py.File(path, 'w') as f:
f.attrs['SCHEMA'] = SCHEMA
f.attrs['SCHEMA_VERSION'] = SCHEMA_VERSION
f.attrs['SCHEMA_VERSION'] = 2
# FIXME: add budget code hash or something
f.attrs['date'] = datetime.datetime.now().isoformat()
if ifo:
f.attrs['ifo'] = ifo.to_yaml()
if plot_style:
f.attrs['plot_style'] = yaml.dump(plot_style)
for key, val in kwargs.items():
if key == 'ifo':
f.attrs['ifo'] = val.to_yaml()
f.attrs[key] = val
f.create_dataset('Freq', data=trace.freq)
tgrp = f.create_group(trace.name)
_write_trace_recursive(tgrp, trace)
def _load_hdf5_v1(f):
def read_recursive(element):
budget = []
for name, item in element.items():
if name == 'Total':
total = item[:]
continue
if isinstance(item, h5py.Group):
trace = read_recursive(item)
trace.name = name
else:
f.attrs[key] = val
f.create_dataset('Freq', data=freq)
tgrp = f.create_group('traces')
_write_trace_recursive(tgrp, traces)
trace = BudgetTrace(
name=name,
style=dict(item.attrs),
psd=item[:],
budget=[],
)
budget.append(trace)
return BudgetTrace(
psd=total,
budget=budget,
)
trace = read_recursive(f['/traces'])
trace.name = 'Total'
trace._freq = f['Freq'][:]
attrs = dict(f.attrs)
if 'ifo' in attrs:
try:
ifo = Struct.from_yaml(attrs['ifo'])
del attrs['ifo']
except yaml.constructor.ConstructorError:
logger.warning("HDF5 load warning: Could not de-serialize 'ifo' YAML attribute.")
ifo = None
trace.style = attrs
trace.ifo = ifo
return trace
def _read_trace_recursive(element):
trace = {}
for name, item in element.items():
if isinstance(item, h5py.Group):
trace[name] = _read_trace_recursive(item)
def _load_hdf5_v2(f):
def read_recursive(element):
psd = element['PSD'][:]
style = yaml.safe_load(element.attrs['style'])
budget = []
if 'budget' in element:
for name, item in element['budget'].items():
trace = read_recursive(item)
trace.name = name
budget.append(trace)
return BudgetTrace(
psd=psd,
budget=budget,
style=style,
)
for name, item in f.items():
if name == 'Freq':
freq = item[:]
else:
trace[name] = item[:], dict(item.attrs.items())
trace = read_recursive(item)
trace.name = name
trace._freq = freq
trace.ifo = None
attrs = dict(f.attrs)
ifo = attrs.get('ifo')
if ifo:
try:
trace.ifo = Struct.from_yaml(ifo)
del attrs['ifo']
except yaml.constructor.ConstructorError:
logger.warning("HDF5 load warning: Could not de-serialize 'ifo' YAML attribute.")
trace.plot_style = yaml.safe_load(attrs['plot_style'])
return trace
def load_hdf5(path):
"""Load GWINC budget data from an HDF5 file.
Returns (freq, traces, attrs). An 'ifo' attr will be
de-serialized from YAML into a Struct object.
Returns BudgetTrace. An 'ifo' attr will be de-serialized from
YAML into a Struct object and assigned as an attribute to the
BudgetTrace.
See HDF5_SCHEMA.
"""
loaders = {
1: _load_hdf5_v1,
2: _load_hdf5_v2,
}
with h5py.File(path, 'r') as f:
# FIXME: check SCHEMA name/version
freq = f['Freq'][:]
traces = _read_trace_recursive(f['/traces'])
attrs = dict(f.attrs)
if 'ifo' in attrs:
try:
attrs['ifo'] = Struct.from_yaml(attrs['ifo'])
except yaml.constructor.ConstructorError:
logger.warning("HDF5 load warning: Could not de-serialize 'ifo' YAML attribute.")
return freq, traces, attrs
version = f.attrs.get('SCHEMA_VERSION', 1)
return loaders[version](f)
gwinc/nb.py
View file @ 169b1afa
import operator
import itertools
import collections
import numpy as np
import scipy.interpolate
from functools import reduce
from . import logger
from .trace import BudgetTrace
def precomp(*precomp_funcs, **precomp_fmaps):
"""BudgetItem.calc decorator to add pre-computed functions
This is intended to decorate BudgetItem.calc() methods with common
functions whose return value are cached for later use. The
functions are supplied with the `freq` array and `ifo` Struct
attributes as arguments, and their return values are provided as
keyword arguments to calc().
For example, if a calc method is defined as:
@precomp(foo=precomp_foo)
@precomp(bar=precomp_bar)
def calc(self, foo, bar):
...
then when the budget is run the precomp functions will be executed
before calc(), roughly the equivalent of:
foo = precomp_foo(self.freq, self.ifo)
bar = precomp_bar(self.freq, self.ifo)
psd = calc(foo=foo, bar=bar)
The execution state of each precomp function is cached so that the
same function is not needlessly executed multiple times.
"""
def decorator(func):
if precomp_funcs:
try:
func._precomp_list.extend(precomp_funcs)
except AttributeError:
func._precomp_list = list(precomp_funcs)
if precomp_fmaps:
try:
func._precomp_mapped.update(precomp_fmaps)
except AttributeError:
func._precomp_mapped = dict(precomp_fmaps)
return func
return decorator
def _precomp_recurse_mapping(func, freq, ifo, _precomp):
"""Recursively execute @precomp decorator functions
Recurses down functions which may themselves have precomp
decorators, building a **kwargs mapping to pass to the wrapped
function call.
"""
for pc_func in itertools.chain(
getattr(func, '_precomp_list', []),
getattr(func, '_precomp_mapped', {}).values(),
):
if pc_func in _precomp:
continue
pc_map = _precomp_recurse_mapping(pc_func, freq, ifo, _precomp=_precomp)
logger.debug("precomp {}".format(pc_func))
_precomp[pc_func] = pc_func(freq, ifo, **pc_map)
kwargs = {
name: _precomp[pc_func] for name, pc_func in getattr(func, '_precomp_mapped', {}).items()
}
return kwargs
def quadsum(data):
......@@ -20,7 +86,6 @@ def quadsum(data):
return np.nansum(data, 0)
class BudgetItem:
"""GWINC BudgetItem class
......@@ -37,6 +102,11 @@ class BudgetItem:
By default any keyword arguments provided are written directly
as attribute variables (as with __init__).
When overloading this method it is recommended to execute the
method from the base class with e.g.:
super().update(**kwargs)
"""
for key, val in kwargs.items():
setattr(self, key, val)
......@@ -50,6 +120,12 @@ class BudgetItem:
"""
return None
def _calc(self, _precomp=None):
"""internal function to call calc with precomp evaluation"""
pcmap = _precomp_recurse_mapping(self.calc, self.freq, self.ifo, _precomp)
logger.debug("calc {}".format(self.name))
return self.calc(**pcmap)
##########
def __init__(self, freq=None, **kwargs):
......@@ -66,10 +142,10 @@ class BudgetItem:
if freq is not None:
assert isinstance(freq, np.ndarray)
self.freq = freq
elif not hasattr(self, 'freq'):
raise AttributeError("Frequency array not provided or defined.")
for key, val in kwargs.items():
setattr(self, key, val)
self._loaded = False
self._precomp = dict()
@property
def name(self):
......@@ -115,7 +191,7 @@ class Calibration(BudgetItem):
e.g. point-by-point product of data and calibration arrays.
"""
cal = self.calc()
cal = self._calc()
assert data.shape == cal.shape, \
"data shape does not match calibration ({} != {})".format(data.shape, cal.shape)
return data * cal
......@@ -134,38 +210,75 @@ class Noise(BudgetItem):
style = {}
"""Trace plot style dictionary"""
def calc_trace(self, calibrations=None, calc=True):
"""Returns noise (PSD, style) tuple.
def _make_trace(self, psd=None, budget=None):
return BudgetTrace(
name=self.name,
style=self.style,
freq=self.freq,
psd=psd,
budget=budget,
)
def calc_trace(self, calibration=1, calc=True, _precomp=None):
"""Calculate noise and return BudgetTrace object
If `calibrations` is not None it is assumed to be a list of
len(self.freq) arrays that will be multiplied to the output
PSD.
`calibration` should either be a scalar or a len(self.freq)
array that will be multiplied to the output PSD of the budget
and all sub noises.
If calc=False, the noise will not be calculated and the PSD
will be None. This is useful for just getting the style.
If `calc` is False, the noise will not be calculated and the
trace PSD will be None. This is useful for just getting the
trace style info.
"""
if _precomp is None:
_precomp = dict()
total = None
if calc:
data = self.calc()
if calibrations:
data *= reduce(
operator.mul,
calibrations,
)
else:
data = None
return data, self.style
total = self._calc(_precomp) * calibration
return self._make_trace(psd=total)
def run(self, **kwargs):
"""Convenience method to load, update, and return calc traces.
"""Run full budget and return BudgetTrace object.
Roughly the equivalent running load(), update(), and
calc_trace() in sequence. Keyword arguments are passed to the
update() method.
NOTE: The load status is cached such that subsequent calls to
this method will not re-execute the load() method.
Equivalent of load(), update(), calc_traces() run in sequence.
Keyword arguments are passed to update().
NOTE: The update() method is only run if keyword arguments
(`kwargs`) are supplied, or if the `ifo` attribute has
changed.
"""
self.load()
self.update(**kwargs)
return self.calc_trace()
if not self._loaded:
self.load()
self._loaded = True
ifo = kwargs.get('ifo', getattr(self, 'ifo'))
if ifo:
if not hasattr(ifo, '_orig_keys'):
logger.debug("new ifo detected")
ifo._orig_keys = tuple(k for k, v in ifo.walk())
ifo_hash = ifo.hash()
kwargs['ifo'] = ifo
else:
ifo_hash = ifo.hash(ifo._orig_keys)
if ifo_hash != getattr(self, '_ifo_hash', 0):
logger.debug("ifo hash change")
kwargs['ifo'] = self.ifo
self._ifo_hash = ifo_hash
if kwargs:
self.update(**kwargs)
# clear precomp cache
self._precomp = dict()
return self.calc_trace(_precomp=self._precomp)
class Budget(Noise):
......@@ -361,95 +474,75 @@ class Budget(Noise):
logger.debug("load {}".format(item))
item.load()
def update(self, **kwargs):
"""Update all noise and cal objects with supplied kwargs."""
def update(self, _precomp=None, **kwargs):
"""Recursively update all noise and cal objects with supplied kwargs.
See BudgetItem.update() for more info.
"""
for key, val in kwargs.items():
setattr(self, key, val)
for name, item in itertools.chain(
self._cal_objs.items(),
self._noise_objs.items()):
logger.debug("update {}".format(item))
item.update(**kwargs)
def calc_noise(self, name):
"""Return calibrated individual noise term.
def calc_noise(self, name, calibration=1, calc=True, _cals=None, _precomp=None):
"""Return calibrated individual noise BudgetTrace.
The noise PSD and calibration transfer functions are
calculated, and the calibrated noise array is returned.
The noise and calibration transfer functions are calculated,
and the calibrated noise BudgetTrace is returned.
`calibration` is an overall calculated calibration to apply to
the noise.
"""
if _precomp is None:
_precomp = dict()
for cal in self._noise_cals[name]:
if _cals:
calibration *= _cals[cal]
else:
obj = self._cal_objs[cal]
calibration *= obj._calc(_precomp)
noise = self._noise_objs[name]
cals = [self._cal_objs[cal].calc() for cal in self._noise_cals[name]]
data = noise.calc_trace(cals)
if isinstance(data, dict):
return data['Total'][0]
else:
return data[0]
def calc(self):
"""Calculate sum of all noises.
"""
data = [self.calc_noise(name) for name in self._noise_objs.keys() if name in self._budget_noises]
return quadsum(data)
def calc_trace(self, calibrations=None, calc=True):
"""Returns a dictionary of noises traces, keyed by noise names.
return noise.calc_trace(
calibration=calibration,
calc=calc,
_precomp=_precomp,
)
Values are (data, style) trace tuples (see Noise.calc_trace).
The key of the budget sum total is 'Total'. The values of sub
budgets are themselves dictionaries returned from
calc_trace() of the sub budget.
def calc_trace(self, calibration=1, calc=True, _precomp=None):
"""Calculate all budget noises and return BudgetTrace object
If `calibrations` is not None it is assumed to be a list of
len(self.freq) array that will be multiplied to the output PSD
of the budget and all sub noises.
`calibration` should either be a scalar or a len(self.freq)
array that will be multiplied to the output PSD of the budget
and all sub noises.
If calc=False, the noise will not be calculated and the PSD
will be None. This is useful for just getting style the
style.
If `calc` is False, the noise will not be calculated and the
trace PSD will be None. This is useful for just getting the
trace style info.
"""
# start by creating an empty OrderedDict used for outputing trace data
# or style info with the following order:
# references
# total
# constituents
d = collections.OrderedDict()
# allocate references
for name, noise in self._noise_objs.items():
if name in self._budget_noises:
continue
d[name] = noise.calc_trace(calc=False)
# allocate total
if self._budget_noises:
d['Total'] = None, self.style
# allocate constituent
for name, noise in self._noise_objs.items():
if name not in self._budget_noises:
continue
d[name] = noise.calc_trace(calc=False)
# if we're not calc'ing, just return the dict now
if not calc:
return d
# calc all calibrations
c = {}
for name, cal in self._cal_objs.items():
c[name] = cal.calc()
# calc all noises
for name, noise in self._noise_objs.items():
cals = [c[cal] for cal in self._noise_cals[name]]
if calibrations:
cals += calibrations
d[name] = noise.calc_trace(
calibrations=cals,
if _precomp is None:
_precomp = dict()
_cals = {}
if calc:
for name, cal in self._cal_objs.items():
_cals[name] = cal._calc(_precomp)
budget = []
for name in self._noise_objs:
trace = self.calc_noise(
name,
calibration=calibration,
calc=calc,
_cals=_cals,
_precomp=_precomp,
)
# calc budget total
constituent_data = []
for name in self._budget_noises:
if isinstance(d[name], dict):
data = d[name]['Total'][0]
else:
data = d[name][0]
constituent_data.append(data)
d['Total'] = quadsum(constituent_data), self.style
return d
budget.append(trace)
total = quadsum([trace.psd for trace in budget])
return self._make_trace(
psd=total, budget=budget
)
gwinc/noise/coatingthermal.py
View file @ 169b1afa
This diff is collapsed.
gwinc/noise/quantum.py
View file @ 169b1afa
......@@ -17,11 +17,12 @@ def sqzOptimalSqueezeAngle(Mifo, eta):
return alpha
def shotrad(f, ifo):
def shotrad(f, ifo, mirror_mass, power):
"""Quantum noise noise strain spectrum
:f: frequency array in Hz
:ifo: gwinc IFO structure
:ifo: gwinc IFO Struct
:power: gwinc power Struct
:returns: strain noise power spectrum at :f:
......@@ -47,7 +48,7 @@ def shotrad(f, ifo):
else:
namespace = globals()
fname = namespace['shotrad' + ifo.Optics.Type]
coeff, Mifo, Msig, Mn = fname(f, ifo)
coeff, Mifo, Msig, Mn = fname(f, ifo, mirror_mass, power)
# check for consistent dimensions
Nfield = Msig.shape[0]
......@@ -290,7 +291,7 @@ def compile_SEC_RES_TF():
print('RES', '=', str(SEC_RES_expr[0]).replace('Matrix', 'np.array'))
def shotradSignalRecycled(f, ifo):
def shotradSignalRecycled(f, ifo, mirror_mass, power):
"""Quantum noise model for signal recycled IFO (see shotrad for more info)
New version July 2016 by JH based on transfer function formalism
......@@ -311,7 +312,7 @@ def shotradSignalRecycled(f, ifo):
L = ifo.Infrastructure.Length # Length of arm cavities [m]
l = ifo.Optics.SRM.CavityLength # SRC Length [m]
T = ifo.Optics.ITM.Transmittance # ITM Transmittance [Power]
m = ifo.Materials.MirrorMass # Mirror mass [kg]
m = mirror_mass # Mirror mass [kg]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bsloss = ifo.Optics.BSLoss # BS Loss [Power]
......@@ -333,7 +334,7 @@ def shotradSignalRecycled(f, ifo):
R = 1 - T - ifo.Optics.Loss # ITM Reflectivity [Power]
P = ifo.gwinc.parm # use precomputed value
P = power.parm # use precomputed value
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
......@@ -460,7 +461,7 @@ def shotradSignalRecycled(f, ifo):
return coeff, Mifo, Msig, Mnoise
def shotradSignalRecycledBnC(f, ifo):
def shotradSignalRecycledBnC(f, ifo, mirror_mass, power):
"""Quantum noise model for signal recycled IFO
See shotrad for more info.
......@@ -492,7 +493,7 @@ def shotradSignalRecycledBnC(f, ifo):
L = ifo.Infrastructure.Length # Length of arm cavities [m]
l = ifo.Optics.SRM.CavityLength # SRC Length [m]
T = ifo.Optics.ITM.Transmittance # ITM Transmittance [Power]
m = ifo.Materials.MirrorMass # Mirror mass [kg]
m = mirror_mass # Mirror mass [kg]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bsloss = ifo.Optics.BSLoss # BS Loss [Power]
......@@ -514,7 +515,7 @@ def shotradSignalRecycledBnC(f, ifo):
gamma_ac = T*c/(4*L) # [KLMTV-PRD2001] Arm cavity half bandwidth [1/s]
epsilon = lambda_arm/(2*gamma_ac*L/c) # [BnC, after 5.2] Loss coefficent for arm cavity
I_0 = ifo.gwinc.pbs # [BnC, Table 1] Power at BS (Power*prfactor) [W]
I_0 = power.pbs # [BnC, Table 1] Power at BS (Power*prfactor) [W]
I_SQL = (m*L**2*gamma_ac**4)/(4*omega_0) # [BnC, 2.14] Power to reach free mass SQL
Kappa = 2*((I_0/I_SQL)*gamma_ac**4)/ \
(Omega**2*(gamma_ac**2+Omega**2)) # [BnC 2.13] Effective Radiation Pressure Coupling
......
gwinc/noise/seismic.py
View file @ 169b1afa
......@@ -6,20 +6,21 @@ import numpy as np
from scipy.interpolate import PchipInterpolator as interp1d
def seismic_suspension_fitered(sus, in_trans):
def seismic_suspension_fitered(sus, sustf, in_trans):
"""Seismic displacement noise for single suspended test mass.
:sus: gwinc suspension structure
:sus: suspension Struct
:sustf: sus transfer function Struct
:in_trans: input translational displacement spectrum
:returns: tuple of displacement noise power spectrum at :f:, and
horizontal and vertical components.
"""
hTable = sus.hTable
vTable = sus.vTable
hTable = sustf.hTable
vTable = sustf.vTable
theta = sus.VHCoupling.theta
theta = sustf.VHCoupling.theta
# horizontal noise total
nh = (abs(hTable)**2) * in_trans**2
......
gwinc/noise/substratethermal.py
View file @ 169b1afa
......@@ -12,52 +12,6 @@ from ..const import BESSEL_ZEROS as zeta
from ..const import J0M as j0m
def substrate_carrierdensity(f, materials, wBeam, exact=False):
"""Substrate thermal displacement noise spectrum from charge carrier density fluctuations
For semiconductor substrates.
:f: frequency array in Hz
:materials: gwinc optic materials structure
:wBeam: beam radius (at 1 / e^2 power)
:exact: whether to use adiabatic approximation or exact calculation (False)
:returns: displacement noise power spectrum at :f:, in meters
"""
H = materials.MassThickness
diffElec = materials.Substrate.ElectronDiffusion
diffHole = materials.Substrate.HoleDiffusion
mElec = materials.Substrate.ElectronEffMass
mHole = materials.Substrate.HoleEffMass
cdDens = materials.Substrate.CarrierDensity
gammaElec = materials.Substrate.ElectronIndexGamma
gammaHole = materials.Substrate.HoleIndexGamma
r0 = wBeam/np.sqrt(2)
omega = 2*pi*f
if exact:
def integrand(k, om, D):
return D * k**3 * exp(-k**2 * wBeam**2/4) / (D**2 * k**4 + om**2)
integralElec = np.array([scipy.integrate.quad(lambda k: integrand(k, om, diffElec), 0, inf)[0] for om in omega])
integralHole = np.array([scipy.integrate.quad(lambda k: integrand(k, om, diffHole), 0, inf)[0] for om in omega])
# From P1400084 Heinert et al. Eq. 15
#psdCD = @(gamma,m,int) 2*(3/pi^7)^(1/3)*kBT*H*gamma^2*m/hbar^2*cdDens^(1/3)*int; %units are meters
# FIXME: why the unused argument here?
def psdCD(gamma, m, int_):
return 2/pi * H * gamma**2 * cdDens * int_
psdElec = psdCD(gammaElec, mElec, integralElec)
psdHole = psdCD(gammaHole, mHole, integralHole)
else:
psdElec = 4*H*gammaElec**2*cdDens*diffElec/(pi*r0**4*omega**2)
psdHole = 4*H*gammaHole**2*cdDens*diffHole/(pi*r0**4*omega**2)
return psdElec + psdHole
def substrate_thermorefractive(f, materials, wBeam, exact=False):
"""Substrate thermal displacement noise spectrum from thermorefractive fluctuations
......@@ -87,7 +41,7 @@ def substrate_thermorefractive(f, materials, wBeam, exact=False):
# From P1400084 Heinert et al. Eq. 15
#psdCD = @(gamma,m,int) 2*(3/pi^7)^(1/3)*kBT*H*gamma^2*m/hbar^2*cdDens^(1/3)*int; %units are meters
psdTR = lambda int_: 2/pi * H * beta**2 * kBT * Temp / (rho*C) * int_;
psdTR = lambda int_: 2/pi * H * beta**2 * kBT * Temp / (rho*C) * int_
psd = psdTR(inte)
psd = 2/pi * H * beta**2 * kBT * Temp / (rho*C) * inte
......
gwinc/noise/suspensionthermal.py
View file @ 169b1afa
......@@ -9,11 +9,12 @@ from ..const import kB
from ..suspension import getJointParams
def suspension_thermal(f, sus):
def suspension_thermal(f, sus, sustf):
"""Suspension thermal noise.
:f: frequency array in Hz
:sus: gwinc suspension structure
:sustf: sus transfer function Struct
:returns: displacement noise power spectrum at :f:
......@@ -25,23 +26,14 @@ def suspension_thermal(f, sus):
pre-calculated and populated into the relevant `sus` fields.
"""
# suspension TFs
hForce = sustf.hForce
vForce = sustf.vForce
# and vertical to beamline coupling angle
theta = sus.VHCoupling.theta
theta = sustf.VHCoupling.theta
noise = np.zeros_like(f)
##########################################################
# Suspension TFs
##########################################################
hForce = sus.hForce
vForce = sus.vForce
##########################################################
# Thermal Noise Calculation
##########################################################
# Here the suspension stage list is reversed so that the last stage
# in the list is the test mass
stages = sus.Stage[::-1]
......
gwinc/plot.py
View file @ 169b1afa
from numpy import sqrt
from collections.abc import Mapping
def plot_noise(
freq,
traces,
def plot_budget(
budget,
ax=None,
**kwargs
):
"""Plot a GWINC noise budget from calculated noises
"""Plot a GWINC BudgetTrace noise budget from calculated noises
If an axis handle is provided it will be used for the plot.
......@@ -22,36 +17,29 @@ def plot_noise(
else:
fig = ax.figure
ylim = kwargs.get('ylim')
for name, trace in traces.items():
if isinstance(trace, Mapping):
trace = trace['Total']
total = budget.asd
ylim = [min(total)/10, max(total)]
style = dict(
color='#000000',
alpha=0.6,
lw=4,
)
style.update(getattr(budget, 'style', {}))
if 'label' in style:
style['label'] = 'Total ' + style['label']
else:
style['label'] = 'Total'
ax.loglog(budget.freq, total, **style)
try:
data = trace[0]
style = dict(**trace[1])
except TypeError:
data = trace
style = {}
# assuming all data is PSD
data = sqrt(data)
if name == 'Total' and not style:
style = dict(
color='#000000',
alpha=0.6,
lw=4,
)
if ylim is None:
ylim = [min(data)/10, max(data)]
for name, trace in budget.items():
style = trace.style
if 'label' not in style:
style['label'] = name
style['label'] = budget.name
if 'linewidth' in style:
style['lw'] = style['linewidth']
elif 'lw' not in style:
style['lw'] = 3
ax.loglog(freq, data, **style)
ax.loglog(budget.freq, trace.asd, **style)
ax.grid(
True,
......@@ -67,14 +55,19 @@ def plot_noise(
)
ax.autoscale(enable=True, axis='y', tight=True)
if ylim:
ax.set_ylim(ylim)
ax.set_xlim(freq[0], freq[-1])
ax.set_ylim(kwargs.get('ylim', ylim))
ax.set_xlim(budget.freq[0], budget.freq[-1])
ax.set_xlabel('Frequency [Hz]')
if 'ylabel' in kwargs:
ax.set_ylabel(kwargs['ylabel'])
if 'title' in kwargs:
ax.set_title(kwargs['title'])
title = kwargs['title']
if 'subtitle' in kwargs:
title += '\n' + kwargs['subtitle']
ax.set_title(title)
return fig
# FIXME: deprecate
plot_noise = plot_budget
gwinc/struct.py
View file @ 169b1afa
......@@ -173,7 +173,9 @@ class Struct(object):
"""
d = {}
for k,v in self.__dict__.items():
for k, v in self.__dict__.items():
if k[0] == '_':
continue
if isinstance(v, Struct):
d[k] = v.to_dict(array=array)
else:
......@@ -222,6 +224,8 @@ class Struct(object):
"""
for k,v in self.__dict__.items():
if k[0] == '_':
continue
if isinstance(v, type(self)):
for sk,sv in v.walk():
yield k+'.'+sk, sv
......@@ -233,6 +237,25 @@ class Struct(object):
except (AttributeError, TypeError):
yield k, v
def hash(self, keys=None):
"""Hash of Struct.walk() data (recursive)
"""
def filter_keys(kv):
k, v = kv
if keys:
return k in keys
else:
return True
def map_tuple(kv):
k, v = kv
if isinstance(v, list):
return k, tuple(v)
else:
return k, v
return hash(tuple(sorted(
map(map_tuple, filter(filter_keys, self.walk()))
)))
def diff(self, other):
"""Return tuple of differences between target IFO.
......