This is gsl-ref.info, produced by makeinfo version 5.1 from gsl-ref.texi. Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 The GSL Team. Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections and no cover texts. A copy of the license is included in the section entitled "GNU Free Documentation License". INFO-DIR-SECTION Software libraries START-INFO-DIR-ENTRY * gsl-ref: (gsl-ref). GNU Scientific Library - Reference END-INFO-DIR-ENTRY  File: gsl-ref.info, Node: The Beta Distribution, Next: The Logistic Distribution, Prev: The t-distribution, Up: Random Number Distributions 20.20 The Beta Distribution =========================== -- Function: double gsl_ran_beta (const gsl_rng * R, double A, double B) This function returns a random variate from the beta distribution. The distribution function is, p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx for 0 <= x <= 1. -- Function: double gsl_ran_beta_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a beta distribution with parameters A and B, using the formula given above. -- Function: double gsl_cdf_beta_P (double X, double A, double B) -- Function: double gsl_cdf_beta_Q (double X, double A, double B) -- Function: double gsl_cdf_beta_Pinv (double P, double A, double B) -- Function: double gsl_cdf_beta_Qinv (double Q, double A, double B) These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the beta distribution with parameters A and B.  File: gsl-ref.info, Node: The Logistic Distribution, Next: The Pareto Distribution, Prev: The Beta Distribution, Up: Random Number Distributions 20.21 The Logistic Distribution =============================== -- Function: double gsl_ran_logistic (const gsl_rng * R, double A) This function returns a random variate from the logistic distribution. The distribution function is, p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx for -\infty < x < +\infty. -- Function: double gsl_ran_logistic_pdf (double X, double A) This function computes the probability density p(x) at X for a logistic distribution with scale parameter A, using the formula given above. -- Function: double gsl_cdf_logistic_P (double X, double A) -- Function: double gsl_cdf_logistic_Q (double X, double A) -- Function: double gsl_cdf_logistic_Pinv (double P, double A) -- Function: double gsl_cdf_logistic_Qinv (double Q, double A) These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the logistic distribution with scale parameter A.  File: gsl-ref.info, Node: The Pareto Distribution, Next: Spherical Vector Distributions, Prev: The Logistic Distribution, Up: Random Number Distributions 20.22 The Pareto Distribution ============================= -- Function: double gsl_ran_pareto (const gsl_rng * R, double A, double B) This function returns a random variate from the Pareto distribution of order A. The distribution function is, p(x) dx = (a/b) / (x/b)^{a+1} dx for x >= b. -- Function: double gsl_ran_pareto_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a Pareto distribution with exponent A and scale B, using the formula given above. -- Function: double gsl_cdf_pareto_P (double X, double A, double B) -- Function: double gsl_cdf_pareto_Q (double X, double A, double B) -- Function: double gsl_cdf_pareto_Pinv (double P, double A, double B) -- Function: double gsl_cdf_pareto_Qinv (double Q, double A, double B) These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Pareto distribution with exponent A and scale B.  File: gsl-ref.info, Node: Spherical Vector Distributions, Next: The Weibull Distribution, Prev: The Pareto Distribution, Up: Random Number Distributions 20.23 Spherical Vector Distributions ==================================== The spherical distributions generate random vectors, located on a spherical surface. They can be used as random directions, for example in the steps of a random walk. -- Function: void gsl_ran_dir_2d (const gsl_rng * R, double * X, double * Y) -- Function: void gsl_ran_dir_2d_trig_method (const gsl_rng * R, double * X, double * Y) This function returns a random direction vector v = (X,Y) in two dimensions. The vector is normalized such that |v|^2 = x^2 + y^2 = 1. The obvious way to do this is to take a uniform random number between 0 and 2\pi and let X and Y be the sine and cosine respectively. Two trig functions would have been expensive in the old days, but with modern hardware implementations, this is sometimes the fastest way to go. This is the case for the Pentium (but not the case for the Sun Sparcstation). One can avoid the trig evaluations by choosing X and Y in the interior of a unit circle (choose them at random from the interior of the enclosing square, and then reject those that are outside the unit circle), and then dividing by \sqrt{x^2 + y^2}. A much cleverer approach, attributed to von Neumann (See Knuth, v2, 3rd ed, p140, exercise 23), requires neither trig nor a square root. In this approach, U and V are chosen at random from the interior of a unit circle, and then x=(u^2-v^2)/(u^2+v^2) and y=2uv/(u^2+v^2). -- Function: void gsl_ran_dir_3d (const gsl_rng * R, double * X, double * Y, double * Z) This function returns a random direction vector v = (X,Y,Z) in three dimensions. The vector is normalized such that |v|^2 = x^2 + y^2 + z^2 = 1. The method employed is due to Robert E. Knop (CACM 13, 326 (1970)), and explained in Knuth, v2, 3rd ed, p136. It uses the surprising fact that the distribution projected along any axis is actually uniform (this is only true for 3 dimensions). -- Function: void gsl_ran_dir_nd (const gsl_rng * R, size_t N, double * X) This function returns a random direction vector v = (x_1,x_2,...,x_n) in N dimensions. The vector is normalized such that |v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1. The method uses the fact that a multivariate Gaussian distribution is spherically symmetric. Each component is generated to have a Gaussian distribution, and then the components are normalized. The method is described by Knuth, v2, 3rd ed, p135-136, and attributed to G. W. Brown, Modern Mathematics for the Engineer (1956).  File: gsl-ref.info, Node: The Weibull Distribution, Next: The Type-1 Gumbel Distribution, Prev: Spherical Vector Distributions, Up: Random Number Distributions 20.24 The Weibull Distribution ============================== -- Function: double gsl_ran_weibull (const gsl_rng * R, double A, double B) This function returns a random variate from the Weibull distribution. The distribution function is, p(x) dx = {b \over a^b} x^{b-1} \exp(-(x/a)^b) dx for x >= 0. -- Function: double gsl_ran_weibull_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a Weibull distribution with scale A and exponent B, using the formula given above. -- Function: double gsl_cdf_weibull_P (double X, double A, double B) -- Function: double gsl_cdf_weibull_Q (double X, double A, double B) -- Function: double gsl_cdf_weibull_Pinv (double P, double A, double B) -- Function: double gsl_cdf_weibull_Qinv (double Q, double A, double B) These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Weibull distribution with scale A and exponent B.  File: gsl-ref.info, Node: The Type-1 Gumbel Distribution, Next: The Type-2 Gumbel Distribution, Prev: The Weibull Distribution, Up: Random Number Distributions 20.25 The Type-1 Gumbel Distribution ==================================== -- Function: double gsl_ran_gumbel1 (const gsl_rng * R, double A, double B) This function returns a random variate from the Type-1 Gumbel distribution. The Type-1 Gumbel distribution function is, p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx for -\infty < x < \infty. -- Function: double gsl_ran_gumbel1_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a Type-1 Gumbel distribution with parameters A and B, using the formula given above. -- Function: double gsl_cdf_gumbel1_P (double X, double A, double B) -- Function: double gsl_cdf_gumbel1_Q (double X, double A, double B) -- Function: double gsl_cdf_gumbel1_Pinv (double P, double A, double B) -- Function: double gsl_cdf_gumbel1_Qinv (double Q, double A, double B) These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-1 Gumbel distribution with parameters A and B.  File: gsl-ref.info, Node: The Type-2 Gumbel Distribution, Next: The Dirichlet Distribution, Prev: The Type-1 Gumbel Distribution, Up: Random Number Distributions 20.26 The Type-2 Gumbel Distribution ==================================== -- Function: double gsl_ran_gumbel2 (const gsl_rng * R, double A, double B) This function returns a random variate from the Type-2 Gumbel distribution. The Type-2 Gumbel distribution function is, p(x) dx = a b x^{-a-1} \exp(-b x^{-a}) dx for 0 < x < \infty. -- Function: double gsl_ran_gumbel2_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a Type-2 Gumbel distribution with parameters A and B, using the formula given above. -- Function: double gsl_cdf_gumbel2_P (double X, double A, double B) -- Function: double gsl_cdf_gumbel2_Q (double X, double A, double B) -- Function: double gsl_cdf_gumbel2_Pinv (double P, double A, double B) -- Function: double gsl_cdf_gumbel2_Qinv (double Q, double A, double B) These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-2 Gumbel distribution with parameters A and B.  File: gsl-ref.info, Node: The Dirichlet Distribution, Next: General Discrete Distributions, Prev: The Type-2 Gumbel Distribution, Up: Random Number Distributions 20.27 The Dirichlet Distribution ================================ -- Function: void gsl_ran_dirichlet (const gsl_rng * R, size_t K, const double ALPHA[], double THETA[]) This function returns an array of K random variates from a Dirichlet distribution of order K-1. The distribution function is p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K = (1/Z) \prod_{i=1}^K \theta_i^{\alpha_i - 1} \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 ... d\theta_K for theta_i >= 0 and alpha_i > 0. The delta function ensures that \sum \theta_i = 1. The normalization factor Z is Z = {\prod_{i=1}^K \Gamma(\alpha_i)} / {\Gamma( \sum_{i=1}^K \alpha_i)} The random variates are generated by sampling K values from gamma distributions with parameters a=alpha_i, b=1, and renormalizing. See A.M. Law, W.D. Kelton, 'Simulation Modeling and Analysis' (1991). -- Function: double gsl_ran_dirichlet_pdf (size_t K, const double ALPHA[], const double THETA[]) This function computes the probability density p(\theta_1, ... , \theta_K) at THETA[K] for a Dirichlet distribution with parameters ALPHA[K], using the formula given above. -- Function: double gsl_ran_dirichlet_lnpdf (size_t K, const double ALPHA[], const double THETA[]) This function computes the logarithm of the probability density p(\theta_1, ... , \theta_K) for a Dirichlet distribution with parameters ALPHA[K].  File: gsl-ref.info, Node: General Discrete Distributions, Next: The Poisson Distribution, Prev: The Dirichlet Distribution, Up: Random Number Distributions 20.28 General Discrete Distributions ==================================== Given K discrete events with different probabilities P[k], produce a random value k consistent with its probability. The obvious way to do this is to preprocess the probability list by generating a cumulative probability array with K+1 elements: C[0] = 0 C[k+1] = C[k]+P[k]. Note that this construction produces C[K]=1. Now choose a uniform deviate u between 0 and 1, and find the value of k such that C[k] <= u < C[k+1]. Although this in principle requires of order \log K steps per random number generation, they are fast steps, and if you use something like \lfloor uK \rfloor as a starting point, you can often do pretty well. But faster methods have been devised. Again, the idea is to preprocess the probability list, and save the result in some form of lookup table; then the individual calls for a random discrete event can go rapidly. An approach invented by G. Marsaglia (Generating discrete random variables in a computer, Comm ACM 6, 37-38 (1963)) is very clever, and readers interested in examples of good algorithm design are directed to this short and well-written paper. Unfortunately, for large K, Marsaglia's lookup table can be quite large. A much better approach is due to Alastair J. Walker (An efficient method for generating discrete random variables with general distributions, ACM Trans on Mathematical Software 3, 253-256 (1977); see also Knuth, v2, 3rd ed, p120-121,139). This requires two lookup tables, one floating point and one integer, but both only of size K. After preprocessing, the random numbers are generated in O(1) time, even for large K. The preprocessing suggested by Walker requires O(K^2) effort, but that is not actually necessary, and the implementation provided here only takes O(K) effort. In general, more preprocessing leads to faster generation of the individual random numbers, but a diminishing return is reached pretty early. Knuth points out that the optimal preprocessing is combinatorially difficult for large K. This method can be used to speed up some of the discrete random number generators below, such as the binomial distribution. To use it for something like the Poisson Distribution, a modification would have to be made, since it only takes a finite set of K outcomes. -- Function: gsl_ran_discrete_t * gsl_ran_discrete_preproc (size_t K, const double * P) This function returns a pointer to a structure that contains the lookup table for the discrete random number generator. The array P[] contains the probabilities of the discrete events; these array elements must all be positive, but they needn't add up to one (so you can think of them more generally as "weights")--the preprocessor will normalize appropriately. This return value is used as an argument for the 'gsl_ran_discrete' function below. -- Function: size_t gsl_ran_discrete (const gsl_rng * R, const gsl_ran_discrete_t * G) After the preprocessor, above, has been called, you use this function to get the discrete random numbers. -- Function: double gsl_ran_discrete_pdf (size_t K, const gsl_ran_discrete_t * G) Returns the probability P[k] of observing the variable K. Since P[k] is not stored as part of the lookup table, it must be recomputed; this computation takes O(K), so if K is large and you care about the original array P[k] used to create the lookup table, then you should just keep this original array P[k] around. -- Function: void gsl_ran_discrete_free (gsl_ran_discrete_t * G) De-allocates the lookup table pointed to by G.  File: gsl-ref.info, Node: The Poisson Distribution, Next: The Bernoulli Distribution, Prev: General Discrete Distributions, Up: Random Number Distributions 20.29 The Poisson Distribution ============================== -- Function: unsigned int gsl_ran_poisson (const gsl_rng * R, double MU) This function returns a random integer from the Poisson distribution with mean MU. The probability distribution for Poisson variates is, p(k) = {\mu^k \over k!} \exp(-\mu) for k >= 0. -- Function: double gsl_ran_poisson_pdf (unsigned int K, double MU) This function computes the probability p(k) of obtaining K from a Poisson distribution with mean MU, using the formula given above. -- Function: double gsl_cdf_poisson_P (unsigned int K, double MU) -- Function: double gsl_cdf_poisson_Q (unsigned int K, double MU) These functions compute the cumulative distribution functions P(k), Q(k) for the Poisson distribution with parameter MU.  File: gsl-ref.info, Node: The Bernoulli Distribution, Next: The Binomial Distribution, Prev: The Poisson Distribution, Up: Random Number Distributions 20.30 The Bernoulli Distribution ================================ -- Function: unsigned int gsl_ran_bernoulli (const gsl_rng * R, double P) This function returns either 0 or 1, the result of a Bernoulli trial with probability P. The probability distribution for a Bernoulli trial is, p(0) = 1 - p p(1) = p -- Function: double gsl_ran_bernoulli_pdf (unsigned int K, double P) This function computes the probability p(k) of obtaining K from a Bernoulli distribution with probability parameter P, using the formula given above.  File: gsl-ref.info, Node: The Binomial Distribution, Next: The Multinomial Distribution, Prev: The Bernoulli Distribution, Up: Random Number Distributions 20.31 The Binomial Distribution =============================== -- Function: unsigned int gsl_ran_binomial (const gsl_rng * R, double P, unsigned int N) This function returns a random integer from the binomial distribution, the number of successes in N independent trials with probability P. The probability distribution for binomial variates is, p(k) = {n! \over k! (n-k)! } p^k (1-p)^{n-k} for 0 <= k <= n. -- Function: double gsl_ran_binomial_pdf (unsigned int K, double P, unsigned int N) This function computes the probability p(k) of obtaining K from a binomial distribution with parameters P and N, using the formula given above. -- Function: double gsl_cdf_binomial_P (unsigned int K, double P, unsigned int N) -- Function: double gsl_cdf_binomial_Q (unsigned int K, double P, unsigned int N) These functions compute the cumulative distribution functions P(k), Q(k) for the binomial distribution with parameters P and N.  File: gsl-ref.info, Node: The Multinomial Distribution, Next: The Negative Binomial Distribution, Prev: The Binomial Distribution, Up: Random Number Distributions 20.32 The Multinomial Distribution ================================== -- Function: void gsl_ran_multinomial (const gsl_rng * R, size_t K, unsigned int N, const double P[], unsigned int N[]) This function computes a random sample N[] from the multinomial distribution formed by N trials from an underlying distribution P[K]. The distribution function for N[] is, P(n_1, n_2, ..., n_K) = (N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K where (n_1, n_2, ..., n_K) are nonnegative integers with sum_{k=1}^K n_k = N, and (p_1, p_2, ..., p_K) is a probability distribution with \sum p_i = 1. If the array P[K] is not normalized then its entries will be treated as weights and normalized appropriately. The arrays N[] and P[] must both be of length K. Random variates are generated using the conditional binomial method (see C.S. Davis, 'The computer generation of multinomial random variates', Comp. Stat. Data Anal. 16 (1993) 205-217 for details). -- Function: double gsl_ran_multinomial_pdf (size_t K, const double P[], const unsigned int N[]) This function computes the probability P(n_1, n_2, ..., n_K) of sampling N[K] from a multinomial distribution with parameters P[K], using the formula given above. -- Function: double gsl_ran_multinomial_lnpdf (size_t K, const double P[], const unsigned int N[]) This function returns the logarithm of the probability for the multinomial distribution P(n_1, n_2, ..., n_K) with parameters P[K].  File: gsl-ref.info, Node: The Negative Binomial Distribution, Next: The Pascal Distribution, Prev: The Multinomial Distribution, Up: Random Number Distributions 20.33 The Negative Binomial Distribution ======================================== -- Function: unsigned int gsl_ran_negative_binomial (const gsl_rng * R, double P, double N) This function returns a random integer from the negative binomial distribution, the number of failures occurring before N successes in independent trials with probability P of success. The probability distribution for negative binomial variates is, p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k Note that n is not required to be an integer. -- Function: double gsl_ran_negative_binomial_pdf (unsigned int K, double P, double N) This function computes the probability p(k) of obtaining K from a negative binomial distribution with parameters P and N, using the formula given above. -- Function: double gsl_cdf_negative_binomial_P (unsigned int K, double P, double N) -- Function: double gsl_cdf_negative_binomial_Q (unsigned int K, double P, double N) These functions compute the cumulative distribution functions P(k), Q(k) for the negative binomial distribution with parameters P and N.  File: gsl-ref.info, Node: The Pascal Distribution, Next: The Geometric Distribution, Prev: The Negative Binomial Distribution, Up: Random Number Distributions 20.34 The Pascal Distribution ============================= -- Function: unsigned int gsl_ran_pascal (const gsl_rng * R, double P, unsigned int N) This function returns a random integer from the Pascal distribution. The Pascal distribution is simply a negative binomial distribution with an integer value of n. p(k) = {(n + k - 1)! \over k! (n - 1)! } p^n (1-p)^k for k >= 0 -- Function: double gsl_ran_pascal_pdf (unsigned int K, double P, unsigned int N) This function computes the probability p(k) of obtaining K from a Pascal distribution with parameters P and N, using the formula given above. -- Function: double gsl_cdf_pascal_P (unsigned int K, double P, unsigned int N) -- Function: double gsl_cdf_pascal_Q (unsigned int K, double P, unsigned int N) These functions compute the cumulative distribution functions P(k), Q(k) for the Pascal distribution with parameters P and N.  File: gsl-ref.info, Node: The Geometric Distribution, Next: The Hypergeometric Distribution, Prev: The Pascal Distribution, Up: Random Number Distributions 20.35 The Geometric Distribution ================================ -- Function: unsigned int gsl_ran_geometric (const gsl_rng * R, double P) This function returns a random integer from the geometric distribution, the number of independent trials with probability P until the first success. The probability distribution for geometric variates is, p(k) = p (1-p)^(k-1) for k >= 1. Note that the distribution begins with k=1 with this definition. There is another convention in which the exponent k-1 is replaced by k. -- Function: double gsl_ran_geometric_pdf (unsigned int K, double P) This function computes the probability p(k) of obtaining K from a geometric distribution with probability parameter P, using the formula given above. -- Function: double gsl_cdf_geometric_P (unsigned int K, double P) -- Function: double gsl_cdf_geometric_Q (unsigned int K, double P) These functions compute the cumulative distribution functions P(k), Q(k) for the geometric distribution with parameter P.  File: gsl-ref.info, Node: The Hypergeometric Distribution, Next: The Logarithmic Distribution, Prev: The Geometric Distribution, Up: Random Number Distributions 20.36 The Hypergeometric Distribution ===================================== -- Function: unsigned int gsl_ran_hypergeometric (const gsl_rng * R, unsigned int N1, unsigned int N2, unsigned int T) This function returns a random integer from the hypergeometric distribution. The probability distribution for hypergeometric random variates is, p(k) = C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t) where C(a,b) = a!/(b!(a-b)!) and t <= n_1 + n_2. The domain of k is max(0,t-n_2), ..., min(t,n_1). If a population contains n_1 elements of "type 1" and n_2 elements of "type 2" then the hypergeometric distribution gives the probability of obtaining k elements of "type 1" in t samples from the population without replacement. -- Function: double gsl_ran_hypergeometric_pdf (unsigned int K, unsigned int N1, unsigned int N2, unsigned int T) This function computes the probability p(k) of obtaining K from a hypergeometric distribution with parameters N1, N2, T, using the formula given above. -- Function: double gsl_cdf_hypergeometric_P (unsigned int K, unsigned int N1, unsigned int N2, unsigned int T) -- Function: double gsl_cdf_hypergeometric_Q (unsigned int K, unsigned int N1, unsigned int N2, unsigned int T) These functions compute the cumulative distribution functions P(k), Q(k) for the hypergeometric distribution with parameters N1, N2 and T.  File: gsl-ref.info, Node: The Logarithmic Distribution, Next: Shuffling and Sampling, Prev: The Hypergeometric Distribution, Up: Random Number Distributions 20.37 The Logarithmic Distribution ================================== -- Function: unsigned int gsl_ran_logarithmic (const gsl_rng * R, double P) This function returns a random integer from the logarithmic distribution. The probability distribution for logarithmic random variates is, p(k) = {-1 \over \log(1-p)} {(p^k \over k)} for k >= 1. -- Function: double gsl_ran_logarithmic_pdf (unsigned int K, double P) This function computes the probability p(k) of obtaining K from a logarithmic distribution with probability parameter P, using the formula given above.  File: gsl-ref.info, Node: Shuffling and Sampling, Next: Random Number Distribution Examples, Prev: The Logarithmic Distribution, Up: Random Number Distributions 20.38 Shuffling and Sampling ============================ The following functions allow the shuffling and sampling of a set of objects. The algorithms rely on a random number generator as a source of randomness and a poor quality generator can lead to correlations in the output. In particular it is important to avoid generators with a short period. For more information see Knuth, v2, 3rd ed, Section 3.4.2, "Random Sampling and Shuffling". -- Function: void gsl_ran_shuffle (const gsl_rng * R, void * BASE, size_t N, size_t SIZE) This function randomly shuffles the order of N objects, each of size SIZE, stored in the array BASE[0..N-1]. The output of the random number generator R is used to produce the permutation. The algorithm generates all possible n! permutations with equal probability, assuming a perfect source of random numbers. The following code shows how to shuffle the numbers from 0 to 51, int a[52]; for (i = 0; i < 52; i++) { a[i] = i; } gsl_ran_shuffle (r, a, 52, sizeof (int)); -- Function: int gsl_ran_choose (const gsl_rng * R, void * DEST, size_t K, void * SRC, size_t N, size_t SIZE) This function fills the array DEST[k] with K objects taken randomly from the N elements of the array SRC[0..N-1]. The objects are each of size SIZE. The output of the random number generator R is used to make the selection. The algorithm ensures all possible samples are equally likely, assuming a perfect source of randomness. The objects are sampled _without_ replacement, thus each object can only appear once in DEST[k]. It is required that K be less than or equal to 'n'. The objects in DEST will be in the same relative order as those in SRC. You will need to call 'gsl_ran_shuffle(r, dest, n, size)' if you want to randomize the order. The following code shows how to select a random sample of three unique numbers from the set 0 to 99, double a[3], b[100]; for (i = 0; i < 100; i++) { b[i] = (double) i; } gsl_ran_choose (r, a, 3, b, 100, sizeof (double)); -- Function: void gsl_ran_sample (const gsl_rng * R, void * DEST, size_t K, void * SRC, size_t N, size_t SIZE) This function is like 'gsl_ran_choose' but samples K items from the original array of N items SRC with replacement, so the same object can appear more than once in the output sequence DEST. There is no requirement that K be less than N in this case.  File: gsl-ref.info, Node: Random Number Distribution Examples, Next: Random Number Distribution References and Further Reading, Prev: Shuffling and Sampling, Up: Random Number Distributions 20.39 Examples ============== The following program demonstrates the use of a random number generator to produce variates from a distribution. It prints 10 samples from the Poisson distribution with a mean of 3. #include #include #include int main (void) { const gsl_rng_type * T; gsl_rng * r; int i, n = 10; double mu = 3.0; /* create a generator chosen by the environment variable GSL_RNG_TYPE */ gsl_rng_env_setup(); T = gsl_rng_default; r = gsl_rng_alloc (T); /* print n random variates chosen from the poisson distribution with mean parameter mu */ for (i = 0; i < n; i++) { unsigned int k = gsl_ran_poisson (r, mu); printf (" %u", k); } printf ("\n"); gsl_rng_free (r); return 0; } If the library and header files are installed under '/usr/local' (the default location) then the program can be compiled with these options, $ gcc -Wall demo.c -lgsl -lgslcblas -lm Here is the output of the program, $ ./a.out 2 5 5 2 1 0 3 4 1 1 The variates depend on the seed used by the generator. The seed for the default generator type 'gsl_rng_default' can be changed with the 'GSL_RNG_SEED' environment variable to produce a different stream of variates, $ GSL_RNG_SEED=123 ./a.out GSL_RNG_SEED=123 4 5 6 3 3 1 4 2 5 5 The following program generates a random walk in two dimensions. #include #include #include int main (void) { int i; double x = 0, y = 0, dx, dy; const gsl_rng_type * T; gsl_rng * r; gsl_rng_env_setup(); T = gsl_rng_default; r = gsl_rng_alloc (T); printf ("%g %g\n", x, y); for (i = 0; i < 10; i++) { gsl_ran_dir_2d (r, &dx, &dy); x += dx; y += dy; printf ("%g %g\n", x, y); } gsl_rng_free (r); return 0; } Here is some output from the program, four 10-step random walks from the origin, The following program computes the upper and lower cumulative distribution functions for the standard normal distribution at x=2. #include #include int main (void) { double P, Q; double x = 2.0; P = gsl_cdf_ugaussian_P (x); printf ("prob(x < %f) = %f\n", x, P); Q = gsl_cdf_ugaussian_Q (x); printf ("prob(x > %f) = %f\n", x, Q); x = gsl_cdf_ugaussian_Pinv (P); printf ("Pinv(%f) = %f\n", P, x); x = gsl_cdf_ugaussian_Qinv (Q); printf ("Qinv(%f) = %f\n", Q, x); return 0; } Here is the output of the program, prob(x < 2.000000) = 0.977250 prob(x > 2.000000) = 0.022750 Pinv(0.977250) = 2.000000 Qinv(0.022750) = 2.000000  File: gsl-ref.info, Node: Random Number Distribution References and Further Reading, Prev: Random Number Distribution Examples, Up: Random Number Distributions 20.40 References and Further Reading ==================================== For an encyclopaedic coverage of the subject readers are advised to consult the book 'Non-Uniform Random Variate Generation' by Luc Devroye. It covers every imaginable distribution and provides hundreds of algorithms. Luc Devroye, 'Non-Uniform Random Variate Generation', Springer-Verlag, ISBN 0-387-96305-7. Available online at . The subject of random variate generation is also reviewed by Knuth, who describes algorithms for all the major distributions. Donald E. Knuth, 'The Art of Computer Programming: Seminumerical Algorithms' (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842. The Particle Data Group provides a short review of techniques for generating distributions of random numbers in the "Monte Carlo" section of its Annual Review of Particle Physics. 'Review of Particle Properties' R.M. Barnett et al., Physical Review D54, 1 (1996) . The Review of Particle Physics is available online in postscript and pdf format. An overview of methods used to compute cumulative distribution functions can be found in 'Statistical Computing' by W.J. Kennedy and J.E. Gentle. Another general reference is 'Elements of Statistical Computing' by R.A. Thisted. William E. Kennedy and James E. Gentle, 'Statistical Computing' (1980), Marcel Dekker, ISBN 0-8247-6898-1. Ronald A. Thisted, 'Elements of Statistical Computing' (1988), Chapman & Hall, ISBN 0-412-01371-1. The cumulative distribution functions for the Gaussian distribution are based on the following papers, 'Rational Chebyshev Approximations Using Linear Equations', W.J. Cody, W. Fraser, J.F. Hart. Numerische Mathematik 12, 242-251 (1968). 'Rational Chebyshev Approximations for the Error Function', W.J. Cody. Mathematics of Computation 23, n107, 631-637 (July 1969).  File: gsl-ref.info, Node: Statistics, Next: Histograms, Prev: Random Number Distributions, Up: Top 21 Statistics ************* This chapter describes the statistical functions in the library. The basic statistical functions include routines to compute the mean, variance and standard deviation. More advanced functions allow you to calculate absolute deviations, skewness, and kurtosis as well as the median and arbitrary percentiles. The algorithms use recurrence relations to compute average quantities in a stable way, without large intermediate values that might overflow. The functions are available in versions for datasets in the standard floating-point and integer types. The versions for double precision floating-point data have the prefix 'gsl_stats' and are declared in the header file 'gsl_statistics_double.h'. The versions for integer data have the prefix 'gsl_stats_int' and are declared in the header file 'gsl_statistics_int.h'. All the functions operate on C arrays with a STRIDE parameter specifying the spacing between elements. * Menu: * Mean and standard deviation and variance:: * Absolute deviation:: * Higher moments (skewness and kurtosis):: * Autocorrelation:: * Covariance:: * Correlation:: * Weighted Samples:: * Maximum and Minimum values:: * Median and Percentiles:: * Example statistical programs:: * Statistics References and Further Reading::  File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Absolute deviation, Up: Statistics 21.1 Mean, Standard Deviation and Variance ========================================== -- Function: double gsl_stats_mean (const double DATA[], size_t STRIDE, size_t N) This function returns the arithmetic mean of DATA, a dataset of length N with stride STRIDE. The arithmetic mean, or "sample mean", is denoted by \Hat\mu and defined as, \Hat\mu = (1/N) \sum x_i where x_i are the elements of the dataset DATA. For samples drawn from a gaussian distribution the variance of \Hat\mu is \sigma^2 / N. -- Function: double gsl_stats_variance (const double DATA[], size_t STRIDE, size_t N) This function returns the estimated, or "sample", variance of DATA, a dataset of length N with stride STRIDE. The estimated variance is denoted by \Hat\sigma^2 and is defined by, \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2 where x_i are the elements of the dataset DATA. Note that the normalization factor of 1/(N-1) results from the derivation of \Hat\sigma^2 as an unbiased estimator of the population variance \sigma^2. For samples drawn from a Gaussian distribution the variance of \Hat\sigma^2 itself is 2 \sigma^4 / N. This function computes the mean via a call to 'gsl_stats_mean'. If you have already computed the mean then you can pass it directly to 'gsl_stats_variance_m'. -- Function: double gsl_stats_variance_m (const double DATA[], size_t STRIDE, size_t N, double MEAN) This function returns the sample variance of DATA relative to the given value of MEAN. The function is computed with \Hat\mu replaced by the value of MEAN that you supply, \Hat\sigma^2 = (1/(N-1)) \sum (x_i - mean)^2 -- Function: double gsl_stats_sd (const double DATA[], size_t STRIDE, size_t N) -- Function: double gsl_stats_sd_m (const double DATA[], size_t STRIDE, size_t N, double MEAN) The standard deviation is defined as the square root of the variance. These functions return the square root of the corresponding variance functions above. -- Function: double gsl_stats_tss (const double DATA[], size_t STRIDE, size_t N) -- Function: double gsl_stats_tss_m (const double DATA[], size_t STRIDE, size_t N, double MEAN) These functions return the total sum of squares (TSS) of DATA about the mean. For 'gsl_stats_tss_m' the user-supplied value of MEAN is used, and for 'gsl_stats_tss' it is computed using 'gsl_stats_mean'. TSS = \sum (x_i - mean)^2 -- Function: double gsl_stats_variance_with_fixed_mean (const double DATA[], size_t STRIDE, size_t N, double MEAN) This function computes an unbiased estimate of the variance of DATA when the population mean MEAN of the underlying distribution is known _a priori_. In this case the estimator for the variance uses the factor 1/N and the sample mean \Hat\mu is replaced by the known population mean \mu, \Hat\sigma^2 = (1/N) \sum (x_i - \mu)^2 -- Function: double gsl_stats_sd_with_fixed_mean (const double DATA[], size_t STRIDE, size_t N, double MEAN) This function calculates the standard deviation of DATA for a fixed population mean MEAN. The result is the square root of the corresponding variance function.  File: gsl-ref.info, Node: Absolute deviation, Next: Higher moments (skewness and kurtosis), Prev: Mean and standard deviation and variance, Up: Statistics 21.2 Absolute deviation ======================= -- Function: double gsl_stats_absdev (const double DATA[], size_t STRIDE, size_t N) This function computes the absolute deviation from the mean of DATA, a dataset of length N with stride STRIDE. The absolute deviation from the mean is defined as, absdev = (1/N) \sum |x_i - \Hat\mu| where x_i are the elements of the dataset DATA. The absolute deviation from the mean provides a more robust measure of the width of a distribution than the variance. This function computes the mean of DATA via a call to 'gsl_stats_mean'. -- Function: double gsl_stats_absdev_m (const double DATA[], size_t STRIDE, size_t N, double MEAN) This function computes the absolute deviation of the dataset DATA relative to the given value of MEAN, absdev = (1/N) \sum |x_i - mean| This function is useful if you have already computed the mean of DATA (and want to avoid recomputing it), or wish to calculate the absolute deviation relative to another value (such as zero, or the median).  File: gsl-ref.info, Node: Higher moments (skewness and kurtosis), Next: Autocorrelation, Prev: Absolute deviation, Up: Statistics 21.3 Higher moments (skewness and kurtosis) =========================================== -- Function: double gsl_stats_skew (const double DATA[], size_t STRIDE, size_t N) This function computes the skewness of DATA, a dataset of length N with stride STRIDE. The skewness is defined as, skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3 where x_i are the elements of the dataset DATA. The skewness measures the asymmetry of the tails of a distribution. The function computes the mean and estimated standard deviation of DATA via calls to 'gsl_stats_mean' and 'gsl_stats_sd'. -- Function: double gsl_stats_skew_m_sd (const double DATA[], size_t STRIDE, size_t N, double MEAN, double SD) This function computes the skewness of the dataset DATA using the given values of the mean MEAN and standard deviation SD, skew = (1/N) \sum ((x_i - mean)/sd)^3 These functions are useful if you have already computed the mean and standard deviation of DATA and want to avoid recomputing them. -- Function: double gsl_stats_kurtosis (const double DATA[], size_t STRIDE, size_t N) This function computes the kurtosis of DATA, a dataset of length N with stride STRIDE. The kurtosis is defined as, kurtosis = ((1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^4) - 3 The kurtosis measures how sharply peaked a distribution is, relative to its width. The kurtosis is normalized to zero for a Gaussian distribution. -- Function: double gsl_stats_kurtosis_m_sd (const double DATA[], size_t STRIDE, size_t N, double MEAN, double SD) This function computes the kurtosis of the dataset DATA using the given values of the mean MEAN and standard deviation SD, kurtosis = ((1/N) \sum ((x_i - mean)/sd)^4) - 3 This function is useful if you have already computed the mean and standard deviation of DATA and want to avoid recomputing them.  File: gsl-ref.info, Node: Autocorrelation, Next: Covariance, Prev: Higher moments (skewness and kurtosis), Up: Statistics 21.4 Autocorrelation ==================== -- Function: double gsl_stats_lag1_autocorrelation (const double DATA[], const size_t STRIDE, const size_t N) This function computes the lag-1 autocorrelation of the dataset DATA. a_1 = {\sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i-1} - \Hat\mu) \over \sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i} - \Hat\mu)} -- Function: double gsl_stats_lag1_autocorrelation_m (const double DATA[], const size_t STRIDE, const size_t N, const double MEAN) This function computes the lag-1 autocorrelation of the dataset DATA using the given value of the mean MEAN.  File: gsl-ref.info, Node: Covariance, Next: Correlation, Prev: Autocorrelation, Up: Statistics 21.5 Covariance =============== -- Function: double gsl_stats_covariance (const double DATA1[], const size_t STRIDE1, const double DATA2[], const size_t STRIDE2, const size_t N) This function computes the covariance of the datasets DATA1 and DATA2 which must both be of the same length N. covar = (1/(n - 1)) \sum_{i = 1}^{n} (x_i - \Hat x) (y_i - \Hat y) -- Function: double gsl_stats_covariance_m (const double DATA1[], const size_t STRIDE1, const double DATA2[], const size_t STRIDE2, const size_t N, const double MEAN1, const double MEAN2) This function computes the covariance of the datasets DATA1 and DATA2 using the given values of the means, MEAN1 and MEAN2. This is useful if you have already computed the means of DATA1 and DATA2 and want to avoid recomputing them.  File: gsl-ref.info, Node: Correlation, Next: Weighted Samples, Prev: Covariance, Up: Statistics 21.6 Correlation ================ -- Function: double gsl_stats_correlation (const double DATA1[], const size_t STRIDE1, const double DATA2[], const size_t STRIDE2, const size_t N) This function efficiently computes the Pearson correlation coefficient between the datasets DATA1 and DATA2 which must both be of the same length N. r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y) = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y) \over \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1) \sum (y_i - \Hat y)^2} } -- Function: double gsl_stats_spearman (const double DATA1[], const size_t STRIDE1, const double DATA2[], const size_t STRIDE2, const size_t N, double WORK[]) This function computes the Spearman rank correlation coefficient between the datasets DATA1 and DATA2 which must both be of the same length N. Additional workspace of size 2*N is required in WORK. The Spearman rank correlation between vectors x and y is equivalent to the Pearson correlation between the ranked vectors x_R and y_R, where ranks are defined to be the average of the positions of an element in the ascending order of the values.  File: gsl-ref.info, Node: Weighted Samples, Next: Maximum and Minimum values, Prev: Correlation, Up: Statistics 21.7 Weighted Samples ===================== The functions described in this section allow the computation of statistics for weighted samples. The functions accept an array of samples, x_i, with associated weights, w_i. Each sample x_i is considered as having been drawn from a Gaussian distribution with variance \sigma_i^2. The sample weight w_i is defined as the reciprocal of this variance, w_i = 1/\sigma_i^2. Setting a weight to zero corresponds to removing a sample from a dataset. -- Function: double gsl_stats_wmean (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) This function returns the weighted mean of the dataset DATA with stride STRIDE and length N, using the set of weights W with stride WSTRIDE and length N. The weighted mean is defined as, \Hat\mu = (\sum w_i x_i) / (\sum w_i) -- Function: double gsl_stats_wvariance (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) This function returns the estimated variance of the dataset DATA with stride STRIDE and length N, using the set of weights W with stride WSTRIDE and length N. The estimated variance of a weighted dataset is calculated as, \Hat\sigma^2 = ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2 Note that this expression reduces to an unweighted variance with the familiar 1/(N-1) factor when there are N equal non-zero weights. -- Function: double gsl_stats_wvariance_m (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN) This function returns the estimated variance of the weighted dataset DATA using the given weighted mean WMEAN. -- Function: double gsl_stats_wsd (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) The standard deviation is defined as the square root of the variance. This function returns the square root of the corresponding variance function 'gsl_stats_wvariance' above. -- Function: double gsl_stats_wsd_m (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN) This function returns the square root of the corresponding variance function 'gsl_stats_wvariance_m' above. -- Function: double gsl_stats_wvariance_with_fixed_mean (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, const double MEAN) This function computes an unbiased estimate of the variance of the weighted dataset DATA when the population mean MEAN of the underlying distribution is known _a priori_. In this case the estimator for the variance replaces the sample mean \Hat\mu by the known population mean \mu, \Hat\sigma^2 = (\sum w_i (x_i - \mu)^2) / (\sum w_i) -- Function: double gsl_stats_wsd_with_fixed_mean (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, const double MEAN) The standard deviation is defined as the square root of the variance. This function returns the square root of the corresponding variance function above. -- Function: double gsl_stats_wtss (const double W[], const size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) -- Function: double gsl_stats_wtss_m (const double W[], const size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN) These functions return the weighted total sum of squares (TSS) of DATA about the weighted mean. For 'gsl_stats_wtss_m' the user-supplied value of WMEAN is used, and for 'gsl_stats_wtss' it is computed using 'gsl_stats_wmean'. TSS = \sum w_i (x_i - wmean)^2 -- Function: double gsl_stats_wabsdev (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) This function computes the weighted absolute deviation from the weighted mean of DATA. The absolute deviation from the mean is defined as, absdev = (\sum w_i |x_i - \Hat\mu|) / (\sum w_i) -- Function: double gsl_stats_wabsdev_m (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN) This function computes the absolute deviation of the weighted dataset DATA about the given weighted mean WMEAN. -- Function: double gsl_stats_wskew (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) This function computes the weighted skewness of the dataset DATA. skew = (\sum w_i ((x_i - \Hat x)/\Hat \sigma)^3) / (\sum w_i) -- Function: double gsl_stats_wskew_m_sd (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN, double WSD) This function computes the weighted skewness of the dataset DATA using the given values of the weighted mean and weighted standard deviation, WMEAN and WSD. -- Function: double gsl_stats_wkurtosis (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) This function computes the weighted kurtosis of the dataset DATA. kurtosis = ((\sum w_i ((x_i - \Hat x)/\Hat \sigma)^4) / (\sum w_i)) - 3 -- Function: double gsl_stats_wkurtosis_m_sd (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN, double WSD) This function computes the weighted kurtosis of the dataset DATA using the given values of the weighted mean and weighted standard deviation, WMEAN and WSD.  File: gsl-ref.info, Node: Maximum and Minimum values, Next: Median and Percentiles, Prev: Weighted Samples, Up: Statistics 21.8 Maximum and Minimum values =============================== The following functions find the maximum and minimum values of a dataset (or their indices). If the data contains 'NaN's then a 'NaN' will be returned, since the maximum or minimum value is undefined. For functions which return an index, the location of the first 'NaN' in the array is returned. -- Function: double gsl_stats_max (const double DATA[], size_t STRIDE, size_t N) This function returns the maximum value in DATA, a dataset of length N with stride STRIDE. The maximum value is defined as the value of the element x_i which satisfies x_i >= x_j for all j. If you want instead to find the element with the largest absolute magnitude you will need to apply 'fabs' or 'abs' to your data before calling this function. -- Function: double gsl_stats_min (const double DATA[], size_t STRIDE, size_t N) This function returns the minimum value in DATA, a dataset of length N with stride STRIDE. The minimum value is defined as the value of the element x_i which satisfies x_i <= x_j for all j. If you want instead to find the element with the smallest absolute magnitude you will need to apply 'fabs' or 'abs' to your data before calling this function. -- Function: void gsl_stats_minmax (double * MIN, double * MAX, const double DATA[], size_t STRIDE, size_t N) This function finds both the minimum and maximum values MIN, MAX in DATA in a single pass. -- Function: size_t gsl_stats_max_index (const double DATA[], size_t STRIDE, size_t N) This function returns the index of the maximum value in DATA, a dataset of length N with stride STRIDE. The maximum value is defined as the value of the element x_i which satisfies x_i >= x_j for all j. When there are several equal maximum elements then the first one is chosen. -- Function: size_t gsl_stats_min_index (const double DATA[], size_t STRIDE, size_t N) This function returns the index of the minimum value in DATA, a dataset of length N with stride STRIDE. The minimum value is defined as the value of the element x_i which satisfies x_i >= x_j for all j. When there are several equal minimum elements then the first one is chosen. -- Function: void gsl_stats_minmax_index (size_t * MIN_INDEX, size_t * MAX_INDEX, const double DATA[], size_t STRIDE, size_t N) This function returns the indexes MIN_INDEX, MAX_INDEX of the minimum and maximum values in DATA in a single pass.  File: gsl-ref.info, Node: Median and Percentiles, Next: Example statistical programs, Prev: Maximum and Minimum values, Up: Statistics 21.9 Median and Percentiles =========================== The median and percentile functions described in this section operate on sorted data. For convenience we use "quantiles", measured on a scale of 0 to 1, instead of percentiles (which use a scale of 0 to 100). -- Function: double gsl_stats_median_from_sorted_data (const double SORTED_DATA[], size_t STRIDE, size_t N) This function returns the median value of SORTED_DATA, a dataset of length N with stride STRIDE. The elements of the array must be in ascending numerical order. There are no checks to see whether the data are sorted, so the function 'gsl_sort' should always be used first. When the dataset has an odd number of elements the median is the value of element (n-1)/2. When the dataset has an even number of elements the median is the mean of the two nearest middle values, elements (n-1)/2 and n/2. Since the algorithm for computing the median involves interpolation this function always returns a floating-point number, even for integer data types. -- Function: double gsl_stats_quantile_from_sorted_data (const double SORTED_DATA[], size_t STRIDE, size_t N, double F) This function returns a quantile value of SORTED_DATA, a double-precision array of length N with stride STRIDE. The elements of the array must be in ascending numerical order. The quantile is determined by the F, a fraction between 0 and 1. For example, to compute the value of the 75th percentile F should have the value 0.75. There are no checks to see whether the data are sorted, so the function 'gsl_sort' should always be used first. The quantile is found by interpolation, using the formula quantile = (1 - \delta) x_i + \delta x_{i+1} where i is 'floor'((n - 1)f) and \delta is (n-1)f - i. Thus the minimum value of the array ('data[0*stride]') is given by F equal to zero, the maximum value ('data[(n-1)*stride]') is given by F equal to one and the median value is given by F equal to 0.5. Since the algorithm for computing quantiles involves interpolation this function always returns a floating-point number, even for integer data types.  File: gsl-ref.info, Node: Example statistical programs, Next: Statistics References and Further Reading, Prev: Median and Percentiles, Up: Statistics 21.10 Examples ============== Here is a basic example of how to use the statistical functions: #include #include int main(void) { double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6}; double mean, variance, largest, smallest; mean = gsl_stats_mean(data, 1, 5); variance = gsl_stats_variance(data, 1, 5); largest = gsl_stats_max(data, 1, 5); smallest = gsl_stats_min(data, 1, 5); printf ("The dataset is %g, %g, %g, %g, %g\n", data[0], data[1], data[2], data[3], data[4]); printf ("The sample mean is %g\n", mean); printf ("The estimated variance is %g\n", variance); printf ("The largest value is %g\n", largest); printf ("The smallest value is %g\n", smallest); return 0; } The program should produce the following output, The dataset is 17.2, 18.1, 16.5, 18.3, 12.6 The sample mean is 16.54 The estimated variance is 5.373 The largest value is 18.3 The smallest value is 12.6 Here is an example using sorted data, #include #include #include int main(void) { double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6}; double median, upperq, lowerq; printf ("Original dataset: %g, %g, %g, %g, %g\n", data[0], data[1], data[2], data[3], data[4]); gsl_sort (data, 1, 5); printf ("Sorted dataset: %g, %g, %g, %g, %g\n", data[0], data[1], data[2], data[3], data[4]); median = gsl_stats_median_from_sorted_data (data, 1, 5); upperq = gsl_stats_quantile_from_sorted_data (data, 1, 5, 0.75); lowerq = gsl_stats_quantile_from_sorted_data (data, 1, 5, 0.25); printf ("The median is %g\n", median); printf ("The upper quartile is %g\n", upperq); printf ("The lower quartile is %g\n", lowerq); return 0; } This program should produce the following output, Original dataset: 17.2, 18.1, 16.5, 18.3, 12.6 Sorted dataset: 12.6, 16.5, 17.2, 18.1, 18.3 The median is 17.2 The upper quartile is 18.1 The lower quartile is 16.5  File: gsl-ref.info, Node: Statistics References and Further Reading, Prev: Example statistical programs, Up: Statistics 21.11 References and Further Reading ==================================== The standard reference for almost any topic in statistics is the multi-volume 'Advanced Theory of Statistics' by Kendall and Stuart. Maurice Kendall, Alan Stuart, and J. Keith Ord. 'The Advanced Theory of Statistics' (multiple volumes) reprinted as 'Kendall's Advanced Theory of Statistics'. Wiley, ISBN 047023380X. Many statistical concepts can be more easily understood by a Bayesian approach. The following book by Gelman, Carlin, Stern and Rubin gives a comprehensive coverage of the subject. Andrew Gelman, John B. Carlin, Hal S. Stern, Donald B. Rubin. 'Bayesian Data Analysis'. Chapman & Hall, ISBN 0412039915. For physicists the Particle Data Group provides useful reviews of Probability and Statistics in the "Mathematical Tools" section of its Annual Review of Particle Physics. 'Review of Particle Properties' R.M. Barnett et al., Physical Review D54, 1 (1996) The Review of Particle Physics is available online at the website .  File: gsl-ref.info, Node: Histograms, Next: N-tuples, Prev: Statistics, Up: Top 22 Histograms ************* This chapter describes functions for creating histograms. Histograms provide a convenient way of summarizing the distribution of a set of data. A histogram consists of a set of "bins" which count the number of events falling into a given range of a continuous variable x. In GSL the bins of a histogram contain floating-point numbers, so they can be used to record both integer and non-integer distributions. The bins can use arbitrary sets of ranges (uniformly spaced bins are the default). Both one and two-dimensional histograms are supported. Once a histogram has been created it can also be converted into a probability distribution function. The library provides efficient routines for selecting random samples from probability distributions. This can be useful for generating simulations based on real data. The functions are declared in the header files 'gsl_histogram.h' and 'gsl_histogram2d.h'. * Menu: * The histogram struct:: * Histogram allocation:: * Copying Histograms:: * Updating and accessing histogram elements:: * Searching histogram ranges:: * Histogram Statistics:: * Histogram Operations:: * Reading and writing histograms:: * Resampling from histograms:: * The histogram probability distribution struct:: * Example programs for histograms:: * Two dimensional histograms:: * The 2D histogram struct:: * 2D Histogram allocation:: * Copying 2D Histograms:: * Updating and accessing 2D histogram elements:: * Searching 2D histogram ranges:: * 2D Histogram Statistics:: * 2D Histogram Operations:: * Reading and writing 2D histograms:: * Resampling from 2D histograms:: * Example programs for 2D histograms::  File: gsl-ref.info, Node: The histogram struct, Next: Histogram allocation, Up: Histograms 22.1 The histogram struct ========================= A histogram is defined by the following struct, -- Data Type: gsl_histogram 'size_t n' This is the number of histogram bins 'double * range' The ranges of the bins are stored in an array of N+1 elements pointed to by RANGE. 'double * bin' The counts for each bin are stored in an array of N elements pointed to by BIN. The bins are floating-point numbers, so you can increment them by non-integer values if necessary. The range for BIN[i] is given by RANGE[i] to RANGE[i+1]. For n bins there are n+1 entries in the array RANGE. Each bin is inclusive at the lower end and exclusive at the upper end. Mathematically this means that the bins are defined by the following inequality, bin[i] corresponds to range[i] <= x < range[i+1] Here is a diagram of the correspondence between ranges and bins on the number-line for x, [ bin[0] )[ bin[1] )[ bin[2] )[ bin[3] )[ bin[4] ) ---|---------|---------|---------|---------|---------|--- x r[0] r[1] r[2] r[3] r[4] r[5] In this picture the values of the RANGE array are denoted by r. On the left-hand side of each bin the square bracket '[' denotes an inclusive lower bound (r <= x), and the round parentheses ')' on the right-hand side denote an exclusive upper bound (x < r). Thus any samples which fall on the upper end of the histogram are excluded. If you want to include this value for the last bin you will need to add an extra bin to your histogram. The 'gsl_histogram' struct and its associated functions are defined in the header file 'gsl_histogram.h'.  File: gsl-ref.info, Node: Histogram allocation, Next: Copying Histograms, Prev: The histogram struct, Up: Histograms 22.2 Histogram allocation ========================= The functions for allocating memory to a histogram follow the style of 'malloc' and 'free'. In addition they also perform their own error checking. If there is insufficient memory available to allocate a histogram then the functions call the error handler (with an error number of 'GSL_ENOMEM') in addition to returning a null pointer. Thus if you use the library error handler to abort your program then it isn't necessary to check every histogram 'alloc'. -- Function: gsl_histogram * gsl_histogram_alloc (size_t N) This function allocates memory for a histogram with N bins, and returns a pointer to a newly created 'gsl_histogram' struct. If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of 'GSL_ENOMEM'. The bins and ranges are not initialized, and should be prepared using one of the range-setting functions below in order to make the histogram ready for use. -- Function: int gsl_histogram_set_ranges (gsl_histogram * H, const double RANGE[], size_t SIZE) This function sets the ranges of the existing histogram H using the array RANGE of size SIZE. The values of the histogram bins are reset to zero. The 'range' array should contain the desired bin limits. The ranges can be arbitrary, subject to the restriction that they are monotonically increasing. The following example shows how to create a histogram with logarithmic bins with ranges [1,10), [10,100) and [100,1000). gsl_histogram * h = gsl_histogram_alloc (3); /* bin[0] covers the range 1 <= x < 10 */ /* bin[1] covers the range 10 <= x < 100 */ /* bin[2] covers the range 100 <= x < 1000 */ double range[4] = { 1.0, 10.0, 100.0, 1000.0 }; gsl_histogram_set_ranges (h, range, 4); Note that the size of the RANGE array should be defined to be one element bigger than the number of bins. The additional element is required for the upper value of the final bin. -- Function: int gsl_histogram_set_ranges_uniform (gsl_histogram * H, double XMIN, double XMAX) This function sets the ranges of the existing histogram H to cover the range XMIN to XMAX uniformly. The values of the histogram bins are reset to zero. The bin ranges are shown in the table below, bin[0] corresponds to xmin <= x < xmin + d bin[1] corresponds to xmin + d <= x < xmin + 2 d ...... bin[n-1] corresponds to xmin + (n-1)d <= x < xmax where d is the bin spacing, d = (xmax-xmin)/n. -- Function: void gsl_histogram_free (gsl_histogram * H) This function frees the histogram H and all of the memory associated with it.  File: gsl-ref.info, Node: Copying Histograms, Next: Updating and accessing histogram elements, Prev: Histogram allocation, Up: Histograms 22.3 Copying Histograms ======================= -- Function: int gsl_histogram_memcpy (gsl_histogram * DEST, const gsl_histogram * SRC) This function copies the histogram SRC into the pre-existing histogram DEST, making DEST into an exact copy of SRC. The two histograms must be of the same size. -- Function: gsl_histogram * gsl_histogram_clone (const gsl_histogram * SRC) This function returns a pointer to a newly created histogram which is an exact copy of the histogram SRC.  File: gsl-ref.info, Node: Updating and accessing histogram elements, Next: Searching histogram ranges, Prev: Copying Histograms, Up: Histograms 22.4 Updating and accessing histogram elements ============================================== There are two ways to access histogram bins, either by specifying an x coordinate or by using the bin-index directly. The functions for accessing the histogram through x coordinates use a binary search to identify the bin which covers the appropriate range. -- Function: int gsl_histogram_increment (gsl_histogram * H, double X) This function updates the histogram H by adding one (1.0) to the bin whose range contains the coordinate X. If X lies in the valid range of the histogram then the function returns zero to indicate success. If X is less than the lower limit of the histogram then the function returns 'GSL_EDOM', and none of bins are modified. Similarly, if the value of X is greater than or equal to the upper limit of the histogram then the function returns 'GSL_EDOM', and none of the bins are modified. The error handler is not called, however, since it is often necessary to compute histograms for a small range of a larger dataset, ignoring the values outside the range of interest. -- Function: int gsl_histogram_accumulate (gsl_histogram * H, double X, double WEIGHT) This function is similar to 'gsl_histogram_increment' but increases the value of the appropriate bin in the histogram H by the floating-point number WEIGHT. -- Function: double gsl_histogram_get (const gsl_histogram * H, size_t I) This function returns the contents of the I-th bin of the histogram H. If I lies outside the valid range of indices for the histogram then the error handler is called with an error code of 'GSL_EDOM' and the function returns 0. -- Function: int gsl_histogram_get_range (const gsl_histogram * H, size_t I, double * LOWER, double * UPPER) This function finds the upper and lower range limits of the I-th bin of the histogram H. If the index I is valid then the corresponding range limits are stored in LOWER and UPPER. The lower limit is inclusive (i.e. events with this coordinate are included in the bin) and the upper limit is exclusive (i.e. events with the coordinate of the upper limit are excluded and fall in the neighboring higher bin, if it exists). The function returns 0 to indicate success. If I lies outside the valid range of indices for the histogram then the error handler is called and the function returns an error code of 'GSL_EDOM'. -- Function: double gsl_histogram_max (const gsl_histogram * H) -- Function: double gsl_histogram_min (const gsl_histogram * H) -- Function: size_t gsl_histogram_bins (const gsl_histogram * H) These functions return the maximum upper and minimum lower range limits and the number of bins of the histogram H. They provide a way of determining these values without accessing the 'gsl_histogram' struct directly. -- Function: void gsl_histogram_reset (gsl_histogram * H) This function resets all the bins in the histogram H to zero.  File: gsl-ref.info, Node: Searching histogram ranges, Next: Histogram Statistics, Prev: Updating and accessing histogram elements, Up: Histograms 22.5 Searching histogram ranges =============================== The following functions are used by the access and update routines to locate the bin which corresponds to a given x coordinate. -- Function: int gsl_histogram_find (const gsl_histogram * H, double X, size_t * I) This function finds and sets the index I to the bin number which covers the coordinate X in the histogram H. The bin is located using a binary search. The search includes an optimization for histograms with uniform range, and will return the correct bin immediately in this case. If X is found in the range of the histogram then the function sets the index I and returns 'GSL_SUCCESS'. If X lies outside the valid range of the histogram then the function returns 'GSL_EDOM' and the error handler is invoked.  File: gsl-ref.info, Node: Histogram Statistics, Next: Histogram Operations, Prev: Searching histogram ranges, Up: Histograms 22.6 Histogram Statistics ========================= -- Function: double gsl_histogram_max_val (const gsl_histogram * H) This function returns the maximum value contained in the histogram bins. -- Function: size_t gsl_histogram_max_bin (const gsl_histogram * H) This function returns the index of the bin containing the maximum value. In the case where several bins contain the same maximum value the smallest index is returned. -- Function: double gsl_histogram_min_val (const gsl_histogram * H) This function returns the minimum value contained in the histogram bins. -- Function: size_t gsl_histogram_min_bin (const gsl_histogram * H) This function returns the index of the bin containing the minimum value. In the case where several bins contain the same maximum value the smallest index is returned. -- Function: double gsl_histogram_mean (const gsl_histogram * H) This function returns the mean of the histogrammed variable, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. The accuracy of the result is limited by the bin width. -- Function: double gsl_histogram_sigma (const gsl_histogram * H) This function returns the standard deviation of the histogrammed variable, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. The accuracy of the result is limited by the bin width. -- Function: double gsl_histogram_sum (const gsl_histogram * H) This function returns the sum of all bin values. Negative bin values are included in the sum.  File: gsl-ref.info, Node: Histogram Operations, Next: Reading and writing histograms, Prev: Histogram Statistics, Up: Histograms 22.7 Histogram Operations ========================= -- Function: int gsl_histogram_equal_bins_p (const gsl_histogram * H1, const gsl_histogram * H2) This function returns 1 if the all of the individual bin ranges of the two histograms are identical, and 0 otherwise. -- Function: int gsl_histogram_add (gsl_histogram * H1, const gsl_histogram * H2) This function adds the contents of the bins in histogram H2 to the corresponding bins of histogram H1, i.e. h'_1(i) = h_1(i) + h_2(i). The two histograms must have identical bin ranges. -- Function: int gsl_histogram_sub (gsl_histogram * H1, const gsl_histogram * H2) This function subtracts the contents of the bins in histogram H2 from the corresponding bins of histogram H1, i.e. h'_1(i) = h_1(i) - h_2(i). The two histograms must have identical bin ranges. -- Function: int gsl_histogram_mul (gsl_histogram * H1, const gsl_histogram * H2) This function multiplies the contents of the bins of histogram H1 by the contents of the corresponding bins in histogram H2, i.e. h'_1(i) = h_1(i) * h_2(i). The two histograms must have identical bin ranges. -- Function: int gsl_histogram_div (gsl_histogram * H1, const gsl_histogram * H2) This function divides the contents of the bins of histogram H1 by the contents of the corresponding bins in histogram H2, i.e. h'_1(i) = h_1(i) / h_2(i). The two histograms must have identical bin ranges. -- Function: int gsl_histogram_scale (gsl_histogram * H, double SCALE) This function multiplies the contents of the bins of histogram H by the constant SCALE, i.e. h'_1(i) = h_1(i) * scale. -- Function: int gsl_histogram_shift (gsl_histogram * H, double OFFSET) This function shifts the contents of the bins of histogram H by the constant OFFSET, i.e. h'_1(i) = h_1(i) + offset.  File: gsl-ref.info, Node: Reading and writing histograms, Next: Resampling from histograms, Prev: Histogram Operations, Up: Histograms 22.8 Reading and writing histograms =================================== The library provides functions for reading and writing histograms to a file as binary data or formatted text. -- Function: int gsl_histogram_fwrite (FILE * STREAM, const gsl_histogram * H) This function writes the ranges and bins of the histogram H to the stream STREAM in binary format. The return value is 0 for success and 'GSL_EFAILED' if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures. -- Function: int gsl_histogram_fread (FILE * STREAM, gsl_histogram * H) This function reads into the histogram H from the open stream STREAM in binary format. The histogram H must be preallocated with the correct size since the function uses the number of bins in H to determine how many bytes to read. The return value is 0 for success and 'GSL_EFAILED' if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture. -- Function: int gsl_histogram_fprintf (FILE * STREAM, const gsl_histogram * H, const char * RANGE_FORMAT, const char * BIN_FORMAT) This function writes the ranges and bins of the histogram H line-by-line to the stream STREAM using the format specifiers RANGE_FORMAT and BIN_FORMAT. These should be one of the '%g', '%e' or '%f' formats for floating point numbers. The function returns 0 for success and 'GSL_EFAILED' if there was a problem writing to the file. The histogram output is formatted in three columns, and the columns are separated by spaces, like this, range[0] range[1] bin[0] range[1] range[2] bin[1] range[2] range[3] bin[2] .... range[n-1] range[n] bin[n-1] The values of the ranges are formatted using RANGE_FORMAT and the value of the bins are formatted using BIN_FORMAT. Each line contains the lower and upper limit of the range of the bins and the value of the bin itself. Since the upper limit of one bin is the lower limit of the next there is duplication of these values between lines but this allows the histogram to be manipulated with line-oriented tools. -- Function: int gsl_histogram_fscanf (FILE * STREAM, gsl_histogram * H) This function reads formatted data from the stream STREAM into the histogram H. The data is assumed to be in the three-column format used by 'gsl_histogram_fprintf'. The histogram H must be preallocated with the correct length since the function uses the size of H to determine how many numbers to read. The function returns 0 for success and 'GSL_EFAILED' if there was a problem reading from the file.  File: gsl-ref.info, Node: Resampling from histograms, Next: The histogram probability distribution struct, Prev: Reading and writing histograms, Up: Histograms 22.9 Resampling from histograms =============================== A histogram made by counting events can be regarded as a measurement of a probability distribution. Allowing for statistical error, the height of each bin represents the probability of an event where the value of x falls in the range of that bin. The probability distribution function has the one-dimensional form p(x)dx where, p(x) = n_i/ (N w_i) In this equation n_i is the number of events in the bin which contains x, w_i is the width of the bin and N is the total number of events. The distribution of events within each bin is assumed to be uniform.  File: gsl-ref.info, Node: The histogram probability distribution struct, Next: Example programs for histograms, Prev: Resampling from histograms, Up: Histograms 22.10 The histogram probability distribution struct =================================================== The probability distribution function for a histogram consists of a set of "bins" which measure the probability of an event falling into a given range of a continuous variable x. A probability distribution function is defined by the following struct, which actually stores the cumulative probability distribution function. This is the natural quantity for generating samples via the inverse transform method, because there is a one-to-one mapping between the cumulative probability distribution and the range [0,1]. It can be shown that by taking a uniform random number in this range and finding its corresponding coordinate in the cumulative probability distribution we obtain samples with the desired probability distribution. -- Data Type: gsl_histogram_pdf 'size_t n' This is the number of bins used to approximate the probability distribution function. 'double * range' The ranges of the bins are stored in an array of N+1 elements pointed to by RANGE. 'double * sum' The cumulative probability for the bins is stored in an array of N elements pointed to by SUM. The following functions allow you to create a 'gsl_histogram_pdf' struct which represents this probability distribution and generate random samples from it. -- Function: gsl_histogram_pdf * gsl_histogram_pdf_alloc (size_t N) This function allocates memory for a probability distribution with N bins and returns a pointer to a newly initialized 'gsl_histogram_pdf' struct. If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of 'GSL_ENOMEM'. -- Function: int gsl_histogram_pdf_init (gsl_histogram_pdf * P, const gsl_histogram * H) This function initializes the probability distribution P with the contents of the histogram H. If any of the bins of H are negative then the error handler is invoked with an error code of 'GSL_EDOM' because a probability distribution cannot contain negative values. -- Function: void gsl_histogram_pdf_free (gsl_histogram_pdf * P) This function frees the probability distribution function P and all of the memory associated with it. -- Function: double gsl_histogram_pdf_sample (const gsl_histogram_pdf * P, double R) This function uses R, a uniform random number between zero and one, to compute a single random sample from the probability distribution P. The algorithm used to compute the sample s is given by the following formula, s = range[i] + delta * (range[i+1] - range[i]) where i is the index which satisfies sum[i] <= r < sum[i+1] and delta is (r - sum[i])/(sum[i+1] - sum[i]).  File: gsl-ref.info, Node: Example programs for histograms, Next: Two dimensional histograms, Prev: The histogram probability distribution struct, Up: Histograms 22.11 Example programs for histograms ===================================== The following program shows how to make a simple histogram of a column of numerical data supplied on 'stdin'. The program takes three arguments, specifying the upper and lower bounds of the histogram and the number of bins. It then reads numbers from 'stdin', one line at a time, and adds them to the histogram. When there is no more data to read it prints out the accumulated histogram using 'gsl_histogram_fprintf'. #include #include #include int main (int argc, char **argv) { double a, b; size_t n; if (argc != 4) { printf ("Usage: gsl-histogram xmin xmax n\n" "Computes a histogram of the data " "on stdin using n bins from xmin " "to xmax\n"); exit (0); } a = atof (argv[1]); b = atof (argv[2]); n = atoi (argv[3]); { double x; gsl_histogram * h = gsl_histogram_alloc (n); gsl_histogram_set_ranges_uniform (h, a, b); while (fscanf (stdin, "%lg", &x) == 1) { gsl_histogram_increment (h, x); } gsl_histogram_fprintf (stdout, h, "%g", "%g"); gsl_histogram_free (h); } exit (0); } Here is an example of the program in use. We generate 10000 random samples from a Cauchy distribution with a width of 30 and histogram them over the range -100 to 100, using 200 bins. $ gsl-randist 0 10000 cauchy 30 | gsl-histogram -100 100 200 > histogram.dat A plot of the resulting histogram shows the familiar shape of the Cauchy distribution and the fluctuations caused by the finite sample size. $ awk '{print $1, $3 ; print $2, $3}' histogram.dat | graph -T X  File: gsl-ref.info, Node: Two dimensional histograms, Next: The 2D histogram struct, Prev: Example programs for histograms, Up: Histograms 22.12 Two dimensional histograms ================================ A two dimensional histogram consists of a set of "bins" which count the number of events falling in a given area of the (x,y) plane. The simplest way to use a two dimensional histogram is to record two-dimensional position information, n(x,y). Another possibility is to form a "joint distribution" by recording related variables. For example a detector might record both the position of an event (x) and the amount of energy it deposited E. These could be histogrammed as the joint distribution n(x,E).  File: gsl-ref.info, Node: The 2D histogram struct, Next: 2D Histogram allocation, Prev: Two dimensional histograms, Up: Histograms 22.13 The 2D histogram struct ============================= Two dimensional histograms are defined by the following struct, -- Data Type: gsl_histogram2d 'size_t nx, ny' This is the number of histogram bins in the x and y directions. 'double * xrange' The ranges of the bins in the x-direction are stored in an array of NX + 1 elements pointed to by XRANGE. 'double * yrange' The ranges of the bins in the y-direction are stored in an array of NY + 1 elements pointed to by YRANGE. 'double * bin' The counts for each bin are stored in an array pointed to by BIN. The bins are floating-point numbers, so you can increment them by non-integer values if necessary. The array BIN stores the two dimensional array of bins in a single block of memory according to the mapping 'bin(i,j)' = 'bin[i * ny + j]'. The range for 'bin(i,j)' is given by 'xrange[i]' to 'xrange[i+1]' in the x-direction and 'yrange[j]' to 'yrange[j+1]' in the y-direction. Each bin is inclusive at the lower end and exclusive at the upper end. Mathematically this means that the bins are defined by the following inequality, bin(i,j) corresponds to xrange[i] <= x < xrange[i+1] and yrange[j] <= y < yrange[j+1] Note that any samples which fall on the upper sides of the histogram are excluded. If you want to include these values for the side bins you will need to add an extra row or column to your histogram. The 'gsl_histogram2d' struct and its associated functions are defined in the header file 'gsl_histogram2d.h'.  File: gsl-ref.info, Node: 2D Histogram allocation, Next: Copying 2D Histograms, Prev: The 2D histogram struct, Up: Histograms 22.14 2D Histogram allocation ============================= The functions for allocating memory to a 2D histogram follow the style of 'malloc' and 'free'. In addition they also perform their own error checking. If there is insufficient memory available to allocate a histogram then the functions call the error handler (with an error number of 'GSL_ENOMEM') in addition to returning a null pointer. Thus if you use the library error handler to abort your program then it isn't necessary to check every 2D histogram 'alloc'. -- Function: gsl_histogram2d * gsl_histogram2d_alloc (size_t NX, size_t NY) This function allocates memory for a two-dimensional histogram with NX bins in the x direction and NY bins in the y direction. The function returns a pointer to a newly created 'gsl_histogram2d' struct. If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of 'GSL_ENOMEM'. The bins and ranges must be initialized with one of the functions below before the histogram is ready for use. -- Function: int gsl_histogram2d_set_ranges (gsl_histogram2d * H, const double XRANGE[], size_t XSIZE, const double YRANGE[], size_t YSIZE) This function sets the ranges of the existing histogram H using the arrays XRANGE and YRANGE of size XSIZE and YSIZE respectively. The values of the histogram bins are reset to zero. -- Function: int gsl_histogram2d_set_ranges_uniform (gsl_histogram2d * H, double XMIN, double XMAX, double YMIN, double YMAX) This function sets the ranges of the existing histogram H to cover the ranges XMIN to XMAX and YMIN to YMAX uniformly. The values of the histogram bins are reset to zero. -- Function: void gsl_histogram2d_free (gsl_histogram2d * H) This function frees the 2D histogram H and all of the memory associated with it.  File: gsl-ref.info, Node: Copying 2D Histograms, Next: Updating and accessing 2D histogram elements, Prev: 2D Histogram allocation, Up: Histograms 22.15 Copying 2D Histograms =========================== -- Function: int gsl_histogram2d_memcpy (gsl_histogram2d * DEST, const gsl_histogram2d * SRC) This function copies the histogram SRC into the pre-existing histogram DEST, making DEST into an exact copy of SRC. The two histograms must be of the same size. -- Function: gsl_histogram2d * gsl_histogram2d_clone (const gsl_histogram2d * SRC) This function returns a pointer to a newly created histogram which is an exact copy of the histogram SRC.  File: gsl-ref.info, Node: Updating and accessing 2D histogram elements, Next: Searching 2D histogram ranges, Prev: Copying 2D Histograms, Up: Histograms 22.16 Updating and accessing 2D histogram elements ================================================== You can access the bins of a two-dimensional histogram either by specifying a pair of (x,y) coordinates or by using the bin indices (i,j) directly. The functions for accessing the histogram through (x,y) coordinates use binary searches in the x and y directions to identify the bin which covers the appropriate range. -- Function: int gsl_histogram2d_increment (gsl_histogram2d * H, double X, double Y) This function updates the histogram H by adding one (1.0) to the bin whose x and y ranges contain the coordinates (X,Y). If the point (x,y) lies inside the valid ranges of the histogram then the function returns zero to indicate success. If (x,y) lies outside the limits of the histogram then the function returns 'GSL_EDOM', and none of the bins are modified. The error handler is not called, since it is often necessary to compute histograms for a small range of a larger dataset, ignoring any coordinates outside the range of interest. -- Function: int gsl_histogram2d_accumulate (gsl_histogram2d * H, double X, double Y, double WEIGHT) This function is similar to 'gsl_histogram2d_increment' but increases the value of the appropriate bin in the histogram H by the floating-point number WEIGHT. -- Function: double gsl_histogram2d_get (const gsl_histogram2d * H, size_t I, size_t J) This function returns the contents of the (I,J)-th bin of the histogram H. If (I,J) lies outside the valid range of indices for the histogram then the error handler is called with an error code of 'GSL_EDOM' and the function returns 0. -- Function: int gsl_histogram2d_get_xrange (const gsl_histogram2d * H, size_t I, double * XLOWER, double * XUPPER) -- Function: int gsl_histogram2d_get_yrange (const gsl_histogram2d * H, size_t J, double * YLOWER, double * YUPPER) These functions find the upper and lower range limits of the I-th and J-th bins in the x and y directions of the histogram H. The range limits are stored in XLOWER and XUPPER or YLOWER and YUPPER. The lower limits are inclusive (i.e. events with these coordinates are included in the bin) and the upper limits are exclusive (i.e. events with the value of the upper limit are not included and fall in the neighboring higher bin, if it exists). The functions return 0 to indicate success. If I or J lies outside the valid range of indices for the histogram then the error handler is called with an error code of 'GSL_EDOM'. -- Function: double gsl_histogram2d_xmax (const gsl_histogram2d * H) -- Function: double gsl_histogram2d_xmin (const gsl_histogram2d * H) -- Function: size_t gsl_histogram2d_nx (const gsl_histogram2d * H) -- Function: double gsl_histogram2d_ymax (const gsl_histogram2d * H) -- Function: double gsl_histogram2d_ymin (const gsl_histogram2d * H) -- Function: size_t gsl_histogram2d_ny (const gsl_histogram2d * H) These functions return the maximum upper and minimum lower range limits and the number of bins for the x and y directions of the histogram H. They provide a way of determining these values without accessing the 'gsl_histogram2d' struct directly. -- Function: void gsl_histogram2d_reset (gsl_histogram2d * H) This function resets all the bins of the histogram H to zero.  File: gsl-ref.info, Node: Searching 2D histogram ranges, Next: 2D Histogram Statistics, Prev: Updating and accessing 2D histogram elements, Up: Histograms 22.17 Searching 2D histogram ranges =================================== The following functions are used by the access and update routines to locate the bin which corresponds to a given (x,y) coordinate. -- Function: int gsl_histogram2d_find (const gsl_histogram2d * H, double X, double Y, size_t * I, size_t * J) This function finds and sets the indices I and J to the bin which covers the coordinates (X,Y). The bin is located using a binary search. The search includes an optimization for histograms with uniform ranges, and will return the correct bin immediately in this case. If (x,y) is found then the function sets the indices (I,J) and returns 'GSL_SUCCESS'. If (x,y) lies outside the valid range of the histogram then the function returns 'GSL_EDOM' and the error handler is invoked.  File: gsl-ref.info, Node: 2D Histogram Statistics, Next: 2D Histogram Operations, Prev: Searching 2D histogram ranges, Up: Histograms 22.18 2D Histogram Statistics ============================= -- Function: double gsl_histogram2d_max_val (const gsl_histogram2d * H) This function returns the maximum value contained in the histogram bins. -- Function: void gsl_histogram2d_max_bin (const gsl_histogram2d * H, size_t * I, size_t * J) This function finds the indices of the bin containing the maximum value in the histogram H and stores the result in (I,J). In the case where several bins contain the same maximum value the first bin found is returned. -- Function: double gsl_histogram2d_min_val (const gsl_histogram2d * H) This function returns the minimum value contained in the histogram bins. -- Function: void gsl_histogram2d_min_bin (const gsl_histogram2d * H, size_t * I, size_t * J) This function finds the indices of the bin containing the minimum value in the histogram H and stores the result in (I,J). In the case where several bins contain the same maximum value the first bin found is returned. -- Function: double gsl_histogram2d_xmean (const gsl_histogram2d * H) This function returns the mean of the histogrammed x variable, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. -- Function: double gsl_histogram2d_ymean (const gsl_histogram2d * H) This function returns the mean of the histogrammed y variable, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. -- Function: double gsl_histogram2d_xsigma (const gsl_histogram2d * H) This function returns the standard deviation of the histogrammed x variable, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. -- Function: double gsl_histogram2d_ysigma (const gsl_histogram2d * H) This function returns the standard deviation of the histogrammed y variable, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. -- Function: double gsl_histogram2d_cov (const gsl_histogram2d * H) This function returns the covariance of the histogrammed x and y variables, where the histogram is regarded as a probability distribution. Negative bin values are ignored for the purposes of this calculation. -- Function: double gsl_histogram2d_sum (const gsl_histogram2d * H) This function returns the sum of all bin values. Negative bin values are included in the sum.  File: gsl-ref.info, Node: 2D Histogram Operations, Next: Reading and writing 2D histograms, Prev: 2D Histogram Statistics, Up: Histograms 22.19 2D Histogram Operations ============================= -- Function: int gsl_histogram2d_equal_bins_p (const gsl_histogram2d * H1, const gsl_histogram2d * H2) This function returns 1 if all the individual bin ranges of the two histograms are identical, and 0 otherwise. -- Function: int gsl_histogram2d_add (gsl_histogram2d * H1, const gsl_histogram2d * H2) This function adds the contents of the bins in histogram H2 to the corresponding bins of histogram H1, i.e. h'_1(i,j) = h_1(i,j) + h_2(i,j). The two histograms must have identical bin ranges. -- Function: int gsl_histogram2d_sub (gsl_histogram2d * H1, const gsl_histogram2d * H2) This function subtracts the contents of the bins in histogram H2 from the corresponding bins of histogram H1, i.e. h'_1(i,j) = h_1(i,j) - h_2(i,j). The two histograms must have identical bin ranges. -- Function: int gsl_histogram2d_mul (gsl_histogram2d * H1, const gsl_histogram2d * H2) This function multiplies the contents of the bins of histogram H1 by the contents of the corresponding bins in histogram H2, i.e. h'_1(i,j) = h_1(i,j) * h_2(i,j). The two histograms must have identical bin ranges. -- Function: int gsl_histogram2d_div (gsl_histogram2d * H1, const gsl_histogram2d * H2) This function divides the contents of the bins of histogram H1 by the contents of the corresponding bins in histogram H2, i.e. h'_1(i,j) = h_1(i,j) / h_2(i,j). The two histograms must have identical bin ranges. -- Function: int gsl_histogram2d_scale (gsl_histogram2d * H, double SCALE) This function multiplies the contents of the bins of histogram H by the constant SCALE, i.e. h'_1(i,j) = h_1(i,j) scale. -- Function: int gsl_histogram2d_shift (gsl_histogram2d * H, double OFFSET) This function shifts the contents of the bins of histogram H by the constant OFFSET, i.e. h'_1(i,j) = h_1(i,j) + offset.  File: gsl-ref.info, Node: Reading and writing 2D histograms, Next: Resampling from 2D histograms, Prev: 2D Histogram Operations, Up: Histograms 22.20 Reading and writing 2D histograms ======================================= The library provides functions for reading and writing two dimensional histograms to a file as binary data or formatted text. -- Function: int gsl_histogram2d_fwrite (FILE * STREAM, const gsl_histogram2d * H) This function writes the ranges and bins of the histogram H to the stream STREAM in binary format. The return value is 0 for success and 'GSL_EFAILED' if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures. -- Function: int gsl_histogram2d_fread (FILE * STREAM, gsl_histogram2d * H) This function reads into the histogram H from the stream STREAM in binary format. The histogram H must be preallocated with the correct size since the function uses the number of x and y bins in H to determine how many bytes to read. The return value is 0 for success and 'GSL_EFAILED' if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture. -- Function: int gsl_histogram2d_fprintf (FILE * STREAM, const gsl_histogram2d * H, const char * RANGE_FORMAT, const char * BIN_FORMAT) This function writes the ranges and bins of the histogram H line-by-line to the stream STREAM using the format specifiers RANGE_FORMAT and BIN_FORMAT. These should be one of the '%g', '%e' or '%f' formats for floating point numbers. The function returns 0 for success and 'GSL_EFAILED' if there was a problem writing to the file. The histogram output is formatted in five columns, and the columns are separated by spaces, like this, xrange[0] xrange[1] yrange[0] yrange[1] bin(0,0) xrange[0] xrange[1] yrange[1] yrange[2] bin(0,1) xrange[0] xrange[1] yrange[2] yrange[3] bin(0,2) .... xrange[0] xrange[1] yrange[ny-1] yrange[ny] bin(0,ny-1) xrange[1] xrange[2] yrange[0] yrange[1] bin(1,0) xrange[1] xrange[2] yrange[1] yrange[2] bin(1,1) xrange[1] xrange[2] yrange[1] yrange[2] bin(1,2) .... xrange[1] xrange[2] yrange[ny-1] yrange[ny] bin(1,ny-1) .... xrange[nx-1] xrange[nx] yrange[0] yrange[1] bin(nx-1,0) xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,1) xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,2) .... xrange[nx-1] xrange[nx] yrange[ny-1] yrange[ny] bin(nx-1,ny-1) Each line contains the lower and upper limits of the bin and the contents of the bin. Since the upper limits of the each bin are the lower limits of the neighboring bins there is duplication of these values but this allows the histogram to be manipulated with line-oriented tools. -- Function: int gsl_histogram2d_fscanf (FILE * STREAM, gsl_histogram2d * H) This function reads formatted data from the stream STREAM into the histogram H. The data is assumed to be in the five-column format used by 'gsl_histogram2d_fprintf'. The histogram H must be preallocated with the correct lengths since the function uses the sizes of H to determine how many numbers to read. The function returns 0 for success and 'GSL_EFAILED' if there was a problem reading from the file.  File: gsl-ref.info, Node: Resampling from 2D histograms, Next: Example programs for 2D histograms, Prev: Reading and writing 2D histograms, Up: Histograms 22.21 Resampling from 2D histograms =================================== As in the one-dimensional case, a two-dimensional histogram made by counting events can be regarded as a measurement of a probability distribution. Allowing for statistical error, the height of each bin represents the probability of an event where (x,y) falls in the range of that bin. For a two-dimensional histogram the probability distribution takes the form p(x,y) dx dy where, p(x,y) = n_{ij}/ (N A_{ij}) In this equation n_{ij} is the number of events in the bin which contains (x,y), A_{ij} is the area of the bin and N is the total number of events. The distribution of events within each bin is assumed to be uniform. -- Data Type: gsl_histogram2d_pdf 'size_t nx, ny' This is the number of histogram bins used to approximate the probability distribution function in the x and y directions. 'double * xrange' The ranges of the bins in the x-direction are stored in an array of NX + 1 elements pointed to by XRANGE. 'double * yrange' The ranges of the bins in the y-direction are stored in an array of NY + 1 pointed to by YRANGE. 'double * sum' The cumulative probability for the bins is stored in an array of NX*NY elements pointed to by SUM. The following functions allow you to create a 'gsl_histogram2d_pdf' struct which represents a two dimensional probability distribution and generate random samples from it. -- Function: gsl_histogram2d_pdf * gsl_histogram2d_pdf_alloc (size_t NX, size_t NY) This function allocates memory for a two-dimensional probability distribution of size NX-by-NY and returns a pointer to a newly initialized 'gsl_histogram2d_pdf' struct. If insufficient memory is available a null pointer is returned and the error handler is invoked with an error code of 'GSL_ENOMEM'. -- Function: int gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * P, const gsl_histogram2d * H) This function initializes the two-dimensional probability distribution calculated P from the histogram H. If any of the bins of H are negative then the error handler is invoked with an error code of 'GSL_EDOM' because a probability distribution cannot contain negative values. -- Function: void gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * P) This function frees the two-dimensional probability distribution function P and all of the memory associated with it. -- Function: int gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * P, double R1, double R2, double * X, double * Y) This function uses two uniform random numbers between zero and one, R1 and R2, to compute a single random sample from the two-dimensional probability distribution P.  File: gsl-ref.info, Node: Example programs for 2D histograms, Prev: Resampling from 2D histograms, Up: Histograms 22.22 Example programs for 2D histograms ======================================== This program demonstrates two features of two-dimensional histograms. First a 10-by-10 two-dimensional histogram is created with x and y running from 0 to 1. Then a few sample points are added to the histogram, at (0.3,0.3) with a height of 1, at (0.8,0.1) with a height of 5 and at (0.7,0.9) with a height of 0.5. This histogram with three events is used to generate a random sample of 1000 simulated events, which are printed out. #include #include #include int main (void) { const gsl_rng_type * T; gsl_rng * r; gsl_histogram2d * h = gsl_histogram2d_alloc (10, 10); gsl_histogram2d_set_ranges_uniform (h, 0.0, 1.0, 0.0, 1.0); gsl_histogram2d_accumulate (h, 0.3, 0.3, 1); gsl_histogram2d_accumulate (h, 0.8, 0.1, 5); gsl_histogram2d_accumulate (h, 0.7, 0.9, 0.5); gsl_rng_env_setup (); T = gsl_rng_default; r = gsl_rng_alloc (T); { int i; gsl_histogram2d_pdf * p = gsl_histogram2d_pdf_alloc (h->nx, h->ny); gsl_histogram2d_pdf_init (p, h); for (i = 0; i < 1000; i++) { double x, y; double u = gsl_rng_uniform (r); double v = gsl_rng_uniform (r); gsl_histogram2d_pdf_sample (p, u, v, &x, &y); printf ("%g %g\n", x, y); } gsl_histogram2d_pdf_free (p); } gsl_histogram2d_free (h); gsl_rng_free (r); return 0; }  File: gsl-ref.info, Node: N-tuples, Next: Monte Carlo Integration, Prev: Histograms, Up: Top 23 N-tuples *********** This chapter describes functions for creating and manipulating "ntuples", sets of values associated with events. The ntuples are stored in files. Their values can be extracted in any combination and "booked" in a histogram using a selection function. The values to be stored are held in a user-defined data structure, and an ntuple is created associating this data structure with a file. The values are then written to the file (normally inside a loop) using the ntuple functions described below. A histogram can be created from ntuple data by providing a selection function and a value function. The selection function specifies whether an event should be included in the subset to be analyzed or not. The value function computes the entry to be added to the histogram for each event. All the ntuple functions are defined in the header file 'gsl_ntuple.h' * Menu: * The ntuple struct:: * Creating ntuples:: * Opening an existing ntuple file:: * Writing ntuples:: * Reading ntuples :: * Closing an ntuple file:: * Histogramming ntuple values:: * Example ntuple programs:: * Ntuple References and Further Reading::  File: gsl-ref.info, Node: The ntuple struct, Next: Creating ntuples, Up: N-tuples 23.1 The ntuple struct ====================== Ntuples are manipulated using the 'gsl_ntuple' struct. This struct contains information on the file where the ntuple data is stored, a pointer to the current ntuple data row and the size of the user-defined ntuple data struct. typedef struct { FILE * file; void * ntuple_data; size_t size; } gsl_ntuple;  File: gsl-ref.info, Node: Creating ntuples, Next: Opening an existing ntuple file, Prev: The ntuple struct, Up: N-tuples 23.2 Creating ntuples ===================== -- Function: gsl_ntuple * gsl_ntuple_create (char * FILENAME, void * NTUPLE_DATA, size_t SIZE) This function creates a new write-only ntuple file FILENAME for ntuples of size SIZE and returns a pointer to the newly created ntuple struct. Any existing file with the same name is truncated to zero length and overwritten. A pointer to memory for the current ntuple row NTUPLE_DATA must be supplied--this is used to copy ntuples in and out of the file.  File: gsl-ref.info, Node: Opening an existing ntuple file, Next: Writing ntuples, Prev: Creating ntuples, Up: N-tuples 23.3 Opening an existing ntuple file ==================================== -- Function: gsl_ntuple * gsl_ntuple_open (char * FILENAME, void * NTUPLE_DATA, size_t SIZE) This function opens an existing ntuple file FILENAME for reading and returns a pointer to a corresponding ntuple struct. The ntuples in the file must have size SIZE. A pointer to memory for the current ntuple row NTUPLE_DATA must be supplied--this is used to copy ntuples in and out of the file.  File: gsl-ref.info, Node: Writing ntuples, Next: Reading ntuples, Prev: Opening an existing ntuple file, Up: N-tuples 23.4 Writing ntuples ==================== -- Function: int gsl_ntuple_write (gsl_ntuple * NTUPLE) This function writes the current ntuple NTUPLE->NTUPLE_DATA of size NTUPLE->SIZE to the corresponding file. -- Function: int gsl_ntuple_bookdata (gsl_ntuple * NTUPLE) This function is a synonym for 'gsl_ntuple_write'.  File: gsl-ref.info, Node: Reading ntuples, Next: Closing an ntuple file, Prev: Writing ntuples, Up: N-tuples 23.5 Reading ntuples ==================== -- Function: int gsl_ntuple_read (gsl_ntuple * NTUPLE) This function reads the current row of the ntuple file for NTUPLE and stores the values in NTUPLE->DATA.  File: gsl-ref.info, Node: Closing an ntuple file, Next: Histogramming ntuple values, Prev: Reading ntuples, Up: N-tuples 23.6 Closing an ntuple file =========================== -- Function: int gsl_ntuple_close (gsl_ntuple * NTUPLE) This function closes the ntuple file NTUPLE and frees its associated allocated memory.  File: gsl-ref.info, Node: Histogramming ntuple values, Next: Example ntuple programs, Prev: Closing an ntuple file, Up: N-tuples 23.7 Histogramming ntuple values ================================ Once an ntuple has been created its contents can be histogrammed in various ways using the function 'gsl_ntuple_project'. Two user-defined functions must be provided, a function to select events and a function to compute scalar values. The selection function and the value function both accept the ntuple row as a first argument and other parameters as a second argument. The "selection function" determines which ntuple rows are selected for histogramming. It is defined by the following struct, typedef struct { int (* function) (void * ntuple_data, void * params); void * params; } gsl_ntuple_select_fn; The struct component FUNCTION should return a non-zero value for each ntuple row that is to be included in the histogram. The "value function" computes scalar values for those ntuple rows selected by the selection function, typedef struct { double (* function) (void * ntuple_data, void * params); void * params; } gsl_ntuple_value_fn; In this case the struct component FUNCTION should return the value to be added to the histogram for the ntuple row. -- Function: int gsl_ntuple_project (gsl_histogram * H, gsl_ntuple * NTUPLE, gsl_ntuple_value_fn * VALUE_FUNC, gsl_ntuple_select_fn * SELECT_FUNC) This function updates the histogram H from the ntuple NTUPLE using the functions VALUE_FUNC and SELECT_FUNC. For each ntuple row where the selection function SELECT_FUNC is non-zero the corresponding value of that row is computed using the function VALUE_FUNC and added to the histogram. Those ntuple rows where SELECT_FUNC returns zero are ignored. New entries are added to the histogram, so subsequent calls can be used to accumulate further data in the same histogram.  File: gsl-ref.info, Node: Example ntuple programs, Next: Ntuple References and Further Reading, Prev: Histogramming ntuple values, Up: N-tuples 23.8 Examples ============= The following example programs demonstrate the use of ntuples in managing a large dataset. The first program creates a set of 10,000 simulated "events", each with 3 associated values (x,y,z). These are generated from a Gaussian distribution with unit variance, for demonstration purposes, and written to the ntuple file 'test.dat'. #include #include #include struct data { double x; double y; double z; }; int main (void) { const gsl_rng_type * T; gsl_rng * r; struct data ntuple_row; int i; gsl_ntuple *ntuple = gsl_ntuple_create ("test.dat", &ntuple_row, sizeof (ntuple_row)); gsl_rng_env_setup (); T = gsl_rng_default; r = gsl_rng_alloc (T); for (i = 0; i < 10000; i++) { ntuple_row.x = gsl_ran_ugaussian (r); ntuple_row.y = gsl_ran_ugaussian (r); ntuple_row.z = gsl_ran_ugaussian (r); gsl_ntuple_write (ntuple); } gsl_ntuple_close (ntuple); gsl_rng_free (r); return 0; } The next program analyses the ntuple data in the file 'test.dat'. The analysis procedure is to compute the squared-magnitude of each event, E^2=x^2+y^2+z^2, and select only those which exceed a lower limit of 1.5. The selected events are then histogrammed using their E^2 values. #include #include #include struct data { double x; double y; double z; }; int sel_func (void *ntuple_data, void *params); double val_func (void *ntuple_data, void *params); int main (void) { struct data ntuple_row; gsl_ntuple *ntuple = gsl_ntuple_open ("test.dat", &ntuple_row, sizeof (ntuple_row)); double lower = 1.5; gsl_ntuple_select_fn S; gsl_ntuple_value_fn V; gsl_histogram *h = gsl_histogram_alloc (100); gsl_histogram_set_ranges_uniform(h, 0.0, 10.0); S.function = &sel_func; S.params = &lower; V.function = &val_func; V.params = 0; gsl_ntuple_project (h, ntuple, &V, &S); gsl_histogram_fprintf (stdout, h, "%f", "%f"); gsl_histogram_free (h); gsl_ntuple_close (ntuple); return 0; } int sel_func (void *ntuple_data, void *params) { struct data * data = (struct data *) ntuple_data; double x, y, z, E2, scale; scale = *(double *) params; x = data->x; y = data->y; z = data->z; E2 = x * x + y * y + z * z; return E2 > scale; } double val_func (void *ntuple_data, void *params) { struct data * data = (struct data *) ntuple_data; double x, y, z; x = data->x; y = data->y; z = data->z; return x * x + y * y + z * z; } The following plot shows the distribution of the selected events. Note the cut-off at the lower bound.  File: gsl-ref.info, Node: Ntuple References and Further Reading, Prev: Example ntuple programs, Up: N-tuples 23.9 References and Further Reading =================================== Further information on the use of ntuples can be found in the documentation for the CERN packages PAW and HBOOK (available online).  File: gsl-ref.info, Node: Monte Carlo Integration, Next: Simulated Annealing, Prev: N-tuples, Up: Top 24 Monte Carlo Integration ************************** This chapter describes routines for multidimensional Monte Carlo integration. These include the traditional Monte Carlo method and adaptive algorithms such as VEGAS and MISER which use importance sampling and stratified sampling techniques. Each algorithm computes an estimate of a multidimensional definite integral of the form, I = \int_xl^xu dx \int_yl^yu dy ... f(x, y, ...) over a hypercubic region ((x_l,x_u), (y_l,y_u), ...) using a fixed number of function calls. The routines also provide a statistical estimate of the error on the result. This error estimate should be taken as a guide rather than as a strict error bound--random sampling of the region may not uncover all the important features of the function, resulting in an underestimate of the error. The functions are defined in separate header files for each routine, 'gsl_monte_plain.h', 'gsl_monte_miser.h' and 'gsl_monte_vegas.h'. * Menu: * Monte Carlo Interface:: * PLAIN Monte Carlo:: * MISER:: * VEGAS:: * Monte Carlo Examples:: * Monte Carlo Integration References and Further Reading::  File: gsl-ref.info, Node: Monte Carlo Interface, Next: PLAIN Monte Carlo, Up: Monte Carlo Integration 24.1 Interface ============== All of the Monte Carlo integration routines use the same general form of interface. There is an allocator to allocate memory for control variables and workspace, a routine to initialize those control variables, the integrator itself, and a function to free the space when done. Each integration function requires a random number generator to be supplied, and returns an estimate of the integral and its standard deviation. The accuracy of the result is determined by the number of function calls specified by the user. If a known level of accuracy is required this can be achieved by calling the integrator several times and averaging the individual results until the desired accuracy is obtained. Random sample points used within the Monte Carlo routines are always chosen strictly within the integration region, so that endpoint singularities are automatically avoided. The function to be integrated has its own datatype, defined in the header file 'gsl_monte.h'. -- Data Type: gsl_monte_function This data type defines a general function with parameters for Monte Carlo integration. 'double (* f) (double * X, size_t DIM, void * PARAMS)' this function should return the value f(x,params) for the argument X and parameters PARAMS, where X is an array of size DIM giving the coordinates of the point where the function is to be evaluated. 'size_t dim' the number of dimensions for X. 'void * params' a pointer to the parameters of the function. Here is an example for a quadratic function in two dimensions, f(x,y) = a x^2 + b x y + c y^2 with a = 3, b = 2, c = 1. The following code defines a 'gsl_monte_function' 'F' which you could pass to an integrator: struct my_f_params { double a; double b; double c; }; double my_f (double x[], size_t dim, void * p) { struct my_f_params * fp = (struct my_f_params *)p; if (dim != 2) { fprintf (stderr, "error: dim != 2"); abort (); } return fp->a * x[0] * x[0] + fp->b * x[0] * x[1] + fp->c * x[1] * x[1]; } gsl_monte_function F; struct my_f_params params = { 3.0, 2.0, 1.0 }; F.f = &my_f; F.dim = 2; F.params = ¶ms; The function f(x) can be evaluated using the following macro, #define GSL_MONTE_FN_EVAL(F,x) (*((F)->f))(x,(F)->dim,(F)->params)  File: gsl-ref.info, Node: PLAIN Monte Carlo, Next: MISER, Prev: Monte Carlo Interface, Up: Monte Carlo Integration 24.2 PLAIN Monte Carlo ====================== The plain Monte Carlo algorithm samples points randomly from the integration region to estimate the integral and its error. Using this algorithm the estimate of the integral E(f; N) for N randomly distributed points x_i is given by, E(f; N) = = V = (V / N) \sum_i^N f(x_i) where V is the volume of the integration region. The error on this estimate \sigma(E;N) is calculated from the estimated variance of the mean, \sigma^2 (E; N) = (V^2 / N^2) \sum_i^N (f(x_i) - )^2. For large N this variance decreases asymptotically as \Var(f)/N, where \Var(f) is the true variance of the function over the integration region. The error estimate itself should decrease as \sigma(f)/\sqrt{N}. The familiar law of errors decreasing as 1/\sqrt{N} applies--to reduce the error by a factor of 10 requires a 100-fold increase in the number of sample points. The functions described in this section are declared in the header file 'gsl_monte_plain.h'. -- Function: gsl_monte_plain_state * gsl_monte_plain_alloc (size_t DIM) This function allocates and initializes a workspace for Monte Carlo integration in DIM dimensions. -- Function: int gsl_monte_plain_init (gsl_monte_plain_state* S) This function initializes a previously allocated integration state. This allows an existing workspace to be reused for different integrations. -- Function: int gsl_monte_plain_integrate (gsl_monte_function * F, const double XL[], const double XU[], size_t DIM, size_t CALLS, gsl_rng * R, gsl_monte_plain_state * S, double * RESULT, double * ABSERR) This routines uses the plain Monte Carlo algorithm to integrate the function F over the DIM-dimensional hypercubic region defined by the lower and upper limits in the arrays XL and XU, each of size DIM. The integration uses a fixed number of function calls CALLS, and obtains random sampling points using the random number generator R. A previously allocated workspace S must be supplied. The result of the integration is returned in RESULT, with an estimated absolute error ABSERR. -- Function: void gsl_monte_plain_free (gsl_monte_plain_state * S) This function frees the memory associated with the integrator state S.  File: gsl-ref.info, Node: MISER, Next: VEGAS, Prev: PLAIN Monte Carlo, Up: Monte Carlo Integration 24.3 MISER ========== The MISER algorithm of Press and Farrar is based on recursive stratified sampling. This technique aims to reduce the overall integration error by concentrating integration points in the regions of highest variance. The idea of stratified sampling begins with the observation that for two disjoint regions a and b with Monte Carlo estimates of the integral E_a(f) and E_b(f) and variances \sigma_a^2(f) and \sigma_b^2(f), the variance \Var(f) of the combined estimate E(f) = (1/2) (E_a(f) + E_b(f)) is given by, \Var(f) = (\sigma_a^2(f) / 4 N_a) + (\sigma_b^2(f) / 4 N_b). It can be shown that this variance is minimized by distributing the points such that, N_a / (N_a + N_b) = \sigma_a / (\sigma_a + \sigma_b). Hence the smallest error estimate is obtained by allocating sample points in proportion to the standard deviation of the function in each sub-region. The MISER algorithm proceeds by bisecting the integration region along one coordinate axis to give two sub-regions at each step. The direction is chosen by examining all d possible bisections and selecting the one which will minimize the combined variance of the two sub-regions. The variance in the sub-regions is estimated by sampling with a fraction of the total number of points available to the current step. The same procedure is then repeated recursively for each of the two half-spaces from the best bisection. The remaining sample points are allocated to the sub-regions using the formula for N_a and N_b. This recursive allocation of integration points continues down to a user-specified depth where each sub-region is integrated using a plain Monte Carlo estimate. These individual values and their error estimates are then combined upwards to give an overall result and an estimate of its error. The functions described in this section are declared in the header file 'gsl_monte_miser.h'. -- Function: gsl_monte_miser_state * gsl_monte_miser_alloc (size_t DIM) This function allocates and initializes a workspace for Monte Carlo integration in DIM dimensions. The workspace is used to maintain the state of the integration. -- Function: int gsl_monte_miser_init (gsl_monte_miser_state* S) This function initializes a previously allocated integration state. This allows an existing workspace to be reused for different integrations. -- Function: int gsl_monte_miser_integrate (gsl_monte_function * F, const double XL[], const double XU[], size_t DIM, size_t CALLS, gsl_rng * R, gsl_monte_miser_state * S, double * RESULT, double * ABSERR) This routines uses the MISER Monte Carlo algorithm to integrate the function F over the DIM-dimensional hypercubic region defined by the lower and upper limits in the arrays XL and XU, each of size DIM. The integration uses a fixed number of function calls CALLS, and obtains random sampling points using the random number generator R. A previously allocated workspace S must be supplied. The result of the integration is returned in RESULT, with an estimated absolute error ABSERR. -- Function: void gsl_monte_miser_free (gsl_monte_miser_state * S) This function frees the memory associated with the integrator state S. The MISER algorithm has several configurable parameters which can be changed using the following two functions.(1) -- Function: void gsl_monte_miser_params_get (const gsl_monte_miser_state * S, gsl_monte_miser_params * PARAMS) This function copies the parameters of the integrator state into the user-supplied PARAMS structure. -- Function: void gsl_monte_miser_params_set (gsl_monte_miser_state * S, const gsl_monte_miser_params * PARAMS) This function sets the integrator parameters based on values provided in the PARAMS structure. Typically the values of the parameters are first read using 'gsl_monte_miser_params_get', the necessary changes are made to the fields of the PARAMS structure, and the values are copied back into the integrator state using 'gsl_monte_miser_params_set'. The functions use the 'gsl_monte_miser_params' structure which contains the following fields: -- Variable: double estimate_frac This parameter specifies the fraction of the currently available number of function calls which are allocated to estimating the variance at each recursive step. The default value is 0.1. -- Variable: size_t min_calls This parameter specifies the minimum number of function calls required for each estimate of the variance. If the number of function calls allocated to the estimate using ESTIMATE_FRAC falls below MIN_CALLS then MIN_CALLS are used instead. This ensures that each estimate maintains a reasonable level of accuracy. The default value of MIN_CALLS is '16 * dim'. -- Variable: size_t min_calls_per_bisection This parameter specifies the minimum number of function calls required to proceed with a bisection step. When a recursive step has fewer calls available than MIN_CALLS_PER_BISECTION it performs a plain Monte Carlo estimate of the current sub-region and terminates its branch of the recursion. The default value of this parameter is '32 * min_calls'. -- Variable: double alpha This parameter controls how the estimated variances for the two sub-regions of a bisection are combined when allocating points. With recursive sampling the overall variance should scale better than 1/N, since the values from the sub-regions will be obtained using a procedure which explicitly minimizes their variance. To accommodate this behavior the MISER algorithm allows the total variance to depend on a scaling parameter \alpha, \Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}. The authors of the original paper describing MISER recommend the value \alpha = 2 as a good choice, obtained from numerical experiments, and this is used as the default value in this implementation. -- Variable: double dither This parameter introduces a random fractional variation of size DITHER into each bisection, which can be used to break the symmetry of integrands which are concentrated near the exact center of the hypercubic integration region. The default value of dither is zero, so no variation is introduced. If needed, a typical value of DITHER is 0.1. ---------- Footnotes ---------- (1) The previous method of accessing these fields directly through the 'gsl_monte_miser_state' struct is now deprecated.  File: gsl-ref.info, Node: VEGAS, Next: Monte Carlo Examples, Prev: MISER, Up: Monte Carlo Integration 24.4 VEGAS ========== The VEGAS algorithm of Lepage is based on importance sampling. It samples points from the probability distribution described by the function |f|, so that the points are concentrated in the regions that make the largest contribution to the integral. In general, if the Monte Carlo integral of f is sampled with points distributed according to a probability distribution described by the function g, we obtain an estimate E_g(f; N), E_g(f; N) = E(f/g; N) with a corresponding variance, \Var_g(f; N) = \Var(f/g; N). If the probability distribution is chosen as g = |f|/I(|f|) then it can be shown that the variance V_g(f; N) vanishes, and the error in the estimate will be zero. In practice it is not possible to sample from the exact distribution g for an arbitrary function, so importance sampling algorithms aim to produce efficient approximations to the desired distribution. The VEGAS algorithm approximates the exact distribution by making a number of passes over the integration region while histogramming the function f. Each histogram is used to define a sampling distribution for the next pass. Asymptotically this procedure converges to the desired distribution. In order to avoid the number of histogram bins growing like K^d the probability distribution is approximated by a separable function: g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ... so that the number of bins required is only Kd. This is equivalent to locating the peaks of the function from the projections of the integrand onto the coordinate axes. The efficiency of VEGAS depends on the validity of this assumption. It is most efficient when the peaks of the integrand are well-localized. If an integrand can be rewritten in a form which is approximately separable this will increase the efficiency of integration with VEGAS. VEGAS incorporates a number of additional features, and combines both stratified sampling and importance sampling. The integration region is divided into a number of "boxes", with each box getting a fixed number of points (the goal is 2). Each box can then have a fractional number of bins, but if the ratio of bins-per-box is less than two, Vegas switches to a kind variance reduction (rather than importance sampling). -- Function: gsl_monte_vegas_state * gsl_monte_vegas_alloc (size_t DIM) This function allocates and initializes a workspace for Monte Carlo integration in DIM dimensions. The workspace is used to maintain the state of the integration. -- Function: int gsl_monte_vegas_init (gsl_monte_vegas_state* S) This function initializes a previously allocated integration state. This allows an existing workspace to be reused for different integrations. -- Function: int gsl_monte_vegas_integrate (gsl_monte_function * F, double XL[], double XU[], size_t DIM, size_t CALLS, gsl_rng * R, gsl_monte_vegas_state * S, double * RESULT, double * ABSERR) This routines uses the VEGAS Monte Carlo algorithm to integrate the function F over the DIM-dimensional hypercubic region defined by the lower and upper limits in the arrays XL and XU, each of size DIM. The integration uses a fixed number of function calls CALLS, and obtains random sampling points using the random number generator R. A previously allocated workspace S must be supplied. The result of the integration is returned in RESULT, with an estimated absolute error ABSERR. The result and its error estimate are based on a weighted average of independent samples. The chi-squared per degree of freedom for the weighted average is returned via the state struct component, S->CHISQ, and must be consistent with 1 for the weighted average to be reliable. -- Function: void gsl_monte_vegas_free (gsl_monte_vegas_state * S) This function frees the memory associated with the integrator state S. The VEGAS algorithm computes a number of independent estimates of the integral internally, according to the 'iterations' parameter described below, and returns their weighted average. Random sampling of the integrand can occasionally produce an estimate where the error is zero, particularly if the function is constant in some regions. An estimate with zero error causes the weighted average to break down and must be handled separately. In the original Fortran implementations of VEGAS the error estimate is made non-zero by substituting a small value (typically '1e-30'). The implementation in GSL differs from this and avoids the use of an arbitrary constant--it either assigns the value a weight which is the average weight of the preceding estimates or discards it according to the following procedure, current estimate has zero error, weighted average has finite error The current estimate is assigned a weight which is the average weight of the preceding estimates. current estimate has finite error, previous estimates had zero error The previous estimates are discarded and the weighted averaging procedure begins with the current estimate. current estimate has zero error, previous estimates had zero error The estimates are averaged using the arithmetic mean, but no error is computed. The convergence of the algorithm can be tested using the overall chi-squared value of the results, which is available from the following function: -- Function: double gsl_monte_vegas_chisq (const gsl_monte_vegas_state * S) This function returns the chi-squared per degree of freedom for the weighted estimate of the integral. The returned value should be close to 1. A value which differs significantly from 1 indicates that the values from different iterations are inconsistent. In this case the weighted error will be under-estimated, and further iterations of the algorithm are needed to obtain reliable results. -- Function: void gsl_monte_vegas_runval (const gsl_monte_vegas_state * S, double * RESULT, double * SIGMA) This function returns the raw (unaveraged) values of the integral RESULT and its error SIGMA from the most recent iteration of the algorithm. The VEGAS algorithm is highly configurable. Several parameters can be changed using the following two functions. -- Function: void gsl_monte_vegas_params_get (const gsl_monte_vegas_state * S, gsl_monte_vegas_params * PARAMS) This function copies the parameters of the integrator state into the user-supplied PARAMS structure. -- Function: void gsl_monte_vegas_params_set (gsl_monte_vegas_state * S, const gsl_monte_vegas_params * PARAMS) This function sets the integrator parameters based on values provided in the PARAMS structure. Typically the values of the parameters are first read using 'gsl_monte_vegas_params_get', the necessary changes are made to the fields of the PARAMS structure, and the values are copied back into the integrator state using 'gsl_monte_vegas_params_set'. The functions use the 'gsl_monte_vegas_params' structure which contains the following fields: -- Variable: double alpha The parameter 'alpha' controls the stiffness of the rebinning algorithm. It is typically set between one and two. A value of zero prevents rebinning of the grid. The default value is 1.5. -- Variable: size_t iterations The number of iterations to perform for each call to the routine. The default value is 5 iterations. -- Variable: int stage Setting this determines the "stage" of the calculation. Normally, 'stage = 0' which begins with a new uniform grid and empty weighted average. Calling VEGAS with 'stage = 1' retains the grid from the previous run but discards the weighted average, so that one can "tune" the grid using a relatively small number of points and then do a large run with 'stage = 1' on the optimized grid. Setting 'stage = 2' keeps the grid and the weighted average from the previous run, but may increase (or decrease) the number of histogram bins in the grid depending on the number of calls available. Choosing 'stage = 3' enters at the main loop, so that nothing is changed, and is equivalent to performing additional iterations in a previous call. -- Variable: int mode The possible choices are 'GSL_VEGAS_MODE_IMPORTANCE', 'GSL_VEGAS_MODE_STRATIFIED', 'GSL_VEGAS_MODE_IMPORTANCE_ONLY'. This determines whether VEGAS will use importance sampling or stratified sampling, or whether it can pick on its own. In low dimensions VEGAS uses strict stratified sampling (more precisely, stratified sampling is chosen if there are fewer than 2 bins per box). -- Variable: int verbose -- Variable: FILE * ostream These parameters set the level of information printed by VEGAS. All information is written to the stream OSTREAM. The default setting of VERBOSE is '-1', which turns off all output. A VERBOSE value of '0' prints summary information about the weighted average and final result, while a value of '1' also displays the grid coordinates. A value of '2' prints information from the rebinning procedure for each iteration. The above fields and the CHISQ value can also be accessed directly in the 'gsl_monte_vegas_state' but such use is deprecated.  File: gsl-ref.info, Node: Monte Carlo Examples, Next: Monte Carlo Integration References and Further Reading, Prev: VEGAS, Up: Monte Carlo Integration 24.5 Examples ============= The example program below uses the Monte Carlo routines to estimate the value of the following 3-dimensional integral from the theory of random walks, I = \int_{-pi}^{+pi} {dk_x/(2 pi)} \int_{-pi}^{+pi} {dk_y/(2 pi)} \int_{-pi}^{+pi} {dk_z/(2 pi)} 1 / (1 - cos(k_x)cos(k_y)cos(k_z)). The analytic value of this integral can be shown to be I = \Gamma(1/4)^4/(4 \pi^3) = 1.393203929685676859.... The integral gives the mean time spent at the origin by a random walk on a body-centered cubic lattice in three dimensions. For simplicity we will compute the integral over the region (0,0,0) to (\pi,\pi,\pi) and multiply by 8 to obtain the full result. The integral is slowly varying in the middle of the region but has integrable singularities at the corners (0,0,0), (0,\pi,\pi), (\pi,0,\pi) and (\pi,\pi,0). The Monte Carlo routines only select points which are strictly within the integration region and so no special measures are needed to avoid these singularities. #include #include #include #include #include #include /* Computation of the integral, I = int (dx dy dz)/(2pi)^3 1/(1-cos(x)cos(y)cos(z)) over (-pi,-pi,-pi) to (+pi, +pi, +pi). The exact answer is Gamma(1/4)^4/(4 pi^3). This example is taken from C.Itzykson, J.M.Drouffe, "Statistical Field Theory - Volume 1", Section 1.1, p21, which cites the original paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74 1800 (1977) */ /* For simplicity we compute the integral over the region (0,0,0) -> (pi,pi,pi) and multiply by 8 */ double exact = 1.3932039296856768591842462603255; double g (double *k, size_t dim, void *params) { double A = 1.0 / (M_PI * M_PI * M_PI); return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2])); } void display_results (char *title, double result, double error) { printf ("%s ==================\n", title); printf ("result = % .6f\n", result); printf ("sigma = % .6f\n", error); printf ("exact = % .6f\n", exact); printf ("error = % .6f = %.2g sigma\n", result - exact, fabs (result - exact) / error); } int main (void) { double res, err; double xl[3] = { 0, 0, 0 }; double xu[3] = { M_PI, M_PI, M_PI }; const gsl_rng_type *T; gsl_rng *r; gsl_monte_function G = { &g, 3, 0 }; size_t calls = 500000; gsl_rng_env_setup (); T = gsl_rng_default; r = gsl_rng_alloc (T); { gsl_monte_plain_state *s = gsl_monte_plain_alloc (3); gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s, &res, &err); gsl_monte_plain_free (s); display_results ("plain", res, err); } { gsl_monte_miser_state *s = gsl_monte_miser_alloc (3); gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s, &res, &err); gsl_monte_miser_free (s); display_results ("miser", res, err); } { gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3); gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s, &res, &err); display_results ("vegas warm-up", res, err); printf ("converging...\n"); do { gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s, &res, &err); printf ("result = % .6f sigma = % .6f " "chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s)); } while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5); display_results ("vegas final", res, err); gsl_monte_vegas_free (s); } gsl_rng_free (r); return 0; } With 500,000 function calls the plain Monte Carlo algorithm achieves a fractional error of 1%. The estimated error 'sigma' is roughly consistent with the actual error-the computed result differs from the true result by about 1.4 standard deviations, plain ================== result = 1.412209 sigma = 0.013436 exact = 1.393204 error = 0.019005 = 1.4 sigma The MISER algorithm reduces the error by a factor of four, and also correctly estimates the error, miser ================== result = 1.391322 sigma = 0.003461 exact = 1.393204 error = -0.001882 = 0.54 sigma In the case of the VEGAS algorithm the program uses an initial warm-up run of 10,000 function calls to prepare, or "warm up", the grid. This is followed by a main run with five iterations of 100,000 function calls. The chi-squared per degree of freedom for the five iterations are checked for consistency with 1, and the run is repeated if the results have not converged. In this case the estimates are consistent on the first pass. vegas warm-up ================== result = 1.392673 sigma = 0.003410 exact = 1.393204 error = -0.000531 = 0.16 sigma converging... result = 1.393281 sigma = 0.000362 chisq/dof = 1.5 vegas final ================== result = 1.393281 sigma = 0.000362 exact = 1.393204 error = 0.000077 = 0.21 sigma If the value of 'chisq' had differed significantly from 1 it would indicate inconsistent results, with a correspondingly underestimated error. The final estimate from VEGAS (using a similar number of function calls) is significantly more accurate than the other two algorithms.  File: gsl-ref.info, Node: Monte Carlo Integration References and Further Reading, Prev: Monte Carlo Examples, Up: Monte Carlo Integration 24.6 References and Further Reading =================================== The MISER algorithm is described in the following article by Press and Farrar, W.H. Press, G.R. Farrar, 'Recursive Stratified Sampling for Multidimensional Monte Carlo Integration', Computers in Physics, v4 (1990), pp190-195. The VEGAS algorithm is described in the following papers, G.P. Lepage, 'A New Algorithm for Adaptive Multidimensional Integration', Journal of Computational Physics 27, 192-203, (1978) G.P. Lepage, 'VEGAS: An Adaptive Multi-dimensional Integration Program', Cornell preprint CLNS 80-447, March 1980  File: gsl-ref.info, Node: Simulated Annealing, Next: Ordinary Differential Equations, Prev: Monte Carlo Integration, Up: Top 25 Simulated Annealing ********************** Stochastic search techniques are used when the structure of a space is not well understood or is not smooth, so that techniques like Newton's method (which requires calculating Jacobian derivative matrices) cannot be used. In particular, these techniques are frequently used to solve combinatorial optimization problems, such as the traveling salesman problem. The goal is to find a point in the space at which a real valued "energy function" (or "cost function") is minimized. Simulated annealing is a minimization technique which has given good results in avoiding local minima; it is based on the idea of taking a random walk through the space at successively lower temperatures, where the probability of taking a step is given by a Boltzmann distribution. The functions described in this chapter are declared in the header file 'gsl_siman.h'. * Menu: * Simulated Annealing algorithm:: * Simulated Annealing functions:: * Examples with Simulated Annealing:: * Simulated Annealing References and Further Reading::  File: gsl-ref.info, Node: Simulated Annealing algorithm, Next: Simulated Annealing functions, Up: Simulated Annealing 25.1 Simulated Annealing algorithm ================================== The simulated annealing algorithm takes random walks through the problem space, looking for points with low energies; in these random walks, the probability of taking a step is determined by the Boltzmann distribution, p = e^{-(E_{i+1} - E_i)/(kT)} if E_{i+1} > E_i, and p = 1 when E_{i+1} <= E_i. In other words, a step will occur if the new energy is lower. If the new energy is higher, the transition can still occur, and its likelihood is proportional to the temperature T and inversely proportional to the energy difference E_{i+1} - E_i. The temperature T is initially set to a high value, and a random walk is carried out at that temperature. Then the temperature is lowered very slightly according to a "cooling schedule", for example: T -> T/mu_T where \mu_T is slightly greater than 1. The slight probability of taking a step that gives higher energy is what allows simulated annealing to frequently get out of local minima.  File: gsl-ref.info, Node: Simulated Annealing functions, Next: Examples with Simulated Annealing, Prev: Simulated Annealing algorithm, Up: Simulated Annealing 25.2 Simulated Annealing functions ================================== -- Function: void gsl_siman_solve (const gsl_rng * R, void * X0_P, gsl_siman_Efunc_t EF, gsl_siman_step_t TAKE_STEP, gsl_siman_metric_t DISTANCE, gsl_siman_print_t PRINT_POSITION, gsl_siman_copy_t COPYFUNC, gsl_siman_copy_construct_t COPY_CONSTRUCTOR, gsl_siman_destroy_t DESTRUCTOR, size_t ELEMENT_SIZE, gsl_siman_params_t PARAMS) This function performs a simulated annealing search through a given space. The space is specified by providing the functions EF and DISTANCE. The simulated annealing steps are generated using the random number generator R and the function TAKE_STEP. The starting configuration of the system should be given by X0_P. The routine offers two modes for updating configurations, a fixed-size mode and a variable-size mode. In the fixed-size mode the configuration is stored as a single block of memory of size ELEMENT_SIZE. Copies of this configuration are created, copied and destroyed internally using the standard library functions 'malloc', 'memcpy' and 'free'. The function pointers COPYFUNC, COPY_CONSTRUCTOR and DESTRUCTOR should be null pointers in fixed-size mode. In the variable-size mode the functions COPYFUNC, COPY_CONSTRUCTOR and DESTRUCTOR are used to create, copy and destroy configurations internally. The variable ELEMENT_SIZE should be zero in the variable-size mode. The PARAMS structure (described below) controls the run by providing the temperature schedule and other tunable parameters to the algorithm. On exit the best result achieved during the search is placed in '*X0_P'. If the annealing process has been successful this should be a good approximation to the optimal point in the space. If the function pointer PRINT_POSITION is not null, a debugging log will be printed to 'stdout' with the following columns: #-iter #-evals temperature position energy best_energy and the output of the function PRINT_POSITION itself. If PRINT_POSITION is null then no information is printed. The simulated annealing routines require several user-specified functions to define the configuration space and energy function. The prototypes for these functions are given below. -- Data Type: gsl_siman_Efunc_t This function type should return the energy of a configuration XP. double (*gsl_siman_Efunc_t) (void *xp) -- Data Type: gsl_siman_step_t This function type should modify the configuration XP using a random step taken from the generator R, up to a maximum distance of STEP_SIZE. void (*gsl_siman_step_t) (const gsl_rng *r, void *xp, double step_size) -- Data Type: gsl_siman_metric_t This function type should return the distance between two configurations XP and YP. double (*gsl_siman_metric_t) (void *xp, void *yp) -- Data Type: gsl_siman_print_t This function type should print the contents of the configuration XP. void (*gsl_siman_print_t) (void *xp) -- Data Type: gsl_siman_copy_t This function type should copy the configuration SOURCE into DEST. void (*gsl_siman_copy_t) (void *source, void *dest) -- Data Type: gsl_siman_copy_construct_t This function type should create a new copy of the configuration XP. void * (*gsl_siman_copy_construct_t) (void *xp) -- Data Type: gsl_siman_destroy_t This function type should destroy the configuration XP, freeing its memory. void (*gsl_siman_destroy_t) (void *xp) -- Data Type: gsl_siman_params_t These are the parameters that control a run of 'gsl_siman_solve'. This structure contains all the information needed to control the search, beyond the energy function, the step function and the initial guess. 'int n_tries' The number of points to try for each step. 'int iters_fixed_T' The number of iterations at each temperature. 'double step_size' The maximum step size in the random walk. 'double k, t_initial, mu_t, t_min' The parameters of the Boltzmann distribution and cooling schedule.  File: gsl-ref.info, Node: Examples with Simulated Annealing, Next: Simulated Annealing References and Further Reading, Prev: Simulated Annealing functions, Up: Simulated Annealing 25.3 Examples ============= The simulated annealing package is clumsy, and it has to be because it is written in C, for C callers, and tries to be polymorphic at the same time. But here we provide some examples which can be pasted into your application with little change and should make things easier. * Menu: * Trivial example:: * Traveling Salesman Problem::  File: gsl-ref.info, Node: Trivial example, Next: Traveling Salesman Problem, Up: Examples with Simulated Annealing 25.3.1 Trivial example ---------------------- The first example, in one dimensional Cartesian space, sets up an energy function which is a damped sine wave; this has many local minima, but only one global minimum, somewhere between 1.0 and 1.5. The initial guess given is 15.5, which is several local minima away from the global minimum. #include #include #include #include /* set up parameters for this simulated annealing run */ /* how many points do we try before stepping */ #define N_TRIES 200 /* how many iterations for each T? */ #define ITERS_FIXED_T 1000 /* max step size in random walk */ #define STEP_SIZE 1.0 /* Boltzmann constant */ #define K 1.0 /* initial temperature */ #define T_INITIAL 0.008 /* damping factor for temperature */ #define MU_T 1.003 #define T_MIN 2.0e-6 gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE, K, T_INITIAL, MU_T, T_MIN}; /* now some functions to test in one dimension */ double E1(void *xp) { double x = * ((double *) xp); return exp(-pow((x-1.0),2.0))*sin(8*x); } double M1(void *xp, void *yp) { double x = *((double *) xp); double y = *((double *) yp); return fabs(x - y); } void S1(const gsl_rng * r, void *xp, double step_size) { double old_x = *((double *) xp); double new_x; double u = gsl_rng_uniform(r); new_x = u * 2 * step_size - step_size + old_x; memcpy(xp, &new_x, sizeof(new_x)); } void P1(void *xp) { printf ("%12g", *((double *) xp)); } int main(int argc, char *argv[]) { const gsl_rng_type * T; gsl_rng * r; double x_initial = 15.5; gsl_rng_env_setup(); T = gsl_rng_default; r = gsl_rng_alloc(T); gsl_siman_solve(r, &x_initial, E1, S1, M1, P1, NULL, NULL, NULL, sizeof(double), params); gsl_rng_free (r); return 0; } Here are a couple of plots that are generated by running 'siman_test' in the following way: $ ./siman_test | awk '!/^#/ {print $1, $4}' | graph -y 1.34 1.4 -W0 -X generation -Y position | plot -Tps > siman-test.eps $ ./siman_test | awk '!/^#/ {print $1, $5}' | graph -y -0.88 -0.83 -W0 -X generation -Y energy | plot -Tps > siman-energy.eps  File: gsl-ref.info, Node: Traveling Salesman Problem, Prev: Trivial example, Up: Examples with Simulated Annealing 25.3.2 Traveling Salesman Problem --------------------------------- The TSP ("Traveling Salesman Problem") is the classic combinatorial optimization problem. I have provided a very simple version of it, based on the coordinates of twelve cities in the southwestern United States. This should maybe be called the "Flying Salesman Problem", since I am using the great-circle distance between cities, rather than the driving distance. Also: I assume the earth is a sphere, so I don't use geoid distances. The 'gsl_siman_solve' routine finds a route which is 3490.62 Kilometers long; this is confirmed by an exhaustive search of all possible routes with the same initial city. The full code can be found in 'siman/siman_tsp.c', but I include here some plots generated in the following way: $ ./siman_tsp > tsp.output $ grep -v "^#" tsp.output | awk '{print $1, $NF}' | graph -y 3300 6500 -W0 -X generation -Y distance -L "TSP - 12 southwest cities" | plot -Tps > 12-cities.eps $ grep initial_city_coord tsp.output | awk '{print $2, $3}' | graph -X "longitude (- means west)" -Y "latitude" -L "TSP - initial-order" -f 0.03 -S 1 0.1 | plot -Tps > initial-route.eps $ grep final_city_coord tsp.output | awk '{print $2, $3}' | graph -X "longitude (- means west)" -Y "latitude" -L "TSP - final-order" -f 0.03 -S 1 0.1 | plot -Tps > final-route.eps This is the output showing the initial order of the cities; longitude is negative, since it is west and I want the plot to look like a map. # initial coordinates of cities (longitude and latitude) ###initial_city_coord: -105.95 35.68 Santa Fe ###initial_city_coord: -112.07 33.54 Phoenix ###initial_city_coord: -106.62 35.12 Albuquerque ###initial_city_coord: -103.2 34.41 Clovis ###initial_city_coord: -107.87 37.29 Durango ###initial_city_coord: -96.77 32.79 Dallas ###initial_city_coord: -105.92 35.77 Tesuque ###initial_city_coord: -107.84 35.15 Grants ###initial_city_coord: -106.28 35.89 Los Alamos ###initial_city_coord: -106.76 32.34 Las Cruces ###initial_city_coord: -108.58 37.35 Cortez ###initial_city_coord: -108.74 35.52 Gallup ###initial_city_coord: -105.95 35.68 Santa Fe The optimal route turns out to be: # final coordinates of cities (longitude and latitude) ###final_city_coord: -105.95 35.68 Santa Fe ###final_city_coord: -103.2 34.41 Clovis ###final_city_coord: -96.77 32.79 Dallas ###final_city_coord: -106.76 32.34 Las Cruces ###final_city_coord: -112.07 33.54 Phoenix ###final_city_coord: -108.74 35.52 Gallup ###final_city_coord: -108.58 37.35 Cortez ###final_city_coord: -107.87 37.29 Durango ###final_city_coord: -107.84 35.15 Grants ###final_city_coord: -106.62 35.12 Albuquerque ###final_city_coord: -106.28 35.89 Los Alamos ###final_city_coord: -105.92 35.77 Tesuque ###final_city_coord: -105.95 35.68 Santa Fe Here's a plot of the cost function (energy) versus generation (point in the calculation at which a new temperature is set) for this problem:  File: gsl-ref.info, Node: Simulated Annealing References and Further Reading, Prev: Examples with Simulated Annealing, Up: Simulated Annealing 25.4 References and Further Reading =================================== Further information is available in the following book, 'Modern Heuristic Techniques for Combinatorial Problems', Colin R. Reeves (ed.), McGraw-Hill, 1995 (ISBN 0-07-709239-2).  File: gsl-ref.info, Node: Ordinary Differential Equations, Next: Interpolation, Prev: Simulated Annealing, Up: Top 26 Ordinary Differential Equations ********************************** This chapter describes functions for solving ordinary differential equation (ODE) initial value problems. The library provides a variety of low-level methods, such as Runge-Kutta and Bulirsch-Stoer routines, and higher-level components for adaptive step-size control. The components can be combined by the user to achieve the desired solution, with full access to any intermediate steps. A driver object can be used as a high level wrapper for easy use of low level functions. These functions are declared in the header file 'gsl_odeiv2.h'. This is a new interface in version 1.15 and uses the prefix 'gsl_odeiv2' for all functions. It is recommended over the previous 'gsl_odeiv' implementation defined in 'gsl_odeiv.h' The old interface has been retained under the original name for backwards compatibility. * Menu: * Defining the ODE System:: * Stepping Functions:: * Adaptive Step-size Control:: * Evolution:: * Driver:: * ODE Example programs:: * ODE References and Further Reading::  File: gsl-ref.info, Node: Defining the ODE System, Next: Stepping Functions, Up: Ordinary Differential Equations 26.1 Defining the ODE System ============================ The routines solve the general n-dimensional first-order system, dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t)) for i = 1, \dots, n. The stepping functions rely on the vector of derivatives f_i and the Jacobian matrix, J_{ij} = df_i(t,y(t)) / dy_j. A system of equations is defined using the 'gsl_odeiv2_system' datatype. -- Data Type: gsl_odeiv2_system This data type defines a general ODE system with arbitrary parameters. 'int (* function) (double t, const double y[], double dydt[], void * params)' This function should store the vector elements f_i(t,y,params) in the array DYDT, for arguments (T,Y) and parameters PARAMS. The function should return 'GSL_SUCCESS' if the calculation was completed successfully. Any other return value indicates an error. A special return value 'GSL_EBADFUNC' causes 'gsl_odeiv2' routines to immediately stop and return. The user must call an appropriate reset function (e.g. 'gsl_odeiv2_driver_reset' or 'gsl_odeiv2_step_reset') before continuing. Use return values distinct from standard GSL error codes to distinguish your function as the source of the error. 'int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);' This function should store the vector of derivative elements in the array DFDT and the Jacobian matrix J_{ij} in the array DFDY, regarded as a row-ordered matrix 'J(i,j) = dfdy[i * dimension + j]' where 'dimension' is the dimension of the system. Not all of the stepper algorithms of 'gsl_odeiv2' make use of the Jacobian matrix, so it may not be necessary to provide this function (the 'jacobian' element of the struct can be replaced by a null pointer for those algorithms). The function should return 'GSL_SUCCESS' if the calculation was completed successfully. Any other return value indicates an error. A special return value 'GSL_EBADFUNC' causes 'gsl_odeiv2' routines to immediately stop and return. The user must call an appropriate reset function (e.g. 'gsl_odeiv2_driver_reset' or 'gsl_odeiv2_step_reset') before continuing. Use return values distinct from standard GSL error codes to distinguish your function as the source of the error. 'size_t dimension;' This is the dimension of the system of equations. 'void * params' This is a pointer to the arbitrary parameters of the system.  File: gsl-ref.info, Node: Stepping Functions, Next: Adaptive Step-size Control, Prev: Defining the ODE System, Up: Ordinary Differential Equations 26.2 Stepping Functions ======================= The lowest level components are the "stepping functions" which advance a solution from time t to t+h for a fixed step-size h and estimate the resulting local error. -- Function: gsl_odeiv2_step * gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T, size_t DIM) This function returns a pointer to a newly allocated instance of a stepping function of type T for a system of DIM dimensions. Please note that if you use a stepper method that requires access to a driver object, it is advisable to use a driver allocation method, which automatically allocates a stepper, too. -- Function: int gsl_odeiv2_step_reset (gsl_odeiv2_step * S) This function resets the stepping function S. It should be used whenever the next use of S will not be a continuation of a previous step. -- Function: void gsl_odeiv2_step_free (gsl_odeiv2_step * S) This function frees all the memory associated with the stepping function S. -- Function: const char * gsl_odeiv2_step_name (const gsl_odeiv2_step * S) This function returns a pointer to the name of the stepping function. For example, printf ("step method is '%s'\n", gsl_odeiv2_step_name (s)); would print something like 'step method is 'rkf45''. -- Function: unsigned int gsl_odeiv2_step_order (const gsl_odeiv2_step * S) This function returns the order of the stepping function on the previous step. The order can vary if the stepping function itself is adaptive. -- Function: int gsl_odeiv2_step_set_driver (gsl_odeiv2_step * S, const gsl_odeiv2_driver * D) This function sets a pointer of the driver object D for stepper S, to allow the stepper to access control (and evolve) object through the driver object. This is a requirement for some steppers, to get the desired error level for internal iteration of stepper. Allocation of a driver object calls this function automatically. -- Function: int gsl_odeiv2_step_apply (gsl_odeiv2_step * S, double T, double H, double Y[], double YERR[], const double DYDT_IN[], double DYDT_OUT[], const gsl_odeiv2_system * SYS) This function applies the stepping function S to the system of equations defined by SYS, using the step-size H to advance the system from time T and state Y to time T+H. The new state of the system is stored in Y on output, with an estimate of the absolute error in each component stored in YERR. If the argument DYDT_IN is not null it should point an array containing the derivatives for the system at time T on input. This is optional as the derivatives will be computed internally if they are not provided, but allows the reuse of existing derivative information. On output the new derivatives of the system at time T+H will be stored in DYDT_OUT if it is not null. The stepping function returns 'GSL_FAILURE' if it is unable to compute the requested step. Also, if the user-supplied functions defined in the system SYS return a status other than 'GSL_SUCCESS' the step will be aborted. In that case, the elements of Y will be restored to their pre-step values and the error code from the user-supplied function will be returned. Failure may be due to a singularity in the system or too large step-size H. In that case the step should be attempted again with a smaller step-size, e.g. H/2. If the driver object is not appropriately set via 'gsl_odeiv2_step_set_driver' for those steppers that need it, the stepping function returns 'GSL_EFAULT'. If the user-supplied functions defined in the system SYS returns 'GSL_EBADFUNC', the function returns immediately with the same return code. In this case the user must call 'gsl_odeiv2_step_reset' before calling this function again. The following algorithms are available, -- Step Type: gsl_odeiv2_step_rk2 Explicit embedded Runge-Kutta (2, 3) method. -- Step Type: gsl_odeiv2_step_rk4 Explicit 4th order (classical) Runge-Kutta. Error estimation is carried out by the step doubling method. For more efficient estimate of the error, use the embedded methods described below. -- Step Type: gsl_odeiv2_step_rkf45 Explicit embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. -- Step Type: gsl_odeiv2_step_rkck Explicit embedded Runge-Kutta Cash-Karp (4, 5) method. -- Step Type: gsl_odeiv2_step_rk8pd Explicit embedded Runge-Kutta Prince-Dormand (8, 9) method. -- Step Type: gsl_odeiv2_step_rk1imp Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method. This algorithm requires the Jacobian and access to the driver object via 'gsl_odeiv2_step_set_driver'. -- Step Type: gsl_odeiv2_step_rk2imp Implicit Gaussian second order Runge-Kutta. Also known as implicit mid-point rule. Error estimation is carried out by the step doubling method. This stepper requires the Jacobian and access to the driver object via 'gsl_odeiv2_step_set_driver'. -- Step Type: gsl_odeiv2_step_rk4imp Implicit Gaussian 4th order Runge-Kutta. Error estimation is carried out by the step doubling method. This algorithm requires the Jacobian and access to the driver object via 'gsl_odeiv2_step_set_driver'. -- Step Type: gsl_odeiv2_step_bsimp Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems. This stepper requires the Jacobian. -- Step Type: gsl_odeiv2_step_msadams A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. This stepper requires the access to the driver object via 'gsl_odeiv2_step_set_driver'. -- Step Type: gsl_odeiv2_step_msbdf A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. This stepper requires the Jacobian and the access to the driver object via 'gsl_odeiv2_step_set_driver'.  File: gsl-ref.info, Node: Adaptive Step-size Control, Next: Evolution, Prev: Stepping Functions, Up: Ordinary Differential Equations 26.3 Adaptive Step-size Control =============================== The control function examines the proposed change to the solution produced by a stepping function and attempts to determine the optimal step-size for a user-specified level of error. -- Function: gsl_odeiv2_control * gsl_odeiv2_control_standard_new (double EPS_ABS, double EPS_REL, double A_Y, double A_DYDT) The standard control object is a four parameter heuristic based on absolute and relative errors EPS_ABS and EPS_REL, and scaling factors A_Y and A_DYDT for the system state y(t) and derivatives y'(t) respectively. The step-size adjustment procedure for this method begins by computing the desired error level D_i for each component, D_i = eps_abs + eps_rel * (a_y |y_i| + a_dydt h |y\prime_i|) and comparing it with the observed error E_i = |yerr_i|. If the observed error E exceeds the desired error level D by more than 10% for any component then the method reduces the step-size by an appropriate factor, h_new = h_old * S * (E/D)^(-1/q) where q is the consistency order of the method (e.g. q=4 for 4(5) embedded RK), and S is a safety factor of 0.9. The ratio E/D is taken to be the maximum of the ratios E_i/D_i. If the observed error E is less than 50% of the desired error level D for the maximum ratio E_i/D_i then the algorithm takes the opportunity to increase the step-size to bring the error in line with the desired level, h_new = h_old * S * (E/D)^(-1/(q+1)) This encompasses all the standard error scaling methods. To avoid uncontrolled changes in the stepsize, the overall scaling factor is limited to the range 1/5 to 5. -- Function: gsl_odeiv2_control * gsl_odeiv2_control_y_new (double EPS_ABS, double EPS_REL) This function creates a new control object which will keep the local error on each step within an absolute error of EPS_ABS and relative error of EPS_REL with respect to the solution y_i(t). This is equivalent to the standard control object with A_Y=1 and A_DYDT=0. -- Function: gsl_odeiv2_control * gsl_odeiv2_control_yp_new (double EPS_ABS, double EPS_REL) This function creates a new control object which will keep the local error on each step within an absolute error of EPS_ABS and relative error of EPS_REL with respect to the derivatives of the solution y'_i(t). This is equivalent to the standard control object with A_Y=0 and A_DYDT=1. -- Function: gsl_odeiv2_control * gsl_odeiv2_control_scaled_new (double EPS_ABS, double EPS_REL, double A_Y, double A_DYDT, const double SCALE_ABS[], size_t DIM) This function creates a new control object which uses the same algorithm as 'gsl_odeiv2_control_standard_new' but with an absolute error which is scaled for each component by the array SCALE_ABS. The formula for D_i for this control object is, D_i = eps_abs * s_i + eps_rel * (a_y |y_i| + a_dydt h |y\prime_i|) where s_i is the i-th component of the array SCALE_ABS. The same error control heuristic is used by the Matlab ODE suite. -- Function: gsl_odeiv2_control * gsl_odeiv2_control_alloc (const gsl_odeiv2_control_type * T) This function returns a pointer to a newly allocated instance of a control function of type T. This function is only needed for defining new types of control functions. For most purposes the standard control functions described above should be sufficient. -- Function: int gsl_odeiv2_control_init (gsl_odeiv2_control * C, double EPS_ABS, double EPS_REL, double A_Y, double A_DYDT) This function initializes the control function C with the parameters EPS_ABS (absolute error), EPS_REL (relative error), A_Y (scaling factor for y) and A_DYDT (scaling factor for derivatives). -- Function: void gsl_odeiv2_control_free (gsl_odeiv2_control * C) This function frees all the memory associated with the control function C. -- Function: int gsl_odeiv2_control_hadjust (gsl_odeiv2_control * C, gsl_odeiv2_step * S, const double Y[], const double YERR[], const double DYDT[], double * H) This function adjusts the step-size H using the control function C, and the current values of Y, YERR and DYDT. The stepping function STEP is also needed to determine the order of the method. If the error in the y-values YERR is found to be too large then the step-size H is reduced and the function returns 'GSL_ODEIV_HADJ_DEC'. If the error is sufficiently small then H may be increased and 'GSL_ODEIV_HADJ_INC' is returned. The function returns 'GSL_ODEIV_HADJ_NIL' if the step-size is unchanged. The goal of the function is to estimate the largest step-size which satisfies the user-specified accuracy requirements for the current point. -- Function: const char * gsl_odeiv2_control_name (const gsl_odeiv2_control * C) This function returns a pointer to the name of the control function. For example, printf ("control method is '%s'\n", gsl_odeiv2_control_name (c)); would print something like 'control method is 'standard'' -- Function: int gsl_odeiv2_control_errlevel (gsl_odeiv2_control * C, const double Y, const double DYDT, const double H, const size_t IND, double * ERRLEV) This function calculates the desired error level of the IND-th component to ERRLEV. It requires the value (Y) and value of the derivative (DYDT) of the component, and the current step size H. -- Function: int gsl_odeiv2_control_set_driver (gsl_odeiv2_control * C, const gsl_odeiv2_driver * D) This function sets a pointer of the driver object D for control object C.  File: gsl-ref.info, Node: Evolution, Next: Driver, Prev: Adaptive Step-size Control, Up: Ordinary Differential Equations 26.4 Evolution ============== The evolution function combines the results of a stepping function and control function to reliably advance the solution forward one step using an acceptable step-size. -- Function: gsl_odeiv2_evolve * gsl_odeiv2_evolve_alloc (size_t DIM) This function returns a pointer to a newly allocated instance of an evolution function for a system of DIM dimensions. -- Function: int gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * E, gsl_odeiv2_control * CON, gsl_odeiv2_step * STEP, const gsl_odeiv2_system * SYS, double * T, double T1, double * H, double Y[]) This function advances the system (E, SYS) from time T and position Y using the stepping function STEP. The new time and position are stored in T and Y on output. The initial step-size is taken as H. The control function CON is applied to check whether the local error estimated by the stepping function STEP using step-size H exceeds the required error tolerance. If the error is too high, the step is retried by calling STEP with a decreased step-size. This process is continued until an acceptable step-size is found. An estimate of the local error for the step can be obtained from the components of the array 'E->yerr[]'. If the user-supplied functions defined in the system SYS returns 'GSL_EBADFUNC', the function returns immediately with the same return code. In this case the user must call 'gsl_odeiv2_step_reset' and 'gsl_odeiv2_evolve_reset' before calling this function again. Otherwise, if the user-supplied functions defined in the system SYS or the stepping function STEP return a status other than 'GSL_SUCCESS', the step is retried with a decreased step-size. If the step-size decreases below machine precision, a status of 'GSL_FAILURE' is returned if the user functions returned 'GSL_SUCCESS'. Otherwise the value returned by user function is returned. If no acceptable step can be made, T and Y will be restored to their pre-step values and H contains the final attempted step-size. If the step is successful the function returns a suggested step-size for the next step in H. The maximum time T1 is guaranteed not to be exceeded by the time-step. On the final time-step the value of T will be set to T1 exactly. -- Function: int gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * E, gsl_odeiv2_control * CON, gsl_odeiv2_step * STEP, const gsl_odeiv2_system * SYS, double * T, const double H, double Y[]) This function advances the ODE-system (E, SYS, CON) from time T and position Y using the stepping function STEP by a specified step size H. If the local error estimated by the stepping function exceeds the desired error level, the step is not taken and the function returns 'GSL_FAILURE'. Otherwise the value returned by user function is returned. -- Function: int gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * E) This function resets the evolution function E. It should be used whenever the next use of E will not be a continuation of a previous step. -- Function: void gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * E) This function frees all the memory associated with the evolution function E. -- Function: int gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * E, const gsl_odeiv2_driver * D) This function sets a pointer of the driver object D for evolve object E. If a system has discontinuous changes in the derivatives at known points, it is advisable to evolve the system between each discontinuity in sequence. For example, if a step-change in an external driving force occurs at times t_a, t_b and t_c then evolution should be carried out over the ranges (t_0,t_a), (t_a,t_b), (t_b,t_c), and (t_c,t_1) separately and not directly over the range (t_0,t_1).  File: gsl-ref.info, Node: Driver, Next: ODE Example programs, Prev: Evolution, Up: Ordinary Differential Equations 26.5 Driver =========== The driver object is a high level wrapper that combines the evolution, control and stepper objects for easy use. -- Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * SYS, const gsl_odeiv2_step_type * T, const double HSTART, const double EPSABS, const double EPSREL) -- Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * SYS, const gsl_odeiv2_step_type * T, const double HSTART, const double EPSABS, const double EPSREL) -- Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * SYS, const gsl_odeiv2_step_type * T, const double HSTART, const double EPSABS, const double EPSREL, const double A_Y, const double A_DYDT) -- Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * SYS, const gsl_odeiv2_step_type * T, const double HSTART, const double EPSABS, const double EPSREL, const double A_Y, const double A_DYDT, const double SCALE_ABS[]) These functions return a pointer to a newly allocated instance of a driver object. The functions automatically allocate and initialise the evolve, control and stepper objects for ODE system SYS using stepper type T. The initial step size is given in HSTART. The rest of the arguments follow the syntax and semantics of the control functions with same name ('gsl_odeiv2_control_*_new'). -- Function: int gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * D, const double HMIN) The function sets a minimum for allowed step size HMIN for driver D. Default value is 0. -- Function: int gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * D, const double HMAX) The function sets a maximum for allowed step size HMAX for driver D. Default value is 'GSL_DBL_MAX'. -- Function: int gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * D, const unsigned long int NMAX) The function sets a maximum for allowed number of steps NMAX for driver D. Default value of 0 sets no limit for steps. -- Function: int gsl_odeiv2_driver_apply (gsl_odeiv2_driver * D, double * T, const double T1, double Y[]) This function evolves the driver system D from T to T1. Initially vector Y should contain the values of dependent variables at point T. If the function is unable to complete the calculation, an error code from 'gsl_odeiv2_evolve_apply' is returned, and T and Y contain the values from last successful step. If maximum number of steps is reached, a value of 'GSL_EMAXITER' is returned. If the step size drops below minimum value, the function returns with 'GSL_ENOPROG'. If the user-supplied functions defined in the system SYS returns 'GSL_EBADFUNC', the function returns immediately with the same return code. In this case the user must call 'gsl_odeiv2_driver_reset' before calling this function again. -- Function: int gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * D, double * T, const double H, const unsigned long int N, double Y[]) This function evolves the driver system D from T with N steps of size H. If the function is unable to complete the calculation, an error code from 'gsl_odeiv2_evolve_apply_fixed_step' is returned, and T and Y contain the values from last successful step. -- Function: int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * D) This function resets the evolution and stepper objects. -- Function: int gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * D, const double HSTART) The routine resets the evolution and stepper objects and sets new initial step size to HSTART. This function can be used e.g. to change the direction of integration. -- Function: int gsl_odeiv2_driver_free (gsl_odeiv2_driver * D) This function frees the driver object, and the related evolution, stepper and control objects.  File: gsl-ref.info, Node: ODE Example programs, Next: ODE References and Further Reading, Prev: Driver, Up: Ordinary Differential Equations 26.6 Examples ============= The following program solves the second-order nonlinear Van der Pol oscillator equation, u''(t) + \mu u'(t) (u(t)^2 - 1) + u(t) = 0 This can be converted into a first order system suitable for use with the routines described in this chapter by introducing a separate variable for the velocity, v = u'(t), u' = v v' = -u + \mu v (1-u^2) The program begins by defining functions for these derivatives and their Jacobian. The main function uses driver level functions to solve the problem. The program evolves the solution from (u, v) = (1, 0) at t=0 to t=100. The step-size h is automatically adjusted by the controller to maintain an absolute accuracy of 10^{-6} in the function values (u, v). The loop in the example prints the solution at the points t_i = 1, 2, \dots, 100. #include #include #include #include int func (double t, const double y[], double f[], void *params) { double mu = *(double *)params; f[0] = y[1]; f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1); return GSL_SUCCESS; } int jac (double t, const double y[], double *dfdy, double dfdt[], void *params) { double mu = *(double *)params; gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2); gsl_matrix * m = &dfdy_mat.matrix; gsl_matrix_set (m, 0, 0, 0.0); gsl_matrix_set (m, 0, 1, 1.0); gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0); gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0)); dfdt[0] = 0.0; dfdt[1] = 0.0; return GSL_SUCCESS; } int main (void) { double mu = 10; gsl_odeiv2_system sys = {func, jac, 2, &mu}; gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0); int i; double t = 0.0, t1 = 100.0; double y[2] = { 1.0, 0.0 }; for (i = 1; i <= 100; i++) { double ti = i * t1 / 100.0; int status = gsl_odeiv2_driver_apply (d, &t, ti, y); if (status != GSL_SUCCESS) { printf ("error, return value=%d\n", status); break; } printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv2_driver_free (d); return 0; } The user can work with the lower level functions directly, as in the following example. In this case an intermediate result is printed after each successful step instead of equidistant time points. int main (void) { const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk8pd; gsl_odeiv2_step * s = gsl_odeiv2_step_alloc (T, 2); gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (1e-6, 0.0); gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (2); double mu = 10; gsl_odeiv2_system sys = {func, jac, 2, &mu}; double t = 0.0, t1 = 100.0; double h = 1e-6; double y[2] = { 1.0, 0.0 }; while (t < t1) { int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &h, y); if (status != GSL_SUCCESS) break; printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv2_evolve_free (e); gsl_odeiv2_control_free (c); gsl_odeiv2_step_free (s); return 0; } For functions with multiple parameters, the appropriate information can be passed in through the PARAMS argument in 'gsl_odeiv2_system' definition (MU in this example) by using a pointer to a struct. It is also possible to work with a non-adaptive integrator, using only the stepping function itself, 'gsl_odeiv2_driver_apply_fixed_step' or 'gsl_odeiv2_evolve_apply_fixed_step'. The following program uses the driver level function, with fourth-order Runge-Kutta stepping function with a fixed stepsize of 0.001. int main (void) { double mu = 10; gsl_odeiv2_system sys = { func, jac, 2, &mu }; gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk4, 1e-3, 1e-8, 1e-8); double t = 0.0; double y[2] = { 1.0, 0.0 }; int i, s; for (i = 0; i < 100; i++) { s = gsl_odeiv2_driver_apply_fixed_step (d, &t, 1e-3, 1000, y); if (s != GSL_SUCCESS) { printf ("error: driver returned %d\n", s); break; } printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv2_driver_free (d); return s; }  File: gsl-ref.info, Node: ODE References and Further Reading, Prev: ODE Example programs, Up: Ordinary Differential Equations 26.7 References and Further Reading =================================== Ascher, U.M., Petzold, L.R., 'Computer Methods for Ordinary Differential and Differential-Algebraic Equations', SIAM, Philadelphia, 1998. Hairer, E., Norsett, S. P., Wanner, G., 'Solving Ordinary Differential Equations I: Nonstiff Problems', Springer, Berlin, 1993. Hairer, E., Wanner, G., 'Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems', Springer, Berlin, 1996. Many of the basic Runge-Kutta formulas can be found in the Handbook of Mathematical Functions, Abramowitz & Stegun (eds.), 'Handbook of Mathematical Functions', Section 25.5. The implicit Bulirsch-Stoer algorithm 'bsimp' is described in the following paper, G. Bader and P. Deuflhard, "A Semi-Implicit Mid-Point Rule for Stiff Systems of Ordinary Differential Equations.", Numer. Math. 41, 373-398, 1983. The Adams and BDF multistep methods 'msadams' and 'msbdf' are based on the following articles, G. D. Byrne and A. C. Hindmarsh, "A Polyalgorithm for the Numerical Solution of Ordinary Differential Equations.", ACM Trans. Math. Software, 1, 71-96, 1975. P. N. Brown, G. D. Byrne and A. C. Hindmarsh, "VODE: A Variable-coefficient ODE Solver.", SIAM J. Sci. Stat. Comput. 10, 1038-1051, 1989. A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker and C. S. Woodward, "SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers.", ACM Trans. Math. Software 31, 363-396, 2005.  File: gsl-ref.info, Node: Interpolation, Next: Numerical Differentiation, Prev: Ordinary Differential Equations, Up: Top 27 Interpolation **************** This chapter describes functions for performing interpolation. The library provides a variety of interpolation methods, including Cubic splines and Akima splines. The interpolation types are interchangeable, allowing different methods to be used without recompiling. Interpolations can be defined for both normal and periodic boundary conditions. Additional functions are available for computing derivatives and integrals of interpolating functions. These interpolation methods produce curves that pass through each datapoint. To interpolate noisy data with a smoothing curve see *note Basis Splines::. The functions described in this section are declared in the header files 'gsl_interp.h' and 'gsl_spline.h'. * Menu: * Introduction to Interpolation:: * Interpolation Functions:: * Interpolation Types:: * Index Look-up and Acceleration:: * Evaluation of Interpolating Functions:: * Higher-level Interface:: * Interpolation Example programs:: * Interpolation References and Further Reading::  File: gsl-ref.info, Node: Introduction to Interpolation, Next: Interpolation Functions, Up: Interpolation 27.1 Introduction ================= Given a set of data points (x_1, y_1) \dots (x_n, y_n) the routines described in this section compute a continuous interpolating function y(x) such that y(x_i) = y_i. The interpolation is piecewise smooth, and its behavior at the end-points is determined by the type of interpolation used.  File: gsl-ref.info, Node: Interpolation Functions, Next: Interpolation Types, Prev: Introduction to Interpolation, Up: Interpolation 27.2 Interpolation Functions ============================ The interpolation function for a given dataset is stored in a 'gsl_interp' object. These are created by the following functions. -- Function: gsl_interp * gsl_interp_alloc (const gsl_interp_type * T, size_t SIZE) This function returns a pointer to a newly allocated interpolation object of type T for SIZE data-points. -- Function: int gsl_interp_init (gsl_interp * INTERP, const double XA[], const double YA[], size_t SIZE) This function initializes the interpolation object INTERP for the data (XA,YA) where XA and YA are arrays of size SIZE. The interpolation object ('gsl_interp') does not save the data arrays XA and YA and only stores the static state computed from the data. The XA data array is always assumed to be strictly ordered, with increasing x values; the behavior for other arrangements is not defined. -- Function: void gsl_interp_free (gsl_interp * INTERP) This function frees the interpolation object INTERP.  File: gsl-ref.info, Node: Interpolation Types, Next: Index Look-up and Acceleration, Prev: Interpolation Functions, Up: Interpolation 27.3 Interpolation Types ======================== The interpolation library provides six interpolation types: -- Interpolation Type: gsl_interp_linear Linear interpolation. This interpolation method does not require any additional memory. -- Interpolation Type: gsl_interp_polynomial Polynomial interpolation. This method should only be used for interpolating small numbers of points because polynomial interpolation introduces large oscillations, even for well-behaved datasets. The number of terms in the interpolating polynomial is equal to the number of points. -- Interpolation Type: gsl_interp_cspline Cubic spline with natural boundary conditions. The resulting curve is piecewise cubic on each interval, with matching first and second derivatives at the supplied data-points. The second derivative is chosen to be zero at the first point and last point. -- Interpolation Type: gsl_interp_cspline_periodic Cubic spline with periodic boundary conditions. The resulting curve is piecewise cubic on each interval, with matching first and second derivatives at the supplied data-points. The derivatives at the first and last points are also matched. Note that the last point in the data must have the same y-value as the first point, otherwise the resulting periodic interpolation will have a discontinuity at the boundary. -- Interpolation Type: gsl_interp_akima Non-rounded Akima spline with natural boundary conditions. This method uses the non-rounded corner algorithm of Wodicka. -- Interpolation Type: gsl_interp_akima_periodic Non-rounded Akima spline with periodic boundary conditions. This method uses the non-rounded corner algorithm of Wodicka. The following related functions are available: -- Function: const char * gsl_interp_name (const gsl_interp * INTERP) This function returns the name of the interpolation type used by INTERP. For example, printf ("interp uses '%s' interpolation.\n", gsl_interp_name (interp)); would print something like, interp uses 'cspline' interpolation. -- Function: unsigned int gsl_interp_min_size (const gsl_interp * INTERP) -- Function: unsigned int gsl_interp_type_min_size (const gsl_interp_type * T) These functions return the minimum number of points required by the interpolation object INTERP or interpolation type T. For example, Akima spline interpolation requires a minimum of 5 points.  File: gsl-ref.info, Node: Index Look-up and Acceleration, Next: Evaluation of Interpolating Functions, Prev: Interpolation Types, Up: Interpolation 27.4 Index Look-up and Acceleration =================================== The state of searches can be stored in a 'gsl_interp_accel' object, which is a kind of iterator for interpolation lookups. It caches the previous value of an index lookup. When the subsequent interpolation point falls in the same interval its index value can be returned immediately. -- Function: size_t gsl_interp_bsearch (const double X_ARRAY[], double X, size_t INDEX_LO, size_t INDEX_HI) This function returns the index i of the array X_ARRAY such that 'x_array[i] <= x < x_array[i+1]'. The index is searched for in the range [INDEX_LO,INDEX_HI]. An inline version of this function is used when 'HAVE_INLINE' is defined. -- Function: gsl_interp_accel * gsl_interp_accel_alloc (void) This function returns a pointer to an accelerator object, which is a kind of iterator for interpolation lookups. It tracks the state of lookups, thus allowing for application of various acceleration strategies. -- Function: size_t gsl_interp_accel_find (gsl_interp_accel * A, const double X_ARRAY[], size_t SIZE, double X) This function performs a lookup action on the data array X_ARRAY of size SIZE, using the given accelerator A. This is how lookups are performed during evaluation of an interpolation. The function returns an index i such that 'x_array[i] <= x < x_array[i+1]'. An inline version of this function is used when 'HAVE_INLINE' is defined. -- Function: int gsl_interp_accel_reset (gsl_interp_accel * ACC); This function reinitializes the accelerator object ACC. It should be used when the cached information is no longer applicable--for example, when switching to a new dataset. -- Function: void gsl_interp_accel_free (gsl_interp_accel* ACC) This function frees the accelerator object ACC.  File: gsl-ref.info, Node: Evaluation of Interpolating Functions, Next: Higher-level Interface, Prev: Index Look-up and Acceleration, Up: Interpolation 27.5 Evaluation of Interpolating Functions ========================================== -- Function: double gsl_interp_eval (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * ACC) -- Function: int gsl_interp_eval_e (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * ACC, double * Y) These functions return the interpolated value of Y for a given point X, using the interpolation object INTERP, data arrays XA and YA and the accelerator ACC. When X is outside the range of XA, the error code 'GSL_EDOM' is returned with a value of 'GSL_NAN' for Y. -- Function: double gsl_interp_eval_deriv (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * ACC) -- Function: int gsl_interp_eval_deriv_e (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * ACC, double * D) These functions return the derivative D of an interpolated function for a given point X, using the interpolation object INTERP, data arrays XA and YA and the accelerator ACC. -- Function: double gsl_interp_eval_deriv2 (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * ACC) -- Function: int gsl_interp_eval_deriv2_e (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * ACC, double * D2) These functions return the second derivative D2 of an interpolated function for a given point X, using the interpolation object INTERP, data arrays XA and YA and the accelerator ACC. -- Function: double gsl_interp_eval_integ (const gsl_interp * INTERP, const double XA[], const double YA[], double A, double B, gsl_interp_accel * ACC) -- Function: int gsl_interp_eval_integ_e (const gsl_interp * INTERP, const double XA[], const double YA[], double A, double B, gsl_interp_accel * ACC, double * RESULT) These functions return the numerical integral RESULT of an interpolated function over the range [A, B], using the interpolation object INTERP, data arrays XA and YA and the accelerator ACC.  File: gsl-ref.info, Node: Higher-level Interface, Next: Interpolation Example programs, Prev: Evaluation of Interpolating Functions, Up: Interpolation 27.6 Higher-level Interface =========================== The functions described in the previous sections required the user to supply pointers to the x and y arrays on each call. The following functions are equivalent to the corresponding 'gsl_interp' functions but maintain a copy of this data in the 'gsl_spline' object. This removes the need to pass both XA and YA as arguments on each evaluation. These functions are defined in the header file 'gsl_spline.h'. -- Function: gsl_spline * gsl_spline_alloc (const gsl_interp_type * T, size_t SIZE) -- Function: int gsl_spline_init (gsl_spline * SPLINE, const double XA[], const double YA[], size_t SIZE) -- Function: void gsl_spline_free (gsl_spline * SPLINE) -- Function: const char * gsl_spline_name (const gsl_spline * SPLINE) -- Function: unsigned int gsl_spline_min_size (const gsl_spline * SPLINE) -- Function: double gsl_spline_eval (const gsl_spline * SPLINE, double X, gsl_interp_accel * ACC) -- Function: int gsl_spline_eval_e (const gsl_spline * SPLINE, double X, gsl_interp_accel * ACC, double * Y) -- Function: double gsl_spline_eval_deriv (const gsl_spline * SPLINE, double X, gsl_interp_accel * ACC) -- Function: int gsl_spline_eval_deriv_e (const gsl_spline * SPLINE, double X, gsl_interp_accel * ACC, double * D) -- Function: double gsl_spline_eval_deriv2 (const gsl_spline * SPLINE, double X, gsl_interp_accel * ACC) -- Function: int gsl_spline_eval_deriv2_e (const gsl_spline * SPLINE, double X, gsl_interp_accel * ACC, double * D2) -- Function: double gsl_spline_eval_integ (const gsl_spline * SPLINE, double A, double B, gsl_interp_accel * ACC) -- Function: int gsl_spline_eval_integ_e (const gsl_spline * SPLINE, double A, double B, gsl_interp_accel * ACC, double * RESULT)  File: gsl-ref.info, Node: Interpolation Example programs, Next: Interpolation References and Further Reading, Prev: Higher-level Interface, Up: Interpolation 27.7 Examples ============= The following program demonstrates the use of the interpolation and spline functions. It computes a cubic spline interpolation of the 10-point dataset (x_i, y_i) where x_i = i + \sin(i)/2 and y_i = i + \cos(i^2) for i = 0 \dots 9. #include #include #include #include #include int main (void) { int i; double xi, yi, x[10], y[10]; printf ("#m=0,S=2\n"); for (i = 0; i < 10; i++) { x[i] = i + 0.5 * sin (i); y[i] = i + cos (i * i); printf ("%g %g\n", x[i], y[i]); } printf ("#m=1,S=0\n"); { gsl_interp_accel *acc = gsl_interp_accel_alloc (); gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, 10); gsl_spline_init (spline, x, y, 10); for (xi = x[0]; xi < x[9]; xi += 0.01) { yi = gsl_spline_eval (spline, xi, acc); printf ("%g %g\n", xi, yi); } gsl_spline_free (spline); gsl_interp_accel_free (acc); } return 0; } The output is designed to be used with the GNU plotutils 'graph' program, $ ./a.out > interp.dat $ graph -T ps < interp.dat > interp.ps The result shows a smooth interpolation of the original points. The interpolation method can be changed simply by varying the first argument of 'gsl_spline_alloc'. The next program demonstrates a periodic cubic spline with 4 data points. Note that the first and last points must be supplied with the same y-value for a periodic spline. #include #include #include #include #include int main (void) { int N = 4; double x[4] = {0.00, 0.10, 0.27, 0.30}; double y[4] = {0.15, 0.70, -0.10, 0.15}; /* Note: y[0] == y[3] for periodic data */ gsl_interp_accel *acc = gsl_interp_accel_alloc (); const gsl_interp_type *t = gsl_interp_cspline_periodic; gsl_spline *spline = gsl_spline_alloc (t, N); int i; double xi, yi; printf ("#m=0,S=5\n"); for (i = 0; i < N; i++) { printf ("%g %g\n", x[i], y[i]); } printf ("#m=1,S=0\n"); gsl_spline_init (spline, x, y, N); for (i = 0; i <= 100; i++) { xi = (1 - i / 100.0) * x[0] + (i / 100.0) * x[N-1]; yi = gsl_spline_eval (spline, xi, acc); printf ("%g %g\n", xi, yi); } gsl_spline_free (spline); gsl_interp_accel_free (acc); return 0; } The output can be plotted with GNU 'graph'. $ ./a.out > interp.dat $ graph -T ps < interp.dat > interp.ps The result shows a periodic interpolation of the original points. The slope of the fitted curve is the same at the beginning and end of the data, and the second derivative is also.  File: gsl-ref.info, Node: Interpolation References and Further Reading, Prev: Interpolation Example programs, Up: Interpolation 27.8 References and Further Reading =================================== Descriptions of the interpolation algorithms and further references can be found in the following books: C.W. Ueberhuber, 'Numerical Computation (Volume 1), Chapter 9 "Interpolation"', Springer (1997), ISBN 3-540-62058-3. D.M. Young, R.T. Gregory 'A Survey of Numerical Mathematics (Volume 1), Chapter 6.8', Dover (1988), ISBN 0-486-65691-8.  File: gsl-ref.info, Node: Numerical Differentiation, Next: Chebyshev Approximations, Prev: Interpolation, Up: Top 28 Numerical Differentiation **************************** The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative. These functions are declared in the header file 'gsl_deriv.h'. * Menu: * Numerical Differentiation functions:: * Numerical Differentiation Examples:: * Numerical Differentiation References::  File: gsl-ref.info, Node: Numerical Differentiation functions, Next: Numerical Differentiation Examples, Up: Numerical Differentiation 28.1 Functions ============== -- Function: int gsl_deriv_central (const gsl_function * F, double X, double H, double * RESULT, double * ABSERR) This function computes the numerical derivative of the function F at the point X using an adaptive central difference algorithm with a step-size of H. The derivative is returned in RESULT and an estimate of its absolute error is returned in ABSERR. The initial value of H is used to estimate an optimal step-size, based on the scaling of the truncation error and round-off error in the derivative calculation. The derivative is computed using a 5-point rule for equally spaced abscissae at x-h, x-h/2, x, x+h/2, x+h, with an error estimate taken from the difference between the 5-point rule and the corresponding 3-point rule x-h, x, x+h. Note that the value of the function at x does not contribute to the derivative calculation, so only 4-points are actually used. -- Function: int gsl_deriv_forward (const gsl_function * F, double X, double H, double * RESULT, double * ABSERR) This function computes the numerical derivative of the function F at the point X using an adaptive forward difference algorithm with a step-size of H. The function is evaluated only at points greater than X, and never at X itself. The derivative is returned in RESULT and an estimate of its absolute error is returned in ABSERR. This function should be used if f(x) has a discontinuity at X, or is undefined for values less than X. The initial value of H is used to estimate an optimal step-size, based on the scaling of the truncation error and round-off error in the derivative calculation. The derivative at x is computed using an "open" 4-point rule for equally spaced abscissae at x+h/4, x+h/2, x+3h/4, x+h, with an error estimate taken from the difference between the 4-point rule and the corresponding 2-point rule x+h/2, x+h. -- Function: int gsl_deriv_backward (const gsl_function * F, double X, double H, double * RESULT, double * ABSERR) This function computes the numerical derivative of the function F at the point X using an adaptive backward difference algorithm with a step-size of H. The function is evaluated only at points less than X, and never at X itself. The derivative is returned in RESULT and an estimate of its absolute error is returned in ABSERR. This function should be used if f(x) has a discontinuity at X, or is undefined for values greater than X. This function is equivalent to calling 'gsl_deriv_forward' with a negative step-size.  File: gsl-ref.info, Node: Numerical Differentiation Examples, Next: Numerical Differentiation References, Prev: Numerical Differentiation functions, Up: Numerical Differentiation 28.2 Examples ============= The following code estimates the derivative of the function f(x) = x^{3/2} at x=2 and at x=0. The function f(x) is undefined for x<0 so the derivative at x=0 is computed using 'gsl_deriv_forward'. #include #include #include double f (double x, void * params) { return pow (x, 1.5); } int main (void) { gsl_function F; double result, abserr; F.function = &f; F.params = 0; printf ("f(x) = x^(3/2)\n"); gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr); printf ("x = 2.0\n"); printf ("f'(x) = %.10f +/- %.10f\n", result, abserr); printf ("exact = %.10f\n\n", 1.5 * sqrt(2.0)); gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr); printf ("x = 0.0\n"); printf ("f'(x) = %.10f +/- %.10f\n", result, abserr); printf ("exact = %.10f\n", 0.0); return 0; } Here is the output of the program, $ ./a.out f(x) = x^(3/2) x = 2.0 f'(x) = 2.1213203120 +/- 0.0000004064 exact = 2.1213203436 x = 0.0 f'(x) = 0.0000000160 +/- 0.0000000339 exact = 0.0000000000  File: gsl-ref.info, Node: Numerical Differentiation References, Prev: Numerical Differentiation Examples, Up: Numerical Differentiation 28.3 References and Further Reading =================================== The algorithms used by these functions are described in the following sources: Abramowitz and Stegun, 'Handbook of Mathematical Functions', Section 25.3.4, and Table 25.5 (Coefficients for Differentiation). S.D. Conte and Carl de Boor, 'Elementary Numerical Analysis: An Algorithmic Approach', McGraw-Hill, 1972.  File: gsl-ref.info, Node: Chebyshev Approximations, Next: Series Acceleration, Prev: Numerical Differentiation, Up: Top 29 Chebyshev Approximations *************************** This chapter describes routines for computing Chebyshev approximations to univariate functions. A Chebyshev approximation is a truncation of the series f(x) = \sum c_n T_n(x), where the Chebyshev polynomials T_n(x) = \cos(n \arccos x) provide an orthogonal basis of polynomials on the interval [-1,1] with the weight function 1 / \sqrt{1-x^2}. The first few Chebyshev polynomials are, T_0(x) = 1, T_1(x) = x, T_2(x) = 2 x^2 - 1. For further information see Abramowitz & Stegun, Chapter 22. The functions described in this chapter are declared in the header file 'gsl_chebyshev.h'. * Menu: * Chebyshev Definitions:: * Creation and Calculation of Chebyshev Series:: * Auxiliary Functions for Chebyshev Series:: * Chebyshev Series Evaluation:: * Derivatives and Integrals:: * Chebyshev Approximation Examples:: * Chebyshev Approximation References and Further Reading::  File: gsl-ref.info, Node: Chebyshev Definitions, Next: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations 29.1 Definitions ================ A Chebyshev series is stored using the following structure, typedef struct { double * c; /* coefficients c[0] .. c[order] */ int order; /* order of expansion */ double a; /* lower interval point */ double b; /* upper interval point */ ... } gsl_cheb_series The approximation is made over the range [a,b] using ORDER+1 terms, including the coefficient c[0]. The series is computed using the following convention, f(x) = (c_0 / 2) + \sum_{n=1} c_n T_n(x) which is needed when accessing the coefficients directly.  File: gsl-ref.info, Node: Creation and Calculation of Chebyshev Series, Next: Auxiliary Functions for Chebyshev Series, Prev: Chebyshev Definitions, Up: Chebyshev Approximations 29.2 Creation and Calculation of Chebyshev Series ================================================= -- Function: gsl_cheb_series * gsl_cheb_alloc (const size_t N) This function allocates space for a Chebyshev series of order N and returns a pointer to a new 'gsl_cheb_series' struct. -- Function: void gsl_cheb_free (gsl_cheb_series * CS) This function frees a previously allocated Chebyshev series CS. -- Function: int gsl_cheb_init (gsl_cheb_series * CS, const gsl_function * F, const double A, const double B) This function computes the Chebyshev approximation CS for the function F over the range (a,b) to the previously specified order. The computation of the Chebyshev approximation is an O(n^2) process, and requires n function evaluations.  File: gsl-ref.info, Node: Auxiliary Functions for Chebyshev Series, Next: Chebyshev Series Evaluation, Prev: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations 29.3 Auxiliary Functions ======================== The following functions provide information about an existing Chebyshev series. -- Function: size_t gsl_cheb_order (const gsl_cheb_series * CS) This function returns the order of Chebyshev series CS. -- Function: size_t gsl_cheb_size (const gsl_cheb_series * CS) -- Function: double * gsl_cheb_coeffs (const gsl_cheb_series * CS) These functions return the size of the Chebyshev coefficient array 'c[]' and a pointer to its location in memory for the Chebyshev series CS.  File: gsl-ref.info, Node: Chebyshev Series Evaluation, Next: Derivatives and Integrals, Prev: Auxiliary Functions for Chebyshev Series, Up: Chebyshev Approximations 29.4 Chebyshev Series Evaluation ================================ -- Function: double gsl_cheb_eval (const gsl_cheb_series * CS, double X) This function evaluates the Chebyshev series CS at a given point X. -- Function: int gsl_cheb_eval_err (const gsl_cheb_series * CS, const double X, double * RESULT, double * ABSERR) This function computes the Chebyshev series CS at a given point X, estimating both the series RESULT and its absolute error ABSERR. The error estimate is made from the first neglected term in the series. -- Function: double gsl_cheb_eval_n (const gsl_cheb_series * CS, size_t ORDER, double X) This function evaluates the Chebyshev series CS at a given point X, to (at most) the given order ORDER. -- Function: int gsl_cheb_eval_n_err (const gsl_cheb_series * CS, const size_t ORDER, const double X, double * RESULT, double * ABSERR) This function evaluates a Chebyshev series CS at a given point X, estimating both the series RESULT and its absolute error ABSERR, to (at most) the given order ORDER. The error estimate is made from the first neglected term in the series.  File: gsl-ref.info, Node: Derivatives and Integrals, Next: Chebyshev Approximation Examples, Prev: Chebyshev Series Evaluation, Up: Chebyshev Approximations 29.5 Derivatives and Integrals ============================== The following functions allow a Chebyshev series to be differentiated or integrated, producing a new Chebyshev series. Note that the error estimate produced by evaluating the derivative series will be underestimated due to the contribution of higher order terms being neglected. -- Function: int gsl_cheb_calc_deriv (gsl_cheb_series * DERIV, const gsl_cheb_series * CS) This function computes the derivative of the series CS, storing the derivative coefficients in the previously allocated DERIV. The two series CS and DERIV must have been allocated with the same order. -- Function: int gsl_cheb_calc_integ (gsl_cheb_series * INTEG, const gsl_cheb_series * CS) This function computes the integral of the series CS, storing the integral coefficients in the previously allocated INTEG. The two series CS and INTEG must have been allocated with the same order. The lower limit of the integration is taken to be the left hand end of the range A.  File: gsl-ref.info, Node: Chebyshev Approximation Examples, Next: Chebyshev Approximation References and Further Reading, Prev: Derivatives and Integrals, Up: Chebyshev Approximations 29.6 Examples ============= The following example program computes Chebyshev approximations to a step function. This is an extremely difficult approximation to make, due to the discontinuity, and was chosen as an example where approximation error is visible. For smooth functions the Chebyshev approximation converges extremely rapidly and errors would not be visible. #include #include #include double f (double x, void *p) { if (x < 0.5) return 0.25; else return 0.75; } int main (void) { int i, n = 10000; gsl_cheb_series *cs = gsl_cheb_alloc (40); gsl_function F; F.function = f; F.params = 0; gsl_cheb_init (cs, &F, 0.0, 1.0); for (i = 0; i < n; i++) { double x = i / (double)n; double r10 = gsl_cheb_eval_n (cs, 10, x); double r40 = gsl_cheb_eval (cs, x); printf ("%g %g %g %g\n", x, GSL_FN_EVAL (&F, x), r10, r40); } gsl_cheb_free (cs); return 0; } The output from the program gives the original function, 10-th order approximation and 40-th order approximation, all sampled at intervals of 0.001 in x.  File: gsl-ref.info, Node: Chebyshev Approximation References and Further Reading, Prev: Chebyshev Approximation Examples, Up: Chebyshev Approximations 29.7 References and Further Reading =================================== The following paper describes the use of Chebyshev series, R. Broucke, "Ten Subroutines for the Manipulation of Chebyshev Series [C1] (Algorithm 446)". 'Communications of the ACM' 16(4), 254-256 (1973)  File: gsl-ref.info, Node: Series Acceleration, Next: Wavelet Transforms, Prev: Chebyshev Approximations, Up: Top 30 Series Acceleration ********************** The functions described in this chapter accelerate the convergence of a series using the Levin u-transform. This method takes a small number of terms from the start of a series and uses a systematic approximation to compute an extrapolated value and an estimate of its error. The u-transform works for both convergent and divergent series, including asymptotic series. These functions are declared in the header file 'gsl_sum.h'. * Menu: * Acceleration functions:: * Acceleration functions without error estimation:: * Example of accelerating a series:: * Series Acceleration References::  File: gsl-ref.info, Node: Acceleration functions, Next: Acceleration functions without error estimation, Up: Series Acceleration 30.1 Acceleration functions =========================== The following functions compute the full Levin u-transform of a series with its error estimate. The error estimate is computed by propagating rounding errors from each term through to the final extrapolation. These functions are intended for summing analytic series where each term is known to high accuracy, and the rounding errors are assumed to originate from finite precision. They are taken to be relative errors of order 'GSL_DBL_EPSILON' for each term. The calculation of the error in the extrapolated value is an O(N^2) process, which is expensive in time and memory. A faster but less reliable method which estimates the error from the convergence of the extrapolated value is described in the next section. For the method described here a full table of intermediate values and derivatives through to O(N) must be computed and stored, but this does give a reliable error estimate. -- Function: gsl_sum_levin_u_workspace * gsl_sum_levin_u_alloc (size_t N) This function allocates a workspace for a Levin u-transform of N terms. The size of the workspace is O(2n^2 + 3n). -- Function: void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * W) This function frees the memory associated with the workspace W. -- Function: int gsl_sum_levin_u_accel (const double * ARRAY, size_t ARRAY_SIZE, gsl_sum_levin_u_workspace * W, double * SUM_ACCEL, double * ABSERR) This function takes the terms of a series in ARRAY of size ARRAY_SIZE and computes the extrapolated limit of the series using a Levin u-transform. Additional working space must be provided in W. The extrapolated sum is stored in SUM_ACCEL, with an estimate of the absolute error stored in ABSERR. The actual term-by-term sum is returned in 'w->sum_plain'. The algorithm calculates the truncation error (the difference between two successive extrapolations) and round-off error (propagated from the individual terms) to choose an optimal number of terms for the extrapolation. All the terms of the series passed in through ARRAY should be non-zero.  File: gsl-ref.info, Node: Acceleration functions without error estimation, Next: Example of accelerating a series, Prev: Acceleration functions, Up: Series Acceleration 30.2 Acceleration functions without error estimation ==================================================== The functions described in this section compute the Levin u-transform of series and attempt to estimate the error from the "truncation error" in the extrapolation, the difference between the final two approximations. Using this method avoids the need to compute an intermediate table of derivatives because the error is estimated from the behavior of the extrapolated value itself. Consequently this algorithm is an O(N) process and only requires O(N) terms of storage. If the series converges sufficiently fast then this procedure can be acceptable. It is appropriate to use this method when there is a need to compute many extrapolations of series with similar convergence properties at high-speed. For example, when numerically integrating a function defined by a parameterized series where the parameter varies only slightly. A reliable error estimate should be computed first using the full algorithm described above in order to verify the consistency of the results. -- Function: gsl_sum_levin_utrunc_workspace * gsl_sum_levin_utrunc_alloc (size_t N) This function allocates a workspace for a Levin u-transform of N terms, without error estimation. The size of the workspace is O(3n). -- Function: void gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * W) This function frees the memory associated with the workspace W. -- Function: int gsl_sum_levin_utrunc_accel (const double * ARRAY, size_t ARRAY_SIZE, gsl_sum_levin_utrunc_workspace * W, double * SUM_ACCEL, double * ABSERR_TRUNC) This function takes the terms of a series in ARRAY of size ARRAY_SIZE and computes the extrapolated limit of the series using a Levin u-transform. Additional working space must be provided in W. The extrapolated sum is stored in SUM_ACCEL. The actual term-by-term sum is returned in 'w->sum_plain'. The algorithm terminates when the difference between two successive extrapolations reaches a minimum or is sufficiently small. The difference between these two values is used as estimate of the error and is stored in ABSERR_TRUNC. To improve the reliability of the algorithm the extrapolated values are replaced by moving averages when calculating the truncation error, smoothing out any fluctuations.  File: gsl-ref.info, Node: Example of accelerating a series, Next: Series Acceleration References, Prev: Acceleration functions without error estimation, Up: Series Acceleration 30.3 Examples ============= The following code calculates an estimate of \zeta(2) = \pi^2 / 6 using the series, \zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + ... After N terms the error in the sum is O(1/N), making direct summation of the series converge slowly. #include #include #include #define N 20 int main (void) { double t[N]; double sum_accel, err; double sum = 0; int n; gsl_sum_levin_u_workspace * w = gsl_sum_levin_u_alloc (N); const double zeta_2 = M_PI * M_PI / 6.0; /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */ for (n = 0; n < N; n++) { double np1 = n + 1.0; t[n] = 1.0 / (np1 * np1); sum += t[n]; } gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err); printf ("term-by-term sum = % .16f using %d terms\n", sum, N); printf ("term-by-term sum = % .16f using %d terms\n", w->sum_plain, w->terms_used); printf ("exact value = % .16f\n", zeta_2); printf ("accelerated sum = % .16f using %d terms\n", sum_accel, w->terms_used); printf ("estimated error = % .16f\n", err); printf ("actual error = % .16f\n", sum_accel - zeta_2); gsl_sum_levin_u_free (w); return 0; } The output below shows that the Levin u-transform is able to obtain an estimate of the sum to 1 part in 10^10 using the first eleven terms of the series. The error estimate returned by the function is also accurate, giving the correct number of significant digits. $ ./a.out term-by-term sum = 1.5961632439130233 using 20 terms term-by-term sum = 1.5759958390005426 using 13 terms exact value = 1.6449340668482264 accelerated sum = 1.6449340668166479 using 13 terms estimated error = 0.0000000000508580 actual error = -0.0000000000315785 Note that a direct summation of this series would require 10^10 terms to achieve the same precision as the accelerated sum does in 13 terms.  File: gsl-ref.info, Node: Series Acceleration References, Prev: Example of accelerating a series, Up: Series Acceleration 30.4 References and Further Reading =================================== The algorithms used by these functions are described in the following papers, T. Fessler, W.F. Ford, D.A. Smith, HURRY: An acceleration algorithm for scalar sequences and series 'ACM Transactions on Mathematical Software', 9(3):346-354, 1983. and Algorithm 602 9(3):355-357, 1983. The theory of the u-transform was presented by Levin, D. Levin, Development of Non-Linear Transformations for Improving Convergence of Sequences, 'Intern. J. Computer Math.' B3:371-388, 1973. A review paper on the Levin Transform is available online, Herbert H. H. Homeier, Scalar Levin-Type Sequence Transformations, .  File: gsl-ref.info, Node: Wavelet Transforms, Next: Discrete Hankel Transforms, Prev: Series Acceleration, Up: Top 31 Wavelet Transforms ********************* This chapter describes functions for performing Discrete Wavelet Transforms (DWTs). The library includes wavelets for real data in both one and two dimensions. The wavelet functions are declared in the header files 'gsl_wavelet.h' and 'gsl_wavelet2d.h'. * Menu: * DWT Definitions:: * DWT Initialization:: * DWT Transform Functions:: * DWT Examples:: * DWT References::  File: gsl-ref.info, Node: DWT Definitions, Next: DWT Initialization, Up: Wavelet Transforms 31.1 Definitions ================ The continuous wavelet transform and its inverse are defined by the relations, w(s,\tau) = \int f(t) * \psi^*_{s,\tau}(t) dt and, f(t) = \int \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau ds where the basis functions \psi_{s,\tau} are obtained by scaling and translation from a single function, referred to as the "mother wavelet". The discrete version of the wavelet transform acts on equally-spaced samples, with fixed scaling and translation steps (s, \tau). The frequency and time axes are sampled "dyadically" on scales of 2^j through a level parameter j. The resulting family of functions {\psi_{j,n}} constitutes an orthonormal basis for square-integrable signals. The discrete wavelet transform is an O(N) algorithm, and is also referred to as the "fast wavelet transform".  File: gsl-ref.info, Node: DWT Initialization, Next: DWT Transform Functions, Prev: DWT Definitions, Up: Wavelet Transforms 31.2 Initialization =================== The 'gsl_wavelet' structure contains the filter coefficients defining the wavelet and any associated offset parameters. -- Function: gsl_wavelet * gsl_wavelet_alloc (const gsl_wavelet_type * T, size_t K) This function allocates and initializes a wavelet object of type T. The parameter K selects the specific member of the wavelet family. A null pointer is returned if insufficient memory is available or if a unsupported member is selected. The following wavelet types are implemented: -- Wavelet: gsl_wavelet_daubechies -- Wavelet: gsl_wavelet_daubechies_centered This is the Daubechies wavelet family of maximum phase with k/2 vanishing moments. The implemented wavelets are k=4, 6, ..., 20, with K even. -- Wavelet: gsl_wavelet_haar -- Wavelet: gsl_wavelet_haar_centered This is the Haar wavelet. The only valid choice of k for the Haar wavelet is k=2. -- Wavelet: gsl_wavelet_bspline -- Wavelet: gsl_wavelet_bspline_centered This is the biorthogonal B-spline wavelet family of order (i,j). The implemented values of k = 100*i + j are 103, 105, 202, 204, 206, 208, 301, 303, 305 307, 309. The centered forms of the wavelets align the coefficients of the various sub-bands on edges. Thus the resulting visualization of the coefficients of the wavelet transform in the phase plane is easier to understand. -- Function: const char * gsl_wavelet_name (const gsl_wavelet * W) This function returns a pointer to the name of the wavelet family for W. -- Function: void gsl_wavelet_free (gsl_wavelet * W) This function frees the wavelet object W. The 'gsl_wavelet_workspace' structure contains scratch space of the same size as the input data and is used to hold intermediate results during the transform. -- Function: gsl_wavelet_workspace * gsl_wavelet_workspace_alloc (size_t N) This function allocates a workspace for the discrete wavelet transform. To perform a one-dimensional transform on N elements, a workspace of size N must be provided. For two-dimensional transforms of N-by-N matrices it is sufficient to allocate a workspace of size N, since the transform operates on individual rows and columns. A null pointer is returned if insufficient memory is available. -- Function: void gsl_wavelet_workspace_free (gsl_wavelet_workspace * WORK) This function frees the allocated workspace WORK.  File: gsl-ref.info, Node: DWT Transform Functions, Next: DWT Examples, Prev: DWT Initialization, Up: Wavelet Transforms 31.3 Transform Functions ======================== This sections describes the actual functions performing the discrete wavelet transform. Note that the transforms use periodic boundary conditions. If the signal is not periodic in the sample length then spurious coefficients will appear at the beginning and end of each level of the transform. * Menu: * DWT in one dimension:: * DWT in two dimension::  File: gsl-ref.info, Node: DWT in one dimension, Next: DWT in two dimension, Up: DWT Transform Functions 31.3.1 Wavelet transforms in one dimension ------------------------------------------ -- Function: int gsl_wavelet_transform (const gsl_wavelet * W, double * DATA, size_t STRIDE, size_t N, gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet_transform_forward (const gsl_wavelet * W, double * DATA, size_t STRIDE, size_t N, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet_transform_inverse (const gsl_wavelet * W, double * DATA, size_t STRIDE, size_t N, gsl_wavelet_workspace * WORK) These functions compute in-place forward and inverse discrete wavelet transforms of length N with stride STRIDE on the array DATA. The length of the transform N is restricted to powers of two. For the 'transform' version of the function the argument DIR can be either 'forward' (+1) or 'backward' (-1). A workspace WORK of length N must be provided. For the forward transform, the elements of the original array are replaced by the discrete wavelet transform f_i -> w_{j,k} in a packed triangular storage layout, where J is the index of the level j = 0 ... J-1 and K is the index of the coefficient within each level, k = 0 ... (2^j)-1. The total number of levels is J = \log_2(n). The output data has the following form, (s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0}, ..., d_{j,k}, ..., d_{J-1,2^{J-1}-1}) where the first element is the smoothing coefficient s_{-1,0}, followed by the detail coefficients d_{j,k} for each level j. The backward transform inverts these coefficients to obtain the original data. These functions return a status of 'GSL_SUCCESS' upon successful completion. 'GSL_EINVAL' is returned if N is not an integer power of 2 or if insufficient workspace is provided.  File: gsl-ref.info, Node: DWT in two dimension, Prev: DWT in one dimension, Up: DWT Transform Functions 31.3.2 Wavelet transforms in two dimension ------------------------------------------ The library provides functions to perform two-dimensional discrete wavelet transforms on square matrices. The matrix dimensions must be an integer power of two. There are two possible orderings of the rows and columns in the two-dimensional wavelet transform, referred to as the "standard" and "non-standard" forms. The "standard" transform performs a complete discrete wavelet transform on the rows of the matrix, followed by a separate complete discrete wavelet transform on the columns of the resulting row-transformed matrix. This procedure uses the same ordering as a two-dimensional Fourier transform. The "non-standard" transform is performed in interleaved passes on the rows and columns of the matrix for each level of the transform. The first level of the transform is applied to the matrix rows, and then to the matrix columns. This procedure is then repeated across the rows and columns of the data for the subsequent levels of the transform, until the full discrete wavelet transform is complete. The non-standard form of the discrete wavelet transform is typically used in image analysis. The functions described in this section are declared in the header file 'gsl_wavelet2d.h'. -- Function: int gsl_wavelet2d_transform (const gsl_wavelet * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2, gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_transform_forward (const gsl_wavelet * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_transform_inverse (const gsl_wavelet * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2, gsl_wavelet_workspace * WORK) These functions compute two-dimensional in-place forward and inverse discrete wavelet transforms in standard form on the array DATA stored in row-major form with dimensions SIZE1 and SIZE2 and physical row length TDA. The dimensions must be equal (square matrix) and are restricted to powers of two. For the 'transform' version of the function the argument DIR can be either 'forward' (+1) or 'backward' (-1). A workspace WORK of the appropriate size must be provided. On exit, the appropriate elements of the array DATA are replaced by their two-dimensional wavelet transform. The functions return a status of 'GSL_SUCCESS' upon successful completion. 'GSL_EINVAL' is returned if SIZE1 and SIZE2 are not equal and integer powers of 2, or if insufficient workspace is provided. -- Function: int gsl_wavelet2d_transform_matrix (const gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK) These functions compute the two-dimensional in-place wavelet transform on a matrix A. -- Function: int gsl_wavelet2d_nstransform (const gsl_wavelet * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2, gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_nstransform_forward (const gsl_wavelet * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_nstransform_inverse (const gsl_wavelet * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2, gsl_wavelet_workspace * WORK) These functions compute the two-dimensional wavelet transform in non-standard form. -- Function: int gsl_wavelet2d_nstransform_matrix (const gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK) -- Function: int gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK) These functions compute the non-standard form of the two-dimensional in-place wavelet transform on a matrix A.  File: gsl-ref.info, Node: DWT Examples, Next: DWT References, Prev: DWT Transform Functions, Up: Wavelet Transforms 31.4 Examples ============= The following program demonstrates the use of the one-dimensional wavelet transform functions. It computes an approximation to an input signal (of length 256) using the 20 largest components of the wavelet transform, while setting the others to zero. #include #include #include #include int main (int argc, char **argv) { int i, n = 256, nc = 20; double *data = malloc (n * sizeof (double)); double *abscoeff = malloc (n * sizeof (double)); size_t *p = malloc (n * sizeof (size_t)); FILE * f; gsl_wavelet *w; gsl_wavelet_workspace *work; w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4); work = gsl_wavelet_workspace_alloc (n); f = fopen (argv[1], "r"); for (i = 0; i < n; i++) { fscanf (f, "%lg", &data[i]); } fclose (f); gsl_wavelet_transform_forward (w, data, 1, n, work); for (i = 0; i < n; i++) { abscoeff[i] = fabs (data[i]); } gsl_sort_index (p, abscoeff, 1, n); for (i = 0; (i + nc) < n; i++) data[p[i]] = 0; gsl_wavelet_transform_inverse (w, data, 1, n, work); for (i = 0; i < n; i++) { printf ("%g\n", data[i]); } gsl_wavelet_free (w); gsl_wavelet_workspace_free (work); free (data); free (abscoeff); free (p); return 0; } The output can be used with the GNU plotutils 'graph' program, $ ./a.out ecg.dat > dwt.dat $ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps  File: gsl-ref.info, Node: DWT References, Prev: DWT Examples, Up: Wavelet Transforms 31.5 References and Further Reading =================================== The mathematical background to wavelet transforms is covered in the original lectures by Daubechies, Ingrid Daubechies. Ten Lectures on Wavelets. 'CBMS-NSF Regional Conference Series in Applied Mathematics' (1992), SIAM, ISBN 0898712742. An easy to read introduction to the subject with an emphasis on the application of the wavelet transform in various branches of science is, Paul S. Addison. 'The Illustrated Wavelet Transform Handbook'. Institute of Physics Publishing (2002), ISBN 0750306920. For extensive coverage of signal analysis by wavelets, wavelet packets and local cosine bases see, S. G. Mallat. 'A wavelet tour of signal processing' (Second edition). Academic Press (1999), ISBN 012466606X. The concept of multiresolution analysis underlying the wavelet transform is described in, S. G. Mallat. Multiresolution Approximations and Wavelet Orthonormal Bases of L^2(R). 'Transactions of the American Mathematical Society', 315(1), 1989, 69-87. S. G. Mallat. A Theory for Multiresolution Signal Decomposition--The Wavelet Representation. 'IEEE Transactions on Pattern Analysis and Machine Intelligence', 11, 1989, 674-693. The coefficients for the individual wavelet families implemented by the library can be found in the following papers, I. Daubechies. Orthonormal Bases of Compactly Supported Wavelets. 'Communications on Pure and Applied Mathematics', 41 (1988) 909-996. A. Cohen, I. Daubechies, and J.-C. Feauveau. Biorthogonal Bases of Compactly Supported Wavelets. 'Communications on Pure and Applied Mathematics', 45 (1992) 485-560. The PhysioNet archive of physiological datasets can be found online at and is described in the following paper, Goldberger et al. PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals. 'Circulation' 101(23):e215-e220 2000.  File: gsl-ref.info, Node: Discrete Hankel Transforms, Next: One dimensional Root-Finding, Prev: Wavelet Transforms, Up: Top 32 Discrete Hankel Transforms ***************************** This chapter describes functions for performing Discrete Hankel Transforms (DHTs). The functions are declared in the header file 'gsl_dht.h'. * Menu: * Discrete Hankel Transform Definition:: * Discrete Hankel Transform Functions:: * Discrete Hankel Transform References::  File: gsl-ref.info, Node: Discrete Hankel Transform Definition, Next: Discrete Hankel Transform Functions, Up: Discrete Hankel Transforms 32.1 Definitions ================ The discrete Hankel transform acts on a vector of sampled data, where the samples are assumed to have been taken at points related to the zeroes of a Bessel function of fixed order; compare this to the case of the discrete Fourier transform, where samples are taken at points related to the zeroes of the sine or cosine function. Specifically, let f(t) be a function on the unit interval and j_(\nu,m) the m-th zero of the Bessel function J_\nu(x). Then the finite \nu-Hankel transform of f(t) is defined to be the set of numbers g_m given by, g_m = \int_0^1 t dt J_\nu(j_(\nu,m)t) f(t), so that, f(t) = \sum_{m=1}^\infty (2 J_\nu(j_(\nu,m)t) / J_(\nu+1)(j_(\nu,m))^2) g_m. Suppose that f is band-limited in the sense that g_m=0 for m > M. Then we have the following fundamental sampling theorem. g_m = (2 / j_(\nu,M)^2) \sum_{k=1}^{M-1} f(j_(\nu,k)/j_(\nu,M)) (J_\nu(j_(\nu,m) j_(\nu,k) / j_(\nu,M)) / J_(\nu+1)(j_(\nu,k))^2). It is this discrete expression which defines the discrete Hankel transform. The kernel in the summation above defines the matrix of the \nu-Hankel transform of size M-1. The coefficients of this matrix, being dependent on \nu and M, must be precomputed and stored; the 'gsl_dht' object encapsulates this data. The allocation function 'gsl_dht_alloc' returns a 'gsl_dht' object which must be properly initialized with 'gsl_dht_init' before it can be used to perform transforms on data sample vectors, for fixed \nu and M, using the 'gsl_dht_apply' function. The implementation allows a scaling of the fundamental interval, for convenience, so that one can assume the function is defined on the interval [0,X], rather than the unit interval. Notice that by assumption f(t) vanishes at the endpoints of the interval, consistent with the inversion formula and the sampling formula given above. Therefore, this transform corresponds to an orthogonal expansion in eigenfunctions of the Dirichlet problem for the Bessel differential equation.  File: gsl-ref.info, Node: Discrete Hankel Transform Functions, Next: Discrete Hankel Transform References, Prev: Discrete Hankel Transform Definition, Up: Discrete Hankel Transforms 32.2 Functions ============== -- Function: gsl_dht * gsl_dht_alloc (size_t SIZE) This function allocates a Discrete Hankel transform object of size SIZE. -- Function: int gsl_dht_init (gsl_dht * T, double NU, double XMAX) This function initializes the transform T for the given values of NU and XMAX. -- Function: gsl_dht * gsl_dht_new (size_t SIZE, double NU, double XMAX) This function allocates a Discrete Hankel transform object of size SIZE and initializes it for the given values of NU and XMAX. -- Function: void gsl_dht_free (gsl_dht * T) This function frees the transform T. -- Function: int gsl_dht_apply (const gsl_dht * T, double * F_IN, double * F_OUT) This function applies the transform T to the array F_IN whose size is equal to the size of the transform. The result is stored in the array F_OUT which must be of the same length. Applying this function to its output gives the original data multiplied by (1/j_(\nu,M))^2, up to numerical errors. -- Function: double gsl_dht_x_sample (const gsl_dht * T, int N) This function returns the value of the N-th sample point in the unit interval, (j_{\nu,n+1}/j_{\nu,M}) X. These are the points where the function f(t) is assumed to be sampled. -- Function: double gsl_dht_k_sample (const gsl_dht * T, int N) This function returns the value of the N-th sample point in "k-space", j_{\nu,n+1}/X.  File: gsl-ref.info, Node: Discrete Hankel Transform References, Prev: Discrete Hankel Transform Functions, Up: Discrete Hankel Transforms 32.3 References and Further Reading =================================== The algorithms used by these functions are described in the following papers, H. Fisk Johnson, Comp. Phys. Comm. 43, 181 (1987). D. Lemoine, J. Chem. Phys. 101, 3936 (1994).  File: gsl-ref.info, Node: One dimensional Root-Finding, Next: One dimensional Minimization, Prev: Discrete Hankel Transforms, Up: Top 33 One dimensional Root-Finding ******************************* This chapter describes routines for finding roots of arbitrary one-dimensional functions. The library provides low level components for a variety of iterative solvers and convergence tests. These can be combined by the user to achieve the desired solution, with full access to the intermediate steps of the iteration. Each class of methods uses the same framework, so that you can switch between solvers at runtime without needing to recompile your program. Each instance of a solver keeps track of its own state, allowing the solvers to be used in multi-threaded programs. The header file 'gsl_roots.h' contains prototypes for the root finding functions and related declarations. * Menu: * Root Finding Overview:: * Root Finding Caveats:: * Initializing the Solver:: * Providing the function to solve:: * Search Bounds and Guesses:: * Root Finding Iteration:: * Search Stopping Parameters:: * Root Bracketing Algorithms:: * Root Finding Algorithms using Derivatives:: * Root Finding Examples:: * Root Finding References and Further Reading::  File: gsl-ref.info, Node: Root Finding Overview, Next: Root Finding Caveats, Up: One dimensional Root-Finding 33.1 Overview ============= One-dimensional root finding algorithms can be divided into two classes, "root bracketing" and "root polishing". Algorithms which proceed by bracketing a root are guaranteed to converge. Bracketing algorithms begin with a bounded region known to contain a root. The size of this bounded region is reduced, iteratively, until it encloses the root to a desired tolerance. This provides a rigorous error estimate for the location of the root. The technique of "root polishing" attempts to improve an initial guess to the root. These algorithms converge only if started "close enough" to a root, and sacrifice a rigorous error bound for speed. By approximating the behavior of a function in the vicinity of a root they attempt to find a higher order improvement of an initial guess. When the behavior of the function is compatible with the algorithm and a good initial guess is available a polishing algorithm can provide rapid convergence. In GSL both types of algorithm are available in similar frameworks. The user provides a high-level driver for the algorithms, and the library provides the individual functions necessary for each of the steps. There are three main phases of the iteration. The steps are, * initialize solver state, S, for algorithm T * update S using the iteration T * test S for convergence, and repeat iteration if necessary The state for bracketing solvers is held in a 'gsl_root_fsolver' struct. The updating procedure uses only function evaluations (not derivatives). The state for root polishing solvers is held in a 'gsl_root_fdfsolver' struct. The updates require both the function and its derivative (hence the name 'fdf') to be supplied by the user.  File: gsl-ref.info, Node: Root Finding Caveats, Next: Initializing the Solver, Prev: Root Finding Overview, Up: One dimensional Root-Finding 33.2 Caveats ============ Note that root finding functions can only search for one root at a time. When there are several roots in the search area, the first root to be found will be returned; however it is difficult to predict which of the roots this will be. _In most cases, no error will be reported if you try to find a root in an area where there is more than one._ Care must be taken when a function may have a multiple root (such as f(x) = (x-x_0)^2 or f(x) = (x-x_0)^3). It is not possible to use root-bracketing algorithms on even-multiplicity roots. For these algorithms the initial interval must contain a zero-crossing, where the function is negative at one end of the interval and positive at the other end. Roots with even-multiplicity do not cross zero, but only touch it instantaneously. Algorithms based on root bracketing will still work for odd-multiplicity roots (e.g. cubic, quintic, ...). Root polishing algorithms generally work with higher multiplicity roots, but at a reduced rate of convergence. In these cases the "Steffenson algorithm" can be used to accelerate the convergence of multiple roots. While it is not absolutely required that f have a root within the search region, numerical root finding functions should not be used haphazardly to check for the _existence_ of roots. There are better ways to do this. Because it is easy to create situations where numerical root finders can fail, it is a bad idea to throw a root finder at a function you do not know much about. In general it is best to examine the function visually by plotting before searching for a root.  File: gsl-ref.info, Node: Initializing the Solver, Next: Providing the function to solve, Prev: Root Finding Caveats, Up: One dimensional Root-Finding 33.3 Initializing the Solver ============================ -- Function: gsl_root_fsolver * gsl_root_fsolver_alloc (const gsl_root_fsolver_type * T) This function returns a pointer to a newly allocated instance of a solver of type T. For example, the following code creates an instance of a bisection solver, const gsl_root_fsolver_type * T = gsl_root_fsolver_bisection; gsl_root_fsolver * s = gsl_root_fsolver_alloc (T); If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of 'GSL_ENOMEM'. -- Function: gsl_root_fdfsolver * gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * T) This function returns a pointer to a newly allocated instance of a derivative-based solver of type T. For example, the following code creates an instance of a Newton-Raphson solver, const gsl_root_fdfsolver_type * T = gsl_root_fdfsolver_newton; gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc (T); If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of 'GSL_ENOMEM'. -- Function: int gsl_root_fsolver_set (gsl_root_fsolver * S, gsl_function * F, double X_LOWER, double X_UPPER) This function initializes, or reinitializes, an existing solver S to use the function F and the initial search interval [X_LOWER, X_UPPER]. -- Function: int gsl_root_fdfsolver_set (gsl_root_fdfsolver * S, gsl_function_fdf * FDF, double ROOT) This function initializes, or reinitializes, an existing solver S to use the function and derivative FDF and the initial guess ROOT. -- Function: void gsl_root_fsolver_free (gsl_root_fsolver * S) -- Function: void gsl_root_fdfsolver_free (gsl_root_fdfsolver * S) These functions free all the memory associated with the solver S. -- Function: const char * gsl_root_fsolver_name (const gsl_root_fsolver * S) -- Function: const char * gsl_root_fdfsolver_name (const gsl_root_fdfsolver * S) These functions return a pointer to the name of the solver. For example, printf ("s is a '%s' solver\n", gsl_root_fsolver_name (s)); would print something like 's is a 'bisection' solver'.  File: gsl-ref.info, Node: Providing the function to solve, Next: Search Bounds and Guesses, Prev: Initializing the Solver, Up: One dimensional Root-Finding 33.4 Providing the function to solve ==================================== You must provide a continuous function of one variable for the root finders to operate on, and, sometimes, its first derivative. In order to allow for general parameters the functions are defined by the following data types: -- Data Type: gsl_function This data type defines a general function with parameters. 'double (* function) (double X, void * PARAMS)' this function should return the value f(x,params) for argument X and parameters PARAMS 'void * params' a pointer to the parameters of the function Here is an example for the general quadratic function, f(x) = a x^2 + b x + c with a = 3, b = 2, c = 1. The following code defines a 'gsl_function' 'F' which you could pass to a root finder as a function pointer: struct my_f_params { double a; double b; double c; }; double my_f (double x, void * p) { struct my_f_params * params = (struct my_f_params *)p; double a = (params->a); double b = (params->b); double c = (params->c); return (a * x + b) * x + c; } gsl_function F; struct my_f_params params = { 3.0, 2.0, 1.0 }; F.function = &my_f; F.params = ¶ms; The function f(x) can be evaluated using the macro 'GSL_FN_EVAL(&F,x)' defined in 'gsl_math.h'. -- Data Type: gsl_function_fdf This data type defines a general function with parameters and its first derivative. 'double (* f) (double X, void * PARAMS)' this function should return the value of f(x,params) for argument X and parameters PARAMS 'double (* df) (double X, void * PARAMS)' this function should return the value of the derivative of F with respect to X, f'(x,params), for argument X and parameters PARAMS 'void (* fdf) (double X, void * PARAMS, double * F, double * DF)' this function should set the values of the function F to f(x,params) and its derivative DF to f'(x,params) for argument X and parameters PARAMS. This function provides an optimization of the separate functions for f(x) and f'(x)--it is always faster to compute the function and its derivative at the same time. 'void * params' a pointer to the parameters of the function Here is an example where f(x) = 2\exp(2x): double my_f (double x, void * params) { return exp (2 * x); } double my_df (double x, void * params) { return 2 * exp (2 * x); } void my_fdf (double x, void * params, double * f, double * df) { double t = exp (2 * x); *f = t; *df = 2 * t; /* uses existing value */ } gsl_function_fdf FDF; FDF.f = &my_f; FDF.df = &my_df; FDF.fdf = &my_fdf; FDF.params = 0; The function f(x) can be evaluated using the macro 'GSL_FN_FDF_EVAL_F(&FDF,x)' and the derivative f'(x) can be evaluated using the macro 'GSL_FN_FDF_EVAL_DF(&FDF,x)'. Both the function y = f(x) and its derivative dy = f'(x) can be evaluated at the same time using the macro 'GSL_FN_FDF_EVAL_F_DF(&FDF,x,y,dy)'. The macro stores f(x) in its Y argument and f'(x) in its DY argument--both of these should be pointers to 'double'.  File: gsl-ref.info, Node: Search Bounds and Guesses, Next: Root Finding Iteration, Prev: Providing the function to solve, Up: One dimensional Root-Finding 33.5 Search Bounds and Guesses ============================== You provide either search bounds or an initial guess; this section explains how search bounds and guesses work and how function arguments control them. A guess is simply an x value which is iterated until it is within the desired precision of a root. It takes the form of a 'double'. Search bounds are the endpoints of an interval which is iterated until the length of the interval is smaller than the requested precision. The interval is defined by two values, the lower limit and the upper limit. Whether the endpoints are intended to be included in the interval or not depends on the context in which the interval is used.  File: gsl-ref.info, Node: Root Finding Iteration, Next: Search Stopping Parameters, Prev: Search Bounds and Guesses, Up: One dimensional Root-Finding 33.6 Iteration ============== The following functions drive the iteration of each algorithm. Each function performs one iteration to update the state of any solver of the corresponding type. The same functions work for all solvers so that different methods can be substituted at runtime without modifications to the code. -- Function: int gsl_root_fsolver_iterate (gsl_root_fsolver * S) -- Function: int gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * S) These functions perform a single iteration of the solver S. If the iteration encounters an unexpected problem then an error code will be returned, 'GSL_EBADFUNC' the iteration encountered a singular point where the function or its derivative evaluated to 'Inf' or 'NaN'. 'GSL_EZERODIV' the derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero. The solver maintains a current best estimate of the root at all times. The bracketing solvers also keep track of the current best interval bounding the root. This information can be accessed with the following auxiliary functions, -- Function: double gsl_root_fsolver_root (const gsl_root_fsolver * S) -- Function: double gsl_root_fdfsolver_root (const gsl_root_fdfsolver * S) These functions return the current estimate of the root for the solver S. -- Function: double gsl_root_fsolver_x_lower (const gsl_root_fsolver * S) -- Function: double gsl_root_fsolver_x_upper (const gsl_root_fsolver * S) These functions return the current bracketing interval for the solver S.  File: gsl-ref.info, Node: Search Stopping Parameters, Next: Root Bracketing Algorithms, Prev: Root Finding Iteration, Up: One dimensional Root-Finding 33.7 Search Stopping Parameters =============================== A root finding procedure should stop when one of the following conditions is true: * A root has been found to within the user-specified precision. * A user-specified maximum number of iterations has been reached. * An error has occurred. The handling of these conditions is under user control. The functions below allow the user to test the precision of the current result in several standard ways. -- Function: int gsl_root_test_interval (double X_LOWER, double X_UPPER, double EPSABS, double EPSREL) This function tests for the convergence of the interval [X_LOWER, X_UPPER] with absolute error EPSABS and relative error EPSREL. The test returns 'GSL_SUCCESS' if the following condition is achieved, |a - b| < epsabs + epsrel min(|a|,|b|) when the interval x = [a,b] does not include the origin. If the interval includes the origin then \min(|a|,|b|) is replaced by zero (which is the minimum value of |x| over the interval). This ensures that the relative error is accurately estimated for roots close to the origin. This condition on the interval also implies that any estimate of the root r in the interval satisfies the same condition with respect to the true root r^*, |r - r^*| < epsabs + epsrel r^* assuming that the true root r^* is contained within the interval. -- Function: int gsl_root_test_delta (double X1, double X0, double EPSABS, double EPSREL) This function tests for the convergence of the sequence ..., X0, X1 with absolute error EPSABS and relative error EPSREL. The test returns 'GSL_SUCCESS' if the following condition is achieved, |x_1 - x_0| < epsabs + epsrel |x_1| and returns 'GSL_CONTINUE' otherwise. -- Function: int gsl_root_test_residual (double F, double EPSABS) This function tests the residual value F against the absolute error bound EPSABS. The test returns 'GSL_SUCCESS' if the following condition is achieved, |f| < epsabs and returns 'GSL_CONTINUE' otherwise. This criterion is suitable for situations where the precise location of the root, x, is unimportant provided a value can be found where the residual, |f(x)|, is small enough.  File: gsl-ref.info, Node: Root Bracketing Algorithms, Next: Root Finding Algorithms using Derivatives, Prev: Search Stopping Parameters, Up: One dimensional Root-Finding 33.8 Root Bracketing Algorithms =============================== The root bracketing algorithms described in this section require an initial interval which is guaranteed to contain a root--if a and b are the endpoints of the interval then f(a) must differ in sign from f(b). This ensures that the function crosses zero at least once in the interval. If a valid initial interval is used then these algorithm cannot fail, provided the function is well-behaved. Note that a bracketing algorithm cannot find roots of even degree, since these do not cross the x-axis. -- Solver: gsl_root_fsolver_bisection The "bisection algorithm" is the simplest method of bracketing the roots of a function. It is the slowest algorithm provided by the library, with linear convergence. On each iteration, the interval is bisected and the value of the function at the midpoint is calculated. The sign of this value is used to determine which half of the interval does not contain a root. That half is discarded to give a new, smaller interval containing the root. This procedure can be continued indefinitely until the interval is sufficiently small. At any time the current estimate of the root is taken as the midpoint of the interval. -- Solver: gsl_root_fsolver_falsepos The "false position algorithm" is a method of finding roots based on linear interpolation. Its convergence is linear, but it is usually faster than bisection. On each iteration a line is drawn between the endpoints (a,f(a)) and (b,f(b)) and the point where this line crosses the x-axis taken as a "midpoint". The value of the function at this point is calculated and its sign is used to determine which side of the interval does not contain a root. That side is discarded to give a new, smaller interval containing the root. This procedure can be continued indefinitely until the interval is sufficiently small. The best estimate of the root is taken from the linear interpolation of the interval on the current iteration. -- Solver: gsl_root_fsolver_brent The "Brent-Dekker method" (referred to here as "Brent's method") combines an interpolation strategy with the bisection algorithm. This produces a fast algorithm which is still robust. On each iteration Brent's method approximates the function using an interpolating curve. On the first iteration this is a linear interpolation of the two endpoints. For subsequent iterations the algorithm uses an inverse quadratic fit to the last three points, for higher accuracy. The intercept of the interpolating curve with the x-axis is taken as a guess for the root. If it lies within the bounds of the current interval then the interpolating point is accepted, and used to generate a smaller interval. If the interpolating point is not accepted then the algorithm falls back to an ordinary bisection step. The best estimate of the root is taken from the most recent interpolation or bisection.  File: gsl-ref.info, Node: Root Finding Algorithms using Derivatives, Next: Root Finding Examples, Prev: Root Bracketing Algorithms, Up: One dimensional Root-Finding 33.9 Root Finding Algorithms using Derivatives ============================================== The root polishing algorithms described in this section require an initial guess for the location of the root. There is no absolute guarantee of convergence--the function must be suitable for this technique and the initial guess must be sufficiently close to the root for it to work. When these conditions are satisfied then convergence is quadratic. These algorithms make use of both the function and its derivative. -- Derivative Solver: gsl_root_fdfsolver_newton Newton's Method is the standard root-polishing algorithm. The algorithm begins with an initial guess for the location of the root. On each iteration, a line tangent to the function f is drawn at that position. The point where this line crosses the x-axis becomes the new guess. The iteration is defined by the following sequence, x_{i+1} = x_i - f(x_i)/f'(x_i) Newton's method converges quadratically for single roots, and linearly for multiple roots. -- Derivative Solver: gsl_root_fdfsolver_secant The "secant method" is a simplified version of Newton's method which does not require the computation of the derivative on every step. On its first iteration the algorithm begins with Newton's method, using the derivative to compute a first step, x_1 = x_0 - f(x_0)/f'(x_0) Subsequent iterations avoid the evaluation of the derivative by replacing it with a numerical estimate, the slope of the line through the previous two points, x_{i+1} = x_i f(x_i) / f'_{est} where f'_{est} = (f(x_i) - f(x_{i-1})/(x_i - x_{i-1}) When the derivative does not change significantly in the vicinity of the root the secant method gives a useful saving. Asymptotically the secant method is faster than Newton's method whenever the cost of evaluating the derivative is more than 0.44 times the cost of evaluating the function itself. As with all methods of computing a numerical derivative the estimate can suffer from cancellation errors if the separation of the points becomes too small. On single roots, the method has a convergence of order (1 + \sqrt 5)/2 (approximately 1.62). It converges linearly for multiple roots. -- Derivative Solver: gsl_root_fdfsolver_steffenson The "Steffenson Method"(1) provides the fastest convergence of all the routines. It combines the basic Newton algorithm with an Aitken "delta-squared" acceleration. If the Newton iterates are x_i then the acceleration procedure generates a new sequence R_i, R_i = x_i - (x_{i+1} - x_i)^2 / (x_{i+2} - 2 x_{i+1} + x_{i}) which converges faster than the original sequence under reasonable conditions. The new sequence requires three terms before it can produce its first value so the method returns accelerated values on the second and subsequent iterations. On the first iteration it returns the ordinary Newton estimate. The Newton iterate is also returned if the denominator of the acceleration term ever becomes zero. As with all acceleration procedures this method can become unstable if the function is not well-behaved. ---------- Footnotes ---------- (1) J.F. Steffensen (1873-1961). The spelling used in the name of the function is slightly incorrect, but has been preserved to avoid incompatibility.  File: gsl-ref.info, Node: Root Finding Examples, Next: Root Finding References and Further Reading, Prev: Root Finding Algorithms using Derivatives, Up: One dimensional Root-Finding 33.10 Examples ============== For any root finding algorithm we need to prepare the function to be solved. For this example we will use the general quadratic equation described earlier. We first need a header file ('demo_fn.h') to define the function parameters, struct quadratic_params { double a, b, c; }; double quadratic (double x, void *params); double quadratic_deriv (double x, void *params); void quadratic_fdf (double x, void *params, double *y, double *dy); We place the function definitions in a separate file ('demo_fn.c'), double quadratic (double x, void *params) { struct quadratic_params *p = (struct quadratic_params *) params; double a = p->a; double b = p->b; double c = p->c; return (a * x + b) * x + c; } double quadratic_deriv (double x, void *params) { struct quadratic_params *p = (struct quadratic_params *) params; double a = p->a; double b = p->b; double c = p->c; return 2.0 * a * x + b; } void quadratic_fdf (double x, void *params, double *y, double *dy) { struct quadratic_params *p = (struct quadratic_params *) params; double a = p->a; double b = p->b; double c = p->c; *y = (a * x + b) * x + c; *dy = 2.0 * a * x + b; } The first program uses the function solver 'gsl_root_fsolver_brent' for Brent's method and the general quadratic defined above to solve the following equation, x^2 - 5 = 0 with solution x = \sqrt 5 = 2.236068... #include #include #include #include #include "demo_fn.h" #include "demo_fn.c" int main (void) { int status; int iter = 0, max_iter = 100; const gsl_root_fsolver_type *T; gsl_root_fsolver *s; double r = 0, r_expected = sqrt (5.0); double x_lo = 0.0, x_hi = 5.0; gsl_function F; struct quadratic_params params = {1.0, 0.0, -5.0}; F.function = &quadratic; F.params = ¶ms; T = gsl_root_fsolver_brent; s = gsl_root_fsolver_alloc (T); gsl_root_fsolver_set (s, &F, x_lo, x_hi); printf ("using %s method\n", gsl_root_fsolver_name (s)); printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "root", "err", "err(est)"); do { iter++; status = gsl_root_fsolver_iterate (s); r = gsl_root_fsolver_root (s); x_lo = gsl_root_fsolver_x_lower (s); x_hi = gsl_root_fsolver_x_upper (s); status = gsl_root_test_interval (x_lo, x_hi, 0, 0.001); if (status == GSL_SUCCESS) printf ("Converged:\n"); printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, x_lo, x_hi, r, r - r_expected, x_hi - x_lo); } while (status == GSL_CONTINUE && iter < max_iter); gsl_root_fsolver_free (s); return status; } Here are the results of the iterations, $ ./a.out using brent method iter [ lower, upper] root err err(est) 1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000 2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000 3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000 4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000 5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300 Converged: 6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666 If the program is modified to use the bisection solver instead of Brent's method, by changing 'gsl_root_fsolver_brent' to 'gsl_root_fsolver_bisection' the slower convergence of the Bisection method can be observed, $ ./a.out using bisection method iter [ lower, upper] root err err(est) 1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000 2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000 3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000 4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000 5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500 6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250 7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625 8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312 9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656 10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828 11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414 Converged: 12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207 The next program solves the same function using a derivative solver instead. #include #include #include #include #include "demo_fn.h" #include "demo_fn.c" int main (void) { int status; int iter = 0, max_iter = 100; const gsl_root_fdfsolver_type *T; gsl_root_fdfsolver *s; double x0, x = 5.0, r_expected = sqrt (5.0); gsl_function_fdf FDF; struct quadratic_params params = {1.0, 0.0, -5.0}; FDF.f = &quadratic; FDF.df = &quadratic_deriv; FDF.fdf = &quadratic_fdf; FDF.params = ¶ms; T = gsl_root_fdfsolver_newton; s = gsl_root_fdfsolver_alloc (T); gsl_root_fdfsolver_set (s, &FDF, x); printf ("using %s method\n", gsl_root_fdfsolver_name (s)); printf ("%-5s %10s %10s %10s\n", "iter", "root", "err", "err(est)"); do { iter++; status = gsl_root_fdfsolver_iterate (s); x0 = x; x = gsl_root_fdfsolver_root (s); status = gsl_root_test_delta (x, x0, 0, 1e-3); if (status == GSL_SUCCESS) printf ("Converged:\n"); printf ("%5d %10.7f %+10.7f %10.7f\n", iter, x, x - r_expected, x - x0); } while (status == GSL_CONTINUE && iter < max_iter); gsl_root_fdfsolver_free (s); return status; } Here are the results for Newton's method, $ ./a.out using newton method iter root err err(est) 1 3.0000000 +0.7639320 -2.0000000 2 2.3333333 +0.0972654 -0.6666667 3 2.2380952 +0.0020273 -0.0952381 Converged: 4 2.2360689 +0.0000009 -0.0020263 Note that the error can be estimated more accurately by taking the difference between the current iterate and next iterate rather than the previous iterate. The other derivative solvers can be investigated by changing 'gsl_root_fdfsolver_newton' to 'gsl_root_fdfsolver_secant' or 'gsl_root_fdfsolver_steffenson'.  File: gsl-ref.info, Node: Root Finding References and Further Reading, Prev: Root Finding Examples, Up: One dimensional Root-Finding 33.11 References and Further Reading ==================================== For information on the Brent-Dekker algorithm see the following two papers, R. P. Brent, "An algorithm with guaranteed convergence for finding a zero of a function", 'Computer Journal', 14 (1971) 422-425 J. C. P. Bus and T. J. Dekker, "Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function", 'ACM Transactions of Mathematical Software', Vol. 1 No. 4 (1975) 330-345  File: gsl-ref.info, Node: One dimensional Minimization, Next: Multidimensional Root-Finding, Prev: One dimensional Root-Finding, Up: Top 34 One dimensional Minimization ******************************* This chapter describes routines for finding minima of arbitrary one-dimensional functions. The library provides low level components for a variety of iterative minimizers and convergence tests. These can be combined by the user to achieve the desired solution, with full access to the intermediate steps of the algorithms. Each class of methods uses the same framework, so that you can switch between minimizers at runtime without needing to recompile your program. Each instance of a minimizer keeps track of its own state, allowing the minimizers to be used in multi-threaded programs. The header file 'gsl_min.h' contains prototypes for the minimization functions and related declarations. To use the minimization algorithms to find the maximum of a function simply invert its sign. * Menu: * Minimization Overview:: * Minimization Caveats:: * Initializing the Minimizer:: * Providing the function to minimize:: * Minimization Iteration:: * Minimization Stopping Parameters:: * Minimization Algorithms:: * Minimization Examples:: * Minimization References and Further Reading::  File: gsl-ref.info, Node: Minimization Overview, Next: Minimization Caveats, Up: One dimensional Minimization 34.1 Overview ============= The minimization algorithms begin with a bounded region known to contain a minimum. The region is described by a lower bound a and an upper bound b, with an estimate of the location of the minimum x. The value of the function at x must be less than the value of the function at the ends of the interval, f(a) > f(x) < f(b) This condition guarantees that a minimum is contained somewhere within the interval. On each iteration a new point x' is selected using one of the available algorithms. If the new point is a better estimate of the minimum, i.e. where f(x') < f(x), then the current estimate of the minimum x is updated. The new point also allows the size of the bounded interval to be reduced, by choosing the most compact set of points which satisfies the constraint f(a) > f(x) < f(b). The interval is reduced until it encloses the true minimum to a desired tolerance. This provides a best estimate of the location of the minimum and a rigorous error estimate. Several bracketing algorithms are available within a single framework. The user provides a high-level driver for the algorithm, and the library provides the individual functions necessary for each of the steps. There are three main phases of the iteration. The steps are, * initialize minimizer state, S, for algorithm T * update S using the iteration T * test S for convergence, and repeat iteration if necessary The state for the minimizers is held in a 'gsl_min_fminimizer' struct. The updating procedure uses only function evaluations (not derivatives).  File: gsl-ref.info, Node: Minimization Caveats, Next: Initializing the Minimizer, Prev: Minimization Overview, Up: One dimensional Minimization 34.2 Caveats ============ Note that minimization functions can only search for one minimum at a time. When there are several minima in the search area, the first minimum to be found will be returned; however it is difficult to predict which of the minima this will be. _In most cases, no error will be reported if you try to find a minimum in an area where there is more than one._ With all minimization algorithms it can be difficult to determine the location of the minimum to full numerical precision. The behavior of the function in the region of the minimum x^* can be approximated by a Taylor expansion, y = f(x^*) + (1/2) f''(x^*) (x - x^*)^2 and the second term of this expansion can be lost when added to the first term at finite precision. This magnifies the error in locating x^*, making it proportional to \sqrt \epsilon (where \epsilon is the relative accuracy of the floating point numbers). For functions with higher order minima, such as x^4, the magnification of the error is correspondingly worse. The best that can be achieved is to converge to the limit of numerical accuracy in the function values, rather than the location of the minimum itself.  File: gsl-ref.info, Node: Initializing the Minimizer, Next: Providing the function to minimize, Prev: Minimization Caveats, Up: One dimensional Minimization 34.3 Initializing the Minimizer =============================== -- Function: gsl_min_fminimizer * gsl_min_fminimizer_alloc (const gsl_min_fminimizer_type * T) This function returns a pointer to a newly allocated instance of a minimizer of type T. For example, the following code creates an instance of a golden section minimizer, const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection; gsl_min_fminimizer * s = gsl_min_fminimizer_alloc (T); If there is insufficient memory to create the minimizer then the function returns a null pointer and the error handler is invoked with an error code of 'GSL_ENOMEM'. -- Function: int gsl_min_fminimizer_set (gsl_min_fminimizer * S, gsl_function * F, double X_MINIMUM, double X_LOWER, double X_UPPER) This function sets, or resets, an existing minimizer S to use the function F and the initial search interval [X_LOWER, X_UPPER], with a guess for the location of the minimum X_MINIMUM. If the interval given does not contain a minimum, then the function returns an error code of 'GSL_EINVAL'. -- Function: int gsl_min_fminimizer_set_with_values (gsl_min_fminimizer * S, gsl_function * F, double X_MINIMUM, double F_MINIMUM, double X_LOWER, double F_LOWER, double X_UPPER, double F_UPPER) This function is equivalent to 'gsl_min_fminimizer_set' but uses the values F_MINIMUM, F_LOWER and F_UPPER instead of computing 'f(x_minimum)', 'f(x_lower)' and 'f(x_upper)'. -- Function: void gsl_min_fminimizer_free (gsl_min_fminimizer * S) This function frees all the memory associated with the minimizer S. -- Function: const char * gsl_min_fminimizer_name (const gsl_min_fminimizer * S) This function returns a pointer to the name of the minimizer. For example, printf ("s is a '%s' minimizer\n", gsl_min_fminimizer_name (s)); would print something like 's is a 'brent' minimizer'.  File: gsl-ref.info, Node: Providing the function to minimize, Next: Minimization Iteration, Prev: Initializing the Minimizer, Up: One dimensional Minimization 34.4 Providing the function to minimize ======================================= You must provide a continuous function of one variable for the minimizers to operate on. In order to allow for general parameters the functions are defined by a 'gsl_function' data type (*note Providing the function to solve::).  File: gsl-ref.info, Node: Minimization Iteration, Next: Minimization Stopping Parameters, Prev: Providing the function to minimize, Up: One dimensional Minimization 34.5 Iteration ============== The following functions drive the iteration of each algorithm. Each function performs one iteration to update the state of any minimizer of the corresponding type. The same functions work for all minimizers so that different methods can be substituted at runtime without modifications to the code. -- Function: int gsl_min_fminimizer_iterate (gsl_min_fminimizer * S) This function performs a single iteration of the minimizer S. If the iteration encounters an unexpected problem then an error code will be returned, 'GSL_EBADFUNC' the iteration encountered a singular point where the function evaluated to 'Inf' or 'NaN'. 'GSL_FAILURE' the algorithm could not improve the current best approximation or bounding interval. The minimizer maintains a current best estimate of the position of the minimum at all times, and the current interval bounding the minimum. This information can be accessed with the following auxiliary functions, -- Function: double gsl_min_fminimizer_x_minimum (const gsl_min_fminimizer * S) This function returns the current estimate of the position of the minimum for the minimizer S. -- Function: double gsl_min_fminimizer_x_upper (const gsl_min_fminimizer * S) -- Function: double gsl_min_fminimizer_x_lower (const gsl_min_fminimizer * S) These functions return the current upper and lower bound of the interval for the minimizer S. -- Function: double gsl_min_fminimizer_f_minimum (const gsl_min_fminimizer * S) -- Function: double gsl_min_fminimizer_f_upper (const gsl_min_fminimizer * S) -- Function: double gsl_min_fminimizer_f_lower (const gsl_min_fminimizer * S) These functions return the value of the function at the current estimate of the minimum and at the upper and lower bounds of the interval for the minimizer S.  File: gsl-ref.info, Node: Minimization Stopping Parameters, Next: Minimization Algorithms, Prev: Minimization Iteration, Up: One dimensional Minimization 34.6 Stopping Parameters ======================== A minimization procedure should stop when one of the following conditions is true: * A minimum has been found to within the user-specified precision. * A user-specified maximum number of iterations has been reached. * An error has occurred. The handling of these conditions is under user control. The function below allows the user to test the precision of the current result. -- Function: int gsl_min_test_interval (double X_LOWER, double X_UPPER, double EPSABS, double EPSREL) This function tests for the convergence of the interval [X_LOWER, X_UPPER] with absolute error EPSABS and relative error EPSREL. The test returns 'GSL_SUCCESS' if the following condition is achieved, |a - b| < epsabs + epsrel min(|a|,|b|) when the interval x = [a,b] does not include the origin. If the interval includes the origin then \min(|a|,|b|) is replaced by zero (which is the minimum value of |x| over the interval). This ensures that the relative error is accurately estimated for minima close to the origin. This condition on the interval also implies that any estimate of the minimum x_m in the interval satisfies the same condition with respect to the true minimum x_m^*, |x_m - x_m^*| < epsabs + epsrel x_m^* assuming that the true minimum x_m^* is contained within the interval.  File: gsl-ref.info, Node: Minimization Algorithms, Next: Minimization Examples, Prev: Minimization Stopping Parameters, Up: One dimensional Minimization 34.7 Minimization Algorithms ============================ The minimization algorithms described in this section require an initial interval which is guaranteed to contain a minimum--if a and b are the endpoints of the interval and x is an estimate of the minimum then f(a) > f(x) < f(b). This ensures that the function has at least one minimum somewhere in the interval. If a valid initial interval is used then these algorithm cannot fail, provided the function is well-behaved. -- Minimizer: gsl_min_fminimizer_goldensection The "golden section algorithm" is the simplest method of bracketing the minimum of a function. It is the slowest algorithm provided by the library, with linear convergence. On each iteration, the algorithm first compares the subintervals from the endpoints to the current minimum. The larger subinterval is divided in a golden section (using the famous ratio (3-\sqrt 5)/2 = 0.3189660...) and the value of the function at this new point is calculated. The new value is used with the constraint f(a') > f(x') < f(b') to a select new interval containing the minimum, by discarding the least useful point. This procedure can be continued indefinitely until the interval is sufficiently small. Choosing the golden section as the bisection ratio can be shown to provide the fastest convergence for this type of algorithm. -- Minimizer: gsl_min_fminimizer_brent The "Brent minimization algorithm" combines a parabolic interpolation with the golden section algorithm. This produces a fast algorithm which is still robust. The outline of the algorithm can be summarized as follows: on each iteration Brent's method approximates the function using an interpolating parabola through three existing points. The minimum of the parabola is taken as a guess for the minimum. If it lies within the bounds of the current interval then the interpolating point is accepted, and used to generate a smaller interval. If the interpolating point is not accepted then the algorithm falls back to an ordinary golden section step. The full details of Brent's method include some additional checks to improve convergence. -- Minimizer: gsl_min_fminimizer_quad_golden This is a variant of Brent's algorithm which uses the safeguarded step-length algorithm of Gill and Murray.  File: gsl-ref.info, Node: Minimization Examples, Next: Minimization References and Further Reading, Prev: Minimization Algorithms, Up: One dimensional Minimization 34.8 Examples ============= The following program uses the Brent algorithm to find the minimum of the function f(x) = \cos(x) + 1, which occurs at x = \pi. The starting interval is (0,6), with an initial guess for the minimum of 2. #include #include #include #include double fn1 (double x, void * params) { return cos(x) + 1.0; } int main (void) { int status; int iter = 0, max_iter = 100; const gsl_min_fminimizer_type *T; gsl_min_fminimizer *s; double m = 2.0, m_expected = M_PI; double a = 0.0, b = 6.0; gsl_function F; F.function = &fn1; F.params = 0; T = gsl_min_fminimizer_brent; s = gsl_min_fminimizer_alloc (T); gsl_min_fminimizer_set (s, &F, m, a, b); printf ("using %s method\n", gsl_min_fminimizer_name (s)); printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "min", "err", "err(est)"); printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, a, b, m, m - m_expected, b - a); do { iter++; status = gsl_min_fminimizer_iterate (s); m = gsl_min_fminimizer_x_minimum (s); a = gsl_min_fminimizer_x_lower (s); b = gsl_min_fminimizer_x_upper (s); status = gsl_min_test_interval (a, b, 0.001, 0.0); if (status == GSL_SUCCESS) printf ("Converged:\n"); printf ("%5d [%.7f, %.7f] " "%.7f %+.7f %.7f\n", iter, a, b, m, m - m_expected, b - a); } while (status == GSL_CONTINUE && iter < max_iter); gsl_min_fminimizer_free (s); return status; } Here are the results of the minimization procedure. $ ./a.out 0 [0.0000000, 6.0000000] 2.0000000 -1.1415927 6.0000000 1 [2.0000000, 6.0000000] 3.2758640 +0.1342713 4.0000000 2 [2.0000000, 3.2831929] 3.2758640 +0.1342713 1.2831929 3 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862 4 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862 5 [2.8689068, 3.2758640] 3.1460585 +0.0044658 0.4069572 6 [3.1346075, 3.2758640] 3.1460585 +0.0044658 0.1412565 7 [3.1346075, 3.1874620] 3.1460585 +0.0044658 0.0528545 8 [3.1346075, 3.1460585] 3.1460585 +0.0044658 0.0114510 9 [3.1346075, 3.1460585] 3.1424060 +0.0008133 0.0114510 10 [3.1346075, 3.1424060] 3.1415885 -0.0000041 0.0077985 Converged: 11 [3.1415885, 3.1424060] 3.1415927 -0.0000000 0.0008175