Mailman 3 November 2011 - SciPy-User

Power Spectral Density in SciPy, not pylab
by Gus Ishere 26 Nov '11

26 Nov '11

I see the psd() function within matplotlib pylab, and I realize that I can use it to get an array of Pxx,freqs. I don't need to make a figure though and I need to take the PSD of many time series (so it would be nice to do it in a vectorized way). Is there a PSD function within SciPy that allows me to select an axis in a n-dimensional array on which to do the PSD estimation, and then return an n-1 dimensional array? Thanks, Gustavo

2 1

ANN: pandas 0.6.0 released
by Wes McKinney 26 Nov '11

26 Nov '11

I'm pleased to announce the pandas 0.6.0 major release. It's been about one month since the last major release. It includes 155 commits and 16 pull requests closing 78 tickets on GitHub. Several new people contributed code to the project for this release. This upgrade is recommended for all users and should not cause any API breakage for 0.5.0 users. There are a lot of miscellaneous new functions and features, many performance enhancements, and a significant number of bugs and corner cases encountered since the 0.5.0 release. See the full release notes below and on GitHub. Some features to look forward to (or help with!) in the next couple releases: - NumPy datetime64 type integration - Enhanced GroupBy, especially for binning time series data - Further performance enhancements to existing functionality Many thanks to all the users who contributed code, bug reports, and suggestions for new features. best, Wes What is it ========== pandas is a Python package providing fast, flexible, and expressive data structures designed to make working with “relational” or “labeled” data both easy and intuitive. It aims to be the fundamental high-level building block for doing practical, real world data analysis in Python. Additionally, it has the broader goal of becoming the most powerful and flexible open source data analysis / manipulation tool available in any language. Links ===== Release Notes: https://github.com/wesm/pandas/blob/master/RELEASE.rst Documentation: http://pandas.sourceforge.net Installers: http://pypi.python.org/pypi/pandas Code Repository: http://github.com/wesm/pandas Mailing List: http://groups.google.com/group/pystatsmodels Blog: http://blog.wesmckinney.com pandas 0.6.0 ============ **Release date:** 11/25/2011 **API Changes** - Arithmetic methods like `sum` will attempt to sum dtype=object values by default instead of excluding them (GH #382) **New features / modules** - Add `melt` function to `pandas.core.reshape` - Add `level` parameter to group by level in Series and DataFrame descriptive statistics (PR #313) - Add `head` and `tail` methods to Series, analogous to to DataFrame (PR #296) - Add `Series.isin` function which checks if each value is contained in a passed sequence (GH #289) - Add `float_format` option to `Series.to_string` - Add `skip_footer` (GH #291) and `converters` (GH #343) options to `read_csv` and `read_table` - Add proper, tested weighted least squares to standard and panel OLS (GH #303) - Add `drop_duplicates` and `duplicated` functions for removing duplicate DataFrame rows and checking for duplicate rows, respectively (GH #319) - Implement logical (boolean) operators &, |, ^ on DataFrame (GH #347) - Add `Series.mad`, mean absolute deviation, matching DataFrame - Add `QuarterEnd` DateOffset (PR #321) - Add matrix multiplication function `dot` to DataFrame (GH #65) - Add `orient` option to `Panel.from_dict` to ease creation of mixed-type Panels (GH #359, #301) - Add `DataFrame.from_dict` with similar `orient` option - Can now pass list of tuples or list of lists to `DataFrame.from_records` for fast conversion to DataFrame (GH #357) - Can pass multiple levels to groupby, e.g. `df.groupby(level=[0, 1])` (GH #103) - Can sort by multiple columns in `DataFrame.sort_index` (GH #92, PR #362) - Add fast `get_value` and `put_value` methods to DataFrame and micro-performance tweaks (GH #360) - Add `cov` instance methods to Series and DataFrame (GH #194, PR #362) - Add bar plot option to `DataFrame.plot` (PR #348) - Add `idxmin` and `idxmax` functions to Series and DataFrame for computing index labels achieving maximum and minimum values (PR #286) - Add `read_clipboard` function for parsing DataFrame from OS clipboard, should work across platforms (GH #300) - Add `nunique` function to Series for counting unique elements (GH #297) - DataFrame constructor will use Series name if no columns passed (GH #373) - Support regular expressions and longer delimiters in read_table/read_csv, but does not handle quoted strings yet (GH #364) - Add `DataFrame.to_html` for formatting DataFrame to HTML (PR #387) - MaskedArray can be passed to DataFrame constructor and masked values will be converted to NaN (PR #396) - Add `DataFrame.boxplot` function (GH #368, others) - Can pass extra args, kwds to DataFrame.apply (GH #376) **Improvements to existing features** - Raise more helpful exception if date parsing fails in DateRange (GH #298) - Vastly improved performance of GroupBy on axes with a MultiIndex (GH #299) - Print level names in hierarchical index in Series repr (GH #305) - Return DataFrame when performing GroupBy on selected column and as_index=False (GH #308) - Can pass vector to `on` argument in `DataFrame.join` (GH #312) - Don't show Series name if it's None in the repr, also omit length for short Series (GH #317) - Show legend by default in `DataFrame.plot`, add `legend` boolean flag (GH #324) - Significantly improved performance of `Series.order`, which also makes np.unique called on a Series faster (GH #327) - Faster cythonized count by level in Series and DataFrame (GH #341) - Raise exception if dateutil 2.0 installed on Python 2.x runtime (GH #346) - Significant GroupBy performance enhancement with multiple keys with many "empty" combinations - New Cython vectorized function `map_infer` speeds up `Series.apply` and `Series.map` significantly when passed elementwise Python function, motivated by PR #355 - Cythonized `cache_readonly`, resulting in substantial micro-performance enhancements throughout the codebase (GH #361) - Special Cython matrix iterator for applying arbitrary reduction operations with 3-5x better performance than `np.apply_along_axis` (GH #309) - Add `raw` option to `DataFrame.apply` for getting better performance when the passed function only requires an ndarray (GH #309) - Improve performance of `MultiIndex.from_tuples` - Can pass multiple levels to `stack` and `unstack` (GH #370) - Can pass multiple values columns to `pivot_table` (GH #381) - Can call `DataFrame.delevel` with standard Index with name set (GH #393) - Use Series name in GroupBy for result index (GH #363) - Refactor Series/DataFrame stat methods to use common set of NaN-friendly function - Handle NumPy scalar integers at C level in Cython conversion routines **Bug fixes** - Fix bug in `DataFrame.to_csv` when writing a DataFrame with an index name (GH #290) - DataFrame should clear its Series caches on consolidation, was causing "stale" Series to be returned in some corner cases (GH #304) - DataFrame constructor failed if a column had a list of tuples (GH #293) - Ensure that `Series.apply` always returns a Series and implement `Series.round` (GH #314) - Support boolean columns in Cythonized groupby functions (GH #315) - `DataFrame.describe` should not fail if there are no numeric columns, instead return categorical describe (GH #323) - Fixed bug which could cause columns to be printed in wrong order in `DataFrame.to_string` if specific list of columns passed (GH #325) - Fix legend plotting failure if DataFrame columns are integers (GH #326) - Shift start date back by one month for Yahoo! Finance API in pandas.io.data (GH #329) - Fix `DataFrame.join` failure on unconsolidated inputs (GH #331) - DataFrame.min/max will no longer fail on mixed-type DataFrame (GH #337) - Fix `read_csv` / `read_table` failure when passing list to index_col that is not in ascending order (GH #349) - Fix failure passing Int64Index to Index.union when both are monotonic - Fix error when passing SparseSeries to (dense) DataFrame constructor - Added missing bang at top of setup.py (GH #352) - Change `is_monotonic` on MultiIndex so it properly compares the tuples - Fix MultiIndex outer join logic (GH #351) - Set index name attribute with single-key groupby (GH #358) - Bug fix in reflexive binary addition in Series and DataFrame for non-commutative operations (like string concatenation) (GH #353) - setupegg.py will invoke Cython (GH #192) - Fix block consolidation bug after inserting column into MultiIndex (GH #366) - Fix bug in join operations between Index and Int64Index (GH #367) - Handle min_periods=0 case in moving window functions (GH #365) - Fixed corner cases in DataFrame.apply/pivot with empty DataFrame (GH #378) - Fixed repr exception when Series name is a tuple - Always return DateRange from `asfreq` (GH #390) - Pass level names to `swaplavel` (GH #379) - Don't lose index names in `MultiIndex.droplevel` (GH #394) - Infer more proper return type in `DataFrame.apply` when no columns or rows depending on whether the passed function is a reduction (GH #389) - Always return NA/NaN from Series.min/max and DataFrame.min/max when all of a row/column/values are NA (GH #384) - Enable partial setting with .ix / advanced indexing (GH #397) - Handle mixed-type DataFrames correctly in unstack, do not lose type information (GH #403) - Fix integer name formatting bug in Index.format and in Series.__repr__ - Handle label types other than string passed to groupby (GH #405) - Fix bug in .ix-based indexing with partial retrieval when a label is not contained in a level - Index name was not being pickled (GH #408) - Level name should be passed to result index in GroupBy.apply (GH #416) Thanks ------ - Craig Austin - Marius Cobzarenco - Joel Cross - Jeff Hammerbacher - Adam Klein - Thomas Kluyver - Jev Kuznetsov - Kieran O'Mahony - Wouter Overmeire - Nathan Pinger - Christian Prinoth - Skipper Seabold - Chang She - Ted Square - Aman Thakral - Chris Uga - Dieter Vandenbussche - carljv - rsamson

1 0

array of matrices
by Gustavo Goretkin 26 Nov '11

26 Nov '11

I want to create an array, call it R, of matrices [M1,M2, M3, ...] . So R is 3D. If v is a (column) vector, then np.dot(M1,v) matrix-multiples v by M1, and produces another (column) vector. I want R to be shaped so that np.dot(R,v) produces an array of column vectors [M1*v, M2*v, M3*v, ...] Specifically, I want to define a function that resembles: def rot_2d(theta): return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) that works even if theta itself is an array. np.dot(rot_2d(np.array([0,np.pi/2,np.pi]),np.array([1,0]) produces an array like [ [1,0], [0,1], [-1,0] ]. Thanks! Gustav

1 0

rvs and broadcasting
by josef.pktd＠gmail.com 26 Nov '11

26 Nov '11

rvs in scipy.stats distributions has a nasty broadcasting if location or scale are arrays and size is not defined for the same shape http://projects.scipy.org/scipy/ticket/1544 (also https://groups.google.com/group/pystatsmodels/browse_thread/thread/e757d73b… ) I was playing with two solutions while I was writing a rvs for the truncated normal. 1) broadcast shape parameters, loc and scale, if they are arrays produce rvs in that shape, and, if in this case size is not the same or 1, then raise a ValueError essentially lower, upper, loc, scale = np.broadcast_arrays(lower, upper, loc, scale) if (np.size(lower) > 1) and (size != (1,)) and (lower.shape != size): raise ValueError('Do you really want this? Then do it yourself.') 2) broadcast shape parameters, loc and scale, for each of these create random variables given by size, the return shape is essentially broadcasted shape concatenated with size, for example assert_equal(truncnorm_rvs(lower*np.arange(4)[:,None], upper, loc=np.arange(5), scale=1, size=(2,3)).shape, (4, 5, 2, 3)) this version is attached. Any opinions about which version should be preferred? (As aside, truncnorm and other distribution with parameter dependent support might also have other problems than just the broadcasting of shape and scale.) Josef

3 3

strange behavior calling odeint from brentq
by Schmidt, Phil 26 Nov '11

26 Nov '11

Hello, I am implementing the shooting method using optimize.brentq() and integrate.odeint(). The following is an outline of my code: def objective(t2, *args): t1 = args[0] x_init = args[1] x_target = args[2] x = odeint(dxdt, x_init, [t1, t2]) return x - x_target t_target = brentq(objective, t1, t2, args=(t1, x_init, x_target)) I have observed that if I place do-nothing statements in the objective function (e.g., print statements or dummy assignments like t1=t1), sometimes I will get different answers for t_target. I have not identified a pattern for when this may or may not occur, but presumably there is some dependency between brentq() and odeint(). I am running Scipy 0.9.0rc3, Python 2.6.5, Windows XP. Can anyone explain why this is happening, and point me to the "right" way to do what I'm attemtping? Thanks, Phil ______________________________________________________________________ This e-mail message may contain privileged and/or confidential information, and is intended to be received only by persons entitled to receive such information. If you have received this e-mail in error, please notify the sender immediately. Please delete it and all attachments from any servers, hard drives or any other media. Other use of this e-mail by you is strictly prohibited. All e-mails and attachments sent and received are subject to monitoring, reading and archival by Watlow. The recipient of this e-mail is solely responsible for checking for the presence of "Viruses" or other "Malware". Watlow accepts no liability for any damage caused by any such code transmitted by or accompanying this e-mail or any attachment. Please note that any views or opinions presented in this email are solely those of the author and do not necessarily represent those of Watlow.

2 1

ANN: Veusz 1.14, a scientific plotting package
by Jeremy Sanders 25 Nov '11

25 Nov '11

Veusz 1.14 ---------- Velvet Ember Under Sky Zenith ----------------------------- http://home.gna.org/veusz/ Copyright (C) 2003-2011 Jeremy Sanders <jeremy(a)jeremysanders.net> and contributors. Licenced under the GPL (version 2 or greater). Veusz is a Qt4 based scientific plotting package. It is written in Python, using PyQt4 for display and user-interfaces, and numpy for handling the numeric data. Veusz is designed to produce publication-ready Postscript/PDF/SVG output. The user interface aims to be simple, consistent and powerful. Veusz provides a GUI, command line, embedding and scripting interface (based on Python) to its plotting facilities. It also allows for manipulation and editing of datasets. Data can be captured from external sources such as Internet sockets or other programs. Changes in 1.14: * Added interactive tutorial * Points in graphs can be colored depending on another dataset and the scale shown in a colorbar widget * Improved CSV import - better data type detection - locale-specific numeric and date formats - single/multiple/none header modes - option to skip lines at top of file - better handling of missing values * Data can be imported from clipboard * Substantially reduced size of output SVG files * In standard data import, descriptor can be left blank to generate dataset names colX * Axis plotting range can be interactively manipulated * If axis is in date-time format, show and allow the min and max values to be in date-time format * ImageFile widget can have image data embedded in document file * Fit widget can update the fit parameters and fit quality to a label widget * Allow editing of 2D datasets in data edit dialog * Add copy and paste dataset command to dataset browser context menu Minor and API changes: * Examples added to help menu * Picker shows date values as dates * Allow descriptor statement in standard data files after a comment character, e.g. "#descriptor x y" * Added some further color maps * Draw key symbols for vector field widget * Import plugin changes - Register classes rather than instances (backward compatibility is retained) - Plugins can return constants and functions (see Constant and Function types) - Add DatasetDateTime for returning date-time datasets * Custom definitions - Add RemoveCustom API to remove custom definitions - AddCustom API can specify order where custom definition is added * C++ code to speed up plotting points of different sizes / colors * Expand files by default in data navigator window * Select created datasets in data edit dialog * Tooltip wrapping used in data navigator window * Grid lines are dropped if they overlap with edge of graph Bug fixes * Fix initial extension in export dialog * Fix crash on hiding pages * Fixed validation for numeric values * Position of grid lines in perpendicular direction for non default positions * Catch errors in example import plugin * Fix crash for non existent key symbols * Fix crash when mismatch of dataset sizes when combining 1D datasets to make 2D dataset Features of package: * X-Y plots (with errorbars) * Line and function plots * Contour plots * Images (with colour mappings and colorbars) * Stepped plots (for histograms) * Bar graphs * Vector field plots * Box plots * Polar plots * Ternary plots * Plotting dates * Fitting functions to data * Stacked plots and arrays of plots * Plot keys * Plot labels * Shapes and arrows on plots * LaTeX-like formatting for text * EPS/PDF/PNG/SVG/EMF export * Scripting interface * Dataset creation/manipulation * Embed Veusz within other programs * Text, CSV, FITS, NPY/NPZ, QDP, binary and user-plugin importing * Data can be captured from external sources * User defined functions, constants and can import external Python functions * Plugin interface to allow user to write or load code to - import data using new formats - make new datasets, optionally linked to existing datasets - arbitrarily manipulate the document * Data picker * Interactive tutorial * Multithreaded rendering Requirements for source install: Python (2.4 or greater required) http://www.python.org/ Qt >= 4.4 (free edition) http://www.trolltech.com/products/qt/ PyQt >= 4.3 (SIP is required to be installed first) http://www.riverbankcomputing.co.uk/pyqt/ http://www.riverbankcomputing.co.uk/sip/ numpy >= 1.0 http://numpy.scipy.org/ Optional: Microsoft Core Fonts (recommended for nice output) http://corefonts.sourceforge.net/ PyFITS >= 1.1 (optional for FITS import) http://www.stsci.edu/resources/software_hardware/pyfits pyemf >= 2.0.0 (optional for EMF export) http://pyemf.sourceforge.net/ PyMinuit >= 1.1.2 (optional improved fitting) http://code.google.com/p/pyminuit/ For EMF and better SVG export, PyQt >= 4.6 or better is required, to fix a bug in the C++ wrapping For documentation on using Veusz, see the "Documents" directory. The manual is in PDF, HTML and text format (generated from docbook). The examples are also useful documentation. Please also see and contribute to the Veusz wiki: http://barmag.net/veusz-wiki/ Issues with the current version: * Some recent versions of PyQt/SIP will causes crashes when exporting SVG files. Update to 4.7.4 (if released) or a recent snapshot to solve this problem. If you enjoy using Veusz, we would love to hear from you. Please join the mailing lists at https://gna.org/mail/?group=veusz to discuss new features or if you'd like to contribute code. The latest code can always be found in the Git repository at https://github.com/jeremysanders/veusz.git.

2 2

Re: [SciPy-User] Problem with ODEINT
by Lofgren, Eric 24 Nov '11

24 Nov '11

I checked the parameters, and they seem to both be reasonable values for the system in question. What's especially confusing to be is if I change the code to something like: SIR = spi.odeint(eq_system, startPop, t_interval, hmin=1e-20) where I've set a smaller step-size, I still get the occasional error like this: lsoda-- warning..internal t (=r1) and h (=r2) are such that in the machine, t + h = t on the next step (h = step size). solver will continue anyway In above, R1 = 0.1000000000000E+01 R2 = 0.5798216400482E-81 lsoda-- at t(=r1) and step size h(=r2), the error test failed repeatedly or with abs(h) = hmin In above, R1 = 0.1000000000000E+01 R2 = 0.5798216400482E-81 Repeated error test failures (internal error). Run with full_output = 1 to get quantitative information. Which suggests the solver is blazing right past the minimum step size I just set. Eric On Nov 24, 2011, at 10:56 AM, <scipy-user-request(a)scipy.org> <scipy-user-request(a)scipy.org> wrote: > Date: Thu, 24 Nov 2011 13:33:57 -0200 > From: Flavio Coelho <fccoelho(a)gmail.com> > Subject: Re: [SciPy-User] Problem with ODEINT > To: SciPy Users List <scipy-user(a)scipy.org> > Message-ID: > <CAAbyzbQP0JV3fyKsH3XiisZAhr-tBTik1Z8HEwi=_VJWgqcnjw(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > Hi, > > I have seen this before with odeint. It is probably related to the > parameterization of your model, which is leading very large derivatives. > This is causing the numerical solver to try to reduce the step size (h) to > try to resolve the dynamics. h is apparently hitting its default lower > limit (something in the vicinity of 1E-85). > > You can try to set a lower hmin, when calling odeint, but it would be best > to go over your parameter values and check if they are reasonable. > > Disclaimer: Ihave not run your script nor checked your parameters, just > took a lok at the error messages. > > good luck, > > Fl?vio

1 0

Problem with ODEINT
by Lofgren, Eric 24 Nov '11

24 Nov '11

I'm attempting to solve a fairly simple compartmental SIR model (http://en.wikipedia.org/wiki/SIR_Model) using SciPy and particularly odeint . I've solved things like this before using this code, and they've always worked, so I'm frankly quite puzzled. The code is as follows: --- #Python implementation of continuous SIR model #Import Necessary Modules import numpy as np import scipy.integrate as spi #Solving the differential equation. Solves over t for initial conditions PopIn def eq_system(startState,t): '''Defining SIR System of Equations''' beta = startState[3] gamma = startState[4] #Creating an array of equations Eqs= np.zeros((3)) Eqs[0]= -beta * startState[0]*startState[1] Eqs[1]= beta * startState[0]*startState[1] - gamma*startState[1] Eqs[2]= gamma*startState[1] return Eqs def model_solve(t): '''Stores all model parameters, runs ODE solver and tracks results''' #Setting up how long the simulation will run t_start = 1 t_end = t t_step = 0.02 t_interval = np.arange(t_start,t_end, t_step) n_steps = (t_end-t_start)/t_step #Setting up initial population state #n_params is the number of parameters (beta and gamma in this case) S0 = 0.99999 I0 = 0.00001 R0 = 0.0 beta = 0.50 gamma = 1/10. n_params = 2 startPop = (S0, I0, R0, beta, gamma) #Create an array the size of the ODE solver output with the parameter values params = np.zeros((n_steps,n_params)) params[:,0] = beta params[:,1] = gamma timer = np.arange(n_steps).reshape(n_steps,1) SIR = spi.odeint(eq_system, startPop, t_interval) #Glue together ODE model output and parameter values in one big array output = hstack((timer,SIR,params)) return output def MonteCarlo_SIR(runs): holder = [ ] for i in range(runs): holder.append(model_solve(100)) results = np.hstack(holder) return results testing = MonteCarlo_SIR(10) print testing print testing.shape --- Ignore for a moment that there's absolutely no reason for a Monte Carlo simulation of a system made up of nothing but constants. And the poorly documented code - this is essentially a running musing right now, and wasn't intended to see the light of day until I hit this problem. When I run this code, sometimes it just works, and I get a nice tidy array of results. But sometimes, errors like the following crop up: lsoda-- warning..internal t (=r1) and h (=r2) are such that in the machine, t + h = t on the next step (h = step size). solver will continue anyway In above, R1 = 0.1000000000000E+01 R2 = 0.0000000000000E+00 intdy-- t (=r1) illegal In above message, R1 = 0.1020000000000E+01 t not in interval tcur - hu (= r1) to tcur (=r2) In above, R1 = 0.1000000000000E+01 R2 = 0.1000000000000E+01 intdy-- t (=r1) illegal In above message, R1 = 0.1040000000000E+01 t not in interval tcur - hu (= r1) to tcur (=r2) In above, R1 = 0.1000000000000E+01 R2 = 0.1000000000000E+01 lsoda-- trouble from intdy. itask = i1, tout = r1 In above message, I1 = 1 In above message, R1 = 0.1040000000000E+01 Illegal input detected (internal error). Run with full_output = 1 to get quantitative information. lsoda-- at t (=r1) and step size h (=r2), the corrector convergence failed repeatedly or with abs(h) = hmin In above, R1 = 0.6556352055692E+01 R2 = 0.2169078563087E-05 Repeated convergence failures (perhaps bad Jacobian or tolerances). Run with full_output = 1 to get quantitative information. lsoda-- warning..internal t (=r1) and h (=r2) are such that in the machine, t + h = t on the next step (h = step size). solver will continue anyway In above, R1 = 0.1000000000000E+01 R2 = 0.5798211925922E-81 lsoda-- warning..internal t (=r1) and h (=r2) are such that in the machine, t + h = t on the next step (h = step size). solver will continue anyway In above, R1 = 0.1000000000000E+01 R2 = 0.8614947121323E-85 ... lsoda-- warning..internal t (=r1) and h (=r2) are such that in the machine, t + h = t on the next step (h = step size). solver will continue anyway In above, R1 = 0.1000000000000E+01 R2 = 0.1722989424265E-83 lsoda-- above warning has been issued i1 times. it will not be issued again for this problem In above message, I1 = 10 I believe those three are are a full tour of the error messages that are cropping up. This definitely occurs a minority of the time, but it's common enough that 4 out of 10 runs of the code above produces at least one error message like that. It seems that the problem is the step size getting small enough that it's beyond the precision of the machine to deal with, but passing an argument of something like hmin = 0.0000001 or something that should be well within range doesn't seem to help. Any idea what's going on? The end goal of this is a somewhat more complex set of equations, and the parameters (i.e. beta and gamma) to be drawn from a distribution, but I'm somewhat concerned about the trouble the solver is seeming to have with what should be a pretty straightforward case. Any suggestions on how to fix this? Or is there another type of ODE solver I should be using? I confess my knowledge of how these things work is...limited. Thanks in advance, Eric

3 2

Integration over Voronoi cells
by Eraldo Pomponi 24 Nov '11

24 Nov '11

Hi folks, I'm working on the integration of a function like: K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1] over Voronoi cells (r is the distance from the point at which is associated the cell). I googled a lot and I found this two useful hints: http://stackoverflow.com/questions/5941113/looking-for-python-package-for-n… http://mathforum.org/kb/message.jspa?messageID=4963570&tstart=0 but I'm still not able to understand how I should do this integration. I have the function that returns the segments ( list of [(x_star,y_start),(x_sop,y,stop)] ) necessary to construct the Voronoi cells associated to a set of points. Could someone suggest how to proceed ? There's also a numerical problem connected with the integration due to the singularity in r==0. Could you suggest which is a reasonable stable integration method available in scipy that could handle the function [1]? Cheers, Eraldo

2 3

numerical integration with square root like singularity
by Daniel Platz 24 Nov '11

24 Nov '11

Hi, I am stuck with a problem in scipy. I want to numerically integrate a function with a square root like singularity at one end of the integration interval. The integral has the form \int_{0}^{1} f(x) * x / sqrt(1-x**2) dx or alternatively \int_{0}^{A} f(x) / sqrt(A-x) dx. Can the quad function in scipy deal with this kind of singularity. I tried to use the points argument of the quad function but I still get warning messages and do not know how much I can trust the results. Alternatively, I was wondering if I can implement a Chebyshev–Gauss quadrature myself to lift the singularity? Or is there a way to do this elegantly using scipy? I would be very glad about a short answer. Best regards, Daniel Platz

3 2