From da96ab184a063cf9010a9d6b3f8cf22f12fbb1fe Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Mon, 29 Mar 2010 13:34:21 +0000 Subject: [PATCH] Added ppimport.py Small bugfix to estimation.py and distributions.py --- pywafo/src/wafo/objects.py | 4 +- pywafo/src/wafo/ppimport.py | 441 ++++++++++++++++++ pywafo/src/wafo/spectrum/models.py | 614 ++++++++++++------------- pywafo/src/wafo/stats/__init__.py | 2 +- pywafo/src/wafo/stats/distributions.py | 19 +- pywafo/src/wafo/stats/estimation.py | 51 +- 6 files changed, 790 insertions(+), 341 deletions(-) create mode 100644 pywafo/src/wafo/ppimport.py diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 4e5b7fa..6cac986 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -1,5 +1,5 @@ -#------------------------------------------------------------------------------- + # Name: module1 # Purpose: # @@ -8,7 +8,7 @@ # Created: 16.09.2008 # Copyright: (c) pab 2008 # Licence: -#------------------------------------------------------------------------------- + #!/usr/bin/env python diff --git a/pywafo/src/wafo/ppimport.py b/pywafo/src/wafo/ppimport.py new file mode 100644 index 0000000..b11fb55 --- /dev/null +++ b/pywafo/src/wafo/ppimport.py @@ -0,0 +1,441 @@ +#!/usr/bin/env python +""" +Postpone module import to future. + +Python versions: 1.5.2 - 2.3.x +Author: Pearu Peterson +Created: March 2003 +$Revision: 922 $ +$Date: 2004-11-27 14:23:27 -0700 (Sat, 27 Nov 2004) $ +""" +__all__ = ['ppimport','ppimport_attr','ppresolve'] + +import os +import sys +import types +import traceback + +DEBUG=0 + +_ppimport_is_enabled = 1 +def enable(): + """ Enable postponed importing.""" + global _ppimport_is_enabled + _ppimport_is_enabled = 1 + +def disable(): + """ Disable postponed importing.""" + global _ppimport_is_enabled + _ppimport_is_enabled = 0 + +class PPImportError(ImportError): + pass + +def _get_so_ext(_cache={}): + so_ext = _cache.get('so_ext') + if so_ext is None: + if sys.platform[:5]=='linux': + so_ext = '.so' + else: + try: + # if possible, avoid expensive get_config_vars call + from distutils.sysconfig import get_config_vars + so_ext = get_config_vars('SO')[0] or '' + except ImportError: + #XXX: implement hooks for .sl, .dll to fully support + # Python 1.5.x + so_ext = '.so' + _cache['so_ext'] = so_ext + return so_ext + +def _get_frame(level=0): + try: + return sys._getframe(level+1) + except AttributeError: + # Python<=2.0 support + frame = sys.exc_info()[2].tb_frame + for i in range(level+1): + frame = frame.f_back + return frame + +def ppimport_attr(module, name): + """ ppimport(module, name) is 'postponed' getattr(module, name) + """ + global _ppimport_is_enabled + if _ppimport_is_enabled and isinstance(module, _ModuleLoader): + return _AttrLoader(module, name) + return getattr(module, name) + +class _AttrLoader(object): + def __init__(self, module, name): + self.__dict__['_ppimport_attr_module'] = module + self.__dict__['_ppimport_attr_name'] = name + + def _ppimport_attr_getter(self): + module = self.__dict__['_ppimport_attr_module'] + if isinstance(module, _ModuleLoader): + # in case pp module was loaded by other means + module = sys.modules[module.__name__] + attr = getattr(module, + self.__dict__['_ppimport_attr_name']) + try: + d = attr.__dict__ + if d is not None: + self.__dict__ = d + except AttributeError: + pass + self.__dict__['_ppimport_attr'] = attr + + return attr + + def __nonzero__(self): + return 1 + + def __getattr__(self, name): + try: + attr = self.__dict__['_ppimport_attr'] + except KeyError: + attr = self._ppimport_attr_getter() + if name=='_ppimport_attr': + return attr + return getattr(attr, name) + + def __repr__(self): + if '_ppimport_attr' in self.__dict__: + return repr(self._ppimport_attr) + module = self.__dict__['_ppimport_attr_module'] + name = self.__dict__['_ppimport_attr_name'] + return "" % (`name`,`module`) + + __str__ = __repr__ + + # For function and class attributes. + def __call__(self, *args, **kwds): + return self._ppimport_attr(*args,**kwds) + + + +def _is_local_module(p_dir,name,suffices): + base = os.path.join(p_dir,name) + for suffix in suffices: + if os.path.isfile(base+suffix): + if p_dir: + return base+suffix + return name+suffix + +def ppimport(name): + """ ppimport(name) -> module or module wrapper + + If name has been imported before, return module. Otherwise + return ModuleLoader instance that transparently postpones + module import until the first attempt to access module name + attributes. + """ + global _ppimport_is_enabled + + level = 1 + parent_frame = p_frame = _get_frame(level) + while '__name__' not in p_frame.f_locals: + level = level + 1 + p_frame = _get_frame(level) + + p_name = p_frame.f_locals['__name__'] + if p_name=='__main__': + p_dir = '' + fullname = name + elif '__path__' in p_frame.f_locals: + # python package + p_path = p_frame.f_locals['__path__'] + p_dir = p_path[0] + fullname = p_name + '.' + name + else: + # python module + p_file = p_frame.f_locals['__file__'] + p_dir = os.path.dirname(p_file) + fullname = p_name + '.' + name + + # module may be imported already + module = sys.modules.get(fullname) + if module is not None: + if _ppimport_is_enabled or isinstance(module, types.ModuleType): + return module + return module._ppimport_importer() + + so_ext = _get_so_ext() + py_exts = ('.py','.pyc','.pyo') + so_exts = (so_ext,'module'+so_ext) + + for d,n,fn,e in [\ + # name is local python module or local extension module + (p_dir, name, fullname, py_exts+so_exts), + # name is local package + (os.path.join(p_dir, name), '__init__', fullname, py_exts), + # name is package in parent directory (scipy specific) + (os.path.join(os.path.dirname(p_dir), name), '__init__', name, py_exts), + ]: + location = _is_local_module(d, n, e) + if location is not None: + fullname = fn + break + + if location is None: + # name is to be looked in python sys.path. + fullname = name + location = 'sys.path' + + # Try once more if module is imported. + # This covers the case when importing from python module + module = sys.modules.get(fullname) + + if module is not None: + if _ppimport_is_enabled or isinstance(module,types.ModuleType): + return module + return module._ppimport_importer() + # It is OK if name does not exists. The ImportError is + # postponed until trying to use the module. + + loader = _ModuleLoader(fullname,location,p_frame=parent_frame) + if _ppimport_is_enabled: + return loader + + return loader._ppimport_importer() + +def _get_frame_code(frame): + filename = frame.f_code.co_filename + lineno = frame.f_lineno + result = '%s in %s:\n' % (filename,frame.f_code.co_name) + if not os.path.isfile(filename): + return result + f = open(filename) + i = 1 + line = f.readline() + while line: + line = f.readline() + i = i + 1 + if (abs(i-lineno)<2): + result += '#%d: %s\n' % (i,line.rstrip()) + if i>lineno+3: + break + f.close() + return result + +def frame_traceback(frame): + if not frame: + return + blocks = [] + f = frame + while f: + blocks.insert(0,_get_frame_code(f)) + f = f.f_back + print '='*50 + print '\n'.join(blocks) + print '='*50 + +class _ModuleLoader(object): + # Don't use it directly. Use ppimport instead. + + def __init__(self,name,location,p_frame=None): + + # set attributes, avoid calling __setattr__ + self.__dict__['__name__'] = name + self.__dict__['__file__'] = location + self.__dict__['_ppimport_p_frame'] = p_frame + + if location != 'sys.path': + from numpy.testing import Tester + self.__dict__['test'] = Tester(os.path.dirname(location)).test + + # install loader + sys.modules[name] = self + + def _ppimport_importer(self): + name = self.__name__ + + try: + module = sys.modules[name] + except KeyError: + raise ImportError,self.__dict__.get('_ppimport_exc_info')[1] + if module is not self: + exc_info = self.__dict__.get('_ppimport_exc_info') + if exc_info is not None: + raise PPImportError,\ + ''.join(traceback.format_exception(*exc_info)) + else: + assert module is self,`(module, self)` + + # uninstall loader + del sys.modules[name] + + if DEBUG: + print 'Executing postponed import for %s' %(name) + try: + module = __import__(name,None,None,['*']) + except Exception,msg: # ImportError: + if DEBUG: + p_frame = self.__dict__.get('_ppimport_p_frame',None) + frame_traceback(p_frame) + self.__dict__['_ppimport_exc_info'] = sys.exc_info() + raise + + assert isinstance(module,types.ModuleType),`module` + + #self.__dict__ = module.__dict__ + self.__dict__.update(module.__dict__) + self.__dict__['_ppimport_module'] = module + + # XXX: Should we check the existence of module.test? Warn? + from numpy.testing import Tester + test = Tester(os.path.dirname(module.__file__)).test + return module + + def __setattr__(self, name, value): + try: + module = self.__dict__['_ppimport_module'] + except KeyError: + module = self._ppimport_importer() + return setattr(module, name, value) + + def __getattr__(self, name): + try: + module = self.__dict__['_ppimport_module'] + except KeyError: + module = self._ppimport_importer() + return getattr(module, name) + + def __repr__(self): + global _ppimport_is_enabled + if not _ppimport_is_enabled: + try: + module = self.__dict__['_ppimport_module'] + except KeyError: + module = self._ppimport_importer() + return module.__repr__() + if '_ppimport_module' in self.__dict__: + status = 'imported' + elif '_ppimport_exc_info' in self.__dict__: + status = 'import error' + else: + status = 'import postponed' + return '' \ + % (`self.__name__`,`self.__file__`, status) + + __str__ = __repr__ + +def ppresolve(a,ignore_failure=None): + """ Return resolved object a. + + a can be module name, postponed module, postponed modules + attribute, string representing module attribute, or any + Python object. + """ + global _ppimport_is_enabled + if _ppimport_is_enabled: + disable() + a = ppresolve(a,ignore_failure=ignore_failure) + enable() + return a + if type(a) is type(''): + ns = a.split('.') + if ignore_failure: + try: + a = ppimport(ns[0]) + except: + return a + else: + a = ppimport(ns[0]) + b = [ns[0]] + del ns[0] + while ns: + if hasattr(a,'_ppimport_importer') or \ + hasattr(a,'_ppimport_module'): + a = getattr(a,'_ppimport_module',a) + if hasattr(a,'_ppimport_attr'): + a = a._ppimport_attr + b.append(ns[0]) + del ns[0] + if ignore_failure and not hasattr(a, b[-1]): + a = '.'.join(ns+b) + b = '.'.join(b) + if b in sys.modules and sys.modules[b] is None: + del sys.modules[b] + return a + a = getattr(a,b[-1]) + if hasattr(a,'_ppimport_importer') or \ + hasattr(a,'_ppimport_module'): + a = getattr(a,'_ppimport_module',a) + if hasattr(a,'_ppimport_attr'): + a = a._ppimport_attr + return a + +def _ppresolve_ignore_failure(a): + return ppresolve(a,ignore_failure=1) + +try: + import pydoc as _pydoc +except ImportError: + _pydoc = None + +if _pydoc is not None: + # Redefine __call__ method of help.__class__ to + # support ppimport. + import new as _new + + _old_pydoc_help_call = _pydoc.help.__class__.__call__ + def _ppimport_pydoc_help_call(self,*args,**kwds): + return _old_pydoc_help_call(self, *map(_ppresolve_ignore_failure,args), + **kwds) + _ppimport_pydoc_help_call.__doc__ = _old_pydoc_help_call.__doc__ + _pydoc.help.__class__.__call__ = _new.instancemethod(_ppimport_pydoc_help_call, + None, + _pydoc.help.__class__) + + _old_pydoc_Doc_document = _pydoc.Doc.document + def _ppimport_pydoc_Doc_document(self,*args,**kwds): + args = (_ppresolve_ignore_failure(args[0]),) + args[1:] + return _old_pydoc_Doc_document(self,*args,**kwds) + _ppimport_pydoc_Doc_document.__doc__ = _old_pydoc_Doc_document.__doc__ + _pydoc.Doc.document = _new.instancemethod(_ppimport_pydoc_Doc_document, + None, + _pydoc.Doc) + + _old_pydoc_describe = _pydoc.describe + def _ppimport_pydoc_describe(object): + return _old_pydoc_describe(_ppresolve_ignore_failure(object)) + _ppimport_pydoc_describe.__doc__ = _old_pydoc_describe.__doc__ + _pydoc.describe = _ppimport_pydoc_describe + +import inspect as _inspect +_old_inspect_getfile = _inspect.getfile +def _ppimport_inspect_getfile(object): + if isinstance(object,_ModuleLoader): + return object.__dict__['__file__'] + return _old_inspect_getfile(_ppresolve_ignore_failure(object)) +_ppimport_inspect_getfile.__doc__ = _old_inspect_getfile.__doc__ +_inspect.getfile = _ppimport_inspect_getfile + +_old_inspect_getdoc = _inspect.getdoc +def _ppimport_inspect_getdoc(object): + return _old_inspect_getdoc(_ppresolve_ignore_failure(object)) +_ppimport_inspect_getdoc.__doc__ = _old_inspect_getdoc.__doc__ +_inspect.getdoc = _ppimport_inspect_getdoc + +_old_inspect_getsource = _inspect.getsource +def _ppimport_inspect_getsource(object): + return _old_inspect_getsource(_ppresolve_ignore_failure(object)) +_ppimport_inspect_getsource.__doc__ = _old_inspect_getsource.__doc__ +_inspect.getsource = _ppimport_inspect_getsource + +import __builtin__ as _builtin +_old_builtin_dir = _builtin.dir +def _ppimport_builtin_dir(*arg): + if not arg: + p_frame = _get_frame(1) + g = p_frame.f_globals + l = p_frame.f_locals + l['_ppimport_old_builtin_dir'] = _old_builtin_dir + r = eval('_ppimport_old_builtin_dir()',g,l) + del r[r.index('_ppimport_old_builtin_dir')] + return r + return _old_builtin_dir(*map(_ppresolve_ignore_failure,arg)) +_ppimport_builtin_dir.__doc__ = _old_builtin_dir.__doc__ +_builtin.dir = _ppimport_builtin_dir diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py index 139c1ad..1e02867 100644 --- a/pywafo/src/wafo/spectrum/models.py +++ b/pywafo/src/wafo/spectrum/models.py @@ -56,14 +56,14 @@ from dispersion_relation import w2k #ppimport.enable() #_wafospectrum = ppimport.ppimport('wafo.spectrum') from core import SpecData1D -sech = lambda x: 1.0/cosh(x) +sech = lambda x: 1.0 / cosh(x) eps = finfo(float).eps -__all__ = ['Bretschneider','Jonswap','Torsethaugen','Wallop','McCormick','OchiHubble', - 'Tmaspec','jonswap_peakfact','jonswap_seastate','spreading', - 'w2k','k2w','phi1'] +__all__ = ['Bretschneider', 'Jonswap', 'Torsethaugen', 'Wallop', 'McCormick', 'OchiHubble', + 'Tmaspec', 'jonswap_peakfact', 'jonswap_seastate', 'spreading', + 'w2k', 'k2w', 'phi1'] # Model spectra @@ -121,17 +121,17 @@ def _gengamspec(wn, N=5, M=4): S = zeros_like(w) ##for w>0 # avoid division by zero - k = flatnonzero(w>0.0) - if k.size>0: - B = N/M - C = (N-1.0)/M + k = flatnonzero(w > 0.0) + if k.size > 0: + B = N / M + C = (N - 1.0) / M # # A = Normalizing factor related to Bretschneider form # A = B**C*M/gamma(C) # S(k) = A*wn(k)**(-N)*exp(-B*wn(k)**(-M)) logwn = log(w.take(k)) - logA = (C*log(B)+log(M)-sp.gammaln(C)) - S.put(k,exp(logA-N*logwn-B*exp(-M*logwn))) + logA = (C * log(B) + log(M) - sp.gammaln(C)) + S.put(k, exp(logA - N * logwn - B * exp(-M * logwn))) return S class ModelSpectrum(object): @@ -160,9 +160,9 @@ class ModelSpectrum(object): if w is None: if wc is None: - wc = 33./self.Tp - w = linspace(0,wc,nw) - S = SpecData1D(self.__call__(w),w) + wc = 33. / self.Tp + w = linspace(0, wc, nw) + S = SpecData1D(self.__call__(w), w) try: h = self.h S.h = h @@ -176,13 +176,13 @@ class ModelSpectrum(object): ''' Check if seastate is valid ''' - if self.Hm0<0: + if self.Hm0 < 0: raise ValueError('Hm0 can not be negative!') - if self.Tp<=0: + if self.Tp <= 0: raise ValueError('Tp must be positve!') - if self.Hm0==0.0: + if self.Hm0 == 0.0: warnings.warn('Hm0 is zero!') self._chk_extra_param() @@ -253,19 +253,19 @@ class Bretschneider(ModelSpectrum): self.chk_seastate() - def __call__(self,wi): + def __call__(self, wi): ''' Return Bretschnieder spectrum ''' w = atleast_1d(wi) - if self.Hm0>0: - wp = 2*pi/self.Tp - wn = w/wp - S = (self.Hm0/4.0)**2/wp * _gengamspec(wn,self.N,self.M) + if self.Hm0 > 0: + wp = 2 * pi / self.Tp + wn = w / wp + S = (self.Hm0 / 4.0) ** 2 / wp * _gengamspec(wn, self.N, self.M) else: S = zeros_like(w) return S -def jonswap_peakfact(Hm0,Tp): +def jonswap_peakfact(Hm0, Tp): ''' Jonswap peakedness factor, gamma, given Hm0 and Tp Parameters @@ -316,17 +316,17 @@ def jonswap_peakfact(Hm0,Tp): -------- jonswap ''' - Hm0,Tp = atleast_1d(Hm0,Tp) + Hm0, Tp = atleast_1d(Hm0, Tp) - x = Tp/sqrt(Hm0) + x = Tp / sqrt(Hm0) gam = ones_like(x) - k1 = flatnonzero(x<=5.14285714285714) - if k1.size>0: # #limiting gamma to [1 7] + k1 = flatnonzero(x <= 5.14285714285714) + if k1.size > 0: # #limiting gamma to [1 7] xk = x.take(k1) - D = 0.036-0.0056*xk # # approx 5.061*Hm0**2/Tp**4*(1-0.287*log(gam)) - gam.put(k1,minimum(exp(3.484*( 1.0-0.1975*D*xk**4.0 ) ),7.0)) # # gamma + D = 0.036 - 0.0056 * xk # # approx 5.061*Hm0**2/Tp**4*(1-0.287*log(gam)) + gam.put(k1, minimum(exp(3.484 * (1.0 - 0.1975 * D * xk ** 4.0)), 7.0)) # # gamma return gam def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'): @@ -401,35 +401,35 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'): # The following formulas are from Lewis and Allos 1990: - zeta = g*fetch/(u10**2) # dimensionless fetch, Table 1 + zeta = g * fetch / (u10 ** 2) # dimensionless fetch, Table 1 #zeta = min(zeta, 2.414655013429281e+004) if method.startswith('h'): - if method[-1]=='3': # Hasselman et.al (1973) - A = 0.076*zeta**(-0.22) - ny= 3.5*zeta**(-0.33) # dimensionless peakfrequency, Table 1 - epsilon1 = 9.91e-8*zeta**1.1 # dimensionless surface variance, Table 1 + if method[-1] == '3': # Hasselman et.al (1973) + A = 0.076 * zeta ** (-0.22) + ny = 3.5 * zeta ** (-0.33) # dimensionless peakfrequency, Table 1 + epsilon1 = 9.91e-8 * zeta ** 1.1 # dimensionless surface variance, Table 1 else: # Hasselman et.al (1976) - A = 0.0662*zeta**(-0.2) - ny = 2.84*zeta**(-0.3) # dimensionless peakfrequency, Table 1 - epsilon1 = 1.6e-7*zeta # dimensionless surface variance, Eq.4 + A = 0.0662 * zeta ** (-0.2) + ny = 2.84 * zeta ** (-0.3) # dimensionless peakfrequency, Table 1 + epsilon1 = 1.6e-7 * zeta # dimensionless surface variance, Eq.4 - sa = 0.07 - sb = 0.09 + sa = 0.07 + sb = 0.09 gam = 3.3 else: - A = 0.074*zeta**(-0.22) # Eq. 10 - ny = 3.57*zeta**(-0.33) # dimensionless peakfrequency, Eq. 11 - epsilon1 = 3.512e-4*A*ny**(-4.)*zeta**(-0.1) # dimensionless surface variance, Eq.12 - sa = 0.05468*ny**(-0.32) # Eq. 13 - sb = 0.078314*ny**(-0.16) # Eq. 14 - gam = maximum(17.54*zeta**(-0.28384),1) # Eq. 15 - - Tp = u10/(ny*g) # Table 1 - Hm0 = 4*sqrt(epsilon1)*u10**2./g # Table 1 - if output[0]=='l': - return Hm0,Tp,gam,sa,sb,A + A = 0.074 * zeta ** (-0.22) # Eq. 10 + ny = 3.57 * zeta ** (-0.33) # dimensionless peakfrequency, Eq. 11 + epsilon1 = 3.512e-4 * A * ny ** (-4.) * zeta ** (-0.1) # dimensionless surface variance, Eq.12 + sa = 0.05468 * ny ** (-0.32) # Eq. 13 + sb = 0.078314 * ny ** (-0.16) # Eq. 14 + gam = maximum(17.54 * zeta ** (-0.28384), 1) # Eq. 15 + + Tp = u10 / (ny * g) # Table 1 + Hm0 = 4 * sqrt(epsilon1) * u10 ** 2. / g # Table 1 + if output[0] == 'l': + return Hm0, Tp, gam, sa, sb, A else: - return dict(Hm0=Hm0,Tp=Tp,gamma=gam,sigmaA=sa,sigmaB=sb,Ag=A) + return dict(Hm0=Hm0, Tp=Tp, gamma=gam, sigmaA=sa, sigmaB=sb, Ag=A) class Jonswap(ModelSpectrum): ''' @@ -527,8 +527,8 @@ class Jonswap(ModelSpectrum): self.method = method self.wnc = wnc - if self.gamma==None or not isfinite(self.gamma) or self.gamma<1: - self.gamma = jonswap_peakfact(Hm0,Tp) + if self.gamma == None or not isfinite(self.gamma) or self.gamma < 1: + self.gamma = jonswap_peakfact(Hm0, Tp) self._preCalculateAg() @@ -539,7 +539,7 @@ class Jonswap(ModelSpectrum): Tp = self.Tp Hm0 = self.Hm0 gam = self.gamma - outsideJonswapRange = Tp>5*sqrt(Hm0) or Tp<3.6*sqrt(Hm0) + outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0) if outsideJonswapRange: txt0 = ''' Hm0,Tp is outside the JONSWAP range. @@ -547,7 +547,7 @@ class Jonswap(ModelSpectrum): ''' warnings.warn(txt0) - if gam<1 or 7 < gam: + if gam < 1 or 7 < gam: txt = ''' The peakedness factor, gamma, is possibly too large. The validity of the spectral density is questionable. @@ -555,30 +555,30 @@ class Jonswap(ModelSpectrum): warnings.warn(txt) - def _localspec(self,wn): + def _localspec(self, wn): Gf = self.peak_e_factor(wn) - return Gf*_gengamspec(wn,self.N, self.M) + return Gf * _gengamspec(wn, self.N, self.M) def _preCalculateAg(self): ''' PRECALCULATEAG Precalculate normalization. ''' - if self.gamma==1: + if self.gamma == 1: self.Ag = 1.0 self.method = 'parametric' elif self.Ag != None: self.method = 'custom' - if self.Ag<=0: + if self.Ag <= 0: raise ValueError('Ag must be larger than 0!') - elif self.method[0]=='i': + elif self.method[0] == 'i': # normalizing by integration self.method = 'integration' - if self.wnc<1.0: + if self.wnc < 1.0: raise ValueError('Normalized cutoff frequency, wnc, must be larger than one!') area1, unused_err1 = integrate.quad(self._localspec, 0, 1) area2, unused_err2 = integrate.quad(self._localspec, 1, self.wnc) area = area1 + area2 - self.Ag = 1.0/area - elif self.method[1]=='p': + self.Ag = 1.0 / area + elif self.method[1] == 'p': self.method = 'parametric' ## # Original normalization ## # NOTE: that Hm0**2/16 generally is not equal to intS(w)dw @@ -587,11 +587,11 @@ class Jonswap(ModelSpectrum): N = self.N M = self.M gammai = self.gamma - parametersOK = (3<=N and N<=50) or (2<=M and M <=9.5) and (1<= gammai and gammai<=20) + parametersOK = (3 <= N and N <= 50) or (2 <= M and M <= 9.5) and (1 <= gammai and gammai <= 20) if parametersOK: - f1NM = 4.1*(N-2*M**0.28+5.3)**(-1.45*M**0.1+0.96) - f2NM = (2.2*M**(-3.3) + 0.57)*N**(-0.58*M**0.37+0.53)-1.04*M**(-1.9)+0.94 - self.Ag = (1+ f1NM*log(gammai)**f2NM)/gammai + f1NM = 4.1 * (N - 2 * M ** 0.28 + 5.3) ** (-1.45 * M ** 0.1 + 0.96) + f2NM = (2.2 * M ** (-3.3) + 0.57) * N ** (-0.58 * M ** 0.37 + 0.53) - 1.04 * M ** (-1.9) + 0.94 + self.Ag = (1 + f1NM * log(gammai) ** f2NM) / gammai ### elseif N == 5 && M == 4, ### options.Ag = (1+1.0*log(gammai).**1.16)/gammai @@ -602,35 +602,35 @@ class Jonswap(ModelSpectrum): else: raise ValueError('Not knowing the normalization because N, M or peakedness parameter is out of bounds!') - if self.sigmaA!=0.07 or self.sigmaB!=0.09: + if self.sigmaA != 0.07 or self.sigmaB != 0.09: warnings.warn('Use integration to calculate Ag when sigmaA~=0.07 or sigmaB~=0.09') - def peak_e_factor(self,wn): + def peak_e_factor(self, wn): ''' PEAKENHANCEMENTFACTOR ''' - w = maximum(atleast_1d(wn),0.0) - sab = where(w>1,self.sigmaB,self.sigmaA) + w = maximum(atleast_1d(wn), 0.0) + sab = where(w > 1, self.sigmaB, self.sigmaA) - wnm12 = 0.5*((w-1.0)/sab)**2.0 - Gf = self.gamma**(exp(-wnm12)) + wnm12 = 0.5 * ((w - 1.0) / sab) ** 2.0 + Gf = self.gamma ** (exp(-wnm12)) return Gf - def __call__(self,wi): + def __call__(self, wi): ''' JONSWAP spectral density ''' w = atleast_1d(wi) - if (self.Hm0>0.0): + if (self.Hm0 > 0.0): - N = self.N - M = self.M - wp = 2*pi/self.Tp - wn = w/wp - Ag = self.Ag + N = self.N + M = self.M + wp = 2 * pi / self.Tp + wn = w / wp + Ag = self.Ag Hm0 = self.Hm0 - Gf = self.peak_e_factor(wn) - S = ((Hm0/4.0)**2/wp*Ag)*Gf*_gengamspec(wn,N,M) + Gf = self.peak_e_factor(wn) + S = ((Hm0 / 4.0) ** 2 / wp * Ag) * Gf * _gengamspec(wn, N, M) else: S = zeros_like(w) return S @@ -674,14 +674,14 @@ def phi1(wi, h, g=9.81): if h == inf: # % special case infinite water depth return ones_like(w) - k1 = w2k(w,0,inf,g=g)[0] - dw1 = 2.0*w/g # % dw/dk|h=inf - k2 = w2k(w,0,h,g=g)[0] + k1 = w2k(w, 0, inf, g=g)[0] + dw1 = 2.0 * w / g # % dw/dk|h=inf + k2 = w2k(w, 0, h, g=g)[0] - k2h = k2*h - dw2 = where(k1!=0,dw1/(tanh(k2h)+k2h/cosh(k2h)**2.0),0) # dw/dk|h=h0 + k2h = k2 * h + dw2 = where(k1 != 0, dw1 / (tanh(k2h) + k2h / cosh(k2h) ** 2.0), 0) # dw/dk|h=h0 - return where(k1!=0,(k1/k2)**3.0*dw2/dw1,0) + return where(k1 != 0, (k1 / k2) ** 3.0 * dw2 / dw1, 0) class Tmaspec(Jonswap): ''' JONSWAP spectrum for finite water depth @@ -755,24 +755,24 @@ class Tmaspec(Jonswap): Zeitschrift. ''' - def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09, - Ag=None, N=5, M=4, method='integration', wnc=6.0, + def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09, + Ag=None, N=5, M=4, method='integration', wnc=6.0, chk_seastate=True, h=42, g=9.81): self.g = g self.h = h - super(Tmaspec, self).__init__(Hm0,Tp,gamma,sigmaA,sigmaB,Ag,N,M,method,wnc,chk_seastate) + super(Tmaspec, self).__init__(Hm0, Tp, gamma, sigmaA, sigmaB, Ag, N, M, method, wnc, chk_seastate) self.type = 'TMA' - def phi(self,w,h=None,g=None): - if h==None: + def phi(self, w, h=None, g=None): + if h == None: h = self.h - if g==None: + if g == None: g = self.g - return phi1(w,h,g) + return phi1(w, h, g) - def __call__(self,w,h=None,g=None): + def __call__(self, w, h=None, g=None): jonswap = super(Tmaspec, self).__call__(w) - return jonswap*self.phi(w,h,g) + return jonswap * self.phi(w, h, g) class Torsethaugen(ModelSpectrum): ''' @@ -863,20 +863,20 @@ class Torsethaugen(ModelSpectrum): self._init_spec() - def __call__(self,w): + def __call__(self, w): ''' TORSETHAUGEN spectral density ''' - return self.wind(w)+self.swell(w) + return self.wind(w) + self.swell(w) def _chk_extra_param(self): Hm0 = self.Hm0 Tp = self.Tp - if Hm0>11 or Hm0>max((Tp/3.6)**2, (Tp-2)*12/11): + if Hm0 > 11 or Hm0 > max((Tp / 3.6) ** 2, (Tp - 2) * 12 / 11): txt0 = '''Hm0 is outside the valid range. The validity of the spectral density is questionable''' warnings.warn(txt0) - if Tp>20 or Tp<3: + if Tp > 20 or Tp < 3: txt1 = '''Tp is outside the valid range. The validity of the spectral density is questionable''' warnings.warn(txt1) @@ -916,33 +916,33 @@ class Torsethaugen(ModelSpectrum): # the parameters on the line below can be dependent on geographical location A10 = 0.7; A1 = 0.5; A20 = 0.6; A2 = 0.3; A3 = 6 - Tf = Af*(Hm0)**(1.0/3.0) - Tl = AL*sqrt(Hm0) # lower limit + Tf = Af * (Hm0) ** (1.0 / 3.0) + Tl = AL * sqrt(Hm0) # lower limit Tu = Au # upper limit #Non-dimensional scales # New call pab April 2005 - El = min(max((Tf-Tp)/(Tf-Tl),0),1) #wind sea - Eu = min(max((Tp-Tf)/(Tu-Tf),0),1) #Swell + El = min(max((Tf - Tp) / (Tf - Tl), 0), 1) #wind sea + Eu = min(max((Tp - Tf) / (Tu - Tf), 0), 1) #Swell - if Tp0: - Tpw = (16*S0*(1-exp(-Hm0/S1))*(0.4)**Nw/( G0w*Hpw**2))**(-1.0/(Nw-1.0)) + if Hpw > 0: + Tpw = (16 * S0 * (1 - exp(-Hm0 / S1)) * (0.4) ** Nw / (G0w * Hpw ** 2)) ** (-1.0 / (Nw - 1.0)) else: Tpw = inf @@ -989,13 +989,13 @@ class Torsethaugen(ModelSpectrum): if monitor: - if (3.6*sqrt(Hm0)<= Tp & Tp<=5*sqrt(Hm0)): + if (3.6 * sqrt(Hm0) <= Tp & Tp <= 5 * sqrt(Hm0)): print(' Jonswap range') print('Hm0 = %g' % Hm0) print('Ns, Ms = %g, %g Nw, Mw = %g, %g' % (Ns, Ms, Nw, Mw)) - print('gammas = %g gammaw = ' % (gammas,gammaw)) - print('Rps = %g Rpw = %g' % (Rps,Rpw)) + print('gammas = %g gammaw = ' % (gammas, gammaw)) + print('Rps = %g Rpw = %g' % (Rps, Rpw)) print('Hps = %g Hpw = %g' % (Hps, Hpw)) print('Tps = %g Tpw = %g' % (Tps, Tpw)) @@ -1003,10 +1003,10 @@ class Torsethaugen(ModelSpectrum): #G0s=Ms/((Ns/Ms)**(-(Ns-1)/Ms)*gamma((Ns-1)/Ms )) #normalizing factor # Wind part - self.wind = Jonswap(Hm0=Hpw, Tp=Tpw, gamma=gammaw, N=Nw, M=Mw, + self.wind = Jonswap(Hm0=Hpw, Tp=Tpw, gamma=gammaw, N=Nw, M=Mw, method=self.method, chk_seastate=False) # Swell part - self.swell = Jonswap(Hm0=Hps, Tp=Tps, gamma=gammas, N=Ns, M=Ms, + self.swell = Jonswap(Hm0=Hps, Tp=Tps, gamma=gammas, N=Ns, M=Ms, method=self.method, chk_seastate=False) class McCormick(Bretschneider): @@ -1060,23 +1060,23 @@ class McCormick(Bretschneider): self.type = 'McCormick' self.Hm0 = Hm0 self.Tp = Tp - if Tz==None: - Tz = 0.8143*Tp + if Tz == None: + Tz = 0.8143 * Tp self.Tz = Tz if chk_seastate: self.chk_seastate() - if M==None and self.Hm0>0: - self._TpdTz = Tp/Tz - M = 1.0/optimize.fminbound(self._localoptfun,0.01,5) + if M == None and self.Hm0 > 0: + self._TpdTz = Tp / Tz + M = 1.0 / optimize.fminbound(self._localoptfun, 0.01, 5) self.M = M - self.N = M+1.0 + self.N = M + 1.0 - def _localoptfun(self,x): + def _localoptfun(self, x): #LOCALOPTFUN Local function to optimize. - y = 1.0+x - return (y**(x)/sp.gamma(y)-self._TpdTz)**2.0 + y = 1.0 + x + return (y ** (x) / sp.gamma(y) - self._TpdTz) ** 2.0 class OchiHubble(ModelSpectrum): ''' OchiHubble bimodal spectral density model. @@ -1141,8 +1141,8 @@ class OchiHubble(ModelSpectrum): self.chk_seastate() self._init_spec() - def __call__(self,w): - return self.wind(w)+self.swell(w) + def __call__(self, w): + return self.wind(w) + self.swell(w) def _init_spec(self): @@ -1156,7 +1156,7 @@ class OchiHubble(ModelSpectrum): [0.65, 0.76], [0.90, 0.44], [0.77, 0.64], - [0.73,0.68], + [0.73, 0.68], [0.92, 0.39]]) wa = array([ [0.7, 1.15], [0.93, 1.5], @@ -1180,31 +1180,31 @@ class OchiHubble(ModelSpectrum): [0.039, 0.000], [0.046, 0.039], [0.046, 0.039]]) - Lpar = array([[3.00, 1.54,-0.062], - [3.00, 2.77,-0.112], - [2.55, 1.82,-0.089], - [2.65, 3.90,-0.085], - [2.60, 0.53,-0.069], - [1.35, 2.48,-0.102], - [4.95, 2.48,-0.102], - [1.80, 2.95,-0.105], - [4.50, 1.95,-0.082], - [6.40, 1.78,-0.069], - [0.70, 1.78,-0.069]]) + Lpar = array([[3.00, 1.54, -0.062], + [3.00, 2.77, -0.112], + [2.55, 1.82, -0.089], + [2.65, 3.90, -0.085], + [2.60, 0.53, -0.069], + [1.35, 2.48, -0.102], + [4.95, 2.48, -0.102], + [1.80, 2.95, -0.105], + [4.50, 1.95, -0.082], + [6.40, 1.78, -0.069], + [0.70, 1.78, -0.069]]) Hm0 = self.Hm0 Lpari = Lpar[self.par] - Li = hstack((Lpari[0],Lpari[1]*exp(Lpari[2]*Hm0))) + Li = hstack((Lpari[0], Lpari[1] * exp(Lpari[2] * Hm0))) - Hm0i = hp[self.par]*Hm0 - Tpi = 2*pi*exp(wb[self.par]*Hm0)/wa[self.par] - Ni = 4*Li+1 - Mi = [4, 4] + Hm0i = hp[self.par] * Hm0 + Tpi = 2 * pi * exp(wb[self.par] * Hm0) / wa[self.par] + Ni = 4 * Li + 1 + Mi = [4, 4] self.swell = Bretschneider(Hm0=Hm0i[0], Tp=Tpi[0], N=Ni[0], M=Mi[0]) - self.wind = Bretschneider(Hm0=Hm0i[1], Tp=Tpi[1], N=Ni[1], M=Mi[1]) + self.wind = Bretschneider(Hm0=Hm0i[1], Tp=Tpi[1], N=Ni[1], M=Mi[1]) def _chk_extra_param(self): - if self.par<0 or 10=1): + if any(X >= 1): raise ValueError('POISSON spreading: X value must be less than 1') return X - def fourier2a(self,r1): + def fourier2a(self, r1): ''' Returns the solution of R1 = sin(A)/A. ''' - A0 = flipud(linspace(0,pi+0.1,1025)) - funA = interp1d(sinc(A0/pi),A0) + A0 = flipud(linspace(0, pi + 0.1, 1025)) + funA = interp1d(sinc(A0 / pi), A0) A0 = funA(r1.ravel()) - A = asarray(A0) + A = asarray(A0) # Newton-Raphson da = ones_like(r1) @@ -1766,18 +1766,18 @@ class Spreading(object): ix = flatnonzero(A) for unused_iy in range(max_count): Ai = A[ix] - da[ix] = (sin(Ai) -Ai*r1[ix])/(cos(Ai)-r1[ix]) + da[ix] = (sin(Ai) - Ai * r1[ix]) / (cos(Ai) - r1[ix]) Ai = Ai - da[ix] # Make sure that the current guess is larger than zero and less than pi #x(ix) = xi + 0.1*(dx(ix) - 9*xi).*(xi<=0) + 0.38*(dx(ix)-6.2*xi +6.2).*(xi>pi) # Make sure that the current guess is larger than zero. - A[ix] = Ai + 0.5*(da[ix] - Ai)*(Ai<=0.0) + A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0) - ix = flatnonzero((abs(da) > sqrt(eps)*abs(A))*(abs(da) > sqrt(eps))) - if ix.size==0: - if any(A>pi): + ix = flatnonzero((abs(da) > sqrt(eps) * abs(A)) * (abs(da) > sqrt(eps))) + if ix.size == 0: + if any(A > pi): raise ValueError('BOX-CAR spreading: The A value must be less than pi') - return A.clip(min=1e-16,max=pi) + return A.clip(min=1e-16, max=pi) warnings.warn('Newton raphson method did not converge.') @@ -1787,52 +1787,52 @@ class Spreading(object): ''' Returns the solution of R1 = besseli(1,K)/besseli(0,K), ''' - K0 = hstack((linspace(0,10,513), linspace(10.00001,100))) - fun0 = lambda x : sp.ive(1,x)/sp.ive(0,x) - funK = interp1d(fun0(K0),K0) + K0 = hstack((linspace(0, 10, 513), linspace(10.00001, 100))) + fun0 = lambda x : sp.ive(1, x) / sp.ive(0, x) + funK = interp1d(fun0(K0), K0) K0 = funK(r1.ravel()) k1 = flatnonzero(isnan(K0)) - if (k1.size>0): + if (k1.size > 0): K0[k1] = 0.0 K0[k1] = K0.max() - ix0 = flatnonzero(r1!=0.0) + ix0 = flatnonzero(r1 != 0.0) K = zeros_like(r1) - fun = lambda x : fun0(x)-r1[ix] + fun = lambda x : fun0(x) - r1[ix] for ix in ix0: - K[ix] = optimize.fsolve(fun,K0[ix]) + K[ix] = optimize.fsolve(fun, K0[ix]) return K - def fourier2b(self,r1): + def fourier2b(self, r1): ''' Returns the solution of R1 = pi/(2*B*sinh(pi/(2*B)). ''' - B0 = hstack((linspace(eps,5,513), linspace(5.0001,100))) - funB = interp1d(self._r1ofsech2(B0),B0) + B0 = hstack((linspace(eps, 5, 513), linspace(5.0001, 100))) + funB = interp1d(self._r1ofsech2(B0), B0) B0 = funB(r1.ravel()) k1 = flatnonzero(isnan(B0)) - if (k1.size>0): + if (k1.size > 0): B0[k1] = 0.0 B0[k1] = max(B0) - ix0 = flatnonzero(r1!=0.0) + ix0 = flatnonzero(r1 != 0.0) B = zeros_like(r1) - fun = lambda x : 0.5*pi/(sinh(.5*pi/x))-x*r1[ix] + fun = lambda x : 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix] for ix in ix0: - B[ix] = abs(optimize.fsolve(fun,B0[ix])) + B[ix] = abs(optimize.fsolve(fun, B0[ix])) return B - def fourier2d(self,r1): + def fourier2d(self, r1): ''' Returns the solution of R1 = exp(-D**2/2). ''' - r = clip(r1,0.,1.0) - return where(r<=0,inf,sqrt(-2.0*log(r))) + r = clip(r1, 0., 1.0) + return where(r <= 0, inf, sqrt(-2.0 * log(r))) - def spread_par_s(self,wn): + def spread_par_s(self, wn): ''' Return spread parameter, S, of COS2S function Parameters @@ -1844,50 +1844,50 @@ class Spreading(object): S : ndarray spread parameter of COS2S functions ''' - if self.method==None: + if self.method == None: # no frequency dependent spreading, # but possible frequency dependent direction - s = atleast_1d(self.s_a) + s = atleast_1d(self.s_a) else: wn_lo = self.wn_lo wn_up = self.wn_up wn_c = self.wn_c - spa = self.s_a - spb = self.s_b - ma = self.m_a - mb = self.m_b + spa = self.s_a + spb = self.s_b + ma = self.m_a + mb = self.m_b # Mitsuyasu et. al and Hasselman et. al parametrization of # frequency dependent spreading - s = where(wn<=wn_c,spa*wn**ma,spb*wn**mb) - s[wn<=wn_lo] = 0.0 + s = where(wn <= wn_c, spa * wn ** ma, spb * wn ** mb) + s[wn <= wn_lo] = 0.0 - k = flatnonzero(wn_up0: - if self.method[0]=='d': + k = flatnonzero(wn_up < wn) + if k.size > 0: + if self.method[0] == 'd': # Donelan et. al. parametrization for B in SECH-2 - s[k] = spb*(wn_up)**mb + s[k] = spb * (wn_up) ** mb # Convert to S-paramater in COS-2S distribution r1 = self.r1ofsech2(s) - s = r1/(1.-r1) + s = r1 / (1. - r1) - elif self.method[0]=='b': + elif self.method[0] == 'b': # Banner parametrization for B in SECH-2 - s3m = spb*(wn_up)**mb - s3p = self._donelan(wn_up) - scale = s3m/s3p #% Scale so that parametrization will be continous - s[k] = scale*self.donelan(wn[k]) + s3m = spb * (wn_up) ** mb + s3p = self._donelan(wn_up) + scale = s3m / s3p #% Scale so that parametrization will be continous + s[k] = scale * self.donelan(wn[k]) r1 = self.r1ofsech2(s) #% Convert to S-paramater in COS-2S distribution - s = r1/(1.-r1) + s = r1 / (1. - r1) else: s[k] = 0.0 - if any(s<0): + if any(s < 0): raise ValueError('The COS2S spread parameter, S(w), value must be larger than 0') return s @@ -1895,16 +1895,16 @@ class Spreading(object): def _donelan(self, wn): ''' High frequency decay of B of sech2 paramater ''' - return 10.0**(-0.4+0.8393*exp(-0.567*log(wn**2))) + return 10.0 ** (-0.4 + 0.8393 * exp(-0.567 * log(wn ** 2))) def _r1ofsech2(self, B): ''' R1OFSECH2 Computes R1 = pi./(2*B.*sinh(pi./(2*B))) ''' realmax = finfo(float).max - tiny = 1./realmax - x = clip(2.*B,tiny,realmax) - xk = pi/x - return where(x<100.,xk/sinh(xk),-2.*xk/(exp(xk)*expm1(-2.*xk))) + tiny = 1. / realmax + x = clip(2. * B, tiny, realmax) + xk = pi / x + return where(x < 100., xk / sinh(xk), -2. * xk / (exp(xk) * expm1(-2. * xk))) @@ -1912,17 +1912,17 @@ def test_some_spectra(): S = Jonswap() w = arange(3.0) - S(w)*phi1(w,30.0) + S(w) * phi1(w, 30.0) S1 = S.tospecdata(w) S1.plot() import pylab as plb - w = plb.linspace(0,2.5) - S = Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11 - plb.plot(w,S(w)) - plb.plot(w,S(w,h=21)) - plb.plot(w,S(w,h=42)) + w = plb.linspace(0, 2.5) + S = Tmaspec(h=10, gamma=1) # Bretschneider spectrum Hm0=7, Tp=11 + plb.plot(w, S(w)) + plb.plot(w, S(w, h=21)) + plb.plot(w, S(w, h=42)) plb.show() plb.close('all') @@ -1930,39 +1930,39 @@ def test_some_spectra(): import pylab as plb #w = plb.linspace(0,3) - w,th = plb.ogrid[0:4,0:6] - k,k2 = w2k(w,th) - k1,k12 = w2k(w,th,h=20) - plb.plot(w,k,w,k2) + w, th = plb.ogrid[0:4, 0:6] + k, k2 = w2k(w, th) + k1, k12 = w2k(w, th, h=20) + plb.plot(w, k, w, k2) plb.show() plb.close('all') - w = plb.linspace(0,2,100) + w = plb.linspace(0, 2, 100) S = Torsethaugen(Hm0=6, Tp=8) - plb.plot(w,S(w),w,S.wind(w),w,S.swell(w)) + plb.plot(w, S(w), w, S.wind(w), w, S.swell(w)) - S1 = Jonswap(Hm0=7, Tp=11,gamma=1) - w = plb.linspace(0,2,100) - plb.plot(w,S1(w)) + S1 = Jonswap(Hm0=7, Tp=11, gamma=1) + w = plb.linspace(0, 2, 100) + plb.plot(w, S1(w)) plb.show() plb.close('all') - Hm0 = plb.arange(1,11) - Tp = plb.linspace(2,16) - T,H = plb.meshgrid(Tp,Hm0) - gam = jonswap_peakfact(H,T) - plb.plot(Tp,gam.T) + Hm0 = plb.arange(1, 11) + Tp = plb.linspace(2, 16) + T, H = plb.meshgrid(Tp, Hm0) + gam = jonswap_peakfact(H, T) + plb.plot(Tp, gam.T) plb.xlabel('Tp [s]') plb.ylabel('Peakedness parameter') - Hm0 = plb.linspace(1,20) + Hm0 = plb.linspace(1, 20) Tp = Hm0 - [T,H] = plb.meshgrid(Tp,Hm0) - gam = jonswap_peakfact(H,T) - v = plb.arange(0,8) - plb.contourf(Tp,Hm0,gam,v) + [T, H] = plb.meshgrid(Tp, Hm0) + gam = jonswap_peakfact(H, T) + v = plb.arange(0, 8) + plb.contourf(Tp, Hm0, gam, v) plb.colorbar() plb.show() plb.close('all') @@ -1970,19 +1970,19 @@ def test_some_spectra(): def test_spreading(): import pylab as plb pi = plb.pi - w = plb.linspace(0,3,257) - theta = plb.linspace(-pi,pi,129) - theta0 = lambda w: w*plb.pi/6.0 - D2 = Spreading('cos2s',theta0=theta0) - d1 =D2(theta,w)[0] + w = plb.linspace(0, 3, 257) + theta = plb.linspace(-pi, pi, 129) + theta0 = lambda w: w * plb.pi / 6.0 + D2 = Spreading('cos2s', theta0=theta0) + d1 = D2(theta, w)[0] t = plb.contour(d1.squeeze()) pi = plb.pi - D = Spreading('wrap_norm',s_a=10.0) + D = Spreading('wrap_norm', s_a=10.0) - w = plb.linspace(0,3,257) - theta = plb.linspace(-pi,pi,129) - d1 = D(theta,w) + w = plb.linspace(0, 3, 257) + theta = plb.linspace(-pi, pi, 129) + d1 = D(theta, w) plb.contour(d1[0]) plb.show() @@ -1994,4 +1994,4 @@ def main(): doctest.testmod() if __name__ == '__main__': - main() \ No newline at end of file + main() diff --git a/pywafo/src/wafo/stats/__init__.py b/pywafo/src/wafo/stats/__init__.py index 7588461..f08000f 100644 --- a/pywafo/src/wafo/stats/__init__.py +++ b/pywafo/src/wafo/stats/__init__.py @@ -6,7 +6,7 @@ Statistics package in WAFO Toolbox. """ from scipy.stats import * from wafo.stats.core import * -#from wafo.stats.distributions import * +from wafo.stats.distributions import * #from wafo.spectrum.core import SpecData1D #import wafo.spectrum.models #import wafo.spectrum.dispersion_relation diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py index a3b4b01..42f50a5 100644 --- a/pywafo/src/wafo/stats/distributions.py +++ b/pywafo/src/wafo/stats/distributions.py @@ -163,7 +163,7 @@ class rv_frozen(object): def __init__(self, dist, *args, **kwds): self.dist = dist loc0, scale0 = map(kwds.get, ['loc', 'scale']) - if isinstance(dist, rv_continuous): + if hasattr(dist, 'fix_loc_scale'): #isinstance(dist, rv_continuous): args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0) self.par = args + (loc0, scale0) else: # rv_discrete @@ -569,7 +569,7 @@ class rv_generic(object): and rv_continuous. """ - def fix_loc_scale(self, args, loc, scale=1): + def _fix_loc_scale(self, args, loc, scale=1): N = len(args) if N > self.numargs: if N == self.numargs + 1 and loc is None: @@ -585,8 +585,8 @@ class rv_generic(object): loc = 0.0 return args, loc, scale - def fix_loc(self, args, loc): - args, loc, scale = self.fix_loc_scale(args, loc) + def _fix_loc(self, args, loc): + args, loc, scale = self._fix_loc_scale(args, loc) return args, loc # These are actually called, and should not be overwritten if you @@ -617,7 +617,7 @@ class rv_generic(object): loc, scale, size, discrete = map(kwds.get, kwd_names, [None] * len(kwd_names)) - args, loc, scale = self.fix_loc_scale(args, loc, scale) + args, loc, scale = self._fix_loc_scale(args, loc, scale) cond = logical_and(self._argcheck(*args), (scale >= 0)) if not all(cond): raise ValueError, "Domain error in arguments." @@ -765,7 +765,8 @@ class rv_continuous(rv_generic): shapes=None, extradoc=None): rv_generic.__init__(self) - + self.fix_loc_scale = self._fix_loc_scale + if badvalue is None: badvalue = nan self.badvalue = badvalue @@ -4518,7 +4519,7 @@ class rv_discrete(rv_generic): shapes=None, extradoc=None): rv_generic.__init__(self) - + self.fix_loc = self._fix_loc if badvalue is None: badvalue = nan self.badvalue = badvalue @@ -4543,9 +4544,9 @@ class rv_discrete(rv_generic): self.pk = take(ravel(self.pk), indx, 0) self.a = self.xk[0] self.b = self.xk[-1] - self.P = make_dict(self.xk, self.pk) + self.P = dict(zip(self.xk, self.pk)) #make_dict(self.xk, self.pk) self.qvals = numpy.cumsum(self.pk, axis=0) - self.F = make_dict(self.xk, self.qvals) + self.F = dict(zip(self.xk, self.qvals)) #make_dict(self.xk, self.qvals) self.Finv = reverse_dict(self.F) self._ppf = new.instancemethod(sgf(_drv_ppf, otypes='d'), self, rv_discrete) diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index d7c6198..b75db37 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -1,26 +1,23 @@ -# Functions to implement several important functions for -# various Continous and Discrete Probability Distributions -# -# Author: Travis Oliphant 2002-2003 +''' +Contains FitDistribution and Profile class, which are +important classes for fitting to various Continous and Discrete Probability Distributions + +Author: Per A. Brodtkorb 2008 +''' from __future__ import division from wafo.plotbackend import plotbackend -from wafo.misc import ecross, findcross, JITImport +from wafo.misc import ecross, findcross from scipy.misc.ppimport import ppimport import numdifftools - +from scipy import special from scipy.linalg import pinv2 from scipy import optimize -# Trick to avoid error due to circular import -#_WAFOCOV = ppimport('wafo.covariance') -_WAFODIST = ppimport('wafo.stats.distributions') -#_WAFODIST = JITImport('wafo.stats.distributions') - from numpy import alltrue, arange, \ ravel, ones, sum, \ zeros, log, sqrt, exp @@ -43,6 +40,14 @@ arr = asarray all = alltrue +def chi2isf(p,df): + return special.chdtri(df, p) +def chi2sf(x, df): + return special.chdtrc(df, x) + +def norm_ppf(q): + return special.ndtri(q) + def valarray(shape, value=nan, typecode=None): """Return an array of all value. """ @@ -87,7 +92,7 @@ class rv_frozen(object): def __init__(self, dist, *args, **kwds): self.dist = dist loc0, scale0 = map(kwds.get, ['loc', 'scale']) - if isinstance(dist, _WAFODIST.rv_continuous): + if hasattr(dist, 'fix_loc_scale'): #isinstance(dist, _WAFODIST.rv_continuous): args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0) self.par = args + (loc0, scale0) else: # rv_discrete @@ -247,7 +252,7 @@ class Profile(object): self.i_free = nonzero(isfree) self.Lmax = Lmax - self.alpha_Lrange = 0.5 * _WAFODIST.chi2.isf(self.alpha, 1) + self.alpha_Lrange = 0.5 * chi2isf(self.alpha, 1) #_WAFODIST.chi2.isf(self.alpha, 1) self.alpha_cross_level = Lmax - self.alpha_Lrange lowLevel = self.alpha_cross_level - self.alpha_Lrange / 7.0 @@ -344,7 +349,8 @@ class Profile(object): pcov = self.fit_dist.par_cov[i_notfixed, :][:, i_notfixed] pvar = sum(numpy.dot(drl, pcov) * drl) - p_crit = _WAFODIST.norm.isf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5 + #p_crit = _WAFODIST.norm.isf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5 + p_crit = -norm_ppf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5 if self.pmin == None: self.pmin = p_opt - 5.0 * p_crit if self.pmax == None: @@ -392,7 +398,7 @@ class Profile(object): if alpha < self.alpha: raise ValueError('Unable to return CI with alpha less than %g' % self.alpha) - cross_level = self.Lmax - 0.5 * _WAFODIST.chi2.isf(alpha, 1) + cross_level = self.Lmax - 0.5 * chi2isf(alpha, 1) #_WAFODIST.chi2.isf(alpha, 1) ind = findcross(self.data, cross_level) N = len(ind) if N == 0: @@ -536,7 +542,7 @@ class FitDistribution(rv_frozen): self.par_cov[:, :] = nan pvar = numpy.diag(self.par_cov) - zcrit = -_WAFODIST.norm.ppf(self.alpha / 2.0) + zcrit = -norm_ppf(self.alpha / 2.0)#_WAFODIST.norm.ppf(self.alpha / 2.0) self.par_lower = self.par - zcrit * sqrt(pvar) self.par_upper = self.par + zcrit * sqrt(pvar) @@ -702,7 +708,7 @@ class FitDistribution(rv_frozen): Other distribution types will introduce deviations in the plot. ''' - bin, limits = numpy.histogram(self.data, normed=True, new=True) + bin, limits = numpy.histogram(self.data, normed=True) #, new=True) limits.shape = (-1, 1) xx = limits.repeat(3, axis=1) xx.shape = (-1,) @@ -793,7 +799,7 @@ class FitDistribution(rv_frozen): C1 = m - sqrt(0.5 * n * v) C2 = sqrt(v / (2.0 * n)) Tn = (T + 0.5 * k * isParUnKnown - C1) / C2 # chi2 with n degrees of freedom - pvalue = _WAFODIST.chi2.sf(Tn, n) + pvalue = chi2sf(Tn, n) #_WAFODIST.chi2.sf(Tn, n) return pvalue def nlogps(self, theta, x): @@ -964,6 +970,7 @@ class FitDistribution(rv_frozen): return - H def main(): + _WAFODIST = ppimport('wafo.stats.distributions') #nbinom(10, 0.75).rvs(3) import matplotlib matplotlib.interactive(True) @@ -986,15 +993,15 @@ def main(): Lx.get_CI(alpha=0.2) # CI for logSF=log(SF) - Lpr = phat.profile(i=1, logSF=log(R), link=phat.dist.link) + Lpr = phat.profile(i=0, logSF=log(R), link=phat.dist.link) Lpr.plot() Lpr.get_CI(alpha=0.075) - dlaplace.stats(0.8, loc=0) + _WAFODIST.dlaplace.stats(0.8, loc=0) # pass - t = planck(0.51000000000000001) + t = _WAFODIST.planck(0.51000000000000001) t.ppf(0.5) - t = zipf(2) + t = _WAFODIST.zipf(2) t.ppf(0.5) import pylab as plb _WAFODIST.rice.rvs(1)