runmat_runtime/builtins/array/indexing/
ind2sub.rs

1//! MATLAB-compatible `ind2sub` builtin with GPU-aware semantics for RunMat.
2
3use runmat_accelerate_api::HostTensorView;
4use runmat_builtins::{Tensor, Value};
5use runmat_macros::runtime_builtin;
6
7use super::common::{build_strides, materialize_value, parse_dims, total_elements};
8use crate::builtins::common::spec::{
9    BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
10    ProviderHook, ReductionNaN, ResidencyPolicy, ScalarType, ShapeRequirements,
11};
12use crate::builtins::common::tensor;
13use crate::{make_cell, register_builtin_fusion_spec, register_builtin_gpu_spec};
14
15#[cfg(feature = "doc_export")]
16use crate::register_builtin_doc_text;
17
18#[cfg(feature = "doc_export")]
19pub const DOC_MD: &str = r#"---
20title: "ind2sub"
21category: "array/indexing"
22keywords: ["ind2sub", "linear index", "subscripts", "column major", "gpu", "nd indexing"]
23summary: "Convert MATLAB column-major linear indices into per-dimension subscript arrays."
24references: []
25gpu_support:
26  elementwise: false
27  reduction: false
28  precisions: ["f32", "f64"]
29  broadcasting: "matlab"
30  notes: "When the active provider implements the `ind2sub` hook (WGPU today), conversions run entirely on the GPU. Other providers fall back to the host and re-upload the results."
31fusion:
32  elementwise: false
33  reduction: false
34  max_inputs: 0
35  constants: "inline"
36requires_feature: null
37tested:
38  unit: "builtins::array::indexing::ind2sub::tests"
39  integration: "builtins::array::indexing::ind2sub::tests::ind2sub_gpu_roundtrip"
40---
41
42# What does the `ind2sub` function do in MATLAB / RunMat?
43`ind2sub(siz, idx)` converts MATLAB's 1-based column-major linear indices back into the individual subscripts for each dimension specified in `siz`.
44
45## How does the `ind2sub` function behave in MATLAB / RunMat?
46- The size vector `siz` supplies one extent per dimension. Each entry must be a positive integer.
47- `idx` can be a scalar or an array of any shape; every element must be a positive integer within `prod(siz)`.
48- The result is always a 1×N cell array (N is `numel(siz)`). Each cell contains a double array that matches the shape of `idx`.
49- Non-integer, complex, NaN, Inf, or out-of-range indices raise MATLAB-compatible errors.
50- Empty inputs produce empty outputs with matching shapes.
51- When called with multiple outputs (`[i,j,k] = ind2sub(...)`) RunMat unpacks the cell array automatically, mirroring MATLAB semantics.
52
53## `ind2sub` Function GPU Execution Behaviour
54When a WGPU-backed provider is active, `ind2sub` executes entirely on the GPU. The shader mirrors MATLAB's validation rules, rejecting non-integer, non-positive, or out-of-range indices with the same diagnostic messages as the CPU path. Providers that do not yet implement the hook fall back to the host implementation; after computing the subscripts RunMat uploads the results back to the active provider so downstream fused kernels continue operating on device-resident data.
55
56## Examples of using the `ind2sub` function in MATLAB / RunMat
57
58### Recovering row and column subscripts from a matrix index
59
60```matlab
61[row, col] = ind2sub([3 4], 8);
62```
63
64Expected output:
65
66```matlab
67row = 2;
68col = 3;
69```
70
71### Extracting multiple matrix indices at once
72
73```matlab
74idx = [7 8 9];
75[rows, cols] = ind2sub([3 5], idx);
76```
77
78Expected output:
79
80```matlab
81rows = [1 2 3];
82cols = [3 3 3];
83```
84
85### Converting indices for a 3-D volume
86
87```matlab
88idx = [3 11];
89[r, c, p] = ind2sub([2 3 4], idx);
90```
91
92Expected output:
93```matlab
94r = [1 1];
95c = [2 3];
96p = [1 2];
97```
98
99### Keeping indices on the GPU
100
101```matlab
102rows = gpuArray([1; 2; 3]);
103cols = gpuArray([4; 4; 4]);
104lin = sub2ind([3 4], rows, cols);
105subs = ind2sub([3 4], lin);  % subs{1} and subs{2} remain gpuArray values
106class(subs{1})
107class(subs{2})
108```
109
110Expected output:
111
112```matlab
113ans =
114    'gpuArray'
115
116ans =
117    'gpuArray'
118```
119
120### Using a single output cell for flexible unpacking
121
122```matlab
123subs = ind2sub(size(magic(4)), 6:9);
124% Access with subs{1}, subs{2}, ...
125```
126
127Expected output:
128
129```matlab
130subs =
131  1×2 cell array
132    {1×4 double}    {1×4 double}
133```
134
135### Validating index ranges before reshaping
136
137```matlab
138idx = 1:numel(A);
139[i, j] = ind2sub(size(A), idx);
140B = accumarray([i(:) j(:)], A(:));
141```
142
143## GPU residency in RunMat (Do I need `gpuArray`?)
144You typically do **not** need to convert tensors manually. When the active provider implements the `ind2sub` hook (WGPU today), the entire conversion stays on the GPU. Otherwise RunMat gathers the inputs, performs validation on the host, and uploads the resulting subscript arrays back to the device so downstream kernels or fusion plans can continue using them without additional `gather` calls.
145
146## FAQ
147
148### What types does `ind2sub` accept for `idx`?
149Numeric and logical inputs are accepted. Logical arrays are treated as double precision (with `true → 1`, `false → 0`); complex values are rejected.
150
151### Can `siz` contain zeros?
152No. Each element of `siz` must be a positive integer, matching MATLAB behaviour.
153
154### How are errors reported for invalid indices?
155Indices that are non-integer, non-positive, or exceed `prod(siz)` raise errors matching MATLAB's wording (e.g., "Index exceeds number of array elements. Index must not exceed …").
156
157### What shape do the output arrays have?
158Each output array matches the shape of `idx`. Scalars produce scalar doubles; vectors remain vectors; higher-dimensional shapes are preserved.
159
160### Does `ind2sub` support GPU arrays?
161Yes. With the WGPU provider the conversion happens entirely on the GPU, including validation. Other providers gather the data to the host, compute the subscripts, and upload them back to the device automatically.
162
163### Can I request fewer outputs than dimensions?
164Yes. In a multi-output context RunMat provides as many outputs as requested in order, just like MATLAB. Any additional dimensions are still available inside the single-output cell if you need them.
165
166### How does `ind2sub` relate to `sub2ind`?
167They are inverse operations: `sub2ind` turns subscripts into linear indices, while `ind2sub` recovers subscripts from those linear indices.
168
169### Does `ind2sub` allocate full copies of the outputs?
170Yes. Each subscript array is materialised as a dense double array matching the shape of `idx`, mirroring MATLAB semantics.
171
172### What happens with empty inputs?
173Empty index arrays produce empty subscript arrays with the same shape.
174
175### Can I use `ind2sub` with more dimensions than two?
176Definitely—`ind2sub` works for any number of dimensions represented in `siz`.
177
178## See Also
179[sub2ind](./sub2ind), [size](../../introspection/size), [find](./find), [gpuArray](../../acceleration/gpu/gpuArray), [gather](../../acceleration/gpu/gather)
180
181## Source & Feedback
182- The full source code for the implementation of the `ind2sub` function is available at: [`crates/runmat-runtime/src/builtins/array/indexing/ind2sub.rs`](https://github.com/runmat-org/runmat/blob/main/crates/runmat-runtime/src/builtins/array/indexing/ind2sub.rs)
183- Found a bug or behavioural difference? Please [open an issue](https://github.com/runmat-org/runmat/issues/new/choose) with details and a minimal repro.
184"#;
185
186pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
187    name: "ind2sub",
188    op_kind: GpuOpKind::Custom("indexing"),
189    supported_precisions: &[ScalarType::F32, ScalarType::F64],
190    broadcast: BroadcastSemantics::Matlab,
191    provider_hooks: &[ProviderHook::Custom("ind2sub")],
192    constant_strategy: ConstantStrategy::InlineLiteral,
193    residency: ResidencyPolicy::NewHandle,
194    nan_mode: ReductionNaN::Include,
195    two_pass_threshold: None,
196    workgroup_size: None,
197    accepts_nan_mode: false,
198    notes: "WGPU provider executes `ind2sub` entirely on-device; other providers fall back to the host implementation and re-upload results to preserve residency.",
199};
200
201register_builtin_gpu_spec!(GPU_SPEC);
202
203pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
204    name: "ind2sub",
205    shape: ShapeRequirements::Any,
206    constant_strategy: ConstantStrategy::InlineLiteral,
207    elementwise: None,
208    reduction: None,
209    emits_nan: false,
210    notes: "Index conversion is eager and does not participate in fusion today.",
211};
212
213register_builtin_fusion_spec!(FUSION_SPEC);
214
215#[cfg(feature = "doc_export")]
216register_builtin_doc_text!("ind2sub", DOC_MD);
217
218#[runtime_builtin(
219    name = "ind2sub",
220    category = "array/indexing",
221    summary = "Convert MATLAB column-major linear indices into per-dimension subscript arrays.",
222    keywords = "ind2sub,linear index,subscripts,column major,gpu indexing",
223    accel = "custom"
224)]
225fn ind2sub_builtin(dims_val: Value, indices_val: Value) -> Result<Value, String> {
226    let (dims_value, dims_was_gpu) = materialize_value(dims_val)?;
227    let dims = parse_dims(&dims_value)?;
228    if dims.is_empty() {
229        return Err("Size vector must have at least one element.".to_string());
230    }
231
232    let total = total_elements(&dims)?;
233    let strides = build_strides(&dims)?;
234
235    if let Some(result) = try_gpu_ind2sub(&dims, &strides, total, &indices_val)? {
236        return Ok(result);
237    }
238
239    let (indices_value, indices_was_gpu) = materialize_value(indices_val)?;
240    let indices_tensor = tensor::value_into_tensor_for("ind2sub", indices_value)?;
241
242    let subscripts = compute_subscripts(&dims, total, &strides, &indices_tensor)?;
243
244    let want_gpu = (dims_was_gpu || indices_was_gpu) && runmat_accelerate_api::provider().is_some();
245
246    let mut outputs: Vec<Value> = Vec::with_capacity(dims.len());
247    for tensor in subscripts {
248        if want_gpu {
249            #[cfg(all(test, feature = "wgpu"))]
250            {
251                if runmat_accelerate_api::provider().is_none() {
252                    let _ = runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
253                        runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
254                    );
255                }
256            }
257            if let Some(provider) = runmat_accelerate_api::provider() {
258                let view = HostTensorView {
259                    data: &tensor.data,
260                    shape: &tensor.shape,
261                };
262                if let Ok(handle) = provider.upload(&view) {
263                    outputs.push(Value::GpuTensor(handle));
264                    continue;
265                }
266            }
267        }
268        outputs.push(tensor::tensor_into_value(tensor));
269    }
270
271    make_cell(outputs, 1, dims.len())
272}
273
274fn try_gpu_ind2sub(
275    dims: &[usize],
276    strides: &[usize],
277    total: usize,
278    indices: &Value,
279) -> Result<Option<Value>, String> {
280    #[cfg(all(test, feature = "wgpu"))]
281    {
282        if let Value::GpuTensor(h) = indices {
283            if h.device_id != 0 {
284                let _ = runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
285                    runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
286                );
287            }
288        }
289    }
290    let provider = match runmat_accelerate_api::provider() {
291        Some(p) => p,
292        None => return Ok(None),
293    };
294    if !provider.supports_ind2sub() {
295        return Ok(None);
296    }
297    let handle = match indices {
298        Value::GpuTensor(handle) => handle,
299        _ => return Ok(None),
300    };
301    if dims.len() != strides.len() {
302        return Err("Size vector must have at least one element.".to_string());
303    }
304    if dims.iter().any(|&d| d > u32::MAX as usize)
305        || strides.iter().any(|&s| s > u32::MAX as usize)
306        || total > u32::MAX as usize
307    {
308        return Ok(None);
309    }
310    let len = if handle.shape.is_empty() {
311        1usize
312    } else {
313        handle.shape.iter().copied().product()
314    };
315    if total == 0 && len > 0 {
316        return Err("Index exceeds number of array elements. Index must not exceed 0.".to_string());
317    }
318    if len > u32::MAX as usize {
319        return Ok(None);
320    }
321    let output_shape = if handle.shape.is_empty() {
322        vec![len, 1]
323    } else {
324        handle.shape.clone()
325    };
326    match provider.ind2sub(dims, strides, handle, total, len, &output_shape) {
327        Ok(handles) => {
328            if handles.len() != dims.len() {
329                return Err(
330                    "ind2sub: provider returned an unexpected number of outputs.".to_string(),
331                );
332            }
333            let values: Vec<Value> = handles.into_iter().map(Value::GpuTensor).collect();
334            make_cell(values, 1, dims.len()).map(Some)
335        }
336        Err(err) => Err(err.to_string()),
337    }
338}
339
340fn compute_subscripts(
341    dims: &[usize],
342    total: usize,
343    strides: &[usize],
344    indices: &Tensor,
345) -> Result<Vec<Tensor>, String> {
346    if strides.len() != dims.len() {
347        return Err("Size vector must have at least one element.".to_string());
348    }
349
350    let len = indices.data.len();
351    let mut outputs: Vec<Vec<f64>> = dims.iter().map(|_| Vec::with_capacity(len)).collect();
352
353    for &value in &indices.data {
354        let idx = coerce_linear_index(value, total)?;
355        let zero_based = idx - 1;
356        for (dim_index, (&dim, &stride)) in dims.iter().zip(strides.iter()).enumerate() {
357            let coord = ((zero_based / stride) % dim) + 1;
358            outputs[dim_index].push(coord as f64);
359        }
360    }
361
362    let output_shape = if indices.shape.is_empty() {
363        vec![len, 1]
364    } else {
365        indices.shape.clone()
366    };
367
368    let mut tensors = Vec::with_capacity(dims.len());
369    for data in outputs {
370        let tensor =
371            Tensor::new(data, output_shape.clone()).map_err(|e| format!("ind2sub: {e}"))?;
372        tensors.push(tensor);
373    }
374    Ok(tensors)
375}
376
377fn coerce_linear_index(value: f64, max_index: usize) -> Result<usize, String> {
378    if !value.is_finite() {
379        return Err("Linear indices must be positive integers.".to_string());
380    }
381    let rounded = value.round();
382    if (rounded - value).abs() > f64::EPSILON {
383        return Err("Linear indices must be positive integers.".to_string());
384    }
385    if rounded < 1.0 {
386        return Err("Linear indices must be positive integers.".to_string());
387    }
388    if rounded > usize::MAX as f64 {
389        return Err("Index exceeds maximum supported size for this platform.".to_string());
390    }
391    let coerced = rounded as usize;
392    if coerced > max_index {
393        return Err(format!(
394            "Index exceeds number of array elements. Index must not exceed {}.",
395            max_index
396        ));
397    }
398    Ok(coerced)
399}
400
401#[cfg(test)]
402mod tests {
403    use super::*;
404    use crate::builtins::common::test_support;
405    use runmat_accelerate_api::HostTensorView;
406    use runmat_builtins::{Tensor, Value};
407
408    fn cell_to_vec(cell: &runmat_builtins::CellArray) -> Vec<Value> {
409        cell.data.iter().map(|ptr| (**ptr).clone()).collect()
410    }
411
412    #[test]
413    fn recovers_matrix_indices() {
414        let dims = Tensor::new(vec![3.0, 4.0], vec![1, 2]).unwrap();
415        let result = ind2sub_builtin(Value::Tensor(dims), Value::Num(8.0)).unwrap();
416        match result {
417            Value::Cell(cell) => {
418                let values = cell_to_vec(&cell);
419                assert_eq!(values.len(), 2);
420                assert_eq!(values[0], Value::Num(2.0));
421                assert_eq!(values[1], Value::Num(3.0));
422            }
423            other => panic!("expected cell output, got {other:?}"),
424        }
425    }
426
427    #[test]
428    fn handles_vector_indices() {
429        let dims = Tensor::new(vec![3.0, 5.0], vec![1, 2]).unwrap();
430        let idx = Tensor::new(vec![7.0, 8.0, 9.0], vec![1, 3]).unwrap();
431        let result =
432            ind2sub_builtin(Value::Tensor(dims), Value::Tensor(idx)).expect("ind2sub result");
433        match result {
434            Value::Cell(cell) => {
435                let values = cell_to_vec(&cell);
436                assert_eq!(values.len(), 2);
437                match &values[0] {
438                    Value::Tensor(t) => {
439                        assert_eq!(t.shape, vec![1, 3]);
440                        assert_eq!(t.data, vec![1.0, 2.0, 3.0]);
441                    }
442                    other => panic!("expected tensor rows, got {other:?}"),
443                }
444                match &values[1] {
445                    Value::Tensor(t) => {
446                        assert_eq!(t.shape, vec![1, 3]);
447                        assert_eq!(t.data, vec![3.0, 3.0, 3.0]);
448                    }
449                    other => panic!("expected tensor cols, got {other:?}"),
450                }
451            }
452            other => panic!("expected cell output, got {other:?}"),
453        }
454    }
455
456    #[test]
457    fn recovers_three_dimensional_indices() {
458        let dims = Tensor::new(vec![2.0, 3.0, 4.0], vec![1, 3]).unwrap();
459        let idx = Tensor::new(vec![3.0, 11.0], vec![1, 2]).unwrap();
460        let result =
461            ind2sub_builtin(Value::Tensor(dims), Value::Tensor(idx)).expect("ind2sub result");
462        if let Value::Cell(cell) = result {
463            let values = cell_to_vec(&cell);
464            assert_eq!(values.len(), 3);
465            assert_eq!(
466                values[0],
467                Value::Tensor(Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap())
468            );
469            assert_eq!(
470                values[1],
471                Value::Tensor(Tensor::new(vec![2.0, 3.0], vec![1, 2]).unwrap())
472            );
473            assert_eq!(
474                values[2],
475                Value::Tensor(Tensor::new(vec![1.0, 2.0], vec![1, 2]).unwrap())
476            );
477        } else {
478            panic!("expected cell output");
479        }
480    }
481
482    #[test]
483    fn errors_on_out_of_range_index() {
484        let dims = Tensor::new(vec![3.0, 4.0], vec![1, 2]).unwrap();
485        let err =
486            ind2sub_builtin(Value::Tensor(dims), Value::Num(13.0)).expect_err("expected failure");
487        assert!(
488            err.contains("Index exceeds number of array elements"),
489            "unexpected error: {err}"
490        );
491    }
492
493    #[test]
494    fn errors_on_zero_index() {
495        let dims = Tensor::new(vec![3.0, 4.0], vec![1, 2]).unwrap();
496        let err =
497            ind2sub_builtin(Value::Tensor(dims), Value::Num(0.0)).expect_err("expected failure");
498        assert!(
499            err.contains("Linear indices must be positive integers"),
500            "unexpected error: {err}"
501        );
502    }
503
504    #[test]
505    fn errors_on_fractional_index() {
506        let dims = Tensor::new(vec![3.0, 4.0], vec![1, 2]).unwrap();
507        let err =
508            ind2sub_builtin(Value::Tensor(dims), Value::Num(2.5)).expect_err("expected failure");
509        assert!(
510            err.contains("Linear indices must be positive integers"),
511            "unexpected error: {err}"
512        );
513    }
514
515    #[test]
516    fn errors_on_invalid_size_elements() {
517        let dims = Tensor::new(vec![3.5, 4.0], vec![1, 2]).unwrap();
518        let err = ind2sub_builtin(Value::Tensor(dims), Value::Num(5.0)).expect_err("expected fail");
519        assert!(
520            err.contains("Size arguments must be positive integers"),
521            "unexpected error: {err}"
522        );
523    }
524
525    #[test]
526    fn ind2sub_gpu_roundtrip() {
527        test_support::with_test_provider(|provider| {
528            let dims = Tensor::new(vec![3.0, 4.0], vec![1, 2]).unwrap();
529            let idx_tensor = Tensor::new(vec![10.0, 11.0], vec![2, 1]).unwrap();
530            let view = HostTensorView {
531                data: &idx_tensor.data,
532                shape: &idx_tensor.shape,
533            };
534            let handle = provider.upload(&view).expect("upload indices");
535            let result = ind2sub_builtin(Value::Tensor(dims), Value::GpuTensor(handle)).unwrap();
536            match result {
537                Value::Cell(cell) => {
538                    let values = cell_to_vec(&cell);
539                    assert_eq!(values.len(), 2);
540                    match &values[0] {
541                        Value::GpuTensor(_) => {}
542                        other => panic!("expected gpu tensor output, got {other:?}"),
543                    }
544                    match &values[1] {
545                        Value::GpuTensor(_) => {}
546                        other => panic!("expected gpu tensor output, got {other:?}"),
547                    }
548                    let rows = test_support::gather(values[0].clone()).expect("gather rows");
549                    assert_eq!(rows.shape, vec![2, 1]);
550                    assert_eq!(rows.data, vec![1.0, 2.0]);
551                    let cols = test_support::gather(values[1].clone()).expect("gather cols");
552                    assert_eq!(cols.shape, vec![2, 1]);
553                    assert_eq!(cols.data, vec![4.0, 4.0]);
554                }
555                other => panic!("expected cell output, got {other:?}"),
556            }
557        });
558    }
559
560    #[test]
561    #[cfg(feature = "doc_export")]
562    fn doc_examples_present() {
563        let blocks = test_support::doc_examples(DOC_MD);
564        assert!(!blocks.is_empty());
565    }
566
567    #[test]
568    #[cfg(feature = "wgpu")]
569    fn ind2sub_wgpu_matches_cpu() {
570        let provider_init = std::panic::catch_unwind(|| {
571            runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
572                runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
573            )
574        });
575        if let Ok(Ok(_)) = provider_init {
576            // provider successfully registered
577        } else {
578            return;
579        }
580
581        let dims_tensor = Tensor::new(vec![3.0, 4.0], vec![1, 2]).unwrap();
582        let idx_tensor = Tensor::new(vec![7.0, 8.0, 9.0], vec![1, 3]).unwrap();
583
584        let cpu = ind2sub_builtin(
585            Value::Tensor(dims_tensor.clone()),
586            Value::Tensor(idx_tensor.clone()),
587        )
588        .expect("cpu ind2sub");
589
590        let provider = runmat_accelerate_api::provider().unwrap();
591        let view = HostTensorView {
592            data: &idx_tensor.data,
593            shape: &idx_tensor.shape,
594        };
595        let handle = provider.upload(&view).expect("upload indices");
596
597        let gpu = ind2sub_builtin(Value::Tensor(dims_tensor), Value::GpuTensor(handle))
598            .expect("gpu ind2sub");
599
600        let cpu_values = match cpu {
601            Value::Cell(cell) => cell_to_vec(&cell),
602            other => panic!("expected cell output, got {other:?}"),
603        };
604        let gpu_values = match gpu {
605            Value::Cell(cell) => cell_to_vec(&cell),
606            other => panic!("expected cell output, got {other:?}"),
607        };
608
609        assert_eq!(cpu_values.len(), gpu_values.len());
610
611        for (cpu_val, gpu_val) in cpu_values.iter().zip(gpu_values.iter()) {
612            let host_cpu = test_support::gather(cpu_val.clone()).expect("gather cpu");
613            let host_gpu = test_support::gather(gpu_val.clone()).expect("gather gpu");
614            assert_eq!(host_cpu.shape, host_gpu.shape);
615            assert_eq!(host_cpu.data, host_gpu.data);
616        }
617    }
618}