Mailman 3 January 2012 - SciPy-User

AUTO: Christopher Peters is out of the office (returning 01/17/2012)
by Christopher Peters Jan. 16, 2012

Jan. 16, 2012

I am out of the office until 01/17/2012. I am out of the office. Please email urgent requests to Mike McDonald. Note: This is an automated response to your message "[SciPy-User] [ANN] global constrained solver with discrete variables" sent on 1/16/2012 3:07:54 PM. This is the only notification you will receive while this person is away.

1 0

[ANN] global constrained solver with discrete variables
by Dmitrey Jan. 16, 2012

Jan. 16, 2012

hi all, I've done support of discrete variables for interalg - free (license: BSD) solver with specifiable accuracy, you can take a look at an example here It is written in Python + NumPy, and I hope it's speed will be essentially increased when PyPy (Python with dynamic compilation) support for NumPy will be done (some parts of code are not vectorized and still use CPython cycles). Also, NumPy funcs like vstack or append produce only copy of data, and it also slows the solver very much (for mature problems). Maybe some bugs still present somewhere - interalg code already became very long, but since it already works, you could be interested in trying to use it right now. Regards, D.

1 0

Python for Scientists - courses in Germany and US
by Mike Müller Jan. 13, 2012

Jan. 13, 2012

Learn NumPy and Much More ========================= Scientists like Python. If you would like to learn more about important libraries for scientific applications, you might be interested in these courses. The course in Germany covers: - Overview of libraries - NumPy - Data storage with text files, Excel, netCDF and HDF5 - matplotlib - Object oriented programming for scientists - Problem solving session The course in the USA covers all this plus: - Extending Python in other languages - Version control - Unit testing More details below. If you have any questions about the courses, please contact me. Mike Python for Scientists and Engineers (Germany) --------------------------------------------- A three-day course covering all the basic tools scientists and engineers need. This course requires basic Python knowledge. Date: 19.01.-21.01.2012 Location: Leipzig, Germany Trainer: Mike Müller Course Language: English Link: http://www.python-academy.com/courses/python_course_scientists.html Python for Scientists and Engineers (USA) ----------------------------------------- This is an extend version of our well-received course for scientists and engineers. Five days of intensive training will give you a solid basis for using Python for scientific an technical problems. The course is hosted by David Beazley (http://www.dabeaz.com). Date: 27.02.-02.03.2012 Location: Chicago, IL, USA Trainer: Mike Müller Course Language: English Link: http://www.dabeaz.com/chicago/science.html

1 0

Scipy.optimize.fmin_slsqp problem
by Hendrik Jan. 12, 2012

Jan. 12, 2012

Hey Folks!! I'm new to the Forum world. I looked for a topic on optimization, but didn't see one, so sorry if I'm duplicating... I run Python 2.6.6.1 I'm using scipy.optimize.fmin_slsqp to solve a non-linear function but it only returns the function value with the initial guess values... Can you please help??? Below this message is my program... I'm struggling with this problem for over 3 months now, asking for advice with no success!! I want to minimize the function "changevars(self.gibbs)" within the constraints "changevars(self.atombalance)". The optimization program only returns the function value with the initial guess values... And i don't know why. When I run the program, I can see the minimum value graphically, but the optimization can't calculate it... Can you please help?? Regards ************************************************************************************************************************************* #!/usr/bin/env python # Calculate mixture equilibriums using minimisation of the Gibbs energy # 201110 Originally by Hendrik Venter # 201111 Significantly reworked by Carl Sandrock from __future__ import division import scipy.optimize import scipy.linalg import math import numpy import atomparser import sys import matplotlib.pyplot as pl from mpl_toolkits.mplot3d import axes3d from matplotlib import cm smallvalue = 1e-10 fig = pl.figure() ax = fig.add_subplot(111, projection='3d') T = 298.15 R = 8.314 RT = R*T mustplot = 1 # 1 for Plotting, 0 for not plotting class compound: """ Basic container for compound properties """ def __init__(self, name, DGf): self.name = name self.DGf = DGf self.parsed = atomparser.parseformula(name) class mixture: """ Container for mixture properties """ def __init__(self, charge): """ Initialise mixture - charge contains tuples of (initialcharge, compound) """ self.N, self.compounds = zip(*charge) self.N = numpy.array(self.N) self.compoundnames = [c.name for c in self.compounds] self.DGf = numpy.array([c.DGf for c in self.compounds]) self.elements = reduce(set.union, (c.parsed.distinctelements() for c in self.compounds)) self.S = numpy.array([c.parsed.counts(self.elements) for c in self.compounds]).T Srank = numpy.count_nonzero(numpy.linalg.svd(self.S, compute_uv=False) > smallvalue) # Calculating the Degrees of freedom self.DOF = self.S.shape[1] - Srank print 'Degrees of Freedom =', self.DOF # Coefficients of the Gibbs function Ncomps = len(self.compounds) # self.A = numpy.tile(numpy.repeat([-1, 1], [1, Ncomps-1]), [Ncomps, 1]) # number of atoms of each element self.atoms = numpy.dot(self.S, self.N) def gibbs(self, N=None): """ Gibbs energy of mixture with N of each compound """ if N is None: N = self.N # TODO: There is every chance that this function is not correct. It needs to be checked. #return sum(N*(self.DGf/(R*T) + numpy.log(numpy.dot(self.A, N)))) logs = numpy.log(sum(N)) return sum(N*(self.DGf/RT + numpy.log(N) - logs)) def atombalance(self, N): """ Atom balance with N of each compound """ #if N is None: N = self.N return numpy.dot(self.S, N) - self.atoms def conversion(self, conversionspec): """ Calculate the composition given a certain conversion. conversionspec is a list of 2-tuples containing a component index or name and a conversion """ #TODO: A and B should only be calculated once #TODO: This does not take into account any existing products in the mixture #TODO: This works only for conversion specified in terms of reagents if len(conversionspec) < self.DOF: raise Exception("Not enough conversions specified.") C = numpy.zeros([len(conversionspec), self.S.shape[1]]) Ic = C.copy() for i, (j, c) in enumerate(conversionspec): if type(j) is str: j = self.compoundnames.index(j) C[i, j] = 1-c Ic[i, j] = 1 A = numpy.vstack([self.S, C]) B = numpy.vstack([self.S, Ic]) # A ni = B nf nf, _, _, _ = scipy.linalg.lstsq(B, numpy.dot(A, self.N)) # assert residuals are neglegable nf[nf<0] = 0 return nf def equilibrium(self, initialconversion): """ Return equilibrium composition as minimum of Gibbs Energy """ # guess initial conditions N0 = self.conversion(initialconversion) logN0 = numpy.log(N0) # This decorator modifies a function to be defined in terms of new variables def changevars(f): def newf(newX): # print 'newX', newX # print 'X', numpy.exp(newX) r = f(numpy.exp(newX)) # print 'f(X)', r return r return newf # Find optimal point in terms of a change of variables logN = scipy.optimize.fmin_slsqp(changevars(self.gibbs), logN0, f_eqcons=changevars(self.atombalance), acc = 1.0E-12) scipy.optimize.fmin N = numpy.exp(logN) print '' print 'Calculated Optimmum Values' print N print '' print 'Calculated Optimum Function value' print m.gibbs(N) return N # Here for a particular mixture: m = mixture([[1., compound('MgO', -568343.0)], [0.5, compound('Al2O3', -1152420.0)], [4., compound('H2O', -237141.0)], [0, compound('Mg(OH)2', -833644.0)], [0, compound('Al(OH)3', -1835750.0)]]) # find Gibbs energy surface for many conversion possibilities Nsteps = 10 convrange = numpy.linspace(0.001, 0.999, Nsteps) gibbssurface = numpy.zeros((Nsteps, Nsteps)) for iMgO, xMgO in enumerate(convrange): for iAl2O3, xAl2O3 in enumerate(convrange): N = m.conversion([('MgO', xMgO), ('Al2O3', xAl2O3/50)]) gibbssurface[iMgO, iAl2O3] = m.gibbs(N) # Try to find equilibrium # Initial conversion guess print '' print "***************Optimizing******************" X0 = [0.5, 0.5] conversionspec = zip(['MgO', 'Al2O3'], X0) # Solve to equilibrium N = m.equilibrium(conversionspec) print '' print '*************Final Values***************' for n, init, name in zip(N, m.N, m.compoundnames): print name, 'initial =', init, 'final =', n # Plot results if mustplot: Xaxis, Yaxis = numpy.meshgrid(convrange, convrange) Zaxis = gibbssurface surf = ax.plot_surface(Xaxis, Yaxis, Zaxis, rstride=1, cstride=1, cmap=cm.jet, linewidth=0, antialiased=False) fig.colorbar(surf, shrink=0.5, aspect=5) ax.set_xlabel('Conversion of MgO') ax.set_ylabel('Conversion of Al2O3') ax.set_zlabel('Gibbs free energy (G/RT)') pl.show()

2 1

Lookfor in SciPy sub-packages?
by klo uo Jan. 12, 2012

Jan. 12, 2012

I just discovered this nifty NumPy utils.py module. Among various useful information retrival functions (as it can be imagined by wrapping pydoc, inspect and re modules) there is also 'lookfor', for which I even quickly made GTK input box ;) However seaching for something within SciPy sub-packages seems pretty impossible to me. I tried various combinations but nothing worked. Can someone provide tip how to return result for 'Butterworth' for example, without loading scipy.signal inside script then passing this module to lookfor function, which seems lake a overkill to this function purpose

1 1

ODE Integration error
by Tim Gallagher Jan. 11, 2012

Jan. 11, 2012

Hi, I'm trying to use the 'vode' solver in scipy.integrate.ode (in version 0.10.0) and I get the following error: DVODE-- RWORK length needed, LENRW (=I1), exceeds LRW (=I2) In above message, I1 = 7328 I2 = 56 I don't see any arguments to anything I can use to determine workspace sizes. I'm sure I did something incorrectly in my calling code, but I don't know how to find out from this error. Any suggestions? Tim

2 1

cs_graph_componets returning wrong results?
by Joohwi Lee Jan. 11, 2012

Jan. 11, 2012

Hi All, To identify connected components of a graph, I am trying to use 'scipy.sparse.cs_graph_components()' function. The adjacency graph looks like array([[1, 0, 0, 0, 1], [0, 1, 1, 1, 1], [0, 0, 1, 1, 1], [0, 0, 0, 1, 1], [0, 0, 0, 0, 1]]) The graph is actually consisting of (1 -> (5 = 2 = 3 = 4) ), where (2,3,4,5) are connected each other and 1 is only connected to 5. I expect cs_graph_components(x) should return (1, [0,0,0,0,0]) which means that there is only one component, but it returns (2, array([0, 1, 1, 1, 0], dtype=int32)). (1-5) & (2-3-4). Is this correct result? - joohwi

2 2

numpydoc, show attributes in autoclass output
by alex arsenovic Jan. 10, 2012

Jan. 10, 2012

i am using the numpydoc sphinx extensions, and would like to modify the way in which the autoclass feature works. right now, i am generating the class's page using autosummary and only the class's methods are shown. i would like to have the class's attributes displayed as well. i have looked at the template in sphinx/source/_template/autosummary/class.rst, and understand what its doing and how to work with it, but the numpydoc appears to be doing something itself, modifying sphinx's autoclass behavior. i poked around and found the cfg option 'numpydoc_show_class_members'. while this does control whether or not the class methods are included when it is set in sphinx/source/conf.py, it doesnt control include the class's attributes, even though it appears as though it should by inspecting its usage in sphinx_ext/docscrape.py. at the moment i have disabled the 'numpydoc_show_class_members', and used the sphinx class.rst template to do what i need, but would be interested to have it working through the numpydoc. thanks! alex

1 0

Trouble using curve_fit
by klo uo Jan. 9, 2012

Jan. 9, 2012

Example provided at the bottom of this page: scipy.optimize.curve_fit<http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fi…>seemed perfect for what I hoped I need, but result was bad. I don't know if I misuse it or I don't have enough data, but as can be seen from attached screen-shot, I did not expect problem for finding obvious exponential function fit. Here is snippet: ======================================== import numpy as np from scipy.optimize import curve_fit def func(x, a, b, c): return a*np.exp(-b*x)+c x = np.array([ 0.8, 11., 12., 16., 32., 37.8, 44., 48., 56., 64., 88., 96., 112., 128., 144., 176., 192. ]) y = np.array([0.94597685600279, 0.95856916599601, 0.96009142950541, 0.96454515552826, 0.97938932735214, 0.98252400815195, 0.98500175787242, 0.98621192462708, 0.98816995007392, 0.98964101933472, 0.99247255046129, 0.99308203517541, 0.99406737810867, 0.99480702681278, 0.99538268958706, 0.99622916581118, 0.99653501465135]) popt, pcov = curve_fit(func, x, y) popt Out[]: array([ 164.57396829 206.83406927 0.98251156]) ======================================== I then found this page: wolfram<http://mathworld.wolfram.com/LeastSquaresFittingExponential.html>and quickly ======================================== def lsq_exp(x, y): # return A and B in y=A*exp(B*x) from numpy import exp, log, sum n = len(x) a_num = (log(y)).sum() * (x**2).sum() - x.sum() * (x*log(y)).sum() b_num = n*(x*log(y)).sum() - x.sum() * (log(y)).sum() den = n*(x**2).sum() - (x.sum())**2 A = exp(a_num/den) b = b_num/den return A,b ======================================== result coefficients: (0.96677003283750318, 0.00021665429178716782) What is wrong with this approach?

2 2

[SciPy-user] optimization with nonlinear constraints
by Dr. Phillip M. Feldman Jan. 9, 2012

Jan. 9, 2012

>From the SciPy documentation, it looks as though there's no way to do optimization with nonlinear constraints. Have I overlooked something? If not, then this seems like a fairly important omission. -- View this message in context: http://old.nabble.com/optimization-with-nonlinear-constraints-tp33101263p33… Sent from the Scipy-User mailing list archive at Nabble.com.

3 2