[BugFixes] Fix bugs in python interface of SampEn2D.

Author: phreer

include/sample_entropy_calculator2d.h

@@ -206,13 +206,14 @@ class SampleEntropyCalculator2DSampling : public SampleEntropyCalculator2D<T> {
                                     unsigned dilation_factor,
                                     unsigned sample_size, unsigned sample_num,
                                     double real_entropy, double real_a_norm,
-                                    double real_b_norm, OutputLevel output_level)
+                                    double real_b_norm, bool random_,
+                                    OutputLevel output_level)
       : SampleEntropyCalculator2D<T>(first, last, m, r, width, height,
                                      moving_step_size, dilation_factor,
                                      output_level),
         _sample_size(sample_size), _sample_num(sample_num),
         _real_entropy(real_entropy), _real_a_norm(real_a_norm),
-        _real_b_norm(real_b_norm) {
+        _real_b_norm(real_b_norm), _random(random_) {
     if (_num_templates < _sample_size) {
       MSG_ERROR(-1, "Sample size (N0) exceeds the number of templates.\n");
     }
@@ -224,13 +225,15 @@ class SampleEntropyCalculator2DSampling : public SampleEntropyCalculator2D<T> {
                                     unsigned dilation_factor,
                                     unsigned sample_size, unsigned sample_num,
                                     double real_entropy, double real_a_norm,
-                                    double real_b_norm, OutputLevel output_level)
+                                    double real_b_norm, bool random_,
+                                    OutputLevel output_level)
       : SampleEntropyCalculator2DSampling(data.cbegin(), data.cend(), m, r,
                                           width, height, moving_step_size,
                                           dilation_factor,
                                           sample_size, sample_num,
                                           real_entropy, real_a_norm,
-                                          real_b_norm, output_level) {}
+                                          real_b_norm, random_,
+                                          output_level) {}
   vector<long long> get_a_vec() {
     if (!_computed)
       ComputeSampleEntropy();
@@ -294,6 +297,7 @@ class SampleEntropyCalculator2DSampling : public SampleEntropyCalculator2D<T> {
   const double _real_entropy;
   const double _real_a_norm;
   const double _real_b_norm;
+  const bool _random;
   vector<long long> _a_vec, _b_vec;
 };
 
@@ -322,6 +326,7 @@ class SampleEntropyCalculator2DSampling : public SampleEntropyCalculator2D<T> {
   using SampleEntropyCalculator2DSampling<T>::_real_entropy;\
   using SampleEntropyCalculator2DSampling<T>::_real_a_norm;\
   using SampleEntropyCalculator2DSampling<T>::_real_b_norm;\
+  using SampleEntropyCalculator2DSampling<T>::_random;\
   using SampleEntropyCalculator2DSampling<T>::_a_vec;\
   using SampleEntropyCalculator2DSampling<T>::_b_vec;\
   using SampleEntropyCalculator2DSampling<T>::_IsMatchedK;\
@@ -389,7 +394,7 @@ template <typename T>
 void SampleEntropyCalculator2DSamplingDirect<T>::_ComputeSampleEntropy() {
   unsigned num_templates = _num_steps_x * _num_steps_y;
   RandomIndicesSamplerWR sampler(num_templates, _sample_size, _sample_num,
-                                 RandomType::SWR_UNIFORM);
+                                 RandomType::SWR_UNIFORM, _random);
   const auto random_indices_array = sampler.GetSampleArrays();
   _a_vec.clear();
   _b_vec.clear();

python/sample_entropy_calculator2d.pxd

@@ -45,6 +45,7 @@ cdef extern from "sample_entropy_calculator2d.h" namespace "sampen":
                                             double real_entropy,
                                             double real_a_norm,
                                             double real_b_norm,
+                                            bool random_,
                                             OutputLevel output_level) except +
     double get_computation_time()
     double get_entropy()

python/sample_entropy_calculator2d.pyx

@@ -62,10 +62,11 @@ cdef class SampEn2DSamplingD:
                 unsigned dilation_factor,
                 unsigned sample_size, unsigned sample_num,
                 double real_entropy, double real_a_norm, double real_b_norm, 
-                OutputLevel level):
+                bool random_, OutputLevel level):
     self.c_ = new SampleEntropyCalculator2DSamplingDirect[int](
         data, r, m, width, height, moving_step_size, dilation_factor,
-        sample_size, sample_num, real_entropy, real_a_norm, real_b_norm, level)
+        sample_size, sample_num, real_entropy, real_a_norm, real_b_norm,
+        random_, level)
 
   def __dealloc__(self):
     if self.c_ != NULL:

python/test_sampen2d.py

@@ -1,3 +1,5 @@
+import os
+import sys
 import sampen2d
 import numpy as np
 from PIL import Image, ImageOps
@@ -8,7 +10,8 @@
 n1 = 20
 
 if __name__ == '__main__':
-  image_path = '../data2d/Original Brodatz/D1.gif'
+  input_dir = os.path.join(os.environ['HOME'], 'workspace/entropy/2d')
+  image_path = os.path.join(input_dir, 'Original_Brodatz/D1.gif')
   image = Image.open(image_path)
   image = ImageOps.grayscale(image)
   image = np.array(image)
@@ -24,6 +27,11 @@
 
   n0 = n * sample_rate
   s = sampen2d.SampEn2DSamplingD(image.flatten(), r_scaled, m, w, h, 1, 1,
-                                 n0, n1, r_e, r_a, r_b, 0)
+                                 n0, n1, r_e, r_a, r_b, True, 0)
+  print('time: ', s.time())
+  print('entropy: ', s.entropy())
+
+  s = sampen2d.SampEn2DSamplingD(image.flatten(), r_scaled, m, w, h, 1, 1,
+                                 n0, n1, r_e, r_a, r_b, True, 0)
   print('time: ', s.time())
   print('entropy: ', s.entropy())

script/python_experiments/exp_utils.py

@@ -84,31 +84,35 @@ def estimate_sampen2d_and_save_statistics(sig, r: float, m: int,
       r_e = s.entropy()
       r_a = s.a_norm()
       r_b = s.b_norm()
+      t = s.time()
     else:
       r_e, r_a, r_b = groundtruth
+      t = 0
     print(f'output_filename [{output_filename}]', file=f)
     print(f'n0 [{n0}] n1 [{n1}] m [{m}] r [{r}] n [{n}]', file=f)
-    print(f'[Direct] sampen2d [{r_e:.6f}] a [{r_a:.6f}] b [{r_b:.6f}]', file=f)
+    print(f'[Direct] sampen2d [{r_e:.6f}] a [{r_a:.6f}] b [{r_b:.6f}] t [{t:.4e}]', file=f)
     list_err_e = np.zeros(n_computations)
     list_err_a = np.zeros(n_computations)
     list_err_b = np.zeros(n_computations)
+    list_time = np.zeros(n_computations)
     for i in range(n_computations):
       s = sampen2d.SampEn2DSamplingD(sig, r_scaled, m, width, height, 1, 1,
-                                     n0, n1, r_e, r_a, r_b, 0)
+                                     n0, n1, r_e, r_a, r_b, True, 0)
       err_e = s.entropy() - r_e
       err_a = s.a_norm() - r_a
       err_b = s.b_norm() - r_b
       list_err_e[i] = err_e
       list_err_a[i] = err_a
       list_err_b[i] = err_b
-      print(f'[SWR] sampen2d [{s.entropy():.6f}] a [{s.a_norm():.6f}] b [{s.b_norm():.6f}]', file=f)
+      list_time[i] = s.time()
+      print(f'[SWR] sampen2d [{s.entropy():.6f}] a [{s.a_norm():.6f}] b [{s.b_norm():.6f}] t [{s.time():.4e}]', file=f)
       print(f'[SWR] err_sampen2d [{err_e :.4e}] err_a [{err_a:.4e}] err_b [{err_b:.4e}]', file=f)
     m_e, std_e, m_abs_err_e, r_m_e, r_std_e, r_m_abs_err_e = get_statistics(list_err_e, r_e)
     m_a, std_a, m_abs_err_a, r_m_a, r_std_a, r_m_abs_err_a = get_statistics(list_err_a, r_a)
     m_b, std_b, m_abs_err_b, r_m_b, r_std_b, r_m_abs_err_b = get_statistics(list_err_b, r_b)
-    print(f'[SWR] m_e [{m_e:.4e}] std_e [{std_e:.4e}] m_abs_err_e [{m_abs_err_e:.4e}]', file=f)    
+    print(f'[SWR] m_e [{m_e:.4e}] std_e [{std_e:.4e}] m_abs_err_e [{m_abs_err_e:.4e}] m_t [{np.mean(list_time):.4e}]', file=f)    
     print(f'[SWR] r_m_e [{r_m_e:.4e}] r_std_e [{r_std_e:.4e}] r_m_abs_err_e [{r_m_abs_err_e:.4e}]', file=f)    
     print(f'[SWR] m_a [{m_a:.4e}] std_a [{std_a:.4e}] m_abs_err_a [{m_abs_err_a:.4e}]', file=f)    
     print(f'[SWR] r_m_a [{r_m_a:.4e}] r_std_a [{r_std_a:.4e}] r_m_abs_err_a [{r_m_abs_err_a:.4e}]', file=f)    
     print(f'[SWR] m_b [{m_b:.4e}] std_b [{std_b:.4e}] m_abs_err_b [{m_abs_err_b:.4e}]', file=f)    
-    print(f'[SWR] r_m_b [{r_m_b:.4e}] r_std_b [{r_std_b:.4e}] r_m_abs_err_b [{r_m_abs_err_b:.4e}]', file=f)    
+    print(f'[SWR] r_m_b [{r_m_b:.4e}] r_std_b [{r_std_b:.4e}] r_m_abs_err_b [{r_m_abs_err_b:.4e}]', file=f)    

script/python_experiments/run_ex_2d.py

@@ -21,13 +21,13 @@
 script_dir = os.path.dirname(sys.argv[0])
 
 def estimate_sampen2d_and_save_task(db_name, record_name, r, m):
-  output_dir = pjoin(script_dir, '../../result/2d_err', db_name)
+  output_dir = pjoin(script_dir, '../../result/2d_err/m%d_r%.2f' % (m, r), db_name)
   output_filename = os.path.splitext(record_name)[0] + '.txt'
   input_path = pjoin(input_dir, db_name, record_name)
   image = Image.open(input_path)
   image = np.array(ImageOps.grayscale(image))
   height, width = image.shape
-  print('start estimating...')
+  print('start %s %s %.2f %d...' % (db_name, record_name, r, m))
   estimate_sampen2d_and_save_statistics(image.flatten(), r, m,
                                         width=width, height=height,
                                         n0=n0, n1=n1, outputdir=output_dir,

src/sampen2d.cpp

@@ -27,8 +27,8 @@ struct Argument {
   unsigned height;
   unsigned moving_step_size;
   unsigned dilation_factor;
-  bool sampling1;
-  bool sampling2;
+  bool sampling;
+  bool random_;
   bool multiscale;
   unsigned multiscale_depth;
   double multiscale_factor;
@@ -72,13 +72,13 @@ char usage[] =
     "    --multiscale-factor <FACTOR>\n"
     "        The scaling factor of resizing when computing multiscale SampEn2D.\n\n"
     "Options:\n"
-    "    --sampling1\n"
-    "        Use sampling method (approach 1).\n"
-    "    --sampling2\n"
-    "        Use sampling method (approach 2).\n"
+    "    --sampling\n"
+    "        Use sampling method to estimate SampEn2D.\n"
     "    --multiscale\n"
     "        Compute multiscale SampEn2D.\n"
-    "    --help\n"
+    "    -r | --random\n"
+    "        Set seed randomly.\n"
+    "    -h | --help\n"
     "        Display this message.\n";
 // clang-format on
 
@@ -99,8 +99,7 @@ void ParseArgument(int argc, char *argv[]) {
   arg.image_filename = parser.getArg("--input");
   arg.output_level =
       static_cast<sampen::OutputLevel>(parser.getArgLong("--output-level", 0));
-  arg.sampling1 = parser.isOption("--sampling1");
-  arg.sampling2 = parser.isOption("--sampling2");
+  arg.sampling = parser.isOption("--sampling");
   arg.multiscale = parser.isOption("--multiscale");
   arg.multiscale_depth = parser.getArgLong("--multiscale-depth", 1);
   arg.multiscale_factor = parser.getArgDouble("--multiscale-factor", 0.5);
@@ -113,9 +112,10 @@ void ParseArgument(int argc, char *argv[]) {
     std::cerr << "Usage: " << argv[0] << usage;
     exit(-1);
   }
-  if (arg.sampling1 || arg.sampling2) {
+  if (arg.sampling) {
     arg.sample_size = parser.getArgLong("--sample-size", 1024);
     arg.sample_num = parser.getArgLong("--sample-num", 20);
+    arg.random_ = parser.isOption("--random") || parser.isOption("-r");
   }
 }
 
@@ -149,7 +149,8 @@ std::vector<double> ComputeMultiscaleSampEn2D(Magick::Image image, double r,
       calculator =
           std::make_shared<SampleEntropyCalculator2DSamplingDirect<int>>(
               data.cbegin(), data.cend(), r, m, width, height, 1, 1,
-              arg.sample_size, arg.sample_num, 0., 0., 0., arg.output_level);
+              arg.sample_size, arg.sample_num, 0., 0., 0., arg.random_,
+              arg.output_level);
     } else {
       calculator = std::make_shared<SampleEntropyCalculator2DDirect<int>>(
           data.cbegin(), data.cend(), r, m, width, height, 1, 1,
@@ -224,11 +225,12 @@ int main(int argc, char *argv[]) {
   double a_norm = sec2dd.get_a_norm();
   double b_norm = sec2dd.get_b_norm();
 
-  if (arg.sampling1) {
+  if (arg.sampling) {
     sampen::SampleEntropyCalculator2DSamplingDirect<int> sec2dds(
         data.cbegin(), data.cend(), arg.r, arg.m, arg.width, arg.height,
         arg.moving_step_size, arg.dilation_factor, arg.sample_size,
-        arg.sample_num, sampen2d, a_norm, b_norm, arg.output_level);
+        arg.sample_num, sampen2d, a_norm, b_norm, arg.random_,
+        arg.output_level);
     std::cout << sec2dds.get_result_str();
   }
   if (arg.multiscale) {

src/utils.cpp

@@ -29,6 +29,7 @@ void ReportVmPeak() {
           MSG_INFO("%s\n", buffer);
         }
       }
+      ifs.close();
     }
 #endif
 }