From 60be54fcc94e4b20b3fd99c4260631f91e3a41da Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Wed, 4 Oct 2023 17:00:45 +0200
Subject: [PATCH] Optimize routines to test point inclusion in closed contour

---
 src/reflections.f90 | 20 +++++++++++++++-----
 src/utils.f90       | 39 ++++++++++++++++++++++++++++++++++++---
 2 files changed, 51 insertions(+), 8 deletions(-)

diff --git a/src/reflections.f90 b/src/reflections.f90
index 7cc157b..3b3c12d 100644
--- a/src/reflections.f90
+++ b/src/reflections.f90
@@ -218,22 +218,32 @@ function inside(xc,yc,n,x,y)
   real(wp_), intent(in) :: x,y
   logical :: inside
   integer, dimension(n) :: jint
-  real(wp_), dimension(n) :: xint
+!   real(wp_), dimension(n) :: xint
   real(wp_), dimension(n+1) :: xclosed,yclosed
   integer :: i,nj
   xclosed(1:n)=xc(1:n)
   yclosed(1:n)=yc(1:n)
   xclosed(n+1)=xc(1)
   yclosed(n+1)=yc(1)
+  ! Count the number of segments for which y lies between
+  ! yclosed(j) and yclosed(j+1)
+  ! and store the indexes j in jint
   call locate_unord(yclosed,n+1,y,jint,n,nj)
   inside=.false.
   if (nj==0) return
+  ! Determine the x coordinate, along the previously identified segments,
+  ! required to match the y coordinate of the inquired point,
+  ! and count the segments to the left of the inquired point.
+  ! An even number means that the point is outside the closed contour.
   do i=1,nj
-    xint(i)=intlinf(yclosed(jint(i)),xclosed(jint(i)), &
-                   yclosed(jint(i)+1),xclosed(jint(i)+1),y)
+    if (x>intlinf(yclosed(jint(i)),xclosed(jint(i)), &
+                  yclosed(jint(i)+1),xclosed(jint(i)+1),y)) &
+      inside=.not.inside
+!     xint(i)=intlinf(yclosed(jint(i)),xclosed(jint(i)), &
+!                    yclosed(jint(i)+1),xclosed(jint(i)+1),y)
   end do
-  call bubble(xint,nj)
-  inside=(mod(locatef(xint,nj,x),2)==1)
+  ! call bubble(xint,nj)
+  ! inside=(mod(locatef(xint,nj,x),2)==1)
 end function inside
 
 
diff --git a/src/utils.f90 b/src/utils.f90
index a2cc142..39a8449 100644
--- a/src/utils.f90
+++ b/src/utils.f90
@@ -93,11 +93,44 @@ contains
     real(wp_), intent(in) :: x
     integer, dimension(m), intent(inout) :: j
     integer :: i
+    logical :: larger_than_last
+
     nj=0
-    do i=1,n-1
-      if (x>a(i).neqv.x>a(i+1)) then
+!     do i=1,n-1
+!       if (x>a(i).neqv.x>a(i+1)) then
+!         nj=nj+1
+!         if (nj<=m) j(nj)=i
+!       end if
+!     end do
+!
+! Alternative formulation
+! Should reduce the number of evaluations of i+1 and x>a(i)
+! * For an array a with only two changes of derivative
+!   (e.g., z of convex R,z contour) and x in range:
+!     * comparisons between floats (x>a(i)): 2*(n-1) --> n+2
+!     * comparisons between booleans (.neqv.): n-1 --> n-1
+!     * comparisons between integers (n<2): 0 --> 1
+!     * sum of integers (i+1, i-1): n --> 2
+!     * booleans assigments: 0 --> 3
+! * For x out of range:
+!     * comparisons between floats (x>a(i)): 2*(n-1) --> n
+!     * comparisons between booleans (.neqv.): n-1 --> n-1
+!     * comparisons between integers (n<2): 0 --> 1
+!     * sum of integers (i+1, i-1): n --> 0
+!     * booleans assigments: 0 --> 1
+! * For a elements alternated above/below x:
+!     * comparisons between floats (x>a(i)): 2*(n-1) --> 2*(n-1)+1
+!     * comparisons between booleans (.neqv.): n-1 --> n-1
+!     * comparisons between integers (n<2): 0 --> 1
+!     * sum of integers (i+1, i-1): n --> n-1
+!     * booleans assigments: 0 --> n
+    if (n<2) return
+    larger_than_last = (x>a(1))
+    do i=2,n
+      if (x>a(i).neqv.larger_than_last) then
+        larger_than_last = (x>a(i))
         nj=nj+1
-        if (nj<=m) j(nj)=i
+        if (nj<=m) j(nj)=i-1
       end if
     end do
   end subroutine locate_unord