Mailman 3 October 2012 - SciPy-User

Someone with OSX 10.7/10.8 please test
by Pauli Virtanen Oct. 30, 2012

Oct. 30, 2012

Hi, If someone with OSX 10.7/10.8 can lend a hand and test this code change, help would be appreciated: https://github.com/scipy/scipy/pull/280 Namely, what needs to be done: git clone git://github.com/scipy/scipy.git git remote add pv git://github.com/pv/scipy-work.git git fetch pv git checkout pv/accelerate-ffc Now, rebuild Scipy from these sources, install to some temporary location, and check import scipy print scipy.__file__ # <- check you get the right one python -c 'import scipy; scipy.test("full", verbose=2)' > test.log 2>&1 Look for failures that look like this: http://projects.scipy.org/scipy/ticket/1618 There shouldn't be any of those. Then compare to the master branch: git checkout origin/master rm -rf build and rebuild, and re-test. Do the failures reappear? If you have github account, you can reply directly to the pull request. (Test logs can be uploaded to pastebin.) Thanks, -- Pauli Virtanen

2 1

Re: [SciPy-User] fsolve return values
by The Helmbolds Oct. 29, 2012

Oct. 29, 2012

Thanks to all who participated. Yes, the "rank-1 method" is billed as "Broyden's rank-a method", altho that does not make it any clearer to me. And, as you pointed out, it won't be clear to many users, either. So I'll have to either say more or say less about it. As far as how frequently the Jacobian gets updated, I guess that depends on the details of the problem and the guessed starting-point. I do have a toy problem (only 3 simultaneous equations in 3 variables) where the number of iterations made is 12, but the number of Jacobian evaluations is 1. So apparently the initial Jacobian was never updated. That did not affect the validity of the solution, however, as it satisfied the 3 equations to at least 8 significant figures. BTW, I'm sure you all are annoyed at the excess question-marks. Does anyone know how to get rid of them when using Yahoo mail on IE 9? I have tried several things, but none seem to be successsful. Gimme some credit for trying, but little for succeeding. Bob H

1 0

leastsq error
by stephen lukacs Oct. 29, 2012

Oct. 29, 2012

hello one and all, i am having a terrible time with optimize.leastsq, even fitting a line, so please help. here is my python code and error >>> import numpy >>> from scipy import optimize >>> def expr_conductance(x, a, b, c): ... return a*0 + b*x + c ... >>> def residual_conductance(p,x,y): ... a, b, c = p ... return y - expr_conductance(x, a, b, c) ... >>> x = [0.99771057137610752, 0.49976145827781415, 0.24821831884394957, 0.12480215949109599, 0.06315070141095365, 0.03065779901355976, 0.015669312142317458, 0.0078799613766362755, 0.0039027338918740067] >>> len(x) 9 >>> y = [2.9954211427522148, 1.9995229165556283, 1.496436637687899, 1.2496043189821919, 1.1263014028219074, 1.0613155980271196, 1.031338624284635, 1.0157599227532725, 1.007805467783748] >>> len(y) 9 >>> [a, b, c], conv = optimize.leastsq(residual_conductance, [0,10000.,5.], args = (x,y)) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 276, in leastsq m = _check_func('leastsq', 'func', func, x0, args, n)[0] File "/usr/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 13, in _check_func res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) File "<stdin>", line 3, in residual_conductance ValueError: operands could not be broadcast together with shapes (9) (90000) >>> i have tried many permutations and looked all over for syntax, but i can not find why this error is there or how to deal with it. i have scipy 0.10.1 and numpy 1.6.1 under python 2.6.6 on a centos 6.3 system. ultimately i want to put make the fit function a*numpy.power(x,(3/2))+b*x+c, which is a bit nonlinear but at this point i can't even get a line to fit. thank you in advance and have a great day. lucas

3 2

Re: [SciPy-User] The QR decomposition values returned by fsolve
by The Helmbolds Oct. 29, 2012

Oct. 29, 2012

I think the following is the source of my confusion. SciPy's docstring for fsolve omits the following information found in the "User Guide for MINPACK" regarding HYBRD and HYBRDJ. The following is not a direct quote, but it's pretty close: The initial value of the Jacobian is not updated until the rank-1 method fails to produce satisfactory progress. I assume that the Jacobian gets updated intermittently, and only when the rank-1 method is not producing satisfactory progress. (So in fact it might never get updated!!) Because fsolve's QR-related outputs (`fjac`, `r`, and `qtf`) are based on the final value of fsolve's internal "approximate Jacobian", they may be quite wide of the mark, unless fsolve "just happens" to return right after the Jacobian has been updated. Accordingly -- unless there is some objection -- in my revision of fsolve's docstring, I'll add to the Notes section something like the following: **Cautionary Note**: According to [the MINPACK User Guide], the initial value of the program's "approximate Jacobian" is estimated (or calculated if `fprime` is supplied by the user), but is updated only when the rank-1 method is not producing satisfactory progress. Because the program's QR-related outputs (`fjac`, `r`, and `qtf`) are based on the program's internal "approximate Jacobian", they should not be used in subsequent analysis unless their validity is confirmed by independent computations. Bob H

3 2

Orthogonal polynomials on the unit circle
by josef.pktd＠gmail.com Oct. 27, 2012

Oct. 27, 2012

http://en.wikipedia.org/wiki/Orthogonal_polynomials_on_the_unit_circle with link to handbook application: goodness of fit for circular data http://onlinelibrary.wiley.com/doi/10.1111/j.1467-842X.2009.00558.x/abstract Are those available anywhere in python land? What's the difference between orthogonal polynomials on the unit circle and periodic polynomials like Fourier series? Josef circular statistics - what's that? It's like TDD, you go in circles

4 8

Re: [SciPy-User] SciPy-User Digest, Vol 110, Issue 46
by The Helmbolds Oct. 27, 2012

Oct. 27, 2012

In response to: > Date: Fri, 26 Oct 2012 16:53:47 -0400 > From: "Moore, Eric (NIH/NIDDK) [F]" <eric.moore2(a)nih.gov> > Subject: Re: [SciPy-User] Request help with fsolve outputs >> -----Original Message----- >> From: The Helmbolds [mailto:helmrp@yahoo.com] >> Sent: Friday, October 26, 2012 11:52 AM >> To: User SciPy >> Subject: [SciPy-User] Request help with fsolve outputs >> >> ----SNIP---- >> >> So please help me out here. What are the fjac and r values that fsolve >> returns (where fjac and r are suppose to be the QR factors for the final >>Jacobian) ? >> How are they related to the above Qend, Rend, and Jend (i.e., the apparently >> correct values from an independent computation)? >> How is the user supposed to use them (i.e., the fjac and r values returned by >> fsolve)? >> Bob H >> > Eric Moore responded: > > I haven't spent any time playing with your example, but I have looked at a > simpler example. It seems that the approximated jacobian can be quite far off > and fsolve can still return correct zeros. How good the approximation ends up > being depends on your choice of initial conditions. I've attached a trivial > example that shows this. > > There is also the epsfcn parameter, but based on the minpack documentation, the > default should be okay for both our examples since we have full machine > precision available in our functions to minimize. > > It would make things much easier for other people if you could post a python > file where you have defined your function and are calling fsolve. I believe > Ralph asked if you could do this the last time you posted this question. > I'm sorry that no one has come in with a simple answer, but this would still > be a good step to make considering your question a little easier for other > people. > > Good luck, > > Eric Eric, thanks for your comments. I gather from them that the descriptive write-up might be improved if it included a statement more or less as follows: If the user's application requires accurate values of the problem's Jacobian and its QR decomposition values evaluated at fsolve's 'root' value, they should not rely on the values for those quantities returned by fsolve, unless confirmed by an independent computation. Bob H

1 0

Re: [SciPy-User] Ubuntu / Scipy / Lapack / BLAS issues
by Tanmoy Laskar Oct. 26, 2012

Oct. 26, 2012

Hi scipy-users! I'm on Ubuntu 12.10 (Quantal), and it appears that my recent upgrade to 12.10 has broken my scipy 0.10.0 installation. A traceback is below. The choking point appears to be lapack in some way. I originally installed scipy through my package manager (sudo apt-get install python-scipy). I would greatly appreciate any help. Please let me know if there is any additional info I could provide to help diagnose the problem. Thanks in advance, Tanmoy In [1]: import numpy In [2]: import scipy In [3]: import scipy.linalg --------------------------------------------------------------------------- ImportError Traceback (most recent call last) /home/tanmoy/Projects/Edo/Reverse_Shocks/GRB/120521C/<ipython-input-3-6362d33d4139> in <module>() ----> 1 import scipy.linalg /usr/local/lib/python2.7/dist-packages/scipy/linalg/__init__.py in <module>() 114 115 from misc import * --> 116 from basic import * 117 from decomp import * 118 from decomp_lu import * /usr/local/lib/python2.7/dist-packages/scipy/linalg/basic.py in <module>() 10 11 from flinalg import get_flinalg_funcs ---> 12 from lapack import get_lapack_funcs 13 from misc import LinAlgError, _datacopied 14 from scipy.linalg import calc_lwork /usr/local/lib/python2.7/dist-packages/scipy/linalg/lapack.py in <module>() 13 14 from scipy.linalg import flapack ---> 15 from scipy.linalg import clapack 16 _use_force_clapack = 1 17 if hasattr(clapack,'empty_module'): ImportError: /usr/local/lib/python2.7/dist-packages/scipy/linalg/clapack.so: undefined symbol: clapack_sgesv

2 2

Request help with fsolve outputs
by The Helmbolds Oct. 26, 2012

Oct. 26, 2012

Please help me out here. I’m trying to rewrite the docstring for the `fsolve.py` routine located on my machine in: C:/users/owner/scipy/scipy/optimize/minpack.py The specific issue I’m having difficulty with is understanding the outputs described in fsolve’s docstring as: 'fjac': the orthogonal matrix, q, produced by the QR factorization of the final approximate Jacobian matrix, stored column wise. 'r': upper triangular matrix produced by QR factorization of same matrix. These are described in SciPy’s minpack/hybrd.f file as: ‘fjac’ is an output n by n array which contains the orthogonal matrix q produced by the qr factorization of the final approximate jacobian. ‘r’ is an output array of length lr which contains the upper triangular matrix produced by the qr factorization of the final approximate jacobian, stored rowwise. For ease in writing, in what follows let’s use the symbols ‘Jend’ for the final approximate Jacobian matrix, and use ‘Q’ and ‘R’ for its QR decomposition matrices. Now consider the problem of finding the solution to the following three nonlinear equations (which we will refer to as 'E'), in three unknowns (u, v, w): 2 * a * u + b * v + d - w * v = 0 b * u + 2 * c * v + e - w * u = 0 -u * v + f = 0 where (a, b, c, d, e, f ) = (2, 3, 7, 8, 9, 2). For inputs to fsolve, we identify (u, v, w) = (x[0], x[1], x[2]). Now fsolve gives the solution array: [uend vend wend] = [ 1.79838825 1.11210691 16.66195357]. With these values, the above three equations E are satisfied to an accuracy of about 9 significant figures. The Jacobian matrix for the three LHS functions in E is: J = np.matrix([[2*a, b-w, -v], [b-w, 2*c, -u], [-v, -u, 0.]]) Note that it’s symmetrical, and if we compute its value using the above fsolve’s ‘end’ solution values we get: Jend = [[ 4. 19.66195357 1.11210691], [ 19.66195357 14. 1.79838825], [ 1.11210691 1.79838825 0. ]] Using SciPy’s linalg package, this Jend has the QR decomposition: Qend = [[-0.28013447 -0.91516674 -0.28981807] [ 0.95679602 -0.24168763 -0.16164302] [ 0.07788487 -0.32257856 0.94333293]] Rend = [[-14.278857 17.08226116 -1.40915124] [ -0. 9.69946027 1.45241144] [ -0. 0. 0.61300558]] and Qend * Rend = Jend to within about 15 significant figures. However, fsolve gives the QR decomposition: qretm = [[-0.64093238 0.75748326 0.1241966 ] [-0.62403598 -0.60841098 0.4903215 ] [-0.44697291 -0.23675978 -0.8626471 ]] rret = [ -7.77806716 30.02199802 -0.819055 -10.74878184 2.00090268 1.02706198] and converting rret to an upper triangular NumPy matrix gives: rretm = [[ -7.77806716 30.02199802 -0.819055 ] [ 0. -10.74878184 2.00090268] [ 0. 0. 1.02706198]] Now qret and rretm bear no obvious relation to Qend and Rend. Although qretm is orthogonal to about 16 significant figures, we find the product: qretm * rretm = [[ 4.98521509 -27.38409295 2.16816676] [ 4.85379376 -12.19513008 -0.2026608 ] [ 3.47658529 -10.87414051 -0.99362993]] which bears no obvious relationship to Jend. The hybrd.f routine in minpack refers to a permutation matrix, p, such that we should have in our notation: p*Jend = qretm*rretm, but fsolve apparently does not return the matrix p, and I don’t see any permutation of Jend that would equal qretm*rretm. The hybrd.f routine does refer to some "scaling" that is going on, but my Fortran is about 40 years too stale for me to interpret it. If we reinterpret rret as meaning the matrix: rretaltm = [[ -7.77806716 30.02199802 -10.74878184] [ 0. -0.819055 2.00090268] [ 0. 0. 1.02706198]] then we get the product: qretm * rretaltm = [[ 4.98521509 -19.86249109 8.53245022] [ 4.85379376 -18.2364849 5.99384603] [ 3.47658529 -13.22510045 3.44468895]] which again bears no obvious relationship to Jend. Using the transpose of qretm in the above product is no help. So please help me out here. What are the fjac and r values that fsolve returns? How are they related to the above Qend, Rend, and Jend? How is the user supposed to use them? Bob H

2 1

tip (maybe): scaling and optimizers
by josef.pktd＠gmail.com Oct. 26, 2012

Oct. 26, 2012

mainly an observation: After figuring out that fmin_slsqp is scale sensitive, I switched to normalizing, rescaling loglikelihood functions in statsmodels. Loglikelihood functions are our main functions for nonlinear optimization. Today I was working by accident on an older branch of statsmodels, and the results I got with fmin_bfgs were awful. After switching to statsmodels master, the results I get with fmin_bfgs are much better (very good: robust and accurate). The impression I got from this and from a discussion with Ian Langmore (on an L1 penalized optimization pull request) is that many scipy optimizers might be scale sensitive in the default settings. Watch the scale of your objective function !? (qualifier: I don't remember if other changes are in statsmodels master and not in my old branch that make optimization more robust.) Josef ------ "anecdotal evidence ain't proof" http://editthis.info/logic/Informal_Fallacies#Anecdotal_Evidence http://sayings.jacomac.de/details.php?id=10 ( http://www.unilang.org/viewtopic.php?f=11&t=38585&start=0&st=0&sk=t&sd=a ) ----

2 2

fmin_slsqp constraint problem
by Neal Becker Oct. 26, 2012

Oct. 26, 2012

I have a ineq constraint: ## constrain poles to be inside unit circle def c(coef, len_z, len_p, dz, dp): p = compose ((coef/opt.scale)[len_z:-1], dp) return np.abs(p) - 1 So this will return a 1D array where each value should satisfy the constraint. fmin_slsqp will not accept this directly: e = fmin_slsqp (obj_fnc, coef*opt.scale, ieqcons=[lambda coef, len_z, len_p, dz, dp: -c(coef, len_z, len_p, dz, dp)], args=(len(lz), len(lp), dz, dp), eqcons=[lambda coef, len_z, len_p, dz, dp: h(coef, len_z, len_p, dz, dp)], full_output=True) Traceback (most recent call last): File "./optimize_pll5.3.2.py", line 519, in <module> run_line (sys.argv) File "./optimize_pll5.3.2.py", line 498, in run_line e = fmin_slsqp (obj_fnc, coef*opt.scale, ieqcons=[lambda coef, len_z, len_p, dz, dp: -c(coef, len_z, len_p, dz, dp)], args=(len(lz), len(lp), dz, dp), eqcons=[lambda coef, len_z, len_p, dz, dp: h(coef, len_z, len_p, dz, dp)], full_output=True) File "/usr/lib64/python2.7/site-packages/scipy/optimize/slsqp.py", line 334, in fmin_slsqp a_ieq[i] = ieqcons_prime[i](x) File "/usr/lib64/python2.7/site-packages/scipy/optimize/optimize.py", line 176, in function_wrapper return function(x, *args) File "/usr/lib64/python2.7/site-packages/scipy/optimize/optimize.py", line 398, in approx_fprime grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon ValueError: setting an array element with a sequence. Any ideas on this?

3 6