Merge branch 'master' into arldatastream

.gitlab-ci.yml

@@ -40,6 +40,7 @@ mac_llvm_testing:
     - mkdir build
     - cd build
     - which cmake
+    - export radix_ENABLE_Fortran=OFF
     - cmake -DDEBUG_OUTPUT=1 -DBUILDNAME=$(uname -s)-LLVM-Debug-${CI_BUILD_REF_NAME} -DCMAKE_BUILD_TYPE=DEBUG -Dradix_ENABLE_TESTS=ON -Dradix_ENABLE_SECONDARY_TESTED_CODE=ON -Dradix_ENABLE_TESTS=ON -DTPL_ENABLE_VTK=ON -Dradix_ENABLE_radixplot=OFF ..
     - ctest -D ExperimentalStart -D ExperimentalBuild -D ExperimentalTest -DExperimentalMemCheck -D ExperimentalSubmit
 
@@ -131,6 +132,7 @@ windows_msvc_testing:
    - mkdir build
    - cd build
    - SET VTK_DIR=c:\vendors\cl\vtk\8.1.0\
+   - SET radix_ENABLE_Fortran=OFF
    - cmake -DBUILD_SHARED_LIBS=ON -DBUILDNAME=Windows-CL-18-Release-%CI_BUILD_REF_NAME% -DCMAKE_BUILD_TYPE=RELEASE -Dradix_ENABLE_SECONDARY_TESTED_CODE=ON -Dradix_ENABLE_TESTS=ON -DTPL_ENABLE_VTK=ON -Dradix_ENABLE_radixplot=OFF -Dradix_ENABLE_radixglls=OFF -G "NMake Makefiles" ..
    - ctest -D ExperimentalStart -D ExperimentalBuild -D ExperimentalTest -D ExperimentalSubmit
 

cmake/CallbackSetupExtraOptions.cmake

@@ -64,8 +64,10 @@ MACRO(TRIBITS_REPOSITORY_SETUP_EXTRA_OPTIONS)
   IF(NOT "$ENV{${PROJECT_NAME}_ENABLE_Fortran}" STREQUAL "" )
      SET(${PROJECT_NAME}_ENABLE_Fortran $ENV{${PROJECT_NAME}_ENABLE_Fortran})
   ELSE()
-     # ensure fortran compiler is set on
-     SET(${PROJECT_NAME}_ENABLE_Fortran OFF CACHE BOOL "" FORCE)
+     # allow fortran compiler to be set on
+     IF(NOT ${PROJECT_NAME}_ENABLE_Fortran)
+        SET(${PROJECT_NAME}_ENABLE_Fortran OFF CACHE BOOL "" FORCE)
+     ENDIF()
   ENDIF()
   # Set up radix cmake directory, used by default option scripts
   SET(radix_CMAKE_DIR "${radix_SOURCE_DIR}/cmake" CACHE PATH "")

radixalgorithm/CMakeLists.txt

@@ -11,6 +11,11 @@ marchingsquares.hh
 marchingsquares.i.hh
 )
 
+IF(${PROJECT_NAME}_ENABLE_Fortran)
+   SET(SOURCE ${SOURCE}
+       mergesort.f90)
+ENDIF()
+
 TRIBITS_ADD_LIBRARY(radixalgorithmlib
    SOURCES ${SOURCE}
    NOINSTALLHEADERS ${HEADERS}

radixalgorithm/marchingsquares.i.hh

 
-#include <queue>
 #include <set>
 #include <vector>
 
@@ -336,7 +335,6 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
 
   //
   // Initialize bit field
-  std::vector<size_t> label_counts(isovalues_size + 1, 0);
   size_t data_size = mData.size();
   for (size_t p_i = 0; p_i < data_size; ++p_i)
   {
@@ -357,12 +355,8 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
 
   typedef std::pair<float, float> Point;
   typedef std::vector<Point> Polygon;
-  struct PolyComparator
-  {
-    bool operator()(const Polygon& a, const Polygon& b)
-    {
-      return a.size() < b.size();
-    }
+  auto poly_comparator = [](const Polygon& a, const Polygon& b) {
+    return (a.size() < b.size());
   };
 
   for (size_t isoi = 0; isoi < isovalues_size; ++isoi)
@@ -378,7 +372,7 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
 
     // reset previous for each contour search
     prev = {0, 0};
-    std::priority_queue<Polygon, std::vector<Polygon>, PolyComparator> pqueue;
+    std::vector<Polygon> pvector;
     while (starting_point(label, start, prev))
     {
       prev = start;
@@ -426,11 +420,28 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
       //
       // Finish connecting the contour by making the last=first
       polygon.push_back(start);
-      pqueue.push(polygon);
+      if (pvector.size() == max_contour_polygons)
+      {
+        // compare against the last element
+        // if it is larger then replace it.
+        if (pvector[pvector.size() - 1].size() < polygon.size())
+        {
+          pvector[pvector.size() - 1] = polygon;
+          std::sort(pvector.begin(), pvector.end(), poly_comparator);
+        }
+        // we ignore as this polygon is not large enough to keep
+      }
+      else
+      {
+        pvector.push_back(polygon);
+        std::sort(pvector.begin(), pvector.end(), poly_comparator);
+      }
 
       radix_tagged_line("clearing label ("
                         << mBit[mColumns * start.first + start.second]
-                        << ") threshold (" << ordered_isovalues[isoi] << ")");
+                        << ") threshold (" << ordered_isovalues[isoi]
+                        << ") starting at [" << start.first << ","
+                        << start.second << "]");
 
       clear_isovalue_label(start.first, start.second,
                            mBit[mColumns * start.first + start.second],
@@ -438,30 +449,28 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
     }
     //
     // save the top N polygons
-    size_t count = std::min(max_contour_polygons, pqueue.size());
+    size_t count = pvector.size();
     contours[contouri].resize(count);
     for (size_t i = 0; i < count; ++i)
     {
-      auto& top      = pqueue.top();
+      auto& polygon  = pvector[i];
       auto& cpolygon = contours[contouri][i];
-      if (top.size() > max_polygon_points)
+      if (polygon.size() > max_polygon_points)
       {
         // std::cout << "Taking every other point given size(" << top.size()
         //          << ") exceeds " << max_polygon_points << std::endl;
         // take every other point
-        size_t polygon_size = top.size() / 2;
+        size_t polygon_size = polygon.size() / 2;
         for (size_t pi = 0; pi < polygon_size; ++pi)
         {
-          cpolygon.push_back(top[pi * 2]);
+          cpolygon.push_back(polygon[pi * 2]);
         }
       }
       else
       {
-        cpolygon.resize(top.size());
-        std::copy(top.begin(), top.end(), cpolygon.begin());
+        cpolygon.resize(polygon.size());
+        std::copy(polygon.begin(), polygon.end(), cpolygon.begin());
       }
-      // remove the top element
-      pqueue.pop();
     }
   }
   return contours;
@@ -472,8 +481,9 @@ void MarchingSquares<data_type>::clear_isovalue_label(
     int row, int col, short label,
     const std::vector<data_type>& ordered_isovalues)
 {
-  if (row < 0 || row == mColumns) return;  // out of bounds
-  if (col < 0 || col == mRows) return;     // out of bounds
+  // we can allow for row=mRows and col=mColumns
+  if (row < 0 || row > mRows) return;     // out of bounds
+  if (col < 0 || col > mColumns) return;  // out of bounds
 
   std::set<size_t> list;
   list.insert(mColumns * row + col);
@@ -488,25 +498,6 @@ void MarchingSquares<data_type>::clear_isovalue_label(
     row = c_i / mColumns;
     // upate the column
     col = c_i % mColumns;
-    // re-initilize bit field
-    data_type value = mData[c_i];
-    if (label > 1)  // not the last isovalue
-    {
-      // radix_tagged_line("Comparing " << value
-      //                               << " >= " << ordered_isovalues[isoi]);
-      if (value >= ordered_isovalues[isoi])
-      {
-        //
-        // save the category (aka there are only number of isovalue categories)
-        // specifically 1-N, 0 is background
-        mBit[c_i] = label - 1;
-        // radix_tagged_line("Found new isovalue (" << ordered_isovalues[isoi]
-        //                                         << ") for [" << row << ","
-        //                                         << col << "]");
-      }
-    }
-    // we didn't find a lower level contour, so zero it out
-    if (mBit[c_i] == label) mBit[c_i] = 0;
     // search neighbors
     for (int direction = 0; direction < 4; ++direction)
     {
@@ -521,10 +512,36 @@ void MarchingSquares<data_type>::clear_isovalue_label(
         // don't add it again
         if (list.find(nc_i) == list.end())
         {
+          radix_tagged_line("\t[" << nr << "," << nc << "]");
           list.insert(nc_i);
         }
       }
     }
+    // re-initilize bit field
+    if (label > 1)  // not the last isovalue
+    {
+      data_type value = mData[c_i];
+      radix_tagged_line("Comparing " << value
+                                     << " >= " << ordered_isovalues[isoi]);
+      if (value >= ordered_isovalues[isoi])
+      {
+        //
+        // save the category (aka there are only number of isovalue categories)
+        // specifically 1-N, 0 is background
+        mBit[c_i] = label - 1;
+        radix_tagged_line("Found new isovalue (" << ordered_isovalues[isoi]
+                                                 << ") for [" << row << ","
+                                                 << col << "]");
+      }
+    }
+    // we didn't find a lower level contour, so zero it out
+    radix_tagged_line("Comparing mBit[" << c_i << "](" << mBit[c_i] << ") to "
+                                        << label);
+    if (mBit[c_i] == label)
+    {
+      radix_tagged_line("\tzeroing [" << row << "," << col << "]");
+      mBit[c_i] = 0;
+    }
     list.erase(it);
   }
 }

radixalgorithm/mergesort.f90

+    subroutine merge_sort(n, a, indices)
+! This routine takes an input array 'a' of dimension (n)
+! and an input array 'indices' of dimension (n)
+! and conducts a merge sort on array 'a' while saving the original order
+! to array 'indices'
+
+! The first pass sorts each pair of consecutive elements [a(i) vs. a(i+1)].
+! The second pass sorts adjacent sequences of length 2 [a(i), a(i+1) vs. a(i+2), a(i+3)]; the sequences are already ordered, due to step 1.
+! The third pass sorts adjacent sequences of length 4 [a(i),...,a(i+3) vs. a(i+4),...,a(i+7)]; the sequences are already ordered, from step 2.
+! This continues for L1 passes, where 2**(L1-1) < n < 2**L1.  Note that n does not need to be a power of 2.
+
+    implicit none
+    !DEC$ ATTRIBUTES DLLEXPORT::merge_sort
+    integer :: i, imax, i0, j, jmax, k, L, L1, m, n
+    real, intent(inout) :: a(n)
+    integer, intent(out) :: indices(n)
+    real, allocatable :: b(:)
+    integer, allocatable :: ai(:)
+    allocate (b(n))
+    allocate (ai(n))
+    do i=1, n
+      ai(i) = i
+    end do
+    L1 = 1
+    m  = 1
+    do while (m < n)    ! Determine L1 so that 2**(L1-1) < n < 2**L1
+      m = m + m
+      L1 = L1 + 1
+    end do
+    L1 = L1 - 1
+    m = 1
+    do L=1, L1
+      k=1
+      do i0=1, n-m+1, m+m
+        i=i0
+        j=i+m
+        imax=j
+        jmax=j+m
+        if(imax > n) imax = n + 1
+        if(jmax > n) jmax = n + 1
+        do while(i < imax  .and.  j < jmax)
+          if(a(i) < a(j)) then
+            b(k) = a(i)
+            indices(k) = ai(i)
+            i = i + 1
+          else
+            b(k) = a(j)
+            indices(k) = ai(j)
+            j = j + 1
+          end if
+          k = k + 1
+        end do
+        do while(i < imax)
+          b(k) = a(i)
+          indices(k) = ai(i)
+          i = i + 1
+          k = k + 1
+        end do
+        do while(j < jmax)
+          b(k) = a(j)
+          indices(k) = ai(j)
+          j = j + 1
+          k = k + 1
+        end do
+      end do
+      m = m + m
+      a = b
+      ai = indices
+    end do
+    if(allocated(b)) deallocate(b)
+    if(allocated(ai)) deallocate(ai)
+    return
+    end subroutine merge_sort
+

radixalgorithm/tests/CMakeLists.txt

@@ -3,3 +3,9 @@ INCLUDE(GoogleTest)
 ADD_GOOGLE_TEST(tstOrdering.cc NP 1)
 ADD_GOOGLE_TEST(tstMarchingSquares.cc NP 1)
 ADD_GOOGLE_TEST(tstChaikins.cc NP 1)
+IF(${PROJECT_NAME}_ENABLE_Fortran)
+TRIBITS_ADD_EXECUTABLE_AND_TEST(tstMergeSortFortran
+    SOURCES tstMergeSort.f90
+    LINKER_LANGUAGE Fortran
+    )
+ENDIF()

radixalgorithm/tests/tstMarchingSquares.cc

@@ -276,7 +276,7 @@ TEST(MarchingSquares, NestedMultiIsovalues)
   {
     std::string line;
     std::getline(stream, line);
-    std::vector<std::string> fields = radix::split_string(",", line);
+    std::vector<std::string> fields = radix::split_string(",", line, true);
     for (size_t i = 0; i < fields.size(); ++i)
     {
       double value = std::atof(fields[i].c_str());
@@ -287,6 +287,7 @@ TEST(MarchingSquares, NestedMultiIsovalues)
   }
   size_t rows    = 256;
   size_t columns = rows;
+  radix_tagged_line("Grid size(" << grid.size() << ")");
   MarchingSquares<double> ms(grid, rows, columns);
 
   std::vector<std::vector<std::vector<std::pair<float, float>>>> contours;

radixalgorithm/tests/tstMergeSort.f90

+    Program run_merge
+    implicit none
+    interface
+      subroutine merge_sort(n, a, indices)
+      integer :: n
+      real, intent(inout) :: a(n)
+      integer, intent(out) :: indices(n)
+      end subroutine
+    end interface
+    integer :: n, i
+    character (len=10) :: d1, d2, d3
+    integer :: time(8)
+    real :: dtime1, dtime2
+    real, allocatable :: a(:)
+    real, allocatable :: blessed(:)
+    integer, allocatable :: indices(:)
+    call DATE_AND_TIME(d1, d2, d3, time)
+    dtime1 = time(6)*60.+time(7)+time(8)/1000.
+    n = 1000
+    allocate (a(n))
+    allocate (blessed(n))
+    allocate (indices(n))
+    ! avoid unitialised value warning
+    indices = 0
+    do i=1, time(7)
+      call RANDOM_NUMBER(dtime2)
+    end do
+    call RANDOM_NUMBER(a)
+    do i=1, n
+      a(i) = a(i) * n / 10.
+    end do
+    ! save this original version of a as blessed version
+    blessed = a
+!    print "(10f12.3)", (a(i,1), i=1, 10), (a(i,2), i=1, 10)
+!    a = a * N / 10.
+    call DATE_AND_TIME(d1, d2, d3, time)
+    dtime2 = time(6)*60.+time(7)+time(8)/1000.
+    dtime1 = dtime2 - dtime1
+    print *, "setup time (s) =", dtime1
+    print *, "First/Last 10 elements before sorting"
+    print "(10f12.3/10i12)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10i12)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
+    call merge_sort(n, a, indices)
+    call DATE_AND_TIME(d1, d2, d3, time)
+    print *, "First/Last 10 elements after sorting"
+    print "(10f12.3/10i12)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10i12)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
+    dtime2 = time(6)*60.+time(7)+time(8)/1000.-dtime2
+    print *, "sort time (s) =", dtime2
+    do i=1, (n-1)
+      ! check order of entire array
+      if(a(i) > a(i+1)) then
+         print *, "Element (",i,") is greater than element (",(i+1),")"
+         stop 99
+      endif
+    end do
+    ! check that the lookup works now
+
+    !print *, "unsorted:(10f3.2)", (blessed(i), i=1, n)
+    !print *, "Indices:(10i5)", (indices(i), i=1, n)
+    !print *, "sorted:(10f3.2)", (a(i), i=1, n)
+    do i=1, n
+      if(blessed(indices(i)) /= a(i)) then
+         print *, "Index (",i,") doesn't not match unsorted(",blessed(indices(i)),") with sorted(",a(i),")"
+         stop 100
+      endif
+    end do
+    stop
+    end program run_merge