Re: The nanflag parameter

[Ștefan-Gabriel Mirea]

Hi Rik,

Thank you for your reply.

On Mon, Jan 16, 2017 at 2:54 AM, Rik <address@hidden> wrote:
> So, my fundamental question for Octave-Maintainers is whether we need a
> nanflag in liboctave, or only in the Octave interpreter?

If the performance deficit turns out to be acceptable, is there a way
to also compute the index array (the second output) using high level
functions? Unluckily, the find function cannot operate along a given
dimension. This question also arises when calling with one output if
the argument is a complex array, as NaN complex values can have
regular imaginary parts and then min/max must return the first one
met. For example, in MATLAB:

>> max([1, NaN + 2i, NaN + 3i], [], 'includenan')

ans =

      NaN + 2.0000i

Otherwise, as far as I can see I'd have to manually go through the
array content. For NDArrays, I'd have to use an index tuple (maybe
somehow simulating that the loop along dim is the innermost, in order
to be able to break it at the first NaN found, unless a full but
cache-friendly loop was preferably), while for Sparse(Complex)Matrix
I'd operate on the three underlying vectors. Is this correct?

As for the two-array-input form, there is always only one output, so a
method similar to your pseudocode is enough for non-complex inputs
including SparseMatrix (by using "or" instead of "any", which
fortunately knows to perform broadcasting).

For two complex input arrays, the current default behaviour is
surprisingly equivalent to "includenan" (because of octave::math::max
with complex parameters, which propagates NaN's). The following
pseudocode, although tricky, is the best I could think of for the
"omitnan" case of two complex inputs:

tmp = __max__(x, y);  # Call old max function
if (! includenan)
  x_extended = x + zeros(size(y));
  y_extended = y + zeros(size(x));

  x_extended_isnan = isnan(x_extended); # or isnan(x) | zeros(size(y));
  non_NaN_operands = ifelse(x_extended_isnan, y_extended, x_extended);

  idx = xor(isnan(x), isnan(y));
  tmp(idx) = non_NaN_operands(idx);
endif

Of course, suppose the equivalent C++ code. It also works on
SparseComplexMatrices.

Stefan