Fast smooth emd review

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(NetEmdSmoothV2)
 export(adaptive_breaks)
 export(area_between_dhist_ecmfs)
 export(as_smoothed_dhist)
@@ -48,6 +49,7 @@ export(min_emd)
 export(min_emd_exhaustive)
 export(min_emd_optimise)
 export(min_emd_optimise_fast)
+export(netEMDSpeedTestSmooth)
 export(netdis)
 export(netdis.plot)
 export(netdis_centred_graphlet_counts)

R/RcppExports.R

@@ -25,3 +25,13 @@ emd_fast_no_smoothing <- function(locations1, values1, locations2, values2) {
     .Call(`_netdist_emd_fast_no_smoothing`, locations1, values1, locations2, values2)
 }
 
+#' @title
+#' Compute EMD
+NULL
+
+#'
+#' @export
+NetEmdSmoothV2 <- function(loc1, val1, binWidth1, loc2, val2, binWidth2) {
+    .Call(`_netdist_NetEmdSmoothV2`, loc1, val1, binWidth1, loc2, val2, binWidth2)
+}
+

R/emd.R

@@ -81,10 +81,49 @@ min_emd_optimise_fast <- function(dhist1, dhist2) {
     min_offset <- soln$minimum
     return(list(min_emd = min_emd, min_offset = min_offset))
   }
-  else {
-    # Fall back on other version if either dhist is smoothed
-    return(min_emd_optimise(dhist1, dhist2))
+  else #if ((dhist1$smoothing_window_width==1) && (dhist2$smoothing_window_width==1))
+  {
+    val1 <- cumsum(dhist1$masses)
+    val2 <- cumsum(dhist2$masses)
+    val1 <- val1/val1[length(val1)]
+    val2 <- val2/val2[length(val2)]
+    loc1=dhist1$locations
+    loc2=dhist2$locations
+    binWidth1=dhist1$smoothing_window_width
+    binWidth2=dhist2$smoothing_window_width
+    count=0
+    # Offset the histograms to make the alignments work
+    loc1Mod <- loc1 - binWidth1/2 
+    loc2Mod <- loc2 - binWidth2/2 
+    # Determine minimum and maximum offset of range in which histograms overlap
+    # (based on sliding histogram 1)
+    min_offset <- min(loc2Mod) - max(loc1Mod) - max(binWidth1,binWidth2)
+    max_offset <- max(loc2Mod) - min(loc1Mod) + max(binWidth1,binWidth2)
+    # Set lower and upper range for optimise algorithm to be somewhat wider than
+    # range defined by the minimum and maximum offset. This guards against a
+    # couple of issues that arise if the optimise range is exactly min_offset 
+    # to max_offset
+    # 1) If lower and upper are equal, the optimise method will throw and error
+    # 2) It seems that optimise is not guaranteed to explore its lower and upper
+    #    bounds, even in the case where one of them is the offset with minimum
+    #    EMD
+    buffer <- 0.1
+    min_offset <- min_offset - buffer
+    max_offset <- max_offset + buffer
+    # Define a single parameter function to minimise emd as a function of offset
+    emd_offset <- function(offset) {
+      temp1<- NetEmdSmoothV2(loc1Mod+offset,val1,binWidth1,loc2Mod,val2,binWidth2)
+      temp1
+    }
+    # Get solution from optimiser
+    soln <- stats::optimise(emd_offset, lower = min_offset, upper = max_offset, 
+                            tol = .Machine$double.eps*1000)
+    # Return mnimum EMD and associated offset
+    min_emd <- soln$objective
+    min_offset <- soln$minimum
+    return(list(min_emd = min_emd, min_offset = min_offset))
   }
+
 }
 
 

R/net_emd_speed_benchmark.R

@@ -32,3 +32,40 @@ netEMDSpeedTest <- function() {
   }
   list(gddBuildTime = gddBuildTime, netEMDtime = netEMDtime)
 }
+
+#' @export
+netEMDSpeedTestSmooth <- function()
+{
+  ##load the data
+  source_dir <- system.file(file.path("extdata", "random"), package = "netdist")
+  print(source_dir)
+  edge_format = "ncol"
+  file_pattern = ""
+  #    source_dir <- system.file(file.path("extdata", "VRPINS"), package = "netdist")
+  #    edge_format = "ncol"
+  #    file_pattern = ".txt"
+  graphs <- read_simple_graphs(source_dir = source_dir, format = edge_format, pattern = file_pattern)
+  n1=names(graphs)
+  lab1=c()
+  gddBuildTime=c()
+  netEMDtime=c()
+  for (i in 1:length(graphs))
+  {
+    for (j in 1:(i))
+    {
+      g1=graphs[[i]]
+      g2=graphs[[j]]
+      lab1=append(lab1,paste(n1[i],n1[j],sep=','))
+      print(paste(n1[i],n1[j],sep=','))
+      fulltimeStart=Sys.time()
+      gdd1=gdd(g1)
+      gdd2=gdd(g2)
+      netEMDStart=Sys.time()
+      net_emd(gdd1,gdd2,smoothing_window_width = 1)
+      endTime=Sys.time()
+      gddBuildTime=append(gddBuildTime,as.double(netEMDStart-fulltimeStart))
+      netEMDtime=append(netEMDtime,as.double(endTime-netEMDStart))
+    }
+  }
+  list(gddBuildTime = gddBuildTime, netEMDtime = netEMDtime)
+}

README.md

@@ -1,5 +1,6 @@
 # Network Comparison
-An R package implementing the Netdis and NetEMD alignment-free network comparison measures.
+An R package implementing the Netdis and NetEMD alignment-free network comparison measures.   
+
 
 ### :warning: BETA: Package under construction (pre-release) :warning:
 Until this package hits release 1.0 anything can change with no notice.

src/Makevars

@@ -1,2 +1,2 @@
 CXX_STD = CXX11
-PKG_CPPFLAGS += -fno-fast-math -msse2 -mfpmath=sse -mstackrealign
+PKG_CPPFLAGS += -fno-fast-math -msse2  -mstackrealign

src/RcppExports.cpp

@@ -5,6 +5,11 @@
 
 using namespace Rcpp;
 
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
 // counts_from_observations
 NumericMatrix counts_from_observations(NumericMatrix features);
 RcppExport SEXP _netdist_counts_from_observations(SEXP featuresSEXP) {
@@ -30,12 +35,29 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// NetEmdSmoothV2
+double NetEmdSmoothV2(NumericVector loc1, NumericVector val1, double binWidth1, NumericVector loc2, NumericVector val2, double binWidth2);
+RcppExport SEXP _netdist_NetEmdSmoothV2(SEXP loc1SEXP, SEXP val1SEXP, SEXP binWidth1SEXP, SEXP loc2SEXP, SEXP val2SEXP, SEXP binWidth2SEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type loc1(loc1SEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type val1(val1SEXP);
+    Rcpp::traits::input_parameter< double >::type binWidth1(binWidth1SEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type loc2(loc2SEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type val2(val2SEXP);
+    Rcpp::traits::input_parameter< double >::type binWidth2(binWidth2SEXP);
+    rcpp_result_gen = Rcpp::wrap(NetEmdSmoothV2(loc1, val1, binWidth1, loc2, val2, binWidth2));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 RcppExport SEXP run_testthat_tests();
 
 static const R_CallMethodDef CallEntries[] = {
     {"_netdist_counts_from_observations", (DL_FUNC) &_netdist_counts_from_observations, 1},
     {"_netdist_emd_fast_no_smoothing", (DL_FUNC) &_netdist_emd_fast_no_smoothing, 4},
+    {"_netdist_NetEmdSmoothV2", (DL_FUNC) &_netdist_NetEmdSmoothV2, 6},
     {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0},
     {NULL, NULL, 0}
 };

src/fastSmoothV2.cpp

@@ -0,0 +1,197 @@
+// Enable C++11
+// [[Rcpp::plugins(cpp11)]]
+#include <Rcpp.h>
+#include <unordered_set>
+#include <algorithm>
+#include <numeric>
+#include <math.h>
+#include "emd_fast_no_smoothing.h" // add_element_kahan()
+#include "fastSmoothV2.h"
+
+using namespace Rcpp;
+
+double bowtie_area(double length, double val1_start, double val1_end,
+                   double val2_start, double val2_end)
+{
+  double midPoint = (val1_start - val2_start) /
+    ((val2_end - val2_start) - (val1_end - val1_start));
+
+  const double midValue = val1_start + midPoint * (val1_end - val1_start);
+
+  midPoint = midPoint * length;
+
+  double topTriangle = 0.5 * midPoint * (midValue - val1_start);
+  double topRectangle = midPoint * (val1_start - val2_start);
+  double bottomTriangle = 0.5 * midPoint * (midValue - val2_start);
+
+  double res = topTriangle + topRectangle - bottomTriangle;
+
+  topTriangle = 0.5 * (length - midPoint) * (val2_end - midValue);
+  topRectangle = 0; // midPoint*(val1_start-val2_start);
+  bottomTriangle = 0.5 * (length - midPoint) * (val1_end - midValue);
+
+  res += topTriangle + topRectangle - bottomTriangle;
+  return res;
+}
+
+// Compute the unsigned area between two line segments
+// assumes that val1_end > val1_start and val2_end > val2_start
+double get_segment(double start, double end, double val1_start,
+                   double val1_end, double val2_start, double val2_end)
+{
+  const double length = end - start;
+
+  double topTriangle;
+  double topRectangle;
+  double bottomTriangle;
+  double midPoint;
+  double midValue;
+  double res = 0;
+
+  bool both_differences_positive = val1_start > val2_start && val1_end >= val2_end;
+  bool both_differences_negative = val1_start <= val2_start && val1_end <= val2_end;
+
+  if (both_differences_positive || both_differences_negative)
+  {
+    // They are in the same order: no bowtie
+    // triangle of seg1
+    topTriangle = 0.5 * length * (val1_end - val1_start);
+    // rectangle between seg1 and seg2
+    topRectangle = length * (val1_start - val2_start);
+    // triangle of seg2 (to be removed)
+    bottomTriangle = 0.5 * length * (val2_end - val2_start);
+
+    const double sign = both_differences_positive?1.0:-1.0;
+    return sign * (topTriangle + topRectangle - bottomTriangle);
+  }
+  else if (val1_start > val2_start) { // bowtie, first case
+    return bowtie_area(length, val1_start, val1_end, val2_start, val2_end);
+  }
+  else { // bowtie, second case
+    return bowtie_area(length, val2_start, val2_end, val1_start, val1_end);
+  }
+}
+
+// cut down and compute segment
+double get_segment_constrained(double seg1L1, double seg1L2,
+                               double seg2L1, double seg2L2,
+                               double seg1V1, double seg1V2,
+                               double seg2V1, double seg2V2)
+{
+  //We have a valid range
+  double valStart1, valEnd1, valStart2, valEnd2;
+  double start,end;
+  double result;
+  start = std::max(seg1L1,seg2L1);
+  end   = std::min(seg1L2,seg2L2);
+  if (start < end) {
+    valStart1 = seg1V1 + (seg1V2 - seg1V1)*(start - seg1L1)/(seg1L2 - seg1L1);
+    valEnd1   = seg1V1 + (seg1V2 - seg1V1)*(end   - seg1L1)/(seg1L2 - seg1L1);
+    valStart2 = seg2V1 + (seg2V2 - seg2V1)*(start - seg2L1)/(seg2L2 - seg2L1);
+    valEnd2   = seg2V1 + (seg2V2 - seg2V1)*(end   - seg2L1)/(seg2L2 - seg2L1);
+    result = get_segment(start, end, valStart1, valEnd1, valStart2, valEnd2);
+    return result;
+  }
+  else {
+    return 0;
+  }
+}
+
+double get_double_segment_constrained(
+  double seg1Loc1, double seg1Loc2, double seg1Loc3,
+  double seg1Val1, double seg1Val2,
+  double seg2Loc1, double seg2Loc2, double seg2Loc3,
+  double seg2Val1, double seg2Val2)
+{
+  double res = 0;
+
+  // compare the linear section with the linear section
+  res += get_segment_constrained(seg1Loc1, seg1Loc2, seg2Loc1, seg2Loc2,
+                                 seg1Val1, seg1Val2, seg2Val1, seg2Val2);
+
+  // compare the linear section with the flat section
+  // This could be easily special cased (saving ~1 if statements ).
+  res += get_segment_constrained(seg1Loc1, seg1Loc2, seg2Loc2, seg2Loc3,
+                                 seg1Val1, seg1Val2, seg2Val2, seg2Val2);
+
+  // compare the flat section with the linear section
+  // This could be easily special cased (saving ~1 if statements ).
+  res += get_segment_constrained(seg1Loc2, seg1Loc3, seg2Loc1, seg2Loc2,
+                                 seg1Val2, seg1Val2, seg2Val1, seg2Val2);
+
+  // compare the flat section with the flat section
+  // This could be easily special cased (saving ~2 if statements ).
+  res += get_segment_constrained(seg1Loc2, seg1Loc3, seg2Loc2, seg2Loc3,
+                                 seg1Val2, seg1Val2, seg2Val2, seg2Val2);
+
+  return res;
+}
+
+//' @title
+//' Compute EMD
+////'
+////' @param loc1 numeric vector.
+////' @param val1 numeric vector.
+////' @param loc2 numeric vector.
+////' @param val2 numeric vector.
+//'
+//' @export
+// [[Rcpp::export]]
+double NetEmdSmoothV2(NumericVector loc1, NumericVector val1, double binWidth1,
+                      NumericVector loc2, NumericVector val2, double binWidth2)
+{
+  OverlappingSegments segs(loc1, loc2, binWidth1, binWidth2);
+
+  double res = 0.0;
+  double compensation = 0.0;
+
+  for (OverlappingSegments::iterator it = segs.begin(), endIt = segs.end();
+       it != endIt; ++it)
+  {
+    // The OverlappingSegments iterator returns pairs of the left
+    // indices of overlapping endpoints: it->first is the index into
+    // loc1, it->second is the index into loc2.  When the smallest
+    // element in one of these vectors is larger than the current
+    // element of the other, an 'index' of -1 is returned.
+
+    // Hist 1
+    // Start of the gradient section in Seg1
+    double curSeg1Loc1 = segs.loc1_left(it->first);
+
+    // End of the gradient section in Seg1
+    double curSeg1Loc2 = segs.loc1_mid(it->first);
+
+    // End of the flat section in Seg1
+    double curSeg1Loc3 = segs.loc1_right(it->first);
+
+    // Start and end values in Seg1: val1 gives the values at *right*
+    // endpoints of the segments.  A value of 0.0 is used before the
+    // first segment.
+    double curSeg1Val1 = (it->first > 0) ? val1[it->first - 1] : 0.0;
+    double curSeg1Val2 = (it->first >= 0) ? val1[it->first] : 0.0;
+
+    // Hist 2
+    // Start of the gradient section in Seg1
+    double curSeg2Loc1 = segs.loc2_left(it->second);
+
+    // End of the gradient section in Seg1
+    double curSeg2Loc2 = segs.loc2_mid(it->second);
+
+    // End of the flat section in Seg1
+    double curSeg2Loc3 = segs.loc2_right(it->second);
+
+    // Start and end values in Seg2: val2 gives the values at *right*
+    // endpoints of the segments.  A value of 0.0 is used before the
+    // first segment.
+    double curSeg2Val1 = (it->second > 0) ? val2[it->second - 1] : 0.0;
+    double curSeg2Val2 = (it->second >= 0) ? val2[it->second] : 0.0;
+
+    double element = get_double_segment_constrained(
+      curSeg1Loc1, curSeg1Loc2, curSeg1Loc3, curSeg1Val1, curSeg1Val2,
+      curSeg2Loc1, curSeg2Loc2, curSeg2Loc3, curSeg2Val1, curSeg2Val2);
+
+    add_element_kahan(res, element, compensation);
+  }
+
+  return res;
+}