From f1b07d74b45a7ff08a8f23263cf748812f5394ac Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Wed, 15 Oct 2014 22:08:58 +0000 Subject: [PATCH] ongoing work --- pywafo/src/wafo/covariance/estimation.py | 2 +- .../wafo/covariance/test/test_covariance.py | 29 ++++++--- pywafo/src/wafo/misc.py | 2 +- pywafo/src/wafo/polynomial.py | 33 +++++----- pywafo/src/wafo/spectrum/core.py | 65 ++++++++++--------- 5 files changed, 74 insertions(+), 57 deletions(-) diff --git a/pywafo/src/wafo/covariance/estimation.py b/pywafo/src/wafo/covariance/estimation.py index 4a2c87c..c538da2 100644 --- a/pywafo/src/wafo/covariance/estimation.py +++ b/pywafo/src/wafo/covariance/estimation.py @@ -76,7 +76,7 @@ class CovarianceEstimatior(object): Lmax = min(300, len(R) - 1) # maximum lag if L is undetermined # finding where ACF is less than 2 st. deviations. sigma = np.sqrt(np.r_[0, R[0] ** 2, - R[0] ** 2 + 2 * np.cumsum(R[1:] ** 2)] / Ncens) + R[0] ** 2 + 2 * np.cumsum(R[1:] ** 2)] / Ncens) lag = Lmax + 2 - (np.abs(R[Lmax::-1]) > 2 * sigma[Lmax::-1]).argmax() if self.window == 'parzen': lag = int(4 * lag / 3) diff --git a/pywafo/src/wafo/covariance/test/test_covariance.py b/pywafo/src/wafo/covariance/test/test_covariance.py index 3b6aa2b..8e34eb3 100644 --- a/pywafo/src/wafo/covariance/test/test_covariance.py +++ b/pywafo/src/wafo/covariance/test/test_covariance.py @@ -1,17 +1,28 @@ -from numpy.testing import (run_module_suite, assert_equal, assert_almost_equal, - assert_array_equal, assert_array_almost_equal) -import numpy as np -from numpy import array, cos, exp, linspace, pi, sin, diff, arange, ones - +from numpy.testing import (run_module_suite, assert_equal, + assert_array_almost_equal, + assert_almost_equal, assert_array_equal) # @UnusedImport @IgnorePep8 import wafo.spectrum.models as sm -from wafo.covariance import CovData1D +# from wafo.covariance import CovData1D + def test_covariance(): Sj = sm.Jonswap() - S = Sj.tospecdata() #Make spec + S = Sj.tospecdata() # Make spec R = S.tocovdata() - - x = R.sim(ns=1000,dt=0.2) + dt = R.sampling_period() + assert_equal(dt, 1.0471975511965976) + S1 = R.tospecdata() + assert_array_almost_equal(S.data[:10], S1.data[:10], 11) + x = R.sim(ns=1000, dt=0.2, iseed=0) + assert_array_almost_equal(x[:10, 0], [0.0, 1.04719755, 2.0943951, + 3.14159265, 4.1887902, 5.23598776, + 6.28318531, 7.33038286, 8.37758041, + 9.42477796], decimal=3) + assert_array_almost_equal(x[:10, 1], [0.22155905, 1.21207066, 1.95670282, + 2.11634902, 1.57967273, 0.2665005, + -0.79630253, -1.31908028, + -2.20056021, -1.84451748], decimal=3) + if __name__ == '__main__': run_module_suite() diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index f256b2f..77ef33a 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -3398,4 +3398,4 @@ def test_hyp2f1(): if __name__ == "__main__": test_docstrings() - # test_hyp2f1() + # test_hyp2f1() \ No newline at end of file diff --git a/pywafo/src/wafo/polynomial.py b/pywafo/src/wafo/polynomial.py index e060fe4..7294157 100644 --- a/pywafo/src/wafo/polynomial.py +++ b/pywafo/src/wafo/polynomial.py @@ -24,7 +24,7 @@ from numpy.fft import fft, ifft from numpy import (zeros, ones, zeros_like, array, asarray, newaxis, arange, logical_or, any, pi, cos, round, diff, all, exp, where, extract, linalg, sign, concatenate, floor, isreal, - conj, remainder, linspace, atleast_1d, hstack, sum) + conj, remainder, linspace, sum) from numpy.lib.polynomial import * # @UnusedWildImport from scipy.misc.common import pade # @UnresolvedImport __all__ = np.lib.polynomial.__all__ @@ -272,8 +272,8 @@ def polydeg(x, y): p = orthofit(x, y, n) ys = orthoval(p, x) # -- Akaike's Information Criterion - aic = 2 * (n + 1) * (1 + (n + 2) / (N - n - 2)) + \ - N * (np.log(2 * pi * sum((ys - y) ** 2) / N) + 1) + aic = (2 * (n + 1) * (1 + (n + 2) / (N - n - 2)) + + N * (np.log(2 * pi * sum((ys - y) ** 2) / N) + 1)) if aic >= AIC: nit += 1 @@ -634,17 +634,17 @@ def poly2hstr(p, variable='x'): ix = 1 for expon in range(order, -1, -1): coef = coefs[order - expon] - #% There is no point in adding a zero term (except if it's the only - #% term, but we'll take care of that later). + # There is no point in adding a zero term (except if it's the only + # term, but we'll take care of that later). if coef == 0: ix += 1 else: - #% Append exponent if necessary. + # Append exponent if necessary. if ix > 1: exponstr = '%.0f' % ix s = '%s**%s' % (s, exponstr) ix = 1 - #% Is it the first term? + # Is it the first term? isfirst = s == '' # We need the coefficient only if it is different from 1 or -1 or @@ -656,7 +656,7 @@ def poly2hstr(p, variable='x'): # We need the variable except in the constant term. needvar = (expon != 0) - #% Add sign, but we don't need a leading plus-sign. + # Add sign, but we don't need a leading plus-sign. if isfirst: if coef < 0: s = '-' # % Unary minus. @@ -666,22 +666,22 @@ def poly2hstr(p, variable='x'): else: s = '%s + ' % s # % Binary plus (addition). - #% Append the coefficient if it is different from one or when it is - #% the constant term. + # Append the coefficient if it is different from one or when it is + # the constant term. if needcoef: coefstr = '%.20g' % abs(coef) s = '%s%s' % (s, coefstr) - #% Append variable if necessary. + # Append variable if necessary. if needvar: - #% Append a multiplication sign if necessary. + # Append a multiplication sign if necessary. if needcoef: if 1 - isfirst: s = '(%s)' % s s = '%s*' % s s = '%s%s' % (s, var) - #% Now treat the special case where the polynomial is zero. + # Now treat the special case where the polynomial is zero. if s == '': s = '0' return s @@ -1170,7 +1170,6 @@ def chebfit(fun, n=10, a=-1, b=1, trace=False): else: f = fun n = len(f) - #raise ValueError('Function must be callable!') # N-1 # c(k) = (2/N) sum w(n) f(n)*cos(pi*k*(2n+1)/(2N)), 0 <= k < N. # n=0 @@ -1618,7 +1617,7 @@ class Cheb1d(object): return self.coeffs[val] def __setitem__(self, key, val): - #ind = self.order - key + # ind = self.order - key if key < 0: raise ValueError("Does not support negative powers.") if key > self.order: @@ -1921,7 +1920,7 @@ def test_polydeg(): x = np.linspace(0, 10, 300) y = np.sin(x ** 3 / 100) ** 2 + 0.05 * np.random.randn(x.size) n = polydeg(x, y) - #n = 2 + # n = 2 p = orthofit(x, y, n) xi = linspace(x.min(), x.max()) ys0 = orthoval(p, x) @@ -1946,4 +1945,4 @@ if __name__ == '__main__': main() else: test_docstrings() - #test_polydeg() + # test_polydeg() diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index bb5e004..deddc8f 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -1,6 +1,6 @@ from __future__ import division from wafo.misc import meshgrid, gravity, cart2polar, polar2cart -from wafo.objects import TimeSeries # mat2timeseries, +from wafo.objects import TimeSeries, mat2timeseries import warnings import os import numpy as np @@ -23,6 +23,7 @@ from wafo.containers import PlotData, now from wafo.misc import sub_dict_select, nextpow2, discretize, JITImport # from wafo.graphutil import cltext from wafo.kdetools import qlevels +from scipy.interpolate.interpolate import interp1d try: from wafo.gaussian import Rind @@ -55,7 +56,7 @@ __all__ = ['SpecData1D', 'SpecData2D', 'plotspec'] def _set_seed(iseed): '''Set seed of random generator''' - if iseed != None: + if iseed is not None: try: random.set_state(iseed) except: @@ -136,7 +137,7 @@ def qtf(w, h=inf, g=9.81): # tmp1 = 0.5*g*k_w./(w.*sqrt(g*h)) # tmp2 = 0.25*w.^2/g -# Wave group velocity + # Wave group velocity c_g = 0.5 * g * (tanh(k_w * h) + k_w * h * (1.0 - tanh(k_w * h) ** 2)) / w h_dii = (0.5 * (0.5 * g * (k_w / w) ** 2. - 0.5 * w ** 2 / g + g * k_w / (w * c_g)) @@ -146,7 +147,7 @@ def qtf(w, h=inf, g=9.81): # k = find(w_1==w_2) # h_d(k) = h_dii - #% The NaN's occur due to division by zero. => Set the isnans to zero + # The NaN's occur due to division by zero. => Set the isnans to zero h_dii = where(isnan(h_dii), 0, h_dii) h_d = where(isnan(h_d), 0, h_d) @@ -197,8 +198,8 @@ def plotspec(specdata, linetype='b-', flag=1): # # S = demospec('dir'); S2 = mkdspec(jonswap,spreading); # plotspec(S,2), hold on -# plotspec(S,3,'g') % Same as previous fig. due to frequency independent spreading -# plotspec(S2,2,'r') % Not the same as previous figs. due to frequency dependent spreading +# plotspec(S,3,'g') # Same as previous fig. due to frequency independent spreading +# plotspec(S2,2,'r') # Not the same as previous figs. due to frequency dependent spreading # plotspec(S2,3,'m') # % transform from angular frequency and radians to frequency and degrees # Sf = ttspec(S,'f','d'); clf @@ -1004,8 +1005,7 @@ class SpecData1D(PlotData): max_sim = 30 tolerance = 5e-4 - L = 200 # %maximum lag size of the window function used in - #%spectral estimate + L = 200 # maximum lag size of the window function used in estimate # ftype = self.freqtype #options are 'f' and 'w' and 'k' # switch ftype # case 'f', @@ -1015,8 +1015,8 @@ class SpecData1D(PlotData): Hm0 = self.characteristic('Hm0') Tm02 = self.characteristic('Tm02') - if not iseed is None: - _set_seed(iseed) # % set the the seed + if iseed is not None: + _set_seed(iseed) # set the the seed n = len(self.data) if ns is None: @@ -1050,25 +1050,27 @@ class SpecData1D(PlotData): SL.data[indZero] = 0 maxS = max(S.data) - # Fs = 2*freq(end)+eps % sampling frequency + # Fs = 2*freq(end)+eps # sampling frequency for ix in xrange(max_sim): x2, x1 = self.sim_nl(ns=np, cases=cases, dt=None, iseed=iseed, - method=method, - fnlimit=fn_limit) - #%x2(:,2:end) = x2(:,2:end) -x1(:,2:end) + method=method, fnlimit=fn_limit, + output='timeseries') + x2.data -= x1.data # x2(:,2:end) = x2(:,2:end) -x1(:,2:end) + S2 = x2.tospecdata(L) + S1 = x1.tospecdata(L) + # TODO: Finish spec.to_linspec # S2 = dat2spec(x2, L) # S1 = dat2spec(x1, L) # %[tf21,fi] = tfe(x2(:,2),x1(:,2),1024,Fs,[],512) # %Hw11 = interp1q(fi,tf21.*conj(tf21),freq) -# if True: -# Hw1 = exp(interp1q(S2.args, log(abs(S1.data / S2.data)), freq)) -# else: -# Geometric mean -# Hw1 = exp( -# (interp1q(S2.args, log(abs(S1.data / S2.data)), -# freq) + log(Hw2)) / 2) + if True: + Hw1 = exp(interp1d(log(abs(S1.data / S2.data)), S2.args)(freq)) + else: + # Geometric mean + Hw1 = exp((interp1d(log(abs(S1.data / S2.data)), S2.args)(freq) + + log(Hw2)) / 2) # end # Hw1 = (interp1q( S2.w,abs(S1.S./S2.S),freq)+Hw2)/2 # plot(freq, abs(Hw11-Hw1),'g') @@ -1085,7 +1087,7 @@ class SpecData1D(PlotData): # end k = nonzero(SL.data < 0)[0] if len(k): # Make sure that the current guess is larger than zero - #%k + # k # Hw1(k) Hw1[k] = min(S1.data[k] * 0.9, S.data[k]) SL.data[k] = max(Hw1[k] * S.data[k], eps) @@ -1120,7 +1122,7 @@ class SpecData1D(PlotData): #%Hw1(end) #%maxS*1e-3 - #%if Hw1(end)*S.>maxS*1e-3, + #%if Hw1[-1]*S.data>maxS*1e-3, #% warning('The Nyquist frequency of the spectrum may be too low') #%end @@ -2639,7 +2641,8 @@ class SpecData1D(PlotData): # function [x2,x,svec,dvec,amp]=spec2nlsdat(spec,np,dt,iseed,method, # truncationLimit) def sim_nl(self, ns=None, cases=1, dt=None, iseed=None, method='random', - fnlimit=1.4142, reltol=1e-3, g=9.81, verbose=False): + fnlimit=1.4142, reltol=1e-3, g=9.81, verbose=False, + output='timeseries'): """ Simulates a Randomized 2nd order non-linear wave X(t) @@ -2894,7 +2897,7 @@ class SpecData1D(PlotData): (f_limit_lo, f_limit_up)) # if nargout>3, -# %compute the sum and frequency effects separately +# #compute the sum and frequency effects separately # [svec, dvec] = disufq((amp.'),w,kw,min(h,10^30),g,nmin,nmax) # svec = svec.' # dvec = dvec.' @@ -2902,7 +2905,7 @@ class SpecData1D(PlotData): # x2s = fft(svec) % 2'nd order sum frequency component # x2d = fft(dvec) % 2'nd order difference frequency component ## -# % 1'st order + 2'nd order component. +# # 1'st order + 2'nd order component. # x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:)) # else if False: @@ -2926,7 +2929,10 @@ class SpecData1D(PlotData): # 1'st order + 2'nd order component. x2[:, 1::] = x[:, 1::] + x2o[0:ns, :].real - + if output=='timeseries': + xx2 = mat2timeseries(x2[:, 1::], x2[:, 0].ravel()) + xx = mat2timeseries(x[:, 1::], x[:, 0].ravel()) + return xx2, xx return x2, x def stats_nl(self, h=None, moments='sk', method='approximate', g=9.81): @@ -4323,6 +4329,7 @@ def test_mm_pdf(): w = np.linspace(0, 4, 256) S1 = Sj.tospecdata(w) # Make spectrum object from numerical values S = sm.SpecData1D(Sj(w), w) # Alternatively do it manually + S0 = S.to_linspec() mm = S.to_mm_pdf() @@ -4332,6 +4339,6 @@ def test_docstrings(): if __name__ == '__main__': - test_docstrings() - # test_mm_pdf() + #test_docstrings() + test_mm_pdf() # main()