view src/may/vector/blockfuncs.yeti @ 454:a926ad176efb

Add negative

module may.vector.blockfuncs;

vec = load may.vector;

load may.vector.type;

import may.bits: VectorBits;

// "internal" vector function to retrieve data for read-only
// purposes without copying
raw =
   (raw' v is ~double[] -> ~double[] = v;
    raw' as vector -> ~double[]);

sum' v = VectorBits#sum(raw v);

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

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

min' v = 
   (dat = raw v;
    var mn = 0;
    for [0..length dat - 1] do i:
        if i == 0 or dat[i] < mn then
            mn := dat[i];
        fi
    done;
    mn);

minindex v =
   (dat = raw v;
    var mn = 0;
    var mi = -1;
    for [0..length dat - 1] do i:
        if i == 0 or dat[i] < mn then
            mn := dat[i];
            mi := i;
        fi
    done;
    mi);

mean v =
    case vec.length v of
        0: 0;
        len: sum' v / len
    esac;

add bb =
   (len = head (sort (map vec.length bb));
    vv = map raw bb;
    out = new double[len];
    for vv do v:
        VectorBits#addTo(out, v, len);
    done;
    vec.vector out);

subtract b1 b2 =
   (v1 = raw b1;
    v2 = raw b2;
    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;
    vec.vector out);

multiply b1 b2 = VectorBits#multiply(raw b1, raw b2);

scaled n v =
    if n == 1 then v
    else VectorBits#scaled(raw v, n);
    fi;

divideBy n v =
    // 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 vec.fromList (map (/ n) (vec.list v));
    fi;

sqr v =
    multiply v v;

rms =
    sqrt . mean . sqr;

abs' =
    vec.fromList . (map abs) . vec.list;

negative =
    vec.fromList . (map (0-)) . vec.list;

sqrt' =
    vec.fromList . (map sqrt) . vec.list;

unityNormalised v = 
   (m = max' (abs' v);
    if m != 0 then
        divideBy m v;
    else
        v;
    fi);

fftshift v =
   (len = vec.length v;
    half = int(len/2 + 0.5); // round up for odd-length sequences
    vec.concat [vec.slice v half len, vec.slice v 0 half]);

ifftshift v =
   (len = vec.length v;
    half = int(len/2); // round down for odd-length sequences
    vec.concat [vec.slice v half len, vec.slice v 0 half]);

{
sum is vector -> number = sum',
mean is vector -> number,
add is list?<vector> -> vector,
subtract is vector -> vector -> vector,
multiply is vector -> vector -> vector, 
divideBy is number -> vector -> vector, 
scaled is number -> vector -> vector,
abs is vector -> vector = abs',
negative is vector -> vector,
sqr is vector -> vector,
sqrt is vector -> vector = sqrt',
rms is vector -> number,
max is vector -> number = max',
min is vector -> number = min',
maxindex is vector -> number,
minindex is vector -> number,
unityNormalised is vector -> vector,
fftshift is vector -> vector,
ifftshift is vector -> vector,
}