Mercurial > hg > may
changeset 559:ab85dfc45e3f
Merge
| author | Chris Cannam |
|---|---|
| date | Sat, 05 Apr 2014 16:02:24 +0100 |
| parents | eeebedb84ca7 (current diff) 6a5c20ceb949 (diff) |
| children | 9a40a9eedba7 |
| files | src/may/stream/test/test_framer.yeti |
| diffstat | 31 files changed, 725 insertions(+), 241 deletions(-) [+] |
line wrap: on
line diff
--- a/.hgsubstate Mon Mar 03 17:03:32 2014 +0000 +++ b/.hgsubstate Sat Apr 05 16:02:24 2014 +0100 @@ -1,1 +1,1 @@ -0a4af51413b7ddd1cc4c058fb732556b849cac21 ext +ea762d69743ec57dd9f10dbfc3ff87ffacf64565 ext
--- a/bin/may Mon Mar 03 17:03:32 2014 +0000 +++ b/bin/may Sat Apr 05 16:02:24 2014 +0100 @@ -30,7 +30,7 @@ YETI_LIBDIR=$YLDIR/../other/yeti fi -CLASSPATH=$YLDIR/may.jar:$JARDIR/ayr.jar:$YETI_LIBDIR/yeti.jar:$YETI_LIBDIR/yeti-lib.jar:$JARDIR/jline-2.11-SNAPSHOT.jar:$JARDIR/jvamp.jar:$JARDIR/yertle.jar:$JARDIR/jtransforms-2.4.jar:$JARDIR/org.jzy3d-0.9.jar:$JARDIR/jogl-all.jar:$JARDIR/gluegen.jar:$JARDIR/gluegen-rt.jar:$JARDIR/opencsv-2.1.jar:$JARDIR/org.convexhull.jar:$CLASSPATH +CLASSPATH=$YLDIR/may.jar:$JARDIR/ayr.jar:$YETI_LIBDIR/yeti.jar:$YETI_LIBDIR/yeti-lib.jar:$JARDIR/jline-2.11-SNAPSHOT.jar:$JARDIR/jvamp.jar:$JARDIR/yertle.jar:$JARDIR/jtransforms-2.4.jar:$JARDIR/jzy3d-swt-0.9.1.jar:$JARDIR/jzy3d-api-0.9.1.jar:$JARDIR/jogl-all.jar:$JARDIR/gluegen.jar:$JARDIR/gluegen-rt.jar:$JARDIR/opencsv-2.1.jar:$JARDIR/org.convexhull.jar:$CLASSPATH YETI_MODULE_SOURCE_PATH=${YETI_LIBDIR}/modules \ LD_LIBRARY_PATH=$SODIR:$LD_LIBRARY_PATH \
--- a/bin/may.bat Mon Mar 03 17:03:32 2014 +0000 +++ b/bin/may.bat Sat Apr 05 16:02:24 2014 +0100 @@ -1,3 +1,3 @@ @echo off set MAY_INIT_MODULES=may.vector;may.matrix;may.complex;may.plot;may.vamp -java -classpath "may.jar;ext\jar\ayr.jar;..\yeti\yeti.jar;ext\jar\jline-2.11-SNAPSHOT.jar;ext\jar\jvamp.jar;ext\jar\yertle.jar;ext\jar\jtransforms-2.4.jar;ext\jar\org.jzy3d-0.9.jar;ext\jar\jogl-all.jar;ext\jar\gluegen.jar;ext\jar\opencsv-2.1.jar;ext\jar\org.convexhull.jar" -Djava.library.path=ext\native\win64 com.particularprograms.ayr -init "%MAY_INIT_MODULES%" +java -classpath "may.jar;ext\jar\ayr.jar;..\yeti\yeti.jar;ext\jar\jline-2.11-SNAPSHOT.jar;ext\jar\jvamp.jar;ext\jar\yertle.jar;ext\jar\jtransforms-2.4.jar;ext\jar\jzy3d-swt-0.9.1.jar;ext\jar\jzy3d-api-0.9.1.jar;ext\jar\jogl-all.jar;ext\jar\gluegen.jar;ext\jar\opencsv-2.1.jar;ext\jar\org.convexhull.jar" -Djava.library.path=ext\native\win64 com.particularprograms.ayr -init "%MAY_INIT_MODULES%"
--- a/bin/yc Mon Mar 03 17:03:32 2014 +0000 +++ b/bin/yc Sat Apr 05 16:02:24 2014 +0100 @@ -34,7 +34,7 @@ shift fi -CLASSPATH=${YLJAR_WITH_COLON}$YETI_LIBDIR/yeti.jar:$YETI_LIBDIR/yeti-lib.jar:$JARDIR/jvamp.jar:$JARDIR/yertle.jar:$JARDIR/jtransforms-2.4.jar:$JARDIR/org.jzy3d-0.9.jar:$JARDIR/jogl-all.jar:$JARDIR/gluegen.jar:$JARDIR/gluegen-rt.jar:$JARDIR/misc/opencsv-2.1.jar:$JARDIR/misc/org.convexhull.jar:$JARDIR/misc/swt.jar:$CLASSPATH +CLASSPATH=${YLJAR_WITH_COLON}$YETI_LIBDIR/yeti.jar:$YETI_LIBDIR/yeti-lib.jar:$JARDIR/jvamp.jar:$JARDIR/yertle.jar:$JARDIR/jtransforms-2.4.jar:$JARDIR/jzy3d-swt-0.9.1.jar:$JARDIR/jzy3d-api-0.9.1.jar:$JARDIR/jogl-all.jar:$JARDIR/gluegen.jar:$JARDIR/gluegen-rt.jar:$JARDIR/misc/opencsv-2.1.jar:$JARDIR/misc/org.convexhull.jar:$JARDIR/misc/swt.jar:$CLASSPATH LD_LIBRARY_PATH=$SODIR:$LD_LIBRARY_PATH \ $JAVA_HOME/bin/java $JAVA_OPTS -classpath "$CLASSPATH" yeti.lang.compiler.yeti "$@"
--- a/bin/yc.bat Mon Mar 03 17:03:32 2014 +0000 +++ b/bin/yc.bat Sat Apr 05 16:02:24 2014 +0100 @@ -1,2 +1,2 @@ -java -classpath "may.jar;..\yeti\yeti.jar;ext\jar\jline-2.11-SNAPSHOT.jar;ext\jar\jvamp.jar;ext\jar\yertle.jar;ext\jar\jtransforms-2.4.jar;ext\jar\org.jzy3d-0.9.jar;ext\jar\jogl-all.jar;ext\jar\gluegen.jar;ext\jar\opencsv-2.1.jar;ext\jar\org.convexhull.jar" -Djava.library.path=ext\native\win64 yeti.lang.compiler.yeti "%~1" +java -classpath "may.jar;..\yeti\yeti.jar;ext\jar\jline-2.11-SNAPSHOT.jar;ext\jar\jvamp.jar;ext\jar\yertle.jar;ext\jar\jtransforms-2.4.jar;ext\jar\jzy3d-swt-0.9.1.jar;ext\jar\jzy3d-api-0.9.1.jar;ext\jar\jogl-all.jar;ext\jar\gluegen.jar;ext\jar\opencsv-2.1.jar;ext\jar\org.convexhull.jar" -Djava.library.path=ext\native\win64 yeti.lang.compiler.yeti "%~1"
--- a/build.xml Mon Mar 03 17:03:32 2014 +0000 +++ b/build.xml Sat Apr 05 16:02:24 2014 +0100 @@ -2,7 +2,7 @@ <property name="jardir" value="${basedir}/ext/jar"/> - <property name="extjars" value="${jardir}/jvamp.jar:${jardir}/yertle.jar:${jardir}/jtransforms-2.4.jar:${jardir}/org.jzy3d-0.9.jar:${jardir}/jogl-all.jar:${jardir}/gluegen.jar:${jardir}/gluegen-rt.jar:${jardir}/opencsv-2.1.jar:${jardir}/org.convexhull.jar"/> + <property name="extjars" value="${jardir}/jvamp.jar:${jardir}/yertle.jar:${jardir}/jtransforms-2.4.jar:${jardir}/jzy3d-swt-0.9.1.jar:${jardir}/jzy3d-api-0.9.1.jar:${jardir}/jogl-all.jar:${jardir}/gluegen.jar:${jardir}/gluegen-rt.jar:${jardir}/opencsv-2.1.jar:${jardir}/org.convexhull.jar"/> <condition property="archtag" value="linux32"> <os family="unix" arch="i386"/>
--- a/src/may/bits/VectorBits.java Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/bits/VectorBits.java Sat Apr 05 16:02:24 2014 +0100 @@ -3,6 +3,15 @@ public class VectorBits { + public static void checkLengths(double[] v1, double[] v2) { + if (v1.length != v2.length) { + throw new IllegalArgumentException + ("Found vector of length " + v2.length + + ", but all so far in this arithmetic operation have had length " + + v1.length); + } + } + public static double sum(double[] v) { double tot = 0.0; int len = v.length; @@ -12,24 +21,37 @@ 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) { + checkLengths(out, 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 void addTo(double[] out, double[] in, int len) { - for (int i = 0; i < len; ++i) { + public static void divideBy(double[] out, double[] in) { + checkLengths(out, in); + for (int i = 0; i < in.length && i < out.length; ++i) { + out[i] /= in[i]; + } + } + + public static void addTo(double[] out, double[] in) { + checkLengths(out, in); + for (int i = 0; i < in.length && i < out.length; ++i) { out[i] += in[i]; } } + public static void subtractFrom(double[] out, double[] in) { + checkLengths(out, 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/complex.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/complex.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -2,7 +2,6 @@ module may.complex; vec = load may.vector; -mm = load may.mathmisc; import java.lang: ClassCastException; @@ -31,13 +30,16 @@ fi, String format() (f v r = - //!!! not right -- want to round to 1dp of exp notation if less than 1e-4 - (n = mm.round (v * 10000); - s = "\(n / 10000)"; + (s = if v == 0.0 then "0.0" + elif abs v >= 1000.0 or abs v < 0.01 then + String#format("%.1E", [v as ~Double]) + else + String#format("%4f", [v as ~Double]) + fi; if r then - (strPad ' ' (8 - strLength s) '') ^ s + (strPad ' ' (9 - strLength s) '') ^ s else - strPad ' ' 7 "\(s)i"; + strPad ' ' 8 "\(s)i"; fi); if imag < 0 then " \(f real true)-\(f (-imag) false)"
--- a/src/may/feature/specdiff.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/feature/specdiff.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -3,7 +3,6 @@ fr = load may.stream.framer; ch = load may.stream.channels; -cm = load may.matrix.complex; load may.feature.feature;
--- a/src/may/mathmisc.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/mathmisc.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -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 Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/matrix.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -40,8 +40,11 @@ module may.matrix; +{ ceil, floor, random } = load may.mathmisc; + vec = load may.vector; -mm = load may.mathmisc; + +load yeti.experimental.json; typedef opaque matrix_t = { size is { rows is number, columns is number }, @@ -361,7 +364,7 @@ }; constMatrix n = generate do row col: n done; -randomMatrix = generate do row col: Math#random() done; +randomMatrix = generate do row col: random () done; identityMatrix = constMatrix 1; transposed m = @@ -503,16 +506,22 @@ done (keys h)); newSparseMatching m1 entries); -sum' m1 m2 = - if (size m1) != (size m2) - then failWith "Matrices are not the same size: \(size m1), \(size m2)"; - elif isSparse? m1 and isSparse? m2 then - sparseSumOrDifference (+) m1 m2; - else - add2 v1 v2 = vec.add [v1,v2]; - denseLinearOp add2 m1 m2; - fi; - +sum' mm = + case mm of + m1::m2::rest: + sum' + (if (size m1) != (size m2) + then failWith "Matrices are not the same size: \(size m1), \(size m2)"; + elif isSparse? m1 and isSparse? m2 then + sparseSumOrDifference (+) m1 m2; + else + add2 v1 v2 = vec.add [v1,v2]; + denseLinearOp add2 m1 m2; + fi :: rest); + [m1]: m1; + _: failWith "Empty argument list"; + esac; + difference m1 m2 = if (size m1) != (size m2) then failWith "Matrices are not the same size: \(size m1), \(size m2)"; @@ -652,7 +661,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))))); @@ -675,11 +684,52 @@ fi; fi; -entryWiseProduct m1 m2 = // or element-wise, or Hadamard product -//!!! todo: faster, sparse version, unit tests +entryWiseProduct mm = + case mm of + m1::m2::rest: + entryWiseProduct + (if (size m1) != (size m2) + then failWith "Matrices are not the same size: \(size m1), \(size m2)"; + else + if isSparse? m1 then + newSparse (size m1) + ((taggerForTypeOf m1) + (map do { i, j, v }: { i, j, v = v * (at' m2 i j) } done + (enumerateSparse m1))) + elif isSparse? m2 then + entryWiseProduct (m2::m1::rest) + else + if isRowMajor? m1 then + fromRows (array (map2 do v1 v2: vec.multiply [v1,v2] done + (asRows m1) (asRows m2))); + else + fromColumns (array (map2 do v1 v2: vec.multiply [v1,v2] done + (asColumns m1) (asColumns m2))); + fi + fi + fi :: rest); + [m1]: m1; + _: failWith "Empty argument list"; + esac; + +entryWiseDivide m1 m2 = if (size m1) != (size m2) then failWith "Matrices are not the same size: \(size m1), \(size m2)"; - else generate do row col: at' m1 row col * at' m2 row col done (size m1); + else + if isSparse? m1 then + newSparse (size m1) + ((taggerForTypeOf m1) + (map do { i, j, v }: { i, j, v = v / (at' m2 i j) } done + (enumerateSparse m1))) + // For m2 to be sparse makes no sense (divide by zero all over + // the shop). + else + if isRowMajor? m1 then + fromRows (array (map2 vec.divide (asRows m1) (asRows m2))); + else + fromColumns (array (map2 vec.divide (asColumns m1) (asColumns m2))); + fi + fi fi; concatAgainstGrain tagger getter counter mm = @@ -691,11 +741,10 @@ concatWithGrain tagger getter counter mm = tagger (array - (concat - (map do m: - n = counter (size m); - map do i: getter i m done [0..n-1] - done mm))); + (concatMap do m: + n = counter (size m); + map do i: getter i m done [0..n-1] + done mm)); sparseConcat direction first mm = (dimension d f = if direction == d then sum (map f mm) else f first fi; @@ -850,18 +899,71 @@ fi fi); +//!!! doc: always dense +repeatedHorizontal n m = + (if n == 1 then m + else + cols = asColumns m; + fromColumns (fold do acc _: acc ++ cols done [] [1..n]) + fi); + +//!!! doc: always dense +repeatedVertical n m = + (if n == 1 then m + else + rows = asRows m; + fromRows (fold do acc _: acc ++ rows done [] [1..n]) + fi); + +//!!! doc: always dense +tiledTo newsize m = + if newsize == size m then + m + elif (height m) == 0 or (width m) == 0 then + zeroMatrixWithTypeOf m newsize; + else + h = ceil (newsize.columns / (width m)); + v = ceil (newsize.rows / (height m)); + if isRowMajor? m then + resizedTo newsize (repeatedHorizontal h (repeatedVertical v m)) + else + resizedTo newsize (repeatedVertical v (repeatedHorizontal h m)) + fi + fi; + minValue m = if width m == 0 or height m == 0 then 0 - else + elif isSparse? m then minv ll = fold min (head ll) (tail ll); - minv (map (.v) (enumerate m)); + minnz = minv (map (.v) (enumerate m)); + if minnz > 0 and nonZeroValues m < (width m * height m) then 0 + else minnz fi; + elif isRowMajor? m then + vec.min (vec.fromList (map vec.min (asRows m))); + else + vec.min (vec.fromList (map vec.min (asColumns m))); fi; maxValue m = if width m == 0 or height m == 0 then 0 - else + elif isSparse? m then maxv ll = fold max (head ll) (tail ll); - maxv (map (.v) (enumerate m)); + maxnz = maxv (map (.v) (enumerate m)); + if maxnz < 0 and nonZeroValues m < (width m * height m) then 0 + else maxnz fi; + elif isRowMajor? m then + vec.max (vec.fromList (map vec.max (asRows m))); + else + vec.max (vec.fromList (map vec.max (asColumns m))); + fi; + +total m = + if isSparse? m then + fold (+) 0 (map (.v) (enumerateSparse m)); + elif isRowMajor? m then + fold (+) 0 (map vec.sum (asRows m)); + else + fold (+) 0 (map vec.sum (asColumns m)); fi; mapRows rf m = @@ -880,21 +982,23 @@ [ "\nColumns \(c0) to \(c1)\n", (map do row: map do v: - //!!! not right -- want to round to 1dp of exp notation if less than 1e-4 - n = mm.round (v * 100000); - strPad ' ' 10 "\(n / 100000)"; + strPad ' ' 10 + (if v == 0 then "0.0" + elif abs v >= 1000.0 or abs v < 0.01 then + String#format("%.2E", [v as ~Double]) + else + String#format("%5f", [v as ~Double]) + fi); done (vec.list row) |> strJoin ""; done (asRows (columnSlice m c0 (c1 + 1))) |> strJoin "\n") ]; - done [0..width m / chunk - 1] |> concat); + done [0..floor(width m / chunk)] |> concat); print' = println . format; eprint' = eprintln . format; -//!!! todo: look at all occurrences of matrix (and complexmatrix) -// construction and make sure we have good apis for those use cases. -// in particular, constructing from literal data (list<list<number>>) -// is much too hard at the moment. +json m = + jsonList (map do r: jsonList (map jsonNum (vec.list r)) done (asRows m)); { size, @@ -922,9 +1026,9 @@ toSparse, toDense, scaled, - resizedTo, minValue, maxValue, + total, asRows, asColumns, sum = sum', @@ -936,6 +1040,11 @@ any = any', product, entryWiseProduct, + entryWiseDivide, + resizedTo, + tiledTo, + repeatedHorizontal, + repeatedVertical, concatHorizontal, concatVertical, rowSlice, @@ -953,6 +1062,7 @@ format, print = print', eprint = eprint', + json, } as { @@ -981,12 +1091,12 @@ toSparse is matrix_t -> matrix_t, toDense is matrix_t -> matrix_t, scaled is number -> matrix_t -> matrix_t, - resizedTo is { rows is number, columns is number } -> matrix_t -> matrix_t, minValue is matrix_t -> number, maxValue is matrix_t -> number, + total is matrix_t -> number, asRows is matrix_t -> list<vec.vector_t>, asColumns is matrix_t -> list<vec.vector_t>, - sum is matrix_t -> matrix_t -> matrix_t, + sum is list?<matrix_t> -> matrix_t, difference is matrix_t -> matrix_t -> matrix_t, abs is matrix_t -> matrix_t, negative is matrix_t -> matrix_t, @@ -994,7 +1104,12 @@ all is (number -> boolean) -> matrix_t -> boolean, any is (number -> boolean) -> matrix_t -> boolean, product is matrix_t -> matrix_t -> matrix_t, - entryWiseProduct is matrix_t -> matrix_t -> matrix_t, + entryWiseProduct is list?<matrix_t> -> matrix_t, + entryWiseDivide is matrix_t -> matrix_t -> matrix_t, + resizedTo is { rows is number, columns is number } -> matrix_t -> matrix_t, + tiledTo is { rows is number, columns is number } -> matrix_t -> matrix_t, + repeatedHorizontal is number -> matrix_t -> matrix_t, + repeatedVertical is number -> matrix_t -> matrix_t, concatHorizontal is list<matrix_t> -> matrix_t, concatVertical is list<matrix_t> -> matrix_t, rowSlice is matrix_t -> number -> number -> matrix_t, @@ -1012,5 +1127,6 @@ format is matrix_t -> string, print is matrix_t -> (), eprint is matrix_t -> (), + json is matrix_t -> json, }
--- a/src/may/matrix/complex.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/matrix/complex.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -5,6 +5,8 @@ vec = load may.vector; cpx = load may.complex; +{ floor } = load may.mathmisc; + typedef opaque complexmatrix_t = { size is { rows is number, columns is number }, real is Some mat.matrix_t | None (), @@ -204,7 +206,7 @@ case p1 of Some m1: case p2 of - Some m2: Some (mat.sum m1 m2); + Some m2: Some (mat.sum [m1, m2]); none: Some m1; esac; none: @@ -239,16 +241,22 @@ none; esac; -sum c1 c2 = - (a = c1.real; - b = c1.imaginary; - c = c2.real; - d = c2.imaginary; - { - size = c1.size, - real = addParts a c, - imaginary = addParts b d, - }); +sum cc = + case cc of + c1::c2::rest: + sum + (a = c1.real; + b = c1.imaginary; + c = c2.real; + d = c2.imaginary; + { + size = c1.size, + real = addParts a c, + imaginary = addParts b d, + } :: rest); + [c1]: c1; + _: failWith "Empty argument list"; + esac; product c1 c2 = (a = c1.real; @@ -305,7 +313,7 @@ (map do row: map cpx.format row |> strJoin ""; done (asRows (columnSlice m c0 (c1 + 1))) |> strJoin "\n") ]; - done [0..width m / chunk - 1] |> concat); + done [0..floor(width m / chunk)] |> concat); print' = println . format; eprint' = eprintln . format; @@ -374,7 +382,7 @@ toSparse is complexmatrix_t -> complexmatrix_t, toDense is complexmatrix_t -> complexmatrix_t, scaled is number -> complexmatrix_t -> complexmatrix_t, - sum is complexmatrix_t -> complexmatrix_t -> complexmatrix_t, + sum is list?<complexmatrix_t> -> complexmatrix_t, difference is complexmatrix_t -> complexmatrix_t -> complexmatrix_t, abs is complexmatrix_t -> mat.matrix_t, product is complexmatrix_t -> complexmatrix_t -> complexmatrix_t,
--- a/src/may/matrix/test/speedtest.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/matrix/test/speedtest.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -127,30 +127,30 @@ println ""; print "CMD + CMD... "; -reportOn (time \(mat.sum cmd cmd)); +reportOn (time \(mat.sum [cmd, cmd])); print "CMD + RMD... "; -reportOn (time \(mat.sum cmd rmd)); +reportOn (time \(mat.sum [cmd, rmd])); print "RMD + CMD... "; -reportOn (time \(mat.sum rmd cmd)); +reportOn (time \(mat.sum [rmd, cmd])); print "RMD + RMD... "; -reportOn (time \(mat.sum rmd rmd)); +reportOn (time \(mat.sum [rmd, rmd])); println ""; print "CMS + CMS... "; -reportOn (time \(mat.sum cms cms)); +reportOn (time \(mat.sum [cms, cms])); print "CMS + RMS... "; -reportOn (time \(mat.sum cms rms)); +reportOn (time \(mat.sum [cms, rms])); print "RMS + CMS... "; -reportOn (time \(mat.sum rms cms)); +reportOn (time \(mat.sum [rms, cms])); print "RMS + RMS... "; -reportOn (time \(mat.sum rms rms)); +reportOn (time \(mat.sum [rms, rms])); println ""; @@ -204,15 +204,15 @@ println ""; print "CMS + CMS... "; -reportOn (time \(mat.sum cms cms)); +reportOn (time \(mat.sum [cms, cms])); print "CMS + RMS... "; -reportOn (time \(mat.sum cms rms)); +reportOn (time \(mat.sum [cms, rms])); print "RMS + CMS... "; -reportOn (time \(mat.sum rms cms)); +reportOn (time \(mat.sum [rms, cms])); print "RMS + RMS... "; -reportOn (time \(mat.sum rms rms)); +reportOn (time \(mat.sum [rms, rms])); ();
--- a/src/may/matrix/test/test_matrix.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/matrix/test/test_matrix.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -227,17 +227,45 @@ compareMatrices (mat.mapColumns (vec.scaled 2) m') m'' ), +"minValue-\(name)": \( + compare (mat.minValue (fromRows [[1,2],[3,4],[5,-1]])) (-1) and + compare (mat.minValue (fromRows [[1,2],[3,0],[5,-1]])) (-1) and + compare (mat.minValue (fromRows [[1,2],[3,0],[5,1]])) 0 and + compare (mat.minValue (fromRows [[],[],[]])) 0 +), + +"maxValue-\(name)": \( + compare (mat.maxValue (fromRows [[1,2],[3,4],[5,-1]])) 5 and + compare (mat.maxValue (fromRows [[1,2],[3,0],[5,-1]])) 5 and + compare (mat.maxValue (fromRows [[-1,-2],[-3,0],[-5,-1]])) 0 and + compare (mat.maxValue (fromRows [[-1,-2],[-3,-4],[-5,-1]])) (-1) and + compare (mat.maxValue (fromRows [[],[],[]])) 0 +), + +"total-\(name)": \( + compare (mat.total (fromRows [[1,2],[3,4],[5,-1]])) 14 and + compare (mat.total (fromRows [[1,2],[3,0],[5,-1]])) 10 and + compare (mat.total (fromRows [[-1,-2],[-3,0],[-5,-1]])) (-12) and + compare (mat.total (fromRows [[-1,-2],[-3,-4],[-5,-1]])) (-16) and + compare (mat.total (fromRows [[],[],[]])) 0 +), + "sum-\(name)": \( compareMatrices - (mat.sum (constMatrix 2 { rows = 3, columns = 4 }) - (constMatrix 1 { rows = 3, columns = 4 })) - (constMatrix 3 { rows = 3, columns = 4 }) + (mat.sum [constMatrix 2 { rows = 3, columns = 4 }, + constMatrix 1 { rows = 3, columns = 4 }]) + (constMatrix 3 { rows = 3, columns = 4 }) and + compareMatrices + (mat.sum [constMatrix 2 { rows = 3, columns = 4 }, + constMatrix 1 { rows = 3, columns = 4 }, + fromRows [[1,2,3,4],[5,6,7,8],[9,10,11,12]]]) + (fromRows [[4,5,6,7],[8,9,10,11],[12,13,14,15]]) ), "sumFail-\(name)": \( try - \() (mat.sum (constMatrix 2 { rows = 3, columns = 4 }) - (constMatrix 1 { rows = 3, columns = 5 })); + \() (mat.sum [constMatrix 2 { rows = 3, columns = 4 }, + constMatrix 1 { rows = 3, columns = 5 }]); false; catch FailureException e: true @@ -257,11 +285,11 @@ { i = 1, j = 0, v = 5 }, { i = 1, j = 1, v = -4 }, // NB this means [1,1] -> 0, sparse zero ]); - tot = mat.sum s t; + tot = mat.sum [s, t]; mat.isSparse? tot and - compareMatrices tot (mat.sum (mat.toDense s) t) and - compareMatrices tot (mat.sum (mat.toDense s) (mat.toDense t)) and - compareMatrices tot (mat.sum s (mat.toDense t)) and + compareMatrices tot (mat.sum [mat.toDense s, t]) and + compareMatrices tot (mat.sum [mat.toDense s, mat.toDense t]) and + compareMatrices tot (mat.sum [s, mat.toDense t]) and compareMatrices tot (mat.newSparseMatrix { rows = 2, columns = 3 } (Rows [ { i = 0, j = 0, v = 1 }, @@ -332,6 +360,32 @@ (fromColumns [[58,139],[64,154]]) ), +"entryWiseProduct-\(name)": \( + compareMatrices + (mat.entryWiseProduct + [fromRows [[1,2,3],[4,5,0]], + fromRows [[6,7,8],[0,1,2]]]) + (fromRows [[6,14,24],[0,5,0]]) and + compareMatrices + (mat.entryWiseProduct + [fromRows [[1,2,3],[4,5,0]], + fromColumns [[6,0],[7,1],[8,2]]]) + (fromRows [[6,14,24],[0,5,0]]) +), + +"entryWiseDivide-\(name)": \( + compareMatrices + (mat.entryWiseDivide + (fromRows [[1,2,3],[4,5,0]]) + (fromRows [[6,7,8],[1,2,3]])) + (fromRows [[1/6,2/7,3/8],[4,5/2,0]]) and + compareMatrices + (mat.entryWiseDivide + (fromRows [[1,2,3],[4,5,0]]) + (fromColumns [[6,1],[7,2],[8,3]])) + (fromRows [[1/6,2/7,3/8],[4,5/2,0]]); +), + "sparseProduct-\(name)": \( s = mat.newSparseMatrix { rows = 2, columns = 3 } (Columns [ { i = 0, j = 0, v = 1 }, @@ -388,6 +442,33 @@ (mat.toSparse (fromColumns [[1,4],[2,5],[3,6]]))) ), +"tiledTo-\(name)": \( + compareMatrices + (mat.tiledTo { rows = 2, columns = 2 } + (fromColumns [[1,4],[2,5],[3,6]])) + (fromColumns [[1,4],[2,5]]) and + compareMatrices + (mat.tiledTo { rows = 3, columns = 4 } + (fromColumns [[1,4],[2,5],[3,6]])) + (fromColumns [[1,4,1],[2,5,2],[3,6,3],[1,4,1]]) and + compareMatrices + (mat.tiledTo { rows = 7, columns = 2 } + (fromColumns [[1,4],[2,5],[3,6]])) + (fromColumns [[1,4,1,4,1,4,1],[2,5,2,5,2,5,2]]) and + compareMatrices + (mat.tiledTo { rows = 1, columns = 1 } + (fromColumns [[1,4],[2,5],[3,6]])) + (fromRows [[1]]) and + compareMatrices + (mat.tiledTo { rows = 2, columns = 3 } + (mat.zeroMatrix { rows = 0, columns = 0 })) + (fromRows [[0,0,0],[0,0,0]]) and + compareMatrices + (mat.tiledTo { rows = 0, columns = 0 } + (mat.zeroMatrix { rows = 4, columns = 3 })) + (fromRows []) +), + "zeroSizeMatrix-\(name)": \( compareMatrices (mat.zeroMatrix { rows = 0, columns = 0 }) @@ -419,6 +500,36 @@ [[1,4],[0,5],[3,6]]; ), +"repeated-horiz-\(name)": \( + compareMatrices + (mat.repeatedHorizontal 2 (fromColumns [[1,4],[0,5]])) + (fromColumns [[1,4],[0,5],[1,4],[0,5]]) and + compareMatrices + (mat.repeatedHorizontal 2 (fromRows [[1,0],[4,5]])) + (fromColumns [[1,4],[0,5],[1,4],[0,5]]) and + compareMatrices + (mat.repeatedHorizontal 0 (fromColumns [[1,4],[0,5]])) + (fromColumns []) and + compareMatrices + (mat.repeatedHorizontal 4 (fromColumns [[]])) + (fromColumns [[],[],[],[]]) +), + +"repeated-vert-\(name)": \( + compareMatrices + (mat.repeatedVertical 2 (fromRows [[1,4],[0,5]])) + (fromRows [[1,4],[0,5],[1,4],[0,5]]) and + compareMatrices + (mat.repeatedVertical 2 (fromColumns [[1,0],[4,5]])) + (fromRows [[1,4],[0,5],[1,4],[0,5]]) and + compareMatrices + (mat.repeatedVertical 0 (fromRows [[1,4],[0,5]])) + (fromRows []) and + compareMatrices + (mat.repeatedVertical 4 (fromRows [[]])) + (fromRows [[],[],[],[]]) +), + "concat-horiz-\(name)": \( compareMatrices (mat.concatHorizontal
--- a/src/may/plot.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/plot.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -19,6 +19,8 @@ import org.jzy3d.plot3d.primitives.axes: ContourAxeBox; import org.jzy3d.contour: DefaultContourColoringPolicy; import org.jzy3d.contour: MapperContourPictureGenerator; +import org.jzy3d.chart.factories.AWTChartComponentFactory; +import org.jzy3d.plot3d.text.drawable: DrawableTextBillboard, DrawableTextBitmap; import javax.imageio: ImageIO; @@ -69,7 +71,7 @@ //!!! doesn't work if either rows or columns is 1 xrange = new Range(0, size.columns - 1); yrange = new Range(0, size.rows - 1); - grid = new OrthonormalGrid(xrange, size.columns, yrange, size.rows); +// grid = new OrthonormalGrid(xrange, size.columns, yrange, size.rows); eprintln "Matrix size: \(size)"; // surface = Builder#buildOrthonormal(grid, mapper); //??? big? // println "Z Bounds: \(surface#getBounds()#getZmin()) -> \(surface#getBounds()#getZmax())"; @@ -128,6 +130,11 @@ */ ()); +addCaption chart colour caption is ~Chart -> ~Color -> string -> () = + (scene = chart#getScene(); + text = new DrawableTextBitmap(caption, new Coord3d(0, 0, 0), colour); + scene#add(text)); + /** * Plot a list of typed structures onto a single chart, and return the * resulting Chart object. @@ -138,11 +145,14 @@ plot structures = (quality = Quality#Nicest; quality#setAnimated(false); - chart = new Chart(quality); + chart = AWTChartComponentFactory#chart(quality, "awt"); var j = 0; for structures do s: colour = distinctColour j; case s of + Caption caption: + addCaption chart (new Color(0, 0, 0)) caption; + j := j - 1; // doesn't count as a distinct structure Grid matrix: plotMatrix chart colour matrix; Contour matrix:
--- a/src/may/plot/chart.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/plot/chart.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -77,7 +77,8 @@ \() ChartLauncher#openStaticChart(chart); fi; if opts.saveTo != "" then - \() chart#screenshot(opts.saveTo); + outfile = new File(opts.saveTo); + \() chart#screenshot(outfile); fi); plotBarChart options values =
--- a/src/may/signal/window.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/signal/window.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -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/convolve.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/convolve.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -19,6 +19,7 @@ // framesize*2, in which each frame contains all channels one // after another (not interleaved) [because we lack a complex // matrix type] + //!!! we now have a complex matrix type, revisit this (forwardTransform = fft.realForward (framesize * 2); do stream: padded = @@ -46,29 +47,36 @@ else parts.fst :: splitInto n parts.snd; fi); -fastConvolvedWith ir framesize s = +fastConvolvedWith ir framesize = // prerequisite: ir and s have same number of channels (framer = zeroPaddedFreqFrames framesize (mat.height ir); - ch = s.channels; + ch = mat.height ir; irfr = framer (syn.precalculated 1 ir); // rate arg is irrelevant here - sigfr = framer s; inverseTransform = fft.realInverse (framesize * 2); - extended = sigfr ++ - map do _: list (cplx.zeros (ch * (framesize + 1))) done - [1..length irfr-1]; - cframes = doFastConvolve irfr extended; - rframes = (mat.zeroMatrix { rows = ch, columns = framesize * 2}) :: - map do fr: - mat.fromRows - (map inverseTransform (splitInto (framesize+1) fr)) - done cframes; - fr.streamed s.sampleRate (framesize * 2) - [ Window win.boxcar, Hop framesize ] - rframes; + do s: + if ch != s.channels then + failWith "Signal stream and IR must have same number of channels (\(s.channels) != \(ch))" + fi; + sigfr = framer s; + extended = sigfr ++ + map do _: list (cplx.zeros (ch * (framesize + 1))) done + [1..length irfr-1]; + cframes = doFastConvolve irfr extended; + rframes = + (mat.toRowMajor (mat.zeroMatrix { rows = ch, columns = framesize * 2})) :: + map do fr: + mat.fromRows + (map inverseTransform (splitInto (framesize+1) fr)) + done cframes; + fr.streamed s.sampleRate (framesize * 2) + [ Window win.boxcar, Hop framesize ] + rframes; + done; ); plainConvolvedWith ir s = // prerequisite: ir and s have same number of channels + //!!! test that and throw (var history = mat.toRowMajor (mat.zeroMatrix { rows = s.channels, columns = mat.width ir - 1 }); s with @@ -127,23 +135,21 @@ fi; nextPowerOfTwo' 1 n); -//!!! todo: partial application so as to preprocess ir (in fast convolution case) -convolvedWith options ir s = //!!! cheap mono thing here +//!!! doc note: this supports partial application (it returns a function that +// expects the signal) so as to pre-process the IR +convolvedWith options ir = //!!! cheap mono thing here //!!! doc note: fast convolution output is padded to next frame boundary - (if mat.height ir != s.channels then - failWith "Signal stream and IR must have same number of channels (\(s.channels) != \(mat.height ir))" - fi; - var type = Fast (); + (var type = Fast (); var framesize = nextPowerOfTwo (mat.width ir); for options \case of - Fast s: if s then type := Fast () else type := Plain () fi; + Fast f: if f then type := Fast () else type := Plain () fi; Framesize n: framesize := nextPowerOfTwo n; esac; case type of Fast (): - fastConvolvedWith ir framesize s; + fastConvolvedWith ir framesize; Plain (): - plainConvolvedWith ir s; + plainConvolvedWith ir; esac); {
--- a/src/may/stream/filter.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/filter.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -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 Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/framer.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -70,7 +70,8 @@ else mat.asRows (mat.concatHorizontal - [mat.zeroMatrix { rows = stream.channels, columns = hop }, + [mat.toRowMajor + (mat.zeroMatrix { rows = stream.channels, columns = hop }), stream.read (framesize - hop)]); fi; overlappingBlockList framesize hop stream framesize initialBuffer; @@ -85,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); @@ -105,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 @@ -122,7 +125,8 @@ (var remaining = frames; var position = 0; channels = mat.height (head frames); // so we don't need to keep a head ptr - var buffered = mat.zeroMatrix { rows = channels, columns = 0 }; + var buffered = mat.toRowMajor + (mat.zeroMatrix { rows = channels, columns = 0 }); { get position () = position, get channels () = channels, @@ -159,9 +163,10 @@ first::rest: (w = mat.width pending; pre = mat.columnSlice pending 0 (w - overlap); - added = mat.sum first - (mat.resizedTo (mat.size first) - (mat.columnSlice pending (w - overlap) w)); + added = mat.sum + [first, + (mat.resizedTo (mat.size first) + (mat.columnSlice pending (w - overlap) w))]; ola rest added (pre::acc)); _: reverse (pending::acc); @@ -170,7 +175,7 @@ first::rest: mat.concatHorizontal (ola rest first []); _: - mat.zeroMatrix { rows = 0, columns = 0 }; + mat.toRowMajor (mat.zeroMatrix { rows = 0, columns = 0 }); esac); streamOverlapping rate { framesize, hop, window } frames = @@ -181,7 +186,8 @@ w = vec.scaled factor (window framesize); channels = mat.height (head frames); // so we don't need to keep a head ptr - var buffered = mat.zeroMatrix { rows = channels, columns = 0 }; + var buffered = + mat.toRowMajor (mat.zeroMatrix { rows = channels, columns = 0 }); syncd = synchronized remaining; @@ -193,7 +199,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/manipulate.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/manipulate.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -41,7 +41,8 @@ if not s.finished? then mat.resizedTo { rows = s.channels, columns = n } (s.read n); else - mat.zeroMatrix { columns = n, rows = s.channels } + mat.toRowMajor + (mat.zeroMatrix { columns = n, rows = s.channels }); fi), }); @@ -49,7 +50,7 @@ (var prepos = 0; zeros n = mat.toRowMajor (prepos := prepos + n; - mat.zeroMatrix { rows = s.channels, columns = n }); + mat.toRowMajor (mat.zeroMatrix { rows = s.channels, columns = n })); delay = if nsamples < 0 then \0 (s.read (-nsamples)); @@ -105,7 +106,9 @@ pos := pos + zc; toAdd := toAdd - zc; mat.concatHorizontal - [r, mat.zeroMatrix { rows = s.channels, columns = zc }]; + [r, + mat.toRowMajor + (mat.zeroMatrix { rows = s.channels, columns = zc })]; fi), close = s.close }); @@ -149,7 +152,7 @@ if sz.columns == 0 then mat.zeroMatrix sz else - mat.sum (mat.resizedTo sz m1) (mat.resizedTo sz m2); + mat.sum [mat.resizedTo sz m1, mat.resizedTo sz m2]; fi); channels = head (sortBy (>) (map (.channels) streams)); { @@ -214,7 +217,8 @@ if s.available == Infinite () then s else var pos = 0; - var cache = mat.zeroMatrix { rows = s.channels, columns = 0 }; + var cache = mat.toRowMajor + (mat.zeroMatrix { rows = s.channels, columns = 0 }); chunks = array []; cachedPartsFor count = (start = pos % (mat.width cache); @@ -262,6 +266,7 @@ duplicated copies s = //!!! doc fact that original s cannot be used independently of this afterward // (so maybe name is misleading?) +//!!! also of course "duplicated" means two of them, while we could have lots array if copies < 2 then map \s [0..copies-1]; else pos = array (map \0 [0..copies-1]); @@ -376,7 +381,7 @@ read n = (m = s.read (n*frac); obtained = maths.ceil ((mat.width m) / frac); - mat.flipped + mat.toRowMajor (mat.fromColumns (map do c: mat.getColumn (c*frac) m done [0..obtained-1]))) };
--- a/src/may/stream/playback.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/playback.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -5,6 +5,7 @@ mat = load may.matrix; af = load may.stream.audiofile; ch = load may.stream.channels; +manip = load may.stream.manipulate; import javax.sound.sampled: AudioSystem, AudioFormat, AudioFormat$Encoding, SourceDataLine; @@ -45,18 +46,26 @@ line.close()); playStream stream = - (line = open { rate = stream.sampleRate, channels = stream.channels }; - blocksize = 10240; - not stream.finished? loop - line.play [(ch.mixedAndInterleavedTo line.channels - (stream.read blocksize))]; - line.close(); - stream.close()); + case stream.available of + Infinite (): + failWith "Refusing to play infinite stream. Use playStreamFor <sec> to limit playback duration"; + _: + line = open { rate = stream.sampleRate, channels = stream.channels }; + blocksize = 10240; + not stream.finished? loop + line.play [(ch.mixedAndInterleavedTo line.channels + (stream.read blocksize))]; + line.close(); + stream.close(); + esac; + +playStreamFor seconds stream = + playStream (manip.withDuration (seconds * stream.sampleRate) stream); playFile filename = playStream (af.open filename); { - open, playMatrix, playStream, playFile, + open, playMatrix, playStream, playStreamFor, playFile, }
--- a/src/may/stream/resample.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/resample.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -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 @@ -297,8 +297,9 @@ pd = phaseData[phase]; for [0..s.channels-1] do ch: +// println "lengths \(vec.length (mat.getRow ch buffer)) and \(vec.length pd.filter)"; 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 @@ -335,7 +336,10 @@ fi) }); -resampledTo = resampledDirectlyTo; +resampledTo targetRate s = + if s.sampleRate == targetRate then s + else resampledDirectlyTo targetRate s; + fi; { kaiserSincFilterFor,
--- a/src/may/stream/syntheticstream.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/syntheticstream.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -26,25 +26,22 @@ close = \(), }); -sinusoid rate freq = - generated rate (sin . (* (2 * pi * freq / rate))); +sinusoid sampleRate freq = + generated sampleRate (sin . (* (2 * pi * freq / sampleRate))); -whiteNoise rate = - generated rate \((Math#random() * 2.0) - 1.0); +whiteNoise sampleRate = + generated sampleRate \((Math#random() * 2.0) - 1.0); -pulseTrain rate freq = - generated rate do n: if n % int (rate / freq) == 0 then 1 else 0 fi done; +silent sampleRate = + generated sampleRate \0; -silent rate = - generated rate \0; - -precalculatedMono rate data = +precalculatedMono sampleRate data = (n = vec.length data; var position = 0; { get position () = position, get channels () = 1, - get sampleRate () = rate, + get sampleRate () = sampleRate, get available () = Known (n - position), get finished? () = not (n > position), read count = ch.deinterleaved 1 @@ -55,14 +52,14 @@ close = \(), }); -precalculated rate data = +precalculated sampleRate data = (n = mat.width data; c = mat.height data; var position = 0; { get position () = position, get channels () = c, - get sampleRate () = rate, + get sampleRate () = sampleRate, get available () = Known (n - position), get finished? () = not (n > position), read count = @@ -73,7 +70,7 @@ close = \(), }); -precalculatedRepeated rate data = +precalculatedRepeated sampleRate data = (n = mat.width data; c = mat.height data; chunks count pos = @@ -89,7 +86,7 @@ { get position () = position, get channels () = c, - get sampleRate () = rate, + get sampleRate () = sampleRate, get available () = Infinite (), get finished? () = false, read count = @@ -99,11 +96,11 @@ close = \(), }); -empty rate channels = // degenerate stream with no data in it, occasionally useful +empty sampleRate channels = // degenerate stream with no data in it, occasionally useful { get position () = 0, get channels () = channels, - get sampleRate () = rate, + get sampleRate () = sampleRate, get available () = Known 0, get finished? () = true, read count = mat.zeroMatrix { rows = channels, columns = 0 }, @@ -117,7 +114,6 @@ precalculatedRepeated is number -> mat.matrix_t -> stream_t, sinusoid is number -> number -> stream_t, whiteNoise is number -> stream_t, - pulseTrain is number -> number -> stream_t, silent is number -> stream_t, empty is number -> number -> stream_t, }
--- a/src/may/stream/test/test_framer.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/test/test_framer.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -5,6 +5,7 @@ vec = load may.vector; mat = load may.matrix; syn = load may.stream.syntheticstream; +win = load may.signal.window; { compare, compareUsing } = load may.test; @@ -12,6 +13,30 @@ testStream n is number -> 'a = syn.precalculatedMono rate (vec.fromList [1..n]); +expectedStream framesize opts n = + (var padded = true; + var hop = framesize; + for opts \case of + Padded p: padded := p; + Hop h: hop := h; + _: (); + esac; + if padded or (hop == framesize) then + testStream n + else + // Unpadded streams can't be perfectly reconstructed from + // windowed frames: the first framesize-hop samples will + // remain windowed (and scaled by the overlap-add) + m = framesize - hop; + factor = hop / (framesize/2); + w = vec.scaled factor (win.hann framesize); + syn.precalculatedMono rate + (vec.concat [ + vec.multiply [vec.fromList [1..m], vec.slice w 0 m], + vec.fromList [m+1..n] + ]) + fi); + compareFrames frames1 frames2 = all id (map2 do f1 f2: compareUsing mat.equal f1 f2 done frames1 (map (mat.newRowVector . vec.fromList) frames2)); @@ -19,7 +44,7 @@ testFramesInvertibleWith framesize opts length expected firstChunkSize = (f = fr.frames framesize opts (testStream length); str = fr.streamed rate framesize opts f; // inversion (framed -> streamed) - ts = testStream length; // newly initialised stream to compare with + ts = expectedStream framesize opts length; (firstChunk = str.read firstChunkSize; compareUsing mat.equal
--- a/src/may/stream/test/test_syntheticstream.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/test/test_syntheticstream.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -37,18 +37,6 @@ compare str.position 6 ), -"pulseTrain": \( - // 2 pulses a second sampled 8 times a second - str = syn.pulseTrain 8 2; - compare str.position 0 and - compare str.channels 1 and - compare str.sampleRate 8 and - compare str.available (Infinite ()) and - compare str.finished? false and - compareApprox epsilon (vec.list (mat.getRow 0 (str.read 9))) [ 1, 0, 0, 0, 1, 0, 0, 0, 1 ] and - compare str.position 9 -), - "silent": \( str = syn.silent 8; compare str.position 0 and
--- a/src/may/stream/test/test_waves.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/test/test_waves.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -5,8 +5,15 @@ mat = load may.matrix; vec = load may.vector; +plt = load may.plot; + { compare, compareUsing, compareMatrices, assert } = load may.test; +compareApprox eps = + compareUsing do v1 v2: all id (map2 do f1 f2: abs (f1 - f2) < eps done v1 v2) done; + +epsilon = 0.000001; + [ "square": \( @@ -21,6 +28,49 @@ 0,-0.942809])) ), +"impulseTrain-naive": \( + // 2 pulses a second sampled 8 times a second (naive pulse train) + str = wav.impulseTrain 8 2; + compare str.position 0 and + compare str.channels 1 and + compare str.sampleRate 8 and + compare str.available (Infinite ()) and + compare str.finished? false and + compareApprox epsilon (vec.list (mat.getRow 0 (str.read 9))) [ 1, 0, 0, 0, 1, 0, 0, 0, 1 ] and + compare str.position 9 +), + +"impulseTrain-naive-frac": \( + // 2.5 pulses a second sampled 10 times a second (still a naive + // pulse train even though the pulse rate is non-integer) + str = wav.impulseTrain 10 2.5; + compare str.position 0 and + compare str.channels 1 and + compare str.sampleRate 10 and + compare str.available (Infinite ()) and + compare str.finished? false and + compareApprox epsilon (vec.list (mat.getRow 0 (str.read 11))) [ 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0 ] and + compare str.position 11 +), + +"impulseTrain-bandlimited": \( + // 3 pulses a second sampled 10 times a second + str = wav.impulseTrain 100 3; + + x = mat.getRow 0 (str.read 200); +false; //!!! todo +/* + compare str.position 0 and + compare str.channels 1 and + compare str.sampleRate 8 and + compare str.available (Infinite ()) and + compare str.finished? false and + compareApprox epsilon (vec.list (mat.getRow 0 (str.read 14))) + [ 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0 ] and + compare str.position 14 +*/ +), + "cycleLengthFor": \( compare (wav.cycleLengthFor 48000 1000) 48 and compare (wav.cycleLengthFor 44100 1000) 441 and
--- a/src/may/stream/waves.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/stream/waves.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -4,60 +4,87 @@ module may.stream.waves; -{ sinusoid } = load may.stream.syntheticstream; +{ sinusoid, generated } = load may.stream.syntheticstream; { sum, scaledBy } = load may.stream.manipulate; { precalculatedRepeated } = load may.stream.syntheticstream; -{ gcd } = load may.mathmisc; +{ resampledTo } = load may.stream.resample; +{ gcd, floor, ceil } = load may.mathmisc; load may.stream.type; -cycleLengthFor rate freq = +cycleLengthFor sampleRate freq = // The number of samples we can cache and repeat from a wavetable, // in order to exactly reproduce a signal of the given frequency // at the given samplerate. - (if rate == int rate and freq == int freq and rate > 0 and freq > 0 then - rate / (gcd rate freq) + (if sampleRate == int sampleRate and + freq == int freq and + sampleRate > 0 and freq > 0 then + sampleRate / (gcd sampleRate freq) else - period = rate / freq; - m = find do i: i * period == int (i * period) done [1..rate]; + period = sampleRate / freq; + m = find do i: i * period == int (i * period) done [1..sampleRate]; if empty? m then 0 else period * (head m) fi fi); -square' rate freq = +square' sampleRate freq = sum (map do n: m = n*2 + 1; scaledBy (1/m) - (sinusoid rate (m * freq)) - done (0 :: [1 .. int (rate/4 / freq) - 1])); + (sinusoid sampleRate (m * freq)) + done (0 :: [1 .. int (sampleRate/4 / freq) - 1])); -saw' rate freq = +saw' sampleRate freq = sum (map do n: scaledBy if n % 2 == 0 then -1/n else 1/n fi - (sinusoid rate (n * freq)) - done [1 .. int (rate/2 / freq)]); + (sinusoid sampleRate (n * freq)) + done [1 .. int (sampleRate/2 / freq)]); -triangle' rate freq = +triangle' sampleRate freq = sum (map do n: m = n*2 + 1; scaledBy if n % 2 == 0 then -1/(m*m) else 1/(m*m) fi - (sinusoid rate (m * freq)) - done (0 :: [1 .. int (rate/4 / freq) - 1])); + (sinusoid sampleRate (m * freq)) + done (0 :: [1 .. int (sampleRate/4 / freq) - 1])); -cached f rate freq = - (n = cycleLengthFor rate freq; +// Bandlimited impulse train: produce naive pulse train at a frequency +// that is an integer fraction of the sample rate (so as to have exact +// sample locations), then resample so as to place the impulses at the +// requested frequency +impulseTrain sampleRate freq = + if freq > sampleRate/2 then + failWith "Can't generate impulse train with frequency > half the sample rate (\(freq) > \(sampleRate)/2)" + else + naiveTrain f = + (spacing = int (sampleRate / f); + generated sampleRate + do n: if n % spacing == 0 then 1 else 0 fi done); + ratio = sampleRate / freq; + if ratio == int ratio then + naiveTrain freq + else + integerRatio = floor ratio; + naiveFreq = sampleRate / integerRatio; + n = naiveTrain (sampleRate / integerRatio); + resampled = resampledTo (sampleRate * naiveFreq / freq) n; + resampled with { sampleRate }; + fi; + fi; + +cached f sampleRate freq = + (n = cycleLengthFor sampleRate freq; if n == 0 then - f rate freq; + f sampleRate freq; else // one could alternatively do - // repeated (precalculated rate ((f rate freq).read n)) + // repeated (precalculated sampleRate ((f sampleRate freq).read n)) // or - // repeated (withDuration n (f rate freq)) - precalculatedRepeated rate ((f rate freq).read n); + // repeated (withDuration n (f sampleRate freq)) + precalculatedRepeated sampleRate ((f sampleRate freq).read n); fi); square = cached square'; @@ -68,6 +95,7 @@ square is number -> number -> stream_t, saw is number -> number -> stream_t, triangle is number -> number -> stream_t, + impulseTrain is number -> number -> stream_t, cycleLengthFor is number -> number -> number, }
--- a/src/may/test.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/test.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -44,6 +44,16 @@ false; fi); +assertException f = + try + \() (f ()); + println "** failed to catch expected exception"; + false; + catch Exception _: + goodCompares := goodCompares + 1; + true; + yrt; + time msg f = (start = System#currentTimeMillis(); result = f (); @@ -62,8 +72,12 @@ println "Test \(name) failed"; name; fi - catch FailureException e: + catch Exception e: println "Test \(name) threw exception: \(e)"; + trace = e#getStackTrace(); + maxLen = 10; + for (take maxLen trace) do e: println " at \(e)" done; + if length trace > maxLen then println " ..." fi; name; yrt; done (sort (keys testHash))); @@ -81,7 +95,7 @@ bad); { - compare, compareUsing, compareMatrices, assert, + compare, compareUsing, compareMatrices, assert, assertException, time, runTests, }
--- a/src/may/vector.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/vector.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -18,7 +18,7 @@ import java.util: Arrays; import may.bits: VectorBits; -{ ceil } = load may.mathmisc; +{ floor, ceil } = load may.mathmisc; /// Return a vector of n zeros. zeros n is number -> ~double[] = @@ -33,6 +33,10 @@ /// Return a vector of length n, containing all ones. ones n = consts 1.0 n; +/// Return a vector of length n, containing random values. +randoms n is number -> ~double[] = + (map \(Math#random()) [1..n]) as ~double[]; + /// Return a vector of the values in the given list. fromList l is list?<number> -> ~double[] = l as ~double[]; @@ -204,24 +208,52 @@ len: sum' v / len esac; +sorted v is ~double[] -> ~double[] = + (v' = copyOf v; + Arrays#sort(v'); + v'); + +quantile q v is number -> ~double[] -> number = + if empty? v then 0 + else + v = sorted v; + dist = q * (length v - 1); + low = floor dist; + high = ceil dist; + if low == high then v[low] + else + v[low] * (high - dist) + v[high] * (dist - low); + fi; + fi; + +median v is ~double[] -> number = + quantile 0.5 v; + +listOp f vv = + case vv of + [v]: v; + v::rest: + (out = copyOf v; + for rest (f out); + out); + _: failWith "Empty argument list"; + esac; + +//!!! doc: throws exception if not all inputs are of the same length 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: throws exception if not all inputs are of the same length 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: throws exception if not all inputs are of the same length +multiply vv is list?<~double[]> -> ~double[] = + listOp do out v: VectorBits#multiplyBy(out, v) done vv; + +//!!! doc: throws exception if not all inputs are of the same length +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 @@ -236,7 +268,9 @@ fi; sqr v = - multiply v v; + (out = copyOf v; + VectorBits#multiplyBy(out, v); + out); rms = sqrt . mean . sqr; @@ -293,6 +327,7 @@ zeros, consts, ones, + randoms, vector v = v, primitive, floats, @@ -308,15 +343,19 @@ slice, resizedTo, reverse, + sorted, repeated, concat = concat', sum = sum', mean, + median, + quantile, add, subtract, multiply, + divide, + scaled, divideBy, - scaled, abs = abs', negative, sqr, @@ -333,6 +372,7 @@ zeros is number -> vector_t, consts is number -> number -> vector_t, ones is number -> vector_t, + randoms is number -> vector_t, vector is ~double[] -> vector_t, primitive is vector_t -> ~double[], floats is vector_t -> ~float[], @@ -348,15 +388,19 @@ slice is vector_t -> number -> number -> vector_t, resizedTo is number -> vector_t -> vector_t, reverse is vector_t -> vector_t, + sorted is vector_t -> vector_t, repeated is vector_t -> number -> vector_t, concat is list?<vector_t> -> vector_t, sum is vector_t -> number, mean is vector_t -> number, + median is vector_t -> number, + quantile is number -> 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, - scaled is number -> vector_t -> vector_t, abs is vector_t -> vector_t, negative is vector_t -> vector_t, sqr is vector_t -> vector_t,
--- a/src/may/vector/test/test_vector.yeti Mon Mar 03 17:03:32 2014 +0000 +++ b/src/may/vector/test/test_vector.yeti Sat Apr 05 16:02:24 2014 +0100 @@ -3,7 +3,7 @@ vec = load may.vector; -{ compare } = load may.test; +{ compare, assertException } = load may.test; [ @@ -163,24 +163,65 @@ compare ((vec.mean . vec.zeros) 0) 0 and compare ((vec.mean . vec.zeros) 5) 0 and compare ((vec.mean . vec.ones) 5) 1 and - compare ((vec.mean . vec.fromList) [1,-2,3,0]) 0.5 + compare ((vec.mean . vec.fromList) [1,-2,3,0]) 0.5 and + compare ((vec.mean . vec.fromList) [1,-2,3,2]) 1 +), + +"median": \( + compare ((vec.median . vec.zeros) 0) 0 and + compare ((vec.median . vec.zeros) 5) 0 and + compare ((vec.median . vec.ones) 5) 1 and + compare ((vec.median . vec.fromList) [-2,1,2,3]) 1.5 and + compare ((vec.median . vec.fromList) [1,-2,3,2]) 1.5 and + compare ((vec.median . vec.fromList) [1,-2,3,0]) 0.5 and + compare ((vec.median . vec.fromList) [1,-2,3,0,6]) 1 +), + +"upperQuartile": \( + compare ((vec.quantile 0.75 . vec.zeros) 0) 0 and + compare ((vec.quantile 0.75 . vec.zeros) 5) 0 and + compare ((vec.quantile 0.75 . vec.ones) 5) 1 and + compare ((vec.quantile 0.75 . vec.fromList) [0,1,2]) 1.5 and + compare ((vec.quantile 0.75 . vec.fromList) [0,1,2,3]) 2.25 and + compare ((vec.quantile 0.75 . vec.fromList) [1,-2,3,2]) 2.25 +), + +"sorted": \( + compare ((vec.list . vec.sorted . vec.zeros) 0) [] and + compare ((vec.list . vec.sorted . vec.zeros) 5) [0,0,0,0,0] and + compare ((vec.list . vec.sorted . vec.ones) 5) [1,1,1,1,1] and + compare ((vec.list . vec.sorted . vec.fromList) [-2,1,2,3]) [-2,1,2,3] and + compare ((vec.list . vec.sorted . vec.fromList) [1,-2,3,2]) [-2,1,2,3] and + compare ((vec.list . vec.sorted . vec.fromList) [1,-2,-2.1,3,0.5,1]) [-2.1,-2,0.5,1,1,3] ), "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.zeros 0, vec.ones 0])) [] and + compare (vec.list (vec.add [vec.consts 3 4, vec.fromList [1,2,3,0] ])) [4,5,6,3] and + assertException \(vec.add [vec.consts 3 3, vec.fromList [1,2,3], vec.fromList [6,7,8,9] ]) ), "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.zeros 0) (vec.ones 0))) [] and + compare (vec.list (vec.subtract (vec.consts 3 4) (vec.fromList [1,2,3,0]))) [2,1,0,3] and + assertException \(vec.subtract (vec.consts (-3) 4) (vec.fromList [1,2,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 0])) [] and + compare (vec.list (vec.multiply [vec.consts (-3) 4, vec.fromList [1,2,3,0]])) [-3,-6,-9,0] and + assertException \(vec.multiply [vec.consts 3 3, vec.fromList [1,2,3], vec.fromList [6,7,8,9]]) +), + +"divide": \( + compare (vec.list (vec.divide (vec.zeros 0) (vec.ones 0))) [] and + compare (vec.list (vec.divide (vec.consts (-3) 3) (vec.fromList [1,2,3]))) [-3,-(3/2),-1] and + assertException \(vec.divide (vec.consts (-3) 4) (vec.fromList [1,2,3])); +), + +"scaled": \( + compare (vec.list (vec.scaled 5 (vec.ones 0))) [] and + compare (vec.list (vec.scaled 5 (vec.fromList [1,2,-3]))) [5,10,-15] ), "divideBy": \( @@ -188,11 +229,6 @@ compare (vec.list (vec.divideBy 5 (vec.fromList [1,2,-3]))) [0.2,0.4,-0.6] ), -"scaled": \( - compare (vec.list (vec.scaled 5 (vec.ones 0))) [] and - compare (vec.list (vec.scaled 5 (vec.fromList [1,2,-3]))) [5,10,-15] -), - "abs": \( compare (vec.list (vec.abs (vec.ones 0))) [] and compare (vec.list (vec.abs (vec.fromList [1,2,-3]))) [1,2,3]