view src/may/vector.yeti @ 486:061dc568cfd8

Minor simplifications in vector

/**
 * Vectors. A May vector is a typesafe, immutable wrapper around a Java
 * primitive array of doubles.
 *
 * Although not as convenient and flexible as a Yeti array<number> or
 * list<number>, a vector can be faster and more compact when dealing
 * with dense data of suitable range and precision such as sampled
 * sequences.
 *
 * Because vectors are immutable, functions that transform them will
 * return the original data for identity transformations
 * (e.g. resizing to the vector's existing length).
 */

module may.vector;

import java.util: Arrays;
import may.bits: VectorBits;

{ ceil } = load may.mathmisc;

/// Return a vector of n zeros.
zeros n is number -> ~double[] =
    new double[n];

/// Return a vector of length n, containing all m.
consts m n is number -> number -> ~double[] =
   (a = zeros n;
    Arrays#fill(a, m);
    a);

/// Return a vector of length n, containing all ones.
ones n = consts 1.0 n;

/// Return a vector of the values in the given list.
fromList l is list?<number> -> ~double[] =
    l as ~double[];

/// Return the given vector as a list.
list' a is ~double[] -> list<number> =
    list a;

/// Return the given vector as a Yeti array.
array' a is ~double[] -> array<number> =
    array a;

/// Return the length of the given vector.
length' v is ~double[] -> number =
    length v;

/// Return true if the given vector is empty (has length 0).
empty?' =
    empty? . list';

/// Return element n in the given vector v. (The function name and
/// argument order are chosen for symmetry with the similar standard
/// library array function.)
at' v n is ~double[] -> number -> number =
    v[n];

/// Return the given vector as a Java primitive float array.
floats a is ~double[] -> ~float[] =
   (len = length' a;
    f = new float[len];
    for [0..len-1] do i:
        f[i] := a[i];
    done;
    f);

/// Return a vector of the values in the given Java primitive float array.
fromFloats ff is ~float[] -> ~double[] =
   (len = length (list ff);
    a = new double[len];
    for [0..len-1] do i:
        a[i] := ff[i];
    done;
    a);

/// Return true if the given vectors are equal, using the standard ==
/// comparator on their elements.
equal v1 v2 is ~double[] -> ~double[] -> boolean =
    list v1 == list v2;

/// Return true if the given vectors are equal, when applying the
/// given numerical comparator to each element.
equalUnder comparator v1 v2 =
    length' v1 == length' v2 and
        all id (map2 comparator (list v1) (list v2));

/// Return another copy of the given vector.
copyOf v is ~double[] -> ~double[] =
    Arrays#copyOf(v, length v);

/// Return the given vector as a primitive array. Modifying the
/// array contents will not affect the original vector.
primitive = copyOf;

/// Return a new vector containing a subset of the elements of the
/// given vector, from index start (inclusive) to index finish
/// (exclusive). (The function name and argument order are chosen for
/// symmetry with the standard library slice and strSlice functions.)
slice v start finish is ~double[] -> number -> number -> ~double[] =
   (len = length v;
    if start == 0 and finish == len then v
    elif start < 0 then slice v 0 finish
    elif start > len then slice v len finish
    elif finish < start then slice v start start
    elif finish > len then slice v start len
    else
        Arrays#copyOfRange(v, start, finish);
    fi);

/// Return a vector of length n, containing the contents of the given
/// vector v. If v is longer than n, the contents will be truncated;
/// if shorter, they will be padded with zeros.
resizedTo n v is number -> ~double[] -> ~double[] =
    if n == length v then v;
    else Arrays#copyOf(v, n);
    fi;

/// Return a vector that is the reverse of the given vector.  Name
/// chosen (in preference to passive "reversed") for symmetry with the
/// standard library list reverse function.
reverse v is ~double[] -> ~double[] =
   (len = length v;
    a = new double[len];
    for [0..len-1] do i:
        a[len-i-1] := v[i];
    done;
    a);

/// Return a single new vector that contains the contents of all the
/// given vectors, in order. (Unlike the standard module list concat
/// function, this one cannot be lazy.)
concat' vv is list?<~double[]> -> ~double[] =
    case vv of
        [v]: v;
        v0::rest:
           (len = sum (map length' vv);
            vout = Arrays#copyOf(v0 is ~double[], len);
            var base = length' v0;
            for rest do v: 
                vlen = length' v;
                System#arraycopy(v, 0, vout, base, vlen);
                base := base + vlen;
            done;
            vout);
        _: zeros 0;
    esac;

/// Return a single new vector that contains the contents of the given
/// vector, repeated n times. The vector will therefore have length n
/// times the length of v.
repeated v n is ~double[] -> number -> ~double[] =
    concat' (map \(v) [1..n]);

sum' v is ~double[] -> number = 
    VectorBits#sum(v);

max' v is ~double[] -> number = 
   (var mx = 0;
    for [0..length v - 1] do i:
        if i == 0 or v[i] > mx then
            mx := v[i];
        fi
    done;
    mx);

maxindex v is ~double[] -> number =
   (var mx = 0;
    var mi = -1;
    for [0..length v - 1] do i:
        if i == 0 or v[i] > mx then
            mx := v[i];
            mi := i;
        fi
    done;
    mi);

min' v is ~double[] -> number = 
   (var mn = 0;
    for [0..length v - 1] do i:
        if i == 0 or v[i] < mn then
            mn := v[i];
        fi
    done;
    mn);

minindex v is ~double[] -> number =
   (var mn = 0;
    var mi = -1;
    for [0..length v - 1] do i:
        if i == 0 or v[i] < mn then
            mn := v[i];
            mi := i;
        fi
    done;
    mi);

mean v is ~double[] -> number =
    case length v of
        0: 0;
        len: sum' v / len
    esac;

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);

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);

multiply b1 b2 is ~double[] -> ~double[] -> ~double[] = 
    VectorBits#multiply(b1, b2);

scaled n v is number -> ~double[] -> ~double[] =
    if n == 1 then v
    else VectorBits#scaled(v, n);
    fi;

divideBy n v is number -> ~double[] -> ~double[] =
    // Not just "scaled (1/n)" -- this way we get exact rationals. In fact
    // the unit test for this function will fail if we use scaled (1/n)
    if n == 1 then v
    else fromList (map (/ n) (list v));
    fi;

sqr v =
    multiply v v;

rms =
    sqrt . mean . sqr;

abs' v is ~double[] -> ~double[] =
    fromList (map abs (list v));

negative v is ~double[] -> ~double[] =
    fromList (map (0-) (list v));

sqrt' v is ~double[] -> ~double[] =
    fromList (map sqrt (list v));

unityNormalised v is ~double[] -> ~double[] = 
   (m = max' (abs' v);
    if m != 0 then
        divideBy m v;
    else
        v;
    fi);

zipped vv is list?<~double[]> -> ~double[] =
    case vv of
    [v]: v;
    first::rest:
       (len = length' first;
        if len != min' (fromList (map length' vv)) then
            failWith "Vectors must all be of the same length";
        fi;
        fromList
           (concat
              (map do i:
                   map do v:
                       at' v i
                       done vv
                   done [0..len-1])));
     _: zeros 0;
    esac;

unzipped n v is number -> ~double[] -> array<~double[]> =
    if n == 1 then array [v]
    else 
        len = length' v;
        vv = array (map \(new double[ceil(len / n)]) [1..n]);
        for [0..len-1] do x:
            vv[x % n][int (x / n)] := at' v x;
        done;
        array vv;
    fi;

typedef opaque vector_t = ~double[];

{
    zeros,
    consts,
    ones,
    vector v = v,
    primitive,
    floats,
    fromFloats,
    fromList,
    list = list',
    array = array',
    length = length',
    empty? = empty?',
    at = at',
    equal,
    equalUnder,
    slice,
    resizedTo,
    reverse,
    repeated,
    concat = concat',
    sum = sum',
    mean,
    add,
    subtract,
    multiply,
    divideBy,
    scaled,
    abs = abs',
    negative,
    sqr,
    sqrt = sqrt',
    rms,
    max = max',
    min = min',
    maxindex,
    minindex,
    unityNormalised,
    zipped,
    unzipped,
} as {
    zeros is number -> vector_t,
    consts is number -> number -> vector_t,
    ones is number -> vector_t,
    vector is ~double[] -> vector_t,
    primitive is vector_t -> ~double[],
    floats is vector_t -> ~float[],
    fromFloats is ~float[] -> vector_t,
    fromList is list?<number> -> vector_t,
    list is vector_t -> list<number>,
    array is vector_t -> array<number>,
    length is vector_t -> number,
    empty? is vector_t -> boolean,
    at is vector_t -> number -> number,
    equal is vector_t -> vector_t -> boolean,
    equalUnder is (number -> number -> boolean) -> vector_t -> vector_t -> boolean,
    slice is vector_t -> number -> number -> vector_t,
    resizedTo is number -> vector_t -> vector_t,
    reverse 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,
    add is list?<vector_t> -> vector_t,
    subtract is vector_t -> vector_t -> vector_t,
    multiply is vector_t -> 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,
    sqrt is vector_t -> vector_t,
    rms is vector_t -> number,
    max is vector_t -> number,
    min is vector_t -> number,
    maxindex is vector_t -> number,
    minindex is vector_t -> number,
    unityNormalised is vector_t -> vector_t,
    zipped is list?<vector_t> -> vector_t,
    unzipped is number -> vector_t -> array<vector_t>,
}