changeset 546:17d5a8986f6f

Make vector multiply take a list of vectors, not just two (in symmetry with add); slightly simplify implementations so as to return the arg unmodified if there is only one
author Chris Cannam
date Mon, 24 Mar 2014 09:59:27 +0000
parents 01863795221c
children 1200e1029ecb
files src/may/bits/VectorBits.java src/may/mathmisc.yeti src/may/matrix.yeti src/may/signal/window.yeti src/may/stream/filter.yeti src/may/stream/framer.yeti src/may/stream/resample.yeti src/may/stream/test/test_framer.yeti src/may/vector.yeti src/may/vector/test/test_vector.yeti
diffstat 10 files changed, 71 insertions(+), 57 deletions(-) [+]
line wrap: on
line diff
--- a/src/may/bits/VectorBits.java	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/bits/VectorBits.java	Mon Mar 24 09:59:27 2014 +0000
@@ -12,36 +12,33 @@
 	return tot;
     }
 
-    public static double[] multiply(double[] v1, double[] v2) {
-	int len = v1.length;
-	if (v2.length < len) {
-	    len = v2.length;
+    public static void multiplyBy(double[] out, double[] in) {
+	for (int i = 0; i < in.length && i < out.length; ++i) {
+	    out[i] *= in[i];
 	}
-	double[] out = new double[len];
-	for (int i = 0; i < len; ++i) {
-	    out[i] = v1[i] * v2[i];
+	for (int i = in.length; i < out.length; ++i) {
+	    out[i] *= 0.0;
 	}
-	return out;
     }
 
-    public static double[] divide(double[] v1, double[] v2) {
-	int len = v1.length;
-	if (v2.length < len) {
-	    len = v2.length;
+    public static void divideBy(double[] out, double[] in) {
+	for (int i = 0; i < in.length && i < out.length; ++i) {
+	    out[i] /= in[i];
 	}
-	double[] out = new double[len];
-	for (int i = 0; i < len; ++i) {
-	    out[i] = v1[i] / v2[i];
-	}
-	return out;
     }
 
-    public static void addTo(double[] out, double[] in, int len) {
-	for (int i = 0; i < len; ++i) {
+    public static void addTo(double[] out, double[] in) {
+	for (int i = 0; i < in.length && i < out.length; ++i) {
 	    out[i] += in[i];
 	}
     }
 
+    public static void subtractFrom(double[] out, double[] in) {
+	for (int i = 0; i < in.length && i < out.length; ++i) {
+	    out[i] -= in[i];
+	}
+    }
+
     public static double[] scaled(double[] v, double factor) {
 	int len = v.length;
 	double[] out = new double[v.length];
--- a/src/may/mathmisc.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/mathmisc.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -65,6 +65,9 @@
         np x 1;
     fi;
 
+random () =
+    Math#random();
+
 {
     gcd,
     factorial,
@@ -76,5 +79,6 @@
     log2,
     isPowerOfTwo,
     nextPowerOfTwo,
+    random,
 }
 
--- a/src/may/matrix.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/matrix.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -653,7 +653,7 @@
     for [0..size.rows - 1] do i:
         row = getRow i m1;
         for [0..size.columns - 1] do j:
-            data[j][i] := vec.sum (vec.multiply row (getColumn j m2));
+            data[j][i] := vec.sum (vec.multiply [row, getColumn j m2]);
         done;
     done;
     newMatrix size (Columns (array (map vec.vector (list data)))));
@@ -689,9 +689,11 @@
             entryWiseProduct m2 m1
         else
             if isRowMajor? m1 then
-                fromRows (array (map2 vec.multiply (asRows m1) (asRows m2)));
+                fromRows (array (map2 do v1 v2: vec.multiply [v1,v2] done
+                   (asRows m1) (asRows m2)));
             else
-                fromColumns (array (map2 vec.multiply (asColumns m1) (asColumns m2)));
+                fromColumns (array (map2 do v1 v2: vec.multiply [v1,v2] done
+                   (asColumns m1) (asColumns m2)));
             fi
         fi
     fi;
--- a/src/may/signal/window.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/signal/window.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -208,19 +208,19 @@
 
 //!!! doc note: recalculates window every time! not fast for multiple blocks, consider windowedFrames then
 windowed windowFunc data =
-    vec.multiply (windowFunc (vec.length data)) data;
+    vec.multiply [windowFunc (vec.length data), data];
 
 windowedRows windowFunc matrix = 
    (w = windowFunc (mat.width matrix);
     mat.fromRows
-       (map (vec.multiply w) (mat.asRows matrix)));
+       (map do v: vec.multiply [v,w] done (mat.asRows matrix)));
 
 windowedFrames windowFunc frames =
     case frames of
         []: frames;
          _: (first = head frames;
              window = windowFunc (vec.length first);
-             map (vec.multiply window) frames);
+             map do v: vec.multiply [v,window] done frames);
     esac;
 
 {
--- a/src/may/stream/filter.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/stream/filter.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -27,7 +27,7 @@
              else win.dirac filterLength;
              fi;
          fi;
-    filter = vec.multiply idealBandpass kw;
+    filter = vec.multiply [idealBandpass, kw];
     filtered = convolve.convolvedWith [Framesize 1024]
        (mat.fromRows (map \filter [1..s.channels]))
         s;
--- a/src/may/stream/framer.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/stream/framer.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -86,7 +86,8 @@
         Hop h: hop := h;
         Padded p: padded := p;
         Window w:
-            windower := mat.mapRows (vec.multiply (w framesize));
+           (window = w framesize;
+            windower := mat.mapRows (do v: vec.multiply [v, window] done));
         FrequencyDomain f:
             if f then
                 transform := mat.mapRows (fft.realForwardMagnitude framesize);
@@ -106,7 +107,8 @@
         Hop h: hop := h;
         Padded p: padded := p;
         Window w:
-            windower := mat.mapRows (vec.multiply (w framesize));
+           (window = w framesize;
+            windower := mat.mapRows (do v: vec.multiply [v, window] done));
         FrequencyDomain f:
         //!!! what if both f and not-f provided in one options list? need reset
             if f then 
@@ -196,7 +198,7 @@
                 reverse acc
             else
                 this = mat.resizedTo { columns = framesize, rows = channels }
-                   (mat.mapRows (vec.multiply w) (head remaining));
+                   (mat.mapRows do v: vec.multiply [v,w] done (head remaining));
                 remaining := tail remaining;
                 framesFor (samples - hop) (this::acc)
             fi;
--- a/src/may/stream/resample.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/stream/resample.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -25,7 +25,7 @@
     length = if length > maxlength then maxlength else length fi;
     kw = win.kaiser (pp with { length });
     sw = win.sinc (n*2) length;
-    vec.multiply sw kw);
+    vec.multiply [sw, kw]);
 
 durationAdjusterFor factor s =
     case s.available of
@@ -298,7 +298,7 @@
 
                 for [0..s.channels-1] do ch:
                     rowdata[ch][i] := vec.sum
-                       (vec.multiply (mat.getRow ch buffer) pd.filter);
+                       (vec.multiply [mat.getRow ch buffer, pd.filter]);
                 done;
             
                 buffer := mat.concatHorizontal
--- a/src/may/stream/test/test_framer.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/stream/test/test_framer.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -32,7 +32,7 @@
         w = vec.scaled factor (win.hann framesize);
         syn.precalculatedMono rate
            (vec.concat [
-                vec.multiply (vec.fromList [1..m]) (vec.slice w 0 m),
+                vec.multiply [vec.fromList [1..m], vec.slice w 0 m],
                 vec.fromList [m+1..n]
             ])
     fi);
--- a/src/may/vector.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/vector.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -204,27 +204,31 @@
         len: sum' v / len
     esac;
 
+listOp f vv =
+    case vv of
+    [v]: v;
+    v::rest:
+       (out = copyOf v;
+        for rest (f out);
+        out);
+    _: failWith "Empty argument list";
+    esac;
+
+//!!! doc: returned vector is same length as first argument (this has changed, formerly it was length of shortest argument)
 add vv is list?<~double[]> -> ~double[] =
-   (len = head (sort (map length' vv));
-    out = new double[len];
-    for vv do v:
-        VectorBits#addTo(out, v, len);
-    done;
-    out);
+    listOp do out v: VectorBits#addTo(out, v) done vv;
 
+//!!! doc: returned vector is same length as first argument (this has changed, formerly it was length of shortest argument)
 subtract v1 v2 is ~double[] -> ~double[] -> ~double[] =
-   (len = if length v1 < length v2 then length v1 else length v2 fi;
-    out = new double[len];
-    for [0..len-1] do i:
-        out[i] := v1[i] - v2[i]
-    done;
-    out);
+    listOp do out v: VectorBits#subtractFrom(out, v) done [v1, v2];
 
-multiply b1 b2 is ~double[] -> ~double[] -> ~double[] = 
-    VectorBits#multiply(b1, b2);
+//!!! doc: returned vector is same length as first argument (this has changed, formerly it was length of shortest argument). If first arg is longer than others, the spare values will become zero
+multiply vv is list?<~double[]> -> ~double[] =
+    listOp do out v: VectorBits#multiplyBy(out, v) done vv;
 
-divide b1 b2 is ~double[] -> ~double[] -> ~double[] = 
-    VectorBits#divide(b1, b2);
+//!!! doc: returned vector is same length as first argument (this has changed, formerly it was length of shortest argument). If first arg is longer than second, the spare values are left untouched (not divided by zero)
+divide v1 v2 is ~double[] -> ~double[] -> ~double[] = 
+    listOp do out v: VectorBits#divideBy(out, v) done [v1, v2];
 
 scaled n v is number -> ~double[] -> ~double[] =
     if n == 1 then v
@@ -239,7 +243,9 @@
     fi;
 
 sqr v =
-    multiply v v;
+   (out = copyOf v;
+    VectorBits#multiplyBy(out, v);
+    out);
 
 rms =
     sqrt . mean . sqr;
@@ -358,7 +364,7 @@
     mean is vector_t -> number,
     add is list?<vector_t> -> vector_t,
     subtract is vector_t -> vector_t -> vector_t,
-    multiply is vector_t -> vector_t -> vector_t, 
+    multiply is list?<vector_t> -> vector_t, 
     divide is vector_t -> vector_t -> vector_t, 
     scaled is number -> vector_t -> vector_t,
     divideBy is number -> vector_t -> vector_t, 
--- a/src/may/vector/test/test_vector.yeti	Fri Mar 21 17:32:48 2014 +0000
+++ b/src/may/vector/test/test_vector.yeti	Mon Mar 24 09:59:27 2014 +0000
@@ -168,24 +168,27 @@
 
 "add": \(
     compare (vec.list (vec.add [vec.zeros 0, vec.ones 5])) [] and
-        compare (vec.list (vec.add [vec.consts 3 4, vec.fromList [1,2,3] ])) [4,5,6] and
-        compare (vec.list (vec.add [vec.consts (-3) 4, vec.fromList [1,2,3] ])) [-2,-1,0] 
+        compare (vec.list (vec.add [vec.consts 3 4, vec.fromList [1,2,3] ])) [4,5,6,3] and
+        compare (vec.list (vec.add [vec.consts (-3) 4, vec.fromList [1,2,3] ])) [-2,-1,0,-3] and
+        compare (vec.list (vec.add [vec.consts 3 3, vec.fromList [1,2,3], vec.fromList [6,7,8,9] ])) [10,12,14] 
 ),
 
 "subtract": \(
     compare (vec.list (vec.subtract (vec.zeros 0) (vec.ones 5))) [] and
-        compare (vec.list (vec.subtract (vec.consts 3 4) (vec.fromList [1,2,3]))) [2,1,0] and
-        compare (vec.list (vec.subtract (vec.consts (-3) 4) (vec.fromList [1,2,3]))) [-4,-5,-6]
+        compare (vec.list (vec.subtract (vec.consts 3 4) (vec.fromList [1,2,3]))) [2,1,0,3] and
+        compare (vec.list (vec.subtract (vec.consts (-3) 4) (vec.fromList [1,2,3]))) [-4,-5,-6,-3]
 ),
 
 "multiply": \(
-    compare (vec.list (vec.multiply (vec.zeros 0) (vec.ones 5))) [] and
-        compare (vec.list (vec.multiply (vec.consts (-3) 4) (vec.fromList [1,2,3]))) [-3,-6,-9]
+    compare (vec.list (vec.multiply [vec.zeros 0, vec.ones 5])) [] and
+        compare (vec.list (vec.multiply [vec.consts (-3) 4, vec.fromList [1,2,3]])) [-3,-6,-9,0] and
+        compare (vec.list (vec.multiply [vec.consts 3 3, vec.fromList [1,2,3], vec.fromList [6,7,8,9]])) [18,42,72]
 ),
 
 "divide": \(
     compare (vec.list (vec.divide (vec.zeros 0) (vec.ones 5))) [] and
-        compare (vec.list (vec.divide (vec.consts (-3) 4) (vec.fromList [1,2,3]))) [-3,-(3/2),-1]
+        compare (vec.list (vec.divide (vec.consts (-3) 3) (vec.fromList [1,2,3]))) [-3,-(3/2),-1] and
+        compare (vec.list (vec.divide (vec.consts (-3) 4) (vec.fromList [1,2,3]))) [-3,-(3/2),-1,-3] // values in first vec beyond end of second are *not* divided by zero, just left alone
 ),
 
 "scaled": \(