Add view_re_im and view_mut_re_im methods

[jturner314]

These methods make it possible to obtain views of the real and imaginary components of elements in an array of complex elements. For example:

use ndarray::prelude::*;
use num_complex::{Complex, Complex64};

let mut arr = array![
    [Complex64::new(1., 2.), Complex64::new(3., 4.)],
    [Complex64::new(5., 6.), Complex64::new(7., 8.)],
    [Complex64::new(9., 10.), Complex64::new(11., 12.)],
];

let Complex { mut re, mut im } = arr.view_mut_re_im();
assert_eq!(re, array![[1., 3.], [5., 7.], [9., 11.]]);
assert_eq!(im, array![[2., 4.], [6., 8.], [10., 12.]]);

re[[0, 1]] = 13.;
im[[2, 0]] = 14.;

assert_eq!(arr[[0, 1]], Complex64::new(13., 4.));
assert_eq!(arr[[2, 0]], Complex64::new(9., 14.));

I've created this PR as a draft because I'm not sure where to put these methods in the code and because I need to add more tests. Also, what do you think of the names of the methods? Perhaps a better design would be split_re_im methods for ArrayView and ArrayViewMut.

Fwiw, NumPy has more general implementation of this type of functionality with arbitrary structs as elements instead of just complex numbers. In principle, we could do the same thing. We'd need a trait to constrain the allowable element types, and if we wanted to support all the structs NumPy can, we'd need to change our strides to be in units of bytes rather than in units of elements. Implementing that would be a lot of work, though, and this simple support for complex numbers is useful today.

[bluss]

Fantastic. I wonder if there should be a constraint on T. We have ensured that Complex is repr(C). Any other constraints on T needed?

The fun return type of the method inspires us to think that what happens if the user inputs an Array<Complex<ArrayView<_, _>>> here? ArrayView are Copy but would it truly be a defined operation?

[jturner314]

Fantastic. I wonder if there should be a constraint on T. We have ensured that Complex is repr(C). Any other constraints on T needed?

I don't think so, but the guidelines aren't very clear on this. The relevant references here are the nomicon and the unsafe code guidelines. The re field is guaranteed to be at offset 0 in the struct, and the im field is guaranteed to be after it, as described by the C17 standard (which is referenced by the unsafe code guidelines):

Within a structure object, the non-bit-field members and the units in which bit-fields reside have addresses that increase in the order in which they are declared. A pointer to a structure object,suitably converted, points to its initial member (or if that member is a bit-field, then to the unit in which it resides), and vice versa. There may be unnamed padding within a structure object, but not at its beginning.

The possibility of padding in Complex<T> is somewhat unclear. Maybe, to be on the safe side, we should add a few asserts (which should be eliminated by the compiler at compile time):

assert_eq!(mem::size_of::<Complex<T>>(), mem::size_of::<T>().checked_mul(2).unwrap());
assert_eq!(mem::align_of::<Complex<T>>(), mem::align_of::<T>());
// Non-empty case only:
let ptr = self.ptr.as_ptr();
assert_eq!(
    unsafe { &(*ptr).im as *const T as usize }
        .checked_sub(ptr as usize)
        .unwrap(),
    mem::size_of::<T>(),
);

The fun return type of the method inspires us to think that what happens if the user inputs an Array<Complex<ArrayView<_, _>>> here? ArrayView are Copy but would it truly be a defined operation?

Yeah, that would be a weird type. Arrays containing arrays are pretty weird. I don't see any soundness issues with it, though. I think the implementation should work for arbitrary T. The return type of Array::<Complex<ArrayView<'a, _, _>, D>::view_re_im would be Complex<ArrayView<'b, ArrayView<'a, _, _>, D>>. In other words, you'd get a Complex struct where the re and im fields would each be a view containing views.

[bluss]

Why shouldn't the interface just be something like the methods fn re(&self) -> ArrayView<f64, ..> and so on? Just to consider an alternative. (I think I would prefer if Complex had fields called real/imag instead of re/im.)

This is an interesting feature, the hard part is just how to make it available to all custom structs. I think that the rustic solution would be a derive proc macro that creates and implements a custom extension trait:

#[derive(NdarrayStructFields)]
struct Foo {
    a: i32,
    b: f32,
}

/* adds and implements the following trait for arrays */
trait NdarrayStructFields_Foo {
    fn view_a(&self) -> ArrayView<i32, ..>;
    fn view_b(&self) -> ArrayView<f32, ..>;
    /* view_mut_a, ... */
}

It would depend on a stable offsetof operation, I guess. A macro could be added as a 90% solution instead of a proc macro. I don't foresee that it's something we implement ourselves if we don't want to (can be external).

[jturner314]

Why shouldn't the interface just be something like the methods fn re(&self) -> ArrayView<f64, ..> and so on?

That would be fine for the immutable case but would be less versatile for the mutable case, since the user wouldn't be able to get mutable views of the real and imaginary parts at the same time. Maybe it would be better to call these methods split_re_im and implement them on ArrayView/ArrayViewMut instead of &ArrayBase/&mut ArrayBase? That would have the additional benefit of preserving the lifetime of the original view in the split views.

This is an interesting feature, the hard part is just how to make it available to all custom structs. I think that the rustic solution would be a derive proc macro that creates and implements a custom extension trait:

Right, that would be my thinking too. The soa_derive crate provides somewhat related functionality for slices/Vecs.

An alternative design would be a more dynamic approach like this:

pub unsafe trait FieldAccess<X> {
    // Returns the offset in the struct of the field with the specified name.
    fn field_offset(&self, name: &str) -> Result<usize, FieldAccessError>;
}

pub enum FieldAccessError {
    // There are no fields with the specified name in the struct.
    MissingField,
    // The specified type of the field is incorrect.
    WrongFieldType,
    // The size of the struct is not evenly divisible by the size of the field.
    //
    // This is a requirement because strides in `ndarray` are in units of the element type,
    // so the stride in units of bytes must be evenly divisible by the size of the element.
    StructSizeNotMultipleOfFieldSize,
}

impl<A, S, D> ArrayBase<S, D>
where
    S: Data<Elem = A>,
    D: Dimension,
{
    pub fn view_field<X>(&self, name: &str) -> Result<ArrayView<'_, X, D>, FieldAccessError>
    where
        A: FieldAccess<X>,
    { ... }
}

(The user's struct would implement FieldAccess<X> for each type X of its fields.) I'm not sure there's much benefit to more dynamic access like this, though.

An external crate could provide all this functionality if we expose a way to construct a view with arbitrary strides (including negative strides).

For now, though, providing this functionality for Complex is useful even if we don't provide it for arbitrary structs.

[bluss]

Yes, using the split idea seems good, should have mostly upsides. Only hard part regardless is discoverability.

[bluss]

I like the approach in this PR, so I'm fine with including this when it's ready.

Complex<T> is a special case of structs that look like packed vectors - like [T; N] for small N. Thus we have the alternative of viewing that kind of data as an array view with an extra inner axis of dimension N. This is something ndarray could implement for various such vectors.

[jturner314]

I've pushed a commit which switches from view_re_im to split_re_im. I still need to write tests.

Unfortunately, rustdoc documents the methods incorrectly on the ArrayBase page of the docs. For instance, it documents

impl<'a, T, D> ArrayView<'a, Complex<T>, D>
where
    D: Dimension,
{
    pub fn split_re_im(self) -> Complex<ArrayView<'a, T, D>> {
        ...
    }
}

as

impl<'a, T, D> ArrayBase<ViewRepr<&'a A>, D>
where
    D: Dimension,
{
    pub fn split_re_im(self) -> Complex<ArrayView<'a, T, D>> {
        ...
    }
}

I spent a little time trying to write a minimal example to demonstrate this issue in order to report it to rustdoc, but haven't finished. For instance, the following example doesn't have the issue:

pub struct ArrayBase<S, D: Dimension>(S, D);

pub struct ViewRepr<T>(T);

pub trait Dimension {}

type ArrayView<'a, A, D> = ArrayBase<ViewRepr<&'a A>, D>;

mod foo {
    use crate::{ArrayView, Dimension};
    use num_complex::Complex;

    impl<'a, A, D> ArrayView<'a, A, D>
    where
        D: Dimension,
    {
        pub fn bar(self) {}
    }

    impl<'a, T, D> ArrayView<'a, Complex<T>, D>
    where
        D: Dimension,
    {
        pub fn foo(self) -> Complex<ArrayView<'a, T, D>> {
            unimplemented!()
        }
    }
}

Anyway, I don't think the rustdoc issue is necessarily a blocker for this PR. The remaining piece is writing tests.