knetfile.{h,c} Random access to remote files
kopen.c Smart stream opening
-kfunc.c Special mathematical functions
-krand.{h,c} Pseudo-random number generator mt19937
-kmin.{h,c} Derivative-free non-linear programming
+kmath.{h,c} Numerical routines
khmm.{h,c} Basic HMM library
ksw.(h,c} Smith-Waterman using SSE2
knhx.{h,c} Newick format parser
+++ /dev/null
-#include <math.h>
-
-
-/* Log gamma function
- * \log{\Gamma(z)}
- * AS245, 2nd algorithm, http://lib.stat.cmu.edu/apstat/245
- */
-double kf_lgamma(double z)
-{
- double x = 0;
- x += 0.1659470187408462e-06 / (z+7);
- x += 0.9934937113930748e-05 / (z+6);
- x -= 0.1385710331296526 / (z+5);
- x += 12.50734324009056 / (z+4);
- x -= 176.6150291498386 / (z+3);
- x += 771.3234287757674 / (z+2);
- x -= 1259.139216722289 / (z+1);
- x += 676.5203681218835 / z;
- x += 0.9999999999995183;
- return log(x) - 5.58106146679532777 - z + (z-0.5) * log(z+6.5);
-}
-
-/* complementary error function
- * \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt
- * AS66, 2nd algorithm, http://lib.stat.cmu.edu/apstat/66
- */
-double kf_erfc(double x)
-{
- const double p0 = 220.2068679123761;
- const double p1 = 221.2135961699311;
- const double p2 = 112.0792914978709;
- const double p3 = 33.912866078383;
- const double p4 = 6.37396220353165;
- const double p5 = .7003830644436881;
- const double p6 = .03526249659989109;
- const double q0 = 440.4137358247522;
- const double q1 = 793.8265125199484;
- const double q2 = 637.3336333788311;
- const double q3 = 296.5642487796737;
- const double q4 = 86.78073220294608;
- const double q5 = 16.06417757920695;
- const double q6 = 1.755667163182642;
- const double q7 = .08838834764831844;
- double expntl, z, p;
- z = fabs(x) * M_SQRT2;
- if (z > 37.) return x > 0.? 0. : 2.;
- expntl = exp(z * z * - .5);
- if (z < 10. / M_SQRT2) // for small z
- p = expntl * ((((((p6 * z + p5) * z + p4) * z + p3) * z + p2) * z + p1) * z + p0)
- / (((((((q7 * z + q6) * z + q5) * z + q4) * z + q3) * z + q2) * z + q1) * z + q0);
- else p = expntl / 2.506628274631001 / (z + 1. / (z + 2. / (z + 3. / (z + 4. / (z + .65)))));
- return x > 0.? 2. * p : 2. * (1. - p);
-}
-
-/* The following computes regularized incomplete gamma functions.
- * Formulas are taken from Wiki, with additional input from Numerical
- * Recipes in C (for modified Lentz's algorithm) and AS245
- * (http://lib.stat.cmu.edu/apstat/245).
- *
- * A good online calculator is available at:
- *
- * http://www.danielsoper.com/statcalc/calc23.aspx
- *
- * It calculates upper incomplete gamma function, which equals
- * kf_gammaq(s,z)*tgamma(s).
- */
-
-#define KF_GAMMA_EPS 1e-14
-#define KF_TINY 1e-290
-
-// regularized lower incomplete gamma function, by series expansion
-static double _kf_gammap(double s, double z)
-{
- double sum, x;
- int k;
- for (k = 1, sum = x = 1.; k < 100; ++k) {
- sum += (x *= z / (s + k));
- if (x / sum < KF_GAMMA_EPS) break;
- }
- return exp(s * log(z) - z - kf_lgamma(s + 1.) + log(sum));
-}
-// regularized upper incomplete gamma function, by continued fraction
-static double _kf_gammaq(double s, double z)
-{
- int j;
- double C, D, f;
- f = 1. + z - s; C = f; D = 0.;
- // Modified Lentz's algorithm for computing continued fraction
- // See Numerical Recipes in C, 2nd edition, section 5.2
- for (j = 1; j < 100; ++j) {
- double a = j * (s - j), b = (j<<1) + 1 + z - s, d;
- D = b + a * D;
- if (D < KF_TINY) D = KF_TINY;
- C = b + a / C;
- if (C < KF_TINY) C = KF_TINY;
- D = 1. / D;
- d = C * D;
- f *= d;
- if (fabs(d - 1.) < KF_GAMMA_EPS) break;
- }
- return exp(s * log(z) - z - kf_lgamma(s) - log(f));
-}
-
-double kf_gammap(double s, double z)
-{
- return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);
-}
-
-double kf_gammaq(double s, double z)
-{
- return z <= 1. || z < s? 1. - _kf_gammap(s, z) : _kf_gammaq(s, z);
-}
-
-/* Regularized incomplete beta function. The method is taken from
- * Numerical Recipe in C, 2nd edition, section 6.4. The following web
- * page calculates the incomplete beta function, which equals
- * kf_betai(a,b,x) * gamma(a) * gamma(b) / gamma(a+b):
- *
- * http://www.danielsoper.com/statcalc/calc36.aspx
- */
-static double kf_betai_aux(double a, double b, double x)
-{
- double C, D, f;
- int j;
- if (x == 0.) return 0.;
- if (x == 1.) return 1.;
- f = 1.; C = f; D = 0.;
- // Modified Lentz's algorithm for computing continued fraction
- for (j = 1; j < 200; ++j) {
- double aa, d;
- int m = j>>1;
- aa = (j&1)? -(a + m) * (a + b + m) * x / ((a + 2*m) * (a + 2*m + 1))
- : m * (b - m) * x / ((a + 2*m - 1) * (a + 2*m));
- D = 1. + aa * D;
- if (D < KF_TINY) D = KF_TINY;
- C = 1. + aa / C;
- if (C < KF_TINY) C = KF_TINY;
- D = 1. / D;
- d = C * D;
- f *= d;
- if (fabs(d - 1.) < KF_GAMMA_EPS) break;
- }
- return exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x)) / a / f;
-}
-double kf_betai(double a, double b, double x)
-{
- return x < (a + 1.) / (a + b + 2.)? kf_betai_aux(a, b, x) : 1. - kf_betai_aux(b, a, 1. - x);
-}
-
-#ifdef KF_MAIN
-#include <stdio.h>
-int main(int argc, char *argv[])
-{
- double x = 5.5, y = 3;
- double a, b;
- printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
- printf("upper-gamma(%lg,%lg): %lg\n", x, y, kf_gammaq(y, x)*tgamma(y));
- a = 2; b = 2; x = 0.5;
- printf("incomplete-beta(%lg,%lg,%lg): %lg\n", a, b, x, kf_betai(a, b, x) / exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b)));
- return 0;
-}
-#endif
+++ /dev/null
-/* The MIT License
-
- Copyright (c) 2008, 2010 by Attractive Chaos <attractor@live.co.uk>
-
- Permission is hereby granted, free of charge, to any person obtaining
- a copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
-
- The above copyright notice and this permission notice shall be
- included in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
- NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
- BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
- ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
- CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
- SOFTWARE.
-*/
-
-/* Hooke-Jeeves algorithm for nonlinear minimization
-
- Based on the pseudocodes by Bell and Pike (CACM 9(9):684-685), and
- the revision by Tomlin and Smith (CACM 12(11):637-638). Both of the
- papers are comments on Kaupe's Algorithm 178 "Direct Search" (ACM
- 6(6):313-314). The original algorithm was designed by Hooke and
- Jeeves (ACM 8:212-229). This program is further revised according to
- Johnson's implementation at Netlib (opt/hooke.c).
-
- Hooke-Jeeves algorithm is very simple and it works quite well on a
- few examples. However, it might fail to converge due to its heuristic
- nature. A possible improvement, as is suggested by Johnson, may be to
- choose a small r at the beginning to quickly approach to the minimum
- and a large r at later step to hit the minimum.
- */
-
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "kmin.h"
-
-static double __kmin_hj_aux(kmin_f func, int n, double *x1, void *data, double fx1, double *dx, int *n_calls)
-{
- int k, j = *n_calls;
- double ftmp;
- for (k = 0; k != n; ++k) {
- x1[k] += dx[k];
- ftmp = func(n, x1, data); ++j;
- if (ftmp < fx1) fx1 = ftmp;
- else { /* search the opposite direction */
- dx[k] = 0.0 - dx[k];
- x1[k] += dx[k] + dx[k];
- ftmp = func(n, x1, data); ++j;
- if (ftmp < fx1) fx1 = ftmp;
- else x1[k] -= dx[k]; /* back to the original x[k] */
- }
- }
- *n_calls = j;
- return fx1; /* here: fx1=f(n,x1) */
-}
-
-double kmin_hj(kmin_f func, int n, double *x, void *data, double r, double eps, int max_calls)
-{
- double fx, fx1, *x1, *dx, radius;
- int k, n_calls = 0;
- x1 = (double*)calloc(n, sizeof(double));
- dx = (double*)calloc(n, sizeof(double));
- for (k = 0; k != n; ++k) { /* initial directions, based on MGJ */
- dx[k] = fabs(x[k]) * r;
- if (dx[k] == 0) dx[k] = r;
- }
- radius = r;
- fx1 = fx = func(n, x, data); ++n_calls;
- for (;;) {
- memcpy(x1, x, n * sizeof(double)); /* x1 = x */
- fx1 = __kmin_hj_aux(func, n, x1, data, fx, dx, &n_calls);
- while (fx1 < fx) {
- for (k = 0; k != n; ++k) {
- double t = x[k];
- dx[k] = x1[k] > x[k]? fabs(dx[k]) : 0.0 - fabs(dx[k]);
- x[k] = x1[k];
- x1[k] = x1[k] + x1[k] - t;
- }
- fx = fx1;
- if (n_calls >= max_calls) break;
- fx1 = func(n, x1, data); ++n_calls;
- fx1 = __kmin_hj_aux(func, n, x1, data, fx1, dx, &n_calls);
- if (fx1 >= fx) break;
- for (k = 0; k != n; ++k)
- if (fabs(x1[k] - x[k]) > .5 * fabs(dx[k])) break;
- if (k == n) break;
- }
- if (radius >= eps) {
- if (n_calls >= max_calls) break;
- radius *= r;
- for (k = 0; k != n; ++k) dx[k] *= r;
- } else break; /* converge */
- }
- free(x1); free(dx);
- return fx1;
-}
+++ /dev/null
-/*
- Copyright (c) 2008, 2010 by Attractive Chaos <attractor@live.co.uk>
-
- Permission is hereby granted, free of charge, to any person obtaining
- a copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
-
- The above copyright notice and this permission notice shall be
- included in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
- NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
- BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
- ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
- CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
- SOFTWARE.
-*/
-
-#ifndef KMIN_H
-#define KMIN_H
-
-#define KMIN_RADIUS 0.5
-#define KMIN_EPS 1e-7
-#define KMIN_MAXCALL 50000
-
-typedef double (*kmin_f)(int, double*, void*);
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
- double kmin_hj(kmin_f func, int n, double *x, void *data, double r, double eps, int max_calls);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
+++ /dev/null
-/*
- 64-bit Mersenne Twister pseudorandom number generator. Adapted from:
-
- http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
-
- which was written by Takuji Nishimura and Makoto Matsumoto and released
- under the 3-clause BSD license.
-*/
-
-#include <stdlib.h>
-#include "krand.h"
-
-#define KR_NN 312
-#define KR_MM 156
-#define KR_UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
-#define KR_LM 0x7FFFFFFFULL /* Least significant 31 bits */
-
-struct _krand_t {
- int mti;
- krint64_t mt[KR_NN];
-};
-
-static void kr_srand0(krint64_t seed, krand_t *kr)
-{
- kr->mt[0] = seed;
- for (kr->mti = 1; kr->mti < KR_NN; ++kr->mti)
- kr->mt[kr->mti] = 6364136223846793005ULL * (kr->mt[kr->mti - 1] ^ (kr->mt[kr->mti - 1] >> 62)) + kr->mti;
-}
-
-krand_t *kr_srand(krint64_t seed)
-{
- krand_t *kr;
- kr = malloc(sizeof(krand_t));
- kr_srand0(seed, kr);
- return kr;
-}
-
-krint64_t kr_rand(krand_t *kr)
-{
- krint64_t x;
- static const krint64_t mag01[2] = { 0, 0xB5026F5AA96619E9ULL };
- if (kr->mti >= KR_NN) {
- int i;
- if (kr->mti == KR_NN + 1) kr_srand0(5489ULL, kr);
- for (i = 0; i < KR_NN - KR_MM; ++i) {
- x = (kr->mt[i] & KR_UM) | (kr->mt[i+1] & KR_LM);
- kr->mt[i] = kr->mt[i + KR_MM] ^ (x>>1) ^ mag01[(int)(x&1)];
- }
- for (; i < KR_NN - 1; ++i) {
- x = (kr->mt[i] & KR_UM) | (kr->mt[i+1] & KR_LM);
- kr->mt[i] = kr->mt[i + (KR_MM - KR_NN)] ^ (x>>1) ^ mag01[(int)(x&1)];
- }
- x = (kr->mt[KR_NN - 1] & KR_UM) | (kr->mt[0] & KR_LM);
- kr->mt[KR_NN - 1] = kr->mt[KR_MM - 1] ^ (x>>1) ^ mag01[(int)(x&1)];
- kr->mti = 0;
- }
- x = kr->mt[kr->mti++];
- x ^= (x >> 29) & 0x5555555555555555ULL;
- x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
- x ^= (x << 37) & 0xFFF7EEE000000000ULL;
- x ^= (x >> 43);
- return x;
-}
-
-#ifdef _KR_MAIN
-int main(int argc, char *argv[])
-{
- long i, N = 200000000;
- krand_t *kr;
- if (argc > 1) N = atol(argv[1]);
- kr = kr_srand(11);
- for (i = 0; i < N; ++i) kr_rand(kr);
-// for (i = 0; i < N; ++i) lrand48();
- free(kr);
- return 0;
-}
-#endif
+++ /dev/null
-#ifndef AC_KRAND_H
-#define AC_KRAND_H
-
-#include <stdint.h>
-
-typedef uint64_t krint64_t;
-
-struct _krand_t;
-typedef struct _krand_t krand_t;
-
-#define kr_drand(_kr) ((kr_rand(_kr) >> 11) * (1.0/9007199254740992.0))
-#define kr_sample(_kr, _k, _cnt) ((*(_cnt))++ < (_k)? *(_cnt) - 1 : kr_rand(_kr) % *(_cnt))
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-krand_t *kr_srand(krint64_t seed);
-krint64_t kr_rand(krand_t *kr);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
projects / klib.git / commitdiff