SLIDEREAD

PYTHON:BATTERIES INCLUDEDUsing Python to SolvePartial Differential EquationsThis article describes two Python modules for solving partial differential equations (PDEs):PyCC is designed as a Matlab-like environment for writing algorithms for solving PDEs, andSyFi creates matrices based on symbolic mathematics, code generation, and the finiteelement method.Our work at the Simula Research Laboratory mostly focuses on computationalapplications in life sciences. Usually,this involves fairly typical partial differential equations such as the incompressible NavierStokes equations, elasticity equations, and parabolicand elliptic PDEs, but these PDEs are typically coupled either with each other or with ordinary differential equations (ODEs). Hence, even though thePDEs themselves are reasonably well understood,the couplings between them make the problems westudy quite challenging.Our design goals are therefore threefold. First,we want to easily define systems of PDEs. Second,we want it to be easy to play with different solutionalgorithms for systems of coupled PDEs. Finally,we want to reuse existing software to avoid reinventing the wheel.We use many good and mature libraries from theWeb, including Dolfin,GiNaC, MayaVi, NumPy,1521-9615/07/ 25.00 IEEECopublished by the IEEE CS and the AIPKENT-ANDRE MARDAL, OLA SKAVHAUG, GLENN T. LINES,GUNNAR A. STAFF, AND ÅSMUND ØDEGÅRDSimula Research Laboratory48THIS ARTICLE HAS BEEN PEER-REVIEWED.PETSc, SciPy, Trilinos, and VTK. In fact, we’remixing these libraries with our own packages: Famms (verification based on the method ofmanufactured solutions), Instant (inlining of C in Python), PyCC (the underlying framework forgluing components together), PySE (a finite differencetoolbox), Swiginac (a Pythoninterface to the symbolic mathematics engineGiNaC), and SyFi (a finite elementtoolbox).Some of these packages are Python modules,whereas the others—thanks to Python’s popularity inscientific computing—are equipped with Python interfaces. By using Python, we don’t have to mix thesepackages at the C level, which is a huge advantage.Solving Systems of PDEsCurrently, our most important application is in cardiac electrophysiology.1 The central model here isthe bidomain model,2 which is a system of two PDEsCOMPUTING IN SCIENCE & ENGINEERING

with the following form: v ( M i v ) ( M i u ) I ion ( v, s ) , t(1)0 (Mi v) ((Mi Me) u).(2)(The domain here is the same for both PDEs—thatis, the heart—but it has two potentials, intra- andextra-cellular, which live inside and outside heartcells.) The primary unknowns here are the transmembrane potential v and the extra cellular potential u. The function Iion(v, s) describes the flow ofions across the cell membrane and can be quitecomplicated; Mi represents intracellular conductivity, and Me is extracellular. The second argumentIion(v, s) is generally a vector of variables, governedby a set of ODEs: s F ( s , t ). t(3)Note that s s(x), so the ODE system is defined ineach point (you can find examples of cell models). Hence, we must solve the systemfor each computational node.The bidomain formulation gives an accurate description of the myocardial tissue’s electrical conduction. Coupled with realistic ODE models of theionic current, we can study a large range of phenomena, including conduction abnormalities due toischmeia or channel myopathies (genetic defects),fibrillation and defibrillation, and drug intervention.Figure 1, for example, shows a geometric model ofthe human heart’s ventricles. The color representsthe transmembrane potential’s magnitude; Figure1a shows normal activation, and Figure 1b showschaotic behavior (which corresponds to a fibrillatory heart with critically reduced pumping ability).We solve the bidomain model in Equations 1through 3 by using an operator-splitting approach,in which we first solve the ODE systems in eachcomputational node at each time step before wesolve the PDE system. Here’s a simple Pythonscript we use for solving this problem:from dolfin import Meshfrom pycc.MatSparse import *import numpyfrom pycc import MatFacfrom pycc import ConjGradfrom pycc.BlockMatrix import *from pycc.Functions import *from pycc.ODESystem import *from pycc.CondGen import *from pycc.IonicODEs import *MAY/JUNE (a)(b)Figure 1. Human heart ventricles. This model compares (a) normalactivity and (b) chaotic behavior.mesh Mesh("Heart.xml.gz")matfac MatFac.MatrixFactory(mesh)M PyCond("Heart.axis")pc.setconductances(3.0e-3, 3e-4)ct ConductivityTensorFunction(pc.conductivity)Ai s(5.0e-3, 1.6e-3)ct ConductivityTensorFunction(pc.conductivity) Construct compound matricesdt 0.1A M dt*AiB dt*AiBt dt*AiC dt*Aie# Create the Block systemAA BlockMatrix((A,B),(Bt,C))prec DiagBlockMatrix((MLPrec(A),MLPrec(C)))v numpy.zeros(A.n, dtype 'd') - 45.0u numpy.zeros(A.n, dtype 'd')x BlockVector(v,u)# Create one ODE systems for each vertexodesys Courtemanche ODESystem()ode solver RKF32(odesys)ionic IonicODEs(A.n, ode solver,odesys)49

)# Solvez numpy.zeros((A.n,), dtype 'd')for i in xrange(0, 10):t i*dtionic.forward(x[0], t, dt)ConjGrad.precondconjgrad(prec, AA,x, BlockVector(M*x[0], z))Although the code seems clean and simple, it’sdue to a powerful combination of C/C /Fortranand Python. The script runs on desktop computerswith meshes that have millions of nodes and cansolve complete problems within minutes or hours.All computer-intensive calculations such as computing matrices, solving linear systems (via algebraic multigrid and the conjugate gradientmethod), and solving ODE systems are done efficiently in C or C .Creating Matrices for Systems of PDEsWe created the tool SyFi to define finite elementsand variational forms, as well as generate C codefor finite element computations. It uses the symbolic mathematics engine GiNaC and its Pythoninterface Swiginac for all its basic mathematical operations. SyFi enables polynomial differentiationand integration on polygonal domains. Furthermore, it uses the computed expressions, such as entries in an element matrix, to generate C code.The following example demonstrates how tocompute an element matrix for the Jacobian of anincompressible power-law fluid’s (nonlinear) stationary Navier-Stokes equations. LetFi ( u u ) N i μ ( u ) u : N i dx ,Twheredef sum(u char,fe):ujs symbolic matrix(1,fe.nbf(),u char)u 0for j in range(0,fe.nbf()):u ujs.op(j)*fe.N(j)u u.evalm()return u, ujsnsd cvar.nsd 3polygon ReferenceTetrahedron()fe VectorCrouzeixRaviart(polygon,1)fe.set size(nsd) # size of basis functions()# create sum u i N iu, ujs sum("u", fe)n symbol("n")mu pow(inner(grad(u), grad(u)),n)for i in range(0,fe.nbf()):# nonlinear power-law diffusion termfi diffusion mu*inner(grad(u),grad(fe.N(i)))# nonlinear convection termuxgradu (u.transpose()*grad(u)).evalm()fi convection inner(uxgradu,fe.N(i), True)fi fi diffusion fi convectionFi polygon.integrate(fi)for j in range(0,fe.nbf()):# differentiate to get the Jacobianuj ujs.op(j)Jij diff(Fi, uj)print "J[%d,%d] %s\n"%(i,j,Jij.printc())u kukNk and (u) u 2n.Then,J ij Fi u j u j T ( u u ) N i μ( u ) u : N i dx ,(4)Here’s the corresponding Python code for computing Equation 4 and generating the C code:from swiginac import *from SyFi import *50Note that both the differentiation and integrationis performed symbolically exactly as we would havedone by hand. This naturally leads to quite efficientcode compared to the traditional way of implementing such integrals—namely, as quadratureloops that involve the evaluation of basis functions,their derivatives, and so on. The printc functiongenerates C code for the expressions; so far, we’veused this system to generate roughly 60,000 lines ofC code for computing various matrices based onvarious finite elements and variational forms.To ease the integration of the generated C code in Python, we developed an inlining tool calledCOMPUTING IN SCIENCE & ENGINEERING

Instant, which lets us generate code, generate thecorresponding wrapper code, compile and link it toan extension module, and then import the moduleon the fly. The following code demonstrates Instantwith a simple example, in which we computey( x ) sin(cos( x )) xsymbolically, generate the corresponding C code, and inline the expression in Python with x asa NumPy array:import swiginac as Sfrom Instant import inline with numpyimport numpy as Nx S.symbol("x")xi S.symbol("x[i]")f S.sin(S.cos(x))dfdx S.diff(f,x)print dfdxstring """void func (int n, double* x, int m,double* y) {if ( n ! m ) {printf("Both arrays should be ofthe same size!");return;}for (int i 0; i n; i ) {y[i] %s;}} """ % dfdx.subs( x xi ).printc()print stringfunc inline with numpy(string, arrays [['n', 'x'], ['m', 'y']])x N.arange(100.0 )y N.zeros(100, dtype 'd')func(x,y)print xprint yWMAY/JUNE e’ve shown that it’s possible tosolve real-life problems in a userfriendly environment by combining Python’s high-level syntax withthe efficiency of compiled languages. However, thisapproach opens up many new possibilities for combining symbolic mathematics and code generation,which is a largely overlooked alternative to traditionalapproaches in finite element simulations. We’re currently implementing fairly advanced finite elementmethods such as the mixed elasticity method,3 andwe also want to simulate human tissue and bloodwith the most realistic models available today.AcknowledgmentsMardal is supported by the Research Council of Norwayunder the grant ES254277.References1. J. Sundnes et al., Computing the Electrical Activity in the HumanHeart, Monographs in Computations Science and Engineering,Springer-Verlag.2. C.S. Henriquez, “Simulating the Electrical Behavior of CardiacTissue Using the Bidomain Model,” Crit. Rev. Biomedical Eng., vol.21, no. 1, 1993, pp. 1–77.3. D.N. Arnold, R.S. Falk, and R. Winther, “Finite Element ExteriorCalculus, Homological Techniques, and Applications,” Acta Numerica, , pp. 1–155.Kent-Andre Mardal is a postdoc at the Simula ResearchLaboratory. His research interests include high-level numerical programming and solution of PDEs. Mardal hasa PhD in scientific computing from the University of Oslo.Skavhaug is a research scientist at the Simula Research Laboratory. His research interests include high-levelnumerical programming, PDEs, and code verification.Skavhaug has a PhD in scientific computing from the University of Oslo. T. Lines is a research scientist at the Simula ResearchLaboratory. His research interests include PDEs and computer simulation of cardiac electrophysiology. Lines has aPhD in scientific computing from the University of Oslo.A. Staff is a consultant at Scandpower PetroleumTechnology. His research interests include software for scientific computing and initial value problems. Staff has aPhD in scientific computing from the University of Oslo.Åsmund Ødegård is an IT manager and part-time research scientist at the Simula Research Laboratory. Hisresearch interests include high-level languages in scientific computing, object-oriented numerics, PDEs, andhigh-level parallelism. Ødegård has a PhD in computational science from the University of Oslo.

PYTHON:BATTERIES INCLUDEDAnalysis of Functional MagneticResonance Imaging in PythonThe authors describe a package for analyzing magnetic resonance imaging (MRI) andfunctional MRI (fMRI) data, which is part of the Neuroimaging in Python (NIPY) project.An international group of leading statisticians, physicists, programmers, andneuroimaging methodologists are developing NIPY for wider use.Magnetic resonance imaging (MRI)measures induced magnetic properties of tissue. It has long been thechosen technique for creatinghigh-resolution anatomical images of the humanbrain. Over the past decade, a new technique calledfunctional MRI (fMRI) has become a powerful andwidely used method for studying human brainfunction. fMRI measures regional blood flowchanges in the brain, which can help researchersidentify the most active brain areas during mentaltasks such as memory and language.Functional MRI AnalysisThe first step of an fMRI analysis—image reconstruction—takes raw data from the scanner and performs a highly customized inverse Fourier transformto create a time series of 3D functional images. A typical next step is to estimate the movement betweenscans via an automated image-matching algorithmand then use that estimate to remove artifacts due to1521-9615/07/ 25.00 IEEECopublished by the IEEE CS and the AIPK. JARROD MILLMANUniversity of California, BerkeleyMATTHEW BRETTMRC Cognitive Brain Science Unit, Cambridge, UK52THIS ARTICLE HAS BEEN PEER-REVIEWED.motion. Researchers commonly relate the activitydetected in the low-resolution functional images toa high-resolution structural image of the same subject. However, to compare between subjects, thefunctional or structural data must be warped tomatch some standard brain, a process that requiressophisticated models of brain anatomy. Finally, statistical techniques are used to determine which brainregions are related to certain tasks or activities.Clearly, fMRI data analysis comes with severalchallenges. First, it has a wide variety of