ModelicaStandardLibrary icon indicating copy to clipboard operation
ModelicaStandardLibrary copied to clipboard

Invalid function Modelica.Math.Matrices.nullSpace calculating max({})

Open rfranke opened this issue 4 years ago • 11 comments

The function ModelicaTest.Math.Matrices3 contains

  import Modelica.Math.Matrices;
  Real N[:,:]= fill(0,0,0);
  Real Xn[0,0];
  Integer n;
algorithm
  (Xn, n) := Matrices.nullSpace(N);

Modelica.Math.Matrices.nullSpace calculates:

  Real sigma[min(size(N, 1), size(N, 2))] "Singular values";
algorithm
  // rank computation
  eps := max(size(N, 1), size(N, 2)) * max(sigma) * Modelica.Constants.eps;

i.e. max(0, 0) * max({}) * Modelica.Constants.eps.

Can this pass any simulation with success?

rfranke avatar Jul 08 '21 18:07 rfranke

The function ModelicaTest.Math.Matrices3 contains

  import Modelica.Math.Matrices;
  Real N[:,:]= fill(0,0,0);
  Real Xn[0,0];
  Integer n;
algorithm
  (Xn, n) := Matrices.nullSpace(N);

Modelica.Math.Matrices.nullSpace calculates:

  Real sigma[min(size(N, 1), size(N, 2))] "Singular values";
algorithm
  // rank computation
  eps := max(size(N, 1), size(N, 2)) * max(sigma) * Modelica.Constants.eps;

i.e. max(0, 0) * max({}) * Modelica.Constants.eps.

Can this pass any simulation with success?

Sort of, max for an empty range is defined as the least value of the type (in this case Real so -Modelica.Constants.inf). https://specification.modelica.org/master/arrays.html#reduction-expressions Multiplying that by 0 and eps should give 0 as result.

However, that is only for a reduction expression and this is max for an array. We should clarify that it handles the empty case in the same manner.

HansOlsson avatar Jul 09 '21 10:07 HansOlsson

The multiplication only gives 0 if a Modelica tool treats Modelica.Constants.inf as large number, such as 1e60. IEEE 754 defines 0 * Infinity = NaN Modelica should be developed towards IEEE 754, i.e. allow ModelicaServices.Machine.inf defined as IEEE 754 Infinity.

rfranke avatar Jul 09 '21 11:07 rfranke

I agree that we should move towards IEEE 754 definitions.

However, I would propose that we view Modelica.Constants.inf as DBL_MAX not as true "infinite" see #2056 for a number of reasons:

  • Code relying on computations involving inf (such as this one) is more likely to continue to work. Not all - see #3390 (and for eps #3357 ).
  • It sort of corresponds to the description (but I'm not consistent - I view the descriptions for eps and small as wrong).

HansOlsson avatar Jul 09 '21 11:07 HansOlsson

It sounds like a mistake to couple the definitions of max and Modelica.Constants.inf. To get sensible properties of max and min, they really need to return the least value, so if there is a IEEE 754 Infinity, one needs to use it. Hence, I think the current formulation of max is good, but nullSpace is a good example of what might break in case Real is extended (?!) to also allow ­+/-Infinity.

henrikt-ma avatar Jul 12 '21 08:07 henrikt-ma

It sounds like a mistake to couple the definitions of max and Modelica.Constants.inf.

That is the current case.

To get sensible properties of max and min, they really need to return the least value, so if there is a IEEE 754 Infinity, one needs to use it.

The second part does not follow, and that's why I argue for using the maximum representable number as in C (DBL_MAX) as it has several advantages:

  • we have existing code, such as this one, that may fail if we use a infinite number.
  • I don't think we have any code in MSL or other standard libraries that relies on it being truly infinite.
  • it is easier to generate as portable code.

Obviously it would be clearer if this and other constants were given clear descriptions and names matching that.

HansOlsson avatar Jul 12 '21 09:07 HansOlsson

It sounds like a mistake to couple the definitions of max and Modelica.Constants.inf.

That is the current case.

They way it is said looks more like a non-normative observation of what Modelica.Constants.inf is than a definition of max in terms of Modelica.Constants.inf. It would be much better it the specification (regardless of the MSL) stated whether Real includes the infinities or not. This would make writing of portable code with max({}) much easier.

To get sensible properties of max and min, they really need to return the least value, so if there is a IEEE 754 Infinity, one needs to use it.

The second part does not follow

Allowing Real to include infinities, but not using it for max({}) would be catastrophic for reasoning and symbolic manipulation. I think it's safe to call these properties essential:

  • max(x, y)x
  • max(x, y)y
  • max({x}) = x
  • max({x, y1, …, yn}) = max(x, max({y1, …, yn}))

With n = 0 and x = -Infinity this means:

  • max({})max(-Infinity, max({})) = max({-Infinity}) = -Infinity

and that's why I argue for using the maximum representable number as in C (DBL_MAX)

My take on this is that it means we need to limit Real to finite values.

henrikt-ma avatar Jul 12 '21 10:07 henrikt-ma

It is hard (appears artificial) to restrict Modelica Real to finite values, because almost all models are executed on IEEE 754 machines.

Here is how Python treats this since version 3.4:

>>> max([], default = 1)
1

"The default argument specifies an object to return if the provided iterable is empty. If the iterable is empty and default is not provided, a ValueError is raised."

Modelica max could be extended accordingly, in order to avoid implicit definitions and problems with IEEE 754. The function nullSpace could write:

max(sigma, default = 1)

or without Modelica language extension:

(if size(sigma, 1) > 0 then max(sigma) else 1)

rfranke avatar Jul 12 '21 11:07 rfranke

I find the Python way extremely ugly. Instead of making max({}) illegal, I'd much rather break backward compatibility by allowing the infinities in Real (which according to my reading of the specification would mean that they should be used by max and min, and that the comment about the relation to Modelica.Constants.inf would become outdated).

henrikt-ma avatar Jul 12 '21 11:07 henrikt-ma

To cut a long discussion short: #3850 would fix the issue.

rfranke avatar Jul 12 '21 17:07 rfranke

It is hard (appears artificial) to restrict Modelica Real to finite values, because almost all models are executed on IEEE 754 machines.

There are two parts of this:

  • One is how values are handled during the execution of algorithms and I understand your statement in that case.
  • Another is the symbolic manipulation, and that's the tricky part.

If Modelica Real variables aren't restricted to finite values then we can have:

  Real a;
  Real x(min=0.1);
equation
  a*x=a; // Solution x=1 unless a=+/-infinite when any finite non-zero x works

  Real a;
  Real x(min=0, max=10);
equation
  a+x=a; // Solution x=0 or a=+/infinite

That doesn't make sense and would likely limit the optimization possibilities in actual models. Saying that non-finite values are allowed sometimes and sometimes not seem quite complicated.

HansOlsson avatar Jul 14 '21 07:07 HansOlsson

That doesn't make sense and would likely limit the optimization possibilities in actual models. Saying that non-finite values are allowed sometimes and sometimes not seem quite complicated.

This is a good argument, and until we make a deliberate choice to allow the infinities, I think it needs to be more clearly specified that they are not allowed. Some sort implicit ban via the definition of max and a reference to Modelica.Constants.inf doesn't count as being clearly specified; the specification should be self-contained on a matter like this, and we need a statement similar to the recommended minimum range of IntegerType (and similar to IntegerType the statement about recommended range should be moved from Lexical Structure to the definition of RealType).

henrikt-ma avatar Jul 15 '21 08:07 henrikt-ma