[Feature Requests] meshgrid

[kojix2]

Hello.

In my case, meshgrid is the most frequently used feature that Numo:: NArray doesn't have.
I would appreciate it if you could consider the implementation.

Thank you. Have a nice day.

[kojix2]

Honestly, I don't know the right way to create a 3D mesh grid with Numo::NArray.

[himotoyoshi]

This is a sample implementation of meshgrid that works for any dimensions and any data types.

require "numo/narray"

class Numo::NArray

  def self.meshgrid (*axes)
    shape = axes.map(&:size)
    axes.map.with_index do |axis, k|
      extended_shape = Array.new(shape.size) { |i| ( i == k ) ? axis.size : 1 }
      axis.class.new(*shape).store(axis.reshape(*extended_shape))
                            # broadcasting to `shape` occurs here
    end
  end

end

a = Numo::Float32.cast([1,2,3])
b = Numo::Int32.cast([10,20,30])
c = Numo::RObject.cast(["a","b","c"])

Numo::NArray.meshgrid(a,b,c)
  # [Numo::SFloat#shape=[3,3,3]
  # [[[1, 1, 1], 
  #   [1, 1, 1], 
  #   [1, 1, 1]], 
  #  [[2, 2, 2], 
  #   [2, 2, 2], 
  #   [2, 2, 2]], 
  #  [[3, 3, 3], 
  #   [3, 3, 3], 
  #   [3, 3, 3]]],
  #  Numo::Int32#shape=[3,3,3]
  # [[[10, 10, 10], 
  #   [20, 20, 20], 
  #   [30, 30, 30]], 
  #  [[10, 10, 10], 
  #   [20, 20, 20], 
  #   [30, 30, 30]], 
  #  [[10, 10, 10], 
  #   [20, 20, 20], 
  #   [30, 30, 30]]],
  #  Numo::RObject#shape=[3,3,3]
  # [[["a", "b", "c"], 
  #   ["a", "b", "c"], 
  #   ["a", "b", "c"]], 
  #  [["a", "b", "c"], 
  #   ["a", "b", "c"], 
  #   ["a", "b", "c"]], 
  #  [["a", "b", "c"], 
  #   ["a", "b", "c"], 
  #   ["a", "b", "c"]]]]

[kojix2]

Nice. But is this exactly the same as numpy's meshgrid?
(Even though the keyword argument is obviously less than numpy's meshgrid)
I think this should be added to NArray if there are no performance issues. meshgrid is used a lot.

[himotoyoshi]

But is this exactly the same as numpy's meshgrid?

Sorry, it is not.

It is just a sample code for your comment as a Proof of Concept in Ruby level (not in C-extension).

Honestly, I don't know the right way to create a 3D mesh grid with Numo::NArray.

More works are needed to realize the same behavior of numpy's meshgrid and also need the input check.

[kojix2]

Good.

There are a lot of methods in the NArray class that have been implemented in Ruby.
Important methods such as hstack and vstack are also implemented in Ruby.
I don't think you have to necessarily implement meshgrid in C.

https://github.com/ruby-numo/numo-narray/blob/master/lib/numo/narray/extra.rb

[himotoyoshi]

This is a sample implementation of meshgrid in Ruby that can accept full numpy's options.

require "numo/narray"

class Numo::NArray

  # full version
  def self.meshgrid (*axes, indexing: "xy", copy: true, sparse: false)
    if sparse != true and copy != true
      raise NotImplementedError, %{can not create coordinate array as a view into the original axis vector}
    end
    case indexing 
    when "xy"
      # no operation
    when "ij"
      axes = axes.map{|axis| axis.seq }
    else
      raise ArgumentError, %{indexing option should be one of "xy" and "ij"}
    end
    shape = axes.map(&:size)
    if sparse
      axes.map.with_index do |axis, k|
        extended_shape = Array.new(shape.size) { |i| ( i == k ) ? axis.size : 1 }
        if copy
          axis.reshape(*extended_shape)
        else
          axis.reshape!(*extended_shape)
        end
      end
    else
      axes.map.with_index do |axis, k|
        extended_shape = Array.new(shape.size) { |i| ( i == k ) ? axis.size : 1 }
        if copy
          axis.class.new(*shape).store(axis.reshape(*extended_shape))
        else
          # Currently, the behavior for here can not be implemented because of the lack of 
          # the ability to generate a view that represent a repetition of the original vector.
        end
      end
    end
  end

end

In Numo::NArray, the behavior specified by the option copy=false and sparse=false can not be implemented because of the lack of the ability to generate a view that represent a repetition of the original vector. On the other hand, the behavior for copy=false and sparse=true can be simulated using reshape! method.

Currently, I think it's better to leave out the copy option from the method because it's confusing. In that case, if sparse=true is specified, it's better to keep the behavior corresponding to copy=false (the one that uses #reshape!) to save memory. This is because it is considered that few people change the contents of the sparse grid arrays returned from meshgrid.

The latter simple verions is here,

class Numo::NArray

  # simple version witout copy option
  def self.meshgrid (*axes, indexing: "xy", sparse: false)
    case indexing 
    when "xy"
      # no operation
    when "ij"
      axes = axes.map{|axis| axis.seq }
    else
      raise ArgumentError, %{indexing option should be one of "xy" and "ij"}
    end
    shape = axes.map(&:size)
    if sparse
      axes.map.with_index do |axis, k|
        extended_shape = Array.new(shape.size) { |i| ( i == k ) ? axis.size : 1 }
        axis.reshape!(*extended_shape)
      end
    else
      axes.map.with_index do |axis, k|
        extended_shape = Array.new(shape.size) { |i| ( i == k ) ? axis.size : 1 }
        axis.class.new(*shape).store(axis.reshape(*extended_shape))
      end
    end
  end

end

[himotoyoshi]

Sorry. I misunderstood the order of the repetitive dimensions. Here's the corrected version.

class Numo::NArray

  # simple version witout copy option
  def self.meshgrid (*axes, indexing: "xy", sparse: false)
    case indexing 
    when "xy"
      # no operation
    when "ij"
      axes = axes.map{|axis| axis.seq }
    else
      raise ArgumentError, %{indexing option should be one of "xy" and "ij"}
    end
    shape = axes.map(&:size)
    
    if sparse
      axes.map.with_index do |axis, k|
        extended_shape = (shape.size-1).downto(0).map { |i| ( i == k ) ? axis.size : 1 }
        axis.reshape!(*extended_shape)
      end
    else
      axes.map.with_index do |axis, k|
        extended_shape = (shape.size-1).downto(0).map { |i| ( i == k ) ? axis.size : 1 }
        axis.class.new(*shape).store(axis.reshape(*extended_shape))
                              # broadcasting to `shape` occurs here
      end
    end
  end

end

[himotoyoshi]

Sorry, It wasn't fixed yet.

class Numo::NArray

  # simple version witout copy option
  def self.meshgrid (*axes, indexing: "xy", sparse: false)
    case indexing 
    when "xy"
      # no operation
    when "ij"
      axes = axes.map{|axis| axis.seq }
    else
      raise ArgumentError, %{indexing option should be one of "xy" and "ij"}
    end
    shape = axes.map(&:size).reverse
    if sparse
      axes.map.with_index do |axis, k|
        extended_shape = (shape.size-1).downto(0).map { |i| ( i == k ) ? axis.size : 1 }
        axis.reshape!(*extended_shape)
      end
    else
      axes.map.with_index do |axis, k|
        extended_shape = (shape.size-1).downto(0).map { |i| ( i == k ) ? axis.size : 1 }
        axis.class.new(*shape).store(axis.reshape(*extended_shape))
                              # broadcasting to `shape` occurs here
      end
    end
  end

end

For example,

require "numo/narray"

a = Numo::Float32.cast([1,2])
b = Numo::Int32.cast([10,20,30])
c = Numo::RObject.cast(["a","b","c","d"])

p Numo::NArray.meshgrid(a,b,c, indexing: 'xy', sparse: false)

It generate the result,

[Numo::SFloat#shape=[4,3,2]
[[[1, 2], 
  [1, 2], 
  [1, 2]], 
 [[1, 2], 
  [1, 2], 
  [1, 2]], 
 [[1, 2], 
  [1, 2], 
  [1, 2]], 
 [[1, 2], 
  [1, 2], 
  [1, 2]]], Numo::Int32#shape=[4,3,2]
[[[10, 10], 
  [20, 20], 
  [30, 30]], 
 [[10, 10], 
  [20, 20], 
  [30, 30]], 
 [[10, 10], 
  [20, 20], 
  [30, 30]], 
 [[10, 10], 
  [20, 20], 
  [30, 30]]], Numo::RObject#shape=[4,3,2]
[[["a", "a"], 
  ["a", "a"], 
  ["a", "a"]], 
 [["b", "b"], 
  ["b", "b"], 
  ["b", "b"]], 
 [["c", "c"], 
  ["c", "c"], 
  ["c", "c"]], 
 [["d", "d"], 
  ["d", "d"], 
  ["d", "d"]]]]

[kojix2]

Another implementation of the mesh grid by @show-o-atakun

• https://github.com/show-o-atakun/meshgrid_numo_narray

The demand for meshgrid is continuing.
I hope that some implementation of meshgrid will be added to NArray.