Indexing arrays by an array of logicals

Language enhancement
1

Imagine that you have two arrays A and B of the same shape, and you would like to set A to 0 wherever B is positive. This can be done by

where (B > 0)
      A = 0
end where

This is great, but in MATLAB, it can be done in a more concise way, which is

A(B > 0) = 0

Here, A is indexed by an array of logicals (masks). Operations are affecting only the positions where the mask is .true.. This seems to me a very convenient feature. It might be something worth considering in Fortran as well.

See here for more information about logical indexing in MATLAB: Matrix Indexing in MATLAB - MATLAB & Simulink

Python has a similar feature:
https://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/AdvancedIndexing.html

4 Likes
2

There are a few ways more ways to achieve the same:

do concurrent (i = 1:size(b), b(i) > 0)
  a(i) = 0 
end do

a(findloc(b > 0,.true.)) = 0

I can imagine a small helper function would make the second version shorter, but perhaps less performant. I guess it doesn’t generalize to arbitrary rank either.

5 Likes
3

If your only purpose is to make the code look nicer and you don’t need elsewhere, nothing prevent you from using

where (b > 0) a = 0

This is standard approval (F2008 7.2.3.1), and to me this is closer to natural language than a(b > 0) = 0.

8 Likes
4

In Fortran, array indices have to be integers; therefore, “an array of logicals” cannot be used to index an array (of any type).

In the proposed assignment

A(B > 0) = 0

what is the type of the expression B > 0 ? Is it an array of integers or is it an array of logical values?

The following short program illustrates the issues. Many ideas may seem “smart” in a narrow context or to offer convenience, but (i) will they work in all contexts allowed by the language, and (ii) will the suggested changes in semantics break existing code?

program WhereWhat
implicit none
integer, parameter :: N = 5
integer :: ix(N), A(N), B(N)
logical :: gt(N)
!
B = [2, 3, 4, 5, 7]
A = B
where (B > 4) A = -A
print 10,'A : ',A
print *,'(B > 4) :      ',(B > 4) ! what is the type of the expression?
!ix = B > 4     ! compile time error, logical array value placed into integer array
gt = B > 4
print 10,'ix = ',ix
print 20,'gt = ',gt
print 10,'A(ix) = ',A(ix)  ! illegal subscripts used, results meaningless
!print *,'A(gt) = ',A(gt) ! compile time error: logical values used as subscripts
10 format(1x,A,T15,5I4)
20 format(1x,A,T15,5L4)
end program
1 Like
5

In MATLAB this would be an array of logical values (represented as zeros and ones).

MATLAB also has a function find which can be used to map logical values (or nonzero elements more generally) to linear indexes. The syntax is:

k = find(X)         % or
[row,col] = find(X)

This could be emulated to some extent in Fortran, but would be more verbose than the where construct, or the use of a logical array as an index, like in MATLAB.

2 Likes
6

Yes, the value of the logical expression B > 0 as represented in Matlab is [1, 1, 1, 0, 0], so we cannot use that expression as an index expression. The 0 indices are less than the lower bound of the array, and the only other value is 1, and the result is not at all what we want. That is the problem: we need to implicitly map the positions that have 1 as the value into a new integer valued index array, [1, 2, 3], and use those values as the indices into A. Furthermore, this would have to be done consistently in all contexts, since in Fortran the value of an expression has an associated type.

1 Like
7

I have written a function that returns a vector index from a logical array

function true_pos(tf) result(ipos)
! return in ipos the positions of the true elements in tf(:)
logical, intent(in) :: tf(:)
integer             :: ipos(count(tf))
integer             :: i,j
j = 0
do i=1,size(tf)
   if (tf(i)) then
      j = j + 1
      ipos(j) = i
   end if
end do
end function true_pos

and would write

A(true_pos(B > 0)) = 0

Logical indexing would be convenient, but given a simple alternative I’m not sure it is worth changing the language. If you are referring to array elements but not setting them you can use the expression

pack(A,B>0)

3 Likes
8

@zaikunzhang thanks for bringing this up. Indeed, Python/NumPy also allows this.

There are always usually some reasons why it wasn’t done in Fortran, as others have pointed out above. It might also be very difficult to make such changes now and keep everything backwards compatible.

However, I still like to have these discussions and to brainstorm such possibilities. If nothing else, one could for example use these to create better compiler error messages, say if you write A(B>0) = 0, the compiler can tell you to write it as where (B>0) A = 0. But I think these discussions can also be useful either to try to improve Fortran, or if somebody is designing a new language, how it could be done.

3 Likes
9

@Beliavsky, I think this would make a nice addition to the standard library. Could you open a new issue for it?

10

Some care is required to make the proposed function work correctly when its argument is not an array with a lower bound of 1, as the following test code demonstrates (it uses Beliavsky’s true_pos function from #7).

program xtpos
implicit none
integer :: ia(2:6) = [1,2,3,4,5], i
!
write(*,'(5(1x,I3))')(i,i=lbound(ia,1),ubound(ia,1))
do i=lbound(ia,1),ubound(ia,1)
   write(*,'(1x,L3)',advance='no')ia(i) > 2
end do
print *
print 10, 'True_Pos',true_pos(ia > 2)
print 10, 'Expected',[4,5,6]
10 format(A8,2x,3(1x,I3))
contains

function true_pos(tf) result(ipos)
! return in ipos the positions of the true elements in tf(:)
logical, intent(in) :: tf(:)
integer             :: ipos(count(tf))
integer             :: i,j
j = 0
do i=1,size(tf)
   if (tf(i)) then
      j = j + 1
      ipos(j) = i
   end if
end do
end function true_pos

end program

The output:

   2   3   4   5   6
   F   F   T   T   T
True_Pos     3   4   5
Expected     4   5   6

A related issue: what should the function return when an array section or an array reference with a vector subscript is the input argument?

2 Likes
11

What? None of you suggested:

a = merge (a,0,b>0)

???

6 Likes
12

Nice! I don’t use merge in my codes, so I haven’t thought of that.

13

When you spend too much time merging pull requests, you forget about the Fortran merge. :slight_smile:

1 Like
14

I generally would not recommend MERGE for any instructions except for those involving scalar arguments. My experience is both the practitioners and compiler implementations don’t quite “get” this intrinsic.

Exhibit A: going by the original post, it should be

a = merge (0, a, b>0)

Exhibit B: compiler implementations can do needless operations.

C:\Temp>a.exe
Block 1: use MERGE
In assign_w_t: performing lhs%n(currently= 1 ) = 1 assignment
In assign_w_t: performing lhs%n(currently= 2 ) = 0 assignment
In assign_w_t: performing lhs%n(currently= 3 ) = 3 assignment
In assign_w_t: performing lhs%n(currently= 4 ) = 4 assignment
In assign_w_t: performing lhs%n(currently= 5 ) = 5 assignment
In assign_w_t: performing lhs%n(currently= 6 ) = 6 assignment
In assign_w_t: performing lhs%n(currently= 7 ) = 7 assignment
In assign_w_t: performing lhs%n(currently= 8 ) = 0 assignment
In assign_w_t: performing lhs%n(currently= 9 ) = 9 assignment
In assign_w_t: performing lhs%n(currently= 10 ) = 10 assignment
a = 1 0 3 4 5 6
7 0 9 10

Block 2: use WHERE
In assign_w_t: performing lhs%n(currently= 2 ) = 0 assignment
In assign_w_t: performing lhs%n(currently= 8 ) = 0 assignment
a = 1 0 3 4 5 6
7 0 9 10

C:\Temp>

Click for code!

Illustrate array operations using a derived type

module w_m
   type :: w_t
      integer :: n = 0
   end type
   generic :: operator(>) => compare_w_t
   generic :: assignment(=) => assign_w_t
   type(w_t), parameter :: ZERO = w_t()
contains
   elemental function compare_w_t( lhs, rhs ) result(r)
      type(w_t), intent(in) :: lhs
      type(w_t), intent(in) :: rhs
      logical :: r
      r = lhs%n > rhs%n 
   end function
   impure elemental subroutine assign_w_t( lhs, rhs )
      type(w_t), intent(inout) :: lhs
      type(w_t), intent(in)    :: rhs
      print *, "In assign_w_t: performing lhs%n(currently=", lhs%n, ") = ", rhs%n, " assignment" 
      lhs%n = rhs%n
   end subroutine 
end module
   use w_m
   type(w_t), dimension(10) :: a, b
   blk1: block
      print *, "Block 1: use MERGE"
      call init()
      a = merge( ZERO, a, b > ZERO )  !<-- use of MERGE intrinsic
      print *, "a = ", a%n
   end block blk1
   print * 
   blk2: block
      print *, "Block 2: use WHERE"
      call init()
      where ( b > ZERO ) a = ZERO   !<-- WHERE statement
      print *, "a = ", a%n
   end block blk2
contains
   subroutine init()
      b(2)%n = 1 ; b(8)%n = 1
      do concurrent ( integer :: i=1:size(a) )
         a(i)%n = i
      end do
   end subroutine 
end

C:\Temp>ifort /standard-semantics a.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.3.0 Build 20210609_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.29.30038.1
Copyright (C) Microsoft Corporation. All rights reserved.

-out:a.exe
-subsystem:console
a.obj

C:\Temp>

3 Likes
15

It is true that MERGE is not well optimized, but it IS the language feature being asked for and if more people used it, compiler writers would work harder to optimize it. I think Cray Fortran already does this to some degree.

One could have said similar things about array assignments in the early 90s.

2 Likes
16

No, in the context of this thread, what is being asked for is some sort of “precis writing” for what the language offers with where ( b > 0 ) a = 0 and that ain’t MERGE.

2 Likes
17

To be fair, merge is one way it could be done (and let’s assume compilers implement this well in the future), but as beautifully shown by you above, the merge syntax is so confusing that anybody can get it wrong. I personally do not use it for the same reason, I just never found it clear enough what it does to warrant using it. where is pretty good. So is the NumPy syntax a[b>0] = 0.

1 Like
18

Thank everyone for the ideas and insights, from which I have learned a lot, particularly the usage of merge proposed by Dr. Fortran @sblionel.

It seems that the idea of logical indexing has gained little popularity during the discussion. There do exist a few languages for numerical computing, namely Python (as mentioned by @certik), Julia, R, and MATLAB, that support logical indexing. So I do not think it is just yet another “smart idea”.

However, I do understand that including logical indexing necessities significant augmentation to the existing indexing mechanism of Fortran, and I do not assume by any means that such augmentation is trivial to implement.

Meanwhile, it will be very interesting to see a snippet illustrating that logical indexing can hurt backward compatibility, namely a piece of code that works now but would become broken if logical indexing were introduced. Mathematically, it is a bit hard for me to imagine how an augmentation/generalization (rather than a modification) can harm backward compatibility.

It may be worth mentioning that A(B>0) = 0 is only a simple illustration for logical indexing. I totally agree that where as suggested by @han190 is a better implementation for this job, because where is closer to the natural language. However, logical indexing can be applied in much broader contexts, for instance, A(B>0) = C(D<0). The true_pos proposed by @Beliavsky is very nice, and it is essentially the same as logical indexing except that it is implemented as a function.

Anyway, even if logical indexing were to be included in the language by a future standard, I would not expect my projects to make use of it. I am quite conservative regarding new features in order to keep my code as portable as possible. I am still refraining from some convenient F2008 features because they are not supported yet by some recent compilers; for example, Absoft 21.0 does not support automatic allocation upon intrinsic assignments. One may complain that the compiler vendors are too slow to support the standards, but one can never blame users for using such compilers.

2 Likes
19

I do not find merge at all confusing but then I use it with high frequency. I find many things confusing and needing a visit to help for the sole reason of infrequent use. Many simple if /else/ endif blocks become one line instead five lines! It saves typing and makes the code easier to scan read IMO.

3 Likes
20

In MoA,
(
(((rav A) = (rav B))/(iota (x red ( rho A)))) Omega (<0 1>) <2 4>
) psi A = 0

So when psi reduced, we just do the assignment.
The rav A = rav B Creates a Boolean vector, when compressed with iota product rho the shape, finds the index, then using Omega, applies gamma prime to get the indices in the two by 4 array , eg. Now we assign 0 to them.