Mailman 3 March 2008 - SciPy-Dev

Python and self-validating algorithms
by Nils Wagner 17 Mar '08

17 Mar '08

Hi all, Is someone aware of a python package for self-validating algorithms (similar to INTLAB for Matlab) ? Cheers, Nils References: http://www.ti3.tu-harburg.de/rump/intlab/ http://www.cs.cas.cz/rohn/matlab/

1 0

scons build scripts in scipy trunk: merge is done
by David Cournapeau 17 Mar '08

17 Mar '08

Hi, Since I could not make svnmerge work, I manually merged my branch to build scipy with scons in the trunk, and the merge is done. Since I did it manually, I may have missed something; I quickly try it on linux, and it built, though. Again, do not use scons to build scipy for anything but testing. thanks, David

2 1

Dealing with stupid compilers: ifort on itanium cannot compile empty source files
by David Cournapeau 15 Mar '08

15 Mar '08

Hi, While helping someone to get scipy compiled with intel compilers on itanium, I encountered a serious problem: ifort cannot compile empty source files on itanium !! The files I got problems with are in ARPACK (two .f files are empty), but there may be other too, I don't know. Putting comments do not seem to work, there must be some statements... What would be the best solution in this case ? Is it ok to remove the files (since they are empty, this cannot be too bad, I guess, except if there is a need to track upstream library) cheers, David

2 1

pbm in numpy.linalg.polyfit and nupy.linalg polyfit
by Julien Le Sommer 13 Mar '08

13 Mar '08

Hi, I just want to mention a possible bug in numpy/lib/polynomial.py with numpy version > 1.0.4 (with python2.4). the modification introduced with revision 3338 of numpy by charris, (see http://www.scipy.org/scipy/numpy/log/trunk/numpy/lib/polynomial.py) has broken the possibility to use polyfit with two-dimensional arrays. IPython gives the following message : ------------------------------------------------------------------------------------------ 308 # scale returned coefficients 309 if scale != 0 : --> 310 c /= vander([scale], order)[0] 311 312 if full : ValueError: shape mismatch: objects cannot be broadcast to a single shape -------------------------------------------------------------------------------------------- with the following script : -------------------------------------------------------------------------------------------- import numpy from numpy.random import * from numpy import linalg #- create some arrays : n = 100 r = 0.5 a = 1.4 # 1D vector T = numpy.array(xrange(n),dtype=numpy.float32)*2/(n-1) -1. X = a*T +r*(rand(n)-0.5) # (n * 2) array aa = numpy.array([[0.7,1.2]]) rr = numpy.array([[0.4,0.6]]) TT = T.reshape((n,1))*numpy.ones((n,2)) XX = aa*TT +rr*(randn(n,2)-0.5) # (n * 3) array aaa = numpy.array([[0.7,1.2,1.]]) rrr = numpy.array([[0.4,0.6,0.2]]) TTT = T.reshape((n,1))*numpy.ones((n,3)) XXX = aaa*TTT +rrr*(randn(n,3)-0.5) #- try numpy.polyfit # works ok (ar,br) = numpy.polyfit(T,X,1) # works ok (ar,br) = numpy.polyfit(T,XX,1) # fails (ar,br) = numpy.polyfit(T,XXX,1) -- Julien Le Sommer, LEGI/MEOM, BP53, 38041 Grenoble cedex 9, France Tel : +33 4 76 82 50 65 Fax : +33 4 76 82 52 71 http://www.legi.hmg.inpg.fr/~lesommer/

3 2

returning string by f2py
by dmitrey 13 Mar '08

13 Mar '08

hi all, I have modified bvls.f file to make it good enough to be connected to scipy (add istop and msg - integer & string description of stop case). However, I always get same msg value that I pass to the routine (via f2py). Could anyone check the modified bvls.f file? http://scipy.org/scipy/scikits/attachment/wiki/OOFiles/bvls.f (Or mb something should be modified in Python code?) Thank you in advance, D.

2 4

Merging branch to build scipy with scons
by David Cournapeau 13 Mar '08

13 Mar '08

Hi, I would like to discuss about the possibility to merge the build_with_scons branch of scipy into the trunk. It would enable building scipy with numscons. The merge should be much easier than for numpy: it only add setup and scons script files, it does not touch any source file, and by default, the new code should not even be seen by the build system, so it is unlikely to cause any unattended consequences (that is, if you do not build scipy with numscons, because it still have some rather subtle issues with fortran linking). cheers, David

3 6

Help with f2py
by Russel Howe 13 Mar '08

13 Mar '08

I am attempting to add support for lsodi from odepack to scipy.integrate.ode. I have almost succeeded, but to be complete I would like to be able to recover values of ydot during integration. The documentation for lsodi recommends a call to intdy to return these values. How do I write a pyf file for this function? From the lsodi documentation: c call intdy (t, k, rwork(21), nyh, dky, iflag) c c the input parameters are.. c c t = value of independent variable where answers are desired c (normally the same as the t last returned by lsodi). c for valid results, t must lie between tcur - hu and tcur. c (see optional outputs for tcur and hu.) c k = integer order of the derivative desired. k must satisfy c 0 .le. k .le. nqcur, where nqcur is the current order c (see optional outputs). the capability corresponding c to k = 0, i.e. computing y(t), is already provided c by lsodi directly. since nqcur .ge. 1, the first c derivative dy/dt is always available with intdy. c rwork(21) = the base address of the history array yh. c nyh = column length of yh, equal to the initial value of neq. t, k, and nyh are simple but how do I pass the base address of an array? Further digging shows the the yh array has dimensions (neq, nqcur) both of which are accessible at call time but unknown at compile time. I use this information to construct an apropriate 2-d array to pass to intdy, but I don't see how to pass the full dimensions to f2py. Attached is my current attempt, plus a small test program. Russel Index: scipy/integrate/setup.py =================================================================== --- scipy/integrate/setup.py (revision 3969) +++ scipy/integrate/setup.py (working copy) @@ -51,7 +51,13 @@ sources=['vode.pyf'], libraries=libs, **newblas) + # lsodi + config.add_extension('lsodi', + sources=['lsodi.pyf'], + libraries=libs, + **newblas) + config.add_data_dir('tests') return config Index: scipy/integrate/ode.py =================================================================== --- scipy/integrate/ode.py (revision 3969) +++ scipy/integrate/ode.py (working copy) @@ -15,18 +15,18 @@ i = 0, ..., len(y0) - 1 -class ode ---------- - A generic interface class to numeric integrators. It has the following methods:: - integrator = ode(f,jac=None) + integrator = ode(f,jac=None,adda=None) integrator = integrator.set_integrator(name,**params) - integrator = integrator.set_initial_value(y0,t0=0.0) + integrator = integrator.set_initial_value(y0,t0=0.0,ydoti=None) integrator = integrator.set_f_params(*args) + integrator = integrator.set_adda_params(*args) integrator = integrator.set_jac_params(*args) y1 = integrator.integrate(t1,step=0,relax=0) + y1 = integrator.implicit_integrate(t1,step=0,relax=0) + ydot = integrator.get_ydot(t1,order=1) flag = integrator.successful() """ @@ -97,6 +97,51 @@ failures, and for this problem one should instead use DVODE on the equivalent real system (in the real and imaginary parts of y). +lsodi +~~~~ +Integrator for linearly implicit systems of first-order odes. lsodi +provides the same methods as vode (adams and bdf). This integrator +accepts the following parameters in set_integrator() method of the ode class: + atol=float|seq + rtol=float + lband=None|int + rband=None|int + method='adams'|'bdf' + with_jacobian=0|1 + nsteps = int + (first|min|max)_step = float + order = int # <=12 for adams, <=5 for bdf +Source: http://www.netlib.org/ode/lsodi.f + +Implicit integration requires two functions do be defined to integrate + d y(t)[i] + a(t,y) * --------- = g(t,y)[i] + dt +where a(t,y) is a square matrix. + +res returns an (neq, ) array and is passed to the ode constructor as f. + def res(t, y, s, ires): + r(t,y,s)=g(t,y)-a(t,y)*s + return r + +adda must modify the provided matrix p. If banded storage is used, +ml and mu provide the upper and lower diagonals, see the lsodi.f +source for full documentation. adda is passed to the ode constructor +with the adda= keyword argument. + def adda(t,y,ml,mu,p,[nrowp]): + p += a(t,y) + return p + +An optional jacobian can be specified which provides the derivatves of +r(t,y,s)=g(t,y)-a(t,y)*s + def jac(t, y, s, ml, mu, nrowp): + d r(t,y,s)[i] + p_{i,j} = -------------- + d y[j] + return p + +To integrate the system to tout use integrator.implicit_integrate(tout) + """ __doc__ += integrator_info @@ -142,11 +187,11 @@ # if myodeint.runner: # IntegratorBase.integrator_classes.append(myodeint) -__all__ = ['ode'] +__all__ = ['ode','lsodi'] __version__ = "$Id$" __docformat__ = "restructuredtext en" -from numpy import asarray, array, zeros, sin, int32, isscalar +from numpy import asarray, array, zeros, sin, int32, isscalar, zeros_like import re, sys #------------------------------------------------------------------------------ @@ -190,7 +235,7 @@ __doc__ += integrator_info - def __init__(self, f, jac=None): + def __init__(self, f, jac=None, adda=None): """ Define equation y' = f(y,t) where (optional) jac = df/dy. @@ -206,20 +251,28 @@ self.stiff = 0 self.f = f self.jac = jac + self.adda = adda self.f_params = () + self.adda_params=() self.jac_params = () self.y = [] - def set_initial_value(self, y, t=0.0): - """Set initial conditions y(t) = y.""" + def set_initial_value(self,y,t=0.0,ydot=None): + """Set initial conditions y(t) = y. In the implicit case, an initial ydot can optionally specified. If ydot is not provided and a(t0,y0) is not singular, lsodi will compute ydot0.""" if isscalar(y): y = [y] n_prev = len(self.y) + if ydot is not None: + self.ydot=array(ydot, float) + if len(y) != len(self.ydot): + raise ValueError, "ydot (len=%d) and y (len=%d) must have same length"%(len(y),len(ydot)) + else: + self.ydot=zeros_like(y) if not n_prev: self.set_integrator('') # find first available integrator self.y = asarray(y, self._integrator.scalar) self.t = t - self._integrator.reset(len(self.y),self.jac is not None) + self._integrator.reset(len(self.y),self.jac is not None, ydot is not None) return self def set_integrator(self, name, **integrator_params): @@ -238,6 +291,7 @@ print 'No integrator name match with %s or is not available.'\ %(`name`) else: + self.integrator_name=name self._integrator = integrator(**integrator_params) if not len(self.y): self.t = 0.0 @@ -258,6 +312,26 @@ self.f_params,self.jac_params) return self.y + def implicit_integrate(self,t,step=0,relax=0): + """Find y=y(t), set y as an initial condition, and return y.""" + if step and self._integrator.supports_implicit_step: + mth = self._integrator.implicit_step + elif relax and self._integrator.supports_implicit_run_relax: + mth = self._integrator.implicit_run_relax + else: + mth = self._integrator.implicit_run + self.y,self.t = mth(self.f, self.adda, self.jac or (lambda :None), + self.y, self.ydot, + self.t, t, + self.f_params, self.adda_params, self.jac_params) + return self.y + + def get_ydot(self, t, order=1): + if self._integrator.supports_intdy: + return self._integrator.intdy(t, order) + else: + raise NotImplementedError, "integrator %s does not support derivative calculation"%self.integrator_name + def successful(self): """Check if integration was successful.""" try: self._integrator @@ -269,6 +343,16 @@ self.f_params = args return self + def set_res_params(self,*args): + """Set extra-parameters for user-supplied function res.""" + self.f_params = args + return self + + def set_adda_params(self,*args): + """Set extra-parameters for user-supplied function adda.""" + self.adda_params = args + return self + def set_jac_params(self,*args): """Set extra parameters for user-supplied function jac.""" self.jac_params = args @@ -293,7 +377,7 @@ integrator_classes = [] scalar = float - def reset(self,n,has_jac): + def reset(self,n,has_jac,has_ydot=False): """Prepare integrator for call: allocate memory, set flags, etc. n - number of equations. has_jac - if user has supplied function for evaluating Jacobian. @@ -305,18 +389,38 @@ defines the stoppage coordinate of the result. """ raise NotImplementedError,\ - 'all integrators must define run(f,jac,t0,t1,y0,f_params,jac_params)' + 'all integrators must define run(f,jac,t0,t1,y0,f_params,jac_params) or' + ' implicit_run(self,res,adda,jac,y0,ydoti,t0,t1,f_params,adda_params,jac_params)' + def implicit_run(self,res,adda,jac,y0,ydoti,t0,t1,f_params,adda_params,jac_params): + """Integrate from t=t0 to t=t1 using y0 as an initial condition. + Return 2-tuple (y1,t1) where y1 is the result and t=t1 + defines the stoppage coordinate of the result. + """ + raise NotImplementedError,\ + 'all integrators must define implicit_run(self,res,adda,jac,y0,ydoti,t0,t1,f_params,adda_params,jac_params)'\ + ' or run(f,jac,t0,t1,y0,f_params,jac_params)' + def step(self,f,jac,y0,t0,t1,f_params,jac_params): """Make one integration step and return (y1,t1).""" raise NotImplementedError,'%s does not support step() method' %\ (self.__class__.__name__) + def implicit_step(self,res,adda,jac,y0,ydoti,t0,t1,f_params,adda_params,jac_params): + """Make one integration step and return (y1,t1).""" + raise NotImplementedError,'%s does not support implicit_step() method' %\ + (self.__class__.__name__) + def run_relax(self,f,jac,y0,t0,t1,f_params,jac_params): """Integrate from t=t0 to t>=t1 and return (y1,t).""" raise NotImplementedError,'%s does not support run_relax() method' %\ (self.__class__.__name__) + def implicit_run_relax(self,res,adda,jac,y0,ydoti,t0,t1,f_params,adda_params,jac_params): + """Integrate from t=t0 to t>=t1 and return (y1,t).""" + raise NotImplementedError,'%s does not support implicit_run_relax() method' %\ + (self.__class__.__name__) + #XXX: __str__ method for getting visual state of the integrator class vode(IntegratorBase): @@ -338,6 +442,9 @@ } supports_run_relax = 1 supports_step = 1 + supports_implicit_run_relax = 0 + supports_implicit_step = 0 + supports_intdy = 0 def __init__(self, method = 'adams', @@ -367,7 +474,7 @@ self.first_step = first_step self.success = 1 - def reset(self,n,has_jac): + def reset(self,n,has_jac,has_ydot=False): # Calculate parameters for Fortran subroutine dvode. if has_jac: if self.mu is None and self.ml is None: @@ -467,6 +574,9 @@ supports_run_relax = 1 supports_step = 1 scalar = complex + supports_implicit_run_relax=0 + supports_implicit_step=0 + supports_intdy=0 def reset(self, n, has_jac): # Calculate parameters for Fortran subroutine dvode. @@ -560,3 +670,248 @@ if zvode.runner: IntegratorBase.integrator_classes.append(zvode) + + +class lsodi(IntegratorBase): + try: + import lsodi as _lsodi + except ImportError: + print sys.exc_value + _lsodi = None + runner = getattr(_lsodi,'lsodi',None) + _intdy = getattr(_lsodi,'intdy',None) + + messages = {2 : 'lsodi was successful.', + -1 : 'excess work done on this call (check all inputs).', + -2 : 'excess accuracy requested (tolerances too small).', + -3 : 'illegal input detected (see printed message).', + -4 : 'repeated error test failures (check all inputs).', + -5 : 'repeated convergence failures (perhaps bad jacobian' + ' supplied or wrong choice of tolerances).', + -6 : 'error weight became zero during problem. (solution' + ' component i vanished, and atol or atol(i) = 0.).', + -7 : 'cannot occur in casual use.', + -8 : 'lsodi was unable to compute the initial dy/dt. In' + ' casual use, this means a(t,y) is initially singular.' + ' Supply ydoti and use istate = 1 on the first call.' + } + supports_run_relax=0 + supports_step=0 + scalar = float + supports_implicit_run_relax=1 + supports_implicit_step=1 + supports_intdy=1 + + def __init__(self, + method="adams", + with_jacobian=0, + banded = 0, + lband=None,uband=None, + rtol=1e-6,atol=[1e-12], + order = 0, + nsteps = 0, + max_step = 0.0, # corresponds to infinite + min_step = 0.0, + first_step = 0.0, # determined by solver + ): + + if re.match(method,r'adams',re.I): self.meth = 1 + elif re.match(method,r'bdf',re.I): self.meth = 2 + else: raise ValueError,'Unknown integration method %s'%(method) + self.with_jacobian = with_jacobian + self.rtol = rtol + self.atol = array(atol) + self.atol = atol + + self.order = order + self.nsteps = nsteps + self.max_step = max_step + self.min_step = min_step + self.first_step = first_step + self.mu=uband + self.ml=lband + self.success = 1 + + def reset(self,n,has_jac,has_ydot=False): + # Calculate parameters for Fortran subroutine lsodi. + self.neq=n + if has_jac: + if self.mu is None and self.ml is None: + miter = 1 + else: + if self.mu is None: self.mu = 0 + if self.ml is None: self.ml = 0 + miter = 4 + else: + if self.mu is None and self.ml is None: + miter = 2 + else: + if self.mu is None: self.mu = 0 + if self.ml is None: self.ml = 0 + miter = 5 + mf = 10*self.meth + miter + if mf in [11,12]: + lrw = 22 + 16*n + n**2 + elif mf in [14,15]: + lrw = 22 + 17*n + (2*ml +mu)*neq + elif mf in [21,22]: + lrw = 22 + 9*n + n**2 + elif mf in [24,25]: + lrw = 22 + 10*n + (2*ml +mu)*neq + else: + raise ValueError,'Unexpected mf=%s'%(mf) + liw = 20 + n + rwork = zeros((lrw,), float) + rwork[4] = self.first_step + rwork[5] = self.max_step + rwork[6] = self.min_step + self.rwork = rwork + iwork = zeros((liw,), int32) + if self.ml is not None: + iwork[0] = self.ml + if self.mu is not None: + iwork[1] = self.mu + iwork[4] = self.order + iwork[5] = self.nsteps + iwork[6] = 2 # mxhnil + self.iwork = iwork + if len(self.atol) == 1: + itol=1 + else: + itol=2 + itask=1 + # don't internally compute ydot if it is provided + if has_ydot: + istate=1 + else: + istate=0 + iopt=1 + self.call_args = [itol,self.rtol,self.atol,itask,istate,iopt,self.rwork,self.iwork,mf] + self.success = 1 + + def implicit_run(self,*args): + #y1,t,istate,self.rwork,self.iwork = self.runner(*(args[:7]+tuple(self.call_args)+args[7:])) + y1,t,istate = self.runner(*(args[:7]+tuple(self.call_args)+args[7:])) + self.call_args[4]=istate + #self.call_args[6]=self.rwork + #self.call_args[7]=self.iwork + if istate <0: + print 'lsodi:',self.messages.get(istate,'Unexpected istate=%s'%istate) + self.success = 0 + #else: + # self.call_args[3] = 2 # upgrade itask from 1 to 2 + return y1,t + + def implicit_step(self,*args): + itask = self.call_args[3] + self.call_args[3] = 2 + r = self.run(*args) + self.call_args[3] = itask + return r + + def implicit_run_relax(self,*args): + itask = self.call_args[3] + self.call_args[3] = 3 + r = self.run(*args) + self.call_args[3] = itask + return r + + def intdy(self, t, order): + """ + Compute derivative of specified order at t. Must be called after a successfull call + to lsodi with tcur-hu<=t<=tcur. + """ + nqcur=self.iwork[13] + yh=self.rwork[20:20+nqcur*self.neq] + yh.resize((self.neq, nqcur)) + dy, iflag = self._intdy(t, order, yh) + if iflag<0: + if iflag==-1: raise ValueError, "t=%s invalid in call to intdy"%t + if iflag==-2: raise ValueError, "order=%s invalid in call to intdy"%order + return dy + + +if lsodi.runner: + IntegratorBase.integrator_classes.append(lsodi) + +def test1(): + def f(t,y): + a = sin(6*t) + return y*y-a+y + + ode_runner = ode(f) + ode_runner.set_integrator('vode') + ode_runner.set_initial_value([0.1,0.11,.1]*10) + + while ode_runner.successful() and ode_runner.t < 50: + y1 = ode_runner.integrate(ode_runner.t+2) + print ode_runner.t,y1[:3] + +def test2(): + # Stiff problem. Requires analytic Jacobian. + def f(t,y): + ydot0 = -0.04*y[0] + 1e4*y[1]*y[2] + ydot2 = 3e7*y[1]*y[1] + ydot1 = -ydot0-ydot2 + return [ydot0,ydot1,ydot2] + def jac(t,y): + jc = [[-0.04,1e4*y[2] ,1e4*y[1]], + [0.04 ,-1e4*y[2]-6e7*y[1],-1e4*y[1]], + [0.0 ,6e7*y[1] ,0.0]] + return jc + r = ode(f,jac).set_integrator('vode', + rtol=1e-4, + atol=[1e-8,1e-14,1e-6], + method='bdf', + ) + r.set_initial_value([1,0,0]) + print 'At t=%s y=%s'%(r.t,r.y) + tout = 0.4 + for i in range(12): + r.integrate(tout) + print 'At t=%s y=%s'%(r.t,r.y) + tout *= 10 + + +def test3(): + def res(t,y,s,ires): + r=array([ -.04*y[0] + 1.0e4*y[1]*y[2] - s[0], + .04*y[0] - 1.0e4*y[1]*y[2] - 3.0e7*y[1]*y[1] - s[1], + y[0] + y[1] + y[2] - 1.0], float) + return r + def adda(t,y,ml,mu,p,nrowp): + p[0,0] += 1.0 + p[1,1] += 1.0 + return p + def jac(t,y,s,ml,mu,nrowp): + p=array([ + [-.040, 1.0e4*y[2], 1.0e4*y[1]], + [ .04, -1.0e4*y[2] - 6.0e7*y[1], -1.0e4*y[1]], + [ 1.0, 1.0, 1.0]], float) + return p + y0=array([1.0,0.0,0.0], float) + ydot0=array([-0.04,0.04,0], float) + r=ode(res, jac=jac, adda=adda).set_integrator('lsodi', + rtol=1e-4, + atol=array([1e-6, 1e-10, 1e-6]), + method='bdf', + with_jacobian=1 + ) + r.set_initial_value(y0, ydot=ydot0) + tout=0.4 + for i in range(12): + r.implicit_integrate(tout) + if not r.successful(): + print "Quitting on error" + break + print "at t=%0.4e, y=%s"%(r.t, r.y) + tout *= 10. + + +if __name__ == "__main__": + print 'Integrators available:',\ + ', '.join(map(lambda c:c.__name__, + IntegratorBase.integrator_classes)) + test1() + test2() + test3() Index: scipy/integrate/lsodi.pyf =================================================================== --- scipy/integrate/lsodi.pyf (revision 0) +++ scipy/integrate/lsodi.pyf (revision 0) @@ -0,0 +1,69 @@ +! -*- f90 -*- +! Note: the context of this file is case sensitive. + +python module lsodi__user__routines + interface lsodi_user_interface + subroutine res(neq,tn,y,s,r,ires) + integer intent(hide),depend(y) :: neq = len(y) + double precision intent(in) :: tn + double precision dimension(neq),intent(in) :: y + double precision dimension(neq),intent(in) :: s + double precision dimension(neq),intent(out):: r + integer intent(in) :: ires + end subroutine res + subroutine adda(neq,t,y,ml,mu,p,nrowp) ! in :lsodiuser: + integer intent(hide),depend(y) :: neq = len(y) + double precision intent(in) :: t + double precision dimension(neq),intent(in) :: y + integer intent(in) :: ml + integer intent(in) :: mu + double precision dimension(nrowp,neq),intent(in,out) :: p + integer optional,check(shape(p,0)==nrowp),depend(p) :: nrowp=shape(p,0) + end subroutine adda + subroutine jac(neq,t,y,s,ml,mu,p,nrowp) + integer intent(hide),depend(y) :: neq = len(y) + double precision intent(in) :: t + double precision dimension(neq),intent(in) :: y + double precision dimension(neq),intent(in) :: s + integer intent(in) :: ml + integer intent(in) :: mu + double precision dimension(nrowp,neq),intent(out) :: p + integer optional,check(shape(p,0)==nrowp),depend(p) :: nrowp=shape(p,0) + end subroutine jac + end interface lsodi_user_interface +end python module lsodi__user__routines + +python module lsodi + interface + subroutine lsodi(res,adda,jac,neq,y,ydoti,t,tout,itol,rtol,atol,itask,istate,iopt,rwork,lrw,iwork,liw,mf) + use lsodi__user__routines + external res + external adda + external jac + integer intent(hide),depend(y) :: neq = len(y) + double precision dimension(neq),intent(in,out) :: y + double precision dimension(neq),intent(in) :: ydoti + double precision intent(in,out) :: t + double precision intent(in) :: tout + integer intent(in) :: itol + double precision intent(in) :: rtol + double precision dimension(neq),intent(in) :: atol + integer intent(in) :: itask + integer intent(in,out) :: istate + integer intent(in) :: iopt + double precision dimension(lrw),intent(in) :: rwork + integer optional,check(len(rwork)>=lrw),depend(rwork),intent(hide) :: lrw=len(rwork) + integer dimension(liw),intent(in) :: iwork + integer optional,check(len(iwork)>=liw),depend(iwork),intent(hide) :: liw=len(iwork) + integer intent(in) :: mf + end subroutine lsodi + subroutine intdy(t, k, yh, nyh, dky, iflag) + double precision intent(in) :: t + integer intent(in) :: k + double precision dimension(nyh,-1),intent(in) :: yh + integer intent(hide),depend(yh) :: nyh=len(yh) + double precision intent(out),dimension(nyh) :: dky + integer intent(out) :: iflag + end subroutine intdy + end interface +end python module lsodi import numpy as N import pylab as P import scipy.integrate as SI """ This equation represents this circuit: ------------------------- | | | | /Rs0 /Rs1 /Rs2 | \ \ \ | / / / / Rload \ \ \ \ | | | / |C0 |C1 |C2 \ --- --- --- | --- --- --- | | | | | ------------------------- Positive current is defined as flowing into a capacitor. Rload * Sum(Qdot) + Q/C + Qdot * Rs = 0 The equation to be integrated is [[Rload+Rs[0] Rload Rload] [Rload Rload+Rs[1] Rload] Qdot = -Q/C [Rload Rload Rload+Rs[2]]] """ Rload=1.0 Rs=N.array([0.1,0.2,0.3]) A=N.ones((3,3))*Rload+N.diagflat(Rs) C=N.array([1.0,2.0,3.0]) def res(t, y, s, ires): return -y/C-N.dot(A, s) def adda(self,t,y,ml,mu,p): return p+A def jac(t,y,ml,mu,p,nrowp): return N.diagflat(1.0/C) t=N.linspace(0,10,300) Q=N.zeros((len(t),3)) Qd=N.zeros((len(t),3)) Q[0]=N.array([1.0,1.0,1.0]) r=SI.ode(res, jac=jac, adda=adda).set_integrator('lsodi', with_jacobian=1) r.set_initial_value(Q[0]) for i, tn in enumerate(t[1:]): r.implicit_integrate(tn) if not r.successful(): print "Quitting on error at i=%s, t=%s"%(i,tn) break Q[i+1]=r.y Qd[i+1]=r.get_ydot(tn) P.clf() P.subplot(211) P.plot(t,Q) P.subplot(212) P.plot(t,Qd) P.show()

3 4

Google Summer of Code Ideas
by Jarrod Millman 13 Mar '08

13 Mar '08

Hello, I have started a Google Summer of Code Ideas page: http://scipy.org/scipy/scipy/wiki/SummerofCodeIdeas Please feel free to add any ideas you have for a summer project especially if you would be interested in mentoring it. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

1 0

[GSoC 2008]Machine learning package in SciPy
by Stéfan van der Walt 12 Mar '08

12 Mar '08

On Tue, Mar 11, 2008 at 11:50 AM, Anton Slesarev <slesarev.anton(a)gmail.com> wrote: > Hi all, > > it might be a good idea to have a machine learning(ML) package in SciPy. As > I understand there are some ML code in SciKits, but it is in raw state? > > There are a lot of machine learning projects, with its own data format, > number of classifiers, feature selection algorithms and benchmarks. But if > you want to compare your own algorithm with some others, you should convert > your data format to input format of every tool you want to use and after > training, you should convert output format of each tools to the single > format to have facility to compare results(for example you want to see > common which features was selected by different tools). > > Now I'm analyzing different ML approaches for the special case of text > classification problem. I couldn't find ML framework appropriate for my > task. I've got two simple requirements for this framework. It should support > sparse data format and has at least svm classifier. For example, Orange [1] > is a vary good data mining project but has poor sparse format support. PyML > [2] has all needed features, but there are problems with installation on > different platforms and code design is not perfect. > > I believe that creation framework, which will be convenient for scientist to > integrate their algorithms to it, is a vary useful challenge. Scientists > often talk about standard machine learning software[3] and may be SciPy will > be appropriate platform for developing such software. > > I can write detailed proposal, but I want to see is it interesting for > someone? Any wishes and recommendations? > > 1. Orange http://magix.fri.uni-lj.si/orange/ > 2. PyML http://pyml.sourceforge.net/ > 3. The Need for Open Source Software in Machine Learning > http://www.jmlr.org/papers/volume8/sonnenburg07a/sonnenburg07a.pdf I also recently learned of Elefant, http://elefant.developer.nicta.com.au/ but haven't had a chance to investigate in more detail. Regards Stéfan

6 8

ANN: EuroSciPy 2008 Conference - Leipzig, Germany
by Travis Vaught 12 Mar '08

12 Mar '08

Greetings, We're pleased to announce the EuroSciPy 2008 Conference to be held in Leipzig, Germany on July 26-27, 2008. http://www.scipy.org/EuroSciPy2008 We are very excited to create a venue for the European community of users of the Python programming language in science. This conference will bring the presentations and collaboration that we've enjoyed at Caltech each year closer to home for many users of SciPy, NumPy and Python generally--with a similar focus and schedule. Call for Participation: ---------------------- If you are a scientist using Python for your computational work, we'd love to have you formally present your results, methods or experiences. To apply to present a talk at this year's EuroSciPy, please submit an abstract of your talk as a PDF, MS Word or plain text file to euroabstracts(a)scipy.org. The deadline for abstract submission is April 30, 2008. Papers and/or presentation slides are acceptable and are due by June 15, 2008. Presentations will be allotted 30 minutes. Registration: ------------ Registration will open April 1, 2008. The registration fee will be 100.00€ for early registrants and will increase to 150.00€ for late registration. Registration will include breakfast, snacks and lunch for Saturday and Sunday. Volunteers Welcome: ------------------ If you're interested in volunteering to help organize things, please email us at info(a)scipy.org.

1 0