Re: Using functions from the math or numpy modules

On May 29, 2008, at 10:22 AM, Etienne Rivard wrote:

> After considering many possible solutions for my problem, I am now  
> back to FiPy.

We win!

> I could find almost every tool I need in FiPy,

Great to know.

> but there is at least one problem remaining.

Quelle dommage!

> The manual of FiPy says that functions such as "exp", etc. are  
> available from "fipy.tools.numerix", I've tried that and as the  
> attached script shows, it works perfectly. However, if "exp" is  
> taken from the "math" or the "numpy" module, FiPy gives an error.

Correct. That's what we mean by "Warning: Generally, things will not  
work as expected if the equivalent function is used from the Numeric  
library." (although we need to update that warning to reflect that we  
use NumPy now).

Clearly more explanation of *why* is in order, though.

The "math" module provides a rudimentary set of functions to operate  
on *single* (floating point) numbers. Everything about FiPy is  
designed around *fields* of numbers, so this module is not useful at  
all and should never be used with FiPy.

As it says in the "numpy" package documentation:

     NumPy
     ==========
     Provides
        1) An array object of arbitrary homogeneous items
        2) Fast mathematical operations over arrays
        3) Linear Algebra, Fourier Transforms, Random Number Generation

NumPy arrays are what FiPy uses internally for the storage of its  
field values, so that's good as far as it goes. However, FiPy depends  
on an additional functionality that NumPy doesn't know anything about.  
FiPy Variable objects are capable of recording functional  
relationships with what's called "lazy evaluation". If you write

     >>> m = 2.
     >>> b = 5.
     >>> x = Variable(3.)
     >>> y = m * x + b

then y contains a functional relationship, as seen by

     >>> y
     ((Variable(value=array(3.0)) * 2.0) + 5.0)

If we actually want its value, we need to do any of `print y`,  
`y.getValue()` or `y()`

     >>> print y
     11.0

The advantage of this is that we can change the value of x, and y will  
update automatically:

     >>> x.setValue(4.)
     >>> print y
     13.0

The "lazy evaluation" part means that nothing will actually be  
calculated until it actually needs to be.

I don't remember the syntax in Matlab, but this is roughly equivalent  
to the distinction between

    y = m * x + b

and

    y := m * x + b

in Mathematica. The first just assigns a value to y; the second  
assigns a dynamic, functional relationship. In FiPy's case, y only  
will see changes in x. By declaring m and b as simple numbers, there  
is no way to change their values and have that change propagate.  
Likewise, if you write

    >>> x = 4.

then you just *replace* the Python variable pointing to the object  
`Variable(3.)` with a new Python variable pointing to the floating  
point number `4.` The expression for y will still work, but you'll no  
longer have an easy way to access and change the value of  
`Variable(3.)`.

In the case of your example, even if you didn't get an error[*] as a  
result of

     def F(y):
         return numpy.exp(5.73*y) - numpy.exp(-11.46*y)

     eq1 = TransientTerm() == ImplicitDiffusionTerm(coeff=0.024) -  
F(u1-u2)
     eq2 = TransientTerm() == ImplicitDiffusionTerm(coeff=0.170) +  
F(u1-u2)

you still would not get the right answer, because eq1 and eq2 would  
have fixed sources with the value of F(u1-u2) established at the  
beginning of the simulation, i.e., F(u1-u2) = F(1. - 0.) for all time,  
which is clearly not what's desired.

[*] I was not expecting to get "TypeError: The coefficient must be a  
rank-0 CellVariable or a scalar value.", but am not altogether  
surprised by it and I don't plan to try to "fix" it. We used to be  
more tolerant of coefficients and sources that had the correct shape,  
but not the correct type, but this caused other problems and so we're  
stricter now. If you really, *really* wanted this fixed source, you  
could obtain it easily enough by writing `CellVariable(mesh=mesh,  
value=F(u1-u2))`.

> I noticed, for example, that the output of "math.exp" is a "float"  
> object and the output of "fipy.tools.numerix.exp" is a  
> "numpy.float64" object.

That's only true in the case of single arguments. If try to pass a  
field (array) of values, math will throw an error and numerix will  
return a "numpy.ndarray", e.g.

     >>> type(math.exp(2.))
     <type 'float'>
     >>> type(numpy.exp(2.))
     <type 'numpy.float64'>
     >>> type(math.exp(numerix.arange(0., 2., 0.1)))
     Traceback (most recent call last):
       File "<stdin>", line 1, in <module>
     TypeError: only length-1 arrays can be converted to Python scalars
     >>> type(numpy.exp(numerix.arange(0., 2., 0.1)))
     <type 'numpy.ndarray'>

> I tried converting so that my function returns this type of object,  
> but I haven't had much success with that method.

You've learned a valuable lesson: type-casting is evil.  8^)

> That's a problem for me because I want to use functions such as  
> "scipy.integrate.quadrature" and "scipy.special.gamma" in my PDEs.

What you need to do is create an "operator Variable" object to perform  
the operation you want. This is not presently documented anywhere, so  
there's no reason you should have known about it unless you were  
*really* familiar with FiPy's internals. We naively never thought that  
FiPy users would need to care about this, but clearly need to  
reevaluate that conclusion.

If we didn't already provide numerix.exp(), then you would need to do:

     >>> v = CellVariable(mesh=mesh, ...)
     >>> exp_v = v._UnaryOperatorVariable(lambda a: numerix.exp(a))

This is precisely what FiPy is doing when you call numerix.exp(v).  
Note that you don't get infinite recursion by defining numerix.exp()  
in terms of numerix.exp() because the arguments are different and  
numerix.exp() is written to do different things with different  
arguments. When you call

     >>> numerix.exp(v).getValue()

you effectively get

     >>> numerix.exp(v.getValue())

Specifically:

     def exp(arr):
         if _isPhysical(arr):
             return arr.exp()
         elif type(arr) is type(array((0))):
             return NUMERIX.exp(arr)
         else:
             return umath.exp(arr)

_isPhysical() is an internal function to identify FiPy's `Variable`  
and `PhysicalField` objects.
NUMERIX is an alias for the underlying array library that FiPy is  
using (always numpy right now, but it used to be Numeric, and it might  
be something else in the future).
umath is like the math package, but better and designed for  
consistency with numpy.

The _UnaryOperatorVariable() method accepts any function that takes a  
single (field) argument and returns a single (field) value. That  
function does not have to be a "lambda" function, but can, if you  
want, be written out as

     >>> def myfunc(a):
     ...     b = a**2
     ...     c = b + 3
     ...     d = sin(c)
     ...     return d

     >>> myfunc_v = v._UnaryOperatorVariable(myfunc)

However, lambda is the clearest and most concise way to do it for  
everything we've ever needed.

So, for gamma, you would do exactly the same thing, e.g,

     >>> gamma_v = v._UnaryOperatorVariable(lambda a:  
scipy.special.gamma(a))

Depending on what you're integrating, calling  
scipy.integrate.quadrature might be more complicated, but should still  
be possible. Tell me a bit more about it and I'll try to explain what  
to do.

In the meantime, we'll think about how best to make all this public  
and properly documented.