efficiency
kent-and
Hi, I just implemented a small toolbox for finite element calculations
based on Lagrangian elements
and compared it to similar code using GiNaC. The difference in
efficiency is remarkable. My sympy
code takes about 20 seconds where corresponding GiNaC code is below
one second.
Any comments on this ?
The code is below:
fem.py:
-----------------------------------------------------------------------------------------
from sympy import *
x,y,z = symbols('xyz')
class ReferenceSimplex:
def __init__(self, nsd):
self.nsd = nsd
coords = []
if nsd <= 3:
coords = symbols('xyz')[:nsd]
else:
coords = []
for d in range(0,nsd):
coords.append(Symbol("x_%d" % d))
self.coords = coords
def integrate(self,f):
coords = self.coords
nsd = self.nsd
limit = 1
for p in coords:
limit -= p
intf = f
for d in range(0,nsd):
p = coords[d]
limit += p
intf = integrate(intf, (p, 0, limit))
return intf
def bernstein_space(order, nsd):
if nsd > 3:
raise RuntimeError("Bernstein only implemented in 1D, 2D, and
3D")
sum = 0
basis = []
coeff = []
if nsd == 1:
b1, b2 = x, 1-x
for o1 in range(0,order+1):
for o2 in range(0,order+1):
if o1 + o2 == order:
aij = Symbol("a_%d_%d" % (o1,o2))
sum += aij*binomial(order,o1)*pow(b1, o1)*pow(b2,
o2)
basis.append(binomial(order,o1)*pow(b1,
o1)*pow(b2, o2))
coeff.append(aij)
if nsd == 2:
b1, b2, b3 = x, y, 1-x-y
for o1 in range(0,order+1):
for o2 in range(0,order+1):
for o3 in range(0,order+1):
if o1 + o2 + o3 == order:
aij = Symbol("a_%d_%d_%d" % (o1,o2,o3))
fac = factorial(order)/
(factorial(o1)*factorial(o2)*factorial(o3))
sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3,
o3)
basis.append(fac*pow(b1, o1)*pow(b2,
o2)*pow(b3, o3))
coeff.append(aij)
if nsd == 3:
b1, b2, b3, b4 = x, y, z, 1-x-y-z
for o1 in range(0,order+1):
for o2 in range(0,order+1):
for o3 in range(0,order+1):
for o4 in range(0,order+1):
if o1 + o2 + o3 + o4 == order:
aij = Symbol("a_%d_%d_%d_%d" %
(o1,o2,o3,o4))
fac = factorial(order)/
(factorial(o1)*factorial(o2)*factorial(o3)*factorial(o4))
sum += aij*fac*pow(b1, o1)*pow(b2,
o2)*pow(b3, o3)*pow(b4, o4)
basis.append(fac*pow(b1, o1)*pow(b2,
o2)*pow(b3, o3)*pow(b4, o4))
coeff.append(aij)
return sum, coeff, basis
def create_point_set(order, nsd):
h = Rational(1,order)
set = []
if nsd == 1:
for i in range(0, order+1):
x = i*h
if x <= 1:
set.append((x,y))
if nsd == 2:
for i in range(0, order+1):
x = i*h
for j in range(0, order+1):
y = j*h
if x + y <= 1:
set.append((x,y))
if nsd == 3:
for i in range(0, order+1):
x = i*h
for j in range(0, order+1):
y = j*h
for k in range(0, order+1):
z = j*h
if x + y + z <= 1:
set.append((x,y,z))
return set
def create_matrix(equations, coeffs):
A = zeronm(len(equations), len(equations))
i = 0; j = 0
for j in range(0, len(coeffs)):
c = coeffs[j]
for i in range(0, len(equations)):
e = equations[i]
d, r = div(e, c)
A[i,j] = d
return A
class Lagrange:
def __init__(self,nsd, order):
self.nsd = nsd
self.order = order
self.compute_basis()
def nbf(self):
return len(self.N)
def compute_basis(self):
order = self.order
nsd = self.nsd
N = []
pol, coeffs, basis = bernstein_space(order, nsd)
points = create_point_set(order, nsd)
equations = []
for p in points:
ex = pol.subs(x, p[0])
if nsd > 1:
ex = ex.subs(y, p[1])
if nsd > 2:
ex = ex.subs(z, p[2])
equations.append(ex )
A = create_matrix(equations, coeffs)
Ainv = A.inv()
b = eye(len(equations))
xx = Ainv*b
for i in range(0,len(equations)):
Ni = pol
for j in range(0,len(coeffs)):
Ni = Ni.subs(coeffs[j], xx[j,i])
N.append(Ni)
self.N = N
sympy_test.py:
---------------------------------------------------------------------
from fem import *
t = ReferenceSimplex(2)
fe = Lagrange(2,2)
u = 0
#compute u = sum_i u_i N_i
us = []
for i in range(0, fe.nbf()):
ui = Symbol("u_%d" % i)
us.append(ui)
u += ui*fe.N[i]
J = zeronm(fe.nbf(), fe.nbf())
for i in range(0, fe.nbf()):
Fi = u*fe.N[i]
for j in range(0, fe.nbf()):
uj = us[j]
integrands = diff(Fi, uj)
J[j,i] = t.integrate(integrands)
print J
syfi_test.py:
--------------------------------------------------------------------------------------------
from swiginac import *
from SyFi import *
t = ReferenceTriangle()
fe = Lagrange(t,2)
u = 0
us = []
for i in range(0, fe.nbf()):
ui = symbol("u_%d" % i)
us.append(ui)
u += ui*fe.N(i)
J = matrix(fe.nbf(), fe.nbf())
for i in range(0, fe.nbf()):
Fi = u*fe.N(i)
for j in range(0, fe.nbf()):
uj = us[j]
integrands = diff(Fi, uj)
J[j,i] = t.integrate(integrands)
print J
Fredrik Johansson
>
>
> Hi, I just implemented a small toolbox for finite element calculations
> based on Lagrangian elements
> and compared it to similar code using GiNaC. The difference in
> efficiency is remarkable. My sympy
> code takes about 20 seconds where corresponding GiNaC code is below
> one second.
>
> Any comments on this ?
From profiling your code, it appears that nearly all the time is spent
in integrate().
This shows quite clearly that we need to optimize the integration code for
polynomials.
Fredrik
Ondrej Certik
> Hi, I just implemented a small toolbox for finite element calculations
> based on Lagrangian elements
> and compared it to similar code using GiNaC. The difference in
> efficiency is remarkable. My sympy
> code takes about 20 seconds where corresponding GiNaC code is below
> one second.
Thanks a lot for writing this. I am aware of SyFi [1] and I wanted to
to write something like you just wrote in sympy. For those who want to
know more about it, I suggest to read the tutorial [2], which explains
everything. I myself am using the libmesh [3] library for finite
elements, which is very good. But the SyFi's approach is very nice.
Did you try it on some real example? How does it compare to the
traditional approach, that libmesh is using? BTW, I am glad that SyFi
is using swiginac[4], that we wrote with Ola Skavhaug. I thought that
when I wrap ginac to python, I'll get an easy to use symbolic library,
but I realized this is not the way.
Can I add your example to the examples directory in SymPy? We use a BSD license.
> Any comments on this ?
As Fredrik has said, we need to improve the integration of
polynomials. See recently discussed issues for that:
http://code.google.com/p/sympy/issues/detail?id=461
http://code.google.com/p/sympy/issues/detail?id=463
But generally, Python is like 200 slower than C++ in my experience, so
this is what has to be expected.
Could you tell us what are the future plans with SyFi and also
Fenics[5] etc? Last time I tried that, it took quite a long time to
compile and it wasn't as versatile as libmesh (particularly I need to
assign boundary conditions and changing conductivity in the body,
which wasn't implemented at that time, but libmesh can do that
easily). Generally I like the fenics ideas, but it wasn't yet ready
for a production use. So I am curious what your own plans are, as I
would like to help if I like it, because good finite elements in
Python is exactly what I need.
Ondrej
[1] http://heim.ifi.uio.no/~kent-and/software/SyFi/doc/SyFi.html
[2] http://heim.ifi.uio.no/~kent-and/software/SyFi/doc/tutorial/tutorial.pdf
[3] http://libmesh.sourceforge.net/
[4] http://swiginac.berlios.de/
[5] http://www.fenics.org/wiki/FEniCS_Project
Fredrik Johansson
> But generally, Python is like 200 slower than C++ in my experience, so
> this is what has to be expected.
I'd say it's closer to 20x slower for most code. (In fact, the
Computer Language Benchmarks Game lists Python as 17x slower than C++
on average.) There gap is larger for number crunching, but symbolic
computing doesn't suffer so badly since operations like comparions and
slicing can be delegated to Python built-in functions.
Did anyone compare sympycore with a C/C++ computer algebra system yet?
IIRC we compared some operation to Maxima and sympycore was 7x or so
slower.
Fredrik
kent-and
On 25 Nov, 22:51, "Ondrej Certik" <ond...@certik.cz> wrote:
> Hi Kent,
>
> > Hi, I just implemented a small toolbox for finite element calculations
> > based on Lagrangian elements
> > and compared it to similar code using GiNaC. The difference in
> > efficiency is remarkable. My sympy
> > code takes about 20 seconds where corresponding GiNaC code is below
> > one second.
>
> Thanks a lot for writing this. I am aware of SyFi [1] and I wanted to
> to write something like you just wrote in sympy. For those who want to
> know more about it, I suggest to read the tutorial [2], which explains
> everything. I myself am using the libmesh [3] library for finite
> elements, which is very good. But the SyFi's approach is very nice.
> Did you try it on some real example? How does it compare to the
> traditional approach, that libmesh is using? BTW, I am glad that SyFi
> is using swiginac[4], that we wrote with Ola Skavhaug. I thought that
> when I wrap ginac to python, I'll get an easy to use symbolic library,
> but I realized this is not the way.
>
equations, Navier-Stokes etc.
But I only use SyFi for generating the C++ code for the finite
elements, element matrices etc.
The generated code is usually very fast, orders of magnitude faster
than traditional quadrature.
Its code and efficiency is similar to FFC. Assembly is typically done
in Dolfin/PyCC and solving
the linear systems is done in Hypre/Trilinos.
> Can I add your example to the examples directory in SymPy? We use a BSD license.
>
> > Any comments on this ?
>
> As Fredrik has said, we need to improve the integration of
> polynomials. See recently discussed issues for that:
>
> But generally, Python is like 200 slower than C++ in my experience, so
> this is what has to be expected.
>
> Could you tell us what are the future plans with SyFi and also
> Fenics[5] etc? Last time I tried that, it took quite a long time to
> compile and it wasn't as versatile as libmesh (particularly I need to
> assign boundary conditions and changing conductivity in the body,
> which wasn't implemented at that time, but libmesh can do that
> easily). Generally I like the fenics ideas, but it wasn't yet ready
> for a production use. So I am curious what your own plans are, as I
> would like to help if I like it, because good finite elements in
> Python is exactly what I need.
> d it wasn't as versatile as libmesh (particularly I need to
> assign boundary conditions and changing conductivity in the body,
> which wasn't implemented at that time, but libmesh can do that
> easily). Generally I like the fenics ideas, but it wasn't yet ready
> for a production use. So I am curious what your own plans are, as I
> would like to help if I like it, because good finite elements in
> Python is exactly what I need.
year. But this
is improving. The future plan of SyFi is to come up with a nice user
interface/language such
that SyFi will work as a compiler translating high-level finite
element python code
to efficient low-level C++ code.
The reason why I am looking at Sympy is that although I am fairly
happy with GiNaC, it has its issues.
Still, I think I'd prefer Sympy (except the efficiency).
Can sympy generate C++ code ?
Kent
kent-and
On 25 Nov, 22:26, "Fredrik Johansson" <fredrik.johans...@gmail.com>
wrote:
Kent
Ondrej Certik
> equations, Navier-Stokes etc.
> But I only use SyFi for generating the C++ code for the finite
> elements, element matrices etc.
> The generated code is usually very fast, orders of magnitude faster
> than traditional quadrature.
> Its code and efficiency is similar to FFC. Assembly is typically done
> in Dolfin/PyCC and solving
> the linear systems is done in Hypre/Trilinos.
Assembly is enough for me, I need to solve the eigenproblem anyway, so
I use blzpack/arpack/primme for that.
I'll try it, this sounds pretty exciting.
> Many of the Fenics projects have been in a state of flux for the last
> year. But this
> is improving. The future plan of SyFi is to come up with a nice user
> interface/language such
> that SyFi will work as a compiler translating high-level finite
> element python code
> to efficient low-level C++ code.
>
> The reason why I am looking at Sympy is that although I am fairly
> happy with GiNaC, it has its issues.
Could you be more specific what issues GiNaC has? The only problem
with ginac that I know of is that it is very difficult to create new
functions,
but as far as I know, SyFi just needs polynomials.
> Still, I think I'd prefer Sympy (except the efficiency).
>
> Can sympy generate C++ code ?
Not yet, but this is fairly easy to do. We can generate a Python code:
In [1]: e = sin(x)**4/y**(-2)
In [2]: e
Out[2]:
2 4
y *sin (x)
In [3]: Basic.set_repr_level(0)
Out[3]: 2
In [4]: e
Out[4]: Mul(Pow(Symbol('y'), Integer(2)), Pow(sin(Symbol('x')), Integer(4)))
So we can generate C++ code as well. But of course only a subset of
sympy features can be exported to C++. Mainly just an expression and a
way to substitute to that
expression. So the above should translate to something like this:
float x;
float y;
e = pow(sin(x), 4)*pow(y,2);
is that correct? If you could give me more info what exactly you need,
I'll implement that.
> > From profiling your code, it appears that nearly all the time is spent
> > in integrate().
> > This shows quite clearly that we need to optimize the integration code for
> > polynomials.
> >
> > Fredrik
>
> Ok thanks, is this easy to do ?
I think it should be fairly easy, that's why I wrote SymPy, so that
those things are easy. We just need to go to:
http://hg.sympy.org/sympy/file/c8da4f2b2324/sympy/integrals/integrals.py
line 122, and we use a general (slow) algorithm to integrate a lot of
functions. So we just check, if we got a polynomial, extract the
coefficients and integrate it. We can use functions from:
http://hg.sympy.org/sympy/file/132066377866/sympy/polynomials/base.py
There is a Polynomial class in there, so my bet is just to implement
.integrate() in that class, to be as efficient as possible and than
just fix ingetrals.py, to try to convert to polynomial if it is a
polynomial.
I'll try it in the evening, if it speeds things up.
Ondrej
Ondrej Certik
So indeed I sped things up by a factor of 10x. See:
http://code.google.com/p/sympy/issues/detail?id=470
for details and a patch. Now it needs to be polished and committed.
Your example now runs for 2.3s. How fast is ginac on this exactly?
Ondrej
kent-and
SyFi/GiNaC runs in 0.3s on my computer.
Kent
kent-and
the institute
I work at in the long run).
The more serious issue (now) is that GiNaC expands everything.
For instance, I'd like the derivative with respect to x of a function
f(x)
to be stored as an object Derivative(f,x), much like GiNaCs integrals
before eval_integ.
Or when doing G(x) = N(x)*M(x), G is not the expanded product of N and
M, but keeps
N and M.
I guess SymPy and GiNaC behave similarly in these matters, but I guess
it is easier to implement
such features in SymPy.
>
> > Still, I think I'd prefer Sympy (except the efficiency).
>
> > Can sympy generate C++ code ?
>
> Not yet, but this is fairly easy to do. We can generate a Python code:
>
> In [1]: e = sin(x)**4/y**(-2)
>
> In [2]: e
> Out[2]:
> 2 4
> y *sin (x)
>
> In [3]: Basic.set_repr_level(0)
> Out[3]: 2
>
> In [4]: e
> Out[4]: Mul(Pow(Symbol('y'), Integer(2)), Pow(sin(Symbol('x')), Integer(4)))
>
> So we can generate C++ code as well. But of course only a subset of
> sympy features can be exported to C++. Mainly just an expression and a
> way to substitute to that
> expression. So the above should translate to something like this:
>
> float x;
> float y;
> e = pow(sin(x), 4)*pow(y,2);
>
> is that correct? If you could give me more info what exactly you need,
> I'll implement that.
>
kent-and
>
> So indeed I sped things up by a factor of 10x. See:
>
> http://code.google.com/p/sympy/issues/detail?id=470
Kent
Ondrej Certik
> > with ginac that I know of is that it is very difficult to create new
> > functions,
> > but as far as I know, SyFi just needs polynomials.
>
> the institute
> I work at in the long run).
>
> The more serious issue (now) is that GiNaC expands everything.
> For instance, I'd like the derivative with respect to x of a function
> f(x)
> to be stored as an object Derivative(f,x), much like GiNaCs integrals
> before eval_integ.
http://code.google.com/p/sympy/issues/detail?id=476
Thanks for noticing. This used to work, but it isn't apparently now in
SymPy. We'll fix that.
> Or when doing G(x) = N(x)*M(x), G is not the expanded product of N and
> M, but keeps
> N and M.
Is this the issue:
http://code.google.com/p/sympycore/issues/detail?id=18
?
So you want:
In [4]: f = x+y
In [5]: g = x+z
In [6]: f*g
Out[6]: (x + y)*(x + z)
In [7]: f+g
Out[7]: y + z + 2*x
So [6] is ok, but you want [7] to be "x+y+x+z" instead? Is it because
of memory requirements (see that issue for more details)?
>
> I guess SymPy and GiNaC behave similarly in these matters, but I guess
> it is easier to implement
> such features in SymPy.
In theory yes. There is the speed issue though with SymPy, it depends
if we can cope with that.
> > float x;
> > float y;
> > e = pow(sin(x), 4)*pow(y,2);
> >
> > is that correct? If you could give me more info what exactly you need,
> > I'll implement that.
> >
>
> This is what I need, although I'd like double instead of float.
I created an issue for that:
http://code.google.com/p/sympy/issues/detail?id=475
feel free to discuss it in there.
> >
> > So indeed I sped things up by a factor of 10x. See:
> >
> > http://code.google.com/p/sympy/issues/detail?id=470
>
> How do I download this ?
you follow the link in the issue to this page:
http://groups.google.com/group/sympy-patches/browse_thread/thread/c51231bb0fd91c71
and download the integration.patch. Then
$ hg clone http://hg.sympy.org/sympy/
$ cd sympy/
$ hg import ~/Desktop/Downloads/integration.patch
applying /home/ondra/Desktop/Downloads/integration.patch
$ time python examples/fem_test.py
1/60 0 -1/360 0 -1/90 -1/360
0 4/45 0 2/45 2/45 -1/90
-1/360 0 1/60 -1/90 0 -1/360
0 2/45 -1/90 4/45 2/45 0
-1/90 2/45 0 2/45 4/45 0
-1/360 -1/90 -1/360 0 0 1/60
real 0m2.486s
user 0m2.460s
sys 0m0.024s
And that's it. BTW, we are still reviewing the patch, so for example
it errorneously modifies your example by doing .expand() in there,
which is not necessary, etc.
So still 8x slower than ginac, but I consider it a very successful
result. Maybe it can still be improved, for example we are converting
to a polynomial and back, many times, this slows things down. The
problem is, I don't know a clear cut way to over come that, see:
http://code.google.com/p/sympy/issues/detail?id=477
If you (and anyone) could give us your thoughts on that issue, it'd be cool.
Ondrej
kent-and
> So you want:
>
> In [4]: f = x+y
>
> In [5]: g = x+z
>
> In [6]: f*g
> Out[6]: (x + y)*(x + z)
>
> In [7]: f+g
> Out[7]: y + z + 2*x
>
> So [6] is ok, but you want [7] to be "x+y+x+z" instead? Is it because
> of memory requirements (see that issue for more details)?
>
>
f = pow(x+1, 12) + pow(x-1, 12)
This is a very efficient representation of f in terms of the number of
operations
needed to evaluate f(x). Once expanded it is not easy to go back to
an
efficient representation without knowing what f was.
BTW: I've tested the patch, great speed-up, great work!
Kent
Ondrej Certik
But I think both ginac and SymPy work in this way already:
In [1]: (x+1)**12 + (x-1)**12
Out[1]:
12 12
(1 + x) + (1 - x)
In [2]: f = (x+1)**12 + (x-1)**12
In [3]: f
Out[3]:
12 12
(1 + x) + (1 - x)
In [4]: f.expand()
Out[4]:
12 2 10 4 8 6
2 + 2*x + 132*x + 132*x + 990*x + 990*x + 1848*x
SymPy's automatic evaluation works in a way to only simplify things,
that can be done in a very fast way. Expanding doesn't belong to that
categhory.
> BTW: I've tested the patch, great speed-up, great work!
Thanks! It will be committed and released soon.
Ondrej