On Wed, 30 Jul 2008 12:43:12 +0200 Johann Rohwer <jr@sun.ac.za> wrote:

On Wednesday, 30 July 2008, Nils Wagner wrote:

...

Hi all,

I am looking for a python tool to solve a set of DAE's, e.g.

m \ddot{x} = 2 \lambda x m \ddot{y} = 2 \lambda y x^2 + y^2 - l^2 = 0

Our group has developed PySUNDIALS (http://pysundials.sourceforge.net), which is a Python package providing bindings for the SUNDIALS suite of solvers (using ctypes). SUNDIALS has a DAE solver (IDA), which should be able to handle your kind of problem. The PySUNDIALS distributrion comes with examples for each of the wrapped solvers.

Johann

Hi Johann, I have installed pysundials. Can you give me an advice how to implement my simple example using pysundials ? idadenx.py seems to be a good starting point. How do I define the set of equations, initial conditions etc. ? Is it necessary to define the Jacobian ? Thank you in advance. Nils

On Wed, 30 Jul 2008 11:28:50 -0400 "Rob Clewley" <rob.clewley@gmail.com> wrote:

Hi,

PyDSTool wraps Hairer and Wanner's Radau 5 code using SWIG. This supports DAEs. An example is in PyDSTool/tests/DAE_example.py, and you get the advantage of the higher level Pointset and Trajectory structures for your solutions, with named coordinates (rather than plain arrays with integer indices), if that appeals to you.

-Rob

Hi Rob, Thank you for your reply. How do I install PyDSTool ? As opposed to other projects like numpy/scipy/matplotlib I can't find a setup.py. /usr/bin/python -i src/PyDSTool/tests/DAE_example.py Traceback (most recent call last): File "src/PyDSTool/tests/DAE_example.py", line 12, in ? from PyDSTool import * ImportError: No module named PyDSTool Nils

On Wed, 30 Jul 2008 11:59:25 -0400 "Rob Clewley" <rob.clewley@gmail.com> wrote:

Erm... yes, there isn't one! I know it's a bit clunky but I have been waiting to add a setup.py until the next major release because I wanted to get some details right about pre-compiling the C and Fortran extensions. Those issues have been ironed out so it shouldn't be long now...

Anyway, it's still very easy. The instructions at http://www.cam.cornell.edu/~rclewley/cgi-bin/moin.cgi/GettingStarted will fill you in but you basically just put the unzipped directory wherever you like and add it to your pythonpath. There is currently no pre-compiling of the extensions, it will all be done when you first invoke the compiler to build your Radau model. Let me know if you have any other questions or if there's something I should add to the documentation to make things more clear.

Cheers Rob

Hi Rob, I think I followed all these instructions. echo $PYTHONPATH /home/nwagner/:/home/nwagner/src/PyDSTool/:/home/nwagner/src/PyDSTool/tests/:/usr/lib/python2.4/site-packages/ nwagner@linux:~> /usr/bin/python -i src/PyDSTool/tests/DAE_example.py Traceback (most recent call last): File "src/PyDSTool/tests/DAE_example.py", line 12, in ? from PyDSTool import * ImportError: No module named PyDSTool Am I missing something ? Cheers, Nils

On Wed, 30 Jul 2008 12:17:17 -0400 "Rob Clewley" <rob.clewley@gmail.com> wrote:

Did you also add a PyDSTool.pth into python's site-packages/ directory?

The file /usr/lib/python2.4/site-packages/PyDSTool.pth consists of two lines /home/nwagner/src/PyDSTool/ /home/nwagner/src/PyDSTool/tests/ Nils

On Wed, 30 Jul 2008 12:36:05 -0400 "Rob Clewley" <rob.clewley@gmail.com> wrote:

As I say in the wiki in the Windows install section (I must add it to the *nix section!) of GettingStarted, my experience is that you'll need to add

/home /home/wagner /home/wagner/src

before these other lines. I don't understand why this is necessary from a technical perspective. It seems daft to me.

Strange, Now it works fine for me, but it's very slow... Thank you again Nils

On Wed, 30 Jul 2008 13:03:24 -0400 "Rob Clewley" <rob.clewley@gmail.com> wrote:

...

Strange, Now it works fine for me, but it's very slow...

I assume you mean the time to be able to run the radau test script?

Exactly,

If so, that's because the current version of the systems needs to recompile everything afresh each time, even though 90% of the code is the same. That's a distutils issue, but the version I'll be releasing in the next few months side-steps that problem.

Will you post an announcement to this list ? I am looking forward to the enhanced version version ! Is PyDSTool available via svn ?

Then, only your vector field will need to be re-compiled if you change it, and the core radau and dopri parts will be compiled once only at the original installation time. Hence the complicated issues getting setup.py working.

I think you'll find that once the model is compiled, it integrates solutions *very* quickly. Radau is very efficient, and there's only a small overhead in simulating models with it using the additional structures from PyDSTool. If this isn't what you find then there's something else wrong!

-Rob

Nils

Please, correct me if I am wrong. Is PyDSTool a package which aims at providing the functionalities of, say, XPP (http://www.math.pitt.edu/~bard/xpp/xpp.html)? That would be great.

In part, yes. It's all here: http://www.cam.cornell.edu/~rclewley/cgi-bin/moin.cgi/ProjectOverview and in the rest of the documentation on the site.

Folks, I have some rate/time data I'm in the process of analyzing. The data represent the natural gas produced from a well over time. The data are fit best by either an exponential model (least squares best fit of log of rate data, just qi and di see below) or a "hyperbolic" model: qt = qi *(1-b*di*t)**(-1.0/b) qt = rate at time t qi = initial rate (calculated) b = decline exponent (rate of change of the decline rate, calculated) di = initial decline rate (calculated) t = time The calculated parameters are found using the optimization routine scipy.optimize.fmin_tnc(). I'm already rejecting some best fit results because either the correlation coefficient is not statistically significant or the di (slope) is not significantly different from 0. The image I've attached shows two things. An exponential best fit is shown by the solid line and a hyperbolic best fit by the dotted line. It is common practice in analyzing this type of data to begin the best fit analysis at an early time point where the data actually start to "decline"; thus, the hyperbolic best fit started at time 2 and the first rate point (rate = about 1350) is ignored. The other thing the graph shows is that there appear to be outliers. My question is, what Python resources are there for testing for outliers and robust statistics? Is RANSAC appropriate for this? Note that when I run ransack.py from the cookbook, I get: Traceback (most recent call last): --Snip other traceback info--- File "C:\Documents and Settings\bnuttall\Desktop\ransac.py", line 132, in fit A = numpy.vstack([data[:,i] for i in self.input_columns]).T AttributeError: 'numpy.ndarray' object has no attribute 'T' I presume the T attribute is supposed to mean the transpose operator and have changed the offending statement(s, there are 4 total) to the form: A = numpy.vstack([data[:,i] for i in self.input_columns]).transpose() Thanks. Brandon Nuttall

Nuttall, Brandon C wrote:

Folks,

I have some rate/time data I'm in the process of analyzing. The data represent the natural gas produced from a well over time. The data are fit best by either an exponential model (least squares best fit of log of rate data, just qi and di see below) or a "hyperbolic" model:

qt = qi *(1-b*di*t)**(-1.0/b)

qt = rate at time t qi = initial rate (calculated) b = decline exponent (rate of change of the decline rate, calculated) di = initial decline rate (calculated) t = time

The calculated parameters are found using the optimization routine scipy.optimize.fmin_tnc(). I'm already rejecting some best fit results because either the correlation coefficient is not statistically significant or the di (slope) is not significantly different from 0.

The image I've attached shows two things. An exponential best fit is shown by the solid line and a hyperbolic best fit by the dotted line. It is common practice in analyzing this type of data to begin the best fit analysis at an early time point where the data actually start to "decline"; thus, the hyperbolic best fit started at time 2 and the first rate point (rate = about 1350) is ignored. The other thing the graph shows is that there appear to be outliers.

My question is, what Python resources are there for testing for outliers and robust statistics?

Is RANSAC appropriate for this? Note that when I run ransack.py from the cookbook, I get:

Traceback (most recent call last):

--Snip other traceback info--- File "C:\Documents and Settings\bnuttall\Desktop\ransac.py", line 132, in fit A = numpy.vstack([data[:,i] for i in self.input_columns]).T AttributeError: 'numpy.ndarray' object has no attribute 'T'

I presume the T attribute is supposed to mean the transpose operator and have changed the offending statement(s, there are 4 total) to the form:

A = numpy.vstack([data[:,i] for i in self.input_columns]).transpose()

Thanks.

Brandon Nuttall

------------------------------------------------------------------------

------------------------------------------------------------------------

_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user

Hi, First, what do you mean by 'outlier'? Please see, for example, http://en.wikipedia.org/wiki/Outlier because outliers fully depend on the (implicit or explicit) model fitted. Note this is really defined for linear models not nonlinear models that you are fitting. Second, NumPy provides the tools to do fit various models but you need to implement your own statistics especially for nonlinear models. (In particular I would strongly suggest looking at likelihood ratio test and Bayesian information criterion for model selection.) Really your best bet is to use to R directly or interface to R with Rpy. That way you also have access to many nonlinear statistical routines that can fit that type of data. If you do not need a parametric model, try fitting splines. These may be more suited to handle the variability near t=60 although that may be more related to data collection and time scale. Regards Bruce