diff --git a/pywafo/src/wafo/demos.py b/pywafo/src/wafo/demos.py new file mode 100644 index 0000000..38d50ab --- /dev/null +++ b/pywafo/src/wafo/demos.py @@ -0,0 +1,49 @@ +''' +Created on 20. jan. 2011 + +@author: pab +''' +import numpy as np +from numpy import exp +from wafo.misc import meshgrid +__all__ = ['peaks', 'humps'] + + +def peaks(x=None, y=None, n=51): + ''' + Return the "well" known MatLab (R) peaks function + evaluated in the [-3,3] x,y range + + Example + ------- + >>> import pylab as plt + >>> x,y,z = peaks() + >>> h = plt.contourf(x,y,z) + + ''' + if x is None: + x = np.linspace(-3, 3, n) + if y is None: + y = np.linspace(-3, 3, n) + + [x1, y1] = meshgrid(x, y) + + z = (3 * (1 - x1) ** 2 * exp(-(x1 ** 2) - (y1 + 1) ** 2) + - 10 * (x1 / 5 - x1 ** 3 - y1 ** 5) * exp(-x1 ** 2 - y1 ** 2) + - 1. / 3 * exp(-(x1 + 1) ** 2 - y1 ** 2)) + + return x1, y1, z + +def humps(x=None): + ''' + Computes a function that has three roots, and some humps. + ''' + if x is None: + y = np.linspace(0, 1) + else: + y = np.asarray(x) + + return 1.0 / ((y - 0.3) ** 2 + 0.01) + 1.0 / ((y - 0.9) ** 2 + 0.04) + 2 * y - 5.2 + +if __name__ == '__main__': + pass \ No newline at end of file diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py index b93be3a..b39ae52 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py @@ -32,7 +32,6 @@ from pylab import * #! from commands used in Chapter 2 #! -pstate = 'off'; #! Section 2.1 Introduction and preliminary analysis #!==================================================== diff --git a/pywafo/src/wafo/graphutil.py b/pywafo/src/wafo/graphutil.py new file mode 100644 index 0000000..519d1a5 --- /dev/null +++ b/pywafo/src/wafo/graphutil.py @@ -0,0 +1,110 @@ +''' +Created on 20. jan. 2011 + +@author: pab + +license BSD +''' +import warnings +import numpy as np +from wafo.plotbackend import plotbackend + +__all__ = ['cltext'] + +def cltext(levels, percent=False, n=4, xs=0.036, ys=0.94, zs=0): + ''' + Places contour level text in the current window + + Parameters + ---------- + levels = vector of contour levels or the corresponding percent which the + contour line encloses + percent = 0 if levels are the actual contour levels (default) + 1 if levels are the corresponding percent which the + contour line encloses + n = maximum N digits of precision (default 4) + Returns + h = handles to the text objects. + CLTEXT creates text objects in the current figure and prints + "Level curves at:" if percent is False and + "Level curves enclosing:" otherwise + and the contour levels or percent. + + NOTE: + -The handles to the lines of text may also be found by + h = findobj(gcf,'gid','CLTEXT','type','text'); + h = findobj(gca,'gid','CLTEXT','type','text'); + -To make the text objects follow the data in the axes set the units + for the text objects 'data' by + set(h,'unit','data') + + Examples: + >>> from wafo.integrate import peaks + >>> import pylab as plt + >>> x,y,z = peaks(); + >>> h = plt.contour(x,y,z) + >>> h = cltext(h.levels) + >>> plt.show() + + data = rndray(1,2000,2); f = kdebin(data,{'kernel','epan','L2',.5,'inc',128}); + contour(f.x{:},f.f,f.cl),cltext(f.pl,1) + + See also + pdfplot + ''' + # TODO : Make it work like legend does (but without the box): include position options etc... + clevels = np.atleast_1d(levels) + _CLTEXT_TAG = 'CLTEXT' + cax = plotbackend.gca() + axpos = cax.get_position() + xint = axpos.intervalx + yint = axpos.intervaly + + xss = xint[0] + xs * (xint[1] - xint[0]) + yss = yint[0] + ys * (yint[1] - yint[0]) + + cf = plotbackend.gcf() # get current figure + #% delete cltext object if it exists + #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + def matchfun(x): + if hasattr(x, 'get_gid'): + return x.get_gid() == _CLTEXT_TAG + return False + h_cltxts = plotbackend.findobj(cf, matchfun); + if len(h_cltxts): + for i in h_cltxts: + try: + cax.texts.remove(i) + except: + warnings.warn('Tried to delete a non-existing CL-text') + + try: + cf.texts.remove(i) + except: + warnings.warn('Tried to delete a non-existing CL-text') + + charHeight = 1 / 33; + delta_y = charHeight; + + if percent: + titletxt = 'Level curves enclosing:'; + else: + titletxt = 'Level curves at:'; + + format = '%0.' + ('%d' % n) + 'g\n' + + cltxt = ''.join([format % level for level in clevels.tolist()]) + + titleProp = dict(gid=_CLTEXT_TAG, horizontalalignment='left', + verticalalignment='center', fontweight='bold', axes=cax) # + ha1 = plotbackend.figtext(xss, yss, titletxt, **titleProp) + + yss -= delta_y; + txtProp = dict(gid=_CLTEXT_TAG, horizontalalignment='left', + verticalalignment='top', axes=cax) + + ha2 = plotbackend.figtext(xss, yss, cltxt, **txtProp) + + return ha1, ha2 +if __name__ == '__main__': + pass \ No newline at end of file diff --git a/pywafo/src/wafo/integrate.py b/pywafo/src/wafo/integrate.py index 4fd5190..f99e79f 100644 --- a/pywafo/src/wafo/integrate.py +++ b/pywafo/src/wafo/integrate.py @@ -2,60 +2,20 @@ from __future__ import division import warnings import copy import numpy as np -from numpy import pi, sqrt, ones, zeros, exp #@UnresolvedImport +from numpy import pi, sqrt, ones, zeros #@UnresolvedImport from scipy import integrate as intg import scipy.special.orthogonal as ort from scipy import special as sp import pylab as plb -from wafo.misc import meshgrid - +from scipy.integrate import simps, trapz +from wafo.misc import is_numlike +from wafo.demos import humps _POINTS_AND_WEIGHTS = {} -__all__ = ['peaks', 'humps', 'is_numlike', 'dea3', 'clencurt', 'romberg', +__all__ = ['dea3', 'clencurt', 'romberg', 'h_roots','j_roots', 'la_roots','p_roots','qrule', 'gaussq', 'richardson', 'quadgr', 'qdemo'] -def peaks(x=None, y=None, n=51): - ''' - Return the "well" known MatLab (R) peaks function - evaluated in the [-3,3] x,y range - - Example - ------- - >>> import pylab as plt - >>> x,y,z = peaks() - >>> h = plt.contourf(x,y,z) - - ''' - if x is None: - x = np.linspace(-3, 3, n) - if y is None: - y = np.linspace(-3, 3, n) - - [x1, y1] = meshgrid(x, y) - - z = (3 * (1 - x1) ** 2 * exp(-(x1 ** 2) - (y1 + 1) ** 2) - - 10 * (x1 / 5 - x1 ** 3 - y1 ** 5) * exp(-x1 ** 2 - y1 ** 2) - - 1. / 3 * exp(-(x1 + 1) ** 2 - y1 ** 2)) - return x1, y1, z - -def humps(x=None): - ''' - Computes a function that has three roots, and some humps. - ''' - if x is None: - y = np.linspace(0, 1) - else: - y = np.asarray(x) - - return 1.0 / ((y - 0.3) ** 2 + 0.01) + 1.0 / ((y - 0.9) ** 2 + 0.04) + 2 * y - 5.2 - -def is_numlike(obj): - 'return true if *obj* looks like a number' - try: - obj + 1 - except TypeError: return False - else: return True def dea3(v0, v1, v2): ''' @@ -530,7 +490,6 @@ def j_roots(n, alpha, beta, method='newton'): a = 2. * j * (j + alfbet) * tmp tmp = tmp + 2 c = 2 * (j - 1 + alpha) * (j - 1 + beta) * tmp - b = (tmp - 1) * (alpha ** 2 - beta ** 2 + tmp * (tmp - 2) * z) L[kp1, :] = (b * L[k0, :] - c * L[km1, :]) / a diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 5ea0724..c467ebe 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -23,7 +23,7 @@ except: floatinfo = finfo(float) -__all__ = ['JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_select', +__all__ = ['is_numlike','JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_select', 'parse_kwargs', 'detrendma', 'ecross', 'findcross', 'findextrema', 'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc', 'findoutliers', 'common_shape', 'argsreduce', @@ -31,6 +31,13 @@ __all__ = ['JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_select', 'discretize', 'discretize2', 'pol2cart', 'cart2pol', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'histgrm', 'num2pistr'] +def is_numlike(obj): + 'return true if *obj* looks like a number' + try: + obj + 1 + except TypeError: return False + else: return True + class JITImport(object): ''' Just In Time Import of module diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 31a264f..5defa28 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -28,7 +28,7 @@ from numpy.fft import fft from numpy.random import randn from scipy.integrate import trapz from pylab import stineman_interp -from matplotlib.mlab import psd +from matplotlib.mlab import psd, detrend_mean import scipy.signal @@ -40,6 +40,9 @@ from wafodata import WafoData from plotbackend import plotbackend import matplotlib from scipy.stats.stats import skew, kurtosis +from scipy.signal.windows import parzen +from scipy import special + matplotlib.interactive(True) _wafocov = JITImport('wafo.covariance') _wafospec = JITImport('wafo.spectrum') @@ -47,6 +50,8 @@ _wafospec = JITImport('wafo.spectrum') __all__ = ['TimeSeries', 'LevelCrossings', 'CyclePairs', 'TurningPoints', 'sensortypeid', 'sensortype'] +def _invchi2(q, df): + return special.chdtri(df, q) class LevelCrossings(WafoData): ''' @@ -856,7 +861,7 @@ class TimeSeries(WafoData): acf.norm = norm return acf - def tospecdata(self, *args, **kwargs): + def tospecdata(self, L=None, tr=None, method='cov', detrend=detrend_mean, window=parzen, noverlap=0, pad_to=None): """ Return power spectral density by Welches average periodogram method. @@ -866,9 +871,6 @@ class TimeSeries(WafoData): if len(data) < NFFT, it will be zero padded to `NFFT` before estimation. Must be even; a power 2 is most efficient. detrend : function - Fs : real, scalar - sampling frequency (samples per time unit). - window : vector of length NFFT or function To create window vectors see numpy.blackman, numpy.hamming, numpy.bartlett, scipy.signal, scipy.signal.get_window etc. @@ -893,11 +895,180 @@ class TimeSeries(WafoData): Bendat & Piersol (1986) Random Data: Analysis and Measurement Procedures, John Wiley & Sons """ - fs = 1. / (2 * self.sampling_period()) - S, f = psd(self.data.ravel(), Fs=fs, *args, **kwargs) + dt = self.sampling_period() + #fs = 1. / (2 * dt) + yy = self.data.ravel() if tr is None else tr.dat2gauss(self.data.ravel()) + yy = detrend(yy) if hasattr(detrend,'__call__') else yy + + S, f = psd(yy, Fs=1./dt, NFFT=nfft, detrend=detrend, window=window, + noverlap=noverlap,pad_to=pad_to, scale_by_freq=True) fact = 2 * 2.0 * pi w = fact * f return _wafospec.SpecData1D(S / fact, w) + def specdata(self, L=None, tr=None, method='cov',dflag='mean', ftype='w'): + ''' + Estimate one-sided spectral density from data. + + Parameters + ---------- + L : scalar integer + maximum lag size of the window function. As L decreases the estimate + becomes smoother and Bw increases. If we want to resolve peaks in + S which is Bf (Hz or rad/sec) apart then Bw < Bf. If no value is given the + lag size is set to be the lag where the auto correlation is less than + 2 standard deviations. (maximum 300) + tr : transformation object + the transformation assuming that x is a sample of a transformed + Gaussian process. If g is None then x is a sample of a Gaussian process (Default) + method : string + defining estimation method. Options are + 'cov' : Frequency smoothing using a parzen window function + on the estimated autocovariance function. (default) + 'psd' : Welch's averaged periodogram method with no overlapping batches + dflag : string + defining detrending performed on the signal before estimation. + 'mean','linear' or 'ma' (= moving average) (default 'mean') + ftype : character + defining frequency type: 'w' or 'f' (default 'w') + + Returns + --------- + spec : SpecData1D object + + + Example + ------- + x = load('sea.dat'); + S = dat2spec(x); + specplot(S) + + See also + -------- + dat2tr, dat2cov + + + References: + ----------- + Georg Lindgren and Holger Rootzen (1986) + "Stationära stokastiska processer", pp 173--176. + + Gareth Janacek and Louise Swift (1993) + "TIME SERIES forecasting, simulation, applications", + pp 75--76 and 261--268 + + Emanuel Parzen (1962), + "Stochastic Processes", HOLDEN-DAY, + pp 66--103 + ''' + + #% Initialize constants + #%~~~~~~~~~~~~~~~~~~~~~ + nugget = 0; #%10^-12; + rate = 2; #% interpolationrate for frequency + tapery = 0; #% taper the data before the analysis + wdef = 1; #% 1=parzen window 2=hanning window, 3= bartlett window + + dt = self.sampling_period() + yy = self.data if tr is None else tr.dat2gauss(self.data) + n = len(yy) + L = min(L,n); + + dflag = dflag.lower() + if dflag=='mean': + yy -= yy.mean() + elif dflag== 'linear': + yy = detrend(yy,1); # % signal toolbox detrend + elif dflag== 'ma': + dL = np.ceil(1200/2/dT); # % approximately 20 min. moving average + yy = detrendma(yy,dL); + dflag = 'mean'; + + + max_L = min(300,n); #% maximum lag if L is undetermined + change_L = L is None + if change_L: + L = min(n-2, int(4./3*max_L+0.5)) + + + if method=='cov' or change_L: + R = self.tocovdata() + if change_L: + #finding where ACF is less than 2 st. deviations. + L = max_L-(np.abs(R.data[max_L::-1])>2*R.stdev[max_L::-1]).argmax() # a better L value + hasattr(windom, '__call__') + if wdef==1: # % modify L so that hanning and Parzen give appr. the same result + L = min(int(4*L/3),n-2) + print('The default L is set to %d' % L) + try: + win = window(2*L-1) + wname = window.__name__ + except: + wname = None + win = window + v = None + Be = None + + nf = rate*2**nextpow2(2*L-2) # Interpolate the spectrum with rate + nfft = 2*nf + +# S = createspec('freq',ftype); +# S.tr = g; +# S.note = ['dat2spec(',inputname(1),'), Method = ' method ]; +# S.norm = 0; % not normalized +# S.L = L; +# S.S = zeros(nf+1,m-1); + + if method=='psd': + S, f = psd(yy, Fs=1./dt, NFFT=nfft, detrend=detrend, window=window, + noverlap=noverlap,pad_to=pad_to, scale_by_freq=True) + + else :# cov method + # add a nugget effect to ensure that round off errors + # do not result in negative spectral estimates + + R.data[:L] = R.data[:L]*win[L::] + R.data[L:] = 0 + + spec = R.tospecdata(rate=2,nugget=nugget) + + if ( ~isempty(p) ), + alpha = (1-p); + #% Confidence interval constants + CI = [v/_invchi2( 1-alpha/2 ,v), v/_invchi2( alpha/2 ,v)]; + + + + + ind = find(Rper<0); + if any(ind) + Rper(ind) = 0; % set negative values to zero + warning('WAFO:DAT2SPEC','negative spectral estimates') + end + + + if wname=='parzen': + v = int(3.71*n/L) # degrees of freedom used in chi^2 distribution + Be = 2*pi*1.33/(L*dT) # % bandwidth (rad/sec) + elif wname=='hanning': + v = int(2.67*n/L); # degrees of freedom used in chi^2 distribution + Be = 2*pi/(L*dT); # % bandwidth (rad/sec) + elif wname=='bartlett': + v = int(3*n/L); # degrees of freedom used in chi^2 distribution + Be = 2*pi*1.33/(L*dT); # bandwidth (rad/sec) + if ftype=='w' + S.w = (0:nf)'/nf*pi/dT; % (rad/s) + S.S = real(Rper(1:(nf+1),1))*dT/pi; % (m^2*s/rad)one sided spectrum + S.Bw = Be; + else % ftype == f + S.f = (0:nf)'/nf/2/dT; % frequency Hz if dT is in seconds + S.S = 2*real(Rper(1:(nf+1),1))*dT; % (m^2*s) one sided spectrum + S.Bw = Be/(2*pi); % bandwidth in Hz + + + + + + def trdata(self, method='nonlinear', **options): ''' Estimate transformation, g, from data. diff --git a/pywafo/src/wafo/source/c_codes/setup.py b/pywafo/src/wafo/source/c_codes/setup.py new file mode 100644 index 0000000..966aded --- /dev/null +++ b/pywafo/src/wafo/source/c_codes/setup.py @@ -0,0 +1,16 @@ +''' +python setup.py build_src build_ext --inplace + +See also http://www.scipy.org/Cookbook/CompilingExtensionsOnWindowsWithMinGW +''' +# File setup.py +def configuration(parent_package='',top_path=None): + from numpy.distutils.misc_util import Configuration + config = Configuration('',parent_package,top_path) + + config.add_extension('c_library', + sources = ['c_library.pyf','c_functions.c']) + return config +if __name__ == "__main__": + from numpy.distutils.core import setup + setup(**configuration(top_path='').todict()) \ No newline at end of file diff --git a/pywafo/src/wafo/source/mvnprd/setup.py b/pywafo/src/wafo/source/mvnprd/setup.py new file mode 100644 index 0000000..6a993bb --- /dev/null +++ b/pywafo/src/wafo/source/mvnprd/setup.py @@ -0,0 +1,30 @@ +''' +python setup.py build_src build_ext --inplace + +See also http://www.scipy.org/Cookbook/CompilingExtensionsOnWindowsWithMinGW +''' + +# File setup.py +def compile_all(): + import os + files = ['mvnprd', 'mvnprodcorrprb'] + compile1_format = 'gfortran -fPIC -c %s.f' + for file in files: + os.system(compile1_format % file) + file_objects = ['%s.o' % file for file in files] + return file_objects + + +def configuration(parent_package='',top_path=None): + from numpy.distutils.misc_util import Configuration + libs = compile_all() + config = Configuration('',parent_package,top_path) + + config.add_extension('mvnprdmod', + libraries = libs, + sources = ['mvnprd_interface.f']) + return config +if __name__ == "__main__": + + from numpy.distutils.core import setup + setup(**configuration(top_path='').todict()) \ No newline at end of file diff --git a/pywafo/src/wafo/spectrum/__init__.py b/pywafo/src/wafo/spectrum/__init__.py index 091a9e0..13b56b1 100644 --- a/pywafo/src/wafo/spectrum/__init__.py +++ b/pywafo/src/wafo/spectrum/__init__.py @@ -5,6 +5,7 @@ Spectrum package in WAFO Toolbox. """ -from core import * #SpecData1D, SpecData2D, cltext +from core import * +#SpecData1D, SpecData2D, cltext import models import dispersion_relation diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index e5600b4..3483d3f 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -19,8 +19,8 @@ from pylab import stineman_interp from dispersion_relation import w2k #, k2w from wafo.wafodata import WafoData, now - from wafo.misc import sub_dict_select, nextpow2, discretize, JITImport, findpeaks #, tranproc +from wafo.graphutil import cltext try: from wafo.gaussian import Rind except ImportError: @@ -41,7 +41,7 @@ from wafo.plotbackend import plotbackend _WAFOCOV = JITImport('wafo.covariance') -__all__ = ['SpecData1D', 'SpecData2D', 'cltext', 'plotspec'] +__all__ = ['SpecData1D', 'SpecData2D', 'plotspec'] def _set_seed(iseed): '''Set seed of random generator''' @@ -235,7 +235,7 @@ def plotspec(specdata, linetype='b-', flag=1): ylbl4_txt = 'Wave Number Spectrum'; else: - raise ValueError('Frequency type unknown') + raise ValueError('Frequency type unknown') if hasattr(specdata, 'norm') and specdata.norm : @@ -607,101 +607,6 @@ def plotspec(specdata, linetype='b-', flag=1): # -def cltext(levels, percent=False, n=4, xs=0.036, ys=0.94, zs=0): - ''' - Places contour level text in the current window - - Parameters - ---------- - levels = vector of contour levels or the corresponding percent which the - contour line encloses - percent = 0 if levels are the actual contour levels (default) - 1 if levels are the corresponding percent which the - contour line encloses - n = maximum N digits of precision (default 4) - Returns - h = handles to the text objects. - CLTEXT creates text objects in the current figure and prints - "Level curves at:" if percent is False and - "Level curves enclosing:" otherwise - and the contour levels or percent. - - NOTE: - -The handles to the lines of text may also be found by - h = findobj(gcf,'gid','CLTEXT','type','text'); - h = findobj(gca,'gid','CLTEXT','type','text'); - -To make the text objects follow the data in the axes set the units - for the text objects 'data' by - set(h,'unit','data') - - Examples: - >>> from wafo.integrate import peaks - >>> import pylab as plt - >>> x,y,z = peaks(); - >>> h = plt.contour(x,y,z) - >>> h = cltext(h.levels) - >>> plt.show() - - data = rndray(1,2000,2); f = kdebin(data,{'kernel','epan','L2',.5,'inc',128}); - contour(f.x{:},f.f,f.cl),cltext(f.pl,1) - - See also - pdfplot - ''' - # TODO : Make it work like legend does (but without the box): include position options etc... - clevels = np.atleast_1d(levels) - _CLTEXT_TAG = 'CLTEXT' - cax = plotbackend.gca() - axpos = cax.get_position() - xint = axpos.intervalx - yint = axpos.intervaly - - xss = xint[0] + xs * (xint[1] - xint[0]) - yss = yint[0] + ys * (yint[1] - yint[0]) - - cf = plotbackend.gcf() # get current figure - #% delete cltext object if it exists - #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - def matchfun(x): - if hasattr(x, 'get_gid'): - return x.get_gid() == _CLTEXT_TAG - return False - h_cltxts = plotbackend.findobj(cf, matchfun); - if len(h_cltxts): - for i in h_cltxts: - try: - cax.texts.remove(i) - except: - warnings.warn('Tried to delete a non-existing CL-text') - - try: - cf.texts.remove(i) - except: - warnings.warn('Tried to delete a non-existing CL-text') - - charHeight = 1 / 33; - delta_y = charHeight; - - if percent: - titletxt = 'Level curves enclosing:'; - else: - titletxt = 'Level curves at:'; - - format = '%0.' + ('%d' % n) + 'g\n' - - cltxt = ''.join([format % level for level in clevels.tolist()]) - - titleProp = dict(gid=_CLTEXT_TAG, horizontalalignment='left', - verticalalignment='center', fontweight='bold', axes=cax) # - ha1 = plotbackend.figtext(xss, yss, titletxt, **titleProp) - - yss -= delta_y; - txtProp = dict(gid=_CLTEXT_TAG, horizontalalignment='left', - verticalalignment='top', axes=cax) - - ha2 = plotbackend.figtext(xss, yss, cltxt, **txtProp) - - return ha1, ha2 class SpecData1D(WafoData): @@ -754,8 +659,7 @@ class SpecData1D(WafoData): self.phi = 0.0 self.v = 0.0 self.norm = False - somekeys = ['angletype', 'phi', 'name', 'h', 'tr', 'freqtype', 'v', - 'type', 'norm'] + somekeys = ['phi', 'name', 'h', 'tr', 'freqtype', 'v','type', 'norm'] self.__dict__.update(sub_dict_select(kwds, somekeys)) diff --git a/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py b/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py new file mode 100644 index 0000000..3b08fcc --- /dev/null +++ b/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py @@ -0,0 +1,32 @@ +''' +Created on 19. juli 2010 + +@author: pab +''' + +from wafo.spectrum.dispersion_relation import w2k,k2w + +def test_k2w(): + ''' + >>> from numpy import arange + >>> k2w(arange(0.01,.5,0.2))[0] + array([ 0.3132092 , 1.43530485, 2.00551739]) + >>> k2w(arange(0.01,.5,0.2),h=20)[0] + array([ 0.13914927, 1.43498213, 2.00551724]) + ''' + + +def test_w2k(): + ''' + >>> w2k(range(4))[0] + array([ 0. , 0.1019368 , 0.4077472 , 0.91743119]) + >>> w2k(range(4),h=20)[0] + array([ 0. , 0.10503601, 0.40774726, 0.91743119]) + ''' + +def main(): + import doctest + doctest.testmod() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/pywafo/src/wafo/spectrum/test/test_models.py b/pywafo/src/wafo/spectrum/test/test_models.py new file mode 100644 index 0000000..44fac36 --- /dev/null +++ b/pywafo/src/wafo/spectrum/test/test_models.py @@ -0,0 +1,68 @@ +import numpy as np +from wafo.spectrum.models import (Bretschneider, Jonswap, OchiHubble, Tmaspec, + Torsethaugen, McCormick, Wallop) + +def test_bretschneider(): + ''' + >>> S = Bretschneider(Hm0=6.5,Tp=10) + >>> S((0,1,2,3)) + array([ 0. , 1.69350993, 0.06352698, 0.00844783]) + ''' +def test_jonswap(): + ''' + + >>> S = Jonswap(Hm0=7, Tp=11,gamma=1) + >>> S((0,1,2,3)) + array([ 0. , 1.42694133, 0.05051648, 0.00669692]) + >>> w = np.linspace(0,5) + >>> S2 = Bretschneider(Hm0=7, Tp=11) + + JONSWAP with gamma=1 should be equal to Bretscneider: + >>> np.all(np.abs(S(w)-S2(w))<1.e-7) + True + ''' + +def test_tmaspec(): + ''' + >>> S = Tmaspec(Hm0=7, Tp=11,gamma=1,h=10) + >>> S((0,1,2,3)) + array([ 0. , 0.70106233, 0.05022433, 0.00669692]) + ''' +def test_torsethaugen(): + ''' + >>> S = Torsethaugen(Hm0=7, Tp=11,gamma=1,h=10) + >>> S((0,1,2,3)) + array([ 0. , 1.19989709, 0.05819794, 0.0093541 ]) + >>> S.wind(range(4)) + array([ 0. , 1.13560528, 0.05529849, 0.00888989]) + >>> S.swell(range(4)) + array([ 0. , 0.0642918 , 0.00289946, 0.00046421]) + ''' + +def test_ochihubble(): + ''' + >>> S = OchiHubble(par=2) + >>> S(range(4)) + array([ 0. , 0.90155636, 0.04185445, 0.00583207]) + ''' +def test_mccormick(): + ''' + >>> S = McCormick(Hm0=6.5,Tp=10) + >>> S(range(4)) + array([ 0. , 1.87865908, 0.15050447, 0.02994663]) + ''' +def test_wallop(): + ''' + >>> S = Wallop(Hm0=6.5, Tp=10) + >>> S(range(4)) + array([ 0.00000000e+00, 9.36921871e-01, 2.76991078e-03, + 7.72996150e-05]) + ''' + + +def main(): + import doctest + doctest.testmod() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/pywafo/src/wafo/spectrum/test/test_specdata1d.py b/pywafo/src/wafo/spectrum/test/test_specdata1d.py new file mode 100644 index 0000000..13f9dec --- /dev/null +++ b/pywafo/src/wafo/spectrum/test/test_specdata1d.py @@ -0,0 +1,200 @@ +import wafo.spectrum.models +from wafo.spectrum import SpecData1D +def test_tocovmatrix(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap() + >>> S = Sj.tospecdata() + >>> acfmat = S.tocov_matrix(nr=3, nt=256, dt=0.1) + >>> acfmat[:2,:] + array([[ 3.06075987, 0. , -1.67750289, 0. ], + [ 3.05246132, -0.16662376, -1.66819445, 0.18634189]]) + ''' +def test_tocovdata(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap() + >>> S = Sj.tospecdata() + >>> Nt = len(S.data)-1 + >>> acf = S.tocovdata(nr=0, nt=Nt) + >>> acf.data[:5] + array([ 3.06093287, 2.23846752, 0.48630084, -1.1336035 , -2.03036854]) + ''' + +def test_to_t_pdf(): + ''' + The density of Tc is computed by: + >>> from wafo.spectrum import models as sm + >>> Sj = sm.Jonswap() + >>> S = Sj.tospecdata() + >>> f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100) + >>> ['%2.3f' % val for val in f.data[:10]] + ['0.000', '0.014', '0.027', '0.040', '0.050', '0.059', '0.067', '0.072', '0.077', '0.081'] + + estimated error bounds + >>> ['%2.4f' % val for val in f.err[:10]] + ['0.0000', '0.0003', '0.0003', '0.0004', '0.0006', '0.0009', '0.0016', '0.0019', '0.0020', '0.0021'] + ''' +def test_sim(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap();S = Sj.tospecdata() + >>> ns =100; dt = .2 + >>> x1 = S.sim(ns,dt=dt) + + >>> import numpy as np + >>> import scipy.stats as st + >>> x2 = S.sim(20000,20) + >>> truth1 = [0,np.sqrt(S.moment(1)[0]),0., 0.] + >>> funs = [np.mean,np.std,st.skew,st.kurtosis] + >>> for fun,trueval in zip(funs,truth1): + ... res = fun(x2[:,1::],axis=0) + ... m = res.mean() + ... sa = res.std() + ... #trueval, m, sa + ... np.abs(m-trueval)>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap();S = Sj.tospecdata() + >>> ns =100; dt = .2 + >>> x1 = S.sim_nl(ns,dt=dt) + + >>> import numpy as np + >>> import scipy.stats as st + >>> x2, x1 = S.sim_nl(ns=20000,cases=40) + >>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk') + >>> truth1[-1] = truth1[-1]-3 + >>> truth1 + [0, 1.7495200310090633, 0.18673120577479801, 0.061988521262417606] + + >>> funs = [np.mean,np.std,st.skew,st.kurtosis] + >>> for fun,trueval in zip(funs,truth1): + ... res = fun(x2[:,1::], axis=0) + ... m = res.mean() + ... sa = res.std() + ... #trueval, m, sa + ... np.abs(m-trueval)<2*sa + True + True + True + True + ''' +def test_stats_nl(): + ''' + >>> import wafo.spectrum.models as sm + >>> Hs = 7. + >>> Sj = sm.Jonswap(Hm0=Hs, Tp=11) + >>> S = Sj.tospecdata() + >>> me, va, sk, ku = S.stats_nl(moments='mvsk') + >>> me; va; sk; ku + 0.0 + 3.0608203389019537 + 0.18673120577479801 + 3.0619885212624176 + ''' +def test_testgaussian(): + ''' + >>> import wafo.spectrum.models as sm + >>> import wafo.transform.models as wtm + >>> import wafo.objects as wo + >>> Hs = 7 + >>> Sj = sm.Jonswap(Hm0=Hs) + >>> S0 = Sj.tospecdata() + >>> ns =100; dt = .2 + >>> x1 = S0.sim(ns, dt=dt) + + >>> S = S0.copy() + >>> me, va, sk, ku = S.stats_nl(moments='mvsk') + >>> S.tr = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4) + >>> ys = wo.mat2timeseries(S.sim(ns=2**13)) + >>> g0, gemp = ys.trdata() + >>> t0 = g0.dist2gauss() + >>> t1 = S0.testgaussian(ns=2**13, t0=t0, cases=50) + >>> sum(t1>t0)<5 + True + ''' +def test_moment(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.moment() + ([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt']) + ''' + +def test_nyquist_freq(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.nyquist_freq() + 3.0 + ''' +def test_sampling_period(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.sampling_period() + 1.0471975511965976 + ''' +def test_normalize(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.moment(2) + ([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt']) + + >>> Sn = S.copy(); Sn.normalize() + + Now the moments should be one + >>> Sn.moment(2) + ([1.0000000000000004, 0.99999999999999967], ['m0', 'm0tt']) + + ''' +def test_characteristic(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.characteristic(1) + (array([ 8.59007646]), array([[ 0.03040216]]), ['Tm01']) + + >>> [ch, R, txt] = S.characteristic([1,2,3]) # fact a vector of integers + >>> ch; R; txt + array([ 8.59007646, 8.03139757, 5.62484314]) + array([[ 0.03040216, 0.02834263, NaN], + [ 0.02834263, 0.0274645 , NaN], + [ NaN, NaN, 0.01500249]]) + ['Tm01', 'Tm02', 'Tm24'] + + >>> S.characteristic('Ss') # fact a string + (array([ 0.04963112]), array([[ 2.63624782e-06]]), ['Ss']) + + >>> S.characteristic(['Hm0','Tm02']) # fact a list of strings + (array([ 4.99833578, 8.03139757]), array([[ 0.05292989, 0.02511371], + [ 0.02511371, 0.0274645 ]]), ['Hm0', 'Tm02']) + ''' +def test_bandwidth(): + ''' + >>> import numpy as np + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=3) + >>> w = np.linspace(0,4,256) + >>> S = SpecData1D(Sj(w),w) #Make spectrum object from numerical values + >>> S.bandwidth([0,1,2,3]) + array([ 0.65354446, 0.3975428 , 0.75688813, 2.00207912]) + ''' +def main(): + import doctest + doctest.testmod() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/pywafo/src/wafo/test/test_objects.py b/pywafo/src/wafo/test/test_objects.py new file mode 100644 index 0000000..5175a44 --- /dev/null +++ b/pywafo/src/wafo/test/test_objects.py @@ -0,0 +1,59 @@ +# -*- coding:utf-8 -*- +""" +Created on 5. aug. 2010 + +@author: pab +""" +import wafo.data +import numpy as np + +def test_timeseries(): + ''' + >>> import wafo.data + >>> import wafo.objects as wo + >>> x = wafo.data.sea() + >>> ts = wo.mat2timeseries(x) + >>> ts.sampling_period() + 0.25 + + Estimate spectrum + >>> S = ts.tospecdata() + >>> S.data[:10] + array([ 0.01350817, 0.0050932 , 0.00182003, 0.00534806, 0.049407 , + 0.25144845, 0.28264622, 0.21398405, 0.18563258, 0.25878918]) + + Estimated covariance function + >>> rf = ts.tocovdata(lag=150) + >>> rf.data[:10] + array([ 0.22368637, 0.20838473, 0.17110733, 0.12237803, 0.07024054, + 0.02064859, -0.02218831, -0.0555993 , -0.07859847, -0.09166187]) + ''' +def test_timeseries_trdata(): + ''' + >>> import wafo.spectrum.models as sm + >>> import wafo.transform.models as tm + >>> from wafo.objects import mat2timeseries + >>> Hs = 7.0 + >>> Sj = sm.Jonswap(Hm0=Hs) + >>> S = Sj.tospecdata() #Make spectrum object from numerical values + >>> S.tr = tm.TrOchi(mean=0, skew=0.16, kurt=0, sigma=Hs/4, ysigma=Hs/4) + >>> xs = S.sim(ns=2**16) + >>> ts = mat2timeseries(xs) + >>> g0, gemp = ts.trdata(monitor=True) # Monitor the development + >>> g1, gemp = ts.trdata(method='m', gvar=0.5 ) # Equal weight on all points + >>> g2, gemp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5]) # Less weight on the ends + >>> S.tr.dist2gauss() + 5.9322684525265501 + >>> np.round(gemp.dist2gauss()) + 6.0 + >>> np.round(g0.dist2gauss()) + 4.0 + >>> np.round(g1.dist2gauss()) + 4.0 + >>> np.round(g2.dist2gauss()) + 4.0 + + ''' +if __name__=='__main__': + import doctest + doctest.testmod() \ No newline at end of file