From 83d7387e99da63c842954759054a773bf49efc8f Mon Sep 17 00:00:00 2001
From: Jameson Graef Rollins <jrollins@finestructure.net>
Date: Fri, 15 Dec 2017 17:12:01 -0800
Subject: [PATCH] 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 | 210 ++++++++++++++---------------------------
1 file changed, 73 insertions(+), 137 deletions(-)
diff --git a/gwinc/test/__main__.py b/gwinc/test/__main__.py
index 117a210..8efaa40 100644
--- a/gwinc/test/__main__.py
+++ b/gwinc/test/__main__.py
@@ -1,165 +1,101 @@
"""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())
--
GitLab