3126.sagews

Hi there!

I am engineer dealing with nonlinear dynamical systems. I made a spkg for Assimulo so that I can reach the Sundials solvers, particularly IDA (see thread on sage-devel).

I am doing  models which depend on the modal-shapes of a structure to describe it's motion.

The model is built symbolically, and the size of the problem grows according to the number of modes used.

Everything is working fine, but too slow. I am looking for ways to speedup the evaluation/substitution on symbolic matrices and vectors.

It seems that fast_float and fast_callable is not available for these cases. In my opinion, this case Fortran seems much simpler than Cython.

My first idea was to translate the symbolic system of equation into Fortran and attach them back via f2py. The idea is to build a translator to compile symbolic entities which are to be evaluated numerically (intensively).

Here I made a simple comparison and illustration of a use case.

In this case simple case you see an improvement of at least 3 orders of magnitude.

Just to mention that Maxima has a f90 translator, but it is very limited. It pretty much does a print() to a file.

How can I build such a translator?

- The solver expects numpy.array at the input and output of user defined residual and jacobian functions.

- Floats must me decorated by the Fortran precision

- Nice if we can pass arrays of parameters instead of long lists of dummy variables

- Observe the 0-based Python to 1-based Fortran.

- Observe the 1**2 instead of 1^2.


 

# bring the system into state-space representation 

M = 2   # number of modes

# position and derivative
u =vector(SR, list(var('y_%d' %i) for i in range(M)) )
ud=vector(SR, list(var('yd_%d' %i) for i in range(M)) )

# velocity and derivative
v =vector(SR, list(var('y_%d' %i) for i in range(M,2*M)) )
vd=vector(SR, list(var('yd_%d' %i) for i in range(M,2*M)) )

# dummy state dependent polynomials (in reality depend on M)
p = u[0]^2+u[0]+1
q = u[1]^2+u[1]+1

# dummy state dependent mass matrix (also depend on M)
m = matrix(SR,M)
m[:,:] = [[0.25*p*q, 0.125*p*q],[0.125*p*q, 0.15*p*q]]

# system of equations (kind-of harmonic oscillator)
# 0 = udot - v
# 0 = M*vdot + c*v + k*u - f(t)

eq1 = ud - v 
eq2 = m * vd  # further terms...

# system as a vector
system = vector(SR,2*M)
system[0:M]= eq1
system[M: ]= eq2

var1SR = list(var('y_%d' %i)  for i in range(2*M))
var2SR = list(var('yd_%d' %i) for i in range(2*M))

# system as a list of fast callables
system_list = list(eq1)+list(eq2)
system_list = [fast_callable(eq,vars=var1SR+var2SR,domain=float) for eq in system_list]

# variables used on substitution
var1 = list('y_%d' %i  for i in range(2*M))
var2 = list('yd_%d' %i for i in range(2*M))

# Residual function needed by the solver (another is needed for Jacobian)
# need numpy.array as input and return

import numpy as NP

def residual(t,*y_and_yd):    
    res=NP.zeros(2*M,float)
    
    #yn =dict(zip(var1,y[:2*M]))
    #ydn=dict(zip(var2,yd[:2*M]))
    
    #state = dict(yn.items() + ydn.items())
    
    # evaluate system at current solver time-step and state
    res[:]=[eq(*y_and_yd) for eq in system_list]
    
    return NP.array(res,float)

# simple test

t=float(0.0)
y_and_yd=NP.ones(4*M,float)

a = residual(t,*y_and_yd); a

array([ 0. , 0. , 3.375, 2.475])

type(a)

<type 'numpy.ndarray'>

system.parent()

Vector space of dimension 4 over Symbolic Ring

I would like to do:

- Automate the translation of the sybolic system to a f2py equivalent.

%fortran
!f90
subroutine residual_fortran(t, y, yd, res, siz )
implicit none

integer, parameter :: DP = kind(1.0d0)

real (DP) :: t
real (DP) :: y(siz)          
real (DP) :: yd(siz)
real (DP) :: res(siz)
integer :: siz

intent (in)  :: t, y, yd, siz
intent (out) :: res

! # Pragmas f2py-specific
!f2py intent (in)  t, y, yd, siz
!f2py intent (out) res

res(1) = yd(1)- y(3)

res(2) = yd(2)- y(4)

res(3) = (y(2)**2 + y(2) + 1)*(0.125_DP*y(1)**2 + 0.125_DP*y(1) + &
& 0.125_DP)*yd(4) + (y(2)**2 + y(2) + 1)*(0.25_DP*y(1)**2 +       &
& 0.25_DP*y(1) + 0.25_DP)*yd(3)

res(4) = (y(2)**2 + y(2) + 1)*(0.125_DP*y(1)**2 + 0.125_DP*y(1) + &
& 0.125_DP)*yd(3) + (y(2)**2 + y(2) + 1)*(0.15_DP*y(1)**2 +       &
& 0.15_DP*y(1) + 0.15_DP)*yd(4)

end subroutine

None

print residual_fortran.__doc__

residual_fortran - Function signature: res = residual_fortran(t,y,yd,[siz]) Required arguments: t : input float y : input rank-1 array('d') with bounds (siz) yd : input rank-1 array('d') with bounds (siz) Optional arguments: siz := len(y) input int Return objects: res : rank-1 array('d') with bounds (siz)

y=NP.ones(2*M,float)
yd=NP.ones(2*M,float)

b= residual_fortran(t,y,yd); b

array([ 0. , 0. , 3.375, 2.475])

type(b)

<type 'numpy.ndarray'>

Comparing the implementations:

timeit('a=residual(t,*y_and_yd)')

625 loops, best of 3: 19 µs per loop

timeit('[eq(*y_and_yd) for eq in system_list]')

625 loops, best of 3: 8.04 µs per loop

timeit('NP.array([eq(*y_and_yd) for eq in system_list],float)')

625 loops, best of 3: 15.7 µs per loop

timeit('a=residual_fortran(t,y,yd)')

625 loops, best of 3: 662 ns per loop