Replaced

sub2index and index2sub with calls to np.ravel_multi_index and np.unravel_index, respectively.

pywafo/.project

@@ -1,6 +1,6 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <projectDescription>
-	<name>pywafo</name>
+	<name>svn_pywafo</name>
 	<comment></comment>
 	<projects>
 	</projects>
@@ -10,16 +10,6 @@
 			<arguments>
 			</arguments>
 		</buildCommand>
-		<buildCommand>
-			<name>org.eclipse.ui.externaltools.ExternalToolBuilder</name>
-			<triggers>auto,full,incremental,</triggers>
-			<arguments>
-				<dictionary>
-					<key>LaunchConfigHandle</key>
-					<value>&lt;project&gt;/.externalToolBuilders/wafo_stats_tests.launch</value>
-				</dictionary>
-			</arguments>
-		</buildCommand>
 	</buildSpec>
 	<natures>
 		<nature>org.python.pydev.pythonNature</nature>

pywafo/.pydevproject

@@ -1,10 +1,5 @@
 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
-<?eclipse-pydev version="1.0"?>
-
-<pydev_project>
-<pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
-<path>/pywafo/src</path>
-</pydev_pathproperty>
+<?eclipse-pydev version="1.0"?><pydev_project>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.7</pydev_property>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>
 </pydev_project>

pywafo/build_all.py

@@ -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()

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()

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.')
+

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[k<prm*M, N[-1+k/betam, wprec], k==Ceiling[prm*M],0,
+# k>prm*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[r<M+1, Do[Do[ AA[j,k]=
+# Limit[Derivative[0,Mod[j-S[r-1],ne[r]]][Psi][k,t],t->x[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()

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):