Just to revisit this one last time, the main reason this was slow is because you were accessing x_DG.dat.data_ro[j,0] and x_DG.dat.data_ro[j,1] inside the innermost, i, loop (twice each!). This came with *a lot* of Python overhead. Since these quantities are independent of i, you could have defined temporary variables outside the i loop: x0 = x_DG.dat.data_ro[j,0] x1 = x_DG.dat.data_ro[j,1] and used these variables within the loop. Just doing this makes the code ~4x faster (75% saving) than what was in your email. Getting rid of the innermost i loop and using array-sized operations, as I did before, then gives the rest of the benefit (an extra factor of 20x or more...) On 12 October 2016 at 21:18, Justin Chang <jychang48@gmail.com> wrote:

Yes that works brilliantly. Thanks Andrew!

On Wed, Oct 12, 2016 at 2:37 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:

...

FWIW, this already gives a *massive* speedup:

for j in range(usize): num = np.zeros(600) den = np.zeros(600)

num[:] = datas[:]/(4*np.power(dataxs[:] - x_DG.dat.data_ro[j, 0], 2) + np.power(datays[:] - x_DG.dat.data_ro[j, 1], 2))

den[:] = 1/(4*np.power(dataxs[:] - x_DG.dat.data_ro[j, 0], 2) + np.power(datays[:] - x_DG.dat.data_ro[j, 1], 2))

u.dat.data[j] = 1e-12*np.exp(sum(num)/sum(den))

There might still be a smarter way to do it.

On 12 October 2016 at 20:23, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:

...

Justin,

Your loops are in standard Python. This is likely to be slow.

To give a simple case, if you want to calculate a = b + c, where these are 1000x1000 NumPy arrays, your current code looks like

for j in range(1000): for i in range(1000): a[i, j] = b[i, j] + c[i, j]

David is suggesting that you replace this with the "NumPy operation"

a[:, :] = b[:, :] + c[:, :]

which will be much faster. In this simple case you could just write a = b + c, but in your real case you'll probably need to use slice notation.

On 12 October 2016 at 19:44, Justin Chang <jychang48@gmail.com> wrote:

...

David,

I tried to follow your steps, but I ended up with something that's even worse performance wise. Clearly I am not doing something right.

Attached is my code. Am I not utilizing numpy operations correctly?

Thanks, Justin

On Wed, Oct 12, 2016 at 5:13 AM, David Ham <David.Ham@imperial.ac.uk> wrote:

...

Hi Justin,

Yes.

1. Make a DG0 "coordinate" field by projecting the coordinates into VectorFunctionSpace(mesh, DG, "0") call it, say X_DG 2. Use numpy operations to collectively interpolate your input data to the points given by X_DG.dat.data_ro 3. Set u.dat.data[:] to those interpolated values.

This will result in going through the python stack O(1) times instead of O(meshsize) which should be waaay faster.

On Wed, 12 Oct 2016 at 11:00 Justin Chang <jychang48@gmail.com> wrote:

...

Miklos I am not sure I follow what you're saying.

Also, my code is pretty slow lol, is there a "more efficient" way of doing this or is this something I will just have to deal with?

On Wed, Oct 12, 2016 at 4:46 AM, Homolya, Miklós < m.homolya14@imperial.ac.uk> wrote:

Maybe you need to call the __init__ is the superclass? ------------------------------ *From:* firedrake-bounces@imperial.ac.uk < firedrake-bounces@imperial.ac.uk> on behalf of Justin Chang < jychang48@gmail.com> *Sent:* 12 October 2016 10:45:11 *To:* firedrake *Subject:* Re: [firedrake] Python expression for loading files

Nevermind, I fixed it by adding this line under __init__

self._user_args = []

Code works now. Thanks again for the help

Justin

On Wed, Oct 12, 2016 at 4:38 AM, Justin Chang <jychang48@gmail.com> wrote:

I tried that, but I now get this error:

Traceback (most recent call last): File "readData.py", line 25, in <module> u.interpolate(createData()) File "/home/justin/Software/firedrake/src/firedrake/firedrake/function.py", line 328, in interpolate return interpolation.interpolate(expression, self, subset=subset) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 36, in interpolate return Interpolator(expr, V, subset=subset).interpolate() File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 56, in __init__ self.callable = make_interpolator(expr, V, subset) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 116, in make_interpolator loops.append(_interpolator(V, f.dat, expr, subset)) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 166, in _interpolator kernel, oriented, coefficients = compile_python_kernel(expr, to_pts, to_element, V, coords) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 310, in compile_python_kernel for _, arg in expression._user_args: AttributeError: 'createData' object has no attribute '_user_args'

Here's my code and text files

Thanks, Justin

On Wed, Oct 12, 2016 at 4:36 AM, Homolya, Miklós < m.homolya14@imperial.ac.uk> wrote:

def value_shape(self): return (1,)

Change (1,) to ()

On Wed, Oct 12, 2016 at 10:35 AM +0100, "Justin Chang" < jychang48@gmail.com> wrote:

How do I make the expression return a scalar?

On Wed, Oct 12, 2016 at 4:28 AM, Homolya, Miklós < m.homolya14@imperial.ac.uk> wrote:

Your function space Q is scalar (), while your expression is a 1D vector (1,). Change one of them to the other.

On Wed, Oct 12, 2016 at 10:05 AM +0100, "Justin Chang" < jychang48@gmail.com> wrote:

Dear all,

I am attempting to create an Expression which uses data from some text files. This is my code:

from firedrake import * import numpy as np

mesh = UnitSquareMesh(100,100) Q = FunctionSpace(mesh,'DG',0)

file_prefix = 'perm600' class createData(Expression): def __init__(self): self.s = np.loadtxt(file_prefix+'_s') self.xs = np.loadtxt(file_prefix+'_xs') self.ys = np.loadtxt(file_prefix+'_ys') def eval(self, val, x): numerator = 0 denominator = 0 for i in range(0,600): numerator += self.s[i]/((self.xs[i]-x[0])** 2+(self.ys[i]-x[1])**2) denominator += 1/((self.xs[i]-x[0])**2+(self.ys[i]-x[1])**2) ssum = numerator/denominator val[0] = 1e-12*exp(ssum) def value_shape(self): return (1,)

u = Function(Q) u.interpolate(createData())

However, I get this error:

Traceback (most recent call last): File "readData.py", line 25, in <module> u.interpolate(createData()) File "/home/justin/Software/firedrake/src/firedrake/firedrake/function.py", line 328, in interpolate return interpolation.interpolate(expression, self, subset=subset) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 36, in interpolate return Interpolator(expr, V, subset=subset).interpolate() File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 56, in __init__ self.callable = make_interpolator(expr, V, subset) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 116, in make_interpolator loops.append(_interpolator(V, f.dat, expr, subset)) File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", line 152, in _interpolator % (len(expr.ufl_shape), len(V.ufl_element().value_shape()))) RuntimeError: Rank mismatch: Expression rank 1, FunctionSpace rank 0

What's going on here? I know you guys generally recommend using UFL for interpolating expressions, but I am not sure how I would implement the above in a UFL way since I have to read data from a textfile.

Any help appreciated :)

Justin

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake