how to handle a parameter with a possible value of infinity?

16 days ago by
What is best practice for handling a parameter which is usually finite but may take the value of infinity?

Think of a form containing a term like "f/C", where is some function, and C is a parameter.
If C has a value of, say, numpy.inf, then ideally the code still works, and the term simply gets cancelled out without causing any other trouble.  This might correspond to an asymptotic solution, for instance.

At the level of numpy that works: 1/numpy.inf is 0.  But applied to FEniCS functions, ffc complains when the script is processed: e.g.
error: ‘inf’ was not declared in this scope
     sp[2] = sp[1] + 0.001 * sp[0] / inf;

It looks like numpy.inf is converted in C++ to "inf" and treated as a variable that is otherwise undefined.

Here is a FEniCS python minimal example:
from dolfin import *
import numpy

mesh = IntervalMesh(100,0,10)
V = FunctionSpace(mesh,"P",2)

fexp = Expression('exp(-x[0])',degree=2)
fcos = Expression('cos(x[0])',degree=2)

#C=0.1  # displays a plot of a perturbed exponential curve
C=numpy.inf   # FFC (JIT) -> error: ‘inf’ was not declared in this scope

f = fexp*(1 + fcos/C)


import matplotlib.pyplot as plt

The general idea of handling inf as a value of a variable still works in principle in C++. i.e. if appropriately defined, e.g. inf=std::numeric_limits<double>::infinity(), then 1/inf is 0 in C++, c.f.
#include <iostream>
#include <limits>

int main()
  double a;

  a = 2;
  std::cout << "finite a=" << a << " inverts to " << 1/a << std::endl;

  a = std::numeric_limits<double>::infinity();
  std::cout << "infinite a=" << a << " inverts to " << 1/a << std::endl;
  return 0;
which outputs
finite a=2 inverts to 0.5
infinite a=inf inverts to 0​

But at the moment FFC appears to simply be converting numpy.inf to a "variable", inf.   Should this be considered a bug in FFC?   Or is there an entirely better way of handling the situation?
Community: FEniCS Project

1 Answer

16 days ago by
Wrapping it as a Constant seems to work.  All of the following possibilities work in the MWE:
C = Constant('inf')​
C = Constant(numpy.inf)
C = Constant(float('inf'))
Good thinking David, that does work. 

I wonder if there is a clear winner in terms of performance.  The alternatives are
1) litter the script with isfinite tests to only include the term if the value of the parameter is finite
2) wrap the parameter in Constant() (your solution)
3) patch FFC to treat inf as C++ inf, not as a "variable" inf

Option 1) makes the code less tidy. might make runtime setup a touch slower but perhaps runtime overall is faster if the term is not included when solving.  One the other hand, tidier code is easier to maintain.
written 16 days ago by Drew Parsons  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »