Add ability to load ifo data from YAML and MATLAB .mat files

This adds a Struct() class that mimics a MATLAB struct by storing values as
class attributes, but also includes methods for converting to/from dicts,
yaml, and mat_struct objects.

A convenience function load_ifo() is included to load an ifo definition
either from a .yaml or .mat file or from an included <ifo>.yaml.

An 'aLIGO.yaml' definition is included.

gwinc/__init__.py

+from ifo import load_ifo
 from util import precompIFO
 from gwinc import noise_calc
 from gwinc import gwinc

gwinc/ifo/__init__.py

+import os
+import re
+import yaml
+import scipy
+import scipy.io
+import numpy as np
+from scipy.io.matlab.mio5_params import mat_struct
+
+
+# HACK: fix loading number in scientific notation
+#
+# https://stackoverflow.com/questions/30458977/yaml-loads-5e-6-as-string-and-not-a-number
+#
+# An apparent bug in python-yaml prevents it from regognizing
+# scientific notation as a float.  The following is a modified version
+# of the parser that recognize scientific notation appropriately.
+#
+loader = yaml.SafeLoader
+loader.add_implicit_resolver(
+    u'tag:yaml.org,2002:float',
+    re.compile(u'''^(?:
+     [-+]?(?:[0-9][0-9_]*)\\.[0-9_]*(?:[eE][-+]?[0-9]+)?
+    |[-+]?(?:[0-9][0-9_]*)(?:[eE][-+]?[0-9]+)
+    |\\.[0-9_]+(?:[eE][-+][0-9]+)?
+    |[-+]?[0-9][0-9_]*(?::[0-5]?[0-9])+\\.[0-9_]*
+    |[-+]?\\.(?:inf|Inf|INF)
+    |\\.(?:nan|NaN|NAN))$''', re.X),
+    list(u'-+0123456789.'))
+
+
+class Struct(object):
+    """Matlab struct-like object
+
+    This is a simple implementation of a MATLAB struct-like object
+    that stores values as attributes of a simple class: and allows assigning to
+    attributes recursively, e.g.:
+
+    >>> s = Struct()
+    >>> s.a = 4
+    >>> s.b = Struct()
+    >>> s.b.c = 8
+
+    Various classmethods allow creating one of these objects from YAML
+    file, a nested dict, or a MATLAB struct object.
+
+    """
+
+    # FIXME: This would be a way to allow setting nested struct
+    # attributes, e.g.:
+    #
+    # >>> s = Struct()
+    # >>> s.a.b.c = 4
+    #
+    # Usage of __getattr__ like this is dangerous and creates
+    # non-intuitive behavior (i.e. an empty struct is returned with
+    # accessing attributes that don't exist).  Is there a way to
+    # accomplish this without that adverse side affect?
+    #
+    # def __getattr__(self, name):
+    #     if name not in self.__dict__:
+    #         self.__dict__[name] = Struct()
+    #     return self.__dict__[name]
+
+    ##########
+
+    def __contains__(self, item):
+        return item in self.__dict__
+
+    def to_dict(self):
+        """Return nested dictionary representation of Struct.
+
+        """
+        d = {}
+        for k,v in self.__dict__.iteritems():
+            if isinstance(v, Struct):
+                d[k] = v.to_dict()
+            else:
+                # handle lists of Structs
+                try:
+                    d[k] = [i.to_dict() for i in v]
+                except (AttributeError, TypeError):
+                    d[k] = v
+        return d
+
+    def to_yaml(self, path=None):
+        """Return YAML representation of Struct as nested dict.
+
+        Or write Struct to YAML file if file 'path' argument
+        specified.
+
+        """
+        y = yaml.dump(self.to_dict(), default_flow_style=False)
+        if path:
+            with open(path, 'w') as f:
+                f.write(y)
+        else:
+            return y
+
+    def __repr__(self):
+        return self.to_yaml().strip('\n')
+
+    def __str__(self):
+        return '<GWINC Struct: {}>'.format(self.__dict__.keys())
+
+    def __iter__(self):
+        return iter(self.__dict__)
+
+    def walk(self):
+        """Iterate over all leaves in the struct tree.
+
+        """
+        for k,v in self.__dict__.iteritems():
+            if type(v) is Struct:
+                for sk,sv in v.walk():
+                    yield k+'.'+sk, sv
+            else:
+                try:
+                    for i,vv in enumerate(v):
+                        for sk,sv in vv.walk():
+                            yield '{}[{}].{}'.format(k,i,sk), sv
+                except (AttributeError, TypeError):
+                    yield k, v
+
+    @classmethod
+    def from_dict(cls, d):
+        """Create Struct from nested dict.
+
+        """
+        c = cls()
+        for k,v in d.iteritems():
+            if type(v) == dict:
+                c.__dict__[k] = Struct.from_dict(v)
+            else:
+                try:
+                    c.__dict__[k] = map(Struct.from_dict, v)
+                except (AttributeError, TypeError):
+                    c.__dict__[k] = v
+        return c
+
+    @classmethod
+    def from_matstruct(cls, s):
+        """Create Struct from scipy.io.matlab mat_struct object.
+
+        """
+        c = cls()
+        # FIXME: handle 'ifo' attribute at top level.  should we
+        # ignore this?
+        try:
+            s = s['ifo']
+        except:
+            pass
+        for k,v in s.__dict__.iteritems():
+            if k in ['_fieldnames']:
+                # skip these fields
+                pass
+            # FIXME: manaually converting Stage list to individual
+            # named elements.  what's the right thing to do here?
+            elif k == 'Stage':
+                for i,stage in enumerate(v):
+                    name = 'Stage{}'.format(i+1)
+                    struct = Struct.from_matstruct(stage)
+                    c.__dict__[name] = struct
+            elif type(v) is mat_struct:
+                c.__dict__[k] = Struct.from_matstruct(v)
+            else:
+                # handle lists of Structs
+                try:
+                    c.__dict__[k] = map(Struct.from_matstruct, v)
+                except:
+                    c.__dict__[k] = v
+        return c
+                
+
+    @classmethod
+    def from_file(cls, path):
+        """Load Struct from .yaml or GWINC .mat file.
+
+        File type will be determined by extension.
+
+        """
+        (root, ext) = os.path.splitext(path)
+        with open(path, 'r') as f:
+            if ext in ['.yaml', '.yml']:
+                d = yaml.load(f, Loader=loader)
+                return cls.from_dict(d)
+            elif ext == '.mat':
+                s = scipy.io.loadmat(f, squeeze_me=True, struct_as_record=False)
+                return cls.from_matstruct(s)
+            else:
+                raise IOError("Unknown file type: {}".format(ext))
+
+
+
+def load_ifo(name_or_path):
+    """Load IFO by name or from file.
+
+    IFO names will correspond to basename of included .yaml IFO
+    definition file.
+
+    When specifying path may be either .yaml or .mat.
+
+    """
+    if os.path.exists(name_or_path):
+        path = name_or_path
+    else:
+        path = os.path.join(os.path.dirname(__file__),
+                            name_or_path+'.yaml')
+    s = Struct.from_file(path)
+    return s
+
+
+##################################################
+
+
+if __name__ == '__main__':
+    import sys
+    ifo = load_ifo(sys.argv[1])
+    # print(ifo.to_yaml())
+    print(ifo.to_dict())

gwinc/ifo/aLIGO.yaml

+# GWINC aLIGO interferometer parameters
+#
+# parameters for quad pendulum suspension updated 3rd May 2006, NAR
+# References:
+# LIGO-T000012-00-D
+# 	* Differentiate between silica and sapphire substrate absorption
+# 	* Change ribbon suspension aspect ratio
+# 	* Change pendulum frequency
+# References:
+# 1. Electro-Optic Handbook, Waynant & Ediger (McGraw-Hill: 1993)
+# 2. LIGO/GEO data/experience
+# 3. Suspension reference design, LIGO-T000012-00
+# 4. Quartz Glass for Optics Data and Properties, Heraeus data sheet,
+#    numbers for suprasil
+# 5. Y.S. Touloukian (ed), Thermophysical Properties of Matter
+#    (IFI/Plenum,1970)
+# 6. Marvin J. Weber (ed) CRC Handbook of laser science and technology,
+#    Vol 4, Pt 2
+# 7. R.S. Krishnan et al.,Thermal Expansion of Crystals, Pergamon Press
+# 8. P. Klocek, Handbook of infrared and optical materials, Marcel Decker,
+#    1991
+# 9. Rai Weiss, electronic log from 5/10/2006
+# 10. Wikipedia online encyclopedia, 2006
+# 11. D.K. Davies, The Generation and Dissipation of Static Charge on
+# dielectrics in a Vacuum, page 29
+# 12. Gretarsson & Harry, Gretarsson thesis
+# 13. Fejer
+# 14. Braginsky
+#
+#
+
+Constants:
+  # Temperature of the Vacuum
+  Temp: 290 # K
+
+Infrastructure:
+  Length: 3995 # m
+  ResidualGas:
+    pressure: 4.0e-7 # Pa
+    mass: 3.35e-27 # kg; Mass of H_2 (ref. 10)
+    polarizability: 7.8e-31 # m^3
+
+TCS:
+  s_cc: 7.024 # Watt^-2
+  s_cs: 7.321 # Watt^-2
+  s_ss: 7.631 # Watt^-2
+  # The hardest part to model is how efficient the TCS system is in
+  # compensating this loss. Thus as a simple Ansatz we define the
+  # TCS efficiency TCSeff as the reduction in effective power that produces
+  # a phase distortion. E.g. TCSeff=0.99 means that the compensated distortion
+  # of 1 Watt absorbed is eqivalent to the uncompensated distortion of 10mWatt.
+  # The above formula thus becomes:
+  # S= s_cc P_coat^2 + 2*s_cs*P_coat*P_subst + s_ss*P_subst^2 * (1-TCSeff)^2
+  #
+  # To avoid iterative calculation we define TCS.SCRloss = S as an input
+  # and calculate TCSeff as an output.
+  # TCS.SRCloss is incorporated as an additional loss in the SRC
+  SRCloss: 0.00
+
+Seismic:
+  Site: 'LHO'                      # LHO or LLO (only used for Newtonian noise)
+  KneeFrequency: 10                # Hz; freq where 'flat' noise rolls off
+  LowFrequencyLevel: 1e-9          # m/rtHz; seismic noise level below f_knee
+  Gamma: .8                        # abruptness of change at f_knee
+  Rho: 1.8e3                       # kg/m^3; density of the ground nearby
+  Beta: 0.5                        # quiet times beta: 0.35-0.60
+  # noisy times beta: 0.15-1.4
+  Omicron: 1                       # Feedforward cancellation factor
+
+Suspension:
+  # 0 for cylindrical suspension
+  #Type: 'Quad'
+  Type: 2
+  # 0: round, 1: ribbons
+  FiberType: 0
+  BreakStress: 750e6            # Pa; ref. K. Strain
+  Temp: 290
+  VHCoupling:
+    theta: 1e-3        # vertical-horizontal x-coupling
+
+  Silica:
+    Rho   : 2200           # Kg/m^3;
+    C     : 772            # J/Kg/K;
+    K     : 1.38           # W/m/kg;
+    Alpha : 3.9e-7         # 1/K;
+    dlnEdT: 1.52e-4        # (1/K), dlnE/dT
+    Phi   : 4.1e-10        # from G Harry e-mail to NAR 27April06 dimensionless units
+    Y     : 72e9           # Pa; Youngs Modulus
+    Dissdepth: 1.5e-2      # from G Harry e-mail to NAR 27April06
+
+  C70Steel:
+    Rho: 7800
+    C: 486
+    K: 49
+    Alpha: 12e-6
+    dlnEdT: -2.5e-4
+    Phi: 2e-4
+    Y: 212e9        # measured by MB for one set of wires
+
+  MaragingSteel:
+    Rho: 7800
+    C: 460
+    K: 20
+    Alpha: 11e-6
+    dlnEdT: 0
+    Phi: 1e-4
+    Y: 187e9
+
+  # ref ---- http://design.caltech.edu/Research/MEMS/siliconprop.html
+  # all properties should be for T ~ 20 K
+  Silicon:
+    Rho: 2330                     # Kg/m^3;  density
+    C: 772                        # J/kg/K   heat capacity
+    K: 4980                       # W/m/K    thermal conductivity
+    Alpha: 1e-9                   # 1/K      thermal expansion coeff
+    # from Gysin, et. al. PRB (2004)  E(T): E0 - B*T*exp(-T0/T)
+    # E0: 167.5e9 Pa   T0: 317 K   B: 15.8e6 Pa/K
+    dlnEdT: 2.5e-10               # (1/K)    dlnE/dT  T=20K
+    Phi: 2e-9                     # Nawrodt (2010)      loss angle  1/Q
+    Y: 150e9                      # Pa       Youngs Modulus
+    Dissdepth: 1.5e-3             # 10x smaller surface loss depth (Nawrodt (2010))
+
+  # Note stage numbering: mirror is at beginning of stack, not end
+  #
+  # last stage length adjusted for d: 10mm and and d_bend = 4mm
+  # (since 602mm is the CoM separation, and d_bend is accounted for
+  # in suspQuad, so including it here would double count)
+  Stage1:
+    Mass: 39.6                    # kg; current numbers May 2006 NAR
+    Length: 0.59                  # m
+    Dilution: .nan                #
+    K: .nan                       # N/m; vertical spring constant
+    WireRadius: .nan              # m
+    Blade: .nan                   # blade thickness
+    NWires: 4
+
+  Stage2:
+    Mass: 39.6
+    Length: 0.341
+    Dilution: 106
+    K: 5200              
+    WireRadius: 310e-6
+    Blade: 4200e-6
+    NWires: 4
+
+  Stage3:
+    Mass: 21.8
+    Length: 0.277
+    Dilution: 80
+    K: 3900
+    WireRadius: 350e-6
+    Blade: 4600e-6
+    NWires: 4
+
+  Stage4:
+    Mass: 22.1
+    Length: 0.416
+    Dilution: 87
+    K: 3400
+    WireRadius: 520e-6
+    Blade: 4300e-6
+    NWires: 2
+
+  Ribbon:
+    Thickness: 115e-6             # m
+    Width: 1150e-6                # m
+
+  Fiber:
+    Radius: 205e-6                # m
+    Blade: 4300e-6
+
+## Optic Material -------------------------------------------------------
+Materials:
+  MassRadius: 0.17                # m; 
+  MassThickness: 0.200            # m; Peter F 8/11/2005
+
+  ## Dielectric coating material parameters----------------------------------
+  Coating:
+    ## high index material: tantala
+    Yhighn: 140e9
+    Sigmahighn: 0.23
+    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
+    Phihighn: 2.3e-4
+    Indexhighn: 2.06539
+
+    ## low index material: silica
+    Ylown: 72e9
+    Sigmalown: 0.17
+    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
+    Philown: 4.0e-5
+    Indexlown: 1.45
+
+  ## Substrate Material parameters--------------------------------------------
+  Substrate:
+    c2 : 7.6e-12                 # Coeff of freq depend. term for bulk mechanical loss, 7.15e-12 for Sup2
+    MechanicalLossExponent: 0.77   # Exponent for freq dependence of silica loss, 0.822 for Sup2
+    Alphas: 5.2e-12              # Surface loss limit (ref. 12)
+    MirrorY: 7.27e10             # N/m^2; Youngs modulus (ref. 4)
+    MirrorSigma: 0.167           # Kg/m^3; Poisson ratio (ref. 4)
+    MassDensity: 2.2e3           # Kg/m^3; (ref. 4)
+    MassAlpha: 3.9e-7            # 1/K; thermal expansion coeff. (ref. 4)
+    MassCM: 739                  # J/Kg/K; specific heat (ref. 4)
+    MassKappa: 1.38              # J/m/s/K; thermal conductivity (ref. 4)
+    RefractiveIndex: 1.45        # mevans 25 Apr 2008
+
+## Laser-------------------------------------------------------------------
+Laser:
+  Wavelength: 1.064e-6 # m
+  Power: 125 # W
+
+## Optics------------------------------------------------------------------
+Optics:
+  Type: 'SignalRecycled'
+  PhotoDetectorEfficiency: 0.95    # photo-detector quantum efficiency
+  Loss: 37.5e-6 # average per mirror power loss
+  BSLoss: 0.5e-3                   # power loss near beamsplitter
+  coupling: 1.0                     # mismatch btwn arms & SRC modes; used to
+
+  #SubstrateAbsorption: 0.5e-4       # 1/m; bulk absorption coef (ref. 2)
+  SubstrateAbsorption: 0.3e-4       # 1/m; 0.3 ppm/cm for Hereaus
+  pcrit: 10                         # W; tolerable heating power (factor 1 ATC)
+  Quadrature:
+    dc: 1.5707963 # pi/2 # demod/detection/homodyne phase
+
+  ITM:
+    BeamRadius: 0.055 # m, 1/e^2 power radius
+    Transmittance: 0.014
+    CoatingThicknessLown: 0.308
+    CoatingThicknessCap: 0.5
+    CoatingAbsorption: 0.5e-6
+    SubstrateAbsorption: 0.3e-4 # 1/m, 0.3 ppm/cm for Hereaus
+  ETM:
+    BeamRadius: 0.062 # m, 1/e^2 power radius
+    Transmittance: 5e-6
+    CoatingThicknessLown: 0.27
+    CoatingThicknessCap: 0.5
+  PRM:
+    Transmittance: 0.03
+  SRM:
+    Transmittance: 0.20
+    CavityLength: 55 # m, ITM to SRM distance
+    Tunephase: 0.0 # SEC tuning
+
+  Curvature: # ROC
+    ITM: 1970
+    ETM: 2192
+
+## Squeezer Parameters------------------------------------------------------
+# Define the squeezing you want:
+#   None: ignore the squeezer settings
+#   Freq Independent: nothing special (no filter cavties)
+#   Freq Dependent = applies the specified filter cavites
+#   Optimal = find the best squeeze angle, assuming no output filtering
+#   OptimalOptimal = optimal squeeze angle, assuming optimal readout phase
+Squeezer:
+  Type: 'None'
+  AmplitudedB: 10         # SQZ amplitude [dB]
+  InjectionLoss: 0.05     # power loss to sqz
+  SQZAngle: 0             # SQZ phase [radians]
+
+  # Parameters for frequency dependent squeezing
+  FilterCavity:
+    fdetune: -14.5  # detuning [Hz]
+    L: 100          # cavity length
+    Ti: 0.12e-3     # input mirror trasmission [Power]
+    Te: 0           # end mirror trasmission
+    Lrt: 100e-6     # round-trip loss in the cavity
+    Rot: 0          # phase rotation after cavity
+
+## Variational Output Parameters--------------------------------------------
+# Define the output filter cavity chain
+#   None = ignore the output filter settings
+#   Chain = apply filter cavity chain
+#   Optimal = find the best readout phase
+OutputFilter:
+  Type: 'None'
+  FilterCavity:
+    fdetune: -30 # detuning [Hz]
+    L: 4000      # cavity length
+    Ti: 10e-3    # input mirror trasmission [Power]
+    Te: 0        # end mirror trasmission
+    Lrt: 100e-6  # round-trip loss in the cavity
+    Rot: 0       # phase rotation after cavity