From 18e3813f6dab61b8d6ac51b2e37ef17954c182cc Mon Sep 17 00:00:00 2001 From: Lee McCuller Date: Mon, 19 Oct 2020 11:26:24 -0400 Subject: [PATCH] fix to external/trf to not require scipy._six --- IIRrational/AAA.py | 8 +- IIRrational/external/scipy_optimize/trf.py | 3 +- IIRrational/file_io/hdf5_io.py | 2 +- IIRrational/pytest.py | 25 ++ IIRrational/representations/zpktf.py | 14 +- test/AAA/test_AAA.py | 2 +- test/AAA/test_AAA_present.py | 373 +++++++++++++++++++++ 7 files changed, 411 insertions(+), 16 deletions(-) create mode 100644 test/AAA/test_AAA_present.py diff --git a/IIRrational/AAA.py b/IIRrational/AAA.py index d7560f3..215e13a 100644 --- a/IIRrational/AAA.py +++ b/IIRrational/AAA.py @@ -149,7 +149,7 @@ def tf_bary_zpk( TF = np.prod([gz / gp for gz, gp in itertools.zip_longest(Gz, Gp, fillvalue = 1)]) TFvals_rel.append(f / TF) TFvals_rel = np.asarray(TFvals_rel) - print(TFvals_rel) + #print(TFvals_rel) gain = np.median(TFvals_rel.real) #this may also get computed using the gain_n/gain_d above, but that fails #when poles or zeros are dropped since one of gain_n or gain_d will be @@ -422,8 +422,8 @@ class rtAAAResults(object): pass self.fit_idx = idx self.fit_dict = self.fit_list[self.fit_idx] - self.p_order = self.fit_dict['p_order'], - self.order = self.fit_dict['order'], + self.p_order = self.fit_dict['p_order'] + self.order = self.fit_dict['order'] self.wvals = self.fit_dict['wvals'] self.zvals = self.zvals_full[:self.p_order] self.fvals = self.fvals_full[:self.p_order] @@ -532,7 +532,7 @@ def order_reduce_zp( for z, p in rpB.r12_list: Q_rank = Q_rank_calc(p, z) - print("rank: ", p, z, Q_rank) + #print("rank: ", p, z, Q_rank) #print(z, p, Q_rank) if Q_rank < Q_rank_cutoff: continue diff --git a/IIRrational/external/scipy_optimize/trf.py b/IIRrational/external/scipy_optimize/trf.py index c6b14fa..febd3f7 100644 --- a/IIRrational/external/scipy_optimize/trf.py +++ b/IIRrational/external/scipy_optimize/trf.py @@ -134,7 +134,6 @@ import declarative from numpy.linalg import norm from scipy.linalg import svd -from scipy._lib.six import string_types from .common import ( solve_lsq_trust_region, @@ -179,7 +178,7 @@ def trf( g = compute_grad(J, f) - jac_scale = isinstance(x_scale, string_types) and x_scale == 'jac' + jac_scale = isinstance(x_scale, str) and x_scale == 'jac' if jac_scale: scale, scale_inv = compute_jac_scale(J) else: diff --git a/IIRrational/file_io/hdf5_io.py b/IIRrational/file_io/hdf5_io.py index e7fc935..c47e119 100644 --- a/IIRrational/file_io/hdf5_io.py +++ b/IIRrational/file_io/hdf5_io.py @@ -14,7 +14,7 @@ def load_hdf5(fname): def write_hdf5(fname, fdict): with h5py.File(fname, 'w') as h5F: - hdf = HDFDeepBunch(h5F, writeable = True, vpath = True) + hdf = HDFDeepBunch(h5F, writeable = True) hdf.update_recursive(fdict) return diff --git a/IIRrational/pytest.py b/IIRrational/pytest.py index f992605..a6179cf 100644 --- a/IIRrational/pytest.py +++ b/IIRrational/pytest.py @@ -48,6 +48,31 @@ def relfile_test(_file_, request, pre = None, post = None, fname = None): return relfile(_file_, testname, fname = fname) class Timer(object): + def __init__(self, N): + self.N = N + + def __iter__(self): + return iter(range(self.N)) + + def __call__(self): + return self.interval / self.N + + def __float__(self): + return self.interval / self.N + + def __str__(self): + time = self() + if time > 10: + return "{:.1f}s".format(time) + elif time > 1: + return "{:.2f}s".format(time) + elif time > .001: + return "{:.1f}ms".format(time * 1e3) + elif time > 1e-6: + return "{:.1f}us".format(time * 1e6) + else: + return "{:.1f}ns".format(time * 1e9) + def __enter__(self): self.start = time.clock() return self diff --git a/IIRrational/representations/zpktf.py b/IIRrational/representations/zpktf.py index e107465..4629b0a 100644 --- a/IIRrational/representations/zpktf.py +++ b/IIRrational/representations/zpktf.py @@ -106,29 +106,27 @@ class ZPKTF(DepBunch): @depB_property def order(self): return max( - len(self.poles) + len(self.poles_overlay), - len(self.zeros) + len(self.zeros_overlay), + len(self.poles), + len(self.zeros), ) @depB_property def order_sos(self): return (max( - len(self.poles) + len(self.poles_overlay), - len(self.zeros) + len(self.zeros_overlay), + len(self.poles), + len(self.zeros), ) + 1)//2 @depB_property def order_relative(self): return ( - (len(self.zeros) + len(self.zeros_overlay)) - - (len(self.poles) + len(self.poles_overlay)) + len(self.zeros) - len(self.poles) ) @depB_property def order_total(self): return ( - len(self.poles) + len(self.poles_overlay) - + len(self.zeros) + len(self.zeros_overlay) + len(self.poles) + len(self.zeros) ) def xfer_eval(self, F_Hz): diff --git a/test/AAA/test_AAA.py b/test/AAA/test_AAA.py index 9e32558..3fc6410 100644 --- a/test/AAA/test_AAA.py +++ b/test/AAA/test_AAA.py @@ -9,7 +9,7 @@ from os import path import scipy.signal from IIRrational.pytest import ( # noqa: F401 - tpath_join, plot, pprint, tpath, tpath_preclear + tpath_join, plot, pprint, tpath, tpath_preclear, Timer ) from IIRrational.representations import asZPKTF from IIRrational.testing import IIRrational_data diff --git a/test/AAA/test_AAA_present.py b/test/AAA/test_AAA_present.py new file mode 100644 index 0000000..c6bc579 --- /dev/null +++ b/test/AAA/test_AAA_present.py @@ -0,0 +1,373 @@ +# -*- coding: utf-8 -*- +""" +""" +from __future__ import division, print_function, unicode_literals, absolute_import + +import pytest +import numpy as np +from os import path +import scipy.signal + +from IIRrational.pytest import ( # noqa: F401 + tpath_join, plot, pprint, tpath, tpath_preclear, Timer +) +from IIRrational.representations import asZPKTF + +from IIRrational.utilities.mpl import mplfigB, generate_stacked_plot_ax, asavefig +from IIRrational.utilities.np import logspaced +from IIRrational.AAA import tfAAA + +asavefig.formats.svg.use = True + + +def test_AAA_view(tpath_join, tpath_preclear, pprint): + ZPK1 = asZPKTF( + ( + ( + -1, -2, + -10, + -2+10j, -2-10j, + ), + ( + -.001, + -.1+.5j, -0.1-.5j, + -.1+5j, -0.1-5j, + -.1+50j, -0.1-50j, + -4, + ), 500 + )) + pprint("ZEROS: ", ZPK1.zeros.fullplane) + pprint("POLES: ", ZPK1.poles.fullplane) + + F_Hz = logspaced(0.5, 120, 200) + axB = fitplot(ZPK1, F_Hz, pprint = pprint) + axB.ax0.legend(framealpha=1) + axB.save(tpath_join('First_demo200pts')) + axB = fitplot(ZPK1, F_Hz, CLG = True, pprint = pprint) + axB.save(tpath_join('First_demo200pts_CLG')) + + ZPK2 = asZPKTF( + ( + ( + -20, + ), + ( + -2 + 2j, -2 - 2j, + ), 10 + )) + axB = fitplot(ZPK2, F_Hz, pprint = pprint) + axB.save(tpath_join('Second_demo200pts')) + axB = fitplot(ZPK1, F_Hz, addZPK=ZPK2, pprint = pprint) + axB.ax0.legend(framealpha=1, loc = "upper right") + axB.save(tpath_join('Third_demo200pts')) + F_Hz = logspaced(0.5, 120, 50) + axB = fitplot(ZPK1, F_Hz, addZPK=ZPK2, pprint = pprint) + axB.ax0.legend(framealpha=1, loc = "upper right") + axB.save(tpath_join('Third_demo50pts')) + return + + +def test_AAA_BNS(tpath_join, tpath_preclear, pprint): + params = dict() + params["m1"] = 30 + params["m2"] = 30 + params["f_min"] = 20 + params["f_max"] = 400 + params["deltaF"] = 1 + params["distance"] = 100e6 + F_Hz, hp, hc = gen_waveform(**params) + + axB = generate_stacked_plot_ax( + ["ax0", "ax1"], + heights_phys_in_default = 3, + heights_phys_in = {}, + height_ratios = {"ax0" : 1, "ax1" : 0.5}, + width_ratios = [1], + xscales = 'log', + xlim = None, + wspacing = .04, + hspace = 0.1, + ) + hp *= 1e21 + axB.ax0.loglog(F_Hz, abs(hp)) + axB.ax1.semilogx(F_Hz, np.angle(hp, deg = True)) + axB.save(tpath_join('BNS')) + + pprint("N points", len(F_Hz)) + + def fitplot( + TF1, + F_Hz, + N = 10, + CLG = False, + addZPK = None, + res_tol = 1e-4, + ): + timer = Timer(N) + with timer: + for _ in timer: + results = tfAAA( + F_Hz = F_Hz, + xfer = TF1, + res_tol = res_tol, + #lf_eager = True, + #degree_max = 20, + #nconv = 1, + #nrel = 10, + #rtype = 'log', + #supports = (1e-2, 1e-1, 4.2e-1, 5.5e-1, 1.5, 2.8, 1, 5e-1, 2), + ) + pprint("Time: {}".format(timer)) + pprint("weights", results.wvals) + pprint('poles', results.poles) + pprint('zeros', results.zeros) + pprint('gain', results.gain) + + TF2 = results(F_Hz) + _, TF3 = scipy.signal.freqs_zpk(results.zeros, results.poles, results.gain, worN = F_Hz) + axB = mplfigB(Nrows = 2) + + axB = generate_stacked_plot_ax( + ["ax0", "ax1"], + heights_phys_in_default = 3, + heights_phys_in = {}, + height_ratios = {"ax0" : 1, "ax1" : 0.5}, + width_ratios = [1], + xscales = 'log', + xlim = None, + wspacing = .04, + hspace = 0.1, + ) + + axB.ax0.loglog(F_Hz, abs(TF1), label = 'Model') + axB.ax0.semilogy(F_Hz, abs(TF2), label = 'Fit Bary. (order {})'.format(results.order)) + axB.ax0.semilogy(F_Hz, abs(TF3), label = 'Bary. ZPK (order {})'.format(ZPKorder(results.zpk))) + axB.ax1.semilogx(F_Hz, np.angle(TF1, deg = True)) + axB.ax1.semilogx(F_Hz, np.angle(TF2, deg = True)) + axB.ax1.semilogx(F_Hz, np.angle(TF3, deg = True)) + axB.ax0.set_ylabel("Magnitude") + axB.ax1.set_ylabel("Phase [deg]") + axB.ax_bottom.set_xlabel("Frequency") + for idx, z in enumerate(results.supports): + kw = dict() + if idx == 0: + kw['label'] = 'supports' + axB.ax0.axvline(z, lw = 0.5, ls = '--', color = 'black', **kw) + axB.ax1.axvline(z, lw = 0.5, ls = '--', color = 'black') + axB.ax0.set_title("AAA Fit, {} Points, max rel. error {:0.0e}, fit time {}".format(len(F_Hz), res_tol, timer)) + axB.ax0.legend(framealpha=1) + return axB + + axB = fitplot(TF1 = hp, F_Hz = F_Hz) + axB.save(tpath_join('BNSfit')) + + axB = fitplot(TF1 = abs(hp)**2, F_Hz = F_Hz) + axB.save(tpath_join('BNSfitPow')) + + +def gen_waveform(**params): + """Generate frequency-domain inspiral waveform + + Returns a tuple of (freq, h_plus^tilde, h_cross^tilde). + + The waveform is generated with the lalsimulation + SimInspiralChooseFDWaveform() function. Keyword arguments are + used to update the default waveform parameters (see DEFAULT_PARAMS + macro). The mass parameters ('m1' and 'm2') should be specified + in solar masses and the 'distance' parameter should be specified + in parsecs**. Waveform approximants may be given as string names + (see `lalsimulation` documentation for more info). + + For example, to generate a 20/20 Msolar BBH waveform: + + >>> hp,hc = waveform.gen_waveform('m1'=20, 'm2'=20) + + **NOTE: The requirement that masses be specified in solar masses + and distances in parsecs is different than that of the underlying + lalsimulation method which expects mass and distance parameters to + be in SI units. + + """ + import lalsimulation + from inspiral_range import waveform + from inspiral_range import const + + iparams = dict(waveform.DEFAULT_PARAMS) + iparams.update(**params) + + # convert to SI units + iparams['distance'] *= const.PC_SI + iparams['m1'] *= const.MSUN_SI + iparams['m2'] *= const.MSUN_SI + iparams['approximant'] = lalsimulation.SimInspiralGetApproximantFromString(iparams['approximant']) + + m = iparams['m1'] + iparams['m2'] + + # calculate delta F based on frequency of inner-most stable + # circular orbit ("fisco") + fisco = (const.c**3)/(const.G*(6**1.5)*2*np.pi*m) + df = 2**(np.max([np.floor(np.log(fisco/4096)/np.log(2)), -6])) + + # FIXME: are these limits reasonable? + if iparams['deltaF'] is None: + iparams['deltaF'] = df + #iparams['f_min'] = 0.1 + #iparams['f_max'] = 10000 + + hp, hc = lalsimulation.SimInspiralChooseFDWaveform(**iparams) + + freq = hp.f0 + np.arange(len(hp.data.data)) * hp.deltaF + print("f0", hp.f0) + select = abs(hp.data.data) > 0 + + return freq[select], hp.data.data[select], hc.data.data[select] + + +def ZPKorder(ZPK): + Z, P, K = ZPK + return max(len(Z), len(P)) + + +def test_AAA_GWINC(tpath_join, tpath_preclear, pprint): + import gwinc + import numpy as np + F_Hz = logspaced(5, 2000, 300) + Budget = gwinc.load_budget('aLIGO') + traces = Budget(F_Hz).run() + fig = gwinc.plot_noise(F_Hz, traces) + ax = fig.add_subplot(1, 1, 1) + ax.set_title("PyGWINC aLIGO noise budget") + ax.set_ylabel("Strain ASD [h/rtHz]") + fig.savefig(tpath_join('GWINC.pdf')) + fig.savefig(tpath_join('GWINC.svg')) + traces2 = dict() + traces3 = dict() + + timer = Timer(3) + with timer: + for N in timer: + for name, (data, other) in traces.items(): + pprint("name", name) + rescale = 1e45 + data2 = rescale * data + select = data2 > 1e-6 + results = tfAAA( + F_Hz = F_Hz[select], + xfer = data2[select], + res_tol = 1e-4, + #lf_eager = True, + #degree_max = 20, + #nconv = 1, + #nrel = 10, + #rtype = 'log', + #supports = (1e-2, 1e-1, 4.2e-1, 5.5e-1, 1.5, 2.8, 1, 5e-1, 2), + ) + + #TF2 = results(F_Hz) + _, TF3 = scipy.signal.freqs_zpk(results.zeros, results.poles, results.gain, worN = F_Hz) + label = other.get('label', name) + other2 = dict(other) + other2['label'] = label + " (order {})".format(ZPKorder(results.zpk)) + if name == 'Total': + traces2[name] = (TF3 / rescale, other) + traces3[name] = (TF3 / rescale, other2) + else: + traces2[name] = (TF3 / rescale, other2) + traces3[name] = (TF3 / rescale, other2) + + fig = gwinc.plot_noise(F_Hz, traces2) + ax = fig.add_subplot(1, 1, 1) + ax.set_title("PyGWINC aLIGO noise budget\nAAA Fit, {} Points, max rel. error {:0.0e}, fit time {}".format(len(F_Hz), 1e-4, timer)) + ax.set_ylabel("Strain ASD [h/rtHz]") + fig.savefig(tpath_join('GWINCfit.pdf')) + fig.savefig(tpath_join('GWINCfit.svg')) + + fig = gwinc.plot_noise(F_Hz, traces3) + ax = fig.add_subplot(1, 1, 1) + ax.set_title("PyGWINC aLIGO noise budget\nAAA Fit, {} Points, max rel. error {:0.0e}, fit time {}".format(len(F_Hz), 1e-4, timer)) + ax.set_ylabel("Strain ASD [h/rtHz]") + fig.savefig(tpath_join('GWINCfitUgly.pdf')) + fig.savefig(tpath_join('GWINCfitUgly.svg')) + + #axB = fitplot(TF1 = hp, F_Hz = F_Hz) + #axB.save(tpath_join('BNSfit')) + + #axB = fitplot(TF1 = abs(hp)**2, F_Hz = F_Hz) + #axB.save(tpath_join('BNSfitPow')) + + +def fitplot( + ZPK, + F_Hz, + N = 10, + CLG = False, + addZPK = None, + res_tol = 1e-4, + pprint = print, +): + TF1 = ZPK.xfer_eval(F_Hz=F_Hz) + order = ZPK.order + if CLG: + TF1 = 1 / (1 - TF1) + if addZPK: + TF1 = TF1 + addZPK.xfer_eval(F_Hz=F_Hz) + order = ZPK.order + addZPK.order + timer = Timer(N) + + with timer: + for _ in timer: + results = tfAAA( + F_Hz = F_Hz, + xfer = TF1, + res_tol = res_tol, + #lf_eager = True, + #degree_max = 20, + #nconv = 1, + #nrel = 10, + #rtype = 'log', + #supports = (1e-2, 1e-1, 4.2e-1, 5.5e-1, 1.5, 2.8, 1, 5e-1, 2), + ) + pprint("Time: {}".format(timer)) + pprint("weights", results.wvals) + pprint('poles', results.poles) + pprint('zeros', results.zeros) + pprint('gain', results.gain) + pprint('poles1', ZPK.poles.fullplane) + pprint('zeros1', ZPK.zeros.fullplane) + pprint('order1', ZPK.order, len(ZPK.zeros), len(ZPK.poles)) + pprint('order2', results.order) + + TF2 = results(F_Hz) + _, TF3 = scipy.signal.freqs_zpk(results.zeros, results.poles, results.gain, worN = F_Hz) + axB = mplfigB(Nrows = 2) + + axB = generate_stacked_plot_ax( + ["ax0", "ax1"], + heights_phys_in_default = 3, + heights_phys_in = {}, + height_ratios = {"ax0" : 1, "ax1" : 0.5}, + width_ratios = [1], + xscales = 'log', + xlim = None, + wspacing = .04, + hspace = 0.1, + ) + + axB.ax0.loglog(F_Hz, abs(TF1), label = 'Model (order {})'.format(order)) + axB.ax0.semilogy(F_Hz, abs(TF2), label = 'Fit Bary. (order {})'.format(results.order)) + axB.ax0.semilogy(F_Hz, abs(TF3), label = 'Bary. ZPK (order {})'.format(ZPKorder(results.zpk))) + axB.ax1.semilogx(F_Hz, np.angle(TF1, deg = True)) + axB.ax1.semilogx(F_Hz, np.angle(TF2, deg = True)) + axB.ax1.semilogx(F_Hz, np.angle(TF3, deg = True)) + axB.ax0.set_ylabel("Magnitude") + axB.ax1.set_ylabel("Phase [deg]") + axB.ax_bottom.set_xlabel("Frequency") + for idx, z in enumerate(results.supports): + kw = dict() + if idx == 0: + kw['label'] = 'supports' + axB.ax0.axvline(z, lw = 0.5, ls = '--', color = 'black', **kw) + axB.ax1.axvline(z, lw = 0.5, ls = '--', color = 'black') + axB.ax0.set_title("AAA Fit, {} Points, max rel. error {:0.0e}, fit time {}".format(len(F_Hz), res_tol, timer)) + axB.ax0.legend(framealpha=1) + return axB -- GitLab