Elementwise operation improvements: return types

[translunar]

Depends upon: Issue #3: Improve SyntaxTree / C / math expression parser

The problem is described in this email thread:

Let's say you have the following matrix operation:
a = (x == y)
where x is a float64 NMatrix, y is an int8 NMatrix. This will currently cause y to be copy-casted to float64, and after the operation a will be a float64 NMatrix as well.

(Note that this is distinct from x.eql?(y) which will return simply true or false.)

The alternative -- which is not implemented but could be done somewhat easily -- is to have a be an NMatrix of chars (dtype = :byte), which would then be 0 or 1.

The reason that I have it this way currently is that I figure if you're doing operations at some precision (e.g., float64), you want as few copy-casts as possible. If you're using something like a = x == y to make a an indicator variable, you don't want to have to copy-cast it to a float64 to use it later in your float64 operations.

But it's also a major space suck. Eight bytes versus one!

The problem doesn't only concern booleans. Let's say you want to take the elementwise square root of an int64 matrix. You should be able to choose the return type: who knows, maybe every value in your matrix is already a square number. It makes sense that you should be able to optionally provide a dtype.

For example, == could be shorthand for the following operation:

a.elementwise(:'==', b, :byte)

which would upcast a or b to the correct type, as always, but then return a :byte NMatrix instead of a :float64. But then the coder could still choose to get back a float64:

a.elementwise(:'==', b, :float64)

and even redefine == to always return the upcasted type.

Why does this depend on issue #3? Well, currently, the dense elementwise template looks kind of like this:

int nm_d_%%TYPE_ABBREV%%_elementwise(const %%TYPE%%* A, const %%TYPE%%* B, %%TYPE%%* C, size_t n, enum NMatrix_Ops op) {
  size_t i;

  switch(op) {
    // ...
  case NM_OP_EQEQ: // ==
    for (i = 0; i < n; ++i) {
      %%TYPE C[i] = A[i] == B[i]%%
    }
    break;
  // ...
  }
}

But the parser can only handle operations where A[i], B[i], and C[i] are all the same type. Go look at blas.c to see the output for Ruby objects, for example:

C[i] = rb_funcall(A[i], rb_intern("=="), 1, B[i]);

Well, actually, that's partially broken, I just realized. The above is what you should get, but instead the output looks like what you'd actually want if you were returning a byte matrix:

C[i] = rb_funcall(A[i], rb_intern("=="), 1, B[i]) == Qtrue;

In conclusion, SyntaxTree needs to be aware of the template types in the function header.

It may be that RDoc includes a C parser up to the job. That's one place to look. If not, we may need to write one from scratch. It's also possible that SyntaxTree could be improved.

But ultimately, we need to give more control to the user.