/* Operations on vectors of floats or doubles. * * Can operate on vectors of doubles, floats, or integers - appropriate * routine is prefixed with D, F, or I. For example, esl_vec_DSet() is * the Set routine for a vector of doubles; esl_vec_ISet() is for integers. * * Contents: * 1. The vectorops API. * 2. Unit tests. * 3. Test driver. * 4. Examples. */ #include "esl_config.h" #include #include #include "easel.h" #include "esl_random.h" #include "esl_vectorops.h" /* Function: esl_vec_{DFIL}Set() * Synopsis: Set all items in vector to scalar value. * * Purpose: Sets all items in to . */ void esl_vec_DSet(double *vec, int n, double value) { int i; for (i = 0; i < n; i++) vec[i] = value; } void esl_vec_FSet(float *vec, int n, float value) { int i; for (i = 0; i < n; i++) vec[i] = value; } void esl_vec_ISet(int *vec, int n, int value) { int i; for (i = 0; i < n; i++) vec[i] = value; } void esl_vec_LSet(int64_t *vec, int n, int64_t value) { int i; for (i = 0; i < n; i++) vec[i] = value; } /* Function: esl_vec_{DFIL}Scale() * Synopsis: Multiply all items in vector by scalar value. * * Purpose: Multiplies all items in by . */ void esl_vec_DScale(double *vec, int n, double scale) { int i; for (i = 0; i < n; i++) vec[i] *= scale; } void esl_vec_FScale(float *vec, int n, float scale) { int i; for (i = 0; i < n; i++) vec[i] *= scale; } void esl_vec_IScale(int *vec, int n, int scale) { int i; for (i = 0; i < n; i++) vec[i] *= scale; } void esl_vec_LScale(int64_t *vec, int n, int64_t scale) { int i; for (i = 0; i < n; i++) vec[i] *= scale; } /* Function: esl_vec_{DFIL}Increment() * Synopsis: Add a scalar to all items in a vector. * Incept: SRE, Mon Mar 21 11:56:57 2005 [St. Louis] * * Purpose: Adds scalar to all items in the -vector . */ void esl_vec_DIncrement(double *v, int n, double x) { int i; for (i = 0; i < n; i++) v[i] += x; } void esl_vec_FIncrement(float *v, int n, float x) { int i; for (i = 0; i < n; i++) v[i] += x; } void esl_vec_IIncrement(int *v, int n, int x) { int i; for (i = 0; i < n; i++) v[i] += x; } void esl_vec_LIncrement(int64_t *v, int n, int64_t x) { int i; for (i = 0; i < n; i++) v[i] += x; } /* Function: esl_vec_{DFIL}Add() * Synopsis: Vector addition of two vectors. * * Purpose: Vector addition. Adds to , leaving * result in . ( is unchanged.). * Both vectors are of size . */ void esl_vec_DAdd(double *vec1, const double *vec2, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i]; } void esl_vec_FAdd(float *vec1, const float *vec2, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i]; } void esl_vec_IAdd(int *vec1, const int *vec2, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i]; } void esl_vec_LAdd(int64_t *vec1, const int64_t *vec2, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i]; } /* Function: esl_vec_{DFIL}AddScaled() * Synopsis: Scale and add it to . * * Purpose: Scales by scalar , and adds that * to . Both vectors are of size . */ void esl_vec_DAddScaled(double *vec1, const double *vec2, double a, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i] * a; } void esl_vec_FAddScaled(float *vec1, const float *vec2, float a, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i] * a; } void esl_vec_IAddScaled(int *vec1, const int *vec2, int a, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i] * a; } void esl_vec_LAddScaled(int64_t *vec1, const int64_t *vec2, int64_t a, int n) { int i; for (i = 0; i < n; i++) vec1[i] += vec2[i] * a; } /* Function: esl_vec_{DFIL}Sum() * Synopsis: Returns $\sum_i x_i$. * * Purpose: Returns the scalar sum of the items in . * * Floating point summations use Kahan compensated * summation, in order to minimize roundoff error * accumulation. Additionally, they are most * accurate if vec[] is sorted in increasing order, from * small to large, so you may consider sorting before * summing it. */ double esl_vec_DSum(const double *vec, int n) { double sum = 0.; double y,t,c; int i; c = 0.0; for (i = 0; i < n; i++) { y = vec[i] - c; t = sum + y; c = (t-sum)-y; sum = t; } return sum; } float esl_vec_FSum(const float *vec, int n) { float sum = 0.; float y,t,c; int i; c = 0.0; for (i = 0; i < n; i++) { y = vec[i] - c; t = sum + y; c = (t-sum)-y; sum = t; } return sum; } int esl_vec_ISum(const int *vec, int n) { int sum = 0; int i; for (i = 0; i < n; i++) sum += vec[i]; return sum; } int64_t esl_vec_LSum(const int64_t *vec, int n) { int64_t sum = 0; int i; for (i = 0; i < n; i++) sum += vec[i]; return sum; } /* Function: esl_vec_{DFIL}Dot() * Synopsis: Return the dot product of two vectors. * * Purpose: Returns the scalar dot product $\cdot$ . * Both vectors are of size . */ double esl_vec_DDot(const double *vec1, const double *vec2, int n) { double result = 0.; int i; for (i = 0; i < n; i++) result += vec1[i] * vec2[i]; return result; } float esl_vec_FDot(const float *vec1, const float *vec2, int n) { float result = 0.; int i; for (i = 0; i < n; i++) result += vec1[i] * vec2[i]; return result; } int esl_vec_IDot(const int *vec1, const int *vec2, int n) { int result = 0; int i; for (i = 0; i < n; i++) result += vec1[i] * vec2[i]; return result; } int64_t esl_vec_LDot(const int64_t *vec1, const int64_t *vec2, int n) { int result = 0; int i; for (i = 0; i < n; i++) result += vec1[i] * vec2[i]; return result; } /* Function: esl_vec_{DFIL}Max() * Synopsis: Return value of the maximum element in a vector. * * Purpose: Returns the maximum value of the values * in . */ double esl_vec_DMax(const double *vec, int n) { double best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] > best) best = vec[i]; return best; } float esl_vec_FMax(const float *vec, int n) { float best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] > best) best = vec[i]; return best; } int esl_vec_IMax(const int *vec, int n) { int best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] > best) best = vec[i]; return best; } int64_t esl_vec_LMax(const int64_t *vec, int n) { int64_t best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] > best) best = vec[i]; return best; } /* Function: esl_vec_{DFIL}Min() * Synopsis: Return value of the minimum element in a vector. * * Purpose: Returns the minimum value of the values * in . */ double esl_vec_DMin(const double *vec, int n) { double best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] < best) best = vec[i]; return best; } float esl_vec_FMin(const float *vec, int n) { float best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] < best) best = vec[i]; return best; } int esl_vec_IMin(const int *vec, int n) { int best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] < best) best = vec[i]; return best; } int64_t esl_vec_LMin(const int64_t *vec, int n) { int64_t best; int i; best = vec[0]; for (i = 1; i < n; i++) if (vec[i] < best) best = vec[i]; return best; } /* Function: esl_vec_{DFIL}ArgMax() * Synopsis: Return index of maximum element in a vector. * * Purpose: Returns the index of the maximum value in the values * in . In case of ties, return the smallest index. * * can be 0 and can be , in which case the * function returns 0. * * Note: Do not change the behavior that the smallest index is * returned in case of ties. Some functions rely on this * behavior: optimal accuracy tracebacks in HMMER for example. */ int esl_vec_DArgMax(const double *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] > vec[best]) best = i; return best; } int esl_vec_FArgMax(const float *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] > vec[best]) best = i; return best; } int esl_vec_IArgMax(const int *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] > vec[best]) best = i; return best; } int esl_vec_LArgMax(const int64_t *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] > vec[best]) best = i; return best; } /* Function: esl_vec_{DFIL}ArgMin() * Synopsis: Return index of minimum element in a vector. * * Purpose: Returns the index of the minimum value in the values * in . */ int esl_vec_DArgMin(const double *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] < vec[best]) best = i; return best; } int esl_vec_FArgMin(const float *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] < vec[best]) best = i; return best; } int esl_vec_IArgMin(const int *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] < vec[best]) best = i; return best; } int esl_vec_LArgMin(const int64_t *vec, int n) { int i; int best = 0; for (i = 1; i < n; i++) if (vec[i] < vec[best]) best = i; return best; } /* Function: esl_vec_{DFILWB}Copy() * Synopsis: Set vector to same values as . * * Purpose: Copies to . is * unchanged. Both vectors are of size . */ void esl_vec_DCopy(const double *src, int n, double *dest) { int i; for (i = 0; i < n; i++) dest[i] = src[i]; } void esl_vec_FCopy(const float *src, int n, float *dest) { int i; for (i = 0; i < n; i++) dest[i] = src[i]; } void esl_vec_ICopy(const int *src, int n, int *dest) { int i; for (i = 0; i < n; i++) dest[i] = src[i]; } void esl_vec_LCopy(const int64_t *src, int n, int64_t *dest) { int i; for (i = 0; i < n; i++) dest[i] = src[i]; } void esl_vec_WCopy(const int16_t *src, int n, int16_t *dest) { int i; for (i = 0; i < n; i++) dest[i] = src[i]; } void esl_vec_BCopy(const int8_t *src, int n, int8_t *dest) { int i; for (i = 0; i < n; i++) dest[i] = src[i]; } /* Function: esl_vec_{DFIL}Swap() * Synopsis: Swap two vectors. * * Purpose: Swaps and , both of size . * * But you're better off swapping the pointers to the vectors, * if that's feasible. */ void esl_vec_DSwap(double *vec1, double *vec2, int n) { int i; double tmp; for (i = 0; i < n; i++) { tmp = vec1[i]; vec1[i] = vec2[i]; vec2[i] = tmp; } } void esl_vec_FSwap(float *vec1, float *vec2, int n) { int i; float tmp; for (i = 0; i < n; i++) { tmp = vec1[i]; vec1[i] = vec2[i]; vec2[i] = tmp; } } void esl_vec_ISwap(int *vec1, int *vec2, int n) { int i; int tmp; for (i = 0; i < n; i++) { tmp = vec1[i]; vec1[i] = vec2[i]; vec2[i] = tmp; } } void esl_vec_LSwap(int64_t *vec1, int64_t *vec2, int n) { int i; int tmp; for (i = 0; i < n; i++) { tmp = vec1[i]; vec1[i] = vec2[i]; vec2[i] = tmp; } } /* Function: esl_vec_{DFILC}Reverse() * Synopsis: Reverse a vector (possibly in place). * * Purpose: Put the values from vector in reversed order in * . Caller provides storage in for at least * values. * * and can be the same, in which case is * reversed in place. * * needs to be used carefully if * is a NUL-terminated string, rather than an array. * If you reverse a string in place (i.e. * ), the trailing NUL will * still be there, and you're fine. If you reverse string * into new storage , you'll need to NUL-terminate * yourself. */ void esl_vec_DReverse(const double *vec, double *rev, int n) { int i; double x; for (i = 0; i < n/2; i++) { x = vec[n-i-1]; rev[n-i-1] = vec[i]; rev[i] = x; } if (n%2) rev[i] = vec[i]; } void esl_vec_FReverse(const float *vec, float *rev, int n) { int i; float x; for (i = 0; i < n/2; i++) { x = vec[n-i-1]; rev[n-i-1] = vec[i]; rev[i] = x; } if (n%2) rev[i] = vec[i]; } void esl_vec_IReverse(const int *vec, int *rev, int n) { int i; int x; for (i = 0; i < n/2; i++) { x = vec[n-i-1]; rev[n-i-1] = vec[i]; rev[i] = x; } if (n%2) rev[i] = vec[i]; } void esl_vec_LReverse(const int64_t *vec, int64_t *rev, int n) { int i; int64_t x; for (i = 0; i < n/2; i++) { x = vec[n-i-1]; rev[n-i-1] = vec[i]; rev[i] = x; } if (n%2) rev[i] = vec[i]; } void esl_vec_CReverse(const char *vec, char *rev, int n) { int i; char x; for (i = 0; i < n/2; i++) { x = vec[n-i-1]; rev[n-i-1] = vec[i]; rev[i] = x; } if (n%2) rev[i] = vec[i]; } /* some static functions to pass to qsort() that the * upcoming Sort() functions will call */ static int qsort_DIncreasing(const void *xp1, const void *xp2) { double x1 = * (double *) xp1; double x2 = * (double *) xp2; if (x1 < x2) return -1; if (x1 > x2) return 1; return 0; } static int qsort_FIncreasing(const void *xp1, const void *xp2) { float x1 = * (float *) xp1; float x2 = * (float *) xp2; if (x1 < x2) return -1; if (x1 > x2) return 1; return 0; } static int qsort_IIncreasing(const void *xp1, const void *xp2) { int x1 = * (int *) xp1; int x2 = * (int *) xp2; if (x1 < x2) return -1; if (x1 > x2) return 1; return 0; } static int qsort_LIncreasing(const void *xp1, const void *xp2) { int64_t x1 = * (int64_t *) xp1; int64_t x2 = * (int64_t *) xp2; if (x1 < x2) return -1; if (x1 > x2) return 1; return 0; } static int qsort_DDecreasing(const void *xp1, const void *xp2) { double x1 = * (double *) xp1; double x2 = * (double *) xp2; if (x1 > x2) return -1; if (x1 < x2) return 1; return 0; } static int qsort_FDecreasing(const void *xp1, const void *xp2) { float x1 = * (float *) xp1; float x2 = * (float *) xp2; if (x1 > x2) return -1; if (x1 < x2) return 1; return 0; } static int qsort_IDecreasing(const void *xp1, const void *xp2) { int x1 = * (int *) xp1; int x2 = * (int *) xp2; if (x1 > x2) return -1; if (x1 < x2) return 1; return 0; } static int qsort_LDecreasing(const void *xp1, const void *xp2) { int64_t x1 = * (int64_t *) xp1; int64_t x2 = * (int64_t *) xp2; if (x1 > x2) return -1; if (x1 < x2) return 1; return 0; } /* Function: esl_vec_{DFIL}SortIncreasing() * Synopsis: Sort vector from smallest to largest. * Incept: SRE, Wed Aug 17 10:44:31 2005 [St. Louis] * * Purpose: Sorts in place, from smallest to largest value. * (That is, is the minimum and is * the maximum.) */ void esl_vec_DSortIncreasing(double *vec, int n) { qsort((void *) vec, n, sizeof(double), qsort_DIncreasing); } void esl_vec_FSortIncreasing(float *vec, int n) { qsort((void *) vec, n, sizeof(float), qsort_FIncreasing); } void esl_vec_ISortIncreasing(int *vec, int n) { qsort((void *) vec, n, sizeof(int), qsort_IIncreasing); } void esl_vec_LSortIncreasing(int64_t *vec, int n) { qsort((void *) vec, n, sizeof(int64_t), qsort_LIncreasing); } /* Function: esl_vec_{DFIL}SortDecreasing() * Synopsis: Sort vector from largest to smallest. * Incept: SRE, Wed Aug 17 10:44:31 2005 [St. Louis] * * Purpose: Sorts in place, from largest to smallest value. * (That is, is the maximum and is * the minimum.) */ void esl_vec_DSortDecreasing(double *vec, int n) { qsort((void *) vec, n, sizeof(double), qsort_DDecreasing); } void esl_vec_FSortDecreasing(float *vec, int n) { qsort((void *) vec, n, sizeof(float), qsort_FDecreasing); } void esl_vec_ISortDecreasing(int *vec, int n) { qsort((void *) vec, n, sizeof(int), qsort_IDecreasing); } void esl_vec_LSortDecreasing(int64_t *vec, int n) { qsort((void *) vec, n, sizeof(int64_t), qsort_LDecreasing); } /* Function: esl_vec_{DFIL}Shuffle() * Synopsis: Shuffle a vector, in place. * * Purpose: Shuffle a vector of items, using the * random number generator . */ int esl_vec_DShuffle(ESL_RANDOMNESS *r, double *v, int n) { double swap; int pos; for ( ; n > 1; n--) { pos = esl_rnd_Roll(r, n); swap = v[pos]; v[pos] = v[n-1]; v[n-1] = swap; } return eslOK; } int esl_vec_FShuffle(ESL_RANDOMNESS *r, float *v, int n) { float swap; int pos; for ( ; n > 1; n--) { pos = esl_rnd_Roll(r, n); swap = v[pos]; v[pos] = v[n-1]; v[n-1] = swap; } return eslOK; } int esl_vec_IShuffle(ESL_RANDOMNESS *r, int *v, int n) { int swap; int pos; for ( ; n > 1; n--) { pos = esl_rnd_Roll(r, n); swap = v[pos]; v[pos] = v[n-1]; v[n-1] = swap; } return eslOK; } int esl_vec_LShuffle(ESL_RANDOMNESS *r, int64_t *v, int n) { int64_t swap; int pos; for ( ; n > 1; n--) { pos = esl_rnd_Roll(r, n); swap = v[pos]; v[pos] = v[n-1]; v[n-1] = swap; } return eslOK; } /* Function: esl_vec_{DFIL}Compare() * Synopsis: Return if two vectors are equal. * Incept: SRE, Mon Nov 6 10:20:28 2006 [Janelia] * * Purpose: Compare to for equality, by * comparing each cognate element pair. Both vectors * are of size . Equality of elements is * defined by being $\leq$ fractional tolerance * for floating point comparisons, and strict equality * for integer comparisons. Return * if the vectors are equal, and if not. * * If , the test always succeeds. In this case, either * and (or both) may be . This * accommodates an occasional convention of leaving empty * vectors . */ int esl_vec_DCompare(const double *vec1, const double *vec2, int n, double tol) { int i; for (i = 0; i < n; i++) if (esl_DCompare(vec1[i], vec2[i], tol) == eslFAIL) return eslFAIL; return eslOK; } int esl_vec_FCompare(const float *vec1, const float *vec2, int n, float tol) { int i; for (i = 0; i < n; i++) if (esl_DCompare(vec1[i], vec2[i], tol) == eslFAIL) return eslFAIL; return eslOK; } int esl_vec_ICompare(const int *vec1, const int *vec2, int n) { int i; for (i = 0; i < n; i++) if (vec1[i] != vec2[i]) return eslFAIL; return eslOK; } int esl_vec_LCompare(const int64_t *vec1, const int64_t *vec2, int n) { int i; for (i = 0; i < n; i++) if (vec1[i] != vec2[i]) return eslFAIL; return eslOK; } /* Function: esl_vec_{DFIL}Dump() * Synopsis: Output vector to a stream as text. * Incept: ER, Thu Jul 21 12:54:56 CDT 2005 [St. Louis] * * Purpose: Given a vector, dump it to stream . * * If