From 43d2eff687a6f94a6a8538ff01ca4180c724a406 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Tue, 21 Dec 2021 19:54:12 +0100
Subject: [PATCH] Do not use a second, 0-based, date index

We are reverting commit e48b7f7e. I think I did not really appreciate
then that the number of days since 1950-1-1 is another date index. It
is clearer to use just one date index, which is in the dbf files. So,
in effect, we replace everywhere the old 0-based date index k by `k -
d_init`. Note that we no longer need argument `d_init` for procedure
`dispatch_snapshot` and we no longer need to write `d_init` to file
`node_id_param.json`. Note also that the value of every node index
changes by the amount `d_init * e_overestim`.
---
 Common/read_snapshot.f90                      | 13 +++++-----
 Inst_eddies/Tests/test_nearby_extr.f90        |  2 +-
 Overlap/Tests/Input/NSO_delta_in.txt          |  2 +-
 .../Input/successive_overlap_periodic_in.txt  |  2 +-
 Overlap/Tests/short_tests.json                | 24 +++++++++----------
 Overlap/Tests/test_get_dispatch_snap.f90      | 11 ++++-----
 Overlap/Tests/test_overlap.f90                |  7 +++---
 Overlap/Tests/test_read_snapshot.f90          |  2 +-
 Overlap/Tests/test_send_recv.f90              |  6 ++---
 Overlap/dispatch_snapshot.f90                 | 13 ++++------
 Overlap/eddy_graph.f90                        | 17 +++++++------
 Overlap/overlap.f90                           |  2 +-
 Overlap/report_graph.py                       | 22 +++++++----------
 13 files changed, 56 insertions(+), 67 deletions(-)

diff --git a/Common/read_snapshot.f90 b/Common/read_snapshot.f90
index 9069e4b1..3239e8b0 100644
--- a/Common/read_snapshot.f90
+++ b/Common/read_snapshot.f90
@@ -41,24 +41,23 @@ contains
 
     !---------------------------------------------------------------------
 
-    if (k == 0) then
+    if (k == d_init) then
        ishape_first = 0
     else
-       ! {k > 0}
-       ishape_first = ishape_last(k) + 1
+       ! {k > d_init}
+       ishape_first = ishape_last(k - d_init) + 1
     end if
 
-    s%number_extr = ishape_last(k + 1) - ishape_first + 1
+    s%number_extr = ishape_last(k - d_init + 1) - ishape_first + 1
     allocate(s%list(s%number_extr))
 
-    do ishape = ishape_first, ishape_last(k + 1)
+    do ishape = ishape_first, ishape_last(k - d_init + 1)
        call read_eddy(e, date, eddy_i, hshp, ishape)
 
        ! Check that all the eddies have the same date index:
-       if (date - d_init /= k) then
+       if (date /= k) then
           print *, "read_snapshot: bad date index"
           print *, "date = ", date
-          print *, "d_init = ", d_init
           print *, "k = ", k
           print *, "ishape = ", ishape
           print *, "eddy_i = ", eddy_i
diff --git a/Inst_eddies/Tests/test_nearby_extr.f90 b/Inst_eddies/Tests/test_nearby_extr.f90
index f2d55824..280415c3 100644
--- a/Inst_eddies/Tests/test_nearby_extr.f90
+++ b/Inst_eddies/Tests/test_nearby_extr.f90
@@ -50,7 +50,7 @@ program test_nearby_extr
 
   call shpc_open(hshp, trim(shpc_dir), rank = 0, pszaccess = "rb")
   call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0)
-  call read_snapshot(s, hshp, nlon, nlat, d_init, k = 0, &
+  call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, &
        corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
        copy = 0, ishape_last = [ishape_last])
   CALL shpc_close(hshp)
diff --git a/Overlap/Tests/Input/NSO_delta_in.txt b/Overlap/Tests/Input/NSO_delta_in.txt
index d2089bdd..3a6620d9 100644
--- a/Overlap/Tests/Input/NSO_delta_in.txt
+++ b/Overlap/Tests/Input/NSO_delta_in.txt
@@ -1,3 +1,3 @@
- &MAIN_NML MAX_DELTA       = 2 /
+ &MAIN_NML MAX_DELTA       = 2, k_test_1 = 20454, k_test_2 = 20455 /
 ,,,,,,1
 ,,,,,,1
diff --git a/Overlap/Tests/Input/successive_overlap_periodic_in.txt b/Overlap/Tests/Input/successive_overlap_periodic_in.txt
index d04199ea..8fa2295e 100644
--- a/Overlap/Tests/Input/successive_overlap_periodic_in.txt
+++ b/Overlap/Tests/Input/successive_overlap_periodic_in.txt
@@ -1,3 +1,3 @@
- &MAIN_NML DIST_LIM = 4,k_test_2 = 0 /
+ &MAIN_NML DIST_LIM = 4,k_test_1 = 20454, k_test_2 = 20454 /
 /
 /
diff --git a/Overlap/Tests/short_tests.json b/Overlap/Tests/short_tests.json
index 023d6ccb..b52e2e5f 100644
--- a/Overlap/Tests/short_tests.json
+++ b/Overlap/Tests/short_tests.json
@@ -47,7 +47,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML k_test_1 = 1 /\n/\n/\n",
+        "input": "&MAIN_NML k_test_1 = 20455, k_test_2 = 20455 /\n/\n/\n",
         "description": "Overlap of a snapshot with itself.",
         "title": "Successive_overlap",
         "command":
@@ -57,7 +57,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML k_test_1 = 1 /\n/\n/\n",
+        "input": "&MAIN_NML k_test_1 = 20455, k_test_2 = 20455 /\n/\n/\n",
         "description": "Same as Successive_overlap, but with larger domain: region 5 instead of region 4.",
         "title": "Successive_overlap_region_5_one_date",
         "command":
@@ -77,7 +77,7 @@
         ]
     },
     {
-	"input": "&MAIN_NML /\n/\n/\n",
+	"input": "&MAIN_NML k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description": "Overlap of different snapshots.",
         "title": "Successive_overlap_different_snapshots",
         "command":
@@ -87,7 +87,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML /\n/\n/\n",
+        "input": "&MAIN_NML k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description":
 	"Overlap of different snapshots. Same as Successive_overlap_different_snapshots, but with a larger region. The identifying numbers of the connected eddies are not the same for all edges.",
         "title": "Successive_overlap_region_2",
@@ -98,7 +98,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML /\n/\n/\n",
+        "input": "&MAIN_NML k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description": "Same as Successive_overlap_region_2, but with a larger region. We get some merging and splitting.",
         "title": "Successive_overlap_region_5",
         "command":
@@ -108,7 +108,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML /\n/\n/\n",
+        "input": "&MAIN_NML k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description": "Global grid, normal 0.25° resolution.",
         "title": "Successive_overlap_global",
         "command":
@@ -118,7 +118,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML max_delta = 3, k_test_2 = 0 /\n/\n/\n",
+        "input": "&MAIN_NML max_delta = 3, k_test_1 = 20454, k_test_2 = 20454 /\n/\n/\n",
         "description": "Overlap of a snapshot with itself. Same as Successive_overlap, except for max_delta.",
         "title": "Non_successive_overlap",
         "command":
@@ -128,7 +128,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML max_delta = 3 /\n/\n/\n",
+        "input": "&MAIN_NML max_delta = 3, k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description":
 	"Overlap of different snapshots. Same as Successive_overlap_different_snapshots, except for max_delta.",
         "title": "NSO_different_snapshots",
@@ -139,7 +139,7 @@
             ]
     },
     {
-        "input": "&MAIN_NML max_delta = 3 /\n/\n/\n",
+        "input": "&MAIN_NML max_delta = 3, k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description":
 	"Same as Successive_overlap_region_5, except for max_delta.",
         "title": "NSO_region_5",
@@ -150,7 +150,7 @@
             ]
     },
     {
-        "input": "&MAIN_NML max_delta = 3 /\n/\n/\n",
+        "input": "&MAIN_NML max_delta = 3, k_test_1 = 20454, k_test_2 = 20455 /\n/\n/\n",
         "description":
 	"Same as Successive_overlap_global, except for max_delta.",
         "title": "NSO_global",
@@ -191,7 +191,7 @@
     },
     {
         "input":
-	"&MAIN_NML CORNER = 5.125, -36.375, NLON = 29, NLAT = 17 /\n",
+	"&MAIN_NML CORNER = 5.125, -36.375, NLON = 29, NLAT = 17, k = 20454 /\n",
         "title": "Get_dispatch_snap",
         "stdout": "test_get_dispatch_snap_stdout.txt",
         "commands": [
@@ -209,7 +209,7 @@
         ]
     },
     {
-        "input": "&MAIN_NML CORNER = 5.125, -36.375, NLON = 29, NLAT = 17, k = 2 /\n",
+        "input": "&MAIN_NML CORNER = 5.125, -36.375, NLON = 29, NLAT = 17, k = 20456 /\n",
         "description": "Receive case for one of the processes in procedure get_snapshot.",
         "title": "Get_dispatch_snap_recv",
         "stdout": "test_get_dispatch_snap_stdout.txt",
diff --git a/Overlap/Tests/test_get_dispatch_snap.f90 b/Overlap/Tests/test_get_dispatch_snap.f90
index 656f32cc..3a0a1c3c 100644
--- a/Overlap/Tests/test_get_dispatch_snap.f90
+++ b/Overlap/Tests/test_get_dispatch_snap.f90
@@ -96,23 +96,22 @@ program test_get_dispatch_snap
   if (rank == 0) call dbf_read_attribute_03(d_init, hshp%extremum, &
        hshp%extr_date, ishape = 0)
   call ezmpi_bcast(d_init, root = 0)
-  k_begin = (rank * n_dates) / n_proc
+  k_begin = d_init + (rank * n_dates) / n_proc
   
   if (rank < n_proc - 1) then
-     k_end = ((rank + 1) * n_dates) / n_proc
+     k_end = d_init + ((rank + 1) * n_dates) / n_proc
   else
-     k_end = n_dates - 1
+     k_end = d_init + n_dates - 1
   end if
   
   call get_snapshot(s, nlon, nlat, ishape_last, corner * deg_to_rad, &
        step * deg_to_rad, copy, hshp, d_init, k, k_end, rank, n_proc, &
        max_delta = 1)
   CALL shpc_close(hshp)
-  call dispatch_snapshot(s, unit_isolated, rank, k_begin, max_delta = 1, &
-       d_init = d_init, k = k)
+  call dispatch_snapshot(s, unit_isolated, rank, k_begin, max_delta = 1, k = k)
 
   if (rank == 0) then
-     call write_snapshot(s, corner, step, nlon, nlat, copy, d = d_init + k)
+     call write_snapshot(s, corner, step, nlon, nlat, copy, d = k)
      close(unit_isolated)
   end if
 
diff --git a/Overlap/Tests/test_overlap.f90 b/Overlap/Tests/test_overlap.f90
index 4a6e2810..d48ff4f7 100644
--- a/Overlap/Tests/test_overlap.f90
+++ b/Overlap/Tests/test_overlap.f90
@@ -68,15 +68,14 @@ program test_overlap
   allocate(flow(max_delta + 1))
   call read_column(ishape_last, file = trim(shpc_dir) // "/ishape_last.txt")
   n_dates = size(ishape_last)
-  call assert(0 <= [k_test_1, k_test_2] .and. [k_test_1, k_test_2] < n_dates, &
-       "test_overlap k_test_1, k_test_2")
+  call assert(d_init <= [k_test_1, k_test_2] .and. [k_test_1, k_test_2] &
+       < d_init + n_dates, "test_overlap k_test_1, k_test_2")
   e_overestim = maxval([ishape_last(1) + 1, ediff1d(ishape_last)])
   call shpc_open(hshp, trim(shpc_dir), rank = 0, pszaccess = "rb")
   call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0)
   open(unit, file = "node_id_param.json", status = "replace", &
        action = "write")
-  write(unit, fmt = *) '{"e_overestim": ', e_overestim, ', "d_init": ', &
-       d_init, '}'
+  write(unit, fmt = *) '{"e_overestim": ', e_overestim, '}'
   close(unit)
   call read_snapshot(flow(1), hshp, nlon, nlat, d_init, k_test_1, corner, &
        step, copy, ishape_last)
diff --git a/Overlap/Tests/test_read_snapshot.f90 b/Overlap/Tests/test_read_snapshot.f90
index efc9ef12..0dc040a2 100644
--- a/Overlap/Tests/test_read_snapshot.f90
+++ b/Overlap/Tests/test_read_snapshot.f90
@@ -68,7 +68,7 @@ program test_read_snapshot
 
   call shpc_open(hshp, trim(shpc_dir), rank = 0, pszaccess = "rb")
   call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0)
-  call read_snapshot(s, hshp, nlon, nlat, d_init, k = 0, &
+  call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, &
        corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
        copy = copy, ishape_last = [ishape_last])
   CALL shpc_close(hshp)
diff --git a/Overlap/Tests/test_send_recv.f90 b/Overlap/Tests/test_send_recv.f90
index ec179173..658e63f2 100644
--- a/Overlap/Tests/test_send_recv.f90
+++ b/Overlap/Tests/test_send_recv.f90
@@ -102,19 +102,19 @@ program test_send_recv
   if (rank == 0) call dbf_read_attribute_03(d_init, hshp%extremum, &
        hshp%extr_date, ishape = 0)
   call ezmpi_bcast(d_init, root = 0)
-  if (rank == 1) call read_snapshot(s, hshp, nlon, nlat, d_init, k = 0, &
+  if (rank == 1) call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, &
        corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
        copy = copy, ishape_last = ishape_last)
   CALL shpc_close(hshp)
 
   if (rank == 1) then
-     call send_snapshot(s, dest = 0, tag = 0)
+     call send_snapshot(s, dest = 0, tag = d_init)
   else
      ! rank == 0
      call MPI_Comm_get_attr(MPI_Comm_world, MPI_TAG_UB, attribute_val, flag)
      call assert(flag, "test_send_recv MPI_Comm_get_attr MPI_TAG_UB")
      print *, "MPI_TAG_UB = ", attribute_val
-     call recv_snapshot(s, nlon, nlat, copy, source = 1, tag = 0)
+     call recv_snapshot(s, nlon, nlat, copy, source = 1, tag = d_init)
      call write_snapshot(s, corner_deg, step_deg, nlon, nlat, copy, d_init)
   end if
 
diff --git a/Overlap/dispatch_snapshot.f90 b/Overlap/dispatch_snapshot.f90
index f3662094..daeb1449 100644
--- a/Overlap/dispatch_snapshot.f90
+++ b/Overlap/dispatch_snapshot.f90
@@ -4,8 +4,7 @@ module dispatch_snapshot_m
 
 contains
 
-  subroutine dispatch_snapshot(s, unit_isolated, rank, k_begin, max_delta, &
-       d_init, k)
+  subroutine dispatch_snapshot(s, unit_isolated, rank, k_begin, max_delta, k)
 
     use derived_types, only: snapshot
     use send_snapshot_m, only: send_snapshot
@@ -15,21 +14,19 @@ contains
     integer, intent(in):: unit_isolated
     ! logical unit for file isolated_nodes_$rank.csv
 
-    integer, intent(in):: rank, k_begin, max_delta, d_init
-    integer, intent(in):: k ! date index (0-based)
+    integer, intent(in):: rank, k_begin, max_delta
+    integer, intent(in):: k ! date index
 
     ! Local:
-    integer i, d
+    integer i
 
     !------------------------------------------------------------------
 
     if (rank == 0 .or. k >= k_begin + max_delta) then
-       d = k + d_init
-       
        do i = 1, s%number_extr
           if (s%list(i)%valid .and. s%list(i)%delta_in == huge(0) &
                .and. s%list(i)%delta_out == huge(0)) &
-               write(unit_isolated, fmt = *) d, i
+               write(unit_isolated, fmt = *) k, i
        end do
     else
        call send_snapshot(s, tag = k, dest = rank - 1)
diff --git a/Overlap/eddy_graph.f90 b/Overlap/eddy_graph.f90
index 028b3751..e75db8ac 100644
--- a/Overlap/eddy_graph.f90
+++ b/Overlap/eddy_graph.f90
@@ -22,7 +22,7 @@ program eddy_graph
 
   integer rank, n_proc, copy, delta, j
   integer:: n_dates = huge(0) ! number of dates to read
-  integer k ! date index ((0-based)
+  integer k ! date index
   integer d_init, k_begin, k_end, k_end_main_loop
   character(len = :), allocatable:: shpc_dir
   integer:: max_delta = 1
@@ -124,8 +124,7 @@ program eddy_graph
   if (rank == 0) then
      open(unit, file = "node_id_param.json", status = "replace", &
           action = "write")
-     write(unit, fmt = *) '{"e_overestim": ', e_overestim, ', "d_init": ', &
-          d_init, '}'
+     write(unit, fmt = *) '{"e_overestim": ', e_overestim, '}'
      close(unit)
   end if
 
@@ -138,13 +137,13 @@ program eddy_graph
   call open_edge_file(rank)
   ! We do not write the title line. That will be handled by eddy_graph.sh.
 
-  k_begin = (rank * n_dates) / n_proc
+  k_begin = d_init + (rank * n_dates) / n_proc
 
   if (rank < n_proc - 1) then
-     k_end = ((rank + 1) * n_dates) / n_proc + max_delta - 1
+     k_end = d_init + ((rank + 1) * n_dates) / n_proc + max_delta - 1
      k_end_main_loop = k_end - max_delta + 1
   else
-     k_end = n_dates - 1
+     k_end = d_init + n_dates - 1
      k_end_main_loop = k_end
   end if
 
@@ -168,7 +167,7 @@ program eddy_graph
 
   do k = k_begin + max_delta + 1, k_end_main_loop
      call dispatch_snapshot(flow(1), unit_isolated, rank, k_begin, max_delta, &
-          d_init, k = k - max_delta - 1)
+          k = k - max_delta - 1)
      flow(:max_delta) = flow(2:)
      call get_snapshot(flow(max_delta + 1), nlon, nlat, ishape_last, corner, &
           step, copy, hshp, d_init, k, k_end, rank, n_proc, max_delta)
@@ -184,7 +183,7 @@ program eddy_graph
   do k = k_end_main_loop + 1, k_end
      ! {rank < n_proc - 1 and k >= k_end - max_delta + 2}
      call dispatch_snapshot(flow(1), unit_isolated, rank, k_begin, max_delta, &
-          d_init, k = k - max_delta - 1)
+          k = k - max_delta - 1)
      flow(:max_delta) = flow(2:)
      call get_snapshot(flow(max_delta + 1), nlon, nlat, ishape_last, corner, &
           step, copy, hshp, d_init, k, k_end, rank, n_proc, max_delta)
@@ -199,7 +198,7 @@ program eddy_graph
 
   do j = 1, max_delta + 1
      call dispatch_snapshot(flow(j), unit_isolated, rank, k_begin, max_delta, &
-          d_init, k = k_end - max_delta - 1 + j)
+          k = k_end - max_delta - 1 + j)
   end do
 
   ! 5. Back matter
diff --git a/Overlap/overlap.f90 b/Overlap/overlap.f90
index b0579f6d..af0658b7 100644
--- a/Overlap/overlap.f90
+++ b/Overlap/overlap.f90
@@ -35,7 +35,7 @@ contains
     integer, intent(in):: e_overestim
     ! over-estimation of the number of eddies at each date
     
-    integer, intent(in):: k ! date index (0-based)
+    integer, intent(in):: k ! date index
     integer, intent(in):: delta ! between 1 and max_delta
 
     integer, intent(in):: j
diff --git a/Overlap/report_graph.py b/Overlap/report_graph.py
index 4003d52a..7d7a1927 100755
--- a/Overlap/report_graph.py
+++ b/Overlap/report_graph.py
@@ -12,22 +12,20 @@ from networkx import nx_agraph
 import os
 import json
 
-def node_to_date_eddy(node_index, e_overestim, d_init):
+def node_to_date_eddy(node_index, e_overestim):
     """Convert from node identification in graph to eddy identification in
     shapefiles.
 
     """
     k = (node_index - 1) // e_overestim
-    date = d_init + k
-    return {"eddy_index": node_index - k * e_overestim, "date_index": k,
-            "date": date}
+    return {"eddy_index": node_index - k * e_overestim, "date_index": k}
 
-def date_eddy_to_node(date, eddy_index, e_overestim, d_init):
+def date_eddy_to_node(date, eddy_index, e_overestim):
     """Convert from eddy identification in shapefiles to node
     identification in graph.
 
     """
-    return (date - d_init) * e_overestim + eddy_index
+    return date * e_overestim + eddy_index
 
 def read_eddy_graph(edgelist, shpc_dir = None):
     if os.access(edgelist, os.R_OK):
@@ -37,7 +35,6 @@ def read_eddy_graph(edgelist, shpc_dir = None):
         node_id_param_file = path.join(dir_edgelist, "node_id_param.json")
         with open(node_id_param_file) as f: node_id_param = json.load(f)
         G.graph["e_overestim"] = node_id_param["e_overestim"]
-        G.graph["d_init"] = node_id_param["d_init"]
 
         if shpc_dir is not None:
             # Read and set attributes of eddies:
@@ -54,7 +51,7 @@ def set_attribute(G, extr_file):
         for shape_rec in s_read:
             n = date_eddy_to_node(shape_rec.record.date,
                                   shape_rec.record.eddy_index,
-                                  G.graph["e_overestim"], G.graph["d_init"])
+                                  G.graph["e_overestim"])
             
             if shape_rec.record.valid == 1 and n in G:
                 G.add_node(n, coordinates = tuple(shape_rec.shape.points[0]),
@@ -65,13 +62,12 @@ def to_eddy_agraph(G):
     my_subgraphs = {}
 
     for n in G:
-        date_eddy = node_to_date_eddy(n, G.graph["e_overestim"],
-                                      G.graph["d_init"])
+        date_eddy = node_to_date_eddy(n, G.graph["e_overestim"])
 
-        if date_eddy["date"] not in my_subgraphs:
-            my_subgraphs[date_eddy["date"]] = []
+        if date_eddy["date_index"] not in my_subgraphs:
+            my_subgraphs[date_eddy["date_index"]] = []
 
-        my_subgraphs[date_eddy["date"]].append(n)
+        my_subgraphs[date_eddy["date_index"]].append(n)
 
     A = nx_agraph.to_agraph(G)
 
-- 
GitLab