USAGE:

    #loadlib "vector.so"

Contents:

Routines for sorting indexing and handling sorted vector elements:

qsort(afRes,afX):

-- quicksort vector afX[]. Return iDim, or -1 on error. Vector afRes[] is arranged in ascending order. Requires |afX|=|afRes|; NaNs are put at low end. return value=size;

countSortLevel(pfSorted):

-- count number of distinct levels in sorted array. For unsorted array return num of not-equal consecutive lines. Tip: call after: nnum=qsort(pfS,pfX); nlev=countlevels(pfS+dimof(pfX)-nnum); NaN are ignored anyway;

rank(pfRanks, pfX):

-- returns rank table pfRanks[] of unsorted pfX[]; return num-of-NaN-entries; pfRanks=pfX is ok; // do not confuse with index(),midranks(),ranks();

index(afIndex, afX):

-- compute index array from data vector. Vector afIndex[] index arragement in ascending order. Requires |afIndex|=|afRes|; Index-stable (sequence is secondary sort criteria); Nan-save: Return value is number of elements NaN with index=0.

midranks(pfMidranks, pfXsorted)

-- returns ranks values pfRanks[] of pfXsorted[]. Ranks range [1..|pfX|], ties return mid-ranks values. Return sum_all ties(n^3-n) for n=tie number of elements with same value.

ranks(pfRanks, pfXsorted)

-- obsolete name of midranks()

percentindex(pfSorted, percent):

-- returns interpolated value into sorted array psSorted[] at given percent level.

findindex(pfIndex, pfX, pfBook, highIdx):

-- find and return closest indeces pfIndex[i] of each element pfX[i] in sorted map pfBook[?]. The found index [0,|pfBook|-1] is rescaled to [0,highIdx]. E.g. highIdx = 100 (percent), 255 (for 8bit images), dim(pfBook)-1 (index same as 0). Return number of elements outside [min,max] of pfBook[] and -1 on error. Requires |pfIndex|=|pfX|.

locateindex(fX,pfBook):

-- locate index j in sorted map given an array pfBook[0..n-1], and given a value fX, returns a value j such that fX is between pfBook[j] and pfBook[j+1]. pfBook[] must be either increasing or descreasing. j=-1 or j=n-1 is returned to indicate that fX is out of range.

findkNN(pfIndex, pfD, pfX, iDim, iIdx, iMode):

-- find the indices of the kNN of the element with index iIdx in pfX.

Arguments:

     pfIndex: vector to hold indices of nearest neighbours
               if dimof(pfIndex)==0, only distances are computed.
               dimof(pfIndex) determines the nr. of neighbours to search
     pfD: vector to hold distances from el. iIdx to all others
               |pfD| must be iN = |pfX|/iDim.
               If iIdx == -1, it must contain a reference vector that
               is the basis for kNN search (see iIdx below)
     pfX: dataset with column dimension iDim.
     iDim: nr. of columns in data set pfX
     iIdx: index of element to ref. position for kNN search.
               if iIdx = -1, the first iDim els of pfD
               are used as reference vector
     iMode: Mode to specify norm for distance evaluation
                0: abs norm
                1: city block norm
                2: euclidiean
                3: minkowski-norm order iMode
     
   The nr. of neighbours seeked is given by dimof(pfIndex)
   Sorting is done by insertion sort, which is efficient for small k.
   If MANY neighbours are needed, better use qsort
    

remap(pfRes, pfSrc, neglBeg, neglEnd):

-- remap vector of length |pfSrc| by interpolation to length |pfRes|. Furthermore a number of elements can be neglected at the beginning (neglBeg>=0) and end (neglEnd>=0) of the original vector pfSrc[]. Return the zoom factor. *

Example image enhancement:

pfDest[] may be significantly smaller than pfSrc in order to accelerate findindex() on further processing.

    static float imgSorted[NX], table[200]; // NX large image, e.g. 512^2
    qsort(imgSorted, X); // find pixel distribution
    remap(table, imgSorted, 0.05*NX, 0); // compact map, trash 5% low-pix
    findindex(Y, X, table, 255); // smoothen pixel distribution,
                                             // enhance low contrast parts

remap2d(pfD,iDWid,pfS,iSWid,iBlock,iMode)

-- remap 2d-vector-picture. A 2d-vector-picture is a rectangular area of vectors, organized in memory as a sequence of all vector-components for all pixels in one row and all rows in a sequence. The structure is determined by the width, height and vector-dimension (which is given by iBlock). Examples are complex-pictures, interleaved RGB-pictures, STFT-Spectra, phase vocoder data... remap2d rescales the 2d-vector-picture in pfS, Width iSWid, height iSHgt = |pfS| / (iBlock*iDWid) to pfD, width iDWid, height iDHgt = |pfD| / (iBlock*iDWid). Both pfS and pfD have subvector dimension iBlock. Operation is depending on Mode:

    iMode & 1 = 0: SAMPLE, no interpolation
    iMode & 1 = 1: AVERAGEING
    iMode & 2 = 0: STRIDE: only process first subvector-entry
    iMode & 2 = 1: BLOCK: process all subvector-entries

    RMAP2_INTERPOLATE: interpolate for upsample / average for downsample
    RMAP2_BLOCK: process each element of subvector, else only first component

remap2dscale(pfD,iDWid,pfS,iSWid,fScaleX,fScaleY,iBlock, iMode)

-- remap2dscale rescales a 2d-vector-picture with subvector dimension iBlock. in pfS, Width iSWid, height iSHgt = |pfS| / (iBlock*iDWid) to pfD, width iDWid, height iDHgt = |pfD| / (iBlock*iDWid). The difference to remap2d is that the X/Y-scaling-factors can/must be specified explicitely. If iDWid<iSWid*fScaleX, the rest is clipped If iDWid>iSWid*fScaleX the rest is set to black. Same clipping happens to the y-dimension. iMode has the same meaning as in remap2d.

subrect( pfD,iDWid,iDXo,iDYo, pfS,iSWid,iSXo,iSYo, iSXb,iSYb,iBlock,iMode)

-- copy rectangular area. subrect copys a rectangular subarea from pfS (Width: iSWid) to pfD (Width: iDWid). The height of the vector-pictures is determined by iDHgt = |pfD|/(iBlock*iDWid) (for S resp.) (...each pixel consists of an iBlock dimensional subvector.). The subrect is specified by its LL-Corner coordinates in S: (iSXo, iSYo) and its width iSXb and height iSYb. A clipping is done if subrect is not inside S or D. Operation is depending on iMode:

   Bit 0: -
   Bit 1: 0: STRIDE (process only first subvector-entry)
   Bit 1: 1: BLOCK (process all subvector-entries)
   

checkorder(pfX):

-- check ordering of vector, return bit mask:

     * bit 0 (&1): descending order occured,
     * bit 1 (&2): constant order occured, and
     * bit 2 (&4): ascending occured

Routines for creating random sequences without doubles or drop outs:

permuteindex(afIndex):

-- compute random index sequence afIndex[] with index between 0 and iD-1. Use as initialisation with nextindex().

nextindex(afIndex):

-- return next index in array and mark used if all items are marked used (=-1) call automatically permuteindex() for next random sequence through the index set 0,1,2,... iD-1.

Routines for histogramming:

histo(pfCount, pfMidpoints, pfX):

-- compute histogram of pfX[] elements. Determine linear bin structure by min and max value of pfX[]. Returns bin counts in pfCount[] and bin mid-points in array pdMidpoints[]. Requires |pfCount|=|pfMidpoints| (calls binning()). Return number of valid data or -1 on error. NaN-save;

binning(pfCount, pfMidpoints, pfX):

-- bin and count vector elements pfX[] in bins defined by mid-points in array pfMidpoints[]. Result in pfCount[]. pfMidpoints[] must be in ascending order; Requires |pfCount=|pfMidpoints|; Return number of valid data or -1 on error. NaN-save;

histoline(pfLine, pfCounts, pfMidpoints):

-- compute polyline for histogram. Given counts in pfCounts[] and ordered bin-midpoints in pfMidpoints[] generate vector containing 4 vertices per bin to draw polyline (using draw_sym). |pfLine| = 8*|pfCounts| =8*|pfXmid|. Returns -1 on error.

entropy(pfCount)

-- compute the entropy S of a distribution in n bins (e.g. result of binning() or histo()). S = -sum(p_i*log_2(p_i)), p_i is the probability p_i=pfCount[i]/sum(pfCount). S ranges between 0 (one single peak) and ln_2(|pfCount|) equals uniform distribution.

Routines for statistic analysis:

quantile(pfRes, pfX):

-- compute required quantile levels. Return |pfRes| uniform spaced quantiles between 0 and 100% (incl.) in pfRes[]. Levels are rounded and interpolated if half between two values (=> median definition); Performs sorting if |pfRes|!=2. Return value =size(pfX);

quantile_i(pfRes, pfX):

-- compute required quantile levels by interpolation. same as quantile() but does interploation at all intermediate leveles.

quantileSet(pfQ, pfX, mode):

-- compute a set of quantile levels of data vector pfX[]. mode bit 0 chooses standard quantile levels (s.below). mode bit 1 turns print option to stdout on/off. If mode out of range 0..3, the standard set level is read out (no computation). Non-standard level set are defined by presetting pfQ[], each between 0 and 1. For standard sets (mode & 1=1) dim(pfQ) is restricted to

    dim(pfQ)=5: min, 25%, 50%, 75%, max
    dim(pfQ)=7: 0%, 5%, 25%, 50%, 75%, 95%, 100%,
    dim(pfQ)=9: 0%, 2%, 5%, 25%, 50%, 75%, 95%, 98% 100%,
    dim(pfQ)=11: 0%, 1%, 2.5%, 5%, 25%, 50%, 75%, 95%, 97.5%, 99%, 100%,
    dim(pfQ)=13: 0%, 1%, 2.5%, 5%, 10%, 25%, 50%, 75%, 90%, 95%, 97.5%, 99%, 100%,
    (round qualtile to next full data entry, or interpolate only at half way between to levels)

mean(pfX):

-- return mean value of vector pfX[]; (see also .stat4)

adev(pfX):

-- return average deviation from mean value of vector pfX[]

var(pfX):

-- return variance of vector pfX[]

std(pfX):

-- return standard deviation value of vector pfX[]

stat4(pfRes, pfX):

-- return mean, std, min, and max of vector pfX[]. This is much faster than calling the corresponding above routines in sequence. NaN save. Returns number of valid elements.

    pfRes[0]: minimum, if |pfRes|>= 1,
    pfRes[1]: maximum, if |pfRes|>= 2,
    pfRes[2]: mean, if |pfRes|>= 3,
    pfRes[3]: standard deviation, if |pfRes|>= 4,
    pfRes[4]: max-min, if |pfRes|>= 5,
    pfRes[5]: max(abs(X)), if |pfRes|>= 6,

dsstat(pfMM, iDimM, pfD):

-- return stat4 results for all columns in dataset pfD. Call stat4 for all column vectors in pfD, collect stat4-results in pfMM. Set iDimM to nr. of columns in pfMM.

moments(pfRes, pfX):

-- return several moments and quantiles levels of vector pfX[], see below. NaN save. Returns number of valid elements. If pfRes[0] is set <0 before call the result will be printed to stdout.

    pfRes[0]: average deviation
    pfRes[1]: mean
    pfRes[2]: standard deviation
    pfRes[3]: skewness
    pfRes[4]: kurtosis
    pfRes[5]: variance
    pfRes[6..n+6]: quantiles in n levels pfRes[] are filled as provided

linreg(pfRes, pfX, pfY)

-- compute linear regression given data points pfX[] and pfY[] find straight line y=ax+b: Requires number of data pnts MIN(dim(pfY),dim(pfX))>=2 and |pfRes|>=2. Returns the correlation coefficient r (Paerson's r) and following parameters:

     pfRes[0]=a = slope,
     pfRes[1]=b = y(x=0),
     pfRes[2]=r (Paerson's r=return value), if |pfRes|>5,
     pfRes[3]=Sum(xy)/(n-1) cross corretion, if |pfRes|>5,
     pfRes[4]=Sum(xx)/(n-1)=Var(x) = Variance x, if |pfRes|>5,
     pfRes[5]=Sum(yy)/(n-1)=Var(y) = Variance y, if |pfRes|>5,
     pfRes[6]=sig_a = probable uncertainty of a, if |pfRes|>6,
     pfRes[7]=sig_b = probable uncertainty of b, if |pfRes|>7,
     pfRes[8]=chi2 = sum of squared regression error, if |pfRes|>8,
     pfRes[9]=r_ab = correlation coefficient(a,b), if |pfRes|>9,
     pfRes[10]=mean(x), if |pfRes|>11,
     pfRes[11]=mean(y), if |pfRes|>11,
     pfRes[12]=mean(x)-sqrt(Var(x)), if |pfRes|>15,
     pfRes[13]=a*pfRes[12]+b, if |pfRes|>15,
     pfRes[14]=mean(x)+sqrt(Var(x)), if |pfRes|>15,
     pfRes[15]=a*pfRes[14]+b, if |pfRes|>=16.

     pfRes[16]=mean(x)-2*sqrt(Var(x)), if |pfRes|>=20,
     pfRes[17..19] same as pfRes[13] but two time std
     ... max 5 blocks of drawable 2-xy-points blocks

linmedreg(pfRes, pfX, pfY)

-- compute robust linear regression y=ax+b by median for given data points pf[] and pfY[]. Superior to linreg() if data set contains some bad outliers. Requires number of data pnts MIN(dim(pfY),dim(pfX))>=2. and |pfRes|>=2 Minimizes the mean absolute deviation MAD = sum(|pfY[i]-a*pfX[i]-b|)/dim(pfX). Return value of MAD and following parameters:

     pfRes[0]=a = slope,
     pfRes[1]=b = y(x=0),
     pfRes[2]=mean(x), if |pfRes|>=4,
     pfRes[3]=mean(y), if |pfRes|>=4,
     pfRes[4]=mean(x)-sqrt(Var(x)), if |pfRes|>=8,
     pfRes[5]=a*pfRes[4]+b, if |pfRes|>=8,
     pfRes[6]=mean(x)+sqrt(Var(x)), if |pfRes|>=8,
     pfRes[7]=a*pfRes[6]+b, if |pfRes|>=8,
     .. [8-11] repeat [4-7] but 2*std(x)..
   

lgamma(x)

-- return the log of the gamma function

gamma(x) = integralint_0^infty(t^x-1\exp(-t))dt;

The factorial n! = gamma(n+1);

gammq(a,x)

-- return the incomplete gamma function for a>0:

Q(a,x) = 1-P(a,x) = gamma(a,x)/gamma(a)

Q(a,x) = 1/gamma(a) * integralint_x^infty(t^a-1\exp(-t))dt

betai(a,b,x)

-- return the incomplete beta function with .a,b>0 and .0<=x<=1

I_x(a,b)=B_x(a,b)/B(a,b)=1/B(a,b)*integralint_0^x(t^a-1(1-t)^b-1)dt

limits are I_0=(a,b)=0 and I_1(a,b)=1 , it raises sharply at about x=a/(a+b)

student_t(t,dof)

-- return student distribution for t with dof degrees of freedom

      pfRet[0] returns uncertainty coefficient of attribte 0, knowing attribute 1
      pfRet[1] returns uncertainty coefficient of attribte 1, knowing attribute 0
      pfRet[2] returns CramersV
      pfRet[3] returns Chi-square statistics
      pfRet[4] returns significance level for chi2 statistics
      pfRet[5] returns number of degrees of freedom
      pfRet[6] return value: entropy of attribute 0 distribution
      pfRet[7] return value: entropy of attribute 1 distribution
      pfRet[8] return value: total entropy
      pfRet[9] return value: iB1
      pfRet[10] return value: iB2
      dimof(pfRet) is minimum 3

covmat(pfC, pfD, pfS, pfWeights, iMode)

-- compute data covariance matrix (pfC, |pfC|=iDim^2) from a subset of elements of pfD (|pfD|=iDim*iN) defined by an index vector pfS (|pfS| = 0 selects all elements), using weights pfWeight[i] for the ith column. iDim = |pfWeights| is assumed. pfWeight may be used to mute certain columns. iMode can by set 1 to substract the sample mean automatically. In this case the mean is returned in pfWeights. (this is done indirectly using C(X-M) = C(X) - C(M) which is only valid if M is the mean of X)

meanvec(pfC, pfD, pfS)

-- compute mean data vector (center of mass) to pfC (iDim = |pfC|) of a subset of rows of data set pfD (|pfD| = iDim*iN). The selection is specified by index vector pfS (if |pfS| = 0 --> take all data)

kerneldensity(pfD,fX,fH,iMode)

-- estimate probability density at fX of a one dimensional data distribution (given by sample vector pfD) using kerneldensity estimation with kernel width fH. Return value prob(x) = 1./iN * sum_i ( K_fH ( fX - pfD[i]) ), where K_fH(.) = 1./fH * K(x/fH) is the fH scaling of kernel K(.), which is selected by iMode. iMode==0 -> triangular kernel K(x) = (1-|x|) if |x|<1. else 0 iMode==1 -> gaussian kernel K(x) = sqrt(1/(2*Pi)) * exp(-0.5*x^2)

kernelde(pfD,pfX,fH,iMode)

-- compute probability density at pfX by kernel density estimation, using the data set pfD, (iDim = |pfX|). Return value is prob(x) = 1/iN * sum_i=0^iN-1 K_fH( x - x_i ). where iN is the nr. of records in the data set, K_fH(x) = 1./fH^d * K(x / fH) is the fH-scaling of the Kernel K(.) d=iDim is the data space dimension and K(.) is the kernel, selected by iMode. iMode==0 -> Triangular kernel K(x) = 1.- |x| if |x|<1., else 0. iMode==1 -> Gaussian kernel K(x) = 1/sqrt(2*Pi)^d exp(-0.5*(x'x)) Use 100+iMode to compute also the gradient at pfX. Return Value is then the density, the gradient is returned in pfX. (currently only 101 implemented)

kernelreg(pfR, pfD, pfS, fSigma, iMode)

-- compute kernel regression for data set pfD. pfS specifies independent and dependent columns by pfS = [iNI, IC_0, ., IC_(iNI-1), iNO, OC_0, ., OC_(iNO-1)]. This means that the iNI cols with indices [IC_0,..,IC_(iNI-1)] are input variables and the iNO cols with indices [OC_0,..,OC_(iNO-1)] are output variables. Example: pfS = [2, 0, 4, 1, 3] -> input cols 0, 4, output column 3. input and output columns should be disjunct (is not tested). Given input cols in pfR (in corresponding columns) kernelreg() computes pfR[j] = Norm*sum_(i=0..iN-1) pfD[j] * exp(-0.5*(dist(i)/fSigma)^2), where j is from (OC_0..OC_(iNO-1)), dist(i) is |pfR-row(pfD,i,iN)|_(input space) and Norm=1/(sum(i=0..iN-1) exp(-0.5*(dist(i)/fSigma)^2)). Kernel regression results are stored into the resp. OC cols in pfR.

Routines for setting a vector:

zero(pfX)

-- set vector to zero. 'DEPRECATED - use pfX:=0

set(pfX,x)

-- set vector to value x. 'DEPRECATED - use pfX:=x

linspace(pfX, x1, x2)

-- return vector pfX[n] of n equally spaced points between x1 and x2 inclusive.

logspace(pfX, e1, e2)

-- return vector of n logarithmically equally spaced points between 10**e1 and 10**e2 inclusive.

gspace(pfX, x1, x2)

-- geometric sequence. Return vector of n logarithmically equally spaced points between x1 and x2 inclusive. I.e. an exponential curve between x1 and x2. (Note difference to logspace() which takes exponents base 10). Return factor pfX[1]/pfX[0] or -1 on error (x1*x2<=0).

rand(pfX,low,high)

-- return random vector with uniform distribution between low and high value.

randn(pfX,std)

-- return random vector with normal distribution, zero mean and standard deviation std.

chirpwave(pfY, li, lf, offs)

-- Return |pfY| dim vector given by pfY[i]=sin(2*pi*i*(1/li+(1/lf-1/li)*i/(|pfY|-1))+offs) which is a sinewave with linear frequency drift, starting with initial wavelength li, ending with final lf, and phase offset offs. Requires |pfY|>1, li>0, and lf>0.

chirpwave_linWL(pfY, li, lf, offs)

-- same as chirpwave() but with linear wavelength shift given by pfY[i]=sin(2*pi*(i+offs)/(li+(lf-li)*i/(|pfY|-1)))

mk_nkernel(pfY,type,width,length)

-- make kernel of length length and type with typical width width. Requires |pfY|>=length. On negative width generate reverse ordered kernel. On negative type generate un-normalized kernel, positive kernel types are normalized to sum(pfY)=1 (which produces on tiny or zero width kernels with one peak element =1). Return value is the sum(pfY) or -1 on error (unknown type). Defined kernel (absolute) shapes |types| are

     KERNEL_00 (i<=width)? 1:0 // flat square kernel of width
     KERNEL_01 exp(-(i/width) ) // exp decay = rolling average
     KERNEL_02 exp(-(i/width)^2) // half Gauss
     KERNEL_11 (i/width) * exp(-(i/width) )
     KERNEL_12 (i/width) * exp(-(i/width)^2)
     KERNEL_21 (i/width)^2 * exp(-(i/width) )
     KERNEL_22 (i/width)^2 * exp(-(i/width)^2)
     KERNEL_31 (i/width)^3 * exp(-(i/width) )
     KERNEL_32 (i/width)^3 * exp(-(i/width)^2)
     KERNEL_33 (i/width)^3 * exp(-(i/width)^3)
     KERNEL_GAUSS
     KERNEL_02m exp(-(i-(|pfY|-1)/2)^2) // mid Gauss bell
     KERNEL_SAW1 max(1-(i/width)), 0) // left edge
     KERNEL_SAW2 max(1-(i/width)^2, 0) // left hill
     KERNEL_SAW3 max(1-(i/width)^3, 0) // left hill
     KERNEL_SAW4 max(1-(i/width)^4, 0) // left hill
     KERNEL_SAW1m max(1-|i-((|pfY|-1)/2)/width|), 0) // mid symm
     KERNEL_SAW2m max(1-(i-((|pfY|-1)/2)/width)^2), 0) // mid hill
     KERNEL_SAW3m max(1-|i-((|pfY|-1)/2)/width)^3|, 0) // mid hill
     KERNEL_SAW4m max(1-(i-((|pfY|-1)/2)/width)^4), 0) // mid hill
     KERNEL_DIRAC
     KERNEL_E kronecker_delta( i, floor(width)) // =base vector
     KERNEL_SINC sin(x)/x, x=i/(width/pi) // sinc left, 1st zero at i=width
     KERNEL_SINCm sin(x)/x, x=(i-(|pfY|-1)/2)/(width/pi) // sinc mid

mk_kernel(pfY,type,width)

-- make kernel with length |pfY|. Same as mk_nkernel(pfY,type,width,size(pfY)).

Setter und getter to half of symmetric square matrix:

Unitary vector routines:

norm0(pfX):

-- return maximum norm of vector pfX[] (=max(|x_i|)). NaN-save.

norm1(pfX):

-- return absolute norm of vector pfX[] (=sum(|x_i|)). NaN-save.

norm2(pfX):

-- return Euclidean norm of vector pfX[] (== norm(pfX)). NaN-save.

norm(pfX):

-- same as norm2(pfX). NaN-save.

sum(pfX):

-- return sum of all elements of vector pfX[]. NaN-save.

prod(pfX):

-- return product of all elements of vector pfX[]. NaN-save.

min(pfX):

-- return minimal element of vector pfX[]. NaN-save.

max(pfX):

-- return maximal element of vector pfX[]. NaN-save.

mini(pfX):

-- return index of minimal element of vector pfX[]. NaN-save.

maxi(pfX):

-- return index of maximal element of vector pfX[]. NaN-save.

Vector operations:

vecvec(float *s,float *x,int n=0,int gap=0,int *best=0,float *dist=0,float *metric=0)

seek for bestmatch position of s in x and return index. Let d=dimof(s). Then, s is compared with d-long blocks of x, starting at positions x+i*(d+gap) and the i that yields the best position is returned. If n>0, only values i<n are considered. If best is non-NULL, the m=dimof(best) best positions are determined in ascending order (i.e., the best position coming first) and returned in best[]. If dist is also non-NULL and of the same dimension as best, dist[] additionally will hold the euclidean distance values (in the same ascending order) for the m positions. If metric is non-NULL and has at least d elements, the squared distance contribution from element at position k in each d-long block will be weighted by a factor metric[k] (which should be nonnegative). To try all consecutive starting positions, use gap = -d+1, which is the lowest allowed value.

sprod(pfX1, pfX2):

-- return scalar (inner) product of vectors pfX1[] and pfX2[].

mprod(pfRik, pfMij, pfMjk, j)

-- return (ixk) matrix product pfRik[] of the (ixj) matrix pfMij[] and (jxk) matrix pfMjk[] all in vector storage (rows=vert x columns=horiz, columns change fastest). j is given, i and k determined from the given vector dimensions. Returns i or -1 on error on dimension mismatch. E.g.: mprod(r24,a23,b34,3) returns the matrix product in matrix r24[]. mprod(r33,a3,b3,1) returns the outer product of two vectors a3*b3'. mprod(s1,a3,b3,3) returns the scalar product of two vectors in s1[0]=a3'*b3.

eprod(pfDest,pfX1,pfX2)

-- return elementwise product in vector pfDest[i]=pfX1[i]*pfX2[i] and return value norm2(). pfDest, pfX1, and pfX2 may point to the same memory. NaN-save.

multadd(pfY, a, pfX, b):

-- return vector Y=a*X+b. Elementwise multiply-add, a,b are scalars. useful for scaling, translation of vectors, e.g. d=multadd(Y,-1,X2,6); multadd(Y,1,X,0); multadd(X,2,X,0)

multaddv(pfY, pfA, pfX, pfB):

-- return vector Y = A *X + B. Elementwise multiply-add, |pfA|,|pfX|, and |pfB| must be either dimension 1 or |pfX|. Useful for copy, addition, substraction, scaling, translation of vectors, e.g.

     c1=1; c0=0; fac=4;
     d=multaddv(Y,&c1,X2,X1); // c1*X2+X1
     multaddv(Y,Y,X,&c0); // same as eprod(Y,Y,X)
     multadd(Y,1,X,-2); // elementwise -2
     multadd(Y,&c0,&c0,&c1); // same as set(Y,1);

add(pfDest, pfX1, pfX2):

-- return addition vector .pfDest[i] = pfX1[i] + pfX2[i] and return its norm. pfDest, pfX1, and pfX2A may point to the same memory. NOT-nan-save, see add_N().

add_N(pfDest, pfX1, pfX2):

-- return addition vector pfDest[i] = pfX1[i] + pfX2[i]. NaN-save (single entries ignores, double kept). -.diff(pfDest, pfX1, pfX2) -- return subtraction vector pfDest[i] = pfX1[i] - pfX2[i] and return its norm. pfDest, pfX1, and pfX2A may point to the same memory. NaN-save. dimof(pfDest)=0 is accepted and leads to no result vector computation, only norm of diff vector is computed.

scale(pfDest, fScale, pfX):

-- returns scaled vector pfDest[i] = fScale * pfX[i] and its norm. pfDest and pfX may point to the same memory. NaN-save.

clip(pfDest, pfX, lowClipValue, highClipValue):

-- returns vector pfDest[i] cliped by minimum lowClipValue and manimum highClipValue. Return number of cliped values. pfDest and pfX may point to the same memory. NaN-values unaffected.

threshold(pfRes,pfX,threshold):

-- threshold operation returns vector with 0 or 1s depending whether element is larger or equal threshold value.

reverse(pfRes, pfX)

-- reverses order of vector elements in array pfX. Reverseing in place with pfRes=pfX is allowed, then return number of swap-operations else return |pfX|. Requires |pfRes|=|pfX|pfX).

transpose(pfNM, pfMN, n)

-- transpose matrix pfMN into pfNM previously with n columns (= n rows afterwards). Return m or -1 on error

      static float m32[6]= {11,12,13, 21,22,23};
      transpose(Y,m32,2); // returns 3
                           // Y={ 11,21, 12,22, 13,23 }

switch_input(pfX, pfVectors, select)

-- switch one (select-ed) of a set of vectors pfVectors[] and copy to output vector pfX[]. Requires |pfVectors| is a multiple of |pfX|. select modulo |pfVectors|/|pfX| choose the vector postion. Return value !=0 indicates wrapping of select.

switch_output(pfVectors, pfX, select)

-- copy pfX to one (select-ed) of a set of vectors pfVectors[].

switch(pfDest, pfSrc, select)

-- switch_input() or switch_output() depending on dimension. pfDest is the desitination.

distmatrix(pfDistMat,pfData,pfWeight,fExponent)

-- compute distance matrix pfDistMat[iRows*iRows] from data table pfData[iRows*iDim] with column weights pfWeights[iDim]. The element i*Row+j is computed as the Minkowski Norm of data iDim-dimensional row i and j. pfDistMat[i*iRow+j]=( Sum_k=0..iDim-1 |pfWeight[k]* (pfData[i*iDim+k]-pfData[j*iDim+k])| ^fExponent )^(1/fExponent). Return the sum of distances between pairs of records; -.distmatrix2(pfDistMat,pfData) :-- the fast version of distmatrix with Euclidean norm (fExponent=2) and unit weighting (all pfWeight[]=1).

distmatrixToEquivalenceTable(pfTab,pfDistMat)

-- compute equivalence table pfTab[iRows] from distance matrix pfDistMat[iRows*iRows] (previous result from e.g. distmatrix()). pfTab[i] is for a distance record 0, for groups of record with zero distance the smallest index i carries the count of "slave" records as gets referenced by them as "master" by shifted negative index -(i+1). return the total number of dependent records (sum of slaves). This function is useful for treating distance driven iterations.

elog(pfDest,fX)

-- returns elementwise logarithm .pfDest[i]=log(pfX[i]). pfDest and pfX may point to the same memory. Return number of elements or -1 on error (negative elements).

elog10(pfDest, fX)

-- returns elementwise logarithm base 10, otherwise same as elog().

eexp(pfDest,pfX)

-- return elementwise exponential value in vector pfDest[i]=exp(pfX[i]). pfDest and pfXmay point to the same memory.

efermi(pfDest,pfX)

-- returns elementwise fermi function pfDest[i]=1/(1+exp(-pfX[i])). pfDest and pfX may point to the same memory; Return number of elements or -1 on error (m==0, dim mismatch).

eabs(pfDest,pfX)

-- return elementwise absolute value in vector pfDest[i]=|pfX[i]| pfDest and pfXmay point to the same memory.

emod(pfDest,pfX,m)

-- returns elementwise modulo pfDest[i]=pfX[i]%m. pfDest and pfX may point to the same memory; Return number of elements or -1 on error (m==0, dim mismatch).

eequal(pfX,f)

-- count elements in pfX[] equal to f

convolve(pfRes, pfX, pfKernel):

-- return convolution pfRes[] of pfKernel[] and pfX[]: pfRes[i]=sum_j(pfX[j+i]*pfKernel[j]). Return number of computed values iN = min(|pfRes|, |pfX|-|pfKernel|+1).

     static float D1Op[2]={ -1, 1}; // 1.derivative
     static float D2Op[3]={ -1, 2, -1}; // 2.derivative
     static float SmoothOp[3]={0.2, 0.6, 0.2};// smoothing
     convolve(Y,X,D1Op);

convolve_0(pfRes, pfX, pfKernel, kernelRefPos)

-- return convolution pfRes[] of pfKernel[] and pfX[]. Same as convolve() but uses zero-padding at edge positions. kernelRefPos gives the target position of the convolution operation, i.e. kernelRefPos=0 is left aligned, kernelRefPos=|pfKernel|/2+1 is mid alligned (odd kernel size) Requires 0<=kernelRefPos<|pfKernel| and |pfRef|=|pfX|. Returns |pfX| or -1 on error.

convolve_e(pfRes, pfX, pfKernel, kernelRefPos)

-- convolution, same as convolve_0(), but with edge extension (missing values take edge value).

convolve_c(pfRes, pfX, pfKernel, kernelRefPos)

-- convolution, same as convolve_0(), but with cyclic extension (missing values take wrapped values).

convolve_(pfRes, pfX, pfKernel, kernelRefPos, mode)

-- switched convolution. If mode=1 same as convolve_e(), on mode=2 convole_c(), else convolve_0().

sparse_copy(pfRes,offRes,incRes, pfX,offX,incX, maxCopy)

-- sparse copy. returns pfRes[] filled with maxCopy elements of pfx[] starting at the offset position offRes/offX with increment incRes/incX. Return value is the number of copy operations. Copying stops beyond vector boundaries. If incRes and/or incX is negative the corresponing offset is taken negative startig at the last element of the vector. A zero or negative maxCopy value is ignored, (incRes=0 && maxCopy>1) is considered as an error. A NaN save version is scopy_N() and scopy_2N().

scopy(pfRes,incRes, pfX,incX, maxCopy)

-- like sparse_copy() but offset arguments dropped since prog_unit allows offsets by calling

    scopy(Y+3, 1, X+3,2, 3);
    scopy(Y,-1, X,1, 0); // same as reverse(Y,X), if |X|=|Y|

sparse_multadd(pfX,stride, a,b, maxNum)

-- sparse multiply add operating in place on max maxNum sparse elements (position: i=0,1*stride,2*stride,...) of vector pfX[]. pfX[i]=a*pfX[i]+b; Operation stops at vector boundaries. If stride is negative the corresponing offset is taken negative starting at the last element of the vector. A zero or negative maxNum value is ignored. Return value is the number of mult-add operations.

sparse_multaddv(pfX,xStride, pfA,aStride, pfB,bStride, maxNum)

-- sparse multiply add operating on max maxNum sparse elements (position: i=0,1,2,...) of vector pfX[] with vector scale and offset. pfX[i*xStride]=pfA[i*aStride] * pfX[i*xStride] + pfB[i*bStride]; Operation stops at any vector boundaries. If a stride is negative the corresponing offset is taken negative starting at the last element of the vector. A zero or negative maxNum value is ignored, Return value is the number of mult-add operations.

sparse_sum(pfX,stride, maxNum)

-- sparse summation on max maxNum sparse elements (position: i=0,1,2,...maxNum-1) of vector pfX[]. Return sum(pfX[i*stride]). Operation stops at vector boundaries. A zero or negative maxNum value is ignored. -.sparse_ssum(pfX,stride, maxNum) sparse square summation on max maxNum sparse elements (position: i=0,1,2,...maxNum-1) of vector pfX[]. Return sum((pfX[i*stride])^2).

Routines for shifting vectors and time series representations:

push(pfX, x):

-- shift vector pfX[] right, return last, insert x at pos 0.

pop(pfX, x):

-- return pos 0, shift vector pfX[] left, insert x at end.

acov(pfACV, pfX):

-- compute data auto covariance function. Return sum_j=i^iD pfX[j]*pfX[j-i] vector of length dim(pfACV)

acor(pfACV, pfX):

-- compute data auto corelation function. Return acov[]/acov[0] normalized.

Routines for complex vector operation:

cabs(pfX,pfC)

-- compute absolute values of complex vector pfC[] stored as Re/Im pairs. Requires |pfX|=|pfC|/2. Return total Norm;

cangle(pfPhi,pfC)

-- compute angle of complex vector pfC[] Re/Im pairs. Requires |pfPhi|=|pfC|/2. All returned components are between -M_PI and M_PI. Return average angle, or -10 on error.

c2polar(pfP,pfC)

-- convert complex Re/Im pairs to abs/angle pairs. Requires |pfP|=|pfC|=muliple of 2. Inplace operation with pfP=pfC is ok. Return norm2 of elements or -1 on error.

polar2c(pfC,pfP)

-- convert abs/angle pairs to complex Re/Im pairs. Requires |pfP|=|pfC|=muliple of 2. Inplace operation with pfP=pfC is ok. Return num of element pairs or -1 on error.

csquare(pfX,pfC)

-- compute square absolute values of complex Re/Im pairs. Requires |pfX|=|pfC|/2. Return total Norm^2;

cprod(pfRes,pfC1,pfC2)

-- compute elementwise complex product of complex Re/Im pairs. |pfRes|,|pfC1|,|pfC2| must be muliple of 2. If dimMatchIsStrict=1 |pfC1|,|pfC2| = |pfRes|,2 (see below). If dimMatchIsStrict=0 perform as many operation as possible. Return num of operations or -1 on error.

Routines for handling NaN:

setNaN(pfX,iL)

-- set iL elements to NaN. iL<=0 is ignored; Return value is NaN itself.

firstNaN(pfX)

-- return first index of element equal NaN. If all are valid return NaN. Use e.g. by: if(firstNaN(X)>=0)...

cntNaN(pfX)

-- return number of elements equal NaN.

size(pfX)

-- returns effective dimension of pfX[] equals dimof(pfX)-cntNaN(pfX).

swapNaN(pfX)

-- swap all NaN-elements of pfX to beginning of vector pfX. Return number of element equal to NaN.

replaceNaN(pfX,y)

-- replace all NaN elements in pfX[] by y Return number of affected elements.

cpNonNaN(pfD, pfS)

-- copy as many as possible elements not-NaN. Return number of copied elements.

scopy_N(pfRes,incRes, pfX,incX, maxCopy)

-- like scopy() but skip NaN; Return number of copy ops. Allows offsets by calling scopy(Y+3, 1, X+3,2, 3); A zero or negative maxCopy value is ignored, (incRes=0 && maxCopy>1) is considered as an error.

scopy_2N(pfR1,pfR2,incR, pfX1,pfX2,incX, maxCopy)

-- similar to scopy_N() but copy pairs of elements if both are valid numbers (skipping any NaN); The stride value incR is for both destination vectors pfR1[]/pfR2[] and incX for both source vectors pfX1[]/pfX2[]. Return number of copy operations. Terminated if any of the four cursors run out of bounds.

add_N(pfDest, pfX1, pfX2)

NaN-save (replace by 0) -- add vector, NaN save version. Replace NaN by 0, return norm(pfDest) and addition vector pfDest[i] = pfX1[i] + pfX2[i]. pfDest, pfX1, and pfX2A may *NOT* point to the same memory, on NaN-elements

Constructor:

extract(int stride, char *srcmask)

creates instance to extract elements from data partitioned into rows of stride elements each (an element can be a vector element or a string token). Srcmask is a mask string that provides a row-wise applied "extract pattern" whose constituents determine how the value at the corresponding position in the row is used:

#

ignore that value

-

same as #

X

extract that value and append it to the result row

[a:b]

find for that value the nearest neigbor n from the range a,a+1,..b and append to the result row the n th unit vector (of dim b-a+1) ("nominal expansion"). Example: [1:3]: value 2.1 will cause a tripel (0 1 0) in the result row.

[n]

if abs(n)>1: same as [1:abs(n)]

[1]

same as 0:1

[0]

map negative values into (1,0), nonnegative ones into (0,1).

[-1]

same as [0].

f

repeat the constituent that follows f times

*

wildcard: represents a suitable number of consecutive #' so that the total number of constituents amounts to d

Methods:

float *elements(float *src, int rows=-1, int first=0)

int *elements(int *src, int rows=-1, int first=0)

byte *elements(byte *src, int rows=-1, int first=0)

return newly allocated vector containing the concatenation of all extracted elements.

float *elements(float *src, int *rowmask)

int *elements(int *src, int *rowmask)

byte *elements(byte *src, int *rowmask)

analogous, but now the subset of rows from which elements will be extracted is selected by the non-zero entries in array rowmask.

char *tokens(char *src, sep=0)

analogous, but now src is considered as a "vector" whose elements are tokens, with stride tokens always forming one row. Return value is a newly allocated string consisting of the concatenation (with separators in between) of the extracted tokens. If sep=0, a blank is used as separator. Currently, there is no way of selecting a subset of rows for this method. Also, there is no nominal expansion, i.e., [n] entries in src are treated as # don't cares.

Routines with special data structure:

copylist(pfY,pfX,pfList,iMax):

-- copy .iMax patterns from pfX[] to pfY[] as defined by list structure, which is encoded in one of the formats below. Default: len=l, srcStride=destStride=len, destOffset=0, d0=0, d1=1...; s? source index in [0..(srcStride-1)]; d? destination index in [-destOffset..(destStride-1-destOffset)]; srcStride/destStride is pattern length to/from (if zero take len); destOffest is an offset for dest indecies; If iMax<=0, copy as many patterns as possible; If iMax negative, warn if dimof(pfX/pfY) is not a multiple of srcStride/destStride. Return number of copied patterns or negative on error.

    List-FORMATS of pfList[] with l=dimof(pfList)
    F1: s0, s1, s2, ... s(l-1)
    F2: -len, srcStride, s0, s1, s2 ... s(len-1),
    F3: -len, srcStride, s0, s1, s2 ... s(len-1), 0, 0, ...
    F4: -len, srcStride, s0, s1 .. s(len-1), destOffset, destStride
    F5: -len, srcStride, s0, s1 .. s(len-1), destOffset, destStride,
                          d0, d1, d2 ... d(len-1), ...

token_cat(pftDest,pftAdd)

-- add a null-terminated sequence of elements (string) to string pftDest[]. Assure separation with global separator. Return null-termination position (or last element), or negative required length, if pftDest ist too short.

token_idx(pfIdx, pfString, fSep):

-- return index set pfIdx[] similar that repeated call to strtoken(). (I.e. indexing element sequences unequal to fSep, evtl stopping at the first zero element if fSep!=0). Note, pfString is not string, but vector formated (no conversion needed) Return value is the number of tokens n found and negative value -n if dimof(pfIdx)<n too small (resize and call again). string-index FORMAT: i0, i1, i2, ... i(n), -1 ... Separators value fSep is stored, see also token_set().

token_set(fSep):

-- set separator, return old. The global separator is used for multiple token copy operations.

token_fetch1(pfRet, pfString):

-- return first token from pfString[] in pfRet[] (0-terminated); (Note that pfSting can be called with offset.) Return token len, or negative on error, i.e. -1=bad=index, else negative required length if out-of-memory (resize and call again). Use global separator, set by token_idx() or token_set().

token_getidx(pfRet, pfString, pfIdx, iElement ):

-- return iElement token from string pfString starting at pfIdx[iElement] (note the difference to token_get(,,iIdx)), or default-string-generation, if token not available. Return token length or negative required length if pfRet is out-of-memory (resize and call again). The reason for double dereference is the ease of usage; missing values, do not result in any out-of-bounds error and need no external treatment.

token_copylist(tDest,tSrc,pfList)

-- copy tokens according to copy-list, see copylist() (only F1/F2-FORMAT) return number of elements copied, or negative required length, or -1 on error. Use separator as set by token_idx() or token_set().

token_setDefaultLabel(iFormat,iOffset)

-- set string format of default text token and additive offset (default=(0,0)):

     0=%i, 1=#%i, 2=row_%i, 3=col_%i, 4=axis_%i, 12=%02i, 13=%03i ... 15=%05
   

copylist_from_extract_mask(pfList,tFmt,iStride):

-- compute copylist from extraction mask string and source stride length iStride. Return copylist length and copylist pfList[] in F1,F2, or F3 FORMAT. If pfList[] is too short, the return value is the negative required length (=-dim for F1, resize here to len+2). The usage of FORMAT F1 is discouraged, since the source stride is not encoded. Therefore allow dimof(pfList)>=2 + returnValue. Return 0 on zero or negative iStride (prompt error) (be careful do not create endless resize-while-loops here) In the format string, use X for to-be-extracted positions, # for to-be-ignored positions. Prefix numbers can be used to abbreviate multiple occurrences, e.g., 5X3#. Special meaning have '*' and '+'. Maximally one wildcard is allowed. '*' expands to skipping, '#'; '+' expands to an including 'X', such that the end of the string matches the end of the iStride position. A zero copylist length is replaced by full copy = '+'.

token_maxLen(pfStr):

-- return longest token length separated by the global separtor stop at element equal to 0

Miscellaneous:

printfv(fmt, pfX, split)

-- printf vector pfX with format ftm (e.g. "%11.5g" or "%7.4f"). Elements are separated by space or newlines, if (split>0) after each split items. If (split==0) no newline is printed. If (split<0) the vector is printed transposed with -split rows.

     printfv("%7.5g", Y, 10); // print Y with 5 sig digits on the terminal

popup(string)

-- write string to a popup alert window (with standard neo interface).

size(pfX)

-- returns number of elements vector pfX[]; Does not count NaN elements. Same as dimof(pfX)-cntNaN

gauss(std):

-- (DEPRECATED) return single random variable with gauss distribution, zero mean and standard deviation =std. Use standard prog_unit grnd(std) function.

poissonbook(pfBook, fMean):

-- return distribution function for Poisson distribution with mean and variance equal fMean. Return probability of too short book. Generate Poisson random variable: USAGE:

      pfBook=resize(20);
      poissonBook(pfBook,1.3);
      iPoisson=locateindex(pfBook,rnd(1));

tic(num)

-- restarts realtime clock .num (modulo tic_toc_MAXCLOCK) and return num of millisec since last restart of the clock num (all clocks are started when lib gets loaded the first time).

toc(num)

-- reads realtime clock num (modulo tic_toc_MAXCLOCK) in milliseconds.

tic_toc_MAXCLOCK

-- number of tic-toc clocks. Note: all prog_units share the same clock set!

pow(x,y)

-- compute power x^y (redundant since "x^y" is valid)

atan2(y,x)

-- compute atan2(y,x)

DEG2RAD

-- factor pi/180 = 0.0174532925199433

RAD2DEG

-- factor 180/pi = 57.2957795130823

MINFLOAT

-- smallest positive float value >0 (compiler constant)

MAXFLOAT

-- largest positive float value.

setDimMatchIsStrict(mode)

-- set .dimMatchIsStrict to 0 or 1, which is honored by some functions (as indicated there); 1=strict, dimension mismatch is considered an error, return -1. 0=tolerant, silently perform operation on the largest common element set (=min(|pfX|,|pfY|) for a simple two vector operation)

dimMatchIsStrict

-- variable to read the current state (see above), which is global to all prog_unit.

objStore(pfX)

-- store in global obj storage. Return identifier. The vector elements are retrievable using objRecall() until progunit-delation/rebuild or explicit removal by objDelete().

objRecall(pfX, fIdent, iMode)

-- recall from obj storage. Retrieve vector-obj and copy as many elements as possible to pfX[]. Return number of elements or negative required length.

       iMode=0: silent anyway
       iMode&1: warn if not present (return value=0)
       iMode&2: warn if pfX[] too small (return value <0)
       iMode&4: pop-up alert, not terminal

objDelete(fIdent)

-- remove from obj storage. Return 1 on sucess, 0 if not present, or negative on internal error

objList(fIdent)

-- list obj storage. If fIdent!=0 ident this element, else list all. Return number of listed elements.

isFingerprint(pfPrint,pfData,fAux)

-- compute fingerprint of vector pfData[] and float fAux and store it in pfPrint[2]. Return 0 if previous fingerprint is unequal to new, else 0; dimof(pfPrint) must be minimal 2; return -1 on error.

svd(m,n)

-- Sigular Value Decomposition class

      svd(m,n) / * constructor (m,n), m>=n * /
      newA<*> / * set new A[0..n*m-1], do SVD and zeroing * /
      x_by_b<**> / * solve A x=b, return x * /
      setZeroClipRatio<1> / * set, default 1e-6 * /
      getW<*> / * get diag of W matrix * /
      getU<*.> / * get col i of U, return W_ii [m] * /
      getV<*.> / * get col i of V, return V_ii [n] * /
      cntZero<0> / * number of zeroed W elements * /
      getDet<0> / * det of W = trace (after zeroing) * /
      getDetOrig<0> / * det of W (before zeroing) * /
      getWmax<0> / * max w value * /
      print(iFlag) / * print to stdout, iFlag select output * /
                                    iFlag & 0x1 print W-matrix
                                    iFlag & 0x2 print U-matrix
                                    iFlag & 0x4 print U-matrix

See also:

      neo -i +/homes/walter/p/nst5/FOLDERS/jwexpl

FILE