atmtools / arts

The Atmospheric Radiative Transfer Simulator
https://www.radiativetransfer.org/
Other
58 stars 28 forks source link

New types for ARTS-3 #612

Open riclarsson opened 1 year ago

riclarsson commented 1 year ago

Hi all,

We should decide the future of the use of a type system in ARTS. There's numerous ideas of how to improve it in a few PRs so I thought we gather those here.

I will present only the problem I have with some of our use of matrices and vectors here, but I know @erikssonpatrick also is considering exposing Complex (and perhaps a few multi-dimensional access to it?) as yet-another-type for the users.

But we should also decide if we want the solution below, or if it would be better for LOS and POS to be the types, and Ppath to contain lists of these.


Vector and Matrix Limitations

The problem: The propagation matrix calculations requires an Agenda call, so its input line-of-sight has to be a Vector today. But Ppath can only give VectorView from the Matrix it holds, so the data has to be copied. VectorView is not and should not be a user type. On the other hand, using ArrayOfVector in Ppath would make the data non-exhaustive (a Vector has a pointer, so it is effectively an Array<> of pointers) which means it would be less efficient to work with Ppath. It must have been decided in the past that this is not a good trade-off to make and that Ppath must have exhaustive data.

The suggested solution: Expose some fixed-size multidimensional arrays to the user. Most notably, Vector2 and Vector3, which can be used as ArrayOfVector2 and ArrayOfVector3 in Ppath. These can then be accessed without the need of copying, and the data will be exhaustive in Array<> form because a VectorX does not need a pointer if implemented as matpack::matpack_constant_data<Numeric, X>.

New user types: Vector2; Vector3; ArrayOfVector2; ArrayOfVector3

erikssonpatrick commented 1 year ago

Hej, As usual, we approach things from different ends. I care more about the user perspective. Will this make ARTS easier to use? That said, there could be competing aspects. Maybe a more clear definition, but makes things more complicated practically. Not many users will look into the path structure. The standard user will mainly meet this when defining sensor_pos and sensor_los. Is it a good idea to convert these WSVs to AoVector3 and AoVector2, respectively? For example, it will be harder to create these WSVs with linspace and similar functions. I think is a quite general situation. Introducing special types gives clearer definitions of variables, but makes it harder to use general methods to create them. What's most important?

riclarsson commented 1 year ago

Yes, we do have a different approaches. I disagree with what I get from your comment, which is that you seem to think that we have so much technical debt in ARTS that we cannot make these improvements to user-space.

It is pretty obvious to me that an error in ws.sensor_pos = [[1, 2]] is better than an error in ws.yCalc(). The former is the place where the error takes place, the latter where the error has consequences. To me, we should try and design ARTS to catch errors as soon as possible, not at the last minute. I therefore am in favor of just using these fixed-sized types wherever the size is fixed.

I don't think we have a user-facing linspace for Matrix, so I do not understand how this could be a concern for the users?

Other than that, it is pretty trivial to do linspace for fixed-sized variables. For example this gives a linearly-spaced sample of VectorX with n array-entries starting at pl, ending at ph with linear-spacing in-between:

template <Index N> using FixedVector = matpack::matpack_constant_data<Numeric, N>;

template <Index N>
void linspace(Array<FixedVector<N>>& arr, FixedVector<N>&& pl, FixedVector<N>&& ph, Index n) {
  ARTS_USER_ERROR_IF(n < 2, "Must have more than 1 entry")
  arr.resize(n);
  FixedVector<N> dp{ph};
  dp -= pl;
  dp /= static_cast<Numeric>(n - 1);
  std::generate(arr.begin(), arr.end() - 1,
                [x = pl, dx = dp]() mutable {
                  FixedVector<N> out{x};
                  x += dx;
                  return out;
                });
  arr.back() = ph;
}

(And since it is a template, you can use it with Vector3 or Vector2 as you please.)

nuitlejour commented 3 weeks ago

I want to offer a bold sugguestion, why not use Julia language to develop a brand new RTM for microwave region? successful example like vSmartMOM.jl, https://github.com/RemoteSensingTools/vSmartMOM.jl