Mailman 3 July 2012 - NumPy-Discussion

swig + numpy + variable length arrays
by Andrew Nelson July 16, 2012

July 16, 2012

Dear list, I am trying to SWIG a C function with the signature: void abelescalcall(long numcoefs, double *coefP, long npoints , double *yP, double *xP); numcoefs corresponds to the number of points in the coefP array. npoints corresponds to the number of points in yP and xP arrays. coefP and xP are inputs (IN_ARRAY1), yP is an output (ARGOUT_ARRAY1). I have been trying to use the following input file: %module reflect %{ #define SWIG_FILE_WITH_INIT #include "myfitfunctions.h" %} %include "numpy.i" %init %{ import_array(); %} %apply (long DIM1, double* IN_ARRAY1){(long numcoefs, double* coefP)} %apply (long DIM1, double* ARGOUT_ARRAY1){(long len2, double* yP)} %apply (long DIM1, double* IN_ARRAY1){(long len3, double *xP)} %include "myfitfunctions.h" %rename (abelescalcall) my_abeles; %inline %{ my_abeles(long numcoefs, double* coefP, long len2, double* yP, long len3, double *xP) { if (len2 != len3) { PyErr_Format(PyExc_ValueError, "Arrays of lengths (%d,%d) given", len2, len3); return 0.0; } abelescalcall(numcoefs, coefP, len2, yP, xP); } %} However, if I look at the wrapper code, and if I try to call my_abeles, then I am asked for 6 parameters. I only want to supply two numpy arrays - coefP and xP. thanks for any help you are able to give. (I spent the whole afternoon trying to get this bloomin thing to work). cheers, Andrew. -- _____________________________________ Dr. Andrew Nelson _____________________________________

2 1

Possible test failure on Debian 6 with Python 2.6
by Chris Ball July 16, 2012

July 16, 2012

Hi, I'm trying to set up various build machines. Some of these are with ShiningPanda.com, which provides a 64-bit Debian 6 machine (as well as Windows 7). This machine has multiple versions of Python installed. Using the build procedure below, I see a test failure with Python 2.6 (and 2.7) but not 2.5 (and 2.4). Any ideas? I guess other people are testing on Debian 6, and if they don't see the test failure with Python 2.6 maybe there's something unusual about the setup of this machine. Build procedure: $ python setup.py build install --prefix=./numpy-install $ cd numpy-install $ export PYTHONPATH=/path/to/numpy-install/lib/python2.6/site-packages $ python ../tools/test-installed-numpy.py --coverage -- --with-xunit Test error: ====================================================================== ERROR: test_multiarray.TestNewBufferProtocol.test_roundtrip ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/slave/jenkins/shiningpanda/jobs/51551925/virtualenvs/d41d8cd9/lib/python2 .6/site-packages/nose/case.py", line 197, in runTest self.test(*self.arg) File "/home/slave/jenkins/workspace/NumPy-Debian6_x64-py26/numpy- install/lib/python2.6/site-packages/numpy/core/tests/test_multiarray.py", line 2537, in test_roundtrip assert_raises(ValueError, self._check_roundtrip, x) File "/home/slave/jenkins/workspace/NumPy-Debian6_x64-py26/numpy- install/lib/python2.6/site-packages/numpy/testing/utils.py", line 1018, in assert_raises return nose.tools.assert_raises(*args,**kwargs) File "/sp/lib/python/cpython-2.6/lib/python2.6/unittest.py", line 336, in failUnlessRaises callableObj(*args, **kwargs) File "/home/slave/jenkins/workspace/NumPy-Debian6_x64-py26/numpy- install/lib/python2.6/site-packages/numpy/core/tests/test_multiarray.py", line 2471, in _check_roundtrip y = np.asarray(x) RuntimeWarning: tp_compare didn't return -1 or -2 for exception [...] Ran 4781 tests in 68.179s FAILED (KNOWNFAIL=4, SKIP=6, errors=1) The full build and test output is here: https://jenkins.shiningpanda.com/scipy/job/NumPy-Debian6_x64-py26/8/console… For comparison, with Python 2.5: https://jenkins.shiningpanda.com/scipy/job/NumPy-Debian6_x64-py25/9/console…

2 1

routine for linear least norms problems with specifiable accuracy
by Dmitrey July 16, 2012

July 16, 2012

hi all, I have wrote a routine to solve dense / sparse problems min {alpha1*||A1 x - b1||_1 + alpha2*||A2 x - b2||^2 + beta1 * ||x||_1 + beta2 * ||x||^2} with specifiable accuracy fTol > 0: abs(f-f*) <= fTol (this parameter is handled by solvers gsubg and maybe amsg2p, latter requires known good enough fOpt estimation). Constraints (box-bound, linear, quadratic) also could be easily connected. This problem is very often encountered in many areas, e.g. machine learning, sparse approximation, see for example http://scikit-learn.org/stable/modules/ ? lastic-net First of all solver large-scale gsubg is recommended. Some hand-tuning of its parameters also could essentially speedup the solver. Also you could be interested in other OpenOpt NSP solvers - ralg and amsg2p (they are medium-scaled although). You can see the source of the routine and its demo result here. You shouldn't expect gsubg will always solve your problem and inform of obtained result with specifiable accuracy - for some very difficult, e.g. extremely ill-conditioned problems it may * fail to solve QP subproblem (default QP solver is cvxopt, you may involve another one, e.g. commercial or free-for-educational cplex) * exit with another stop criterion, e.g. maxIter has been reached, or maxShoots have been exceeded (usually latter means you have reached solution, but it cannot be guaranteed in the case) First of all I have created the routine to demonstrate gsubg abilities; I haven't decided yet commit or not commit the routine to OpenOpt, with or without special class for this problem; in either case you can very easily create problems like this one in FuncDesigner (without having to write a routine for derivatives) to solve them by gsubg or another NSP solver; however, IIRC FuncDesigner dot() doesn't work with sparse matrices yet

2 2

Help building NumPy on Windows
by Chris Ball July 16, 2012

July 16, 2012

Hi, I'm having some trouble building numpy on a 64-bit Windows 7 machine. I'm probably accidentally missing a step following the build process described at http://scipy.org/Installing_SciPy/Windows; it would be great if someone could spot what! Here's what I did: 1. installed python 2.7 from python.org 2. installed mingw32 from the link above (results in gcc 3.4.5) 3. added c:\mingw\bin to front of path 4. ran "python setup.py build --compiler=mingw32 install --prefix=numpy-install" I got the following error: [...] C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\system_info.py:1409: UserWarning: Lapack (http://www.netlib.org/lapack/) sources not found. Directories to search for the sources can be specified in the numpy/distutils/site.cfg file (section [lapack_src]) or by setting the LAPACK_SRC environment variable. warnings.warn(LapackSrcNotFoundError.__doc__) customize GnuFCompiler gnu: no Fortran 90 compiler found gnu: no Fortran 90 compiler found customize GnuFCompiler using config C compiler: gcc -mno-cygwin -O2 -Wall -Wstrict-prototypes compile options: '-DNPY_MINGW_USE_CUSTOM_MSVCR -D__MSVCRT_VERSION__=0x0900 - Inumpy\core\src\private -Inumpy\core\src -Inumpy\core -Inumpy\core\src\npymath - Inumpy\core\src\multiarray -Inumpy\core\src\umath -Inumpy\core\src\npysort - Inumpy\core\include -Ic:\Python27\include -Ic:\Python27\PC -c' gcc -mno-cygwin -O2 -Wall -Wstrict-prototypes -DNPY_MINGW_USE_CUSTOM_MSVCR - D__MSVCRT_VERSION__=0x0900 -Inumpy\core\src\private -Inumpy\core\src - Inumpy\core -Inumpy\core\src\npymath -Inumpy\core\src\multiarray - Inumpy\core\src\umath -Inumpy\core\src\npysort -Inumpy\core\include - Ic:\Python27\include -Ic:\Python27\PC -c _configtest.c -o _configtest.o Found executable c:\mingw\bin\gcc.exe g++ -mno-cygwin _configtest.o -lmsvcr90 -o _configtest.exe Found executable c:\mingw\bin\g++.exe c:\mingw\bin\..\lib\gcc\mingw32\3.4.5\..\..\..\..\mingw32\bin\ld.exe: cannot find -lmsvcr90 collect2: ld returned 1 exit status failure. removing: _configtest.exe.manifest _configtest.c _configtest.o Traceback (most recent call last): File "setup.py", line 214, in <module> setup_package() File "setup.py", line 207, in setup_package configuration=configuration ) File "C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\core.py", line 186, in setup return old_setup(**new_attr) File "c:\Python27\lib\distutils\core.py", line 152, in setup dist.run_commands() File "c:\Python27\lib\distutils\dist.py", line 953, in run_commands self.run_command(cmd) File "c:\Python27\lib\distutils\dist.py", line 972, in run_command cmd_obj.run() File "C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\command\build.py", line 37, in run old_build.run(self) File "c:\Python27\lib\distutils\command\build.py", line 127, in run self.run_command(cmd_name) File "c:\Python27\lib\distutils\cmd.py", line 326, in run_command self.distribution.run_command(command) File "c:\Python27\lib\distutils\dist.py", line 972, in run_command cmd_obj.run() File "C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\command\build_src.py", line 152, in run self.build_sources() File "C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\command\build_src.py", line 163, in build_sources self.build_library_sources(*libname_info) File "C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\command\build_src.py", line 298, in build_library_sources sources = self.generate_sources(sources, (lib_name, build_info)) File "C:\Users\ceball\npslavetest\Windows7_x64_mingw32_3.4.5- py27\build\numpy\distutils\command\build_src.py", line 385, in generate_sources source = func(extension, build_dir) File "numpy\core\setup.py", line 648, in get_mathlib_info raise RuntimeError("Broken toolchain: cannot link a simple C program") RuntimeError: Broken toolchain: cannot link a simple C program Can anyone see what I've missed from this? I'm not sure what the situation is for supporting building with Microsoft compilers, but on another Windows 7 machine (also 64 bits) with MS Visual Studio 9.0 installed, the build completes but I get test failures. You can see the full output of this here: https://jenkins.shiningpanda.com/scipy/job/NumPy-Windows7_x64- py27/26/consoleFull

2 1

Istalling Numpy and Scipy on preinstalled Python 2.6 on Mac
by Naser Nikandish July 16, 2012

July 16, 2012

Hi, I need to install numpy and scipy on preinstalled Python 2.6 on my Mac Lion. Is there anyway to do it? I am aware that Lion OS comes with Python 2.7 as well. But I have to install it on Python 2.6. I really appreciate any help. Cheers

3 2

Select-based median (in light of code freeze)
by Mike Ressler July 15, 2012

July 15, 2012

Hi, A couple of years ago there was a flurry of work partially at my instigation at SciPy 2009 to build a better median function based on a select algorithm rather than a sort algorithm. It seemed that it had progressed quite far, but the code in lib/function_base.py still uses a sort. Has the select work been abandoned or could it be merged in with a little spit-and-polish? (Or has it been and I'm just too naive to have seen it?) It seemed that the work had gotten quite close to release - I would hate to see it lost. (Not to mention I am greedy for the speed-up.) Regards, Mike -- mike.ressler(a)alum.mit.edu

4 5

ufunc and nditer flags (was Re: Code Freeze for NumPy 1.7)
by Nathaniel Smith July 15, 2012

July 15, 2012

On Sun, Jul 15, 2012 at 6:18 PM, jay bourque <jay.bourque(a)continuum.io> wrote: > Just added PR #359. The purpose is to allow the nditer object operand and > iter flags to be set for a ufunc to provide better control over how an array > is iterated over by a ufunc and how the ufunc uses the operands passed to > it. One specific motivation for this is to be able to specify an input > operand to a ufunc as being read/write instead of read only. Huh. My first gut reaction to this is that it's an argument *against* merging this change, because ufuncs *shouldn't* be writing to their inputs. Maybe I'm wrong, but... obviously there is more context here than we've heard so far. Can you explain what you're actually trying to accomplish? -n

3 2

Extracting sub-fields from an array as a view (PR 350)
by Travis Oliphant July 15, 2012

July 15, 2012

In https://github.com/numpy/numpy/pull/350/files , javius provides a patch to allow field extraction from a structured array to return a view instead of a copy. Generally, this is consistent with the desire to have NumPy return views whenever it can. The same idea underlies the change to the diagonal method. Suppose 'myarr' is a structured array with fields ['lat', 'long', 'meas1', 'meas2', 'meas3', 'meas4']. Currently, myarr[['lat', 'long', 'mesa3']] will return a copy of the data in the underlying array. The proposal is to have this return a view, but do it in a two-stage approach so that a first version returns a copy with the WARN_ON_WRITE flag set introduced in NumPy 1.7. A later version will remove the flag (and the copy). What are thoughts on this proposal and which version of NumPy it should go in? -Travis

4 3

Matrix rank default tolerance - is it too low?
by Matthew Brett July 14, 2012

July 14, 2012

Hi, I noticed that numpy.linalg.matrix_rank sometimes gives full rank for matrices that are numerically rank deficient: If I repeatedly make random matrices, then set the first column to be equal to the sum of the second and third columns: def make_deficient(): X = np.random.normal(size=(40, 10)) deficient_X = X.copy() deficient_X[:, 0] = deficient_X[:, 1] + deficient_X[:, 2] return deficient_X then the current numpy.linalg.matrix_rank algorithm returns full rank (10) in about 8 percent of cases (see appended script). I think this is a tolerance problem. The ``matrix_rank`` algorithm does this by default: S = spl.svd(M, compute_uv=False) tol = S.max() * np.finfo(S.dtype).eps return np.sum(S > tol) I guess we'd we want the lowest tolerance that nearly always or always identifies numerically rank deficient matrices. I suppose one way of looking at whether the tolerance is in the right range is to compare the calculated tolerance (``tol``) to the minimum singular value (``S.min()``) because S.min() in our case should be very small and indicate the rank deficiency. The mean value of tol / S.min() for the current algorithm, across many iterations, is about 2.8. We might hope this value would be higher than 1, but not much higher, otherwise we might be rejecting too many columns. Our current algorithm for tolerance is the same as the 2-norm of M * eps. We're citing Golub and Van Loan for this, but now I look at our copy (p 261, last para) - they seem to be suggesting using u * |M| where u = (p 61, section 2.4.2) eps / 2. (see [1]). I think the Golub and Van Loan suggestion corresponds to: tol = np.linalg.norm(M, np.inf) * np.finfo(M.dtype).eps / 2 This tolerance gives full rank for these rank-deficient matrices in about 39 percent of cases (tol / S.min() ratio of 1.7) We see on p 56 (section 2.3.2) that: m, n = M.shape 1 / sqrt(n) . |M|_{inf} <= |M|_2 So we can get an upper bound on |M|_{inf} with |M|_2 * sqrt(n). Setting: tol = S.max() * np.finfo(M.dtype).eps / 2 * np.sqrt(n) gives about 0.5 percent error (tol / S.min() of 4.4) Using the Mathworks threshold [2]: tol = S.max() * np.finfo(M.dtype).eps * max((m, n)) There are no false negatives (0 percent rank 10), but tol / S.min() is around 110 - so conservative, in this case. So - summary - I'm worrying our current threshold is too small, letting through many rank-deficient matrices without detection. I may have misread Golub and Van Loan, but maybe we aren't doing what they suggest. Maybe what we could use is either the MATLAB threshold or something like: tol = S.max() * np.finfo(M.dtype).eps * np.sqrt(n) - so 2 * the upper bound for the inf norm = 2 * |M|_2 * sqrt(n) . This gives 0 percent misses and tol / S.min() of 8.7. What do y'all think? Best, Matthew [1] http://matthew-brett.github.com/pydagogue/floating_error.html#machine-epsil… [2] http://www.mathworks.com/help/techdoc/ref/rank.html Output from script: Percent undetected current: 9.8, tol / S.min(): 2.762 Percent undetected inf norm: 39.1, tol / S.min(): 1.667 Percent undetected upper bound inf norm: 0.5, tol / S.min(): 4.367 Percent undetected upper bound inf norm * 2: 0.0, tol / S.min(): 8.734 Percent undetected MATLAB: 0.0, tol / S.min(): 110.477 <script> import numpy as np import scipy.linalg as npl M = 40 N = 10 def make_deficient(): X = np.random.normal(size=(M, N)) deficient_X = X.copy() if M > N: # Make a column deficient deficient_X[:, 0] = deficient_X[:, 1] + deficient_X[:, 2] else: # Make a row deficient deficient_X[0] = deficient_X[1] + deficient_X[2] return deficient_X matrices = [] ranks = [] ranks_inf = [] ranks_ub_inf = [] ranks_ub2_inf = [] ranks_mlab = [] tols = np.zeros((1000, 6)) for i in range(1000): m = make_deficient() matrices.append(m) # The SVD tolerances S = npl.svd(m, compute_uv=False) S0 = S.max() # u in Golub and Van Loan == numpy eps / 2 eps = np.finfo(m.dtype).eps u = eps / 2 # Current numpy matrix_rank algorithm ranks.append(np.linalg.matrix_rank(m)) # Which is the same as: tol_s0 = S0 * eps # ranks.append(np.linalg.matrix_rank(m, tol=tol_s0)) # Golub and Van Loan suggestion tol_inf = npl.norm(m, np.inf) * u ranks_inf.append(np.linalg.matrix_rank(m, tol=tol_inf)) # Upper bound of |X|_{inf} tol_ub_inf = tol_s0 * np.sqrt(N) / 2 ranks_ub_inf.append(np.linalg.matrix_rank(m, tol=tol_ub_inf)) # Times 2 fudge tol_ub2_inf = tol_s0 * np.sqrt(N) ranks_ub2_inf.append(np.linalg.matrix_rank(m, tol=tol_ub2_inf)) # MATLAB algorithm tol_mlab = tol_s0 * max(m.shape) ranks_mlab.append(np.linalg.matrix_rank(m, tol=tol_mlab)) # Collect tols tols[i] = tol_s0, tol_inf, tol_ub_inf, tol_ub2_inf, tol_mlab, S.min() rel_tols = tols / tols[:, -1][:, None] fmt = 'Percent undetected %s: %3.1f, tol / S.min(): %2.3f' max_rank = min(M, N) for name, ranks, mrt in zip( ('current', 'inf norm', 'upper bound inf norm', 'upper bound inf norm * 2', 'MATLAB'), (ranks, ranks_inf, ranks_ub_inf, ranks_ub2_inf, ranks_mlab), rel_tols.mean(axis=0)[:5]): pcnt = np.sum(np.array(ranks) == max_rank) / 1000. * 100 print fmt % (name, pcnt, mrt) </script>

5 18

New Python tool for searching maximum stable set of a graph
by Dmitrey July 14, 2012

July 14, 2012

Hi all, In the OpenOpt software (BSD-licensed, http://openopt.org ) we have implemented new class - STAB - searching for maximum stable set of a graph. networkx graphs are used as input arguments. Unlike networkx maximum_independent_set() we focus on searching for exact solution (this is NP-Hard problem). interalg or OpenOpt MILP solvers are used, some GUI features and stop criterion (e.g. maxTime, maxCPUTime, fEnough) can be used. Optional arguments are includedNodes and excludedNodes - nodes that have to be present/absent in solution. See http://openopt.org/STAB for details. Future plans (probably very long-term although) include TSP and some other graph problems. ------------------------- Regards, Dmitrey.

1 0