numerical_recipes.h File Reference

#include <math.h>
#include "xmipp_memory.h"
#include "xmipp_macros.h"
#include <algorithm>
#include <vector>

Include dependency graph for numerical_recipes.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

#define	TINY   1.0e-20;
void	nrerror (const char error_text[])
double	ran1 (int *idum)
double	gasdev (int *idum)
double	tdev (double nu, int *idum)
void	ksone (double data[], int n, double(*func)(double), double *d, double *prob)
double	probks (double alam)
void	four1 (double data[], int nn, int isign)
void	realft (double data[], int n, int isign)
void	fourn (double data[], int nn[], int ndim, int isign)
void	indexx (int n, double arrin[], int indx[])
double	bessj0 (double x)
double	bessj3_5 (double x)
double	bessj1_5 (double x)
double	bessi0 (double x)
double	bessi1 (double x)
double	bessi0_5 (double x)
double	bessi1_5 (double x)
double	bessi2 (double x)
double	bessi3 (double x)
double	bessi2_5 (double x)
double	bessi3_5 (double x)
double	bessi4 (double x)
double	gammln (double xx)
double	gammp (double a, double x)
double	betacf (double a, double b, double x)
double	betai (double a, double b, double x)
double	sinc (double x)
size_t	fact (int num)
size_t	binom (int n, int k)
void	svdcmp (double *a, int m, int n, double *w, double *v)
void	svbksb (double *u, double *w, double *v, int m, int n, double *b, double *x)
void	convlv (double *data, int n, double *respns, int m, int isign, double *ans)
void	realft (double *data, int n, int isign)
void	twofft (double *data1, double *data2, double *fft1, double *fft2, int n)
void	four1 (double *data, int nn, int isign)
void	powell (double *p, double *xi, int n, double ftol, int &iter, double &fret, double(*func)(double *, void *), void *prm, bool show)
int	nnls (double *a, int m, int n, double *b, double *x, double *rnorm, double *w, double *zz, int *index)
int	nnlsWght (int N, int M, double *A, double *b, double *weight)
void	grobfd (int, int, double *, double *, void(*)(int, int, double *, double *, void *), void *)
void	grcnfd (int, int, double *, double *, void(*)(int, int, double *, double *, void *), void *)
void	cfsqp (int, int, int, int, int, int, int, int, int, int *, int, int, int, int *, double, double, double, double, double *, double *, double *, double *, double *, double *, void(*)(int, int, double *, double *, void *), void(*)(int, int, double *, double *, void *), void(*)(int, int, double *, double *, void(*)(int, int, double *, double *, void *), void *), void(*)(int, int, double *, double *, void(*)(int, int, double *, double *, void *), void *), void *)
void	wtn (double a[], unsigned long nn[], int ndim, int isign, void(*wtstep)(double [], unsigned long, int))
void	pwtset (int n)
void	pwt (double a[], unsigned long n, int isign)
template<class T >
void	ludcmp (T *a, int n, int *indx, T *d)
template<class T >
void	lubksb (T *a, int n, int *indx, T b[])
template<class T >
void	gaussj (T *a, int n, T *b, int m)
void	choldc (double *a, int n, double *p)
void	cholsl (double *a, int n, double *p, double *b, double *x)
void	polint (double *xa, double *ya, int n, double x, double &y, double &dy)

Macro Definition Documentation

TINY

#define TINY   1.0e-20;

Definition at line 167 of file numerical_recipes.h.

Function Documentation

bessi0()

double bessi0

(

double

x

)

Definition at line 286 of file numerical_recipes.cpp.

{
    double y, ax, ans;
    if ((ax = fabs(x)) < 3.75)
    {
        y = x / 3.75;
        y *= y;
        ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
                                          + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
    }
    else
    {
        y = 3.75 / ax;
        ans = (exp(ax) / sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
                                      + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
                                                                + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
                                                                                            + y * 0.392377e-2))))))));
    }
    return ans;
}

bessi0_5()

double bessi0_5

(

double

x

)

Definition at line 547 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : sqrt(2 / (PI*x))*sinh(x);
}

bessi1()

double bessi1

(

double

x

)

Definition at line 308 of file numerical_recipes.cpp.

{
    double ax, ans;
    double y;
    if ((ax = fabs(x)) < 3.75)
    {
        y = x / 3.75;
        y *= y;
        ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
                               + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
    }
    else
    {
        y = 3.75 / ax;
        ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
                                  - y * 0.420059e-2));
        ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2
                                + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
        ans *= (exp(ax) / sqrt(ax));
    }
    return x < 0.0 ? -ans : ans;
}

bessi1_5()

double bessi1_5

(

double

x

)

Definition at line 551 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : sqrt(2 / (PI*x))*(cosh(x) - sinh(x) / x);
}

bessi2()

double bessi2

(

double

x

)

Definition at line 555 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : bessi0(x) - ((2*1) / x) * bessi1(x);
}

bessi2_5()

double bessi2_5

(

double

x

)

Definition at line 559 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : bessi0_5(x) - ((2*1.5) / x) * bessi1_5(x);
}

bessi3()

double bessi3

(

double

x

)

Definition at line 563 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : bessi1(x) - ((2*2) / x) * bessi2(x);
}

bessi3_5()

double bessi3_5

(

double

x

)

Definition at line 567 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : bessi1_5(x) - ((2*2.5) / x) * bessi2_5(x);
}

bessi4()

double bessi4

(

double

x

)

Definition at line 571 of file numerical_recipes.cpp.

{
    return (x == 0) ? 0 : bessi2(x) - ((2*3) / x) * bessi3(x);
}

bessj0()

double bessj0

(

double

x

)

Definition at line 250 of file numerical_recipes.cpp.

{
    double ax, z;
    double xx, y, ans, ans1, ans2;

    if ((ax = fabs(x)) < 8.0)
    {
        y = x * x;
        ans1 = 57568490574.0 + y * (-13362590354.0 +
                                    y * (651619640.7
                                         + y * (-11214424.18 +
                                                y * (77392.33017 +
                                                     y * (-184.9052456)))));
        ans2 = 57568490411.0 + y * (1029532985.0 +
                                    y * (9494680.718
                                         + y * (59272.64853 +
                                                y * (267.8532712 +
                                                     y * 1.0))));
        ans = ans1 / ans2;
    }
    else
    {
        z = 8.0 / ax;
        y = z * z;
        xx = ax - 0.785398164;
        ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                          + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
        ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                                       + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                                                                      - y * 0.934935152e-7)));
        ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
    }
    return ans;
}

bessj1_5()

double bessj1_5

(

double

x

)

Definition at line 575 of file numerical_recipes.cpp.

{
    double rj, ry, rjp, ryp;
    bessjy(x, 1.5, &rj, &ry, &rjp, &ryp);
    return rj;
}

bessj3_5()

double bessj3_5

(

double

x

)

Definition at line 581 of file numerical_recipes.cpp.

{
    double rj, ry, rjp, ryp;
    bessjy(x, 3.5, &rj, &ry, &rjp, &ryp);
    return rj;
}

betacf()

double betacf	(	double	a,
		double	b,
		double	x
	)

Definition at line 630 of file numerical_recipes.cpp.

{
    double qap, qam, qab, em, tem, d;
    double bz, bm = 1.0, bp, bpp;
    double az = 1.0, am = 1.0, ap, app, aold;
    int m;

    qab = a + b;
    qap = a + 1.0;
    qam = a - 1.0;
    bz = 1.0 - qab * x / qap;
    for (m = 1;m <= ITMAX;m++)
    {
        em = (double) m;
        tem = em + em;
        d = em * (b - em) * x / ((qam + tem) * (a + tem));
        ap = az + d * am;
        bp = bz + d * bm;
        d = -(a + em) * (qab + em) * x / ((qap + tem) * (a + tem));
        app = ap + d * az;
        bpp = bp + d * bz;
        aold = az;
        am = ap / bpp;
        bm = bp / bpp;
        az = app / bpp;
        bz = 1.0;
        if (fabs(az - aold) < (EPS*fabs(az)))
            return az;
    }
    nrerror("a or b too big, or ITMAX too small in BETACF");
    return 0;
}

betai()

double betai	(	double	a,
		double	b,
		double	x
	)

Definition at line 612 of file numerical_recipes.cpp.

{
    double bt;
    if (x < 0.0 || x > 1.0)
        nrerror("Bad x in routine BETAI");
    if (x == 0.0 || x == 1.0)
        bt = 0.0;
    else
        bt = exp(gammln(a + b) - gammln(a) - gammln(b) + a * log(x) + b * log(1.0 - x));
    if (x < (a + 1.0) / (a + b + 2.0))
        return bt*betacf(a, b, x) / a;
    else
        return 1.0 -bt*betacf(b, a, 1.0 - x) / b;

}

binom()

size_t binom ( int  n, int  k  )

inline

Definition at line 103 of file numerical_recipes.h.

{
    size_t factor=1;
    for (int i = n; i > (n-k); --i)
        factor *= i;
    return factor/fact(k);
}

cfsqp()

void cfsqp	(	int	,
		int	,
		int	,
		int	,
		int	,
		int	,
		int	,
		int	,
		int	,
		int *	,
		int	,
		int	,
		int	,
		int *	,
		double	,
		double	,
		double	,
		double	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		void(*)(int, int, double *, double *, void *)	,
		void(*)(int, int, double *, double *, void *)	,
		void(*)(int, int, double *, double *, void(*)(int, int, double *, double *, void *), void *)	,
		void(*)(int, int, double *, double *, void(*)(int, int, double *, double *, void *), void *)	,
		void *
	)

choldc()

void choldc	(	double *	a,
		int	n,
		double *	p
	)

Definition at line 9347 of file numerical_recipes.cpp.

{
    int i, j, k;
    double sum;

    for (i = 1;i <= n;i++)
    {
        for (j = i;j <= n;j++)
        {
            for (sum = a[i*n+j], k = i - 1;k >= 1;k--)
                sum -= a[i*n+k] * a[j*n+k];
            if (i == j)
            {
                if (sum <= 0.0)
                    nrerror("choldc failed");
                p[i] = sqrt(sum);
            }
            else
                a[j*n+i] = sum / p[i];
        }
    }
}

cholsl()

void cholsl	(	double *	a,
		int	n,
		double *	p,
		double *	b,
		double *	x
	)

Definition at line 9370 of file numerical_recipes.cpp.

{
    int i, k;
    double sum;

    for (i = 1;i <= n;i++)
    {
        for (sum = b[i], k = i - 1;k >= 1;k--)
            sum -= a[i*n+k] * x[k];
        x[i] = sum / p[i];
    }
    for (i = n;i >= 1;i--)
    {
        for (sum = x[i], k = i + 1;k <= n;k++)
            sum -= a[k*n+i] * x[k];
        x[i] = sum / p[i];
    }
}

convlv()

void convlv	(	double *	data,
		int	n,
		double *	respns,
		int	m,
		int	isign,
		double *	ans
	)

fact()

size_t fact ( int  num )

inline

Definition at line 91 of file numerical_recipes.h.

{
    size_t value = 1;
    if (num !=1 && num != 0)
    {
        for (int i = num; i > 0; --i)
            value *= i;
    }
    return value;

}

four1() [1/2]

void four1	(	double	data[],
		int	nn,
		int	isign
	)

four1() [2/2]

void four1	(	double *	data,
		int	nn,
		int	isign
	)

fourn()

void fourn	(	double	data[],
		int	nn[],
		int	ndim,
		int	isign
	)

gammln()

double gammln

(

double

xx

)

Definition at line 589 of file numerical_recipes.cpp.

{
    double x, tmp, ser;
    static double cof[6] =
        {
            76.18009173, -86.50532033, 24.01409822,
            -1.231739516, 0.120858003e-2, -0.536382e-5
        };
    int j;

    x = xx - 1.0;
    tmp = x + 5.5;
    tmp -= (x + 0.5) * log(tmp);
    ser = 1.0;
    for (j = 0;j <= 5;j++)
    {
        x += 1.0;
        ser += cof[j] / x;
    }
    return -tmp + log(2.50662827465*ser);
}

gammp()

double gammp	(	double	a,
		double	x
	)

Definition at line 9328 of file numerical_recipes.cpp.

{
    double gamser, gammcf, gln;

    if (x < 0.0 || a <= 0.0)
        nrerror("Invalid arguments in routine gammp");
    if (x < (a + 1.0))
    {
        gser(&gamser, a, x, &gln);
        return gamser;
    }
    else
    {
        gcf(&gammcf, a, x, &gln);
        return 1.0 -gammcf;
    }
}

gasdev()

double gasdev

(

int *

idum

)

Definition at line 106 of file numerical_recipes.cpp.

{
    static int iset = 0;
    static double gset;
    double fac, r, v1, v2;

    if (iset == 0)
    {
        do
        {
            v1 = 2.0 * ran1(idum) - 1.0;
            v2 = 2.0 * ran1(idum) - 1.0;
            r = v1 * v1 + v2 * v2;
        }
        while (r >= 1.0);
        fac = sqrt(-2.0 * log(r) / r);
        gset = v1 * fac;
        iset = 1;
        return v2*fac;
    }
    else
    {
        iset = 0;
        return gset;
    }
}

gaussj()

template<class T >

void gaussj	(	T *	a,
		int	n,
		T *	b,
		int	m
	)

Definition at line 266 of file numerical_recipes.h.

{
    T temp;
    int i, icol=0, irow=0, j, k, l, ll;
    T big, dum;
    double pivinv;

    std::vector<int> buffer(3*n);
    auto *indxc= buffer.data()-1;
    auto *indxr= indxc + n;
    auto *ipiv= indxr + n;
    for (j = 1;j <= n;j++)
        ipiv[j] = 0;
    for (i = 1;i <= n;i++)
    {
        big = (T)0;
        for (j = 1;j <= n;j++)
            if (ipiv[j] != 1)
                for (k = 1;k <= n;k++)
                {
                    if (ipiv[k] == 0)
                    {
                        if (fabs((double)a[j*n+k]) >= (double) big)
                        {
                            big = ABS(a[j*n+k]);
                            irow = j;
                            icol = k;
                        }
                    }
                    else if (ipiv[k] > 1)
                        nrerror("GAUSSJ: Singular Matrix-1");
                }
        ++(ipiv[icol]);
        if (irow != icol)
        {
            for (l = 1;l <= n;l++)
                SWAP(a[irow*n+l], a[icol*n+l], temp)
                for (l = 1;l <= m;l++)
                    SWAP(b[irow*n+l], b[icol*n+l], temp)
                }
        indxr[i] = irow;
        indxc[i] = icol;
        if (a[icol*n+icol] == 0.0)
            nrerror("GAUSSJ: Singular Matrix-2");
        pivinv = 1.0f / a[icol*n+icol];
        a[icol*n+icol] = (T)1;
        for (l = 1;l <= n;l++)
            a[icol*n+l] = (T)(pivinv * a[icol*n+l]);
        for (l = 1;l <= m;l++)
            b[icol*n+l] = (T)(pivinv * b[icol*n+l]);
        for (ll = 1;ll <= n;ll++)
            if (ll != icol)
            {
                dum = a[ll*n+icol];
                a[ll*n+icol] = (T)0;
                for (l = 1;l <= n;l++)
                    a[ll*n+l] -= a[icol*n+l] * dum;
                for (l = 1;l <= m;l++)
                    b[ll*n+l] -= b[icol*n+l] * dum;
            }
    }
    for (l = n;l >= 1;l--)
    {
        if (indxr[l] != indxc[l])
            for (k = 1;k <= n;k++)
                SWAP(a[k*n+indxr[l]], a[k*n+indxc[l]], temp);
    }
}

grcnfd()

void grcnfd	(	int	,
		int	,
		double *	,
		double *	,
		void(*)(int, int, double *, double *, void *)	,
		void *
	)

grobfd()

void grobfd	(	int	,
		int	,
		double *	,
		double *	,
		void(*)(int, int, double *, double *, void *)	,
		void *
	)

indexx()

void indexx	(	int	n,
		double	arrin[],
		int	indx[]
	)

Definition at line 206 of file numerical_recipes.cpp.

{
    int l, j, ir, indxt, i;
    double q;

    for (j = 1;j <= n;j++)
        indx[j] = j;
    l = (n >> 1) + 1;
    ir = n;
    for (;;)
    {
        if (l > 1)
            q = arrin[(indxt=indx[--l])];
        else
        {
            q = arrin[(indxt=indx[ir])];
            indx[ir] = indx[1];
            if (--ir == 1)
            {
                indx[1] = indxt;
                return;
            }
        }
        i = l;
        j = l << 1;
        while (j <= ir)
        {
            if (j < ir && arrin[indx[j]] < arrin[indx[j+1]])
                j++;
            if (q < arrin[indx[j]])
            {
                indx[i] = indx[j];
                j += (i = j);
            }
            else
                j = ir + 1;
        }
        indx[i] = indxt;
    }
}

ksone()

void ksone	(	double	data[],
		int	n,
		double(*)(double)	func,
		double *	d,
		double *	prob
	)

Definition at line 165 of file numerical_recipes.cpp.

{
    std::sort(data, data + n);
    double fn, ff, en, dt, fo=0.;
    en = (double)n;
    *d = 0.;
    for (int j=1; j<=n; j++)
    {
        fn = j / en;
        ff = (*func)(data[j]);
        dt = XMIPP_MAX(fabs(fo - ff), fabs(fn - ff));
        if (dt> *d)
            *d = dt;
        fo = fn;
    }
    *prob = probks(sqrt(en)*(*d));
}

lubksb()

template<class T >

void lubksb	(	T *	a,
		int	n,
		int *	indx,
		T	b[]
	)

Definition at line 237 of file numerical_recipes.h.

{
    int i, ii = 0, ip, j;
    T sum;

    for (i = 1;i <= n;i++)
    {
        ip = indx[i];
        sum = b[ip];
        b[ip] = b[i];
        if (ii)
            for (j = ii;j <= i - 1;j++)
                sum -= a[i*n+j] * b[j];
        else if (sum)
            ii = i;
        b[i] = sum;
    }
    for (i = n;i >= 1;i--)
    {
        sum = b[i];
        for (j = i + 1;j <= n;j++)
            sum -= a[i*n+j] * b[j];
        b[i] = sum / a[i*n+i];
    }
}

ludcmp()

template<class T >

void ludcmp	(	T *	a,
		int	n,
		int *	indx,
		T *	d
	)

Definition at line 170 of file numerical_recipes.h.

{
    int i, imax=0, j, k;
    T big, dum, sum, temp;

    std::vector<T> buffer(n);
    auto *vv= buffer.data()-1;
    *d = (T)1.0;
    for (i = 1;i <= n;i++)
    {
        big = (T)0.0;
        for (j = 1;j <= n;j++)
            if ((temp = (T)fabs((double)a[i*n+j])) > big)
                big = temp;
        if (big == (T)0.0)
            nrerror("Singular matrix in routine LUDCMP");
        vv[i] = (T)1.0 / big;
    }
    for (j = 1;j <= n;j++)
    {
        for (i = 1;i < j;i++)
        {
            sum = a[i*n+j];
            for (k = 1;k < i;k++)
                sum -= a[i*n+k] * a[k*n+j];
            a[i*n+j] = sum;
        }
        big = (T)0.0;
        for (i = j;i <= n;i++)
        {
            sum = a[i*n+j];
            for (k = 1;k < j;k++)
                sum -= a[i*n+k] * a[k*n+j];
            a[i*n+j] = sum;
            if ((dum = vv[i] * (T)fabs((double)sum)) >= big)
            {
                big = dum;
                imax = i;
            }
        }
        if (j != imax)
        {
            for (k = 1;k <= n;k++)
            {
                dum = a[imax*n+k];
                a[imax*n+k] = a[j*n+k];
                a[j*n+k] = dum;
            }
            *d = -(*d);
            vv[imax] = vv[j];
        }
        indx[j] = imax;
        if (a[j*n+j] == 0.0)
            a[j*n+j] = (T) TINY;
        if (j != n)
        {
            dum = (T)1.0 / (a[j*n+j]);
            for (i = j + 1;i <= n;i++)
                a[i*n+j] *= dum;
        }
    }
}

nnls()

int nnls	(	double *	a,
		int	m,
		int	n,
		double *	b,
		double *	x,
		double *	rnorm,
		double *	w,
		double *	zz,
		int *	index
	)

Definition at line 1165 of file numerical_recipes.cpp.

{
    int pfeas, ret = 0, iz, jz, iz1, iz2, npp1, *index;
    double d1, d2, sm, up, ss, *w, *zz;
    int iter, k, j = 0, l, itmax, izmax = 0, nsetp, ii, jj = 0, ip;
    double temp, wmax, t, alpha, asave, dummy, unorm, ztest, cc;


    /* Check the parameters and data */
    if (m <= 0 || n <= 0 || a == NULL || b == NULL || x == NULL)
        return(2);
    /* Allocate memory for working space, if required */
    if (wp != NULL)
        w = wp;
    else
        w = (double*)calloc(n, sizeof(double));
    if (zzp != NULL)
        zz = zzp;
    else
        zz = (double*)calloc(m, sizeof(double));
    if (indexp != NULL)
        index = indexp;
    else
        index = (int*)calloc(n, sizeof(int));
    if (w == NULL || zz == NULL || index == NULL)
        return(2);

    /* Initialize the arrays INDEX[] and X[] */
    for (k = 0; k < n; k++)
    {
        x[k] = 0.;
        index[k] = k;
    }
    iz2 = n - 1;
    iz1 = 0;
    nsetp = 0;
    npp1 = 0;

    /* Main loop; quit if all coeffs are already in the solution or */
    /* if M cols of A have been triangularized */
    iter = 0;
    itmax = n * 3;
    while (iz1 <= iz2 && nsetp < m)
    {
        /* Compute components of the dual (negative gradient) vector W[] */
        for (iz = iz1; iz <= iz2; iz++)
        {
            j = index[iz];
            sm = 0.;
            for (l = npp1; l < m; l++)
                sm += a[j*m+l] * b[l];
            w[j] = sm;
        }

        while (1)
        {
            /* Find largest positive W[j] */
            for (wmax = 0., iz = iz1; iz <= iz2; iz++)
            {
                j = index[iz];
                if (w[j] > wmax)
                {
                    wmax = w[j];
                    izmax = iz;
                }
            }

            /* Terminate if wmax<=0.; */
            /* it indicates satisfaction of the Kuhn-Tucker conditions */
            if (wmax <= 0.0)
                break;
            iz = izmax;
            j = index[iz];

            /* The sign of W[j] is ok for j to be moved to set P. */
            /* Begin the transformation and check new diagonal element to avoid */
            /* near linear dependence. */
            asave = a[j*m+npp1];
            _nnls_h12(1, npp1, npp1 + 1, m, &a[j*m+0], 1, &up, &dummy, 1, 1, 0);
            unorm = 0.;
            if (nsetp != 0)
                for (l = 0; l < nsetp; l++)
                {
                    d1 = a[j*m+l];
                    unorm += d1 * d1;
                }
            unorm = sqrt(unorm);
            d2 = unorm + (d1 = a[j*m+npp1], fabs(d1)) * 0.01;
            if ((d2 - unorm) > 0.)
            {
                /* Col j is sufficiently independent. Copy B into ZZ, update ZZ */
                /* and solve for ztest ( = proposed new value for X[j] ) */
                for (l = 0; l < m; l++)
                    zz[l] = b[l];
                _nnls_h12(2, npp1, npp1 + 1, m, &a[j*m+0], 1, &up, zz, 1, 1, 1);
                ztest = zz[npp1] / a[j*m+npp1];
                /* See if ztest is positive */
                if (ztest > 0.)
                    break;
            }

            /* Reject j as a candidate to be moved from set Z to set P. Restore */
            /* A[npp1,j], set W[j]=0., and loop back to test dual coeffs again */
            a[j*m+npp1] = asave;
            w[j] = 0.;
        } /* while(1) */
        if (wmax <= 0.0)
            break;

        /* Index j=INDEX[iz] has been selected to be moved from set Z to set P. */
        /* Update B and indices, apply householder transformations to cols in */
        /* new set Z, zero subdiagonal elts in col j, set W[j]=0. */
        for (l = 0; l < m; ++l)
            b[l] = zz[l];
        index[iz] = index[iz1];
        index[iz1] = j;
        iz1++;
        nsetp = npp1 + 1;
        npp1++;
        if (iz1 <= iz2)
            for (jz = iz1; jz <= iz2; jz++)
            {
                jj = index[jz];
                _nnls_h12(2, nsetp - 1, npp1, m, &a[j*m+0], 1, &up,
                          &a[jj*m+0], 1, m, 1);
            }
        if (nsetp != m)
            for (l = npp1; l < m; l++)
                a[j*m+l] = 0.;
        w[j] = 0.;
        /* Solve the triangular system; store the solution temporarily in Z[] */
        for (l = 0; l < nsetp; l++)
        {
            ip = nsetp - (l + 1);
            if (l != 0)
                for (ii = 0; ii <= ip; ii++)
                    zz[ii] -= a[jj*m+ii] * zz[ip+1];
            jj = index[ip];
            zz[ip] /= a[jj*m+ip];
        }

        /* Secondary loop begins here */
        while (++iter < itmax)
        {
            /* See if all new constrained coeffs are feasible; if not, compute alpha */
            for (alpha = 2.0, ip = 0; ip < nsetp; ip++)
            {
                l = index[ip];
                if (zz[ip] <= 0.)
                {
                    t = -x[l] / (zz[ip] - x[l]);
                    if (alpha > t)
                    {
                        alpha = t;
                        jj = ip - 1;
                    }
                }
            }

            /* If all new constrained coeffs are feasible then still alpha==2. */
            /* If so, then exit from the secondary loop to main loop */
            if (alpha == 2.0)
                break;
            /* Use alpha (0.<alpha<1.) to interpolate between old X and new ZZ */
            for (ip = 0; ip < nsetp; ip++)
            {
                l = index[ip];
                x[l] += alpha * (zz[ip] - x[l]);
            }

            /* Modify A and B and the INDEX arrays to move coefficient i */
            /* from set P to set Z. */
            k = index[jj+1];
            pfeas = 1;
            do
            {
                x[k] = 0.;
                if (jj != (nsetp - 1))
                {
                    jj++;
                    for (j = jj + 1; j < nsetp; j++)
                    {
                        ii = index[j];
                        index[j-1] = ii;
                        _nnls_g1(a[ii*m+j-1], a[ii*m+j], &cc, &ss, &a[ii*m+j-1]);
                        for (a[ii*m+j] = 0., l = 0; l < n; l++)
                            if (l != ii)
                            {
                                /* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
                                temp = a[l*m+j-1];
                                a[l*m+j-1] = cc * temp + ss * a[l*m+j];
                                a[l*m+j] = -ss * temp + cc * a[l*m+j];
                            }
                        /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
                        temp = b[j-1];
                        b[j-1] = cc * temp + ss * b[j];
                        b[j] = -ss * temp + cc * b[j];
                    }
                }
                npp1 = nsetp - 1;
                nsetp--;
                iz1--;
                index[iz1] = k;

                /* See if the remaining coeffs in set P are feasible; they should */
                /* be because of the way alpha was determined. If any are */
                /* infeasible it is due to round-off error. Any that are */
                /* nonpositive will be set to zero and moved from set P to set Z */
                for (jj = 0; jj < nsetp; jj++)
                {
                    k = index[jj];
                    if (x[k] <= 0.)
                    {
                        pfeas = 0;
                        break;
                    }
                }
            }
            while (pfeas == 0);

            /* Copy B[] into zz[], then solve again and loop back */
            for (k = 0; k < m; k++)
                zz[k] = b[k];
            for (l = 0; l < nsetp; l++)
            {
                ip = nsetp - (l + 1);
                if (l != 0)
                    for (ii = 0; ii <= ip; ii++)
                        zz[ii] -= a[jj*m+ii] * zz[ip+1];
                jj = index[ip];
                zz[ip] /= a[jj*m+ip];
            }
        } /* end of secondary loop */
        if (iter > itmax)
        {
            ret = 1;
            break;
        }
        for (ip = 0; ip < nsetp; ip++)
        {
            k = index[ip];
            x[k] = zz[ip];
        }
    } /* end of main loop */
    /* Compute the norm of the final residual vector */
    sm = 0.;
    if (npp1 < m)
        for (k = npp1; k < m; k++)
            sm += (b[k] * b[k]);
    else
        for (j = 0; j < n; j++)
            w[j] = 0.;
    *rnorm = sqrt(sm);
    /* Free working space, if it was allocated here */
    if (wp == NULL)
        free(w);
    if (zzp == NULL)
        free(zz);
    if (indexp == NULL)
        free(index);
    return(ret);
} /* nnls_ */

nnlsWght()

int nnlsWght	(	int	N,
		int	M,
		double *	A,
		double *	b,
		double *	weight
	)

Definition at line 1455 of file numerical_recipes.cpp.

{
    int n, m;
    double *w;

    /* Check the arguments */
    if (N < 1 || M < 1 || A == NULL || b == NULL || weight == NULL)
        return(1);

    /* Allocate memory */
    w = (double*)malloc(M * sizeof(double));
    if (w == NULL)
        return(2);

    /* Check that weights are not zero and get the square roots of them to w[] */
    for (m = 0; m < M; m++)
    {
        if (weight[m] <= 1.0e-20)
            w[m] = 0.0;
        else
            w[m] = sqrt(weight[m]);
    }

    /* Multiply rows of matrix A and elements of vector b with weights*/
    for (m = 0; m < M; m++)
    {
        for (n = 0; n < N; n++)
        {
            A[n*M+m] *= w[m];
        }
        b[m] *= w[m];
    }

    free(w);
    return(0);
}

nrerror()

void nrerror

(

const char

error_text[]

)

Definition at line 35 of file numerical_recipes.cpp.

{
    fprintf(stderr, "Numerical Recipes run-time error...\n");
    fprintf(stderr, "%s\n", error_text);
    fprintf(stderr, "...now exiting to system...\n");
    exit(1);
}

polint()

void polint	(	double *	xa,
		double *	ya,
		int	n,
		double	x,
		double &	y,
		double &	dy
	)

Definition at line 9390 of file numerical_recipes.cpp.

{
    int i, m, ns = 1;
    double den, dif, dift, ho, hp, w;
    dif = fabs(x - xa[1]);
    std::vector<double> buffer(2*n);
    auto *c= buffer.data()-1;
    auto *d= c + n;
    for (i = 1;i <= n;i++)
    {
        if ((dift = fabs(x - xa[i])) < dif)
        {
            ns = i;
            dif = dift;
        }
        c[i] = ya[i];
        d[i] = ya[i];
    }
    y = ya[ns--];
    for (m = 1;m < n;m++)
    {
        for (i = 1;i <= n - m;i++)
        {
            ho = xa[i] - x;
            hp = xa[i+m] - x;
            w = c[i+1] - d[i];
            if ((den = ho - hp) == 0.0)
            {
                nrerror("error in routine polint\n");
            }
            den = w / den;
            d[i] = hp * den;
            c[i] = ho * den;
        }
        y += (dy = (2 * ns < (n - m) ? c[ns+1] : d[ns--]));
    }
}

powell()

void powell	(	double *	p,
		double *	xi,
		int	n,
		double	ftol,
		int &	iter,
		double &	fret,
		double(*)(double *, void *)	func,
		void *	prm,
		bool	show
	)

Definition at line 877 of file numerical_recipes.cpp.

{
    int i, ibig, j;
    double t, fptt, fp, del;
    std::vector<double> buffer(3*n);
    auto *pt= buffer.data()-1;
    auto *ptt= pt + n;
    auto *xit= ptt + n;
    bool   different_from_0;

    fret = (*func)(p,prm);
    for (j = 1;j <= n;j++)
        pt[j] = p[j];

    for (iter = 1;;(iter)++)
    {
        /* By coss ----- */
        if (show)
        {
            std::cout << iter << " (" << p[1];
            for (int co = 2; co <= n; co++)
                std::cout << "," << p[co];
            std::cout << ")--->" << fret << std::endl;
        }
        /* ------------- */

        fp = fret;
        ibig = 0;
        del = 0.0;
        for (i = 1;i <= n;i++)
        {
            different_from_0 = false; // CO
            for (j = 1;j <= n;j++)
            {
                xit[j] = xi[j*n+i];
                if (xit[j] != 0)
                    different_from_0 = true;
            }
            if (different_from_0)
            {
                fptt = fret;
                linmin(p, xit, n, fret, func, prm);
                if (fabs(fptt - fret) > del)
                {
                    del = fabs(fptt - fret);
                    ibig = i;
                }
                /* By coss ----- */
                if (show)
                {
                    std::cout << "   (";
                    if (i == 1)
                        std::cout << "***";
                    std::cout << p[1];
                    for (int co = 2; co <= n; co++)
                    {
                        std::cout << ",";
                        if (co == i)
                            std::cout << "***";
                        std::cout << p[co];
                    }
                    std::cout << ")--->" << fret << std::endl;
                }
                /* ------------- */
            }
        }
        if (2.0*fabs(fp - fret) <= ftol*(fabs(fp) + fabs(fret)) || n==1) return;
        if (iter == ITMAX)
            nrerror("Too many iterations in routine POWELL");
        for (j = 1;j <= n;j++)
        {
            ptt[j] = 2.0 * p[j] - pt[j];
            xit[j] = p[j] - pt[j];
            pt[j] = p[j];
        }
        fptt = (*func)(ptt,prm);
        if (fptt < fp)
        {
#define SQR(a) ((a)*(a))
            t = 2.0 * (fp - 2.0 * fret + fptt) * SQR(fp - fret - del) - del * SQR(fp - fptt);
            if (t < 0.0)
            {
                linmin(p, xit, n, fret, func, prm);
                for (j = 1;j <= n;j++)
                    xi[j*n+ibig] = xit[j];
            }
        }
    }
}

probks()

double probks

(

double

alam

)

Definition at line 184 of file numerical_recipes.cpp.

{
    int j;
    double a2, fac=2.0, sum=0.0, term, termbf=0.0;
    double EPS1=0.001, EPS2=1.0e-8;

    a2 = -2.0 * alam * alam;
    for (j = 1; j<= 100; j++)
    {
        term = fac * exp(a2*j*j);
        sum += term;
        if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum)
            return sum;
        fac = -fac;
        termbf = fabs(term);
    }
    return 1.0;

}

pwt()

void pwt	(	double	a[],
		unsigned long	n,
		int	isign
	)

pwtset()

void pwtset

(

int

n

)

ran1()

double ran1

(

int *

idum

)

Definition at line 58 of file numerical_recipes.cpp.

{
    static long ix1, ix2, ix3;
    static double r[98];
    double temp;
    static int iff = 0;
    int j;

    if (*idum < 0 || iff == 0)
    {
        iff = 1;
        ix1 = (IC1 - (*idum)) % M1;
        ix1 = (IA1 * ix1 + IC1) % M1;
        ix2 = ix1 % M2;
        ix1 = (IA1 * ix1 + IC1) % M1;
        ix3 = ix1 % M3;
        for (j = 1;j <= 97;j++)
        {
            ix1 = (IA1 * ix1 + IC1) % M1;
            ix2 = (IA2 * ix2 + IC2) % M2;
            r[j] = (ix1 + ix2 * RM2) * RM1;
        }
        *idum = 1;
    }
    ix1 = (IA1 * ix1 + IC1) % M1;
    ix2 = (IA2 * ix2 + IC2) % M2;
    ix3 = (IA3 * ix3 + IC3) % M3;
    j = 1 + ((97 * ix3) / M3);
    if (j > 97 || j < 1)
        nrerror("RAN1: This cannot happen.");
    temp = r[j];
    r[j] = (ix1 + ix2 * RM2) * RM1;
    return temp;
}

realft() [1/2]

void realft	(	double	data[],
		int	n,
		int	isign
	)

realft() [2/2]

void realft	(	double *	data,
		int	n,
		int	isign
	)

sinc()

double sinc ( double  x )

inline

Definition at line 81 of file numerical_recipes.h.

{
    if (fabs(x)<0.0001)
        return 1;
    else
    {
        double arg=PI*x;
        return sin(arg)/arg;
    }
}

svbksb()

void svbksb	(	double *	u,
		double *	w,
		double *	v,
		int	m,
		int	n,
		double *	b,
		double *	x
	)

Definition at line 1791 of file numerical_recipes.cpp.

{
    int jj, j, i;
    double s;

    std::vector<double> buffer(n);
    auto *tmp= buffer.data()-1;
    for (j = 1;j <= n;j++)
    {
        s = 0.0;
        if (w[j])
        {
            for (i = 1;i <= m;i++)
                s += u[i*n+j] * b[i];
            s /= w[j];
        }
        tmp[j] = s;
    }
    for (j = 1;j <= n;j++)
    {
        s = 0.0;
        for (jj = 1;jj <= n;jj++)
            s += v[j*n+jj] * tmp[jj];
        x[j] = s;
    }
}

svdcmp()

void svdcmp	(	double *	a,
		int	m,
		int	n,
		double *	w,
		double *	v
	)

Definition at line 1507 of file numerical_recipes.cpp.

{
    double Norm, Scale;
    double c, f, g, h, s;
    double x, y, z;
    long i, its, j, jj, k, l = 0L, nm = 0L;
    bool    Flag;
    int     MaxIterations = SVDMAXITER;

    std::vector<double> buffer(Columns*Columns);
    auto *rv1= buffer.data();
    g = Scale = Norm = 0.0;
    for (i = 0L; (i < Columns); i++)
    {
        l = i + 1L;
        rv1[i] = Scale * g;
        g = s = Scale = 0.0;
        if (i < Lines)
        {
            for (k = i; (k < Lines); k++)
            {
                Scale += fabs(U[k * Columns + i]);
            }
            if (Scale != 0.0)
            {
                for (k = i; (k < Lines); k++)
                {
                    U[k * Columns + i] /= Scale;
                    s += U[k * Columns + i] * U[k * Columns + i];
                }
                f = U[i * Columns + i];
                g = (0.0 <= f) ? (-sqrt(s)) : (sqrt(s));
                h = f * g - s;
                U[i * Columns + i] = f - g;
                for (j = l; (j < Columns); j++)
                {
                    for (s = 0.0, k = i; (k < Lines); k++)
                    {
                        s += U[k * Columns + i] * U[k * Columns + j];
                    }
                    f = s / h;
                    for (k = i; (k < Lines); k++)
                    {
                        U[k * Columns + j] += f * U[k * Columns + i];
                    }
                }
                for (k = i; (k < Lines); k++)
                {
                    U[k * Columns + i] *= Scale;
                }
            }
        }
        W[i] = Scale * g;
        g = s = Scale = 0.0;
        if ((i < Lines) && (i != (Columns - 1L)))
        {
            for (k = l; (k < Columns); k++)
            {
                Scale += fabs(U[i * Columns + k]);
            }
            if (Scale != 0.0)
            {
                for (k = l; (k < Columns); k++)
                {
                    U[i * Columns + k] /= Scale;
                    s += U[i * Columns + k] * U[i * Columns + k];
                }
                f = U[i * Columns + l];
                g = (0.0 <= f) ? (-sqrt(s)) : (sqrt(s));
                h = f * g - s;
                U[i * Columns + l] = f - g;
                for (k = l; (k < Columns); k++)
                {
                    rv1[k] = U[i * Columns + k] / h;
                }
                for (j = l; (j < Lines); j++)
                {
                    for (s = 0.0, k = l; (k < Columns); k++)
                    {
                        s += U[j * Columns + k] * U[i * Columns + k];
                    }
                    for (k = l; (k < Columns); k++)
                    {
                        U[j * Columns + k] += s * rv1[k];
                    }
                }
                for (k = l; (k < Columns); k++)
                {
                    U[i * Columns + k] *= Scale;
                }
            }
        }
        Norm = ((fabs(W[i]) + fabs(rv1[i])) < Norm) ? (Norm) : (fabs(W[i]) + fabs(rv1[i]));
    }
    for (i = Columns - 1L; (0L <= i); i--)
    {
        if (i < (Columns - 1L))
        {
            if (g != 0.0)
            {
                for (j = l; (j < Columns); j++)
                {
                    V[j * Columns + i] = U[i * Columns + j] / (U[i * Columns + l] * g);
                }
                for (j = l; (j < Columns); j++)
                {
                    for (s = 0.0, k = l; (k < Columns); k++)
                    {
                        s += U[i * Columns + k] * V[k * Columns + j];
                    }
                    for (k = l; (k < Columns); k++)
                    {
                        if (s != 0.0)
                        {
                            V[k * Columns + j] += s * V[k * Columns + i];
                        }
                    }
                }
            }
            for (j = l; (j < Columns); j++)
            {
                V[i * Columns + j] = V[j * Columns + i] = 0.0;
            }
        }
        V[i * Columns + i] = 1.0;
        g = rv1[i];
        l = i;
    }
    for (i = (Lines < Columns) ? (Lines - 1L) : (Columns - 1L); (0L <= i); i--)
    {
        l = i + 1L;
        g = W[i];
        for (j = l; (j < Columns); j++)
        {
            U[i * Columns + j] = 0.0;
        }
        if (g != 0.0)
        {
            g = 1.0 / g;
            for (j = l; (j < Columns); j++)
            {
                for (s = 0.0, k = l; (k < Lines); k++)
                {
                    s += U[k * Columns + i] * U[k * Columns + j];
                }
                f = s * g / U[i * Columns + i];
                for (k = i; (k < Lines); k++)
                {
                    if (f != 0.0)
                    {
                        U[k * Columns + j] += f * U[k * Columns + i];
                    }
                }
            }
            for (j = i; (j < Lines); j++)
            {
                U[j * Columns + i] *= g;
            }
        }
        else
        {
            for (j = i; (j < Lines); j++)
            {
                U[j * Columns + i] = 0.0;
            }
        }
        U[i * Columns + i] += 1.0;
    }
    for (k = Columns - 1L; (0L <= k); k--)
    {
        for (its = 1L; (its <= MaxIterations); its++)
        {
            Flag = true;
            for (l = k; (0L <= l); l--)
            {
                nm = l - 1L;
                if ((fabs(rv1[l]) + Norm) == Norm)
                {
                    Flag = false;
                    break;
                }
                if ((fabs(W[nm]) + Norm) == Norm)
                {
                    break;
                }
            }
            if (Flag)
            {
                c = 0.0;
                s = 1.0;
                for (i = l; (i <= k); i++)
                {
                    f = s * rv1[i];
                    rv1[i] *= c;
                    if ((fabs(f) + Norm) == Norm)
                    {
                        break;
                    }
                    g = W[i];
                    h = Pythag(f, g);
                    W[i] = h;
                    h = 1.0 / h;
                    c = g * h;
                    s = -f * h;
                    for (j = 0L; (j < Lines); j++)
                    {
                        y = U[j * Columns + nm];
                        z = U[j * Columns + i];
                        U[j * Columns + nm] = y * c + z * s;
                        U[j * Columns + i] = z * c - y * s;
                    }
                }
            }
            z = W[k];
            if (l == k)
            {
                if (z < 0.0)
                {
                    W[k] = -z;
                    for (j = 0L; (j < Columns); j++)
                    {
                        V[j * Columns + k] = -V[j * Columns + k];
                    }
                }
                break;
            }
            if (its == MaxIterations) return;
            x = W[l];
            nm = k - 1L;
            y = W[nm];
            g = rv1[nm];
            h = rv1[k];
            f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
            g = Pythag(f, 1.0);
            f = ((x - z) * (x + z) + h * ((y / (f + ((0.0 <= f) ? (fabs(g))
                                                : (-fabs(g))))) - h)) / x;
            c = s = 1.0;
            for (j = l; (j <= nm); j++)
            {
                i = j + 1L;
                g = rv1[i];
                y = W[i];
                h = s * g;
                g = c * g;
                z = Pythag(f, h);
                rv1[j] = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g * c - x * s;
                h = y * s;
                y *= c;
                for (jj = 0L; (jj < Columns); jj++)
                {
                    x = V[jj * Columns + j];
                    z = V[jj * Columns + i];
                    V[jj * Columns + j] = x * c + z * s;
                    V[jj * Columns + i] = z * c - x * s;
                }
                z = Pythag(f, h);
                W[j] = z;
                if (z != 0.0)
                {
                    z = 1.0 / z;
                    c = f * z;
                    s = h * z;
                }
                f = c * g + s * y;
                x = c * y - s * g;
                for (jj = 0L; (jj < Lines); jj++)
                {
                    y = U[jj * Columns + j];
                    z = U[jj * Columns + i];
                    U[jj * Columns + j] = y * c + z * s;
                    U[jj * Columns + i] = z * c - y * s;
                }
            }
            rv1[l] = 0.0;
            rv1[k] = f;
            W[k] = x;
        }
    }
}

tdev()

double tdev	(	double	nu,
		int *	idum
	)

Definition at line 136 of file numerical_recipes.cpp.

{
    static int iset = 0;
    static double gset;
    double fac, r, v1, v2;

    if (iset == 0)
    {
        do
        {
            v1 = 2.0 * ran1(idum) - 1.0;
            v2 = 2.0 * ran1(idum) - 1.0;
            r = v1 * v1 + v2 * v2;
        }
        while (r >= 1.0);
        fac = sqrt(nu*(pow(r,-2.0/nu) -1.0)/r);
        gset = v1 * fac;
        iset = 1;
        return v2*fac;
    }
    else
    {
        iset = 0;
        return gset;
    }
}

twofft()

void twofft	(	double *	data1,
		double *	data2,
		double *	fft1,
		double *	fft2,
		int	n
	)

wtn()

void wtn	(	double	a[],
		unsigned long	nn[],
		int	ndim,
		int	isign,
		void(*)(double [], unsigned long, int)	wtstep
	)