overhaul tests to use new gwinc_matlab functions

This directly compares the output of gwinc() and gwinc_matlab().  Any
noises that differ by more than 1% are noted and the return value of the
test is non-zero.  Differences can be plotted with the -p option.

gwinc/test/__main__.py

 """Beginnings of a test suite for pygwinc
 
-todo:
-
- * unittest
- * write test to compare matgwinc default IFOs and pygwinc default
-
 """
 import os
+import sys
 import signal
-import scipy
-import scipy.io
+import pickle
 import numpy as np
 import matplotlib.pyplot as plt
+import argparse
 
 import logging
-logging.basicConfig(level=logging.INFO)
-
-from .. import load_ifo, noise_calc
-from ..plot import plot_noise
-from ..ifo import Struct
-# NOTE: gwinc needs to be imported before matlab for some very strange
-# reason:
-#
-# Traceback (most recent call last):
-#   File "./ifo2txt.py", line 9, in <module>
-#     import gwinc
-#   File "/home/jrollins/ligo/src/pygwinc/gwinc/__init__.py", line 9, in <module>
-#     from . import plot
-#   File "/home/jrollins/ligo/src/pygwinc/gwinc/plot.py", line 2, in <module>
-#     import matplotlib.pyplot as plt
-#   File "/usr/lib/python2.7/dist-packages/matplotlib/pyplot.py", line 29, in <module>
-#     import matplotlib.colorbar
-#   File "/usr/lib/python2.7/dist-packages/matplotlib/colorbar.py", line 32, in <module>
-#     import matplotlib.artist as martist
-#   File "/usr/lib/python2.7/dist-packages/matplotlib/artist.py", line 15, in <module>
-#     from .transforms import (Bbox, IdentityTransform, TransformedBbox,
-#   File "/usr/lib/python2.7/dist-packages/matplotlib/transforms.py", line 39, in <module>
-#     from matplotlib._path import (affine_transform, count_bboxes_overlapping_bbox,
-# ImportError: /opt/matlab/R2016b/extern/engines/python/dist/matlab/engine/glnxa64/../../../../../../../sys/os/glnxa64/libstdc++.so.6: version `CXXABI_1.3.9' not found (required by /usr/lib/python2.7/dist-packages/matplotlib/_path.x86_64-linux-gnu.so)
-#
-# no idea what that's about
-import matlab.engine
-
-
-GWINC_MAT_SCRIPT = '''
-source = SourceModel();
-source.BlackHole.Mass1 = 20; % [MSol]
-source.BlackHole.Mass2 = 20; % [MSol]
-
-flo = 5;
-fhi = 6000;
-npoints = 3000;
-f = logspace(log10(flo), log10(fhi), npoints);
-
-ifo = {ifoname}();
-ifo = precompIFO(ifo);
-[score, noises, ifo] = gwinc(f, [], ifo, source, 0);
-'''
-
-def run_gwinc(eng, ifoname, matpath):
-    # This doesn't work I think because the ifo.Suspension.Stages is not a
-    # "scalar structure":
-    # ifo = eng.IFOModel()
-    # ifo = eng.precompIFO(ifo)
-    script = GWINC_MAT_SCRIPT.format(ifoname=ifoname)
-    print()
-    print(script)
-    print()
-    eng.eval(script, nargout=0)
-    eng.save(matpath, 'score', 'noises', 'ifo', nargout=0)
+logging.basicConfig(level=logging.WARNING)
 
+from .. import load_ifo, gwinc
 
-def main():
 
-    ifoname = 'IFOModel'
-    gwincpath = os.path.join(os.path.dirname(__file__),
-                             'gwinc')
-    logging.info('gwincpath: {}'.format(gwincpath))
-    matpath = os.path.join(os.path.dirname(__file__),
-                           'gwinc.mat')
-    logging.info('matpath: {}'.format(matpath))
+IFO = 'aLIGO'
+FLO = 5
+FHI = 6000
+NPOINTS = 3000
 
-    if not os.path.exists(matpath):
-        logging.info("starting matlab engine...")
-        eng = matlab.engine.start_matlab()
-        eng.addpath(gwincpath)
-        logging.info("running gwinc matlab script...")
 
-        run_gwinc(eng, ifoname, matpath)
-        eng.exit()
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--plot', '-p', action='store_true', help='plot differences')
+    args = parser.parse_args()
+
+    ifo = load_ifo(IFO)
 
+    freq = np.logspace(np.log10(FLO), np.log10(FHI), NPOINTS)
+
+    matdata = os.path.join(os.path.dirname(__file__), 'matlab.pkl')
+    if os.path.exists(matdata):
+        logging.info("using existing {}...".format(matdata))
+        with open(matdata, 'r') as f:
+            mdata = pickle.load(f)
     else:
-        logging.info("Using existing {}...".format(matpath))
-
-    logging.info("loading {}...".format(matpath))
-    mat_dict = scipy.io.loadmat(matpath,
-                                squeeze_me=True,
-                                struct_as_record=False)
-
-    ifo = Struct.from_matstruct(mat_dict['ifo'])
-    matnoises = Struct.from_matstruct(mat_dict['noises']).to_dict()
-    f = matnoises['Freq']
-
-    # FIXME: pygwinc expects some newer attributes, for newer
-    # suspensionthermal calc that requires stage dilution, among other
-    # things.  Should make sure it's fully backwards compatible.
-    ifo.Seismic.Omicron = 1
-    ifo.Suspension.FiberType = 0
-    ifo.Suspension.Stage[0].Dilution = np.NaN
-    ifo.Suspension.Stage[1].Dilution = 106
-    ifo.Suspension.Stage[2].Dilution = 80
-    ifo.Suspension.Stage[3].Dilution = 87
-
-    logging.info("calculating noise...")
-    noises = noise_calc(ifo, f)
-
-    # plt.figure()
-    # plot_noise(noises)
-    # plt.title('pygwinc with matgwinc ifo')
-    # plt.show(block=False)
-
-    # mapping from pygwinc names to matgwinc names
-    MAPPING = {
-        'Quantum Vacuum': 'Quantum',
-        'Seismic': 'Seismic',
-        'Newtonian Gravity': 'Newtonian',
-        'Suspension Thermal': 'SuspThermal',
-        'Coating Brownian': 'MirrorThermal',
-        # 'Coating Thermo-Optic': 
-        # 'ITM Thermo-Refractive':
-        # 'ITM Carrier Density':
-        'Substrate Brownian': 'SubBrown',
-        'Substrate Thermo-Elastic': 'SubTE',
-        'Excess Gas': 'ResGas',
-        }
+        logging.info("calculating matlab noise...")
+        gwincpath = os.path.join(os.path.dirname(__file__), 'gwinc')
+        from ..gwinc_matlab import gwinc_matlab
+        score, noises, ifo = gwinc_matlab(freq, ifo, gwincpath=gwincpath)
+        mdata = dict(score=score, noises=noises, ifo=ifo)
+        with open(matdata, 'w') as f:
+            pickle.dump(mdata, f)
+
+
+    score, noises, ifo = gwinc(freq, ifo)
+
+    mnoises = mdata['noises']
 
     diffs = {}
     for name, noise in noises.iteritems():
+        if name == 'Freq':
+            continue
         try:
-            matnoise = matnoises[MAPPING[name]]
+            mnoise = mnoises[name]
         except KeyError:
             continue
-        logging.info("plot compare: {}".format(name))
-        # FIXME: for CoatTO
-        if isinstance(matnoise, dict):
-            logging.info("  skip")
+        # logging.info("compare: {}".format(name))
+
+        diff = np.sqrt(mnoise) - np.sqrt(noise)
+        frac = abs(diff / np.sqrt(noise))
+
+        if max(frac) < 0.01:
             continue
-        diff = matnoise - noise
-        logging.info("  max: {}, min: {}".format(max(diff), min(diff)))
-        diffs[name] = diff
 
-    fig, ax = plt.subplots(len(diffs), sharex=True)
-    for i, name in enumerate(diffs):
-        diff = diffs[name]
-        ax[i].semilogx(f, diff, label=name)
-        ax[i].legend(loc='lower right')
-    ax[0].set_title('difference: (mat_gwinc - pygwinc)')
-    plt.xlabel('frequency [Hz]')
-    plt.show()
+        logging.warning("Excessive difference: {}".format(name))
+        logging.warning("  max: {:e}, min: {:e}".format(max(frac), min(frac)))
+
+        diffs[name] = frac
+
+
+    if args.plot:
+        spec = (len(diffs), 2)
+        for i, name in enumerate(diffs):
+            ax = plt.subplot2grid(spec, (i, 0))
+            ax.loglog(freq, np.sqrt(noises[name]), label='pygwinc')
+            ax.loglog(freq, np.sqrt(mnoises[name]), label='matlab')
+            ax.grid()
+            ax.legend(loc='upper right')
+            # ax.set_title(name)
+            ax.set_ylabel(name)
+    
+            ax = plt.subplot2grid(spec, (i, 1))
+            ax.loglog(freq, diffs[name], label=name)
+            ax.grid()
+            # ax.set_title(name)
+    
+        plt.suptitle("noises that differ by more than 1% [(mat-py)/py]")
+        plt.show()
+
+
+    if len(diffs) > 0:
+        return 1
+    return 0
 
 
 if __name__ == '__main__':
     signal.signal(signal.SIGINT, signal.SIG_DFL)
-    main()
+    sys.exit(main())