Mailman 3 March 2013 - SciPy-User

ANN: Pyzo distro 2013b
by Almar Klein March 19, 2013

March 19, 2013

Dear all, We're pleased to announce release 2013b of Pyzo distro, a Python distribution for scientific computing based on Python 3. For more information, and to give it a try, visit http://www.pyzo.org. The most notable changes since the last release: - Pyzo is also build for Mac. - Python 3.3 is now used. - Added an executable to launch IPython notebook. - New versions of several packages, including the IDE. More information: Pyzo is available for Windows, Mac and Linux. The distribution is portable, thus providing a way to install a scientific Python stack on computers without the need for admin rights. Naturally, Pyzo distro is complient with the scipy-stack, and comes with additional packages such as scikit-image and scikit-learn. With Pyzo we want to make scientific computing in Python easier accessible. We especially hope to make it easier for newcomers (such as Matlab converts) to join our awesome community. Pyzo uses IEP as the default front-end IDE, and IPython (with notebook) is also available. Happy coding, Almar

4 10

scipy.stats beginner: help with fitting distribution
by Arvind Thiagarajan March 18, 2013

March 18, 2013

Hi all, I'm a beginner with SciPy so this may be a basic question. I am trying to fit a Rice distribution to some data using scipy.stats. However, I first tried some test code which doesn't seem to give me a very good fit. I tried the following code: >>> b = [0.3,] >>> samples = rice.rvs(b, loc=0, scale=1, size=1000) >>> rice.fit(samples) (0.0012012190480231357, -0.0023216862043629813, 1.024758538166374) Though the loc and scale seem ok (close to 0 and 1 respectively), I was hoping the first return value would be closer to 0.3 since the shape parameter I supplied while generating random samples was 0.3 (in the call to rvs). I initially thought 1000 samples was perhaps too few, but get similarly poor results with 10,000 samples as well. Concluded that I must be doing something wrong, or that I'm misinterpreting the usage for the rvs() or fit() functions. Is the fit function supposed to try to approximately find the same shape parameter I passed in to rvs? Or should I be interpreting the fit result differently? Wondering if anyone has any ideas? Thanks a lot, Arvind

2 3

OpenOpt Suite release 0.45
by Dmitrey March 15, 2013

March 15, 2013

Hi all, > I'm glad to inform you about new OpenOpt Suite release 0.45 (2013-March-15): * Essential improvements for FuncDesigner interval analysis (thus affect interalg) * Temporary walkaround for a serious bug in FuncDesigner automatic differentiation kernel due to a bug in some versions of Python or NumPy, may affect optimization problems, including (MI)LP, (MI)NLP, TSP etc * Some other minor bugfixes and improvements > > --------------------------- > Regards, D. > http://openopt.org/Dmitrey

1 0

Problem with ndimage.interpolation.zoom
by Eric Emsellem March 14, 2013

March 14, 2013

Hi, I have a problem with the zoom function in the ndimage module of scipy. I am trying to "expand/zoom" an array by a factor of e.g., 3 or 4. I noticed that when the input array has a shape with odd number of e.g. lines/rows etc, something odd happens. When you take a random 3d array of 2x2x2 it works, and with a times-3 zoom, it tri-plicates each value in the output array. However, with the following example you'll see that the border values are replicated only twice and the central value 4 times... What is happening and how can I just get an array with just n times each input value? # Try: test = random.random(9) test3 = test.reshape((3,3)) test3z = scipy.ndimage.interpolation.zoom(test3, 3, mode='constant', prefilter=False, order=0) # and look at the output test3z... the border values are only there twice. Not three times... The central value is there 4 times... Thanks for any input. Eric

2 1

Installing Polymode
by Manders, Mark March 14, 2013

March 14, 2013

I was trying to install polymode into my Python 2.7 (32bit), for some reason the following errors come up. It can't seem to find the numpy 'lapack' libraries even though I have numpy & scipy installed. I also have boost_python installed but it can't find anything from that either. Z:\polymode>setup.py install W: Using scipy hankel ratio - this may be inaccurate! lapack_opt_info: lapack_mkl_info: mkl_info: libraries mkl,vml,guide not found in C:\Python27\lib libraries mkl,vml,guide not found in C:\ libraries mkl,vml,guide not found in C:\Python27\libs NOT AVAILABLE NOT AVAILABLE atlas_threads_info: Setting PTATLAS=ATLAS libraries ptf77blas,ptcblas,atlas not found in C:\Python27\lib libraries lapack_atlas not found in C:\Python27\lib libraries ptf77blas,ptcblas,atlas not found in C:\ libraries lapack_atlas not found in C:\ libraries ptf77blas,ptcblas,atlas not found in C:\Python27\libs libraries lapack_atlas not found in C:\Python27\libs numpy.distutils.system_info.atlas_threads_info NOT AVAILABLE atlas_info: libraries f77blas,cblas,atlas not found in C:\Python27\lib libraries lapack_atlas not found in C:\Python27\lib libraries f77blas,cblas,atlas not found in C:\ libraries lapack_atlas not found in C:\ libraries f77blas,cblas,atlas not found in C:\Python27\libs libraries lapack_atlas not found in C:\Python27\libs numpy.distutils.system_info.atlas_info NOT AVAILABLE C:\Python27\lib\site-packages\numpy\distutils\system_info.py:1340: UserWarning: Atlas (http://math-atlas.sourceforge.net/) libraries not found. Directories to search for the libraries can be specified in the numpy/distutils/site.cfg file (section [atlas]) or by setting the ATLAS environment variable. warnings.warn(AtlasNotFoundError.__doc__) lapack_info: libraries lapack not found in C:\Python27\lib libraries lapack not found in C:\ libraries lapack not found in C:\Python27\libs NOT AVAILABLE C:\Python27\lib\site-packages\numpy\distutils\system_info.py:1351: UserWarning: Lapack (http://www.netlib.org/lapack/) libraries not found. Directories to search for the libraries can be specified in the numpy/distutils/site.cfg file (section [lapack]) or by setting the LAPACK environment variable. warnings.warn(LapackNotFoundError.__doc__) lapack_src_info: NOT AVAILABLE C:\Python27\lib\site-packages\numpy\distutils\system_info.py:1354: 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__) NOT AVAILABLE boost_python_info: NOT AVAILABLE Traceback (most recent call last): File "Z:\polymode\setup.py", line 69, in <module> setup_package() File "Z:\polymode\setup.py", line 63, in setup_package cmdclass = {'doc' : generate_api_docs} File "C:\Python27\lib\site-packages\numpy\distutils\core.py", line 152, in set up config = configuration() File "Z:\polymode\setup.py", line 44, in configuration config.add_subpackage(package_name) File "C:\Python27\lib\site-packages\numpy\distutils\misc_util.py", line 1002, in add_subpackage caller_level = 2) File "C:\Python27\lib\site-packages\numpy\distutils\misc_util.py", line 971, i n get_subpackage caller_level = caller_level + 1) File "C:\Python27\lib\site-packages\numpy\distutils\misc_util.py", line 908, i n _get_configuration_from_setup_py config = setup_module.configuration(*args) File "Polymode\setup.py", line 10, in configuration config.add_subpackage('mathlink') File "C:\Python27\lib\site-packages\numpy\distutils\misc_util.py", line 1002, in add_subpackage caller_level = 2) File "C:\Python27\lib\site-packages\numpy\distutils\misc_util.py", line 971, i n get_subpackage caller_level = caller_level + 1) File "C:\Python27\lib\site-packages\numpy\distutils\misc_util.py", line 908, i n _get_configuration_from_setup_py config = setup_module.configuration(*args) File "Polymode\mathlink\setup.py", line 71, in configuration raise NotFoundError,'no lapack/blas resources found' numpy.distutils.system_info.NotFoundError: no lapack/blas resources found

4 8

ANN: IEP 3.2 - the Interactive Editor for Python
by Almar Klein March 13, 2013

March 13, 2013

Dear all, We're pleased to announce version 3.2 of the Interactive Editor for Python. IEP is a cross-platform Python IDE focused on interactivity and introspection, which makes it very suitable for scientific computing. Its practical design is aimed at simplicity and efficiency. Binaries are available for Windows, Linux, and Mac. A selection of the changes since 3.1: - This is the first release for which all binaries are build with Python 3.3 and PySide? <https://code.google.com/p/iep/w/edit/PySide>. - New file browser tool, which replaces the old file browser and project manager tools. It combines the power of both in one simple interface. It also has functionality for peeking inside python modules. Since its design uses a generalization of the file system, implementing alternative file systems (like zip files or remote machines) should be quite easy in the future. - IEP now also comes with two fonts: 'DejaVu?<https://code.google.com/p/iep/w/edit/DejaVu> Sans Mono' (the default font) and Adobes new 'Source code pro'. - IEP now supports multiple languages. Translations for Dutch, French, Spanish and Catalan are available. Website: http://code.google.com/p/iep/ Discussion group: http://groups.google.com/group/iep_ Release notes: https://code.google.com/p/iep/wiki/Release<https://code.google.com/p/iep/wiki/Release?ts=1363206546&updated=Release#ve…> Happy coding! Almar

1 0

fast N-d interpolators Intergrid and Barypol
by denis March 12, 2013

March 12, 2013

Folks, two small fast interpolators: [Intergrid](http://denis-bz.github.com/docs/intergrid.html): interpolate in an N-d box grid, uniform or non-uniform. This is just a wrapper for scipy.ndimage.map_coordinates and numpy.interp. [Barypol](http://denis-bz.github.com/docs/barypol.html): interpolate in a uniform N-d box grid, using d + 1 corners of a simplex (triangle, tetrahedron ...) around each query point. >From Munos and Moore, "Variable Resolution Discretization in Optimal Control", 1999, 24p; see the pictures and example on pp. 4-5. It's implemented in a C++ header barypol.h, with a Cython wrapper (both of which more knowledgeable people could certainly improve). In 4d, 5d, 6d, Intergrid does around 3M, 2M, .8M interpolations / second. Barypol is ~ 5 times faster, but not as smooth. (Both will of course become cache-bound for large grids.) See interpol/test/*.log and http://github.com/denis-bz/interpol . Comments are welcome, testcases most welcome. cheers -- denis

1 0

HMM in Python
by Subhabrata Banerjee March 11, 2013

March 11, 2013

Dear Group, I do not know whether it is an out of box question but as you are great scientific computing people so I thought to post it here. I am looking a simple example work out for forward probability of Hidden Markov Model(HMM) on the problem of sequence labelling. If any one can kindly help me out. Regards, Subhabrata.

2 2

Python for Card Game Modeling
by Lorenzo Isella March 11, 2013

March 11, 2013

Dear All, Unfortunately, I cannot come up with a script or a concrete example at this stage. I am trying to investigate statistically a card game (not a game like poker, but something more like Magic). There are various kinds of cards (agents, missions, victory points) each obeying a given set of rules. The purpose of the statistical analysis is to come up with a measure of the strength of each card. Any idea if "out there" there is any Python (better if using SciPy/Numpy) package/script which can help me? Any suggestion is welcome Lorenzo

1 0

maximising a function
by mdekauwe March 11, 2013

March 11, 2013

Hi, I am trying to use the scipy library to optimise 4 parameters so that they return the maximum for a given function. To do this I am using the COBYLA function and setting a few constraints, essentially that these parameter values are positive and that when summed they are less than a given value (N_p). Can I check what I have is the correct implementation (see below)? The second issue is that part of the function depends on the solution to a quadratic. If I don't resolve positive roots I want to be able to return some value to let the minimiser know that these are bad values. I can't work out how to do this? I tried returning a very large value but this didn't see to work, so then I tried to use NaNs. What is the correct way to do this? thanks, Code (this will run). #!/usr/bin/env python """ Optimise the distribution of nitrogen within the photosynthetic system to maximise photosynthesis for a given PAR and CO2 concentration. Reference: ========= * Medlyn (1996) The Optimal Allocation of Nitrogen within the C3 Photsynthetic System at Elevated CO2. Australian Journal of Plant Physiology, 23, 593-603. """ import sys import numpy as np import scipy.optimize class PhotosynthesisModel(object): """ Photosynthesis model based on Evans (1989) References: ----------- * Evans, 1989 """ def __init__(self, alpha=0.425, beta=0.7, g_m=0.4, theta=0.7, K_cat=24.0, K_s=1.25E-04, MAXIMISE=False): """ Parameters ---------- alpha : float quantum yield of electron transport (mol mol-1) beta : float constant of proportionality between atmospheric and intercellular CO2 concentrations (-) g_m : float Conductance to CO2 transfer between intercellular spaces and sites of carboxylation (mol m-2 s-1 bar-1) theta : float curvature of the light response (-) K_cat : float specific activitt of Rubisco (mol CO2 mol-1 Rubisco s-1) K_s : float constant of proportionality between N_s and Jmax (g N m2 s umol-1) MAXIMISE : logical if we want to minimise a function we need to return f(x), if we want to maximise a function we return -f(x) """ self.alpha = alpha self.beta = beta self.g_m = g_m self.theta = theta self.K_cat = K_cat self.K_s = K_s self.sign = 1.0 if MAXIMISE: self.sign = -1.0 else: self.sign = 1.0 def calc_photosynthesis(self, N_pools=None, par=None, Km=None, Rd=None, gamma_star=None, Ca=None, N_p=None, error=False): """ Leaf temperature is assumed to be constant = 25 deg C. Parameters ---------- N_pools : list of floats list containing the 5 N pools: N_c - amount of N in chlorophyll (mol m-2) N_s - amount of N in soluble protein other than Rubisco (mol m-2) N_e - amount of N in electron transport components (mol m-2) N_r - amount of N in Rubisco (mol m-2) N_other - amount of leaf N not involved in photosynthesis (mol m-2) par : float incident PAR (umol m-2 s-1) Km : float effective Michalis-Menten coefficent of Rubisco (umol mol-1) Rd : float rate of dark respiration (umol m2 s-1) gamma_star : float leaf CO2 compensation point (umol mol-1) Ca : float atmospheric CO2 concentration (umol mol-1) N_p : float total amount of leaf N to be distrubuted (mol m-2) Returns: -------- An : float Net leaf assimilation rate [umol m-2 s-1] """ # unpack...we need to pass as a list for the optimisation step (N_e, N_c, N_r, N_other) = N_pools # Leaf absorptance depends on the chlorophyll protein complexes (N_c) absorptance = self.calc_absorptance(N_c) # Max rate of electron transport (umol m-2 s-1) Jmax = self.calc_jmax(N_e, N_c) # Max rate of carboxylation velocity (umol m-2 s-1) Vcmax = self.calc_vcmax(N_r) N_s = self.calc_n_alloc_soluble_protein(Jmax) # store values so I can get them all later self.N_pool_store = np.array([N_e, N_c, N_r, N_s, N_other]) # rate of electron transport, a saturating function of absorbed PAR J, error = self.quadratic(a=self.theta, b=-(self.alpha * absorptance * par + Jmax), c=self.alpha * absorptance * par * Jmax) if error: return None, error # CO2 concentration in the intercellular air spaces Ci = self.beta * Ca # Photosynthesis when Rubisco is limiting a = 1. / self.g_m b = (Rd - Vcmax) / self.g_m - Ci - Km c = Vcmax * (Ci - gamma_star) - Rd * (Ci + Km) Wc, error = self.quadratic(a=a, b=b, c=c) if error: return None, error # Photosynthesis when electron transport is limiting a = -4. / self.g_m b = (Rd - J / 4.0) / self.g_m - Ci - 2.0 * gamma_star c = J / 4.0 * (Ci - gamma_star) - Rd * (Ci + 2.0 * gamma_star) Wj, error = self.quadratic(a=a, b=b, c=c) if error: return None, error Aleaf = np.minimum(Wc, Wj) Cc = Ci - (Aleaf / self.g_m); An = (1.0 - (gamma_star / Cc)) * Aleaf - Rd return An, error def assim(self, Cc, a1, a2): """Photosynthesis with the limitation defined by the variables passed as a1 and a2, i.e. if we are calculating vcmax or jmax limited. Parameters: ---------- Cc : float CO2 concentration at the site of of carboxylation gamma_star : float CO2 compensation point in the abscence of mitochondrial respiration a1 : float variable depends on whether the calculation is light or rubisco limited. a2 : float variable depends on whether the calculation is light or rubisco limited. Returns: ------- assimilation_rate : float assimilation rate assuming either light or rubisco limitation. """ return a1 * (Cc / (Cc + a2)) def quadratic(self, a=None, b=None, c=None, error=False): """ minimilist quadratic solution as root for J solution should always be positive, so I have excluded other quadratic solution steps. If roots for J are negative or equal to zero then the data is bogus Parameters: ---------- a : float co-efficient b : float co-efficient c : float co-efficient Returns: ------- val : float positive root """ d = b**2 - 4.0 * a * c # discriminant if np.any(d <= 0.0) or np.any(np.isnan(d)): error = True root1 = -999.9 root2 = -999.9 else: root1 = (-b + np.sqrt(d)) / (2.0 * a) root2 = (-b - np.sqrt(d)) / (2.0 * a) return np.minimum(root1, root2), error def calc_jmax(self, N_e, N_c, a_j=15870.0, b_j=2775.0): """ Evans (1989( found a linear reln between the rate of electron transport and the total amount of nitrogen in the thylakoids per unit chlorophyll Parameters: ---------- N_e : float amount of N in electron transport components (mol m-2) N_c : float amount of N in chlorophyll (mol m-2) a_j : float umol (mol N)-1 s-1 b_j : float umol (mol N)-1 s-1 Returns: -------- Jmax : float Max rate of electron transport (umol m-2 s-1) """ Jmax = a_j * N_e + b_j * N_c return Jmax def calc_vcmax(self, N_r): """ Calculate Rubisco activity Parameters: ---------- N_r : float amount of N in Rubisco (mol m-2) Returns: -------- Vcmax : float Max rate of carboxylation velocity (umol m-2 s-1) """ conv = 7000.0 / 44.0 # mol N to mol Rubiso return self.K_cat * conv * N_r def calc_absorptance(self, N_c): """ Calculate absorptance in the leaf based on N in the chlorophyll protein complexes Parameters: ---------- N_c : float amount of N in chlorophyll (mol m-2) Returns: -------- absorptance : float Leaf absorptance (-) """ return (25.0 * N_c) / (25.0 * N_c + 0.076) def calc_n_alloc_soluble_protein(self, Jmax): """ Calculate N allocated to soluble protein Parameters: -------- Jmax : float Max rate of electron transport (umol m-2 s-1) Returns: -------- N_s : float amount of N in soluble protein other than Rubisco (mol m-2) """ return self.K_s * Jmax def optimise_func(self, x0, par=None, Km=None, Rd=None, gamma_star=None, Ca=None, N_p=None): """ Minimisation and maximisation are equivalent, to get the global maximium we just need to return -f(x) instead of f(x) """ An, error = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star, Ca, N_p) # I want only positive values for xopt if np.any(x0<0.0): return np.nan (N_e, N_c, N_r, N_other) = x0 # Leaf absorptance depends on the chlorophyll protein complexes (N_c) absorptance = self.calc_absorptance(N_c) # Max rate of electron transport (umol m-2 s-1) Jmax = self.calc_jmax(N_e, N_c) # Max rate of carboxylation velocity (umol m-2 s-1) Vcmax = self.calc_vcmax(N_r) N_s = self.calc_n_alloc_soluble_protein(Jmax) # xopt can't be > N_p if (N_s + N_e + N_r + N_c + N_other) > N_p: return np.nan An, err = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star, Ca, N_p) if err: return np.nan else: return self.sign * np.mean(An) def constraint1(self, x, *args): """ N_s + N_e + N_r + N_c < N_p returns a positive number if within bound and 0.0 it is exactly on the edge of the bound """ (par, Km, Rd, gamma_star, Ca, N_p) = args # unpack...we need to pass as a list for the optimisation step (N_e, N_c, N_r, N_other) = x # Leaf absorptance depends on the chlorophyll protein complexes (N_c) absorptance = self.calc_absorptance(N_c) # Max rate of electron transport (umol m-2 s-1) Jmax = self.calc_jmax(N_e, N_c) # Max rate of carboxylation velocity (umol m-2 s-1) Vcmax = self.calc_vcmax(N_r) N_s = self.calc_n_alloc_soluble_protein(Jmax) return N_p - (N_s + N_e + N_r + N_c + N_other) def constraint2(self, x, *args): """ N_e > 0.0 returns a positive number if within bound and 0.0 it is exactly on the edge of the bound """ return x[0] def constraint3(self, x, *args): """ N_c > 0.0 returns a positive number if within bound and 0.0 it is exactly on the edge of the bound """ return x[1] def constraint4(self, x, *args): """ N_r > 0.0 returns a positive number if within bound and 0.0 it is exactly on the edge of the bound """ return x[2] def constraint5(self, x, *args): """ N_other > 0.0 returns a positive number if within bound and 0.0 it is exactly on the edge of the bound """ return x[3] def optimise_test(): # default values Km = 544. Rd = 0.5 gamma_star = 38.6 par = np.arange(10.0, 1500.0, 20.0) Ca = np.ones(len(par)) * 350.0 P = PhotosynthesisModel(MAXIMISE=True) func = P.optimise_func # short name for func # initial guess # Testing, just equally divide up the total N available N_p = 0.09 x0 = np.array([N_p/5.0, N_p/5.0, N_p/5.0, N_p/5.0]) args = (par, Km, Rd, gamma_star, Ca, N_p) cons = ({'type': 'ineq', 'fun': P.constraint1, 'args': args}, {'type': 'ineq', 'fun': P.constraint2, 'args': args}, {'type': 'ineq', 'fun': P.constraint3, 'args': args}, {'type': 'ineq', 'fun': P.constraint4, 'args': args}, {'type': 'ineq', 'fun': P.constraint5, 'args': args}) result = scipy.optimize.minimize(func, x0, method="COBYLA", constraints=cons, args=args) if result["success"]: print result.x if __name__ == "__main__": optimise_test() -- View this message in context: http://scipy-user.10969.n7.nabble.com/maximising-a-function-tp17980.html Sent from the Scipy-User mailing list archive at Nabble.com.

1 1