May

/**

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

{ floor, 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 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[];

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

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[] =

    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[] =

    listOp do out v: VectorBits#subtractFrom(out, v) done [v1, v2];

//!!! 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

    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 =

   (out = copyOf v;

    VectorBits#multiplyBy(out, v);

    out);

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

euclideanDistance v1 v2 is ~double[] -> ~double[] -> number =

    VectorBits#euclideanDistance(v1, v2);

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,

    randoms,

    vector v = v,

    primitive,

    floats,

    fromFloats,

    fromList,

    list = list',

    array = array',

    length = length',

    empty? = empty?',

    at = at',

    equal,

    equalUnder,

    slice,

    resizedTo,

    reverse,

    sorted,

    repeated,

    concat = concat',

    sum = sum',

    mean,

    median,

    quantile,

    add,

    subtract,

    multiply,

    divide,

    scaled,

    divideBy,

    abs = abs',

    negative,

    sqr,

    sqrt = sqrt',

    rms,

    max = max',

    min = min',

    maxindex,

    minindex,

    unityNormalised,

    euclideanDistance,

    zipped,

    unzipped,

} as {

    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[],

    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,

    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 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, 

    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,

    euclideanDistance is vector_t -> vector_t -> number,

    zipped is list?<vector_t> -> vector_t,

    unzipped is number -> vector_t -> array<vector_t>,

}