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