ricker continuous wavelet transform and correct gaussian1d convolution.

java/tech/v3/datatype/Convolve1D.java

@@ -3,6 +3,7 @@
 import clojure.lang.RT;
 import java.util.function.BiFunction;
 import java.util.function.DoubleUnaryOperator;
+import java.util.Arrays;
 
 
 
@@ -15,8 +16,98 @@ public static enum Mode
   };
   public static enum EdgeMode {
     Clamp,
+    Nearest, //Synonym for Clamp
+    Reflect,
+    Constant,
+    //Mirror, - This one makes no sense
+    Wrap,
     Zero,
   }
+  public static class Edging {
+    public final EdgeMode mode;
+    public final double constant;
+    public Edging(EdgeMode mode) {
+      this(mode,0.0);
+    }
+    public Edging(EdgeMode _mode, double _constant) {
+      mode = _mode;
+      constant = _constant;
+    }
+
+    public double[] apply(double[] data, int nWin, Mode mode) {
+      final int extraLen;
+      final int nWinDec = nWin - 1;
+      switch(mode) {
+      case Full: extraLen = nWinDec * 2; break;
+	//Convolving a even sized window means we have to distribute 1 extra on one
+	//side of the input or another.
+      case Same: extraLen = (nWin % 2) == 0 ? nWinDec : nWin; break;
+      case Valid: extraLen = 0; break;
+      default: throw new RuntimeException("Unrecognized mode.");
+      }
+
+      return apply(data, data.length + extraLen);
+    }
+
+    public double[] apply(double[] data, int new_len) {
+      int len = data.length;
+      if (len == new_len)
+	return data;
+      int left_overflow = (new_len - len) / 2;
+      int right_overflow = new_len - left_overflow - len;
+      int max_overflow = Math.max(left_overflow, right_overflow);
+      int end_start = len + left_overflow;
+      double[] retval = new double[new_len];
+      System.arraycopy(data,0,retval,left_overflow,len);
+      switch(mode) {
+      case Constant:
+      case Zero:
+	if(constant != 0.0) {
+	  Arrays.fill(retval, 0, left_overflow, constant);
+	  Arrays.fill(retval, left_overflow+len, new_len - 1, constant);
+	}
+	break;
+      case Nearest:
+      case Clamp:
+	double left = data[0];
+	double right = data[len-1];
+	Arrays.fill(retval, 0, left_overflow, left);
+	Arrays.fill(retval, left_overflow+len, new_len - 1, right);
+	break;
+
+      case Reflect:
+	/* abcddcba|abcd|dcbaabcd */
+	for(int idx = 0; idx < max_overflow; ++idx) {
+	  boolean even = ((idx / len) % 2) == 0;
+	  int ary_pos = idx % len;
+	  if (! even) ary_pos = len - ary_pos - 1;
+
+	  if (idx < left_overflow)
+	    retval[left_overflow - idx - 1] = data[ary_pos];
+
+	  int write_pos = end_start + idx;
+	  //Odd values give off-by-one
+	  if (write_pos < new_len)
+	    retval[write_pos] = data[len-ary_pos-1];
+	}
+	break;
+      case Wrap:
+	/* abcdabcd|abcd|abcdabcd */
+	for(int idx = 0; idx < max_overflow; ++idx) {
+	  int ary_pos = idx % len;
+	  if (idx < left_overflow)
+	    retval[left_overflow - idx - 1] = data[len - ary_pos - 1];
+
+	  int write_pos = end_start + idx;
+	  //Odd values give off-by-one
+	  if (write_pos < new_len)
+	    retval[write_pos] = data[ary_pos];
+	}
+	break;
+      }
+      return retval;
+    }
+  }
   //Precondition is that ndata is <= nwindow
   //The difference between convolve and correlate is convolve reverses
   //the window.  Callers can reverse the window before passing it in.
@@ -29,12 +120,10 @@ public static enum EdgeMode {
      int stepsize,
      DoubleUnaryOperator in_finalize,
      Mode mode,
-     EdgeMode edgeMode) {
+     Edging edging) {
 
     final int nData = data.length;
     final int nWin = window.length;
-    if (nData < nWin)
-      throw new RuntimeException("Data size must be <= window size");
     final int nWinDec = nWin - 1;
     final int nResult;
     switch (mode) {
@@ -52,24 +141,7 @@ public static enum EdgeMode {
     final DoubleUnaryOperator finalize = in_finalize;
 
 
-    final int windowOffset;
-    switch(mode) {
-    case Full: windowOffset = nWinDec; break;
-    case Same: windowOffset = (nWinDec / 2); break;
-    case Valid: windowOffset = 0; break;
-    default: throw new RuntimeException("Unrecognized mode.");
-    }
-
-
-    final double left_edge;
-    final double right_edge;
-    switch(edgeMode) {
-    case Clamp: left_edge=data[0]; right_edge=data[nData-1]; break;
-    case Zero: left_edge=0.0; right_edge=0.0; break;
-    default: throw new RuntimeException("Unrecognized edge mode");
-    }
-
-
+    double[] conv_data = edging.apply(data, nWin, mode);
     double[] retval = new double[nResult];
     parallelizer.apply( new Long(nResult), new BiFunction<Long,Long,Object> () {
 	public Object apply(Long sidx, Long glen) {
@@ -79,18 +151,13 @@ public Object apply(Long sidx, Long glen) {
 
 	  //This loop is an interesting case for a blog post on the vectorization API.
 	  for(int idx = start_idx; idx < end_idx; ++idx ) {
-	    int data_start_idx = (idx*stepsize) - windowOffset;
+	    int data_start_idx = (idx*stepsize);
 	    double sum = 0.0;
 	    for (int win_idx = 0; win_idx < nWin; ++win_idx) {
 	      int data_idx = data_start_idx + win_idx;
-	      double data_val;
-	      if (data_idx < 0)
-		data_val = left_edge;
-	      else if (data_idx >= nData)
-		data_val = right_edge;
-	      else
-		data_val = data[data_idx];
-	      sum += data_val * window[win_idx];
+	      //System.out.println(String.format("idx %d, win_idx %d, data_idx %d, data_len %d",
+	      //			       idx, win_idx, data_idx, conv_len));
+	      sum += conv_data[data_idx] * window[win_idx];
 	    }
 	    retval[idx] = finalize.applyAsDouble(sum);
 	  }

project.clj

@@ -23,7 +23,7 @@
                              [com.taoensso/nippy "3.1.1"]
                              [ch.qos.logback/logback-classic "1.1.3"]
                              [com.clojure-goes-fast/clj-memory-meter "0.1.0"]
-                             [techascent/tech.viz "6.00-beta-16"]]
+                             [techascent/tech.viz "6.00-beta-16-1"]]
               :test-paths ["neanderthal" "test"]}
              :jdk-16 {:jvm-opts ["--add-modules" "jdk.incubator.foreign"
                                  "-Dforeign.restricted=permit"