Compare revisions

Christopher Wipf · Jameson Rollins · Jameson Rollins · Jameson Rollins · Jameson Rollins · Jameson Rollins
--- a/HDF5_SCHEMATA
+++ b/HDF5_SCHEMATA
@@ -2,7 +2,9 @@


 This file describes a schemata for HDF5 storage of noise trace data
-and plot styling GWINC noise budgets.
+and plot styling GWINC noise budgets.  Python functions for writing
+budget data to and reading budget data from this format are included
+in the `gwinc.io` module.

 HDF5 is a heirarchical, structured data storage format [0].  Content
 is organized into a heirarchical folder-like structure, with two
@@ -72,8 +74,9 @@ and an int schema version number:
 /.attrs['SCHEMA_VERSION'] = 1
 ```

-Top-level attributes are generally used for specifying overall plot styling, but the
-following root attributes are typically defined:
+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:
 ```
 /.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
 /.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>

--- a/IFO.md
+++ b/IFO.md
@@ -2,8 +2,12 @@

 CI-generated plots and data for all IFOs included in pygwinc.

+Ranges in "Mpc" are for binary neutron stars (BNS) using the
+[inspiral_range](gwinc/inspiral-range>) package.
+
 ![IFO compare](https://gwinc.docs.ligo.org/pygwinc/ifo/all_compare.png)

+
 ## aLIGO

 * [ifo.yaml](gwinc/ifo/aLIGO/ifo.yaml)

--- a/README.md
+++ b/README.md
@@ -83,6 +83,23 @@ the path to the budget module/package:
 $ 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:
+```shell
+$ python3 -m gwinc -i Aplus
+GWINC interactive shell
+
+The 'ifo' Struct and 'traces' 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 [1]: 
+```
+
 See command help for more info:
 ```shell
 $ python3 -m gwinc -h

--- a/gwinc/__init__.py
+++ b/gwinc/__init__.py
@@ -10,6 +10,9 @@ from .struct import Struct
 from .plot import plot_noise


+logger = logging.getLogger('gwinc')
+
+
 def _load_module(name_or_path):
    """Load module from name or path.

@@ -61,15 +64,15 @@ def load_budget(name_or_path):
        bname, ext = os.path.splitext(os.path.basename(path))

        if ext in Struct.STRUCT_EXT:
-            logging.info("loading struct {}...".format(path))
+            logger.info("loading struct {}...".format(path))
            ifo = Struct.from_file(path)
            bname = 'aLIGO'
            modname = 'gwinc.ifo.aLIGO'
-            logging.info("loading budget {}...".format(modname))
+            logger.info("loading budget {}...".format(modname))

        else:
            modname = path
-            logging.info("loading module path {}...".format(modname))
+            logger.info("loading module path {}...".format(modname))

    else:
        if name_or_path not in IFOS:
@@ -79,7 +82,7 @@ def load_budget(name_or_path):
            ))
        bname = name_or_path
        modname = 'gwinc.ifo.'+name_or_path
-        logging.info("loading module {}...".format(modname))
+        logger.info("loading module {}...".format(modname))

    mod, modpath = _load_module(modname)

@@ -126,13 +129,13 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
    finesse = ifo.gwinc.finesse
    prfactor = ifo.gwinc.prfactor
    if ifo.Laser.Power * prfactor != pbs:
-        logging.warning("Thermal lensing limits input power to {} W".format(pbs/prfactor))
+        logger.warning("Thermal lensing limits input power to {} W".format(pbs/prfactor))

    # report astrophysical scores if so desired
    score = None
    if source:
-        logging.warning("Source score calculation currently not supported.  See `inspiral-range` package for similar functionality:")
-        logging.warning("https://git.ligo.org/gwinc/inspiral-range")
+        logger.warning("Source score calculation currently not supported.  See `inspiral-range` package for similar functionality:")
+        logger.warning("https://git.ligo.org/gwinc/inspiral-range")
        # score = int731(freq, noises['Total'], ifo, source)
        # score.Omega = intStoch(freq, noises['Total'], 0, ifo, source)

@@ -140,31 +143,31 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
    # output graphics

    if plot:
-        logging.info('Laser Power:            %7.2f Watt' % ifo.Laser.Power)
-        logging.info('SRM Detuning:           %7.2f degree' % (ifo.Optics.SRM.Tunephase*180/np.pi))
-        logging.info('SRM transmission:       %9.4f' % ifo.Optics.SRM.Transmittance)
-        logging.info('ITM transmission:       %9.4f' % ifo.Optics.ITM.Transmittance)
-        logging.info('PRM transmission:       %9.4f' % ifo.Optics.PRM.Transmittance)
-        logging.info('Finesse:                %7.2f' % finesse)
-        logging.info('Power Recycling Gain:   %7.2f' % prfactor)
-        logging.info('Arm Power:              %7.2f kW' % (parm/1000))
-        logging.info('Power on BS:            %7.2f W' % pbs)
+        logger.info('Laser Power:            %7.2f Watt' % ifo.Laser.Power)
+        logger.info('SRM Detuning:           %7.2f degree' % (ifo.Optics.SRM.Tunephase*180/np.pi))
+        logger.info('SRM transmission:       %9.4f' % ifo.Optics.SRM.Transmittance)
+        logger.info('ITM transmission:       %9.4f' % ifo.Optics.ITM.Transmittance)
+        logger.info('PRM transmission:       %9.4f' % ifo.Optics.PRM.Transmittance)
+        logger.info('Finesse:                %7.2f' % finesse)
+        logger.info('Power Recycling Gain:   %7.2f' % prfactor)
+        logger.info('Arm Power:              %7.2f kW' % (parm/1000))
+        logger.info('Power on BS:            %7.2f W' % pbs)

        # coating and substrate thermal load on the ITM
        PowAbsITM = (pbs/2) * \
                    np.hstack([(finesse*2/np.pi) * ifo.Optics.ITM.CoatingAbsorption,
                               (2 * ifo.Materials.MassThickness) * ifo.Optics.ITM.SubstrateAbsorption])

-        logging.info('Thermal load on ITM:    %8.3f W' % sum(PowAbsITM))
-        logging.info('Thermal load on BS:     %8.3f W' %
+        logger.info('Thermal load on ITM:    %8.3f W' % sum(PowAbsITM))
+        logger.info('Thermal load on BS:     %8.3f W' %
                     (ifo.Materials.MassThickness*ifo.Optics.SubstrateAbsorption*pbs))
        if (ifo.Laser.Power*prfactor != pbs):
-            logging.info('Lensing limited input power: %7.2f W' % (pbs/prfactor))
+            logger.info('Lensing limited input power: %7.2f W' % (pbs/prfactor))

        if score:
-            logging.info('BNS Inspiral Range:     ' + str(score.effr0ns) + ' Mpc/ z = ' + str(score.zHorizonNS))
-            logging.info('BBH Inspiral Range:     ' + str(score.effr0bh) + ' Mpc/ z = ' + str(score.zHorizonBH))
-            logging.info('Stochastic Omega: %4.1g Universes' % score.Omega)
+            logger.info('BNS Inspiral Range:     ' + str(score.effr0ns) + ' Mpc/ z = ' + str(score.zHorizonNS))
+            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)


--- a/gwinc/__main__.py
+++ b/gwinc/__main__.py
 from __future__ import print_function
 import os
 import signal
+import logging
 import argparse
 import numpy as np

-import logging
-logging.basicConfig(
-    format='%(message)s',
-    level=os.getenv('LOG_LEVEL', logging.WARNING),
-)
+from . import IFOS, load_budget, plot_noise, logger

-from . import IFOS, load_budget, plot_noise
+logger.setLevel(os.getenv('LOG_LEVEL', 'WARNING').upper())
+formatter = logging.Formatter('%(name)s: %(message)s')
+handler = logging.StreamHandler()
+handler.setFormatter(formatter)
+logger.addHandler(handler)

 ##################################################

-description = """Plot GWINC noise budget for specified IFO.
+description = """GWINC noise budget tool

-Available included IFOs: {}
+IFOs can be specified by name of included canonical budget (see
+below), or by path to a budget module (.py), description file
+(.yaml/.mat/.m), or HDF5 data file (.hdf5/.h5).  Available included
+IFOs are:

-""".format(', '.join(["'{}'".format(ifo) for ifo in IFOS]))
+  {}
+
+""".format(', '.join(["{}".format(ifo) for ifo in IFOS]))
 # for ifo in available_ifos():
 #     description += "  '{}'\n".format(ifo)
 description += """

-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' 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.
-
-Individual IFO parameters can be overriden with the --ifo option:
+By default the noise budget of the specified IFO will be loaded and
+plotted with an interactive plotter.  Individual IFO parameters can be
+overriden with the --ifo option:

  gwinc --ifo Optics.SRM.Tunephase=3.14 ...

-If the inspiral_range package is installed, various figures of merit
+If the --save option is specified the plot will be saved directly to a
+file (without display) (various file formats are supported, indicated
+by file extension).  If the requested extension is 'hdf5' or 'h5' then
+the noise traces and IFO parameters will be saved to an HDF5 file.
+
+If the inspiral_range package is available, various figures of merit
 can be calculated for the resultant spectrum with the --fom option,
 e.g.:

  gwinc --fom horizon ...
  gwinc --fom range:m1=20,m2=20 ...

-See documentation for inspiral_range package for details.
-
+See the inspiral_range package documentation for details.
 """

 IFO = 'aLIGO'
 FREQ = '5:3000:6000'
+FOM = 'range:m1=1.4,m2=1.4'
 DATA_SAVE_FORMATS = ['.hdf5', '.h5']

 parser = argparse.ArgumentParser(
@@ -57,10 +60,10 @@ parser = argparse.ArgumentParser(
    description=description,
    formatter_class=argparse.RawDescriptionHelpFormatter)
 parser.add_argument(
-    '--freq', '-f', metavar='FSPEC',
-    help="frequency array specification in Hz, either as 'flo:fhi' or 'flo:npoints:fhi' [{}]".format(FREQ))
+    '--freq', '-f', metavar='FLO:[NPOINTS:]FHI',
+    help="frequency array specification in Hz [{}]".format(FREQ))
 parser.add_argument(
-    '--ifo', '-o',
+    '--ifo', '-o', metavar='PARAM=VAL',
    #nargs='+', action='extend',
    action='append',
    help="override budget IFO parameter (may be specified multiple times)")
@@ -68,14 +71,14 @@ parser.add_argument(
    '--title', '-t',
    help="plot title")
 parser.add_argument(
-    '--fom',
-    help="calculate inspiral range for resultant spectrum ('func[:param=val,param=val]')")
+    '--fom', metavar='FUNC[:PARAM=VAL,...]',
+    help="use inspiral_range.FUNC to calculate range figure-of-merit on resultant spectrum")
 group = parser.add_mutually_exclusive_group()
 group.add_argument(
    '--interactive', '-i', action='store_true',
-    help="interactive plot with interactive shell")
+    help="launch interactive shell after budget processing")
 group.add_argument(
-    '--save', '-s', action='append',
+    '--save', '-s', metavar='PATH', action='append',
    help="save plot (.png/.pdf/.svg) or budget traces (.hdf5/.h5) to file (may be specified multiple times)")
 group.add_argument(
    '--yaml', '-y', action='store_true',
@@ -85,13 +88,13 @@ group.add_argument(
    help="print IFO as text table to stdout and exit (budget not calculated)")
 group.add_argument(
    '--diff', '-d', metavar='IFO',
-    help="show differences table between another IFO description and exit (budget not calculated)")
-group.add_argument(
+    help="show difference table between IFO and another IFO description (name or path) and exit (budget not calculated)")
+parser.add_argument(
    '--no-plot', '-np', action='store_false', dest='plot',
-    help="supress plotting")
+    help="suppress plotting")
 parser.add_argument(
    'IFO',
-    help="IFO name, or path to budget module (.py), description file (.yaml/.mat/.m), or HDF5 data file (.hdf5/.h5)")
+    help="IFO name or path")


 def freq_from_spec(spec):
@@ -114,13 +117,17 @@ def main():
    # initial arg processing

    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")
+        if args.ifo:
+            parser.exit(2, "Can not override ifo parameters when loading traces from file.\n")
        from .io import load_hdf5
        Budget = None
        freq, traces, attrs = load_hdf5(args.IFO)
-        ifo = getattr(attrs, 'IFO', None)
+        ifo = attrs.get('ifo')
+        # FIXME: deprecate 'IFO'
+        ifo = attrs.get('IFO', ifo)
        plot_style = attrs
-        if args.freq:
-            logging.warning("ignoring frequency specification for frequencies defined in HDF5...")

    else:
        Budget = load_budget(args.IFO)
@@ -135,14 +142,6 @@ def main():
        plot_style = getattr(Budget, 'plot_style', {})
        traces = None

-    out_data_files = set()
-    out_plot_files = set()
-    if args.save:
-        for path in args.save:
-            if os.path.splitext(path)[1] in DATA_SAVE_FORMATS:
-                out_data_files.add(path)
-        out_plot_files = set(args.save) - out_data_files
-
    if args.ifo:
        for paramval in args.ifo:
            param, val = paramval.split('=', 1)
@@ -150,22 +149,24 @@ def main():

    if args.yaml:
        if not ifo:
-            parser.exit(2, "no IFO structure available.")
+            parser.exit(2, "No 'ifo' structure available.\n")
        print(ifo.to_yaml(), end='')
        return
    if args.text:
        if not ifo:
-            parser.exit(2, "no IFO structure available.")
+            parser.exit(2, "No 'ifo' structure available.\n")
        print(ifo.to_txt(), end='')
        return
    if args.diff:
        if not ifo:
-            parser.exit(2, "no IFO structure available.")
-        fmt = '{:30} {:>20} {:>20}'
+            parser.exit(2, "No 'ifo' structure available.\n")
        Budget = load_budget(args.diff)
        diffs = ifo.diff(Budget.ifo)
        if diffs:
+            w = max([len(d[0]) for d in diffs])
+            fmt = '{{:{}}} {{:>20}} {{:>20}}'.format(w)
            print(fmt.format('', args.IFO, args.diff))
+            print(fmt.format('', '-----', '-----'))
            for p in diffs:
                k = str(p[0])
                v = repr(p[1])
@@ -173,15 +174,17 @@ def main():
                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)
+    out_data_files = set()
+    out_plot_files = set()
+    if args.save:
+        args.plot = False
+        for path in args.save:
+            if os.path.splitext(path)[1] in DATA_SAVE_FORMATS:
+                out_data_files.add(path)
+        out_plot_files = set(args.save) - out_data_files

-    if args.plot:
-        if args.save:
+    if args.plot or out_plot_files:
+        if out_plot_files:
            # FIXME: this silliness seems to be the only way to have
            # matplotlib usable on systems without a display.  There must
            # be a better way.  'AGG' is a backend that works without
@@ -194,11 +197,22 @@ def main():
        try:
            from matplotlib import pyplot as plt
        except RuntimeError:
-            logging.warning("no display, plotting disabled.")
+            logger.warning("no display, plotting disabled.")
            args.plot = False

-    if args.fom:
+    try:
        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
+    except ModuleNotFoundError:
+        if args.fom:
+            raise
+        logger.warning("inspiral_range package not available, figure of merit will not be calculated...")
+    if args.fom:
        try:
            range_func, fargs = args.fom.split(':')
        except ValueError:
@@ -208,32 +222,43 @@ def main():
        for param in fargs.split(','):
            if not param:
                continue
-            p,v = param.split('=')
+            p, v = param.split('=')
            if not v:
                raise ValueError('missing parameter value "{}"'.format(p))
-            if p != 'approximant':
+            try:
                v = float(v)
+            except ValueError:
+                pass
            range_params[p] = v

    ##########
    # main calculations

    if not traces:
-        logging.info("calculating budget...")
+        logger.info("calculating budget...")
        traces = Budget(freq=freq, ifo=ifo).run()

-    # logging.info('recycling factor: {: >0.3f}'.format(ifo.gwinc.prfactor))
-    # logging.info('BS power:         {: >0.3f} W'.format(ifo.gwinc.pbs))
-    # logging.info('arm finesse:      {: >0.3f}'.format(ifo.gwinc.finesse))
-    # logging.info('arm power:        {: >0.3f} kW'.format(ifo.gwinc.parm/1000))
+    # 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))
+
+    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)
+    else:
+        plot_style['title'] = "GWINC Noise Budget: {}".format(args.IFO)

    if args.fom:
-        logging.info("calculating inspiral {}...".format(range_func))
+        logger.info("calculating inspiral {}...".format(range_func))
        H = inspiral_range.CBCWaveform(freq, **range_params)
-        logging.debug("params: {}".format(H.params))
+        logger.debug("params: {}".format(H.params))
        fom = eval('inspiral_range.{}'.format(range_func))(freq, traces['Total'][0], H=H)
-        logging.info("{}({}) = {:.2f} Mpc".format(range_func, fargs, fom))
-        fom_title = '{func} {m1}/{m2} Msol: {fom:.2f} Mpc'.format(
+        logger.info("{}({}) = {:.2f} Mpc".format(range_func, fargs, fom))
+        fom_title = 'inspiral {func} {m1}/{m2} $\mathrm{{M}}_\odot$: {fom:.0f} Mpc'.format(
            func=range_func,
            m1=H.params['m1'],
            m2=H.params['m2'],
@@ -246,16 +271,26 @@ def main():

    # interactive shell plotting
    if args.interactive:
-        from IPython.terminal.embed import InteractiveShellEmbed
-        ipshell = InteractiveShellEmbed(
-            banner1='''
-GWINC interactive plotter
+        banner = """GWINC interactive shell

-You may interact with plot using "plt." methods, e.g.:
+The 'ifo' Struct and 'traces' data 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:
+
+In [.]: plot_noise(freq, traces, **plot_style)
+"""
+        banner += """
+You may interact with the plot using the 'plt' functions, e.g.:

->>> plt.title("foo")
->>> plt.savefig("foo.pdf")
-''',
+In [.]: plt.title("foo")
+In [.]: plt.savefig("foo.pdf")
+"""
+        from IPython.terminal.embed import InteractiveShellEmbed
+        ipshell = InteractiveShellEmbed(
+            banner1=banner,
            user_ns={
                'freq': freq,
                'traces': traces,
@@ -265,8 +300,9 @@ You may interact with plot using "plt." methods, e.g.:
            },
        )
        ipshell.enable_pylab()
-        ipshell.ex("fig = plot_noise(freq, traces, **plot_style)")
-        ipshell.ex("plt.title('{}')".format(plot_style['title']))
+        if args.plot:
+            ipshell.ex("fig = plot_noise(freq, traces, **plot_style)")
+            ipshell.ex("plt.title('{}')".format(plot_style['title']))
        ipshell()

    ##########
@@ -276,19 +312,19 @@ You may interact with plot using "plt." methods, e.g.:
    if out_data_files:
        from .io import save_hdf5
        attrs = dict(plot_style)
-        attrs['IFO'] = ifo.to_yaml()
        for path in out_data_files:
-            logging.info("saving budget traces {}...".format(path))
+            logger.info("saving budget traces: {}".format(path))
            save_hdf5(
                path=path,
                freq=freq,
                traces=traces,
+                ifo=ifo,
                **attrs,
            )

    # standard plotting
    if args.plot or out_plot_files:
-        logging.info("plotting noises...")
+        logger.debug("plotting noises...")
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        plot_noise(
@@ -300,6 +336,7 @@ You may interact with plot using "plt." methods, e.g.:
        fig.tight_layout()
        if out_plot_files:
            for path in out_plot_files:
+                logger.info("saving budget plot: {}".format(path))
                fig.savefig(path)
        else:
            plt.show()

--- a/gwinc/gwinc_matlab.py
+++ b/gwinc/gwinc_matlab.py
 import os
 import copy
 import tempfile
-import logging

 import scipy.io

+from . import logger
 from . import const
 from . import suspension
 from .struct import Struct
@@ -33,10 +33,10 @@ class Matlab:
        if not os.path.exists(os.path.join(gwincpath, 'gwinc.m')):
            raise IOError("Invalid MATLAB GWINC path: '{}'".format(gwincpath))

-        logging.info("starting MATLAB engine...")
+        logger.info("starting MATLAB engine...")
        MATLAB_ENGINE = matlab.engine.start_matlab()

-        logging.info("gwinc path: {}".format(gwincpath))
+        logger.info("gwinc path: {}".format(gwincpath))
        MATLAB_ENGINE.addpath(gwincpath)



--- a/gwinc/ifo/Aplus/ifo.yaml
+++ b/gwinc/ifo/Aplus/ifo.yaml
@@ -193,19 +193,19 @@ Materials:
  ## Dielectric coating material parameters----------------------------------
  Coating:
    ## high index material: tantala
-    Yhighn: 124e9                 # LMA (Granata at LVC) 2017 (was 140)
-    Sigmahighn: 0.28              # LMA (Granata at LVC) 2017 (was 0.23)
+    Yhighn: 120e9                 # Ta2O5-TiO2 from 2020 LMA https://iopscience.iop.org/article/10.1088/1361-6382/ab77e9
+    Sigmahighn: 0.29              # 2020 LMA
    CVhighn: 2.1e6                # Crooks et al, Fejer et al
    Alphahighn: 3.6e-6            # 3.6e-6 Fejer et al, 5e-6 from Braginsky
    Betahighn: 1.4e-5             # dn/dT, value Gretarrson (G070161)
    ThermalDiffusivityhighn: 33   # Fejer et al
-    Indexhighn: 2.06539
+    Indexhighn: 2.09              # 2020 LMA
    Phihighn: 9.0e-5              # tantala mechanical loss
    Phihighn_slope: 0.1

    ## low index material: silica
-    Ylown: 72e9
-    Sigmalown: 0.17
+    Ylown: 70e9                   # 2020 LMA
+    Sigmalown: 0.19               # 2020 LMA
    CVlown: 1.6412e6              # Crooks et al, Fejer et al
    Alphalown: 5.1e-7             # Fejer et al
    Betalown: 8e-6                # dn/dT,  (ref. 14)

--- a/gwinc/ifo/__main__.py
+++ b/gwinc/ifo/__main__.py
@@ -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

@@ -38,13 +36,20 @@ def main():

    for name, budget in budgets.items():
        data = budget.calc()
-        BNS_range = inspiral_range.range(freq, data)
-        label = '{name:<{pad}} {bns:>6.0f} Mpc'.format(
+        try:
+            import inspiral_range
+            label_range = ' {:>6.0f} Mpc'.format(
+                inspiral_range.range(freq, data),
+            )
+        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, np.sqrt(data), label=label, lw=2)

    ax.grid(
        True,
@@ -65,8 +70,7 @@ def main():

    ax.set_xlabel('Frequency [Hz]')
    ax.set_ylabel(PLOT_STYLE['ylabel'])
-    ax.set_title("""PyGWINC reference IFO strain comparison
-(with BNS range)""")
+    ax.set_title("PyGWINC reference IFO strain comparison")

    if args.save:
        fig.savefig(args.save)

--- a/gwinc/ifo/aLIGO/ifo.yaml
+++ b/gwinc/ifo/aLIGO/ifo.yaml
@@ -194,25 +194,25 @@ Materials:
  ## Dielectric coating material parameters----------------------------------
  Coating:
    ## high index material: tantala
-    Yhighn: 124e9                 # LMA (Granata at LVC) 2017 (was 140)
-    Sigmahighn: 0.28              # LMA (Granata at LVC) 2017 (was 0.23)
+    Yhighn: 120e9                 # Ta2O5-TiO2 from 2020 LMA https://iopscience.iop.org/article/10.1088/1361-6382/ab77e9
+    Sigmahighn: 0.29              # 2020 LMA
    CVhighn: 2.1e6                # Crooks et al, Fejer et al
    Alphahighn: 3.6e-6            # 3.6e-6 Fejer et al, 5e-6 from Braginsky
    Betahighn: 1.4e-5             # dn/dT, value Gretarrson (G070161)
    ThermalDiffusivityhighn: 33   # Fejer et al
-    Indexhighn: 2.06539
-    Phihighn: 3.6e-4              # loss angle at 100Hz (Gras 2018)
+    Indexhighn: 2.09              # 2020 LMA
+    Phihighn: 3.89e-4             # loss angle at 100Hz (Gras 2020)
    Phihighn_slope: 0.1

    ## low index material: silica
-    Ylown: 72e9
-    Sigmalown: 0.17
+    Ylown: 70e9                   # 2020 LMA
+    Sigmalown: 0.19               # 2020 LMA
    CVlown: 1.6412e6              # Crooks et al, Fejer et al
    Alphalown: 5.1e-7             # Fejer et al
    Betalown: 8e-6                # dn/dT,  (ref. 14)
    ThermalDiffusivitylown: 1.38  # Fejer et al
    Indexlown: 1.45
-    Philown: 4e-5                 # loss angle at 100Hz
+    Philown: 2.3e-5               # loss angle at 100Hz 2020 LMA
    Philown_slope: 0              # G1600641 and arXiv:1712.05701 suggest
                                  # slopes between 0 and 0.3, depending on
                                  # deposition method. Slawek's analysis in
@@ -268,6 +268,6 @@ Optics:
    CavityLength: 55              # m, ITM to SRM distance
    Tunephase: 0.0                # SEC tuning

-  Curvature:                      # ROC
-    ITM: 1970
-    ETM: 2192
+  Curvature:                      # RoCs per E080511, E080512
+    ITM: 1934
+    ETM: 2245
--- a/gwinc/ifo/noises.py
+++ b/gwinc/ifo/noises.py
 import numpy as np
 from numpy import pi, sin, exp, sqrt

+from .. import logger
 from .. import const
 from ..struct import Struct
 from .. import nb
@@ -30,7 +31,7 @@ def arm_cavity(ifo):
    if (g1 * g2 * (1 - g1 * g2)) <= 0:
        raise Exception('Unstable arm cavity g-factors.  Change ifo.Optics.Curvature')
    elif gcav < 1e-3:
-        logging.warning('Nearly unstable arm cavity g-factors.  Reconsider ifo.Optics.Curvature')
+        logger.warning('Nearly unstable arm cavity g-factors.  Reconsider ifo.Optics.Curvature')

    ws = sqrt(L * ifo.Laser.Wavelength / pi)
    w1 = ws * sqrt(abs(g2) / gcav)
@@ -50,7 +51,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)
@@ -64,7 +64,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
@@ -78,15 +78,15 @@ 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:
-        logging.warning('P_BS exceeds BS Thermal limit!')
+        logger.warning('P_BS exceeds BS Thermal limit!')

    return pbs, parm, finesse, prfactor, Tpr

@@ -140,7 +140,6 @@ 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
@@ -178,7 +177,7 @@ class StandardQuantumLimit(nb.Noise):
    style = dict(
        label="Standard Quantum Limit",
        color="#000000",
-        linestyle=":", # Dotted.
+        linestyle=":",
    )

    def calc(self):
@@ -339,7 +338,7 @@ class ITMThermoRefractive(nb.Noise):

    def calc(self):
        finesse = ifo_power(self.ifo)[2]
-        gPhase =  finesse * 2/np.pi
+        gPhase = finesse * 2/np.pi
        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
        n = noise.substratethermal.substrate_thermorefractive(
            self.freq, self.ifo.Materials, wBeam_ITM)

--- a/gwinc/io.py
+++ b/gwinc/io.py
 import h5py
+import yaml
 import datetime

+from . import logger
+from . import Struct
+

 SCHEMA = 'GWINC noise budget'
 SCHEMA_VERSION = 1
@@ -19,7 +23,13 @@ def _write_trace_recursive(grp, traces):


 def save_hdf5(path, freq, traces, **kwargs):
-    """Save GWINC traces dict to HDF5 file.
+    """Save GWINC budget data to an HDF5 file.
+
+    The `freq` argument should be the frequency array, and `traces`
+    should be the traces (recursive) dictionary.  Keyword arguments
+    are stored in the HDF5 top level 'attrs' key-value store.  If an
+    'ifo' keyword arg is supplied, it is assumed to be a Struct and
+    will be serialized to YAML for storage.

    See HDF5_SCHEMA.

@@ -30,7 +40,10 @@ def save_hdf5(path, freq, traces, **kwargs):
        # FIXME: add budget code hash or something
        f.attrs['date'] = datetime.datetime.now().isoformat()
        for key, val in kwargs.items():
-            f.attrs[key] = val
+            if key == 'ifo':
+                f.attrs['ifo'] = val.to_yaml()
+            else:
+                f.attrs[key] = val
        f.create_dataset('Freq', data=freq)
        tgrp = f.create_group('traces')
        _write_trace_recursive(tgrp, traces)
@@ -47,14 +60,22 @@ def _read_trace_recursive(element):


 def load_hdf5(path):
-    """Load traces from HDF5 file.
+    """Load GWINC budget data from an HDF5 file.

-    Returns a recursive traces dictionary.  See HDF5_SCHEMA.
+    Returns (freq, traces, attrs).  An 'ifo' attr will be
+    de-serialized from YAML into a Struct object.
+
+    See HDF5_SCHEMA.

    """
    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.items())
+        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
--- a/gwinc/nb.py
+++ b/gwinc/nb.py
-import os
-import logging
 import operator
 import itertools
-import importlib
-import importlib.util
 import collections
 import numpy as np
 import scipy.interpolate
 from functools import reduce

+from . import logger
+

 def quadsum(data):
    """Calculate quadrature sum of list of data arrays.
@@ -220,7 +218,7 @@ class Budget(Noise):
    references = []
    """List of reference noise classes, or (ref, cal) tuples"""

-    def __init__(self, *args, noises=None, **kwargs):
+    def __init__(self, freq=None, noises=None, **kwargs):
        """Initialize Budget object.

        See BudgetItem for base initialization arguments.
@@ -230,15 +228,20 @@ class Budget(Noise):
        be used to filter the noises initialized in this budget.

        """
-        super().__init__(*args, **kwargs)
-        # store args and kwargs for later use
-        self.args = args
+        super().__init__(freq, **kwargs)
+        # store kwargs for later use
        self.kwargs = kwargs
-        # FIXME: special casing the IFO here, in case it's defined as
-        # a class attribute rather than passed at initialization.  we
-        # do this because we're not defining a standard way to extract
-        # IFO variables that get passed around in a reasonable way.
-        # how can we clarify this?
+        # record the frequency array as a kwarg if it's definied as a
+        # class attribute
+        if freq is not None:
+            self.kwargs['freq'] = freq
+        else:
+            self.kwargs['freq'] = getattr(self, 'freq', None)
+        # FIXME: special casing the ifo kwarg here, in case it's
+        # defined as a class attribute rather than passed at
+        # initialization.  we do this because we're not defining a
+        # standard way to extract IFO variables that get passed around
+        # in a reasonable way.  how can we clarify this?
        if 'ifo' not in kwargs and getattr(self, 'ifo', None):
            self.kwargs['ifo'] = getattr(self, 'ifo', None)
        # all noise objects keyed by name
@@ -276,13 +279,12 @@ class Budget(Noise):
            noise = nc
            cals = []
        noise_obj = noise(
-            *self.args,
            **self.kwargs
        )
        name = noise_obj.name
        if noise_filt and name not in noise_filt:
            return
-        logging.debug("init {}".format(noise_obj))
+        logger.debug("init {}".format(noise_obj))
        self._noise_objs[name] = noise_obj
        for cal in cals:
            self.__add_calibration(cal, [name])
@@ -290,12 +292,11 @@ class Budget(Noise):

    def __add_calibration(self, cal, noises):
        cal_obj = cal(
-            *self.args,
            **self.kwargs
        )
        name = cal_obj.name
        if name not in self._cal_objs:
-            logging.debug("init {}".format(cal_obj))
+            logger.debug("init {}".format(cal_obj))
            self._cal_objs[name] = cal_obj
        # register noises for this calibration
        for noise in noises:
@@ -357,7 +358,7 @@ class Budget(Noise):
        for name, item in itertools.chain(
                self._cal_objs.items(),
                self._noise_objs.items()):
-            logging.debug("load {}".format(item))
+            logger.debug("load {}".format(item))
            item.load()

    def update(self, **kwargs):
@@ -365,7 +366,7 @@ class Budget(Noise):
        for name, item in itertools.chain(
                self._cal_objs.items(),
                self._noise_objs.items()):
-            logging.debug("update {}".format(item))
+            logger.debug("update {}".format(item))
            item.update(**kwargs)

    def calc_noise(self, name):
@@ -433,12 +434,14 @@ class Budget(Noise):
        # calc all calibrations
        c = {}
        for name, cal in self._cal_objs.items():
+            logger.debug("calc {}".format(name))
            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
+            logger.debug("calc {}".format(name))
            d[name] = noise.calc_trace(
                calibrations=cals,
            )

--- a/gwinc/noise/quantum.py
+++ b/gwinc/noise/quantum.py
@@ -4,11 +4,12 @@
 from __future__ import division
 import numpy as np
 from numpy import pi, sqrt, arctan, sin, cos, exp, log10, conj
-import logging

+from .. import logger
 from .. import const
 from ..struct import Struct

+
 def sqzOptimalSqueezeAngle(Mifo, eta):
    vHD = np.array([[sin(eta), cos(eta)]])
    H = getProdTF(vHD, Mifo)[0]
@@ -53,10 +54,10 @@ def shotrad(f, ifo):
    Nfreq = len(f)
    if any(np.array([Mifo.shape[0], Mifo.shape[1], Mn.shape[0]]) != Nfield) or \
       any(np.array([Mifo.shape[2], Msig.shape[2], Mn.shape[2]]) != Nfreq):
-        logging.debug(Mifo.shape)
-        logging.debug(Msig.shape)
-        logging.debug(Mn.shape)
-        logging.debug(Nfield, Nfreq)
+        logger.debug(Mifo.shape)
+        logger.debug(Msig.shape)
+        logger.debug(Mn.shape)
+        logger.debug(Nfield, Nfreq)
        raise Exception('Inconsistent matrix sizes returned by %s' % str(fname))

    # deal with non-standard number of fields
@@ -93,7 +94,7 @@ def shotrad(f, ifo):
        raise Exception("Unknown Readout Type")

    if eta_orig is not None:
-        # logging.warn((
+        # logger.warn((
        #     'Quadrature.dc is redundant with '
        #     'Squeezer.Readout and is deprecated.'
        # ))
@@ -129,13 +130,13 @@ def shotrad(f, ifo):
        pass

    elif sqzType == 'Freq Independent':
-        logging.debug('You are injecting %g dB of frequency independent squeezing' % SQZ_DB)
+        logger.debug('You are injecting %g dB of frequency independent squeezing' % SQZ_DB)

    elif sqzType == 'Optimal':
        # compute optimal squeezing angle
        alpha = sqzOptimalSqueezeAngle(Mifo, eta)

-        logging.debug('You are injecting %g dB of squeezing with optimal frequency dependent squeezing angle' % SQZ_DB)
+        logger.debug('You are injecting %g dB of squeezing with optimal frequency dependent squeezing angle' % SQZ_DB)

    elif sqzType == 'OptimalOptimal':
        # compute optimal squeezing angle, assuming optimal readout phase
@@ -144,10 +145,10 @@ def shotrad(f, ifo):
        MsigPD = Msig * sqrt(1 - lambda_PD)
        alpha = sqzOptimalSqueezeAngle(Mifo, [], [R, lambda_in], MsigPD, MnPD)

-        logging.debug('You are injecting %g dB of squeezing with optimal FD squeezing angle, for optimal readout phase' % SQZ_DB)
+        logger.debug('You are injecting %g dB of squeezing with optimal FD squeezing angle, for optimal readout phase' % SQZ_DB)

    elif sqzType == 'Freq Dependent':
-        logging.debug('You are injecting %g dB of squeezing with frequency dependent squeezing angle' % SQZ_DB)
+        logger.debug('You are injecting %g dB of squeezing with frequency dependent squeezing angle' % SQZ_DB)

    else:
        raise Exception('ifo.Squeezer.Type must be None, Freq Independent, Optimal, or Frequency Dependent, not "%s"' % sqzType)
@@ -176,7 +177,7 @@ def shotrad(f, ifo):
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    # Inject squeezed field into the IFO via some filter cavities
    if sqzType == 'Freq Dependent' and 'FilterCavity' in ifo.Squeezer:
-        logging.debug('  Applying %d input filter cavities' % np.atleast_1d(ifo.Squeezer.FilterCavity).size)
+        logger.debug('  Applying %d input filter cavities' % np.atleast_1d(ifo.Squeezer.FilterCavity).size)
        Mr, Msqz = sqzFilterCavityChain(f, np.atleast_1d(ifo.Squeezer.FilterCavity), Msqz)

    #####################################################
@@ -198,14 +199,14 @@ def shotrad(f, ifo):
            pass

        elif ifo.OutputFilter.Type == 'Chain':
-            logging.debug('  Applying %d output filter cavities' % np.atleast_1d(ifo.OutputFilter.FilterCavity).size)
+            logger.debug('  Applying %d output filter cavities' % np.atleast_1d(ifo.OutputFilter.FilterCavity).size)

            Mr, Mnoise = sqzFilterCavityChain(f, np.atleast_1d(ifo.OutputFilter.FilterCavity), Mnoise)
            Msig = getProdTF(Mr, Msig)
            #  Mnoise = getProdTF(Mn, Mnoise);

        elif ifo.OutputFilter.Type == 'Optimal':
-            logging.debug('  Optimal output filtering!')
+            logger.debug('  Optimal output filtering!')

            # compute optimal angle, including upcoming PD losses
            MnPD = sqzInjectionLoss(Mnoise, lambda_PD)

--- a/gwinc/plot.py
+++ b/gwinc/plot.py
@@ -27,14 +27,16 @@ def plot_noise(
    for name, trace in traces.items():
        if isinstance(trace, Mapping):
            trace = trace['Total']
+
        try:
-            data, style = trace
-        except:
+            data = trace[0]
+            style = dict(**trace[1])
+        except TypeError:
            data = trace
            style = {}
        # assuming all data is PSD
        data = sqrt(data)
-        if name == 'Total':
+        if name == 'Total' and not style:
            style = dict(
                color='#000000',
                alpha=0.6,
@@ -48,6 +50,7 @@ def plot_noise(
            style['lw'] = style['linewidth']
        elif 'lw' not in style:
            style['lw'] = 3
+
        ax.loglog(freq, data, **style)

    ax.grid(

--- a/gwinc/struct.py
+++ b/gwinc/struct.py
@@ -174,7 +174,7 @@ class Struct(object):
        """
        d = {}
        for k,v in self.__dict__.items():
-            if isinstance(v, type(self)):
+            if isinstance(v, Struct):
                d[k] = v.to_dict(array=array)
            else:
                if isinstance(v, list):
@@ -188,7 +188,9 @@ class Struct(object):
                    except AttributeError:
                        if array:
                            v = np.array(v)
-                elif isinstance(v, int):
+                # FIXME: there must be a better way to just match all
+                # numeric scalar types
+                elif isinstance(v, (int, float, np.int, np.float)):
                    v = float(v)
                d[k] = v
        return d
@@ -206,12 +208,12 @@ class Struct(object):
        else:
            return y

-    # def __repr__(self):
-    #     return self.to_yaml().strip('\n')
-
    def __str__(self):
        return '<GWINC Struct: {}>'.format(list(self.__dict__.keys()))

+    def __repr__(self):
+        return self.__str__()
+
    def __iter__(self):
        return iter(self.__dict__)


--- a/gwinc/suspension.py
+++ b/gwinc/suspension.py
 from __future__ import division
 from numpy import pi, sqrt, sin, cos, tan, real, imag
 import numpy as np
-import logging

+from . import logger
 from .const import g


@@ -122,7 +122,7 @@ def getJointParams(sus, n):

    WireMaterialUpper = sus[wireMatUpper]
    WireMaterialLower = sus[wireMatLower]
-    logging.debug('stage {} wires: {}, {}'.format(n, wireMatUpper, wireMatLower))
+    logger.debug('stage {} wires: {}, {}'.format(n, wireMatUpper, wireMatLower))

    # support blade (upper joint only)
    if 'BladeMaterial' in stage:
@@ -132,7 +132,7 @@ def getJointParams(sus, n):

    BladeMaterial = sus[bladeMat]

-    logging.debug('stage {} blades: {}'.format(n, bladeMat))
+    logger.debug('stage {} blades: {}'.format(n, bladeMat))

    ULparams = ((0, TempUpper, WireMaterialUpper, BladeMaterial),
                (1, TempLower, WireMaterialLower, None))

--- a/gwinc/test/__main__.py
+++ b/gwinc/test/__main__.py
@@ -17,7 +17,7 @@ from ..io import load_hdf5

 logging.basicConfig(
    format='%(message)s',
-    level=os.getenv('LOG_LEVEL', logging.INFO))
+    level=os.getenv('LOG_LEVEL', 'INFO').upper())

 try:
    import inspiral_range
@@ -44,7 +44,7 @@ def git_find_upstream_name():
        ).stdout
    except subprocess.CalledProcessError as e:
        logging.error(e.stderr.split('\n')[0])
-    for remote in remotes.split('\n'):
+    for remote in remotes.strip().split('\n'):
        name, url, fp = remote.split()
        if 'gwinc/pygwinc.git' in url:
            return name
No results found