Class: CArray::BincountND

Inherits:
Object
  • Object
show all
Defined in:
lib/carray/bincount_nd.rb

Overview

Joint counts of M discrete integer variables — the discrete sibling of Histogram, where a value is its own bin index and there are no edges. Built by CArray#bincount_nd rather than constructed directly.

For a plain 1-D discrete count use CArray#bincount; for continuous data use CArray#histogram.

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(lengths:, fiber_shape: [], weights_dtype: nil) ⇒ BincountND

Returns a new instance of BincountND.

Allocates a new N-D discrete bincount accumulator.

Parameters:

  • lengths (Array<Integer>) —

    per-dimension label ranges; each must be >= 1.

  • fiber_shape (Array<Integer>) (defaults to: []) —

    shape of the leading axes.

  • weights_dtype (Symbol, nil) (defaults to: nil) —

    data_type for weighted accumulators; nil for pure counts (int64).

Raises:

  • (ArgumentError) 


80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
 
# File 'lib/carray/bincount_nd.rb', line 80

def initialize (lengths:, fiber_shape: [], weights_dtype: nil)
  @lengths = lengths.map(&:to_i)
  raise ArgumentError, "lengths must be a non-empty list" if @lengths.empty?
  @lengths.each_with_index do |l, k|
    raise ArgumentError, "lengths[#{k}] must be >= 1" if l < 1
  end
  @m = @lengths.size
  @fiber_shape = fiber_shape.map(&:to_i).freeze
  @weighted = !weights_dtype.nil?
  @counts_dtype = @weighted ? weights_dtype : :int64
  ext_dims = @lengths.map { |l| l + 1 }       # +1: upper overflow cell
  ext_shape = @fiber_shape + ext_dims
  @full_counts = CArray.public_send(@counts_dtype, *ext_shape).fill(0)
  @sample_axis  = nil
  @channel_axis = nil
end

Instance Attribute Details

#fiber_shape ⇒ Object (readonly)

Returns the value of attribute fiber_shape.



98
99
100
 
# File 'lib/carray/bincount_nd.rb', line 98

def fiber_shape
  @fiber_shape
end

#full_counts ⇒ Object (readonly)

Returns the value of attribute full_counts.



98
99
100
 
# File 'lib/carray/bincount_nd.rb', line 98

def full_counts
  @full_counts
end

#lengths ⇒ Object (readonly)

Returns the value of attribute lengths.



98
99
100
 
# File 'lib/carray/bincount_nd.rb', line 98

def lengths
  @lengths
end

#m ⇒ Object (readonly)

Returns the value of attribute m.



98
99
100
 
# File 'lib/carray/bincount_nd.rb', line 98

def m
  @m
end

Instance Method Details

#+(other) ⇒ BincountND

Returns a new BincountND whose counts are the element-wise sum of self and other. Both operands must share lengths, fiber_shape, and weighted state.

Parameters:

Returns:

Raises:

  • (ArgumentError) —

    when structure does not match.

Raises:

  • (ArgumentError) 


290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
 
# File 'lib/carray/bincount_nd.rb', line 290

def + (other)
  raise ArgumentError, "type mismatch" unless other.is_a?(BincountND)
  raise ArgumentError, "M mismatch" unless @m == other.m
  raise ArgumentError, "lengths mismatch" unless @lengths == other.lengths
  raise ArgumentError, "fiber_shape mismatch" unless @fiber_shape == other.fiber_shape
  raise ArgumentError, "weighted/unweighted mismatch" unless @weighted == other.weighted?

  result = self.class.send(:new,
                           lengths: @lengths,
                           fiber_shape: @fiber_shape,
                           weights_dtype: @weighted ? @counts_dtype : nil)
  rf = result.instance_variable_get(:@full_counts)
  rf[] = @full_counts + other.full_counts
  result.instance_variable_set(:@sample_axis, @sample_axis)
  result.instance_variable_set(:@channel_axis, @channel_axis)
  result
end

#add(chunk, axis: nil, weights: nil) ⇒ self

Accumulates chunk (per-sample discrete labels) into self. Locks the sample/channel axes on the first call. Labels must be non-negative; labels >= lengths[k] fold into the upper overflow cell of dim k.

Parameters:

  • chunk (CArray) —

    integer labels with shape fiber_shape + (A, M).

  • axis (Array(Integer, Integer), Integer, nil) (defaults to: nil) —

    [sample, channel] axis pair.

  • weights (CArray, nil) (defaults to: nil) —

    per-sample weights (required iff weighted accumulator).

Returns:

  • (self)

Raises:

  • (ArgumentError) —

    on shape / axis / label / weighted mismatch.

Raises:

  • (ArgumentError) 


161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
 
# File 'lib/carray/bincount_nd.rb', line 161

def add (chunk, axis: nil, weights: nil)
  # Keep the labels in their native integer type (no int64 coercion): an
  # int32 label array stays int32 through the ravel, and `bincount` picks
  # a uint32 output when the table fits.  Forcing int64 would materialise
  # a cast of the whole chunk.
  chunk = CArray.wrap_readonly(chunk)

  # M=1 convenience: accept chunks without the trailing channel axis.
  if @m == 1 && chunk.ndim == @fiber_shape.size + 1
    chunk = chunk.reshape(*(chunk.shape + [1]))
    if axis.is_a?(Integer)
      ax = CArray.normalize_axis(axis, chunk.ndim - 1, "add axis")
      axis = [ax, chunk.ndim - 1]
    end
  end

  ax = axis || [-2, -1]
  ax = [ax] if ax.is_a?(Integer)
  raise ArgumentError, "axis must be [sample, channel]" unless ax.is_a?(Array) && ax.size == 2
  sample_ax  = CArray.normalize_axis(ax[0], chunk.ndim, "sample axis")
  channel_ax = CArray.normalize_axis(ax[1], chunk.ndim, "channel axis")
  raise ArgumentError, "same axis used twice" if sample_ax == channel_ax

  if @sample_axis.nil?
    @sample_axis = sample_ax
    @channel_axis = channel_ax
  elsif @sample_axis != sample_ax || @channel_axis != channel_ax
    raise ArgumentError,
          "axis mismatch (locked at [#{@sample_axis}, #{@channel_axis}], got [#{sample_ax}, #{channel_ax}])"
  end

  expected_ndim = @fiber_shape.size + 2
  unless chunk.ndim == expected_ndim
    raise ArgumentError,
          "chunk.ndim=#{chunk.ndim} expected #{expected_ndim} " \
          "(fiber #{@fiber_shape.inspect} + sample + channel)"
  end
  unless chunk.shape[channel_ax] == @m
    raise ArgumentError, "channel axis length #{chunk.shape[channel_ax]} != M=#{@m}"
  end
  chunk_fiber = chunk.shape.dup
  [sample_ax, channel_ax].sort.reverse.each { |p| chunk_fiber.delete_at(p) }
  unless chunk_fiber == @fiber_shape
    raise ArgumentError,
          "fiber shape mismatch: chunk yields #{chunk_fiber.inspect}, expected #{@fiber_shape.inspect}"
  end

  return self if chunk.shape[sample_ax] == 0

  if weights
    raise ArgumentError, "weights given but accumulator is unweighted" unless @weighted
    weights = CArray.wrap_readonly(weights, @counts_dtype)
    expected_w_shape = chunk.shape.dup
    expected_w_shape.delete_at(channel_ax)
    unless weights.shape == expected_w_shape
      raise ArgumentError,
            "weights shape #{weights.shape.inspect} expected #{expected_w_shape.inspect}"
    end
  elsif @weighted
    raise ArgumentError, "weights required (accumulator is weighted)"
  end

  # --- ravel + bincount --------------------------------------------
  # Each label is its own bin: clamp to the upper overflow cell and ravel
  # the M channels into one flat index, then let the dedicated bincount
  # kernel scatter.  Discrete "binning" is a cheap, vectorisable ravel, so
  # this beats a hand-fused scalar kernel (bench: a fused C kernel was
  # ~4.4 vs ~1.6 ns/sample).  With fibers we loop one small ravel+bincount
  # per fiber so each fiber's counts slice stays L1-resident, rather than
  # one giant bincount over the whole F*total_ext table (which is cache-
  # cold and ~1.7x slower).  See devel/bench_bincount_nd_gate.rb.
  ext_sizes   = @lengths.map { |l| l + 1 }      # +1: upper overflow cell
  strides_ext = ext_sizes.each_with_index.map { |_, k| ext_sizes[(k + 1)..].inject(1, :*) }
  total_ext   = ext_sizes.inject(:*)
  widen       = total_ext > 0x7fffffff          # int64 flat for big joint tables

  # canonical [fiber..., sample, channel] view (channel last); weights to
  # [fiber..., sample].  Skip the transpose when the layout is already
  # canonical (the usual case) — a transpose view would force `reshape`
  # below to materialise a full copy.
  fiber_axes = (0...chunk.ndim).to_a - [sample_ax, channel_ax]
  perm   = fiber_axes + [sample_ax, channel_ax]
  tchunk = perm == (0...chunk.ndim).to_a ? chunk : chunk.transpose(*perm)
  tweights = nil
  if weights
    shift  = ->(p) { p < channel_ax ? p : p - 1 }
    w_perm = fiber_axes.map(&shift) + [shift.call(sample_ax)]
    tweights = w_perm == (0...weights.ndim).to_a ? weights : weights.transpose(*w_perm)
  end

  # One pass for the negative-label check (masked-aware).  `min` returns
  # UNDEF when every sample is masked: that is a well-defined no-op (all
  # samples dropped -> counts unchanged), so bail before the label-range
  # checks below (`chunk.min` / `b.max` would otherwise hit UNDEF and the
  # FLAT path arithmetic would raise on it).
  mn = chunk.min
  return self if mn == UNDEF
  raise ArgumentError, "bincount_nd: negative label" if mn < 0

  if @fiber_shape.empty?
    # Flat: the ravel is cheap + vectorisable, so the separate
    # vectorised ravel + tuned `bincount` beats any fused kernel.
    # Clamp a channel only when it actually overflows (decided once).
    ravel = nil
    (0...@m).each do |k|
      b = tchunk[nil, k]
      b = b.clip(0, @lengths[k]) if b.max > @lengths[k] - 1
      b = b.int64 if widen
      term = strides_ext[k] == 1 ? b : b * strides_ext[k]
      ravel = ravel.nil? ? term : ravel + term
    end
    chunk_counts = ravel.bincount(weights: tweights, length: total_ext)
    @full_counts[] = @full_counts + chunk_counts.reshape(*@full_counts.shape)
  else
    # Fiber: a dedicated C kernel counts each fiber into its own
    # L1-resident counts slice in one pass (no per-fiber Ruby loop, no
    # giant cache-cold bincount, no int coercion).  Clamp is inline in C.
    tchunk.send(:bincount_nd_count_ki, @full_counts, tweights)
  end
  self
end

#counts ⇒ CArray

Returns the in-range counts view with shape fiber_shape + (L_0, ..., L_{M-1}), excluding the upper overflow cell.

Returns:



105
106
107
108
 
# File 'lib/carray/bincount_nd.rb', line 105

def counts
  idx = [nil] * @fiber_shape.size + @lengths.map { |l| 0...l }
  @full_counts[*idx]
end

#overflow(axis: nil) ⇒ CArray

Upper-overflow marginal on dim axis (= samples whose dim-axis label was >= length[axis]); other dims marginalised. shape = fiber_shape. For M=1, axis: may be omitted.

Returns the upper-overflow marginal on dimension axis (samples whose dim-axis label was >= lengths[axis]); other dimensions are marginalised. For 1-D accumulators axis may be omitted.

Parameters:

  • axis (Integer, nil) (defaults to: nil) —

    dimension to marginalise.

Returns:

Raises:

  • (ArgumentError) —

    when axis is required but omitted.

Raises:

  • (ArgumentError) 


121
122
123
124
125
126
127
128
129
 
# File 'lib/carray/bincount_nd.rb', line 121

def overflow (axis: nil)
  raise ArgumentError, "axis: keyword required (M=#{@m})" if axis.nil? && @m > 1
  ax = axis.nil? ? 0 : CArray.normalize_axis(axis, @m, "overflow")
  base = [nil] * @fiber_shape.size
  bin_idx = (0...@m).map { |k| k == ax ? @lengths[k] : nil }   # overflow cell on ax
  slice = @full_counts[*(base + bin_idx)]
  (@m - 1).times { slice = slice.accumulate(axis: slice.ndim - 1) }
  slice
end

#overflow_total ⇒ CArray

Returns the per-fiber count of samples that overflowed on any dimension.

Returns:



143
144
145
 
# File 'lib/carray/bincount_nd.rb', line 143

def overflow_total
  sum_along_bin_axes(@full_counts) - sum_along_bin_axes(counts)
end

#total ⇒ CArray

Returns the per-fiber sample total (in-range plus overflow) with shape fiber_shape.

Returns:



135
136
137
 
# File 'lib/carray/bincount_nd.rb', line 135

def total
  sum_along_bin_axes(@full_counts)
end