From c93efca2669c5048e410b35d814df5ead9e1d725 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 9 Jul 2015 15:58:12 +0000 Subject: [PATCH] Replaced sub2index and index2sub with calls to np.ravel_multi_index and np.unravel_index, respectively. --- pywafo/.project | 12 +- pywafo/.pydevproject | 7 +- pywafo/build_all.py | 60 ---------- pywafo/setup.py | 165 +++++++++++++++++---------- pywafo/src/wafo/f2py_tools.py | 2 + pywafo/src/wafo/integrate.py | 203 +++++++++++++++++++++++++++++++++- pywafo/src/wafo/misc.py | 55 +-------- 7 files changed, 313 insertions(+), 191 deletions(-) delete mode 100644 pywafo/build_all.py diff --git a/pywafo/.project b/pywafo/.project index 731ae23..723695c 100644 --- a/pywafo/.project +++ b/pywafo/.project @@ -1,6 +1,6 @@ - pywafo + svn_pywafo @@ -10,16 +10,6 @@ - - org.eclipse.ui.externaltools.ExternalToolBuilder - auto,full,incremental, - - - LaunchConfigHandle - <project>/.externalToolBuilders/wafo_stats_tests.launch - - - org.python.pydev.pythonNature diff --git a/pywafo/.pydevproject b/pywafo/.pydevproject index 78946ac..d844cda 100644 --- a/pywafo/.pydevproject +++ b/pywafo/.pydevproject @@ -1,10 +1,5 @@ - - - - -/pywafo/src - + python 2.7 Default diff --git a/pywafo/build_all.py b/pywafo/build_all.py deleted file mode 100644 index 04f70f8..0000000 --- a/pywafo/build_all.py +++ /dev/null @@ -1,60 +0,0 @@ -""" -builds all extensions - -In notepad (or your favourite text editor) create a new text file and enter the following in the file: -[build] -compiler = mingw32 - -Save the file as "C:\Python25\Lib\distutils\distutils.cfg". -This will tell python to use the MinGW compiler when compiling extensions - - -See also http://www.scipy.org/Cookbook/CompilingExtensionsOnWindowsWithMinGW -""" -import os -import shutil - -def compile_all(): - wd = os.getcwd() - pkg_name = 'wafo' - root_dir = os.path.join(wd,'src',pkg_name) - root_src = os.path.join(root_dir, 'source') - buildscript = 'build_all.py' - - # if Linux - if os.name == 'posix': - build_call = 'python %s' % buildscript - lib_ext = '.so' - - # On Windows - elif os.name == 'nt': - build_call = 'python.exe %s' % buildscript - lib_ext = '.pyd' - - # give an Error for other OS-es - else: - raise UserWarning, \ - 'Untested platform:', os.name - - for root, dirs, files in os.walk(root_src): - dir1 = [dir for dir in dirs if not os.path.exists(os.path.join(root,dir,buildscript))] - for dir in dir1: - dirs.remove(dir) ## don't visit directories without buildscript - if buildscript in files: - print('Building: ', root) - #buildfile ='python.exe %s' % os.path.join(root,buildscript) - os.chdir(root) - t = os.system(build_call) - print(t) - - for file in os.listdir('.'): - if file.endswith(lib_ext): - dest_file = os.path.join(root_dir, file) - if os.path.exists(dest_file): - os.remove(dest_file) - shutil.copy(os.path.join(root,file), root_dir) - - os.chdir(wd) - -if __name__=='__main__': - compile_all() diff --git a/pywafo/setup.py b/pywafo/setup.py index e100976..f94bb86 100644 --- a/pywafo/setup.py +++ b/pywafo/setup.py @@ -16,38 +16,35 @@ python setup.py sdist bdist_wininst upload --show-response #!/usr/bin/env python import os +import shutil import sys import subprocess import re import warnings - +from Cython.Build import cythonize MAJOR = 0 MINOR = 1 MICRO = 2 ISRELEASED = False VERSION = '%d.%d.%d' % (MAJOR, MINOR, MICRO) - - -#sys.argv.append("develop") -#sys.argv.append("install") +# sys.argv.append("build_src") +# sys.argv.append("build_ext") +# sys.argv.append("--inplace") +# sys.argv.append("develop") +# sys.argv.append("install") DISTUTILS_DEBUG = True -pkg_name = 'wafo' -root_dir = os.path.join('src',pkg_name) +PKG_NAME = 'wafo' +ROOT_DIR = os.path.join('src',PKG_NAME) # make sure we import from this package, not an installed one: -sys.path.insert(0, root_dir) +sys.path.insert(0, ROOT_DIR) import info -#import wafo - -if True: #__file__ == 'setupegg.py': - # http://peak.telecommunity.com/DevCenter/setuptools - from setuptools import setup, Extension, find_packages -else: - from distutils.core import setup +from setuptools import find_packages # setup, Extension +from numpy.distutils.core import setup, Extension # as FExtension -# Return the svn version as a string, raise a ValueError otherwise def svn_version(): + '''Return the svn version as a string, raise a ValueError otherwise''' from numpy.compat import asstr env = os.environ.copy() @@ -84,59 +81,107 @@ short_version='%(version)s' version='%(version)s' release=%(isrelease)s """ - fid = open(os.path.join(root_dir,filename), 'w') + fid = open(os.path.join(ROOT_DIR,filename), 'w') try: fid.write(cnt % {'version': VERSION, 'isrelease': str(ISRELEASED)}) finally: fid.close() -if __name__=='__main__': - write_version_py() + +def get_library_extension(): + '''Return extension of an executable library''' + if os.name == 'posix': # executable library on Linux has extension .so + lib_ext = '.so' + elif os.name == 'nt': # extension on Windows is .pyd + lib_ext = '.pyd' + else: + raise UserWarning('Platform not supported:', os.name) + return lib_ext + + +def compile_all(): + wd = os.getcwd() + root_dir = os.path.join(wd,'src',PKG_NAME) + root_src = os.path.join(root_dir, 'source') + buildscript = 'build_all.py' + lib_ext = get_library_extension() + + if os.name == 'nt': # On Windows + build_call = 'python.exe %s' % buildscript + else: + build_call = 'python %s' % buildscript + for root, dirs, files in os.walk(root_src): + dir1 = [dir for dir in dirs + if not os.path.exists(os.path.join(root, dir, buildscript))] + for dir in dir1: + dirs.remove(dir) # don't visit directories without buildscript + if buildscript in files: + print('Building: ', root) + os.chdir(root) + t = os.system(build_call) + print(t) + + for file in os.listdir('.'): + if file.endswith(lib_ext): + dest_file = os.path.join(root_dir, file) + if os.path.exists(dest_file): + os.remove(dest_file) + shutil.copy(os.path.join(root, file), root_dir) + os.chdir(wd) + + +def setup_package(): + write_version_py() + join = os.path.join packages = find_packages('src') for p in packages: print(p) - package_paths =[p.replace(pkg_name+'.','').replace(pkg_name,'').replace('.',os.path.sep) for p in packages] - test_paths = [os.path.join(pkg_path,'test') for pkg_path in package_paths - if os.path.exists(os.path.join(root_dir,pkg_path,'test'))] - testscripts = [os.path.join(subtst, f) for subtst in test_paths - for f in os.listdir(os.path.join(root_dir, subtst)) + def convert_package2path(p): + return p.replace(PKG_NAME + '.', + '').replace(PKG_NAME, '').replace('.', os.path.sep) + package_paths = [convert_package2path(p) for p in packages] + test_paths = [join(pkg_path, 'test') for pkg_path in package_paths + if os.path.exists(join(ROOT_DIR, pkg_path, 'test'))] + testscripts = [join(subtst, f) for subtst in test_paths + for f in os.listdir(join(ROOT_DIR, subtst)) if not (f.startswith('.') or f.endswith('~') or f.endswith('.old') or f.endswith('.bak'))] - - #subpackages = ('spectrum','data','transform','covariance') - #subpackagesfull = [os.path.join(pkg_name,f) for f in subpackages] - - #subtests = [os.path.join(subpkg,'test') for subpkg in subpackages] - - #testscripts = [os.path.join(subtst, f) for subtst in subtests - # for f in os.listdir(os.path.join(root_dir, subtst)) - # if not (f.startswith('.') or f.endswith('~') or - # f.endswith('.old') or f.endswith('.bak'))] + datadir = 'data' - datafiles = [os.path.join(datadir, f) for f in os.listdir(os.path.join(root_dir, datadir)) + datafiles = [join(datadir, f) for f in os.listdir(join(ROOT_DIR, datadir)) if not (f.startswith('.') or f.endswith('~') or f.endswith('.old') or f.endswith('.bak') or f.endswith('.py') or f.endswith('test') )] + if 'build_ext' in sys.argv: + compile_all() + lib_ext = get_library_extension() + libs = [f for f in os.listdir(join(ROOT_DIR)) if f.endswith(lib_ext)] - # executable library on Linux has extension .so - if os.name == 'posix': - lib_ext = '.so' - - # extension on Windows is .pyd - elif os.name == 'nt': - lib_ext = '.pyd' - - # give an Error for other OS-es - else: - raise UserWarning, \ - 'Platform not supported:', os.name - - libs = [f for f in os.listdir(os.path.join(root_dir)) if f.endswith(lib_ext) ] -# libs = [f for f in os.listdir(os.path.join(root_dir)) if f.endswith('.pyd') ] - - packagedata = testscripts + datafiles + libs #['c_library.pyd'] #,'disufq1.c','diffsumfunq.pyd','diffsumfunq.pyf','findrfc.c','rfc.pyd','rfc.pyf'] - + packagedata = testscripts + datafiles + libs + +# ext_module_list = cythonize(join(ROOT_DIR, "primes.pyx")) +# +# for ext_module in ext_module_list: +# if not isinstance(ext_module, Extension): +# ext_module.__class__ = Extension + +# for name, src_files in [('mvn',('mvn.pyf', 'mvndst.f')), +# ('c_library',('c_library.pyf', 'c_functions.c'))]: +# sources = [join(ROOT_DIR, 'source', name, f) for f in src_files] +# ext_module_list.append(Extension(name='%s.%s' % (PKG_NAME, name), +# sources=sources)) + +# sources = [join(ROOT_DIR, 'source', 'mreg', 'cov2mmpdfreg_intfc.f'), ] +# libs = [join(ROOT_DIR, 'source', 'mreg', f) +# for f in ['dsvdc', 'mregmodule', 'intfcmod'] ] +# ext_module_list.append(Extension(name='wafo.covmod', sources=sources, +# libraries=libs)) + +# mvn_sources = [join(root_mvn, 'source', 'mvn', 'mvn.pyf'), +# join(root_mvn, 'source', 'mvn','mvndst.f')] +# ext_module_list.append(Extension(name='wafo.mvn', sources=mvn_sources)) + setup( version = VERSION, author='WAFO-group', @@ -146,10 +191,11 @@ if __name__=='__main__': install_requires = ['numpy>=1.4','numdifftools>=0.2'], license = "GPL", url='http://code.google.com/p/pywafo/', - name = pkg_name, + name = PKG_NAME, package_dir = {'': 'src'}, packages = packages, package_data = {'': packagedata}, + # ext_modules = ext_module_list, classifiers=[ 'Development Status :: 4 - Beta', 'Intended Audience :: Education', @@ -159,11 +205,8 @@ if __name__=='__main__': 'Programming Language :: Python :: 2.6', 'Topic :: Scientific/Engineering :: Mathematics', ], - #packages = [package_name,] + list(subpackagesfull), - #package_data = {package_name: packagedata}, - #package_data = {'': ['wafo.cfg']}, - #scripts = [os.path.join('bin', f) - # for f in os.listdir('bin') - # if not (f.startswith('.') or f.endswith('~') or - # f.endswith('.old') or f.endswith('.bak'))], ) + + +if __name__=='__main__': + setup_package() \ No newline at end of file diff --git a/pywafo/src/wafo/f2py_tools.py b/pywafo/src/wafo/f2py_tools.py index 45bcb5a..77ce15e 100644 --- a/pywafo/src/wafo/f2py_tools.py +++ b/pywafo/src/wafo/f2py_tools.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- """ f2py c_library.pyf c_functions.c -c @@ -61,3 +62,4 @@ def f2py_call_str(): except NameError: raise UserWarning('Couldn\'t locate f2py. ' 'Should be part of NumPy installation.') + diff --git a/pywafo/src/wafo/integrate.py b/pywafo/src/wafo/integrate.py index 1e5b090..5669a45 100644 --- a/pywafo/src/wafo/integrate.py +++ b/pywafo/src/wafo/integrate.py @@ -8,6 +8,7 @@ from scipy import special as sp from wafo.plotbackend import plotbackend as plt from scipy.integrate import simps, trapz from wafo.demos import humps +from pychebfun import Chebfun _EPS = np.finfo(float).eps _POINTS_AND_WEIGHTS = {} @@ -1423,8 +1424,208 @@ def test_docstrings(): doctest.testmod() +def levin_integrate(): + ''' An oscillatory integral + Sheehan Olver, December 2010 + + + (Chebfun example quad/LevinIntegrate.m) + + This example computes the highly oscillatory integral of + + f * exp( 1i * w * g ), + + over (0,1) using the Levin method [1]. This method computes the integral + by rewriting it as an ODE + + u' + 1i * w * g' u = f, + + so that the indefinite integral of f * exp( 1i * w * g ) is + + u * exp( 1i * w * g ). + + + + We use as an example + + f = 1 / ( x + 2 ); + g = cos( x - 2 ); + w = 100000; + + # + References: + + [1] Levin, D., Procedures for computing one and two-dimensional integrals + of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538 + ''' + exp = np.exp + domain=[0, 1] + x = Chebfun.identity(domain=domain) + f = 1./(x+2) + g = np.cos(x-2) + D = np.diff(domain) + + + # Here is are plots of this integrand, with w = 100, in complex space + w = 100; + line_opts = dict(line_width=1.6) + font_opts = dict(font_size= 14) + # + + intg = f*exp(1j*w*g) + xs, ys, xi, yi, d = intg.plot_data(1000) + #intg.plot(with_interpolation_points=True) + #xi = np.linspace(0, 1, 1024) +# plt.plot(xs, ys) # , **line_opts) +# plt.plot(xi, yi, 'r.') +# #axis equal +# plt.title('Complex plot of integrand') #,**font_opts) +# plt.show('hold') + ## + # and of just the real part +# intgr = np.real(intg) +# xs, ys, xi, yi, d = intgr.plot_data(1000) + #intgr.plot() +# plt.plot(xs, np.real(ys)) # , **line_opts) +# plt.plot(xi, np.real(yi), 'r.') + #axis equal +# plt.title('Real part of integrand') #,**font_opts) +# plt.show('hold') + + ## + # The Levin method will be accurate for large and small w, and the time + # taken is independent of w. Here we take a reasonably large value of w. + w = 1000; + intg = f*exp(1j*w*g) + val0 = np.sum(intg) + # val1 = sum(intg) + print(val0) + ## + # Start timing + #tic + + ## + # Construct the operator L + L = D + 1j*w*np.diag(g.differentiate()) + + ## + # From asymptotic analysis, we know that there exists a solution to the + # equation which is non-oscillatory, though we do not know what initial + # condition it satisfies. Thus we find a particular solution to this + # equation with no boundary conditions. + + u = L / f + + ## + # Because L is a differential operator with derivative order 1, \ expects + # it to be given a boundary condition, which is why the warning message is + # displayed. However, this doesn't cause any problems: though there are, + # in fact, a family of solutions to the ODE without boundary conditions + # due to the kernel + # + # exp(- 1i * w * g), + # + # it does not actually matter which particular solution is computed. + # Non-uniqueness is also not an issue: \ in matlab is least squares, hence + # does not require uniqueness. The existence of a non-oscillatory solution + # ensures that \ converges to a u with length independent of w. + # + # One could prevent the warning by applying a boundary condition consistent + # with the rest of the system, that is + # L.lbc = {L(1,:),f(0)}; + + ## + # Now we evaluate the antiderivative at the endpoints to obtain the + # integral. + + u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0)) + + #toc + + + ## + # Here is a way to compute the integral using Clenshaw--Curtis quadrature. + # As w becomes large, this takes an increasingly long time as the + # oscillations must be resolved. + + #tic + sum( f*exp(1j*w*g) ) + #toc + + +# aLevinTQ[omega_,a_,b_,f_,g_,nu_,wprec_,prm_,test_,basis_,np_]:= +# +# Module[{r,betam,A,AA,BB,S,F,w=N[omega, wprec]}, +# M=Length[nu]-1; +# PB[k_,t_]:=If[basis==1,t^k,ChebyshevT[k,t]]; +# +# ff[t_]:=((b-a)/2)*f[(b-a)*t/2+(a+b)/2]; +# +# gg[t_]:=g[(b-a)*t/2+(a+b)/2]; +# dgg[t_]:=Derivative[1][gg][t]; +# +# If[test==0, betam=Min[Abs[dgg[-1]*w], Abs[dgg[1]*w]]; +# While[prm*M/betam >=1, betam=2*betam]]; +# If[test>0,x[k_]:=N[Cos[k*Pi/M], wprec],x[k_]:= +# Which[kprm*M, N[1-(M-k)/betam, wprec]]]; +# +# Psi[k_,t_]:=Derivative[0,1][PB][k,t]+I*w*dgg[t]*PB[k,t]; +# +# ne[j_]:=nu[[j+1]]; S[-1]=0; S[j_]:=Sum[ne[i],{i,0,j}]; +# nn=S[M]-1; +# A=ConstantArray[0,{nn+1,nn+1}]; +# F=ConstantArray[0,nn+1]; r=0; +# While[rx[r]], +# {k,0,S[M]-1}],{j,S[r-1],S[r]-1}]; +# +# Do[BB[j]=Limit[Derivative[Mod[j-S[r-1],ne[r]]][ff][t], +# t->x[r]],{j,S[r-1],S[r]-1}]; +# Do[F[[j]]=BB[j-1],{j,S[r-1]+1,S[r]}]; +# Do[Do[A[[j,k]]=AA[j-1,k-1],{k,1,S[M]}],{j,S[r-1]+1,S[r]}]; +# r=r+1;]; (*sv=SingularValueList[N[A,wprec]]; +# con=sv[[1]]/sv[[-1]]; Print["cond2(A)= ",N[con,3]];*) +# LS=Block[{MaxExtraPrecision=0},LeastSquares[N[A, wprec],F]]; +# vvv[t_]:=Sum[LS[[k+1]]*PB[k,t], {k,0,nn}]; +# NR=vvv[1]*Exp[I*w*gg[1]]-vvv[-1]*Exp[I*w*gg[-1]]; +# Print["n=", np+ii+2s-2, ", Result= ", N[NR, wprec/2+5]]; +# If[ii==0,PR=NR];]; +# (* End of subroutine aLevinTQ /A.I. Hascelik, July 2013/ *) +# +# def main_levin(): +# a=1; b=2; +# omega=100; +# prm=1/2; +# f[t_]=Exp[4t] +# g[t_]=t+Exp[4t]*Gamma[t] +# +# dg[t_]:=Derivative[1][g][t]; +# +# prec=16 +# wprec=2*prec +# delta = min(abs(omega*dg(a)), abs(omega*dg(b))) +# alpha = min(abs(omega*g(a)), abs(omega*g(b))) +# s=1; #(*if s>1, the integral is computed by Q_s^L*) +# test= 1 if delta<10 or alpha <=10 or s>1 else 0 +# +# m = 1 if s>1 else np.floor(prec/max(np.log10(beta+1),1)+2) +# nc = 2*m+1 #(*or np=2m, number of collocation points*) +# basis=1; # (*take basis=0 for the Chebysev basis*) +# for ii in range(0, 2, 2): +# nu = np.ones((nc+ii,)) # ConstantArray[1,nc+ii]; +# nu[0] = s +# nu[-1] = s +# #nu[[1]]=s; +# #nu[[-1]]=s; +# aLevinTQ[omega,a,b,f,g,nu,wprec,prm,test,basis,nc], +# #{ii,0,2,2}]; +# Print["Error= ",Abs[NR-PR]]; + + if __name__ == '__main__': - test_docstrings() + levin_integrate() + # test_docstrings() # qdemo(np.exp, 0, 3, plot_error=True) # plt.show('hold') # main() diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 6fa00cb..c6f1be3 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -240,29 +240,7 @@ def args_flat(*args): raise ValueError('Number of arguments must be 1 or 3!') -def _check_and_adjust_shape(shape, nsub=None): - s = np.atleast_1d(shape) - ndim = len(s) - if ndim < 1: - raise ValueError('Shape vector must have at least 1 element.') - ndim = len(s) - if nsub is None: - nsub = ndim - if ndim <= nsub: # add trailing singleton dimensions - s = np.hstack([s, np.ones(nsub - ndim, dtype=int)]) - else: # Adjust for linear indexing on last element - s = np.hstack([s[:nsub - 1], np.prod(s[nsub - 1:])]) - return s - - -def _sub2index_factor(shape, order='C'): - '''Return multiplier needed for calculating linear index from subscripts. - ''' - step = 1 if order == 'F' else -1 # C order - return np.hstack([1, np.cumprod(shape[::step][:-1])])[::step] - - -def index2sub(shape, index, nsub=None, order='C'): +def index2sub(shape, index, order='C'): ''' Returns Multiple subscripts from linear index. @@ -272,8 +250,6 @@ def index2sub(shape, index, nsub=None, order='C'): shape of array index : linear index into array - nsub : int optional - Number of subscripts returned. default nsub=len(shape) order : {'C','F'}, optional The order of the linear index. 'C' means C (row-major) order. @@ -300,18 +276,7 @@ def index2sub(shape, index, nsub=None, order='C'): -------- sub2index ''' - ndx = np.atleast_1d(index) - s = _check_and_adjust_shape(shape, nsub) - k = _sub2index_factor(s, order) - n = len(s) - step = -1 if order == 'F' else 1 # C order - subscripts = [0, ] * n - for i in range(n)[::step]: - vi = np.remainder(ndx, k[i]) - subscript = np.array((ndx - vi) / k[i], dtype=int) - subscripts[i] = subscript - ndx = vi - return tuple(subscripts) + return np.unravel_index(index, shape, order=order) def sub2index(shape, *subscripts, **kwds): @@ -350,21 +315,7 @@ def sub2index(shape, *subscripts, **kwds): -------- index2sub ''' - nsub = len(subscripts) - s = _check_and_adjust_shape(shape, nsub) - k = _sub2index_factor(s, **kwds) - - ndx = 0 - s0 = np.shape(subscripts[0]) - for i, subscript in enumerate(subscripts): - np.testing.assert_equal( - s0, np.shape(subscript), - 'The subscripts vectors must all be of the same shape.') - if (np.any(subscript < 0)) or (np.any(s[i] <= subscript)): - raise IndexError('Out of range subscript.') - - ndx = ndx + k[i] * subscript - return ndx + return np.ravel_multi_index(subscripts, shape, **kwds) def is_numlike(obj):